Long-lived entanglement of two multilevel atoms in a waveguide

We study the presence of nontrivial bound states of two multilevel quantum emitters and the photons propagating in a linear waveguide. We characterize the conditions for the existence of such states and determine their general properties, focusing in particular on the entanglement between the two emitters, that increases with the number of excitations. We discuss the relevance of the results for entanglement preservation and generation by spontaneous relaxation processes.


Introduction
The physics of quantum systems confined in one-dimensional (1D) geometries has recently attracted a lot of attention [1,2], and is motivating interesting theoretical and experimental research. The behavior of an excited atom coupled to a field is among the peculiarities of such systems: although decay by spontaneous emission occurs in free (three-dimensional) space, boundary conditions and artificial dimensional reduction alter the picture, enhancing or inhibiting (and sometimes hindering) decay. These effects have been extensively studied and observed in cavity-QED settings [3][4][5][6][7][8][9][10][11][12][13][14], where the spectrum of the electromagnetic field is discrete. It is much less trivial that similar phenomena occur in effectively 1D unbound systems, in which the field spectrum is continuous and photons are free to propagate in 1D space. Dimensional reduction can be implemented in a range of experimental platforms, that include cold atoms in tightly focused fields [15][16][17], photonic crystals [18][19][20][21], optical fibers [22,23], quantum dots in photonic nanowires [24,25], and superconducting qubits in integrated circuit waveguides [26][27][28][29]. Theoretical studies focused on the interplay between the spectral features of the field and the structure of the emitters [30][31][32][33][34][35][36][37].
The vacuum of quasi-1D fields and their coupling with quantum emitters (real or artificial atoms) can be engineered by properly adjusting the distance between the emitters and a perfect mirror at one end of the system [28,38,39]. However, the interplay between absorption, stimulated and spontaneous emission provides a quantum emitter with mirror-like properties [29,40,41]. Hence, a pair of emitters can confine the field in the region between them, yielding nontrivial bound states above the threshold for photon propagation, that can be exploited for their robust entanglement features [42][43][44][45][46][47].
The objective of this article is to study the possible stable configurations of two multilevel atoms placed in a 1D cavity. See figure 1. In analogy with the 2-level case, each atom behaves both as an emitter and a mirror, confining the photon field and giving rise to a bound state endowed with highly nontrivial entanglement between the two atoms and between the atoms and the field. The discussed effects are non-perturbative and enable entanglement generation by relaxation. We will first adopt a method of resolution that applies to the case of two identical atoms that can be modeled as truncated harmonic oscillators. The technique consists in solving the problem for full-fledged harmonic oscillators, endowed with infinite number of levels. Due to the rotatingwave form of the Hamiltonian and the ensuing conservation law, the evolution in a given sector will involve only a finite number of atomic excitations. Finally, we will consider a perturbative method that enables to extend the analysis to the case of general level spacings and couplings, yielding the conditions for which the decay rate vanishes at second order in the atom-field coupling constant.

Harmonic-oscillator model
We consider a pair of emitters consisting of distinguishable harmonic oscillators A and B, with the same characteristic frequency ω 1 , placed in a linear waveguide at a distance d. We assume that the photons coupled to the oscillators belong to a single nondegenerate transverse mode of the waveguide, with dispersion relation ω(k). Hence, at the zeroth order in the coupling constant, the oscillator frequency must be larger than the low-energy cutoff to enable propagation along the guide, and smaller than the other mode cutoffs to justify the assumption of coupling to a single mode. These conditions can be typically realized in a linear rectangular waveguide, where the dispersion relation of the lowest-energy mode reads w = + ( ) ( ) k k M 2 21 2 , with M inversely proportional to the longer side of the guide cross-section. However, we will keep the discussion as general as possible. In the dipolar and rotating-wave approximations, the Hamiltonian reads ò ò where b(k) and b † (k) are the photon field operators in longitudinal momentum space, satisfying the canonical commutation relation , while b J and † b J (J=A, B) are the canonical harmonic oscillator operators, satisfying JK . The real coupling function g(k) naturally decouples at high frequencies, such that ò w (see e.g. [47] for photon waveguides). The excited (number) states are created by acting on the vacuum ñ = ñ Ä ñ The Hamiltonian(2) commutes with the total number of excitations and does not mix different sectors, belonging to different values of  . Due to this conservation law and the robustness of our approximations [47], our analysis applies equally well to a pair of harmonic oscillators and to a pair of (N+1)-level atoms (with equally spaced levels) in a waveguide. See figure 1.

