Atomic waveguide QED with atomic dimers

Quantum emitters coupled to a waveguide is a paradigm of quantum optics, whose essential properties are described by waveguide quantum electrodynamics (QED). We study the possibility of observing the typical features of the conventional waveguide QED scenario in a system where the role of the waveguide is played by a one-dimensional subwavelength atomic array. For the role of emitters, we propose to use anti-symmetric states of atomic dimers - a pair of closely spaced atoms - as effective two-level systems, which significantly reduces the effect of free-space spontaneous emission. We solve the dynamics of the system both when the dimer frequency lies inside and when it lies outside the band of modes of the array. Along with well-known phenomena of collective emission into the guided modes and waveguide mediated long-range dimer-dimer interactions, we uncover significant non-Markovian corrections which arise from both the finiteness of the array and through retardation effects.


I. INTRODUCTION
The interaction between quantum emitters and a structured reservoir leads to a rich phenomenology. It has recently attracted renewed interest [1] due to experimental progress in novel experimental platforms such as superconducting waveguide quantum electrodynamics (QED) [2][3][4][5], cold atoms near nanofibres [6,7] or photonic crystal waveguides [8][9][10], and quantum dots [11][12][13]. The dynamics of emitters coupled to such nonconventional reservoirs is expected to exhibit several distinctive features such as collective super and subradiant emission into the reservoir [9,14], long-range dipoledipole interactions mediated by the reservoir [15][16][17][18], and non-Markovian effects [19][20][21]. The observation of many of these phenomena is, however, challenging even with the unprecedented level of control achieved today in several experiments. These difficulties arise from both imperfections in the fabrication of these devices as well as from the unavoidable absorption of photons in the waveguide (see [1] and references therein).
Recently, Masson and Asenjo-Garcia have proposed the concept of atomic waveguide QED [40]. Inspired by waveguide QED, where optical emitters are coupled to * david.castells@mpq.mpg.de a waveguide, they considered the case in which atoms are coupled to a one-dimensional subwavelength atomic array which plays the role of the waveguide. The motivation is that an atomic waveguide in free space is a conceptually simple, clean optical medium that allows, in principle, to eliminate any intrinsic internal losses or imperfections which affect conventional waveguides or photonic crystals, and may feature very low disorder. The observation of characteristic waveguide QED phenomena in an atomic waveguide setup poses, however, several fundamental challenges [40]. On the one hand, efficient coupling between external "impurity" atoms and the atomic waveguide is hindered by free space decay due to the presence of superradiant (bright) modes of the array. On the other hand, the dynamics of the emitters shows signs of non-Markovian effects, which ultimately spoil some of the interesting phenomena of conventional waveguide QED. Some of these difficulties may be overcome by reducing the interatomic separation of the waveguide and by placing the impurity atoms extremely close to the array [40,41], but it is challenging to achieve the required deep subwavelength regime experimentally.
In this article, we propose a modification of this setup to mitigate these problems. Specifically, we propose (i) to use atomic dimers -a pair of closely spaced atoms -as effective two-level emitters coupled to the atomic waveguide, and (ii) to control the dimers' linewidth with a Raman transition. We show that a dimer behaves as an effective two-level system formed by its ground state and its anti-symmetric state. The anti-symmetric state features a reduced coupling to free-space modes and, at the same time, an increased coupling to the array's guided modes. Additionally, we show that by controlling the decay rate of the dimer atoms via a Raman transition, it is possible to recover a Markovian regime for the dynamics of dimers coupled to an atomic waveguide. We derive simple models for our setup that predict collective emission from the impurity dimers into an array's subradiant mode, and coherent long-range interactions be-tween impurity dimers mediated by the array. We verify both observations numerically for the case of an atomic array with interatomic separation of quarter-wavelength, a regime where the simpler case of single atoms coupled to an atomic waveguide is hampered by free-space decay and non-Markovian effects [40]. We also study the effects and different non-Markovian behaviors arising from the finiteness of the chain and the reduced group velocity at the band edge. Our results show a clear advantage of using dimers over atoms and highlight a promising route towards observing non-Markovian waveguide QED physics, which has received considerable attention recently [21,[42][43][44].
This article is organised as follows. The theoretical model, as well as details on the dimer-array coupling and on the Raman transition used to control the dimer's linewidth are described in Sec. II. The physics of dimers coupled to an atomic array is described in Sec. III (Sec. IV) for the case of dimer's frequency lying inside (outside) the band of guided modes of the array. We discuss possible generalisations of the case presented here and draw our conclusions in Sec. V. We leave additional details on the derivation to the appendices.

