supplement.pdf

# Collective photon emission in solid state environments: Concatenating non-
markovian and markovian dynamics

Devashish Pandey  Department of Electrical and Photonics Engineering,
Technical University of Denmark, Kgs. Lyngby, Denmark  Martijn Wubs
Department of Electrical and Photonics Engineering, Technical University of
Denmark, Kgs. Lyngby, Denmark  NanoPhoton - Center for Nanophotonics,
Technical University of Denmark, Kgs. Lyngby, Denmark<EMAIL_ADDRESS>
###### Abstract

Collective light emission and multi-qubit dynamics of solid-state quantum
emitters are affected both by their coupling to the light field and to lattice
vibrations. The effect of phonons on quantum emitters is twofold: polaron
formation is described by ultrafast non-markovian dynamics, while slower
dephasing is well described by exponential decay. Both temperature-dependent
processes will affect collective emission and entanglement, but they are
usually not modeled, probably due to a lack of efficient methods especially
for more than two emitters. So here we propose and compare two methods: the
first method concatenates the fast and slow phonon dynamics, and the second is
the polaron method. For a single quantum emitter, we show that the dynamical
equations are identical in both methods, while predictions for two or more
emitters also agree very well. Both of our methods incorporate non-markovian
dynamics due to phonons demonstrating the temperature sensitivity of the
collective photon emission. Utilizing a simplified markovian model instead may
not be accurate enough especially for quantum information applications: for
example, we show how the markovian model may considerably overestimate the
two-emitter concurrence, except at very low temperatures. Our concatenation
and polaron methods can be applied to an arbitrary number and type of quantum
emitters, and beyond the bulk GaAs environment that we consider here.
Especially the concatenation method can take phonon effects into account at
the same computational cost as modelling the emitter-photon interaction alone.
Finally, we present approximate analytical expressions for the collective
emission spectrum for $N$ emitters on a one-dimensional chain.

††journal: opticajournal

##  1 Introduction

With quantum information processing and sensing applications in mind, in
recent years, there has been a surge of research in generating single-photon
emitters. Various strategies have been explored to fabricate arrays of them
[1, 2, 3, 4], which may exhibit collective emission (sub- and superradiance)
with potential applications in photovoltaics [5, 6], energy harvesting
complexes [7, 8], superradiant lasers [9, 10], and quantum memories [11].
Superradiance at very low temperatures has already been demonstrated in solid-
state systems such as quantum dots and nitrogen-vacancy centers in diamond
[12, 13, 14]. Phonon-induced decoherence make the observation of room-
temperature superradiance usually very challenging, so its recent observation
in perovskites is remarkable [15].

Considerable work has already been dedicated to modeling the collective
effects with emitters interacting with only the photon bath [16, 17, 18, 19,
20, 21, 22, 23, 24, 8, 25]. Usually, the markovian approximation is excellent
for the photon bath, as spectral features of the bath do not vary much within
the linewidths of the emitters. In the strong-coupling regime of cavity QED
and in some other carefully engineered situations [26, 27, 28] the markovian
approximation for the photon bath is no longer valid. Such situations are not
considered here.

On the other hand, solid-state emitters also interact with phonon baths, for
which the markovian approximation is typically less well justified. Modeling
collective photon emission in solid-state environments, especially taking into
account the non-markovian phonon bath, is still in its early stages with only
a limited number of works addressing this issue [6, 7, 29, 30, 31]. One of the
main goals of our work is to propose novel and computationally efficient
methods for collective emission in solid-sate environments, to obtain simple
analytical results that give maximum physical insight with few approximations.

Already for a single solid-state emitter, modeling the emission
characteristics can become a complex affair, as it depends on both emitter-
photon and emitter-phonon interactions, but master equations that account for
both types of interactions have already been utilized to explain and stimulate
fascinating experiments [32, 33, 34, 35]. A weak system-bath coupling
justifies the common second-order perturbation approximation, resulting in
computationally convenient master equations of the Redfield or Lindblad types.
These equations become less accurate not only under stronger system-bath
coupling but also at higher temperatures, as they do not capture the non-
markovian dynamics of the quantum emitter due to interactions with phonons
[32, 33, 36].

The non-markovian dephasing dynamics of a quantum emitter in a bulk phonon
bath is typically orders of magnitude faster than its photon emission
dynamics. This can be exploited in calculations of the dynamics to obtain both
physical insight and computational efficiency, for example, in the initial-
slip scheme and its augmented form the concatenation scheme [37, 38, 39, 40,
41]. In the initial-slip scheme, the effect of the ultrafast initial non-
markovian phonon dynamics is accounted for only by modifying the initial state
for a markovian master equation of the quantum system [39, 40]. The
concatenation scheme, on the other hand, concatenates the initial non-
markovian dynamics obtained from the second-order perturbation theory with the
subsequent markovian dynamics, giving a more accurate description of the
system dynamics at all relevant time scales [41].

Here we adopt such a concatenation approach, but we generalize it in three
important ways: first, we replace the second-order perturbation theory with
the exact fast first-order phonon dynamics. Second, we also take the slower
second-order phonon dynamics into account, in the markovian approximation,
following Refs. [42, 43, 44]. Third, we make a generalization from one emitter
to $N$ interacting emitters.

We compare multi-emitter spectra in the resulting concatenation approach with
spectra obtained with a second powerful approach to deal with non-markovian
phonon dynamics, namely the well-known polaron approach [45, 46, 33, 47, 36,
43, 48, 49, 6, 50, 7, 51] that we here also generalize to $N$ emitters. The
polaron master equation is obtained after transforming the Hamiltonian into
the so-called polaron frame, resulting in non-additive interactions between
photon and phonon baths [49, 51]. The polaron approach is the more accurate
yet more complex approach of the two. Therefore, one of our main aims is to
study how well the multi-emitter spectra agree in the two approaches.

In Sec. 2, we describe the interaction of the emitter(s) with the phonon and
photon baths using the concatenation scheme and the polaron method for three
different scenarios corresponding to one, two and more than two emitters.
First, we describe a single emitter and compare the dynamics in our two
methods in Sec. 22.1. In Sec. 22.2, we extend this to two emitters interacting
via the dipole-dipole interaction. Here, we also present the consequences of
neglecting the non-markovian phonon effects in calculations of concurrence, an
important measure of two-qubit entanglement in quantum information processing.
In Sec. 22.3, we generalize the polaron method and our concatenation scheme
for $N$ emitters in a 1D chain. We conclude in Sec. 3. Furthermore, we support
our main results by detailed calculations in the Supplementary Material (SM).

##  2 Emitter(s) interacting with photon and phonon baths

###  2.1 Scenario I : $N=1$

Let us now introduce our model Hamiltonian and some assumptions, first for a
single emitter, before discussing single-emitter and multi-emitter spectra. An
emitter coupled to both the phonon and photon bath can be described by the
Hamiltonian

$H^{\rm 1e}=H_{0}+H_{\rm e-pt}+H_{({\rm e-pn})_{1}}+H_{({\rm
e-pn})_{2}}+H_{\rm pt}+H_{\rm pn},$ (1)

where the superscript ‘1e’ indicates that here we are first considering the
Hamiltonian of a single emitter. Here $H_{\rm
pt}=\sum_{q}\omega_{q}a_{q}^{\dagger}a_{q}$ and $H_{\rm
pn}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$ are the photon and phonon bath
Hamiltonians, respectively and we work in the units of $\hbar=1$.
$H_{0}=\omega_{0}\sigma^{\dagger}\sigma$ is the system Hamiltonian describing
a two-level emitter with resonance frequency $\omega_{0}$, while $\sigma$,
($\sigma^{\dagger}$) is the lowering (raising) atomic operator. $H_{\rm
e-pt}=\sum_{q}(h_{q}a_{q}+h_{q}^{*}a_{q}^{\dagger})\sigma^{\dagger}+\mbox{h.c.}$,
describes the emitter-photon coupling modeled within the dipole approximation,
with $h_{q}$ the coupling between the emitter and the photon field of mode $q$
and $a_{q}$ ($a_{q}^{\dagger}$) is the corresponding photon annihilation
(creation) operator. The mode index $q$ labels both the wavevector and
corresponding polarization indices of transverse plane waves. Likewise, the
$b_{k}$ ($b_{k}^{\dagger}$) are phonon annihilation (creation) operators.

To be specific, we will use the material parameters that describe the emitter-
acoustic phonon coupling in bulk GaAs [52, 53, 44, 50]. The dominant electron-
phonon coupling in quantum dots is deformation-potential coupling to
longitudinal acoustic (LA) phonons, and we will disregard other couplings
[54]. This results in a three-dimensional (3D) phonon spectral density $J_{\rm
pn}(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})=\alpha\omega^{3}\mbox{exp}(-\omega^{2}/\omega_{\rm
c}^{2})$ where $g_{k}$ is the emitter-phonon coupling constant, $\alpha$ is
the deformation potential coupling constant and $\omega_{\rm c}$ is the cutoff
frequency which is given by the ratio of the speed of sound and the size of
the QD [36]. For the light-matter coupling, for simplicity we will also
consider bulk GaAs. We neglect photon loss, which is an excellent
approximation at inter-emitter distance scales of micrometers or below on
which emission can be collective.

Our model Hamiltonian defined in (1) contains two types of emitter-phonon
interactions that both will affect emission spectra. In the presence of
phonons, an excited emitter undergoes a fast non-markovian dephasing on the
picosecond scale, followed by a slower dephasing responsible for the
temperature-dependent broadening of the zero-phonon line (ZPL) [55, 42, 43]
(see supplementary material (SM) Sec. 1). The linear emitter-phonon coupling
$H_{{(\rm
e-pn)}_{1}}=\sigma^{\dagger}\sigma\sum_{k}g_{k}b_{k}^{\dagger}+\mbox{h.c.}$ is
responsible for the fast dephasing. The excitation of the emitter leads to a
sudden change in the charge configuration of the lattice, resulting in a fast
decay of coherence of the emitter until the lattice energy is renormalized to
a new equilibrium. At this point, a polaron is formed in the ionic lattice, at
the characteristic time $\tau=\tau_{\rm P}$ [56]. For GaAs, we estimate
$\tau_{\rm P}\approx 2$ ps (see SM Sec. 1). The slower dephasing is modeled by
a quadratic coupling between the emitter and phonons $H_{{(\rm
e-pn)}_{2}}=\sigma^{\dagger}\sigma\sum_{k,k^{\prime}}f_{k,k^{\prime}}(b_{k}^{\dagger}+b_{k})(b_{k^{\prime}}^{\dagger}+b_{k^{\prime}})$
where $f_{k,k^{\prime}}$ describes phonon-assisted virtual transitions to
higher-lying states of the emitter that cause an inelastic scattering of
phonons from mode $k$ to $k^{\prime}$ [42, 43]. This leads to markovian pure
dephasing of the excited state of the emitter with a dephasing time
$1/\gamma_{\rm pd}$ [42].

Having specified our complete Hamiltonian, we now first calculate single-
emitter dynamics in the concatenation approach, which is expected to be
accurate because the initial polaron formation is much faster than spontaneous
emission and markovian phonon pure dephasing. The initial non-markovian
dynamics is therefore governed by the Hamiltonian $H_{\rm NM}^{\rm
1e}=H_{0}+H_{{(\rm e-pn)_{1}}}+H_{\rm pn}$, which is the exactly solvable
independent boson model (see SM Sec. 2). We then use the exact unitary
operator associated with $H_{\rm NM}^{\rm 1e}$, to calculate the two-time
field correlation [19] using the non-markovian evolution (see [36] and SM Sec.
2),

$g^{(1)}(\tau)=\gamma\langle\sigma^{\dagger}\sigma(\tau)\rangle_{\rm
NME}=\gamma\rho_{\rm
ee}\mathcal{C}(\tau)e^{-i\omega_{0}^{\prime}\tau},\;\;\tau\leq\tau_{\rm P}.$
(2)

