Polariton dynamics in one-dimensional arrays of atoms coupled to waveguides

Photons strongly coupled to material systems constitute a novel system for realizing non-linear optics at the level of individual photons and studying the dynamics of non-equilibrium quantum many-body system. We give a simple physical polariton-picture of the dynamics of photons coupled to a one-dimensional array of two-level atoms. This picture allows a fully analytical description of the dynamics in terms of polariton scattering inside the medium and reflections of the polaritons from the edge of the array. We show that inelastic collisions, previously identified in small systems, also occur in infinite systems and are related to the existence of multiple bands in the dispersion relation. The developed theory constitutes an effective field theory for the dynamics, which can be used for studies of non-linear optics and many-body dynamics. As a specific example we map the system to the Lieb–Liniger model and show that a so-called Tonks–Girardeau gas of photons is a stable eigenstate of the system in the limit of many emitters.

Here, we develop a concise fully analytical description of the dynamics of photons coupled to a regular array of multiple two-level emitters, as illustrated in Fig. 1 a).This provides an easy to understand picture of the essential dynamics of these systems.Our theory provides a full description of the input and output of the system [Fig. 1 (b) and (d)] as well as the propagation and interaction of the excitations [Fig. 1 (c)].Furthermore we identify an inelastic scattering mechanism where the photons exchange energy, as illustrated in Fig. 1 d).Such } Figure 1.a) Equidistant two level systems with distance d are coupled to a one-dimensional waveguide with decay rates, to the right (left) ΓR (ΓL), and to the side ΓS.The latter leading to photons being lost entirely from the system.b)&d) Input/output is described by assuming the photons to individually couple to the polariton modes of a semi-infinite array.c) The scattering of polaritons is described inside the chain ignoring any influence from the boundaries.The scattering may contain inelastic scattering to different photons energies while preserving the total energy.inelastic collisions have previously been identified for fewemitter systems [24].Here we show that such effect also exists in bulk systems and elucidate the physics behind them.Our results set the stage for explorations of more advanced many-body phenomena in these systems.We discuss this for the particular example of Tonks-Giradeau gasses [46,47].System-The system of interest consists of leftand right-propagating waveguide modes with field operaters E L (z) and E R (z), respectively, with and c is the group velocity of light in the waveguide.The two level emitter states are described by Pauli operators, e.g.σ + µ = |e µ g| µ exciting the µ-th emitter, all with the same distance d to the neighbours.We consider the typical one-dimensional photon-emitter-interaction Hamiltonian under the rotating wave approximation, and Here, H tot has been transformed to the interaction picture absorbing the excitation energy of the emitters.All photon frequencies ω are thus considered relative to the resonance frequency.By considering potentially deviating couplings Γ L = Γ R we incorporate an arbitrary level of chirality in the atom chain [48].
Polaritons-A polariton is a quasiparticle emerging from the coupling of a photon to a two level system.These polaritons can most easily be described after integrating out the photon degree of freedom.Then the dynamics of re-emission and absorption from emitter to emitter is contained in the effective Hamiltonian [18,20,26] with k 0 being the resonance wave number.The rate Γ S of photons ejected out of the waveguide leads to an overall depletion of the polariton number.While this constant depletion rate can be a considerable limitation for experimental realizations, it does not qualitatively influence the calculations to come apart from introducing a constant loss rate.It will thus be omitted from this point onward.
For an infinite chain the eigenstates of the system are given by Bloch's theorem and the effective Hamiltonian (3) then leads to the dispersion relation plotted in Fig. 2.
Input-We first study the coupling of photons into and out of the atom array as illustrated in Fig. 1 b) and  d).Using the well-known input-output formalism [36] we have where we without loss of generality chose a rightpropagating input field approaching a semi-infinite array.
Taking the spatial Laplace transforms of Eq. ( 5) and the Fourier transform in time we arrive at the input relation for an incoming field The transmission probability of photons entering the medium or equivalently the probability of polariton modes in the chain being ejected at the boundary.These depend on k, since the frequency is dependent on k.For k close to k0 the polaritons are more detuned giving a higher probability to be transmitted.For k0d close to 0 or π we obtain a linear scaling of the reflection probability leading to a N −3 scaling of the longest living subradiant states (see text).For increasing chirality, the transmission |t(k)| 2 /v k increases monotonously due to the suppression of backscattering from each emitter.

