Pronounced non-Markovian features in multiply-excited, multiple-emitter waveguide-QED: Retardation-induced anomalous population trappin
Alexander Carmele, Nikolett Nemet, Victor Canela, Scott Parkins
PPronounced non-Markovian features in multiply-excited, multiple-emitterwaveguide-QED: Retardation-induced anomalous population trapping
Alexander Carmele, Nikolett Nemet,
1, 2
Victor Canela,
1, 2 and Scott Parkins
1, 2 Department of Physics, University of Auckland, Private Bag 92019, Auckland, New Zealand Dodd-Walls Centre for Photonic and Quantum Technologies, New Zealand (Dated: February 5, 2020)The Markovian approximation is widely applied in the field of quantum optics due to the weakfrequency dependence of the vacuum field amplitude, and in consequence non-Markovian effects aretypically regarded to play a minor role in the optical electron-photon interaction. Here, we givean example where non-Markovianity changes the qualitative behavior of a quantum optical system,rendering the Markovian approximation quantitatively and qualitatively insufficient. Namely, westudy a multiple-emitter, multiple-excitation waveguide quantum-electrodynamic (waveguide-QED)system and include propagation time delay. In particular, we demonstrate anomalous populationtrapping as a result of the retardation in the excitation exchange between the waveguide and threeinitially excited emitters. Allowing for local phases in the emitter-waveguide coupling, this pop-ulation trapping cannot be recovered using a Markovian treatment, proving the essential role ofnon-Markovian dynamics in the process. Furthermore, this time-delayed excitation exchange allowsfor a novel steady state, in which one emitter decays entirely to its ground state while the other tworemain partially excited.
I. INTRODUCTION
One-dimensional (1D) waveguide-QED systems are at-tractive platforms for engineering light-matter interac-tions and studying collective behavior in the ongoingefforts to construct scalable quantum networks [1–12].Such systems are realized in photonic-like systems includ-ing photonic crystal waveguides [13–19], optical fibers[20–24], or metal and graphene plasmonic waveguides[25–28]. Due to their one-dimensional structure, long-distance interactions become significant [3, 5, 29]. As aresult of these interactions mediated by left- and right-moving quantized electromagnetic fields, strongly entan-gled dynamics and collective, cooperative effects relatedto Dicke sub- and superradiance emerge [1, 6, 12, 17, 22–24, 30–37].In the framework of standard quantum optics, thesesystems are widely explored in the Markovian, single-emitter or single-excitation limit [9, 31, 38–41]. Such lim-its can be described by a variety of theoretical methodsincluding real-space approaches [5, 38, 42, 43], a Green’sfunction approach [44–47], Lindblad master equations[48, 49], input-output theory [50–54], and the Lippmann-Schwinger equation [55–57]. Already in these regimes,exciting features have been predicted. For example,strong photon-photon interactions can in principle beengineered, allowing for quantum computation protocolsusing flying qubits (propagating photons) and multilevelatoms [5, 11, 56, 58, 59]. Furthermore, bound statesin the continuum are addressed via a joint two-photonpulse, showing that excitation trapping via multiple-photon scattering can occur without band-edge effectsor cavities [7, 43, 55, 60].Beyond the single-excitation and/or single-emitterlimit, the Markovian approximation becomes question-able and the aforementioned methods problematic [61– 68]. In this work, we employ the matrix-product staterepresentation to study exactly this regime, the multiple-excitation and multiple-emitter limit. We focus, inparticular, on the three-emitter and three-photon case,treating the emitters as two-level systems, which coupleto the left- and right-moving photons and thereby inter-act with each other, subject to time delays associatedwith the propagation time of photons between emitters[58, 69, 70]. We choose throughout the paper the triply-excited state as the initial state and compare the re-laxation dynamics in the Markovian and non-Markovian (a) (b)
FIG. 1. Scheme of the simulated waveguide QED system. (a)The system consists of three identical emitter with transitionfrequency ω which couple to left- and right moving quantizedlight fields via the decay constant √ γ . (b) Due to the delay(in the scheme two time steps τ = τ = 2∆ = 2 γ − / τ = 4∆) a closed loop is formed between the first andthird emitter interacting with their respective past bins. Theinteraction strongly depends on the phases ω τ and ω τ . a r X i v : . [ qu a n t - ph ] F e b cases. To compare both scenarios on the same footing,we employ the quantum stochastic Schr¨odinger equationapproach [61, 71, 72] and numerically solve the modelusing a matrix-product-state algorithm [59, 73–77] as analternative to the t-DMRG method in position space [78].We report on striking differences between the Markovianand non-Markovian description. First, we find that in thecase of non-Markovian excitation exchange, the triply-excited initial state allows for population trapping, instrong contrast to the Markovian description. Second,time-delayed excitation exchange allows for anomalouspopulation trapping, in which one emitter relaxes com-pletely into its ground state while the two other emit-ters form a singly-excited dark state together with thewaveguide field in between. No local phase combinationin the Markovian case allows for such anomalous popula-tion trapping, rendering the non-Markovian descriptionqualitatively and quantitatively different from a Marko-vian treatment. II. MODEL
To demonstrate the importance of retardation-inducedeffects and the underlying non-Markovian dynamics, wechoose a system consisting of three identical emitterswith transition frequency ω . All three emitters interactwith left- ( l ( † ) ω ) and right-moving photons ( r ( † ) ω ) in a one-dimensional waveguide, as depicted in Fig. 1(a). To fo-cus on the retardation-induced effects, we neglect out-of-plane losses which inevitably enforce a fully thermalized,trivial steady-state in the ground state, and render thenon-Markovian effects a transient, nevertheless impor-tant feature for waveguide-based counting experiments.The Hamiltonian governing the free evolution of the com-bined, one-dimensional waveguide photon-emitter systemreads: H / (cid:126) = ω (cid:88) i =1 σ i + (cid:90) dω ω (cid:0) r † ω r ω + l † ω l ω (cid:1) , (1)where the emitters are treated as two-level systems, with | (cid:105) as the ground state and | (cid:105) as the excited state, andwith σ ijn := | i (cid:105) nn (cid:104) j | , the flip operator of the n-th emitter.The interaction Hamiltonian describes the emitters inter-acting with right and left moving photons at the emitters’positions: H I = (cid:126) g (cid:88) i =1 σ i (cid:90) dω (cid:16) r † ω e iωx i /c + l † ω e − iωx i /c (cid:17) + H.c. , where we have assumed a frequency-independent cou-pling of the emitters to the quantized light field. The po-sition of the second emitter is chosen as x = 0, leading to x = − d / − cτ / x = d / cτ / c the speed of light in thewaveguide. After transforming into the interaction pic-ture with respect to the free evolution Hamiltonian, and applying a time-independent phase shift to the left- andright-moving photonic field, the transformed Hamilto-nian reads H I ( t ) = (cid:126) g (cid:82) ( H r,ωI + H l,ωI ) dω , where H r,ωI = r † ω ( t ) (cid:0) σ + σ e − i ω τ + σ e − i ω τ (cid:1) + H.c. ,H l,ωI = l † ω ( t ) (cid:0) σ e − i ω τ + σ e − i ω τ + σ (cid:1) + H.c. , (2)with r † ω ( t ) = r † ω (0) exp[ i ( ω − ω ) t ] and l † ω ( t ) = l † ω (0) exp[ i ( ω − ω ) t ], and τ = ( τ + τ ) /
2, cf. App. A.In the following, the left- and right-moving excitationsare treated collectively: R † ( t ) = (cid:90) dωr † ω ( t ) , L † ( t ) = (cid:90) dωl † ω ( t ) . (3)Given these definitions, the non-Markovian interactionHamiltonian reads: H NM I ( t ) / (cid:126) = g (cid:0) σ (cid:0) R † ( t ) + e iω τ L † ( t − τ ) (cid:1) + H.c. (cid:1) + g (cid:16) σ R † ( t − τ / e i ω τ + H.c. (cid:17) (4)+ g (cid:16) σ L † ( t − τ / e i ω τ + H.c. (cid:17) + g (cid:0) σ (cid:0) R † ( t − τ ) e iω τ + L † ( t ) (cid:1) + H.c. (cid:1) . In the following, we compare the Markovian with thenon-Markovian case. The Markovian case neglects retar-dation effects between the excitation exchange, thereforein the Markovian approximation we set R ( † ) ( t − t (cid:48) ) ≈ R ( † ) ( t ) and L ( † ) ( t − t (cid:48) ) ≈ L ( † ) ( t ). In this approxima-tion, only the local phases but not the retardation inthe amplitude are taken into account. Consequently, theMarkovian interaction Hamiltonian reads: H M I ( t ) / (cid:126) = g (cid:104) L † ( t ) (cid:16) σ e iω τ + σ e i ω τ + σ (cid:17) + H.c. (cid:105) + g (cid:104) R † ( t ) (cid:16) σ + σ e i ω τ + e iω τ σ (cid:17) + H.c. (cid:105) , (5)where the emitters interact with time-local collectiveright- and left-moving fields and no time delay is presentin the interaction. We solve for the system’s dynamics inboth cases using the time-discrete Schr¨odinger equationwith the time-step size ∆ up to time N ∆ in N steps,cf. Fig. 1(b): | ψ ( n ) (cid:105) = U NM/M ( n, n − | ψ ( n − (cid:105) (6)= exp (cid:34) − i (cid:126) (cid:90) n ∆( n − H NM/M I ( t (cid:48) ) dt (cid:48) (cid:35) | ψ ( n − (cid:105) , where ∆ is small enough to minimize the error in theSuzuki-Trotter expansion [59, 73–77], and the evolutionis taken either in the Markovian (M) or in the non-Markovian limit (NM). Here, the wavefunction is inMPS form: | ψ ( n ) (cid:105) = (cid:88) s,l ··· l N r ··· r N L [ l ] R [ r ] · · · S [ s ] L [ l n ] R [ r n ] [ l n +1 ] [ r n +1 ] . . . | l , r · · · s, l n , r n , l n +1 , r n +1 · · · l N , r N (cid:105) , (7) FIG. 2. Algorithm to compute the next time step via MPSin case of a Markovian dynamics, i.e. we set R ( † ) ( n − m ) ≈ R ( † ) ( n ) and L ( † ) ( n − m ) ≈ L ( † ) ( n ) and m ∆ = τ . where we assume vacuum input states for m > n withidentity tensor L [ l m ] = [ l m ] and R [ r m ] = [ r m ] in the timebins for the right- and left-moving field corresponding toa vacuum input state from left and right, and the indices s count the degrees of freedom for the emitter system,and r i , l i the number of excitations in the time bin for theleft- and right-moving field. In the following, we assumedim[ s ] = 2 = 8, and choose a time-step size to guaranteethat the dimension of the reservoir excitation does notexceed dim[ r i ] = dim[ l i ] = 3. However, all results havealso been calculated with dim[ r i ] = dim[ l i ] = 4 to proveconvergence.To efficiently simulate the multiple-emitter, multiple-excitation case, we employ the matrix-product-statetechnique described in [59, 73–77], and choose a collec-tive basis for the flip operators of the emitters to al-low for entangled initial states: | ijk (cid:105) = | ( i − +( j − + ( k − (cid:105) , which leads to, e.g., σ ≡| (cid:105)(cid:104) | + | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + | (cid:105)(cid:104) | , or | (cid:105) = | (cid:105) the triply-excited state. We discuss in the corresponding sectionsthe simulation protocol in detail. III. MARKOVIAN LIMIT: NO TIME DELAY
We start our investigation in the Markovian limit,and calculate the system’s dynamics with U M (Eq.(6)) and the initial state | ψ (0) (cid:105) = | (cid:105) until thesteady state is reached. The Markovian case allowsfor a master-equation treatment with g = √ πγ [61–68, 79]. Tracing out the left- and right-movingphotons leads to a collective jump operator, J := √ γ (cid:0) σ exp[ iϕ ] + σ + σ exp[ iϕ ] (cid:1) . The phases ϕ i can be chosen individually via local unitary transforma-tions, or they arise from the spatial position withouttaking the finite distance into account in the evolution[6, 70]. In the following, we nevertheless solve the dy- FIG. 3. The impact of different phase choices ϕ :=( ϕ , ϕ , ϕ ) = ([0 , π ] , , [0 , π ]) in the atom-waveguide cou-plings in the Markovian limit with photon operators: R ( t − t (cid:48) ) = R ( t ) and L ( t − t (cid:48) ) = L ( t ), on the integrated reservoirpopulation in the steady state. If the system is initialized inthe triply-excited state (green line), for all choices of phases,all excitation is radiated into the reservoir. If the system isinitialized in a superposition of doubly-(orange line) or singly-excited states (red line), the only case where all excitation isradiated into the reservoir is when all phases are a multipleof 2 π . namics using the quantum stochastic Schr¨odinger equa-tion ignoring time-delay effects to give a Markovian evo-lution.In Fig. 2 the simulation protocol is depicted. In thetime-discrete basis, the time-local Hamiltonian in thematrix-product-operator MPO(n) form acts only on thepresent reservoir bins n and the system state s due to theMarkovian approximation. After applying the evolutionoperator via contracting the physical indices (in the cor-responding color code, red for the system bin, green forright moving field, blue for left moving field) in step (a),an entangled system-reservoir matrix is created. To writethe MPS in the canonical form, a Schmidt value decom-position is performed, and the entangled system-reservoirstate is expressed as a matrix product (c), after whichthe next MPO(n+1) can be applied. In this manner, thesteady state can be calculated step by step in an efficientand less memory-consuming way. Due to the Marko-vian approximation, i.e. we set R ( † ) ( n − m ) ≈ R ( † ) ( n )and L ( † ) ( n − m ) ≈ L ( † ) ( n ) and m ∆ = τ , no rear-ranging of the MPS is necessary. The simulations aredone, in the Markovian and non-Markovian case, untilthe steady state is reached. In the MPS representation,(Eq. (7)), the steady state is reached when the applica-tion of the MPO leads to the identical tensor combinationfor every time step after the Schmidt value decomposi-tion and swapping procedure, i.e. the MPO and subse-quent re-arranging only acts as an index-shifting operator l n , r n → l n +1 , r n +1 .In Fig. 3, the phase dependence of the integrated reser- FIG. 4. Algorithm to compute the next time step via MPS in case of retardation of the left- (blue) and right-moving quantizedfield (green). The emitter 2 (middle) interacts with a time delay of 2∆ with the right-moving field, taking emitter 1 as point ofreference for the right-moving field, and emitter 2 also interacts with the emitted left-moving field of emitter 3, taking emitter3 as reference point for the left-moving field, resulting in a delay of 4∆. The dashed squares couple to the system s at time n ∆.Since the canonical form of the MPS needs to be maintained, swapping procedures need to be applied where the orthogonalitycenter is swapped from n − voir excitation I in the steady state is plotted for differ-ent initial states: I = (cid:80) N − fn =0 (cid:10) R † ( n ) R ( n ) (cid:11) + (cid:10) L † ( n ) L ( n ) (cid:11) with N as the number of time steps to reach the steadystate and f = ( τ + τ ) / ∆. The phases are permu-tated by changing ϕ , ϕ between 0 to 2 π . If the systemis initialized in the triply-excited state | ψ (0) (cid:105) = | (cid:105) , thephases have no impact at all on the steady-state valuesand all excitation will eventually be radiated into thereservoirs on the left of emitter one and on the right ofemitter three, leading for all phase permutations to theintegrated reservoir occupation of 3, cf. Fig. 3 (greenline).In contrast to the triply-excited case, the steady statesof the emitters initially in a superposition of singly-( | (cid:105) + | (cid:105) + | (cid:105) , red line) and doubly-excited states( | (cid:105) + | (cid:105) + | (cid:105) , orange line) are strongly influenced bythe choice of phases. For those initial states, only if thephase difference vanishes, ϕ = 2 π = ϕ , is all radia-tion emitted into the reservoir. For all other phase com-binations, population trapping occurs, and all emittershave a finite probability to be found in the excited state.Population trapping is created due to the fact that the in-dividual decay of the emitters allows to populate a darkstate of the Hamiltonian. In the two-emitter case, thestandard light-matter Hamiltonian can be written in thecollective basis as: H sl − m = (cid:126) g (cid:88) i =1 , ( σ i + σ i ) (8) ≡ (cid:126) g [ | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + h.c.]and therefore with | D (cid:105) = ( | (cid:105) − | (cid:105) ) / √ H sl − m | D (cid:105) = 0. Therefore, with cor-responding phase differences, such dark states can bedriven via individual decays and lead to dark state pop-ulation, or population trapping, e.g. [3, 6, 23, 24, 44,73, 80, 81]. We conclude that, within the Markoviantreatment, we find that either all emitters relax intotheir ground state or none. For systems initialized inthe triply-excited state, excitation trapping cannot beachieved. And for the superradiant, symmetric singly-excited and doubly-excited initial states, the emitters un-dergo complete decay only in the case of vanishing phasedifference. We show now that including retardation andback-action effects changes this picture completely. IV. NON-MARKOVIAN DYNAMICS:SYMMETRIC TIME DELAY.
In the non-Markovian case, the MPO is not only act-ing on the present reservoir bins n but also on the pastleft- and right-moving reservoir bins n − m for m >
0. InFig. 4, the simulation protocol is schematically explainedfor the case when emitter 2 (middle) interacts with a timedelay of τ / τ / n of theright-moving (green) and n − τ = 6∆. In Fig. 4, thedotted squares are those with which the MPO at timestep n interacts, and the orthogonality center is initiallyin n − P opu l a ti on Emitter 1Emitter 2Emitter 3 t/ τ γτ=0.5 γτ=1.25 γτ=6.25 FIG. 5. The dynamics of the emitter populations for differentfeedback lengths with phase ω τ / π for a system initiallyin the triply-excited state. Already short feedback times ( τ =2ns, i.e. γτ = 0 .
5) lead to population trapping in contrastto the Markovian case. For longer feedback ( τ = 25ns, i.e. γτ = 6 . the dotted squares need to be arranged in the MPS nextto each other to avoid a memory-consuming contractionof the MPS without gain of information. To ascertainthe normalization, the swapping procedures start fromleft to right, cf. Fig. 4(a-b). The swapping guaranteesthat the essential entanglement in between the reservoirbins is preserved. Every swap creates a new orthogonal-ity center, but in (c) the n − s and creates a well-defined orthogonality center before theMPO is applied. After applying the MPO, the resevoirbins are arranged in the previous, canonical form but theorthogonality center is left at n − τ = τ . In thecase of quantum coherent feedback, the delay between theexcitation exchanges introduces a corresponding phase[59, 73–75, 82–87]. In the following, we assume a tran-sition frequency of the emitters to yield: ω τ / π ,i.e. exp[ ± iω τ ] = 1 = exp[ ± iω τ / (cid:10) σ (cid:11) = (cid:80) i =1 |(cid:104) i − | ψ ( t ) (cid:105)| . Even for short feedback in compari-son to the decay time, i.e., γτ = 0 .
5, population trappingis observed (upper panel). Emitter one (black dotted)and three (orange solid line) exhibit the same dynamicsdue to the symmetry of the system. Both start to de- t/ τ P opu l a ti on Emitter 1Emitter 2Emitter 3
FIG. 6. The dynamics of the emitter populations if all emit-ters are initially in their excited state: | ψ (0) (cid:105) = | (cid:105) for a phasechoice of ω ( τ + τ ) = 3 π , and τ = 2 τ . Emitter 3 (greenline) decays completely into its ground state while emitter 1(black line) and 2 (orange line) remain in a partially excitedstate. This steady state is impossible to reach in the Marko-vian treatment if no additional interactions are included. cay exponentially as expected before the first excitationwith a neighboring emitter takes place t ∈ [0 , τ / γτ = 6 .
25, aregime where the feedback phase ceases to have a stronginfluence, we observe an interesting oscillatory behaviorin the emitter populations due to the feedback and finiteexcitation, which settles eventually to a small but finitesteady-state value.This example shows that allowing even for a shortretardation and back-action time, the dynamics of theemitter populations changes qualitatively and quantita-tively. Population trapping from an initial triply-excitedemitter state cannot be recovered just with local phasesin the Hamiltonian. This impossibility is lifted due toa time-delayed coherent feedback mechanism. We em-phasize that for very long delay the need for a particularchoice of ω τ / π is partially lifted, and it takes di-vergingly long for the emitter population to decay. How-ever, for long delay times γτ (cid:29)
1, the population is alsotrapped in the reservoir between the emitters and theabsolute stored population in the emitter system is ex-ponentially small.
V. NON-MARKOVIAN DYNAMICS:ASYMMETRIC TIME DELAY.
Until now, we have discussed symmetric time delay τ = τ , or | x | = | x | . Asymmetric time delay pro-vides a further example how the full non-Markovian andquantized description of many-excitation dynamics inwaveguide-QED deviates qualitatively from the Marko-vian treatment. As shown above, in the Markovian treat-ment either all emitters remain partially excited or noneof them do. For the triply-excited state, no populationtrapping occurs, and for other initial states the systemis not able to reach a steady state with only one emitterin the ground state and the other emitters partially ex-cited. We show now that the non-Markovian descriptionwith asymmetric time delay allows for another exampleof anomalous population trapping, where one emitter de-cays completely into its ground state whereas the otheremitters have a finite probability to be found in the ex-cited state.In Fig. 6, we choose a phase ω τ = 3 π and delay times γτ = 1 between the left (1) and middle emitter (2),and γτ = 0 . t < τ , all emitters radiate unperturbed intothe reservoir. For τ < t < τ , the left emitter (blackline) continues to radiate unperturbed whereas the mid-dle and right emitters start to interact with the emit-ted photons. Due to symmetry, both the right and mid-dle emitters exhibit the same decay behavior for t < τ .This picture changes for larger times, as now the middleemitter’s field starts to constructively interfere with theright-moving photons from emitter one. Emitter threeinteracts with its own past emission and decays faster,while emitters one and two start to form a superpositionstate. After several roundtrip times, γt (cid:38)
15, emitterthree has decayed, and no emission takes place.Interestingly, a necessary condition for this feature tohappen is asymmetric feedback. A symmetric feedback τ = τ exhibits, as in the Markovian case, only fi-nite population in all emitters, or none. This effect de-pends only on the destructive and constructive interfer-ence between left- and right-moving photons. For differ-ent ϕ = ω τ , a different positioning needs to be chosen.Quantity γτ determines the extent of population trap-ping between emitter one and two, but not the qualita-tive effect. VI. CONCLUSION
We have investigated a waveguide-QED system con-sisting of three emitters initialized in the triply-excited state, which interact via left- and right-moving photons.We compared the Markovian and the non-Markoviancase, i.e., without and with time delay in propagationbetween them. In the Markovian case, only a local phaseis taken into account but no delayed amplitude in there-emission events. We recovered the well-known results,that the triply-excited state decays, independent ofphase choice, while the doubly- and singly-superradiantsuperposition state shows population trapping for anynon-vanishing phase differences. In strong contrast, anon-Markovian excitation exchange results in populationtrapping even if the system is initialized in the triply-excited state. Furthermore, quantum feedback allowsfor states in which two emitters form a superpositionstate together with a part of the reservoir, whereas thethird emitter relaxes entirely into the ground state;a state that is not possible to realize in a Markoviansetup if only local phases in the jump operators, and noadditional interactions, are assumed. These examplesprove the significance of time delay in many-emitter,many-excitation systems and the possibility of entirelynew physics beyond the Markovian regime in thesteady-state and 1D (or β →
1) limit considered here,whereby the emitters radiate purely into the (detectable)waveguide modes, which is a regime already in reach ofvarious waveguide-QED platforms [88, 89].AC gratefully acknowledges support from the DeutscheForschungsgemeinschaft (DFG) through the projectB1 of the SFB 910, and from the European UnionsHorizon 2020 research and innovation program under theSONAR grant agreement no. [734690]. The calculationswere performed using the ITensor Library [90].
Appendix A: Hamiltonian of the multiple-emitterwaveguide-QED system
The free evolution Hamiltonian of the combined one-dimensional waveguide photons and emitters systemreads: H / (cid:126) = (cid:88) i =1 , , ω i σ i + (cid:90) dω ω (cid:0) r † ω r ω + l † ω l ω (cid:1) , (A1)where the emitters are treated as two-level systems with | (cid:105) as the ground state and | (cid:105) as the excited state, and σ ijn = | i (cid:105) nn (cid:104) j | for the n-th emitter, with σ n = | (cid:105) nn (cid:104) | the de-excitation operator of the n-th emitter. The in-teraction Hamiltonian consists of the emitter interactingwith right- and left-moving photons at the emitter’s po-sition: H I / (cid:126) = (cid:88) i =1 , , σ i (cid:90) dω g i ( ω ) (cid:16) r † ω e iωx i /c + l † ω e − iωx i /c (cid:17) + (cid:88) i =1 , , σ i (cid:90) dω g ∗ i ( ω ) (cid:16) r ω e − iωx i /c + l ω e iωx i /c (cid:17) . (A2)The positions of the atoms are assumed to be centeredaround x = 0, so in the case of three atoms, the mid-dle atom is located at x = 0, the first (from left)atom at x = − d / − cτ /
2, and the third atomat x = d / cτ /
2. Choosing a rotating frame cor-responding to the middle atoms frequency ω and thereservoir modes of the left and right-moving photons, thetotal Hamiltonian reads: H I / (cid:126) = δ σ + δ σ (A3)+ (cid:88) i =1 , , (cid:18) σ i (cid:90) dω g ∗ i ( ω ) r † ω e iωτ i / e − i ( ω − ω ) t + H.c. (cid:19) + (cid:88) i =1 , , (cid:18) σ i (cid:90) dω g ∗ i ( ω ) l † ω e − iωτ i / e − i ( ω − ω ) t + H.c. (cid:19) , with δ = ω − ω and δ = ω − ω . For convenience,we transform this Hamiltonian. The left atom interactswithout delay with the right-moving field. Also, we wantthe right atom to interact with the left-moving field with-out delay: To achieve this, we apply unitary transforma-tions: U = exp (cid:20) − iτ (cid:90) dωωr † ω r ω (cid:21) → U r r † ω U † r = r † ω e − iωτ / U = exp (cid:20) − iτ (cid:90) dωωl † ω l ω (cid:21) → U l l † ω U † l = l † ω e − iωτ / . Now, the interaction Hamiltonian reads: H I ( t ) / (cid:126) = (A4)= δ σ + δ σ + (cid:90) dω (cid:18) σ g ∗ ( ω ) (cid:16) r † ω + l † ω e − iω ( τ + τ ) / (cid:17) e − i ( ω − ω ) t (A5)+ σ g ∗ ( ω ) (cid:16) r † ω e − iωτ / + l † ω e − iωτ / (cid:17) e − i ( ω − ω ) t + σ g ∗ ( ω ) (cid:16) r † ω e − iω ( τ + τ ) / + l † ω (cid:17) e − i ( ω − ω ) t + H.c. (cid:19) . In the following, the left and right-moving excitations aretreated collectively, and we assume the coupling elementsto be constants with respect to frequency, i.e., g i ( ω ) = (cid:112) γ i / (2 π ). New operators are introduced: R † ( t ) = (cid:90) dωr † ω e i ( ω − ω ) t / √ π, (A6) L † ( t ) = (cid:90) dωl † ω e i ( ω − ω ) t / √ π. (A7)Given these definitions, the interaction Hamiltonianreads: H I ( t ) / (cid:126) = δ σ + δ σ + √ γ (cid:0) σ (cid:0) R † ( t ) + e iω τ L † ( t − τ ) (cid:1) + H.c. (cid:1) (A8)+ √ γ (cid:16) σ R † ( t − τ / e iω τ / + H.c. (cid:17) (A9)+ √ γ (cid:16) σ L † ( t − τ / e iω τ / + H.c. (cid:17) + √ γ (cid:0) σ (cid:0) R † ( t − τ ) e iω τ + L † ( t ) (cid:1) + H.c. (cid:1) , where τ = ( τ + τ ) /
2. In the following, we use thenotation for the collective states corresponding to | ijk (cid:105) = | ( i − + ( j − + ( k − (cid:105) , which leads to, e.g., σ = | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + | (cid:105)(cid:104) | + | (cid:105)(cid:104) | .To solve the corresponding Schr¨odinger equation, weswitch to a time-discrete evolution picture, and integratefrom 0 to ∆ as the first time step from | ψ (0) (cid:105) . | ψ ( n ) (cid:105) = U ( n, n − | ψ ( n − (cid:105) = exp (cid:34) − i (cid:126) (cid:90) n ∆( n − H I ( t (cid:48) ) dt (cid:48) (cid:35) | ψ ( n − (cid:105) (A10)We can now introduce in the discrete time-bin basis, thecollective bath operators as:∆ R † ( n ) = (cid:90) n ∆( n − (cid:90) dωr † ω e i ( ω − ω ) t dt √ π ∆ , (A11)∆ L † ( n ) = (cid:90) n ∆( n − (cid:90) dωl † ω e i ( ω − ω ) t dt √ π ∆ . (A12)Given the time bin dynamics, we can now study dif-ferent cases. Due to the richness of the model, we fo-cus in the following only on the identical emitter case: γ i = γ and ω = ω = ω = ω . In this case, wehave as free parameters only: γ, τ , τ , and the feed-back phase ω ∆ / φ = ω ( τ + τ ) / π , in its feedback interaction themiddle atom is subjected to the phase: φ = φ/ (1+ τ /τ )and φ = φ/ (1+ τ /τ ). For the symmetric case it followsobviously: φ , = φ/ [1] M. C. Booth, M. Atat¨ure, G. Di Giuseppe, B. E. A. Saleh,A. V. Sergienko, and M. C. Teich, Physical Review A , 023815 (2002). [2] H. J. Kimble, Nature , 1023 (2008).[3] D. Petrosyan and M. Fleischhauer, Physical Review Let-ters , 170501 (2008). [4] M. Hafezi, D. E. Chang, V. Gritsev, E. Demler, andM. D. Lukin, Physical Review A , 013822 (2012).[5] H. Zheng and H. U. Baranger, Physical Review Letters , 113601 (2013).[6] S. Zeeb, C. Noh, A. Parkins, and H. Carmichael, Physi-cal Review A , 023829 (2015).[7] G. Calaj´o, F. Ciccarello, D. Chang, and P. Rabl, PhysicalReview A , 033833 (2016).[8] I. S¨ollner, S. Mahmoodian, S. L. Hansen, L. Midolo,A. Javadi, G. Kirˇsansk˙e, T. Pregnolato, H. El-Ella, E. H.Lee, J. D. Song, et al. , Nature nanotechnology , 775(2015).[9] D. Roy, C. M. Wilson, and O. Firstenberg, Reviews ofModern Physics , 021001 (2017).[10] H. Pichler, S. Choi, P. Zoller, and M. D. Lukin, Proc.Natl. Acad. Sci. , 11362 (2017).[11] P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeu-tel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller,Nature , 473 (2017).[12] M. Bello, G. Platero, J. I. Cirac, and A. Gonz´alez-Tudela, Science Advances (2019).[13] A. Mekis, J. Chen, I. Kurland, S. Fan, P. R. Villeneuve,and J. Joannopoulos, Physical Review Letters , 3787(1996).[14] P. Colman, C. Husko, S. Combri´e, I. Sagnes, C. Wong,and A. De Rossi, Nature Photonics , 862 (2010).[15] A. Goban, C.-L. Hung, J. D. Hood, S.-P. Yu, J. A. Muniz,O. Painter, and H. J. Kimble, Physical Review Letters , 063601 (2015).[16] A. Gonz´alez-Tudela, C.-L. Hung, D. E. Chang, J. I.Cirac, and H. J. Kimble, Nature Photonics , 320 (2015).[17] T. Aoki, A. Parkins, D. Alton, C. Regal, B. Dayan, E. Os-tby, K. J. Vahala, and H. J. Kimble, Physical ReviewLetters , 083601 (2009).[18] S. Mahmoodian, K. Prindal-Nielsen, I. S¨ollner, S. Stobbe,and L. Peter, Optical Materials Express , 43 (2017).[19] J.-H. Kim, S. Aghaeimeibodi, C. J. K. Richardson, R. P.Leavitt, and E. Waks, Nano Letters , 4734 (2018).[20] E. Vetsch, D. Reitz, G. Sagu´e, R. Schmidt, S. Dawkins,and A. Rauschenbeutel, Physical Review Letters ,203603 (2010).[21] A. Goban, K. Choi, D. Alton, D. Ding, C. Lacroˆute,M. Pototschnig, T. Thiele, N. Stern, and H. Kimble,Physical Review Letters , 033603 (2012).[22] P. Solano, P. Barberis-Blostein, F. K. Fatemi, L. A.Orozco, and S. L. Rolston, Nature Communications ,1857 (2017).[23] D. H. White, S. Kato, N. N´emet, S. Parkins, and T. Aoki,Physical Review Letters , 253603 (2019).[24] S. Kato, N. N´emet, K. Senga, S. Mizukami, X. Huang,S. Parkins, and T. Aoki, Nature communications ,1160 (2019).[25] A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno,L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal,Physical Review Letters , 020501 (2011).[26] J. Christensen, A. Manjavacas, S. Thongrattanasiri,F. H. L. Koppens, and F. J. Garc´ı´a de Abajo, ACSnano , 431 (2011).[27] E. Berm´udez-Urena, C. Gonzalez-Ballestero, M. Geisel-mann, R. Marty, I. P. Radko, T. Holmgaard,Y. Alaverdyan, E. Moreno, F. J. Garc´ıa-Vidal, S. I.Bozhevolnyi, et al. , Nature Communications , 7883(2015). [28] M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B.Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang,and O. Painter, Nature , 692 (2019).[29] E. Shahmoon, P. Griˇsins, H. P. Stimming, I. Mazets, andG. Kurizki, Optica , 725 (2016).[30] A. Serafini, S. Mancini, and S. Bose, Phys. Rev. Lett. , 010503 (2006).[31] V. I. Yudson and P. Reineker, Phys. Rev. A , 052713(2008).[32] D. E. Chang, L. Jiang, A. V. Gorshkov, and H. J. Kim-ble, New Journal of Physics , 063003 (2012).[33] A. F. Van Loo, A. Fedorov, K. Lalumiere, B. C. Sanders,A. Blais, and A. Wallraff, Science , 1494 (2013),1407.6747.[34] C. Gonzalez-Ballestero, F. J. Garc´ıa-Vidal, andE. Moreno, New Journal of Physics , 073015 (2013).[35] P. Facchi, M. S. Kim, S. Pascazio, F. V. Pepe, D. Po-marico, and T. Tufarelli, Phys. Rev. A , 043839(2016).[36] I. M. Mirza and J. C. Schotland, Phys. Rev. A , 012309(2016).[37] G.-Z. Song, E. Munro, W. Nie, L.-C. Kwek, F.-G. Deng,and G.-L. Long, Physical Review A , 023814 (2018).[38] A. H. Kiilerich and K. Mølmer, Phys. Rev. Lett. ,123604 (2019).[39] E. V. Stolyarov, Physical Review A , 023857 (2019).[40] F. Dinc, A. M. Bra´nczyk, and I. Ercan, arXiv preprintarXiv:1809.05164 (2018).[41] E. V. Denning, J. Iles-Smith, and J. Mork, arXivpreprint arXiv:1907.13423 (2019).[42] Y.-L. L. Fang, H. U. Baranger, et al. , Physical Review A , 053845 (2015).[43] G. Calaj´o, Y.-L. L. Fang, H. U. Baranger, F. Ciccarello, et al. , Physical Review Letters , 073601 (2019).[44] S. Hughes and G. S. Agarwal, Physical Review Letters , 063601 (2017).[45] Z. Liao, X. Zeng, H. Nha, and M. S. Zubairy, PhysicaScripta , 063004 (2016).[46] T. F. See, C. Noh, and D. G. Angelakis, Physical ReviewA , 053845 (2017).[47] M. P. Schneider, T. Sproll, C. Stawiarski, P. Schmitteck-ert, and K. Busch, Physical Review A , 013828 (2016).[48] T. Shi, D. E. Chang, and J. I. Cirac, Physical Review A , 053834 (2015).[49] B. Q. Baragiola, R. L. Cook, A. M. Bra´nczyk, andJ. Combes, Physical Review A , 013811 (2012).[50] S. Xu and S. Fan, Physical Review A , 043845 (2015).[51] S. Fan, S¸. E. Kocaba¸s, and J.-T. Shen, Physical ReviewA , 063821 (2010).[52] K. A. Fischer, R. Trivedi, V. Ramasesh, I. Siddiqi, andJ. Vuˇckovi´c, Quantum , 69 (2018).[53] T. Caneva, M. T. Manzoni, T. Shi, J. S. Douglas, J. I.Cirac, and D. E. Chang, New Journal of Physics ,113001 (2015).[54] D. Witthaut and A. S. Sørensen, New Journal of Physics , 043052 (2010).[55] H. Zheng, D. J. Gauthier, and H. U. Baranger, PhysicalReview A , 063816 (2010).[56] H. Zheng, D. J. Gauthier, and H. U. Baranger, PhysicalReview A , 043832 (2012).[57] J.-T. Shen and S. Fan, Physical Review A , 062709(2007).[58] V. Paulisch, H. Kimble, and A. Gonz´alez-Tudela, NewJournal of Physics , 043041 (2016). [59] H. Pichler and P. Zoller, Physical Review Letters ,093601 (2016).[60] T. Shi, Y.-H. Wu, A. Gonz´alez-Tudela, and J. I. Cirac,Physical Review X , 021027 (2016).[61] C. Gardiner, P. Zoller, and P. Zoller, Quantum noise:a handbook of Markovian and non-Markovian quantumstochastic methods with applications to quantum optics ,Vol. 56 (Springer Science & Business Media, 2004).[62] H.-P. Breuer, F. Petruccione, et al. , The theory of openquantum systems (Oxford University Press on Demand,2002).[63] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini,Reviews of Modern Physics , 021002 (2016).[64] I. De Vega and D. Alonso, Reviews of Modern Physics , 015001 (2017).[65] A. Carmele and S. Reitzenstein, Nanophotonics , 655(2019).[66] M. Strauß, A. Carmele, J. Schleibner, M. Hohn,C. Schneider, S. H¨ofling, J. Wolters, and S. Reitzenstein,Physical review letters , 107401 (2019).[67] S. Bounouar, M. Strauß, A. Carmele, P. Schnauber,A. Thoma, M. Gschrey, J.-H. Schulze, A. Strittmatter,S. Rodt, A. Knorr, and S. Reitzenstein, Phys. Rev. Lett. , 233601 (2017).[68] Y. Su, D. Bimberg, A. Knorr, and A. Carmele, Phys.Rev. Lett. , 113604 (2013).[69] K. Sinha, P. Meystre, E. A. Goldschmidt, F. K.Fatemi, S. L. Rolston, and P. Solano, arXiv preprintarXiv:1907.12017 (2019).[70] T. Ramos, B. Vermersch, P. Hauke, H. Pichler, andP. Zoller, Physical Review A , 062104 (2016).[71] C. Gardiner, A. Parkins, and P. Zoller, Physical ReviewA , 4363 (1992).[72] H. J. Carmichael, Statistical methods in quantum optics2: Non-classical fields (Springer Science & Business Me-dia, 2009).[73] P.-O. Guimond, H. Pichler, A. Rauschenbeutel, andP. Zoller, Phys. Rev. A , 033829 (2016). [74] Y. Lu, N. L. Naumann, J. Cerrillo, Q. Zhao, A. Knorr,and A. Carmele, Physical Review A , 063840 (2017).[75] L. Droenner, N. L. Naumann, E. Sch¨oll, A. Knorr, andA. Carmele, Physical Review A , 023840 (2019).[76] E. S´anchez-Burillo, L. Mart´ın-Moreno, J. J. Garc´ıa-Ripoll, and D. Zueco, Phys. Rev. Lett. , 013601(2019).[77] C. Sch¨on, K. Hammerer, M. M. Wolf, J. I. Cirac, andE. Solano, Physical Review A , 032311 (2007).[78] C. Cascio, J. C. Halimeh, I. P. McCulloch, A. Recati,and I. de Vega, Phys. Rev. A , 013845 (2019).[79] S. Kreinberg, T. Grbeˇsi´c, M. Strauß, A. Carmele, M. Em-merling, C. Schneider, S. H¨ofling, X. Porte, and S. Re-itzenstein, Light: Science & Applications , 1 (2018).[80] N. N´emet, A. Carmele, S. Parkins, and A. Knorr, Phys.Rev. A , 023805 (2019).[81] A. Carmele, B. Vogell, K. Stannigel, and P. Zoller, NewJournal of Physics , 063042 (2014).[82] N. N´emet and S. Parkins, Physical Review A , 023809(2016).[83] A. Carmele, J. Kabuss, F. Schulze, S. Reitzenstein, andA. Knorr, Physical Review Letters , 013601 (2013).[84] J. Kabuss, F. Katsch, A. Knorr, and A. Carmele, JOSAB , C10 (2016).[85] S. Whalen, A. Grimsmo, and H. Carmichael, QuantumScience and Technology , 044008 (2017).[86] J. Kabuss, D. O. Krimer, S. Rotter, K. Stannigel,A. Knorr, and A. Carmele, Phys. Rev. A , 053801(2015).[87] A. L. Grimsmo, Physical Review Letters , 060402(2015).[88] L. Scarpelli, B. Lang, F. Masia, D. M. Beggs, E. A. Mul-jarov, A. B. Young, R. Oulton, M. Kamp, S. H¨ofling,C. Schneider, and W. Langbein, Phys. Rev. B ,035311 (2019).[89] M. Arcari, I. S¨ollner, A. Javadi, S. Lindskov Hansen,S. Mahmoodian, J. Liu, H. Thyrrestrup, E. H. Lee, J. D.Song, S. Stobbe, and P. Lodahl, Phys. Rev. Lett.113