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

W. Alvarez-Giron<EMAIL_ADDRESS>Instituto de Investigaciones en
Matemáticas Aplicadas y en Sistemas, Universidad Nacional Autónoma de México,
Ciudad Universitaria, 04510, DF, México.  P. Solano<EMAIL_ADDRESS>Departamento de Física, Facultad de Ciencias Físicas y Matemáticas,
Universidad de Concepción, Concepción, Chile  CIFAR Azrieli Global Scholars
program, CIFAR, Toronto, Canada.  K. Sinha<EMAIL_ADDRESS>School of
Electrical, Computer and Energy Engineering, Arizona State University, Tempe,
AZ 85287-5706, USA  P. Barberis-Blostein<EMAIL_ADDRESS>Instituto de
Investigaciones en Matemáticas Aplicadas y en Sistemas, Universidad Nacional
Autónoma de México, Ciudad Universitaria, 04510, DF, México.

###### Abstract

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.

Introduction.— Coupling quantum emitters via waveguides is essential for
building large-scale and interconnected quantum systems, which find promising
applications in quantum information processing and distributed quantum sensing
Kimble (2008); Duan and Monroe (2010); Monroe  _et al._ (2014a). Platforms
based on waveguide QED enable tunable, efficient, and long-ranged interactions
between quantum emitters in the optical and microwave regime Goban  _et al._
(2015); Solano  _et al._ (2017a, b); Kim  _et al._ (2018); Wen  _et al._
(2019); Mirhosseini  _et al._ (2019); Sheremet  _et al._ (2021); Magnard  _et
al._ (2020); Tiranov  _et al._ (2023). For example, state-of-the-art
experiments allow one to tune out dispersive dipole-dipole interactions by
strategically positioning the atoms along a waveguide Martín-Cano  _et al._
(2011); van Loo  _et al._ (2013); Pichler  _et al._ (2015) while highly
reducing coupling to non-guided modes Lecamp  _et al._ (2007); Astafiev  _et
al._ (2010); Arcari  _et al._ (2014); Zang  _et al._ (2016); Scarpelli  _et
al._ (2019).

However, long-ranged emitter-emitter coupling in waveguide QED systems
necessitates careful consideration of time-delayed interactions. When the time
taken by a photon to propagate between two distant emitters becomes comparable
to their characteristic lifetimes, the system exhibits surprisingly rich
delay-induced non-Markovian dynamics Zheng and Baranger (2013a); Tufarelli
_et al._ (2014); Guimond  _et al._ (2016); Zhang  _et al._ (2020a); Del Ángel
_et al._ (2022); Solano  _et al._ (2023). Some examples of such dynamics
include collective spontaneous emission rates exceeding those of Dicke
superradiance and formation of highly delocalized atom-photon bound states
Sinha  _et al._ (2020a, b); Dinç and Brańczyk (2019); Calajó  _et al._ (2019);
Facchi  _et al._ (2019); Guo  _et al._ (2020); Yao and Hughes (2009); Crowder
_et al._ (2020); Lee  _et al._ (2023). In addition, time-delayed feedback can
assist in preparing photonic cluster states Pichler  _et al._ (2017) and
single photon sources with improved coherence and indistinguishability Crowder
_et al._ (2023).

Such novel effects invite one to revisit canonical quantum optical phenomena
in the context of waveguide-coupled emitters in a quantum network. For
example, a striking result in collective quantum optics is the spontaneous
emergence of correlations between atoms in a fully inverted ensemble, which
leads to the well-known effect of superfluorescence Bonifacio and Lugiato
(1975); Glauber and Haake (1978); Vrehen  _et al._ (1980); Clemens  _et al._
(2003)111 We use the term superfluorescence Bonifacio and Lugiato (1975) as a
particular case of superradiance, wherein the collective decay rate of an
atomic ensemble is enhanced beyond individual spontaneous emission rate of a
single atom Dicke (1954); Gross and Haroche (1982).. Such atom-atom
correlations emerge without external driving fields, post-selection, or
initial state preparation Gross and Haroche (1982). Thus, the question arises:
can we spontaneously generate and engineer entanglement between distant
emitters in the presence of delay?

Figure 1: Schematic of the system: Two excited atoms with a transition
frequency $\omega_{0}$ near a waveguide at positions $z_{1,2}=\pm d/2$. The
emission rate of each atom into the waveguide is $\gamma$, the propagation
delay of the field between the two atoms is $\tau=d/v$, with $v$ as the
velocity of the EM field in the waveguide.

In this work, we demonstrate that one can spontaneously generate and stabilize
quantum entanglement between two initially excited emitters in the presence of
retardation. We investigate the non-Markovian dynamics of the emitters coupled
via a waveguide, illustrating that for interatomic distances of an integer
multiple of half their resonant wavelength, $\lambda/2$, the system ends in a
dark state where atom-atom entanglement suddenly emerges. We refer to the
phenomena as delayed-induced subfluorescence, as 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. Moreover, the probability of ending in such a
dark state increases linearly for small interatomic separation, reaching its
maximum when the atoms are separated by a distance $d\approx 0.895v/\gamma$,
with $v$ as the speed of light in the waveguide and $\gamma$ as the individual
atomic spontaneous emission rate. For interatomic separations comparable to
the coherence length associated with a spontaneously emitted photon
$(v/\gamma)$ memory effects of the electromagnetic environment become
prominent, necessitating a non-Markovian treatment of the system dynamics.

Theoretical model.— We consider two two-level atoms, with transition frequency
$\omega_{0}$ between the ground $\left|g\right\rangle$ and excited
$\left|e\right\rangle$ states, coupled through a waveguide. The free
Hamiltonian of the total atoms+field system is given by
$H_{0}=(\hbar\omega_{0}/2)\sum_{r=1}^{2}\hat{\sigma}_{z}^{(r)}+\sum_{\alpha,\eta}\hbar\omega_{\alpha}\hat{a}^{\dagger}_{\alpha,\eta}\hat{a}_{\alpha,\eta}$,
where $\omega_{\alpha}$ represents the frequency of the field modes, and
$\eta$ corresponds to the direction of propagation. The interaction
Hamiltonian in the interaction picture is (having made the electric-dipole and
rotating-wave approximations):

$\displaystyle\hat{H}_{\text{int}}=-i\hbar\sum_{r=1}^{2}\sum_{\alpha,\eta}g_{\alpha}\hat{\sigma}_{+}^{(r)}\hat{a}_{\alpha,\eta}e^{-i\Delta_{\alpha}t}e^{ik_{\alpha\eta}z_{r}}+\text{H.c.}\,,$
(1)

where $g_{\alpha}$ is the coupling between atoms and the guided mode $\alpha$
with frequency $\omega_{\alpha}$,
$\Delta_{\alpha}=\omega_{\alpha}-\omega_{0}$, and
$k_{\alpha\eta}=\eta\omega_{\alpha}/v$, with $\eta=\pm 1$ specifying the
direction of propagation. The position of the atoms is $z_{1,\,2}=\pm d/2$,
and the raising and lowering operators for the $r$-th atom are defined as
$\hat{\sigma}_{+}^{(r)}=\left(\hat{\sigma}_{-}^{(r)}\right)^{\dagger}=\left|e\right\rangle_{r}\left\langle
g\right|_{r}$.

We consider the initial state with both atoms excited and the field in the
vacuum state, $\left|e_{1}e_{2},\\{0\\}\right\rangle$. 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$:

$\displaystyle\left|\psi(t)\right\rangle=\bigg{\\{}a(t)\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{+}^{(2)}+\sum_{r=1}^{2}\sum_{\alpha,\eta}b_{\alpha\eta}^{(r)}(t)\hat{\sigma}_{+}^{(r)}\hat{a}_{\alpha,\eta}^{\dagger}+\mathop{\sum\sum}_{\alpha,\eta\neq\beta,\eta^{\prime}}\frac{c_{\alpha\eta,\beta\eta^{\prime}}(t)}{2}\hat{a}_{\alpha,\eta}^{\dagger}\hat{a}_{\beta,\eta^{\prime}}^{\dagger}+\sum_{\alpha,\eta}\frac{c_{\alpha\eta}(t)}{\sqrt{2}}\hat{a}_{\alpha,\eta}^{\dagger}\hat{a}_{\alpha,\eta}^{\dagger}\bigg{\\}}\left|g_{1}g_{2},\\{0\\}\right\rangle\,,$
(2)

where $\left|g_{1}g_{2},\\{0\\}\right\rangle$ is the ground state of the
system. The complex coefficients $a(t)$, $b_{\alpha\eta}^{(r)}(t)$, and
$c_{\alpha\eta,\beta\eta^{\prime}}(t)$ correspond to the probability
amplitudes of having an excitation in both the atoms, $r^{\mathrm{th}}$ atom
and field mode $\left\\{\alpha,\eta\right\\}$, and field modes
$\left\\{\alpha,\eta\right\\}$ and $\left\\{\beta,\eta^{\prime}\right\\}$,
respectively. The coefficient $c_{\alpha\eta}(t)$ represents the probability
amplitude of exciting two photons in mode $\left\\{\alpha,\eta\right\\}$.