Bound states
The  = 0 sector contains only the vacuum ñ |0 . The existence of nontrivial atom-photon bound states in the  = 1 sector, occurring for discrete values of the interatomic distance and having no counterpart in more than one dimension, was proven in the simpler case of two-level atoms [47]. We shall prove that the presence of these bound states is fundamental in determining the properties of highly-excited sectors.
Since the Hamiltonian is quadratic in the field operators, it can be diagonalized by a proper linear combination of the bosonic field operators. In particular, consider the generic combination

bosonic annihilation operator and satisfies the canonical commutation relation
for some real E. Then it is immediate to see that f † b creates bound eigenstates of H from the vacuum. Moreover, applying it N times to the vacuum will create eigenstates of H, belonging to the sector  = N , with eigenvalues NE.
By (5), equation (7) is equivalent to If the system admits a solution with  w < ( ) E M k , then the amplitude f(k) is square integrable and the state is also normalizable, provided condition(3) holds. These solutions correspond to bound states below the threshold for photon propagation: therefore, they occur even in the case of a single excited emitter coupled to the waveguide.
The solutions above threshold (E>M) are more interesting and nontrivial, since in such conditions a single excited emitter would spontaneously decay by photon emission. In fact, this happens in most situations also for a pair of emitters. The relevant properties of these solutions can be determined by a number of non-perturbative and general arguments, that can be eventually specialized to specific dispersion relations and to the case of small coupling. Normalizability of f(k) requires that the poles k i in(10), solutions to E=ω(k), are compensated by zeros in the numerator: thus, i must vanish at all poles. Let us exclude the former possibility, that does not depend on the atomic state, and focus on the latter: since for all i, in order to obtain a nontrivial solution, the poles must be constrained by the conditions (k i −k j )d=2πn, with n an integer [47,48]. Plugging these results into equations (8)-(9) must yield a real solution E, which does not depend on the choice of the pole k i . This sets a strong limitation to the possibility of bound states above threshold. A solution can in principle be found when both ω(k) and g 2 (k) are symmetric in k, and ω is increasing in | | k , in which case = -=k k k 1 2 and the energy w = (¯) E k of the bound state must satisfy the coupled equations The conditions on the existence of a bound state in this case imply f f = - This quantity, which is an exact result as well as(12)-(13), differs from unity by an order l º   g 2 2 in the perturbative regime, but can become small for w (¯) k close to the threshold M [47].
In order to clarify the effect of a finite threshold for propagation, let us consider the case of dipolar coupling to a massive photon (e.g. a guided mode), with the form factor derived from electrodynamics: can generally exist even in the case ω 1 <M, if the distance d is large enough to make w (¯) k of the same order as the interaction integral in equation (12). Actually, in the plot in panel (b), corresponding to the intermediate value ε=0.5 M, resonant bound states with ω 1 <M, corresponding to n=1, 2, 3, can be observed for sufficiently large d. Finally, the case ε=M shows the possibility that the curves extend to very small values of ω 1 , at which the rotating-wave approximation is no longer valid, also crossing the line ω 1 =0.
Once the solution curves, that relate characteristic frequencies and distances for each n, have been determined, it is possible to compute the atomic probability p at by equation (14). In figure 3, this quantity has been represented, as a function of ω 1 /M, in the case n=1 for three different coupling constants. It is worth observing that the value of p at is close to one only in the perturbative regime, namely for a sufficiently small coupling constant and a large enough characteristic frequency, which requires smaller interatomic distances to form the bound state.