and the Pauli operators in momentum space
Here, we have a factor f (k where k is the degenerate wavenumber for which ω k = ω k .
To obtain a probability density for being coupled into the atom array as a polariton, we have to correct for the group velocity The transmission probability |t(k)| 2 /v k becomes smallest for k → nπ, n ∈ Z (i.e. the points closest to resonance), where the photon is entirely reflected if Γ L = Γ R (so that v k /v k = 1) expanding the well-known mirror-like behavior for k 0 d = π [41,50,51] to general k 0 , see Fig. 2. On the other hand, for a completely chiral setup with Γ L = 0 we get r(k) = 0 and the photonic wavefunction enters the medium in its entirety.
Output-For the output scenario we start with the output relation [36] with E in (z N , t) = 0 for a semi-infinite atom array ending at position z N .Doing a similar transformation as above, we arrive at (see [49]) analogous to the input scenario, only rescaled by the group velocity.Also, we have an additional phase factor depending on the length L = d(N − 1) of the atom array.From here we can calculate the probability for a polariton to leave the array as a photon The output probability thus exactly matches the input probability in accordance with Helmholtz reciprocity.As a direct application, this result allows us to reproduce the previously reported N −3 scaling of the subradiant states for Γ R = Γ L [7,9]: The lowest energy eigenstates for an array with N atoms are found at k = πξ/N d with ξ ∈ N + .Wavepackets localized around a certain k will reach the boundaries at a rate v k /N d.Altogether, this leads to a total ejection rate of γ = v k (1−|r(k)| 2 )/N d and expanding the transmission t(k) to lowest order in k reproduces the results from Ref. [9].
Polariton scattering-The only ingredient missing for a complete picture of the two-photon dynamics is the scattering event shown in Fig. 1 c).To describe the scattering of two polaritons we change to center of mass and relative coordinates with absolute momentum K = (k 1 + k 2 )/2 and relative momentum coordinate q = (k 1 − k 2 )/2.The two-excitation momentum basis states can then be expressed as where |0 is the polariton vacuum state where all atoms are in the ground states |g .
To solve the scattering problem we apply the effective Hamiltonian (3) on |q, K to obtain [20] where the total energy ω q,K = ω k1 + ω k2 and a(q, p) = sin qd/[cos qd − cos pd].
To get rid of the second and third line in Eq. ( 12) and turn it into an eigen-equation we make a scattering ansatz |Ψ q,K = |q, K + t 1 | − q, K + t 2 | − q , K , where q is a degenerate momentum number fulfilling ω q ,K = ω q,K with |q | = |q| (see below).Requiring |k 0 − K, K and |k 0 + K, K to vanish from H sys |Ψ q,K , then give two equations determining t 1 and t 2 as a function of q, K, k 0 and Γ R , Γ L .The complete solutions to these equations are given in the supplementary material [49].
In contrast to scattering in free space, the scattering amplitudes depend not only on the relative momentum For energy regions with four-fold degeneracy (blue region) we have inelastic scattering with |t1| 2 ≤ 1.In d) instead of a minimum at qd = 0 we have a tiny local maximum in energy which secures a broad region of inelastic scattering.In the yellow regions we require complex q to solve the eigenequation (12) and the imaginary part is plotted as dashed-dotted lines.b)&e): imaginary part of the degenerate q into which the scattering takes place.Whenever this imaginary part is non-zero we have scattering into meta-stable resonance states with finite extension, while for Im{q d} = 0 we have a combination of elastic and inelastic scattering.The maxima of q in the density plots are actually weakly ascending poles resembling strongly localized resonance states [20].c)&f): The transmission probability |t1| 2 ≤ 1 for elastic scattering into modes with −q as a function of q and K.When Im{q d} = 0 the scattering only goes to localized resonance states, so that |t1| 2 = 1 is required to conserve probability [20], whereas |t1| 2 < 1 for Im{q d} = 0.
q but also on their mean momentum K.The nature of the scattering depends heavily on the degeneracy of the total energy ω q,K shown in Fig. 3 a) and d).By symmetry the opposite relative momenta q and −q always have the same energy ω q,K = ω −q,K , but in some regions of phase space [marked with blue in Fig. 3

a) and d)]
there is a four-fold degeneracy with real q, −q, q , and −q fulfilling ω q,K = ω q ,K .In this case the scattering has an elastic component with a probability |t 1 | 2 to scatter into | − q and an inelastic component scattering into | − q with a probability |t 2 | 2 v q ,K /v q,K ; we show in the supplementary material [49] that the continuity equation /v q,K = 1 holds.Here v q,K = ∂ q ω q,K is the relative group velocity and the probabilities are symmetric under q → −q or q ↔ q .
Physically an inelastic collision means that the photons will be going out with different individual momenta and thus also different energies.Experimentally such events could thus easily be measured and represent a peculiar form of four wave mixing, which could, e.g., be used for completely frequency conversion of two single photons by picking t 1 → 0, e.g. the dark blue areas in Fig. 3 f).In most cases this change of energy results in a sign flip of the energy (relative to resonance ω = 0) such that one of the photons switches to a different branch of the dispersion relation, see Fig. 2. In case of a non-chiral setup this is always the case as the two branches are either purely convex or concave, and we cannot have momentum and energy conservation without branch-hopping.For strong chirality, however, this convexity can be broken and inelastic scattering can take place within the branches, see the red dots in Fig. 2 which are part of the region of inelastic scattering inside the triangle of Fig. 3 e).
As opposed to the situation above, for some combinations of q, K [marked with yellow in Fig. 3 a) and d)] the two-polariton dispersion relation only shows a two fold degeneracy ω q,K = ω −q,K .This means that we cannot directly solve the eigenequation ( 12) with real relative momenta, and thus inelastic scattering.In this case, however, there exist a different degeneracy ω q,K = ω q ,K , where q becomes complex, see Fig. 3 b)&e).The wave number q has a positive imaginary part, which leads to a localized, normalizable wave function.Physical such states represent scattering resonances, where the two polaritons bind together for a finite time [20].In the regions, where we have scattering resonances, we only have a single outgoing states | − q and we therefore find |t 1 | 2 = 1 as required to conserve probability.
Whether inelastic scattering or scattering into resonance states takes place depends on the setup parameters k 0 , Γ R /Γ L and the incoming momenta K, q.In general there is a rather complex interplay between these different effects, which is visualized in Fig. 3 c) & f).
Many particle limit-An interesting application of the above results is to use them as a starting point for studying more complex many-body dynamics in a long chain of atoms.A particular example of such systems is the Lieb-Lininger model describing massive particles with short range interactions in one dimension [52].Here the dynamics is governed by a scattering length a.In the limit a → 0 the ground state of this model is a Tonks-Girardeu gas describing impenetrable particles and fermionization of bosons [46,47].Below we show how our result predict that a similar behavior can also be produced in our setup for a non-chiral system Γ R = Γ L = Γ 0 .
To resemble the situation in the Lieb-Liniger model, we are interested in a situation where we only have elastic scattering.The most promising approach is thus to stay as close to k → 0 (corresponding to K, q → 0) as possible, where we are dominated by elastic scattering, at the boundaries.Near K, q ≈ 0 we have |r(K1)| 2 |r(k2)| 2 ≈ 1 so that polaritons are trapped inside the ensemble facilitating the creation of a Tonks-Girardeu gas.c)&d) Exact real and imaginary part of t1, respectively, in our setup at K = 0 with k0d = 0.05π (red), k0d = 0.5π (purple), and k0d = π (blue).The dotted lines are the result of the short range interaction assumed in the Lieb-Liniger model.While this provides the best description for small k0 the broadest region of t1 = −1 is found for k0d → π (although still being restricted to very small q).All plots are for a non-chiral setup ΓR = ΓL.see 4 a).By expanding ω q quadratically at the origin we find that the polaritons in this region can be described as massive particles with an effective mass m e = 2 (1 − cos k 0 d) 2 /d 2 Γ 0 sin k 0 d.Furthermore in this region we have near perfect reflection at the boundaries, as shown in 4 b).This ensures that the polaritons will be trapped inside the system enabling the study of long term many-body dynamics.
To fully match the dynamics of the Lieb-Liniger model we also need to ensure that the scattering dynamics match the model.For K = 0 the solution is especially simple and only involves scattering into −q with an amplitude and t 2 = 0 independent of the degree of chirality.By expanding the amplitude (13) for small k we find the scattering length allowing a complete mapping to the Lieb-Lininger model.As opposed to the Lieb-Liniger model, however, the model we use here is fundamentally a lattice model.This difference is explored in Fig. 4 c)&d) where we compare the scattering amplitude with the Lieb-Lininger model.As seen in the figure the resemblance is best for small k 0 d.On the other hand, to maximize the region of phase space with t 1 = −1 consistent with Tonks-Girardeu-like behavior, it may be desirable to work in a regime with k 0 d → π.
In total the mapping to the Lieb-Liniger model presented above indicate that it is indeed possible to realize a Tonks-Giradeau gas in these system.The analytical theory derived here provides a firm starting point for a more detailed analysis of this fascinating possibility, but a complete analysis is beyond the scope of this paper.
Conclusion-We have described the dynamics of polaritons in atom arrays in the limit of many two-level emitters.Our concise analytical model allows for an precise grasp of the rich physics in the system, featuring an interplay of elastic and inelastic scattering as well as scattering into scattering resonances.The results we obtain allow a better understanding of previously obtained results for reflection on resonance [41] and sub-radiance [8,9].At the same time the analytical description of the scattering process paves the way for understanding complex many-body effects.We have outlined a first step towards this by mapping the system to the Lieb-Liniger model and the Tonks-Girardeu gas.This is, however, only one out of the many possibilities.In particular the analytical approach brought forward in this paper may be extended to more involved setups including three-level emitters with deviating chiral coupling [21] and atom arrays in two or three dimensions.Such extensions may widen the range of many body dynamics, which may be realized with these systems.Furthermore, the theory developed here may also be applied to enhance photonic quantum technologies [44,45].
where we used in the second line Laplace transform into frequency space together with a substitution ∂ k ω k = v k , and the third line is achieved via the residue theorem.In case of Γ R = Γ L = Γ 0 we have k (ω) = −k(ω) and which simplifies the infinitesimal energy range of emitted photons to Here it makes sense to compare to the N −3 decay rate of subradiant states: Assuming a chain length of L = N d and the lowest energy states k = ξπ/N d with ξ ∈ N + we obtain the emission rate of waves travelling back and forth from one end of a finite chain to the other as which is in agreement with results found by diagonalization [9].