In this work, our emphasis will be on the case of perfectly coupled atoms;
under this condition the modes $\alpha$ are the guided modes of the waveguide
and $\gamma$ is the decay rate of the emitters into those modes. Defining
$B_{\alpha\eta}^{(r)}=b_{\alpha\eta}^{(r)}(t)e^{-i\Delta_{\alpha}t}$ and
$C_{\alpha\eta}=c_{\alpha\eta}(t)e^{-2i\Delta_{\alpha}t}$, and formally
integrating the equation for $c_{\alpha\eta,\beta\eta^{\prime}}(t)$, the
Schrödinger equation yields the following system of delay-differential
equations for the excitation amplitudes:

$\displaystyle\dot{a}(t)=$
$\displaystyle-\sum_{\alpha,\eta}\sum_{s=1}^{2}g_{\alpha}B_{\alpha\eta}^{(s)}(t)e^{-ik_{\alpha\eta}z_{s}}\,,$
(3) $\displaystyle\dot{B}_{\alpha\eta}^{(r)}(t)=$
$\displaystyle-\left(i\Delta_{\alpha}+\frac{\gamma}{2}\right)B_{\alpha\eta}^{(r)}(t)-\sqrt{2}g_{\alpha}C_{\alpha\eta}(t)e^{ik_{\alpha\eta}z_{r}}+$
$\displaystyle
g_{\alpha}^{\ast}a(t)e^{ik_{\alpha\eta}z_{r}}-\frac{\gamma}{2}e^{i\phi}e^{-i\Delta_{\alpha}t}B_{\alpha\eta}^{(s)}(t-\tau)\Theta(t-\tau)$
(4) $\displaystyle\dot{C}_{\alpha\eta}(t)=$
$\displaystyle-2i\Delta_{\alpha}C_{\alpha\eta}(t)+\sqrt{2}g_{\alpha}^{\ast}\sum_{s=1}^{2}B_{\alpha\eta}^{(s)}(t)e^{-ik_{\alpha\eta}z_{s}}\,.$
(5)

We present a detailed derivation in the Appendix. We solve the system of
equations via Laplace transforms considering the initial conditions $a(0)=1$
and $B_{\alpha\eta}(0)=C_{\alpha\eta}(0)=0$.

System dynamics.— The Laplace transform of $a(t)$ is (see Appendix)

$\displaystyle\tilde{a}(s)$ $\displaystyle\approx$
$\displaystyle\left\\{s-\sum_{\alpha,\eta}|g_{\alpha}|^{2}\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}\cos(\omega_{\alpha}\tau)}{\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}+s+i\Delta_{\alpha}+\frac{\gamma}{2}\right]\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}-s-i\Delta_{\alpha}-\frac{\gamma}{2}\right]}\right\\}^{-1}\,,$
(6)

where $\phi=\omega_{0}\tau$ is the resonance field propagation phase.

Figure 2:  Probability $P^{(1)}(t)$ of having only one of the atoms excited
for a delay of $\gamma\tau=0.895$ in comparison with that for coincident and
independent atoms ($\gamma\tau=0$ and $\gamma\tau\to\infty$, respectively). We
note that for a propagation phase of $\phi=n\pi$, there is a formation of a
steady atom-photon bound state as indicated in Eq. (Delay-induced spontaneous
dark state generation from two distant excited atoms).

We write the sum over modes on the RHS of Eq. (6) as an integral over
frequencies by introducing the density of modes $\rho(\omega_{\alpha})$.
Additionally, we use the Wigner-Weisskopf approximation and set
$\rho(\omega_{\alpha}),\,g_{\alpha}\approx\rho(\omega_{0}),\,g_{0}$, evaluated
at the atomic resonance frequency. With these considerations we obtain
$\tilde{a}(s)=1/(s+\gamma)$ for all $\tau$. Taking the inverse Laplace
transform we get $a(t)=e^{-\gamma t}$, which gives the time-dependent
probability of having two atomic excitations

$\displaystyle P^{(2)}(t)=\left\lvert a(t)\right\rvert^{2}=e^{-2\gamma t}\,.$
(7)

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

$\displaystyle
P^{(1)}(t)\approx\rho(\omega_{0})\sum_{r=1,2}\sum_{\eta}\int_{0}^{\infty}d\omega_{\alpha}\,\left\lvert
b_{\alpha\eta}^{(r)}(t)\right\rvert^{2},$ (8)

where the time-dependent solutions $b_{\alpha\eta}^{(1,2)}(t)$ are derived in
the Appendix. In the limits of two coincident atoms with $\tau=0$, and two
infinitely distant atoms with $\tau\to\infty$, we obtain the expected
Markovian dynamics Lehmberg (1970); Alvarez-Giron and Barberis-Blostein
(2020):

$P^{(1)}(t)=\begin{cases}2\gamma t\,e^{-2\gamma t}&\text{for $\tau=0$}\\\
2e^{-\gamma t}(1-e^{-\gamma t})&\text{for $\tau\to\infty$}\end{cases}\,.$ (9)

Figure 2 shows $P^{(1)}(t)$ for different values of the delay $\tau$ and the
propagation phase $\phi$. For times $t<\tau$, the atoms decay independently.
At the time $t=\tau$, the field emitted by one atom reaches the other,
modifying their decay dynamics depending on the value of the propagation phase
$\phi$. Their instantaneous decay rate, given by the negative slope of
$P^{(1)}(t)$, could momentarily exceed that of standard superfluorescence
($\tau=0$). The delay condition $\gamma\tau\approx 0.375$ and $\phi=n\pi$
maximizes the instantaneous decay rate. This phenomenon is a signature of
superduperradiance, reported in Ref. Sinha  _et al._ (2020a). 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 $\phi=n\pi$
for any delay $\tau$ between the atoms:

$\displaystyle\left|\psi(t\to\infty)\right\rangle=\mathop{\sum\sum}_{\alpha,\eta\neq\beta,\eta^{\prime}}\frac{c_{\alpha\eta,\beta\eta^{\prime}}(t\to\infty)}{2}\left|g_{1}g_{2},\\{1_{\alpha,\eta}1_{\beta,\eta^{\prime}}\\}\right\rangle$
$\displaystyle+\left|e_{1}g_{2}\right\rangle\sum_{\alpha,\eta}b_{\alpha\eta}^{(1)}\left|1_{\alpha,\eta}\right\rangle+\left|g_{1}e_{2}\right\rangle\sum_{\alpha,\eta}b_{\alpha\eta}^{(2)}\left|1_{\alpha,\eta}\right\rangle$
$\displaystyle\,,$ (10)

where we have removed the explicit time dependence of the steady-state
coefficients. The steady-state amplitude for mode
$\left\\{\alpha,\eta\right\\}$ is

$b_{\alpha\eta}^{\left(\begin{subarray}{c}1\\\ 2\end{subarray}\right)}=\mp\eta
e^{\pm i\eta
n\pi/2}\frac{g_{\alpha}^{*}}{1+\frac{\gamma\tau}{2}}\frac{i\sin\left(\Delta_{\alpha}\tau/2\right)}{\gamma-i\Delta_{\alpha}}\,,$
(11)

and $c_{\alpha\eta,\beta\eta^{\prime}}$ is given in the Appendix. The last two
terms in Eq. (Delay-induced spontaneous dark state generation from two distant
excited atoms) 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 Calajó  _et al._ (2019). Using the Born rule and Eq.
(11), we obtain that its probability is $2P^{(1)}\left(t\to\infty\right)$,
where:

$\displaystyle
P^{(1)}\left(t\to\infty\right)=\frac{\sinh\left(\frac{\gamma\tau}{2}\right)}{(1+\frac{\gamma\tau}{2})^{2}e^{\frac{\gamma\tau}{2}}}\,.$
(12)

We note that the probability of ending up in a BIC state is maximum for
$\gamma\tau\approx 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 approximation to solve for the dynamics of the system. Using Eq.
(Delay-induced spontaneous dark state generation from two distant excited
atoms) and Eq. (11) we obtain the following reduced density matrix for the
atomic subsystem

$\hat{\rho}_{A}^{\pm}(t\to\infty)=P^{(1)}\left|\Psi^{\pm}\right\rangle\left\langle\Psi^{\pm}\right|+\left(1-P^{(1)}\right)\left|g_{1}g_{2}\right\rangle\left\langle
g_{1}g_{2}\right|\,,$ (13)

