$\displaystyle
C_{\text{eff}}(\Delta)=\frac{\textbf{G}^{\top}\textbf{G}}{\kappa\gamma_{\text{eff}}(\Delta)}.$
(85)

As mentioned in Sec. A, for a single emitter the cooperativity is independent
on the emitter properties (dipole moment, decay rate) and simply depends on
the cavity finesse. For $\mathcal{N}$ uncoupled emitters the same holds true
and only a linear increase with $\mathcal{N}$ will be obtained. For coupled
emitters, an interesting decoupling of the dipole moment in the direction of
the cavity (quantified by the term $\textbf{G}^{\top}\textbf{G}$) from the
collective radiative properties $\gamma_{\text{eff}}(\Delta)$ can be achieved.
As $\gamma_{\text{eff}}$ is not a natural constant of the ensemble, but
strongly dependent on the relative positioning and phase of individual
emitters, one can reach subradiant states with $\gamma_{\text{eff}}\ll\gamma$.
By proper design of the cavity transverse field amplitude profile, the
numerator can at the same time be maximized, resulting in a scaling up of
$C_{\text{eff}}$ well above the independent emitter case
$\mathcal{N}g^{2}/(\kappa\gamma)$. This is illustrated in Fig. 11a where a
scaling with $\mathcal{N}^{4}$ for a 1D chain is shown possible. In Figs. 11b
and c, we consider scans of the cavity transmission for a system of four
coupled emitters. The system will have three subradiant and one superradiant
state with resonances given by Eq. (14) (the superradiant state is located at
$\omega_{0}+2\Omega\cos(\pi/5)$). We adress the system either
antisymmetrically with $\mathbf{G}=(g,-g,g,-g)^{\top}$ [Fig. 11b] or
symmetrically with $\mathbf{G}=(g,g,g,g)^{\top}$ [Fig. 11c]. While
antisymmetric excitation overlaps with two subradiant states, symmetric
excitation overlaps with one subradiant and one superradiant state. In Ref.
Plankensteiner _et al._ (2017) it is shown that phased excitation could be
realized by using the higher-order TEM-modes of the optical cavity.
Alternatively, one could use symmetric addressing around the optimal points
utilized experimentally in Ref. Rui _et al._ (2020) or symmetric addressing
combined with a magnetic field gradient as in Ref. Plankensteiner _et al._
(2015).

###  B Hybrid cavities with 2D subwavelength mirrors

Figure 12: (a) A hybrid cavity comprised of an atomic array (left mirror) and
a standard, flat frequency mirror (right). The cavity mode $a$ interacts with
two external continua $b(\omega)$, $c(\omega)$ and the surface-confined mode
on the array $d$. (b) Cavity response for $\zeta_{0}=10$. The gray dashed line
is the expected Lorentzian profile with linewidth
$\kappa_{\text{flat}}=\omega_{\text{FSR}}/(2\pi\zeta_{0}^{2})$ in the
Markovian regime where $\gamma_{d}=\omega_{\text{FSR}}$. For narrower mirror
linewidths $\gamma_{d}=\omega_{\text{FSR}}/10$ (blue line) and
$\gamma_{d}=\omega_{\text{FSR}}/100$ (red line) the cavity shows an asymmetric
transmission profile with a strongly reduced linewidth at the level of
$\gamma_{d}/\zeta_{0}$. (c) Transmission profile for a double-sided hybrid
cavity, where also the right mirror is replaced by an atomic array (same
values for $\gamma_{d}$ as in (b)) shows a symmetric response with an even
narrower linewidth.

We now aim to provide a theoretical framework for hybrid cavities where
subwavelength arrays are to be used as end mirrors as suggested in Refs.
Shahmoon _et al._ (2017); Cernotik _et al._ (2019); Rui _et al._ (2020). The
standard approach for cavity loss and input-output relations taken in Sec. IV
is based on the assumption that the tunneling rate of photons between a
cavity-confined mode $a$ and the outside continuum of modes ($b(\omega)$ to
the left and $c(\omega)$ to the right) is flat around the cavity resonance.
This assumption is strongly modified in the case of subwavelength reflective
arrays and the presence of a surface-confined resonance (mode $d$) has to
explicitly be taken into account (see illustration in Fig. 12a). The roadmap
to construct an input-output quantum theory for such hybrid cavities follows
the steps described in Ref. Cernotik _et al._ (2019) where it has been applied
to photonic crystal mirrors. The procedure consists of three steps: i) the
derivation of equations of motion for cavity field operators based on a
coupled modes model where the photon-exchange processes between modes $a$ and
$d$ are included, ii) the derivation of the cavity amplitude transmission from
a classical transfer-matrix approach based on the expression of the
reflectivity of the array $r(\omega)$ derived in Sec. C and iii) the
extraction of the phenomenologically introduced parameters by matching the
predictions of the two theories. We apply this procedure here to a hybrid
cavity made of one flat mirror and one subwavelength array and point out
expected improvements in the cavity finesse where the hybrid design makes use
of two emitter array mirrors.

Coupled modes theory \- We consider the situation depicted in Fig. 12a where
the cavity mode $a$ is coupled to the surface-confined mode $d$ at some
complex rate $\mathcal{G}$ and the photon tunneling rates $U,V,W\in\mathbb{R}$
are frequency independent. Tunneling will then give rise to loss rates
$\kappa_{L}=\pi V^{2}$, $\kappa_{R}=\pi W^{2}$, $\gamma_{d}=\pi U^{2}$ (see
Appendix A1). The dynamics is then described by the Langevin equations

$\displaystyle\dot{a}$
$\displaystyle=-(i\omega_{a}+\kappa)a-\mathcal{G}d+\sqrt{2\kappa_{L}}b_{\mathrm{in}}+\sqrt{2\kappa_{R}}c_{\mathrm{in}},$
(86a) $\displaystyle\dot{d}$
$\displaystyle=-(i\omega_{d}+\gamma_{d})d-\mathcal{G}a+\sqrt{2\gamma_{d}}b_{\mathrm{in}},$
(86b)

showing that the mode $a$ is subject to two types of noises:
$b_{\mathrm{in}}$, entering through the left and $c_{\mathrm{in}}$, entering
through the right mirror. In a first approximation, we will consider that the
surface-confined mode reacts to the same noise term $b_{\mathrm{in}}$
consisting of modes on the left of the cavity but extra independent noise
terms stemming for example from nonradiative losses could be added. The total
decay rate is $\kappa=\kappa_{L}+\kappa_{R}$ and the input fields have the
usual nonvanishing correlation functions
$\braket{b_{\mathrm{in}}(t)b_{\mathrm{in}}^{\dagger}(t^{\prime})}=\delta(t-t^{\prime})$
[and similarly for $c_{\mathrm{in}}(t)$]. The associated output fields follow
the input-output relations
$b_{\mathrm{out}}=b_{\mathrm{in}}-\sqrt{2\kappa_{L}}a-\sqrt{2\gamma_{d}}d,$
and $c_{\mathrm{out}}=c_{\mathrm{in}}-\sqrt{2\kappa_{R}}a$.
We are mainly interested here in deriving the free parameters in the coupled
modes model, i.e., $\mathcal{G}$, $\omega_{a}$ and $\kappa_{L}$. For more in-
depth discussions highlighting the non-Markovianity of such hybrid cavities we
refer to Ref. Cernotik _et al._ (2019) where a full analysis of photonic
crystal mirror cavities is provided.

The transmission coefficient
$t(\omega)=\braket{c_{\mathrm{out}}(\omega)}/\braket{b_{\mathrm{in}}(\omega)}$
is easily expressed as
$t(\omega)=-\sqrt{2\kappa_{R}}\braket{a(\omega)}/\braket{b_{\mathrm{in}}(\omega)}$
in the case of sole driving through the left cavity port and can be
analytically derived from Eqs. (86) in the frequency domain

$t(\omega)=\frac{\sqrt{2\kappa_{R}}[\epsilon_{d}(\omega)\mathcal{G}\sqrt{2\gamma_{d}}-\sqrt{2\kappa_{L}}]}{\epsilon_{a}^{-1}(\omega)-\mathcal{G}^{2}\epsilon_{d}(\omega)},$
(87)

where we introduce the susceptibilities
$\epsilon_{a}^{-1}(\omega)=\kappa+i(\omega_{a}-\omega)$ and
$\epsilon_{d}^{-1}(\omega)=\gamma_{d}+i(\omega_{d}-\omega)$. This is the main
result of the coupled modes theory, which at this point still retains complete
generality; in the following we will apply it to the particular scenario of a
subwavelength reflective array.

Transfer matrix results \- The transfer matrix approach, on the other hand,
consists in solving the classical one-dimensional wave propagation in a one-
dimensional setup with two mirrors parametrized by the polarizabilities
$\zeta_{0}$ (right mirror) and $\zeta_{L}(\omega)$ (left, subwavelength array
mirror). In linear response theory, one can find the transmission function of
the setup for any incoming plane wave at a given frequency $\omega$. We use
the parametrization
$\zeta_{L}(\omega)=-ir(\omega)/(1+r(\omega))=\gamma_{d}/(\omega_{d}-\omega)$
with $r(\omega)=-i\gamma_{d}/((\omega-\omega_{d})+i\gamma_{d})$ for the
subwavelength array [cf. Eq. (55)] and choose a fixed $\zeta_{0}\gg 1$
(corresponding to a close to unity reflectivity) for the right mirror. The
classical transmission coefficient (see Appendix A4) then reads

$\tilde{t}(\omega)=\frac{1}{(1-i\zeta_{0})[1-i\zeta_{L}(\omega)]e^{-i\theta}+\zeta_{0}\zeta_{L}(\omega)e^{i\theta}},$
(88)

where $\theta=\omega\ell/c$. For $\zeta_{L}(\omega)$ and $\zeta_{0}$ infinite,
the expression above reaches unit absolute value at resonances
$\omega_{m}=m\times\omega_{\text{FSR}}$ for any positive integer $m$. The
resonances are separated by the free spectral range
$\omega_{\text{FSR}}=c\pi/\ell$. Let us first fix a given resonance number $m$
such that $\omega_{m}$ lies in the neighborhood of the mirror resonance and
see that for finite $\zeta_{0}$ and assuming $\zeta_{L}(\omega)$ is flat in
frequency and equal to $\zeta_{0}$, the transmission can reach unity at a
shifted $\omega^{\prime}_{m}=\omega_{m}+\omega_{\text{FSR}}/(\pi\zeta_{0})$.
The linewidth of such a resonance is then
$\kappa_{\text{flat}}=\omega_{\text{FSR}}/(2\pi\zeta_{0}^{2})$.
When the two mirrors have unequal reflectivities, the transmission of the
optical resonator is always less than unity. Therefore, for variable
susceptibility $\zeta_{L}$ we ask now that a resonance $\bar{\omega}$ should
be reached for which $|\zeta_{L}(\bar{\omega})|=\zeta_{0}$ such that the
transmission is unity. There are two solutions for this equation and we pick
the one where the mirror resonance sits to the left of the cavity resonance.
The reason is that the reflectivity at $\omega_{d}$ is exactly unity which
means the cavity transmission reaches zero. With the choice
$\bar{\omega}=\omega_{d}+\gamma_{d}/\zeta_{0}$ the zero of the cavity sits on
the left of the cavity resonance. We also fix $\bar{\omega}=\omega_{m}$, which
in practice means that one first determines $\gamma_{d},\zeta_{0}$ and then
adjusts the cavity length.
For varying $\gamma_{d}$, Fig. 12b shows a scan of the cavity resonance. For
$\gamma_{d}\gg\kappa_{\text{flat}}$, the expected Lorentzian response is
obtained (gray dashed line) while for decreasing $\gamma_{d}$ at the level of
$\kappa_{\text{flat}}$, a very narrow (around $\gamma_{d}/\zeta_{0}$)
asymmetric Fano profile is obtained. The zero of the hybrid cavity
transmission is at $\omega_{d}$ while unity is reached at $\omega_{m}$.
Notice that the double hybrid cavity has an even narrower linewidth as
illustrated in Fig. 12c for the same regimes as provided for Fig. 12b. In
addition, the transmission profile is symmetric and has the advantage of
presenting two Fano resonances, situated symmetrically with respect to the
cavity resonance.

Extraction of parameters \- Comparison of the two expressions in Eq. (88) and
Eq. (87) can allow for the identification of $\mathcal{G}$, $\omega_{a}$ and
$\kappa_{L}$. In a first step, we ask that the zero of the transmission at the
mirror resonance vanishes $t(\omega_{d})=0$ corresponding to the point at
which the mirror reflectivity is unity. This leads immediately to the
identification $\mathcal{G}=\sqrt{\kappa_{L}\gamma_{d}}$. The derived coupling
between the surface resonance and the cavity-confined mode is purely
dissipative, in stark contrast with the situation treated in Ref. Cernotik _et
al._ (2019), where an additional real photon exchange process between modes
$a$ and $d$ occurs. In the next step, we ask that unit transmission
$|t(\omega)|^{2}=1$ is reached at the resonance
$\omega_{m}=\omega_{d}+\gamma_{d}/\zeta_{0}$. From Eq. (87), after a few
calculations, we obtain a solution for the left mirror loss as
$\kappa_{L}=\kappa_{R}+\zeta_{0}(\omega_{a}-\omega_{m})+i|(\omega_{a}-\omega_{m})-\zeta_{0}\kappa_{R}|$.
As the loss rate is defined in the real domain, we ask for the imaginary part
to vanish which gives $\omega_{a}=\omega_{m}+\zeta_{0}\kappa_{R}$. This also
leads to $\kappa_{L}=\kappa_{R}(1+\zeta_{0}^{2})$. Further fitting of the
transmission profiles also identifies
$\kappa_{R}=(\kappa_{\text{flat}}/2)\gamma_{d}\zeta_{0}/(\omega_{\text{FSR}}+\gamma_{d}\zeta_{0})$.
With the proper definition of the parameters appearing in Eqs. (86), one can
proceed in solving various problems involving hybrid cavities with quantum-
emitter ensembles (eventually in the direction of Fano cavity lasing as in
Ref. Mork _et al._ (2014)) or with movable mirrors. In the direction of
quantum optomechanics, Ref. Cernotik _et al._ (2019) has shown that the Fano
profile of photonic crystal mirror cavities can be constructively utilized to
lead to cooling in the absence of heating. The mechanism is based on sideband
resolution by fitting a Stokes, heating sideband inside the Fano resonance.
For subwavelength emitter arrays the same can be realized by tailoring the
relationship between the cavity length and the array resonance. In addition,
the double-sided hybrid cavity provides a scenario in which both sidebands can
be inhibited as it presents symmetrically placed Fano dips.

###  C Superradiant lasers

The Tavis-Cummings Hamiltonian introduced in Sec. A together with the
collective spin algebra from Sec. D are sufficient for a basic understanding
of the main advantages presented by superradiant lasers (as treated in Refs.
Bohnet _et al._ (2012, 2014); Norcia _et al._ (2016)). First, we will assume a
minimal model for a laser comprised of a cavity containing an incoherently
pumped gain medium. Under certain conditions (implying the existence of a
pumping threshold) a non-zero intracavity field amplitude can be obtained,
which signifies a coherent-light output, i.e., lasing. The distinction between
the good and bad cavity regimes (behavior quantified by the cavity losses)
will then be seen to characterize the difference between standard and
superradiant lasers. Finally, we will make the connection between superradiant
lasing and superradiance as described by Eq. (21) in Sec. D by deriving an
effective master equation for the gain medium after eliminating the cavity
field.
The minimal model assumes $\mathcal{N}$ identical emitters undergoing
(independent) decay with the standard Lindblad term (as defined in Eq. (8))
$\mathcal{L}_{\gamma}$ (each emitter has a decay rate $\gamma$ with
corresponding collapse operator $\sigma_{i}^{\dagger}$). Additionally, the
incoherent pump assumes a Lindblad term $\mathcal{L}_{\gamma_{\text{p}}}$
(each emitter is pumped at a rate $\gamma_{\text{p}}$ with corresponding
collapse operator $\sigma_{i}^{\dagger}$ which brings population to the
excited state). The pump can for example be realized as illustrated in Fig.
13a by excitation into an intermediate level $\ket{i}$ followed by fast decay
to the lasing level $\ket{e}$. For $\gamma_{\text{p}}>\gamma$ population
inversion can be produced, however without exciting the emitter dipoles
$\braket{\sigma_{j}}=0$ as all processes are incoherent. To produce a non-zero
average dipole moment, which would correspond to the generation of nonthermal
light, the coupling to the optical resonator is crucial. From the Hamiltonian
in Eq. (60) (with $g_{j}=g$ and $\omega_{0}=\omega_{c}$) supplemented with the
two Lindblad terms (and in terms of collective operators defined in Sec. D),
we can derive the equations of motion:

$\displaystyle\braket{\dot{a}}$
$\displaystyle=-\kappa\braket{a}-ig\braket{S},$ (89a)
$\displaystyle\braket{\dot{S}}$
$\displaystyle=-(\gamma+\gamma_{\text{p}})\braket{S}+2ig\braket{aS_{z}},$
(89b) $\displaystyle\braket{\dot{S_{z}}}$
$\displaystyle=-2(\gamma+\gamma_{\text{p}})\braket{S_{z}}+ig\braket{a^{\dagger}S-aS^{\dagger}}+\mathcal{N}(\gamma_{\text{p}}-\gamma).$
(89c)