Here $\rho_{\rm ee}=\langle e|\rho(0)|e\rangle$ is the population of the
excited state at the initial time $t=0$,
$\mathcal{C}(\tau)=\mathcal{C}_{\infty}(T)\mbox{exp}(\phi(\tau))$ is the
phonon correlation function with
$\mathcal{C}_{\infty}(T)=\mbox{exp}(-\int_{0}^{\infty}d\omega J_{\rm
pn}(\omega)\coth(\beta\omega/2)/\omega^{2})$ its steady-state value, where
$\beta=1/k_{\rm B}T$ is the inverse temperature and $k_{B\rm}$ the Boltzmann
constant. Notice that we have made the temperature dependence of this steady-
state value explicit, which we will use henceforth. $\mathcal{C}_{\infty}(T)$
is referred to as the Franck-Condon factor and describes the overlap of the
displaced phonon state and the phonon ground state. Furthermore,
$\phi(\tau)=\int_{0}^{\infty}d\omega J_{\rm
pn}(\omega)\left(\cos(\omega\tau)\coth(\beta\omega/2)-i\sin(\omega\tau))\right/\omega^{2}$
is the phonon propagator and $\omega_{0}^{\prime}=\omega_{0}-\int d\omega
J_{\rm pn}(\omega)/\omega$, where the reorganization energy $\int d\omega
J_{\rm pn}(\omega)/\omega$ is also known as the polaron shift [57].

The subsequent dynamics at $\tau>\tau_{\rm P}$ is modeled within the markovian
approximation and encompasses both spontaneous decay and the quadratic
emitter-phonon coupling, as described by the Hamiltonian $H_{\rm MK}^{\rm
1e}=H_{0}+H_{\rm e-pt}+H_{{(\rm e-pn)}_{2}}+H_{\rm pt}+H_{\rm pn}$. This leads
in the usual way to a master equation of the Lindblad form,
${d\rho}/{dt}=\gamma\mathcal{D}_{\sigma}(\rho)+2\gamma_{\rm
pd}(T)\mathcal{D}_{\sigma^{\dagger}\sigma}(\rho)$, where
$\mathcal{D}_{\sigma}(\rho)=(\sigma\rho\sigma^{\dagger}-1/2\\{\sigma^{\dagger}\sigma,\rho\\}_{+})$
describes the spontaneous decay with a decay rate $\gamma$, and where
$\\{\\}_{+}$ is the anticommutator. The dissipator
$\mathcal{D}_{\sigma^{\dagger}\sigma}(\rho)$ describes the temperature-
dependent pure dephasing of the excited emitter state at the rate $\gamma_{\rm
pd}(T)$ due to the quadratic phonon interaction, which has been calculated for
GaAs quantum dots in [43, 44] (see also SM Sec. 1). Our concatenation approach
is consistent as long as markovian phonon pure dephasing and spontaneous decay
are negligible on the time scale $\tau_{\rm P}$, i.e.
$\mbox{exp}(-\Gamma\tau_{\rm P})\approx 1$, where
$\Gamma\equiv{\gamma}/{2}+\gamma_{\rm pd}(T)$. For times $\tau>\tau_{\rm P}$
we can apply the quantum regression theorem (QRT) [58] and find from our
concatenation approach at all times the field correlation (see SM Sec. 2 for
details)

$g^{(1)}(\tau)=\gamma\rho_{\rm
ee}\mathcal{C}(\tau)\mbox{e}^{-(\Gamma+i\omega_{0}^{\prime})\tau}.\;\;\;$ (3)

This first-order correlation function can be measured experimentally [55, 59].

Interestingly, for a single quantum emitter, polaron theory predicts exactly
the same two-time correlation as in (3) [60, 48, 43, 36, 61]. We will
introduce the polaron theory below in Sec. 2.2, but in brief, the equivalence
between both approaches is related to the fact that also in polaron theory, it
is assumed that the polaron formation happens on a faster time scale than all
other processes. We think it is simpler and offers more physical insight to
use our concatenation approach to arrive at Eq. (3). Regardless, both methods
agree on markovian and non-markovian effects due to the emitter-phonon
coupling.

For later use, we also briefly review the corresponding known single-emitter
frequency spectrum that can be obtained from the first-order correlation
function by utilizing the optical Wiener-Khinchin theorem [58]. In the
frequency spectrum $S(\omega)$, the fast initial phonon-assisted non-markovian
decay of the correlations results in broad meV-range sidebands, while the
subsequent slow markovian decay due to photons and phonons results in a finite
linewidth of the ZPL due to spontaneous emission and a further broadening of
the ZPL due to pure dephasing [62, 57, 42, 43, 50]. In SM Sec. 1 we depict the
emission spectrum using (3), demonstrating the temperature dependence of
photon emission into the ZPL and sidebands.

The total spectrum $S(\omega)$ can be written as $S_{\rm ZPL}(\omega)+S_{\rm
SB}(\omega)$, i.e. the sum of the contributions from the ZPL and from the
sidebands. The contribution of only the ZPL to the total spectrum can be
singled out by making in Eq. (3) the replacement
$\mathcal{C}(\tau)\rightarrow\mathcal{C}_{\infty}(T)$ [48, 6]. In this
approximation, the concatenation scheme reduces to the initial-slip scheme
where $\mathcal{C}_{\infty}(T)$ provides the initial condition for the
evolution of the markovian master equation. The approximation results in an
analytical expression for the spectrum near the zero-phonon line of the form
$S_{\rm ZPL}^{\rm 1e}(\omega)=\mathcal{C}_{\infty}(T)\gamma\rho_{\rm
ee}\Gamma/((\omega-\omega_{0}^{\prime})^{2}+\Gamma^{2})$. Thus
$\mathcal{C}_{\infty}(T)$ also measures the proportion of the emission that is
channeled into the ZPL [48, 61]. Materials and operating temperatures $T$ such
that $\mathcal{C}_{\infty}(T)$ is closer to unity are more interesting for
quantum information processing applications. To summarize this section, we
studied the dynamics of a single emitter interacting with both phonon and
photon baths, and we found an exact agreement between our concatenation and
polaron methods.

###  2.2 Scenario II : $N=2$

Next, we extend our approach to two identical emitters. Frequency tuning real
solid-state emitters such that they become optically identical is an art in
itself [30, 63]. For definiteness, we assume the two emitters to be embedded
in bulk GaAs and separated by a distance $r_{12}$, both interacting with the
photon and phonon baths. This can be described by the Hamiltonian

$H^{\rm 2e}=\sum_{n=1}^{2}\left(H_{0}^{(n)}+H_{\rm e-pt}^{(n)}+H_{({\rm
e-pn})_{1}}^{(n)}+H_{({\rm e-pn})_{2}}^{(n)}\right)+H_{\rm pt}+H_{\rm pn},$
(4)

where we use the same nomenclature as in Sec. 22.1 except that now the
coupling constants $h_{q}^{(n)}$, $g_{k}^{(n)}$ and $f_{k,k^{\prime}}^{(n)}$
depend on the position of the $n$-th emitter.

To solve the two-emitter dynamics, we first use our concatenation method. From
the Heisenberg equations of motion corresponding to the Hamiltonian (4), we
find the first-order correlation function
$g^{(1)}(\tau)=\sum_{n,m=1}^{2}\gamma_{nm}\langle\sigma_{n}^{\dagger}\sigma_{m}(\tau)\rangle$
where $\gamma_{nm}$ are the inter-emitter decay rates when $n\neq m$ and
single-emitter decay rates associated with the $n$-th emitter otherwise [19].
Similar to Sec. 22.1, the initial dynamics will again be dominated by the
linear emitter-phonon coupling, now given by the Hamiltonian $H_{\rm NM}^{\rm
2e}=\sum_{n=1}^{2}\left(H_{0}^{(n)}+H_{\rm e-pn}^{(n)}\right)+H_{\rm pn}$.
Then we can express the total correlation as (see SM Sec. 3),

$g^{(1)}(\tau)=\mathcal{C}(\tau)\sum_{n,m=1}^{2}\gamma_{nm}\rho_{nm}e^{-i\omega_{0}^{\prime}\tau},\;\tau\leq\tau_{\rm
P},$ (5)

where $\rho_{nm}=\langle n|\rho(0)|m\rangle$, with $|n\rangle$ representing
the state in the single excitation subspace.

After the non-markovian phonon correlations have reached a steady state, we
next consider in our concatenation method only the markovian effects due to
the photon and photon baths, given by the Hamiltonian $H_{\rm MK}^{\rm
2e}=\sum_{n=1}^{2}(H_{0}^{(n)}+H_{\rm e-pt}^{(n)}+H_{({\rm
e-pn})_{2}}^{(n)})+H_{\rm pt}+H_{\rm pn}$. Similar to the single-emitter case,
the emitter-photon and the quadratic phonon interactions can be accurately
described by a markovian master equation, which for the two identical emitters
can be written as
${d\rho}/{dt}=-i[H_{s},\rho]+\gamma_{+}\mathcal{D}_{\sigma_{+}}(\rho)+\gamma_{-}\mathcal{D}_{\sigma_{-}}(\rho)+2\gamma_{\rm
pd}(T)(\mathcal{D}_{\sigma_{+}^{\dagger}\sigma_{+}}(\rho)+\mathcal{D}_{\sigma_{-}^{\dagger}\sigma_{-}}(\rho))$.
This master equation is in diagonal form and describes both collective photon
emission and dephasing due to local phonon baths [19, 23]. The local phonon
correlations due to an excited emitter become uncorrelated at $\tau>\tau_{\rm
P}$ since $\mathcal{C}(\tau)\rightarrow\mathcal{C}_{\infty}(T)$. This gives a
phonon correlation length $L_{\rm P}=c\tau_{\rm P}$, where $c$ is the speed of
sound. Since we work in the single-excitation regime and we chose the inter-
emitter separations $r_{12}>L_{\rm P}$, the polaron formation occurs on a
shorter time than the travel time for phonons to adjacent emitters. Therefore
we assume uncorrelated inter-emitter phonon correlations implying independent
phonon baths. The diagonalization of the master equation was achieved by using
the symmetric operator $\sigma_{+}=(\sigma_{1}+\sigma_{2})/\sqrt{2}$ and its
antisymmetric counterpart $\sigma_{-}=(\sigma_{1}-\sigma_{2})/\sqrt{2}$. The
collective emission is characterized by the superradiant (subradiant) decay
rates $\Gamma_{+}=\gamma_{+}/2+\gamma_{\rm pd}(T)$
($\Gamma_{-}=\gamma_{-}/2+\gamma_{\rm pd}(T)$), where
$\gamma_{\pm}=\gamma\pm\gamma_{\rm col}$ since we assumed the emitters to be
identical and in a bulk environment, they have identical decay rates, i.e.,
$\gamma_{11}=\gamma_{22}\equiv\gamma$, while
$\gamma_{12}=\gamma_{21}\equiv\gamma_{\rm col}$ is the inter-emitter decay
rate. The factor $\mathcal{C}(\tau)$ in the correlation function (5)
incorporates the coherent evolution of emitters influenced by a linear
coupling of phonon modes at a picosecond timescale where the markovian effects
due to phonons and photons are irrelevant. This initial evolution does not
change emitter populations but leads to phonon-induced decoherence. On the
other hand, the effect of the photon bath at $\tau>\tau_{\rm P}$ in this two-
emitter configuration induces not only decay but also results in frequency
renormalization, called the Lamb shift, which is reflected in the coherent
part of the markovian master equation via the Hamiltonian $H_{s}=\omega_{\rm
col}(\sigma_{+}^{\dagger}\sigma_{+}-\sigma_{-}^{\dagger}\sigma_{-})$, where
$\omega_{\rm col}$ is the collective Lamb shift. For bulk GaAs, the inter-
emitter decay rates $\gamma_{\rm col}$ and collective Lamb shifts $\omega_{\rm
col}$ depend strongly on the distance $r_{12}$ between the embedded emitters;
in an H-aggregate configuration, where the parallel dipoles are orthogonal to
$\mathbf{r_{12}}$, they are given by $\gamma_{\rm
col}/\gamma\propto\sin\vartheta/\vartheta+\cos\vartheta/\vartheta^{2}-\sin\vartheta/\vartheta^{3}$
and $\omega_{\rm
col}/\gamma\propto-\cos\vartheta/\vartheta+\sin\vartheta/\vartheta^{2}+\cos\vartheta/\vartheta^{3}$.
Here, $\vartheta=\mbox{n}q_{0}^{\prime}r_{12}$ with $q_{0}^{\prime}$ the
wavevector associated to the phonon-renormalized emitter resonance frequency
$\omega_{0}^{\prime}$, and $\rm n$ is the refractive index of GaAs [64, 19].
The eigenstates associated with these two-emitter configurations in the
diagonal representation are $\omega_{0}^{\prime}\pm\omega_{\rm col}$.

