Emergent causality and the N-photon scattering matrix in waveguide QED

The simplest model that describes a waveguide-QED setup consists of a one-dimensional bosonic medium and a scatterer. Using units such that ${\hslash }=1$, it reads

$$\begin{eqnarray}&&H={H}_{0}+{H}_{\mathrm{sc}}+\int ({g}_{k}{G}^{\dagger }{a}_{k}+{g}_{k}^{* }{{Ga}}_{k}^{\dagger }){\rm{d}}{k}.\end{eqnarray} \tag{ 2 }$$

The first term stands for the free-Hamiltonian of the photons

$$\begin{eqnarray}&&{H}_{0}=\int {\omega }_{k}\,{a}_{k}^{\dagger }{a}_{k}{\rm{d}}{k},\end{eqnarray} \tag{ 3 }$$

with frequency ${\omega }_{k}$ for momentum k, which is created (annihilated) by the corresponding Fock operator a_k (${a}_{k}^{\dagger }$), satisfying $[{a}_{k},{a}_{{k}^{{\prime} }}^{\dagger }]=\delta (k-{k}^{{\prime} })$. The last two terms are the Hamiltonian ${H}_{\mathrm{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

$$\begin{eqnarray}&&{g}_{x}=\displaystyle \frac{1}{\sqrt{2\pi }}\int {\rm{d}}{k}\,{{\rm{e}}}^{{\rm{i}}{k}{x}}{g}_{k}\end{eqnarray} \tag{ 4 }$$

have a finite support centered around ${x}_{\mathrm{sc}}=0$. The model (2) is not exactly solvable in general. For instance, if the scatterer is a two-level system, ${H}_{\mathrm{sc}}\propto {\sigma }_{z}$ and $G={\sigma }_{x}$ the model is the celebrated spin-boson model [39], which results in a nontrivial ground state with localized photonic excitations around the scatterer.

The discussion below assumes a single photonic band ${\omega }_{k}\in [{\omega }_{\min },{\omega }_{\max }]$ and typically a chiral medium $k\geqslant 0$, ${\partial }_{k}{\omega }_{k}\geqslant 0$. This is a rather standard simplification which does not affect the generality and applicability of our results. We generalize our results to nonchiral media in appendix E. The structure for S⁰ is essentially identical in that case, so the conclusions of the paper also hold for nonchiral media. Anyway, chiral waveguides can be realized experimentally, for instance with photonic crystals [40, 41], or nanofibers coupled to nanoparticles or atoms [14, 42]. Besides, we can consider more generic dispersion relations by introducing additional degrees of freedom in the photons (band index, etc) and keeping track of those quantum numbers in a trivial extension of our results.

Our model is lossless. Losses are negligible in several experimental platforms, such as superconducting transmission lines interacting with superconducting qubits [21, 22]. Besides, the one-dimensional approximation is valid for all the implementations considered in the references we included above.

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 ${\psi }_{\bar{k}\bar{x}}{(t)}^{\dagger }$ for a wave packet as

$$\begin{eqnarray}&&{\psi }_{\bar{k}\bar{x}}{(t-{t}_{0})}^{\dagger }=\int {{\rm{e}}}^{{\rm{i}}{k}\bar{x}-{\rm{i}}{\omega }_{k}(t-{t}_{0})}{\phi }_{\bar{k}}(k){a}_{k}^{\dagger }\,{\rm{d}}{k}.\end{eqnarray} \tag{ 5 }$$

The wavefunction ${\phi }_{\bar{k}}(p)=\phi (p-\bar{k})\in {{ \mathcal L }}^{2}$ is normalized and centered around the average momentum $\bar{k}$. The exponential factor ${{\rm{e}}}^{{\rm{i}}{k}\bar{x}}$ ensures the wave packet is centered around $\bar{x}$ in position space at time $t={t}_{0}$.

As wave packets we will use both Gaussian

$$\begin{eqnarray}&&{\phi }_{\bar{k}}(k)=\displaystyle \frac{1}{\sqrt[4]{2\pi }\sqrt{\sigma }}\exp [-{(k-\bar{k})}^{2}/4{\sigma }^{2}],\end{eqnarray} \tag{ 6 }$$

and Lorentzian envelopes

$$\begin{eqnarray}&&{\phi }_{\bar{k}}(k)=\sqrt{\displaystyle \frac{\sigma }{\pi }}\displaystyle \frac{1}{k-\bar{k}+{\rm{i}}\sigma }.\end{eqnarray} \tag{ 7 }$$

These wave functions are only approximately localized in the sense that the probability of finding a photon decays exponentially far away from the center $\bar{x}$. The width σ in momentum space implies a localization length $1/\sigma$ in position space.

Figure 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 ${\bar{k}}_{1}$ or ${\bar{k}}_{2}$ determines the group velocity at which the photons move ${v}_{g}(k)={\partial }_{k}{\omega }_{k}$. The wave packets may be distorted due both to the dispersive nature of the medium and the interaction with the scatterer.

Zoom In Zoom Out Reset image size

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 $| {{\rm{\Omega }}}_{\mu }\rangle$ are indistinguishable from the vacuum state $| \mathrm{vac}\rangle$ far away from the scatterer. Mathematically this occurs when

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\bar{x}\to \pm \infty }\langle {{\rm{\Omega }}}_{\mu }| O(\bar{x},{\rm{\Delta }})| {{\rm{\Omega }}}_{\mu }\rangle =\langle \mathrm{vac}| {O}(\bar{x},{\rm{\Delta }})| \mathrm{vac}\rangle ,\end{eqnarray} \tag{ 8 }$$

where $O(\bar{x},{\rm{\Delta }})$ is an operator with compact support in the finite interval $\bar{x}-{\rm{\Delta }}/2\lt x\lt \bar{x}+{\rm{\Delta }}/2$ and the vacuum state $| \mathrm{vac}\rangle$ is such that ${a}_{k}| \mathrm{vac}\rangle =0$ $\forall k$.

Besides, the free particle states must satisfy the asymptotic condition [43]:

$$\begin{eqnarray}&&\parallel U(t)| {\rm{\Psi }}\rangle -{U}^{0}(t)| {{\rm{\Psi }}}_{\mathrm{in}/\mathrm{out}}\rangle \parallel \mathop{\longrightarrow }\limits^{t\to \mp \infty }0,\end{eqnarray} \tag{ 9 }$$

with U(t) the evolution operator of the full Hamiltonian (2) and ${U}^{0}(t)={{\rm{e}}}^{-{{\rm{i}}{H}}_{0}t}$ the free-evolution operator.

The scattering operator S relates the amplitude of the output and input fields through

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}\rangle =S\,| {{\rm{\Psi }}}_{\mathrm{in}}\rangle ,\end{eqnarray} \tag{ 10 }$$

which, using (9), has the formal expression:

$$\begin{eqnarray}&&S=\mathop{\mathrm{lim}}\limits_{{t}_{\pm }\to \pm \infty }{U}_{I}({t}_{+},{t}_{-}).\end{eqnarray} \tag{ 11 }$$