The above set of equations is generally not solvable. However, an expansion
around averages $a=\alpha+\delta a$, $S=s+\delta S$ and $S_{z}=s_{z}+\delta
S_{z}$ where all $\delta a,\delta S$ and $\delta S_{z}$ are zero-averaged
quantum noise terms, can lead to a simplification. The observation (which one
can eventually infer from numerical simulations) is that $\braket{\delta
a\delta S}\ll\alpha s$ (and similar for all other two-operator correlations)
in the limit that $\mathcal{N}$ is large. In such a case, a much simpler set
of equations can be obtained for averages

$\displaystyle\dot{\alpha}$ $\displaystyle=-\kappa\alpha-igs,$ (90a)
$\displaystyle\dot{s}$ $\displaystyle=-(\gamma+\gamma_{\text{p}})s+2ig\alpha
s_{z},$ (90b) $\displaystyle\dot{s}_{z}$
$\displaystyle=-2(\gamma+\gamma_{\text{p}})s_{z}+ig(\alpha^{*}s-\alpha
s^{*})+\mathcal{N}(\gamma_{\text{p}}-\gamma).$ (90c)

Figure 13: (a) A gain medium is represented by an incoherently pumped (via
intermediate state $\ket{i}$) three-level system with lasing occuring on the
$\ket{e}$-$\ket{g}$-transition. In a standard laser, stimulated emission
produces a large intracavity field amplitude exhibiting coherence past
threshold. In a superradiant laser, the intracavity field amplitude is instead
kept small and coherence is stored in the collective atomic dipole thus
protecting the laser linewidth from thermal effects. (b) Behavior of
intracavity photon number $|\alpha|^{2}$ and coherence $|s|^{2}$ in the
standard, good cavity regime ($\kappa=0.1g$) and (c) in the superradiant, bad
cavity regime ($\kappa=40g$). Other parameters are $\gamma=10^{-2}g$ and
$\mathcal{N}=200$. Notice that the upper threshold
$\gamma_{\text{p}}^{+}\approx\mathcal{N}g^{2}/\kappa$ is very different in
both regimes.

Let us now assume that a steady state can be reached (by setting all
derivatives to zero) in which case the solution of the above equations yields
a population difference $s_{z}=1/(2\mathcal{C}_{\text{p}})$, with
$\mathcal{C}_{\text{p}}=g^{2}/[\kappa(\gamma+\gamma_{\text{p}})]$, the cavity-
field intensity (photon number)

$\displaystyle|\alpha|^{2}$
$\displaystyle=\frac{1}{\kappa}\left[\frac{\mathcal{N}(\gamma_{\text{p}}-\gamma)}{2}-\frac{\gamma+\gamma_{\text{p}}}{2\mathcal{C}_{\text{p}}}\right],$
(91)

and the dipole coherence $|s|/|\alpha|=\kappa/g$. The expression above is
positive only above a lasing threshold (signifying the existence of a lasing
steady state), which can be derived by finding the required values for the
pumping strength $\gamma_{\text{p}}$. Requiring that in the lasing regime one
has $|\alpha|^{2}>0$, yields a quadratic equation with upper and lower
thresholds given by

$\displaystyle\gamma_{\text{p}}^{\pm}=\frac{g^{2}\mathcal{N}}{2\kappa}-\gamma\pm\sqrt{\frac{g^{2}\mathcal{N}}{\kappa}\left(\frac{g^{2}\mathcal{N}}{4\kappa}-2\gamma\right)}.$
(92)

If $\mathcal{N}g^{2}/\kappa$ is large compared to $\gamma$, the thresholds can
be approximated as $\gamma_{\text{p}}^{-}\approx\gamma$,
$\gamma_{\text{p}}^{+}\approx\mathcal{N}g^{2}/\kappa$. Note that the lower
threshold corresponds to achieving population inversion while the upper
threshold implies that a too large incoherent pumping rate $\gamma_{\text{p}}$
prevents the formation of coherence in the system.
While the analysis above is valid in any regime, the great advantage brought
on by superradiant lasers is that the coherence is stored in the gain medium
which means that the laser linewidth is unaffected by common problems in
optical resonators operating at large intracavity powers (such as thermal
vibrations of the end mirrors). This can be seen in the ratio
$|s|/|\alpha|=\kappa/g$ which indicates that for $\kappa\gg g$ the intracavity
field can be negligible while the gain medium coherence is very large [see
Figs. 13b,c].
Let us now provide an alternative understanding of the superradiant lasing
problem, by connecting to Sec. D which we do by focusing on sketching the
derivation of an effective master equation for the $\mathcal{N}$ emitters in
the case of a lossy cavity. We assume at time $t=0$ the emitters to be fully
inverted, while the cavity mode is in the $\ket{\text{vac}}$ vacuum state. We
neglect spontaneous emission into free space as the cavity will provide a
dominant channel of energy loss by imposing that its decay rate $\kappa$ is
much larger than $g\sqrt{\mathcal{N}}$. The evolution of the system is then
described by

$\dot{\rho}_{\text{tot}}(t)=i\left[\rho_{\text{tot}}(t),\mathcal{H}\right]+\mathcal{L}_{\text{cav}}[\rho_{\text{tot}}].$
(93)

We notice first that the collapse operator $a$ only couples to the total spin
ladder operators $S$ and $S^{\dagger}$ which insures that the dynamics only
takes place within the symmetric manifold with maximal spin length
$\mathcal{N}/2$. The second observation is that, owing to the quick loss rate
$\kappa$ which is much larger than the coherent light-matter exchange rate
$g\sqrt{\mathcal{N}}$, the cavity field is permanently in a very low occupancy
state. This allows us to apply the standard master equation procedure (see
Appendix A1) where we assume the crude factorization
$\rho_{\text{tot}}(t)=\rho(t)\otimes\ket{\text{vac}}\bra{\text{vac}}$ at all
times. The next step is in evaluating two-time correlation terms with the new
time-dependent field operator $\mathcal{F}(t)=ga(t)$ (equivalent to the
expression in Eq. (5)). The situation is now much simpler than in the free-
space spontaneous emission case, as one can use the two-time correlations
$\braket{a(t^{\prime})a^{\dagger}(t^{\prime\prime})}=e^{-\kappa(t^{\prime}-t^{\prime\prime})}$
to find an effective damping rate $g^{2}/\kappa$. The dynamics can then be
reduced to the atomic subspace and incorporated in a master equation

$\displaystyle\dot{\rho}(t)=i[\rho(t),\mathcal{H}_{0}]+\mathcal{L}_{e}^{{}_{\text{cav}}}[\rho],$
(94)

with $\mathcal{H}_{0}=\omega_{0}\left[S_{z}+(\mathcal{N}/2)\mathbbm{1}\right]$
and the cavity-induced collective spontaneous emission

$\displaystyle\mathcal{L}_{e}^{{}_{\text{cav}}}[\rho]=\frac{g^{2}}{\kappa}\left(2S\rho
S^{\dagger}-S^{\dagger}S\rho-\rho S^{\dagger}S\right).$ (95)

The resulting dynamics is of course identical to the one described in Sec. D
showing the emergence of quick bursts of light as superradiant pulses but with
a rate established by the cavity $g^{2}/\kappa$. The effect is the spontaneous
synchronization of dipoles by the cavity to give rise to the superradiant
lasing regime. Here, the lossy cavity acts mainly as a communication bus to
drive the synchronization between the atoms and the average photon number
inside the cavity is typically well below one. For example, Ref. Bohnet _et
al._ (2012) reports the realization of a superradiant laser with an average
photon number of less than $0.2$ and a single-atom cooperativity
$\mathcal{C}=7.7\cdot 10^{-3}$, while the lasing regime achieved in Ref Norcia
_et al._ (2016) uses a MHz linewidth optical clock transition as its gain
medium and is envisioned to act as an active atomic clock insensitive to
fluctuations in reference cavity length (a fundamental limitation in
conventional lasers).

###  D Further remarks

Composite systems made up of flat, standard mirrors and two-dimensional single
emitter-thick regular arrays can find a multitude of applications in the
direction of quantum technologies. When used in a standard cavity QED
scenario, they can act as quick phase switchers with applications in precision
spectroscopy of quantum network characterization Plankensteiner _et al._
(2017) or in hybrid quantum optomechanical setups with enhanced photon-phonon
couplings Dantan _et al._ (2014). In the strong reflectivity limit, subradiant
arrays have been interfaced with two-dimensional semiconductor monolayers to
show enhanced quantum nonlinear optical properties Wild _et al._ (2018). The
identification of dark collective states between the two mirrors has been
proposed to allow for the preparation of Bell superpositions states between
the two subradiant layers Guimond _et al._ (2019), in a double-mirror hybrid
setup as discussed above. As such subradiant mirrors are of extremely small
mass, their zero-point motion is much larger than that of traditional
dielectric mirrors used in standard optomechanics: this opens new
opportunities for quantum optomechanics at the single photon-phonon level
Shahmoon _et al._ (2019, 2020).

##  VI Quantum optics with molecules

A relatively recent novel direction of research with great promise in the
direction of quantum technologies is the engineering of cavity-dressed
materials. In the field of molecular polaritonics, embedding organic
semiconductors in optical cavities allows for the design of novel materials
with enhanced properties such as exciton and charge transport Orgiu _et al._
(2015); Schachenmayer _et al._ (2015); Feist and Garcia-Vidal (2015);
Hagenmüller _et al._ (2017, 2018), superconductive behavior Sentef _et al._
(2018); Thomas _et al._ (2019) or modified chemical reactivity Hutchison _et
al._ (2012); Galego _et al._ (2017a); Herrera and Spano (2016); Martinez-
Martinez _et al._ (2018); Kampschulte and Denschlag (2018) due to light-
modified energy potential surfaces. This approach is also extended to study,
control and design phase transitions in quantum materials by quantum light
Kiffner _et al._ (2019); Wang _et al._ (2019b); Rohn _et al._ (2020); Sentef
_et al._ (2020); Ashida _et al._ (2020) etc.
Quantum emitters widely utilized in such experiments are of a much more
complex nature than the two-level system approximation. We exemplify here how
the previously introduced methods (master equation, quantum Langevin
equations) and concepts can be extended to additionally include couplings
between electronic transitions and vibrations or phonon modes. While this
approach is in principle amenable to a wide range of solid-state emitters
where phononic couplings play a crucial role, such as quantum dots or vacancy
centers in diamond, we focus here on the specific case of molecules. In
particular, organic molecules have emerged as tunable and efficient light-
matter interfaces due to their relatively large dipole moments, wide range of
transition frequencies and narrow linewidths at cryogenic temperatures
Toninelli _et al._ (2021). As such, single-molecule impurities embedded in
solid-state host matrices hold promise as single-photon sources Basché _et
al._ (1992); De Martini _et al._ (1996); Lounis and Moerner (2000); Pazzagli
_et al._ (2018) or nonlinear optical elements Maser _et al._ (2016). There
also is interest in exploiting the strong inherent coupling between electronic
transitions and the molecular vibrational degrees for the realization of
quantum optomechanical effects at the single molecule level Roelli _et al._
(2016); Neuman _et al._ (2019); Benz _et al._ (2016).
In this section, we will mainly follow the approach introduced in Ref. Reitz
_et al._ (2019) by considering the time evolution of a polaron operator (i.e.
an “effective”, vibrationally-dressed electronic dipole operator). We proceed
with a first-principle derivation of the Holstein Hamiltonian for electron-
vibron interactions and derive the standard Franck-Condon physics and
absorption/emission properties of molecules. We then illustrate how one can
circumvent a major downfall of molecular systems, i.e., they do not have
closed transitions, by employing the Purcell effect: a cavity-dressed molecule
can then behave as an ideal quantum emitter. Next we analyze molecular dimers
as good candidates for the observation and exploration of controllable near-
field couplings. Finally we introduce and analyze cooperative processes
occurring between near-field coupled disordered molecules leading to what is
known as the Förster resonance energy transfer.

###  A Optical response of a single molecule

Figure 14: Displaced oscillator model used to derive the Holstein Hamiltonian.
The Morse potential energy landscapes (solid gray) of electronic ground and
excited state are approximated as harmonic oscillators (dashed black) and
their eigenstates are illustrated in green. The position mismatch between
ground and excited state potential energy surfaces $R_{\text{ge}}$ gives rise
to vibronic coupling.

To include the electron-vibration (vibronic) coupling into our formalism, we
start with a first-principle derivation for a single electron coupled to a
single nuclear coordinate of reduced mass $\mu$. Notice that, for the simplest
case of a homonuclear diatomic molecule made of nuclei each with mass $m$, a
single vibrational mode exists corresponding to relative motion with effective
mass $\mu=m/2$; the equilibria coordinates then correspond to the bond length
in the ground and excited states. Under the Born-Oppenheimer approximation,
the ground and excited state potentials are obtained by solving the
Schrödinger equation for each given fixed internuclear distance. We assume
that, along the nuclear coordinate, the equilibrium positions for ground
(coordinate $R_{\text{g}}$, state vector $\ket{g}$) and excited (coordinate
$R_{\text{e}}$, state vector $\ket{e}$) electronic orbitals are slightly
shifted with respect to each other (as illustrated in Fig. 14). We can then
write the total Hamiltonian
$\mathcal{H}_{\text{mol}}=\mathcal{H}_{\text{el}}\otimes\mathcal{H}_{\text{vib}}$
as

$\mathcal{H}_{\text{mol}}=\mathcal{V}_{e}(\hat{R},\hat{P})\otimes\sigma^{\dagger}\sigma+\mathcal{V}_{g}(\hat{R},\hat{P})\otimes\sigma\sigma^{\dagger}.$

where the potential surfaces are described by the quantized position and
momentum operator satisfying $[\hat{R},\hat{P}]=i$. The dynamics takes place
now in an infinite-dimensional Hilbert space: in the electronic subspace the
algebra is that of a spin $1/2$ and is described by projectors
$\sigma^{\dagger}\sigma$ (into the excited state) and $\sigma\sigma^{\dagger}$
(into the ground state) while in the motional subspace one can for example
describe the dynamics in terms of plane waves. However, a great simplification
can be obtained by performing a harmonic expansion of the potential surfaces
around the minima

$\displaystyle\mathcal{V}_{e}(\hat{R},\hat{P})$
$\displaystyle=\omega_{0}+\frac{\hat{P}^{2}}{2\mu}+\frac{1}{2}\mu\nu^{2}\left(\hat{R}-R_{e}\right)^{2},$
(96a) $\displaystyle\mathcal{V}_{g}(\hat{R},\hat{P})$
$\displaystyle=\frac{\hat{P}^{2}}{2\mu}+\frac{1}{2}\mu\nu^{2}\left(\hat{R}-R_{g}\right)^{2}.$
(96b)

Finally, we can introduce small oscillations around the equilibria
$\hat{Q}=\hat{R}-R_{\text{g}}$ and subsequently
$\hat{R}-R_{\text{e}}=\hat{Q}+R_{\text{g}}-R_{\text{e}}=:\hat{Q}-R_{\text{ge}}$
leading to

$\displaystyle\mathcal{H}_{\text{mol}}$
$\displaystyle=\left(\omega_{0}+\frac{1}{2}\mu\nu^{2}R_{\text{ge}}^{2}\right)\mathds{1}_{\text{vib}}\otimes\sigma^{\dagger}\sigma$
(97)
$\displaystyle+\left(\frac{\hat{P}^{2}}{2\mu}+\frac{1}{2}\mu\nu^{2}\hat{Q}^{2}\right)\otimes\mathds{1}_{\text{el}}-\mu\nu^{2}R_{\text{ge}}\hat{Q}\otimes\sigma^{\dagger}\sigma.$
(98)

The last term is a renormalization of the bare electronic transition frequency
energy which will naturally go away when diagonalizing the Hamiltonian via the
polaron transformation resulting in $\omega_{0}$ as the natural electronic
transition frequency. We can now rewrite the momentum and position operators
in terms of bosonic operators $\hat{Q}=r_{\text{zpm}}(b^{\dagger}+b)$,
$\hat{P}=ip_{\text{zpm}}(b^{\dagger}-b)$ by introducing the zero-point motion
$r_{\text{zpm}}=1/\sqrt{2\mu\nu}$ and $p_{\text{zpm}}=\sqrt{\mu\nu/2}$.
Reexpressing the terms above via the factor $\sqrt{S}=\mu\nu
R_{\text{ge}}r_{\text{zpm}}$ ($S$ is the Huang-Rhys factor) yields the
Holstein-Hamiltonian Holstein (1959)

$\mathcal{H}_{\text{mol}}=(\omega_{0}+S\nu)\sigma^{\dagger}\sigma+\nu
b^{\dagger}b-\sqrt{S}\nu(b^{\dagger}+b)\sigma^{\dagger}\sigma,$ (99)

A molecular box illustration is shown in Fig. 15a as a minimal model for a
molecule.
One can bring this Hamiltonian into diagonal form via the polaron
transformation
$U_{\text{pol}}^{\dagger}=({\cal{D}^{\dagger}})^{\sigma^{\dagger}\sigma}=\ket{g}\bra{g}+{\cal{D}^{\dagger}}\ket{e}\bra{e}$
where the displacement is defined as
${\cal{D}}=\exp(-i\sqrt{2S}p)=\exp[\sqrt{S}(b^{\dagger}-b)]$ with the
dimensionless momentum quadrature $p=i(b^{\dagger}-b)/\sqrt{2}$ as generator.
The position quadrature is defined analogously as
$q=(b^{\dagger}+b)/\sqrt{2}$. The displacement operator creates a coherent
state with amplitude $\sqrt{S}$ when applied to the vibrational ground state
${\cal{D}}\ket{0}_{\text{vib}}=\ket{\sqrt{S}}_{\text{vib}}$. Note that the
vibrational creation and annihilation operators transform in the polaron
picture as