We once again concatenate the non-markovian dynamics of (5) with the markovian
dynamics that is calculated by the QRT, to arrive at the first-order
correlation function at all times (see SM Sec. 3 for details)

$g^{(1)}(\tau)=\mathcal{C}(\tau)\sum_{l\in(+,-)}\gamma_{l}\rho_{ll}\mbox{e}^{-(\Gamma_{l}+i\omega_{l}^{\prime})\tau},\;\;\;\forall\tau,$
(6)

where $\omega_{l}^{\prime}=\omega_{0}^{\prime}\pm\omega_{\rm col}$ and
$\Gamma_{l}=\Gamma_{\pm}$. The calculation of (6) is made in the collective
basis where $\rho_{ll}=\langle\pm|\rho(0)|\pm\rangle$, with
$|+\rangle=(|e,g\rangle+|g,e\rangle)/\sqrt{2}$ is the collective symmetric and
$|-\rangle=(|e,g\rangle-|g,e\rangle)/\sqrt{2}$ the collective antisymmetric
state. Since the collective effects manifest themselves in the $\mu$eV range
around the ZPL, while phonon sidebands are spread across the meV range, we can
suppress these sidebands with a similar procedure as before for the single
emitter, which transforms the concatenation scheme to the simpler initial-slip
scheme [39, 41]. This aids in deducing an analytical expression for the
emission spectrum near the ZPL in the two-emitter configuration,

$S_{\rm ZPL}^{\rm
2e}(\omega)=\mathcal{C}_{\infty}(T)\sum_{l\in(+,-)}\gamma_{l}\rho_{ll}(t)\Gamma_{l}/\left((\omega-\omega_{l}^{\prime})^{2}+\Gamma^{2}_{l}\right),$
(7)

where we sum over the symmetric and antisymmetric eigenstates. The spectrum
(7) clearly will show the effects of collective emission, while also
reflecting non-markovian behavior of the phonon bath.

Now, turning to the polaron method, we subject the two-emitter Hamiltonian
$H^{\rm 2e}$ to the unitary transformation $\tilde{H}^{\rm
2e}=\mbox{exp}(S)H^{\rm 2e}\mbox{exp}(-S)$ where
$\mbox{exp}(S)=\mbox{exp}(\sum_{n=1}^{2}S_{n})$, and
$S_{n}=(g_{k}^{(n)}b^{\dagger}-g_{k}^{*(n)}b)/\omega_{k}$, to set the stage to
derive a polaron master equation (PE) (see SM Sec. 4 for derivation). From
$\tilde{H}^{\rm 2e}$ we can calculate the first-order correlation function as
$g^{(1)}(\tau)=\sum_{n,m=1}^{2}\mathcal{C}_{nm}(\tau)\gamma_{nm}\langle\sigma_{n}^{\dagger}\sigma_{m}(\tau)\rangle_{\rm
PE}$, where the two-time correlation function is evaluated with the help of
the polaron master equation. Here $\mathcal{C}_{nm}(\tau)$ is the phonon bath
correlation function defined as $\mathcal{C}_{nm}(\tau)=\mathcal{C}(\tau)$
when $n=m$ and $\mathcal{C}_{\infty}(T)$ when $n\neq m$ [6]. Expressing
$g^{(1)}(\tau)$ in the diagonal form using the symmetric ($\sigma_{+}$) and
antisymmetric operators ($\sigma_{-}$) as before, while neglecting the
sidebands, we obtain
$g^{(1)}(\tau)=\mathcal{C}_{\infty}(T)\sum_{l\in(+,-)}\gamma_{l}\langle\sigma_{l}^{\dagger}\sigma_{l}(\tau)\rangle_{\rm
DPE}$, where now the correlation is evaluated using the diagonal form of the
two-emitter polaron master equation (DPE)

$\frac{d\rho}{dt}=-i[H_{\rm
P},\rho]+\Upsilon_{+}\mathcal{D}_{\sigma_{+}}(\rho)+\Upsilon_{-}\mathcal{D}_{\sigma_{-}}(\rho)\\\
+2\gamma_{\rm
pd}(T)\left(\mathcal{D}_{\sigma_{+}^{\dagger}\sigma_{+}}(\rho)+\mathcal{D}_{\sigma_{-}^{\dagger}\sigma_{-}}(\rho)\right).$
(8)

Here $\Upsilon_{\pm}=\gamma\pm\mathcal{C}_{\infty}(T)\gamma_{\rm col}$ and
$\Omega_{\rm col}=\mathcal{C}_{\infty}(T)\omega_{\rm col}$ are the phonon-
renormalized temperature-dependent collective decay rate and Lamb shift
respectively (see SM Sec. 5 for details) and
$H_{P}=\sum_{l\in(+,-)}\Omega_{l}\sigma_{l}^{\dagger}\sigma_{l}$ where
$\Omega_{\pm}^{\prime}=\omega_{0}^{\prime}\pm\Omega_{\rm col}$ are the
eigenenergies. Here, one can notice that the single emitter decay rate is not
influenced by the phonon bath, which is a consequence of the markovian
approximation in the photon bath. On the other hand, the inter-emitter decay
rates do exhibit a non-additive phonon influence, see SM Sec. 5. Thus we can
evaluate the first-order correlation function through the QRT utilizing (8),
which results in the ZPL contribution to the spectrum as

$S_{\rm ZPL}^{\rm
2e}(\omega)=\mathcal{C}_{\infty}(T)\sum_{l\in(+,-)}\gamma_{l}\rho_{ll}(t)\digamma_{l}\big{/}\left((\omega-\Omega_{l}^{\prime})^{2}+\digamma^{2}_{l}\right),$
(9)

where $\digamma_{l}=\Upsilon_{l}/2+\gamma_{\rm pd}(T)$. This concludes the
derivation of the two-emitter spectra in the concatenation and the polaron
approaches, where in both cases we provided analytical approximations for the
ZPL parts.

####  2.2.1 Comparison of concatenation and polaron methods

In this section, we investigate the differences in the two approaches
(concatenation scheme and polaron method) in the ZPL spectra by using the
combined initial state of the two atoms as $|\psi(0)\rangle=|e,g\rangle$, a
localized excitation which equivalently can be written as
$|\psi(0)\rangle=1/\sqrt{2}(|+\rangle+|-\rangle)$. This choice of initial
state will result in a contribution from both terms in the emission spectrum
in (7) ((9)), showing up as two distinct peaks, a higher-energy superradiant
and a lower-energy subradiant peak, separated by the collective Lamb shift,
$\omega_{+}^{\prime}-\omega_{-}^{\prime}=2\omega_{\rm col}$
($\Omega_{+}^{\prime}-\Omega_{-}^{\prime}=2\Omega_{\rm col}$). The population
dynamics corresponding to this initial state manifests as bi-exponential
decay, as illustrated in SM Sec. 6.

Figure 1: Temperature-dependent emission spectra for two identical emitters
prepared in the initial state $|\psi(0)\rangle=|e,g\rangle$, separated by a
distance of $\lambda/25$, calculated using the concatenation approach (in
blue) and the polaron approach (in red). Simulation parameters: refractive
index $\rm n=3.5$; $\lambda=940$ nm; $\gamma=850$ MHz; $\alpha=0.025$;
$\omega_{\rm c}=1.49\rm ps^{-1}$ and the size of the quantum dots is $4.5$ nm.

In Fig. 1 we depict the ZPL emission spectra based on (7) ((9)) for four
different temperatures and $r_{12}=\lambda/25$ (see SM Sec. 6 for spectrum
with $r_{12}=\lambda/15$). We choose these short distances for bulk GaAs with
the aim to see a clear manifestation of superradiant and subradiant emission;
to observe similar spectra for longer distances, waveguide geometries can be
used instead [63]. The widths of the peaks in Fig. 1 correspond to
superradiant and subradiant decay rates, given by $\Gamma_{+}$ and
$\Gamma_{-}$ ($\digamma_{+}$ and $\digamma_{-}$). A very close agreement
between the two approaches can be observed. The agreement is best at the
lowest and highest temperatures shown, while at intermediate temperatures, the
agreement is still good. The decay rate enhancement $\Gamma_{+}/\Gamma_{-}$
can be used to quantify the collective emission [63], and it has the value of
$13.1\approx\digamma_{+}/\digamma_{-}$, at $T=1$ K and approaches unity at
high temperatures, as shown in Fig. 1. The emission spectra obtained with the
two methods are different, because the collective decay rate $\Upsilon_{\rm
col}(T)$ and collective Lamb shift $\Omega_{\rm col}(T)$ in the polaron
approach are temperature dependent due to phonon renormalization, while such a
renormalization is not accounted for in the concatenation approach (see SM
Sec. 5 and Ref. [6]). This results in slightly different widths of the
emission peaks and their respective resonance positions.

We will now quantify how much the ZPL spectra obtained with two methods agree.
Fig. 2(a) depicts the ratio of the polaron-renormalized collective emission
rates ($\Omega_{\rm col}$ and $\Upsilon_{\rm col}$) and the collective rates
obtained from the concatenation scheme ($\omega_{\rm col}$ and $\gamma_{\rm
col}$), which is equal to $\mathcal{C}_{\infty}$(T), which therefore measures
the variation in the resonance position and the width of the resonance peaks
in the two methods. Fig. 2(b) shows very similar overall decay rates in both
methods: since $\mathcal{C}_{\infty}(T)$ differs negligibly from unity at very
low temperatures, while at the highest temperatures the markovian phonon pure
dephasing dominates, with the same value for $\gamma_{\rm pd}(T)$ in both
methods, we find only small differences between the overall decay rates
$\Gamma_{\pm}$ and $\digamma_{\pm}$ at all temperatures. We can also define a
relative difference between the two methods as $\Delta\equiv 100(\rm CS-
PM)/CS$, where $\mbox{CS}=\Gamma_{+}+\Gamma_{-}+2\omega_{\rm col}$ constitutes
parameters deduced from the concatenation scheme while
$\mbox{PM}=\digamma_{+}+\digamma_{-}+2\Omega_{\rm col}$ constitutes the
polaron method. CS and PM both describe the sum of the widths of the super and
subradiant emission plus the energy difference between their associated
resonance energies. In Fig. 2(c), we plot $\Delta$ for three different emitter
separations. In all three cases, $\Delta$ starts at a low value at low
temperatures, culminating in its highest value at intermediate temperatures,
and then decreases monotonously for higher temperatures, in agreement with our
previous qualitative observations in Fig. 1. Also, a larger separation between
emitters results in smaller $\Delta$ (see SM Sec. 6). In summary, the
deviations between the ZPL spectra obtained with the concatenation and polaron
methods will for most purposes be negligible at all temperatures.