Input
For the input scenario in the perfectly chiral case we can start with a semi-infinite emitter chain with z µ ∈ [0, ∞)  which is very close to the degeneracy condition (B1), leads to a quadratic equation with coefficients proportional to (B2) (the factor being cos[(k 0 + K)d] − cos[(k 0 − K)d]) and vanishes when Eq. (B3) holds.Thus for degenerate k the term (B8) always vanishes and the continuity equation is fulfilled at all degrees of chirality r = Γ R /Γ L .Whenever k is complex, the scattering involves meta-stable states with finite extension and |t 1 | = 1.

Figure 2 .
Figure 2. Left: Dispersion relation for different values of k0d (dashed/solid lines) and chirality, ΓR/ΓL = 1 (thick black curves) or 4 (thin red curves).The symbols mark inelastic scattering between incoming (open symbols) into outgoing (filled symbols) wave numbers.Right:The transmission probability of photons entering the medium or equivalently the probability of polariton modes in the chain being ejected at the boundary.These depend on k, since the frequency is dependent on k.For k close to k0 the polaritons are more detuned giving a higher probability to be transmitted.For k0d close to 0 or π we obtain a linear scaling of the reflection probability leading to a N −3 scaling of the longest living subradiant states (see text).For increasing chirality, the transmission |t(k)| 2 /v k increases monotonously due to the suppression of backscattering from each emitter.