$\displaystyle U_{\text{pol}}^{\dagger}bU_{\text{pol}}$
$\displaystyle=b+\sqrt{S}\sigma^{\dagger}\sigma,$ (100a) $\displaystyle
U_{\text{pol}}^{\dagger}b^{\dagger}U_{\text{pol}}$
$\displaystyle=b^{\dagger}+\sqrt{S}\sigma^{\dagger}\sigma,$ (100b)

and furthermore the following commutation relations are fulfilled
$[p,{\cal{D}^{\dagger}}]=0$,
$[q,{\cal{D}^{\dagger}}]=-\sqrt{2S}{\cal{D}^{\dagger}}$. The Hamiltonian in
the diagonal basis then becomes
$\tilde{\mathcal{H}}_{\text{mol}}=\omega_{0}\sigma^{\dagger}\sigma+\nu
b^{\dagger}b$ and has eigenvectors $\\{\ket{g,n},\ket{e,n}\\}$ while the
eigenvectors in the original basis can be obtained by an inverse polaron
transform and are given by $\\{\ket{g,n}$, ${\cal{D}}\ket{e,n}\\}$.

Figure 15: (a) The dynamics of a molecule with a single nuclear coordinate can
be seen as the coupling of an electronic transition operator $\sigma$ to a
single vibrational mode $b$ with a strength given by the Huang-Rhys factor
$S$. (b) For many vibrational modes, a molecule can be constructed by coupling
an electronic transition operator $\sigma$ to $n$ vibrational modes with
strengths $S_{k}$, frequencies $\nu_{k}$ and linewidths $\Gamma_{k}$ described
by the (Lorentzian) spectral density function $J(\omega)$. (c) Energy
eigenstates of the electron-vibron system denoted by $\ket{g,n}$ and
$\ket{e,n}$. The zero phonon line is the $\ket{g,0}$ to $\ket{e,0}$
transition. For small $S<1$ the width of the arrows indicate the strength of
the optical transitions. For a molecule with many vibrational levels,
spontaneous emission outside the zero phonon line at summed rate
$\gamma^{\prime}$ can be larger than into it at rate $\gamma_{00}$. (d) The
enhancement of spontaneous emission via resonant coupling to an optical cavity
can lead to a great modification of the strength of the zero phonon line thus
turning a molecule into a closed quantum system.

Franck-Condon physics \- Let us now ask what is the optical response of the
molecule to a weak laser excitation modelled via
$\mathcal{H}_{\ell}=i\eta(\sigma^{\dagger}e^{-i\omega_{\ell}t}-\sigma
e^{i\omega_{\ell}t})$. In an interaction picture, reached after performing a
transformation of the Hamiltonian with
$e^{i\omega_{\ell}\sigma^{\dagger}\sigma t}$, one has

$\displaystyle\tilde{\mathcal{H}}=(\omega_{0}-\omega_{\ell})\sigma^{\dagger}\sigma+\nu
b^{\dagger}b+i\eta\left(\sigma^{\dagger}{\cal{D}^{\dagger}}-\sigma{\cal{D}}\right).$
(101)

Starting with the molecule in the absolute ground state $\ket{g,0}$, the
result of the drive is to excite the electron to state $\ket{e}$ while
exciting motion on the nuclear coordinate to a coherent motional state. One
can immediately compute the probability of absorption

$\displaystyle
P_{\text{abs}}=|\braket{e,n}{\tilde{\mathcal{H}}}{g,0}|^{2}=\braket{e,n}{\sigma^{\dagger}{\cal{D}^{\dagger}}}{g,0}|^{2}=f^{n}_{S},$
(102)

showing a Poissonian distribution $f^{n}_{S}=e^{-S}S^{n}/n!$ of the
occupancies of number states $\ket{n}$ with coefficients known as Franck-
Condon factors. In particular, the so-called zero-phonon line (transition
corresponding to $n=0$) is reduced by a factor $e^{-S}$ at zero temperature.
The Huang-Rhys factor $S$ should be interpreted as the average number of
vibrational quanta created upon excitation or deexciation of the molecule. The
probability of the emission process is similarly computed between $\ket{e,0}$
and final states $\ket{g,n}$ and yields the same Poissonian distribution in
$n$. This is a known result in molecular spectroscopy where the absorption and
emission spectra are mirror images of each other. For small $S<1$ the physics
of absorption and emission in a molecule is illustrated in Fig. 15c where the
width of the arrows indicate the strength of the optical transitions between
eigenstates of the free Hamiltonian (without the Holstein interaction).
The Holstein Hamiltonian in Eq. (99) can be generalized to a more complex
molecular box (is illustrated in Fig. 15b) involving $n_{v}$ independent
nuclear coordinates with frequencies $\nu_{k}$ and coupling strengths
$\sqrt{S_{k}}$. The total Hamiltonian is then expressed as

$\displaystyle\mathcal{H}_{\text{mol}}=\tilde{\omega}_{0}\sigma^{\dagger}\sigma+\sum_{k=1}^{n_{v}}\nu_{k}b_{k}^{\dagger}b_{k}-\sum_{k=1}^{n_{v}}\sqrt{S}_{k}\nu_{k}(b_{k}^{\dagger}+b_{k})\sigma^{\dagger}\sigma,$
(103)

where the accumulated frequency shift is
$\tilde{\omega}_{0}=\omega_{0}+\sum_{k=1}^{n_{v}}S_{k}\nu_{k}$.
Diagonalization of this Hamiltonian is straightforward as it involves a
collective displacement for all nuclear coordinates
$U_{\text{pol}}^{\dagger}=(\prod_{k}\mathcal{D}_{k}^{\dagger})^{\sigma^{\dagger}\sigma}$
with $\mathcal{D}_{k}=\exp(-i\sqrt{2S_{k}}p_{k})$. The expectation value of
the displacement operator
$\braket{{\cal{D}}_{k}}=\braket{{\cal{D}}_{k}^{\dagger}}=e^{-S_{k}\braket{p_{k}^{2}}}$
(assuming the vibrational mode to be in a thermal state) can be expressed as
$e^{-S_{k}(\bar{n}_{k}+1/2)}$ (expanding the displacement operator and making
use of the fact that Gaussian states are characterized by their second-order
moments and odd moments are vanishing) where the average thermal occupancy is
given by $2\bar{n}_{k}=\coth(\beta\nu_{k}/2)-1$ with $\beta=1/(k_{B}T)$.

Branching ratio manipulation \- The branching ratio quantifies the rate of
emission into the zero-phonon line versus all other lines. For ideal quantum
emitters this is automatically unity in the absence of coupling to any other
channels of de-excitation. For an isolated molecule this is instead given by
$\alpha=\braket{{\cal{D}}}^{2}=\prod_{k}\braket{{\cal{D}}_{k}}^{2}=e^{-\sum_{k}S_{k}}$
(at zero temperature) stemming from loss of excitation into other vibrational
levels. Complex molecules have many vibrations: even for small Huang-Rhys
factors the sum in the exponential gives a considerable reduction of the
oscillator strength of the zero-phonon line. The Purcell effect (introduced in
Sec. A), stemming from the resonant interaction of the zero-phonon line
transition with a lossy confined optical resonance has been experimentally
proven to turn molecular emitters into almost perfect two-level quantum
emitters Wang _et al._ (2019a). Assuming that the cavity linewidth is narrow
compared to the vibrational frequency $\kappa\ll\nu$ and only couples to the
zero-phonon line transition, we can define by $g_{00}=g\braket{{\cal{D}}}$ the
effective reduced coupling between the cavity and the zero-phonon line with a
corresponding cooperativity $\mathcal{C}_{00}=g_{00}^{2}/(\kappa\gamma)$. This
will lead to an enlarged zero-phonon transition linewidth
$\tilde{\gamma}_{00}=\gamma_{00}(1+\mathcal{C}_{00})$ as illustrated in Fig.
15 while other transitions are unperturbed and sum up to
$\gamma^{\prime}=\gamma-\gamma_{00}$. The cavity-modified branching ratio can
then be expressed as

$\displaystyle\alpha_{\text{cav}}=\frac{\tilde{\gamma}_{00}}{\tilde{\gamma}_{00}+\gamma^{\prime}}=\frac{(1+\mathcal{C}_{00})e^{-\sum_{k}S_{k}}}{1+\mathcal{C}_{00}e^{-\sum_{k}S_{k}}},$
(104)

which reproduces the bare molecule result for $\mathcal{C}_{00}=0$ and goes to
unity in the limit that $\mathcal{C}_{00}\to\infty$. Experimentally, values of
$\mathcal{C}_{00}$ close to $100$ have been reached Wang _et al._ (2019a),
showing that almost unity branching ratios can be achieved.

Quantum Langevin equations for polarons \- We will now proceed by deriving a
quantum Langevin equation for a polaron operator
$\tilde{\sigma}=\sigma{\cal{D}^{\dagger}}$ (i.e., a vibrationally dressed
electronic dipole operator). This will allow us to include all coherent and
incoherent effects at the level of the equations of motion. The collective
displacement operator ${\cal{D}}^{\dagger}$ can refer to a large set of
molecular vibrational modes
${\cal{D}^{\dagger}}=\prod_{k=1}^{n_{v}}{\cal{D}}_{k}^{\dagger}$. The method
can also be extended to include coupling to an (eventually infinite) set of
phonon modes of the host environment Reitz _et al._ (2020); Clear _et al._
(2020). We will assume a Brownian noise dissipation model (see Appendix A5)
for the molecular vibrational modes. Starting with the Holstein Hamiltonian
one can derive the following equations of motion for vibrations

$\displaystyle\dot{q}_{k}$ $\displaystyle=\nu_{k}p_{k},$ (105a)
$\displaystyle\dot{p}_{k}$
$\displaystyle=-2\Gamma_{k}p_{k}-\nu_{k}q_{k}-\sqrt{2S_{k}}\nu_{k}\sigma^{\dagger}\sigma+\xi_{k}.$
(105b)

In addition to the environment induced damping at rates $2\Gamma_{k}$ and
associated thermal input noise $\xi_{k}$ (see Appendix A5 for noise
properties), vibrations are driven by their coupling to the electronic degree
of freedom. For the electronic transition, the polaron operator can be
computed from
$\dot{\tilde{\sigma}}=\dot{\sigma}\mathcal{D}^{\dagger}+\sigma\dot{\mathcal{D}}^{\dagger}$
where the full Langevin equation for $\sigma$ is derived from the Holstein
Hamiltonian and the Linbdlad form radiative emission with the transformation
indicated by the relation in Eq. (63). We obtain then

$\displaystyle\dot{\tilde{\sigma}}=$
$\displaystyle-[\gamma-i(\omega_{\ell}-\omega_{0})]\tilde{\sigma}-2i\sum_{k}\sqrt{2S_{k}}\Gamma_{k}p_{k}\tilde{\sigma}$
(106)
$\displaystyle+\eta{\cal{D}^{\dagger}}+\sqrt{2\gamma}\tilde{\sigma}_{\text{in}}+i\sum_{k}\sqrt{2S_{k}}\xi_{k}\tilde{\sigma}.$

The equation is derived under the assumption of weak driving $\eta\ll\gamma$
and therefore small occupancy of the electronic excited state
$\braket{\sigma^{\dagger}\sigma}\ll 1$ and with special attention devoted to
the rules for taking the time derivative of any exponential operator
$e^{A(t)}$ with $[\dot{A}(t),A(t)]\neq 0$. Moreover, some of the terms cancel
(see Ref. Reitz _et al._ (2019) for full details) and an extremely simple
equation of motion for the polaron operator is obtained

$\displaystyle\dot{\tilde{\sigma}}=-[\gamma-i(\omega_{\ell}-\omega_{0})]\tilde{\sigma}+\sqrt{2\gamma}\tilde{\Sigma}_{\text{in}}.$
(107)

Here we have combined the input fields affecting the electronic transition
into the expression
${\Sigma}_{\text{in}}=\eta/\sqrt{2\gamma}+\sigma_{\text{in}}$ and the
displaced input field (drive plus noise) is
$\tilde{\Sigma}_{\text{in}}={\Sigma}_{\text{in}}{\cal{D}^{\dagger}}$. A
solution of the bare electronic dipole operators can be obtained by a formal
integration as a sum $\sigma(t)=\sigma_{\text{tr}}(t)+\sigma_{\text{ss}}(t)$
between transient and steady state solutions

$\displaystyle\sigma_{\text{tr}}(t)$
$\displaystyle=\sigma(0){\cal{D}}(t){\cal{D}^{\dagger}}(0)e^{-\left[\gamma-i(\omega_{\ell}-\omega_{0})\right]t},$
(108a) $\displaystyle\sigma_{\text{ss}}(t)$
$\displaystyle=\sqrt{2\gamma}\int_{0}^{t}dt^{\prime}e^{-\left[\gamma-i(\omega_{\ell}-\omega_{0})\right](t-t^{\prime})}\Sigma_{\text{{in}}}(t^{\prime}){\cal{D}}(t){\cal{D}^{\dagger}}(t^{\prime}).$
(108b)

To derive the optical response of the electronic degree of freedom in the
presence of its coupling to vibrations one then has to evaluate two-time
correlation functions of the displacement operator. This is factorizable as we
consider that all vibrational modes are independent from each other, i.e.,
$\braket{\mathcal{D}(t){\cal{D}^{\dagger}}(t^{\prime})}=\prod_{k}\braket{\mathcal{D}_{k}(t){\cal{D}^{\dagger}}_{k}(t^{\prime})}$.
By expanding the displacement operators and applying the Isserlis’ theorem (or
Wick’s probability theorem for multivariate normal distributions) one obtains

$\displaystyle\braket{\mathcal{D}_{k}(t){\cal{D}}_{k}^{\dagger}(t^{\prime})}=e^{-2S_{k}(\braket{p_{k}^{2}}-\braket{p_{k}(t)p_{k}(t^{\prime})})},$
(109)

where we made use of the fact that the momentum variance is stationary, i.e.,
$\braket{p_{k}^{2}}=\braket{p_{k}(t)^{2}}=1/2+\bar{n}_{k}$. This can be easily
computed by using the properties for the Brownian noise model (see Appendix
A5)

$\displaystyle\braket{p_{k}(t)p_{k}(t^{\prime})}=\left[\left(\bar{n}_{k}+\frac{1}{2}\right)\cos(\nu_{k}\tau)-\frac{i}{2}\sin(\nu_{k}\tau)\right]e^{-\Gamma_{k}\tau},$
(110)

with $\tau=t-t^{\prime}$. While expressions can be derived for any temperature
of the bath, in the following we will, for the sake of simplicity, focus on
the zero temperature case for all vibrational modes. This is typically a very
good approximation as molecular vibrations typically lie in the range of THz
frequencies and are therefore barely occupied as the condition
$\hbar\nu_{k}\ll k_{B}T$ is typically fulfilled. The displacement correlation
function in this case simplifies to

$\displaystyle\braket{\mathcal{D}_{k}(t){\cal{D}}_{k}^{\dagger}(t^{\prime})}=e^{-S_{k}}e^{S_{k}\exp\left[-(\Gamma_{k}+i\nu_{k})\tau\right]},$
(111)

which we will make use of in the next subsection to derive analytical
expressions for absorption and emission lines and linewidths.

Absorption and emission \- The absorption and emission spectra of a driven
molecule can be readily obtained from Eq. (108). Considering large times $t\gg
1/\gamma$, the first term in Eq. (108) goes towards zero and the average
dipole moment can be expressed as (taking the average over both electronic and
vibrational degrees of freedom)

$\displaystyle\braket{\sigma(t)}=\eta\int_{0}^{t}dt^{\prime}e^{-(\gamma-i\Delta)(t-t^{\prime})}\braket{{\cal{D}}(t){\cal{D}^{\dagger}}(t^{\prime})},$
(112)

with the laser detuning $\Delta=\omega_{\ell}-\omega_{0}$. The population of
the electronic excited state evolves according to
$\partial_{t}\braket{{\sigma^{\dagger}\sigma}}=-2\gamma\braket{\sigma^{\dagger}\sigma}+2\eta\Re\braket{\sigma(t)}$
meaning that the steady state absorption profile can be directly computed from
the expectation value of the coherence $\braket{\sigma(t)}$. We define
$S_{\text{abs}}(\Delta)=\braket{\sigma^{\dagger}\sigma}$ as the absorption
profile for a varying laser frequency with expression

$\displaystyle
S_{\text{abs}}(\Delta_{\ell})=2\eta\int_{-\infty}^{t}dt^{\prime}e^{-2\gamma(t-t^{\prime})}\Re\braket{\sigma(t^{\prime})}.$
(113)

To evaluate this integral, we expand the displacement correlation function
from Eq. (111) in a power series

$\displaystyle\braket{{\cal{D}}_{k}(t){\cal{D}}_{k}^{\dagger}(t^{\prime})}=e^{-S_{k}}\sum_{n_{k}=0}^{\infty}\frac{S_{k}^{n_{k}}}{n_{k}!}e^{-n_{k}(\Gamma_{k}+i\nu_{k})(t-t^{\prime})}.$
(114)

The absorption spectrum can then be readily expressed as a sum over all
indices $\\{n_{k}\\}=n_{1},\ldots,n_{n_{v}}$ weighted by the Franck-Condon
factors $f_{n_{k}}^{S_{k}}$