Figure 2: (a) Temperature dependence of the collective Lamb shift (collective
decay rate) obtained from the polaron method normalised with the collective
Lamb shift (collective decay rate), obtained within the concatenation scheme,
for two emitters separated by distance $0.05\lambda$. The ratio is equivalent
to the linear coupling factor $\mathcal{C}_{\infty}(T)$. (b) Magnitudes of the
superradiant decay rate $\digamma_{+}$ and $\Gamma_{+}$ obtained from the
polaron theory and our concatenating scheme for different temperatures. (c)
The percentage difference between the two schemes $\Delta$ for emitter
separations $\lambda/25$ (black line), $\lambda/20$ (red) and $\lambda/15$
(blue).

####  2.2.2 Consequence of neglecting non-markovian phonon dynamics

Until now we have compared collective emission in two models that both include
markovian as well as non-markovian dynamics due to phonon interactions. In our
modeling, following Refs. [42, 43, 44], the pure-dephasing rate originated
from the markovian approximation to the second-order phonon interaction, while
the first-order phonon interaction was responsible for the non-markovian fast
initial dynamics. We note in passing that it would take a one-dimensional
phonon bath to obtain from the first-order phonon interaction a nonzero pure
dephasing rate in the markovian approximation, while this approximation for a
higher-dimensional bath leads to a spurious decoherence-free subspace [65, 66,
67, 68, 69].

As remarked in the Introduction, the influence of phonons on collective
emission is typically neglected altogether. Or better, if phonon effects are
taken into account, then the state of the art is to only include them through
a pure-dephasing rate in a master equation. For example, in the remarkable
recent work by Tiranov et al. [63], where two distant quantum dots in a
photonic crystal waveguide were shown to emit light collectively, a markovian
master equation was utilized where the emission and emitter interaction were
described by the electromagnetic Green function, while a phenomenological
dephasing rate accounted for the interaction of the quantum dots with phonons.
Our models only reduce to such a markovian master equation by neglecting the
first-order interaction with phonons, while retaining the second-order
interaction. It is therefore not obvious at the outset that this unusual
approximation procedure gives accurate results.

Figure 3: (Top panel) Emission spectrum of (7) with and without the linear
phonon coupling, shown in solid and dashed lines, respectively. Two quantum
dots are prepared in the superradiant eigenstate $|+\rangle$ (red curves) and
the subradiant state $|-\rangle$ (black), at a temperature of (a) 4K and (b)
25K. The solid and dashed curves coincide in panel (a). Panels (c) and (d)
depict the corresponding concurrence dynamics.

Therefore, in Fig. 3, we compare sub- and superradiant spectra of our
concatenation model disregarding the sidebands, with the corresponding spectra
in the markovian model, as before for the two identical quantum dots in bulk
GaAs. In more detail, we compare the approximate spectra of (7) for the
concatenation model with its markovian approximation, which amounts to setting
$\mathcal{C}_{\infty}(T)$ in (7) to unity. We see that at 4K, the non-
markovian dynamics can well be neglected in the $\mu$eV-range around the ZPL,
while for 25K, the effect of neglecting the non-markovian dynamics is to
overestimate the amplitudes of the spectra, while the shapes of the spectra
are the same with or without the non-markovian dynamics in this frequency
range. Experimental spectra are often shown in arbitrary units, which thereby
become insensitive to non-markovian effects around the frequency range where
collective effects are probed.

So does this confirm that phonon-induced non-markovian dynamics is negligible?
The answer is negative: in the first place, because on a larger (meV) range
sidebands, the shapes of the spectra become different with or without non-
markovian dynamics as will be shown in Fig. 4 below. In the second place,
collective emission may be used as a signature of entanglement [63, 70], but
with quantum information applications in mind, the actual interest may be in
the entanglement dynamics itself.

Therefore, as the measure of two-qubit entanglement, we study concurrence,
defined as
$\mathcal{E}(\rho)=\mbox{max}\\{0,\sqrt{\alpha_{1}}-\sqrt{\alpha_{2}}-\sqrt{\alpha_{3}}-\sqrt{\alpha_{4}}\\}$
where $\alpha_{i}$ are the eigenvalues in descending order of their magnitude,
of the matrix $\rho\sigma_{\rm y}^{(1)}\otimes\sigma_{\rm
y}^{(2)}\rho^{*}\sigma_{\rm y}^{(1)}\otimes\sigma_{\rm y}^{(2)}$, with
$\sigma_{\rm y}^{(i)}$ the Pauli matrix related to the $i$-th qubit [71, 67].
In our concatenation model, the system density matrix at the time
$\tau=\tau_{\rm P}$ is given by

$\rho(\tau_{\rm P})=\sum_{n=0}^{2}|n\rangle\langle
n|\rho_{nn}+\mathcal{C}_{\infty}^{2}(T)\sum_{n\neq m=1}^{2}|n\rangle\langle
m|\rho_{nm},$ (10)

where $|0\rangle=|g,g\rangle$ is the collective ground state (see SM Sec. 7
for the derivation), which then serves as the initial state for the markovian
master equation for time $\tau>\tau_{\rm P}$. Notice the quadratic dependence
of the coherences on $\mathcal{C}_{\infty}(T)$, which also entails that the
concurrence will depend sensitively on any fast initial non-markovian
dephasing. In the fully markovian model, again $\mathcal{C}_{\infty}(T)$ is
taken to be unity.

In Fig. 3(c,d), we depict the two-emitter entanglement dynamics for the
initial states $|+\rangle$ and $|-\rangle$ corresponding to the temperature-
dependent spectra in the first two panels. We compare the concatenation model
of (6) with the fully markovian model. We can see that discarding the phonon-
induced non-markovian effects results in overestimating the concurrence. It is
a small effect at 4K, which is the operating temperature of Ref. [63], but a
large nonnegligible effect already at 25K. So, although the full markovian
model for the dynamics may suffice to describe the shape of the collective
emission spectra close to the ZPL, it fails to describe the entanglement
dynamics except at very low temperatures, where according to our models both
the markovian and the non-markovian dephasing can be strongly suppressed.

Finally, we conclude this section by comparing the full spectra $S(\omega)$,
i.e. also including the sidebands, obtained by numerical simulations in Fig. 4
pertaining to different methods discussed in this section. We find that the
concatenation scheme and the polaron method give very similar spectra, where
the largest (but still small) differences occur around the ZPL, the frequency
interval of Fig. 1, where the collective effects manifest themselves. The
figure also illustrates that the concatenation and the initial-slip schemes
indeed coincide around the ZPL. On the other hand, employing the markovian
master equation not only fails to reproduce the sidebands but also may
overestimates the amplitude of the spectra around the ZPL as clearly
demonstrated in the top right panel of Fig. 3.

Figure 4: (a) Comparison of emission spectra computed with the concatenation
model (solid blue line), the polaron method (dashed red line), the initial
slip scheme (dashed yellow line) and the markovian master equation (dashed
green line) at T = 4 K for two emitters prepared at an initial state
$|\psi(0)\rangle=|e,g\rangle$ with emitters separated by distance
$\lambda/25$. The range of frequencies in the main panel is in the meV range,
depicting the asymmetry in sidebands while the y-axis is on a logarithmic
scale. The inset shows a zoomed image of the same plot in the $\mu$eV range,
depicting the super- and subradiant splitting with the y-axis on a linear
scale, for a good comparison with the plots in Fig. 1. (b) A similar plot as
in panel (a) but for T = 25 K, resulting in symmetrical sidebands.

###  2.3  Scenario III: $N>2$

A one-dimensional chain of emitters placed over a photonic structure, for
instance, a nanofiber or a photonic crystal waveguide, is a promising platform
for observing and engineering collective effects [72, 73, 20, 74]. In this
section, we model collective emission in a 1D chain of $N$ equidistant
identical emitters in bulk GaAs, generalising the nomenclature and ideas
developed in Secs. 22.1 and 22.2. Thus, we first model the fast phonon
dynamics for $N$ emitters interacting with independent phonon baths given by
the Hamiltonian, $H_{\rm NM}^{\rm Ne}=\sum_{n=1}^{N}(H_{0}^{(n)}+H_{\rm
e-pn}^{(n)})+H_{\rm pn}$ and obtain an exact time evolution, for time
$\tau\leq\tau_{\rm P}$. From this exact solution, we obtain the first-order
correlation function that captures the non-markovian phonon dynamics as
follows (see SM Sec. 3 for details),

$g^{(1)}(\tau)=\mathcal{C}(\tau)\sum_{n,m=1}^{N}\gamma_{nm}\rho_{nm}e^{-i\omega_{0}^{\prime}\tau},\;\tau\leq\tau_{\rm
P},$ (11)

where $\rho_{nm}=\langle n|\rho(0)|m\rangle$, with the states,
$|n\rangle\equiv|g,g,...,e_{n},...,g\rangle$ being the localised states, where
only the $n$-th emitter is excited, together spanning the single-excitation
manifold, while $\gamma_{nm}$ are the inter-emitter decay rates. For
$\tau>\tau_{\rm P}$ we can once again use the Markovian description described
by the Hamiltonian $H_{\rm MK}^{\rm Ne}=\sum_{n=1}^{N}(H_{0}^{(n)}+H_{\rm
e-pt}^{(n)}+H_{({\rm e-pn})_{2}}^{(n)})+H_{\rm pt}+H_{\rm pn}$, analogous to
Sec. 22.2. From this Hamiltonian, one obtains the first-order correlation
function of the form
$g^{(1)}(\tau)=\sum_{n,m=1}^{N}\gamma_{nm}\langle\sigma_{n}^{\dagger}\sigma_{m}(\tau)\rangle_{\rm
ME}$, where the inter-emitter decay rates $\gamma_{nm}$ satisfy
$\gamma_{n,n+1}=\gamma_{\rm col}$ and $\gamma_{n,m}=\gamma_{m,n}$. The
correlation $\langle\cdot\cdot\cdot\rangle_{\rm ME}$ will be evaluated with
the aid of the QRT utilizing the multi-emitter markovian master equation (ME),
which we will derive in the nearest-neighbour interaction approximation in
order to arrive at some analytical results. (It is not necessary to make the
approximation for our numerical schemes to work, though.) As the name already
indicates, in the nearest-neighbour interaction approximation only the dipole-
dipole interaction $\omega_{nm}$ with nearest neighbours are considered, i.e.
$\omega_{nm}=\omega_{\rm col}$ if $|n-m|=1$ and zero otherwise, where $n$ and
$m$ label the $n$-th and $m$-th emitters, respectively [23]. Such a master
equation written in terms of atomic operators $\sigma_{n}$ is given by

$\frac{d\rho}{dt}=-i\sum_{n=1}^{N}\bigg{[}\omega_{\rm
col}(\sigma_{n}^{\dagger}\sigma_{n+1}+\sigma_{n+1}^{\dagger}\sigma_{n}),\rho\bigg{]}\\\
+\sum_{n,m=1}^{N}\gamma_{nm}\mathcal{D}_{\sigma_{nm}}(\rho)+2\gamma_{\rm
pd}(T)\sum_{n=1}^{N}\mathcal{D}_{\sigma_{n}^{\dagger}\sigma_{n}}(\rho),$ (12)

where the dissipator
$\mathcal{D}_{\sigma_{nm}}(\rho)=\sigma_{n}\rho\sigma_{m}^{\dagger}-\frac{1}{2}\\{\sigma_{n}^{\dagger}\sigma_{m},\rho\\}_{+}$.
For analytical as well as numerical calculations, it is convenient to use the
diagonal form of such a master equation. For $N=2$, obtaining a simultaneous
exact diagonalisation of the coherent and dissipative parts of the master
equation described in Sec. 22.2 was rather intuitive by using the symmetric
and antisymmetric combinations of the atomic operators to obtain the
collective operators. For $N>2$, generally, the simultaneous exact
diagonalisation of the the coherent and dissipative parts of the master
equation is not possible, as outlined in SM Sec. 3. We will resort to an
approximate treatment for both methods to obtain analytical results.