with $+$ corresponding to the case when $n$ is odd and $-$ when it is even and
$\left|\Psi^{\pm}\right\rangle=(\left|e_{1}g_{2}\right\rangle\pm\left|g_{1}e_{2}\right\rangle)/\sqrt{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 $\phi=2n\pi$
$(\phi=(2n+1)\pi)$, the dark state that appears in the late-time limit
corresponds to the antisymmetric state $\left|\Psi^{-}\right\rangle$
(symmetric state $\left|\Psi^{+}\right\rangle$).

(a)

(b)

Figure 3: (a) Concurrence as a function of time. The atoms do not have initial
correlations $\mathcal{C}(0)=0$, after some time $t_{\mathrm{SBE}}$ there is a
sudden birth of entanglement”. For $\phi=(m+1/2)\pi$ $\mathcal{C}(t)=0$ for
all time. If $\phi=m\pi$ concurrence reaches a stationary value
$\mathcal{C}_{\text{ss}}$, given by Eq. (12). (b) Correlation between atomic
dipoles. The atoms decay independently for times $t<\tau$. At $t=\tau$, we
observe an emergence of quantum correlations between the atoms. When
$\phi=m\pi$, the atoms reach a dark steady state for all delays $\tau$. The
maximum concurrence is $\approx 0.141$ at $\gamma\tau\approx 0.895$.

Delayed interactions create quantum correlations between the atoms. We
quantify it using the concurrence $\mathcal{C}(t)$ of the reduced density
matrix of the atoms
$\hat{\rho}_{A}(t)=Tr_{F}(|\psi(t)\rangle\langle\psi(t)\lvert)$ Wootters
(1998). For $\tau\to\infty$, the concurrence is zero throughout the evolution
since the initially uncorrelated atoms evolve independently. Remarkably, when
$\tau=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
$\tau$, we numerically study the dynamics of concurrence, shown in Fig. 3 (a).
We begin studying its behavior for $\gamma\tau\approx 0.895$ (corresponding to
the maximum value of $P^{(1)}\left(t\to\infty\right)$) and $\gamma\tau\approx
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_{\text{SBE}}$, when there is a sudden birth of entanglement
(SBE) López  _et al._ (2008); Mazzola  _et al._ (2009); Ficek and Tanaś
(2003); Tanaś and Ficek (2004); Ficek and Tanaś (2007, 2008); Ashrafi and
Naderi (2013). For $\phi=n\pi$, the concurrence increases until reaching a
stationary value, whereas for $\phi=(n+\frac{1}{2})\pi$, it always remains
zero. In general, for other values of $\phi$, the concurrence suddenly departs
from zero and slowly decays after reaching a maximum value. For
$\gamma\tau\approx 0.895$, we obtain the maximum value for the concurrence.
Fig. 3 (b) shows the emergence of atom-atom correlations defined by
$\mathrm{Tr}\left[\hat{\rho}_{A}(t)\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{-}^{(2)}\right]$.
We note that the atom-atom correlations develop as soon as the atoms ‘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 anti-symmetric 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.

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 atom-photon bound state, with sudden birth of entanglement between
the atoms. Such atom-atom entanglement appears as a result of the non-
Markovian time-delayed feedback of the spontaneously radiated EM field acting
on the emitters. Furthermore, 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 circuits Zheng  _et al._ (2013); Monroe  _et
al._ (2014b); Javadi  _et al._ (2015); Tian  _et al._ (2021) and quantum
networks Yan  _et al._ (2015); Coles  _et al._ (2016); Chang  _et al._ (2020);
Zhong  _et al._ (2021) that benefit from preparing and manipulating long-live
dark-states. In this context, there have been several proposals to generate a
steady entangled state; compared to delayed induced subfluorescence, these
schemes necessitate extra degrees of control, such as external driving fields
Paulisch  _et al._ (2016); Zanner  _et al._ (2022); Gonzalez-Tudela  _et al._
(2011); Zheng and Baranger (2013b), initial entanglement Sinha  _et al._
(2020a), or chiral emission in front of a mirror Zhang  _et al._ (2020b); Ask
and Johansson (2022).

It is not only the generation Horodecki  _et al._ (2009) but also the
stabilization of entangled states that are the key to more efficient 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 Plenio and Knight (1998), the environment can
also be proposed as a generator Kraus  _et al._ (2008); Verstraete  _et al._
(2009) and stabilizer of entanglement Tacchino  _et al._ (2018). Our results
highlight a novel way of environment-assisted entanglement generation and
stabilization via non-Markovian time-delayed feedback.

Adding delay to the most straightforward collective system of two two-level
atoms leads to new phenomenology, breaking the Born and Markov approximations,
non-trivially modifying the dynamics, and spontaneously creating steady-state
quantum correlations. 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 Clemens  _et al._ (2003);
Masson and Asenjo-Garcia (2022); Masson  _et al._ (2020). 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.

Our reductionist approach, under ideal conditions, allows for a clear
interpretation of the effects of delay and how it induces novel quantum optics
phenomena. Moreover, the results contribute to the almost 70 years old
discussion on superradiance and superfluorescence, phenomena at the core of
quantum optics with unceasing interest Timothy Noe II  _et al._ (2012); Rainò
_et al._ (2018); Masson  _et al._ (2020); Ferioli  _et al._ (2021); Huang  _et
al._ (2022); Masson and Asenjo-Garcia (2022).

Acknowledgments.— We acknowledge helpful discussions with Saikat Guha. P.S. is
a CIFAR Azrieli Global Scholar in the Quantum Information Science Program.
This work was support by DGAPA-UNAM under grant IG120518 from Mexico, as well
as CONICYT-PAI grant 77190033, and FONDECYT grant N∘ 11200192 from Chile. K.S.
was supported in part by the John Templeton Foundation Award No. 62422.

##  I Appendix

##  II Equations of motion: Derivation

Using Schrödinger equation with the Hamiltonian and the state ansatz given in
the main text, we get the following differential equations for the probability
amplitudes

$\displaystyle\dot{a}(t)$ $\displaystyle=$
$\displaystyle-\sum_{r=1}^{2}\sum_{\alpha,\eta}g_{a}e^{-i\Delta_{\alpha}t}b_{\alpha\eta}^{(r)}(t)e^{-ik_{\alpha\eta}z_{r}}\,,$
(14) $\displaystyle\dot{b}_{\alpha\eta}^{(r)}(t)$ $\displaystyle=$
$\displaystyle
g_{\alpha}^{*}e^{i\Delta_{\alpha}t}a(t)e^{ik_{\alpha\eta}z_{r}}-\sum_{\beta,\eta^{\prime}}g_{\beta}e^{-i\Delta_{\beta}t}c_{\alpha\eta,\beta\eta^{\prime}}(t)e^{ik_{\beta\eta^{\prime}}z_{r}}-\sqrt{2}g_{\alpha}e^{-i\Delta_{\alpha}t}c_{\alpha\eta}(t)e^{ik_{\alpha\eta}z_{r}}\,,$
(15) $\displaystyle\dot{c}_{\alpha\eta,\beta\eta^{\prime}}(t)$
$\displaystyle=$ $\displaystyle
2g_{\alpha}^{*}e^{i\Delta_{\alpha}t}\sum_{r=1}^{2}b_{\beta\eta^{\prime}}^{(r)}(t)e^{-ik_{\alpha\eta}z_{r}}+2g_{\beta}^{*}e^{i\Delta_{\beta}t}\sum_{r=1}^{2}b_{\alpha\eta}^{(r)}(t)e^{-ik_{\beta\eta^{\prime}}z_{r}}\,,$
(16) $\displaystyle\dot{c}_{\alpha\eta}(t)$ $\displaystyle=$
$\displaystyle\sqrt{2}g_{\alpha}^{*}e^{i\Delta_{\alpha}t}\sum_{r=1}^{2}b_{\alpha\eta}^{(r)}(t)e^{-ik_{\alpha\eta}z_{r}}\,,$
(17)

where we consider atomic positions such that $z_{1}=-z_{2}$.

We integrate (16) with $c_{\alpha\eta,\beta\eta^{\prime}}(0)=0$, and
substitute in (15), obtaining

$\displaystyle\small\dot{b}_{\alpha\eta}^{(r)}(t)$ $\displaystyle=$
$\displaystyle
g_{\alpha}^{*}e^{i\Delta_{\alpha}t}a(t)e^{ik_{\alpha\eta}z_{r}}-2g_{\alpha}^{*}\sum_{s=1}^{2}\int_{0}^{t}dTe^{i\Delta_{\alpha}(t-T)}e^{-ik_{\alpha\eta}z_{s}}\sum_{\beta,\eta^{\prime}}g_{\beta}e^{-i\Delta_{\beta}t}b^{(s)}_{\beta\eta^{\prime}}(t-T)e^{ik_{\beta\eta^{\prime}}z_{r}}$
$\displaystyle-\sqrt{2}g_{\alpha}e^{-i\Delta_{\alpha}t}c_{\alpha\eta}(t)e^{ik_{\alpha\eta}z_{r}}-2\sum_{s=1}^{2}\int_{0}^{t}dT\,b_{\alpha\eta}^{(s)}(t-T)\sum_{\beta,\eta^{\prime}}|g_{\beta}|^{2}e^{-i\Delta_{\beta}T}e^{ik_{\beta\eta^{\prime}}(z_{r}-z_{s})}\,.$

For $t\neq 0$, the second term of the RHS can be approximated as zero assuming
$g_{\beta}$ constant and $b_{\beta\eta^{\prime}}(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 $\rho(w_{\alpha})$
Loudon (1983). For guided modes $k_{\alpha\eta}=\eta\frac{\omega_{\alpha}}{v}$
with $v$ being the phase velocity of the field inside the waveguide, and
$\eta=\pm 1$ labels forward or backward propagation direction of the field
along the waveguide. Thus $\sum_{\alpha,\eta}\to\sum_{\eta=\pm
1}\int_{0}^{\infty}d\omega_{\alpha}\,\rho(\omega_{\alpha})$. 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 $\omega_{\alpha}\approx\omega_{0}$ in the evaluation of
$\rho(\omega_{\alpha})$ and $g_{\alpha}$ functions.

Taking into account all the considerations above we obtain

$\displaystyle\sum_{\beta,\eta^{\prime}}|g_{\beta}|^{2}e^{-i\omega_{\beta}T}e^{ik_{\beta\eta^{\prime}}(z_{r}-z_{s})}$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega_{\beta}\,\rho(\omega_{\beta})|g_{\beta}|^{2}\bigg{\\{}e^{-i\omega_{\beta}(T-\tau_{rs})}+e^{i\omega_{\beta}(T+\tau_{rs})}\bigg{\\}}$
$\displaystyle\approx$
$\displaystyle\rho(\omega_{0})|g_{0}|^{2}\int_{0}^{\infty}d\omega_{\beta}\bigg{\\{}e^{-i\omega_{\beta}(T-\tau_{rs})}+e^{i\omega_{\beta}(T+\tau_{rs})}\bigg{\\}}\,,$

where $\tau_{rs}=\tfrac{z_{r}-z_{s}}{v}$. We make use of Sokhotski-Plemelj
theorem

$\displaystyle\int_{0}^{\infty}d\omega_{\beta}\,e^{-i\omega_{\beta}a}=-i\text{PV}\left(\frac{1}{a}\right)+\pi\delta(a)\,,$
(18)

where PV refers to Cauchy principal value. Absorbing the contribution of the
principal value (which corresponds to the Lamb shifts) into the atomic
transition frequency Gross and Haroche (1982) we obtain

$\displaystyle\sum_{\beta,\eta^{\prime}}|g_{\beta}|^{2}e^{-i\Delta_{\beta}T}e^{ik_{\beta\eta}(z_{r}-z_{s})}$
$\displaystyle=$ $\displaystyle\frac{\gamma
e^{i\omega_{0}T}}{2}\frac{\delta(T-\tau_{rs})+\delta(T+\tau_{rs})}{2}\,,$ (19)

where the single atom decay rate to guided modes is defined as $\gamma\equiv
4\pi\rho(\omega_{0})|g_{0}|^{2}$.

Using (19) we obtain the equation of motion for
$\dot{b}_{\alpha\eta}^{(r)}(t)$

$\displaystyle\small\dot{b}_{\alpha\eta}^{(r)}(t)$ $\displaystyle=$
$\displaystyle
g_{\alpha}^{*}e^{i\Delta_{\alpha}t}a(t)e^{ik_{\alpha\eta}z_{r}}-\sqrt{2}g_{\alpha}e^{-i\Delta_{\alpha}t}c_{\alpha\eta}(t)e^{ik_{\alpha\eta}z_{r}}-\frac{\gamma}{2}b_{\alpha\eta}^{(r)}(t)-\frac{\gamma}{2}e^{i\phi}b_{\alpha\eta}^{(s)}(t-\tau)\Theta(t-\tau)\,,\qquad
r\neq s\,,$

with $\tau=|\tau_{12}|=\frac{|z_{1}-z_{2}|}{v}$, $\phi=\omega_{0}\tau$ and
$\Theta$ is the Heaviside step function.

##  III Equations of motion: Solution

We use the Laplace transform to solve the delayed differential equations of
motion (3)-(5) in the main text. As an intermediate step, we define the
variables

$\displaystyle\tilde{B}_{\alpha\eta}^{(\pm)}(s)=\tilde{B}_{\alpha\eta}^{(1)}(s)e^{-ik_{\alpha\eta}z_{1}}\pm\tilde{B}_{\alpha\eta}^{(2)}(s)e^{-ik_{\alpha\eta}z_{2}}\,.$
(20)

Thus

$\displaystyle s\,\tilde{a}(s)-1$
$\displaystyle=-\sum_{\alpha,\eta}g_{\alpha}\tilde{B}_{\alpha\eta}^{(+)}(s)$
(21)
$\displaystyle\left[s+i\Delta_{\alpha}+\frac{\gamma}{2}+\frac{4|g_{\alpha}|^{2}}{s+2i\Delta_{\alpha}}\right]\tilde{B}_{\alpha\eta}^{(+)}(s)$
$\displaystyle=2g_{\alpha}^{*}\tilde{a}(s)-\frac{\gamma}{2}e^{i\phi}e^{-i\Delta_{\alpha}\tau-s\tau}\left\\{\cos(k_{\alpha\eta}d)\tilde{B}_{\alpha\eta}^{(+)}(s)+i\sin(k_{\alpha\eta}d)\tilde{B}_{\alpha\eta}^{(-)}(s)\right\\}$
(22)
$\displaystyle\left[s+i\Delta_{\alpha}+\frac{\gamma}{2}\right]\tilde{B}_{\alpha\eta}^{(-)}(s)$
$\displaystyle=\frac{\gamma}{2}e^{i\phi}e^{-i\Delta_{\alpha}\tau-s\tau}\left\\{\cos(k_{\alpha\eta}d)\tilde{B}_{\alpha\eta}^{(-)}(s)+i\sin(k_{\alpha\eta}d)\tilde{B}_{\alpha\eta}^{(+)}(s)\right\\}\,,$
(23)

where we use that

$\displaystyle\tilde{C}_{\alpha\eta}(s)=\sqrt{2}g_{\alpha}^{*}\frac{\tilde{B}_{\alpha\eta}^{(+)}(s)}{s+2i\Delta_{\alpha}}\,,$

with $a(0)=1$, $B_{\alpha\eta}^{(1,2)}(0)=0$ and $C_{\alpha\eta}(0)=0$ as
initial conditions.

The solutions of the previous equations are

$\displaystyle\tilde{B}_{\alpha\eta}^{(+)}(s)$
$\displaystyle=-2g_{\alpha}^{*}\tilde{a}(s)\left[s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)\right]\left[\left\\{\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\sin(k_{\alpha\eta}d)\right\\}^{2}\right.$
$\displaystyle\left.+\left\\{s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)\right\\}\left\\{s+i\Delta_{\alpha}+\frac{\gamma}{2}+\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)+\frac{4\lvert
g_{\alpha}\rvert^{2}}{s+2i\Delta_{\alpha}}\right\\}\right]^{-1}$ (24)
$\displaystyle\tilde{B}_{\alpha\eta}^{(-)}(s)$
$\displaystyle=-2g_{\alpha}^{*}\tilde{a}(s)\left[i\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\sin(k_{\alpha\eta}d)\right]\left[\left\\{\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\sin(k_{\alpha\eta}d)\right\\}^{2}\right.$
$\displaystyle\left.+\left\\{s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)\right\\}\left\\{s+i\Delta_{\alpha}+\frac{\gamma}{2}+\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)+\frac{4\lvert
g_{\alpha}\rvert^{2}}{s+2i\Delta_{\alpha}}\right\\}\right]\,.$ (25)

