How to describe collective decay of uncoupled modes in the input-output formalism

Naively, the Markov approximation prevents both information and energy back flow into a system. While true for a system as a whole, this can breakdown for individual subsystems. In the absence of a coupling term in a system Hamiltonian, different modes in an open quantum system will still develop strong correlations if they radiate into the same quantum field–the electromagnetic continuum of states Dicke (1954). Collective decay, including both Dicke sub-radiance and superradiance, have been observed experimentally in atomic ensembles via near field interactions Gross and Haroche (1982); Kalachev and Kröll (2006); Kalachev (2007); Goban et al. (2015); Cidrim et al. (2020) and far field interactions van Enk et al. (1999); Solano et al. (2017); Asenjo-Garcia et al. (2017); Manzoni et al. (2018); Sinha et al. (2020a); Masson et al. (2020a, b) (e.g. a bi-directional transmission channel such as an optical fiber, allowing for a propagating mediating field), and can give rise to non-Markovian dynamics Sinha et al. (2020b). Even without a bi-directional near or far-field interaction initially uncoupled degrees of freedom in a system can spontaneously radiatively couple Cardimona and Stroud (1983); Tscherbul and Brumer (2014), generating coherence between system degrees of freedom Altenmüller (1995).

In this paper, we derive an extension to the input-output formalism Gardiner and Collett (1985) for open quantum systems with multiple modes radiating into the same continuum such that they decay collectively. Building upon early work in quantum network theory Yurke and Denker (1984), input-output theory is a powerful theoretical toolkit where quantum devices are modeled as networks of modes with quantized inputs and outputs corresponding to continuum degrees of freedom Gardiner and Collett (1985) with extensions to fermionic systems Bi Sun and Milburn (1999); Search et al. (2002); Gardiner (2004) and applications in squeezed-mode generation and parametric amplification Collett and Gardiner (1984); Wu et al. (1986), cascaded quantum systems Carmichael (1993); Gardiner (1993), superconducting circuits Wilson et al. (2011); Korotkov (2011), finite-temperature optomechanics Barchielli and Vacchini (2015), single photon detectors Young et al. (2018); Propp and van Enk (2019a), general quantum amplifiers Caves (1982); Propp and van Enk (2019b); Epstein et al. (2021), quantum transduction Randles and van Enk (2022), and a wide range topics in quantum optics (for a thorough review see Ref. Combes et al. (2017))¹¹1The input-output formalism is also in the same spirit as the LSZ reduction, paving the way to natural connections to high-energy physics .. While adding more continua is straightforward (the effect is generally additive Giacosa et al. (2020)), adding additional system degrees of freedom coupled to the same continua requires modification of the standard quantum Langevin equation.

The main result of our paper is the generalized quantum Langevin equation

$\displaystyle\dot{\hat{c}}_{i}(t)$

$\displaystyle=$

$\displaystyle-i\omega_{i}\hat{c}_{i}(t)-i\!\sum\limits_{j}g_{ij}\hat{c}_{j}(t)-\!\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\hat{c}_{j}(t)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(t)$

describing the Heisenberg-picture evolution of a collection of bosonic system mode annihilation operators $\hat{c}_{i}(t)$ with mode resonance frequencies $\omega_{i}$ and decay rates $\gamma_{i}$ to a shared continuum described by a single noise operator $\hat{b}_{\rm in}(t)$ (Fig. 1). For generality, we have also included inter-mode couplings $g_{ij}$ which are assumed w.l.o.g. to be real. As we will see from (I), it makes a great deal of difference whether a set of modes decay locally (each to their own continuum) or globally (to a single, shared continuum). We will also contrast (I) with the standard quantum Langevin equation, where (naively) applying the results of Ref. Gardiner and Collett (1985) yields

$\displaystyle\dot{\hat{c}}_{i}(t)$

$\displaystyle=$

$\displaystyle-i\omega_{i}\hat{c}_{i}(t)-i\!\sum\limits_{j\neq i}g_{ij}\hat{c}_{j}(t)-\frac{\gamma_{i}}{2}\hat{c}_{j}(t)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(t)$

which is missing the correction terms due to spontaneous radiative coupling. For systems where multiple bosonic modes decay to a single continua, (I) and (I) are not equivalent. Examples include multiple modes in optical cavities Suh et al. (2004); Sundaresan et al. (2015); Link et al. (2022), structured continua with multiple pseudomodes Hughes et al. (2018), quantum dots Livache et al. (2019), multimode laser dynamics Vladimirov et al. (2014), and, for fermions, multiple P-orbitals in an atom or molecule Higginbottom et al. (2016). Crucially, we are not studying cascaded quantum systems where one quantum systems drives another in series Carmichael (1993); for our purposes, we are considering multiple modes decaying in parallel.

In addition to descriptions of multiple modes decaying to the same continuum using other formalisms (dressed states Cardimona and Stroud (1983); Tscherbul and Brumer (2014), collision models Ciccarello (2017), and multimode cavity QED Sundaresan et al. (2015); Kollár et al. (2015); McKay et al. (2015)), there have been a few descriptions of systems with spontaneous radiative coupling in the input-output formalism. In studying fermionic systems, Search et al. Search et al. (2002) include the correct modification to the quantum Langevin equation, and (correctly) interpret it as a indirect coupling between modes induced by the shared continuum of states (that is, spontaneous radiative coupling). Additionally, an equation similar to the generalized quantum Langevin equation (I) appears in the appendix of Ref. Yuan et al. (2020) without derivation or discussion. Therefore, the primary aim of this paper is to carefully derive and discuss the form of the generalized quantum Langevin equation and the physics contained therein. We organize the paper as follows. In section II, we rederive the input-output theory for multi-mode systems decaying to a shared global continuum, arriving at a modified quantum Langevin equation (I). We then study uncoupled parallel networks, emphasizing the strong correlations between discrete mode amplitudes. In section III, we further specialize to two two-mode examples: a) two fermionic qubits and b) two bosonic cavity modes (which can be used to make continuous variable qubits Munro et al. (2000)). For both cases, we observe collective decay phenomena including whole system subradiance and subsystem population revival, the latter of which we explore in more detail in section IV. Finally in section V, we propose applications for reservoir engineering in quantum information experiments and for modeling of non-Markovian open quantum systems in the style of the method of pseudomodes Garraway (1997a, b); Dalton et al. (2001); Mazzola et al. (2008); Pleasance and Petruccione (2021), before briefly revisiting general interpretation in the conclusions.