Assuming that each emitter interacts with its own local bulk phonon
environment, we can describe the problem in a collective basis in an
$N$-dimensional single-excitation manifold. The states $|n\rangle$ describing
localized single-photon excitations can equivalently be written in terms of
linear combinations of the collective states, $|l\rangle$ i.e.
$|n\rangle=\sum_{l=1}^{N}M_{nl}|l\rangle$, where
$M_{nl}=\sqrt{2/(N+1)}\sin\left(\pi nl/(N+1)\right)$ are orthonormal
coefficients [75, 23]. The collective states
$|l\rangle=\sigma_{l}^{\dagger}|0\rangle$ where
$|0\rangle\equiv|g,g,...,g\rangle$ as the collective ground state and
$\sigma_{l}^{\dagger}$ is the collective raising operator which excites the
collective ground state to the collective state $|l\rangle$. The atomic
operators similarly undergo the same transformation, i.e.
$\sigma_{n}=\sum_{l=1}^{N}M_{nl}|\sigma_{l}\rangle$. Thus, the $N$-emitter
master equation can be expressed in terms of the collective operators
$\sigma_{l}$ (details in SM Sec. 3), leading to a diagonal form of the master
equation corresponding to $N$ independent photon decay channels, each one
associated with one of the collective states. The resulting eigenstates are
given by $\omega_{l}=2\omega_{\rm col}\cos\left(\pi l/(N+1)\right)$, with the
corresponding collective decay rates
$\gamma_{l}=\sum_{n,m=1}^{N}\gamma_{nm}M_{nl}M_{ml}$. In the Dicke limit, i.e.
with zero separation between the emitters, $\gamma_{nm}=\gamma$ implying that
$\gamma_{1}\approx\gamma N$ (superradiant decay) and $\gamma_{N}\approx 0$
(subradiant decay) [75].

The total first-order correlation function in terms of the collective states
can be expressed as
$g^{(1)}(\tau)=\mathcal{C}_{\infty}(T)\sum_{l=1}^{N}\gamma_{l}\langle\sigma_{l}^{\dagger}\sigma_{l}(\tau)\rangle_{\rm
DME}$, where the correlation $\langle\cdot\cdot\cdot\rangle_{\rm DME}$ is now
evaluated using the QRT from the diagonalised multi-emitter master equation

$\frac{d\rho}{dt}=\sum_{l=1}^{N}\left(-i\left[\omega_{l}\sigma_{l}^{\dagger}\sigma_{l},\rho\right]+\gamma_{l}\mathcal{D}_{\sigma_{l}}(\rho)+2\gamma_{\rm
pd}(T)\mathcal{D}_{\sigma_{l}^{\dagger}\sigma_{l}}(\rho)\right).$ (13)

Concatenating the diagonal form of (11) and the markovian dynamics evaluated
via (13), we obtain the first-order correlation function at all times in the
collective basis

$g^{(1)}(\tau)=\mathcal{C}(\tau)\sum_{l=1}^{N}\gamma_{l}\rho_{ll}\mbox{e}^{-(\Gamma_{l}+i\omega_{l}^{\prime})\tau},\;\;\;\forall\tau,$
(14)

where $\rho_{ll}=\langle l|\rho(0)|l\rangle$ and
$\omega_{l}^{\prime}=\omega_{0}^{\prime}+2\omega_{\rm col}\cos\left(\pi
l/(N+1)\right)$. The emission spectrum for $N$ emitters neglecting the
sidebands evaluated from (14) can therefore be given as

$S_{\rm ZPL}^{\rm
Ne}(\omega)=\mathcal{C}_{\infty}(T)\sum_{l=1}^{N}\gamma_{l}\rho_{ll}\Gamma_{l}/\left((\omega-\omega_{l}^{\prime})^{2}+\Gamma^{2}_{l}\right),$
(15)

where the total collective decay rate is $\Gamma_{l}=\gamma_{l}/2+\gamma_{\rm
pd}(T)$.

In the polaron method we obtain an $N$-emitter polaron master equation (PE) by
subjecting the Hamiltonian describing an $N$-emitter phonon and photon
interaction via $H^{\rm Ne}=\sum_{n=1}^{N}\big{(}H_{0}^{(n)}+H_{\rm
e-pt}^{(n)}+H_{({\rm e-pn})_{1}}^{(n)}+H_{({\rm
e-pn})_{2}}^{(n)}\big{)}+H_{\rm pt}+H_{\rm pn}$, to the unitary polaron
transformation $\tilde{H}^{\rm Ne}=\mbox{exp}({S})H^{\rm Ne}\mbox{exp}(-S)$
where now $\mbox{exp}(S)=\mbox{exp}(\sum_{n=1}^{N}S_{n})$ with $S_{n}$ defined
in Sec. 22.2. We can thereby obtain the first-order correlation function by
utilizing the definition of the phonon bath correlations in Sec. 22.2 and
suppressing the sidebands as
$g^{(1)}(\tau)=\mathcal{C}_{\infty}(T)\sum_{n,m=1}^{N}\gamma_{nm}\langle\sigma_{n}^{\dagger}\sigma_{m}(\tau)\rangle_{\rm
PE}$, where the $N$-emitter polaron master equation (PE) is given by (see SM
Sec. 4 for details),

$\frac{d\rho}{dt}=-i\sum_{n=1}^{N}\bigg{[}\omega_{0}^{\prime}\sigma_{n}^{\dagger}\sigma_{n}+\Omega_{\rm
col}\left(\sigma_{n}^{\dagger}\sigma_{n+1}+\sigma_{n+1}^{\dagger}\sigma_{n}\right),\rho\bigg{]}\\\
+\sum_{n,m=1}^{N}\Upsilon_{nm}\mathcal{D}_{\sigma_{n,m}}(\rho)+2\gamma_{\rm
pd}(T)\sum_{n=1}^{N}\mathcal{D}_{\sigma_{n}^{\dagger}\sigma_{n}}(\rho),$ (16)

and where the inter-emitter decay rate $\Upsilon_{nm}=\gamma$ when $n=m$ while
$\Upsilon_{nm}=\gamma_{nm}\mathcal{C}_{\infty}(T)$ when $n\neq m$ (see SM Sec.
5 for details). In the diagonal representation the first-order correlation
takes the form
$g^{(1)}(\tau)=\mathcal{C}_{\infty}(T)\sum_{l=1}^{N}\gamma_{l}\langle\sigma_{l}^{\dagger}\sigma_{l}(\tau)\rangle_{\rm
DPE}$, where the correlations are computed using the multi-emitter diagonal
form of the polaron master equation (DPE)

$\frac{d\rho}{dt}=\sum_{l=1}^{N}\left(-i\left[\Omega_{l}^{\prime}\sigma_{l}^{\dagger}\sigma_{l},\rho\right]+\Upsilon_{l}\mathcal{D}_{\sigma_{l}}(\rho)+2\gamma_{\rm
pd}(T)\mathcal{D}_{\sigma_{l}^{\dagger}\sigma_{l}}(\rho)\right),$ (17)

which was obtained from the same diagonalization scheme as was used to
diagonalize (12), and where
$\Upsilon_{l}=\sum_{n,m=1}^{N}\Upsilon_{nm}M_{nl}M_{ml}$, and
$\Omega_{l}^{\prime}=\omega_{0}^{\prime}+2\Omega_{\rm col}\cos\left(\pi
l/(N+1)\right)$. Thus, the first-order correlation function takes the form

$g^{(1)}(\tau)=\mathcal{C}_{\infty}(T)\sum_{l=1}^{N}\gamma_{l}\rho_{ll}\mbox{e}^{-(\digamma_{l}+i\Omega_{l}^{\prime})\tau},\;\;\;\forall\tau,$
(18)

where $\digamma_{l}=\Upsilon_{l}/2+\gamma_{\rm pd}(T)$. From (18), we obtain
the emission spectrum

$S_{\rm ZPL}^{\rm
Ne}(\omega)=\mathcal{C}_{\infty}(T)\sum_{l=1}^{N}\gamma_{l}\rho_{ll}\digamma_{l}/\left((\omega-\Omega_{l}^{\prime})^{2}+\digamma^{2}_{l}\right).$
(19)

Figure 5: Emission spectrum for eight equidistant emitters on a line with each
emitter separated by a distance of $\lambda/25$ from the other, prepared in
the initial state $|\psi(0)\rangle=|e_{1},g,g,g,...\rangle$, manifesting
collective emission at different temperatures using our concatenation approach
(in blue) and the polaron approach (in red) for different temperatures.

Choosing $N=8$ and using as our initial state the localized single-photon
state $|\psi(0)\rangle=|e_{1},g,g,g,...\rangle$ results in the emission
spectrum depicted in Fig. 5, from the concatenation scheme (polaron method)
utilizing (15) ((18)). The two approaches show appreciable agreement for the
lowest and highest temperatures, with a maximum deviation at intermediate
temperatures that is still small, just as we found for two emitters in Sec.
22.2.

Although for $N=8$ there are eight eigenfrequencies $\omega_{l}$ (and
$\Omega_{l}$), in Fig. 5(a) only six distinct resonance peaks can be
distinguished since the peaks corresponding to the three highest eigenenergies
have merged together. These multi-peak spectral characteristics can be
regarded as signatures of collective phenomena. A complete absence of
collective effects would result in only one emission peak corresponding to
single-emitter decay. A decay-rate enhancement $\Gamma_{\rm max}/\Gamma_{\rm
min}=946\approx\digamma_{\rm max}/\digamma_{\rm min}$ at T = 1K is deduced,
which reduces to 6.85 already at 25K and approaches unity as the temperature
increases beyond 100K. This decay rate enhancement quantifies the collective
emission and is an indication of entanglement between the multiple emitters.

##  3 Conclusions

We presented two methods to describe collective light emission in photonic
media, which both account for markovian as well as non-markovian dephasing due
to interactions with phonons. We proposed and compared the concatenation
approach and the polaron method, and deduced simple analytical results for an
arbitrary number of emitters interacting with a common photon bath and
independent phonon baths. The methods give very similar results. While both
methods are computationally efficient, the concatenation approach is the
simplest of the two that takes into account the phonon effects without much
additional complexity compared to all-optical master equations.

Our proposed concatenation scheme concatenates the exact solution of the
initial non-Markovian dynamics due to phonons with the subsequent markovian
dynamics due to both photon and phonon baths as decribed by a multi-emitter
markovian master equation. Our analysis was carried out for emitters embedded
in a bulk GaAs environment. Utilising the bulk phonon spectral density for
GaAs allows us to identify an ultrafast cutoff time scale $\tau_{\rm P}$ on
which the phonon correlations reach a steady state, the key feature that made
our concatenation scheme possible. In the second, markovian part of the
dynamics, the effects of the photon and phonon baths are additive, in the
sense that they are described by independent exponential decay rates.

The other method, the polaron method, is a state-of-the-art technique that has
been utilised to study a single emitter interacting with phonons and photons.
We extended it for a multi-emitter 1D chain configuration exhibiting
collective effects by deriving a polaron master equation for an arbitrary
number of emitters. In general, in the polaron method the optical decay rates
are modified by the phonon bath, as an example of the so-called non-additive
effects of the two baths. However, here we considered the common situations
where local optical densities of states vary little across the emitter
linewidths (flat spectral density approximation), in which case the non-
additivity of single-emitter decay rates disappears, resulting in single-
emitter decay rates that are identical to the ones in the concatenation
scheme. On the other hand, inter-emitter decay rates (and hence also
collective emission rates) do exhibit the non-additive effects of the phonon
bath in the polaron method, whether one makes the flat photon spectral density
approximation or not, as a consequence of the assumption that the phonon baths
of different emitters are independent.

