Emergent Causality and the N-photon Scattering Matrix in Waveguide QED

In this work we discuss the emergence of approximate causality in a general setup from waveguide QED -i.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-matrix in momentum space for linear dispersion relations. In this case, we compute the decay of the fluorescence (photon-photon correlations) caused by this S-matrix.


I. INTRODUCTION
Causality is expected to hold in every circumstance. The causality principle states that two experiments which are space-like separated, such that no signal travelling at the speed of light can connect them, must provide uncorrelated results [1]. In Quantum Field Theory (QFT), strict causality imposes that two operators A(x, t) and B(y, t ) acting on two space-like separated points (x, t) and (y, t ), must commute, where c is the speed of light (we restrict ourselves to 1 + 1 dimensions). Another consequence of causality in QFT appears in the study of scattering events or collisions: scattering matrices describing causally disconnected events must "cluster", or decompose into a product of independent scattering matrices [2]. In fact, all acceptable QFT interactions must result in S-matrices fulfilling such a decomposition [3]. Nonrelativistic quantum mechanics is an effective theory which allows signals to propagate arbitrarily fast, but which may give rise to different forms of emergent approximate causality. The typical examples are lowenergy models in solid state, where quasiparticle excitations have a maximum group velocity. In this case, there exists an approximate light cone, outside of which the correlations between operators are exponentially suppressed. This emergent causality was rigorously demonstrated by Lieb and Robinson [4] for spin-models on lattices with bounded interactions that decay rapidly with the distance. Lieb-Robinson bounds not only imply causality in the information-theoretical sense [5], but lead to important results in the static properties of manybody Hamiltonians, such as the clustering of correlations and the area law in gapped models [6,7].
The main result in this work is the structure of the Nphoton S-matrix in 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 (RWA), the Markovian limit, or weak light-matter coupling. To derive the S-matrix decomposition 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 ex-amples. The first one is a numerical study of scattering in the ultrastrong coupling limit [33,35], 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 with a general scatterer, which admits exact calculations. Here, we find the shape of the S-matrix from general principles, including the inelastic processes. We recover the nontrivial form computed by Xu and Fan for a particular case in [36] and find the natural generalization of the standard cluster decomposition. The paper has the following organization. Sect. II presents the nonrelativistic Hamiltonian that models the interaction between propagating photons and quantum impurities, the concept of wave packet, a review of the scattering theory needed, and two conditions necessary for the validity of our results. Sect. III summarizes our formal theory arriving to the general N -photon scattering compatible with causality. Sect. IV presents the examples applying the theory. We close this work with further comments and outlooks. Intermediate lemmas, theorems, and technical issues are discussed in the appendices.

II. MODEL AND SCATTERING THEORY
A. Waveguide QED model The simplest model that describes a waveguide-QED setup consists of a one-dimensional bosonic medium and a scatterer. Using units such that = 1, it reads The first term stands for the free-Hamiltonian of the photons with frequency ω k for momentum k, which is created (annihilated) by the corresponding Fock operator a k (a † k ), satisfying [a k , a † k ] = δ(k − k ). The last two terms are the Hamiltonian H sc of the finite-dimensional system, which is the scatterer, and the dipolar interaction term described by the bounded operators G and the coupling strengths g k . We assume that the coupling strengths in position space have a finite support centered around x sc = 0. The model (2) is not exactly solvable in general. For instance, if the scatterer is a two-level system, H sc ∝ σ z and G = σ x the model is the celebrated spin-boson model [37], which results in a nontrivial ground state with localized photonic excitations around the scatterer.  Figure 1. Two incoming photons with average momentak1 (red) andk2 (green), initially centered around distant points x1 andx2 (l → ∞), scatter against a general quantum object. The scatterer-field can have several bound states (localized and not propagating). In the figure, the scatterer-field is in one of those bound states |Ων (gray region). If the first incoming photon leaves the scattering region in another localized eigenstate |Ω λ the second photon meets the interaction region in a different state that the found by the first wave packet. If this occurs [see main text] the scattering matrix cannot be just a product, it must differentiate the order in which both events happen.
The discussion below assumes a single photonic band ω k ∈ [ω min , ω max ] and typically a chiral medium k ≥ 0, ∂ k ω k ≥ 0. This is a rather standard simplification which does not affect the generality and applicability of our results. More generic dispersion relations and nonchiral medium can be taken into account by introducing additional degrees of freedom in the photons (chirality, band index, etc.) and keeping track of those quantum numbers in a trivial extension of our results.

B. Localized wave packets
In order to talk about causality, we introduce a set of localized wave packets to which an approximate position can be ascribed. As we will see below, approximate localization becomes essential in the discussion, allowing us to discuss the order in which photons interact with the scatterer.
Let us introduce the creation operator ψkx(t) † for a wave packet as The wavefunction φk(p) = φ(p −k) ∈ L 2 is normalized and centered around the average momentumk. The exponential factor e ikx ensures the wave packet is centered aroundx in position space at time t = t 0 .
As wave packets we will use both Gaussian and Lorentzian envelopes These wave functions are only approximately localized in the sense that the probability of finding a photon decays exponentially far away from the centerx. The width σ in momentum space implies a localization length 1/σ in position space. Note that our definition of the wave packet lacks factors such as √ ω k or ω 1/2 k which typically appear when transforming back to position space from a linear bosonic problem that was diagonalized in frequency space. This is a convenient definition that avoids divergences when computing things such as the number of photons. The choice of prefactors is ultimately irrelevant when we take the limit σ → 0 in many of the argumentations below. Fig. 1 illustrates the collision of two approximately localized wave packets against a quantum impurity in a chiral medium. The average momentum of the wave packets k 1 ork 2 determines the group velocity at which the photons move v g (k) = ∂ k ω k . The wave packets may be distorted due both to the dispersive nature of the medium and the interaction with the scatterer.