We now carefully consider the case of a quantum network of multiple discrete modes decaying to a single bosonic continuum (Fig. 1), following closely (but not identically) along the lines of Gardiner and Collet’s seminal paper Gardiner and Collett (1985). Crucially, we will allow independent system operators to couple to the same Markovian continuum in parallel. We work in natural units ($\hbar=1$) and in $1$D, noting that the results of $1$D input-output theory hold in $3$D Gardiner (1986); van Enk (2004).

Discrete modes are operators with sharply-peaked spectra in the Heisenberg picture. The Hamiltonian describing a network of arbitrary discrete modes interacting with each other as well as a single bosonic continuum can be divided into three terms $\hat{H}=\hat{H}_{\rm sys}+\hat{H}_{\rm b}+\hat{H}_{\rm int}$. The first term is a system Hamiltonian

$\displaystyle\hat{H}_{\rm sys}$

$\displaystyle=$

$\displaystyle\sum_{i}\omega_{i}\hat{c}^{\dagger}_{i}\hat{c}_{i}+\sum\limits_{\{i,j\}}g_{ij}\left(\hat{c}^{\dagger}_{i}\hat{c}_{j}+\hat{c}^{\dagger}_{j}\hat{c}_{i}\right)$

(3)

where $\omega_{i}$ is the resonance frequency of the $i$th discrete mode, $g_{ij}$ is the coupling between the $i$th and $j$th discrete mode (w.l.o.g. assumed to be real), and $\hat{c}_{i}$ and $\hat{c}^{\dagger}_{i}$ are the lowering and raising operators for the $i$th discrete mode. At this stage, the system discrete mode operators $\hat{c}_{i}$ are arbitrary and may satisfy bosonic, fermionic, or anyonic commutation relations. Similarly, we leave the couplings $g_{ij}$ arbitrary to incorporate a general network geometry, except to define $g_{ii}=0\,\forall i$. (Any self-interaction can be absorbed into the discrete mode energy $\omega_{i}$.)

The second term is a bosonic Hamiltonian

$\displaystyle\hat{H}_{\rm b}$

$\displaystyle=$

$\displaystyle\int d\omega\omega\hat{b}^{\dagger}(\omega)\hat{b}(\omega)$

(4)

where $\hat{b}(\omega)$ and $\hat{b}^{\dagger}(\omega)$ are the annihilation and creation operators for the monochromatic electromagnetic mode with frequency $\omega$ within the bosonic $1$D continuum of modes, satisfying the canonical commutation relation $[\hat{b}(\omega),\hat{b}^{\dagger}(\omega^{\prime})]=\delta(\omega-\omega^{\prime})$.

The third term is the interaction Hamiltonian

$\displaystyle\hat{H}_{\rm int}$

$\displaystyle=$

$\displaystyle i\int d\omega\sum_{i}\kappa_{i}(\omega)\left(\hat{b}^{\dagger}(\omega)\hat{c}_{i}-\hat{c}^{\dagger}_{i}\hat{b}(\omega)\right)$

(5)

where $\kappa_{i}(\omega)$ is a frequency-dependent coupling describing exchange of excitations between the discrete mode network and the continuum at a single location $x=0$ in the rotating wave approximation (which is valid in the weak-coupling limit). If one relaxes the assumption that the discrete modes occupy the same position, there are additional factors $e^{x_{i}\omega_{i}/c}$ in the interaction Hamiltonian but the derivation below remains effectively unchanged.

We now move to the Heisenberg picture. Here the equations of motion describing continuum and discrete mode time-dependent operator evolution are

	$\displaystyle\dot{\hat{b}}(t,\omega)$	$\displaystyle=$	$\displaystyle-i\omega\hat{b}(t,\omega)+\sum_{j}\kappa_{j}(\omega)\hat{c}_{j}(t)$		(6)
	$\displaystyle\dot{\hat{c}}_{i}(t)$	$\displaystyle=$	$\displaystyle-i[\hat{c}_{i}(t),\hat{H}_{\rm sys}]-\int d\omega\kappa_{i}(\omega)[\hat{c}_{i}(t),\hat{c}^{\dagger}_{i}(t)]\hat{b}(t,\omega)$
		$\displaystyle\rightarrow$	$\displaystyle-i\omega\hat{c}_{i}(t)-i\sum_{j}g_{ij}\hat{c}_{j}(t)-\int d\omega\kappa_{i}(\omega)\hat{b}(t,\omega)$

where in the last line we have assumed the discrete mode operator $\hat{c}_{i}(t)$ to be bosonic and satisfy the canonical bosonic commutation relation for discrete modes $[\hat{c}_{i}(t),\hat{c}^{\dagger}_{i}(t)]=1$. This will also describe fermionic (and anyonic) discrete modes in the weak-field limit so that there is (at most) a single excitation in each mode. Detailed study of strong-field fermionic discrete modes are beyond the scope of the present paper, but we give an analogous derivation of collective decay for fermionic modes in the Appendix A.

The first line of (6) has the formal solution

$\displaystyle\hat{b}(t,\omega)$

$\displaystyle=$

$\displaystyle e^{-i\omega(t-t_{0})}\hat{b}(t_{0},\omega)+\int_{t-t_{0}}^{t}dt^{\prime}e^{-i\omega(t-t^{\prime})}\sum_{j}\kappa_{j}(\omega)\hat{c}_{j}(t^{\prime})$