II. SYSTEM DESCRIPTION
We consider a one-dimensional atomic array of N atoms with resonance frequency ω 0 and lattice spacing d, and n impurity atoms with resonance frequency ω imp 0 placed at a distance h from the chain [see Fig. 1(a)]. We assume all atoms to be polarised along the z-axis, which coincides with the direction of the atomic array. The effective dynamics of an atomic ensemble coupled to a continuum of quantized electromagnetic modes was first derived in Ref. [45]. In the electric-dipole and rotatingwave approximation, and in the Markovian regime, the photonic environment can be eliminated yielding an effective non-Hermitian Hamiltonian ( = 1) that, together with stochastic quantum jump operators, describes the dynamics of the system [46][47][48]. Here, Γ 0 is the emission into free space of a single atom (we assume that the chain and impurity atoms have the same decay rate Γ 0 ), k 0 = 2π/λ 0 = ω 0 /c, and G 0 (r i , r j ) is the free space Green's tensor describing the field at atom i generated by atom j at an energy ω 0 , For the terms i = j, we use −3πG zz 0 (r i , r i ) = −i/2, implicitly including the Lamb shift in the energy of the emitters. In the absence of a driving term and within the single excitation sector, as we consider here, the dynamics of the system is completely characterised by the non-hermitian Hamiltonian Eq. (1). In this case, in fact, a quantum jump prepares the system in the collective ground state which does not evolve under the action of Eq. (1).
We can separate Eq. (1) into three terms including the chain, the impurity, and impurity-chain interaction parts of the Hamiltonian, The chain Hamiltonian can be written aŝ In the single excitation regime,b k can be taken to be bosonic annihilation operators. Closed expressions for modes energy J k and decay rate Γ k are derived in Ref. [31] for the case of an infinitely long array. In the following, we are interested in the experimentally relevant case of an array with open boundary conditions. In this case, an analytical exact expression for the single excitation eigenmodes of Eq. (4) is not available. However, when the atomic array is sufficiently long (N 1), the eigenmodes ofĤ chain can still be understood as spin waves of a definite quasi-momentum k, where the value k corresponds to the point in reciprocal space where the eigenmode wavefunction is peaked [22,23,31]. In this limit, an accurate ansatz can be provided for the single where k ν d = πν/(N + 1) with ν = 1, 2, . . . , N . Accordingly, J k and Γ k in Eq. (5) do not have analytically expressions, but are well approximated by the expressions for an infinite chain whenever N 1. In the following, unless otherwise specified, we will always refer to the case of finite chains with open boundary conditions, as this is the most relevant case for the experimental realisation of atomic waveguide QED. For d < λ 0 /2, there exist single excitation eigenstates ofĤ chain with a quasi-momentum k > k 0 lying outside the light cone of free space electromagnetic modes. These collective states, which are intuitively understood as excitations propagating along the array, have been shown to decay at a rate ∼ Γ 0 /N 3 due to scattering of the field through the ends of the array [31,49], and are referred to as the collective subradiant modes of the atomic array. The dimer is aligned with the center, and the single atom with the closest site to the center of an array with N = 100 and d = h = λ0/4. The dimer anti-symmetric state interacts destructively (constructively) with the low(high)-k modes. Hence, in comparison with a single atom, the coupling of the dimer to the guided modes is stronger, while the coupling to the superradiant modes is highly suppressed. The finiteness of the chain adds standing-wave envelopes to the Bloch states, which explains the difference between even and odd modes represented with dark and light lines, respectively. (c) Atomic level scheme of the setup proposed in this work. The two-level impurity atom is realized with a threelevel Λ atom, in which the transition |g1 → |e is driven by a laser field with detuning ∆. The detuning and spontaneous emission rate are much larger than any other energy in the system and, thus, |e α can be eliminated resulting in an effective two-level atom |g α 1 -|g α 2 (see text for details). The decay rate from |e α to |g α 1 is assumed to be much slower than Γ0. The chain modes are diagonalized in k-space, with energy J k and free-space decay Γ k .
The second term in Eq. (3), which includes the impurity atoms and interactions among them, readŝ Here we defined (g ij − iγ ij /2) = −3πΓ 0 G zz 0 (r i , r j ), and replaced the Pauli matrix of the impurity atoms with bosonic operators (σ i ge →ĉ i ) as we restrict to the oneexcitation dynamics [23]. According to the definition above for the Green's tensor, when i = j we have The last term in Eq. (3) contains the interactions between chain atoms and impurity atoms. Using the definitions above, the interaction Hamiltonian between the impurity atoms and the eigenmodes of the chain readŝ Expanding Eq. (2) in spherical coordinates, the coupling of a single impurity atom at position r i to a chain mode with quasi-momentum k reads [see Appendix A] where ξ k (z i ) is given in Eq. (6) and Here, κ m (k) = (k/k 0 + mλ 0 /d), and H (1) 0 is the Hankel function of the first kind and zeroth order. Note that for k > k 0 , Eq. (10) is purely real, which means that the coupling of an impurity to those Bloch modes is coherent. The coupling Eq. (10) between a single impurity and a mode k in the array is plotted in Fig. 1(b). For the case of an infinite array, the coupling can be obtained from Eq. (10) after the substitution ξ k (z i ) → e −ikzi / √ N . Since achieving small d is increasingly challenging experimentally, we consider here the case of d = λ 0 /4. For this parameter choice, the free-space decay of the most superradiant state, Γ k=0 , is comparable to the width of the band J k in Eq. (5). The dynamics of an emitter with energy lying within the band of the array is thus dominated by the dissipative resonant interaction with the array's superradiant modes. We overcome this problem using as impurity, instead of single atomic emitters, the collective excitation of two neighbouring atoms, which we name atomic dimer.