Here, ${U}_{I}({t}_{+},{t}_{-})={{\rm{e}}}^{{{\rm{i}}{H}}_{0}{t}_{-}}\,{{\rm{e}}}^{-{\rm{i}}{H}({t}_{+}-{t}_{-})}\,{{\rm{e}}}^{-{{\rm{i}}{H}}_{0}{t}_{+}}$ is the evolution operator in the interaction picture. Using again equation (9) leads to $| {{\rm{\Psi }}}_{\mathrm{in}/\mathrm{out}}\rangle ={U}_{0}^{\dagger }({t}_{-/+})U({t}_{-/+})| {\rm{\Psi }}\rangle \equiv | {\rm{\Psi }}({t}_{-/+}){\rangle }_{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:

$$\begin{eqnarray}&&A\equiv \langle {{\rm{\Omega }}}_{\mu }| {\psi }^{\mathrm{out}}({t}_{+})S{\psi }^{\mathrm{in}}{({t}_{-})}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle \end{eqnarray} \tag{ 12 }$$

with ${\psi }^{\mathrm{in}}{({t}_{-})}^{\dagger }={\psi }_{\bar{k}\bar{x}}{({t}_{-})}^{\dagger }$ and an analogous definition for ${\psi }^{\mathrm{out}}({t}_{+})$ and the photon mean position $\bar{x}$ 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.

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 equation (8) (see appendix A) for the unique ground state of Hamiltonian (2), which reads

$$\begin{eqnarray}&&\langle {{\rm{\Omega }}}_{0}| {\psi }_{\bar{k}\bar{x}}^{\dagger }{\psi }_{\bar{k}\bar{x}}| {{\rm{\Omega }}}_{0}\rangle \leqslant { \mathcal O }(| \bar{x}{| }^{-n}),\,| \bar{x}| \to \infty \end{eqnarray} \tag{ 13 }$$

provided that: (i) for all k, $| {g}_{k}/{\omega }_{k}| \lt \infty$ and (ii) that the correlators ${C}_{{kp}}=\langle {{\rm{\Omega }}}_{0}| {a}_{k}^{\dagger }{a}_{p}| {{\rm{\Omega }}}_{0}\rangle$ 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 [34, 35, 44, 45], 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 [35].

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, $\{| {{\rm{\Omega }}}_{\mu }\rangle \}$, which are localized in the sense of equation (8). Notice that with this assumption (2) has a well defined theory (see appendix B.3). This condition allows the expression of the elements of S in the momentum basis:

$$\begin{eqnarray}&&{({S}_{{\bf{pk}}})}_{\mu \nu }=\left\langle {{\rm{\Omega }}}_{\mu }| \displaystyle \prod _{i}{a}_{{p}_{i}}S\displaystyle \prod _{j}{a}_{{k}_{j}}^{\dagger }| {{\rm{\Omega }}}_{\nu }\right\rangle .\end{eqnarray} \tag{ 14 }$$

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}_{\mathrm{sc}}\propto {\sigma }_{z}$ and $G={\sigma }_{x}$) exactly conserves the number of excitations within the RWA, 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 (see [35, 37, 46] and section 4.1).

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 ({\partial }_{k}{\omega }_{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(x-y,t-{t}^{{\prime} })=| \bar{x}-\bar{y}| -c| t-t^{\prime} |$ and prove in appendix B that

$$\begin{eqnarray}&&\parallel [{\psi }_{\bar{k}\bar{x}}(t),{\psi }_{\bar{p}\bar{y}}{(t^{\prime} )}^{\dagger }]\parallel ={ \mathcal O }\left(\displaystyle \frac{1}{| D{| }^{n}}\right)+{ \mathcal O }\left(\displaystyle \frac{1}{| {D}_{0}{| }^{n-1}}\right),\end{eqnarray} \tag{ 15 }$$

with $D\equiv d(x-y,t-{t}^{{\prime} })$ and ${D}_{0}\equiv \min \{d(\bar{x},t),d(\bar{x},{t}_{0}),d(\bar{y},t),d(\bar{y},{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₀. In the Heisenberg picture, the phases ${\rm{i}}{k}(\bar{x}-\bar{y})-{\rm{i}}{\omega }_{k}(t-{t}^{{\prime} })$ can be bounded by the distance $d(x-y,t-{t}^{{\prime} })$. Using the Riemann–Lebesgue lemma ($\int {{\rm{e}}}^{{\rm{i}}{k}{z}}f(k)\,{\rm{d}}{k}\to 0$, as $z\to \infty$) 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₀, 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, (an)harmonic oscillators, master equations, and spin-boson lattices [6, 7, 47–51]. It is important to remark that the approximate causality in equation (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.

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⁰ and an interacting part T, both in momentum space

$$\begin{eqnarray}&&{S}_{{\bf{pk}}}={S}_{{\bf{pk}}}^{0}+{{\rm{i}}{T}}_{{\bf{pk}}}.\end{eqnarray} \tag{ 16 }$$

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 $| {\bar{x}}_{i}-{\bar{x}}_{j}| \to \infty$. This, together with energy conservation, imposes the constraint ${{\rm{i}}{T}}_{{\bf{pk}}}={{\rm{i}}{C}}_{{\bf{pk}}}\delta ({E}_{{\bf{p}}}-{E}_{{\bf{k}}})$ [1]. In this limit the only term contributing to the scattering amplitude is the free part, S⁰. In QFT (typically) occurs momentum conservation which implies that

$$\begin{eqnarray}&&{S}_{{\bf{pk}}}^{0}=\displaystyle \frac{1}{N!}\left(\displaystyle \prod _{n=1}^{N}{S}_{{p}_{n}{k}_{n}}+\mathrm{permutations}[{k}_{n}\leftrightarrow {k}_{m},{p}_{n}\leftrightarrow {p}_{m}]\right),\end{eqnarray} \tag{ 17 }$$

with ${S}_{{p}_{n}{k}_{n}}\propto \delta ({\omega }_{{p}_{n}}-{\omega }_{{k}_{n}})$ the one-photon S-matrix. In order to clarify the notation $\mathrm{permutations}[{k}_{n}\leftrightarrow {k}_{m},{p}_{n}\leftrightarrow {p}_{m}]$, let us write the two-photon S matrix:

$$\begin{eqnarray}&&{S}_{{\bf{pk}}}^{0}=\displaystyle \frac{1}{2}({S}_{{p}_{1}{k}_{1}}{S}_{{p}_{2}{k}_{2}}+{S}_{{p}_{2}{k}_{1}}{S}_{{p}_{1}{k}_{2}}+{S}_{{p}_{2}{k}_{2}}{S}_{{p}_{1}{k}_{1}}+{S}_{{p}_{1}{k}_{2}}{S}_{{p}_{2}{k}_{1}})={S}_{{p}_{1}{k}_{1}}{S}_{{p}_{2}{k}_{2}}+{S}_{{p}_{2}{k}_{1}}{S}_{{p}_{1}{k}_{2}}.\end{eqnarray} \tag{ 18 }$$

This is nothing but the cluster decomposition. Fourier transforming ${S}_{{\bf{pk}}}^{0}$, this structure also holds

$$\begin{eqnarray}&&{S}_{{\bf{yx}}}^{0}=\displaystyle \frac{1}{N!}\left(\displaystyle \prod _{n={1}^{N}}{S}_{{y}_{n}{x}_{n}}+\mathrm{permutations}[{x}_{n}\leftrightarrow {x}_{m},{y}_{n}\leftrightarrow {y}_{m}]\right).\end{eqnarray} \tag{ 19 }$$

This shall be relevant in the following section, where we will work in position space.

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⁰ may not have the structure given by equation (17).

To understand how causality fixes the form of S⁰ we refer to our figure 1 where two well separated wave packets interact with a scatterer. The scattering amplitude is,

$$\begin{eqnarray*}&&A=\langle {{\rm{\Omega }}}_{\mu }| \displaystyle \prod _{m=1}^{2}{\psi }_{{\bar{p}}_{m}{\bar{y}}_{m}}^{\mathrm{out}}({t}_{+})\displaystyle \prod _{n=1}^{2}{\psi }_{{\bar{k}}_{n}{\bar{x}}_{n}}^{\mathrm{in}}{({t}_{-})}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle .\end{eqnarray*}$$

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 equation (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,

$$\begin{eqnarray}&&A\simeq \langle {{\rm{\Omega }}}_{\mu }| {\psi }_{{\bar{p}}_{2}{\bar{y}}_{2}}^{\mathrm{out}}({t}_{+}){\psi }_{{\bar{k}}_{2}{\bar{y}}_{2}}^{\mathrm{in}}{({t}_{-})}^{\dagger }{\psi }_{{\bar{p}}_{1}{\bar{x}}_{1}}^{\mathrm{out}}({t}_{+}){\psi }_{{\bar{k}}_{1}{\bar{x}}_{1}}^{\mathrm{in}}{({t}_{-})}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle .\end{eqnarray} \tag{ 20 }$$

Let us know insert the identity between the operators ${\psi }_{{\bar{k}}_{2}{\bar{x}}_{2}}^{\mathrm{in}}{({t}_{-})}^{\dagger }$ and ${\psi }_{{\bar{p}}_{1}{\bar{y}}_{1}}^{\mathrm{out}}({t}_{+})$. Recalling the conditions discussed in section 2.4, namely the localized nature for the ground states together with the fact that there is not particle creation, just $\{| {{\rm{\Omega }}}_{\lambda }\rangle \}{}_{\lambda =0}^{M-1}$ will contribute to the identity. The final result is:

$$\begin{eqnarray}&&{A}_{12}=\displaystyle \sum _{\lambda =0}^{M-1}{A}_{1,\nu \to \lambda }{A}_{2,\lambda \to \mu },\end{eqnarray} \tag{ 21 }$$

with ${A}_{1,\nu \to \lambda }=\langle {{\rm{\Omega }}}_{\lambda }| {\psi }_{{\bar{p}}_{1}{\bar{y}}_{1}}^{\mathrm{out}}({t}_{+}){\psi }_{{\bar{k}}_{1}{\bar{x}}_{1}}^{\mathrm{in}}{({t}_{-})}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle$ and similarly for ${A}_{2,\lambda \to \mu }$. We can generalize this expression to N photons, with initial average positions ${\bar{x}}_{1}\gt {\bar{x}}_{2}\,\gt ...\gt \,{\bar{x}}_{N}$ and asymptotic ground states ${\lambda }_{0}:= \nu$ and ${\lambda }_{N}:= \mu$

$$\begin{eqnarray}&&A=\displaystyle \sum _{{\lambda }_{1},\ldots ,{\lambda }_{N-1}=0}^{M-1}\displaystyle \prod _{n=1}^{N}{A}_{n,{\lambda }_{N+1-n}\to {\lambda }_{N-n}},\end{eqnarray} \tag{ 22 }$$

with

$$\begin{eqnarray}&&{A}_{n,{\lambda }_{N+1-n}\to {\lambda }_{N-n}}=\langle {{\rm{\Omega }}}_{{\lambda }_{N-n}}| {\psi }_{{\bar{p}}_{n}{\bar{y}}_{n}}^{\mathrm{out}}({t}_{+}){\psi }_{{\bar{k}}_{n}{\bar{x}}_{n}}^{\mathrm{in}}{({t}_{-})}^{\dagger }| {{\rm{\Omega }}}_{{\lambda }_{N+1-n}}\rangle .\end{eqnarray} \tag{ 23 }$$

The sketched constructive demonstration (a complete demonstration is given in appendix 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={{\rm{\Pi }}}_{n}{A}_{n}$. In this case, the $S$-matrix in momentum space recovers the typical structure in QFT (see equation (17)). However, when inelastic-scattering events occur, the sum in (22) leads to a particular structure for the free part of the scattering matrix S⁰ that we discuss now.

We now find the structure for S⁰ in position space compatible with the amplitude (22). For the sake of simplicity, we work with chiral waveguides and a monotonously growing dispersion relation, ${\partial }_{k}{\omega }_{k}\geqslant 0$. Therefore, we can order the events using step functions, eliminating unphysical contributions (e.g. the wave packet ${\psi }_{{\bar{k}}_{2}{\bar{x}}_{2}}$ arriving before than ${\psi }_{{\bar{k}}_{1}{\bar{x}}_{1}}$, see figure 1). Some algebra, fully described in appendix D yields that S⁰ has the following structure

$$\begin{eqnarray}&&{({S}_{{\bf{yx}}}^{0})}_{\mu \nu }=\displaystyle \sum _{{\lambda }_{1}...{\lambda }_{N-1}=0}^{M-1}\displaystyle \prod _{n=1}^{N}{({S}_{{y}_{n}{x}_{n}})}_{{\lambda }_{N+1-n}{\lambda }_{N-n}}\displaystyle \prod _{m=1}^{N-1}\theta ({y}_{m+1}-{y}_{m})\,+\mathrm{permutations}[{x}_{n}\leftrightarrow {x}_{m},{y}_{n}\leftrightarrow {y}_{m}].\end{eqnarray} \tag{ 24 }$$

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 (19). However, strikingly, for $M\gt 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}_{{\bf{pk}}}^{0})}_{\mu \nu }$ 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 section 4.2. The final result is that ${({S}_{{\bf{pk}}}^{0})}_{\mu \nu }$ cannot be written as a product of one-photon S-matrices. This has been recently pointed out in the particular example of a Λ atom coupled to a waveguide within the RWA and Markovian approximations with point-like coupling by Xu and Fan [38] using the input–output formalism [27].

As we said in the introduction, the generalization to nonchiral waveguides is straightworward. We explain the details in appendix E.

Let us consider a system described by the following Hamiltonian

$$\begin{eqnarray}&&H={\rm{\Delta }}{\sigma }^{+}{\sigma }^{-}+\epsilon \displaystyle \sum _{x}{a}_{x}^{\dagger }{a}_{x}-J\displaystyle \sum _{x}({a}_{x}^{\dagger }{a}_{x+1}+{a}_{x+1}^{\dagger }{a}_{x})+g({\sigma }^{-}+{\sigma }^{+})({a}_{0}+{a}_{0}^{\dagger }).\end{eqnarray} \tag{ 25 }$$

The scatterer is a two-level system described by the ladder operators ${\sigma }^{\pm }$ 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}^{\dagger }]={\delta }_{{xy}}$. We work in units such that the lattice spacing is the unit of length.

The lattice model (second and third terms of the Hamiltonian) is diagonalized in momentum space, giving raise to a cosine-shaped dispersion relation, ${\omega }_{k}=\epsilon -2J\cos k$. Some recent implementations of such cosine-shaped dispersion are superconducting coupled resonators and photonic crystals (see [23, 40, 41] respectively).

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 RWA part, $g({\sigma }^{+}{a}_{0}+{\sigma }^{-}{a}_{0}^{\dagger })$, and the so-called counter-rotating terms, $g({\sigma }^{-}{a}_{0}+{\sigma }^{+}{a}_{0}^{\dagger })$. The latter can be neglected if g is small enough compared to the other energies of the full system. This is known as the 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 $| \mathrm{vac}\rangle$ with ${\sigma }^{-}| \mathrm{vac}\rangle ={a}_{x}| \mathrm{vac}\rangle =0$ $\forall 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 (25). 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 [24, 52–54] for which the following simulations are of interest. An important remark is that, despite the fact that the number of excitations $\hat{N}={\sum }_{x}{a}_{x}^{\dagger }{a}_{x}+{\sigma }^{+}{\sigma }^{-}$ is not a good quantum number, i.e. $[H,\hat{N}]\ne 0$, numerical simulations indicate that the total number of flying photons is asymptotically conserved throughout the simulation [35, 37]. Therefore, the second condition needed for proving our results is fulfilled (see section 2.4).

We have studied this model using the matrix-product-state variational ansatz, a celebrated method for describing the low-energy sector of one-dimensional many-body systems [55–58], which has been recently adapted to the photonic world in [35, 37, 59]. Using this ansatz, we computed the nontrivial minimum-energy state [35], which consists of a photonic cloud exponentially localized around the qubit, see figure 2. This result confirms our theoretical predictions from equation (13) and implies that the minimum-energy state $| {{\rm{\Omega }}}_{0}\rangle$ can be approximated by the vacuum far away from the qubit.

Zoom In Zoom Out Reset image size

According to the previous result, we can generate free wave packets, such as input and output states of equations (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 equation (C1), with $| {{\rm{\Omega }}}_{\nu }\rangle =| {{\rm{\Omega }}}_{0}\rangle$. Both wave packets will be Gaussians, equation (6), with mean momentum $\bar{k}$ and width σ. The numerical simulations show that the scattering is elastic for the chosen parameters ($\epsilon =1$, $J=1/\pi$, ${\rm{\Delta }}=\epsilon =1$, and g = 0.3) [35].

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, ${\phi }_{{p}_{1},{p}_{2}}(t)=\langle {{\rm{\Omega }}}_{0}| {a}_{{p}_{1}}{a}_{{p}_{2}}| {\rm{\Psi }}(t)\rangle$. 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

$$\begin{eqnarray}&&F=\int {{\rm{d}}{p}}_{1}{{\rm{d}}{p}}_{2}\,| {\phi }_{{p}_{1},{p}_{2}}({t}_{+}){| }^{2},\end{eqnarray} \tag{ 26 }$$

with p₁ and p₂ such that ${\omega }_{{p}_{1}}+{\omega }_{{p}_{2}}=2({\omega }_{\bar{k}}\pm {\sigma }_{\omega })$ and ${\omega }_{{p}_{1}},{\omega }_{{p}_{2}}\rlap{/}{\in }({\omega }_{\bar{k}}-{\sigma }_{\omega },{\omega }_{\bar{k}}+{\sigma }_{\omega })$, being ${\sigma }_{\omega }$ the width of the input wave packets in energy space. Figure 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 ${\phi }_{{p}_{1},{p}_{2}}({t}_{+})\ne 0$ even though $| {p}_{1}| \ne \bar{k}$ or $| {p}_{2}| \ne \bar{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 [60]. 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 ${\phi }_{{p}_{1},{p}_{2}}({t}_{+})$ goes to zero if $| {p}_{1}| \ne \bar{k}$ or $| {p}_{2}| \ne \bar{k}$. All this is a consequence of the cluster decomposition, see equation (22) and theorem 4 in appendix C.

Zoom In Zoom Out Reset image size

We set ${\omega }_{k}=c| k|$ in H₀. The scatterer and interaction are described by

$$\begin{eqnarray}&&{H}_{\mathrm{sc}}=\displaystyle \sum _{\nu =0}^{M-1}{E}_{\nu }| {{\rm{\Omega }}}_{\nu }\rangle \langle {{\rm{\Omega }}}_{\nu }| +\displaystyle \sum _{J=0}^{M^{\prime} -1}{\tilde{E}}_{J}| J\rangle \langle J| ,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{int}}=\displaystyle \sum _{J=0}^{M^{\prime} -1}\displaystyle \sum _{\nu =0}^{M-1}{g}_{J,\nu }(| J\rangle \langle {{\rm{\Omega }}}_{\nu }| {a}_{0}+{\rm{H}}.{\rm{c}}.),\end{eqnarray} \tag{ 28 }$$

where $\{| {{\rm{\Omega }}}_{\nu }\rangle \}$ and $\{| J\rangle \}$ are the ground and decaying states of the scatterer, respectively, $\{{E}_{\nu }\}$ and $\{{\tilde{E}}_{J}\}$ are their energies, M and $M^{\prime}$ is the number of ground and excited states, respectively, and ${g}_{J,\nu }$ is the coupling strength corresponding to the transition $| {{\rm{\Omega }}}_{\nu }\rangle \leftrightarrow | J\rangle$ (see figure 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^{\prime} =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\gt 0$, which simplifies the final expressions, so we can start from equation (24). Before writing down the two-photon S⁰-matrix in momentum space, we need the one-photon scattering matrix. Imposing energy conservation, it has to be

$$\begin{eqnarray}&&{({S}_{{pk}})}_{\mu \nu }={t}_{\mu \nu }(k)\delta (p+{E}_{\mu }-k-{E}_{\nu }),\end{eqnarray} \tag{ 29 }$$

with k and p the incident and outgoing momenta, respectively, and $| {{\rm{\Omega }}}_{\nu }\rangle$ and $| {{\rm{\Omega }}}_{\mu }\rangle$ the initial and final ground states. The factor ${t}_{\mu \nu }(k)$ is the so-called transmission amplitude. The Dirac delta guarantees energy conservation. Then, the two-photon S⁰-matrix, equation (24) in momentum space is

$$\begin{eqnarray}{({S}_{{\bf{pk}}}^{0})}_{\mu \nu } & = & \displaystyle \frac{1}{{(2\pi )}^{2}}\iint {({S}_{{\bf{yx}}}^{0})}_{\mu \nu }{{\rm{e}}}^{-{\rm{i}}{{\bf{p}}}^{T}{\bf{y}}+{\rm{i}}{{\bf{x}}}^{T}{\bf{k}}}{{\rm{d}}}^{2}{\bf{y}}{{\rm{d}}}^{2}{\bf{x}}\\ & = & \displaystyle \frac{{\rm{i}}}{2\pi }\displaystyle \sum _{n,m=1}^{2}\displaystyle \sum _{\lambda =0}^{M-1}\displaystyle \frac{{t}_{\mu \lambda }({k}_{n}){t}_{\lambda \nu }({k}_{n^{\prime} })}{{p}_{m}+{E}_{\mu }-{k}_{n}-{E}_{\lambda }+{\rm{i}}{0}^{+}}\,\delta ({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{1}-{k}_{2}-{E}_{\nu }).\end{eqnarray} \tag{ 30 }$$

Here, ${n}^{{\prime} }\ne n$, e.g., ${n}^{{\prime} }=2$ if n = 1. The computation is detailed in appendix F. This structure has recently been found by Xu and Fan for a Λ atom (M = 2, $M^{\prime} =1$) within the RWA and Markovian approximations [38]. At first sight (30) may look striking. The matrix S⁰ is not the product of two Dirac-delta functions conserving the single-photon energy, as discussed in section 3.2. The mathematical origin of the structure can be traced back to its form in position space, equation (24). The Heaviside functions set the order in which the different wave packets impinge on the scatterer. The product of Dirac-delta functions is recovered if M = 1 (see appendix F). Besides, equation (30) is also remarkable because presents the generalization of the cluster decomposition for the S-matrix (see equations (16) and (17)) when inelastic processes occur in the scattering.

Zoom In Zoom Out Reset image size

A consequence of (30) is that S⁰ contributes to the fluorescence F, equation (26). This seems to contradict our previous arguments, since S⁰ is built from causally disconnected one-photon events (they do not overlap in the scatterer). To solve the apparent paradox we recall that (30) is a matrix element in momentum space (delocalized photons). For wave packets (5), the scattering amplitude is the integral of these wave packets with (30). 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⁰ and T (see equation (16)). Energy conservation imposes that ${({T}_{{p}_{1}{p}_{2}{k}_{1}{k}_{2}})}_{\mu \nu }={({C}_{{p}_{1}{p}_{2}{k}_{1}{k}_{2}})}_{\mu \nu }\,\delta ({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{1}-{k}_{2}-{E}_{\nu })$. 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}_{{p}_{1}{p}_{2}{k}_{1}{k}_{2}})}_{\mu \nu }$ has simple poles with imaginary parts $\{{\gamma }_{n}^{C}\}$. Similarly, we expect that divergences of ${t}_{\mu \nu }(k)$ come from simple poles with imaginary parts $\{{\gamma }_{n}^{t}\}$. As far as we know, this structure has been found for all S-matrices in waveguide QED [27, 38, 61, 62].

Let us write down the input state in momentum space

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{in}}\rangle =\int {{\rm{d}}{k}}_{1}{{\rm{d}}{k}}_{2}\,{\phi }_{1}({k}_{1}){\phi }_{2}({k}_{2}){{\rm{e}}}^{{{\rm{i}}{k}}_{2}l}{a}_{{k}_{1}}^{\dagger }{a}_{{k}_{2}}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle .\end{eqnarray} \tag{ 31 }$$

The functions ${\phi }_{1}(k)$ and ${\phi }_{2}(k)$ are localized far away the scattering region in position space. The exponential factor ${{\rm{e}}}^{{{\rm{i}}{k}}_{2}l}$ ensures the separation between both wave packets is l. The output state reads

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}\rangle =S| {{\rm{\Psi }}}_{\mathrm{in}}\rangle =\displaystyle \sum _{\mu }\int {{\rm{d}}{p}}_{1}{{\rm{d}}{p}}_{2}\,{\phi }_{\mu }^{\mathrm{out}}({p}_{1},{p}_{2}){a}_{{p}_{1}}^{\dagger }{a}_{{p}_{2}}^{\dagger }| {{\rm{\Omega }}}_{\mu }\rangle \end{eqnarray} \tag{ 32 }$$

with the two-photon wave function ${\phi }_{\mu }^{\mathrm{out}}({p}_{1},{p}_{2})$

$$\begin{eqnarray}&&{\phi }_{\mu }^{\mathrm{out}}({p}_{1},{p}_{2})\propto \displaystyle \sum _{n=1}^{2}\displaystyle \sum _{m=1}^{2}\displaystyle \int {{\rm{d}}{k}}_{n}\left(\displaystyle \frac{{\rm{i}}}{2\pi }\displaystyle \sum _{\lambda }\displaystyle \frac{{t}_{\mu \lambda }({k}_{n}){t}_{\lambda \nu }({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{n}-{E}_{\nu })}{{p}_{m}+{E}_{\mu }-{k}_{n}-{E}_{\lambda }+{\rm{i}}{0}^{+}}\right.+{\rm{i}}{({\tilde{C}}_{{p}_{1}{p}_{2}{k}_{n}})}_{\mu \nu }\phantom{\rule[-0000]{}{3.7ex}})\\ &&\quad \times \,({\phi }_{1}({k}_{n}){{\rm{e}}}^{{\rm{i}}({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{n}-{E}_{\nu })l}{\phi }_{2}({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{n}-{E}_{\nu })+{\phi }_{1}({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{n}-{E}_{\nu }){{\rm{e}}}^{{{\rm{i}}{k}}_{n}l}{\phi }_{2}({k}_{n})),\end{eqnarray} \tag{ 33 }$$

being ${({\tilde{C}}_{{p}_{1}{p}_{2}{k}_{n}})}_{\mu \nu }=\int {{\rm{d}}{k}}_{\bar{n}}{({C}_{{p}_{1}{p}_{2}{k}_{n}{k}_{\bar{n}}})}_{\mu \nu }\delta ({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{n}-{k}_{\bar{n}}-{E}_{\nu })$, with $\bar{n}\ne n$. Even though this expression is cumbersome, we can clearly identify the contribution of S⁰ and T. We solve this integral by means of the residue theorem. Each pole ${\gamma }_{n}^{t}$ and ${\gamma }_{n}^{C}$, together with the exponentials ${{\rm{e}}}^{{{\rm{i}}{k}}_{n}l}$ and ${{\rm{e}}}^{{\rm{i}}({p}_{1}+{p}_{2}+{E}_{\mu }-{k}_{n}-{E}_{\lambda })l}$, gives an exponentially decaying term, ${{\rm{e}}}^{-| {\gamma }_{n}^{t}| l}$ or ${{\rm{e}}}^{-| {\gamma }_{n}^{C}| l}$. We choose Lorentzian envelopes for the wave packets. They have a pole at $\bar{k}-{\rm{i}}\sigma$ (see equation (7)). In consequence, the wave packets will give a term proportional to ${{\rm{e}}}^{-\sigma l}$. Lastly, the imaginary part of the pole of the first term vanishes, $\sim {\rm{i}}{0}^{+}$, so it gives a nondecaying term, ${{\rm{e}}}^{-{0}^{+}l}=1$. The real part of this denominator imposes the single-photon-energy conservation. Thus, it results in the amplitude for the single-photon events, ${\sum }_{\lambda }{A}_{1,\nu \to \lambda }{A}_{2,\lambda \to \mu }$. Therefore, nor S⁰ neither T contains fluorescent terms as the separation between the wave packets grows. The technical details are in appendix G.

As a final application, one can find experimentally the poles of the one- and two-photon scattering matrices $\{{\gamma }_{n}^{t}\}$ and $\{{\gamma }_{n}^{C}\}$ by measuring the decay of F with the distance.

We first prove causal relations in a free theory. In order to do so, we work with localized wave packets ${\psi }_{\bar{k}\bar{x}}(t)$, equation (5). Actual calculations are done with Gaussian wave packets, equation (6). The following two lemmas are used in the demonstration of the theorem.

Lemma 2. Let the dispersion relation ${\omega }_{k}$ have an upper bounded group velocity ${v}_{k}={\partial }_{k}{\omega }_{k}$:

$$\begin{eqnarray}&&| {v}_{k}| \leqslant c.\end{eqnarray} \tag{ B1 }$$

Then, the function $f(k)={kx}-{\omega }_{k}t$ only has stationary points if the distance to the light cone is nonnegative. In other words

$$\begin{eqnarray}&&{d}_{c}(x,t)=| x| -c| t| \gt 0\iff | {f}^{{\prime} }(k)| \gt 0,\forall k.\end{eqnarray} \tag{ B2 }$$

Proof. Solving the equation $f^{\prime} (k)=x-{\partial }_{k}{\omega }_{k}t\,=\,0$ leads to the condition $\tfrac{x}{t}={v}_{k}$ or $| x/t| =| {v}_{k}| \leqslant c$. Then, provided $f^{\prime} (k)=0$, it follows $| x| \leqslant c| t| \Rightarrow {d}_{c}(x,t)\leqslant 0$, which shows (B2). □

Lemma 3. Assume that ${\omega }_{k}$ is n-times differentiable and that every derivative $| {\omega }_{k}^{(r\leqslant n)}|$ is upper bounded by an mth order polynomial in $| k|$. Then the following integral bound applies

$$\begin{eqnarray*}&&\left|\int {{\rm{e}}}^{{\rm{i}}{k}{x}-\tfrac{1}{{\sigma }^{2}}{(k-{k}_{0})}^{2}-{\rm{i}}{\omega }_{k}t}p(k){\rm{d}}{k}\right|=\max ({\sigma }^{m+n+r},1)\max ({t}^{n},1){ \mathcal O }\left(\displaystyle \frac{1}{| x{| }^{n}}\right),\end{eqnarray*}$$

where p(k) is a polynomial of degree $r.$

Proof. Result 5.1 from [64] states that the integral $I(x)={\int }_{a}^{b}{{\rm{e}}}^{{\rm{i}}{k}{x}}q(k){\rm{d}}{k}$ may be integrated by parts n times, obtaining

$$\begin{eqnarray}&&I(x)=\displaystyle \sum _{s=0}^{n-1}{\left(\displaystyle \frac{{\rm{i}}}{x}\right)}^{s+1}[{{\rm{e}}}^{{\rm{i}}{a}{x}}{q}^{(s)}(a)-{{\rm{e}}}^{{\rm{i}}{b}{x}}{q}^{(s)}(b)]+{\epsilon }_{n}(x),\end{eqnarray} \tag{ B3 }$$

where the error term satisfies

$$\begin{eqnarray}&&{\epsilon }_{n}(x)={\left(\displaystyle \frac{{\rm{i}}}{x}\right)}^{n}\int {{\rm{e}}}^{{\rm{i}}{k}{x}}{q}^{(n)}(k){\rm{d}}{k}=o({x}^{-n})\end{eqnarray} \tag{ B4 }$$

provided that q(k) is n-times differentiable and that ${q}^{(n)}\in {L}^{1}$. Based on the conditions of the lemma, this is satisfied since $q(k)={{\rm{e}}}^{-\tfrac{1}{{\sigma }^{2}}{(k-{k}_{0})}^{2}-{\rm{i}}{\omega }_{k}t}p(k)$. The limits of the integral may be easily extended to $\pm \infty ,$ as explained in result 5.2 from [64]. Since ${x}^{-s}{q}^{(s)}(a)\to 0$ when $a\to \pm \infty ,\forall x,$ we obtain

$$\begin{eqnarray}&&I(x)=\int {{\rm{e}}}^{{\rm{i}}{k}{x}}q(k)\,{\rm{d}}{k}={\left(\displaystyle \frac{{\rm{i}}}{x}\right)}^{n}\int {{\rm{e}}}^{{\rm{i}}{k}{x}}{q}^{(n)}(k){\rm{d}}{k}.\end{eqnarray} \tag{ B5 }$$

Moreover, ${q}^{(n)}$, resulting from a product of derivatives of ${\omega }_{k}t,$ $-{k}^{2}/{\sigma }^{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 ${({\partial }_{k}{\omega }_{k})}^{n}$, and the monomials in $| k|$, which produce another prefactor ${\sigma }^{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)\to 0\to 0$ when $a\to \pm \infty ,\forall 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}=\int {\rm{d}}{k}\,{\omega }_{k}{a}_{k}^{\dagger }{a}_{k}$. Let ${\psi }_{\bar{k}\bar{x}}(t)$ and ${\psi }_{\bar{p}\bar{y}}(t^{\prime} )$ 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}| =| {\partial }_{k}{\omega }_{k}| \leqslant 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:

$$\begin{eqnarray}&&| {\partial }_{k}^{(r\leqslant n)}{\omega }_{k}| \leqslant {a}_{r}+{(| k| /{b}_{r})}^{m},\,0\lt {a}_{r},{b}_{r}\lt +\infty .\end{eqnarray} \tag{ B6 }$$

The commutator between these wave packets is small whenever they are outside of their respective light cones, that is, whenever $d=| \bar{y}-\bar{x}| -c| t^{\prime} -t| \gg 0$,

$$\begin{eqnarray}&&\parallel [{\psi }_{\bar{k}\bar{x}}(t),{\psi }_{\bar{p}\bar{y}}{(t^{\prime} )}^{\dagger }]\parallel ={ \mathcal O }\left(\displaystyle \frac{1}{| d{| }^{n}}\right),\,d\to \infty .\end{eqnarray} \tag{ B7 }$$

Proof. Let us assume that the model evolves freely according to the free Hamiltonian ${H}_{0}=\int {\rm{d}}{k}\,{\omega }_{k}{a}_{k}^{\dagger }{a}_{k}$. In this case, our wave packet operators have the simple form

$$\begin{eqnarray}&&{\psi }_{\bar{k}\bar{x}}(t)=\int {{\rm{e}}}^{{\rm{i}}{k}\bar{x}-{\rm{i}}{\omega }_{k}t}{\phi }_{\bar{k}}{(k)}^{* }{a}_{k}(0){\rm{d}}{k},\end{eqnarray} \tag{ B8 }$$

and analogously for ${\psi }_{\bar{p}\bar{y}}(t^{\prime} )$. The commutator between operators reads

$$\begin{eqnarray}&&I:= [{\psi }_{\bar{k}\bar{x}}(t),{\psi }_{\bar{p}\bar{y}}{(t^{\prime} )}^{\dagger }]=\int {{\rm{e}}}^{{\rm{i}}{k}(\bar{x}-\bar{y})-{\rm{i}}{\omega }_{k}(t-t^{\prime} )}{\phi }_{\bar{k}}(k){\phi }_{\bar{p}}{(k)}^{* }{\rm{d}}{k}.\end{eqnarray} \tag{ B9 }$$

Let $d={d}_{c}(\bar{x}-\bar{y},t-t^{\prime} )=| \bar{x}-\bar{y}| -c| t-t^{\prime} | \gt 0$, using lemma 2 we know that the exponent has no stationary point. Assuming w.l.o.g. $\bar{x}\gt \bar{y}$, $t\gt {t}^{{\prime} }$ (other combinations are analogous) and writing ${\tilde{\omega }}_{k}={\omega }_{k}-{ck},$ we obtain

$$\begin{eqnarray*}&&I=\displaystyle \int {{\rm{e}}}^{{\rm{i}}{k}(\bar{x}-\bar{y})-{\rm{i}}{\omega }_{k}(t-t^{\prime} )}{\phi }_{\bar{k}}{(k)}^{* }{\phi }_{\bar{p}}(k){\rm{d}}{k}=\displaystyle \int {{\rm{e}}}^{{\rm{i}}{k}{{d}}_{c}(\bar{x}-\bar{y},t-t^{\prime} )-{\rm{i}}{\tilde{\omega }}_{k}(t-t^{\prime} )}{\phi }_{\bar{k}}{(k)}^{* }{\phi }_{\bar{p}}(k){\rm{d}}{k}.\end{eqnarray*}$$

The exponent ${\tilde{\omega }}_{k}={\omega }_{k}-{ck}$ is n-times differentiable and is upper bounded in modulus by a polynomial of degree $m\geqslant 1$. Lemma 3 therefore allows us to bound the commutator by a term ${ \mathcal O }({d}^{-n})$. □

Note that for a linear dispersion, ${\omega }_{k}=c| k|$, we can rewrite this integral as a function of the distance between world lines from equation (B2), $d=(\bar{x}-\bar{y})-c(t-t^{\prime} )$. Introducing ${k}_{\pm }=(\bar{k}\pm \bar{p})/2$ and using our Gaussian wave packets (6), we obtain

$$\begin{eqnarray}&&| I| =\exp \left[-\displaystyle \frac{{k}_{-}^{2}}{{\sigma }^{2}}-\displaystyle \frac{{d}^{2}{\sigma }^{2}}{4}\right].\end{eqnarray} \tag{ B10 }$$

This bound is better than the one we have found but it is compatible with lemma 3 and theorem 2.

Causal relation (B7) can be extended to the full model (2).

Theorem 3. Let H be the light–matter Hamiltonian given by equation (2). We assume the conditions of theorem 2: differentiable, polynomially bounded functions ${\omega }_{k}$ and g_k, with degrees $n\geqslant 2$. Then, all wave packets outside the light cone of the scatterer evolve approximately with the free Hamiltonian, H₀. More precisely, if $(\bar{x},{t}_{1})$ and $(\bar{x},{t}_{0})$ are two points outside the light cone

$$\begin{eqnarray}&&{\psi }_{\bar{k}\bar{x}}({t}_{1})={U}_{0}{({t}_{1},{t}_{0})}^{\dagger }{\psi }_{\bar{k}\bar{x}}({t}_{0}){U}_{0}({t}_{1},{t}_{0})+{ \mathcal O }\left(\displaystyle \frac{1}{| {d}_{\min }{| }^{n-1}}\right),\end{eqnarray} \tag{ B11 }$$

where ${d}_{\min }=\min \{d(\bar{x},{t}_{1}),d(\bar{x},{t}_{0})\}\gg 0$ and

$$\begin{eqnarray}&&{U}_{0}(t,{t}_{0})=\exp (-{\rm{i}}(t-{t}_{0}){H}_{0})\end{eqnarray} \tag{ B12 }$$

is the free-evolution operator for the photons at time t₀.

Proof. We start by building the Heisenberg equations for the operators

$$\begin{eqnarray}&&{\partial }_{t}{a}_{k}(t)=-{\rm{i}}{\omega }_{k}{a}_{k}(t)-{{\rm{i}}{g}}_{k}G(t).\end{eqnarray} \tag{ B13 }$$

Making the change of variables ${a}_{k}(t)={{\rm{e}}}^{-{\rm{i}}{\omega }_{k}t}{b}_{k}(t)$, we have

$$\begin{eqnarray}&&{\partial }_{t}{b}_{k}(t)=-{{\rm{i}}{g}}_{k}G(t){{\rm{e}}}^{{\rm{i}}{\omega }_{k}t},\end{eqnarray} \tag{ B14 }$$

so that the wave packet operators evolved from some initial time t_s are

$$\begin{eqnarray}&&{\psi }_{\bar{k}\bar{x}}(t)=\int {{\rm{e}}}^{{\rm{i}}{k}\bar{x}-{\rm{i}}{\omega }_{k}t}\left[{b}_{k}({t}_{s})-{\rm{i}}{\int }_{{t}_{s}}^{t}{g}_{k}G(\tau ){{\rm{e}}}^{+{\rm{i}}{\omega }_{k}\tau }{\rm{d}}\tau \right]{\phi }_{\bar{k}}(k){\rm{d}}{k}\end{eqnarray} \tag{ B15 }$$

$$\begin{eqnarray}&&=\,{U}_{0}(t,{t}_{s}){\psi }_{\bar{k}\bar{x}}({t}_{s}){U}_{0}{(t,{t}_{s})}^{\dagger }-{\rm{i}}{\int }_{{t}_{s}}^{t}\left[\int {{\rm{e}}}^{{\rm{i}}{k}\bar{x}-{\rm{i}}{c}(t-\tau )}{g}_{k}{\phi }_{\bar{k}}(k){\rm{d}}{k}\right]\,G(\tau ){\rm{d}}\tau \end{eqnarray} \tag{ B16 }$$

$$\begin{eqnarray}&&=\,{U}_{0}(t,{t}_{s}){\psi }_{\bar{k}\bar{x}}({t}_{s}){U}_{0}{(t,{t}_{s})}^{\dagger }-{\rm{i}}{\int }_{0}^{t-{t}_{s}}\left[\int {{\rm{e}}}^{{\rm{i}}{k}\bar{x}-{\rm{i}}{c}\tau ^{\prime} }{g}_{k}{\phi }_{\bar{k}}(k){\rm{d}}{k}\right]\,G(\tau ){\rm{d}}\tau ^{\prime} .\end{eqnarray} \tag{ B17 }$$

The first part corresponds to free evolution, while the second part is an error term $\varepsilon (t)$, which can be bounded. We will assume without loss of generality $\parallel G\parallel =1$, with $| | \cdot | |$ the Hilbert–Schmidt norm, and $| {t}_{1}| \gt | {t}_{0}|$. We have to choose the integration limits t and t_s so that $\mathrm{sign}(\tau ^{\prime} )=\mathrm{sign}(x)$. If $x\gt 0$ then ${t}_{1}\gt {t}_{0}\gt 0$ and $(t,{t}_{s})=({t}_{1},{t}_{0})$ is a good choice. If $x\lt 0$ then $0\gt {t}_{0}\gt {t}_{1}$ and again $(t,{t}_{s})=({t}_{1},{t}_{0})$ is also a valid choice ($\tau ^{\prime} \lt 0$). This means we can introduce $\tau ^{\prime\prime} =\mathrm{sign}(x)\tau ^{\prime} \geqslant 0$ and bound

$$\begin{eqnarray}&&| \varepsilon ({t}_{1})| \leqslant {\int }_{0}^{| {t}_{1}-{t}_{0}| }\left|\int {{\rm{e}}}^{{\rm{i}}\mathrm{sign}(\bar{{x}}){{kd}}_{c}(| \bar{{x}}| ,\tau ^{\prime\prime} )}q(k){\rm{d}}{k}\right|{\rm{d}}\tau ^{\prime\prime} \leqslant {\int }_{0}^{| {t}_{1}-{t}_{0}| }{ \mathcal O }\left(\displaystyle \frac{1}{{{d}}_{c}{(| \bar{{x}}| ,\tau ^{\prime\prime} )}^{n}}\right){\rm{d}}\tau ^{\prime\prime} \end{eqnarray} \tag{ B18 }$$

$$\begin{eqnarray}&&\leqslant { \mathcal O }\left({\left.\displaystyle \frac{1}{c(n-1)}\displaystyle \frac{1}{{(| \bar{x}| -c\tau )}^{n-1}}\right|}_{\tau =0}^{\tau =| {t}_{1}-{t}_{0}| }\right)\ \leqslant { \mathcal O }\left(\displaystyle \frac{1}{{d}_{c}{(| \bar{x}| ,| {t}_{1}-{t}_{0}| )}^{n-1}}\right).\end{eqnarray} \tag{ B19 }$$

Here we have taken into account that ${d}_{c}(| \bar{x}| ,\tau ^{\prime\prime} )\geqslant {d}_{c}(| \bar{x}| ,| {t}_{1}-{t}_{0}| )\gt 0$ in the domain of integration. We can now use the fact that ${d}_{c}(| \bar{x}| ,| {t}_{1}-{t}_{0}| )\geqslant {d}_{c}(| \bar{x}| ,| {t}_{1}| )\geqslant \min \{{d}_{c}(\bar{x},{t}_{1}),{d}_{c}(\bar{x},{t}_{0})\}$, obtaining the expression in the theorem. □

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 section 2.3 and equation (9)). However, using theorems 1 and 3 we have that, given a state $| {\rm{\Psi }}\rangle \equiv {\psi }_{\bar{k},\bar{x}}{({t}_{0})}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle$, then

$$\begin{eqnarray}U({t}_{\pm })| {\rm{\Psi }}\rangle & = & U({t}_{\pm }){\psi }_{\bar{k},\bar{x}}({t}_{0})U{({t}_{\pm })}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle \\ & = & {U}_{0}({t}_{\pm }){\psi }_{\bar{k},\bar{x}}({t}_{0}){U}_{0}{({t}_{\pm })}^{\dagger }| {{\rm{\Omega }}}_{\nu }\rangle \equiv {U}_{0}({t}_{\pm })| {{\rm{\Psi }}}_{\mathrm{in}}\rangle .\end{eqnarray} \tag{ B20 }$$

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}_{\pm }\to \pm \infty$) from the scatterer and, by means of theorem 1 and the conditions presented in 2.4 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.

We first need to compute the one photon amplitude. Let the one-photon input state be,

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{in}}^{1}\rangle ={\psi }_{{\bar{k}}_{1},{\bar{x}}_{1}}^{\mathrm{in}\,\dagger }| {{\rm{\Omega }}}_{\nu }\rangle ,\end{eqnarray} \tag{ D1 }$$

with the creation operator ${\psi }_{{\bar{k}}_{1},{\bar{x}}_{1}}^{\mathrm{in}\,\dagger }$ given by equation (5), removing the time dependence. For simplicity, we absorb the factor ${{\rm{e}}}^{{\rm{i}}{k}{\bar{x}}_{1}}$ into the wave packet: ${\phi }_{{\bar{k}}_{1},{\bar{x}}_{1}}(k)={{\rm{e}}}^{{\rm{i}}{k}{\bar{x}}_{1}}{\phi }_{{\bar{k}}_{1}}(k)$. In position space, the output state will read

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}^{1}\rangle =S| {{\rm{\Psi }}}_{\mathrm{in}}^{1}\rangle =\displaystyle \sum _{\mu =1}^{M}\int {\rm{d}}{y}{\rm{d}}{x}\,{({S}_{{yx}})}_{\mu \nu }{\phi }_{{\bar{k}}_{1},{\bar{x}}_{1}}(x)| y,{{\rm{\Omega }}}_{\mu }\rangle .\end{eqnarray} \tag{ D2 }$$

Defining

$$\begin{eqnarray}&&{\phi }_{1,\mu \nu }(y)=\int {\rm{d}}{x}\,{({S}_{{yx}})}_{\mu \nu }{\phi }_{{\bar{k}}_{1},{\bar{x}}_{1}}(x)\end{eqnarray} \tag{ D3 }$$

and

$$\begin{eqnarray}&&| {\xi }_{\mathrm{out}}^{1}{\rangle }_{1,\mu \nu }=\int {\rm{d}}{y}\,{\phi }_{1,\mu \nu }(y)| y;{{\rm{\Omega }}}_{\mu }\rangle ,\end{eqnarray} \tag{ D4 }$$

being $| y;{{\rm{\Omega }}}_{\mu }\rangle ={a}_{y}^{\dagger }| {{\rm{\Omega }}}_{\mu }\rangle$ the state with a photon at y and the scatterer in the ground state $| {{\rm{\Omega }}}_{\mu }\rangle$, the output state (D2) can be rewritten as

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}^{1}\rangle =\displaystyle \sum _{\mu =1}^{M}| {\xi }_{\mathrm{out}}^{1}{\rangle }_{1,\mu \nu }.\end{eqnarray} \tag{ D5 }$$

