Bound states and entanglement generation in waveguide quantum electrodynamics

We investigate the behavior of two quantum emitters (two-level atoms) embedded in a linear waveguide, in a quasi-one-dimensional configuration. Since the atoms can emit, absorb and reflect radiation, the pair can spontaneously relax towards an entangled bound state, under conditions in which a single atom would instead decay. We analyze the properties of these bound states, which occur for resonant values of the interatomic distance, and discuss their relevance with respect to entanglement generation. The stability of such states close to the resonance is studied, as well as the properties of non resonant bound states, whose energy is below the threshold for photon propagation.


I. INTRODUCTION
An excited atom in free space unavoidably decays towards its ground state through spontaneous emission. Boundary conditions and artificial dimensional reduction can drastically modify the picture, providing situations in which the decay is enhanced, inhibited or even completely hindered [1][2][3][4][5][6][7][8][9]. While confinement in optical cavities has long been a common way to study the effects of geometry [10][11][12], one-dimensional systems have recently emerged as another promising stage for the observation of interesting quantum electrodynamics (QED) phenomena. Nowadays a variety of quantum emitters (atoms for brevity) can be coupled to quasi one-dimensional fields such as waveguides, optical fibers and microwave transmission lines [13][14][15][16][17][18][19]. Alternatively, the effective reduction to one dimension can be obtained by tightly focusing photons [20][21][22]. These impressive experimental advances have opened the way to unexplored nonperturbative regimes of QED, and have motivated work on the interaction between atoms and waveguides in different geometries [23][24][25][26][27][28][29][30].
In this context, an interesting problem is the study of atoms in semi-infinite linear waveguides, where one end of the guide behaves as a perfect mirror [31][32][33]. For selected values of the atom-mirror distance a nontrivial bound state exists, in which the probability of atomic excitation is finite, even when photons emitted through spontaneous decay can propagate in the guide [34,35]. The optical path between the atom and the mirror is crucial for the existence of this kind of resonance. It is worth noting that even a single atom exhibits a mirror-like behavior in one dimension [21,22,[36][37][38][39]. One may thus consider the interaction of two atoms, mediated by the exchange of photons propagating in one dimension, and exploit the dual behavior of each atom as both an emitter and a mirror. Such interaction can give rise to stable configurations in which the atoms display significant en-tanglement, while the field is confined between the atoms and does not propagate [39][40][41]. Besides the fundamental interest of few-body QED in quasi 1D geometries, where non-Markovian effects easily come into play [42], such a system is thus interesting from the point of view of generating entanglement, an important resource in Quantum Information, by relaxation. Indeed, if a bound state exists in which the two atoms are entangled, an initially factorized atomic state can spontaneously relax towards a state with finite entanglement. Relaxation occurs after an initial transient in which photon exchange builds up quantum correlations. Differently from other methods of entanglement generation in waveguide-QED [43], this process would not require a continuous pumping of energy into the system, and would ideally provide a constant entanglement in time after the initial transient.
In this paper we show how the properties of bound and quasi-bound entangled states in waveguide-QED can be studied in great depth and generality by exploiting the resolvent formalism [44,46]. Studying the resolvent, one notices the presence of a number of poles in the so-called complex-energy plane. Each pole can be associated with a (generally unstable) state, and the imaginary part of a pole is proportional to the inverse lifetime of the state. This allows us to immediately identify a favorable situation for entanglement generation by relaxation: we need one of these poles to be a long-lived entangled state (i.e., the pole must have a negligible imaginary part), while the remaining poles must be fast-decaying states. Under such conditions, a separable atomic state would quickly relax onto an entangled metastable state. While the metastable state will eventually decay due to losses and imperfections, our analysis allows to clearly identify the relevant timescales of the problem. Thus, we can give a clear indication of what degree of losses and imperfections a given system is able to tolerate while still allowing the generation of long-lived entanglement. Importantly, our formalism automatically takes into account a number arXiv:1604.01300v2 [quant-ph] 21 Oct 2016 Figure 1. Two two-level atoms are placed at relative distance d in a one-dimensional waveguide, with propagation direction along the x axis. Both atoms possess the same internal structure (for brevity we only sketch the levels of emitter A) and interact through the mediation of waveguide photons. The waveguide is characterized by its one-dimensional photon dispersion relation ω(k), with k being the photon momentum. We first focus on the TE1,0 mode of an infinite waveguide of rectangular cross section, and then generalize to a wide class of one-dimensional dispersion relations in Sec. V.
of physical effect that are often neglected: these include non-markovian effects, time-delay (due to the finite propagation speed of photons), threshold effects (due to the presence of either high-or low-frequency cutoffs in the dispersion relation). Even though, for definiteness, we focus for the most part on the dispersion relation typical of rectangular waveguides, we will also outline how a wide class of physically relevant dispersion relations can be tackled within the same framework.
Our paper is organized as follows. In section II we introduce the Hamiltonian of our model, and illustrate how a suitable choice of the inter-atomic distance gives rise to entangled bound states above the frequency threshold for photon propagation. Section III is devoted to the study of poles in the complex-energy plane, which allows us to extract crucial information relevant to the entanglementby-relaxation protocol. In section IV we extend our analysis to off-resonant bound states, whose energy is below the low-frequency cutoff of the waveguide. We outline in section V how our study can be generalized in a straightforward manner to any dispersion relation that satisfies appropriate conditions. Finally, we draw our conclusions in section VI.

II. THE MODEL
We describe the dynamics of two two-level atoms A and B, situated in an infinite waveguide of rectangular cross section, with sides L y < L z , see Fig. 1. When longitudinal propagation occurs with long wavelength compared to the transverse size, interaction between atoms and field can be reduced to a coupling with the lowestcutoff-energy TE 1,0 mode, in which the electric field vibrates along the z direction and has a sine modulation in the y direction [48]. In this situation, the electromagnetic field is effectively scalar and massive. The inter-acting atoms and photons are described, in dipolar and rotating wave approximations, by the Hamiltonian where ω 0 is the bare energy separation between the atomic ground |g and first-excited states |e , λ is the coupling constant (see Appendix B), d is the A-B distance, ω(k) is the photon dispersion relation, and b(k) (b † (k)) is the annihilation (creation) field operator, satisfying the canonical commutation relation [b(k), b † (k )] = δ(k − k ). Henceforth, we will focus on the dispersion ω(k) = √ k 2 + M 2 of the TE 1,0 mode in the waveguide, characterized by a mass M ∝ L −1 y . However, as discussed in section V, our approach is applicable to a wide class of one-dimensional dispersion relations. The effective mass M provides a natural cutoff to the coupling. The Hamiltonian (1) commutes with the excitation number where N at = |e A e A | + |e B e B | is the atomic excitation number. The N = 0 sector is 1-dimensional and is spanned by the bare ground state |g A , g B ; vac . We shall focus instead on the dynamics in the N = 1 sector, where the states read (3) where |ϕ := dk ϕ(k)b † (k)|vac is a one-photon state, and |c A | 2 + |c B | 2 + dk|ϕ(k)| 2 = 1.
In the small-coupling regime, an isolated excited atom with ω 0 M would decay to the ground state. We shall demonstrate that, when two atoms are considered, a resonance effect emerges, yielding a bound state. Using the expansion (3) the eigenvalue equation, H|ψ = E|ψ , reads The field amplitude ϕ(k) has two simple poles at k = ±k = ± √ E 2 − M 2 . Thus, when E > M , the integrals in (4)- (5) are finite only if c A +c B e ±ikd = 0, yieldingkd = nπ for positive integers n. This implies that a bound state can exist only for discrete values of the interatomic distance d. Moreover, in the first component of such an eigenstate (3), the atoms are in a maximally entangled (singlet or triplet) state, namely c A = (−1) n+1 c B . To determine the distances at which the bound state exists, let us first compute the energy eigenvalue, which after the resonance condition is the solution of and if the wavenumberk is real The properties of states with E < M , to which an imaginary wavenumber can be associated, will be discussed in Section IV.
To complete the characterization of the bound state, we shall analyze the atomic populations and the field energy density. The former can be immediately obtained using the normalization condition on the states (3) as Where we use the shorthands c B | 2 as the probability associated to the N at = 1 sector, one gets Notice that, despite being apparently of order λ 2 , the correction to unity is given by the ratio between powers of two small quantities, namely the effective coupling constant λ/M , and the wavenumbers ratiok/M . The resulting number can be of order one, even at small coupling constants. Observe that the probability vanishes likek 3 /n at very smallk: this behavior is physically motivated by the fact that, as the energy approaches the cutoff, the distance between the atoms must increase to infinity in all bound states. Let us finally analyze the energy density of the electromagnetic fields. Neglecting the exponentially suppressed contribution of the square-root cuts, the energy density turns out to be related to the Fourier transform of the photon amplitude, as , and E n (x) 0 outside. Thus, the field is confined between the two atoms, and modulated with periodicity π/k, with nodes at the positions of the atoms which act as mirrors. This explains the occurrence of such bound states for discrete values (8) of the interatomic distance.
Moreover, the structure of the bound state is where s = (−1) n+1 and |Ψ ± = (|e A , g B ± |g A , e B )/ √ 2 are (maximally entangled) Bell states. This is a key feature which enables entanglement generation by atomphoton interaction. Indeed, suppose that d = d n : a factorized initial state, say |ψ(0) = |e A , g B ⊗ |vac , can be expanded into a "stable" and a "decaying" part as with ψ ⊥ n |ψ n = 0. After a transient of the order of |ψ ⊥ n 's lifetime (see discussion in the following), the atomic density matrix ρ at (t) := Tr field |ψ(t) ψ(t)| approaches (15) in which the atoms have a finite probability, determined by (10), to be maximally entangled. In Figure 2 we display the atomic entanglement in the asymptotic state, as measured by the concurrence [45]. However, one could also measure the photon state and obtain, with a finite probability, a maximally entangled atomic state. The strategy is therefore the following: one prepares a factorized state, and measures whether a photon is emitted. If (after a few lifetimes) no photon has been observed, the atomic state is projected over the maximally entangled Bell state |Ψ s . This can be achieved with higher probabilities for larger values of ω 0 . In realistic scenarios this simplified picture is challenged by the presence of losses, such that it is no longer possible to prepare an exact Bell state. Nevertheless, if losses occur on sufficiently long timescales (as compared to the decay rate of the fast pole -see section III below), and provided the detector efficiency is high enough, it remains possible to achieve high fidelity with a Bell state.

III. TIME EVOLUTION AND BOUND STATE STABILITY
Let us now study the general evolution of an initial state in the atomic sector N at = 1. We will use the resolvent formalism [44,46] to illustrate that the system  relaxes towards the bound state, and to quantify the robustness of the bound state against small variations in the model parameters (such as the A-B distance). We remark that the usefulness of the resolvent formalism goes beyond the analysis of stable states, in that it provides crucial information on the relevant timescales of the problem. Indeed, the entanglement-by-relaxation protocol described in the previous section relies on the fast decay of the unstable Bell state. The analysis of the resolvent enables to determine the lifetime of this unstable state, which must be much shorter than the typical timescales of waveguide or atomic losses, as well as the inevitably finite lifetime of the bound state (due for example to imperfect control of the A-B distance). Whenever these conditions are met, the effectiveness of the protocol is guaranteed and a long-lived entangled state may be prepared by relaxation.
The resolvent G(z) = (z − H) −1 , with z the complex energy variable, has singularities only on the real axis (on the first Riemann sheet) and the study of additional singularities (on the other Riemann sheets) yields crucial information about the dynamical stability of the system: in particular, a pole with a non-vanishing imaginary component signals a decay process. The resolvent approach yields results that are consistent with those obtained from the analysis of the Laplace transform of the time evolution [41].
For λ = 0, the free resolvent G 0 (z) = (z − H 0 ) −1 has a pole on the real axis, at z = ω 0 , corresponding to the excited states of atoms A or B. When interaction is turned on, this singularity splits into two simple poles, which generally migrate into the second Riemann sheet. We shall see from a non-perturbative analysis that, under resonance conditions, one of the poles falls on the real axis (and is therefore very long-lived), while the other one has a very short lifetime. Let G(z) and G 0 (z) be the restrictions to the N at = 1 sector of the interacting and free resolvent, respectively. In the basis {|e A , g B , |g A , e B } one gets G 0 (z) = 1 z − ω 0 1 0 0 1 (16) and is called self energy. The resolvent G(z) is analytic in the whole complex energy plane, except at points on the real axis that belong to the spectrum of the Hamiltonian H. In particular it exhibits simple poles at the eigenvalues of H and cuts along its continuous spectrum [44,46]. However, it can happen that some complex poles show up on the second Riemann sheet, through the analytic continuation of the resolvent G II (z) from the upper half-plane to the lower half-plane under the cut [46,47]. These poles physically correspond to unstable states with energy and decay rates given by their real and imaginary part, respectively.
The particular form of the interaction Hamiltonian V in (1) enables one to exactly evaluate the self energy: . (20) Due to the bare energy degeneracy and the symmetric structure of the self energy, the propagator can be diagonalized as where with spectral densities The self-energy functions Σ ± (z) are analytic in the cut complex energy plane C \ [M, +∞) and have a purely imaginary discontinuity across the cut proportional to the spectral density: During the continuation process into the second Riemann sheet through the cut, the self energy (22) will thus get an additional term Note that the new term has in general a nonvanishing imaginary part and is the analytical continuation of the discontinuity of the self-energy function across the cut. Now, a pole of G(z) on the second sheet must satisfy the equation for s = +1 or s = −1, where Σ II s (z) is the branch (25) of the self energy in the second sheet. By plugging (25) and (23) into (27) we get (28) It is evident from (28) how the energetic degeneracy at λ = 0 is lifted by interactions. Notice the presence of an imaginary component, detecting decay.
The last ingredient we need in order to get a closed expression for the complex energy poles is the evaluation of Σ s (z) in (28). Thus, let us rewrite (22) as an integral over k: The integrand function can be analytically continued to the complex k plane using the principal determination of the square root, which has nonnegative real part for all values of its argument, and is characterized by a branch cut for k 2 + M 2 < 0 , that is Two first-order poles, symmetric with respect to the origin of the k plane, are also present whenever Re(z) > 0: By deforming the integration contours as in Figure 3 and applying Jordan's theorem, Σ s is split in two terms coming from the upper branch cut and from one of the two poles (see Figure 3). Specifically, when Im(z) > 0, the pole k 0 (z) lies in the upper half plane, and the integral involves the residue Instead, when Im(z) < 0, the deformed contour in the upper plane encircles −k 0 (z), where the residue yields Finally, the integrals along the cut read where the contribution from e −χd , that is not amenable to an explicit closed form in terms of simple functions, is nevertheless suppressed like e −M d and can be neglected for large values of M d.
We are now able to recognize the real resonant poles discussed in the first part of the paper as special solutions of Eq. (28). Indeed, assuming that the complex energy pole (26) is far from the branching point z = M and that its imaginary part is almost vanishing, one can decouple the real and imaginary parts of (28) and obtain from Eqs. (32)-(35) where k p = E 2 p − M 2 . Hence, we find that the poles in the second Riemann sheet have a cyclic behavior with respect to d. This result is in agreement with the one obtained in Ref. [40] in the Markovian approximation. In particular, when d = d n , as defined in Eq. (8), the real part of the pole equations is solved by E p = k2 + M 2 . In this case, one of the poles corresponds to the entangled bound state, and has vanishing imaginary part, while the other signals an unstable state with associated decay rate Even if, strictly speaking, bound states only occur for discrete values of d, it can be readily checked that while the energy shift is linear, the decay rate of the stable pole is quadratic for d → d n implying that the state |ψ n remains very long-lived close to resonance. Eq. (39) quantifies the robustness of the bound states against variations of the parameter d. Note how Eqs. (38) and (39) provide essential information on the feasibility and effectiveness of the entanglement generation protocol: firstly, it is necessary that the condition γ (s) p γ (u) p is satisfied, which is equivalent to the conditionk 2 (d−d n ) 2 1. Secondly, γ (u) p must be much larger than any decay rate associated with loss processes (e.g. waveguide losses). Even though approximate analytical expressions such as Eqs. (38) and (39) are extremely valuable, we emphasize that our methodology is capable of capturing the exact behaviour of the poles against variations in the model parameters. To illustrate this, in Fig. 4 we show the trajectories of the poles (26) in the complex energy plane, obtained by fixing M and d and varying the bare excitation energy ω 0 . On the one hand, we are thus able to assess quantitatively the robustness of bound states against variations in ω 0 . On the other hand, Fig. 4 demonstrates how our methodology allows one to interpolate seamlessly between perturbative and non-perturbative regimes.

IV. OFF-RESONANT BOUND STATES
Let us briefly discuss the behavior of bound states with E < M . In this case, the atoms are not expected to decay. Nonetheless, they interact by coupling to the evanescent modes of the waveguide. Scrutiny of Eqs. shows that there are bound states for all d, whose energy satisfies where we have again neglected the O(e −M d ) contributions from branch-cut integration in the complex k plane. If the coupling is small and the excitation energy ω 0 is far from the threshold M for photon emission, the above equations reduce to an effective Hamiltonian eigenvalue equation in the N at = 1 sector. The eigenvalues read with c A = c B for the plus sign (ground state) and c A = −c B for the minus sign. These bound states are not associated to any resonance. It is also possible to check that the electromagnetic energy density falls like exp(− √ M 2 − E 2 |x|) away from the atoms. Since the eigenstates of the effective Hamiltonian are Bell states, the evolution of an initially factorized state is characterized by oscillations between two orthogonal maximally entangled states with period 2π/β(ω). Compared to entanglement by relaxation, this mechanism yields unit concurrence [41]. On the other hand, the process can be very slow, since the energy splitting is exponentially suppressed with the interatomic distance, and requires the fine tuning of an optimal time to stop the interactions, which is not required in the spontaneous entanglement process described in Section II.
The case E → M is more interesting, since the physics becomes nonperturbative. Equations (4)-(5) admit a singlet and a triplet solution. For the singlet case, the bound state with E = M is obtained at a finite excitation energy: in which, due to a cancellation of divergences, the correction to the bare energy is still perturbative in λ 2 . This singlet solution approximates the dark eigenstate |Ψ − ⊗ |vac occurring at d = 0. The triplet, instead, survives as a real eigenstate even for ω 0 ≥ M . However, since c A = c B implies that the integrals over the field become divergent in this limit, the population in the N at = 1 sector is suppressed to fulfill normalization, and the contribution of this pole to the expansion (14) can be safely neglected.

V. EXTENSION TO GENERIC DISPERSION RELATIONS
While we have worked out in detail the case of a rectangular waveguide, we emphasize that our methods can be applied to a generic dispersion relation ω(k). We start by noticing that Eqs. (4)-(6) lead in full generality to the implicit condition which must be satisfied by the bound state energy E. For the existence of a resonant (i.e. above threshold) bound state, it is evident that also the conditionkd = nπ, n ∈ N must hold, wherek > 0 satisfies E = ω(±k). If this were not the case, the right hand side of Eq. (45) would diverge. We assume that suchk exists and is unique. This is the case, for example, when ω(k) is an increasing lower-bounded function of |k|. Moreover, the possibility of non-resonant eigenstates below threshold follows as in the case of a rectangular waveguide. Moving on to the complex energy plane, the analysis of poles proceeds along the same lines as in section III, albeit the existence of compact analytical expressions will rely on the specific functional form of ω(k). The pole contribution to the self energies can be generalized by replacing the denominators in Eqs. (33)- (34) with ω(k 0 )ω (k 0 ), and the square root in the exponentials with k 0 . In the perturbative regime, this change does not affect formally the ratio of the decay rates of stable and unstable poles close to a resonance, namely [see Eqs. (38) where d n = nπ/k. We can see that quantitative differences between models arise in the inversion of the dispersion relation as a function of the energy. The quantity in Eq. (46) gives a clear indication of the potential of a given model to generate entanglement by relaxation. While losses would inevitably degrade the quality of the achievable entangled state, Eq. (46) may be seen as posing a fundamental limit to the entangling capabilities of a given system, a limit which would persist even in an idealized lossless scenario.

VI. CONCLUSIONS AND OUTLOOK
We analyzed stable and unstable states of a pair of atoms in a waveguide, finding that an entangled bound state exists for discrete values of the interatomic distance. This implies that an initially factorized atomic state can spontaneously relax towards a long-lived entangled state. By analyzing the poles of the resolvent operator, we have shown how to quantify the robustness of the entangled bound state to small variations in the model parameters, and how to identify the timescales that are crucial for the preparation of an entangled state by relaxation. While it has been pointed out that quantum computation may be achievable in waveguide-QED trough effective photon-photon interactions [49], focusing on the atomic degrees of freedom may also hold significant potential for applications in Quantum Information [50]. Further investigation will thus be devoted to the analysis of manyatom systems [51][52][53][54], in which photon-mediated interactions could possibly produce stable configurations such as W states or cluster states. electric and magnetic fields on the surface S read with ∂/∂n denoting the normal derivative with respect to the surface. Transverse electric (TE) modes are characterized by E x = 0 everywhere in the guide and obtained by imposing ∂B x /∂n = 0 on the surface. On the other hand, transverse magnetic (TM) modes have B x = 0 identically. If the waveguide is rectangular, the boundary conditions for TE modes reduce to which limits the form of the longitudinal magnetic field to the real part of with m, n ∈ N 2 \{(0, 0)} and B 0 a constant. The integers m and n label the mode TE m,n . The dispersion relation with respect to the longitudinal momentum has the same form as a massive relativistic particle, with ω m,n (0) = v mπy , where the mass term is called the cutoff frequency of the mode, and v = (µ ) −1/2 is the phase velocity in the waveguide, assumed isotropic and nondispersive with magnetic permeability µ and dielectric constant . Since L y < L z , the TE 1,0 mode has the lowest cutoff frequency. It can be proved [48] that ω 1,0 (0) is also lower than the cutoffs of all TM modes. Thus, at sufficiently low energy the contribution of the higher energy modes can be neglected, and propagation occurs effectively in one dimension. The TE 1,0 mode is characterized by the following behavior of the fields E z = i ω 1,0 (k)L y B 0 π sin πy L y e i(kx−ω1,0(k)t) , (A7) with the other three components vanishing. These fields can be derived from the (transverse) vector potential A z = L y B 0 π sin πy L y e i(kx−ω1,0(k)t) .
This stucture can be simplified if one assumes that the dominant contribution to the integrals comes from the poles of F (k) ∼ A + (k − k 0 ) −1 + A − (k + k 0 ) −1 . Neglecting the corrections yielded by square-root branch-cut integration, one obtains which is used to compute the energy density for the resonant states.