$\displaystyle
S_{\text{abs}}(\Delta)=\frac{\eta^{2}}{\gamma}\sum_{\\{n_{k}\\}=0}^{\infty}\frac{f^{n_{1}}_{S_{1}}\ldots
f^{n_{n_{v}}}_{S_{n_{v}}}\left[\gamma+\Gamma_{\\{n_{k}\\}}\right]}{\left[\gamma+\Gamma_{\\{n_{k}\\}}\right]^{2}+\left[\Delta-\nu_{\\{n_{k}\\}}\right]^{2}},$
(115)

where the denominator indicates a series of blue-shifted absorption sidebands
with increased linewidths
$\Gamma_{\\{n_{k}\\}}=\sum_{\\{n_{k}\\}}n_{k}\Gamma_{k}$ at positions
$\omega_{\ell}=\omega+\sum_{k}\nu_{\\{n_{k}\\}}$ with
$\nu_{\\{n_{k}\\}}=\sum_{\\{n_{k}\\}}n_{k}\nu_{k}$.
For the definition of the emission spectrum one can consider the transient
dynamics of a molecule which is partially excited at time $t=0$:
$\braket{\sigma^{\dagger}\sigma(0)}=p_{0}$. The emission spectrum is then
defined as the Fourier transform
$S_{\text{em}}(\omega)=2\Re\int_{0}^{\infty}d\tau\braket{\mathcal{T}\left\\{\sigma^{\dagger}(\tau)\sigma(0)\right\\}}e^{-i\omega\tau}$
(where $\mathcal{T}$ denotes time ordering) and can be computed as (normalized
by $p_{0}$)

$\displaystyle
S_{\text{em}}(\omega)=\sum_{\\{n_{k}\\}=0}^{\infty}\frac{2f^{n_{1}}_{S_{1}}\ldots
f^{n_{n_{v}}}_{S_{n_{v}}}\left[\gamma+\Gamma_{\\{n_{k}\\}}\right]}{\left[\gamma+\Gamma_{\\{n_{k}\\}}\right]^{2}+\left[\omega-\omega_{0}+\nu_{\\{n_{k}\\}}\right]^{2}}.$
(116)

The denominator indicates a series of red-shifted Stokes lines at frequencies
$\omega=\omega_{0}-\nu_{\\{n_{k}\\}}$.

###  B Near-field coupled molecules

Let us now consider two electronically coupled molecules. We start by
discussing the vibronic dimer, a well-known model to describe the interplay
between electronic and vibrational interactions in molecular aggregates which
has been studied both theoretically Witkowski and Moffitt (1960); Riley and
Krausz (1990); Eisfeld _et al._ (2005) and experimentally Lippitz _et al._
(2004); Diehl _et al._ (2014). As opposed to the pure electronic dimer
discussed in Sec. B, this model gives rise to a more complex energy landscape
with states possessing both electronic and vibrational character. Effects
showing quantum coherence in electronically excited dimers have been
experimentally investigated in a combination of single molecule spectroscopy
and quantum chemistry techniques Diehl _et al._ (2014). We then move on to
describe the dipole-dipole mediated unidirectional energy transfer process
(FRET) between two molecules.

Vibronic dimer model \- The vibronic dimer model consists of two molecules
(indexed here by $1$ and $2$), each with a single electronic ground and
excited state and each coupled to a single harmonic nuclear coordinate
$q_{1/2}$ with Huang-Rhys factor $S$. The molecules are coupled to each other
electronically (e.g. via dipole-dipole interaction) with strength $\Omega$.
Assuming that the excited state energies of the two molecules are identical
and using the single excitation manifold with states
$\ket{e}_{1}\otimes\ket{g}_{2}$ and $\ket{g}_{1}\otimes\ket{e}_{2}$, the
Hamiltonian of the dimer can be represented by

Figure 16: Potential energy surfaces $\epsilon_{\pm}$ of the excited state
manifold of the vibronic dimer model as a function of the (dimensionless)
relative motion coordinate $Q_{-}$ for $\Omega=0$ (red dashed),
$2\Omega=\Delta\epsilon/4$ (red solid) and $2\Omega=\Delta\epsilon$ (orange).

$\displaystyle\mathcal{H}_{\text{dim}}\\!=\\!\begin{pmatrix}-\sqrt{2S}\nu
q_{1}&\Omega\\\ \Omega&-\sqrt{2S}\nu
q_{2}\end{pmatrix}+\frac{\nu}{2}\sum_{i=1,2}\left(p_{i}^{2}+q_{i}^{2}\right)\mathbb{1},$
(117)

where $\mathbb{1}$ is the $2\times 2$ unity matrix. Introducing symmetrized
coordinates $Q_{\pm}=(q_{1}\pm q_{2})/\sqrt{2}$ and electronic operators
$\sigma_{\pm}=(\sigma_{1}\pm\sigma_{2})/\sqrt{2}$, the Hamiltonian transforms
to

$\displaystyle\tilde{\mathcal{H}}_{\text{dim}}$
$\displaystyle=\begin{pmatrix}\Omega&-\sqrt{S}\nu Q_{-}\\\ -\sqrt{S}\nu
Q_{-}&-\Omega\end{pmatrix}-\sqrt{S}\nu Q_{+}\mathbb{1}$ (118)
$\displaystyle+\frac{\nu}{2}\left(P_{+}^{2}+P_{-}^{2}+Q_{+}^{2}+Q_{-}^{2}\right)\mathbb{1},$

from which one can see that only the coordinate $Q_{-}$ couples between
electronic states in the collective basis. The coordinate $Q_{+}$ instead
leads to a constant shift and can eventually be eliminated by a proper choice
of the origin. The Hamiltonian in Eq. (118) can be easily diagonalized leading
to the eigenvalues for the excited state potential energy surfaces

$\displaystyle\epsilon_{\pm}=\frac{\nu}{2}\left(Q_{+}^{2}+Q_{-}^{2}\right)-\sqrt{S}\nu
Q_{+}\pm\sqrt{S\nu^{2}Q_{-}^{2}+\Omega^{2}}.$ (119)

The potential energy surfaces are plotted in Fig. 16 for various coupling
strengths $\Omega$. One can see that the coupling lifts the degeneracy of the
potential surfaces at $Q_{-}=0$ and leads to a splitting of $2\Omega$.
Depending on $\Omega$, the lower energy surface either has two minima or just
a single minimum. An important quantity is the reorganization energy
$\Delta\epsilon=2S\nu$ defined as the vertical energy separation between the
potential curves at each minimum for zero coupling $\Omega=0$ [as illustrated
in Fig. 16]. One finds that a double minimum of the lower energy surface
located at $Q_{-}^{\pm}=\pm\sqrt{S[1-\Omega^{2}/(S^{2}\nu^{2})]}$ exists when
the level splitting $2\Omega$ is smaller than the reorganisation energy
$|2\Omega|<\Delta\epsilon$.

Förster resonance energy transfer \- Förster resonance energy transfer (FRET)
is the main excitation transfer mechanism in photosynthetic light harvesting
complexes, where energy is transfered with high efficiency from a light-
harvesting antenna through a network of molecules (chromophores) to a reaction
centre Förster (1948); Blankenship (2014). Although a lot of progress has been
made in understanding and unravelling the reasons for the high efficiency of
the photosynthetic energy transfer occuring in nature over the last years,
there is still ongoing debate e.g. about the role of quantum coherence or
entanglement in the energy transfer process Engel _et al._ (2007); Streltsov
_et al._ (2017). In the sense of this tutorial, FRET can be seen as a truly
cooperative process since it involves a complex interplay between electronic
and nuclear degrees of freedom as well as both coherent and incoherent
couplings. Generally speaking, the requirement for FRET is a spectral overlap
between the emission profile of the donor molecule and the absorption profile
of the acceptor, giving rise to a resonant transfer process typically
accompanied by quick vibrational relaxation May and Kühn (2011).
Here we exemplify the application of the QLEs formalism to derive a
perturbative expression for the FRET rate between two near-field coupled
molecules $D$ and $A$ with energy mismatch $\Delta=\omega_{D}-\omega_{A}$,
first by considering a simple configuration where each molecule has a single
vibrational mode $b_{D}$ and $b_{A}$ (with corresponding vibrational
frequencies $\nu_{D}$ and $\nu_{A}$ and Huang-Rhys factors $S_{D}$ and
$S_{A}$). We will assume an initially excited donor molecule which can undergo
a resonant exchange with the acceptor’s excited state vibrational manifold. In
a standard scenario, a unidirectional process emerges as the acceptor quickly
relaxes to its vibrational ground state $\Gamma_{D}\gg\Omega$. We start with
the equations of motion for the polaron-transformed dipole operators
$\tilde{\sigma}_{D}=\mathcal{D}_{D}^{\dagger}\sigma_{D}$ and
$\tilde{\sigma}_{A}=\mathcal{D}_{A}^{\dagger}\sigma_{A}$ with the displacement
operators $\mathcal{D}_{D}=e^{\sqrt{S_{D}}(b_{D}^{\dagger}-b_{D})}$ and
$\mathcal{D}_{A}=e^{\sqrt{S_{A}}(b_{A}^{\dagger}-b_{A})}$ which shows coherent
coupling induced by the dipole-dipole near-field exchange of virtual photons

$\displaystyle\dot{\tilde{\sigma}}_{D}=-\left(\gamma_{D}+i\omega_{D}\right)\tilde{\sigma}_{D}-i\Omega\tilde{\sigma}_{A}\mathcal{D}_{A}\mathcal{D}_{D}^{\dagger}+\sqrt{2\gamma}\tilde{\sigma}_{D}^{\text{in}},$
(120a)
$\displaystyle\dot{\tilde{\sigma}}_{A}=-\left(\gamma_{A}+i\omega_{A}\right)\tilde{\sigma}_{A}-i\Omega\tilde{\sigma}_{D}\mathcal{D}_{D}\mathcal{D}_{A}^{\dagger}+\sqrt{2\gamma}\tilde{\sigma}_{A}^{\text{in}}.$
(120b)

To derive a perturbative expression for the energy transfer rate, we consider
a scenario with an initial full excitation of $D$ and no excitation of $A$.
Then from above one can derive an equation of motion for the acceptor’s
population

$\displaystyle\dot{P}_{A}=-2\gamma_{A}P_{A}+2\Omega\Im\braket{\sigma_{A}^{\dagger}\sigma_{D}}+\sqrt{2\gamma_{A}}\braket{\sigma_{A}^{\dagger}\sigma_{A}^{\text{in}}+\sigma_{A}^{\dagger,\text{in}}\sigma_{A}},$
(121)

which shows how the two-particle correlations
$\braket{\sigma_{A}^{\dagger}\sigma_{D}}$ appear as source terms for the
acceptor population. The procedure to compute
$2\Omega\Im\braket{\sigma_{A}^{\dagger}\sigma_{D}}$ involves then formal
integration of the equation of motion for the acceptor polaron operator (as
detailed in the Appendix A6) under the assumption that
$\Gamma_{D},\Gamma_{A}\gg\gamma_{D},\gamma_{A}$. Finally, one arrives at an
expression for the energy transfer rate showing proportionality to the donor
population
$2\Omega\braket{\sigma_{A}^{\dagger}\sigma_{D}}=\kappa_{\text{ET}}P_{D}(t)$

$\displaystyle\kappa_{\text{ET}}=\sum_{n_{D},n_{A}}\frac{2\Omega^{2}f_{S_{D}}^{n_{D}}f_{S_{A}}^{n_{A}}(n_{D}\Gamma_{D}+n_{A}\Gamma_{A})}{(n_{D}\Gamma_{D}+n_{A}\Gamma_{A})^{2}+[\Delta-
n_{D}\nu_{D}-n_{A}\nu_{A}]^{2}},$ (122)

where, for simplicity, we have set $\gamma_{D}=\gamma_{A}$. The denominator of
the above expression asks that the resonance condition
$\omega_{D}-n_{D}\nu_{D}=\omega_{A}+n_{A}\nu_{A}$ is fulfilled, i.e.,
fluorescence lines of the donor $\ket{e_{D},0}\to\ket{g_{D},n_{D}}$ have to
overlap with absorption lines of the acceptor
$\ket{g_{A},0}\to\ket{e_{A},n_{A}}$, both of which are weighted by the
respective Franck-Condon factors.
In the case of many vibrational modes $n$ for donor and acceptor, we can
generalize the result by writing general displacements
$\mathcal{D}_{A}=\prod_{k=1}^{n_{v}}\mathcal{D}_{A}^{k}$ and
$\mathcal{D}_{D}=\prod_{k=1}^{n_{v}}\mathcal{D}_{D}^{k}$ for all vibrational
modes. The rate is computed by assuming multiple paths of energy transfer
between the two molecules involving all vibrational modes. We assume an
initially electronically excited state with no vibrations present
$\ket{e_{D};0_{1},0_{2}...0_{n_{v}}}$ of molecule $D$ and ground state without
vibrations $\ket{g_{A};0_{1},0_{2}...0_{n_{v}}}$ for molecule $A$. The
emission of molecule $D$ leads it into state
$\ket{g_{D};m_{1},...m_{k},...m_{n_{v}}}$ and resonant interactions can occur
with state $\ket{e_{A};l_{1},...l_{k^{\prime}},...l_{n_{v}}}$ of molecule $A$.
Summing over all these processes leads to an analytical cumbersome expression
for $\kappa_{\text{ET}}$ (listed in Appendix A6) seen simply as a
generalization of the energy transfer rate from above. Such an analytical
result is the discrete version of the well established integral formulation
May and Kühn (2011) describing the overlap between the emission spectrum of
molecule $D$ and absorption spectrum of molecule $A$ as illustrated in Fig.
17a for a donor-acceptor pair with $8$ vibrational modes and spectral density
$J(\omega)=\sum_{k}2S_{k}\nu_{k}^{2}\Gamma_{k}/[\Gamma_{k}^{2}+(\omega-\nu_{k})^{2}]$
(assumed to be identical for both molecules) stemming from the coupling of the
molecular vibrations to some external phonon bath allowing for vibrational
relaxation at rates $\Gamma_{k}$. The process is unidirectional as shown in
Fig. 17b as $\omega_{D}-\omega_{A}$ dictates the direction of the energy flow.

Figure 17: Energy transfer in donor-acceptor configurations. (a) Overlap
between acceptor emission spectrum and donor absorption spectrum leads to the
FRET process. The results assume a given vibrational spectral energy
$J(\omega)$ (assumed here to be identical for both molecules) shown in the
inset, for $8$ vibrational modes considered. (b) Numerical analysis of the
energy transfer rate (for $\Omega=30\gamma$) shows a unidirectional flow from
the higher energy molecule (donor) to the lower energy one (acceptor)
optimized around a detuning corresponding to the overlap of the donor
electronic state to a high density of vibrational states in the acceptor
excited state manifold.

###  C Further remarks

We have covered here the electron-vibron interaction and described the optical
response to a weak classical optical drive. We have however hinted that
coupling to confined optical modes can lead to strong modifications on the
matter side, such as the described model of a molecule turned into an
efficient quantum emitter Wang _et al._ (2019a). Indeed, many experimental and
theoretical efforts have been recently targeted into the direction of cavity
charge and energy transport enhancement Orgiu _et al._ (2015); Schachenmayer
_et al._ (2015); Feist and Garcia-Vidal (2015); Hagenmüller _et al._ (2017,
2018), Förster resonance energy transfer (FRET) enhancement Zhong _et al._
(2016, 2017); Garcia-Vidal and Feist (2017); Du _et al._ (2018); Reitz _et
al._ (2018), modified chemical reactivity Hutchison _et al._ (2012); Galego
_et al._ (2017a); Herrera and Spano (2016); Martinez-Martinez _et al._ (2018);
Kampschulte and Denschlag (2018), polariton dynamics Kansanen _et al._ (2019);
Schwartz _et al._ ; Canaguier-Durand _et al._ (2015), etc. On the theory side
a large effort is aimed at translating standard cavity QED concepts (as
introduced in Sec. A) with two-level systems into the realm of molecules. To
this end, theoretical efforts consider a generalized light-electronic-
vibrations problem modeled as a Holstein-Tavis-Cummings Hamiltonian.
Investigations aim at providing an understanding of the vibrationally induced
cavity polariton asymmetry Canaguier-Durand _et al._ (2015); Neuman and
Aizpurua (2018), vibrationally dressed polaritons Zeb _et al._ (2018), dark
vibronic polaritons Herrera and Spano (2017, 2018), developing a cavity Born-
Oppenheimer theory Flick _et al._ (2017a, b) or deriving relevant simplified
models for large scale numerics in the mesoscopic limit del Pino _et al._
(2018).
We also remark that molecules are natural quantum mechanical platforms
interfacing electronic and vibrational degrees of freedom via an interaction
that resembles the radiation pressure Hamiltonian (via a boson-spin
replacement) utilized in quantum optomechanics; as vibration at terahertz
frequencies are practically in the ground state even at room temperature and
the strength of the coherent coupling can be comparable to the vibrational
frequency, molecules could be a good platform for studying quantum state
transfer between light and motion Benz _et al._ (2016). Moreover, molecules
are good single photon emitters as they are characterized by large emission
rates and they allow for highly efficient collection schemes, thus utilizable
in applications for photon-photon entanglement generation and optical quantum
computing Lettow _et al._ (2010).

##  VII Conclusions