where $\hat{b}(t_{0},\omega)$ describes initial value of the continuum field at the time $t_{0}$. Inserting (II) into the second line of (6) yields an equation of motion for discrete mode evolution with the continuum field integrated out

	$\displaystyle\dot{\hat{c}}_{i}(t)$	$\displaystyle=$	$\displaystyle-i\omega_{i}\hat{c}_{i}(t)-i\sum\limits_{j}g_{ij}\hat{c}_{j}(t)-\int d\omega\kappa_{i}(\omega)e^{-i\omega(t-t_{0})}\hat{b}(t_{0},\omega)$
		$\displaystyle-$	$\displaystyle\int d\omega\kappa_{i}(\omega)\int_{t-t_{0}}^{t}dt^{\prime}e^{-i\omega(t-t^{\prime})}\sum_{j}\kappa_{j}(\omega)\hat{c}_{j}(t^{\prime}).$

We now invoke the first Markov approximation $\kappa_{i}(\omega)=\sqrt{\frac{\gamma_{i}}{2\pi}}$ for each discrete mode with one key modification; in addition to assuming the decay rates $\gamma_{i}$ to be constant over each discrete mode’s frequency $\omega_{i}$, we must also assume that they are also constant over each other discrete mode’s frequency $\omega_{j}$ in order to pull them out of the frequency integral in (II). To be clear, this is a slightly stronger assumption than is made in Ref. Gardiner and Collett (1985). However, it is subsumed by the standard Markov approximation when the energy spacing $\omega_{i}-\omega_{j}$ is of the same order as the decay rates $\gamma_{i}$ (that is, the spectral half-widths of the discrete modes).

We also define an input field

$\displaystyle\hat{b}_{\rm in}(t)$

$\displaystyle=$

$\displaystyle\frac{1}{\sqrt{2\pi}}\int d\omega e^{-i\omega(t-t_{0})}\hat{b}(t_{0},\omega)$

(9)

satisfying the canonical commutation relation $[\hat{b}_{\rm in}(t),\hat{b}_{\rm in}^{\dagger}(t^{\prime})]=\delta(t-t^{\prime})$, and make the usual arguments for

$\displaystyle\int_{t-t_{0}}^{t}dt^{\prime}f(t^{\prime})\delta(t-t^{\prime})$

$\displaystyle=$

$\displaystyle\frac{1}{2}f(t).$

(10)

Now, we can rewrite (II) in a familiar form

$\displaystyle\dot{\hat{c}}_{i}(t)$

$\displaystyle=$

$\displaystyle-i\omega_{i}\hat{c}_{i}(t)-i\!\sum\limits_{j}g_{ij}\hat{c}_{j}(t)-\!\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\hat{c}_{j}(t)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(t)$

where the single decay term in the traditional Langevin equation (I) has been replaced with a sum over all of the modes $j$ (including the mode $i$ itself). This generalized Langevin equation for multiple bosonic modes is generic, and remains unchanged if the decay rates $\gamma_{i}$ and frequencies $\omega_{i}$ are made time-dependent, making the equation amenable to study for experiments involving specific temporal modes of light Titulaer and Glauber (1966); Kiilerich and Mølmer (2020); Propp and van Enk (2020), as well as realistic noisy systems in the lab²²2Note however, that the correlations between modes in 13 do not remain unchanged when $\gamma_{i}$ and $\omega_{i}$ are made time-dependent; there are additional convolution terms that disrupt the exact correlations..

We note that the first, second, and fourth term in(II) are standard in the quantum Langevin treatment: the first term corresponds to the free evolution of the discrete mode, as its phase rotates with frequency $\omega_{i}$, the second term describes the interaction between different discrete modes (recall that we define $g_{ii}=0\,\forall i$), and the fourth term contains the input to the network: a mean-zero random variable, once a partial trace is taken over the environment Gardiner and Zoller (2000). It is the third term in the generalized quantum Langevin equation that is of particular interest here; in addition to the standard damping term $-\frac{\gamma_{i}}{2}\hat{c}_{i}$ one expects in the quantum Langevin equation, the sum contains additional damping terms proportional to other discrete mode operators. Although these terms relate modes operators to one another like a standard coupling term, they are neither contained in the Hamiltonian, nor do they appear in (II) like a physical coupling $g_{ij}$ (they’re missing a $-i$). Indeed, these additional terms are absent in the original work by Gardiner and Collet Gardiner and Collett (1985), and might be surprising for the following reason:

Taking $g_{ij}=0\,\forall i$, the generalized quantum Langevin equation (II) describes an uncoupled network of discrete modes. Although the discrete modes all couple to the same continuum, that continuum is Markovian. As such, one might think that no energy or information can pass between modes. However, the third term in (II) remains in this limit, and encodes quantum correlations between discrete modes (that is, collective decay³³3We note that the collective decay here described is distinct from the collective coupling discussed in Refs. Latune et al. (2019a, b, c, 2020). There, each mode is indistinguishable to the continuum (they have the same energy and decay rates) resulting in a a factorized interaction Hamiltonian and a collective all-to-all coupling. This results in non-trivial finite-temperature thermodynamic effects in addition to non-Markovian evolution. Notably, the collective decay here only induces pair-wise interactions and correlations between modes (13). Whether this also can lead to interesting quantum thermodynamics remains to be seen.) that are inevitable whenever multiple modes couple to the same continuum. To see this, consider the Fourier transform of (II)

$\displaystyle-i\Delta_{i}\hat{c}_{i}(\omega)$

$\displaystyle=$

$\displaystyle-\!\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\hat{c}_{j}(\omega)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(\omega)$