Figure 3 .
Figure 3.The scattering of two polaritons for different relative and absolute momentum q and K with a)-c) k0d = π/2 and ΓR/ΓL = 1, d)−f) k0d = 0.9 π and ΓR/ΓL = 4. a)&d): the two-polariton dispersion relation ωq,K (solid lines) for a fixed Kd = −1 marked by dashed lines in b), c), e), and f).For energy regions with four-fold degeneracy (blue region) we have inelastic scattering with |t1| 2 ≤ 1.In d) instead of a minimum at qd = 0 we have a tiny local maximum in energy which secures a broad region of inelastic scattering.In the yellow regions we require complex q to solve the eigenequation(12) and the imaginary part is plotted as dashed-dotted lines.b)&e): imaginary part of the degenerate q into which the scattering takes place.Whenever this imaginary part is non-zero we have scattering into meta-stable resonance states with finite extension, while for Im{q d} = 0 we have a combination of elastic and inelastic scattering.The maxima of q in the density plots are actually weakly ascending poles resembling strongly localized resonance states[20].c)&f): The transmission probability |t1| 2 ≤ 1 for elastic scattering into modes with −q as a function of q and K.When Im{q d} = 0 the scattering only goes to localized resonance states, so that |t1| 2 = 1 is required to conserve probability[20], whereas |t1| 2 < 1 for Im{q d} = 0.