The probability amplitude will read

$$\begin{eqnarray}&&{A}_{1,\nu \to \mu }=\langle {{\rm{\Omega }}}_{\mu }| {\psi }_{{\bar{p}}_{1},{\bar{y}}_{1}}^{\mathrm{out}}\,S\,{\psi }_{{\bar{k}}_{1},{\bar{x}}_{1}}^{\mathrm{in}\,\dagger }| {{\rm{\Omega }}}_{\nu }\rangle =\int {\rm{d}}{y}\,{\phi }_{{\bar{p}}_{1},{\bar{y}}_{1}}{(y)}^{* }{\phi }_{1,\mu \nu }(y).\end{eqnarray} \tag{ D6 }$$

If the wave packets are monochromatic with momenta k₁ and p₁, respectively, this amplitude is

$$\begin{eqnarray}&&{A}_{1,\nu \to \mu }={({S}_{{p}_{1}{k}_{1}})}_{\mu \nu }.\end{eqnarray} \tag{ D7 }$$

The two-photon wave packet, as sketched in figure 1, is

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{in}}^{2}\rangle ={\psi }_{{\bar{k}}_{1},{\bar{x}}_{1}}^{\mathrm{in}\,\dagger }{\psi }_{{\bar{k}}_{2},{\bar{x}}_{2}}^{\mathrm{in}\,\dagger }| {{\rm{\Omega }}}_{\nu }\rangle .\end{eqnarray} \tag{ D8 }$$