A. Atomic dimers
We consider a dimer to be formed by two neighbouring impurity atoms at positions r i = (h, 0, z i ) T and r i+1 = (h, 0, z i+1 ) T . For convenience, we label the two atoms forming a dimer "a" and "b" with positions r a i = r i and r b i = r i+1 , respectively. We represent the collective single excitation of an atomic dimer by the bosonic operator and its Hermitian conjugate, whereĉ a † i+1 . For λ = −1 (λ = 1), Eq. (11) creates an antisymmetric (symmetric) excitation which predominantly couples to the sub(super)-radiant array modes exploiting their short(long)-wavelength nature. The coupling to the chain modes of the states of a dimer centered at ρ i = z · (r a i + r b i )/2, with atoms aligned to the array atoms as depicted in Fig. 1 This expression is similar to Eq. (10) with an additional factor of √ 2 and a quarter-sine-wave envelope. The anti-symmetric configuration (λ = −1) eliminates the undesired dissipative coupling at small k due to a destructive interference of the interaction of each atom with the chain modes, as can be seen in Fig. 1(b). In the case of the dimer, coupling to half of the modes becomes zero at the center of the chain due to the particular symmetry of the collective dimer states, as indicated in Eq. (A6) and Eq. (A7). Unless otherwise specified, we set h = d.
The dimer states in Eq. (11) are the eigenstates of Eq. (1) with two atoms. Their eigenvalues are for a separation r ab between the two atoms. For r ab < λ 0 /2, as is the case for subwavelength arrays in one dimension, the anti-symmetric dimer state has a reduced linewidth as compared to the single atom. In particular, for d = λ 0 /4, the decay rate Γ 0− Γ 0 /4. We have considered alternative setups, such as a 2x2 quadrupole in an anti-symmetric configuration, but the reduction of freespace decay is hindered by the reduction of its effective coupling to the chain modes [see Appendix A 1]. We have also investigated the case in which the dimer atoms are separated by a smaller distance, for which we predict only a small improvement in performance [see Appendix A 2]. When the emitter's resonance energy is resonant with the subradiant part of the band, or it lies outside of the band, it is possible to adiabatically eliminate both the chain's superradiant part as well as the symmetric state of the dimer. The dynamics of the system can thus be modelled as an effective two-level system -formed by the dimer's ground and anti-symmetric states -coupled coherently to a set of subradiant modes [see Appendix B]. In the rest of the text we often refer to the resulting effective two-level system simply as dimer. Note that the expressions Eq. (12) and Eq. (13) are derived here for the case of a single dimer coupled to the array. The same derivation applies to the case of multiple dimers if these are sufficiently distant, such that the interaction in Eq. (7) between atoms belonging to different dimers is negligible.
Despite the reduced effective two-level system decay rate Γ 0− , the dimer's linewidth is comparable to the bandwidth (BW) of the array for d = λ 0 /4, leading to strong non-Markovian effects. To enter the Markovian regime, we can further reduce the dimer's linewidth using a Raman scheme as depicted in Fig. 1(c), and as we now describe.

B. Raman Transition
In the Raman scheme, the dimer atoms are initialized in additional metastable levels |g a,b 1 , and driven into the excited state of the dipole transition included in Eq. (1), |e i =σ i eg |g i 2 , by a laser with Rabi frequency Ω, as illustrated in Fig. 1(c). The driving frequency, ω R , is detuned with respect to the energy difference between the two states by ∆. The dimer atoms are excited without affecting the chain atoms by placing the chain on a node of the laser field. To describe the dynamics of the system including the Raman transition on the dimer atoms, we use the Hamiltonian in Eq. (3) and include the energy of the levels |g a,b 1 , and the interaction terms due to the driving laser, Ω 2 |g α 1 e α | e iω R t + h.c. . To remove the time dependence of the new interaction term, we move to a frame rotating with ω R , yieldinĝ We use the evolution under the Hamiltonian in Eq. (14), withĤ being the real-space expression of Eq. (3), for all numerical results presented in this paper, which we use to benchmark the predictions obtained with analytical derivations and additional approximations. For ∆, Γ 0 Ω, g ab , γ ab , g i k , γ i k , the states |e a,b remain weakly populated at all times, which we eliminate to second order in perturbation theory [see Appendix C]. The resulting dynamics can be approximated by an effective Hamiltonian that includes effective two-level impurity atoms |g α 1 -|g α 2 . For this, we define new creation operators for the dimer states, The population of all other states is negligible. The loss of total population is due to the dimer's free space decay. The unshaded region corresponds to the Markovian regime, in which the population is predicted by Fermi's golden rule (FGR). The shaded region corresponds to the non-Markovian regime. (b) Purcell enhancement extracted from the population dynamics. Γ1D is obtained by fitting an exponential decay to the corresponding region of the dimer's population evolution, and Γ tot 0 , which includes all sources of decay, by fitting the evolution of the total population in the one-excitation subspace. The solid lines show the prediction using FGR. Note the difference in performance between using a dimer or an atom. (b) and (c) The black arrow indicates the quasi-momentum of the chain mode with which the dimer in (a) is at resonance. The deviation of the points at higher k from the predictions is due to the non-Markovian dynamics close to the band edge. (c) Time until the dynamics enter the shaded region in (a), measured from the numerical evolution as the time at which the population at the dimer differs from the Markovian prediction by a 6%. The solid line shows the prediction using the length of the chain and the group velocity at k to calculate the time at which the excitation in the corresponding guided mode reaches the dimer after being reflected at the ends of the chain. (d) Dispersion relation of the chain. The shaded regions indicate the energies for which the emission rate of a symmetric state of n = 3 dimers into the chain shows n Γ1D superradiance for dimer-dimer separations of L = 5, 7, 9, and 11d. Since the possible spacing between the impurities in our setup is a multiple of the array spacing, d, only a discrete set of energies give rise to such superradiance. The dashed lines indicate the predicted energy at which superradiance is observed, which corresponds to k = (L − d)/L · π/d. The inset displays one of the plots used to obtain the green shaded regions. For L = 9, it shows the ratio between the decay rate into the chain modes extracted from the numerical evolution with n = 3 and the expected decay rate of a single dimer into the chain modes using FGR. Parameters for all the plots are N = 500, d = λ0/4, ∆ = 8 Γ0, and Ω = 0.2 Γ0.
We can further simplify the resulting Hamiltonian by eliminating the dimer symmetric state and the superradiant modes of the array. Using ∆ Γ 0 , the effective coupling between chain and dimer becomes coherent. Although we derive the effective Hamiltonian considering one dimer, the extension to multiple dimers in the case in which their interaction through free space is negligible is straightforward. Hence, for ∆ Γ 0 Ω, g ab , γ ab , g i k , γ i k and distant dimers, we observe the dynamics of Eq. (14) to be well approximated by the effective Hamiltonian with Note that g k 2 /Γ 0− is indepen-dent of Ω and ∆ as long as ∆ Γ 0 Ω, g ab , γ ab , g i k , γ i k is fulfilled. Obtaining smaller Γ 0− , however, makes the Markovian regime accessible, since Γ 0− /BW ∼ (Ω/∆) 2 , while g k can be modified with other system parameters, such as tuning the resonance energy of the emitter with respect to the band, as we discuss in section III. In conventional waveguide QED, different physics arises when the emitter energy is inside the band of guided modes compared to when it lies in the band-gap. In the following, we consider these two regimes separately in the context of atomic dimers coupled to an atomic waveguide. For simplicity, except when the finiteness of the chain is explicitly relevant, we use the infinite array form of ξ k (ρ i ).