and

$\displaystyle\tilde{a}(s)$
$\displaystyle=\left[s-2\sum_{\alpha,\eta}|g_{\alpha}|^{2}\left[s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)\right]\left[\left\\{\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\sin(k_{\alpha}d)\right\\}^{2}\right.\right.$
$\displaystyle\left.\left.+\left\\{s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-(s+i\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)\right\\}\left\\{s+i\Delta_{\alpha}+\frac{\gamma}{2}+\frac{\gamma}{2}e^{i\phi}e^{-(s+\Delta_{\alpha})\tau}\cos(k_{\alpha\eta}d)+\frac{4\lvert
g_{\alpha}\rvert^{2}}{s+2i\Delta_{\alpha}}\right\\}\right]\right]^{-1}\,.$
(26)

In order to find the inverse Laplace transform of Eqs. (III), (III) and (III)
we use that $|g_{\alpha}|^{2}/\gamma\sim 10^{-4}$ for EM modes in the visible
range Barnes  _et al._ (2020) to neglect the contribution of the term
$4|g_{\alpha}|^{2}/(s+i\Delta_{\alpha})$ in the denominator. This corresponds
to neglecting the probability of exciting two equal $\alpha$ modes traveling
in the same direction compared to the probability of exciting two different
modes. Considering this, we get Eq. (6) in the main text and

