Non-Markovian collective emission from macroscopically separated emitters

TL;DR

This paper investigates how two quantum emitters interact with a waveguide in a non-Markovian regime, revealing enhanced or suppressed collective emission due to time-delayed feedback effects.

Contribution

It introduces a model of non-Markovian collective emission from separated emitters, highlighting the role of retarded backaction and feedback in modifying decay rates.

Findings

Enhanced superradiance observed due to feedback

Inhibition of emission beyond Dicke limits

Non-Markovian dynamics significantly alter decay behavior

Abstract

We study the collective radiative decay of a system of two two-level emitters coupled to a one-dimensional waveguide in a regime where their separation is comparable to the coherence length of a spontaneously emitted photon. The electromagnetic field propagating in the cavity-like geometry formed by the emitters exerts a retarded backaction on the system leading to strongly non-Markovian dynamics. The collective spontaneous emission rate of the emitters exhibits an enhancement or inhibition beyond the usual Dicke super- and sub-radiance due to a self-consistent coherent time-delayed feedback.

Tables

Table 1. Table 1: List of experimental parameters achieved up to date in different platforms.

Emitter	${\omega }_0/2\pi$(THz)	${\gamma }_0/2\pi$(MHz)	${\gamma }_{1\text{D}}/\gamma$	$v_g/c$	d (m)	$\eta$
Atoms Kato2019	380	6	0.1	0.7-0.1	2	$\sim 1$
Quantum dots (QD) Kim2018	230	$17\times {10}^3$	0.7	0.08	$15\times {10}^{-6}$	$\sim 0.07$
Superconducting circuits (SC) vanloo2013	$6\times {10}^{-3}$	20	0.9	0.5	$18\times {10}^{-3}$	$\sim 0.02$

Equations

$$H_{\mathrm{int}}=$$

$$|\Psi \left(t\right)⟩=\sum _{m=1}^2{c}_m\left(t\right){\hat{\sigma }}_{+}^{(m)}|g,g,\left\{0\right\}⟩+$$

$${\int }_0^{\infty }\mathrm{d}\omega \left[c_a\left(\omega ,t\right){\hat{a}}^{†}\left(\omega \right)+c_b\left(\omega ,t\right){\hat{b}}^{†}\left(\omega \right)\right]|g,g,\left\{0\right\}⟩,$$

$${\dot{c}}_m\left(t\right)=-\frac{\gamma }{2}\left[c_m\left(t\right)+\beta c_n\left(t-d/v_g\right)\Theta \left(t-d/v_g\right)e^{i{\varphi }_p}\right].$$

$${\tilde{c}}_{\sup }\left(s\right)=$$

$${\tilde{c}}_{\mathrm{sub}}\left(s\right)=$$

$${\tilde{c}}_{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}\left(\tilde{s}\right)$$

$${\gamma }^{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}\approx \frac{1\pm \beta }{1\mp \beta \eta /2}\gamma .$$

$$c_{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}(t)=\frac{1}{\sqrt{2}}\sum _{n\in \mathbb{Z}}{\alpha }_n^{(\pm )}{e}^{-{\gamma }_n^{(\pm )}t/2},$$

$${\gamma }_{\mathrm{inst}}=\gamma \left[1-\frac{W_0\left(-\frac{\eta }{2}\beta e^{\eta /2}\right)}{\eta /2}\right].$$

$${\mathcal{C}}_{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}=$$

$$\left[2\pm \left\{e^{-\eta {\gamma }_n^{(\pm )}/(2\gamma )}+e^{-\eta {{\gamma }_m^{(\pm )}}^{\ast }/(2\gamma )}\right\}.\right]$$

$$c_{\begin{array}{c}\sup \\ \sup \end{array}}(t)$$

$$=\frac{1}{2\pi i}{\int }_{-\infty -iϵ}^{+\infty -iϵ}\frac{\mathrm{d}\tilde{z}}{\sqrt{2}}\frac{e^{i\tilde{z}\gamma t}}{\tilde{z}-\frac{i}{2}\mp \frac{i}{2}\beta e^{-i\eta \tilde{z}}}.$$

$$\tilde{z}-\frac{i}{2}\mp \frac{i}{2}\beta e^{-i\eta \tilde{z}}=0⟹\left(i\eta \bar{z}\right)e^{i\eta \bar{z}}=\mp \frac{\eta }{2}\beta e^{\eta /2}⟹{\tilde{z}}_n^{(\pm )}=\frac{i}{2}\left[1-\frac{W_n\left(\mp \beta \frac{\eta }{2}{e}^{\eta /2}\right)}{\eta /2}\right],$$