A. Markovian regime
When the resonance energy of a dimer lies within the sub-radiant region of the chain's band, there is a coherent transfer of population between the dimer and the guided modes of the chain. In the Markovian regime [see the white region in Fig. 2(a)], the transfer of population can be modeled as a plane wave emitted into the resonant chain mode k, where Γ 1D is the effective decay rate of the dimer excitation into the chain. A large Purcell factor P = Γ 1D /Γ 0− corresponds to the desired regime of predominant decay of the emitters into the chain modes. The decay rate Γ 1D obtained from numerical simulations agrees with the prediction using Fermi's golden rule (FGR), In Fig. 2(b), we plot P as a function of k as extracted from the numerical evolution, and compare it to the prediction using FGR. We compare the case of a dimer and of a single atom, showing that the former allows for a substantial improvement in P with respect to the latter. The coupling g k increases with k, while ∂ k J k decreases and becomes zero at the band edge. Hence, larger k are favorable and lead to a divergence in the Γ 1D predicted by Eq. (17). The Markovian assumption, however, breaks down at large Γ 1D before reaching such divergence, as observed from the deviation between the model and the numerical results in Fig. 2(b,c). The value of P at a certain k depends on the parameters d and h, but not on Ω or ∆. The value of k at which the Markovian approximation breaks down, however, depends on d, h, and also on Ω/∆: smaller Ω/∆ reduces both the coupling strength and the linewidth of the dimer. Hence, the evolution stays Markovian at higher k, for which P is larger. This improvement in P is achieved at the cost of slower dynamics. The presence of a 1D bath allows for different dimers along the chain to have a finite probability to interact with the photon emitted by the originally excited dimer. Such interaction leads to a constructive interference if n dimers are prepared in a symmetric statê i− and placed in an atomic mirror configuration: |ρ i − ρ j |k = 2πq, with q ∈ Z. In this state, the decay rate is enhanced by Γ sym 1D = nΓ 1D , while emission into free space is unaltered. This type of superradiance is also observed in conventional waveguide QED [14]. In Fig. 2(d), we show this feature with three dimers. Specifically, we plot the band J k and indicate which energies correspond to values of k that satisfy the atomic mirror condition. The shaded regions indicate the resonance energy of the dimers for which we observe a decay rate into the chain a factor of three larger than the FGR prediction, Eq. (17). We extract these regions from plots like the one that we show as an inset, in which we obtain Γ 1D as a function of k. This collective emission can be exploited to achieve larger Γ 1D /Γ 0− .

B. Non-Markovian regime
The breakdown of the Markovian regime due to the diverging density of states at the edge of the band discussed before is also observable in waveguide QED setups [50,51]. In our system, we observe another type of non-Markovianity due to the finite chain length. Note that the results discussed above are independent of the number of chain atoms, N , except that the decay rate Γ k of the subradiant chain modes is slightly larger for smaller N . However, the length of the chain determines the time for which the dynamics stays Markovian [see the shaded region in Fig. 2(a), and Fig. 2(c)]. For a dimer aligned with the center of the chain, the outgoing plane waves return to the dimer after a characteristic time-scale where v g is the mode's group velocity, due to reflection at the ends of the chain. We understand these non-Markovian effects as the retarded back-action via the reflected electric field of the emitter. Mathematically, these non-Markovian effects originate from the discrete spectrum of the atomic waveguide [52], as the dynamics of the emitter at long times is able to resolve the energy difference between two chain modes.
The evolution of the chain and dimer population in the non-Markovian regime is highly dependent on the resonance energy of the dimer, which makes the system highly tunable. We can distinguish two particular cases (i) when the dimer is resonant with an anti-symmetric chain mode and (ii) when its energy is resonant with a symmetric chain mode. The phase of the reflected wave at the dimer's position differs by a factor of π in the two cases. In the first case, the reflected wave is in phase with the dimer and leads to an enhanced emission [note the kink in Fig. 2(a)]. In the second case, the reflected wave accumulates a difference in dynamical phase. Since the coupling rate of the dimer with a symmetric mode,ν, at the center of the chain is zero, the dimer energy effectively lies between two anti-symmetric modes. For large N , the energy difference with the two modes is approximately equal, ∆E J kν±1 − J kν , and τ π/ J kν±1 − J kν , accumulating a phase difference ϕ = ∆E τ π. The reflected wave is, thus, out of phase and leads to an increase of population in the dimer [see Fig. 6]. The interaction with the reflected field can also be understood as the dimer interacting with its mirror image [53,54], in which the non-Markovian effects are a form of retarded Dicke super or subradiance, with the emitters having (i) parallel or (ii) opposite polarization, respectively [21]. The slow propagation of the guided modes, especially close to the band edge, enhances the retardation effects responsible of the non-Markovian behavior.