We deduced approximate analytical expressions for the emission spectrum around
the zero-phonon line, after employing a known approximate simultaneous
diagonalization scheme for the coherent and dissipative parts of the master
equation in both methods. We also made a rigorous quantitative comparison of
the two methods for different temperatures and found that they agree very
well. Thus using the polaron method instead of the simpler concatenation
method seems like overkill, at least for the flat photon spectral densities
considered here.

Decoherence due to the environment of solid-state emitters is often modeled as
a single markovian dephasing rate, which from the perspective of our
Hamiltonian model amounts to neglecting the linear phonon interaction (non-
markovian dephasing) while keeping the quadratic interaction (markovian
dephasing). We studied the accuracy of this procedure for various
temperatures. As a main result, we found that with a purely markovian model
for phonon effects one may considerably overestimate the amount of quantum
entanglement (measured as concurrence) between quantum emitters. However, for
quantum dots in bulk GaAs at temperatures down to a few degrees Kelvin, the
markovian model gave almost the same entanglement dynamics as our
concatenation model. Indeed, in the limit of zero temperature, both models
reduce to the same well-known all-optical multi-qubit master equations.

For simplicity, we considered bulk GaAs as the photonic environment. For
appreciable collective emission, we put the emitters at strongly subwavelength
distances, much shorter than in the recent experiments in photonic crystal
waveguides [63]. This illustrates that waveguide structures are better suited
for observing collective emission at longer distances. Our formalism can be
generalized to such inhomogeneous lossless media, by replacing the transverse
plane waves with the optical eigenmodes of those media. In most inhomogeneous
media that do not involve emitters resonantly coupled to optical cavities, the
flat photon spectral density approximation will still hold.

Finally, in the regime where we work it is natural to assume independent
phonon baths to model the collective effects in the presence of phonons.
Scenarios where correlated phonon baths as in Refs. [66, 76, 68, 77, 78, 79]
show their marks would introduce another source of non-markovian phonon
behaviour, and in the future it would be interesting to study their effects on
collective light emission as well.

Funding M.W. and D.P. acknowledge the Independent Research Fund Denmark
Natural Sciences (Project No. 0135-00403B). M.W. acknowledges the support by
the Danish National Research Foundation through NanoPhoton Center for
Nanophotonics, Grant No. DNRF147, and Center for Nanostructured Graphene,
Grant No. DNRF103. Acknowledgments We thank Mads Anders Jørgensen, Jakob
Hummelgaard, Tobias Egebjerg, Nicolas Stenger and Sanshui Xiao for stimulating
discussions. Disclosures The authors declare no conflicts of interest. Data
availability Data underlying the results presented in this paper are not
publicly available but may be obtained from the authors upon reasonable
request.

## References

  * [1] C. Lagoin, S. Suffit, K. Baldwin, L. Pfeiffer, and F. Dubin, “Dual-density waves with neutral and charged dipolar excitons of GaAs bilayers,” Nature Materials 22, 170–174 (2023).
  * [2] C. Toninelli, I. Gerhardt, A. S. Clark, A. Reserbat-Plantey, S. Götzinger, Z. Ristanović, M. Colautti, P. Lombardi, K. D. Major, I. Deperasińska, W. H. Pernice, F. H. Koppens, B. Kozankiewicz, A. Gourdon, V. Sandoghdar, and M. Orrit, “Single organic molecules for photonic quantum technologies,” Nature Materials 20, 1615–1628 (2021).
  * [3] K. Parto, S. I. Azzam, K. Banerjee, and G. Moody, “Defect and strain engineering of monolayer WSe2 enables site-controlled single-photon emission up to 150 K,” Nature Communications 12, 3585 (2021).
  * [4] C. Fournier, A. Plaud, S. Roux, A. Pierret, M. Rosticher, K. Watanabe, T. Taniguchi, S. Buil, X. Quélin, J. Barjon, J. P. Hermier, and A. Delteil, “Position-controlled quantum emitters with reproducible emission wavelength in hexagonal boron nitride,” Nature Communications 12, 3779 (2021).
  * [5] A. Fruchtman, R. Gómez-Bombarelli, B. W. Lovett, and E. M. Gauger, “Photocell optimization using dark state protection,” Physical Review Letters 117, 203603 (2016).
  * [6] D. M. Rouse, E. M. Gauger, and B. W. Lovett, “Optimal power generation using dark states in dimers strongly coupled to their environment,” New Journal of Physics 21, 063025 (2019).
  * [7] S. Tomasi, D. M. Rouse, E. M. Gauger, B. W. Lovett, and I. Kassal, “Environmentally improved coherent light harvesting,” Journal of Physical Chemistry Letters 12, 6143–6151 (2021).
  * [8] M. Moreno-Cardoner, R. Holzinger, and H. Ritsch, “Efficient nano-photonic antennas based on dark states in quantum emitter rings,” Optics Express 30, 10779 (2022).
  * [9] J. G. Bohnet, Z. Chen, J. M. Weiner, D. Meiser, M. J. Holland, and J. K. Thompson, “A steady-state superradiant laser with less than one intracavity photon,” Nature 484, 78–81 (2012).
  * [10] I. E. Protsenko, A. V. Uskov, E. C. André, J. Mørk, and M. Wubs, “Quantum Langevin approach for superradiant nanolasers,” New Journal of Physics 23, 063010 (2021).
  * [11] A. I. Lvovsky, B. C. Sanders, and W. Tittel, “Optical quantum memory,” Nature Photonics 3, 706–714 (2009).
  * [12] M. Scheibner, T. Schmidt, L. Worschech, A. Forchel, G. Bacher, T. Passow, and D. Hommel, “Superradiance of quantum dots,” Nature Physics 3, 106–110 (2007).
  * [13] K. Cong, Q. Zhang, Y. Wang, G. T. Noe, A. Belyanin, and J. Kono, “Dicke superradiance in solids [Invited],” Journal of the Optical Society of America B 33, C80 (2016).
  * [14] A. Angerer, K. Streltsov, T. Astner, S. Putz, H. Sumiya, S. Onoda, J. Isoya, W. J. Munro, K. Nemoto, J. Schmiedmayer, and J. Majer, “Superradiant emission from colour centres in diamond,” Nature Physics 14, 1168–1172 (2018).
  * [15] M. Biliroglu, G. Findik, J. Mendes, D. Seyitliyev, L. Lei, Q. Dong, Y. Mehta, V. V. Temnov, F. So, and K. Gundogdu, “Room-temperature superfluorescence in hybrid perovskites and its origins,” Nature Photonics 16, 324–329 (2022).
  * [16] R. H. Lehmberg, “Radiation from an N-atom system. I. General formalism,” Physical Review A 2, 883–888 (1970).
  * [17] T. Gruner and D. G. Welsch, “Green-function approach to the radiation-field quantization for homogeneous and inhomogeneous Kramers-Kronig dielectrics,” Physical Review A - Atomic, Molecular, and Optical Physics 53, 1818–1829 (1996).
  * [18] M. Wubs, L. G. Suttorp, and A. Lagendijk, “Multiple-scattering approach to interatomic interactions and superradiance in inhomogeneous dielectrics,” Physical Review A - Atomic, Molecular, and Optical Physics 70, 053823 (2004).
  * [19] Z. Ficek and S. Swain, _Quantum Interference and Coherence: Theory and Experiments_ , vol. 100 (Springer, Berlin, 2005).
  * [20] A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, “Exponential improvement in photon storage fidelities using subradiance & "selective radiance" in atomic arrays,” Physical Review X 7, 031024 (2017).
  * [21] Y. X. Zhang and K. Mølmer, “Subradiant emission from regular atomic arrays: Universal scaling of decay rates from the generalized Bloch theorem,” Physical Review Letters 125, 253601 (2020).
  * [22] J. Feist, A. I. Fernández-Domínguez, and F. J. García-Vidal, “Macroscopic QED for quantum nanophotonics: Emitter-centered modes as a minimal basis for multiemitter problems,” Nanophotonics 10, 477–489 (2020).
  * [23] M. Reitz, C. Sommer, and C. Genes, “Cooperative quantum phenomena in light-matter platforms,” PRX Quantum 3, 010201 (2022).
  * [24] M. Sánchez-Barquilla, F. J. García-Vidal, A. I. Fernández-Domínguez, and J. Feist, “Few-mode field quantization for multiple emitters,” Nanophotonics 11, 4363–4374 (2022).
  * [25] M. A. Jørgensen and M. Wubs, “Quantifying the breakdown of the rotating-wave approximation in single-photon superradiance,” Journal of Physics B: Atomic, Molecular and Optical Physics 55, 195401 (2022).
  * [26] N. Vats and S. John, “Non-Markovian quantum fluctuations and superradiance near a photonic band edge,” Physical Review A - Atomic, Molecular, and Optical Physics 58, 4168–4185 (1998).
  * [27] I. S. Nikolaev, P. Lodahl, A. F. Van Driel, A. F. Koenderink, and W. L. Vos, “Strongly nonexponential time-resolved fluorescence of quantum-dot ensembles in three-dimensional photonic crystals,” Physical Review B - Condensed Matter and Materials Physics 75, 115302 (2007).
  * [28] K. H. Madsen, S. Ates, T. Lund-Hansen, A. Löffler, S. Reitzenstein, A. Forchel, and P. Lodahl, “Observation of non-markovian dynamics of a single quantum dot in a micropillar cavity,” Physical Review Letters 106, 233601 (2011).
  * [29] R. Holzinger, S. A. Oh, M. Reitz, H. Ritsch, and C. Genes, “Cooperative subwavelength molecular quantum emitter arrays,” Physical Review Research 4, 033116 (2022).
  * [30] Z. X. Koong, M. Cygorek, E. Scerri, T. S. Santana, S. I. Park, J. D. Song, E. M. Gauger, and B. D. Gerardot, “Coherence in cooperative photon emission from indistinguishable quantum emitters,” Science Advances 8, eabm8171 (2022).
  * [31] J. Wiercinski, E. M. Gauger, and M. Cygorek, “Phonon coupling versus pure dephasing in the photon statistics of cooperative emitters,” Physical Review Research 5 (2023).
  * [32] A. J. Ramsay, A. V. Gopal, E. M. Gauger, A. Nazir, B. W. Lovett, A. M. Fox, and M. S. Skolnick, “Damping of exciton Rabi rotations by acoustic phonons in optically excited InGaAs/GaAs quantum dots,” Physical Review Letters 104, 017402 (2010).
  * [33] D. P. McCutcheon and A. Nazir, “Quantum dot Rabi rotations beyond the weak exciton-phonon coupling regime,” New Journal of Physics 12, 113042 (2010).
  * [34] S. Hughes and H. J. Carmichael, “Phonon-mediated population inversion in a semiconductor quantum-dot cavity system,” New Journal of Physics 15, 053039 (2013).
  * [35] J. H. Quilter, A. J. Brash, F. Liu, M. Glässl, A. M. Barth, V. M. Axt, A. J. Ramsay, M. S. Skolnick, and A. M. Fox, “Phonon-assisted population inversion of a single InGaAs/GaAs quantum dot by pulsed laser excitation,” Physical Review Letters 114, 137401 (2015).
  * [36] A. Nazir and D. P. McCutcheon, “Modelling exciton-phonon interactions in optically driven quantum dots,” Journal of Physics Condensed Matter 28, 103002 (2016).
  * [37] U. Geigenmüller, U. M. Titulaer, and B. U. Felderhof, “Systematic elimination of fast variables in linear systems,” Physica A 119, 41–52 (1983).
  * [38] F. Haake and R. Reibold, “Strong damping and low-temperature anomalies for the harmonic oscillator,” Physical Review A 32, 2462–2475 (1985).
  * [39] A. Suárez, R. Silbey, and I. Oppenheim, “Memory effects in the relaxation of quantum open systems,” Journal of Chemical Physics 97, 5101–5107 (1992).
  * [40] T. Yu, L. Diósi, N. Gisin, and W. T. Strunz, “Post-markov master equation for the dynamics of open quantum systems,” Physics Letters A 265, 331–336 (2000).
  * [41] Y. C. Cheng and R. J. Silbey, “Markovian approximation in the relaxation of open quantum systems,” Journal of Physical Chemistry B 109, 21399–21405 (2005).
  * [42] E. A. Muljarov and R. Zimmermann, “Dephasing in quantum dots: Quadratic coupling to acoustic phonons,” Physical Review Letters 93, 237401 (2004).
  * [43] A. Reigue, J. Iles-Smith, F. Lux, L. Monniello, M. Bernard, F. Margaillan, A. Lemaitre, A. Martinez, D. P. McCutcheon, J. Mørk, R. Hostein, and V. Voliotis, “Probing electron-phonon interaction through two-photon interference in resonantly driven semiconductor quantum dots,” Physical Review Letters 118, 233602 (2017).
  * [44] P. Tighineanu, C. L. Dreeßen, C. Flindt, P. Lodahl, and A. S. Sørensen, “Phonon decoherence of quantum dots in photonic structures: Broadening of the zero-phonon line and the role of dimensionality,” Physical Review Letters 120, 257401 (2018).
  * [45] M. Wubs, K. Saito, S. Kohler, P. Hänggi, and Y. Kayanuma, “Gauging a quantum heat bath with dissipative Landau-Zener transitions,” Physical Review Letters 97, 200404 (2006).
  * [46] A. Nazir, “Correlation-dependent coherent to incoherent transitions in resonant energy transfer dynamics,” Physical Review Letters 103, 146404 (2009).
  * [47] D. P. McCutcheon and A. Nazir, “Coherent and incoherent dynamics in excitonic energy transfer: Correlated fluctuations and off-resonance effects,” Physical Review B 83, 165101 (2011).
  * [48] J. Iles-Smith, D. P. McCutcheon, J. Mørk, and A. Nazir, “Limits to coherent scattering and photon coalescence from solid-state quantum emitters,” Physical Review B 95, 201305 (2017).
  * [49] H. Maguire, J. Iles-Smith, and A. Nazir, “Environmental nonadditivity and Franck-Condon physics in nonequilibrium quantum systems,” Physical Review Letters 123, 093601 (2019).
  * [50] E. V. Denning, J. Iles-Smith, N. Gregersen, and J. Mork, “Phonon effects in quantum dot single-photon sources,” Optical Materials Express 10, 222 (2020).
  * [51] D. M. Rouse, E. M. Gauger, and B. W. Lovett, “Analytic expression for the optical exciton transition rates in the polaron frame,” Physical Review B 105, 014302 (2022).
  * [52] J. S. Blakemore, “Semiconducting and other major properties of gallium arsenide,” Journal of Applied Physics 53, R123–R181 (1982).
  * [53] M. Cardona and N. E. Christensen, “Acoustic deformation potentials and heterostructure band offsets in semiconductors,” Physical Review B 35, 6182–6194 (1987).
  * [54] B. Krummheuer, V. M. Axt, and T. Kuhn, “Theory of pure dephasing and the resulting absorption line shape in semiconductor quantum dots,” Physical Review B 65, 1–12 (2002).
  * [55] P. Borri, W. Langbein, S. Schneider, U. Woggon, R. L. Sellin, D. Ouyang, and D. Bimberg, “Ultralong dephasing time in InGaAs quantum dots,” Physical Review Letters 87, 157401 (2001).
  * [56] U. Hohenester, “Quantum control of polaron states in semiconductor quantum dots,” Journal of Physics B 40, S315 (2007).
  * [57] G. D. Mahan, _Many-Particle Physics_ (Springer US, Boston, 2000).
  * [58] H. J. Carmichael and M. O. Scully, “Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations,” Physics Today 53, 78–80 (2000).
  * [59] A. J. Brash, J. Iles-Smith, C. L. Phillips, D. P. McCutcheon, J. O’Hara, E. Clarke, B. Royall, L. R. Wilson, J. Mørk, M. S. Skolnick, A. M. Fox, and A. Nazir, “Light scattering from solid-state quantum emitters: beyond the atomic picture,” Physical Review Letters 123, 167403 (2019).
  * [60] Y. Chassagneux, A. Jeantet, T. Claude, and C. Voisin, “Effect of phonon-bath dimensionality on the spectral tuning of single-photon emitters in the Purcell regime,” Physical Review B 97, 205124 (2018).
  * [61] J. Iles-Smith, D. P. McCutcheon, A. Nazir, and J. Mørk, “Phonon scattering inhibits simultaneous near-unity efficiency and indistinguishability in semiconductor single-photon sources,” Nature Photonics 11, 521–526 (2017).
  * [62] C. B. Duke and G. D. Mahan, “Phonon-broadened impurity spectra. I. Density of states,” Physical Review 139, A1965 (1965).
  * [63] A. Tiranov, V. Angelopoulou, C. J. van Diepen, B. Schrinski, O. A. Dall’Alba Sandberg, Y. Wang, L. Midolo, S. Scholz, A. D. Wieck, A. Ludwig, A. S. Sørensen, and P. Lodahl, “Collective super- and subradiant dynamics between distant optical quantum emitters,” Science 379, 389–393 (2023).
  * [64] J. B. Trebbia, Q. Deplano, P. Tamarat, and B. Lounis, “Tailoring the superradiant and subradiant nature of two coherently coupled quantum emitters,” Nature Communications 13, 2962 (2022).
  * [65] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger, “Dynamics of the dissipative two-state system,” Reviews of Modern Physics 59, 1–85 (1987).
  * [66] G. Massimo Palma, K. A. Suominen, and A. K. Ekert, “Quantum computers and dissipation,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 452, 567–584 (1996).
  * [67] R. Doll, M. Wubs, P. Hänggi, and S. Kohler, “Limitation of entanglement due to spatial qubit separation,” Europhysics Letters 76, 547–553 (2006).
  * [68] R. Doll, M. Wubs, P. Hänggi, and S. Kohler, “Incomplete pure dephasing of N-qubit entangled W states,” Physical Review B 76, 045317 (2007).
  * [69] R. Doll, P. Hänggi, S. Kohler, and M. Wubs, “Fast initial qubit dephasing and the influence of substrate dimensions on error correction rates,” European Physical Journal B 68, 523–527 (2009).
  * [70] F. Lohof, D. Schumayer, D. A. W. Hutchinson, and C. Gies, “Signatures of superradiance as a witness to multipartite entanglement,” Phys. Rev. Lett. 131, 63601 (2023).
  * [71] W. K. Wootters, “Entanglement of formation of an arbitrary state of two qubits,” Physical Review Letters 80, 2245–2248 (1998).
  * [72] E. Vetsch, D. Reitz, G. Sagué, R. Schmidt, S. T. Dawkins, and A. Rauschenbeutel, “Optical interface created by laser-cooled atoms trapped in the evanescent field surrounding an optical nanofiber,” Physical Review Letters 104, 203603 (2010).
  * [73] J. D. Hood, A. Goban, A. Asenjo-Garcia, M. Lu, S. P. Yu, D. E. Chang, and H. J. Kimble, “Atom-atom interactions around the band edge of a photonic crystal waveguide,” PNAS 113, 10507–10512 (2016).
  * [74] K. Sinha, P. Meystre, E. A. Goldschmidt, F. K. Fatemi, S. L. Rolston, and P. Solano, “Non-Markovian collective emission from macroscopically separated emitters,” Physical Review Letters 124, 043603 (2020).
  * [75] D. Plankensteiner, L. Ostermann, H. Ritsch, and C. Genes, “Selective protected state preparation of coupled dissipative quantum emitters,” Scientific Reports 5, 16231 (2015).
  * [76] J. H. Reina, L. Quiroga, and N. F. Johnson, “Decoherence of quantum registers,” Physical Review A 65, 032326 (2002).
  * [77] A. A. Valido, D. Alonso, and S. Kohler, “Gaussian entanglement induced by an extended thermal environment,” Physical Review A 88, 042303 (2013).
  * [78] C. Addis, P. Haikka, S. McEndoo, C. Macchiavello, and S. Maniscalco, “Two-qubit non-Markovianity induced by a common environment,” Physical Review A 87, 052109 (2013).
  * [79] J. Jeske and J. H. Cole, “Derivation of Markovian master equations for spatially correlated decoherence,” Physical Review A 87, 052138 (2013).