Quantum systems, comprised of many particles, inherently coupled via common
reservoirs (the electromagnetic vacuum, optical resonators, optical
waveguides, etc.) can exhibit cooperative behavior with some specific
properties scaling more favorably than what would be expected from a
collection of uncoupled particles. In order to analytically and numerically
understand the cooperativity of light-matter platforms, we have detailed two
alternative approaches based on either the time evolution of the density
operator (master equation approach) or of system operators (stochastic quantum
Langevin equations). Based on this toolbox, this tutorial introduced a set of
equations successfully applied to chains, rings and arrays of quantum emitters
or ensembles of molecules within the confined mode of optical resonators.
A first set of applications that we have tackled are based on effects stemming
from the dipole-dipole interaction occurring in dense quantum emitter
ensembles. This interaction allows for the hopping of excitations in one or
two-dimensional regular arrays and for the engineering of energy bands with
special properties such as band gaps or Dirac points. Thus, such systems are
ideal platforms for the investigation of topological quantum photonics with
built-in nonlinearities. Their subradiant properties allow for applications in
quantum metrology, in the design of high-fidelity photon-storage platforms and
in the generation of entangled Bell states of two or more quantum emitter
arrays (hinting towards quantum networks of spatially distant subwavelength
arrays). Their enhanced reflectivity properties render them useful as
extremely light mechanical resonators for quantum nano-optomechanical
applications or as metasurfaces with enhanced optical nonlinearities at the
level of a single photon. Also, subradiant rings illuminated by incoherent
light can provide a natural gain medium via supported waveguide resonances
leading to thresholdless nanolasers.
A second set of applications that we have discussed are based on the
enhancement of the photon-emitter coupling in cavity QED. We have demonstrated
that the light-matter cooperativity can be greatly enhanced by addressing
collective subradiant states of one or two-dimensional arrays which are
characterized by narrow antiresonances in transmission. A hybrid cavity
approach, where flat mirrors are replaced with frequency-dependent mirrors
(such as subwavelength arrays or photonic crystal mirrors) can lead to narrow
cavity resonances useful for example in resolved sideband optomechanics. We
provided an input-output theory, valid outside Markovian regimes with wide
applicability. In the same context of cavity QED, we provided a simple theory
of superradiant lasers with potential applications in enhancing the stability
of atomic clocks.
As molecular systems, quantum dots, vacancy centers in solid-state hosts are
increasingly important for quantum technology applications, it is crucial to
understand and to tailor the electron-vibration interaction. We provided here
a first-principle derivation of the Holstein Hamiltonian for electron-vibron
interactions in simple molecular systems and extended the QLEs approach to
polaron physics. This theory leads to very simple analytical insight in the
understanding of molecular spectral lines, their absorption and emission
profiles and processes such as FRET migration of energy in donor-acceptor
configurations. We also made the analytical connection between the molecular
dimer model, widely used in quantum chemistry for studying vibronic quantum
coherence and the dipole-dipole interaction Hamiltonian.

## Acknowledgments

We acknowledge useful discussion with A. Dantan. We acknowledge financial
support from the Max Planck Society and from the German Federal Ministry of
Education and Research, co-funded by the European Commission (project RouTe),
project number 13N14839 within the research program “Photonik Forschung
Deutschland”. This work was funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project-ID 429529648 – TRR 306 QuCoLiMa
(“Quantum Cooperativity of Light and Matter”). M. R. acknowledges financial
support from the International Max Planck Research School - Physics of Light
(IMPRS-PL).