IV. EMITTERS IN THE BAND-GAP
If the dimer state resonance energy is located at the band-gap and the detuning with the band edge, strength to the corresponding chain modes, g π d , the emission into the chain is blocked. However, an atom-photon bound state with an exponentially decaying tail is formed [15]. The width of the tail scales as 1/ √ δ and, for sufficiently small δ, an overlap between distant atoms can be obtained [17,18]. Adiabatically eliminating the chain modes in Eq. (15), we obtain the following expression for the effective long-range coupling between distant dimers For small δ, the major contributions to the sum over k above are concentrated around the band edge. We thus approximate the band at the edge of the Brillouin zone as J π d (1−x) = J π d −A d x 2 , and g k g π d , since g k varies slowly close to k = π/d [see Fig. 1(b)]. Likewise, the decay rate of the most subradiant modes can be approximated by Γ π d (1−x) γ N x 2 [31], with γ N /Γ 0 1/N . For subwavelength arrays, γ N /A d < 1/N . For compactness, we usẽ With these approximations and in the continuum limit, Σ k → N d 2π dk, we obtain a closed form for the effective coupling between two dimers mediated by the guided modes of the array, with where we identify the length scale of the interactions discussed above, l = Ã d /δ. These expressions are valid as long as δ g ij eff , such that elimination of the chain modes is justified. For a sufficiently large g ij eff in relation to the free-space decay, Eq. (19) predicts Rabi oscillations between the two dimers with a fidelity dictated by the ratio g ij eff /Γ 0− . This ratio does not depend on the Raman transition parameters ∆ and Ω. However, a small Ω/∆ allows to reduce g ij eff and, thus, one can still fulfill δ g ij eff with a smaller detuning δ. Fixing = g ij eff /δ, we can rewrite Eq. (20) as Note that, as expected from the continuum limit, Eq. (21) is independent of N , since |g − k | ∼ 1/ √ N , with the exception of small corrections due to a finite Γ k (N ), that vanishes at large N . The dependency on ∆/Ω indicates that the effective coupling can be made arbitrarily large at the expense of slower dynamics. We consider a system of two dimers separated by a distance L. Initializing the system with the excitation in one of the dimers, the chain mediates Rabi oscillations between the two otherwise non-interacting dimers, as shown in Fig. 3(a) for a distance L = 14 d. We aim at maximising the fidelity of preparing the first dimer again in the excited state after a single Rabi cycle. The total error of the protocol, defined as the population lost to free space after one full cycle, comes from (i) the free space decay rate Γ 0− of the dimers and (ii) the finite linewidth Γ k of the guided modes. There is also a small dephasing contribution in the exponential envelope, which we neglect in our discussion. We measure the error from the numerical evolution as one minus the first relative maximum of population at the initial dimer. Another source to the measured error is (iii) the transfer of population to the chain modes. When Γ 0− Re[g ij eff ], the first contribution is approximated as πΓ 0− / Re[g ij eff ] ∼ − 1 /3 , which is minimized for large values of . The third contribution, instead, grows with , as the maximum population transferred to the chain during the dynamics can be shown to be upper bound by a function proportional to [see Appendix E 1]. The second contribution to the error is independent of , and has a value γ N /(2A d ) ∼ 1/N , which sets a lower bound to the error independent of the ratio Ω/∆. In Fig. 3(b), we plot the error for dimers interacting with chains of different lengths. We compare the numerical results with the predictions by Eq. (18) and Eq. (19), which both include the first and second sources of error described above. We discuss four scenarios labelled (1-4) in Fig. 3(b). Case (1) corresponds to the regime δ g k , for which the only relevant source of error is Γ 0− . Case (2) corresponds to the regime δ ∼ g k , in which a smaller error is obtained in expense of an exchange of population with the chain modes. Case (3) shows that placing the dimer closer to the chain as compared to the lattice spacing d reduces the error thanks to the larger effective coupling rate g π d /Γ 0− [40]. The improvement is remarkable already at small N . Because of the discrete nature of the modes, the prediction in the continuum limit, Eq. (19), leads to an overestimation of the resonance frequency of the modes near the band edge [see Appendix E 2] and, hence, to the disagreement with the case of a finite chain, especially at smaller N . A larger effective coupling between distant dimers can, thus, be obtained by taking δ < 0, while staying off-resonant with the chain modes. For this, we define a new detuning between the dimer and the highest-energy array mode, The condition δ 2 g ij eff can again be made arbitrarily small by tuning Ω/∆. For small δ 2 , however, Γ k also becomes a dominant source of error. Minimizing the error for the simplified model including only the interaction with the highest-energy chain mode, we obtain an optimal δ 2 [see Appendix E 3], for which we predict an error that scales as 1/N , as in case (4) in Fig. 3(b). The error deviates from the prediction at larger N , as the energy spacing between chain modes is reduced with N and, thus, the contribution of further chain modes becomes non-negligible, for which one should go back to use Eq. (18). Note that, although barely captured in the plots, the population in the chain modes becomes nonnegligible for smaller δ/g ij eff , as in cases (2) and (4) with larger N .