references See pages - of supplement.pdf


# Two Hamiltonian Cycles

Vaidy Sivaraman  Department of Mathematics and Statistics, Mississippi State
University, Mississippi State, MS 39762, USA<EMAIL_ADDRESS>and  Thomas
Zaslavsky  Dept. of Mathematical Sciences, Binghamton University (SUNY),
Binghamton, NY 13902-6000, USA<EMAIL_ADDRESS>
###### Abstract.

If the line graph of a graph $G$ decomposes into Hamiltonian cycles, what is
$G$? We answer this question for decomposition into two cycles.

###### Key words and phrases:

Hamilton decomposition, line graph

###### 2010 Mathematics Subject Classification:

Primary 05C45; Secondary 05C76

## Prelude

Two problems that are known to be NP-complete are: Is a cubic graph
Hamiltonian [2]? Does a graph $L$ decompose into two Hamiltonian cycles [8]?
An example for $L$ is the line graph of a cubic graph $G$: $L(G)$ has a
Hamilton decomposition if and only if $G$, the root graph, is Hamiltonian;
this was discovered by Kotzig [4] and later Martin [6]. (We say a graph is
_Hamilton decomposable_ if its edge set partitions into Hamiltonian cycles.)
The proofs are by polynomial-time constructions. Since Hamiltonicity of a
cubic graph is NP-complete, it follows that Hamiltonian decomposability of a
4-regular graph remains NP-complete when restricted to line graphs of cubic
graphs.

An opposite question seems not to have been asked. (It was briefly considered
by the second author for an undergraduate graph theory examination.)

######  Problem 1.

The line graph $L(G)$ of a simple graph $G$ is decomposable into two
Hamiltonian cycles. What is $G$?

We answer this question by characterizing all root graphs $G$. One would like
to describe all Hamilton decomposable line graphs. We offer some thoughts at
the end of this paper.

## Presto

######  Theorem 2.

The graph $G$ of Problem 1 is either $K_{1,5}$, or the first subdivision of a
$4$-regular graph $G^{\prime}$ that decomposes into two Hamiltonian cycles, or
a Hamiltonian cubic graph.

The first subdivision $SG$ of a graph $G$ is obtained by subdividing every
edge into a path of length 2. The 4-regular graph $G^{\prime}$ is not assumed
to be simple.

The first case gives line graph $K_{5}$, which is obviously a solution to
Problem 1. In the second case $G=SG^{\prime}$ is not itself Hamiltonian; we
shall explain how its line graph decomposes in terms of a known Hamilton
decomposition of $G^{\prime}$. If $G^{\prime}$ has $n$ vertices, there are
$2^{n}$ different ways to decompose $L(G)$ using that decomposition of
$G^{\prime}$. In the third case we assume a Hamiltonian cycle is known. If
there are $2n$ vertices in $G$, there are $2^{n}$ different ways to decompose
$L(G)$ using that Hamiltonian cycle.

The reductions from the line graph and its Hamiltonian cycles to the root
graph and its Hamiltonian cycle(s) can be performed in polynomial time, since
a graph is reconstructible from its line graph in polynomial time [5].

## Moderato

We describe the construction of the Hamilton decomposition of our 4-regular
$L(G)$ in each case other than the familiar $K_{5}$. We give the cubic
construction, from [4, 6], for completeness.