$$\frac{1}{\tilde{z}-\frac{i}{2}\mp \frac{i}{2}\beta e^{-i\eta \tilde{z}}}$$

$${\alpha }_n^{(\pm )}=\underset{\tilde{z}\to {\tilde{z}}_n^{(\pm )}}{\lim }\frac{\tilde{z}-{\tilde{z}}_n^{(\pm )}}{\tilde{z}-\frac{i}{2}\mp \frac{i}{2}\beta e^{-i\eta \tilde{z}}}$$

$$\frac{1}{\tilde{z}-\frac{i}{2}\mp \frac{i}{2}\beta e^{-i\eta \tilde{z}}}$$

$$c_{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}(t)=\frac{1}{\sqrt{2}}\sum _{n\in \mathbb{Z}}{\alpha }_n^{(\pm )}{e}^{-{\gamma }_n^{(\pm )}t/2},$$

$${\alpha }_n^{(\pm )}$$

$${\gamma }_n^{(\pm )}$$

$$c_{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}(t)\approx \frac{1}{2\pi i}{\int }_{-i\infty +ϵ}^{+i\infty +ϵ}\frac{\mathrm{d}\tilde{s}}{\sqrt{2}}\frac{e^{\gamma \tilde{s}t}}{\tilde{s}+1/2}\left[\sum _{j=0}^{\infty }\frac{(\pm \beta {)}^j{e}^{-j\eta \tilde{s}}}{{\left(\tilde{s}+1/2\right)}^j}\right].$$

$$c_{\begin{array}{c}\sup \\ \mathrm{sub}\end{array}}(t)\approx \frac{1}{\sqrt{2}}\sum _{j=0}^{\infty }$$

$$c_a\left(\omega ,t\right)$$

$$c_b\left(\omega ,t\right)$$

$$c_a^{\sup }\left(\omega ,t\right)$$

$$c_a^{\mathrm{sub}}\left(\omega ,t\right)$$

$$I_{\sup }(x,t)/I_0=$$

$$\int \mathrm{d}{k}_2\left\{\hat{a}\left(k_2\right)e^{i{k}_{2}x}+\hat{b}\left(k_2\right)e^{-i{k}_{2}x}\right\}{e}^{i{\omega }_{2}t}]|\Psi \left(t\right)⟩$$

$$=$$

$$+e^{i\left(k_1+k_2\right)x}{c}_b^{\ast }\left({\omega }_1,t\right)c_a\left({\omega }_2,t\right)+e^{-i\left(k_1+k_2\right)x}{c}_a^{\ast }\left({\omega }_1,t\right)c_b\left({\omega }_2,t\right)]e^{-i\left({\omega }_1-{\omega }_2\right)t}$$

$$=$$

$$I_{\sup }(x,t)/I_0=$$

$$=$$

$${+\underset{\text{2a}}{\underbrace{\left\{\Theta \left(t-{\tau }_2\right)-\Theta \left(-{\tau }_2\right)\right\}{e}^{-\left(i{\omega }_0+\frac{{\gamma }_n^{(+)}}{2}\right)\left(t-{\tau }_2\right)}}}+\underset{\text{2b}}{\underbrace{\left\{\Theta \left(t+{\tau }_2\right)-\Theta \left({\tau }_2\right)\right\}{e}^{-\left(i{\omega }_0+\frac{{\gamma }_n^{(+)}}{2}\right)\left(t+{\tau }_2\right)}}}]|}^2,$$

$$I_{\mathrm{sub}}(x,t)/I_0=$$

Full text

*Introduction.—*Long-distance interactions are a central tenet of many quantum systems and processes, including large-scale quantum networks, distributed quantum sensing and information processing Kimble08 ; Schoelkopf08 ; Pichler17 ; Pichler16 ; Ge2018 . When the separations between emitters become comparable to the coherence length of the photons mediating their interaction, interference properties of the electromagnetic (EM) field can be modified due to retardation. In such cases, the backaction of the EM field on the emitters leads to a coherent time-delayed feedback on the system dynamics Whalen17 ; Grimsmo15 , thus rendering it non-Markovian BPbook ; Breuer16 ; deVega17 ; Fleming12 .