## References

  * Deutsch (2020) I. H. Deutsch, “Harnessing the power of the second quantum revolution,” PRX Quantum 1, 020101 (2020).
  * Facchinetti _et al._ (2016) G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, “Storing light with subradiant correlations in arrays of atoms,” Phys. Rev. Lett. 117, 243601 (2016).
  * Asenjo-Garcia _et al._ (2017) A. Asenjo-Garcia, M. Moreno-Cardoner, A. Albrecht, H. J. Kimble, and D. E. Chang, “Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays,” Phys. Rev. X 7, 031024 (2017).
  * Perczel _et al._ (2017) J. Perczel, J. Borregaard, D. E. Chang, H. Pichler, S. F. Yelin, P. Zoller, and M. D. Lukin, “Topological quantum optics in two-dimensional atomic arrays,” Phys. Rev. Lett. 119, 023603 (2017).
  * Bettles _et al._ (2017) R. J. Bettles, J. Minar, C. S. Adams, I. Lesanovsky, and B. Olmos, “Topological properties of a dense atomic lattice gas,” Phys. Rev. A 96, 041603 (2017).
  * Perczel _et al._ (2020) J. Perczel, J. Borregaard, D. E. Chang, S. F. Yelin, and M. D. Lukin, “Topological quantum optics using atomlike emitter arrays coupled to photonic crystals,” Phys. Rev. Lett. 124, 083603 (2020).
  * Glaetzle _et al._ (2017) A. W. Glaetzle, K. Ender, D. S. Wild, S. Choi, H. Pichler, M. D. Lukin, and P. Zoller, “Quantum spin lenses in atomic arrays,” Phys. Rev. X 7, 031049 (2017).
  * Guimond _et al._ (2019) P.-O. Guimond, A. Grankin, D. V. Vasilyev, B. Vermersch, and P. Zoller, “Subradiant Bell states in distant atomic arrays,” Phys. Rev. Lett. 122, 093601 (2019).
  * Chang _et al._ (2018) D. E. Chang, J. S. Douglas, A. González-Tudela, C.-L. Hung, and H. J. Kimble, “Colloquium: Quantum matter built from nanoscopic lattices of atoms and photons,” Rev. Mod. Phys. 90, 031002 (2018).
  * Ostermann _et al._ (2013) L. Ostermann, H. Ritsch, and C. Genes, “Protected state enhanced quantum metrology with interacting two-level ensembles,” Phys. Rev. Lett. 111, 123601 (2013).
  * Manzoni _et al._ (2018) M. T. Manzoni, M. Moreno-Cardoner, A. Asenjo-Garcia, J. V. Porto, A. V. Gorshkov, and D. E. Chang, “Optimization of photon storage fidelity in ordered atomic arrays,” New Journal of Physics 20, 083048 (2018).
  * Zhang and Mølmer (2020) Y.-X. Zhang and K. Mølmer, “Subradiant emission from regular atomic arrays: Universal scaling of decay rates from the generalized Bloch theorem,” Phys. Rev. Lett. 125, 253601 (2020).
  * Rui _et al._ (2020) J. Rui, D. Wei, A. Rubio-Abadal, S. Hollerith, J. Zeiher, D. M. Stamper-Kurn, C. Gross, and I. Bloch, “A subradiant optical mirror formed by a single structured atomic layer,” Nature 583, 369 (2020).
  * Shahmoon _et al._ (2019) E. Shahmoon, M. D. Lukin, and S. F. Yelin, “Chapter one - collective motion of an atom array under laser illumination,” (Academic Press, 2019) p. 1.
  * Shahmoon _et al._ (2020) E. Shahmoon, M. D. Lukin, and S. F. Yelin, “Quantum optomechanics of a two-dimensional atomic array,” Phys. Rev. A 101, 063833 (2020).
  * Plankensteiner _et al._ (2017) D. Plankensteiner, C. Sommer, H. Ritsch, and C. Genes, “Cavity antiresonance spectroscopy of dipole coupled subradiant arrays,” Phys. Rev. Lett. 119, 093601 (2017).
  * Plankensteiner _et al._ (2019) D. Plankensteiner, C. Sommer, M. Reitz, H. Ritsch, and C. Genes, “Enhanced collective Purcell effect of coupled quantum emitter systems,” Phys. Rev. A 99, 043843 (2019).
  * Bettles _et al._ (2020) R. J. Bettles, M. D. Lee, S. A. Gardiner, and J. Ruostekoski, “Quantum and nonlinear effects in light transmitted through planar atomic arrays,” Communications Physics 3, 141 (2020).
  * Haroche and Kleppner (1989) S. Haroche and D. Kleppner, “Cavity Quantum Electrodynamics,” Phys. Today 42, 24 (1989).
  * Berman (1994) P. Berman, _Cavity Quantum Electrodynamics_ , Advances in atomic, molecular, and optical physics (Academic Press, 1994).
  * Walther _et al._ (2006) H. Walther, B. T. Varcoe, B. Englert, and T. Becker, “Cavity Quantum Electrodynamics,” Rep. Prog. Phys. 69, 1325 (2006).
  * Haroche and Raimond (2013) S. Haroche and J.-M. Raimond, _Exploring the Quantum: Atoms, Cavities, and Photons_ (Oxford University Press, 2013).
  * Dicke (1954) R. H. Dicke, “Coherence in spontaneous radiation processes,” Phys. Rev. 93, 99 (1954).
  * Gross and Haroche (1982) M. Gross and S. Haroche, “Superradiance: An essay on the theory of collective spontaneous emission,” Physics Reports 93, 301 (1982).
  * Bohnet _et al._ (2012) J. Bohnet, Z. Chen, J. Weiner, D. Meiser, M. Holland, and J. K. Thompson, “A steady-state superradiant laser with less than one intracavity photon,” Nature 484, 78 (2012).
  * Norcia _et al._ (2016) M. A. Norcia, M. N. Winchester, J. R. K. Cline, and J. K. Thompson, “Superradiance on the millihertz linewidth strontium clock transition,” Science Advances 2 (2016).
  * Wild _et al._ (2018) D. S. Wild, E. Shahmoon, S. F. Yelin, and M. D. Lukin, “Quantum nonlinear optics in atomically thin materials,” Phys. Rev. Lett. 121, 123606 (2018).
  * Cernotik _et al._ (2019) O. Cernotik, A. Dantan, and C. Genes, “Cavity quantum electrodynamics with frequency-dependent reflectors,” Phys. Rev. Lett. 122, 243601 (2019).
  * Orgiu _et al._ (2015) E. Orgiu, J. George, J. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samori, and T. W. Ebbesen, “Conductivity in organic semiconductors hybridized with the vacuum field,” Nat. Mater. 14, 1123 (2015).
  * Schachenmayer _et al._ (2015) J. Schachenmayer, C. Genes, E. Tignone, and G. Pupillo, “Cavity-enhanced transport of excitons,” Phys. Rev. Lett. 114, 196403 (2015).
  * Feist and Garcia-Vidal (2015) J. Feist and F. J. Garcia-Vidal, “Extraordinary exciton conductance induced by strong coupling,” Phys. Rev. Lett. 114, 196402 (2015).
  * Hagenmüller _et al._ (2017) D. Hagenmüller, J. Schachenmayer, S. Schütz, C. Genes, and G. Pupillo, “Cavity-enhanced transport of charge,” Phys. Rev. Lett. 119, 223601 (2017).
  * Hagenmüller _et al._ (2018) D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genes, and G. Pupillo, “Cavity-assisted mesoscopic transport of fermions: Coherent and dissipative dynamics,” Phys. Rev. B 97, 205303 (2018).
  * Zhong _et al._ (2016) X. Zhong, T. Chervy, S. Wang, J. George, A. Thomas, J. Hutchinson, E. Devaux, C. Genet, and T. W. Ebbesen, “Non-radiative energy transfer mediated by hybrid light-matter states,” Angew. Chem. 55, 6202 (2016).
  * Zhong _et al._ (2017) X. Zhong, T. Chervy, L. Zhang, A. Thomas, J. George, C. Genet, J. Hutchinson, and T. W. Ebbesen, “Energy transfer between spatially separated entangled molecules,” Angew. Chem. 56, 9034 (2017).
  * Hutchison _et al._ (2012) J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, “Modifying chemical landscapes by coupling to vacuum fields,” Angew. Chem. Int. Ed. 51, 1592 (2012).
  * Galego _et al._ (2017a) J. Galego, F. Garcia-Vidal, and J. Feist, “Suppressing photochemical reactions with quantized light fields,” Nat. Commun. 7, 13841 (2017a).
  * Herrera and Spano (2016) F. Herrera and F. C. Spano, “Cavity-controlled chemistry in molecular ensembles,” Phys. Rev. Lett. 116, 238301 (2016).
  * Galego _et al._ (2017b) J. Galego, F. J. Garcia-Vidal, and J. Feist, “Many-molecule reaction triggered by a single photon in polaritonic chemistry,” Phys. Rev. Lett. 119, 136001 (2017b).
  * Martinez-Martinez _et al._ (2018) L. A. Martinez-Martinez, R. F. Ribeiro, J. Campos-González-Angulo, and J. Yuen-Zhou, “Can ultrastrong coupling change ground-state chemical reactions?” ACS Photonics 5, 167 (2018).
  * Lettow _et al._ (2010) R. Lettow, Y. L. A. Rezus, A. Renn, G. Zumofen, E. Ikonen, S. Götzinger, and V. Sandoghdar, “Quantum interference of tunably indistinguishable photons from remote organic molecules,” Phys. Rev. Lett. 104, 123605 (2010).
  * Wang _et al._ (2019a) D. Wang, H. Kelkar, D. Martin-Cano, D. Rattenbacher, A. Shkarin, T. Utikal, S. Götzinger, and V. Sandoghdar, “Turning a molecule into a coherent two-level quantum system,” Nature Physics 15, 483 (2019a).
  * Förster (1948) T. Förster, “Zwischenmolekulare Energiewanderung und Fluoreszenz,” Ann. Phys. 437, 55 (1948).
  * Blankenship (2014) R. E. Blankenship, _Molecular Mechanisms of Photosynthesis, 2nd Edition_ (John Wiley & Sons, Inc., Hoboken, NJ, 2014).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione, _The theory of open quantum systems_ (Oxford University Press, 2002).
  * Weiss (1999) U. Weiss, _Quantum Dissipative Systems_, Series in modern condensed matter physics (World Scientific, 1999).
  * Scully and Zubairy (1997) M. Scully and M. Zubairy, _Quantum Optics_ (Cambridge University Press, 1997).
  * Gardiner and Zoller (2004) C. Gardiner and P. Zoller, _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ , Vol. 56 (Springer Science & Business Media, 2004).
  * Vogel and Welsch (2006) W. Vogel and D. Welsch, _Quantum Optics_ (Wiley, 2006).
  * Walls and Milburn (2012) D. Walls and G. Milburn, _Quantum Optics_, Springer Study Edition (Springer Berlin Heidelberg, 2012).
  * Holstein (1959) T. Holstein, “Studies of polaron motion: Part I. The molecular-crystal model,” Annals of Physics 8, 325 (1959).
  * Reitz _et al._ (2019) M. Reitz, C. Sommer, and C. Genes, “Langevin approach to quantum optics with molecules,” Phys. Rev. Lett. 122, 203602 (2019).
  * Reitz _et al._ (2020) M. Reitz, C. Sommer, B. Gurlek, V. Sandoghdar, D. Martin-Cano, and C. Genes, “Molecule-photon interactions in phononic environments,” Phys. Rev. Research 2, 033270 (2020).
  * Diehl _et al._ (2014) F. P. Diehl, C. Roos, A. Duymaz, B. Lunkenheimer, A. Köhn, and T. Basche, “Emergence of coherence through variation of intermolecular distances in a series of molecular dimers,” J. Phys. Chem. Lett 5, 262 (2014).
  * Loudon (1973) R. Loudon, _The Quantum Theory of Light_, Oxford (Clarendon Press, 1973).
  * Lehmberg (1970) R. Lehmberg, “Radiation from an N-atom system. I. General formalism,” Phys. Rev. A 2, 883 (1970).
  * Flick _et al._ (2017a) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry,”  114, 3026 (2017a).
  * Hettich _et al._ (2002) C. Hettich, C. Schmitt, J. Zitzmann, S. Kühn, I. Gerhardt, and V. Sandoghdar, “Nanometer resolution and coherent optical dipole coupling of two individual molecules,” Science 298, 385 (2002).
  * Plankensteiner _et al._ (2015) D. Plankensteiner, L. Ostermann, H. Ritsch, and C. Genes, “Selective protected state preparation of coupled dissipative quantum emitters,” Sci. Rep. 5, 16231 (2015).
  * Hebenstreit _et al._ (2017) M. Hebenstreit, B. Kraus, L. Ostermann, and H. Ritsch, “Subradiance via entanglement in atoms with several independent decay channels,” Phys. Rev. Lett. 118, 143602 (2017).
  * Plenio (2005) M. B. Plenio, “Logarithmic negativity: a full entanglement monotone that is not convex,” Phys. Rev. Lett. 95, 090503 (2005).
  * Alaee _et al._ (2020) R. Alaee, B. Gurlek, M. Albooyeh, D. Martín-Cano, and V. Sandoghdar, “Quantum metamaterials with magnetic response at optical frequencies,” Phys. Rev. Lett. 125, 063601 (2020).
  * Ballantine and Ruostekoski (2020a) K. E. Ballantine and J. Ruostekoski, “Optical magnetism and Huygens’ surfaces in arrays of atoms induced by cooperative responses,” Phys. Rev. Lett. 125, 143604 (2020a).
  * Bekenstein _et al._ (2020) R. Bekenstein, I. Pikovski, H. Pichler, E. Shahmoon, S. F. Yelin, and M. D. Lukin, “Quantum metasurfaces with atom arrays,” Nature Physics 16, 676–681 (2020).
  * Grankin _et al._ (2018) A. Grankin, P. O. Guimond, D. V. Vasilyev, B. Vermersch, and P. Zoller, “Free-space photonic quantum link and chiral quantum optics,” Phys. Rev. A 98, 043825 (2018).
  * Ostermann _et al._ (2014) L. Ostermann, D. Plankensteiner, H. Ritsch, and C. Genes, “Protected subspace Ramsey spectroscopy,” Phys. Rev. A 90, 053823 (2014).
  * Holzinger _et al._ (2020) R. Holzinger, D. Plankensteiner, L. Ostermann, and H. Ritsch, “Nanoscale coherent light source,” Phys. Rev. Lett. 124, 253603 (2020).
  * Bettles _et al._ (2016) R. J. Bettles, S. A. Gardiner, and C. S. Adams, “Enhanced optical cross section via collective coupling of atomic dipoles in a 2D array,” Phys. Rev. Lett. 116, 103602 (2016).
  * Shahmoon _et al._ (2017) E. Shahmoon, D. S. Wild, M. D. Lukin, and S. F. Yelin, “Cooperative resonances in light scattering from two-dimensional atomic arrays,” Phys. Rev. Lett. 118, 113601 (2017).
  * Su _et al._ (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, “Solitons in polyacetylene,” Phys. Rev. Lett. 42, 1698 (1979).
  * Asbóth _et al._ (2016) J. K. Asbóth, L. Oroszlány, and A. Pályi, _A Short Course on Topological Insulators_ (Springer, New York, 2016).
  * Ramsey (1990) N. Ramsey, _Molecular beams_ (Oxford University Press, New York, 1990).
  * Cogdell _et al._ (2006) R. J. Cogdell, A. Gall, and J. Köhler, “The architecture and function of the light-harvesting apparatus of purple bacteria: from single molecules to in vivo membranes,” Quarterly Reviews of Biophysics 39, 227 (2006).
  * Bourne Worster _et al._ (2019) S. Bourne Worster, C. Stross, F. M. W. C. Vaughan, N. Linden, and F. R. Manby, “Structure and efficiency in bacterial photosynthetic light harvesting,” The Journal of Physical Chemistry Letters 10, 7383 (2019).
  * Caycedo-Soler _et al._ (2017) F. Caycedo-Soler, C. A. Schroeder, C. Autenrieth, A. Pick, R. Ghosh, S. F. Huelga, and M. B. Plenio, “Quantum redirection of antenna absorption to photosynthetic reaction centers,” The Journal of Physical Chemistry Letters 8, 6015 (2017).
  * Glicenstein _et al._ (2020) A. Glicenstein, G. Ferioli, N. Sibalic, L. Brossard, I. Ferrier-Barbut, and A. Browaeys, “Collective shift in resonant light scattering by a one-dimensional atomic chain,” Phys. Rev. Lett. 124, 253602 (2020).
  * Moreno-Cardoner _et al._ (2019) M. Moreno-Cardoner, D. Plankensteiner, L. Ostermann, D. E. Chang, and H. Ritsch, “Subradiance-enhanced excitation transfer between dipole-coupled nanorings of quantum emitters,” Phys. Rev. A 100, 023806 (2019).
  * Needham _et al._ (2019) J. A. Needham, I. Lesanovsky, and B. Olmos, “Subradiance-protected excitation transport,” New Journal of Physics 21, 073061 (2019).
  * Ballantine and Ruostekoski (2020b) K. E. Ballantine and J. Ruostekoski, “Subradiance-protected excitation spreading in the generation of collimated photon emission from an atomic array,” Phys. Rev. Research 2, 023086 (2020b).
  * Tavis and Cummings (1968) M. Tavis and F. W. Cummings, “Exact solution for an N-molecule—radiation-field Hamiltonian,” Phys. Rev. 170, 379 (1968).
  * Jaynes and Cummings (1963) E. T. Jaynes and F. W. Cummings, “Comparison of quantum and semiclassical radiation theories with application to the beam maser,” Proceedings of the IEEE 51, 89 (1963).
  * Braak (2011) D. Braak, “Integrability of the Rabi Model,” Phys. Rev. Lett. 107, 100401 (2011).
  * Rabi (1936) I. I. Rabi, “On the process of space quantization,” Phys. Rev. 49, 324 (1936).
  * Emary and Brandes (2003) C. Emary and T. Brandes, “Quantum chaos triggered by precursors of a quantum phase transition: The Dicke model,” Phys. Rev. Lett. 90, 044101 (2003).
  * Braak (2013) D. Braak, “Solution of the Dicke model for N = 3,” Journal of Physics B: Atomic, Molecular and Optical Physics 46, 224007 (2013).
  * Hepp and Lieb (1973) K. Hepp and E. H. Lieb, “On the superradiant phase transition for molecules in a quantized radiation field: the Dicke maser model,” Annals of Physics 76, 360 (1973).
  * Brankov _et al._ (1975) I. Brankov, V. Zagrebnov, and I. Tonchev, “Asymptotically exact solution of the generalized Dicke model,” Theor Math Phys 22, 13 (1975).
  * Bogolubov and Plechko (1976) N. Bogolubov and V. Plechko, “A class of exactly soluble many-body Hamiltonians with the interaction of substance and boson field,” Physica A: Statistical Mechanics and its Applications 82, 163 (1976).
  * Wang and Hioe (1973) Y. K. Wang and F. T. Hioe, “Phase transition in the Dicke model of superradiance,” Phys. Rev. A 7, 831 (1973).
  * Jacobs (2010) K. Jacobs, _Stochastic Processes for Physicists: Understanding Noisy Systems_ (Cambridge University Press, New York, 2010).
  * Sommer _et al._ (2021) C. Sommer, M. Reitz, F. Mineo, and C. Genes, “Molecular polaritonics in dense mesoscopic disordered ensembles,” Phys. Rev. Research 3, 033141 (2021).
  * Miroshnichenko _et al._ (2010) A. E. Miroshnichenko, S. Flach, and Y. S. Kivshar, “Fano resonances in nanoscale structures,” Rev. Mod. Phys. 82, 2257 (2010).
  * Chang-Hasnain and Yang (2012) C. J. Chang-Hasnain and W. Yang, “High-contrast gratings for integrated optoelectronics,” Adv. Opt. Photon. 4, 379 (2012).
  * Zhou _et al._ (2014) W. Zhou, D. Zhao, Y.-C. Shuai, H. Yang, S. Chuwongin, A. Chadha, J.-H. Seo, K. X. Wang, V. Liu, Z. Ma, and S. Fan, “Progress in 2D photonic crystal Fano resonance photonics,” Progress in Quantum Electronics 38, 1–74 (2014).
  * Limonov _et al._ (2017) M. F. Limonov, M. V. Rybin, A. N. Poddubny, and Y. S. Kivshar, “Fano resonances in photonics,” Nature Photonics 11, 543 (2017).
  * Zeytinoğlu _et al._ (2017) S. Zeytinoğlu, C. Roth, S. Huber, and A. Imamoğlu, “Atomically thin semiconductors as nonlinear mirrors,” Physical Review A 96, 031801 (2017).
  * Back _et al._ (2018) P. Back, S. Zeytinoglu, A. Ijaz, M. Kroner, and A. Imamoğlu, “Realization of an electrically tunable narrow-bandwidth atomically thin mirror using monolayer MoSe2,” Phys. Rev. Lett. 120, 037401 (2018).
  * Scuri _et al._ (2018) G. Scuri, Y. Zhou, A. A. High, D. S. Wild, C. Shu, K. De Greve, L. A. Jauregui, T. Taniguchi, K. Watanabe, P. Kim, M. D. Lukin, and H. Park, “Large excitonic reflectivity of monolayer MoSe2 encapsulated in hexagonal boron nitride,” Phys. Rev. Lett. 120, 037402 (2018).
  * Mork _et al._ (2014) J. Mork, Y. Chen, and M. Heuck, “Photonic Crystal Fano Laser: Terahertz modulation and ultrashort pulse generation,” Phys. Rev. Lett. 113, 163901 (2014).
  * Bohnet _et al._ (2014) J. G. Bohnet, Z. Chen, J. M. Weiner, K. C. Cox, and J. K. Thompson, “Linear-response theory for superradiant lasers,” Phys. Rev. A 89, 013806 (2014).
  * Dantan _et al._ (2014) A. Dantan, B. Nair, G. Pupillo, and C. Genes, “Hybrid cavity mechanics with doped systems,” Phys. Rev. A 90, 033820 (2014).
  * Sentef _et al._ (2018) M. A. Sentef, M. Ruggenthaler, and A. Rubio, “Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity,” Sci. Adv. 4, 6969 (2018).
  * Thomas _et al._ (2019) A. Thomas, E. Devaux, K. Nagarajan, T. Chervy, M. Seidel, D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genet, G. Pupillo, and T. W. Ebbesen, “Exploring superconductivity under strong coupling with the vacuum electromagnetic field,” arXiv:1911.01459 (2019).
  * Kampschulte and Denschlag (2018) T. Kampschulte and J. H. Denschlag, “Cavity-controlled formation of ultracold molecules,” New Journal of Physics 20, 123015 (2018).
  * Kiffner _et al._ (2019) M. Kiffner, J. R. Coulthard, F. Schlawin, A. Ardavan, and D. Jaksch, “Manipulating quantum materials with quantum light,” Phys. Rev. B 99, 085116 (2019).
  * Wang _et al._ (2019b) X. Wang, E. Ronca, and M. A. Sentef, “Cavity quantum electrodynamical Chern insulator: Towards light-induced quantized anomalous Hall effect in graphene,” Phys. Rev. B 99, 235156 (2019b).
  * Rohn _et al._ (2020) J. Rohn, M. Hörmann, C. Genes, and K. P. Schmidt, “Ising model in a light-induced quantized transverse field,” Phys. Rev. Research 2, 023131 (2020).
  * Sentef _et al._ (2020) M. A. Sentef, J. Li, F. Künzel, and M. Eckstein, “Quantum to classical crossover of Floquet engineering in correlated quantum systems,” Phys. Rev. Research 2, 033033 (2020).
  * Ashida _et al._ (2020) Y. Ashida, A. Imamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, “Quantum electrodynamic control of matter: Cavity-enhanced ferroelectric phase transition,” Phys. Rev. X 10, 041027 (2020).
  * Toninelli _et al._ (2021) 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. L. Koppens, B. Kozankiewicz, A. Gourdon, V. Sandoghdar, and M. Orrit, “Single organic molecules for photonic quantum technologies,” Nature Materials (2021).
  * Basché _et al._ (1992) T. Basché, W. E. Moerner, M. Orrit, and H. Talon, “Photon antibunching in the fluorescence of a single dye molecule trapped in a solid,” Phys. Rev. Lett. 69, 1516 (1992).
  * De Martini _et al._ (1996) F. De Martini, G. Di Giuseppe, and M. Marrocco, “Single-mode generation of quantum photon states by excited single molecules in a microcavity trap,” Phys. Rev. Lett. 76, 900 (1996).
  * Lounis and Moerner (2000) B. Lounis and W. E. Moerner, “Single photons on demand from a single molecule at room temperature,” Nature 407, 491 (2000).
  * Pazzagli _et al._ (2018) S. Pazzagli, P. Lombardi, D. Martella, M. Colautti, B. Tiribilli, F. S. Cataliotti, and C. Toninelli, “Self-assembled nanocrystals of polycyclic aromatic hydrocarbons show photostable single-photon emission,” ACS Nano 12, 4295 (2018).
  * Maser _et al._ (2016) A. Maser, B. Gmeiner, T. Utikal, S. Götzinger, and V. Sandoghdar, “Few-photon coherent nonlinear optics with a single molecule,” Nature Photonics 10, 450 (2016).
  * Roelli _et al._ (2016) P. Roelli, C. Galland, N. Piro, and T. J. Kippenberg, “Molecular cavity optomechanics as a theory of plasmon-enhanced Raman scattering,” Nat. Nanotechnol. 11, 164 (2016).
  * Neuman _et al._ (2019) T. Neuman, R. Esteban, G. Giedke, M. K. Schmidt, and J. Aizpurua, “Quantum description of surface-enhanced resonant Raman scattering within a hybrid-optomechanical model,” Phys. Rev. A 100, 043422 (2019).
  * Benz _et al._ (2016) F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg, “Single-molecule optomechanics in “picocavities”,” Science 354, 726 (2016).
  * Clear _et al._ (2020) C. Clear, R. C. Schofield, K. D. Major, J. Iles-Smith, A. S. Clark, and D. P. S. McCutcheon, “Phonon-induced optical dephasing in single organic molecules,” Phys. Rev. Lett. 124, 153602 (2020).
  * Witkowski and Moffitt (1960) A. Witkowski and W. Moffitt, “Electronic spectra of dimers: Derivation of the fundamental vibronic equation,” J. Chem. Phys. 33, 872–875 (1960).
  * Riley and Krausz (1990) M. J. Riley and E. R. Krausz, “The symmetric dimer: Effects of second-order vibronic coupling,” Chemical Physics 148, 229 (1990).
  * Eisfeld _et al._ (2005) A. Eisfeld, L. Braun, W. T. Strunz, J. S. Briggs, J. Beck, and V. Engel, “Vibronic energies and spectra of molecular dimers,” J. Chem. Phys. 122, 134103 (2005).
  * Lippitz _et al._ (2004) M. Lippitz, C. G. Hübner, T. Christ, H. Eichner, P. Bordat, A. Herrmann, K. Müllen, and T. Basché, “Coherent electronic coupling versus localization in individual molecular dimers,” Phys. Rev. Lett. 92, 103001 (2004).
  * Engel _et al._ (2007) G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mančal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, “Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems,” Nature 446, 782 (2007).
  * Streltsov _et al._ (2017) A. Streltsov, G. Adesso, and M. B. Plenio, “Colloquium: Quantum coherence as a resource,” Rev. Mod. Phys. 89, 041003 (2017).
  * May and Kühn (2011) V. May and O. Kühn, _Charge and energy transfer dynamics in molecular systems_ (Wiley-VCH Verlag GmbH, New York, 2011).
  * Garcia-Vidal and Feist (2017) F. J. Garcia-Vidal and J. Feist, “Long-distance operator for energy transfer,” Science 357, 1357 (2017).
  * Du _et al._ (2018) M. Du, L. A. Martínez-Martínez, R. F. Ribeiro, Z. Hu, V. M. Menon, and J. Yuen-Zhou, “Theory for polariton-assisted remote energy transfer,” Chem. Sci. 9, 6659 (2018).
  * Reitz _et al._ (2018) M. Reitz, F. Mineo, and C. Genes, “Energy transfer and correlations in cavity-embedded donor-acceptor configurations,” Scientific Reports 8, 9050 (2018).
  * Kansanen _et al._ (2019) K. S. U. Kansanen, A. Asikainen, J. J. Toppari, G. Groenhof, and T. T. Heikkilä, “Theory for the stationary polariton response in the presence of vibrations,” Phys. Rev. B 100, 245426 (2019).
  * (131) T. Schwartz, J. A. Hutchinson, J. Léonard, C. Genet, S. Haacke, and T. W. Ebbesen, “Polariton dynamics under strong Light–Molecule coupling,” ChemPhysChem 14, 125.
  * Canaguier-Durand _et al._ (2015) A. Canaguier-Durand, C. Genet, A. Lambrecht, T. W. Ebbesen, and S. Reynaud, “Non-Markovian polariton dynamics in organic strong coupling,” The European Physical Journal D 69, 24 (2015).
  * Neuman and Aizpurua (2018) T. Neuman and J. Aizpurua, “Origin of the asymmetric light emission from molecular exciton polaritons,” Optica 5, 1247 (2018).
  * Zeb _et al._ (2018) M. A. Zeb, P. G. Kirton, and J. Keeling, “Exact states and spectra of vibrationally dressed polaritons,” ACS Photonics 5, 249 (2018).
  * Herrera and Spano (2017) F. Herrera and F. C. Spano, “Dark vibronic polaritons and the spectroscopy of organic microcavities,” Phys. Rev. Lett. 118, 223601 (2017).
  * Herrera and Spano (2018) F. Herrera and F. C. Spano, “Theory of nanoscale organic cavities: The essential role of vibration-photon dressed states,” ACS Photonics 5, 65 (2018).
  * Flick _et al._ (2017b) J. Flick, H. Appel, M. Ruggenthaler, and A. Rubio, “Cavity Born–Oppenheimer approximation for correlated electron–nuclear-photon systems,” Journal of Chemical Theory and Computation 13, 1616 (2017b).
  * del Pino _et al._ (2018) J. del Pino, F. A. Y. N. Schröder, A. W. Chin, J. Feist, and F. J. Garcia-Vidal, “Tensor network simulation of non-Markovian dynamics in organic polaritons,” Phys. Rev. Lett. 121, 227401 (2018).
  * Haake (1973) F. Haake, “Statistical treatment of open systems by generalized master equations,” in _Springer Tracts in Modern Physics: Ergebnisse der exakten Naturwissenschaftenc; Volume 66_, edited by G. Höhler (Springer Berlin Heidelberg, Berlin, Heidelberg, 1973) p. 98.
  * Novotny and Hecht (2006) L. Novotny and B. Hecht, _Principles of Nano-Optics_ (Cambridge University Press, 2006).
  * Genes _et al._ (2008) C. Genes, D. Vitali, P. Tombesi, S. Gigan, and M. Aspelmeyer, “Ground-state cooling of a micromechanical oscillator: Comparing cold damping and cavity-assisted cooling schemes,” Phys. Rev. A 77, 033804 (2008).

## Appendix

###  A1 Master equations for open-system dynamics

We provide here a short review of the master equation approach in open-system
dynamics closely following Ref. Breuer and Petruccione (2002) but also
extensively covered in, among others, Refs. Loudon (1973); Weiss (1999);
Scully and Zubairy (1997); Gardiner and Zoller (2004); Vogel and Welsch
(2006); Walls and Milburn (2012). In particular we illustrate how this
approach deals with the process of spontaneous emission of a single and many
quantum emitters, with photon loss of an optical cavity mode and how it can be
used to show the occurrence of cavity Dicke superradiant decay.

General formalism

Let us consider a system $S$ which is weakly coupled to a bath (or reservoir)
$B$. The Hamiltonian of the total system is given by

$\displaystyle\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{\text{int}}=\mathcal{H}_{S}\otimes\mathbb{1}_{B}+\mathbb{1}_{S}\otimes\mathcal{H}_{B}+\mathcal{H}_{\text{int}},$
(A1)

where $\mathcal{H}_{S}$ and $\mathcal{H}_{B}$ denote the free Hamiltonians of
the system and bath, respectively, and $\mathcal{H}_{\text{int}}$ accounts for
the interaction between the two. The symbols $\mathbb{1}_{B}$ and
$\mathbb{1}_{S}$ represent the unity in the $B$ and $S$ Hilbert spaces. In the
Schrödinger picture, the density operator satisfies
$\dot{\rho}=-i[\mathcal{H},\rho]$. We remove the free Hamiltonian by going
into an interaction picture with a unitary transformation
$\rho_{I}(t)=\mathcal{U}_{0}^{\dagger}(t)\rho(t)\mathcal{U}_{0}$ where
$\mathcal{U}_{0}(t)=e^{-i\mathcal{H}_{0}t}$. The von Neumann equation in the
interaction picture then takes the form

$\displaystyle\frac{d}{dt}\rho(t)=-i[\mathcal{H}_{\text{int}}(t),\rho(t)],$
(A2)

where the interaction Hamiltonian is now time-dependent
$\mathcal{H}_{\text{int}}(t)=\mathcal{U}^{\dagger}_{0}(t)\mathcal{H}_{\text{int}}\mathcal{U}_{0}(t)$
and we removed the index $I$ which served to indicate the interaction picture.
The formal solution of this is given by

$\displaystyle\rho(t)=\rho(0)-i\int_{0}^{t}ds\,[\mathcal{H}_{\text{int}}(s),\rho(s)].$
(A3)

Inserting this back into Eq. (A2) and taking the trace over the bath
$\rho_{S}=\mathrm{Tr}_{B}[\rho]$, one finds

$\displaystyle\frac{d}{dt}\rho_{S}(t)=-\int_{0}^{t}ds\,\mathrm{Tr}_{B}[\mathcal{H}_{\text{int}}(t),[\mathcal{H}_{\text{int}}(s),\rho(s)]],$
(A4)

with $\mathrm{Tr}_{B}[\mathcal{H}_{\text{int}},\rho(0)]=0$. The above
expression still contains the density operator of the total system on its
right-hand side. The density matrix of the total system should show deviations
on the order of $\mathcal{H}_{\text{int}}$ from an uncorrelated state
$\rho(t)=\rho_{S}(t)\otimes\rho_{B}+\mathcal{O}(\mathcal{H}_{\text{int}})$.
Assuming the coupling between system and bath to be weak (Born approximation),
one can neglect terms higher than second order in $\mathcal{H}_{\text{int}}$
and obtain an integro-differential equation for the reduced system operator

$\displaystyle\frac{d}{dt}\rho_{S}(t)=-\int_{0}^{t}ds\,\mathrm{Tr}_{B}[\mathcal{H}_{\text{int}}(t),[\mathcal{H}_{\text{int}}(s),\rho_{S}(s)\otimes\rho_{B}]].$
(A5)

A more detailed discussion of the Born approximation can also be found in Ref.
Haake (1973). This equation of motion can be brought in time-local form by
replacing $\rho_{S}(s)$ with $\rho_{S}(t)$ which is then called the Redfield
equation. Finally, to obtain a Markovian master equation, one first performs a
change of variables $s\to t-s$. Under the assumption that the bath
correlations decay quickly as compared to the time scale over which the state
of the system varies (as induced by the bath), one may replace the upper bound
of the integral by infinity and obtain

$\displaystyle\frac{d}{dt}\rho_{S}(t)=-\int_{0}^{\infty}ds\,\mathrm{Tr}_{B}[\mathcal{H}_{\text{int}}(t),[\mathcal{H}_{\text{int}}(t-s),\rho_{S}(t)\otimes\rho_{B}]].$
(A6)

Spontaneous emission of a single atom

Let us consider an interaction of the form
$\mathcal{H}_{\text{int}}=\sigma^{\dagger}\mathcal{F}(t)+\mathcal{F}^{\dagger}(t)\sigma$.
For spontaneous emission of a single atom, the term driving the emitter
expresses as
$\mathcal{F}(t)=\sum_{\mathbf{k}}g_{\mathbf{k}}a_{\mathbf{k}}e^{i(\omega_{0}-\omega_{k})t}$,
where the coupling is given by

$\displaystyle
g_{\mathbf{k}}=\sqrt{\frac{\omega_{k}}{2\epsilon_{0}\mathcal{V}}}\bm{\epsilon}_{\mathbf{k}}^{(\lambda)}\cdot\mathbf{d}_{\text{eg}}.$
(A7)

Note that in our notation, the sum over $k$ vectors implicitly contains a sum
over the two orthogonal polarizations $\bm{\epsilon}_{\mathbf{k}}^{(1)}$ and
$\bm{\epsilon}_{\mathbf{k}}^{(2)}$. Inserting the interaction Hamiltonian into
Eq. (A6), using the cyclic property of the trace (e.g.,
$\mathrm{Tr}_{B}[\mathcal{F}(t)\rho_{B}(0)\mathcal{F}^{\dagger}(t-s)]=\mathrm{Tr}_{B}[\mathcal{F}^{\dagger}(t-s)\mathcal{F}(t)\rho_{B}(0)]=\braket{\mathcal{F}^{\dagger}(t-s)\mathcal{F}(t)}$)
and assuming a bath $\rho_{B}=\ket{\mathrm{vac}}\bra{\mathrm{vac}}$, the only
remaining terms are

$\displaystyle\dot{\rho}(t)=\int_{0}^{\infty}ds\left\\{\left[\sigma\rho(t)\sigma^{\dagger}-\sigma^{\dagger}\sigma\rho(t)\right]\braket{\mathcal{F}(t)\mathcal{F}^{\dagger}(t-s)}+\mathrm{h.c.}\right\\}.$
(A8)

In consistence with the main text of the tutorial, we now denote the system
density matrix simply by $\rho$. One is left with evaluating the term
$\braket{\mathcal{F}(t)\mathcal{F}^{\dagger}(t-s)}=\sum_{\mathbf{k}}|g_{\mathbf{k}}|^{2}e^{i(\omega_{0}-\omega_{k})s}$
since
$\braket{a_{\mathbf{k}}a_{\mathbf{k}^{\prime}}^{\dagger}}=\delta_{\mathbf{k}{\mathbf{k}^{\prime}}}$.
Identifying the minimal volume in $k$-space which is occupied by two modes of
orthogonal polarization and same $k$ vector as $(2\pi/\ell)^{3}$, we can
replace the sum by an integral and go into spherical coordinates

$\displaystyle\frac{1}{\mathcal{V}}\sum_{\mathbf{k}}\to\int\frac{d^{3}k}{(2\pi)^{3}}=\frac{1}{(2\pi
c)^{3}}\int_{0}^{\infty}d\omega_{k}\omega_{k}^{2}\int d\Omega_{k}.$ (A9)

Additionally, one still has to account for the sum over the two orthogonal
polarizations which can be performed as
$\sum_{\lambda}(\bm{\epsilon}_{\mathbf{k}}^{(\lambda)}\cdot\mathbf{d}_{\text{eg}})^{2}=d_{\text{eg}}^{2}-{(\mathbf{k}\cdot\mathbf{d}_{\text{eg}})^{2}}/k^{2}$.
For the integration over $s$, one can make use of the Sokhotski-Plemelj
theorem

$\displaystyle\int d\omega_{k}\int_{0}^{\infty}dse^{\pm
i(\omega_{k}-\omega_{0})s}=\int
d\omega_{k}\left[\pi\delta(\omega_{k}-\omega_{0})\pm
i\mathcal{P}\left(\frac{1}{\omega_{k}-\omega_{0}}\right)\right],$ (A10)

where $\mathcal{P}$ denotes the Cauchy principal value. While the real part of
the above expression gives rise to decay, the imaginary part corresponds to
the Lamb shift which will eventually give rise to a small frequency
renormalization of the atom. Neglecting this small correction and only
considering the decay, one can finally arrive at a master equation in Lindblad
form

$\displaystyle\dot{\rho}(t)=\gamma\left(2\sigma\rho\sigma^{\dagger}-\\{\sigma^{\dagger}\sigma,\rho\\}\right),$
(A11)

where $\gamma=(\omega_{0}^{3}d_{\text{eg}}^{2})/(6\pi c^{3}\epsilon_{0})$ is
the spontaneous emission rate.

The master equation for collective dynamics of closely spaced quantum-emitter
ensembles

Let us now consider the case of $\mathcal{N}$ identical emitters located at
positions $\mathbf{R}_{j}$ mutually coupled to the electromagnetic vacuum
described by the interaction
$\mathcal{H}_{\text{int}}=\sum_{j}(\sigma_{j}^{\dagger}\mathcal{F}_{j}(t)+\mathcal{F}_{j}^{\dagger}(t)\sigma_{j})$
with
$\mathcal{F}_{j}(t)=\sum_{\mathbf{k}}g_{\mathbf{k}}a_{\mathbf{k}}e^{i\mathbf{k}\mathbf{R}_{j}}e^{-i(\omega_{k}-\omega_{0})t}$,
where we assume equal orientation of all transition dipoles
$\mathbf{d}_{\text{eg}}^{j}=\mathbf{d}_{\text{eg}}$. In this case, the master
equation takes the form

$\displaystyle\dot{\rho}(t)=\sum_{j,j^{\prime}=1}^{\mathcal{N}}\int_{0}^{\infty}ds\left\\{\sigma_{j}\rho(t)\sigma_{j^{\prime}}^{\dagger}\left(\mathcal{C}_{j^{\prime}j}(s)+\mathcal{C}_{jj^{\prime}}^{*}(s)\right)-\sigma_{j}^{\dagger}\sigma_{j^{\prime}}\rho(t)\mathcal{C}_{jj^{\prime}}(s)-\rho(t)\sigma_{j}^{\dagger}\sigma_{j^{\prime}}\mathcal{C}_{j^{\prime}j}^{*}(s)\right\\},$
(A12)

where we denote the correlations as

$\displaystyle\mathcal{C}_{jj^{\prime}}(s)=\braket{\mathcal{F}_{j}(t)\mathcal{F}_{j^{\prime}}^{\dagger}(t-s)}=\sum_{\mathbf{k}}|g_{\mathbf{k}}|^{2}e^{i\mathbf{k}(\mathbf{R}_{j}-\mathbf{R}_{j^{\prime}})}e^{i(\omega_{0}-\omega_{k})s},$
(A13)

and used furthermore that
$\braket{\mathcal{F}_{j}(t)\mathcal{F}_{j^{\prime}}^{\dagger}(t-s)}=\braket{\mathcal{F}_{j^{\prime}}(t-s)\mathcal{F}_{j}^{\dagger}(t)}^{*}$.
After summing over the two orthogonal polarizations we now again perform the
continuum limit as in the previous paragraph and obtain for the correlation
function Eq. (A13)

$\displaystyle\braket{\mathcal{F}_{j}(t)\mathcal{F}_{j^{\prime}}^{\dagger}(t-s)}=\frac{d_{\text{eg}}^{2}}{4\pi^{2}c^{3}\epsilon_{0}}\int
d\omega_{k}\omega_{k}^{3}F(kR_{jj^{\prime}})e^{i(\omega_{0}-\omega_{k})s},$
(A14)

where the function $F(kR)$ is defined as the solid angle integral
$F(kR)=1/(4\pi)\int
d\Omega_{k}(1-(\mathbf{d}_{\text{eg}}\cdot\mathbf{k})^{2})/(d_{\text{eg}}^{2}k^{2})\cdot
e^{i\mathbf{k}\mathbf{R}}$. This integral can be calculated as

$\displaystyle
F(kR)=\frac{1}{4\pi}\left(1+\frac{\left(\mathbf{e}_{\mathbf{d}}\cdot\nabla_{\mathbf{R}}\right)^{2}}{k^{2}}\right)\int_{0}^{\pi}d\theta_{k}e^{ikR\cos\theta_{k}}\int_{0}^{2\pi}d\phi_{k}=\left[1+\frac{(\mathbf{e}_{\mathbf{d}}\cdot\nabla_{\mathbf{R}})^{2}}{k^{2}}\right]\frac{\sin(kR)}{kR}.$
(A15)

Finally, resolving the integral over $s$ by means of Eq. (A10) results in

$\displaystyle\int_{0}^{\infty}ds\braket{\mathcal{F}_{j}(t)\mathcal{F}_{j^{\prime}}^{\dagger}(t-s)}=\frac{3\gamma}{2}F(k_{0}R_{jj^{\prime}})-i\frac{3\gamma}{2}G(k_{0}R_{jj^{\prime}})=:\Gamma_{jj^{\prime}}+i\Omega_{jj^{\prime}},$
(A16)

where

$\displaystyle G(k_{0}R)=\frac{1}{\pi\omega_{0}^{3}}\mathcal{P}\int
d\omega_{k}\frac{\omega_{k}^{3}}{\omega_{k}-\omega_{0}}F(kR).$ (A17)

All together, this yields the master equation Eq. (6) of the main text.

The master equation for cavity-induced Dicke superradiance

For $\mathcal{N}$ identical emitters placed within the volume of an optical
resonator, the interaction Hamiltonian is given by
$\mathcal{H}_{\text{int}}=S^{\dagger}\mathcal{F}(t)+\mathcal{F}^{\dagger}(t)S$
with $\mathcal{F}(t)=ga(t)$. Using that the correlation of the field operator
is given by (assuming the cavity field to be in the vacuum field due to its
fast decay)

$\displaystyle\braket{\mathcal{F}(t)\mathcal{F}^{\dagger}(t-s)}=g^{2}e^{-\kappa
s},$ (A18)

the term appearing in the master equation is

$\displaystyle\dot{\rho}(t)=\frac{g^{2}}{\kappa}\left(2S\rho
S^{\dagger}-\\{S^{\dagger}S,\rho\\}\right).$ (A19)

Damped cavity mode

The tunneling of photons between a confined cavity mode $a$ to the outside
continuum (modes on the left $b(\omega)$ and modes on the right $c(\omega)$)
can be phenomenologically described via the interaction Hamiltonian
$\mathcal{H}_{\text{int}}=a^{\dagger}\mathcal{F}(t)+\mathcal{F}^{\dagger}(t)a$
where $\mathcal{F}(t)=\int
d\omega[V(\omega)b(\omega)+W(\omega)c(\omega)]e^{-i(\omega-\omega_{c})t}$ with
frequency-dependent rates $V(\omega)$ and $W(\omega)$. Evaluating the
correlation (approximating the outside field as
$\rho_{B}=\rho_{B}^{(L)}\otimes\rho_{B}^{(R)}$ with
$\rho_{B}^{(L,R)}=\ket{\mathrm{vac}}_{L,R}\bra{\mathrm{vac}}_{L,R}$)

$\displaystyle\braket{\mathcal{F}(t)\mathcal{F}^{\dagger}(t-s)}=\int
d\omega[V(\omega)^{2}+W(\omega)^{2}]e^{-i(\omega-\omega_{c})s},$ (A20)

and subsequently the integral over $s$, gives rise to loss rates via left and
right mirrors $\kappa_{L}=\pi V(\omega_{c})^{2}$, $\kappa_{R}=\pi
W(\omega_{c})^{2}$ and a master equation in Lindblad form

$\displaystyle\dot{\rho}(t)=\kappa\left(2a\rho
a^{\dagger}-\\{a^{\dagger}a,\rho\\}\right),$ (A21)

where $\kappa=\kappa_{L}+\kappa_{R}$ sums the loss rates via the left and
right mirror.

###  A2 Functional dependence of $F(kR)$, $G(kR)$

The expression for $F(kR)$ in Eq. (A15) can be brought into a more convenient
form widely used in the literature by expressing the nabla operator in
spherical coordinates

$\displaystyle\nabla_{\mathbf{R}}=\partial_{R}\mathbf{e}_{R}+\frac{1}{R}\partial_{\theta}\mathbf{e}_{\theta}+\frac{1}{R\sin\theta}\partial_{\phi}\mathbf{e}_{\phi}.$
(A22)

Assuming the dipoles to be aligned along the $z$ direction
$\mathbf{e}_{\mathbf{d}}=\mathbf{e}_{z}$ with
$\mathbf{e}_{z}=\cos\theta\mathbf{e}_{R}-\sin\theta\mathbf{e}_{\theta}$, we
can express

$\displaystyle\mathbf{e}_{\mathbf{d}}\cdot\nabla_{\mathbf{R}}=\cos\theta\,\partial_{R}-\frac{\sin\theta}{R}\partial_{\theta},$
(A23)

Calculating the derivatives, one obtains

$\displaystyle
F(kR)=(1-\cos^{2}\theta)\frac{\sin(kR)}{kR}+(1-3\cos^{2}\theta)\left[\frac{\cos(kR)}{(kR)^{2}}-\frac{\sin(kR)}{(kR)^{3}}\right].$
(A24)

The calculation of $G(kR)$ involves evaluation of the principal value integral
Eq. (A17). We just specify the result here for brevity

$\displaystyle
G(kR)=(1-\cos^{2}\theta)\frac{\cos(kR)}{kR}-(1-3\cos^{2}\theta)\left[\frac{\sin(kR)}{(kR)^{2}}+\frac{\cos(kR)}{(kR)^{3}}\right].$
(A25)

###  A3 Lattice sum of spherical waves

The sum over spherical waves of a 2D array with atoms periodically arranged at
positions $\mathbf{R}_{j}$ in the $x$-$y$ plane can be expressed along the $z$
direction by going into the continuum as

$\displaystyle\mathcal{S}=\sum_{j}\frac{e^{ik|z\mathbf{e}_{z}-\mathbf{R}_{j}|}}{|z\mathbf{e}_{z}-\mathbf{R}_{j}|}=\sum_{j}\int
d\mathbf{R}_{\parallel}\frac{e^{ik\sqrt{\mathbf{R}_{\parallel}^{2}+z^{2}}}}{\sqrt{\mathbf{R}_{\parallel}^{2}+z^{2}}}\delta(\mathbf{R}_{\parallel}-\mathbf{R}_{j}),$
(A26)

where we denote the in-plane components by $\mathbf{R}_{\parallel}=(x,y)$. The
spherical wave can be expanded by means of a Weyl expansion as Novotny and
Hecht (2006)

$\displaystyle\mathcal{S}=\frac{i}{2\pi}\sum_{j}\int
d\mathbf{R}_{\parallel}\int
d\mathbf{k}_{\parallel}e^{-i\mathbf{k}_{\parallel}\cdot\mathbf{R}_{\parallel}}\frac{e^{ik_{z}|z|}}{k_{z}}\delta(\mathbf{R}_{\parallel}-\mathbf{R}_{j}),$
(A27)

with in-plane momentum $\mathbf{k}_{\parallel}=(k_{x},k_{y})$ and
$k_{z}=\sqrt{k^{2}-|\mathbf{k}_{\parallel}|^{2}}$. Poisson’s summation formula
allows one to turn the sum over the real-space lattice into a sum over the
reciprocal lattice

$\displaystyle\sum_{j}\int
d\mathbf{R}_{\parallel}e^{-i\mathbf{k}_{\parallel}\cdot\mathbf{R}_{\parallel}}\delta(\mathbf{R}_{\parallel}-\mathbf{R}_{j})=\frac{(2\pi)^{2}}{a^{2}}\sum_{m,n}\delta(\mathbf{k}_{\parallel}-\mathbf{q}_{m,n}),$
(A28)

with reciprocal lattice vectors given by
$\mathbf{q}_{m,n}=(q_{m},q_{n})=\frac{2\pi}{a}(m,n)$ for a square lattice with
lattice constant $a$. This yields for the sum

$\displaystyle\mathcal{S}=\frac{2\pi
i}{a^{2}}\sum_{m,n}\frac{e^{i\sqrt{k^{2}-|\mathbf{q}_{m,n}|^{2}}|z|}}{\sqrt{k^{2}-|\mathbf{q}_{m,n}|^{2}}}.$
(A29)

###  A4 One-dimensional transfer-matrix theory

Figure 18: In- and outgoing waves for an arbitrary scatterer with reflectivity
and transmission coefficients $(r,t)$.

For an arbitrary scatterer (mirror, atom, etc.) at a fixed position, we assume
the modes on the left as $Ae^{-ikx}$ (left propagating) and $Be^{ikx}$ (right
propagating) and on the right as $Ce^{-ikx}$ (left propagating) and $De^{ikx}$
(right propagating), as illustrated in Fig. 18. We assume the scatterer to
have a complex reflectivity $r$ and transmittivity $t$. The two are connected
as $t=1+r$ and in the absence of absorption one has $|t|^{2}+|r|^{2}=1$. The
outgoing fields can be related to the incoming fields as

$\displaystyle A$ $\displaystyle=rB+tC,$ (A30a) $\displaystyle D$
$\displaystyle=tB+rC.$ (A30b)

Note that in the case where one has no incoming wave from the right ($C=0$),
this simply gives $D=tB$ and $A=rB$. The fields on the left can be connected
to the fields on the right via

$\displaystyle\begin{pmatrix}A\\\
B\end{pmatrix}=\frac{1}{t}\begin{pmatrix}t^{2}-r^{2}&r\\\
-r&1\end{pmatrix}\begin{pmatrix}C\\\ D\end{pmatrix}.$ (A31)

For the two-mirror arrangement considered in Sec. B, each mirror can be
characterized by its polarizability $\zeta_{j}$, $j=R,L$, where the
polarizability generally depends on frequency $\zeta_{j}=\zeta_{j}(\omega)$.
The transmission and reflection coefficients for each mirror are given by
$t_{j}=1/(1-i\zeta_{j})$ and $r_{j}=i\zeta_{j}/(1-i\zeta_{j})$, respectively
and the transfer matrix of each mirror reads

$\displaystyle\mathbf{T}_{j}$
$\displaystyle=\begin{pmatrix}1+i\zeta_{j}&i\zeta_{j}\\\
-i\zeta_{j}&1-i\zeta_{j}\end{pmatrix},$ (A32)

while the matrix for free-space propagation is given by
$\mathbf{T}_{f}=\mathrm{diag}(e^{i\theta},e^{-i\theta})$ where
$\theta=\omega\ell/c$. The transfer matrix of the whole system is then
obtained by multiplying the three individual matrices

$\displaystyle\mathbf{T}$
$\displaystyle=\mathbf{T}_{R}\mathbf{T}_{f}\mathbf{T}_{L}$ (A33)
$\displaystyle=\begin{pmatrix}(1+i\zeta_{R})(1+i\zeta_{L})e^{i\theta}+\zeta_{R}\zeta_{L}e^{-i\theta}&i(1-i\zeta_{L})\zeta_{R}e^{-i\theta}+i\zeta_{L}(1+i\zeta_{R})e^{i\theta}\\\
-i(1+i\zeta_{L})\zeta_{R}e^{i\theta}-i\zeta_{L}(1-i\zeta_{R})e^{-i\theta}&(1-i\zeta_{R})(1-i\zeta_{L})e^{-i\theta}+\zeta_{R}\zeta_{L}e^{i\theta}\end{pmatrix},$

from which the cavity transmission coefficient can be obtained as
$t=1/T_{22}$.

###  A5 Brownian noise model for vibrational relaxation

A particular case of dissipation arises in the case of a quantum harmonic
oscillator with free energy $\mathcal{H}=\nu(p^{2}+q^{2})/2$ interacting with
a heat bath of (mutually independent) harmonic oscillators. In the Brownian
motion model

$\displaystyle\dot{q}$ $\displaystyle=\nu p,$ (A34a) $\displaystyle\dot{p}$
$\displaystyle=-2\Gamma p-\nu q+\xi,$ (A34b)

the decay $2\Gamma$ stems from the correlations of the stochastic input noise
$\xi$ which is only affecting the momentum quadrature. This model is a
consequence of displacement-displacement interactions between the oscillator
and an infinite surrounding bath
$\mathcal{H}_{\text{int}}=-q\sum_{k}\alpha_{k}{Q}_{k}$ with coupling
coefficients $\alpha_{k}$ Weiss (1999). The correlations of the input noise
$\xi$ at a given bath temperature in the time domain read Genes _et al._
(2008)

$\braket{\xi(t)\xi(t^{\prime})}=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega(t-t^{\prime})}S_{\text{th}}(\omega).$
(A35)