By definition, the output state is

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}^{2}\rangle =S| {{\rm{\Psi }}}_{\mathrm{in}}^{2}\rangle .\end{eqnarray} \tag{ D9 }$$

Here, we are interested in he limit of well separated incident photons. Thus, only the linear part of the scattering matrix S⁰ is considered. We introduce the identity operator

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}^{2}\rangle ={\mathbb{I}}S{\mathbb{I}}| {{\rm{\Psi }}}_{\mathrm{in}}^{2}\rangle ,\end{eqnarray} \tag{ D10 }$$

with

$$\begin{eqnarray}&&{\mathbb{I}}=\displaystyle \frac{1}{2}\displaystyle \sum _{\mu =1}^{M}\int {{\rm{d}}{x}}_{1}{{\rm{d}}{x}}_{2}\,| {x}_{1}{x}_{2};{{\rm{\Omega }}}_{\nu }\rangle \langle {x}_{1}{x}_{2};{{\rm{\Omega }}}_{\nu }| ,\end{eqnarray} \tag{ D11 }$$

being $| {x}_{1}{x}_{2};{{\rm{\Omega }}}_{\nu }\rangle ={a}_{{x}_{1}}^{\dagger }{a}_{{x}_{2}}^{\dagger }| {{\rm{\Omega }}}_{\mu }\rangle$ the symmetrized state with two photons at x₁ at x₂ and the scatterer at $| {{\rm{\Omega }}}_{\nu }\rangle$.