where we have let $g_{ij}\rightarrow 0$ and defined a frequency-dependent detuning $\Delta_{i}=\omega-\omega_{i}$. Since the sum in the right hand side runs over all $j$ (including the mode $i$), the only $i$-dependence on the right hand side is in the common factor $\sqrt{\gamma_{i}}$. Dividing both sides by this leaves the right hand side completely $i$-independent, so that we can conclude

$\displaystyle\frac{\Delta_{i}\hat{c}_{i}(\omega)}{\sqrt{\gamma_{i}}}$

$\displaystyle=$

$\displaystyle\frac{\Delta_{j}\hat{c}_{j}(\omega)}{\sqrt{\gamma_{j}}},\,\forall i,j.$

(13)

Taking the expectation value of (13), we see that the amplitude of a discrete mode is entirely suppressed when the input field is on-resonance with another discrete mode⁴⁴4Indeed, we see in (13) that the spectrum of one mode $\hat{c}_{i}(\omega)$ completely determines every other Heisenberg operator’s spectrum $\hat{c}_{j}(\omega)$ except at the singular frequency, $\omega_{j}$. Since the same is true for their Hermitian conjugates, we can conclude that the electric field amplitudes $\hat{E}_{i}\sim\hat{c}_{i}+\hat{c}_{i}^{\dagger}$ of each mode are also maximally correlated. However, this does not hold for higher order moments of the electric field or mode operators; for example, in (II) the evolution of the number operator has a (non-trivially) mode-dependent driving term $-\sqrt{\gamma_{i}}\left({\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)+\hat{c}^{\dagger}_{i}(t)\hat{b}_{\rm in}(t)}\right)$.: if $\Delta_{i}=0$, then $\left\langle\hat{c}_{j}(\omega)\right\rangle=0\,\forall j\neq i$. This effect becomes more sharply peaked as the decay rates (and thus spectral widths of the discrete modes) decrease. Again, this is surprising! On-resonance is in general not a sufficient condition for perfect absorption by a qubit (or transmission through a quantum network Propp and van Enk (2019a)), yet it is sufficient for the complete suppression of all other discrete mode amplitudes! Furthermore, when two discrete modes become exactly degenerate, one of them decouples from the continuum entirely as discussed in Mazzola et al. (2009). This is a result of the unitarity constraint: if two discrete modes are exactly degenerate, they must couple to different continua with different quantum numbers, and will be discussed in more detail in section IV where we study the subsystem dynamics and memory kernel.

Using the correlations in (13), we note that we can exactly solve for the behavior of each mode individually in the spectral domain

$\displaystyle-i\Delta_{i}\hat{c}_{i}(\omega)$

$\displaystyle=$

$\displaystyle-\left(\!\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\frac{\Delta_{i}}{\Delta_{j}}\right)\hat{c}_{i}(\omega)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(\omega).$

for the case of no direct coupling ($g_{ij}=0\,\forall ij$) (as well as the special case of $g_{ij}\propto\sqrt{\gamma_{i}}\,\forall i$⁵⁵5This corresponds to two modes coupled indirectly through their evanescent fields.). More generally, if there is a direct coupling between $N$ discrete modes $\hat{c}_{i}(t)$, one can always rediagonalize the system into a set of $N$ uncoupled modes $\tilde{\hat{c}}_{i}(t)$, leading again to (II) but with modified decay rates and resonances. Similarly, one can imagine doing the same thing here: treating the extra decay terms as complex couplings and diagonalizing the system to remove them. This is undesirable for our analysis, firstly because the new modes are no longer eigenmodes of the Hamiltonian (unless they are degenerate in energy), and secondly because our analysis leaves the resonance structure and decay rates of the discrete modes individually manipulable, even allowing for complete decoupling from the environment⁶⁶6For examples of the diagonalization process, see Refs. Cardimona and Stroud (1983); Altenmüller (1995); there diagonalizing indeed makes more sense as the coupling to the continuum is always present, so that it is not possible to access the “bare” atomic states. Here, we have ..

Returning to the derived generalized quantum Langevin equation (I), we note the following: though we have solved for the evolution of the mode amplitude operators $\hat{c}_{i}(t)$, this is inconvenient for studying the full evolution of the mode number operators themselves $\hat{N}_{i}(t)=\hat{c}^{\dagger}_{i}(t)\hat{c}_{i}(t)$ ⁷⁷7To get the right open system dynamics for $\hat{N}_{i}(t)$ at finite temperature from $\hat{c}_{i}(t)$ alone, one must use the Stratonovich calculus (which includes the integration endpoints) when evaluating the formal integration so that $\left\langle\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)+\hat{c}^{\dagger}_{i}(t)\hat{b}_{\rm in}(t)\right\rangle$ in (II) contains implicitly the second-order noise term $\gamma_{i}N_{th}$ Gardiner and Zoller (2000). That calculating $\hat{c}^{\dagger}_{i}(t)\hat{c}_{i}(t)$ from the equations of motion for $\hat{c}_{i}(t)$ as in (I) gives the right zero temperature dynamics for $N_{i}(t)$ is a convenient result of unitarity of the combined system-environment dynamics; this would not be true of an analysis using the adjoint representation, for instance., nor the behavior of their expectation values at finite temperature. This is due to the failure of the quantum Langevin equation to give explicitly the correct noise terms for higher-order operators ⁸⁸8Indeed, it is straightforward to see from Eq. 2.12 of Ref. Gardiner and Collett (1985) that the quantum Langevin approach can only ever give the lowest order noise terms explicitly; higher order noise terms are contained implicitly and must be extracted by use of the Stratanovich calculus Gardiner and Zoller (2000).; in general, one must use the quantum regression theorem to derive the equations of motion that the expectation value of $\hat{N}_{i}(t)$ follows Lax (1966); Swain (1981); Gardiner and Zoller (2000); Steck (2018); Blocher and Mølmer (2019). However, here we use a different trick⁹⁹9Note that the trick we use in (15) is in the style of the Ito calculus Gardiner and Zoller (2000), so that we can later separately evaluate system and environment expectation values as usual $\left\langle\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)\right\rangle=\left\langle\hat{b}^{\dagger}_{\rm in}(t)\right\rangle\left\langle\hat{c}_{i}(t)\right\rangle$. to get the correct finite-temperature noise term in the equations of motion for $N_{i}(t)$:

We can rewrite the generalized quantum Langevin equation (II) as an infinitesimal expansion in dt

$\displaystyle\hat{c}_{i}(t+dt)$

$\displaystyle=$

$\displaystyle\hat{c}_{i}(t)+dt\left(-i\omega_{i}\hat{c}_{i}(t)-i\!\sum\limits_{j}g_{ij}\hat{c}_{j}(t)-\!\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\hat{c}_{j}(t)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(t)\right)$

(15)

and expand the infinitesimally evolved number operator

	$\displaystyle\hat{N}_{i}(t+dt)\equiv\hat{c}_{i}^{\dagger}(t+dt)\hat{c}_{i}(t+dt)$	$\displaystyle=$	$\displaystyle\hat{c}_{i}^{\dagger}(t)\hat{c}_{i}(t)+-i\,dt\!\sum\limits_{j\neq i}g_{ij}\left(\hat{c}_{i}^{\dagger}(t)\hat{c}_{j}(t)-\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t)\right)$
		$\displaystyle-$	$\displaystyle\!dt\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\left(\hat{c}_{i}^{\dagger}(t)\hat{c}_{j}(t)+\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t)\right)-dt\sqrt{\gamma_{i}}\left(\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)+\hat{c}_{i}^{\dagger}(t)\hat{b}_{\rm in}(t)\right)$
		$\displaystyle+$	$\displaystyle dt^{2}\gamma_{i}\hat{b}_{\rm in}^{\dagger}(t)\hat{b}_{\rm in}(t)+\mathscr{O}(dt^{2})$

where the reader will notice we have kept a single term of order $dt^{2}$. While this term is second order in the evolution of the quantum operator, it reduces to first order in calculation of the expectation value as, under the Markov approximation, $dt\left\langle\hat{b}_{\rm in}^{\dagger}(t)\hat{b}_{\rm in}(t)\right\rangle=N_{\rm th}$, with $N_{\rm th}$ the average number of excitations at the resonance energy which are incident upon the system from the continuum. From this expression, we can derive Heisenberg equations of motion for the expectation value¹⁰¹⁰10The associated finite-temperature Heisenberg equations of motion for the number operators themselves (as opposed to their expectation values) are intrinsically stochastic and cannot be calculated deterministically. of the number operator

	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\left\langle{\dot{\hat{N}}_{i}(t)}\right\rangle$	$\displaystyle=$	$\displaystyle-\gamma_{i}\left\langle\hat{N}_{i}(t)\right\rangle-i\sum_{j\neq i}g_{ij}\left(\left\langle\hat{c}_{i}^{\dagger}(t)\hat{c}_{j}(t)-\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t)\right\rangle\right)-\sqrt{\gamma_{i}}\sum_{j\neq i}\frac{\sqrt{\gamma_{j}}}{2}\left(\left\langle\hat{c}_{i}^{\dagger}(t)\hat{c}_{j}(t)+\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t)\right\rangle\right)$
		$\displaystyle-$	$\displaystyle\sqrt{\gamma_{i}}\left(\left\langle\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)+\hat{c}^{\dagger}_{i}(t)\hat{b}_{\rm in}(t)\right\rangle\right)+\gamma_{i}N_{th}$
		$\displaystyle\equiv$	$\displaystyle-\gamma_{i}\left\langle\hat{N}_{i}(t)\right\rangle-\sum_{j\neq i}g_{ij}\left\langle\hat{Y}_{ij}(t)\right\rangle-\sum_{j\neq i}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\left\langle\hat{O}_{ij}(t)\right\rangle-\sqrt{\gamma_{i}}\left(\left\langle\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)+\hat{c}^{\dagger}_{i}(t)\hat{b}_{\rm in}(t)\right\rangle\right)+\gamma_{i}N_{th}$
	.				(17)

In the first line of (II), the first term corresponds to decay, the second term correspond to inter-mode coupling (which we are not interested in and will set to zero shortly), the third term corresponds to coherent driving from the continuum (or stochastically, a mean-zero random variable that will disappear upon taking an expectation value), the fourth term describes the correlations the modes develop through their connection to a shared continuum, and the last term is incoherent thermal driving as discussed in the paragraph above. In the second line we define two new Hermitian operators: $\hat{O}_{ij}(t)\equiv\hat{c}_{i}^{\dagger}(t)\hat{c}_{j}(t)+\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t)$ and $\hat{Y}_{ij}(t)\equiv i(\hat{c}_{i}^{\dagger}(t)\hat{c}_{j}(t)-\hat{c}_{j}^{\dagger}(t)\hat{c}_{i}(t))$. These obey their own Heisenberg equations of motion with expectation values

	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\left\langle\dot{\hat{O}}_{ij}(t)\right\rangle=(\omega_{i}-\omega_{j})\left\langle\hat{Y}_{ij}(t)\right\rangle$	$\displaystyle-$	$\displaystyle\sum_{k}\left(g_{ik}\left\langle\hat{Y}_{jk}(t)\right\rangle+g_{jk}\left\langle\hat{Y}_{ik}(t)\right\rangle\right)-\sum_{k}\left(\frac{\sqrt{\gamma_{i}\gamma_{k}}}{2}\left\langle\hat{O}_{jk}(t)\right\rangle+\frac{\sqrt{\gamma_{j}\gamma_{k}}}{2}\left\langle\hat{O}_{ik}(t)\right\rangle\right)$
		$\displaystyle-$	$\displaystyle\sqrt{\gamma_{i}}\left(\left\langle\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{j}(t)+\hat{c}^{\dagger}_{j}(t)\hat{b}_{\rm in}(t)\right\rangle\right)-\sqrt{\gamma_{j}}\left(\left\langle\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)+\hat{c}^{\dagger}_{i}(t)\hat{b}_{\rm in}(t)\right\rangle\right)+2\sqrt{\gamma_{i}\gamma_{j}}N_{th}$
	$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\left\langle\dot{\hat{Y}}_{ij}(t)\right\rangle=(\omega_{i}-\omega_{j})\left\langle\hat{O}_{ij}(t)\right\rangle$	$\displaystyle+$	$\displaystyle\sum_{k}\left(g_{ik}\left\langle\hat{O}_{jk}(t)\right\rangle-g_{jk}\left\langle\hat{O}_{ik}(t)\right\rangle\right)+\sum_{k}\left(\frac{\sqrt{\gamma_{i}\gamma_{k}}}{2}\left\langle\hat{Y}_{jk}(t)\right\rangle-\frac{\sqrt{\gamma_{j}\gamma_{k}}}{2}\left\langle\hat{Y}_{ik}(t)\right\rangle\right)$
		$\displaystyle-$	$\displaystyle i\sqrt{\gamma_{i}}\left({\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{j}(t)-\hat{c}^{\dagger}_{j}(t)\hat{b}_{\rm in}(t)}\right)+i\sqrt{\gamma_{j}}\left({\hat{b}^{\dagger}_{\rm in}(t)\hat{c}_{i}(t)-\hat{c}^{\dagger}_{i}(t)\hat{b}_{\rm in}(t)}\right).$