$\displaystyle\tilde{B}_{\alpha}^{\left(\begin{subarray}{c}1\\\
2\end{subarray}\right)}(s)$ $\displaystyle\approx$ $\displaystyle-
g_{\alpha}^{*}\tilde{a}(s)e^{\pm
ik_{\alpha\eta}d/2}\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}e^{\mp
ik_{\alpha\eta}d}}{\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}+s+i\Delta_{\alpha}+\frac{\gamma}{2}\right]\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}-s-i\Delta_{\alpha}-\frac{\gamma}{2}\right]}\,.$
(27)

##  IV Solving the integrals for $\tilde{a}(s)$

Eq. (6) with the sum approximated as an integral is

$\displaystyle\tilde{a}(s)$ $\displaystyle\simeq$
$\displaystyle\left\\{s-\frac{\gamma}{\pi}\tilde{a}(s)\int_{-\omega_{0}}^{\infty}\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}-\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}\cos(\Delta_{\alpha}\tau+\phi)}{\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}+s+i\Delta_{\alpha}+\frac{\gamma}{2}\right]\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}-s-i\Delta_{\alpha}-\frac{\gamma}{2}\right]}d\Delta_{\alpha}\right\\}^{-1}\,$
(28) $\displaystyle=$
$\displaystyle\left\\{s-\frac{\gamma}{\pi}\tilde{a}(s)I(s,\gamma,\tau,\phi)\right\\}^{-1}\,,$

where we extend the lower limit in the integral to $-\infty$ as an
approximation. We rewrite the denominator inside the integral using the poles
for the $\Delta_{\alpha}$ variable, determined by the characteristic equation

$\displaystyle\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}-\sigma\left(s+i\Delta_{\alpha}+\frac{\gamma}{2}\right)=0\qquad\Longrightarrow\qquad\Delta^{(\sigma)}_{k}=i\left\\{s+\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right\\}\,,$ (29)

with $W_{k}$ denoting the $k$th branch of Lambert W function Corless  _et al._
(1996) and $\sigma=\pm 1$. For simplicity, we introduce
$r=\frac{\gamma\tau}{2}e^{\gamma\tau/2+i\phi}$. Using partial fraction
decomposition Rao and Ahmed (1968) we obtain

$\displaystyle\frac{1}{\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}+s+i\Delta_{\alpha}+\frac{\gamma}{2}\right]\left[\frac{\gamma}{2}e^{i\phi}e^{-s\tau}e^{-i\Delta_{\alpha}\tau}-s-i\Delta_{\alpha}-\frac{\gamma}{2}\right]}=\frac{\tau}{2}\sum_{k=-\infty}^{\infty}\sum_{\sigma=\pm
1}\frac{1}{W_{k}(\sigma r)+W_{k}^{2}(\sigma
r)}\frac{i}{\Delta_{\alpha}-\Delta_{k}^{(\sigma)}}\,.$ (30)

Using the Cauchy’s integral formula we get

$\displaystyle\int_{-\infty}^{\infty}\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}}{\Delta_{\alpha}-\Delta_{k}^{(\sigma)}}d\Delta_{\alpha}=2\pi
i\frac{W_{k}{(\sigma r)}}{\tau}\,,$ (31)
$\displaystyle\int_{-\infty}^{\infty}\frac{e^{-i\Delta_{\alpha}\tau}\cos(\Delta_{\alpha}\tau+\phi)}{\Delta_{\alpha}-\Delta_{k}^{(\sigma)}}d\Delta_{\alpha}=2\pi
ie^{i\phi}\,.$ (32)

For the second integral we use that
$\textrm{Im}\left(\Delta_{k}^{(\sigma)}\right)>0$ because we can take
$\textrm{Re}(s)$ as large as we need, obtaining

$\displaystyle
I(s,\gamma,\tau,\phi)=-\pi\sum_{k=-\infty}^{\infty}\sum_{\sigma=\pm
1}\frac{1}{1+W_{k}(\sigma r)}\left(1+\frac{\tau e^{i\phi}}{2W_{k}(\sigma
r)}\right)\,.$ (33)

Using the following identities of the Lambert functionsCohen (2020)

$\displaystyle\sum_{k=-\infty}^{\infty}\frac{1}{1+W_{k}(z)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,,$ (34)
$\displaystyle\sum_{k=-\infty}^{\infty}\frac{1}{W_{k}(z)+W_{k}^{2}(z)}$
$\displaystyle=$ $\displaystyle\frac{1}{z}\,,$ (35)

we obtain $I(s,\gamma,\tau,\phi)=-\pi$.

##  V Inverse Laplace transform of $\tilde{B}_{\alpha\eta}^{(1,2)}$

To obtain $B_{\alpha\eta}^{(1,2)}(t)$ we apply inverse Laplace transform of
Eq. (27). First, similar to what is done in Appendix IV, we rewrite the result
using the poles of the denominator but for $s$ variable. Defining

$\displaystyle
s_{k,\alpha}^{(\sigma)}=-\frac{\gamma}{2}-i\Delta_{\alpha}+\frac{1}{\tau}W_{k}(\sigma
r)\,,$ (36)

and using that $k_{\alpha\eta}d=\eta(\Delta_{\alpha}\tau+\phi)$ and
$\tilde{a}(s)=1/(s+\gamma)$, we obtain

$\displaystyle B_{\alpha\eta}^{\left(\begin{subarray}{c}1\\\
2\end{subarray}\right)}(t)=-\frac{g_{\alpha}^{*}\tau}{2}\sum_{\sigma=\pm
1}\sum_{k=\infty}^{\infty}\frac{e^{\pm ik_{\alpha\eta}d/2}}{W_{k}(\sigma
r)+W_{k}^{2}(\sigma
r)}\bigg{[}\mathcal{L}^{-1}\left(\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}}{(s+\gamma)(s-s_{k,\alpha}^{(\sigma)})}\right)-\frac{\gamma}{2}e^{i\phi\mp
k_{\alpha\eta}d}\mathcal{L}^{-1}\left(\frac{e^{-s\tau}}{(s+\gamma)(s-s_{k,\alpha}^{(\sigma)})}\right)\bigg{]}$
(37)

where $\mathcal{L}^{-1}$ denotes inverse Laplace transform. Using that

$\displaystyle\mathcal{L}^{-1}\left(\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}}{(s+\gamma)(s-s_{k,\alpha}^{(\sigma)})}\right)$
$\displaystyle=$
$\displaystyle\frac{\frac{\gamma}{2}-i\Delta_{\alpha}}{s_{k,\alpha}^{(\sigma)}+\gamma}e^{-\gamma
t}+\frac{\frac{\gamma}{2}+i\Delta_{\alpha}+s_{k,\alpha}^{(\sigma)}}{s_{k,\alpha}^{(\sigma)}+\gamma}e^{s_{k,\alpha}^{(\sigma)}t}$
(38)
$\displaystyle\mathcal{L}^{-1}\left(\frac{e^{-s\tau}}{(s+\gamma)(s-s_{k,\alpha}^{(\sigma)})}\right)$
$\displaystyle=$
$\displaystyle\left(e^{s_{k,\alpha}^{(\sigma)}(t-\tau)}-e^{-\gamma(t-\tau)}\right)\frac{\Theta(t-\tau)}{s_{k,\alpha}^{(\sigma)}+\gamma}$
(39)

and taking $b_{\alpha\eta}^{(r)}=B_{\alpha\eta}^{(r)}(t)e^{i\Delta_{\alpha}t}$
we get

$\displaystyle b_{\alpha\eta}^{\left(\begin{subarray}{c}1\\\
2\end{subarray}\right)}(t)$ $\displaystyle=$
$\displaystyle-\frac{g_{\alpha}^{*}\tau}{2}e^{i\Delta_{\alpha}t}e^{\pm
i\eta(\Delta_{\alpha}\tau+\phi)/2}\sum_{\sigma=\pm
1}\sum_{k=-\infty}^{\infty}\frac{1}{W_{k}(\sigma r)+W_{k}^{2}(\sigma
r)}\frac{e^{-\gamma
t}}{\frac{\gamma}{2}-i\Delta_{\alpha}+\frac{1}{\tau}W_{k}(\sigma r)}\times$
$\displaystyle\bigg{[}\frac{\gamma}{2}-i\Delta_{\alpha}+\frac{W_{k}(\sigma
r)}{\tau}e^{\left[\frac{\gamma}{2}-i\Delta_{\alpha}+\frac{1}{\tau}W_{k}(\sigma
r)\right]t}+\frac{\gamma}{2}e^{i(1-\eta)\phi-i(1+\eta)\Delta_{\alpha}\tau+\gamma\tau}\Theta(t-\tau)\bigg{\\{}1-e^{\left[\frac{\gamma}{2}-i\Delta_{\alpha}+\frac{1}{\tau}W_{k}(\sigma
r)\right](t-\tau)}\bigg{\\}}\bigg{]}\,,$