C. Scattering operator
In the typical scattering geometry, the interaction occurs in a finite region. Besides, it is assumed that asymptotically far away from that region the field is a linear combination of free-particle states (generated via creation operators on the non-interacting vacuum) even in the presence of the scatterer-waveguide interaction.
A sufficient condition for this is that both the ground state and any non-propagating excited state accessible by scattering |Ω µ are indistinguishable from the vacuum state |vac far away from the scatterer. Mathematically this occurs when where O(x, ∆) is an operator with compact support in the finite intervalx−∆/2 < x <x+∆/2 and the vacuum state |vac is such that a k |vac = 0 ∀k. Besides, the free particle states must satisfy the asymptotic condition [38]: with U (t) the evolution operator of the full Hamiltonian (2) and U 0 (t) = e −iH0t the free-evolution operator.
The scattering operator S relates the amplitude of the output and input fields through which, using (9), has the formal expression: Here, U I (t + , t − ) = e iH0t− e −iH(t+−t−) e −iH0t+ is the evolution operator in the interaction picture. Using again Eq. (9) leads to |Ψ in/out = U † 0 (t −/+ )U (t −/+ )|Ψ ≡ |Ψ(t −/+ ) I , which shows that the input and output fields are represented in the interaction picture.
Related quantities are the scattering amplitudes. For example, the single-photon amplitude is defined as: with ψ in (t − ) † = ψkx(t − ) † and an analogous definition for ψ out (t + ) and the photon mean positionx being well separated from the scatterer. One of the goals of this work is to find the most general form for the amplitude A compatible with causality, thus providing a more clear understanding of the structure of the scattering matrix.

D. Sufficient conditions for having a well-defined scattering theory
Given a general Hamiltonian (2), it is not generally known whether the condition (8) is satisfied. Thus, the existence of scattering states must be assumed. In this work, we provide a further evidence of the validity of this assumptions by demonstrating a limited version of Eq. (8) (see App. A) for the unique ground state of Hamiltonian (2), which reads provided that (i) for all k, |g k /ω k | < ∞ and (ii) that the correlators C kp = Ω 0 |a † k a p |Ω 0 are n-differentiable functions.
Unfortunately, this result is insufficient for treating the most general case. It is well known that the Hamiltonian (2) may support excited eigenstates which are localized around the scattering center [32,33,39,40], which in the literature are usually referred as ground states. Two paradigmatic examples of scatterer with multiple ground states are the three-level Λ atom, with two electronic ground state, and a two-level system coupled to a cavity array in the ultrastrong coupling regime [33].
However, we have been unable to find a general proof that (8) is satisfied (and thus that input and output states can be defined) for non propagating excited states that appear in these systems. In order to make any progress, and as usual in the literature, we have instead assumed a plausible first condition: the Hamiltonian (2) has a finite set of ground states, {|Ω µ }, which are localized in the sense of Eq. (8). Notice that with this assumption (2) has a well defined theory [See. App. B 3]. This condition allows the expression of the elements of S in the momentum basis: In this paper we will also assume a second condition: the N-photon scattering process conserves the number of flying photons in the input and output states. We only provide results for the sector of the scattering matrix that conserves the number of excitations, excluding us from considering other scattering channels, such as downconversion processes. Notice, however, that a large number of systems fulfill this condition. For instance, the unbiased spin-boson model (where H sc ∝ σ z and G = σ x ) exactly conserves the number of excitations within the rotating-wave-approximation, which is valid when the coupling strength is much smaller than the photon energy. But even in the ultrastrong coupling regime, when counter-rotating terms are important, numerical simulations have shown that the scattering process conserves the number of flying excitations within numerical uncertainties [cf. Refs. [33,35,41] and Sect. IV A].