Entanglement
We can now study the entanglement properties of the bound states, that occur, at fixed characteristic frequency ω 1 , whenever the distance d between A and B is such that the system (12)-(13) is satisfied for some n. The bound state can be expanded as  where we have introduced the normalized states of the emitters It is also convenient to introduce the antiresonant combination associated with a state that is orthogonal to both the states associated with b res and b f . The state ñ |N , characterized by a fixed number of total excitations, has components in the emitter excitation sectors ranging from  = 0 at , correlated with the presence of N photons, to  = N at , with no photon in the guide. The state, defined by the above expressions, can be exactly determined for any value of the parameters, after numerical computation of the atomic probability as in (14) for a given resonance. An interesting and relevant case is represented by the perturbative regime, in which  p 1 at and the bound state is dominated by the contribution from  = N at , namely In this regime, the bound state is thus characterized by a large projection on the state with maximal entanglement between atoms. Notice that the state y ñ -| ( ) AB N 1 , that is dominant in -ñ |N 1 , appears in the O(Nλ 2 ) term of equation (21).
Two comments are in order. First, the state belongs to the  = N sector and never leaves it, under the action of the Hamiltonian(2). Thus, although the preceding analysis has been done for two harmonic oscillators, it is still valid for two (N+1)-level systems in the waveguide. This point will be formalized in the final part of the article. Second, the reduced density matrix of the two emitters where H c , that commutes with b f , accounts for the continuous spectrum and the spontaneous decay of excited states orthogonal to ñ |N . Here we have assumed that there is only one resonant bound state in the spectrum, which is valid in the small coupling regime. In the one-excitation sector (namely two 2-level atoms), it has been demonstrated [47] that the presence of the bound state enhances the decay rate of the orthogonal states towards the configuration in which the emitters are both in their ground states and a photon propagates in the waveguide.
If the emitters are placed at infinite distance from one another, an initial state close to y ñ Ä ñ | | ( ) vac AB N with N>0 would rapidly decay to an orthogonal state with a smaller number of emitter excitations. By contrast, in a waveguide, with d and ω 1 close to the resonance conditions, such a state will be left almost invariant by the Hamiltonian evolution, with a dressing (slight, in the perturbative regime) due to the imperfect superposition of the initial and the bound state (i.e. < p 1 at ). Therefore, on timescales smaller than the waveguide losses, entanglement is preserved without imposing constraints or external control.
Another interesting application is related to the decay of the unstable component of an arbitrary initial state. An initial state r r = Ä ñá | | vac vac , which has a nonvanishing projection on ñ |N , will relax towards an excited atomic state. The probability of having maximally entangled atoms is given by Therefore, the atoms, starting from a general density matrix, which can be e.g. factorized, eventually reach a highly entangled state with N excitation with finite probability. The typical relaxation time of the reduced atomic density matrix towards its asymptotic state is determined by the lifetime of the component y ñ -| † ( ) b AB N anti 1 , and coincides with the lifetime of a single atomic excitation created in the antiresonant state [47]. The described entanglement generation strategy is similar to the Hamiltonian generation of entanglement, in which an initial (factorized) state is let to coherently evolve until it reaches an entangled state. Such a procedure can be applied in the case of the bound states well below the threshold for photon propagation. However, the Hamiltonian nature of the evolution in the AB Hilbert space yields oscillations, which implies that the evolution must be stopped at a proper time to obtain the desired state. This drawback is absent in entanglement generation by relaxation, in which the final state is approached asymptotically. This technique does not require energy pumping into the system, since a constant entanglement is reached after an initial transient [47].
Let us go back to the discussion of the entanglement properties of y ñ | ( ) AB N . The emitter state in the sector  = 1 at (two 2-level atoms) has been extensively studied, and is particularly interesting since it corresponds to one of the two Bell states, according to the sign of -( ) 1 n (see equations (13)- (17)). It is thus maximally entangled in the  = 1 at subspace. The bound state = ñ |N 2 , relative to the pair of 3-level atoms in figure 1, reads as  ¥ N , where Γ is the Euler gamma function. Strictly speaking, this quantity measures entanglement between A and its environment (B + field). However, since the state of the field is quasi factorized at small coupling, it is also an approximate measure of entanglement between the two emitters A and B. On one hand, purity(28) scales more slowly than the minimal value (N+1) −1 in the sector, corresponding to maximally mixed reduced density matrices. On the other hand, this result is consistent with the minimal purity for states whose reduced density matrices are effectively approximated by the superposition of O(N 1/2 ) states.
It is also possible to determine the entanglement properties of coherent and incoherent superpositions of the bound states ñ |N . For example, one can consider the 'pseudothermal' state AB E E b b th whose reduced density matrix is the thermal average of (26), yielding the purity å å with Nth the average excitation number. Another interesting case is the coherent state . Since it is a product of coherent states, , there is no entanglement between the two atoms. In particular, in the small-coupling limit, the atom density matrix r a AB is dominated by the projection on , has no effect on the existence and properties of state ñ |N . On the other hand, the bosonic commutation relations satisfied by all operators, which is reflected in the bosonic character of the b f operator, are fundamental in the derivation of the commutator equation (7), and in the subsequent reduction of the problem to the determination of f A , f B , f(k), and E. If the oscillators are truncated atN , then the existence of a bound state in the sector  = >N N does not follow automatically from the presence of a bound state in the lowest-excited  = 1 sector. Summarizing, in bound states, the excited levels of each oscillator must be in a sufficient number to absorb all the photons in the state. And, of course, considerations on the pseudothermal and coherent states are valid, only approximately, if the average occupation number is not close toN .