Non-Markovian open system dynamics can have a variety of physical origins such as structured bath spectral densities, strong system-bath couplings, low temperatures, or initial system-bath correlations among others Breuer16 ; deVega17 ; HPZ92 ; HPZ93 ; Vasile14 ; Groblacher15 . The effects of non-Markovianity have been investigated in collective atomic states in the context of structured reservoirs Thanopulos19 ; Vats98 ; John97 ; AGT2017PRA ; AGT2017PRL and in the strong-coupling regime Daniele19 . Furthermore, delay-induced non-Markovian dynamics has been previously shown in the context of spontaneous emission of single atoms Giessen96 ; Dorner02 ; Tufarelli13 ; Tufarelli14 ; Cook87 ; Carmele13 ; Beige02 ; Guimond16 ; Guimond17 , bound states in continuum (BIC) of the EM field Hsu16 ; Fong17 ; Calajo19 ; Facchi19 ; Dinc18 , and entanglement generation in emitters coupled to waveguides Pichler16 ; Carlos13 .

Additional insights into non-Markovian effects in such a regime can be gained from studying the simple, yet rich quantum optics phenomenon of collective spontaneous emission of two two-level emitters. Cooperative effects in spontaneous emission have an extensive historical background Dicke ; Haroche ; Rehler71 ; Eberly72 , and have been experimentally observed across a range of physical systems Skribanowitz ; Gross76 ; Pavolini85 ; Devoe96 ; Scheibner07 ; Mlynek14 ; Solano2017 ; Kim2018 ; Chen2018 . While the influence of retardation on these effects has been previously studied in Refs. Milonni74 ; Arecchi70 , the non-Markovian dynamics emerging in macroscopically delocalized collective systems is yet unexplored.

In this letter we study the collective radiative dynamics of a pair of macroscopically separated emitters and show that it exhibits non-Markovian features caused by a self-consistent coherent time-delayed feedback. We specifically consider here the emitters are prepared in a super- or sub-radiant electronic state, and present an exact analytical solution of the dynamics of the collective spontaneous emission. We demonstrate that the retarded backaction of the EM field on the emitters can lead to a further enhanced (inhibited) spontaneous emission rate for superradiant (subradiant) states beyond the usual Dicke superradiance (subradiance) Dicke ; Haroche .

We consider two two-level emitters coupled to a waveguide are separated by a distance $d$ comparable to the coherence length ${\sim v_{g}/\gamma}$ of a spontaneously emitted photon, with $v_{g}$ being the group velocity of the field and $\gamma$ the spontaneous emission rate of the individual emitters (see Fig. 1). To gain an intuitive understanding of the non-Markovian nature of this system, consider the following apparent “superradiance paradox:” Assume that the distance $d$ between two emitters prepared in a superradiant state is smaller than the coherence length of an independently emitted photon, but larger than that of a superradiant photon, $v_{g}/\gamma>d>v_{g}/(2\gamma)$. Given that superradiance is an interference effect, one would expect to observe superradiant emission if there is no way to distinguish which atom emitted the field FicekTanas . Now if the emitters radiate collectively, with an emission rate $2\gamma$, then the coherence length of the emitted photons ($v_{g}/(2\gamma)$) is too short to allow for the fields radiated by the two emitters to interfere, suggesting that they should have emitted independently. On the other hand, if we assume that they do emit independently, then the coherence length of the emitted photons ($v_{g}/\gamma$) is long enough that there should be interference and as a result the emitters should emit at the superradiant rate of $2\gamma$ instead. This seeming paradox points to the failure of the Markov approximation: the conventional notion of an exponential decay defining the photon coherence length is no longer valid, and it is necessary instead to consider a full non-Markovian treatment of the system dynamics.

Formal development.— The total Hamiltonian for the emitters+ field system is $H=H_{E}+H_{F}+H_{\mathrm{int}}$, where $H_{E}=\hbar\omega_{0}\sum_{m=1,2}\hat{\sigma}_{+}^{(m)}\hat{\sigma}_{-}^{(m)}$ is the free Hamiltonian for the emitters of resonance frequency $\omega_{0}$, with $\hat{\sigma}_{+}^{(m)}$ is the creation operator of an excitation in the $m^{th}$ emitter; $H_{F}=\hbar\omega\int_{0}^{\infty}\mathrm{d}\omega\left[\hat{a}^{\dagger}\left(\omega\right)\hat{a}\left(\omega\right)+\hat{b}^{\dagger}\left(\omega\right)\hat{b}\left(\omega\right)\right]$ is the free Hamiltonian for the field, with $\hat{a}\left(\omega\right)$ and $\hat{b}\left(\omega\right)$ the annihilation operators for the right- and left-propagating field modes of the waveguide respectively; and $H_{\mathrm{int}}$ is the emitter-field interaction Hamiltonian.