V. DISCUSSION AND OUTLOOK
In this work, we have proposed a setup to achieve a coherent and Markovian interaction between an emitter and the subradiant modes of an atomic chain, a system which parallels conventional waveguide QED of atoms coupled to waveguides. Our proposal is based on two main ingredients: (i) the use of ground and antisymmetric dimer states as an effective two-level system and (ii) the use of a Raman transition to control the linewidth of the dimer. The first method exploits the particular symmetry of the dimer state to improve the coherent coupling to the chain's guided modes by decoupling the emitter from the highly radiating modes of the chain. The second method allows to reduce the dimer's linewidth as compared to the chain's bandwidth, thus achieving the regime of Markovian dynamics. Accordingly, we observe similar dynamics as in conventional waveguide QED both in the in-band and bandgap regimes for the case of an atomic chain with interatomic separation of quarter wavelength. Along with the well-known Markovian dynamics, we also observe non-Markovian effects due to the finiteness of the chain and retardation effects introduced by the slow group velocity at the band edge.
In this work, we focused on a simple linear geometry for both the chain and the dimers, which parallels the conventional setup of atoms coupled to a linear waveguide. Several generalisations can also be considered. In particular, for the band-gap regime of emitter-chain coupling, the fidelity for coherent excitation exchange between the dimer and the chain is ultimately limited by the intrinsic decay of the chain's dark modes. This error can be reduced by considering emitters coupled to an atomic ring -the atomic equivalent of a ring-resonator-, where subradiant modes are expected to have an exponential suppression of the decay rate as ∼ Γ 0 exp(−N ) [31,55]. Furthermore, the idea of exploiting the dimer interference to reduce the coupling to the array's bright modes might be extended to 2D and 3D lattices. In higher dimensional lattices, the use of dimers (or the corresponding generalization) is expected to be particularly advantageous due to the scaling with the array's size of the superradiant modes' linewidth [23]. A viable way to improve the chain-mediated coherent coupling between emitters is to consider an emitter-chain separation smaller than the interatomic separation of the chain. A reduced emitterchain separation, which could be achieved by engineering the trapping potential using optical superlattices [56][57][58][59][60], is experimentally more challenging but leads to large improvements in the fidelity as shown in Fig. 3. Finally, while we only considered d = λ 0 /4, our methods apply equally well to smaller interatomic separations, where the system's dynamics would benefit from the additional reduction in the decay to free-space.
This work paves the way toward observing and exploiting the rich phenomenology of waveguide QED in a clean, atom-based, setup. The additional non-Markovian ef-fects due to the finiteness of the array are difficult to observe in standard waveguide QED and are a distinguishing feature of this platform. In this appendix we outline the derivations of a closed expression for the coupling between a single atom emitter and a chain mode k, as in Eq. (10), and of the extension to an atomic dimer to obtain Eq. (12). We, then, go on to discuss alternative setups.
The coupling between an impurity atom at r i and a Bloch mode of the chain in Eq. (8) reads We express Eq. (2) in cylindrical coordinates by using where ρ i > ρ j , and J m and H Plugging Eq. (A3) into Eq. (A1) yields Eq. (10).
For an even N , setting z = 0 at the center of the chain, r j = djẑ with j = −N + 1 2 , −N + 3 2 , . . . , N − 1 2 , and for dimer's atoms in position r ± i = [hρ + (ρ i ± ρ 0 )ẑ], the coupling to the symmetric (λ = 1) and anti-symmetric (λ = −1) state reads Assuming that the dimer is located far from the edges of the chain, we can extend the sum to infinite j without affecting its total value. Shifting the origin to ρ iẑ , the sum above writes where (−)ρ 0 and j stand for r − ρ iẑ and r j , respectively. For ρ 0 = d/2, i.e., the emitters are aligned with the chain atoms as in the text, we obtain Eq. (12).
Repeating the above treatment with the finite-chain Ansatz in Eq. (6) to replace e −ikrj , we derive the following definitions in the final result for λ = −1, and for λ = 1, (A7)

Alternative setup: a 2x2 quadruplet
Let us now consider a 2 × 2 plaquette of atoms and compare it to the case of a dimer considered in the main text. In general, we consider an effective two-level impurity with levels |0 and |Ψ − . For the dimer, we use |0 = |gg , and |Ψ − = (|eg − |ge )/ √ 2 is the antisymmetric state of the two atoms. For the case of a plaquette, |0 = |gggg is the ground state of the four plaquette atoms, and is the anti-symmetric state of the plaquette in both x and z directions, where the general state of the plaquette is given by |ν A ν B ν C ν D , with ν = g, e. In Fig. 4, we compute the evolution of the system initialised in |Ψ − in the presence of an atomic chain for both the dimer (left panel) and plaquette (right panel) configurations. We tune the atomic transition of the single impurity atoms differently for a dimer and a plaquette, such that in both cases the collective state |Ψ − has the same energy with respect to the band of subradiant modes of the array. Fig. 4 shows that, while the 2 × 2 plaquette allows to attain a smaller free-space decay as compared to the dimer, it also leads to a reduction of the coupling with the chain modes and, thus, to a longer time scale for the transfer of population to the chain. After normalizing the time-axis of the two plots by the time scale of the decay rate of the respective impurity state, one can observe no qualitative difference between the evolution of the two systems.

Alternative setup: ρ0 < d/2
For ρ 0 < d/2, we can repeat the steps above to obtain and for the coupling rate of the anti-symmetric eigenstate of the dimer, and replacing the sine with a cosine for the symmetric state. From Eq. (13), we know that γ ab can be further reduced by reducing ρ 0 . Eq. (A9), however, is also reduced in value for ρ 0 < 0.5. A plot of the ratio |g i− k | 2 /Γ 0−a ratio we want to optimize as we see in section III and IV -as a function of ρ 0 shows a slow monotonic increase when moving towards smaller ρ 0 . The values of the ratio at ρ 0 = d/2 and ρ 0 → 0 for d = λ 0 /4 differ by a factor of two.