where we defined the thermal noise spectrum $S_{\text{th}}(\omega)$

$S_{\text{th}}(\omega)=\frac{2\Gamma\omega}{\nu}\left[\coth{\frac{\hbar\omega}{2k_{\text{B}}T}}+1\right].$
(A36)

While generally, the noise is not delta-correlated, one can approximate Eq.
(A35) in the limit of large temperatures as (using that in this limit
$\bar{n}+1/2\approx k_{B}T/(\hbar\nu)$ with the average thermal occupancy
$\bar{n}=\left[\mathrm{exp}(\hbar\nu/k_{B}T)-1\right]^{-1}$ and
$\coth(x)\approx 1/x$ for $x\ll 1$)

$\braket{\xi(t)\xi(t^{\prime})}\approx
2\Gamma(2\bar{n}+1)\delta(t-t^{\prime})+i\frac{2\Gamma}{\nu}\delta^{\prime}(t-t^{\prime}),$
(A37)

where $\delta^{\prime}(t-t^{\prime})$ is the time derivative of the delta
function. From the spectrum evaluated at $\pm\nu$ one can obtain the cooling
and heating contributions: $S_{\text{th}}(\nu)=4\Gamma(\bar{n}+1)$ and
$S_{\text{th}}(-\nu)=4\Gamma\bar{n}$. Generally, from these properties one can
always estimate the damping rate as
$2\Gamma=(S_{\text{th}}(\nu)-S_{\text{th}}(-\nu))/2$ and the equilibrium
thermal occupancy as
$\bar{n}=S_{\text{th}}(-\nu)/[S_{\text{th}}(\nu)+S_{\text{th}}(-\nu)]$. In the
case of small temperatures $k_{B}T/(\hbar\nu)\ll 1$, the thermal spectrum
becomes very asymmetric
$S_{\text{th}}(\omega)=\left[4\omega\Gamma/\nu\right]\theta(\omega)$ as it
vanishes for negative frequencies (we introduced the Heaviside function
$\theta(\omega)$). Generally, one can proceed with a Fourier domain analysis
of the steady state of the equations above which leads to
$p(\omega)=\epsilon(\omega)\xi(\omega)$ and $q(\omega)=i\nu p(\omega)/\omega$
with the mechanical susceptibility describing the response of the phonon bath