Figure 4 .
Figure 4. : a) Real part of the transmission coefficient for elastic scattering t1 for k0d = 0.95π, ΓR/ΓL = 1.In the center, we observe a broad region of pure elastic scattering with t1 = −1.b) Product of reflection probabilities |r(k1)| 2 |r(k2)| 2at the boundaries.Near K, q ≈ 0 we have |r(K1)| 2 |r(k2)| 2 ≈ 1 so that polaritons are trapped inside the ensemble facilitating the creation of a Tonks-Girardeu gas.c)&d) Exact real and imaginary part of t1, respectively, in our setup at K = 0 with k0d = 0.05π (red), k0d = 0.5π (purple), and k0d = π (blue).The dotted lines are the result of the short range interaction assumed in the Lieb-Liniger model.While this provides the best description for small k0 the broadest region of t1 = −1 is found for k0d → π (although still being restricted to very small q).All plots are for a non-chiral setup ΓR = ΓL.
A11) and directly change to momentum space, σ − k(t) = √ d µ e −ikzµ σ − µ (t), for which we need to know the time evolution of the annihilation operators σ − k .We use the full Lindblad master equation of the reduced system[18] c k = σ − k + r(k)σ − k we have ∂ t c k = −iω k c k .For this operator the input equation becomes∂ t c k (t) = − i Γ R dE in (z 1 , t) ∞ µ=0 e ik0zµ [e −ikzµ + r(k)e ik zµ ] − iω k c k (t) = − i Γ R dE in (z 1 , t)[f (k − k 0 ) + r(k)f (k − k 0 )] − iω k c k (t), (A14)and isolation of c k (t) leads after back and forth transformation between time and frequency (analoguously to the output calculation) toc k (t) = −i Γ R dE in (z 1 , ω k )[f (k − k 0 ) + r(k)f (k − k 0 )]e −iω k t ,(A15)The entrance probability for a photon into the chain can be calculated just like in the output case as|c k (t = 0)| 2 dk = Γ R |f (k − k 0 ) + r(k)f (k − k 0 )| 2 |E in (z 1 , ω k )| 2 ∂k ∂ω dω = 1 d F (k, k 0 )|E in (z 1 , ω k )| 2 dω, (A16) with F (k, k 0 ) = (1 − |r(k)| 2) if the setup is completely non-chiral.After sufficiently long time (after which the wave package in question is localized far enough away from the boundaries) the σ − k (t) in c k (t) do no longer meaningfully contribute and we can map c k → σ − k .