Appendix B: Effective model for dimer-chain interaction
In this appendix, we derive the effective Hamiltonian for the interaction betweeen the anti-symmetric state of a dimer and the subradiant modes of the chain. According to Eq. (3), the evolution of the dimer and of the chain within the single excitation subspace readṡ where a λ are defined in Eq. (11). Here, we stop using the hat notation for the operators and consider a single dimer for simplicity. The extension to multiple dimers is straightforward.
As we can see in Fig. 1, interference suppresses the coupling between the anti-symmetric state of the dimer and the most dissipative modes of the chain. The same applies between the symmetric state and the subradiant modes. For a dimer anti-symmetric state energy, E − , resonant with the subradiant array modes, J k>k0 , and d = λ 0 /4, the symmetric state energy, E + , is resonant with the superradiant modes, J k<k0 . In such case, We set the zero of energy at E − and define ∆ + = (E + − E − ) and ∆ k = (J k − E − ). Under these conditions and assuming that the system is initialized in the dimer anti-symmetric state, a − , we adiabatically eliminate the symmetric dimer state and the superradiant modes of the chain by assuming that they acquire a negligible population,ȧ + =ḃ k<k0 = 0, solve the corresponding equations in Eq. (B1), and substitute the resulting expressions for a + and b k<k0 in the expressions forȧ − andḃ k>k0 . We obtain effective equations of motion of the form of Eq. (B1) including only the anti-symmetric dimer state and the subradiant modes of the chain by doing the substitutions with whereg λ * k = g λ * k − i 2 γ λ * k = ξ λ * k |g λ k | − i 2 |γ λ k | . These corrections are negligible at the edge of the chain's band, as verified numerically by comparing the evolution under the full Hamiltonian and the effective Hamiltonian derived in this appendix. This is true because the coupling rate of the dimer anti-symmetric state with the most superradiant chain modes around k = 0, with which the dimer is resonant due to their large linewidth, is close to zero. Since the edge of the band is our region of interest, we can hence model the system with a two-level impurity consisting of the ground and anti-symmetric state of the dimer interacting coherently with the guided (k > k 0 ) modes of the array. For lower impurity energies closer to J k0 , the dimer becomes resonant with additional, less broad, superradiant modes with which the coupling is non-zero [see Fig. 2]. In this case, Eq. (B2) introduces finite corrections to the effective dynamics. This additional coupling to superradiant channels is further suppressed after we introduce a Raman transition in the dimer atoms due to the reduction of the dimer states' linewidth, as we see in Appendix C.