We proceed by making the electric-dipole and rotating-wave approximations (RWA) and expressing the emitters-field interaction Hamiltonian in the interaction picture with respect to the total free Hamiltonian $H_{E}+H_{F}$ as

[TABLE]

where $g(\omega)$ is the atom-field coupling strength footnoteop ; footnoteg . To isolate the non-Markovian behavior arising from the retardation effects from that due to a structured reservoir we assume a flat spectral density of the field modes around the resonance of the emitters such that $g\left(\omega\right)\approx g\left(\omega_{0}\right)$.

Assuming that the total emitters plus field system is initially prepared in the single-excitation manifold, and considering that in the RWA the Hamiltonian preserves the total number of excitations, the state at time $t>0$ is

[TABLE]

where $c_{m}$ and $c_{a,b}(\omega)$ are the excitation amplitudes for the $m^{\mathrm{th}}$ emitter and the guided field modes with frequency $\omega$ respectively, and $\left|g,g,\left\{0\right\}\right\rangle$ is the ground state of the total system, with $\left|\left\{0\right\}\right\rangle$ the field vacuum state. Tracing out the field modes, evolution of emitter excitation amplitudes is given by Asl03

[TABLE]

for $m\neq n$, where $\phi_{p}\equiv k_{0}d=2p\pi$ is the field phase difference upon propagation, which we assume to be an integer multiple of $2\pi$, $\gamma\equiv\gamma_{1D}+\gamma_{3D}$ is the total spontaneous emission rate, and $\gamma_{1D}=\beta\gamma\equiv 4\pi\left|g\left(\omega_{0}\right)\right|^{2}$ is the spontaneous emission rate into the waveguide, with $\beta$ the coupling efficiency of the emitters to the waveguide. The second term in Eq. (3) represents the retarded backaction of the other emitter via the field with a delay $t=d/v_{g}$.

For emitters initially in the super- or sub-radiant states $\left|\Psi_{\begin{subarray}{c}\mathrm{sup}\\ \mathrm{sub}\end{subarray}}\right\rangle\equiv\frac{1}{\sqrt{2}}\left(\left|eg\right\rangle\pm\left|ge\right\rangle\right)$ footnote1 , one can write the Laplace transformed coefficients $\tilde{c}_{m}\left(s\right)\equiv\int_{0}^{\infty}\mathrm{d}t\,e^{-st}c_{m}\left(t\right)$ as

[TABLE]

where $\tilde{s}\equiv s/\gamma$ and $\eta\equiv d\gamma/v_{g}$ is the separation between the emitters normalized by the photon coherence length. Here $\tilde{c}_{\mathrm{sup}}=\tilde{c}^{\mathrm{sup}}_{1}=\tilde{c}^{\mathrm{sup}}_{2}$ and $\tilde{c}_{\mathrm{sub}}=\tilde{c}^{\mathrm{sub}}_{1}=-\tilde{c}^{\mathrm{sub}}_{2}$ are the Laplace space probability amplitudes for the super- and sub-radiant cases respectively footnote2 .

Consider next the case where the emitters are slightly separated, $\eta\ll 1$. Up to linear terms in $\eta$

[TABLE]

which yields an effective spontaneous emission rate

[TABLE]

For a small but finite delay $0<\eta\ll 1$, this can potentially exceed the usual Dicke superradiant emission rate of $2\gamma$ for $\beta=1$. Also, for a subradiant state with an imperfect coupling ($\beta<1$), the effective decay for slightly separated emitters can be slower than that for coincident ones. This surprising enhancement and inhibition of the collective spontaneous emission can be attributed to a stimulated emission as the correlated field emitted from one of the emitters interferes with that from the other Cray82 . The separation dependence of the collective emission rate in addition to the phase difference, demonstrates the influence of retardation on the interference.

We now consider the general case of arbitrarily separated emitters, for which we present an exact analytical solution of the equations of motion (3) based on a well-developed mathematical treatment of delay differential equations (see Corless96 and the Supplemental Material (SM) SM for details). The general expression for the excitation amplitudes of the emitters is

[TABLE]

