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

TL;DR

This paper explores how approximate causality emerges in waveguide QED systems and influences the structure of the N-photon scattering matrix, with implications for understanding photon interactions and correlations.

Contribution

It proves a Lieb-Robinson-type bound for continuous models and characterizes the N-photon scattering matrix structure in waveguide QED.

Findings

Verified cluster decomposition in two-photon scattering

Analyzed ground state and correlations in ultrastrong coupling regime

Derived S-matrix structure for inelastic and Raman scattering

Abstract

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.

Figures

Equations

$$[A(x,t),B(y,t^{\prime })]=0\text{if}|x-y|-c|t-t^{\prime }|>0,$$

$$H=H_0+H_{\text{sc}}+\int (g_k{G}^{†}{a}_k+g_k^{\ast }G{a}_k^{†})dk.$$

$$H_0=\int {\omega }_k{a}_k^{†}{a}_{k}dk,$$

$$g_x=\frac{1}{\sqrt{2\pi }}\int dk{e}^{ikx}{g}_k$$

$${\psi }_{\bar{k}\bar{x}}(t-t_0{)}^{†}=\int e^{ik\bar{x}-i{\omega }_k(t-t_0)}{\varphi }_{\bar{k}}(k)a_k^{†}dk.$$

$${\varphi }_{\bar{k}}(k)=\frac{1}{\sqrt[4]{2\pi }\sqrt{\sigma }}\exp \left[-(k-\bar{k}{)}^2/4{\sigma }^2\right],$$

$${\varphi }_{\bar{k}}(k)=\sqrt{\frac{\sigma }{\pi }}\frac{1}{k-\bar{k}+i\sigma }.$$

$$\underset{\bar{x}\to \pm \infty }{\lim }⟨{\Omega }_{\mu }⟩O(\bar{x},\Delta ){\Omega }_{\mu }=⟨\mathrm{vac}⟩O(\bar{x},\Delta )\mathrm{vac},$$

$$\Vert U(t)|\Psi ⟩-U^0(t)|{\Psi }_{in/out}⟩\Vert \stackrel{t\to \mp \infty }{⟶}0,$$

$$|{\Psi }_{\mathrm{out}}⟩=S|{\Psi }_{\mathrm{in}}⟩,$$

$$S=\underset{t_{\pm }\to \pm \infty }{\lim }{U}_I(t_{+},t_{-}).$$

$$A\equiv ⟨{\Omega }_{\mu }|{\psi }^{\mathrm{out}}(t_{+})S{\psi }^{\mathrm{in}}(t_{-}{)}^{†}|{\Omega }_{\nu }⟩$$

$$⟨{\Omega }_0⟩{\psi }_{\bar{k}\bar{x}}^{†}{\psi }_{\bar{k}\bar{x}}{\Omega }_0\le \mathcal{O}(|\bar{x}{|}^{-n}),\text{leavevmode}|\bar{x}|\to \infty$$

$$(S_{\mathbf{pk}}{)}_{\mu \nu }=⟨{\Omega }_{\mu }⟩\prod _i{a}_{p_i}S\prod _j{a}_{k_j}^{†}{\Omega }_{\nu }.$$

$$\Vert [{\psi }_{\bar{k}\bar{x}}(t),{\psi }_{\bar{p}\bar{y}}(t^{\prime }{)}^{†}]\Vert =\mathcal{O}\left(\frac{1}{|D{|}^n}\right)+\mathcal{O}\left(\frac{1}{|D_0{|}^{n-1}}\right),$$

$$S_{\text{p}\text{k}}=S_{\text{p}\text{k}}^0+i{T}_{\text{p}\text{k}}.$$

$$S_{\text{p}\text{k}}^0=\frac{1}{N!}\prod _{n=1}^N{S}_{p_n{k}_n}+\text{permutations}[k_n\leftrightarrow k_m,p_n\leftrightarrow p_m],$$

$$S_{\text{y}\text{x}}^0=\frac{1}{N!}\prod _{n=1^N}{S}_{y_n{x}_n}+\text{permutations}[x_n\leftrightarrow x_m,y_n\leftrightarrow y_m].$$

$$A=⟨{\Omega }_{\mu }|\prod _{m=1}^2{\psi }_{{\bar{p}}_m{\bar{y}}_m}^{\mathrm{out}}(t_{+})\prod _{n=1}^2{\psi }_{{\bar{k}}_n{\bar{x}}_n}^{\mathrm{in}}(t_{-}{)}^{†}|{\Omega }_{\nu }⟩.$$

$$A\simeq ⟨{\Omega }_{\mu }|{\psi }_{{\bar{p}}_2{\bar{y}}_2}^{\mathrm{out}}(t_{+}){\psi }_{{\bar{k}}_2{\bar{y}}_2}^{\mathrm{in}}(t_{-}{)}^{†}{\psi }_{{\bar{p}}_1{\bar{x}}_1}^{\mathrm{out}}(t_{+}){\psi }_{{\bar{k}}_1{\bar{x}}_1}^{\mathrm{in}}(t_{-}{)}^{†}|{\Omega }_{\nu }⟩.$$