Appendix C: Effective Hamiltonian with a Raman transition
In this appendix, we extend the effective description developed in Appendix B to the case of impurity atoms driven on a Raman transition (see section II B). We start from Eq. (14) and separate it into a bare,Ĥ 0 , and an interacting,V , part aŝ and for impurity atoms at position (h, 0, (C3) where we use the definition ofg * k in Appendix B.
We shift the energy by ∆ − ω imp 0 , such that the excited states |e α evolve fast, and define the projectorŝ We then calculatePV (QĤ 0Q ) −1VP and use that, for ∆, Γ 0 Ω, g ab , γ ab , g i k , γ i k , and to second order in perturbation [61,62], The resulting effective Hamiltonian in the dimer eigenstate basis and after undoing the previous energy shift reads, In the regime ∆ Γ 0 Ω, g ab , γ ab , g i k , γ i k , the dimer energy and the dimer-chain coupling can be approximated as in Eq. (15). For the dimer energy resonant with the subradiant region of the chain's band, we can eliminate the symmetric dimer state and the superradiant chain modes in the same way than in Appendix B. Note that with the Raman transition, the symmetric state is also resonant with the subradiant modes. Nevertheless, its coupling is smaller and its linewidth larger than the one of the anti-symmetric state, and we verify that it can be safely neglected when studying the dynamics of the latter. In contrast to Appendix B, the reduced dimer linewidth due to the Raman transition maintains the correction due to the superradiant chain modes small at energies close to J k0 . The additional ∼ g 2 k term in Eq. (C7) introduces small shifts in the energy of the chain modes. We account for them by correcting the k axes in Fig. 2 with the momentum extracted from the chain mode be-ing excited, which might not coincide exactly with the k expected from J k . In Fig. 5, we compare the evolution of a dimer anti-symmetric state resonant with the band of subradiant chain modes computed both with the full Hamiltonian in Eq. (14), and the effective Hamiltonian in Eq. (15), and verify that the approximations discussed here are reasonable. In this appendix, we show the non-Markovian dynamics of a dimer relaxing into an atomic chain when the dimer is resonant with a symmetric array modes. At long times, comparable to one over the energy difference between the discrete guided modes of the array, we observe non-Markovian effects in the dynamics of a dimer with energy inside the chain's band. We understand these effects as the retarded back-action via the electric field of the emitter, which is reflected back at the ends of the chain. For a dimer aligned with the center of the chain, and with its anti-symmetric state energy resonant with a chain mode with even symmetry, such as k · d/π = 0.922 in Fig. 6, the coupling rate to that mode is zero. Therefore, the dimer's energy effectively lies between k · d/π = 0.920 and k · d/π = 0.924. The reflected wave is, thus, out of phase with the dimer, as discussed in the text, and leads to an increase of the dimer's population, as we see when the evolution enters the shaded area in Fig. 6. Appendix E: Effective dimer-dimer long range interaction in the band-gap In this appendix, we provide additional detail to the case of a multiple dimers interacting with an atomic chain when the dimers' energy lies within the chain's bandgap. In particular, we derive the effective model for the dimers-chain evolution in this regime, we explain the origin of the discrepancy between the case of finite and infinite long chain, and we derive the optimal error's scaling for the excitation swapping between dimers mediated by the chain.
The evolution of the operators a i− and b k under the Hamiltonian in Eq. (15) iṡ Moving to a rotating frame,ã i− (t) = a i− (t) e iω imp 0 t and b k (t) = b k (t) e iJ k t , solving forb k (t), and substituting into the equation forȧ i− (t), we obtain an effective equation for the dimer's population, We assume that, after summing over momenta, the time integral only contributes for a small correlation time τ c (Markov approximation). Since g k is approximately constant close to the band edge, the region that is closest in resonance with the dimers, τ c is short. Assuming that the dimer operator evolves over time scales much longer than τ c , we approximateã † i− (t ) ã † i− (t). For t τ c and ω at − J π/d > 0, the equation for the evolution of the dimer population readṡ from which we can infer the effective dimer-dimer interaction in Eq. (18). The biggest contributions to g ij eff come from the k near the band edge. To obtain a closed form of the effective coupling between dimers through the guided modes of the chain, we approximate g k g π d , and do a quadratic band approximation of J k around the band edge. Expanding the analytical form of the dispersion relation for an infinite chain [31] around k = π/d and truncating to second order, we obtain J π , such that the discrete k µ = π d 1 − µ N +1 , the decay rate of the most subradiant modes scales as Γ µ γ N µ 2 /N 2 [31]. Thus, Γ π d (1−x) γ N x 2 , and in the continuous band approximation, where δ = ω imp 0 − J π d . For δ A d , the integration limits can be extended to infinity without affecting the solution. We introduce η by substituting (1 − x) with (η − x). Using the convolution theorem (f * g)(η) = F −1 {F(f ) · F(g)}, and the results of the Fourier transform F(e ikρ x )(ν) = δ(kρ − ν), and and Eq. (19) and Eq. (20) follow trivially. Note that the exponential envelope in real space introduces a length scale for the interactions: l = A d /δ.

Error due to populating the chain modes
In the regime in which δ g π d is not satisfied, but = g π d /δ is still small, the transition probability at time, t, from an initial dimer state with energy J π d + δ to the chain modes, to first order [63] (E7) In the continuum limit, Σ k → N d 2π dk, we extract an expression for the absolute maximum of population at the chain Note that this probability is an upper bound both because it is the maximum of population at the chain, and because the integral in the continuum limit overestimates the value of the sum over k for a small δ compared to J π d − J π d (1−1/N ) , as discussed in section IV and Appendix E 2. Since now we compute 1/δ 2 , the mismatch between the results from using a continuum or a discrete band scales faster than in the calculation of g ij eff . Finally, the error defined in the text is not well suited to capture error due to populating the chain modes, as the time scale of the oscillations between chain and dimers is much shorter than the one of the oscillations between dimers. The error at the dimer's maximum is, therefore, most likely measured at a moment in which the population in the chain is zero, and the probability that such point of time exists close to the maximum of dimer population is high.
2. Prediction mismatch due to the finite N In Fig. 3, we observe a difference between the effective dimer-dimer interaction as described by Eq. (18) and Eq. (19) for a dimer' energy in the band-gap, as depicted in Fig. 7 (left). This disagreement is due to the discreteness of the chain modes, for which the highest-energy mode has quasi-momentum k µ=1 [see Appendix E], with k µ=1 π d (1 − 1/N ). Since the biggest contributions to g ij eff come from the modes closest in energy to the band edge, and J kµ=1 < J π d , the larger detuning between the dimer and the mode with k µ=1 (δ 2 in Fig. 7) explains the smaller coupling rates predicted by Eq. (18). In other words, the approximation in Eq. (E5) may not be appropriate below a certain N .
We can better understand the mismatch between the discrete and infinite chain predictions by looking at the integrand in Eq. (E5) as we do in Fig. 7 (right). For simplicity, we focus the comparison on the most dominant term of the discrete sum, the one with k µ=1 , for which the difference in total value between the two sides of Eq. (E5) is captured by the non-shaded areas in Fig. 7. For instance, the ratio between the non-shaded part of the area under f (x) and the full integral is rather small for = 2 × 10 −3 and N = 500, and the offset between the corresponding marker and dashed line in Fig. 3 is also small. However, the steepness of the function at small x can lead to large underestimations of the area under the function. This becomes critical at smaller N . For = 2 × 10 −3 and N = 100, the non-shaded part amounts to multiple times the shaded part of the area under f (x). The value of the mismatch between the discrete and continuum predictions also depends on δ, as smaller δ are better able to resolve the detuning of the mode k µ=1 with the band edge. In other words, f (x) becomes steeper, and the amount of underestimation of the integral using the discrete sum is larger, as shown in Fig. 7 for = 10 −1 . This, again, can be verified by observing Fig. 3. The predictions with Eq. (19) are, thus, in general overly optimistic. However, by studying the discrete spectrum, we can optimize the predicted coupling rates, as we do in the following section.
3. Optimal error for δ < 0 Analyzing the system from the picture of a discrete spectrum of array modes with momenta k µ , as defined above and sketched in Fig. 7 (left), instead of a continuous band, we can improve the long-range coupling rate between dimers. In this case, the expression for the effective dimer-dimer coupling reads Defining a detuning δ 2 from the highest energy mode of the array (µ = 1), .

(E11)
To stay off-resonant with the mode k µ=1 , we need Thus, we approximate Eq. (E11) with the first term of the sum, where for simplicity we have also used ξ − k1 (ρ i )ξ − * k1 (ρ j ) 1, valid for a long chain and dimers located near its center.
As discussed in section IV, the dimer-dimer long-range coupling has three main sources of error. One is due to Γ 0− , given by the decay of the dimer population after one full Rabi cycle, π Γ 0− /Re[g ij eff ]. A second source is due to the small but finite linewidth of the array modes that mediate the interaction, Γ k , and write 2π Im[g ij eff ]/Re[g ij eff ]. Thirdly, the finite population transfer to the array mode if δ 2 g π d is not satisfied. The maximum of population in the chain is 4g 2 π d /(δ 2 2 + 4g 2 π d ). The error in Fig. 3 is due to the free-space decay of the dimers and that of the array modes and reads By minimizing the error, we obtain an optimal value of δ 2 for which we predict an error Since g − π d scales as 1/ √ N , the optimal error scales as 1/N .