where $\alpha_{n}^{(\pm)}\equiv\left[1+W_{n}\left(\mp\frac{\eta}{2}\beta e^{\eta/2}\right)\right]^{-1}$ and the effective decay rate $\gamma_{n}^{(\pm)}\equiv\gamma\left[1-2W_{n}\left(\mp\frac{\eta}{2}\beta e^{\eta/2}\right)/\eta\right]$, with $W_{n}(x)$ the $n^{\mathrm{th}}$ branch of the Lambert $W$-function, which is commonly used to describe systems that exhibit time-delayed feedback Corless96 ; Asl03 . We now discuss the consequences of this analytical solution, which is the main result of this work.

*Results.—*Consider first the dynamics of a superradiant initial state. From Eq. (8) and the properties of the Lambert-$W$ function one finds that the superradiant solution has imaginary exponents for $\eta>\eta_{c}$, where we have introduced the normalized critical distance $\eta_{c}\equiv 2W_{0}\left(1/(e\beta)\right)$ Asl03 . Thus for $\eta\geq\eta_{c}$, the atomic excitation amplitudes exhibit oscillations as the atoms decay to their ground state. These can be understood in terms of a field wavepacket bouncing back and forth between the emitters Milonni74 ; SM . For $\beta=1$ this occurs for separations $d>0.56v_{g}/\gamma$, as shown in Fig 2.

For separations $\eta<\eta_{c}$ the emitters radiate independently until a time $\gamma t=\eta$ and collectively afterwards, with an instantaneous decay rate given by

[TABLE]

For a given value of $\beta$, this rate reaches a maximum $\gamma_{\mathrm{inst}}^{\mathrm{max}}$ when the normalized emitter separation equals its critical value $\eta=\eta_{c}$, with ${\gamma_{\mathrm{inst}}^{\mathrm{max}}/\gamma=1-\frac{W_{0}\left(-1/e\right)}{W_{0}\left(1/(e\beta)\right)}}$, as shown in Fig. 3. In the absence of losses and for perfect emitter-waveguide coupling efficiency $(\beta=1)$ the maximum instantaneous spontaneous emission rate is $\gamma_{\mathrm{inst}}^{\mathrm{max}}/\gamma\approx 4.59$, in stark contrast with superradiant emission in Markovian systems.

In the case of a subradiant initial state, for $\beta=1$, the steady state of the dynamics corresponds to a BIC Calajo19 . The probability of reaching the BIC starting initially in the subradiant state of the atoms is given by ${\left|\left\langle\Psi\left(t\rightarrow\infty\right)|\Psi_{\mathrm{BIC}}\right\rangle\right|^{2}=1/\left(1+\eta/2\right)}$ SPIE . The total probability of the emitters being excited in the steady state as $\left|c_{1,2}^{\mathrm{sub}}(\infty)\right|^{2}\rightarrow\frac{1}{2\left({1+\eta/2}\right)^{2}}$, see Fig. 2 (b). We also note that for an initial subradiant state with a delay of $\eta\approx 0.8$, it is possible to achieve a maximal emitter-field steady state entanglement SM .

It is also instructive to explore the cooperative nature of the atom-field dynamics from the perspective of the emitted field intensity $I\left(x,t\right)\propto\left\langle\Psi\left(t\right)|\hat{E}^{\dagger}\left(x,t\right)\hat{E}\left(x,t\right)|\Psi\left(t\right)\right\rangle$, where the electric field operator is $\hat{E}\left(x,t\right)\propto\int_{0}^{\infty}\mathrm{d}k\left[\hat{a}\left(k\right)e^{ikx}+\hat{b}\left(k\right)e^{-ikx}\right]e^{i\omega t}$ and $|\Psi(t)\rangle$ the state of the system (see SM for details). Fig. 4 shows that the fields emitted by the two emitters in the superradiant (subradiant) case interfere constructively (destructively) when the light cones of the two emitters reach each other. Thereafter, depending on their relative phase they produce an interference pattern that is either constructive, leading to a collective ‘superduperradiant’ burst with an instantaneous emission rate greater than $2\gamma$, or destructive, resulting in the a perfect reflection of the field into the optical cavity created by the two atoms. The non-exponential decay of the emitters is an unambiguous signature of the non-Markovian evolution of the system, a result of the self-consistent backaction of the EM field bath, which is accounted for, by a departure from the usual Lindblad dynamics. We further quantify the non-Markovianity of the system in the Supplemental Material [68], which shows that the system is non-Markovian for any value of $\eta$, approaching a Markovian behavior for $\eta\rightarrow 0$.