Introducing (D11) in (D10) and considering (D8) and (24) we get

$$\begin{eqnarray}| {{\rm{\Psi }}}_{\mathrm{out}}^{2}\rangle & = & \displaystyle \frac{1}{4}\displaystyle \int {{\rm{d}}{y}}_{1}{{\rm{d}}{y}}_{2}{{\rm{d}}{x}}_{1}{{\rm{d}}{x}}_{2}\displaystyle \sum _{\mu ,\lambda =1}^{M}\displaystyle \sum _{n,m=1}^{2}{({S}_{{y}_{n}{x}_{m}})}_{\mu \lambda }{({S}_{{y}_{{n}^{{\prime} }}{x}_{{m}^{{\prime} }}})}_{\lambda \nu }\theta ({y}_{{n}^{{\prime} }}-{y}_{n})({\phi }_{{\bar{k}}_{1},{\bar{x}}_{1}}({x}_{1}){\phi }_{{\bar{k}}_{2},{\bar{x}}_{2}}({x}_{2})\\ & & +\,{\phi }_{{\bar{k}}_{1},{\bar{x}}_{1}}({x}_{2}){\phi }_{{\bar{k}}_{2},{\bar{x}}_{2}}({x}_{1}))| {y}_{1}{y}_{2};{{\rm{\Omega }}}_{\mu }\rangle .\end{eqnarray} \tag{ D12 }$$