where $W_{k}$ denotes the $k$-branch of Lambert-W function.

##  VI Stationary values

In the late-time limit $\gamma t\to\infty$ we have

$\displaystyle b_{\alpha\eta}^{\left(\begin{subarray}{c}1\\\
2\end{subarray}\right)}(t)\approx-\frac{g_{\alpha}^{*}\tau}{2}\sum_{\sigma=\pm
1}\sum_{k=\infty}^{\infty}\frac{e^{\pm
i\eta(\Delta_{\alpha}\tau+\phi)/2}}{W_{k}(\sigma r)+W_{k}^{2}(\sigma
r)}\frac{e^{-\left[\frac{1}{2}-\frac{1}{\gamma\tau}W_{k}(\sigma
r)\right]\gamma
t}}{\frac{\gamma}{2}-i\Delta_{\alpha}+\frac{1}{\tau}W_{k}(\sigma
r)}\bigg{[}\frac{W_{k}(\sigma
r)}{\tau}-\frac{\gamma}{2}e^{i(1-\eta)\phi}e^{-i\eta\Delta_{\alpha}\tau}e^{-\left[\frac{1}{2}-\frac{1}{\gamma\tau}W_{k}(\sigma
r)\right]\gamma\tau}\bigg{]}\,.$

The above expression does not decay to zero if

$\displaystyle\Re\left[\frac{1}{2}-\frac{1}{\gamma\tau}W_{k}\left(\sigma\frac{\gamma\tau}{2}e^{\gamma\tau/2}e^{i\phi}\right)\right]=0\,.$
(41)

Thus, the terms of the expression that will be non-zero in the long time limit
are $k=0$, and $\\{\sigma=+1,\,\phi=2\pi n\\}$ or
$\\{\sigma=-1,\,\phi=(2n+1)\pi\\}$, with $n\in\mathbb{N}^{0}$, because

$\displaystyle
W_{0}\left(\frac{\gamma\tau}{2}e^{\gamma\tau/2}\right)=\frac{\gamma\tau}{2}\,.$

Taking into account the above, we get a stationary value given by

$\displaystyle b_{\alpha\eta,\,\text{ss}}^{\left(\begin{subarray}{c}1\\\
2\end{subarray}\right)}=-\frac{g_{\alpha}^{*}}{2}\frac{e^{\pm
i\eta(\Delta_{\alpha}\tau+n\pi)/2}}{1+\frac{\gamma\tau}{2}}\frac{1-e^{\mp
i\eta\Delta_{\alpha}\tau}}{\gamma-i\Delta_{\alpha}}\,.$ (42)

Then, we obtain

$\displaystyle P^{(1)}\left(t\to\infty\right)$ $\displaystyle=$
$\displaystyle\frac{\rho(\omega_{0})|g_{0}|^{2}}{2(1+\frac{\gamma\tau}{2})^{2}}\sum_{\eta=\pm
1}\int_{-\infty}^{\infty}d\Delta_{\alpha}\left\lvert\frac{1-e^{-i\eta\Delta_{\alpha}\tau}}{\gamma-i\Delta_{\alpha}}\right\rvert^{2}\,,$

which leads to Eq. (12) in the main text.

For correlations between atomic dipoles, we find the stationary value

$\displaystyle\mathrm{Tr}\left[\hat{\rho}_{A}(t\to\infty)\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{-}^{(2)}\right]$
$\displaystyle=$
$\displaystyle\frac{\rho(\omega_{0})|g_{0}|^{2}\cos(n\pi)}{4(1+\frac{\gamma\tau}{2})^{2}}\sum_{\eta=\pm
1}\int_{-\infty}^{\infty}d\Delta_{\alpha}\frac{\left(1-e^{-i\eta\Delta_{\alpha}\tau}\right)^{2}}{\left\lvert\gamma-i\Delta_{\alpha}\right\rvert^{2}}=-\frac{\cos(n\pi)}{2}\frac{\sinh\left(\frac{\gamma\tau}{2}\right)}{(1+\frac{\gamma\tau}{2})^{2}e^{\frac{\gamma\tau}{2}}}\,.$
(43)

Therefore
$\mathrm{Tr}\left[\hat{\rho}_{A}(t\to\infty)\hat{\sigma}_{+}^{(1)}\hat{\sigma}_{-}^{(2)}\right]=P^{(1)}(t\to\infty)/2$.

##  VII Solution for $c_{\alpha\eta,\beta\eta^{\prime}}(t)$

Applying Laplace transform to equation (16) and defining
$C_{\alpha\eta,\beta\eta^{\prime}}(t)=c_{\alpha\eta,\beta\eta^{\prime}}(t)e^{-i(\Delta_{\alpha}+\Delta_{\beta})t}$
we get

$\displaystyle\tilde{C}_{\alpha\eta,\beta\eta^{\prime}}(s)$ $\displaystyle=$
$\displaystyle
2g_{\alpha}^{*}\sum_{r=1}^{2}\frac{\tilde{B}_{\beta\eta^{\prime}}^{(r)}(s)e^{-ik_{\alpha\eta}z_{r}}}{s+i(\Delta_{\alpha}+\Delta_{\beta})}+2g_{\beta}^{*}\sum_{r=1}^{2}\frac{\tilde{B}_{\alpha\eta}^{(r)}(s)e^{-ik_{\beta\eta^{\prime}}z_{r}}}{s+i(\Delta_{\alpha}+\Delta_{\beta})}\,.$
(44)

To find the function in the time domain, we substitute Eq. (27) and use the
poles (36). Then, applying the inverse Laplace transform we arrive to the
expression

$\displaystyle c_{\alpha\eta,\beta\eta^{\prime}}(t)$ $\displaystyle=$
$\displaystyle 2\tau g_{\alpha}^{*}g_{\beta}^{*}\sum_{\sigma=\pm
1}\sum_{k=-\infty}^{\infty}\frac{e^{i(\Delta_{\alpha}+\Delta_{\beta})t}}{W_{k}(\sigma
r)+W_{k}^{2}(\sigma
r)}\bigg{[}\cos\left[(k_{\alpha\eta}-k_{\beta\eta^{\prime}})d/2\right]\left\\{M_{\alpha,\beta}^{(k,\sigma)}(t)+M_{\beta,\alpha}^{(k,\sigma)}(t)\right\\}$
(45)
$\displaystyle-\frac{\gamma}{2}e^{i\phi}\cos\left[(k_{\alpha\eta}+k_{\beta\eta^{\prime}})d/2\right]\left\\{N_{\alpha,\beta}^{(k,\sigma)}(t)e^{-i\Delta_{\alpha}\tau}+N_{\beta,\alpha}^{(k,\sigma)}(t)e^{-i\Delta_{\beta}\tau}\right\\}\bigg{]}\,,$

where we have defined the functions