$\displaystyle\epsilon(\omega)=\frac{i\omega}{\omega^{2}-\nu^{2}+2i\Gamma\omega}.$
(A38)

One can also make use of the fact that the noise is always delta-correlated in
Fourier domain
$\braket{\xi(\omega)\xi\omega^{\prime})}=S_{\text{th}}(\omega)\delta(\omega+\omega^{\prime})$.
One can then proceed by computing the time-domain correlations (using that
$\epsilon(-\omega)=\epsilon^{*}(\omega)$)

$\displaystyle\braket{p(t)p(t^{\prime})}$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega
e^{-i\omega(t-t^{\prime})}|\epsilon(\omega)|^{2}S_{\text{th}}(\omega),$ (A39a)
$\displaystyle\braket{q(t)q(t^{\prime})}$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\omega
e^{-i\omega(t-t^{\prime})}\frac{\nu}{\omega}|\epsilon(\omega)|^{2}S_{\text{th}}(\omega).$
(A39b)

Generally, these expressions can be solved by a contour integral in the
complex plane. We note that the integral of the momentum correlation function
diverges and generally requires the introduction of a cutoff frequency
$\Lambda>\nu$ to keep it finite. However, for $\Gamma\ll\nu$, the expression
for $|\epsilon(\omega)|^{2}S_{\text{th}}(\omega)$ is sharply peaked around
$\pm\nu$. In this case, one can expand around the poles $\omega=\nu+\delta$
and $\omega=\nu-\delta$ with $|\delta|\ll|\nu|$ and only keep leading order
terms in $\delta$. The integrals in Eq. (A39) can then be approximated as the
Fourier transforms of Lorentzian lineshapes which is much easier to solve and
results in (for $t>t^{\prime}$)

$\displaystyle\braket{p(t)p(t^{\prime})}=\left[\left(\bar{n}+\frac{1}{2}\right)\cos(\nu\tau)-\frac{i}{2}\sin(\nu\tau)\right]e^{-\Gamma|\tau|},$
(A40)

with $\tau=t-t^{\prime}$.

###  A6 Derivation of the FRET rate

Let us evaluate the term $2\Omega\Im\braket{\sigma_{A}^{\dagger}\sigma_{D}}$
which is the one that give rise to energy transfer between the two molecules.
Formal integration of the equation of motion for the acceptor gives

$\displaystyle\sigma_{A}^{\dagger}(t)=\sigma_{A}(0)e^{-(\gamma_{A}-i\omega_{A})t}\mathcal{D}_{A}(0)\mathcal{D}_{A}^{\dagger}(t)+\int_{0}^{t}dt^{\prime}e^{-(\gamma_{A}-i\omega_{A})(t-t^{\prime})}\left[i\Omega\sigma_{D}^{\dagger}(t^{\prime})+\sqrt{2\gamma}\sigma_{A}^{\dagger,\text{in}}(t^{\prime})\right]\mathcal{D}_{A}(t^{\prime})\mathcal{D}_{A}^{\dagger}(t).$
(A41)

Due to the quick vibrational relaxation, we can make a great simplification by
neglecting the backaction of the acceptor’s dipole moment onto the donor and
assume free evolution for the donor

$\displaystyle\dot{\tilde{\sigma}}_{D}\approx-\left(\gamma_{D}+i\omega_{D}\right)\tilde{\sigma}_{D}+\sqrt{2\gamma}\tilde{\sigma}_{D}^{\text{in}}.$
(A42)

Integrating this, one can calculate the term
$2\Omega\Im\braket{\sigma_{A}^{\dagger}\sigma_{D}}$ which requires the
correlation function (assuming $t\geq t^{\prime}\geq 0$ and $t\gg
1/\Gamma_{D}$)

$\displaystyle\braket{{\cal{D}}_{D}(0){\cal{D}}_{D}^{\dagger}(t^{\prime}){\cal{D}}_{D}(t){\cal{D}}_{D}^{\dagger}(0)}=e^{-S_{D}}e^{S_{D}e^{-(\Gamma_{D}-i\nu_{D})(t-t^{\prime})}}e^{-S_{D}e^{-(\Gamma_{D}+i\nu_{D})t^{\prime}}}e^{S_{D}e^{-(\Gamma_{D}-i\nu_{D})t^{\prime}}},$
(A43)

where the four-operator correlation function can be reduced to a two-operator
correlation function by commuting ${\cal{D}}_{D}^{\dagger}(0)$ with
${\cal{D}}_{D}(t)$ and ${\cal{D}}_{D}^{\dagger}(t^{\prime})$ and using that
${\cal{D}}_{D}(0){\cal{D}}^{\dagger}_{D}(0)=\mathds{1}$. Under the assumption
of $\Gamma_{D},\Gamma_{A}\gg\gamma_{D},\gamma_{A}$, one can finally arrive at
an expression for the energy transfer rate showing proportionality to the
donor population
$2\Omega\braket{\sigma_{A}^{\dagger}\sigma_{D}}=\kappa_{\text{ET}}P_{D}(t)$
(assuming $\gamma_{D}=\gamma_{A}$):

$\displaystyle\kappa_{\text{ET}}=\sum_{n_{D},n_{A}}\frac{2\Omega^{2}f_{S_{D}}^{n_{D}}f_{S_{A}}^{n_{A}}(n_{D}\Gamma_{D}+n_{A}\Gamma_{A})}{(n_{D}\Gamma_{D}+n_{A}\Gamma_{A})^{2}+[\Delta-
n_{D}\nu_{D}-n_{A}\nu_{A}]^{2}}.$ (A44)

The denominator of the above expression asks that the resonance condition
$\omega_{D}-n_{D}\nu_{D}=\omega_{A}+n_{A}\nu_{A}$ is fulfilled, i.e.,
fluorescence lines of the donor $\ket{e_{D},0}\to\ket{g_{D},n_{D}}$ have to
overlap with absorption lines of the acceptor
$\ket{g_{A},0}\to\ket{e_{A},n_{A}}$, both of which are weighted by the
respective Franck-Condon factors. In the case of many vibrational modes $n$
for donor and acceptor, we can generalize the result by writing general
displacements $\mathcal{D}_{A}=\prod_{k=1}^{n}\mathcal{D}_{A}^{k}$ and
$\mathcal{D}_{D}=\prod_{k=1}^{n}\mathcal{D}_{D}^{k}$ for all vibrational
modes. The equations of motion can then be expressed in the same form as in
Eqs. (120)

$\displaystyle\dot{\tilde{\sigma}}_{D}=-\left(\gamma+i\omega_{D}\right)\tilde{\sigma}_{D}-i\Omega\tilde{\sigma}_{A}\mathcal{D}_{A}\mathcal{D}_{D}^{\dagger}+\sqrt{2\gamma}\tilde{\sigma}_{D}^{\text{in}},$
(A45a)
$\displaystyle\dot{\tilde{\sigma}}_{A}=-\left(\gamma+i\omega_{A}\right)\tilde{\sigma}_{A}-i\Omega\tilde{\sigma}_{D}\mathcal{D}_{D}\mathcal{D}_{A}^{\dagger}+\sqrt{2\gamma}\tilde{\sigma}_{A}^{\text{in}},$
(A45b)

except that now the displacement is produced by all modes. We assume that the
two molecules have the same vibrational properties
$\braket{\mathcal{D}_{A}^{k}(t^{\prime})\mathcal{D}_{A}^{k,\dagger}(t)}=\braket{\mathcal{D}_{D}^{k}(t^{\prime})\mathcal{D}_{D}^{k,\dagger}(t)}$.
This leads to a generalized energy transfer rate

$\displaystyle\kappa_{\text{ET}}=\sum_{\\{m_{k}=0\\}}^{\infty}\sum_{\\{l_{k}=0\\}}^{\infty}2\prod_{k=1}^{n}e^{-2S_{k}}\frac{S_{k}^{(m_{k}+l_{k})}}{m_{k}!l_{k}!}\frac{\sum_{k=1}^{n}(m_{k}+l_{k})\Gamma_{k}\Omega^{2}}{[\sum_{k=1}^{n}(m_{k}+l_{k})\Gamma_{k}]^{2}+[\omega_{D}-\omega_{A}-\sum_{k=1}^{n}(m_{k}+l_{k})\nu_{k}]^{2}}.$
(A46)

Here the sums go over all indices $\\{m_{k}\\}=m_{1},\ldots,m_{n}$ and
$\\{l_{k}\\}=l_{1},\ldots,l_{n}$ referring to the possible occupancies of all
$n$ vibrational modes of the two involved molecules.