The second auxiliary equation (II) has the form of a discrete site particle number current density¹¹¹¹11This is clear from how we have defined $\hat{Y}_{ij}\equiv i(\hat{c}_{i}^{\dagger}\hat{c}_{j}-\hat{c}_{j}^{\dagger}\hat{c}_{i})$ in terms of a difference of two terms, the first moving an excitation from the $i$th mode to the $j$th mode and the second doing the inverse., where the coherences between discrete modes (characterized by the first auxiliary equation (II) for $\hat{O}_{ij}(t)$) drive energy exchange between modes. We pause to note that collective decay is a unique feature of open quantum systems; closed systems require direct couplings to enable energy and information exchange, whereas open systems can accomplish this through their decay processes.

In the absence of a coherent input field $\left(\left\langle\hat{b}_{\rm in}(t)\right\rangle=0\right)$ and at zero temperature ($N_{\rm th}=0$), the set of equations (II–II) describing the evolution of the number operator and the two auxiliary operators form a closed system of equations, which can be written as a Gramm matrix. This can be solved using the standard methods of matrix eigenvalue decomposition, yielding (in the simplest case of two discrete modes) an exact solution determined from initial conditions at $t=0$. Note that in the absence of initial coherences $\left\langle\hat{O}_{ij}(0)\right\rangle=\left\langle\hat{Y}_{ij}(0)\right\rangle=0$. For the case study of two non-interacting bosonic modes, we solve the equations of motion for the number operator, auxiliary operators, and discrete mode amplitude operators numerically using the computational package Mathematica, as detailed in Appendix B.

Before moving on to our numerical study of the two-mode system, we perform one more calculation to better understand the dynamics of the single-mode subsystem.

We can make use of the strong correlations between modes to integrate out one mode analogously to as we did for the continuum in our derivation of the quantum Langevin equation in (II). In the absence of direct coupling ($g_{ij}=0$), we perform this operation for an $N$-mode parallel network as follows:

In the quantum Langevin equation in the absence of coupling, the time evolution of a mode $\hat{c}_{i}(t)$ depends on the other $N-1$ modes (and itself) via a term $-\sum_{j}\frac{\sqrt{\gamma_{i}\gamma_{j}}}{2}\hat{c}_{j}(t)$. However, in the frequency domain, we can make use of the strong correlations between the discrete mode amplitude operators in (13), resulting in the spectral solution for the single modes alone (II). Using the method of residues, we can move back to the time domain and arrive at a time-convolution-full (as opposed to time-convolutionless) Langevin equation for a single mode

$\displaystyle\!\dot{\hat{c}}_{i}(t)$

$\displaystyle=$

$\displaystyle-i\omega_{i}\hat{c}_{i}(t)-(\tilde{\gamma}_{i}*\hat{c}_{i})(t)-\sqrt{\gamma_{i}}\hat{b}_{\rm in}(t)$

(20)

where we have defined a time-dependent memory kernel (the generalization of a time-dependent decay rate)

$\displaystyle\!\!\!\!\!\!\!\!\ \tilde{\gamma}_{i}(t)\equiv\frac{\gamma_{i}}{2}2\pi\delta(t)+\sum\limits_{j\neq i}\frac{\gamma_{j}}{2}\left(2\pi\delta(t)+i\pi(\omega_{j}-\omega_{i})e^{i\omega_{j}t}{\rm Sign}[t]\right).$

We note that the memory kernel (III) is highly time non-local, except for when the resonance energies become degenerate. In this case, all that remains is a modification of the decay rate. This effect can be observed in spectral transmittance through a two-sided network comprised of two uncoupled modes in parallel, as studied in Ref. Propp and van Enk (2019a). When the two modes become exactly degenerate, a frequency of maximum transmission and a frequency of perfect reflection annihilate, leaving one state with a modified decay rate (and thus spectral width) and the other discrete mode perfectly decoupled from the continuum as shown in Fig. 2. Indeed, this is necessary to preserve unitary, and is precisely the dark state generation observed and discussed in Ref. Mazzola et al. (2009).