$$A_{12}=\sum _{\lambda =0}^{M-1}{A}_{1,\nu \to \lambda }{A}_{2,\lambda \to \mu },$$

$$A=\sum _{{\lambda }_1,\dots ,{\lambda }_{N-1}=0}^{M-1}\prod _{n=1}^N{A}_{n,{\lambda }_{N+1-n}\to {\lambda }_{N-n}},$$

$$A_{n,{\lambda }_{N+1-n}\to {\lambda }_{N-n}}$$

$$=⟨{\Omega }_{{\lambda }_{N-n}}|{\psi }_{{\bar{p}}_n{\bar{y}}_n}^{\mathrm{out}}(t_{+}){\psi }_{{\bar{k}}_n{\bar{x}}_n}^{\mathrm{in}}(t_{-}{)}^{†}|{\Omega }_{{\lambda }_{N+1-n}}⟩.$$

$$(S_{\text{y}\text{x}}^0{)}_{\mu \nu }=\sum _{{\lambda }_1\dots {\lambda }_{N-1}=0}^{M-1}\prod _{n=1}^N(S_{y_n{x}_n}{)}_{{\lambda }_{n-1}{\lambda }_n}\prod _{m=1}^{N-1}\theta (y_{m+1}-y_m)+\text{permutations}[x_n\leftrightarrow x_m,y_n\leftrightarrow y_m],$$

$$H$$

$$+g({\sigma }^{-}+{\sigma }^{+})(a_0+a_0^{†}).$$

$$F=\int d{p}_{1}d{p}_2|{\varphi }_{p_1,p_2}(t_{+}){|}^2,$$

$$H_{\mathrm{sc}}=\sum _{\nu =0}^{M-1}{E}_{\nu }|{\Omega }_{\nu }⟩⟨{\Omega }_{\nu }|+\sum _{J=0}^{M^{\prime }-1}{\tilde{E}}_J|J⟩⟨J|,$$

$$H_{\mathrm{int}}=\sum _{J=0}^{M^{\prime }-1}\sum _{\nu =0}^{M-1}{g}_{J,\nu }(|J⟩⟨{\Omega }_{\nu }|a_0+H.c.),$$

$$(S_{pk}{)}_{\mu \nu }=t_{\mu \nu }(k)\delta (p+E_{\mu }-k-E_{\nu }),$$

$$(S_{\mathbf{pk}}^0{)}_{\mu \nu }$$

Full text

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 Weinberg (1996). In Quantum Field Theory (QFT), strict causality imposes that two operators $A(x,t)$ and $B(y,t^{\prime})$ acting on two space-like separated points $(x,t)$ and $(y,t^{\prime})$, must commute,

[TABLE]

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 Wichmann and Crichton (1963). In fact, all acceptable QFT interactions must result in $S$-matrices fulfilling such a decomposition Weinberg (1995).

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 low-energy 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 Lieb and Robinson (1972) 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 Bravyi et al. (2006), but lead to important results in the static properties of many-body Hamiltonians, such as the clustering of correlations and the area law in gapped models Nachtergaele and Sims (2006); Hastings (2007).

In this work we demonstrate the existence and explore the consequences of emergent causality in the nonrelativistic framework of waveguide QED Lodahl et al. (2015); Zhou et al. (2008a); Zheng et al. (2013); Lu et al. (2014); Roy et al. (2017). Theses systems consist of photons propagating in low-dimensional environments —waveguides, photonic crystals, etc—, interacting with local quantum systems. Such models do not satisfy Lorentz or translational invariance, they are typically dispersive, and the photon-matter interaction may become highly non-perturbative. Experimental implementations include dielectrics Mitsch et al. (2014); Yu et al. (2014), cavity arrays Vučković and Yamamoto (2003), metals Chang et al. (2007), diamond structures Sipahigil et al. (2016); Bhaskar et al. (2016), and superconductors Astafiev et al. (2010); Hoi et al. (2011); Liu and Houck (2017); Forn-Díaz et al. (2017) interacting with atoms, molecules, quantum dots, color centers in diamond or superconducting qubits. The focus of waveguide QED is set on quantum processes involving few photons and scatterers. In this regard, it is not surprising that there exists an extensive theoretical literature for waveguide-QED systems Roy et al. (2017), which develops a variety of analytical and numerical methods for the study of the $N$-photon $S$-matrix Shen and Fan (2005); Zhou et al. (2008b); Fan et al. (2010); Roy (2013); Zheng and Baranger (2013); Xu et al. (2013); Laakso and Pletyukhov (2014); Gonzalez-Ballestero et al. (2014); Shi et al. (2015); Longo et al. (2010); Sánchez-Burillo et al. (2014); Şükrü Ekin Kocabaş (2016).