with ${n}^{{\prime} }\ne n$ and ${m}^{{\prime} }\ne m$. Now, we have to compute integrals as

$$\begin{eqnarray}&&C=\int {{\rm{d}}{x}}_{1}{{\rm{d}}{x}}_{2}\,\displaystyle \sum _{n,m}{({S}_{{y}_{n}{x}_{m}})}_{\mu \lambda }{({S}_{{y}_{{n}^{{\prime} }}{x}_{{m}^{{\prime} }}})}_{\lambda \nu }{\phi }_{{\bar{k}}_{i},{\bar{x}}_{i}}({x}_{1}){\phi }_{{\bar{k}}_{j},{\bar{x}}_{j}}({x}_{2})\theta ({y}_{{n}^{{\prime} }}-{y}_{n}).\end{eqnarray} \tag{ D13 }$$

Using equation (D3)

$$\begin{eqnarray}&&C=\displaystyle \sum _{n=1}^{2}({\phi }_{i,\mu \lambda }({y}_{n}){\phi }_{j,\lambda \nu }({y}_{{n}^{{\prime} }})+{\phi }_{j,\mu \lambda }({y}_{n}){\phi }_{i,\lambda \nu }({y}_{{n}^{{\prime} }}))\theta ({y}_{{n}^{{\prime} }}-{y}_{n}).\end{eqnarray} \tag{ D14 }$$