######  Example 3.

This example is more general than is needed for Theorem 2. Suppose $G$ is the
first subdivision of a $2h$-regular graph $G^{\prime}$, $h\geq 2$, that
decomposes into $h$ Hamiltonian cycles $H_{1}^{\prime},\ldots,H_{h}^{\prime}$.
We convert $L(SH_{i}^{\prime})$ into a Hamiltonian cycle $H_{i}$ of $L(G)$ as
follows.

A divalent vertex in $G$, incident with edges, say, $c$ and $d$ derived by
subdividing $c^{\prime}$ in $G^{\prime}$, gives an edge $cd$ that we put in
$H_{i}$ if $c^{\prime}\in H_{i}^{\prime}$.

Consider a $2h$-valent vertex $u$ in $G$. At $u$ there are two edges of $G$
belonging to each $SH_{i}^{\prime}$, call them $e_{i}$ and $f_{i}$. The edges
$e_{1},f_{1},\ldots,e_{h},f_{h}$ become vertices of a $2h$-clique $K_{2h}(u)$
in $L(G)$. The graph $K_{2h}(u)$ decomposes into $h$ Hamiltonian paths
$H_{i}(u)$ with distinct endpoints that can be chosen arbitrarily [3, proof of
Theorem 2.3.3]; we choose them to have respective endpoints $e_{i},f_{i}$ for
each $i=1,\ldots,h$. Replace the edge $e_{i}f_{i}$ in $L(SH_{i}^{\prime})$ by
the path $H_{i}(u)$. Doing this for all $2h$-valent vertices in $G$ gives $h$
edge-disjoint 2-factors of $L(G)$: $H_{i}$ derived from $L(SH_{i}^{\prime})$
for each $i$. Each $H_{i}$ is 2-regular by construction and it is connected
because $H_{i}^{\prime}$ is connected; therefore it is a Hamiltonian cycle of
$L(G)$. We have a decomposition because each edge of $L(G)$ has been assigned
to exactly one $H_{i}$.

In [7] Muthusamy and Paulraja proved half the cases of the similar result that
$L(G^{\prime})$ itself decomposes into Hamiltonian cycles, and indeed more
generally that if $G^{\prime}$ is any Hamilton decomposable graph then
$L(G^{\prime})$ is also Hamilton decomposable. Example 3 differs in that we
subdivide $G^{\prime}$ first, which makes the root graph no longer Hamilton
decomposable.

######  Example 4.

Suppose $G$ is a cubic graph with a Hamiltonian cycle $H$. Let $M$ be the
complementary matching: $M=E(G)\setminus E(H)$. Orient the edges of $M$
arbitrarily. Expand $H_{0}=L(H)$ into two Hamiltonian cycles of $L(G)$, as
follows.

Begin with two copies of $H_{0}$, labeled $H_{1}$ and $H_{2}$. Consider an
edge $ef$ of $L(H)$ corresponding to two edges of $H$ with common vertex $u$
and let $m=uu^{\prime}$ be the matching edge at $u$. The vertex $u^{\prime}$
belongs to two edges of $H$, say $e^{\prime}$ and $f^{\prime}$, forming the
edge $e^{\prime}f^{\prime}$ in $H_{0}$. If $u$ is the tail of $uu^{\prime}$,
replace the edge $ef$ of $H_{1}$ by the path $emf$ and replace the edge
$e^{\prime}f^{\prime}$ in $H_{2}$ by $e^{\prime}mf^{\prime}$, but in $H_{2}$
retain the edge $ef$ and in $H_{1}$ retain the edge $e^{\prime}f^{\prime}$.
Doing this for all matching edges makes $H_{1}$ and $H_{2}$ into edge-disjoint
Hamiltonian cycles of $L(G)$.

The decompositions depend on the choice of Hamilton decomposition in Example 3
and of Hamiltonian cycle in Example 4. Thus, in general there may be a great
many Hamilton decompositions of $L(G)$ besides the $2^{n}$ mentioned earlier.

## Andante

We prove that the examples in Theorem 2 are the only ones. We assume given a
line graph $L=L(G)$, but not the root graph $G$, and a decomposition of $L$
into Hamiltonian cycles $H_{1}$ and $H_{2}$. A root graph can be calculated
quickly from $L$ [5], thus the vertex cliques in $L$ are known. In particular,
we know when $G=K_{1,5}$ because then $L$ is a 5-clique. If it is not, then we
know when $G$ is $(4,2)$-biregular by the existence of 4-cliques in $L$, and
$G$ is cubic in the remaining case. Thus, we can consider each type of root
graph separately. We leave the case $L=K_{5}$ to the reader.

_Suppose $G$ is $(4,2)$-biregular._ We can reconstruct $G^{\prime\prime}$, a
graph isomorphic to $G^{\prime}$, directly from the 4-cliques in $L$. The
vertices $k_{i}$ of $G^{\prime\prime}$ are the 4-cliques in $L$, and there is
an edge $k_{i}k_{j}$ of $G^{\prime\prime}$ for each edge between a vertex of
$L$ in $k_{i}$ and one in $k_{j}$. We know $G^{\prime\prime}\cong G^{\prime}$
because each vertex of $G^{\prime\prime}$ is the vertex clique of $L=L(G)$
that corresponds to a quadrivalent vertex in $SG^{\prime}$, that is, a vertex
of $G^{\prime}$.

The graph $G^{\prime\prime}$ enables us to deduce the two Hamiltonian cycles
$H_{i}^{\prime\prime}$ of $G^{\prime\prime}$ from the $H_{i}$ in $L$. Since a
vertex in $G^{\prime\prime}$, considered as a 4-clique in $L$, is
quadrivalent, each Hamiltonian cycle $H_{i}$ of $L$ must enter and leave that
4-clique exactly once using separate edges from those used by the other one.
Thus, $H_{i}$ acts as a divalent subgraph $H_{i}^{\prime\prime}$ of
$G^{\prime\prime}$. This $H_{i}^{\prime\prime}$ is connected because $H_{i}$
is, hence it is a Hamiltonian cycle of $G^{\prime\prime}$. That proves the
characterization of $G$ in the $(4,2)$-biregular case of $L$.

Note that, if we begin with $G^{\prime}$ and its Hamiltonian cycles
$H_{i}^{\prime}$, construct $L=L(SG^{\prime})$ and its Hamiltonian cycles
$H_{i}^{\prime}$ as in Example 3, and then construct $H_{i}^{\prime\prime}$ in
$G^{\prime\prime}$, then the natural isomorphism of $G^{\prime\prime}$ to $G$
carries each $H_{i}^{\prime\prime}$ to $H_{i}$.

_Suppose $G$ is cubic._ We reconstruct the Hamiltonian cycle in $G$ from $L$
by examining the vertex triangles in $L$. No such triangle can have all three
edges in a single one of $H_{1}$ and $H_{2}$, so there are two kinds of vertex
triangle: type 1 has two edges from $H_{1}$ and type 2 has two edges from
$H_{2}$. Suppose $e_{1},e=e_{2},e_{3}$ are the vertices of a type 1 vertex
triangle in $L$, derived from a vertex $v$ of $G$, let $u_{1},w,u_{3}$ be the
other endpoints of $e_{1},e,e_{3}$ in $G$, and say the path $e_{1}ee_{3}$ is
in $H_{1}$, so $H_{2}$ contains the path $e_{1}e_{3}$. The trace of $H_{1}$ in
$G$ follows the re-entrant trajectory $u_{1}vwvu_{3}$ (in some direction) and
that of $H_{2}$ follows the simple trajectory $u_{1}vu_{2}$. If we remove the
re-entrance from the trace of $H_{1}$, then $H_{1}$ and $H_{2}$ are identical.

The other endpoint of $e$ in $G$, namely $w$, forms its own vertex triangle in
$L$ with vertices $f_{1},e,f_{3}$. Let $x_{1},x_{3}$ be the other endpoints of
$f_{1},f_{3}$ in $G$. In $L$, $e$ is incident with the four edges
$ee_{1},ee_{3},ef_{1},ef_{3}$. The first two belong to $H_{1}$, so the second
pair must belong to $H_{2}$. Therefore, the trace of $H_{2}$ in $G$ contains
the re-entrant trajectory $x_{1}wuwx_{3}$ and that of $H_{1}$ contains the
path $x_{1}wx_{3}$. The re-entrant edge here is the same edge $e=uw$ as in the
trace $u_{1}vwvu_{3}$ of $H_{1}$ from the $v$-vertex triangle. This implies
that the re-entrant edges are a well-defined set and a 1-factor of $G$. The
complementary 2-factor of $G$ is the simplified trace of both $H_{1}$ and
$H_{2}$ (that is, the trace after removing re-entrant edges), thus it is
connected, so it is a Hamiltonian cycle of $G$. This proves the
characterization of $G$ in the 4-regular case of $L$.

It is again easy to see that if we begin with cubic $G$ having a Hamiltonian
cycle $H$, apply the construction of $L(G)$ and its Hamilton decomposition
$H_{1}$ and $H_{2}$ into as in Example 4, and return to $G$ by the preceding
construction of a new Hamiltonian cycle derived from $H_{1}$ and $H_{2}$, the
new Hamiltonian cycle of $G$ is nothing but $H$ again.

## Coda

One naturally thinks about characterizing all line graphs that are Hamilton
decomposable, say into $h$ Hamiltonian cycles. The root graph is then
$(k+1,2h-k+1)$-biregular for some $k\leq h$. The star, the case $k=0$,
generalizes to $K_{2h+1}$, whose line graph $K_{2h+1}$ decomposes into $h$
Hamiltonian cycles. The 4-regular solution generalizes to the case $k=1$,
where we have a $(2h,2)$-biregular root graph as in Example 3. In the case
$k=h$ we have the line graph of an $h+1$-regular graph $G$, generalizing our
cubic root graph. If $G$ can be proved Hamilton decomposable (for odd $h$) or
Hamilton decomposable aside from a 1-factor (for even $h$), then this case is
solved by [7]. If $G$ can be proved merely to be Hamiltonian or to have a
Hamiltonian 3-factor, this same case is solved by [1]. Aside from that, there
is the virtual certainty that each value of $k$ introduces a new family of
solutions.

## References

  * [1] Darryn Bryant, Sara Herke, Barbara Maenhaut, Benjamin R. Smith, Hamilton decompositions of line graphs. Submitted. arXiv:2012.00988.
  * [2] M.R. Garey, D.S. Johnson, R. Endre Tarjan, The planar Hamiltonian circuit problem is NP-complete. SIAM J. Comput. 5(4) (1976), 704–714. MR 444516.
  * [3] Nora Hartsfield and Gerhard Ringel, _Pearls in Graph Theory_. Academic Press, Boston, 1990. Rev. ed., 1994. Repr. Dover Publications, Mineola, N.Y., 2003. MR 1282717. Zbl 703.05001.
  * [4] Anton Kotzig, Aus der Theorie der endlichen regulären Graphen dritten und vierten Grades. (In Slovak.) Časopis Pěst. Mat. 82 (1957), 76–92. MR 90815.
  * [5] Philippe G.H. Lehot, An optimal algorithm to detect a line graph and output its root graph. J. Assoc. Comput. Mach. 21 (1974), 569–575. MR 347690.
  * [6] Pierre Martin, Cycles hamiltoniens dans les graphes 4-réguliers 4-connexes. Aequationes Math. 14(1-2) (1976), 37–40. MR 414442.
  * [7] A. Muthusamy and P. Paulraja, Hamiltonian cycle decomposition of line graphs and a conjecture of Bermond. J. Combin. Theory Ser. B 64(1) (1995), 1–16. MR 1328291. Zbl 831.05056.
  * [8] B. Péroche, NP-completeness of some problems of partitioning and covering in graphs. Discrete Appl. Math. 8(2) (1984), 195–208. MR 743024.