From the highly time-non-local nature of the memory kernel in (III), it is tempting to conclude from this alone the general behavior is non-Markovian; after all, a time-non-local memory kernel for the evolution of quantum states is a hallmark of non-Markovianity Chruściński and Kossakowski (2010). However, this memory kernel is in a Heisenberg equation of motion for a system operator, and not a master equation for density matrix evolution. Modern definitions of non-Markovian dynamics for open quantum systems rely on properties of density matrices (i.e. a monotonically decreasing trace distance, for a thorough review see Ref. Breuer et al. (2016) and are geared towards Schrödinger or interaction picture analyses (although there is already some work on Heisenberg picture analyses, see Ref. de Vega and Alonso (2017); Li et al. (2018)). The task of formally connecting the existing state-centric measures of non-Markovianity to state-independent Heisenberg picture studies is left for future work.

To observe memory effects (including energy back flow Breuer et al. (2009); Guarnieri et al. (2016)) we must look at a specific example, which is the task we turn to now.

We will now focus on the case we are most interested in: two uncoupled ($g_{ij}=0$) discrete modes (e.g. cavity modes, qubits in the weak-field limit) decaying to a continuum of states with no initial excitations or coherences in the external field $\left(\left\langle\hat{b}_{\rm in}(t)\right\rangle=0\right)$ under the first Markov approximation $\left(\frac{d\gamma_{i}}{d\omega}=0\,\forall i\right)$.

First, we study two modes each initialized to the Fock state with a single excitation (for fermions, the excited state). When the two modes are distantly spaced in energy as in Fig. 3, we observe only a slight sub-radiant deviation between the cases where the modes decay globally to a shared continuum, or locally to their own continua. As the two modes approach each other in energy, this subradiance becomes more dramatic and the lifetimes of both modes are extended substantially (Fig. 4). Similarly, the effect on the total occupancy of both modes is subradiant, resulting in non-exponential decay of the full system as the energy spacing decreases (Fig. 5). When the degeneracy becomes exact, one mode decouples from the continuum completely as made clear in (III) and observed experimentally as a dark state in Pellegrino et al. (2014).

We also consider the case where one mode is initialized to a Fock state with a single excitation, and the other mode is initialized in the ground state (Fig. 6). In this case, we observe the same subradiance for the qubit that starts in the excited state. However, the qubit starting in the ground state gains a substantial population before decaying. This is evidence of genuine energy transfer between the two modes, despite both there being no direct coupling between them and their individual decays to the shared global continuum being Markovian.

In the three cases studied so far, the expectation value of the number operator never exceeds unity for either mode. Hence, the numerical results plotted previously will apply to fermionic qubits as well as bosonic modes. Now, we will restrict our study to bosons and study what happens when we initialize each mode in a coherent state.

For an initial coherent state, there is an additional operator of interest: the electric field amplitude. While a measurement of the single-mode electric field operator $\hat{E}_{i}\sim\hat{c}_{i}+\hat{c}_{i}^{\dagger}$ is always zero for Fock states (and thus not of interest), here it is non-zero and exhibits phase-sensitivity: similarly to the photon counting measurements in Figs. 7 and 8, for two in-phase initial coherent states we observe that at short times one mode is sub-radiant while the other is super-radiant with both exhibiting subradiance at longer times (Fig. 9), and for two out-of-phase initial coherent states we observe the faster decaying mode drives the slower decaying modes electric field amplitude (Fig. 10).

First, we notice a change between initial super-and sub-radiant decay depending on the relative phase between the two initial coherent states: when the modes start in phase with each as in Fig. 7, we observe initial subradiance for one mode and superradiance for the other, eventually resulting in population revival in the faster decaying state and subradiance for both modes at longer times; when the modes start out of phase with each other as in Fig. 8, the result is subradiance for both modes as in the Fock state case, but with the faster decaying mode driving the population of the other mode. It is the case where both states are sub-radiant that leads to a population increase in the slower decaying mode. This is easily explained as a result of wave-interference; since the coherent states are initially out of phase, they initially interfere destructively in the continuum and population from the faster decaying mode end up in the slower decaying mode despite their being a) no direct coupling and b) only an mutual driving of the same continuum field—a driving that would be irreversible and incoherent for a single mode alone. Similarly for the in-phase case, initial constructive interference in the continuum results in the faster decaying state exhibiting superradiance (Fig. 9).

Figs. 7–8 can also be understood in terms of particles and non-Markovianity; there is a non-zero probability for an excitation leaking out of a mode to scatter off of an effective structured continuum and return to its original source. This gives rise to energy back-flow, from which we can conclude that the system is evolving in a genuinely non-Markovian fashion Breuer et al. (2009); Guarnieri et al. (2016).

1) Reservoir engineering is a well-known method for modifying the decay and decoherence of a system by “engineering” a desired system-environment interaction Poyatos et al. (1996). By modifying enviornmental parameters in addition to (or instead of) system parameters, one can induce a plethora of system behaviors such as noise squeezing Rabl et al. (2004) (and relatedly, amplification Metelmann and Clerk (2014)), non-reciprocal photo-transmission Metelmann and Clerk (2015), and generation of entanglement between two subsystems Krauter et al. (2011); Muschik et al. (2012). A more pedestrian but nonetheless very important application is the stabilization of quantum states Sarlette et al. (2011). As we have discussed, adding a second auxiliary mode of similar energy will extend the lifetime of a first mode coupled to the same continuum, even if that auxiliary mode is not initially populated (Fig. 6). This is a rudimentary and simple to implement form of reservoir engineering; the presence of a second auxiliary mode will act as an environmental memory to the first, slowing down the dissipation of information and energy into a continuum (or photon loss channel ¹²¹²12For example, it is natural to describe a photo detector as a quantum network of discrete states interacting with both a target input continuum and one or more loss channels, also modeled as continua. Photon loss contributes to both inefficiency Propp and van Enk (2019a) and POVM impurity Propp and van Enk (2020), and can be mitigated with auxiliary mode reservoir engineering.).