Noting that in the presence of delay the instantaneous collective decay rate can exceed that of standard Dicke superradiance, one might wonder if the total collective emission into the waveguide also gets enhanced. An important figure of merit to quantify the collective nature of the system in this regard is its cooperativity $\mathcal{C}\equiv\gamma_{\mathrm{in}}/\gamma_{3D}$ PabloReview , such that $\gamma_{\mathrm{in}}=\lim_{t\rightarrow\infty}\int_{0}^{\infty}\mathrm{d}\omega\,\left[\left|c_{a}\left(\omega,t\right)\right|^{2}+\left|c_{b}\left(\omega,t\right)\right|^{2}\right]$ is the fraction of the field emitted into the waveguide and $\gamma_{3D}=\gamma(1-\beta)$ is the fraction of the field that escapes out to the non-guided modes footnote4 . This can be evaluated as SM

[TABLE]

For $\eta>0$, the cooperativity for a superradiant state is reduced compared to that of coincident emitters $(\eta=0)$ as the total collective emission into the guided modes decreases with the emitter separation. In contrast, for an anti-symmetric state we find an enhanced emission into the waveguide as $\eta$ is increased. This is due to the emission of the field into guided modes by the individual emitters until $\gamma t=\eta$, before they start acting collectively (see Fig. 4 (b)). Given that cooperativity is an important figure of merit in quantum information applications, this result illustrates that retardation effects need to be carefully considered in quantum network protocols based on long distance emitters Kurizki15 .

Summary and outlook.— We have shown that the collective radiative decay of two emitters coupled to a one-dimensional waveguide is subject to non-Markovian modifications due to the time-delayed backaction of the electromagnetic field upon the emitters. When prepared in a superradiant initial state they can exhibit time-dependent decay rates that can instantaneously surpass the standard Dicke superradiance rate. The system also allows for long-lived subradiant states characterized by a bound state in the field trapped in the region between the emitters. These effects can be understood as a combination of Dicke super- or subradiance and a retardation of the field wavepacket where the electromagnetic field senses its boundary conditions with a significant delay.

A key parameter for characterizing the dynamics is the emitter separation relative to the photon coherence length $\eta\equiv d\gamma/v_{g}$. It captures the combined physical origin of non-Markovian behavior, as an appreciable value of $\eta$ can be achieved by increasing the emitter separation $d$, but also by increasing the system-environment coupling as in Carlos13 or by exploiting slow group velocities achievable in the presence of a band gap or near a band edge AGT2017PRL . Importantly, as $\eta$ is increased the system dynamics requires keeping track of field correlation functions of increasing order. We note that the non-Markovianity in this case arises explicitly due to retarded backaction effects, despite having a flat spectral density for the bath.

Experimental observations of these effects could be realized across a number of platforms, including quantum dots in photonic waveguides Kim2018 , atoms near optical nano-fibers Solano2017 ; Johnson2019 ; Kato2019 , and superconducting qubits coupled by coplanar waveguides vanloo2013 ; Mirhosseini2018 . Table I in the SM SM summarizes experimental parameters accessible so far. For a system of atoms coupled to nanofibers, values of $\eta\sim 1$ have already been realized Kato2019 . Given the rapid experimental progress in all these platforms, the retarded collective effects studied here can become relevant in a near future.

With the emerging possibility of preparing collective dipoles subject to internal retardation effects and observing their associated complex dynamics in sharp contrast with the more familiar case of emitters confined in sub-wavelength regions, our work adds a new intricacy that has been little explored in the past. Given that the enhancement in the retarded collective decay of two emitters relies on a pairwise time-delayed feedback, it will be interesting to determine the scaling of these effects with the number of emitters. We also note that similar dynamics can arise in a system of linear oscillators JTH2015 , indicating that such collective retarded dissipation should be observable in classical systems as well. It would be then interesting to extend the present dynamics from the single-excitation case considered here to multiple excitations, where one can observe genuinely quantum non-Markovian effects, such as the phenomenon of superfluorescence Schuurmans with retardation, where all the emitters decay collectively from a fully excited state.

Acknowledgements.—P.S. and K.S. are grateful to Luis A. Orozco and William D. Philips for initiating the ideas for this work. K.S. would like to thank Alejandro González-Tudela, Daniel Malz, Darrick E. Chang, and Hyok Sang Han for fruitful discussions, and Guiseppe Calajó for his BIC contribution to the manuscript.

