In this work we discuss the emergence of approximate causality in a general setup from waveguide QEDi.e. a one-dimensional propagating field interacting with a scatterer. We prove that this emergent causality translates into a structure for the N-photon scattering matrix. Our work builds on the derivation of a Lieb-Robinson-type bound for continuous models and for all coupling strengths, as well as on several intermediate results, of which we highlight: (i) the asymptotic independence of space-like separated wave packets, (ii) the proper definition of input and output scattering states, and (iii) the characterization of the ground state and correlations in the model. We illustrate our formal results by analyzing the two-photon scattering from a quantum impurity in the ultrastrong coupling regime, verifying the cluster decomposition and ground-state nature. Besides, we generalize the cluster decomposition if inelastic or Raman scattering occurs, finding the structure of the S-matrixin momentum space for linear dispersion relations. In this case, we compute the decay of the fluorescence (photon-photon correlations) caused by this S-matrix. models do not satisfy Lorentz or translational invariance, they are typically dispersive, and the photon-matter interaction may become highly non-perturbative. Experimental implementations include dielectrics [14,15], cavity arrays [16], metals [17], diamond structures [18,19], and superconductors [20][21][22][23][24] interacting with atoms, molecules, quantum dots, color centers in diamond or superconducting qubits. The focus of waveguide QED is set on quantum processes involving few photons and scatterers. In this regard, it is not surprising that there exists an extensive theoretical literature for waveguide-QED systems [13], which develops a variety of analytical and numerical methods for the study of the N-photon S-matrix [25][26][27][28][29][30][31][32][33][34][35][36].The main result in this work is the structure of the N-photon S-matrixin waveguide QED, rigorously deduced from emergent causality constraints. Our result builds on a general model of light-matter interactions, without any approximations such as the rotating-wave approximation (RWA), the Markovian limit, or weak light-matter coupling. To derive the S-matrixdecomposition we are assisted by several intermediate and important results, of which we remark: (i) the freedom of wave packets far away from the scatterer, (ii) Lieb-Robinson-like independence relations and approximate light-cones for propagating wave packets, (iii) a characterization of the ground state correlation properties, and (iv) a proper definition and derivation of scattering input and output states.We illustrate our results with two representative examples. The first one is a numerical study of scattering in the ultrastrong coupling limit [35,37], where we demonstrate the clustering decomposition and the nature of the ground state predicted by our intermediate results. The second is an analytical study of a non-dispersive medium interacting ...