Following the sketch drawn in figure 1, if ${x}_{m}\lt {x}_{{m}^{{\prime} }}$, then ${\phi }_{1}({x}_{m}){\phi }_{2}({x}_{{m}^{{\prime} }})$ is zero, so ${\phi }_{1,\mu \nu }({y}_{n}){\phi }_{2,\mu \nu }({y}_{{n}^{{\prime} }})$ is zero if ${y}_{n}\lt {y}_{{n}^{{\prime} }}$. Therefore, choosing i = 1 and j = 2, the integral C reads

$$\begin{eqnarray}&&C=\displaystyle \sum _{n=1}^{2}{\phi }_{2,\mu \lambda }({y}_{n}){\phi }_{1,\lambda \nu }({y}_{{n}^{{\prime} }}).\end{eqnarray} \tag{ D15 }$$

One can easily show that the same expression holds if we take i = 2 and j = 1. The output state, equation (D12), then reads

$$\begin{eqnarray}&&| {{\rm{\Psi }}}_{\mathrm{out}}^{2}\rangle =\displaystyle \frac{1}{2}\int {{\rm{d}}{y}}_{1}{{\rm{d}}{y}}_{2}\,\displaystyle \sum _{\mu ,\lambda =1}^{M}({\phi }_{2,\mu \lambda }({y}_{1}){\phi }_{1,\lambda \nu }({y}_{2})+{\phi }_{2,\mu \lambda }({y}_{2}){\phi }_{1,\lambda \nu }({y}_{1}))| {y}_{1}{y}_{2};{{\rm{\Omega }}}_{\mu }\rangle .\end{eqnarray} \tag{ D16 }$$