I Solution to the dynamics using Lambert $W$-functions

From taking an inverse Laplace transform of equations (4) and (5) in the main text, we can write the time dependent atomic excitation amplitudes as

[TABLE]

The pole of the denominator is determined by the characteristic equation

[TABLE]

where we have introduced $\bar{z}\equiv\tilde{z}-i/2$ in the intermediate step. Here $W(z)$ is the Lambert $W$-function, or more precisely a set of functions $W_{n}(z)$ comprising the $n$ branches of the inverse relation of the function $f(z)=ze^{z}$, where $z$ is a complex number. In other words $z=f^{-1}\left(ze^{z}\right)\equiv W\left(ze^{z}\right)$, $W(z)=f^{-1}(z)$, and $W_{n}(x)$ is its $n$-th branch Corless96 .

For a given value of $\beta$ there is a critical distance between the atoms defined by $\beta\frac{\eta_{c}}{2}e^{\eta_{c}/2}=1/e$, or $\eta_{c}=2W_{0}\left(1/e\beta\right),$ such that for $\eta>\eta_{c}$ the dominant pole contributions from $\tilde{z}_{0}$ acquire a real component, leading to oscillatory dynamics. We can now write the Laurent series expansion of the denominator of the integrand in (S12)

[TABLE]

such that

[TABLE]

where we have used the property of the $W$-function that $W(z_{0})e^{W(z_{0})}=z_{0}$ 111It can be seen from the definition of the $W$-function that given $W(ze^{z})=z$, if we substitute $y\equiv ze^{z}$ in this definition, one obtains $W(y)=z$. Substituting $W(y)=z$ back in $y\equiv ze^{z}$ yields $y=W(y)e^{W(y)}$.. This gives

[TABLE]

such that we can write the inverse Laplace transform of (S12) as

[TABLE]

where we have defined

[TABLE]

as given in the main text.

We remark that the Lambert $W$-function dependence of the dynamics on the separation is a characteristic of time-delayed systems with self-consistent retarded backaction Corless96 that has been studied in detail in control theory literature and in a variety of other physical problems ranging from instrument design Inst to the AdS/CFT correspondence Ads .

We also note that while the real part of $\gamma^{(\pm)}$ leads to non-Markovian effects in collective spontaneous emission, the imaginary part can similarly lead to delay-induced collective effects in the van der Waals shifts between the emitters Sinha18 ; Fuchs18 .

II Alternative solution in terms of wavepacket oscillations

Consider Eqs. (4) and (5) from the main text, noting that $\left|\left(\tilde{s}+1/2\right)/\left(\beta/2e^{-\eta\tilde{s}}\right)\right|<1$ such that the denominator in Eqs. (4) and (5) in the main text be expanded as the following series

[TABLE]

The dynamics in this limit has been previously studied Milonni74 ; Arecchi70 ; Dorner02 , leading to the atomic excitation amplitudes.

Taking the inverse Laplace transform in the above equation, we obtain

[TABLE]

where the index $j$ physically corresponds to the number of round trips of the photon wavepacket between the two emitters. The multiple reflections of the photon wavepacket modify the time evolution of the spontaneous decay, making it non-exponential, as given in Milonni74 . It is helpful to consider this physical picture and expansion in the large $\eta$ limit, as it allows one to see how the wavepackets bounce between the atoms. However, the solution in terms of $W$-functions allows one to understand the collective dynamics of the emitters more effectively, allowing one to – (1) determine the critical separation $\eta_{c}$ up to which the emitters decay monotonically, as opposed to exhibiting oscillatory dynamics, (2) evaluate the instantaneous spontaneous emission rate for the atomic decay after the onset of collective dynamics, and (3) calculate the probabilities of atomic excitation in the steady state when a BIC state is formed.

III Field Intensity

One can write the dynamics of the field amplitudes from equations (1) and (2) in the main text as follows

[TABLE]

we note that for an initial superradiant (subradiant) state, the fields emitted into the left and right going modes have same (opposite) phases such that

[TABLE]

where we have used Equation (9) from the main text for the emitter dynamics. The superscripts refer to the initial state of the emitters being in a symmetric (sup) or anti-symmetric (sub) superposition.

