Delay-induced spontaneous dark state generation from two distant excited atoms

We investigate the collective non-Markovian dynamics of two fully excited two-level atoms coupled to a one-dimensional waveguide in the presence of delay. We demonstrate that analogous to the well-known superfluorescence phenomena, where an inverted atomic ensemble synchronizes to enhance its emission, there is a `subfluorescence' effect that synchronizes the atoms into an entangled dark state depending on the interatomic separation. Our results are pertinent to long-distance quantum networks, presenting a mechanism for spontaneous entanglement generation between distant quantum emitters.


I. INTRODUCTION
Large-scale interconnected quantum systems offer promising applications in quantum information processing and distributed quantum sensing [1][2][3].As the experimental capabilities for directly interconnecting distant quantum nodes grow [4], we have yet to develop the theoretical toolbox to efficiently analyze systems of multiple qubits collectively exchanging photons over large distances, where delay effects cannot be neglected.Such many-body systems exhibit rich non-Markovian and nonlinear dynamics, arising from delayed and multiphoton interactions.Given the complexity associated with such systems, even the simplest scenarios, e.g., the spontaneous decay of two fully excited distant atoms, remains unexplored.
The most studied cases of coupled emitters consider nearby atoms, neglecting the delay time the field needs to propagate between them.This scenario is a landmark of quantum optics across platforms, such as atoms in free space [5,6], inside a leaky cavity [7], and near waveguides [8].In all of these cases, ensembles of initially fully excited two-level atoms decay into their ground state.As atoms decay, correlations among them spontaneously and momentarily emerge, leading to the well-known collective phenomenon of superfluorescence [6,[9][10][11] [12].Such atom-atom correlations, which are absent in the initial state, emerge without external driving fields or post-selection [5].
Recent studies of the effects emerging from delayed interactions considered either many atoms with a single photon [41,45,[53][54][55], many photons with a single atom [44,[56][57][58], or numerical exploration of threephotons and three-atoms [50], inviting us to revisit the FIG. 1. Schematic of the system: two excited atoms with a transition frequency ω0 near a waveguide at positions z1,2 = ±d/2.The emission rate of each atom into the waveguide is γ and the propagation delay of the field between the two atoms is τ = d/v, with v as the velocity of the EM field in the waveguide.

arXiv:2303.06559v2 [quant-ph] 5 Jun 2024
theoretical description of other canonical quantum optical phenomena, such as superfluorescence, in the context of waveguide-coupled emitters in a large-scale quantum network.
In this paper, we show that a system of twoatoms, with two-excitations, and delayed-interactions (see Fig. 1) exhibits novel non-Markovian collective dynamics.In particular, for fully inverted atoms, the instantaneous decay can be faster than in superfluorescence and steady-state entanglement can suddenly emerge.In the latter case, the electromagnetic (EM) field can have two photons trapped between the two atoms, an effect that can not be captured by Markovian dynamics.We refer to this phenomenon, which leads to the spontaneous generation of a dark state, as delayed-induced subfluorescence -the subradiant counterpart of the well-known superfluorescence effect.The resulting steady state is a delocalized hybrid atom-photon state, requiring a description beyond the usual Born and Markov approximations and qualitatively different from the case without delay, where the system always ends in the ground state.The rich and unique dynamics in such a seemingly simple system demonstrates non-Markovian delay as a mechanism for trapping light, entanglement generation and stabilization in long-distance interacting quantum nodes.

II. THEORETICAL MODEL
We consider two two-level atoms, with transition frequency ω 0 between the ground |g⟩ and excited |e⟩ states, coupled through a waveguide.In this work, our emphasis is on the case of perfectly coupled atoms, the free Hamiltonian of the total atoms+field system is given by The interaction Hamiltonian in the interaction picture is (having made the electric-dipole and rotating-wave approximations): where g α is the coupling between atoms and the guided mode α with frequency ω α , ∆ α = ω α − ω 0 , and k αη = ηω α /v, with η = ±1 specifying the direction of propagation.The position of the atoms is z 1, 2 = ±d/2, σ(r) z = |e⟩ r ⟨e| r −|g⟩ r ⟨g| r and the raising and lowering operators for the r-th atom are defined as σ(r) We consider the initial state with both atoms excited and the field in the vacuum state, |e 1 e 2 , {0}⟩.As a consequence of the rotating-wave approximation, the total number of atomic and field excitations are conserved, suggesting the following ansätz for the state of the system at time t: where |g 1 g 2 , {0}⟩ is the ground state of the system.The complex coefficients a(t), b (r) αη (t), and c αη,βη ′ (t) correspond to the probability amplitudes of having an excitation in both the atoms, r th atom and field mode {α, η}, and field modes {α, η} and {β, η ′ }, respectively.The coefficients c αη (t) represent the probability amplitude of exciting two photons in modes {α, η}.The Hamiltonian in Eq. (1) models two emitters interacting via the quantized electromagnetic field; their distance, codified in the phase e ikαηzr , has the effect of introducing a delay in the equations of motion for the quantum state coefficients.This delayed interaction between the emitters results from the finite propagation speed of light.
αη (t)e −i∆αt and C αη = c αη (t)e −2i∆αt , and formally integrating the equation for c αη,βη ′ (t), the Schrödinger equation yields the following system of delay-differential equations for the excitation amplitudes: where γ is the decay rate of the emitters into the guided modes.We present a detailed derivation in Appendix A.

III. SYSTEM DYNAMICS
The Laplace transform of a(t) is (see Appendix B) where ϕ = ω 0 τ is the resonance field propagation phase.We write the sum over modes on the right-hand side (RHS) of Eq. ( 6) as an integral over frequencies by introducing the density of modes ρ(ω α ).Additionally, we use the Wigner-Weisskopf approximation and set ρ(ω α ), g α ≈ ρ(ω 0 ), g 0 , evaluated at the atomic resonance frequency, a good approximation in the presence of delay effects in waveguide QED [40].With these considerations we obtain ã(s) = 1/(s + γ) for all τ .Taking the inverse Laplace transform we get a(t) = e −γt , which gives the time-dependent probability of having two atomic excitations We remark that the probability of having two atomic excitations is independent of the delay between atoms.The probability of having only one of the atoms excited is (t) 2. Probability P (1) (t) of having only one of the atoms excited for a delay of γτ = 0.895 compared to that of coincident and independent atoms (γτ = 0 and γτ → ∞, respectively).We note that for a propagation phase of ϕ = nπ, there is a formation of a steady atom-photon bound state as indicated in Eq. (10).
where the time-dependent solutions b (1,2) αη (t) are derived in Appendix B and given by (B21).In the limits of two coincident atoms with τ = 0, and two infinitely distant atoms with τ → ∞, we obtain the expected Markovian dynamics [59,60]: Figure 2 shows P (1) (t) for different values of the delay τ and the propagation phase ϕ.For times t < τ , the atoms decay independently.At time t = τ , the field emitted by one atom reaches the other, modifying their decay dynamics depending on the value of the propagation phase ϕ.Their instantaneous decay rate, given by the negative slope of P (1) (t), could momentarily exceed that of standard superfluorescence (τ = 0).The delay condition γτ ≈ 0.375 and ϕ = nπ maximizes the instantaneous decay rate.This phenomenon is a signature of superduperradiance, reported in Ref. [41].Notably, in this case the atom-atom coherence that is necessary to modify the instantaneous decay emerges spontaneously.
In the late-time limit, the following steady state appears when ϕ = nπ for any delay τ between the atoms (see Appendix C): where we have removed the explicit time dependence of the steady-state coefficients.The steady-state amplitude for the mode {α, η} is b and c αη,βη ′ is given in Eq. (C4)in Appendix C. The last two terms in Eq. ( 10) represent a bound state in the continuum (BIC) that corresponds to having one shared excitation between the atoms and one propagating photon mode in between them [44].Using the Born rule and Eq. ( 11), we obtain that its probability is 2P (1) (t → ∞), where: We note that the probability of ending up in a BIC state is maximum for γτ ≈ 0.895 with a probability of 0.282.
The fact that the atoms+field state is non-separable shows that one cannot use the Born-Markov approximation to solve for the dynamics of the system.Using Eq. ( 10) and Eq. ( 11) we obtain the following reduced density matrix for the atomic subsystem ) with + corresponding to the case when n is odd and − when it is even and |Ψ ± ⟩ = (|e 1 g 2 ⟩±|g 1 e 2 ⟩)/ √ 2 are single excitation Bell states.As the initially inverted atoms evolve into a superposition of radiative and non-radiative states upon delayed-interactions, the atoms decay into a superposition of ground state and an entangled dark steady state.For a field propagation phase ϕ = 2nπ (ϕ = (2n + 1)π), the dark state that appears in the latetime limit corresponds to the antisymmetric state |Ψ − ⟩ (symmetric state |Ψ + ⟩).
Delayed interactions create quantum correlations between the atoms.We quantify it using the concurrence C(t) of the reduced density matrix of the atoms ρA (t) = T r F (|ψ(t)⟩⟨ψ(t)|) [61].For τ → ∞, the concurrence is zero throughout the evolution since the initially uncorrelated atoms evolve independently.Remarkably, when τ = 0, the concurrence is also zero throughout the evolution, even when atoms transition to a superradiant behavior.This case exemplifies that entanglement is not necessary for superradiance.For intermediate values of τ , we numerically study the dynamics of concurrence, shown in Fig. 3 (a).We begin studying its behavior for γτ ≈ 0.895 (corresponding to the maximum value of P (1) (t → ∞)) and γτ ≈ 0.375 (corresponding to the largest instantaneous decay rate).Since the system lacks initial correlations, the concurrence is zero from t = 0 to a certain time, t SBE , when there is a sudden birth of entanglement (SBE) [62][63][64][65][66][67][68].For ϕ = nπ, the concurrence increases until reaching a stationary value, whereas for ϕ = (n+ 1  2 )π, it always remains zero.In general, for other values of ϕ, the concurrence suddenly departs from zero and slowly decays after reaching a maximum value.For γτ ≈ 0.895, we obtain the maximum value for the concurrence.Figure 3 (b) shows the emergence of atom-atom correlations defined by Tr ρA (t)σ (1) + σ(2) − .We note that the atom-atom correlations develop as soon as the atoms  12).(b) Correlation between atomic dipoles.The atoms decay independently for times t < τ .At t = τ , we observe an emergence of quantum correlations between the atoms.When ϕ = mπ, the atoms reach a dark steady state for all delays τ .The maximum concurrence is ≈ 0.141 at γτ ≈ 0.895.
"see" each other, while the concurrence takes longer to emerge.Delayed atom-atom interactions couple the fully excited state of the atoms to both single excitation symmetric and antisymmetric states.After the build-up of correlations, the most radiative state collectively decays, while the non-radiative atomic state remains.The system, therefore, spontaneously evolves into an entangled dark state in the presence of retardation.Such appearance of quantum correlations is a striking example of environment-assisted spontaneous entanglement generation.

IV. FIELD DYNAMICS
We now focus on the dynamics of the field sourced by the non-Markovian behavior of the atoms.Field operators are often described through mode expansion, as shown in Eq. (D1), or Green's functions [69].Its expectation values involve integrating over all the modes or knowing the Green's function for the field, which for waveguides is usually given by a mode expansion [70].When the system is Markovian, an input-output relation simplifies the calculation of field observables, and the sum over field modes is replaced by a sum over atomic operators: Ê+ (r) = Ê+ Input (r) + i f (r, r i )σ i − , where the positive part of the field operator, Ê+ (r), at r, is proportional to the input field, Ê+ Input (r), plus the field sourced by the atomic operators σi − .Here r i is the position of atom i, and f (r, r i ) is a complex function that depends on the modes of the system [14,71].
The probability of measuring two photons at the same time at position z is proportional to the second-order correlation function [72] g (2) where we take the maximum over the distance between the atoms, in the normalization constant, to compare the second-order correlation function for different time lags.By substituting the input-output relation for the Markovian case into g (2) and replacing it with the stationary solution, Eq. ( 10), we see that the dynamics does not lead to a two-photon bound state.However, the system under consideration is non-Markovian, and one must calculate the field observables from the full expression for the field operator.
Figure 4 shows the second-order correlation function in the stationary regime (t → ∞) (see Appendix D for details), for different distances between the atoms.It can be seen that two photons are trapped between the atoms, a result that is distinctly non-Markovian, demonstrating the rich interplay of non-Markovianity and nonlinear atom-photon interactions.Although the maximum probability of having the atoms excited is given at a distance γτ ≈ 0.895, the probability of having two photons excited at a particular z is greater at γτ ≈ 0.375, where the instantaneous decay is maximum.This is explained by the fact that the probability of having two photons trapped in a particular position increases when the probability of the atomic excitation diminishes, because the energy from the atomic excitation should go to the field.
A salient question pertains to the maximum probability of having two photons between the atoms.This probability depends on the value of the autocorrelation function and the distance between the two atoms.Figure 5 depicts γτ g (2) (z in ), which is proportional to the probability of having two photons in the waveguide region separating the two atoms as a function of the delay.The distance between the atoms where the probability of having two photons is maximum is γτ ≈ 1.0.Note that we chose a normalization that does not depend on the transverse profile (see Appendix D for details).To obtain the actual probability that two photon detection happen at the same time and position, it is useful to calculate the field energy in the guided field modes, in the limit t → ∞ (see Appendix C) where P (1) (t → ∞) is the probability of having one atomic excitation in the long time limit (see Eq. ( 12)).This result states that the total field energy corresponds to the energy of two photons with frequency ω 0 (the total energy) minus the energy accumulated in the atoms that has not been emitted into the waveguide.The normalization constant, dz⟨ Ê(+) (z) Ê(−) (z)⟩, is given by the total energy around the fiber, E total (see Eq. ( 15)) multiplied by the square of the mode profile at the atomic frequency transition ω 0 .

V. SUMMARY AND OUTLOOK
We have analyzed the spontaneous decay of two fully inverted atoms coupled through a waveguide in the presence of retardation effects.A remarkable result is the spontaneous creation of a steady delocalized atomphoton bound state, with sudden birth of entanglement between the atoms.Furthermore, we demonstrate that such a delay can create two-photon bound states, wherein one can have two photons in the waveguide region between the two atoms.Such states appear as a result of the non-Markovian time-delayed feedback of the spontaneously radiated EM field acting on the emitters.Addi- tionally, the collective decay of the two atoms can be momentarily enhanced beyond standard superfluorescence and subsequently inhibited, demonstrating the rich non-Markovian dynamics of such a system.
Such delay-induced spontaneous steady state entanglement generation can have implications in the rapidly growing field of waveguide QED, a field with promising applications in quantum information processing [3,[73][74][75][76][77][78][79] that benefit from preparing and manipulating longlived dark-states.In this context, there have been several proposals to generate a steady entangled state; compared to delay-induced subfluorescence, these schemes necessitate extra degrees of control, such as external driving fields [36,[80][81][82], initial entanglement [41], or chiral emission in front of a mirror [39,83].
It is not only the generation [84] but also the stabilization of entangled states that are critical to developing efficient quantum devices.In contrast to the idea that the interaction of quantum systems with their environment leads to decoherence and can degrade the entanglement between the components of a quantum system [85], the environment can also be proposed as a generator [86,87] and stabilizer of entanglement [88].Our results demonstrate that non-Markovian time-delayed feedback can be a mechanism for environment-assisted entanglement generation and stabilization.
Adding delay to the most straightforward collective system of two two-level atoms leads to different phenomenology, breaking the Born and Markov approximations, non-trivially modifying the dynamics, and spontaneously creating steady-state quantum correlations between the atoms and multiphoton bound states.Our results further the understanding of the significant role delay plays in quantum optics and present the outset of studying more complex phenomena involving many-body interactions [6,89,90].However, studying this scenario is challenging as the complexity of the problem increases with the number of emitters, and known methods based on master equation approaches fail because neither the Markov nor the Born approximations are valid.Furthermore, although we neglect dispersive dipole-dipole interactions and coupling to non-guided modes to highlight the consequences of delay-induced non-Markovianity, all these effects can coexist, creating richer dynamics that will increase in complexity as the system scales up.

VI. ACKNOWLEDGMENTS
We acknowledge helpful discussions with Saikat Guha and Elizabeth Goldschmidt.P.S. is a CIFAR Azrieli Global Scholar in the Quantum Information Science Program.This work was supported by DGAPA-UNAM under grant IG101421 and IG101324 from Mexico, as well as CONICYT-PAI grant 77190033, and FONDE-CYT grant N • 11200192 from Chile.K.S. acknowledges funding support from the John Templeton Foundation Award No. 62422, the National Science Foundation under Grant No. PHY-2309341 and the Army Research Office through Grant No. W911NF2410080.

Appendix A: Equations of motion: Derivation
Using Schrödinger equation with the Hamiltonian (1) and the state ansätz (2), we get the following differential equations for the probability amplitudes where we consider atomic positions such that z 1 = −z 2 .
We integrate (A3) with c αη,βη ′ (0) = 0, and substitute in (A2), obtaining ḃ(r) For t ̸ = 0, the second term of the RHS can be approximated as zero assuming g β constant and b βη ′ (t) evolving slowly, which is consistent with the Wigner-Weisskopf approximation.We transform the sum over frequencies to integrals by using the densities of modes ρ(w α ) [91].For guided modes k αη = η ωα v being the phase velocity of the field inside the waveguide, and η = ±1 labels forward or backward propagation direction of the field along the waveguide.Thus α,η → η=±1 ∞ 0 dω α ρ(ω α ).To study the non-Markovian effects due only to delay and not to a structured reservoir, we assume a flat spectral density of field modes around the resonance of the emitters such that ω α ≈ ω 0 in the evaluation of ρ(ω α ) and g α functions.
Taking into account all the considerations above we obtain where τ rs = zr−zs v .We make use of the Sokhotski-Plemelj theorem where PV refers to the Cauchy principal value.Absorbing the contribution of the principal value (which corresponds to the Lamb shifts) into the atomic transition frequency [5] we obtain where the single atom decay rate to guided modes is defined as Using (A6) we obtain the equation of motion for ḃ(r) , ϕ = ω 0 τ and Θ is the Heaviside step function.

Appendix B: Equations of motion: Solution
We use the Laplace transform to solve the delayed differential equations of motion (3)- (5).As an intermediate step, we define the variables where we use that αη (0) = 0 and C αη (0) = 0 as initial conditions.The solutions of the previous equations are

(B7)
In order to find the inverse Laplace transform of Eqs.(B5), (B6) and (B7) we use |g α | 2 /γ ∼ 10 −4 for EM modes in the visible range [92] to neglect the contribution of the term 4|g α | 2 /(s + i∆ α ) in the denominator.This corresponds to neglecting the probability of exciting two equal α modes traveling in the same direction compared to the probability of exciting two different modes.Considering this, we get Eq.( 6) and where we extend the lower limit in the integral to −∞ as an approximation.We rewrite the denominator inside the integral using the poles for the variable ∆ α , determined by the characteristic equation, with W k denoting the kth branch of Lambert W function [93] and σ = ±1.For simplicity, we introduce r = γτ 2 e γτ /2+iϕ .Using partial fraction decomposition [94] we obtain Using Cauchy's integral formula we get For the second integral we use that Im ∆ (σ) k > 0 because we can take Re(s) as large as we need, obtaining Using the following identities of the Lambert functions [95] ∞ we obtain I(s, γ, τ, ϕ) = −π.

Inverse Laplace transform of B(1,2)
αη To obtain B (1,2) αη (t) we apply the inverse Laplace transform of Eq. (B8).First, similar to what is done in Appendix B 1, we rewrite the result using the poles of the denominator but for the variable s.Defining and using that k αη d = η(∆ α τ + ϕ) and ã(s) = 1/(s + γ), we obtain where L −1 denotes the inverse Laplace transform.Using that and taking b where W k represents the k-th branch of the Lambert-W function.
To find the function in the time domain, we substitute Eq. (B8) and use the poles (B17).Then, applying the inverse Laplace transform, we arrive at the expression where we have defined the functions These functions are not invariant to the interchange of labels α, β.

Appendix C: Stationary values
In the late-time limit γt → ∞ we have for (B21) the value b ( The above expression does not decay to zero if Thus, the terms of the expression that will be non-zero in the long time limit are k = 0, and {σ = +1, ϕ = 2πn} or {σ = −1, ϕ = (2n + 1)π}, with n ∈ N 0 , because Taking into account the above, we get a stationary value given by b Then, we obtain which leads to Eq. (12).For correlations between atomic dipoles, we find the stationary value Tr ρA (t → ∞)σ Finally, the function c αη,βη ′ (t) from (B23) has a simplified expression in the long-time limit γt → ∞ and ϕ = nπ.In this scenario, the value is equal to 2 r=1 e −i(η+uη ′ )ϕzr e −i(η∆α+uη ′ ∆ β )zrτ  where n 2 (r) is the refractive index of the cylindrical waveguide.
To calculate the correlation function G (2) (z, z) = ⟨ Ê(−) (z, t) Ê(−) (z, t) Ê(+) (z, t) Ê(+) (z, t)⟩, for the position z along the waveguide and in the long-time limit, we use the Wigner-Weisskopf approximation (see Appendix A) to evaluate E α ≈ E 0 and e α ≈ e 0 , and the state (10).We get the following, G (2) FIG. 3. (a)Concurrence as a function of time.The atoms do not have initial correlations C(0) = 0, after some time t SBE there is a sudden birth of entanglement".For ϕ = (m + 1/2)π C(t) = 0 for all time.If ϕ = mπ concurrence reaches a stationary value Css, given by Eq. (12).(b) Correlation between atomic dipoles.The atoms decay independently for times t < τ .At t = τ , we observe an emergence of quantum correlations between the atoms.When ϕ = mπ, the atoms reach a dark steady state for all delays τ .The maximum concurrence is ≈ 0.141 at γτ ≈ 0.895.

FIG. 4 .
FIG.4.Second order correlation function as a function of the distance between the atoms.The vertical lines serve to show the location of the atoms.

FIG. 5 .
FIG.5.Probability of detecting two photons in the same position, between the atoms, as a function of the delay, for the steady state case ϕ = mπ.