III. CAUSALITY AND THE N -PHOTON SCATTERING MATRIX
A. Approximate causality We are describing waveguide QED using nonrelativistic models for which strict causality (1) does not apply. However, as a foundational result we have been able to prove that the waveguide-QED model (2) supports an approximate form of causality. This form states that there exists an approximate light cone, defined by the maximum group velocity, c = max(∂ k ω k ). Two wave-packet operators which are outside their respective cones and far away from the scatterer approximately commute.
To be precise, we define the distance d( , d(ȳ, t 0 )} the distance between the packets and the minimum distance between them and the scatterer respectively. The power n stands because we use that the dispersion relation is n-times differentiable. A sketch of the proof is as follows. First, we prove (15) for free fields, i.e. for wave packets moving under H 0 . In the Heisenberg picture, the phases ik(x −ȳ) − iω k (t − t ) can be bounded by the distance d(x − y, t − t ). Using the Riemann-Lebesgue lemma ( e ikz f (k) dk → 0, as z → ∞) we find the power law decay, |D| −n . Causality is thereby linked to the cancellation or averaging of fast oscillations in the unitary dynamics. Applying a similar technique to the interaction term in (2) allows us to prove that packets away the influence of the scatterer evolve freely, producing the second algebraic decay term |D 0 | 1−n . This leads the second decay |D 0 | 1−n . If their evolution can be approximated by the evolution under H 0 , what we found for the commutator of free-evolving packets holds also in the interacting part.
This result is analogous to Lieb-Robinson-type bounds that were initially developed for a lattice of locally interacting spins [4], and which were later generalized to finite-dimensional models, anharmonic oscillators, master equations, and spin-boson lattices [6,[42][43][44][45][46]. It is important to remark that the approximate causality in Eq. (15) is not obtained for the free theory, but for the full waveguide-QED model. As a consequence, it can be used to derive important results on the photon-scatterer interaction.

B. Causality and the scattering matrix
Causality imposes restrictions on the S-matrix [3], among which is the cluster decomposition that we summarize here. For now, let us consider the case of a unique ground state and split the S-matrix into a free part S 0 and an interacting part T , both in momentum space The interacting part T accounts for processes in which two or more photons coincide and interact simultaneously with the scatterer. Causality is then invoked to argue that they cannot influence each other if the input events are space-like separated. Thus, T does not contribute to the scattering amplitude as wave packets fall apart |x i −x j | → ∞. This, together with energy conservation, imposes the constraint iT pk = iC pk δ(E p − E k ) [1]. In this limit the only term contributing to the scattering amplitude is the free part, S 0 . In QFT (typically) occurs momentum conservation which implies that with S pnkn ∝ δ(ω pn − ω kn ) the one-photon S-matrix. This is nothing but the cluster decomposition. Fourier transforming S 0 pk , this structure also holds This shall be relevant in the following section, where we will work in position space.

C. Generalized cluster decomposition
Our goal is to explain how approximate causality (15) implies a cluster decomposition for the S-matrix. We will also show that in waveguide QED the photon momenta need not be conserved and that S 0 may not have the structure given by Eq. (17).
To understand how causality fixes the form of S 0 we refer to our Fig. 1 where two well separated wave packets interact with a scatterer. The scattering amplitude is, Note that for a sufficiently large separation of the wave packets, the output state of the first packet must be causally disconnected. This implies that the input operator for the first wave packet must commute with the output operator for the second packet [see Eq. (15)]. Notice that the second output and the first input will not commute in general. We can then approximate, at any degree of accuracy, the above amplitude as, Let us know insert the identity between the operators ψ in k2x2 (t − ) † and ψ out p1ȳ1 (t + ). Recalling the conditions discussed in Sect. II D, namely the localized nature for the ground states together with the fact that there is not particle creation, just {|Ω λ } M −1 λ=0 will contribute to the identity. The final result is: with A 1,ν→λ = Ω λ | ψ out p1ȳ1 (t + )ψ in k1x1 (t − ) † |Ω ν and similarly for A 2,λ→µ . We can generalize this expression to N photons, with initial average positionsx 1 >x 2 > · · · > x N and asymptotic ground states λ 0 := ν and λ N : with The sketched constructive demonstration (a complete demonstration is given in App. C) has confirmed that causality imposes that the amplitude can be built from single photon events whenever those are well separated. Inelastic processes yield the sum over intermediate states.
If only one ground state is considered, the amplitude is the product A = Π n A n . In this case, the S-matrix in momentum space recovers the typical structure in QFT (see Eq. (17)). However, when inelastic-scattering events occur, the sum in (21) leads to a particular structure for the free part of the scattering matrix S 0 that we discuss now. We now find the structure for S 0 in position space compatible with the amplitude (21). For the sake of simplicity, we work with chiral waveguides and a monotonously growing group velocity, ∂ k ω k ≥ 0. Therefore, we can order the events using step functions, eliminating unphysical contributions (e.g. the wave packet ψk 2x2 arriving before than ψk 1x1 , see Fig. 1). Some algebra, fully described in Appendix D yields that S 0 has the following structure The sum over intermediate states and the Heaviside functions are a direct consequence of causality, since they order the different wave packets and keep track of the state of the scatterer for each arrival. Nevertheless, if the ground state is unique (M = 1), the step functions cancel out and we recover the structure described by (18). However, strikingly, for M > 1 this S-matrix cannot be written as a product of one-photon scattering matrices, up to permutations, due to the Heaviside functions. In order to shed light on this, it is convenient to move to momentum space. Although (S 0 pk ) µν cannot be analytically calculated for a general dispersion relation, a mathematical expression can be found for a linear one. This calculation will be presented in Sect. IV B The final result is that (S 0 pk ) µν cannot be written as a product of one-photon S-matrices. This has been recently pointed out in the particular example of a Λ atom by Xu and Fan [36].

IV. APPLICATIONS
The set of previous theorems and conditions create a framework that describes many useful problems and experiments in waveguide QED. We are now going to illustrate two particular problems which are amenable to numerical and analytical treatment, and which highlight the main features of all the results.
The first problem consists of a two-level system that is ultrastrongly coupled to a photonic crystal. The scattering dynamics has to be computed numerically. The simulations fully conform to our our framework, showing the fast decay of photon-qubit dressing with the distance, the independence of space-like separated wave packets, and the decomposition of the two-photon scattering amplitude as a product (for the chosen parameters, the onephoton scattering is elastic).
The second problem consists of a general scatterer with several ground states that is coupled to a non-dispersive medium and it serves to illustrate the breakdown of the S-matrix decomposition in momentum space

A. Ultrastrong scattering
Let us consider a system described by the following Hamiltonian The scatterer is a two-level system described by the ladder operators σ ± and the level splitting ∆. The lattice tight-binding Hamiltonian, describes an array of identical cavities with frequency , cavity-cavity coupling J, and bosonic modes [a x , a † y ] = δ xy . The lattice model is diagonalized in momentum space, giving raise to a cosine-shaped dispersion relation, ω k = − 2J cos k. The scatterer-waveguide interaction, which is described by the last term, is point-like and g is the coupling constant.
The light-matter interaction term can be expressed as a sum of the rotating-wave part, g(σ + a 0 + σ − a † 0 ), and the so-called counter-rotating terms, g(σ − a 0 + σ + a † 0 ). The latter can be neglected if g is small enough compared to the other energies of the full system. This is known as the rotating-wave approximation (RWA). It is well known that the RWA simplifies the problem because (i) the new effective model conserves the number of excitations and (ii) the ground state is the trivial vacuum |vac with σ − |vac = a x |vac = 0 ∀x. However, when the coupling strength is large enough -the so-called ultrastrong coupling regime-, the RWA fails to describe the dynamics and one has to use the full Rabi model (24). This regime not only represents an interesting and challenging problem where we can test our theoretical framework, but it describes a family of current experiments [22,[47][48][49] for which the following simulations are of interest. An important remark is that, despite the fact that the number of excitationsN = x a † x a x + σ + σ − is not a good quantum number, i.e. [H,N ] = 0, numerical simulations indicate that the total number of flying photons is asymptotically conserved throughout the simulation [33,35]. Therefore, the second condition needed for proving our results is fulfilled (see Sect. II D).
We have studied this model using the matrix-productstate variational ansatz, a celebrated method for describing the low-energy sector of one-dimensional many-body systems [50][51][52][53], which has been recently adapted to the photonic world in [33,35,54]. Using this ansatz, we computed the nontrivial minimum-energy state [33], which consists of a photonic cloud exponentially localized around the qubit, see Fig. 2. This result confirms our theoretical predictions from Eq. (13) and implies that the minimum-energy state |Ω 0 can be approximated by the vacuum far away from the qubit.
According to the previous result, we can generate free wave packets, such as input and output states of Eqs. (C1) and (C2)) by inserting photons far away from the scatterer. We have used the MPS ansatz to study the evolution of input states which consist of a pair of photons, see Eq. (C1), with |Ω ν = |Ω 0 . Both wave packets will be Gaussians, Eq. (6), with mean momentumk and width σ. The numerical simulations show that the scattering is elastic for the chosen parameters ( = 1, J = 1/π, ∆ = = 1, and g = 0.3) [33].
We have also demonstrated numerically that the correlation between output photons vanish as the separation between the input wave packet increases. Our study aimed at computing the two-photon wave function in momentum space, φ p1,p2 (t) = Ω 0 | a p1 a p2 |Ψ(t) . This was used to compute the fluorescence F at time t + , the number of output photons whose energy and momentum differ from the input wave packets. More precisely F = dp 1 dp 2 |φ p1,p2 (t + )| 2 , with p 1 and p 2 such that ω p1 + ω p2 = 2(ωk ± σ ω ) and ω p1 , ω p2 ∈ (ωk − σ ω , ωk + σ ω ), being σ ω the width of the input wave packets in energy space. Fig. 3(g) shows F as function of the distance between the incident wave packets. When the wave packets are close enough the fluorescence maximizes and the output wave function shows a nontrivial structure, with φ p1,p2 (t + ) = 0 even though |p 1 | =k or |p 2 | =k (see panels (a) and (c)). The wave function has also a rich structure in position space, with antibunching in the reflection component and superbunching in the transmission one (see panels (b) and (d)). This structure was already found in the RWA [55]. For long distances, the fluorescence F vanishes [see panels (e) and (f)]. In these cases, the output state is clearly uncorrelated: in position space it is formed by two well-defined wave packets and φ p1,p2 (t + ) goes to zero if |p 1 | =k or |p 2 | =k. All this is a consequence of the cluster decomposition, see Eq.  We set ω k = c|k| in H 0 . The scatterer and interaction are described by where {|Ω ν } and {|J } are the ground and decaying states of the scatterer, respectively, {E ν } and {Ẽ J } are their energies, M and M is the number of ground and excited states, respectively, and g J,ν is the coupling strength corresponding to the transition |Ω ν ↔ |J (see Fig. 4). This is a prototypical situation in waveguide QED. E.g., if there are two ground states, M = 2, and the decaying state is unique, M = 1, the scatterer is a Λ atom. From now on, we work in units such that c = 1.
We further assume chiral waveguides: the scatterer only couples to k > 0, which simplifies the final expressions, so we can start from Eq. (23). Before writing down the two-photon S 0 -matrix in momentum space, we need the one-photon scattering matrix. Imposing energy conservation, it has to be with k and p the incident and outgoing momenta, respectively, and |Ω ν and |Ω µ the initial and final ground states. The factor t µν (k) is the so-called transmission amplitude. The Dirac delta guarantees energy conservation. Then, the two-photon S 0 -matrix, Eq. (23) in momentum space is Here, n = n, e.g., n = 2 if n = 1. The computation is detailed in Appendix E. This structure has recently been found by Xu and Fan for a Λ atom (M = 2, M = 1) within the RWA and Markovian approximations [36]. At first sight (29) may look striking. The matrix S 0 is not the product of two Dirac-delta functions conserving the single-photon energy, as discussed in Sect. III B. The mathematical origin of the structure can be traced back to its form in position space, Eq. (29) is also remarkable because presents the natural generalization of the cluster decomposition for the S-matrix [Cf. Eqs. (16) and (17)] when inelastic processes occur in the scattering.
A consequence of (29) is that S 0 contributes to the fluorescence F , Eq. (25). This seems to contradict our previous arguments, since S 0 is built from causally disconnected one-photon events (they do not overlap in the scatterer). To solve the apparent paradox we recall that (29) is a matrix element in momentum space (delocalized photons). For wave packets (5), the scattering amplitude is the integral of these wave packets with (29). In doing so we find that the fluorescence decays to zero as the separation grows, thus solving the puzzle.
In what follows the fluorescence decay is discussed within the full S-matrix, i.e. we consider the contributions to F from S 0 and T (see Eq. (16)). Energy conservation imposes that (T p1p2k1k2 ) µν = (C p1p2k1k2 ) µν δ(p 1 + Since the contribution of T vanishes as the photon-photon separation increases, C must be sufficiently smooth, at least smoother than a Dirac delta [1]. Then, we assume that (C p1p2k1k2 ) µν has simple poles with imaginary parts {γ C n }. Similarly, we expect that divergences of t µν (k) come from simple poles with imaginary parts {γ t n }. As far as we know, this structure has been found for all S-matrices in waveguide QED [25,36,56,57].
Let us write down the input state in momentum space The functions φ 1 (k) and φ 2 (k) are localized far away the scattering region in position space. The exponential factor e ik2l ensures the separation between both wave packets is l. The output state reads |Ψ out = S |Ψ in = µ dp 1 dp 2 φ out µ (p 1 , p 2 )a † p1 a † p2 |Ω µ (31) with the two-photon wave function φ out µ (p 1 , p 2 ) being (C p1p2kn ) µν = dkn(C p1p2knkn ) µν δ(p 1 + p 2 + E µ − k n − kn − E ν ), withn = n. Even though this expression is cumbersome, we can clearly identify the contribution of S 0 and T . We solve this integral by means of the residue theorem. Each pole γ t n and γ C n , together with the exponentials e iknl and e i(p1+p2+Eµ−kn−E λ )l , gives an exponentially decaying term, e −|γ t n |l or e −|γ C n |l . We choose Lorentzian envelopes for the wave packets. They have a pole atk − iσ [see Eq. (7)]. In consequence, the wave packets will give a term proportional to e −σl . Lastly, the imaginary part of the pole of the first term vanishes, ∼ i0 + , so it gives a nondecaying term, e −0 + l = 1. The real part of this denominator imposes the single-photonenergy conservation. Thus, it results in the amplitude for the single-photon events, λ A 1,ν→λ A 2,λ→µ . Therefore, nor S 0 neither T contains fluorescent terms as the separation between the wave packets grows. The technical details are in App. F.
As a final application, one can find experimentally the poles of the one-and two-photon scattering matrices {γ t n } and {γ C n } by measuring the decay of F with the distance.

V. FINAL COMMENTS
Our work represents a significant evolution over the field-theoretical methods [12] that have been so successfully adapted to the study of waveguide QED. Developing an extensive set of theorems shown in the appendices, we have completed a program that derives the properties of the N -photon S-matrix from the emergent causal structure of a nonrelativistic photonic system. This, together with the fact that the ground states of the Hamiltonian are trivial far away from the scatterer and the asymptotic independence of input and output wave packets, allows us to build a consistent scattering theory. Among the consequences of this framework, we have explained how the existence of Raman (inelastic) processes modifies the usual form of the cluster decomposition to produce a structure that includes the particular example developed in [36].
Our formal results also provide insight in the outcome of simulations for problems where no analytical derivation is possible, such as a qubit ultrastrongly coupled to a waveguide [33,35]. As a second example, we have considered a non-dispersive media ω k = c|k|, where we found the general form for the scattering matrix in momentum space (independent of the scatterer and the coupling to the waveguide), which has been recently calculated for a Λ atom [36] as a particular case. On top of that, we have clarified how fluorescence decays in a general scattering experiment.
Throughout the previous discussion we have focused our attention to scattering processes which involve the same number of flying photons both at the input and the output [See Sect. II D], but this is just a convenient restriction that can be lifted. One may incorporate more scattering channels for the photons using extra indices to keep track of the photon-annihilation and creation processes, which results in a slightly more involved version of Theorem 4. In particular, we can incorporate photoncreation events (see e.g. [54]). Finally, our program can be extended to treat other systems, deriving a cluster decomposition for the scattering of spin waves in quantummagnetism models or for fermionic excitations in manybody systems.
Appendix A: The ground state of the light-matter interaction In this appendix we demonstrate that the ground state converges to the trivial vacuum far away from the scatterer, Eq. (13). The next lemma is neccessary to proof the main theorem. Lemma 1. Given the waveguide-QED model (2), we have the following bounds for the expectation values on its minimum-energy state |Ω 0 , Proof. Let us assume that |Ω 0 is the minimum-energy state of H as given by Eq. (2), and thus (H−E 0 ) |Ω 0 = 0. The energy of the unnormalized state created by any operator O must be larger or equal to that of the ground state, χ|(H − E 0 )|χ ≥ 0. Using (A2) we conclude with the useful relation Let us take O = a k . The previous statement leads to or equivalently Using Cauchy-Schwatz, this translates into the upper bound (A7) Once the diagonal elements of the correlation matrix are bounded the nondiagonal can also be bounded. The correlation matrix is positive C ≥ 0 with C kp = Ω 0 |a † k a p |Ω 0 . A property of positive matrices is [58] |C kp | ≤ |C kk ||C pp | (A8) which implies (A1).
With this lemma at hand we state: Theorem 1. Let us define ψ † kxx as the operator (5) removing the time-dependent part, where φk(k) is infinitely differentiable with a finite support K centered aroundk. Then, the expected value of ψ † kx ψkx in the minimumenergy state fulfills where we choose x sc = 0. Moreover, if we can assume that a † k a p is an n-times differentiable function of k and p, the bound will be improved Proof. Let us compute the expectation value of the number operator for a wave packet N := Ω 0 |ψ † kx ψkx|Ω 0 , We can rewrite N as the Fourier transform of another We are now going to assume that φk(k) is a test function with compact support K of size |K| centered aroundk, and infinitely differentiable. We will also assume that within its support |g k /ω k | 2 GG † ≤ C φ for some constant C φ . Then we can bound Assuming that Ω 0 |a † k a p |Ω 0 is n-times differentiable and using the Riemann-Lebesgue theorem, we have then that at long distances. We first prove causal relations in a free theory. In order to do so, we work with localized wave packets ψkx(t), Eq. (5). Actual calculations are done with Gaussian wave packets, Eq. (6). The following two lemmas are used in the demonstration of the theorem. Lemma 2. Let the dispersion relation ω k have an upper bounded group velocity v k = ∂ k ω k : Then, the function f (k) = kx − ω k t only has stationary points if the distance to the light cone is nonnegative. In other words Then, provided f (k) = 0, it follows |x| ≤ c|t| ⇒ d c (x, t) ≤ 0, which shows (B2).
Lemma 3. Assume that ω k is n-times differentiable and that every derivative |ω (r≤n) k | is upper bounded by an mth order polynomial in |k|. Then the following integral bound applies where p(k) is a polynomial of degree r.
Proof. Result 5.1 from Ref. [59] states that the integral I(x) = b a e ikx q(k) dk may be integrated by parts n times, obtaining where the error term satisfies provided that q(k) is n-times differentiable and that q (n) ∈ L 1 . Based on the conditions of the lemma, this is satisfied since q(k) = e − 1 σ 2 (k−k0) 2 −iω k t p(k). The limits of the integral may be easily extended to ±∞, as explained in Result 5.2 from [59]. Since x −s q (s) (a) → 0 when a → ±∞, ∀x, we obtain Moreover, q (n) , resulting from a product of derivatives of ω k t, −k 2 /σ 2 and the polynomial p(k) of degree r, is bounded by a polynomial of at most (m + n + r)-th order in |k|. Such a polynomial is integrable together with the Gaussian wave packet giving a constant prefactor. In estimating this factor, we can take the worst-case scenario for the terms in t, which appears at most n times together with (∂ k ω k ) n , and the monomials in |k|, which produce another prefactor σ m+n+r .
Note that it would suffice to consider q(k) as a test function or even a Schwartz function since in this case all the differentiability requisities are fullfilled and x −s q (s) (a) → 0 → 0 when a → ±∞, ∀x still holds, because these functions and their derivatives are rapidly decreasing.
With these lemmas at hand we can prove Theorem 2. Let the Hamiltonian be given just by the photonic part, H 0 = dk ω k a † k a k . Let ψkx(t) and ψpȳ(t ) denote two localized wave packets of the form (6). We will assume that (i) the absolute value for the group velocity of these wave packets is upper bounded by a constant c within the domain of the wave packets (|v k | = |∂ k ω k | ≤ c) and (ii) the dispersion relation is n-times differentiable and that each derivative is upper bounded by a polynomial of at most order m: The commutator between these wave packets is small whenever they are outside of their respective light cones, that is, whenever d = |ȳ −x| − c|t − t| 0, Proof. Let us assume that the model evolves freely according to the free Hamiltonian H 0 = dk ω k a † k a k . In this case, our wave packet operators have the simple form and analogously for ψpȳ(t ). The commutator between operators reads using Lemma 2 we know that the exponent has no stationary point. Assuming w.l.o.g.x >ȳ, t > t (other combinations are analogous) and writingω k = ω k − ck, we obtain The exponentω k = ω k − ck is n-times differentiable and is upper bounded in modulus by a polynomial of degree m ≥ 1. Lemma 3 therefore allows us to bound the commutator by a term O(d −n ).
Note that for a linear dispersion, ω k = c|k|, we can rewrite this integral as a function of the distance between world lines from Eq. (B2), d = (x −ȳ) − c(t − t ). Introducing k ± = (k ±p)/2 and using our Gaussian wave packets (6), we obtain This bound is better than the one we have found but it is compatible with Lemma 3 and Theorem 2.

Full model causality
Causal relation (B8) can be extended to the full model (2). Theorem 3. Let H be the light-matter Hamiltonian given by Eq. (2). We assume the conditions of Theorem 2: differentiable, polynomially bounded functions ω k and g k , with degrees n ≥ 2. Then, all wave packets outside the light cone of the scatterer evolve approximately with the free Hamiltonian, H 0 . More precisely, if (x, t 1 ) and (x, t 0 ) are two points outside the light cone is the free-evolution operator for the photons at time t 0 .
Proof. We start by building the Heisenberg equations for the operators Making the change of variables a k (t) = e −iω k t b k (t), we have so that the wave packet operators evolved from some initial time t s are The first part corresponds to free evolution, while the second part is an error term ε(t), which can be bounded. We will assume without loss of generality G = 1, with || · || the Hilbert-Schmidt norm, and |t 1 | > |t 0 |. We have to choose the integration limits t and t s so that sign(τ ) = sign(x). If x > 0 then t 1 > t 0 > 0 and (t, t s ) = (t 1 , t 0 ) is a good choice. If x < 0 then 0 > t 0 > t 1 and again (t, t s ) = (t 1 , t 0 ) is also a valid choice (τ < 0). This means we can introduce τ = sign(x)τ ≥ 0 and bound (B20) Here we have taken into account that d c (|x|, τ ) ≥ d c (|x|, |t 1 − t 0 |) > 0 in the domain of integration. We can now use the fact that d c (|x|, obtaining the expression in the theorem.

Asymptotic Condition
One important limitation of Theorem 3 is that it is focused on the operators, not on the states themselves. This is a key point. For having a well defined scattering theory, the asymptotic condition must holds [See Sect II C and Eq. (9)]. However, using Theorems 3 and 1 we have that, given a state |Ψ ≡ ψk ,x (t 0 ) † |Ω ν , then The first equality is up to a global phase. In the second line, we have used Theorem 3. In the last line, we can introduce input (output) states since the wave packets are well separated (t ± → ±∞) from the scatterer and, by means of Theorem 1 and the conditions presented in II D they are well defined free particle states.
This last result warrants that, under rather general conditions, the light-matter Hamiltonian (2) gives a physical scattering theory.
with |x n −x m | → ∞ ∀n = m. Thus, the scattering amplitude of going to with |ȳ n −ȳ m | → ∞ ∀n = m, is reduced to a product of single-photon events: being λ 0 = µ and λ N = ν, with the wave packet operators given in the Heisenberg picture for t = t ± → ±∞.
The proof is based directly on causality. Therefore, we find convenient to discuss it here.
Proof. The proof is done for the two-photon scattering. The generalization for N photons is straightforward. The scattering operator S is nothing but the evolution operator in the interaction picture, cf. Eq. (11). This permits to write the scattering amplitudes as, In the second equality we have dropped an irrelevant global phase. Here, ψ † in and ψ † out are operators creating wave packets localized far away from the scatterer. Because of Theorem 1, they are well defined N -photon wave packets.
Using Eqs. (C1) and (C2) the amplitude is given by As |x 1 −x 2 | can be arbitrarily large, we can always choose a time t 1 such that ψ out p1ȳ1 (t) † |Ω µ is well separated from the scatterer for t > t 1 , so ψ out . Besides, t 1 is such that the second wave packet is still far away from the scatterer. Therefore ψ in Finally, we insert the identity between the operators ψ in k2x2 (t − ) † and ψ out p1ȳ1 (t + ). Assuming there is not particle creation and just the ground states {|Ω λ } M −1 λ=0 will contribute to the identity, M −1 λ=0 |Ω λ Ω λ |, and we arrive to (C3).
This comes because ψ in k2x2 (t − ) † and ψ out p1ȳ1 (t + ) asymptotically commute but not ψ in k1x1 (t − ) † and ψ out p2ȳ2 (t + ). This is a clear signature of causality, saying which one is arriving first. Lastly, notice that if the ground state is unique, |Ω λn = |Ω 0 , this ordering is not important as the amplitude is simply the product of single-photon scattering amplitudes.
Appendix D: Scattering amplitude from Eq. (23) In this appendix, we prove that (23) is consistent with the amplitude factorization from Theorem 4, Eq. (C3). We do it in the two-photon subspace.
Before, we need the one-photon amplitude as an intermediate result.

One photon
We first need to compute the one photon amplitude. Let the one-photon input state be, with the creation operator ψ in † k1,x1 given by Eq. (5), removing the time dependence. For simplicity, we absorb the factor e ikx1 into the wave packet: φk 1 ,x1 (k) = e ikx1 φk 1 (k). In position space, the output state will read (D2) Defining and |ξ 1 out 1,µν = dy φ 1,µν (y)|y; Ω µ , being |y; Ω µ = a † y |Ω µ the state with a photon at y and the scatterer in the ground state |Ω µ , the output state (D2) can be rewritten as The probability amplitude will read If the wave packets are monochromatic with momenta k 1 and p 1 , respectively, this amplitude is

Two photons
The two-photon wave packet, as sketched in Fig. 1, is By definition, the output state is Here, we are interested in he limit of well separated incident photons. Thus, only the linear part of the scattering matrix S 0 is considered. We introduce the identity operator with being |x 1 x 2 ; Ω ν = a † x1 a † x2 |Ω µ the symmetrized state with two photons at x 1 at x 2 and the scatterer at |Ω ν .
Finally, the probability amplitude of going to the output state ψ out † p1,ȳ1 ψ out † p2,ȳ2 |Ω µ will be the overlap between this state and (D16). Using (D6) as expected.
In the calculations, we have set Ω µ |ψ out pi,ȳi S ψ in † kjxj |Ω ν = 0 for i = j, since we assume that both incident wave packets are far away.
A final comment is in order. Without the step functions in (23), the unphysical amplitude A 2,ν→λ A 1,λ→µ would appear in the final probability amplitude.
Appendix E: S 0 in momentum space Here, we show S 0 in momentum space follows Eq. (29). After that, we prove the Dirac-delta structure is recovered if the ground state is unique.
Let us write (S 0 p1p2k1k2 ) µν as the Fourier transform of (S 0 y1y2x1x2 ) µν Due to the form of (S 0 y1y2x1x2 ) µν , (23), we have to compute integrals as Notice that (S yx ) µν is the Fourier transform of (S pk ) µν , Eq. (28). Therefore, Considering this in (E1), we get with n = n and m = m. The Fourier transform of the step function is Therefore, integrating Eq. (E4) first in y 1 and later in y 2 , we get which is the expression given in the main text, Eq. (29). This result has been recently reported for a Λ atom by Xu and Fan in [36]. Here, we show this is completely general due to our ansatz (Eq. (23)). Lastly, we prove that Eq. (29) is formed by two Dirac-delta functions if M = 1. To do so, we use the following identity with P the principal value. Applying this identity to Eq. (29) we get, t(k n )t(k n ) −iπδ(p m − k n ) + P 1 p m − k n δ(p 1 + p 2 − k 1 − k 2 ).
In these expressions, the wave packets φk n (k) are Lorentzian functions [see Eq. (7)]. The out state is com-puted by means of Eq. (10) |Ψ out = S |Ψ in = ISI |Ψ in . (F3) With I the identity operator in the two-photon sector: I = 1/2 dp 1 dp 2 µ a † p1 a † p2 |Ω µ Ω µ | a p1 a p2 . The scattering matrix S in momentum space is (S p1p2k1k2 ) µν = (S 0 p1p2k1k2 ) µν + i(T p1p2k1k2 ) µν , with (S 0 p1p2k1k2 ) µν given by Eq. (29) and (T p1p2k1k2 ) µν = (C p1p2k1k2 ) µν δ(p 1 + p 2 + E µ − k 1 − k 2 − E ν ) yielding |Ψ out = dp 1 dp 2 Which is nothing but Eq. (32) that we have rewritten here for the discussion. As said in Sect. IV B, we assume that t µν (k) and (C p1p2knkn ) µν have simple poles with imaginary parts {γ t n } and {γ C n } respectively. Then, this integral is solved by taking complex contours and applying the residue theorem. In order to integrate the term proportional to e i(p1+p2+Eµ−kn−E λ )l , we take the contour shown in Fig. 5(a) so that the exponential factor does not diverge. For the same reason, for that proportional to e iknl we take the contour of Fig. 5(b). As t and C have first-order poles, when integrating each pole, we just have to evaluate the rest of the function at the pole. Then, t and C give terms proportional to e −|γ t n |l and e −|γ C n |l , respectively. Now we consider the contribution to the integral of the wave packets, φk n (k). We choose Lorentzian functions, with a simple pole at k =k n − iσ (see Eq. (7)). In consequence, we have a term proportional to e −σl . Lastly, the denominator in the first term has a pole with zero imaginary part. Therefore, its contribution does not decay with l. Importantly enough, this pole enforces singlephoton energy conservation giving single-photon amplitudes, λ A 1,ν→λ A 2,λ→µ .
Finally, let us mention that we do not need to impose that that t µν (k) and (C p1p2knkn ) µν have simple poles. Higher order poles, by virtue of the Cauchy Integral formula for the derivatives, also would yield exponential decay. (32). We show the poles coming from the Lorentzian, ±iσ, those coming from one of the transmission amplitudes or from C, ±iΓ, and those with vanishing imaginary part. The real parts are arbitrary.