The intensity of the field emitted by the atoms as a function of position and time can be evaluated as $I\left(x,t\right)=\frac{\epsilon_{0}c}{2}\left\langle\Psi\left(t\right)|\hat{E}^{\dagger}\left(x,t\right)\hat{E}\left(x,t\right)|\Psi\left(t\right)\right\rangle$, where $\hat{E}\left(x,t\right)=\int_{0}^{\infty}\mathrm{d}k\,\mathcal{E}_{k}\left[\hat{a}\left(k\right)e^{ikx}+\hat{b}\left(k\right)e^{-ikx}\right]e^{i\omega t}$ is the electric field operator at position $x$ and time $t$. We assume $\mathcal{E}_{k}\approx\mathcal{E}_{k_{0}}$ to be constant for all $k$. More explicitly, for the superradiant state we obtain

[TABLE]

where we have used the fact that the field amplitudes $c_{a}(\omega,t)=c_{b}(\omega,t)$ for an initial superradiant state, which also leads to a symmetric distribution around $x=0$. The normalization constant $I_{0}=\frac{\epsilon_{0}c\mathcal{E}_{k_{0}}^{2}}{2}$. We thus obtain the intensity as

[TABLE]

where $\tau_{\begin{subarray}{c}1\\ 2\end{subarray}}\equiv x\pm d/2$, the terms $n$a and $n$b correspond to the fields emitted by atom $n=1,2$ into the left and right going modes, a and b, respectively. The terms $n$a and $n$b together make up the light cones emanating from the $n$th emitter. It can be seen that the two light cones interfere constructively with each other going outwards from the atoms and contribute to a beyond superradiant burst, as can be seen from Fig. 4 of the main text.

Similarly for the subradiant state

[TABLE]

We note that as an important difference from the superradiant case, for $n=0$, $\gamma_{n}^{(-)}=0$, which leads the integral $\int_{0}^{\infty}\mathrm{d}k\frac{e^{-i\omega_{0}t}-e^{-i\omega t}}{i\left(\omega-\omega_{0}\right)}$ to diverge due to a pole contribution on the real axis (which corresponds to a divergent self-energy term otherwise). We therefore take only the principal value of the integral to obtain

[TABLE]

Again, similar to the superradiant case, the four terms correspond to fields emitted by the two atoms into the left and right propagating field modes, which in this case interfere destructively with each other outside of the atomic cavity, and form a standing wave inside of the atomic cavity, as illustrated in Fig.4 in the main text.

IV Cooperativity

We remark that the cooperativity for the system can be calculated as $\mathcal{C}=\gamma_{\mathrm{in}}/\gamma_{3D}$ as discussed in the main text. The total emission into the waveguide can be calculated as

[TABLE]

Substituting (S24) in the above this can be simplified to Eq. (10) in the main text.

V Non-Markovianity measure via coherence

We calculate the coherence measure of non-Markovianity as defined in Chanda16 . We first calculate the $l_{1}$-norm of coherence as

[TABLE]

where $\rho(t)$ refers to the system density matrix corresponding to the two emitters given by

[TABLE]

where, for an intial superradiant state

[TABLE]

For an initial subradiant state

[TABLE]

This yields the $l_{1}$-norm of the coherence for the super- and sub-radiant states as

[TABLE]

Further considering a general initial state in the single-excitation subspace

[TABLE]

we find the time evolved state as

[TABLE]

The $l_{1}$-norm of the coherence for the above state is given as

[TABLE]

For a given initial state, one can further define the coherence measure of non-Markovianity as

[TABLE]

In the above expression, we calculate the non-Markovianity measure for a restricted initial state space of the system Chanda16 , considering the initial state space of single emitter excitation. Since the integral is defined over only positive values of $d\mathcal{C}/dt$, a non-zero value of $\mathcal{N}$ immediately implies that the system is non-Markovian.

We can also obtain a non-Markovianity measure for the super- and sub-radiant states as

[TABLE]

We now optimize the coherence measure of non-Markovianity over the initial state parameters to obtain $\mathcal{N}=\mathcal{N}_{\mathrm{sub}}$. This shows that according to the coherence measure of non-Markovianity, an initially subradiant states remains the most non-Markovian throughout the evolution.

VI Correlation between the emitters and the field

We calculate the correlations between the emitters and the EM field via the linear entropy of the reduced system density matrix defined as Bose2000

[TABLE]

Evaluating the above for a super-and sub-radiant state

[TABLE]

We plot the above linear entropy of entanglement in Fig. S3.

VII Possible implementations

Table 1 we lists some potential systems that study the system of two emitters coupled to waveguides experimentally, and the parameter values that have been realized in these platforms.