$\displaystyle M_{\alpha,\beta}^{(k,\sigma)}(t)$ $\displaystyle=$
$\displaystyle\mathcal{L}^{-1}\left(\frac{s+i\Delta_{\alpha}+\frac{\gamma}{2}}{\left[s+i(\Delta_{\alpha}+\Delta_{\beta}\right]\left[s+\gamma\right]\left[s-s_{k,\alpha}^{(\sigma)}\right]}\right)$
(46) $\displaystyle=$
$\displaystyle\frac{\gamma-2i\Delta_{\alpha}}{\left[\Delta_{\alpha}+i\left(\frac{\gamma}{2}+\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]\left[\Delta_{\alpha}+\Delta_{\beta}+i\gamma\right]}\frac{e^{-\gamma
t}}{2}-\frac{\gamma-2i\Delta_{\beta}}{\left[\Delta_{\beta}+i\left(\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]\left[\Delta_{\alpha}+\Delta_{\beta}+i\gamma\right]}\frac{e^{-i(\Delta_{\alpha}+\Delta_{\beta})t}}{2}$
$\displaystyle+\frac{W_{k}(\sigma
r)}{\tau}\frac{e^{-\left[\frac{1}{2}-\frac{1}{\gamma\tau}W_{k}(\sigma
r)\right]\gamma
t-i\Delta_{\alpha}t}}{\left[\Delta_{\alpha}+i\left(\frac{\gamma}{2}+\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]\left[\Delta_{\beta}+i\left(\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]}\,,$ $\displaystyle N_{\alpha,\beta}^{(k,\sigma)}(t)$
$\displaystyle=$
$\displaystyle\mathcal{L}^{-1}\left(\frac{e^{-s\tau}}{\left[s+i(\Delta_{\alpha}+\Delta_{\beta}\right]\left[s+\gamma\right]\left[s-s_{k,\alpha}^{(\sigma)}\right]}\right)$
(47) $\displaystyle=$
$\displaystyle-\Bigg{(}\frac{e^{-\gamma(t-\tau)}}{\left[\Delta_{\alpha}+i\left(\frac{\gamma}{2}+\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]\left[\Delta_{\alpha}+\Delta_{\beta}+i\gamma\right]}+\frac{e^{-i(\Delta_{\alpha}+\Delta_{\beta})(t-\tau)}}{\left[\Delta_{\beta}+i\left(\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]\left[\Delta_{\alpha}+\Delta_{\beta}+i\gamma\right]}$
$\displaystyle-\frac{e^{-\left[\frac{1}{2}-\frac{1}{\gamma\tau}W_{k}(\sigma
r)\right]\gamma(t-\tau)-i\Delta_{\alpha}(t-\tau)}}{\left[\Delta_{\alpha}+i\left(\frac{\gamma}{2}+\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]\left[\Delta_{\beta}+i\left(\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right)\right]}\Bigg{)}\Theta(t-\tau)\,.$

These functions are not invariant to the interchange of the labels
$\alpha,\beta$.

In the long-time limit $\gamma t\to\infty$ and $\phi=n\pi$ the solution takes
the form

$\displaystyle c_{\alpha\eta,\beta\eta^{\prime}}(t\to\infty)$ $\displaystyle=$
$\displaystyle
2g_{\alpha}^{*}g_{\beta}^{*}\bigg{[}\frac{\cos\left[(k_{\alpha\eta}-k_{\beta\eta^{\prime}})d/2\right]-e^{in\pi}\cos\left[(k_{\alpha\eta}+k_{\beta\eta^{\prime}})d/2\right]}{1+\frac{\gamma\tau}{2}}\left\\{\frac{e^{i\Delta_{\alpha}t}}{\Delta_{\alpha}(\Delta_{\beta}+i\gamma)}+\frac{e^{i\Delta_{\beta}t}}{\Delta_{\beta}(\Delta_{\alpha}+i\gamma)}\right\\}$
$\displaystyle-\frac{\tau}{2\left[\Delta_{\alpha}+\Delta_{\beta}+i\gamma\right]}\sum_{\sigma=\pm
1}\sum_{k=-\infty}^{\infty}\frac{1}{W_{k}(\sigma r)+W_{k}^{2}(\sigma
r)}\left\\{\frac{L_{\alpha\eta,\beta\eta^{\prime}}}{\Delta_{\alpha}+i\left(\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right)}+\frac{L_{\beta\eta^{\prime},\alpha\eta}}{\Delta_{\beta}+i\left(\frac{\gamma}{2}-\frac{1}{\tau}W_{k}(\sigma
r)\right)}\right\\}\bigg{]}\,,$

with

$\displaystyle
L_{\alpha\eta,\beta\eta^{\prime}}=(\gamma-2i\Delta_{\alpha})\cos\left[(k_{\alpha\eta}-k_{\beta\eta^{\prime}})d/2\right]-\gamma
e^{in\pi+i\Delta_{\alpha}\tau}\cos\left[(k_{\alpha\eta}+k_{\beta\eta^{\prime}})d/2\right]\,.$
(49)

## References

  * Kimble (2008) H. J. Kimble, Nature 453, 1023 (2008).
  * Duan and Monroe (2010) L.-M. Duan and C. Monroe, Rev. Mod. Phys. 82, 1209 (2010).
  * Monroe  _et al._ (2014a) C. Monroe, R. Raussendorf, A. Ruthven, K. R. Brown, P. Maunz, L.-M. Duan, and J. Kim, Phys. Rev. A 89, 022317 (2014a).
  * Goban  _et al._ (2015) A. Goban, C.-L. Hung, J. D. Hood, S.-P. Yu, J. A. Muniz, O. Painter, and H. J. Kimble, Phys. Rev. Lett. 115, 063601 (2015).
  * Solano  _et al._ (2017a) P. Solano, P. Barberis-Blostein, F. K. Fatemi, L. A. Orozco, and S. L. Rolston, Nat. Commun. 8, 1857 (2017a).
  * Solano  _et al._ (2017b) P. Solano, J. A. Grover, J. E. Hoffman, S. Ravets, F. K. Fatemi, L. A. Orozco, S. L. Rolston, and A. Optical Nanofibers:, Adv. At. Mol. Opt. Phys. 66, 439 (2017b).
  * Kim  _et al._ (2018) J.-H. Kim, S. Aghaeimeibodi, C. J. K. Richardson, R. P. Leavitt, and E. Waks, Nano Letters 18, 4734 (2018).
  * Wen  _et al._ (2019) P. Y. Wen, K.-T. Lin, A. F. Kockum, B. Suri, H. Ian, J. C. Chen, S. Y. Mao, C. C. Chiu, P. Delsing, F. Nori, G.-D. Lin, and I.-C. Hoi, Phys. Rev. Lett 123, 233602 (2019).
  * Mirhosseini  _et al._ (2019) M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B. Dieterle, A. J. Keller, A. Asenjo-Garcia, D. E. Chang, and O. Painter, Nature 569, 692 (2019).
  * Sheremet  _et al._ (2021) A. S. Sheremet, M. I. Petrov, I. V. Iorsh, A. V. Poshakinskiy, and A. N. Poddubny, arXiv:2103.06824 (2021).
  * Magnard  _et al._ (2020) P. Magnard, S. Storz, P. Kurpiers, J. Schär, F. Marxer, J. Lütolf, T. Walter, J.-C. Besse, M. Gabureac, K. Reuer, A. Akin, B. Royer, A. Blais, and A. Wallraff, Phys. Rev. Lett. 125, 260502 (2020).
  * Tiranov  _et al._ (2023) A. Tiranov, V. Angelopoulou, C. J. van Diepen, B. Schrinski, O. A. D. Sandberg, Y. Wang, L. Midolo, S. Scholz, A. D. Wieck, A. Ludwig, A. S. Sørensen, and P. Lodahl, Science 379, 389 (2023).
  * Martín-Cano  _et al._ (2011) D. Martín-Cano, A. González-Tudela, L. Martín-Moreno, F. J. García-Vidal, C. Tejedor, and E. Moreno, Phys. Rev. B 84, 235306 (2011).
  * van Loo  _et al._ (2013) A. F. van Loo, A. Fedorov, K. Lalumière, B. C. Sanders, A. Blais, and A. Wallraff, Science 342, 1494 (2013).
  * Pichler  _et al._ (2015) H. Pichler, T. Ramos, A. J. Daley, and P. Zoller, Phys. Rev. A 91, 042116 (2015).
  * Lecamp  _et al._ (2007) G. Lecamp, P. Lalanne, and J. P. Hugonin, Phys. Rev. Lett. 99, 023902 (2007).
  * Astafiev  _et al._ (2010) O. Astafiev, A. M. Zagoskin, A. Abdumalikov Jr, Y. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai, Science 327, 840 (2010).
  * Arcari  _et al._ (2014) M. Arcari, I. Söllner, 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, 093603 (2014).
  * Zang  _et al._ (2016) X. Zang, J. Yang, R. Faggiani, C. Gill, P. G. Petrov, J.-P. Hugonin, K. Vynck, S. Bernon, P. Bouyer, V. Boyer, and P. Lalanne, Phys. Rev. Appl. 5, 024003 (2016).
  * Scarpelli  _et al._ (2019) L. Scarpelli, B. Lang, F. Masia, D. M. Beggs, E. A. Muljarov, A. B. Young, R. Oulton, M. Kamp, S. Höfling, C. Schneider, and W. Langbein, Phys. Rev. B 100, 035311 (2019).
  * Zheng and Baranger (2013a) H. Zheng and H. U. Baranger, Phys. Rev. Lett. 110, 113601 (2013a).
  * Tufarelli  _et al._ (2014) T. Tufarelli, M. S. Kim, and F. Ciccarello, Phys. Rev. A 90, 012113 (2014).
  * Guimond  _et al._ (2016) P.-O. Guimond, A. Roulet, H. N. Le, and V. Scarani, Phys. Rev. A 93, 023808 (2016).
  * Zhang  _et al._ (2020a) B. Zhang, S. You, and M. Lu, Phys. Rev. A 101, 032335 (2020a).
  * Del Ángel  _et al._ (2022) A. Del Ángel, P. Solano, and P. Barberis-Blostein, arXiv:2212.01972 (2022).
  * Solano  _et al._ (2023) P. Solano, P. Barberis-Blostein, and K. Sinha, Phys. Rev. A 107, 023723 (2023).
  * Sinha  _et al._ (2020a) K. Sinha, P. Meystre, E. A. Goldschmidt, F. K. Fatemi, S. L. Rolston, and P. Solano, Phys. Rev. Lett. 124, 043603 (2020a).
  * Sinha  _et al._ (2020b) K. Sinha, A. González-Tudela, Y. Lu, and P. Solano, Phys. Rev. A 102, 043718 (2020b).
  * Dinç and Brańczyk (2019) F. Dinç and A. M. Brańczyk, Phys. Rev. Research 1, 032042(R) (2019).
  * Calajó  _et al._ (2019) G. Calajó, Y.-L. L. Fang, H. U. Baranger, and F. Ciccarello, Phys. Rev. Lett. 122, 073601 (2019).
  * Facchi  _et al._ (2019) P. Facchi, D. Lonigro, S. Pascazio, F. V. Pepe, and D. Pomarico, Phys. Rev. A 100, 023834 (2019).
  * Guo  _et al._ (2020) L. Guo, A. F. Kockum, F. Marquardt, and G. Johansson, Phys. Rev. Research 2, 043014 (2020).
  * Yao and Hughes (2009) P. Yao and S. Hughes, Opt. Express 17, 11505 (2009).
  * Crowder  _et al._ (2020) G. Crowder, H. Carmichael, and S. Hughes, Phys. Rev. A 101, 023807 (2020).
  * Lee  _et al._ (2023) A. Lee, H. S. Han, F. K. Fatemi, S. L. Rolston, and K. Sinha, Phys. Rev. A 107, 013701 (2023).
  * Pichler  _et al._ (2017) H. Pichler, S. Choi, P. Zoller, and M. D. Lukin, Proceedings of the National Academy of Sciences 114, 11362 (2017).
  * Crowder  _et al._ (2023) G. Crowder, L. Ramunno, and S. Hughes, arXiv:2302.08093 (2023).
  * Bonifacio and Lugiato (1975) R. Bonifacio and L. A. Lugiato, Phys. Rev. A 11, 1507 (1975).
  * Glauber and Haake (1978) R. Glauber and F. Haake, Physics Letters A 68, 29 (1978).
  * Vrehen  _et al._ (1980) Q. Vrehen, M. Schuurmans, and D. Polder, Nature 285, 70 (1980).
  * Clemens  _et al._ (2003) J. P. Clemens, L. Horvath, B. C. Sanders, and H. J. Carmichael, Phys. Rev. A 68, 023809 (2003).
  * Dicke (1954) R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Gross and Haroche (1982) M. Gross and S. Haroche, Physics Reports 93, 301 (1982).
  * Lehmberg (1970) R. H. Lehmberg, Phys. Rev. A 2, 889 (1970).
  * Alvarez-Giron and Barberis-Blostein (2020) W. Alvarez-Giron and P. Barberis-Blostein, Journal of Physics A: Mathematical and Theoretical 53, 435301 (2020).
  * Wootters (1998) W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).
  * López  _et al._ (2008) C. E. López, G. Romero, F. Lastra, E. Solano, and J. C. Retamal, Phys. Rev. Lett. 101, 080503 (2008).
  * Mazzola  _et al._ (2009) L. Mazzola, S. Maniscalco, J. Piilo, K.-A. Suominen, and B. M. Garraway, Phys. Rev. A 79, 042302 (2009).
  * Ficek and Tanaś (2003) Z. Ficek and R. Tanaś, Journal of Modern Optics 50, 2765 (2003).
  * Tanaś and Ficek (2004) R. Tanaś and Z. Ficek, Journal of Optics B: Quantum and Semiclassical Optics 6, S90 (2004).
  * Ficek and Tanaś (2007) Z. Ficek and R. Tanaś, in  _Quantum-Atom Optics Downunder_ (Optica Publishing Group, 2007) p. QME11.
  * Ficek and Tanaś (2008) Z. Ficek and R. Tanaś, Phys. Rev. A 77, 054301 (2008).
  * Ashrafi and Naderi (2013) M. Ashrafi and M. Naderi, Journal of Modern Optics 60, 331 (2013).
  * Zheng  _et al._ (2013) H. Zheng, D. J. Gauthier, and H. U. Baranger, Phys. Rev. Lett. 111, 090502 (2013).
  * Monroe  _et al._ (2014b) C. Monroe, R. Raussendorf, A. Ruthven, K. R. Brown, P. Maunz, L.-M. Duan, and J. Kim, Phys. Rev. A 89, 022317 (2014b).
  * Javadi  _et al._ (2015) A. Javadi, I. Söllner, M. Arcari, S. L. Hansen, L. Midolo, S. Mahmoodian, G. Kiršanskė, T. Pregnolato, E. H. Lee, J. D. Song, S. Stobbe, and P. Lodahl, Nature Communications 6, 8655 (2015).
  * Tian  _et al._ (2021) Z. Tian, P. Zhang, and X.-W. Chen, Phys. Rev. Appl. 15, 054043 (2021).
  * Yan  _et al._ (2015) W.-B. Yan, B. Liu, L. Zhou, and H. Fan, Europhysics Letters 111, 64005 (2015).
  * Coles  _et al._ (2016) R. J. Coles, D. M. Price, J. E. Dixon, B. Royall, E. Clarke, P. Kok, M. S. Skolnick, A. M. Fox, and M. N. Makhonin, Nature Communications 7, 11183 (2016).
  * Chang  _et al._ (2020) H.-S. Chang, Y. P. Zhong, A. Bienfait, M.-H. Chou, C. R. Conner, E. Dumur, J. Grebel, G. A. Peairs, R. G. Povey, K. J. Satzinger, and A. N. Cleland, Phys. Rev. Lett. 124, 240502 (2020).
  * Zhong  _et al._ (2021) Y. Zhong, H.-S. Chang, A. Bienfait, É. Dumur, M.-H. Chou, C. R. Conner, J. Grebel, R. G. Povey, H. Yan, D. I. Schuster, and A. N. Cleland, Nature 590, 571 (2021).
  * Paulisch  _et al._ (2016) V. Paulisch, H. J. Kimble, and A. González-Tudela, New Journal of Physics 18, 043041 (2016).
  * Zanner  _et al._ (2022) M. Zanner, T. Orell, C. M. F. Schneider, R. Albert, S. Oleschko, M. L. Juan, M. Silveri, and G. Kirchmair, Nature Physics 18, 538 (2022).
  * Gonzalez-Tudela  _et al._ (2011) A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, Phys. Rev. Lett. 106, 020501 (2011).
  * Zheng and Baranger (2013b) H. Zheng and H. U. Baranger, Phys. Rev. Lett. 110, 113601 (2013b).
  * Zhang  _et al._ (2020b) B. Zhang, S. You, and M. Lu, Phys. Rev. A 101, 032335 (2020b).
  * Ask and Johansson (2022) A. Ask and G. Johansson, Phys. Rev. Lett. 128, 083603 (2022).
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * Plenio and Knight (1998) M. B. Plenio and P. L. Knight, Rev. Mod. Phys. 70, 101 (1998).
  * Kraus  _et al._ (2008) B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nature Physics 5, 633 (2009).
  * Tacchino  _et al._ (2018) F. Tacchino, A. Auffèves, M. F. Santos, and D. Gerace, Phys. Rev. Lett. 120, 063604 (2018).
  * Masson and Asenjo-Garcia (2022) S. J. Masson and A. Asenjo-Garcia, Nature Communications 13, 1 (2022).
  * Masson  _et al._ (2020) S. J. Masson, I. Ferrier-Barbut, L. A. Orozco, A. Browaeys, and A. Asenjo-Garcia, Phys. Rev. Lett. 125, 263601 (2020).
  * Timothy Noe II  _et al._ (2012) G. Timothy Noe II, J.-H. Kim, J. Lee, Y. Wang, A. K. Wójcik, S. A. McGill, D. H. Reitze, A. A. Belyanin, and J. Kono, Nature Physics 8, 219 (2012).
  * Rainò  _et al._ (2018) G. Rainò, M. A. Becker, M. I. Bodnarchuk, R. F. Mahrt, M. V. Kovalenko, and T. Stöferle, Nature 563, 671 (2018).
  * Ferioli  _et al._ (2021) G. Ferioli, A. Glicenstein, F. Robicheaux, R. T. Sutherland, A. Browaeys, and I. Ferrier-Barbut, Phys. Rev. Lett. 127, 243602 (2021).
  * Huang  _et al._ (2022) K. Huang, K. K. Green, L. Huang, H. Hallen, G. Han, and S. F. Lim, Nature Photonics 16, 737 (2022).
  * Loudon (1983) R. Loudon,  _The Quantum Theory of Light_ , 2nd ed. (Clarendon Press, Oxford, 1983).
  * Barnes  _et al._ (2020) W. L. Barnes, S. A. R. Horsley, and W. L. Vos, Journal of Optics 22, 073501 (2020).
  * Corless  _et al._ (1996) R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, Advances in Computational Mathematics 5, 329 (1996).
  * Rao and Ahmed (1968) K. R. Rao and N. Ahmed, IEEE Transactions on Education 11, 152 (1968).
  * Cohen (2020) H. Cohen, arXiv:2012.11698 (2020).