Stable states in the general case
We now relax the assumptions of equally-spaced levels and oscillator-like interaction. Unlike in the previous sections, where the analysis could be perfomed non-perturbatively, we will work at second order in the atomfield coupling.
Let us consider a Hamiltonian describing two identical (N+1)-level atoms A and B with the same spectra (with ω 0 =0), placed at a distance d from each other and coupled with the guided field. We assume that the interaction with the field can induce transitions between neighboring levels ñ | j I ( j=0, 1, K, N with I=A, B), with an arbitrary couplings: , generally decays into N photons. We have shown in the previous sections that, in the harmonic-oscillator-like case (ω j =j ω 1 and = ( ) ( ) g k j g k j 1 ) an exact eigenstate with nonvanishing amplitude in the  = N sector can exist. In all other cases, one can nonetheless identify the most stable combinations (33) by analyzing the decay rate at second order in the atom-photon coupling [50]. This quantity can be immediately computed from the action of the interaction Hamiltonian on the general state, yielding The integral can be computed by determining the set of momenta Let us assume that dispersion relation and form factors are symmetric, and that the argument of the delta function in equation (35) has only two symmetric solutions ¯(| | ) k c 1 2 . Then, the decay rate reads Let us consider two resonance cases, with =  (¯) kd cos 1. If = g g 1 2 1 2 , one recovers the truncated harmonic oscillators analyzed in sections 3-4, with * =  c 1 2 1 . In the interesting case of equal couplings g 1 /g 2 =1, stationarity and normalization lead to the solution * =  c 1 3 1 , yielding a uniform minimizer. In both cases, the decay rate is also vanishing at the minimum.
We finally comment on the robustness with respect to small variations of the level spacings. The minimal decay rate can be written as a function of the bare energy levels as , with *  c the minimizing sequence of coefficients. Generally, the minimal rate varies linearly with small displacements of the energy levels. However, a relevant exception is represented by the cases, like those mentioned before for N=2, in which w G =  ( ) 0 N . Since Ḡ N is nonnegative, if it is also a smooth function of w  , then its gradient with respect to w  must vanish as well. In such cases, the decay rate increases quadratically in the variation of the bare energy levels.

Conclusion
We have investigated the existence and properties of the stable states of a pair of N-level atoms in a waveguide, discussing the possibility to generate and preserve robust entanglement between the two atoms. As a case study, we have considered identical distinguishable emitters with uniform level spacing. The generalization to the presence of (small or large) asymmetries in the excitation energy and in the coupling, or uneven levels, has been considered in section 5. Among possible interesting applications, we mention the use of a single two-level emitter as a dynamical probe of the state of the field, along the guidelines discussed in [28]. These ideas can be generalized, by using (entangled) atomic multilevel pairs, paving the way to unprecedented possibilities and possibly super-resolution.