The main result in this work is the structure of the $N$-photon $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 examples. The first one is a numerical study of scattering in the ultrastrong coupling limit Sánchez-Burillo et al. (2014, 2015), 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 Xu and Fan (2016) 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

II.1 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 $\hbar=1$, it reads

[TABLE]

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

[TABLE]

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_{\text{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

[TABLE]

have a finite support centered around $x_{\rm sc}=0$. The model (2) is not exactly solvable in general. For instance, if the scatterer is a two-level system, $H_{\text{sc}}\propto\sigma_{z}$ and $G=\sigma_{x}$ the model is the celebrated spin-boson model  Leggett et al. (1987), 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_{\text{min}},\omega_{\text{max}}]$ and typically a chiral medium $k\geq 0$, $\partial_{k}\omega_{k}\geq 0$. This is a rather standard simplification which does not affect the generality and applicability of our results. More generic dispersion relations and non-chiral 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.

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

[TABLE]

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 $e^{ik\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

[TABLE]

and Lorentzian envelopes

[TABLE]

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 $\sigma$ in momentum space implies a localization length $1/\sigma$ in position space. Note that our definition of the wave packet lacks factors such as $\sqrt{\omega_{k}}$ or $\omega_{k}^{1/2}$ 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 $\sigma\to 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 $\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.

II.3 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 $\ket{\Omega_{\mu}}$ are indistinguishable from the vacuum state $\ket{\rm{vac}}$ far away from the scatterer. Mathematically this occurs when

[TABLE]

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

Besides, the free particle states must satisfy the asymptotic condition Taylor (1972):

[TABLE]

with $U(t)$ the evolution operator of the full Hamiltonian (2) and $U^{0}(t)=e^{-iH_{0}t}$ the free-evolution operator.

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

[TABLE]

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

[TABLE]

Here, $U_{I}(t_{+},t_{-})=e^{iH_{0}t_{-}}\,e^{-iH(t_{+}-t_{-})}\,e^{-iH_{0}t_{+}}$ is the evolution operator in the interaction picture. Using again Eq. (9) leads to $|\Psi_{\rm in/out}\rangle=U_{0}^{\dagger}(t_{-/+})U(t_{-/+})|\Psi\rangle\equiv|\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:

[TABLE]

with $\psi^{\rm in}(t_{-})^{\dagger}=\psi_{\bar{k}\bar{x}}(t_{-})^{\dagger}$ and an analogous definition for $\psi^{\rm 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.

II.4 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

[TABLE]

provided that (i) for all $k$, $|g_{k}/\omega_{k}|<\infty$ and (ii) that the correlators $C_{kp}=\braket{\Omega_{0}}{a_{k}^{\dagger}a_{p}}{\Omega_{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 Longo et al. (2010); Sánchez-Burillo et al. (2014); Calajó et al. (2016); Shi et al. (2016), which in the literature are usually referred as ground states. Two paradigmatic examples of scatterer with multiple ground states are the three-level $\Lambda$ atom, with two electronic ground state, and a two-level system coupled to a cavity array in the ultrastrong coupling regime Sánchez-Burillo et al. (2014).

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, $\{\ket{\Omega_{\mu}}\}$, 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:

[TABLE]

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_{\rm sc}\propto\sigma_{z}$ and $G=\sigma_{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. Sánchez-Burillo et al. (2014, 2015); Shi et al. (2017) and Sect. IV.1].

III Causality and the $N$-photon scattering matrix

III.1 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(\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 App. B that

[TABLE]

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_{0}$. In the Heisenberg picture, the phases $ik(\bar{x}-\bar{y})-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}^{ikz}f(k)\,dk\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_{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 Lieb and Robinson (1972), and which were later generalized to finite-dimensional models, anharmonic oscillators, master equations, and spin-boson lattices Hastings and Koma (2006); Nachtergaele and Sims (2006); Nachtergaele et al. (2009); Poulin (2010); Barthel and Kliesch (2012); Jünemann et al. (2013). 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.

III.2 Causality and the scattering matrix

Causality imposes restrictions on the $S$-matrix  Weinberg (1995), 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

[TABLE]

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 $iT_{\textbf{p}\textbf{k}}=iC_{\textbf{p}\textbf{k}}\delta(E_{\textbf{p}}-E_{\textbf{k}})$ Weinberg (1996). 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

[TABLE]

with $S_{p_{n}k_{n}}\propto\delta(\omega_{p_{n}}-\omega_{k_{n}})$ the one-photon $S$-matrix. This is nothing but the cluster decomposition. Fourier transforming $S^{0}_{\bf pk}$, this structure also holds

[TABLE]

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

III.3 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,

[TABLE]

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,

[TABLE]

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

[TABLE]

with $A_{1,\nu\to\lambda}=\bra{\Omega_{\lambda}}\psi_{\bar{p}_{1}\bar{y}_{1}}^{\rm out}(t_{+})\psi_{\bar{k}_{1}\bar{x}_{1}}^{\rm in}(t_{-})^{\dagger}\ket{\Omega_{\nu}}$ and similarly for $A_{2,\lambda\to\mu}$. We can generalize this expression to $N$ photons, with initial average positions $\bar{x}_{1}>\bar{x}_{2}>\dots>\bar{x}_{N}$ and asymptotic ground states $\lambda_{0}:=\nu$ and $\lambda_{N}:=\mu$

[TABLE]

with

[TABLE]

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=\Pi_{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, $\partial_{k}\omega_{k}\geq 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 Fig. 1). Some algebra, fully described in Appendix D yields that $S^{0}$ has the following structure

[TABLE]

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_{\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 Sect. IV.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 $\Lambda$ atom by Xu and Fan Xu and Fan (2016).

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 one-photon 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

IV.1 Ultrastrong scattering

Let us consider a system described by the following Hamiltonian

[TABLE]

The scatterer is a two-level system described by the ladder operators $\sigma^{\pm}$ and the level splitting $\Delta$. The lattice tight-binding Hamiltonian, describes an array of identical cavities with frequency $\epsilon$, cavity-cavity coupling $J$, and bosonic modes $[a_{x},a_{y}^{\dagger}]=\delta_{xy}$.

The lattice model is diagonalized in momentum space, giving raise to a cosine-shaped dispersion relation, $\omega_{k}=\epsilon-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(\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 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 $|{\rm vac}\rangle$ with $\sigma^{-}|{\rm vac}\rangle=a_{x}|{\rm 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 (IV.1). 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 Niemczyk et al. (2010); Forn-Díaz et al. (2010, 2017); Yoshihara et al. (2017) 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}]\neq 0$, numerical simulations indicate that the total number of flying photons is asymptotically conserved throughout the simulation Sánchez-Burillo et al. (2014, 2015). Therefore, the second condition needed for proving our results is fulfilled (see Sect. II.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 Vidal (2003, 2004); García-Ripoll (2006); Verstraete et al. (2008), which has been recently adapted to the photonic world in Sánchez-Burillo et al. (2014, 2015, 2016a). Using this ansatz, we computed the nontrivial minimum-energy state Sánchez-Burillo et al. (2014), 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 $\ket{\Omega_{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. (68) and (69)) 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. (68), with $\ket{\Omega_{\nu}}=\ket{\Omega_{0}}$. Both wave packets will be Gaussians, Eq. (6), with mean momentum $\bar{k}$ and width $\sigma$. The numerical simulations show that the scattering is elastic for the chosen parameters ($\epsilon=1$, $J=1/\pi$, $\Delta=\epsilon=1$, and $g=0.3$) Sánchez-Burillo et al. (2014).

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)=\bra{\Omega_{0}}a_{p_{1}}a_{p_{2}}\ket{\Psi(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

[TABLE]

with $p_{1}$ and $p_{2}$ such that $\omega_{p_{1}}+\omega_{p_{2}}=2(\omega_{\bar{k}}\pm\sigma_{\omega})$ and $\omega_{p_{1}},\omega_{p_{2}}\not\in(\omega_{\bar{k}}-\sigma_{\omega},\omega_{\bar{k}}+\sigma_{\omega})$, being $\sigma_{\omega}$ 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 $\phi_{p_{1},p_{2}}(t_{+})\neq 0$ even though $|p_{1}|\neq\bar{k}$ or $|p_{2}|\neq\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 Shen and Fan (2007). 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}|\neq\bar{k}$ or $|p_{2}|\neq\bar{k}$. All this is a consequence of the cluster decomposition, see Eq. (21) and Th. 4 in App. C.

IV.2 Inelastic scattering and linear dispersion relation: the cluster decomposition revisited

We set $\omega_{k}=c|k|$ in $H_{0}$. The scatterer and interaction are described by

[TABLE]

where $\{\ket{\Omega_{\nu}}\}$ and $\{\ket{J}\}$ 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 $\ket{\Omega_{\nu}}\leftrightarrow\ket{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^{\prime}=1$, the scatterer is a $\Lambda$ 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

[TABLE]

with $k$ and $p$ the incident and outgoing momenta, respectively, and $|\Omega_{\nu}\rangle$ and $|\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^{0}$-matrix, Eq. (23) in momentum space is

[TABLE]

Here, $n^{\prime}\neq n$, e.g., $n^{\prime}=2$ if $n=1$. The computation is detailed in Appendix E. This structure has recently been found by Xu and Fan for a $\Lambda$ atom ($M=2$, $M^{\prime}=1$) within the RWA and Markovian approximations Xu and Fan (2016). 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.2. The mathematical origin of the structure can be traced back to its form in position space, Eq. (23). 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 App. E]. Besides, 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_{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 Weinberg (1996). Then, we assume that $(C_{p_{1}p_{2}k_{1}k_{2}})_{\mu\nu}$ has simple poles with imaginary parts $\{\gamma^{C}_{n}\}$. Similarly, we expect that divergences of $t_{\mu\nu}(k)$ come from simple poles with imaginary parts $\{\gamma^{t}_{n}\}$. As far as we know, this structure has been found for all $S$-matrices in waveguide QED Fan et al. (2010); Rephaeli et al. (2011); Xu and Fan (2016); Sánchez-Burillo et al. (2016b).

Let us write down the input state in momentum space

[TABLE]

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

[TABLE]

with the two-photon wave function $\phi_{\mu}^{\rm out}(p_{1},p_{2})$

[TABLE]

being $(\tilde{C}_{p_{1}p_{2}k_{n}})_{\mu\nu}=\int dk_{\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}\neq 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 $\gamma_{n}^{t}$ and $\gamma_{n}^{C}$, together with the exponentials $e^{ik_{n}l}$ and $e^{i(p_{1}+p_{2}+E_{\mu}-k_{n}-E_{\lambda})l}$, gives an exponentially decaying term, $e^{-|\gamma_{n}^{t}|l}$ or $e^{-|\gamma_{n}^{C}|l}$. We choose Lorentzian envelopes for the wave packets. They have a pole at $\bar{k}-i\sigma$ [see Eq. (7)]. In consequence, the wave packets will give a term proportional to $e^{-\sigma l}$. Lastly, the imaginary part of the pole of the first term vanishes, $\sim i0^{+}$, so it gives a nondecaying term, $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^{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 $\{\gamma_{n}^{t}\}$ and $\{\gamma_{n}^{C}\}$ by measuring the decay of $F$ with the distance.

V Final Comments

Our work represents a significant evolution over the field-theoretical methods Roy et al. (2017) 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 Xu and Fan (2016).

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 Sánchez-Burillo et al. (2014, 2015). As a second example, we have considered a non-dispersive media $\omega_{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 $\Lambda$ atom Xu and Fan (2016) 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.4], 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 photon-creation events (see e.g. Sánchez-Burillo et al. (2016a)). Finally, our program can be extended to treat other systems, deriving a cluster decomposition for the scattering of spin waves in quantum-magnetism models or for fermionic excitations in many-body systems.

Acknowledgements.

We acknowledge support by the Spanish Ministerio de Economía y Competitividad within projects MAT2014-53432-C5-1-R and FIS2015-70856-P, the Gobierno de Aragón (FENOL group), and CAM Research Network QUITEMAD+ S2013/ICE-2801.

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 $\ket{\Omega_{0}}$,

[TABLE]

Proof.

Let us assume that $\ket{\Omega_{0}}$ is the minimum-energy state of $H$ as given by Eq. (2), and thus $(H-E_{0})\ket{\Omega_{0}}=0$. The energy of the unnormalized state

[TABLE]

created by any operator $O$ must be larger or equal to that of the ground state, $\braket{\chi}{(H-E_{0})}{\chi}\geq 0$. Using (34)

[TABLE]

we conclude with the useful relation

[TABLE]

Let us take $O=a_{k}$. The previous statement leads to

[TABLE]

or equivalently

[TABLE]

Using Cauchy-Schwatz, this translates into the upper bound

[TABLE]

Once the diagonal elements of the correlation matrix are bounded the nondiagonal can also be bounded. The correlation matrix is positive $C\geq 0$ with $C_{kp}=\braket{\Omega_{0}}{a_{k}^{\dagger}a_{p}}{\Omega_{0}}$. A property of positive matrices is Horn and Johnson (1985)

[TABLE]

which implies (33). ∎

With this lemma at hand we state:

Theorem 1**.**

Let us define $\psi_{\bar{k}\bar{x}x}^{\dagger}$ as the operator (5) removing the time-dependent part, where $\phi_{\bar{k}}(k)$ is infinitely differentiable with a finite support $K$ centered around $\bar{k}$. Then, the expected value of $\psi_{\bar{k}\bar{x}}^{\dagger}\psi_{\bar{k}\bar{x}}$ in the minimum-energy state fulfills

[TABLE]

where we choose $x_{\text{sc}}=0$. Moreover, if we can assume that $\braket{a_{k}^{\dagger}a_{p}}$ is an $n$-times differentiable function of $k$ and $p$, the bound will be improved

[TABLE]

Proof.

Let us compute the expectation value of the number operator for a wave packet $N:=\braket{\Omega_{0}}{\psi_{\bar{k}\bar{x}}^{\dagger}\psi_{\bar{k}\bar{x}}}{\Omega_{0}}$,

[TABLE]

We can rewrite $N$ as the Fourier transform of another function $N=\int e^{iu\bar{x}}F(u)\,du$, where

[TABLE]

We are now going to assume that $\phi_{\bar{k}}(k)$ is a test function with compact support $K$ of size $|K|$ centered around $\bar{k}$, and infinitely differentiable. We will also assume that within its support $|g_{k}/\omega_{k}|^{2}\langle GG^{\dagger}\rangle\leq C_{\phi}$ for some constant $C_{\phi}$. Then we can bound

[TABLE]

Assuming that $\langle\Omega_{0}|a_{k}^{\dagger}a_{p}|\Omega_{0}\rangle$ is $n$-times differentiable and using the Riemann-Lebesgue theorem, we have then that

[TABLE]

at long distances. ∎

Appendix B Approximate causality

B.1 Free-field causality

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)$, 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 $\omega_{k}$ have an upper bounded group velocity $v_{k}=\partial_{k}\omega_{k}$:

[TABLE]

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

[TABLE]

Proof.

Solving the equation $f^{\prime}(k)=x-\partial_{k}\omega_{k}t=0$ leads to the condition $\frac{x}{t}=v_{k}$ or $|x/t|=|v_{k}|\leq c$. Then, provided $f^{\prime}(k)=0$, it follows $|x|\leq c|t|\Rightarrow d_{c}(x,t)\leq 0$, which shows (48). ∎

Lemma 3**.**

Assume that $\omega_{k}$ is $n$-times differentiable and that every derivative $|\omega_{k}^{(r\leq n)}|$ is upper bounded by an $m$-th order polynomial in $|k|$. Then the following integral bound applies

[TABLE]

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

Proof.

Result 5.1 from Ref. Olver (2014) states that the integral $I(x)=\int_{a}^{b}e^{ikx}q(k)\,dk$ may be integrated by parts $n$ times, obtaining

[TABLE]

where the error term satisfies

[TABLE]

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)=e^{-\frac{1}{\sigma^{2}}{(k-k_{0})^{2}}-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 Olver (2014). Since $x^{-s}q^{(s)}(a)\to 0$ when $a\to\pm\infty,\,\forall x,$ we obtain

[TABLE]

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 dk\,\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}|\leq 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$:

[TABLE]

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$,

[TABLE]

Proof.

Let us assume that the model evolves freely according to the free Hamiltonian $H_{0}=\int dk\;\omega_{k}a_{k}^{\dagger}a_{k}$. In this case, our wave packet operators have the simple form

[TABLE]

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

[TABLE]

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

[TABLE]

The exponent $\tilde{\omega}_{k}=\omega_{k}-ck$ is $n$-times differentiable and is upper bounded in modulus by a polynomial of degree $m\geq 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 Eq. (48), $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

[TABLE]

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

B.2 Full model causality

Causal relation (54) 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 $\omega_{k}$ and $g_{k}$, with degrees $n\geq 2$. Then, all wave packets outside the light cone of the scatterer evolve approximately with the free Hamiltonian, $H_{0}$. More precisely, if $(\bar{x},t_{1})$ and $(\bar{x},t_{0})$ are two points outside the light cone

[TABLE]

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

[TABLE]

is the free-evolution operator for the photons at time $t_{0}$.

Proof.

We start by building the Heisenberg equations for the operators

[TABLE]

Making the change of variables $a_{k}(t)=e^{-i\omega_{k}t}b_{k}(t)$, we have

[TABLE]

so that the wave packet operators evolved from some initial time $t_{s}$ are

[TABLE]

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 $\|{G}\|=1$, with $||\cdot||$ the Hilbert-Schmidt norm, and $|t_{1}|>|t_{0}|$. We have to choose the integration limits $t$ and $t_{s}$ so that $\rm{sign}(\tau^{\prime})=\rm{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 ($\tau^{\prime}<0$). This means we can introduce $\tau^{\prime\prime}=\rm{sign}(x)\tau^{\prime}\geq 0$ and bound

[TABLE]

Here we have taken into account that $d_{c}(|\bar{x}|,\tau^{\prime\prime})\geq d_{c}(|\bar{x}|,|t_{1}-t_{0}|)>0$ in the domain of integration. We can now use the fact that $d_{c}(|\bar{x}|,|t_{1}-t_{0}|)\geq d_{c}(|\bar{x}|,|t_{1}|)\geq\min\{d_{c}(\bar{x},t_{1}),d_{c}(\bar{x},t_{0})\}$, obtaining the expression in the theorem.

∎

B.3 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.3 and Eq. (9)]. However, using Theorems 3 and 1 we have that, given a state $|\Psi\rangle\equiv\psi_{\bar{k},\bar{x}}(t_{0})^{\dagger}|\Omega_{\nu}\rangle$, then

[TABLE]

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 II.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.

Appendix C Scattering amplitude decomposition

Theorem 4**.**

Let us suppose the input state is

[TABLE]

with $|\bar{x}_{n}-\bar{x}_{m}|\to\infty$ $\forall n\neq m$. Thus, the scattering amplitude of going to

[TABLE]

with $|\bar{y}_{n}-\bar{y}_{m}|\to\infty$ $\forall n\neq m$, is reduced to a product of single-photon events:

[TABLE]

being $\lambda_{0}=\mu$ and $\lambda_{N}=\nu$, with the wave packet operators given in the Heisenberg picture for $t=t_{\pm}\to\pm\infty$.

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,

[TABLE]

In the second equality we have dropped an irrelevant global phase. Here, $\psi_{\rm in}^{\dagger}$ and $\psi_{\rm out}^{\dagger}$ 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. (68) and (69) the amplitude is given by

[TABLE]

As $|\bar{x}_{1}-\bar{x}_{2}|$ can be arbitrarily large, we can always choose a time $t_{1}$ such that $\psi_{\bar{p}_{1}\bar{y}_{1}}^{\text{out}}(t)^{\dagger}\ket{\Omega_{\mu}}$ is well separated from the scatterer for $t>t_{1}$, so $\psi_{\bar{p}_{1}\bar{y}_{1}}^{\text{out}}(t)\cong U_{0}(t,t_{1})^{\dagger}\psi^{\text{out}}_{\bar{p}_{1}\bar{y}_{1}}(t_{1})U_{0}(t,t_{1})$. Besides, $t_{1}$ is such that the second wave packet is still far away from the scatterer. Therefore $\psi_{\bar{k}_{2}\bar{x}_{2}}^{\text{in}}(t^{\prime})\cong U_{0}(t^{\prime},t_{1})^{\dagger}\psi_{\bar{k}_{2}\bar{x}_{2}}^{\text{in}}(t)U_{0}(t^{\prime},t_{1})$, for $t^{\prime}<t_{1}$. Using Theorem 2, $[\psi_{\bar{p}_{1}\bar{y}_{1}}^{\text{out}}(t_{+}),\psi_{\bar{k}_{2}\bar{y}_{2}}^{\text{in}}(t_{-})^{\dagger}]\to 0$ and Eq. (72), the amplitude equals to

[TABLE]

Finally, we insert the identity between the operators $\psi_{\bar{k}_{2}\bar{x}_{2}}^{\rm in}(t_{-})^{\dagger}$ and $\psi_{\bar{p}_{1}\bar{y}_{1}}^{\rm out}(t_{+})$. Assuming there is not particle creation and just the ground states $\{\ket{\Omega_{\lambda}}\}_{\lambda=0}^{M-1}$ will contribute to the identity, $\sum_{\lambda=0}^{M-1}\ket{\Omega_{\lambda}}\bra{\Omega_{\lambda}}$, and we arrive to (70).

This comes because $\psi_{\bar{k}_{2}\bar{x}_{2}}^{\rm in}(t_{-})^{\dagger}$ and $\psi_{\bar{p}_{1}\bar{y}_{1}}^{\rm out}(t_{+})$ asymptotically commute but not $\psi_{\bar{k}_{1}\bar{x}_{1}}^{\rm in}(t_{-})^{\dagger}$ and $\psi_{\bar{p}_{2}\bar{y}_{2}}^{\rm out}(t_{+})$. This is a clear signature of causality, saying which one is arriving first. Lastly, notice that if the ground state is unique, $\ket{\Omega_{\lambda_{n}}}=\ket{\Omega_{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. (70). We do it in the two-photon subspace.

Before, we need the one-photon amplitude as an intermediate result.

D.1 One photon

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

[TABLE]

with the creation operator $\psi_{\bar{k}_{1},\bar{x}_{1}}^{\text{in}\;\dagger}$ given by Eq. (5), removing the time dependence. For simplicity, we absorb the factor $e^{ik\bar{x}_{1}}$ into the wave packet: $\phi_{\bar{k}_{1},\bar{x}_{1}}(k)=e^{ik\bar{x}_{1}}\phi_{\bar{k}_{1}}(k)$. In position space, the output state will read

[TABLE]

Defining

[TABLE]

and

[TABLE]

being $\ket{y;\Omega_{\mu}}=a_{y}^{\dagger}\ket{\Omega_{\mu}}$ the state with a photon at $y$ and the scatterer in the ground state $\ket{\Omega_{\mu}}$, the output state (75) can be rewritten as

[TABLE]

The probability amplitude will read

[TABLE]

If the wave packets are monochromatic with momenta $k_{1}$ and $p_{1}$, respectively, this amplitude is

[TABLE]

D.2 Two photons

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

[TABLE]

By definition, the output state is

[TABLE]

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

[TABLE]

with

[TABLE]

being $|x_{1}x_{2};\Omega_{\nu}\rangle=a_{x_{1}}^{\dagger}a_{x_{2}}^{\dagger}\ket{\Omega_{\mu}}$ the symmetrized state with two photons at $x_{1}$ at $x_{2}$ and the scatterer at $\ket{\Omega_{\nu}}$.

Introducing (84) in (83) and considering (81) and (23) we get

[TABLE]

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

[TABLE]

Using Eq. (76)

[TABLE]

Following the sketch drawn in Fig. 1, if $x_{m}<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}<y_{n^{\prime}}$. Therefore, choosing $i=1$ and $j=2$, the integral $C$ reads

[TABLE]

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

[TABLE]

Finally, the probability amplitude of going to the output state $\psi_{\bar{p}_{1},\bar{y}_{1}}^{\text{out}\;\dagger}\psi_{\bar{p}_{2},\bar{y}_{2}}^{\text{out}\;\dagger}|\Omega_{\mu}\rangle$ will be the overlap between this state and (89). Using (D.1)

[TABLE]

as expected. In the calculations, we have set $\langle\Omega_{\mu}|\psi_{\bar{p}_{i},\bar{y}_{i}}^{\text{out}}\;S\;\psi_{\bar{k}_{j}\bar{x}_{j}}^{\text{in}\;\dagger}|\Omega_{\nu}\rangle=0$ for $i\neq 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,\nu\to\lambda}A_{1,\lambda\to\mu}$ 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_{p_{1}p_{2}k_{1}k_{2}}^{0})_{\mu\nu}$ as the Fourier transform of $(S_{y_{1}y_{2}x_{1}x_{2}}^{0})_{\mu\nu}$

[TABLE]

Due to the form of $(S_{y_{1}y_{2}x_{1}x_{2}}^{0})_{\mu\nu}$, (23), we have to compute integrals as

[TABLE]

Notice that $(S_{yx})_{\mu\nu}$ is the Fourier transform of $(S_{pk})_{\mu\nu}$, Eq. (28). Therefore,

[TABLE]

Considering this in (91), we get

[TABLE]

with $n^{\prime}\neq n$ and $m^{\prime}\neq m$. The Fourier transform of the step function is

[TABLE]

Therefore, integrating Eq. (E) first in $y_{1}$ and later in $y_{2}$, we get

[TABLE]

which is the expression given in the main text, Eq. (29). This result has been recently reported for a $\Lambda$ atom by Xu and Fan in Xu and Fan (2016). 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

[TABLE]

with $\mathcal{P}$ the principal value. Applying this identity to Eq. (29) we get,

[TABLE]

Now, we sum over $n$ and $m$

[TABLE]

Applying the constraint imposed by the global Dirac delta to $p_{2}$ to the second row, it is straightforward to see that they cancel each other, arriving to

[TABLE]

which is the usual expression in translational invariant (momentum conserving) QFT for the cluster decomposition, which also holds in waveguide QED if the ground state is unique.

Appendix F Fluorescence decay

In this appendix, we calculate how the correlations and thus the fluorescence decay as the distance $l$ between the packets grows (See Figs. 1 and 3).

The input state (68) in momentum space is given by,

[TABLE]

with

[TABLE]

In these expressions, the wave packets $\phi_{\bar{k}_{n}}(k)$ are Lorentzian functions [see Eq. (7)]. The out state is computed by means of Eq. (10)

[TABLE]

With $\mathbb{I}$ the identity operator in the two-photon sector: $\mathbb{I}=1/2\int dp_{1}dp_{2}\sum_{\mu}a_{p_{1}}^{\dagger}a_{p_{2}}^{\dagger}\ket{\Omega_{\mu}}\bra{\Omega_{\mu}}a_{p_{1}}a_{p_{2}}$. The scattering matrix $S$ in momentum space is $(S_{p_{1}p_{2}k_{1}k_{2}})_{\mu\nu}=(S^{0}_{p_{1}p_{2}k_{1}k_{2}})_{\mu\nu}+i(T_{p_{1}p_{2}k_{1}k_{2}})_{\mu\nu}$, with $(S^{0}_{p_{1}p_{2}k_{1}k_{2}})_{\mu\nu}$ given by Eq. (29) and $(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})$ yielding

[TABLE]

with

[TABLE]

Which is nothing but Eq. (32) that we have rewritten here for the discussion. As said in Sect. IV.2, we assume that $t_{\mu\nu}(k)$ and $(C_{p_{1}p_{2}k_{n}k_{\bar{n}}})_{\mu\nu}$ have simple poles with imaginary parts $\{\gamma_{n}^{t}\}$ and $\{\gamma_{n}^{C}\}$ 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(p_{1}+p_{2}+E_{\mu}-k_{n}-E_{\lambda})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^{ik_{n}l}$ 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^{-|\gamma_{n}^{t}|l}$ and $e^{-|\gamma_{n}^{C}|l}$, respectively.

Now we consider the contribution to the integral of the wave packets, $\phi_{\bar{k}_{n}}(k)$. We choose Lorentzian functions, with a simple pole at $k=\bar{k}_{n}-i\sigma$ (see Eq. (7)). In consequence, we have a term proportional to $e^{-\sigma 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 single-photon energy conservation giving single-photon amplitudes, $\sum_{\lambda}A_{1,\nu\to\lambda}A_{2,\lambda\to\mu}$.

Finally, let us mention that we do not need to impose that that $t_{\mu\nu}(k)$ and $(C_{p_{1}p_{2}k_{n}k_{\bar{n}}})_{\mu\nu}$ have simple poles. Higher order poles, by virtue of the Cauchy Integral formula for the derivatives, also would yield exponential decay.