Finally, the probability amplitude of going to the output state ${\psi }_{{\bar{p}}_{1},{\bar{y}}_{1}}^{\mathrm{out}\,\dagger }{\psi }_{{\bar{p}}_{2},{\bar{y}}_{2}}^{\mathrm{out}\,\dagger }| {{\rm{\Omega }}}_{\mu }\rangle$ will be the overlap between this state and (D16). Using (D6)

$$\begin{eqnarray}&&{A}_{\mathrm{in}\to \mathrm{out}}=\langle {{\rm{\Omega }}}_{\mu }| {\psi }_{{\bar{p}}_{1},{\bar{y}}_{1}}^{\mathrm{out}}{\psi }_{{\bar{p}}_{2},{\bar{y}}_{2}}^{\mathrm{out}}S{\psi }_{{\bar{k}}_{1},{\bar{x}}_{1}}^{\mathrm{in}\,\dagger }{\psi }_{{\bar{k}}_{2},{\bar{x}}_{2}}^{\mathrm{out}\,\dagger }| {{\rm{\Omega }}}_{\nu }\rangle =\displaystyle \sum _{\lambda =0}^{M-1}{A}_{1,\nu \to \lambda }{A}_{2,\lambda \to \mu },\end{eqnarray} \tag{ D17 }$$

as expected. In the calculations, we have set $\langle {{\rm{\Omega }}}_{\mu }| {\psi }_{{\bar{p}}_{i},{\bar{y}}_{i}}^{\mathrm{out}}\,S\,{\psi }_{{\bar{k}}_{j}{\bar{x}}_{j}}^{\mathrm{in}\,\dagger }| {{\rm{\Omega }}}_{\nu }\rangle =0$ for $i\ne j$, since we assume that both incident wave packets are far away.

A final comment is in order. Without the step functions in (24), the unphysical amplitude ${A}_{2,\nu \to \lambda }{A}_{1,\lambda \to \mu }$ would appear in the final probability amplitude.