2) Conversely, another application of this phenomenon is purely theoretical. The pseudomode representation is a well-studied and useful calculational tool for understanding the dynamics of non-Markovian open quantum systems, especially in the strong-coupling regime Garraway (1997a, b); Dalton et al. (2001); Mazzola et al. (2008); Pleasance and Petruccione (2021). In essence, non-Markovianities in the system-environment are absorbed into coupled, potentially fictitious auxiliary system modes through Fano diagonalization. In addition to coupling to the quantum system in question, each pseudomode couples to a number of local Markovian continua. For the simplest example, consider the case of an atom interacting with a continuum through a simple Fabry-Perot frequency filter. Here, the effect in the spectral domain is a Lorentzian suppression of decay to the continuum Siegman (1986). This filter-induced non-Markovianity is naturally described by an intermediary mode coupling to the atom directly, and then decaying to the flat Markovian continuum (free space).

More generally, the standard pseudomode method relies on a decomposition of the Fourier transformed decay rate or memory kernel into a sum of approximately Lorentzian and Lorentzian-like distributions. The closer to a Lorentzian profile, the more exact the approximation of the environment as an intermediating pseudomode network becomes. As currently formulated, the Fano-diagonalization process that extracts the pseudomodes from the system-environment coupling is only capable of extracting pseudomodes that couple to each other, to the quantum system of interest, and to their own local continua. Allowing pseudomodes to couple to shared global continua as well could more accurately describe environments with spectra from from Lorentzianity. In the absence of an analytic Fano-diagonalization-like process to extract pseudomodes coupled to shared global continua, machine learning may provide a natural method for decomposing open quantum systems into traditional pseudomode networks augmented by collectively decaying pseudomodes ¹³¹³13It bears mentioning that there are intersections between the two applications proposed here; non-Markovianity can be a resource for reservoir engineering as explored in Refs. Poggi et al. (2017); Mirkin et al. (2019), and it is natural to consider the inverse machine learning application: tailoring non-Markovian systems to perform specific reservoir engineering tasks..

Having studied the mathematical of this phenomenon and explored possible applications, we briefly return our attention to the question of interpretation. At first glance, it is somewhat puzzling that we have applied the first Markov approximation whenever possible in the derivation in the theory section, and yet ended up with population revival and memory effects. To this three remarks are warranted: firstly, population revival within a subsystem does not violate Markovianity for the larger system, and the decay of the total population—while non-exponential—is nonetheless monotonic and Markovian (Fig. 5). Secondly, the origin of this effect is the last term of (II): in pulling both $\kappa_{i}(\omega)$ and $\kappa_{j}(\omega)$ outside the integral, we invoke not only the Markov approximation for each mode individually but also for each pair of modes. In this way, making the Markov approximation for each state individually conspires to hide a broader assumption giving rise to this collective decay phenomenon, memory effects for the single-mode subsystems, and non-exponential subradiance for the total system. Thirdly, we would like to reiterate that collective decay is well-known to generate non-Markovianity in a number of systems Fleming and Hu (2012); Breuer et al. (2016); de Vega and Alonso (2017); Sinha et al. (2020c). Furthermore, it is known that multi-mode optical cavities can spontaneously generate coherences Sundaresan et al. (2015) and, for an atom coupled to a multi-mode optical cavity, the subsystem dynamics are non-Markovian Link et al. (2022).

This phenomenon also can be understood as a breakdown of both the immediate irreversibility usually associated with Markovian decay for each mode, and the separability of each mode’s decay into separate processes in $1$D ¹⁴¹⁴14The question of how closely atoms must be spaced to exhibit collective decay was discussed extensively in a recent publication along a similar line to our own (though there the focus is systems with many-body interactions in free space), see Ref. Masson and Asenjo-Garcia (2021),.; the correlations between discrete mode amplitudes prevent the collective decay process from being decomposable into separate decays. Since the modes interact with the same quantum field (the electromagnetic continuum of states) at the same (approximate) physical location, there is simply no time-scale fast enough for energy to dissipate from one mode without reaching the other modes. Thus, the second state acts as environmental memory to the first. To consider a classical analogy, consider the two charged oscillators: classically the oscillation of one will drive the electromagnetic field at that position, which inevitably drives the other as well. For the quantum case, this means there is simply no regime where the Markov approximation can be justified or interpreted in the usual manner. Invoking it mathematically results in the memory effects observed here. That the induced correlated decay disappears far from degeneracy is consistent with this interpretation (the light emitted from one mode is unlikely to drive another far from resonance), as is the observation that this effect disappears for exactly degenerate states (in $1$D, there simply are no other degrees of freedom and so one mode must decouple to preserve unitarity).

In conclusion, we have derived and discussed the features of the generalized quantum Langevin equation (I), amenable to study of systems with multiple modes decaying to the same continuum. For initially uncoupled modes, spontaneous radiative coupling induces pair-wise correlations between discrete mode amplitudes. For the case study of a two-mode system, this dramatically extends the lifetimes of both modes creating non-exponential total population decay, and also induces population revival within the single-mode subsystems. The lifetime extension makes this phenomenon a natural candidate for reservoir engineering, and the simple non-Markovian structure induced suggest natural incorporation of this phenomenon into previously developed methods for modeling non-Markovian open quantum systems.

Collective decay phenomena are ubiquitous in quantum optics; by explicating their study within the input-output formalism, we hope to support their future study.

This work would not have been possible without a number of individuals. This paper is a direct result of a series of productive conversations with Graeme Pleasance; I am very grateful for our chance meeting at the cancelled 2020 APS March Meeting, as well as his expertise on open quantum systems and pseudomodes. I also gratefully acknowledge Steven van Enk, both for his useful comments on this project as well as for first teaching me how to modify input-output theory to incorporate the correlations between modes coupled to the same continuum. I would also like to thank Pablo Solano, Carlton C. Caves, Justin Dressel, Francisco Elohim Becerra, Manuel Munoz, and Philip Blocher for their valuable comments.

This work was supported by National Science Foundation Grant No. PHY-1630114, and is based upon work partially supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator.