††thanks: Present Address: Alpine Quantum Technologies, 6020 Innsbruck,
Austria

# Cooperative quantum phenomena in light-matter platforms

Michael Reitz  Max Planck Institute for the Science of Light, Staudtstraße 2,
D-91058 Erlangen, Germany  Department of Physics, University of Erlangen-
Nuremberg, Staudtstraße 7, D-91058 Erlangen, Germany  Christian Sommer  Max
Planck Institute for the Science of Light, Staudtstraße 2, D-91058 Erlangen,
Germany  Claudiu Genes  Max Planck Institute for the Science of Light,
Staudtstraße 2, D-91058 Erlangen, Germany  Department of Physics, University
of Erlangen-Nuremberg, Staudtstraße 7, D-91058 Erlangen, Germany

###### Abstract

Quantum cooperativity is evident in light-matter platforms where quantum-
emitter ensembles are interfaced with confined optical modes and are coupled
via the ubiquitous electromagnetic quantum vacuum. Cooperative effects can
find applications, among other areas, in topological quantum optics, in
quantum metrology or in quantum information. This tutorial provides a set of
theoretical tools to tackle the behavior responsible for the onset of
cooperativity by extending open quantum system dynamics methods, such as the
master equation and quantum Langevin equations, to electron-photon
interactions in strongly coupled and correlated quantum-emitter ensembles. The
methods are illustrated on a wide range of current research topics such as the
design of nanoscale coherent-light sources, highly reflective quantum
metasurfaces or low intracavity power superradiant lasers. The analytical
approaches are developed for ensembles of identical two-level quantum emitters
and then extended to more complex systems where frequency disorder or vibronic
couplings are taken into account. The relevance of the approach ranges from
atoms in optical lattices to quantum dots or molecular systems in solid-state
environments.

###### pacs:

42.50.Nn, 42.50.Ct, 42.50.-p

###### Contents

  1. I Introduction
  2. II Cooperativity of light and matter
    1. A Collective radiative emission
    2. B Superradiance and subradiance
    3. C The single excitation subspace
    4. D The collective Bloch sphere: Dicke superradiance
  3. III Subwavelength quantum emitter arrays
    1. A Band structure and topology of 1D chains
    2. B Applications of quantum-emitter rings and chains
    3. C Optical response of 2D subwavelength mirrors
    4. D Further remarks
  4. IV Cooperativity in cavity QED
    1. A Cavity QED with coupled quantum emitters
    2. B Input-output formalism for operators
    3. C Cavity QED with disordered ensembles
  5. V Applications in cavity QED
    1. A Antiresonance spectroscopy with 1D or 2D arrays
    2. B Hybrid cavities with 2D subwavelength mirrors
    3. C Superradiant lasers
    4. D Further remarks
  6. VI Quantum optics with molecules
    1. A Optical response of a single molecule
    2. B Near-field coupled molecules
    3. C Further remarks
  7. VII Conclusions
  8. A1 Master equations for open-system dynamics
  9. A2 Functional dependence of $F(kR)$, $G(kR)$
  10. A3 Lattice sum of spherical waves
  11. A4 One-dimensional transfer-matrix theory
  12. A5 Brownian noise model for vibrational relaxation
  13. A6 Derivation of the FRET rate

##  I Introduction

Some of the most intriguing phenomena in nature, both in the classical and
quantum domains, are a product of cooperative effects, i.e., they cannot be
understood by sole consideration of the individual constituents as they arise
from the interplay among them. While at the fundamental science level, the
understanding of many of those problems in the quantum domain already poses a
great intellectual challenge, there is an ever increasing interest to build,
control and harness complex cooperative platforms for emerging quantum
technologies Deutsch (2020).
Light-matter platforms provide an optimal playground for the observation and
exploitation of quantum cooperative effects. Quantum light, either multimode,
as naturally arising in the quantum electromagnetic vacuum or single mode, as
confined in the small volume of an optical resonator, can induce strong
interactions among quantum emitters (QEs). Cooperativity then occurs e.g. in
free space under high-density conditions and manifests itself in a strongly
modified material response owed to a continuous scattering and rescattering of
photons between the matter constituents. Two main aspects brought on by the
common coupling of an ensemble of $\mathcal{N}$ quantum emitters to an
electromagnetic environment are dipole-dipole interactions, stemming from a
virtual exchange of photons and collective radiative emission, stemming from
the loss of excitation into the infinite number of the electromagnetic vacuum
modes. While the former is included as a coherent effect, the latter is an
incoherent one, observable either as an increase (superradiance) or a decrease
(subradiance) in the collective spontaneous emission rate compared to that of
an isolated emitter.
Free-space near-field dipole-dipole interactions are utilized in structured
subwavelength arrays of quantum emitters where they allow for the hopping of
surface excitations. Such arrays are ideal platforms for achieving strong
light-matter interactions and high fidelities for photon storage capabilities
Facchinetti _et al._ (2016); Asenjo-Garcia _et al._ (2017). Furthermore, these
platforms can allow for the propagation of photons or excitations to be
protected against disorder and scattering caused by defects, in an approach
dubbed topological quantum optics Perczel _et al._ (2017); Bettles _et al._
(2017); Perczel _et al._ (2020). An advantage over linear topological photonic
systems is the intrinsic nonlinearity of the quantum emitters which could lead
to a rich many-body physics dynamics. Individual addressing of a single qubit
emitter has been suggested via quantum spin lenses Glaetzle _et al._ (2017)
and an approach to quantum networking with composite quantum systems comprised
of many atomic arrays has been proposed building on theoretical results
showing the production of a Bell entangled superposition quantum state for two
distant arrays Guimond _et al._ (2019).
Free-space subradiance properties are an important resource for applications
ranging from quantum-information processing Chang _et al._ (2018) to metrology
Ostermann _et al._ (2013); Facchinetti _et al._ (2016); Manzoni _et al._
(2018). Such cooperative dissipative effects have been extensively studied,
especially in the direction of engineering robust many-particle quantum states
characterized by extremely long lifetimes Asenjo-Garcia _et al._ (2017); Zhang
and Mølmer (2020). Subradiance has been experimentally employed to show near-
unity ultrathin reflectors Rui _et al._ (2020) with potential applications in
nano-optomechanics Shahmoon _et al._ (2019, 2020), antiresonance spectroscopy
Plankensteiner _et al._ (2017, 2019) and nonlinear quantum optics Bettles _et
al._ (2020).
Cavity quantum electrodynamics (cavity QED) Haroche and Kleppner (1989);
Berman (1994); Walther _et al._ (2006); Haroche and Raimond (2013) can
additionally mediate and enhance emitter-emitter interactions by interfacing
quantum emitters with confined optical resonances. The famous model of Dicke
superradiance Dicke (1954); Gross and Haroche (1982), showing the quick burst
of spontaneous emission from $\mathcal{N}$ indistiguishable quantum emitters
can be simulated in the context of cavity QED and finds application in the
operation of superradiant lasers Bohnet _et al._ (2012); Norcia _et al._
(2016). The combination of subwavelength arrays with plain mirrors allows for
the design of hybrid cavities with one or two frequency-dependent end mirrors,
superior in performance to ones made up by plain, frequency-insensitive
mirrors Wild _et al._ (2018); Cernotik _et al._ (2019).
While many light-matter platforms are designed at the level of quasi-pure
quantum emitters, i.e., considering only electron-photon interactions, a
variety of extremely promising directions in quantum engineering are utilizing
more complex quantum emitters such as organic molecules, quantum dots or color
centers. In particular at the level of organic molecules, the strong
confinement of light modes in micro-cavities gives rise to a novel research
direction into vacuum-dressed or cavity-dressed materials, i.e., vacuum-
hybridized materials with enhanced properties. At the level of mesoscopic
systems, changes in charge conductivity Orgiu _et al._ (2015); Schachenmayer
_et al._ (2015); Feist and Garcia-Vidal (2015); Hagenmüller _et al._ (2017,
2018), energy transfer rates Zhong _et al._ (2016, 2017), chemical reactivity
Hutchison _et al._ (2012); Galego _et al._ (2017a); Herrera and Spano (2016);
Galego _et al._ (2017b); Martinez-Martinez _et al._ (2018) have been
experimentally observed and theoretically studied. At the level of single
molecules, the focus is in producing reliable single quantum emitters as
single photon sources with applications in entanglement generation or in
optical quantum computing Lettow _et al._ (2010). Recent results have shown
that a cavity-dressed molecule can act as an almost ideal quantum emitter,
exhibiting a closed electronic transition Wang _et al._ (2019a). For two or
more closely spaced molecules, cooperativity is naturally manifested in the
process of Förster resonance energy transfer based on energy exchange via
near-field dipole-dipole interactions followed by quick vibrational relaxation
Förster (1948); Blankenship (2014), a key process in photosynthetic light
harvesting.
Many of the aforementioned applications can be understood within the formalism
of open quantum system dynamics extended to the more complex problem of
correlated matter, such as occurring in a coupled quantum emitter ensemble. To
this end, this tutorial utilizes two competing, interconnected approaches, one
at the level of the density operator time evolution, i.e., the master equation
(ME) and the other following time dynamics of system operators, i.e., the
quantum Langevin equations (QLEs) approach. As the two formalisms are standard
textbook methods Breuer and Petruccione (2002); Weiss (1999); Scully and
Zubairy (1997); Gardiner and Zoller (2004); Vogel and Welsch (2006); Walls and
Milburn (2012), this tutorial proceeds with more complex aspects of light-
matter interactions such as the emergence of a cooperative ME for coupled
quantum-emitter systems in free space, introduced in Sec. A together with the
consequential occurrence of subradiance and superradiance, which are tackled
in Sec. B. Particular aspects such as the single excitation subspace and Dicke
superradiance Dicke (1954) are discussed in Sec. C and Sec. D, respectively.
The master equation allows for the derivation of equations of motion for
$\mathcal{N}$ mutually coupled emitters, utilized in understanding on how
energy dispersion relations and energy band gaps in one-dimensional (1D)
arrays can be engineered (in Sec. A). Also in 1D chain and ring
configurations, applications in quantum metrology, quantum information and
lasing are presented in Sec. B. The optical response of 2D subwavelength
arrays is derived in Sec. C to show the perfect reflection of incoming light
around certain collective resonances.
Fundamental aspects of cavity QED with correlated quantum emitters are
introduced in Sec. A followed by the input-output theory for operators
detailed in Sec. B. The question of frequency disorder, relevant in the case
of more complex quantum-emitter ensembles affected, for example, by
inhomogeneous broadening, is tackled in Sec. C. Applications of the formalism
introduced in Sec. IV are described in Sec. A on 1D arrays towards
antiresonance spectroscopy applications and in 2D highly reflective arrays
employed in hybrid cavities exhibiting Fano-like narrow resonances. Finally,
for three-dimensional mesoscopic ensembles of invertable two-level systems,
the theory of Dicke superradiance is applied to the characterization of
superradiant lasers inside lossy cavities.
The transition to increase the complexity of the considered quantum emitter
systems is undertaken in Sec. A by providing a derivation for electron-vibron
interactions in the form of a Holstein Hamiltonian Holstein (1959), and
emphasizing the application of QLEs at the level of polarons to characterize
molecular emission and absorption Reitz _et al._ (2019, 2020). The effect of
cooperative behavior as mediated by dipole-dipole interactions is tackled in
Sec. B to analytically characterize the FRET (Förster resonance energy
transfer) process as well as for the description of the physics of molecular
dimers Diehl _et al._ (2014), a standard model used in quantum chemistry.

##  II Cooperativity of light and matter

We provide in this section a quick review of collective effects in an ensemble
of free space standing $\mathcal{N}$ quantum emitters with emphasis on the
dipole-dipole interactions and collective decay mediated by the
electromagnetic vacuum. The emergent properties of subradiance, superradiance
and dipole-dipole-induced collective energy shifts are then introduced under
the weak excitation assumption (weak driving). Eventually, dynamics in the
full extended Hilbert space of the $\mathcal{N}$ emitters is treated by
introducing the collective Bloch sphere and the resulting Dicke superradiance
regime. This formalism is the starting point for the understanding of
applications exemplified in Sec. III and Sec. V.

###  A Collective radiative emission

The standard system we consider in the next sections is that of an ensemble of
$\mathcal{N}$ identical quantum emitters modeled for simplicity as
hydrogenlike atoms with nuclei fixed in random positions $\mathbf{R}_{j}$,
where $j=1,\ldots,\mathcal{N}$. The (single) electronic degree of freedom for
each emitter is described by quantized momentum $\hat{\mathbf{p}}_{j}$ and
position $\hat{\mathbf{r}}_{j}$ (relative to its respective nucleus) operators
appearing in the Hamiltonian
$h_{j}=\hat{\mathbf{p}}_{j}^{2}/(2\mu)+V(\hat{\mathbf{r}}_{j})$ consisting of
the kinetic energy (where the reduced mass is denoted by $\mu$) and an
electrostatic potential $V(\hat{\mathbf{r}}_{j})$. Out of the infinite set of
eigenvectors of $h_{j}$, we pick the lowest energy one (and set its energy to
zero) $\ket{g_{j}}$ and assume that all the physics described in the following
involves transitions to only one excited state $\ket{e_{j}}$ at energy
$\omega_{0}$ (we set $\hbar=1$). This amounts to a two-level system
approximation where the unity of the Hilbert space is a sum of only two
projectors $\mathds{1}_{j}=\ket{g_{j}}\bra{g_{j}}+\ket{e_{j}}\bra{e_{j}}$. To
quantify transitions between the two levels, we introduce the standard ladder
(Pauli) operators $\sigma_{j}=\ket{g_{j}}\bra{e_{j}}$,
$\sigma^{\dagger}_{j}=\ket{e_{j}}\bra{g_{j}}$. Their commutator
$[\sigma_{j}^{\dagger},\sigma_{j}]=\sigma_{z}^{j}$ is the population
difference operator. The free Hamiltonian can then be written as
${h}_{j}=\omega_{0}\sigma_{j}^{\dagger}\sigma_{j}$ or alternatively as
${h}_{j}=\omega_{0}\sigma_{z}^{j}/2$ (up to a constant energy shift).
To describe the interaction of the emitters with the electromagnetic vacuum,
one introduces a fictitious perfectly reflecting box of volume
$\mathcal{V}=\ell^{3}$ and follows a standard quantization procedure for the
electromagnetic field imposing periodic boundary conditions Breuer and
Petruccione (2002); Weiss (1999); Scully and Zubairy (1997); Gardiner and
Zoller (2004); Vogel and Welsch (2006); Walls and Milburn (2012); Loudon
(1973). This leads to a plane-wave expansion of the electric field operator

$\hat{\mathbf{E}}(\mathbf{R})=\textstyle\sum_{\mathbf{k}}\mathcal{E}_{k}\left(a_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}}+a^{\dagger}_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}\right){\bm{\epsilon}}_{\mathbf{k}},$
(1)

where the allowed $k$ vectors are multiples of $2\pi/\ell$ on each Cartesian
direction and the index $\mathbf{k}$ runs over all possible $k$ vectors and
also over the two orthogonal polarizations with unit vectors
${\bm{\epsilon}}_{\mathbf{k}}$. The bosonic operators follow the commutation
$[a_{\mathbf{k}},a_{\mathbf{k^{\prime}}}^{\dagger}]=\delta_{\mathbf{k}\mathbf{k}^{\prime}}$
and their action is to create and destroy excitations in a given field mode.
The zero-point electric field amplitude for a given mode with frequency
$\omega_{k}=ck$ is defined as
$\mathcal{E}_{k}=\sqrt{\omega_{k}/(2\epsilon_{0}\mathcal{V})}$ (where $c$ is
the speed of light and $\epsilon_{0}$ is the vacuum permittivity). The
Hamiltonian for the field inside the box can then be written as a sum over an
infinite number of bosonic modes
$\mathcal{H}_{\text{vac}}=\sum_{\mathbf{k}}\omega_{k}(a_{\mathbf{k}}^{\dagger}a_{\mathbf{k}}+1/2)$.

The coupling of light and matter occurs within the formalism of the minimal
coupling Hamiltonian Breuer and Petruccione (2002); Weiss (1999); Scully and
Zubairy (1997); Gardiner and Zoller (2004); Vogel and Welsch (2006); Walls and
Milburn (2012); Loudon (1973)

$\mathcal{H}=\sum_{j=1}^{\mathcal{N}}\frac{[\hat{\mathbf{p}}_{j}-e\hat{\mathbf{A}}(\mathbf{R}_{j}+\hat{\mathbf{r}}_{j})]^{2}}{2\mu}+V(\hat{\mathbf{r}}_{j})+\mathcal{H}_{\text{vac}},$
(2)

describing the total quantum system of charges and electromagnetic vacuum
modes. Notice that the canonical momentum of the electron is shifted to
$\hat{\mathbf{p}}_{j}=\mu\dot{\hat{\mathbf{r}}}_{j}+e\hat{\mathbf{A}}(\mathbf{R}_{j}+\hat{\mathbf{r}}_{j})$
where $\hat{\mathbf{A}}$ is the vector potential operator of the
electromagnetic field. The observation that the size of the electronic orbital
is much smaller than the wavelength $\lambda_{0}=2\pi c/\omega_{0}$ associated
with the optical transition between the two electronic orbitals allows for a
simplified form of the Hamiltonian

$\mathcal{H}=\sum_{j=1}^{\mathcal{N}}h_{j}+\sum_{j=1}^{\mathcal{N}}\hat{\mathbf{d}}_{j}\cdot\hat{\mathbf{E}}(\mathbf{R}_{j})+\mathcal{H}_{\text{vac}}.$
(3)

In this so-called dipole approximation, the light-matter interaction
Hamiltonian
$\mathcal{H}_{\text{int}}=\textstyle{\sum_{j=1}^{\mathcal{N}}\hat{\mathbf{d}}_{j}\cdot\hat{\mathbf{E}}(\mathbf{R}_{j})}$
involves only the dipole moment operator
$\hat{\mathbf{d}}_{j}=-e\hat{\mathbf{r}}_{j}$ and the electric field operator
evaluated at the position of the nuclei. Under the two-level assumption, the
dipole operator is written as
$\hat{\mathbf{d}}_{j}=\mathbf{d}_{\text{eg}}\sigma_{j}+\text{h.c.}$, where the
transition dipole moment is computed between the two states
$\mathbf{d}_{\text{eg}}=\bra{g_{j}}\hat{\mathbf{d}}_{j}\ket{e_{j}}$ (assuming
identical emitters). Moreover, one performs an additional rotating-wave
approximation (RWA) where fast oscillating terms $\sigma_{j}a_{\mathbf{k}}$
and $\sigma_{j}^{\dagger}a_{\mathbf{k}}^{\dagger}$ (oscillating under free
evolution with $\omega_{0}+\omega_{k})$ are neglected. Under these
approximations, the light-matter interaction part of the Hamiltonian can be
written as

$\mathcal{H}_{\text{int}}=\sum_{\mathbf{k}}\sum_{j=1}^{\mathcal{N}}\left(g_{\mathbf{k}}\sigma_{j}^{\dagger}a_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}_{j}}+g^{*}_{\mathbf{k}}a^{\dagger}_{\mathbf{k}}\sigma_{j}e^{-i\mathbf{k}\cdot\mathbf{R}_{j}}\right),$
(4)

describing photon-emitter energy exchanges at rates
$g_{\mathbf{k}}=\mathcal{E}_{k}{\bm{\epsilon}}_{\mathbf{k}}\cdot\mathbf{d}_{\text{eg}}$.
Notice that this interaction Hamiltonian conserves the number of excitations
in the system.

Figure 1: Schematics of the procedure followed to obtain the open-system
dynamics formulation for a system of $\mathcal{N}$ quantum emitters (denoted
by $\sigma_{j}$ and $\sigma_{j^{\prime}}$) subject to the quantum
electromagnetic vacuum. The key point is that, for a pair of emitters,
coupling to the common reservoir modes (denoted by $a_{\mathbf{k}}$ and
$a_{\mathbf{k}^{\prime}}$) gives rise to separation-dependent dipole-dipole
interactions at rate $\Omega_{jj^{\prime}}$ as well as to collective decay
rates $\gamma_{jj^{\prime}}$.

The Hamiltonian $\mathcal{H}$ governs unitary evolution in an infinite-
dimensional Hilbert space (owing to the infinite number of electromagnetic
modes). One could then, in principle, use a deterministic Schrödinger equation
or equivalently the von Neumann equation for the density operator
$\rho_{\text{tot}}(t)$ to compute the state of the system at any time.
However, this is an extremely complex task requiring an immense computational
resource even for small size systems. Instead, in order to drastically reduce
complexity, an open-system dynamics approach is usually followed, which
consists in reducing the system of interest to the Hilbert space of dimension
$2^{\mathcal{N}}$ of the matter part only (as schematically illustrated in
Fig. 1). A main observation allowing such reduction of complexity is that the
time-dependent light-matter interaction Hamiltonian can be simply expressed in
the interaction picture as
$\mathcal{H}_{\text{int}}=\sum_{j=1}^{\mathcal{N}}[\sigma_{j}\mathcal{F}_{j}^{\dagger}(t)+\sigma_{j}^{\dagger}\mathcal{F}_{j}(t)]$
describing an excitation exchange between the emitters and the field, with the
time-dependent operator

$\mathcal{F}_{j}(t)=\sum_{\mathbf{k}}g_{\mathbf{k}}a_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{R}_{j}}e^{-i(\omega_{k}-\omega_{0})t},$
(5)

acting solely in the Hilbert space of the photon modes.
The open-system dynamics procedure is then based on a weak-coupling assumption
and implies that the electromagnetic degrees of freedom are traced over to
obtain an equation of motion for
$\rho(t)=\text{Tr}_{\text{em}}\left[\rho_{\text{tot}}(t)\right]$. The steps
for this procedure are summarized in Appendix A1 closely following the
standard textbook derivation in Ref. Breuer and Petruccione (2002) but also
widely covered in, among others, Refs. Weiss (1999); Scully and Zubairy
(1997); Gardiner and Zoller (2004); Vogel and Welsch (2006); Walls and Milburn
(2012). In the final master equation all system properties are derived from
two-time correlations of the operators $\mathcal{F}_{j}(t)$. For $\mathcal{N}$
quantum emitters in the electromagnetic vacuum at zero temperature, the
resulting collective quantum master equation reads

$\dot{\rho}(t)=i\left[\rho(t),\sum_{j=1}^{\mathcal{N}}h_{j}+\mathcal{H}_{dd}\right]+\mathcal{L}_{e}[\rho].$
(6)

The last term in the equation above describes irreversible loss of excitation
into the electromagnetic vacuum and is expressed in the form of a
superoperator (an operator acting on density operators) defined as Lehmberg
(1970)

$\mathcal{L}_{e}[\rho]=\sum_{j,j^{\prime}=1}^{\mathcal{N}}\gamma_{jj^{\prime}}\left[2\sigma_{j}\rho(t)\sigma^{\dagger}_{j^{\prime}}-\sigma^{\dagger}_{j}\sigma_{j^{\prime}}\rho(t)-\rho(t)\sigma^{\dagger}_{j}\sigma_{j^{\prime}}\right].$
(7)

Notice that this is not in standard, diagonal Lindblad form Breuer and
Petruccione (2002); Weiss (1999); Scully and Zubairy (1997); Gardiner and
Zoller (2004); Vogel and Welsch (2006); Walls and Milburn (2012) which, for a
collapse operator $\mathcal{O}$ and collapse rate $\gamma_{\mathcal{O}}$, is
defined as

$\mathcal{L}_{\gamma}[\rho]=\gamma_{\mathcal{O}}\left[2\mathcal{O}\rho(t)\mathcal{O}^{\dagger}-\mathcal{O}^{\dagger}\mathcal{O}\rho(t)-\rho(t)\mathcal{O}^{\dagger}\mathcal{O}\right],$
(8)

and describes decay at rate $\gamma_{\mathcal{O}}$ through a single channel
with operator $\mathcal{O}$.
The decay rates are $\gamma_{jj^{\prime}}=(3\gamma/2)F(k_{0}R_{jj^{\prime}})$
where the relative distance vector is
$R_{jj^{\prime}}=|\mathbf{R}_{j}-\mathbf{R}_{j^{\prime}}|$ and the function of
position is defined as

$\displaystyle
F(kR)=\left[1+\frac{(\mathbf{e}_{\mathbf{d}}\cdot\nabla_{\mathbf{R}})^{2}}{k^{2}}\right]\frac{\sin(kR)}{kR},$
(9)

with unit vector $\mathbf{e}_{\mathbf{d}}$ in the direction of
$\mathbf{d}_{\text{eg}}$. The rate
$\gamma=\omega^{3}_{0}d_{\text{eg}}^{2}/(6\pi c^{3}\epsilon_{0})$ is the
spontaneous emission rate of a single emitter (note that the decay rate of the
excited-state population is given by $2\gamma$). The collective radiative
decay properties are given by the oscillatory behavior of $F(kR)$ as
illustrated in Fig. 2b. Notice that, for emitters separated by much more than
$\lambda_{0}$, the decay becomes purely diagonal as expected for
noninteracting, independent emitters.
The coherent term $\mathcal{H}_{dd}$ in Eq. (6) instead describes a dipole-
dipole interaction characterized by a virtual exchange of excitation via the
vacuum modes without loss of photons

$\mathcal{H}_{dd}=\sum_{j,j^{\prime}:j\neq
j^{\prime}}^{\mathcal{N}}\Omega_{jj^{\prime}}\sigma^{\dagger}_{j}\sigma_{j^{\prime}}.$
(10)

The dipole-dipole exchange rate
$\Omega_{jj^{\prime}}=-(3\gamma/2)G(k_{0}R_{jj^{\prime}})$ can be obtained
from the same function that characterizes the collective decay rates via the
definition

$\displaystyle G(k_{0}R)=\frac{c}{\pi\omega_{0}^{3}}\mathcal{P}\int
dk\frac{(ck)^{3}}{ck-\omega_{0}}F(kR),$ (11)

where $\mathcal{P}$ denotes the Cauchy principal value. The explicit
functional dependence of $F(kR)$ and $G(kR)$ on distance is detailed in
Appendix A2. The behavior illustrated in Fig. 2b shows that the dipole-dipole
interaction ceases at large distances, as expected, but diverges at close
separations. This is, however, only an artefact of the initial assumptions
that the dipole-electric field interaction is valid at any interparticle
distance. This is of course not true, as for separations on the order of the
size of the orbitals, one ends up with a fundamentally quantum many-body
problem where the tunneling of electrons between neighboring emitters has to
be taken into account (leading to molecule formation, hybridization of
orbitals, etc). More involved models, based, for example, on a quantum
electrodynamics density-functional formalism, can then be employed Flick _et
al._ (2017a).

###  B Superradiance and subradiance

Figure 2: (a) Decay and eigenstates in the collective two-emitter basis.
Superradiant and subradiant decay channels emerge through the formation of
symmetric and antisymmetric superposition states. (b) Dipole-dipole coupling
$\Omega_{12}$ and collective decay rate $\gamma_{12}$ scaling with normalized
interatomic distance $k_{0}a$ for parallel ($\parallel$) and perpendicular
($\perp$) orientation of dipoles (with respect to the connection vector
$\mathbf{R}_{12}$). (c) Possible experimental realization of strong dipole-
dipole couplings for two identical chromophores connected by an insulating
bridge of length much smaller than $\lambda_{0}$. The green arrows represent
the orientation of the transition dipoles. Adapted from Ref. Diehl _et al._
(2014).

The superoperator in Eq. (7) describes nontrivial mutual decay characterized
by the matrix $\mathbf{\Gamma}$ made up by the elements
$\gamma_{jj^{\prime}}$. However, this term is not in standard diagonal form
(characterized by a single collapse operator) as it is not comprised of
$\mathcal{N}$ independent decay channels. One can perform a basis
transformation with a matrix T (such that $\textbf{T}^{-1}=\textbf{T}^{\top}$)
which diagonalizes $\bm{\Gamma}$ such that
$\text{diag}\left(\tilde{\gamma}_{1},...,\tilde{\gamma}_{\mathcal{N}}\right)=\textbf{T}^{\top}\bm{\Gamma}\textbf{T}$
where $\tilde{\gamma}_{k}$ is the $k$th eigenvalue of the decay matrix. In
order to achieve the Lindblad form defined in Eq. (8) one can define a set of
collapse operators $\Pi_{k}=\sum_{j}T_{jk}\sigma_{j}$ such that

$\displaystyle\mathcal{L}_{e}[\rho]$
$\displaystyle=\sum_{k=1}^{\mathcal{N}}\tilde{\gamma}_{k}\left(2\Pi_{k}\rho\Pi_{k}^{\dagger}-\Pi_{k}^{\dagger}\Pi_{k}\rho-\rho\Pi_{k}^{\dagger}\Pi_{k}\right),$
(12)

now describes $\mathcal{N}$ independent decay channels each with an associated
collapse operator $\Pi_{k}$ and associated loss rate $\tilde{\gamma}_{k}$.
Notice that the preparation of a collective quantum superposition that decays
at one of the rates $\tilde{\gamma}_{k}$ would require the application of the
Hermitian conjugate of the corresponding collapse operator $\Pi_{k}^{\dagger}$
to the collective ground state
$\ket{G}=\ket{g}_{1}\otimes\ket{g}_{2}\otimes...\otimes\ket{g}_{\mathcal{N}}$
of the system. In practice this is, however, not a straightforward task.

Two coupled emitters \- Let us first make the connection between collective
radiative effects and the symmetry of quantum superpositions by considering
the simplest example of two quantum emitters separated by distance $a$.
Diagonalization of the decay matrix leads to superradiant/subradiant decay
channels characterized by $\gamma_{\pm}=\gamma\pm\bar{\gamma}$ and at the same
time renders the Hamiltonian in diagonal form with eigenenergies
$\omega_{\pm}=\omega_{0}\pm\Omega$. The eigenstates are
symmetric/antisymmetric superpositions
$\ket{\pm}=(\ket{e}_{1}\otimes\ket{g}_{2}\pm\ket{g}_{1}\otimes\ket{e}_{2})/\sqrt{2}$.
This is illustrated in Fig. 2a in the collective basis where the other two
states are the fully excited one $\ket{E}=\ket{e}_{1}\otimes\ket{e}_{2}$ and
the ground state $\ket{G}=\ket{g}_{1}\otimes\ket{g}_{2}$. The splitting
between the two levels in the collective basis $2\Omega=2\Omega_{12}(a)$ as
well as the magnitude and the sign of the mutual decay rate
$\bar{\gamma}=\gamma_{12}(a)$ depend strongly on the particular choice of
dipole orientations as well as on separation (as illustrated in Fig. 2b).
However, for distances below half a wavelength $a<\lambda_{0}/2$, the
antisymmetric (symmetric) states always have a subradiant (superradiant)
character. One can then understand the connection between state symmetry and
radiative properties in terms of an destructive (constructive) interference of
radiative paths.
The strong scaling of the near-field dipole-dipole coupling with the
interparticle distance renders such a simple system valuable for experimental
applications in superresolution imaging. For example, experimentally it has
been shown that optical resolution of fluorescent molecules can be achieved at
distances as small as $12\,\mathrm{nm}$ Hettich _et al._ (2002). In addition,
for strongly coupled emitters, chemical or mechanical means can be employed to
correct energy shifts and render the closely spaced emitters indistinguishable
such that a source of indistinguishable photons can be achieved. Finally, we
remark that subwavelength separations (even at the level of less than
$10\,\mathrm{nm}$) could be achieved in assembled molecular dimers Diehl _et
al._ (2014) where two chromophores are coupled via an insulating bridge
(illustrated in Fig. 2c). The connection between the dipole-dipole exchange
for two-level systems and a vibronic dimer model for chromophores is detailed
in Sec. B.

###  C The single excitation subspace

While the simultaneous diagonalization of the Lindblad term and the
Hamiltonian is generally not possible for $\mathcal{N}>2$, one can still get
some intuition in the nature of cooperative decay in the particular case of
$\mathcal{N}$ equally spaced quantum emitters in a 1D chain configuration. We
start by analyzing the eigenstates and eigenvalues of the dipole-dipole
Hamiltonian under the nearest-neighbor approximation and only up to a single
excitation shared in the whole system (reducing the Hilbert space from
$2^{\mathcal{N}}$ to only contain $\mathcal{N}+1$ states). We then follow with
an analytical derivation of the dissipation rates of such collective states.
As a consequence of the fact that the Lindblad term (i.e., the decay matrix
$\bm{\Gamma}$) and the Hamiltonian do not commute, dissipation is not diagonal
in the collective basis.

Eigenstates of the dipole-dipole Hamiltonian \- Let us first inspect the
Hamiltonian in the analytically solvable case where a nearest-neighbor (NN)
approximation is performed. This is justified as the dipole-dipole
interactions scales as $R_{jj^{\prime}}^{-3}$ for distances smaller than a
wavelength and thus the nearest neighbor coupling is almost an order of
magnitude larger than that of the next to nearest neighbor one. The resulting
Hamiltonian is then in the form of a tridiagonal Toeplitz matrix

$\displaystyle\mathcal{H}=\omega_{0}\sum_{j}\sigma_{j}^{\dagger}\sigma_{j}+\sum_{j}\Omega(\sigma_{j}^{\dagger}\sigma_{j+1}+\sigma_{j+1}^{\dagger}\sigma_{j}),$
(13)

where $\Omega$ is the coupling between two neighbors. The free Hamiltonian has
eigenenergies of degeneracy
$C_{n}^{\mathcal{N}}=\mathcal{N}!/[(\mathcal{N}-n)!n!]$ for a given level of
energy $n\omega_{0}$ where $n$ ranges from $0$ for the ground state to
$\mathcal{N}$ for the highest excited state. We illustrate this in Fig. 3a
where we plot the eigenvalues of the Hamiltonian up to the second excitation
manifold. The reduced single-excitation manifold on which we focus is
characterized by the particle basis $\ket{j}=\sigma_{j}^{\dagger}\ket{G}$ with
$j=1,\ldots,\mathcal{N}$ and the collective ground state $\ket{G}$. The
resulting eigenenergies of the Hamiltonian are

$\epsilon_{k}=\omega_{0}+2\Omega\cos{\left[\frac{\pi
k}{\mathcal{N}+1}\right]},$ (14)

for an index $k$ running from $1$ to $\mathcal{N}$ (and the trivial energy of
the ground state is $0$). The corresponding eigenstates of the Hamiltonian are
then described in the collective exciton basis

$\displaystyle\tilde{\ket{k}}=\sum_{j}\sqrt{\frac{2}{\mathcal{N}+1}}\sin{\left(\frac{\pi
kj}{\mathcal{N}+1}\right)}\ket{j}=\sum_{j}f_{jk}\ket{j}.$ (15)

Figure 3: (a) Energy scheme of the collective $\mathcal{N}$-emitter system.
The $C_{n}^{\mathcal{N}}$-fold degeneracy of a given $n$-excitation manifold
is lifted by the dipole-dipole interactions and contains states with
superradiant (red) and subradiant (blue) character. (b) Decay rates of most
superradiant $\gamma_{\text{max}}$ and most subradiant state
$\gamma_{\text{min}}$ for an equidistant chain with lattice constant
$a=\lambda_{0}/4$ as a function of $\mathcal{N}$ in the single-excitation
manifold. The decay rate of the most subradiant state shows a scaling with
$\mathcal{N}^{-3}$ (light blue circles). From the diagonalization of the decay
matrix only, one would obtain an exponential scaling of the subradiant state
with $\mathcal{N}$ (dark blue stars). (c) Entanglement of six-emitter
equidistant chain with separation $a=0.1\lambda_{0}$ which is initialized in
the most subradiant state $\tilde{\ket{\mathcal{N}}}$. The logarithmic
negativity of each emitter with respect to the $\mathcal{N}-1$ other emitters
is plotted for different points in time.

Notice that the reverse transformation is straightforward
$\ket{j}=\sum_{k}f_{jk}\tilde{\ket{k}}$. The transformation also holds at the
level of operators $\tilde{\sigma}_{k}=\sum_{j}f_{jk}\sigma_{j}$ (in the
single excitation manifold) which sees the diagonal form of the Hamiltonian as
$\mathcal{H}=\sum_{k}\epsilon_{k}\tilde{\sigma}_{k}^{\dagger}\tilde{\sigma}_{k}$.

Collective dissipation \- To compute the decay rates of the Hamiltonian
eigenstates, we use the master equation to derive the equation of motion for
the population component $\rho_{kk}=\bra{\tilde{k}}\rho\ket{\tilde{k}}$. We
arrive at

$\dot{\rho}_{kk}=-\sum_{k^{\prime}}\left\\{\sum_{jj^{\prime}}\gamma_{jj^{\prime}}f_{jk}f_{j^{\prime}k^{\prime}}\right\\}(\rho_{kk^{\prime}}+\rho_{k^{\prime}k}),$
(16)

which shows that the eigenstates of the Hamiltonian couple within the first
excitation manifold in addition to the decay to the ground state. The diagonal
elements can then be estimated by setting $k=k^{\prime}$ to derive
$\tilde{\gamma}_{k}=\sum_{jj^{\prime}}\gamma_{jj^{\prime}}f_{jk}f_{j^{\prime}k}$
and more explicitly

$\displaystyle\tilde{\gamma}_{k}=\frac{2}{\mathcal{N}+1}\sum_{jj^{\prime}}\gamma_{jj^{\prime}}\sin{\left(\frac{\pi
kj}{\mathcal{N}+1}\right)}\sin{\left(\frac{\pi
kj^{\prime}}{\mathcal{N}+1}\right)}.$ (17)

From this expression, one can derive a scaling of subradiant states with
roughly $\mathcal{N}^{-3}$ (see lower panel of Fig. 3b for the scaling of the
most subradiant state) which will be utilized for cavity antiresonance
spectroscopy applications in Sec. A. More generally, scalings with
$\mathcal{N}^{-s}$ up to $s=5$ can be reached Zhang and Mølmer (2020). Also,
superradiance emerges, which for small distances $a$ is characterized by a
rate proportional to $\mathcal{N}$ but eventually saturates with increasing
$\mathcal{N}$ (around $\mathcal{N}=\lambda_{0}/a$, as illustrated in Fig. 3b).
Notice that generally, for $k\neq k^{\prime}$, the dissipative couplings
$\textstyle\sum_{jj^{\prime}}\gamma_{jj^{\prime}}f_{jk}f_{j^{\prime}k^{\prime}}$
are not vanishing, which is the direct consequence of the fact that the
Lindblad term and Hamiltonian do not commute.
Up to here we have only considered the dissipative properties of the
eigenstates of the Hamiltonian. Instead, one can directly inspect the Lindblad
term by a diagonalization of the decay matrix in Eq. (12). Numerical results
in Fig. 3b indicate that in this case subradiant states are present which are
characterized by decay rates scaling exponentially down with $\mathcal{N}$. As
these states are not eigenstates of the Hamiltonian, they are not easily
addressable. The targeting scheme would then involve some excitation scheme
that applies the collective operator $\Pi^{\dagger}_{k}$ to the ground state
$\ket{G}$. In Sec. B we will show how a combination of a magnetic field
gradient together with a pulsed-laser excitation could instead be used to
perform such an action.
In Sec. A we provide an alternative treatment in terms of excitations
propagating on the 1D support by hopping between neighboring sites via the
dipole-dipole exchange. The discrete index $k$ then describes the
quasimomentum of an individual collective mode and its location with respect
to the light cone distinguishes between superradiant versus subradiant modes.

Entanglement properties \- The collective eigenstates of the dipole-dipole
Hamiltonian (both super- and subradiant ones) commonly feature non-classical
correlations Plankensteiner _et al._ (2015); Hebenstreit _et al._ (2017),
rendering them as an interesting resource for quantum information processing;
highly subradiant states are of course even more useful due to the increased
lifetime of correlations. To this end we analyze the logarithmic negativity
Plenio (2005), which is an entanglement monotone. For a bipartite system
consisting of the subsystems $A$ and $B$, it is defined as

$\displaystyle E_{N}(\rho)$
$\displaystyle=\log_{2}\left(|\rho^{T_{A}}|\right),$ (18)

where $\rho^{T_{A}}$ denotes the partial transpose with respect to the
subsystem $A$ and $|\cdot|$ is the tracenorm. In Fig. 3c, we show the time
evolution of the state with $k=\mathcal{N}$. At distinct time points, we
compute the logarithmic negativity for each emitter (i.e., we choose our
bipartite system to consist of the $i$th emitter and the rest of the chain).
One can see, that the amount of bipartite entanglement is significantly larger
in the center of the chain, even in the initial state. Over time, this
behavior is retained and correlation is only slowly lost due to excitation
loss of the chain. Even at $t=100\gamma^{-1}$ there still is considerable
entanglement in the system.

###  D The collective Bloch sphere: Dicke superradiance

Let us now go beyond the single-excitation manifold and consider a famous
example introduced by Dicke Dicke (1954) which shows the generation of a
superradiant pulse from an ensemble of indistinguishable quantum emitters. The
model assumes an idealized case of $\mathcal{N}$ quantum emitters within a
very small volume and neglects their dipole-dipole interactions. While this is
per se an unrealistic assumption (a densely packed ensemble of QEs exhibits
large dipole-dipole frequency shifts), we indicate later in Sec. C how this
model can be realized in the context of cavity QED and is relevant for the
physics of lasing in what is known as bad cavity superradiant lasers Bohnet
_et al._ (2012).
In the following we will make use of the Bloch-sphere illustration in Fig. 4a,
where $\mathcal{N}$ two-level systems can be described by a collective spin of
length $\mathcal{N}/2$. This follows from $1/2\otimes 1/2\otimes...\otimes
1/2=0\oplus 1\oplus...\oplus\mathcal{N}/2$ for $\mathcal{N}$ even and
$1/2\otimes 1/2\otimes...\otimes 1/2=1/2\oplus
3/2\oplus...\oplus\mathcal{N}/2$ for $\mathcal{N}$ odd (the tensor product
space of dimensions $2^{\mathcal{N}}$ of $\mathcal{N}$ spins $1/2$ can be
written as a direct sum of spaces of different sizes). We can then introduce
collective spin operators $S_{z}=\textstyle\sum_{j}\sigma_{z}^{j}/2$ and
$S=\textstyle\sum_{j}\sigma_{j}$ and the total spin
$\mathbf{S}^{2}=S_{z}^{2}+(S^{\dagger}S+SS^{\dagger})/2$. As in the standard
angular momentum algebra, it is then possible to find a collective basis,
known as the Dicke basis of $\mathbf{S}^{2}$ and $S_{z}$, indexed by two
quantum numbers $\ket{s,m_{s}}$

$\displaystyle\mathbf{S}^{2}\ket{s,m_{s}}$
$\displaystyle=s(s+1)\ket{s,m_{s}},$ (19a) $\displaystyle S_{z}\ket{s,m_{s}}$
$\displaystyle=m_{s}\ket{s,m_{s}}.$ (19b)

Figure 4:  (a) Cross section of collective $\mathcal{N}$-emitter Bloch sphere
(for an odd number of spins) showing an onionlike structure with Bloch vector
(red) pointing towards the north pole representing the fully excited state
$\ket{E}$. The lines of latitude (blue dotted lines) correspond to the
$\ket{\mathcal{N}/2,m_{\mathcal{N}/2}}$ states. (b) Dicke superradiance for
$\mathcal{N}=100$ emitters initialized in the fully excited state $\ket{E}$:
The top shows time dynamics of $\braket{S_{z}}$ while bottom shows the
normalized decay rate
$\gamma_{\text{sup}}=-\partial_{t}\braket{S_{z}}/(\braket{{S}_{z}}+\mathcal{N}/2)$.
The inset shows the time evolution of $\braket{S^{\dagger}S}$ with
superradiant ‘burst’ reaching a maximum when $\braket{S_{z}}=m_{s}$ reaches
zero. The dashed gray lines show the result for pure exponential decay of
independent, uncoupled emitters.

The spin quantum number $s$ can take either integer or half-integer values and
is bounded by $0,1/2\leq s\leq\mathcal{N}/2$ where the zero holds for even and
the $1/2$ holds for odd values of $\mathcal{N}$. The so-called inversion
quantum number $m_{s}$ which measures the projection of the collective spin
onto the $z$ axis can take values from $-s\leq m_{s}\leq s$ and equals
$-\mathcal{N}/2$ in the ground state $\ket{G}=\ket{g}^{\otimes\mathcal{N}}$
and $\mathcal{N}/2$ in the fully excited state
$\ket{E}=\ket{e}^{\otimes\mathcal{N}}$. The ladder operators $S$ and
$S^{\dagger}$ act on the Dicke states as

$\displaystyle S^{\dagger}\ket{s,m_{s}}$
$\displaystyle=\sqrt{(s-m_{s})(s+m_{s}+1)}\ket{s,m_{s}+1},$ (20a)
$\displaystyle S\ket{s,m_{s}}$
$\displaystyle=\sqrt{(s+m_{s})(s-m_{s}+1)}\ket{s,m_{s}-1}.$ (20b)

In a closed system interacting with the electromagnetic field, the
interactions are mediated by $S$ and $S^{\dagger}$. Thereby the selection
rules for optical transitions are given by $\Delta m_{s}=\pm 1$, $\Delta s=0$
and the Hilbert space splits into non-interacting subspaces defined by the
quantum number $s$ with dimension $2s+1$ (illustrated in Fig. 4a).
Let us now describe the phenomenon of Dicke superradiance in the ideal case
where all mutual decay rates are equal and maximal
$\gamma_{jj^{\prime}}=\gamma$. A trivial observation is now that the Lindblad
decay term in Eq. (7) assumes a very simple form with a single superradiant
decay channel and collective collapse operator

$\displaystyle\mathcal{L}_{e}[\rho]$ $\displaystyle=\gamma\left[2S\rho
S^{\dagger}-S^{\dagger}S\rho-\rho S^{\dagger}S\right].$ (21)

We will follow the time evolution under such a Lindblad term for an initially
fully excited ensemble characterized by the state vector
$\ket{\mathcal{N}/2,\mathcal{N}/2}$. Notice that the action of the collapse
operator cannot take the system out of the $s=\mathcal{N}/2$ symmetric
manifold containing $\mathcal{N}+1$ Dicke states. We can derive an equation of
motion for the population difference operator average
$\braket{\dot{S}_{z}}=-2\gamma\braket{S^{\dagger}S}$ which simply states that
the loss rate is proportional to the emitted intensity. For a single quantum
emitter one has $\braket{\sigma^{\dagger}\sigma}=\braket{\sigma_{z}}/2+1/2$,
implying that the spontaneous emission always follows an exponential law. The
population difference of a collective Dicke state $\ket{s,m_{s}}$ decays
instead at a state-dependent rate
$2\gamma\braket{S^{\dagger}S}=2\gamma(s+m_{s})(s-m_{s}+1)$. For the initially
fully excited state, the decay rate is $2\gamma\mathcal{N}$ which is the same
as expected for $\mathcal{N}$ independently decaying emitters. However,
particle-particle correlations start building up during the evolution and by
the time the $\ket{\mathcal{N}/2,0}$ state (for even $\mathcal{N}$) is
reached, the emission is superradiant and scales approximately as
$\gamma\mathcal{N}^{2}/2$ (see Fig. 4b). Notice that the crucial effect of
correlations can be distinguished by rewriting
$\braket{S^{\dagger}S}=\braket{\sum_{i}\sigma_{i}^{\dagger}\sigma_{i}}+\braket{\sum_{i\neq
j}\sigma_{i}^{\dagger}\sigma_{j}}$. Using that the sum over populations is
given by $\braket{\sum_{i}\sigma_{i}^{\dagger}\sigma_{i}}=s+m_{s}$, the
dipole-dipole correlation between emitters $i$ and $j$ can be estimated as
$\braket{\sigma_{i}^{\dagger}\sigma_{j}}=(s^{2}-m_{s}^{2})/(\mathcal{N}(\mathcal{N}-1))$,
reaching a maximum of approximately $1/4$ for $m_{s}=0$ and becoming zero for
$m_{s}=\pm s$.

##  III Subwavelength quantum emitter arrays

The interplay between dipole-dipole interactions and collective radiance in
quantum-emitter ensembles leads to a multitude of applications of 1D and 2D
subwavelength arrays such as in nonlinear quantum optics Bettles _et al._
(2020); Wild _et al._ (2018), nano-optomechanics Shahmoon _et al._ (2020), the
design of quantum metamaterials with magnetic response at optical frequencies
Alaee _et al._ (2020); Ballantine and Ruostekoski (2020a), as platforms for
quantum information processing Glaetzle _et al._ (2017); Guimond _et al._
(2019); Bekenstein _et al._ (2020) or as chiral light-matter interfaces
Grankin _et al._ (2018). These applications are based on the fact that such
structures can support collective surface resonances that can interact in a
controllable fashion with impinging fields.

We will start by studying the dispersion relations of the surface modes on 1D
platforms by means of a Bloch ansatz, showing the occurrence of band gaps and
Dirac points. This direction has recently emerged showing the promise of
subwavelength arrays for topological quantum optics implementations Perczel
_et al._ (2017); Bettles _et al._ (2017). We then proceed by providing
illustrations on 1D emitter chains and rings aimed at showing the usefulness
of subradiance as a resource for: i) improved frequency sensitivity (as shown
in Refs. Ostermann _et al._ (2013, 2014)), ii) robust quantum memories
Plankensteiner _et al._ (2015) and iii) the design of nanoscale coherent light
sources as recently proposed in Ref. Holzinger _et al._ (2020). Finally, we
derive the optical response of two dimensional subwavelength arrays around
certain confined surface-mode resonances Bettles _et al._ (2016); Shahmoon _et
al._ (2017); Bettles _et al._ (2020) and describe a regime recently
experimentally tackled showing close to unity reflectivity for arrays of
optically trapped atoms Rui _et al._ (2020).

###  A Band structure and topology of 1D chains

Figure 5: (a) Energy band structure for a simplified model of a chain of
nearest neighbor coupled identical emitters, for two different alignments of
dipoles (perpendicular and parallel to the chain). (b) Top: band structure
considering the exact expression of near-field coupling $\Omega_{jj^{\prime}}$
and in the presence of collective radiative decay $\gamma_{jj^{\prime}}$.
Bottom: corresponding decay rate as a function of the quasi momentum $q$. The
light cone characterized by $\omega=c|q|$, here plotted at the two level
resonance $\omega_{0}=ck_{0}$ is indicated by the dashed black lines.
Parameters are $a/\lambda_{0}=0.2$ with $\mathcal{N}=200$ as in Asenjo-Garcia
_et al._ (2017). (c) Band structure for a chain of equidistant emitters with
two different alternating species. A band gap of size
$|\omega_{1}-\omega_{2}|$ around
$\omega_{\text{av}}=(\omega_{1}+\omega_{2})/2$ is opened. (d) Additionally,
when the interaction between the emitters has an alternating sign, we obtain a
Dirac point at $q=0$ (for the case of degenerate frequencies
$\omega_{1}=\omega_{2}$).

.

Up to this point, we have solved for collective resonances and their
associated collective radiation rates by a direct diagonalization of the
Hamiltonian with near-field coupling terms. One can however take a solid-state
approach instead where the array provides a crystalline structure for the
quasiexcitations propagating on its surface. As the simplest example, let us
first revisit the one-dimensional equidistant chain of $\mathcal{N}$ emitters
with identical frequencies as already discussed in Sec. C described by the
Hamiltonian

$\mathcal{H}=\omega_{0}\sum_{j}\sigma^{\dagger}_{j}\sigma_{j}+\sum_{j^{\prime}\neq
j}\Omega_{jj^{\prime}}\sigma^{\dagger}_{j}\sigma_{j^{\prime}}.$ (22)

We first assume nearest neighbor (NN) coupling
$\Omega_{jj^{\prime}}=\Omega\delta_{j,j^{\prime}\pm 1}$ which insures
translational symmetry with periodicity $a$ (equivalently stated, the unit
cell contains one site only). We impose periodic boundary conditions (PBCs) to
simulate the mesoscopic case by asking that at the edge
$\Omega_{1\mathcal{N}}=\Omega$. We then ask the question: what kind of
excitations can propagate in this chain, i.e., what kind of dispersion
relations such excitations will exhibit. To this end, we start with the
equations of motion for the expectation values $\beta_{j}=\braket{\sigma_{j}}$
which can be obtained from
$\braket{\dot{\sigma}_{j}}=\mathrm{Tr}[\sigma_{j}\dot{\rho}]$ where
$\dot{\rho}=-i[\mathcal{H},\rho]$ (neglecting radiative emission in a first
step)

$\dot{\beta}_{j}=-i\omega_{0}\beta_{j}-i\Omega\beta_{j-1}-i\Omega\beta_{j+1}.$
(23)

We furthermore assume a weak excitation limit
$\braket{{\sigma}_{z}^{j}}\approx-1$. In a compact matrix formulation, we
rewrite $\dot{\mathbf{v}}=M\mathbf{v}$ where
$\mathbf{v}=(\beta_{1},\cdots,\beta_{\mathcal{N}})^{\top}$ and the drift
matrix is expressed as a circulant Toeplitz matrix

$M=-i\left(\begin{array}[]{cccccc}\omega_{0}&\Omega&0&\cdots&0&\Omega\\\
\Omega&\omega_{0}&\Omega&\cdots&0&0\\\ 0&\Omega&\omega_{0}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&\cdots&\omega_{0}&\Omega\\\
\Omega&0&0&\cdots&\Omega&\omega_{0}\end{array}\right).$ (24)

The direct diagonalization of the matrix above gives the expected
$\mathcal{N}$ collective modes listed in Sec. C with a state index from $1$ to
$\mathcal{N}$. However, we now instead look for a dispersion relationship with
a quasimomentum $q$, to which end we plug the ansatz
$\beta_{j}=A_{q}e^{i(qaj-\omega t)}$ into Eq. (23) which straightforwardly
leads to the following dispersion relation

$\omega(q)=\omega_{0}+2\Omega\cos(qa).$ (25)

From the application of the PBCs we have $\beta_{1}=\beta_{\mathcal{N}+1}$
which leads to $qa\mathcal{N}=2\pi m$ where $m$ is an integer. This indicates
that the allowed quasimomenta are $q=2\pi m/(\mathcal{N}a)$ and we can fix the
first Brillouin zone to $m\in\\{-\mathcal{N}/2,\mathcal{N}/2\\}$ corresponding
to $q$ from $-\pi/a$ (left propagating wave) to $\pi/a$ (right propagating
wave). The resulting dispersion relations for two arrangements of dipoles
(parallel and perpendicular to the chain axis) are indicated in Fig. 5a.

Beyond NN with dissipation \- The results presented above describe a
simplified model where only nearest-neighbor interactions are employed and
spontaneous emission is disregarded. We now proceed by analyzing the full
equations of motion including dissipation which are obtained from Eq. (6) and
read

$\dot{\beta}_{j}=-i(\omega_{0}-i\gamma)\beta_{j}-\sum_{j^{\prime}\neq
j}(i\Omega_{jj^{\prime}}+\gamma_{jj^{\prime}})\beta_{j^{\prime}}.$ (26)

These equations take into account the exact behavior of
$\Omega_{jj^{\prime}}$, $\gamma_{jj^{\prime}}$ (as defined in Sec. A and App.
A2) and lead to the following dispersion relation

$\displaystyle\omega(q)\\!-\\!i\gamma(q)\\!$ $\displaystyle=$
$\displaystyle\\!\omega_{0}\\!-\\!i\gamma+\\!\sum_{k\neq
j}(\Omega_{jk}\\!-\\!i\gamma_{jk})e^{iqa(k-j)}.$ (27)

The result is numerically illustrated in Fig. 5b and compared with the
approximation analytically obtained for NN coupling in Fig. 5a. The
corresponding collective decay is plotted in the lower panel of Fig. 5b. The
vertical lines indicate the location of the light cone, where the
quasimomentum of the excitation propagating on the surface equals the wave
vector of the photon $k_{0}=2\pi/\lambda_{0}$, which for the chosen distance
$a=\lambda_{0}/5$ leads to $q_{0}a/\pi=0.4$. Notice that the interpretation in
terms of the location with respect to the light cone is straightforward as
waves with quasimomentum larger (in absolute value) than that of the photon
cannot escape the chain and therefore are subradiant (in agreement with
results presented in Ref. Asenjo-Garcia _et al._ (2017)).

Emerging band gaps \- Let us now consider a slightly more complex structure
with alternating emitter frequencies $\omega_{1}$ and $\omega_{2}$ and
identical NN couplings at rate $\Omega$ (as illustrated in Fig. 5c). One then
immediately notices that this implies a lattice with a double unit cell. We
denote the two types of emitters by an upper index such that their amplitudes
are $(\beta^{(1)}_{j},\beta^{(2)}_{j})$ and the site index $j$ runs from $1$
to $\mathcal{N}/2$ (such that we keep the chain length at $\mathcal{N}a$). The
dipole-dipole interaction now couples the two kinds of emitters with each
other within the unit cell and also between neighboring unit cells leading to
the following set of coupled equations

$\displaystyle\dot{\beta}^{(1)}_{j}$
$\displaystyle=-i\omega_{1}\beta^{(1)}_{j}-i\Omega\beta^{(2)}_{j-1}-i\Omega\beta^{(2)}_{j},$
(28a) $\displaystyle\dot{\beta}^{(2)}_{j}$
$\displaystyle=-i\omega_{2}\beta^{(2)}_{j}-i\Omega\beta^{(1)}_{j+1}-i\Omega\beta^{(1)}_{j}.$
(28b)

We now ask for propagating waves with the following profile
$\beta^{(1)}_{j}=A_{q}e^{i(2qaj-\omega t)}$ and
$\beta^{(2)}_{j}=B_{q}e^{i(2qaj-\omega t)}$ and plug in this ansatz into the
above equations to result in the following eigenvalue problem

$\left(\begin{array}[]{cc}\omega_{1}-\omega(q)&\Omega\left(1+e^{-i2qa}\right)\\\
\Omega\left(1+e^{i2qa}\right)&\omega_{2}-\omega(q)\end{array}\right)\left(\begin{array}[]{c}A_{q}\\\
B_{q}\end{array}\right)=0.$ (29)

Two types of eigenvalues $\omega_{\pm}(q)$ result from diagonalization of the
above matrix corresponding to two distinct energy bands

$\omega_{\pm}(q)=\frac{\omega_{1}+\omega_{2}}{2}\pm\sqrt{\left(\frac{\omega_{1}-\omega_{2}}{2}\right)^{2}+4\Omega^{2}\cos^{2}(qa)}.$
(30)

The particularity of the system is that the two bands exhibit a band gap
$|\omega_{1}-\omega_{2}|$ at the edge of the first Brillouin zone (as depicted
in Fig. 5c). The PBC now impose that
$\beta^{(1)}_{1}=\beta^{(1)}_{\mathcal{N}/2+1}$ which leads to the condition
$q=\pi m/(\mathcal{N}a)$. The first Brillouin zone is now defined by
$m\in\\{-\mathcal{N}/2,\mathcal{N}/2\\}$ corresponding to $q$ varying between
$-\pi/(2a)$ to $\pi/(2a)$.

Dirac points \- Further complexity in the band structure can be engineering by
assuming alternating coupling strengths (for example from $\Omega$ to
$-\Omega$ for consecutive pairs as in Fig. 5d). Noticing that this situation
is again characterized by a double unit cell, we follow the same steps as
above to ask for two kind of propagating waves and reach the following
eigenvalue problem

$\left(\begin{array}[]{cc}\omega_{1}-\omega(q)&\Omega\left(1-e^{-i2qa}\right)\\\
\Omega\left(1-e^{i2qa}\right)&\omega_{2}-\omega(q)\end{array}\right)\left(\begin{array}[]{c}A_{q}\\\
B_{q}\end{array}\right)=0,$ (31)

which leads to the following dispersion relation for the two branches

$\omega_{\pm}(q)=\frac{\omega_{1}+\omega_{2}}{2}\pm\sqrt{\left(\frac{\omega_{1}-\omega_{2}}{2}\right)^{2}+4\Omega^{2}\sin^{2}(qa)}.$
(32)

The quasimomentum varies as above between $-\pi/(2a)$ to $\pi/(2a)$. The
resulting dispersion curve plotted in Fig. 5d shows the emergence of an
avoided crossing in the center of the Brillouin zone which exhibits a Dirac
point with linear dispersion relation at frequency degeneracy
$\omega_{1}=\omega_{2}$.

Berry phase \- A particularly interesting case occurs for identical frequency
emitters (transition frequency $\omega_{0}$) with alternating interaction
strengths $\Omega_{1}$ and $\Omega_{2}$. This can be mapped onto the Su-
Schrieffer-Heeger model (SSH) Su _et al._ (1979), which for electrons
describes the emergence of bulk insulating phases distinguished by topological
invariants and which in nature occurs in polyacetylene molecules.

Figure 6: (a) Chain with alternating dipole-dipole couplings forming the
geometry of the Su-Schrieffer-Heeger (SSH) model. (b) Two distinct bands
emerge, separated by a band gap $2|\Omega_{1}-\Omega_{2}|$ and with shape
independent of the relation $\Omega_{1}\lessgtr\Omega_{2}$. The vertical
dashed line shows the corresponding energies and eigenvectors corresponding to
a given parameter $q$. The arrow in the lower band shows the path chosen for
the integration to derive the accumulated phase from one end to the other one
of first Brillouin zone. (c) The accumulated phase on the lower energy branch
shows considerable differences between the case of $\Omega_{1}>\Omega_{2}$ and
that of $\Omega_{1}<\Omega_{2}$.

.

As before we write the equations of motion

$\displaystyle\dot{\beta}^{(1)}_{j}$
$\displaystyle=-i\omega_{0}\beta^{(1)}_{j}-i\Omega_{2}\beta^{(2)}_{j-1}-i\Omega_{1}\beta^{(2)}_{j},$
(33a) $\displaystyle\dot{\beta}^{(2)}_{j}$
$\displaystyle=-i\omega_{0}\beta^{(2)}_{j}-i\Omega_{2}\beta^{(1)}_{j+1}-i\Omega_{1}\beta^{(1)}_{j}.$
(33b)

and testing for the two kinds of propagating waves results in the following
eigenvalue problem

$\left(\begin{array}[]{cc}\omega_{0}-\omega&\Omega_{1}+\Omega_{2}e^{-i2qa}\\\
\Omega_{1}+\Omega_{2}e^{i2qa}&\omega_{0}-\omega\end{array}\right)\left(\begin{array}[]{c}A_{q}\\\
B_{q}\end{array}\right)=0.$ (34)

Here we will focus both on the energy dispersion curves but also on the
characteristics of their corresponding eigenstates. The eigenvalues are given
by

$\omega_{\pm}(q)=\omega_{0}\pm\sqrt{(\Omega_{1}-\Omega_{2})^{2}+4\Omega_{1}\Omega_{2}\cos^{2}(qa)}.$
(35)

While the band structure is insensitive to the exchange of $\Omega_{1}$ and
$\Omega_{2}$, the eigenvectors are not. The normalized eigenstates are
analytically expressed in vector form as

$\ket{\varphi_{\pm}(q)}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\pm
e^{-i\phi(q)}\\\ 1\end{array}\right),$ (36)

where the $q$-dependent phase is given by the following relation

$e^{-i\phi(q)}=\frac{\Omega_{1}+\Omega_{2}e^{-i2qa}}{\sqrt{\Omega_{1}^{2}+\Omega_{2}^{2}+2\Omega_{1}\Omega_{2}\cos{2qa}}}.$
(37)

The eigenvalue problem can now be written formally as
$\mathcal{H}(q)\ket{\varphi_{n}(q)}=\omega_{n}(q)\ket{\varphi_{n}(q)}$ (where
$n$ stands for $\pm$). We now proceed by assuming that a path $q(t)$ is taken,
with $q(t=0)=-\pi/(2a)$ and $q(\tau)=\pi/(2a)$; this is illustrated in Fig.
6b, particularized to the lower energy band. By moving adiabatically slow,
tunneling to the orthogonal eigenstate is not allowed. Generally, we can then
write the state of the system at any time $t$ as
$\ket{\psi(t)}=e^{-i\theta(t)}\ket{\varphi_{n}(q(t))}$ where the time-
dependent phase could be computed directly by the application of the time-
ordered evolution operator with a time-dependent Hamiltonian to the initial
state. However, a more elegant solution comes from simply writing the
Schrödinger equation
$\mathcal{H}(q(t))\ket{\psi(t)}=i\partial_{t}\ket{\psi(t)}$ explicitly, to
arrive at an equation for the phase

$\left(\partial_{t}\theta(t)\right)\ket{\varphi_{n}(q(t))}=\omega_{n}(q(t))\ket{\varphi_{n}(q(t))}-i\partial_{t}\ket{\varphi_{n}(q(t))}.$
(38)

We sandwich the equation above with $\bra{\varphi_{n}(q(t))}$ and integrate to
obtain two distinct contributions. The integral of the energy over the whole
band vanishes as the final and initial point are equal in energy. The second
contribution is path dependent and reads

$\theta(\tau)=-i\int^{\tau}_{0}\braket{\varphi_{n}(q(t))}{\partial_{t}}{\varphi_{n}(q(t))}dt.$
(39)

The integration over time can be turned into a path integration such that

$\displaystyle\theta=-i\int^{q(\tau)}_{q(0)}\braket{\varphi_{n}(q)}{\partial_{q}}{\varphi_{n}(q)}dq.$
(40)

This describes the path dependent Berry phase where
$A_{n}(q)=\braket{\varphi_{n}(q)}{\partial_{q}}{\varphi_{n}(q)}$ can be
identified as the Berry potential. The Berry phase is gauge invariant for
closed integration loops. In general, the temporal aspect of the adiabatic
deformation is a more or less fictitious process but it is helpful to unravel
the topological structure. Integrating the Berry phase for our Hamiltonian
$\mathcal{H}(q)$ over the Brillouin zone $(-\pi/(2a),\pi/(2a))$ for the lower
band results in

$\nu=\frac{1}{\pi}\oint
A_{-}(q)dq=\left\\{\begin{array}[]{cr}0&\Omega_{1}>\Omega_{2}\\\
1&\Omega_{1}<\Omega_{2}\\\
\mathrm{undef.}&\Omega_{1}=\Omega_{2}\end{array}\right.$ (41)

where the matrix element is easily computed from the vector expression of the
eigenstates in Eq. (36) to result in $A_{-}(q)=(1/2)d\phi(q)/dq$. For
differentiation we use the following expression for
$\phi(q)=\mathrm{arctan2}(\Omega_{1}+\Omega_{2}\cos(2qa),\Omega_{2}\sin(2qa))$
and we note that $\nu$ is similar to a winding number. The result of the
integral is a topological invariant showing that the two cases
$\Omega_{1}>\Omega_{2}$ and $\Omega_{1}<\Omega_{2}$ with the same band
structure outcome are topologically different which is expressed by the Berry
phase for a closed parameter path (see Fig. 6c). In order to go from one
topological phase $\nu$ in the bulk to another here given by smoothly varying
$\Omega_{1}$ and $\Omega_{2}$, one needs to cross a point
$(\Omega_{1}=\Omega_{2})$ where the band gap is zero which would violate
adiabaticity. This shows that the Hamiltonians for the two insulating phases
occurring for $\Omega_{1}>\Omega_{2}$ and $\Omega_{1}<\Omega_{2}$ are not
adiabatically equivalent which is manifested in the difference of the
topological invariant. In the case of open boundary conditions one finds that
the number of edge states forms a topological invariant as well Asbóth _et
al._ (2016).

###  B Applications of quantum-emitter rings and chains

Subradiance can be exploited towards applications in quantum metrology for
sensitive frequency detection as well as in quantum information for the
engineering of robust quantum memories. Moreover, subradiance in symmetric
arrangements, such as rings, can be utilized to design nanoscale light sources
acting as thresholdless nanolasers. To describe these effects, we make use
here of the collective master equation in the single excitation regime
introduced in Sec. C and of the Bloch sphere representation introduced in Sec.
D.

Figure 7: (a) Ramsey metrology with phased excitation (shown on the individual
spin Bloch sphere) for a chain of six emitters. Minimum sensitivity is plotted
against interrogation time $\tau$ showing improved sensitivity for $m=3$
(blue, dotted curve) as opposed to the independent case (black, dashed) or
symmetric illumination (red, full line). (b) Phased excitation achieved via
magnetic field gradient. For two coupled emitters, adiabatic population
transfer to the robust state is achieved via detuned illumination of the
symmetric state and subsequent transfer to the robust antisymmetric state via
the magnetic field. The population reaches values close to unity and decays
much slower than $e^{-2\gamma t}$. (c) Subwavelength structure
($a=\lambda_{0}/2$) of a ring resonator made of five quantum emitters
surrounding an incoherently pumped central emitter. The second-order
correlation function at zero delay $g^{(2)}(0)$ is presented as a function of
the pumping rate $\gamma_{\text{p}}$.

.

Quantum metrology \- Ramsey interferometry is routinely used in quantum
metrology for the most sensitive measurements of optical clock frequencies.
The conventional method of separated oscillatory fields Ramsey (1990) assumes
an ensemble of atoms initially in the ground state illuminated by a laser of
frequency $\omega_{\ell}$ in two consecutive steps separated by time $\tau$.
The two pulses are assumed to be instantaneous and tuned as $\pi/2$ pulses
that can be visualized as $\pi/2$ rotations around the $y$ axis on the
collective Bloch sphere of radius $\mathcal{N}/2$ (introduced in Sec. D).
After the interrogation time $\tau$ the population difference
$\braket{S_{z}}(\omega_{\ell})$ is monitored as a function of the laser
frequency. Its behavior is sinusoidal where the argument is the total
accumulated phase $(\omega_{0}-\omega_{\ell})\tau$ stemming from the mismatch
of the (drifting, variable) laser frequency $\omega_{\ell}$ and the constant
frequency separation of the emitter ground and excited states $\omega_{0}$.
Analysis of the monitored population difference curve indicates a minimal
sensitivity:

$\delta\omega=\min\left[\frac{\Delta
S^{z}(\omega,\tau)}{\left|\partial_{\omega}\langle
S^{z}\rangle(\omega,\tau)\right|}\right],$ (42)

where the minimization is performed with respect to
$\omega=\omega_{0}-\omega_{\ell}$. For weakly decaying emitters
($\gamma\ll\tau^{-1}$), the sensitivity is simply $1/(\tau\sqrt{\mathcal{N}})$
and can obviously be optimized by using longer interrogation times and more
emitters. However, longer interrogation times bring decay into play, while
operation at higher densities achieved by increasing $\mathcal{N}$
indefinitely are marked by the onset of cooperative effects such as dipole-
dipole shifts and super/subradiance.
For independently decaying emitters, the application of the master equation
introduced in Sec. II leads to a simple analytical solution for both
$\braket{S_{z}}$ and $\Delta S^{z}$ allowing one to estimate
$\left[\delta\omega\right]_{\text{indep}}=e^{\gamma\tau}/(\tau\sqrt{\mathcal{N}})$.
Further optimization with respect to the interrogation time gives an optimal
$\tau_{\text{opt}}=1/\gamma$ and optimal sensitivity
$e\gamma/\sqrt{\mathcal{N}}$, which shows that the main impediment of Ramsey
interferometry subject to radiative loss is the limited interrogation times
available. This expression hints towards a deterioration of the sensitivity in
the case of equally illuminated dense ensembles where symmetric collective
states are addressed which are typically characterized by superradiant
behavior.
To protect against such detrimental effects, Refs. Ostermann _et al._ (2013,
2014) introduced an alternative procedure which uses an additional step in the
Ramsey sequence aimed at hiding the collective states into decoherence-free
subspaces. To this end, one complements the $\pi/2$ pulse with a phase
distribution pulse, which for a particular atom $j$ is represented by a
rotation around the $z$-direction with the angle $\varphi_{j}^{(m)}=2\pi
m(j-1)/\mathcal{N}$, where $m=1,\ldots,[\mathcal{N}/2]$ and $[\mathcal{N}/2]$
is the integer before $\mathcal{N}/2$. Rotations in the single emitter
subspace are defined as
$\mathcal{R}^{(j)}_{\mu}[\varphi]=\exp\left(i\varphi\,\sigma^{j}_{\mu}/2\right)$
where $\mu\in\\{x,y,z\\}$ and $\sigma^{j}_{x}=\sigma_{j}+\sigma_{j}^{\dagger}$
and $\sigma^{j}_{y}=i(\sigma_{j}-\sigma_{j}^{\dagger})$. The first generalized
Ramsey pulse operator of such an asymmetric Ramsey technique is then

$\mathcal{R}_{1}=\,\bigotimes_{j}\mathcal{R}^{(j)}_{z}\left[\varphi_{j}^{(m)}\right]\cdot\mathcal{R}^{(j)}_{y}\left[\frac{\pi}{2}\right].$
(43)

The action of this phase distribution pulse is illustrated in Fig. 7a on the
single particle Bloch sphere. In the second step, after the free evolution
where robustness is now expected owing to the folding of the collective state
into the subradiant part of the Bloch sphere, at time $\tau$ the phase spread
is (instantaneously) reversed and a $\pi/2$ pulse follows leading to the
second generalized Ramsey pulse

$\mathcal{R}_{2}=\bigotimes_{j}\,\mathcal{R}^{(j)}_{y}\left[\frac{\pi}{2}\right]\cdot\mathcal{R}^{(j)}_{z}\left[-\varphi_{j}^{(m)}\right].$
(44)

Finally, detection takes place as before and the sensitivity is minimized with
respect to frequency to obtain results as shown in Fig. 7a. Here, the
sensitivity for symmetric illumination is worse than that for independent
emitters while the interrogation times and consequently the minimum frequency
sensitivity can be considerably increased for phased excitations with $m\neq
0$. Further results and analytical considerations can be found in Refs.
Ostermann _et al._ (2013, 2014).

Quantum memories \- We now restrict the discussion to the single-excitation
manifold of a 1D emitter chain. Here, we aim at targeting collective
subradiant states as they show both robustness against decoherence while
exhibiting multipartite entanglement. As pointed out in Sec. C, eigenstates of
the dipole-dipole Hamiltonian with lower energy typically exhibit subradiance.
To access such states, one can proceed by first selecting them with individual
addressing by tailoring the laser light amplitude to fit the shape of the
collective state one wishes to address. For example, for a collective state
$\tilde{\ket{k}}$ introduced in Sec. C one can provide geometrical matching
with the following driving Hamiltonian

$\mathcal{H}_{k}=\eta\sum_{j=1}^{\mathcal{N}}\sin{\left(\frac{\pi
kj}{\mathcal{N}+1}\right)}\left(\sigma
e^{i\omega_{\ell}t}+\sigma^{\dagger}e^{-i\omega_{\ell}t}\right).$ (45)

Moreover, imposing the condition for resonance by setting
$\omega_{\ell}=\epsilon_{k}$ ensures that states far enough from the desired
targeted one are only weakly populated. In Ref. Plankensteiner _et al._ (2015)
it is shown that enhanced lifetimes much larger than $(2\gamma)^{-1}$ can be
reached by such tailored excitation.
While tailored phase excitation with subwavelength resolution might pose great
challenges, alternative methods could be envisioned: for example, symmetric
addressing could be combined with the application of a magnetic field
gradient. The effect of magnetic field gradient applied along the direction of
the emitter chain is to progressively shift the excited state by a quantity
$\Delta_{j}=\Delta_{\text{B}}(j-1)$ from the first emitter with $j=1$ to the
end of the chain $j=\mathcal{N}$. During the duration $\tau$ of an applied
laser pulse, this amounts to a rotation around the $z$ axis of the Bloch
vector of each emitter by the angle $\Delta_{j}\tau$. For a conveniently
chosen duration $\tau$, the effective phase difference between neighboring
emitters can be controlled: for example with the choice
$\Delta_{\text{B}}\tau=\pi$ a completely antisymmetric superposition can be
constructed. For example, for two coupled emitters, Fig. 7b shows results for
pulsed adiabatic transfer of population from the ground state to the
antisymmetric state by a simultaneous off-resonant drive of the symmetric
state and the action of the magnetic field gradient. The population then ends
up in the protected antisymmetric state which then decays much slower than the
timescale defined by $(2\gamma)^{-1}$. More results and analytical
calculations for larger systems sizes can be found in Ref. Plankensteiner _et
al._ (2015).

Nanoscale coherent-light sources \- Subwavelength spaced ensembles of quantum
emitters in ring-like configurations resemble the structure of certain
biological light-harvesting complexes (LHCs), which have been shown to act as
an extremely efficient system of antennae in the photosynthetic process
Cogdell _et al._ (2006). While the complexity of such biological systems is
incredibly hard to tackle even at the computational level, a number of
theoretical works have proposed simplified phenomenological models where, for
example, the combined effects of disorder, vibronic coupling, electron-phonon
couplings, are included as time-dependent frequency shifts (in the
Hamiltonian) or as dephasing (as an additional Lindblad term) Bourne Worster
_et al._ (2019); Caycedo-Soler _et al._ (2017).
Inspired by such naturally occurring systems, a recent theoretical proposal
has shown the possibility of designing a thresholdless laser, i.e., a coherent
light source Holzinger _et al._ (2020) in the configuration depicted in Fig.
7d. Here, a central emitter acts as a gain medium and is pumped by incoherent
light at some rate $\gamma_{\text{p}}$. The role of the optical resonator (as
present in standard lasing systems) is then taken by the optical modes defined
by the geometry of the ring emitters (indexed by $1,\ldots,\mathcal{N}$). In
order to test the properties of the system (emission rate, coherence of light,
etc.) one can proceed to solve the master equation under the assumption that
the incoherent pumping can be modeled by a Lindblad term with a
$\sigma_{\text{p}}^{\dagger}$ collapse operator (that, opposed to the
spontaneous emission case, takes population from the ground state to the
excited state).
A first observation is that the physics of the system can be reduced to solely
the interaction of the symmetric mode
$\ket{\psi_{\mathrm{sym}}}=(1/\sqrt{\mathcal{N}})\sum_{j}\ket{j}=\sigma^{\dagger}_{\mathrm{sym}}\ket{G}$
of the ring resonator with the central atom. The reduced Hamiltonian is then

$\mathcal{H}_{\text{ring}}=\Omega_{\mathrm{sym}}\sigma^{\dagger}_{\mathrm{sym}}\sigma_{\mathrm{sym}}+\sqrt{\mathcal{N}}\Omega(\sigma^{\dagger}_{\mathrm{sym}}\sigma+\mathrm{h.c.}),$
(46)

where $\Omega_{\text{sym}}=\sum_{j=2}^{\mathcal{N}}\Omega_{1j}$ is the dipole
energy shift of the symmetric state and the decay of the system is governed by
$\mathcal{L}[\rho]=\mathcal{L}_{\gamma}[\rho]+\mathcal{L}_{\gamma_{\text{p}}}[\rho]+\mathcal{L}_{\mathrm{sym}}[\rho]$.
For favorable geometries, it is possible to obtain a subradiant decay rate for
the symmetric mode of the ring resonator. Additionally, such configurations
can reduce the dissipative coupling of the central atom, strengthening the
lossless transport of population from the central atom to the ring. We can now
analyze the emitted light properties by using the fact that the far field is
proportional to the sum of the dipole operators such that we can define the
normalized second-order correlation function at zero time delay

$g^{(2)}(0)=\frac{\sum_{ijkl}\braket{\sigma^{\dagger}_{i}\sigma^{\dagger}_{j}\sigma_{k}\sigma_{l}}}{\left|\sum_{mn}\braket{\sigma^{\dagger}_{m}\sigma_{n}}\right|^{2}}.$
(47)

Numerical simulations, plotted in Fig. 7c show a close to unity $g^{(2)}(0)$,
indicating a coherent state as an output of such a thresholdless laser. For
higher pumping rates, antibunching of light is expected. Further results and
considerations can be found in Refs. Holzinger _et al._ (2020).

###  C Optical response of 2D subwavelength mirrors

Figure 8: (a) Two dimensional quantum-emitter array illustrated as a regular
pattern of near-field coupled dipoles with lattice constant $a$. (b) The array
is placed in the $x$-$y$ plane and is illuminated by a plane wave with wave
number $k$. (c) Effective decay rate $\gamma_{\text{eff}}$ and frequency shift
$\Omega_{\text{eff}}$ of an array consisting of $20\times 20$ atoms as a
function of $a/\lambda_{0}$. The dotted black curve shows the approximation
$\gamma_{\text{eff}}=\gamma\frac{3}{4\pi}(\lambda_{0}/a)^{2}$ valid for
lattice constants $a<\lambda_{0}$. (d) Reflection $|r|^{2}$ and transmission
coefficient $|t|^{2}=|1+r|^{2}$ of an atomic array as a function of the
lattice constant $a$ (for resonant illumination of the emitters). (e)
Calculated intensity profile $I(y,z)$ (in arbitrary units) in the $y$-$z$
plane for an array of $12\times 12$ atoms (indicated by red dots) with lattice
constant $a=0.8\lambda_{0}$ which is resonantly illuminated by a plane wave.
(f) Targeted addressing of super- and subradiant states of a 2D array with
$a=0.2\lambda_{0}$ as a function of the number of emitters in the array
$\mathcal{N}$. As indicated in the inset, subradiant states are targeted by
illumination with antisymmetric phases $f_{nm}=(-1)^{n+m}$ in a checkerboard
pattern while superradiant states are obtained from uniform, symmetric
illumination. The effective decay rate of the subradiant states scales
proportional to $\mathcal{N}^{-5.5}$ as shown by the dashed line.

Let us now consider the situation depicted in Figs. 8a,b where a 2D periodic
array of quantum emitters is positioned in the $x$-$y$ plane. To derive the
optical response of such a structure, we consider excitation in the form of a
plane wave with wavenumber $k=\omega_{\ell}/c=2\pi/\lambda$ impinging at
normal incidence, along the $z$ axis.
The electric field (source field) at some position $\mathbf{R}$ emitted by a
collection of emitters, each located at $\mathbf{R}_{j}$, takes the following
expression in the far-field limit Plankensteiner _et al._ (2019); Shahmoon _et
al._ (2017) (only positive frequency component and in the direction of the
emitter dipole moment)

$\displaystyle\mathcal{E}_{\text{dip}}(\mathbf{R},t)\\!=\\!-\frac{3\gamma}{2d_{\text{eg}}}\sum_{j}\beta_{j}(t)\left(\frac{e^{ik|\mathbf{R}-\mathbf{R}_{j}|}}{k|\mathbf{R}-\mathbf{R}_{j}|}\right)e^{-i\omega_{\ell}t}.$
(48)

The total amplitude is given by the sum of dipole-radiated fields, which, in
the linear regime, are proportional to the expectation value of the individual
particle dipole operators $\beta_{j}$ (in a frame rotating at the laser
frequency). We aim at describing regimes where plane waves propagating through
the 2D structure either get reflected or transmitted only in the $z$ direction
while scattering in other directions is inhibited. To this end, we assume a
constant phase illumination (denoting this regime as a symmetric driving case)
where the electric field amplitude does not depend on $x$ and $y$. The total
field can then be written as the sum of the incident field and the field
radiated by the emitters

$\displaystyle\mathcal{E}(z)=\mathcal{E}_{\text{in}}e^{ikz}+\mathcal{E}_{\text{dip}}(z).$
(49)

Let us first look at the case of just a single emitter where the radiated
field is a dipole pattern with the $z$-direction amplitude

$\displaystyle\mathcal{E}(z)=\mathcal{E}_{\text{in}}\left[e^{ikz}+\frac{\alpha_{\text{p}}\pi}{\epsilon_{0}\lambda^{2}}\frac{e^{ik|z|}}{|z|}\right],$
(50)

falling off with increasing distance $z$. We have introduced the single-atom
polarizability $\alpha_{\text{p}}=-d_{\text{eg}}\beta/\mathcal{E}_{\text{in}}$
which can be expressed in terms of the resonant wavelength, linewidth and
$\omega_{\ell}$ as
$\alpha_{\text{p}}=-(3\epsilon_{0}\lambda_{0}^{3})/(4\pi^{2})\gamma/[(\omega_{\ell}-\omega_{0})+i\gamma]$.
The extinction cross section of a quantum emitter $\sigma_{\text{ext}}$ (i.e.,
the effective area seen by an impinging photon) can be related to the
polarizability via
$\sigma_{\text{ext}}(\omega_{\ell})=k\Im(\alpha_{\text{p}})/\epsilon_{0}$
which in the case of resonant illumination simply becomes
$\sigma_{\text{ext}}(\omega_{\ell}=\omega_{0})=3\lambda_{0}^{2}/(2\pi)$.
Notice the interesting aspect that the cross section is much larger
(proportional to $\lambda_{0}^{2}$) than the square of the actual size of the
electronic orbital which can be a milion times smaller.
For many emitters, the driving for the symmetric illumination case can be
included as
$\mathcal{H}_{\ell}=\sum_{j=1}^{\mathcal{N}}\eta(\sigma_{j}^{\dagger}e^{-i\omega_{\ell}t}+\sigma_{j}e^{i\omega_{\ell}t})$
with the Rabi frequency $\eta=d_{\text{eg}}\mathcal{E}_{\text{in}}$.
Transforming into a frame rotating at the laser frequency $\omega_{\ell}$
yields the equations of motion for the dipole amplitudes

$\displaystyle\dot{\beta}_{j}=i(\omega_{\ell}-\omega_{0})\beta_{j}-\sum_{j^{\prime}=1}^{\mathcal{N}}\mathcal{M}_{jj^{\prime}}\beta_{j^{\prime}}-i\eta.$
(51)

Here, the collective effects coming from the coherent and incoherent emitter-
emitter interactions are contained in the off-diagonal elements of the matrix
$\mathcal{M}_{jj^{\prime}}=(i\Omega_{jj^{\prime}}+\gamma_{jj^{\prime}})$
(while the diagonal elements $\mathcal{M}_{jj}=\gamma_{jj}$ give the
independent decay dynamics). The assumption of constant drive $\eta$ over the
whole array combined with the assumption that the array is quasi-infinite
leads to a trivial solution where $\beta_{j}=\beta$ for any
$j=1,...\mathcal{N}$ can be easily estimated from the equation above as
$\beta=\eta/[(\omega_{\ell}-\omega_{0})+i\sum_{j^{\prime}=1}^{\mathcal{N}}\mathcal{M}_{1j^{\prime}}]$.
This in practice means that only a symmetric superposition of all emitters is
excited (only the symmetric surface mode is activated). Notice that the sum
$\sum_{j^{\prime}=1}^{\mathcal{N}}\mathcal{M}_{jj^{\prime}}$ does not depend
on the index $j$ which is why we set it to $1$ above and the result can be
cast into the form
$i\sum_{j^{\prime}=1}^{\mathcal{N}}\mathcal{M}_{1j^{\prime}}=-\Omega_{\text{eff}}+i\gamma_{\text{eff}}$.
The term $\Omega_{\text{eff}}$ leads to a shift of the collective resonance of
the array (as also observed in 1D configurations Glicenstein _et al._ (2020)),
while $\gamma_{\text{eff}}$ describes the effective decay rate of the array.
Also notice that, as we will also detail more in the end of the subsection,
the solution $\beta_{j}=\beta$ holds solely under constant illumination
conditions. Generally, for a laser drive of the form $\mathcal{E}(z)f(x,y)$,
i.e., showing phase or spatial imprinting in the transverse direction,
different surface modes of the array would be driven (depending on the overlap
integral between the function $f(x,y)$ and the transverse profile of the
surface mode.
The sum over the spatial distribution along the $z$-direction can be estimated
by a plane-wave expansion as (for derivation see Appendix A3)

$\displaystyle\sum_{j}\frac{e^{ik|z\mathbf{e}_{z}-\mathbf{R}_{j}|}}{|z\mathbf{e}_{z}-\mathbf{R}_{j}|}=\frac{2\pi
i}{a^{2}}\sum_{m,n=-\infty}^{\infty}\frac{e^{ik_{mn}|z|}}{k_{mn}},$ (52)

where $k_{mn}=\sqrt{k^{2}-q_{m}^{2}-q_{n}^{2}}$ with $q_{m}=(2\pi m/a)$
($q_{n}$ analogously). Here, $q_{m}$ and $q_{n}$ represent the quasimomenta of
the surface modes on the array propagating along the $x$\- and $y$-direction
(Appendix A3). The wave number $k_{mn}$ is real if $(a/\lambda)^{2}\geq
m^{2}+n^{2}$. For a subwavelength lattice, this can then obviously only be
fulfilled with the choice $m,n=0$: this means that only a mode with a $k$
vector equal in amplitude to the impinging laser can propagate in the
$z$-direction. This is exactly the symmetric mode which is also the only one
that constant illumination can activate. Notice that, even if the illumination
phase would slightly vary over the array, any surface mode which is
accidentally excited would radiate in directions other than $z$. One can
finally express the source field radiated by the dipoles in the far-field
limit as

$\displaystyle\mathcal{E}_{\text{dip}}(z)=i\pi\left(\frac{\lambda}{a}\right)^{2}\frac{\alpha_{\text{p}}^{\text{eff}}}{\epsilon_{0}\lambda_{0}^{3}}\mathcal{E}_{\text{in}}e^{ik|z|},$
(53)

with a renormalized effective polarizability summed over the whole ensemble
response as
$\alpha_{\text{p}}^{\text{eff}}=-d_{\text{eg}}\sum_{j}\beta_{j}/(\mathcal{N}\mathcal{E}_{\text{in}})$
and expressed as

$\displaystyle\alpha_{\text{p}}^{\text{eff}}=-\frac{3\epsilon_{0}\lambda_{0}^{3}}{4\pi^{2}}\frac{\gamma}{(\omega_{\ell}-\omega_{0}-\Omega_{\text{eff}})+i\gamma_{\text{eff}}}.$
(54)

We illustrate in Fig. 5c the behavior of the collective rates
$\Omega_{\text{eff}},\gamma_{\text{eff}}$ as a function of the lattice
constant $a/\lambda_{0}$. We note that similarly to the two-particle
interactions discussed in Sec. A, the real and imaginary parts are not
independent but can be connected by a Kramers-Kronig relation. Regions with
$\gamma_{\text{eff}}<\gamma$ (negative values in Fig. 8c) correspond to
collective subradiant behavior and are of particular interest. One can see
that a special operation point occurs e.g. at $a/\lambda_{0}\approx 0.8$ where
the collective frequency shift vanishes while a pronounced subradiant behavior
remains, facilitating the experimental realization of subradiant optical
mirrors with cold atoms in an optical lattice Rui _et al._ (2020).
Notice that the effective cross section per emitter can be considerably
increased as, on resonance, one has
$\sigma_{\text{ext}}^{\text{eff}}(\omega_{\ell})=\sigma_{\text{ext}}(\omega_{\ell})\gamma/\gamma_{\text{eff}}$.
For $a<\lambda_{0}$ one can furthermore approximate the effective decay rate
as $\gamma_{\text{eff}}=3\gamma(\lambda_{0}/a)^{2}/(4\pi)$ [see Fig. 8c]
Shahmoon _et al._ (2017), showing a decrease by a factor of roughly $2.68$
around the optimal operation point at $a/\lambda_{0}\approx 0.8$. This can
also be connected to an increase in the overall reflectivity of the array. To
derive this, we write
$\mathcal{E}_{\text{dip}}(z)=r(\omega_{\ell})\mathcal{E}_{\text{in}}e^{ik|z|}$
where the complex reflectivity amplitude reads (considering small detunings
around the resonance)

$\displaystyle
r(\omega_{\ell})=-i\frac{\gamma_{\text{eff}}}{(\omega_{\ell}-\omega_{0}-\Omega_{\text{eff}})+i\gamma_{\text{eff}}},$
(55)

while the transmission amplitude is obtained as $t=1+r$. In Fig. 8d we plot
the absolute square of these quantities as a function of the lattice constant
for resonant illumination $\omega_{\ell}=\omega_{0}$. One can see that, while
generally the reflectivity of the array is high over a broad range of
separations, for certain values of $a/\lambda_{0}\approx(0.2,0.8)$ the atomic
dipoles can even act as a perfect mirror and reflect the entire input field
with unit efficiency Shahmoon _et al._ (2017). Moreover, for the whole region
where the approximation
$\gamma_{\text{eff}}=3\gamma(\lambda_{0}/a)^{2}/(4\pi)$ holds, the mirror
shows no losses (quantified by the scattering outside the $z$ axis mode),
i.e., $|r|^{2}+|t|^{2}=1$. This is of course only valid in the absence of any
other channels of nonradiative decay (at rate $\gamma_{\text{nr}}$), in which
case the denominator of $r$ acquires an extra term
$i(\gamma_{\text{eff}}+\gamma_{\text{nr}})$. The resulting intensity
distribution of the total electric field is shown in Fig. 8(e), revealing that
the emitters can indeed shut off the transmission of an incident plane wave
for $z>0$.
Let us stress that the expressions above are only valid under symmetric
illumination conditions, which is the relevant experimental situation as
tackled in Ref. Rui _et al._ (2020). As optical lattices have interatomic
distances at the level of $\lambda_{0}/2$ or larger, the simplest operation is
around a point where $\Omega_{\text{eff}}$ vanishes (at $a=0.8\lambda_{0}$)
and unit reflectivity is reached as soon as the laser is resonant to the
emitter transition $\omega_{\ell}=\omega_{0}$. However, the subradiance of
such a symmetrically excited collective state is only a factor of around
$4\pi\times 0.8^{2}/3\approx 2.68$ (at $a=0.8\lambda_{0}$) smaller than
$\gamma$. To fully exploit the scalability of subradiance with $\mathcal{N}$,
one could instead assume antisymmetric, or phased illumination conditions. For
example, Fig. 8f shows that extremely narrow resonances scaling as
$\mathcal{N}^{-5.5}$ can be reached for very dense arrays at
$a=0.1\lambda_{0}$ with antisymmetric phases in a checkerboard pattern where
$f_{nm}=(-1)^{n+m}$ (for $n=1,...,\sqrt{\mathcal{N}}$ and
$m=1,...,\sqrt{\mathcal{N}}$) (deviations from the scaling are due to
imperfect addressing of subradiant states). Let us only sketch how the
derivation above will change in such a case. First, Eq. (49) is changed to
include the $x$ and $y$ spatial dependence of the incoming field in some
function $f(x,y)$ with the periodicity equal to the lattice constant $a$.
Then, the combination of surface modes which are excited by such an
illumination is computed from the steady state of Eqs. (51) where each
$\beta_{nm}$ is driven by an amplitude $\eta f_{nm}$. Finally, replacing the
new solution for $\beta_{nm}$ into Eq. (48), the sum in Eq. (52) needs to be
recomputed. Finally, the specific modes which can propagate into the far field
in the $z$ direction will depend on the specific phased illumination pattern
chosen.

###  D Further remarks

One- and two-dimensional ensembles of coupled quantum emitters arranged in
regular patterns provide an ideal platform for achieving strong light-matter
interactions and high fidelities for photon storage capabilities Facchinetti
_et al._ (2016); Asenjo-Garcia _et al._ (2017). Their subradiance properties
are an important resource for applications ranging from quantum-information
processing Chang _et al._ (2018) and metrology Ostermann _et al._ (2013);
Facchinetti _et al._ (2016); Manzoni _et al._ (2018) to excitation transfer
Moreno-Cardoner _et al._ (2019); Needham _et al._ (2019); Ballantine and
Ruostekoski (2020b) and it is envisioned that one can build quantum matter in
a bottom-up approach from nanoscopic lattices of atoms and photons Chang _et
al._ (2018). A single subradiant array has been shown to act as an ideal
quantum memory with efficient storage and retrieval Manzoni _et al._ (2018)
and it has been suggested Rui _et al._ (2020) that this geometry could lead to
vast improvements in the error bound of quantum memories. Regular or honeycomb
lattices of three-level systems in a $V$-configuration, where time-reversal
symmetry is broken by the application of a magnetic field, have been proposed
as platforms for studying topological phenomena with strongly interacting
photons Perczel _et al._ (2017); Bettles _et al._ (2017) and further
improvements have been proposed in the form of interfacing the arrays with
two-dimensional photonic crystals Perczel _et al._ (2020). A crucial aspect of
these proposals that distinguish them from linear topological photonic systems
is the intrinsic nonlinearity of the quantum emitters which could lead to a
rich many-body physics dynamics on such subradiant lattices. Furthermore,
composite quantum systems comprised of many atomic arrays could find
applications in quantum networking: at the level of two distant layers,
nonlocal entangled Bell superposition states have been shown to exist Guimond
_et al._ (2019). While quantum information processing at the level of atomic
layers require qubit encoding on delocalized spin states over the whole array,
quantum spin lenses have been recently introduced, where incoming flying qubit
photons can be mapped and stored in single atoms Glaetzle _et al._ (2017).

##  IV Cooperativity in cavity QED

Light-matter interactions can be greatly enhanced by using optical elements
which confine electromagnetic fields in very small volumes. This is obvious
from the scaling of the zero-point electric field amplitude
$\mathcal{E}_{k}=\sqrt{\omega_{k}/(2\epsilon_{0}\mathcal{V})}$ which indicates
that stronger field amplitudes per photon mode are achieved for smaller mode
volumes. This led to the development of cavity QED as a subfield of quantum
optics specializing in the description of coherent light-matter interactions
in cavities Haroche and Kleppner (1989); Berman (1994); Walther _et al._
(2006); Haroche and Raimond (2013). This section provides fundamental concepts
and tools of cavity QED with coupled quantum emitter systems, which are then
utilized to describe applications in Sec. V.
We start by introducing a master equation approach to intracavity light-matter
interactions encompassing both collective emitter loss and cavity photon loss.
We then move on to a quantum Langevin equations approach supplemented with an
input-ouput formalism allowing access to the correlations of cavity output
field operators. Finally, to address inhomogeneous emitter ensembles as often
present in experimental setups, we tackle the question of strong frequency
disorder effects on light-matter interactions.

###  A Cavity QED with coupled quantum emitters

The simplest example of an optical cavity is the co-planar design known as a
Fabry-Pérot cavity comprised of two highly reflective parallel mirrors. For a
distance $\ell$ between the mirrors and assuming at first perfect
reflectivity, such a setup defines resonances conditioned by
$\ell=n\lambda/2$. The fundamental mode has, therefore, a wavelength of
$2\ell$ with frequency $\pi c/\ell$ and higher harmonics are multiples of this
mode. The quantity $\omega_{\text{FSR}}=\pi c/\ell$ is known as the free
spectral range and gives the frequency difference between consecutive optical
resonances of the structure. Assuming two mirrors positioned at $z=0$ and
$z=\ell$ and a transverse area $\mathcal{S}$ of the supported optical
resonances, one can proceed in quantizing the field inside the optical
resonator by associating bosonic operators $a_{n}$ to all resonances and
writing the total Hamiltonian as
$\mathcal{H}_{c}=\sum_{n}\omega_{n}a_{n}^{\dagger}a_{n}$. The electric field
operator can then be decomposed as

$\hat{\mathbf{E}}(z)=\textstyle\sum_{n}\mathcal{E}_{n}\left(a_{n}+a_{n}^{\dagger}\right)\sin{(k_{n}z)}{\bm{\epsilon}}_{n},$
(56)

where the sum runs over all allowed wave vectors $k_{n}=n\pi/\ell$ and the
zero-point electric field amplitude
$\mathcal{E}_{n}=\sqrt{\omega_{n}/(\epsilon_{0}\ell\mathcal{S})}$ and all
polarizations ${\bm{\epsilon}}_{n}$ are assumed orthogonal to the
$z$-direction. For all situations we will consider in the following, a single
cavity mode suffices. We denote the mode by the operator $a$ at frequency
$\omega_{c}$ such that the free Hamiltonian becomes
$\mathcal{H}_{c}=\omega_{c}a^{\dagger}a$.

ME for a driven, lossy cavity mode \- Perfect mirrors define infinitely sharp
resonances; in order to allow the in-coupling of light into the optical
resonator side mirrors with slightly less than unit reflectivity are used. For
a double-sided optical resonator, the localized resonances can then couple to
the infinite number of modes to the left and right of the cavity, leading to a
loss of the intracavity photons. This loss can be described within the master
equation formalism by following a phenomenological model where photons from
mode $a$ can tunnel to the left $b(\omega)$ and right $c(\omega)$ continuum of
free radiation modes via an excitation exchange Hamiltonian. The standard
open-system dynamics approach detailed in Appendix A1 leads to a Lindblad form
of the cavity photon decay

$\mathcal{L}_{\kappa}[\rho]=\kappa\left[2a\rho
a^{\dagger}-a^{\dagger}a\rho-\rho a^{\dagger}a\right],$ (57)

where $\kappa=\kappa_{R}+\kappa_{L}$ encompasses total losses via both the
right-side and left-side mirrors.
To allow for the driving of the cavity mode let us now assume a continuous-
wave laser with power $\mathcal{P}$ entering the cavity from the left side.
This can be included in the following Hamiltonian

$\mathcal{H}_{\ell}=i\eta(a^{\dagger}e^{-i\omega_{\ell}t}-ae^{i\omega_{\ell}t}),$
(58)

where the drive amplitude is
$\eta=\sqrt{2\kappa_{L}\mathcal{P}/\omega_{\ell}}$. The equation of motion for
the expectation value of the cavity operator can be derived from the master
equation (with Hamiltonian $\mathcal{H}_{c}+\mathcal{H}_{\ell}$ and Lindblad
term $\mathcal{L}_{\kappa}[\rho]$) leading to

$\braket{\dot{a}}=-\kappa\braket{a}-i(\omega_{c}-\omega_{\ell})\braket{a}+\eta.$
(59)

The steady state of the equation above (reached by requiring that
$\braket{\dot{a}}=0$) describes the expected Lorentzian response
$\braket{a}_{\text{ss}}=\eta/(\kappa+i(\omega_{c}-\omega_{\ell}))$ of the
cavity field exhibiting a linewidth $\kappa$ and resonance frequency
$\omega_{c}$.

Tavis-Cummings Hamiltonian \- Let us consider now $\mathcal{N}$ nonidentical
quantum emitters with transition frequencies $\omega_{j}$ placed inside the
optical cavity around a single optical mode of interest. Assuming the field to
be varying only along the cavity axis, for emitter $j$ positioned at $z_{j}$
within the cavity volume, the dipolar interaction is characterized by a
position-dependent coupling strength
$g_{j}=\mathcal{E}d_{\text{eg}}\sin(kz_{j})$. The light-matter Hamiltonian
describing excitation exchange between the emitters and the cavity mode can
then be cast in the standard Tavis-Cummings Tavis and Cummings (1968) form

$\displaystyle\mathcal{H}_{\text{TC}}=\sum_{j=1}^{\mathcal{N}}\omega_{j}\sigma^{\dagger}_{j}\sigma_{j}+\omega_{c}a^{\dagger}a+\sum_{j=1}^{\mathcal{N}}\left(g_{j}a\sigma^{\dagger}_{j}+g_{j}^{*}\sigma_{j}a^{\dagger}\right).$
(60)

This interaction is a particular form of the free-space case in Eq. (3) where
the coupling to the cavity mode $a$ is enhanced while the coupling to all
other modes is inhibited. At the single emitter level, the interaction is
known as the Jaynes-Cummings Hamiltonian Jaynes and Cummings (1963). Notice
that such Hamiltonian conserves the excitation number and can be exactly
solved Tavis and Cummings (1968). The excitation nonconserving case (known
also as counter-rotating terms) for $\mathcal{N}=1$ is the exactly integrable
Braak (2011), Rabi Hamiltonian Rabi (1936). The extension of the Rabi model to
$\mathcal{N}>1$ for identical emitters $\omega_{j}=\omega_{0}$ and identical
couplings $g_{j}=g$ is a simplified case describing exchange between a single
collective operator $\sum_{j}\sigma_{j}$ and the cavity mode. This is known as
the Dicke model, generally nonintegrable Emary and Brandes (2003); Braak
(2013) except for when $\omega_{0}=0$, and in the thermodynamic limit with
$\mathcal{N}\rightarrow\infty$ Hepp and Lieb (1973); Brankov _et al._ (1975);
Bogolubov and Plechko (1976) where phase transitions have been identified Wang
and Hioe (1973).

Linear optical response \- For dense, interacting ensembles, the dipole-dipole
interactions $\mathcal{H}_{dd}$ defined in Eq. (10) and the collective decay
terms as defined in Eq. (7) play an important role. To follow the system
evolution, we then write the master equation with the total Hamiltonian
$\mathcal{H}_{\text{TC}}+\mathcal{H}_{dd}+\mathcal{H}_{\ell}$ and Lindblad
term $\mathcal{L}_{\kappa}[\rho]+\mathcal{L}_{e}[\rho]$. From here, one can
derive equations of motion for single operator amplitudes. The emerging set of
coupled equations is not closed, i.e., expectation values of single operators
are coupled to two-operator correlations which again couple to three or more
operator correlations and so on. However, in a first approximation, assuming
weak driving conditions ($\eta\ll\gamma$) some factorization of operators
under the approximation $\braket{\sigma_{z}^{j}}\approx-1$ can occur. This
approximation is sufficient to describe the linear optical response as used
throughout Sec. A. In Sec. C instead, we will go beyond this approximation to
describe population-inverted systems exhibiting lasing behavior.
Let us denote the expectation values by $\alpha=\braket{a}$ and
$\beta_{j}=\braket{\sigma_{j}}$ and derive

$\displaystyle\dot{\alpha}$
$\displaystyle=-(\kappa+i\omega_{c})\alpha-i\sum_{j=1}^{\mathcal{N}}g_{j}\beta_{j}+\eta
e^{-i\omega_{\ell}t}$ (61a) $\displaystyle\dot{\beta}_{j}$
$\displaystyle=-i\omega_{j}\beta_{j}-ig_{j}\alpha-\sum_{j^{\prime}=1}^{\mathcal{N}}\mathcal{M}_{jj^{\prime}}\beta_{j^{\prime}}.$
(61b)

In Sec. B these equations will be extended to the quantum regime within the
quantum Langevin equations formalism.

Regimes of cavity QED \- Let us quickly review a few standard concepts of
cavity QED with non-interacting ensembles by setting $\Omega_{jj^{\prime}}=0$
and $\gamma_{jj^{\prime}}=\gamma\delta_{jj^{\prime}}$. In steady state, the
normalized transmission of the cavity $t=\kappa\alpha/\eta$ can be easily
computed from Eqs. (61) to lead to

$\displaystyle
t=\kappa\left[\kappa+i(\omega_{c}-\omega_{\ell})+\sum_{j=1}^{\mathcal{N}}\frac{|g_{j}|^{2}}{\gamma+i(\omega_{j}-\omega_{\ell})}\right]^{-1}.$
(62)

The cavity transmission can then show different regimes either characterized
by a single Lorentzian peak (empty cavity or weak light-emitter coupling), a
peak with a narrow dip (the Purcell, or antiresonant regime for intermediate
couplings) or a polaritonic regime, where two cavity resonances are resolved
in transmission. These regimes are dependent on the magnitude of the loss
rates relative to the coupling strength. Mathematically, they can be easily
deduced from the eigenproblem defined by the evolution matrix of Eqs. (61).
For $\mathcal{N}=1$, diagonalization of Eqs. (61) leads to normal mode
splitting when $g>|\kappa-\gamma|/2$ with new frequencies
$\omega_{\pm}=\omega_{0}\pm\sqrt{g^{2}+(\kappa-\gamma)^{2}/4}$ (on resonance
$\omega_{0}=\omega_{c}$). The difference $\omega_{+}-\omega_{-}$ is known as
the vacuum Rabi splitting (VRS). The newly formed quantum states in such a
strong coupling regime are hybrid light-matter ones, i.e., polaritons
$(\ket{g1}\pm\ket{e0})/\sqrt{2}$ and can be read in the cavity transmission
$t$ as two well-separated peaks. For $\mathcal{N}>1$, the collective strong
coupling regime can be defined when $g_{\mathcal{N}}>|\kappa-\gamma|/2$ with a
collective coupling rate $g_{\mathcal{N}}=\sqrt{\sum_{j}{|g_{j}|^{2}}}$. The
loss rates of the polaritons are equal in this regime $(\kappa+\gamma)/2$. One
can also identify a weak coupling regime for which
$\gamma<g_{\mathcal{N}}<\kappa$ but with a strong cooperativity
$\mathcal{C}_{\mathcal{N}}=g_{\mathcal{N}}^{2}/(\kappa\gamma)\gg 1$. For many
emitters equally coupled to the cavity mode, the cooperativity shows a
collective enhancement proportional to $\mathcal{N}$. For a single emitter,
this regime (known as the Purcell regime) leads to a renormalization of the
intrinsic spontaneous emission rate to $\gamma(1+\mathcal{C})$ understood as
the addition of a new decay channel via loss through the cavity. A scan of $t$
as a function of laser frequency shows in such a case a dip in the cavity
transmission around $\omega_{\ell}=\omega_{0}$ known as a cavity
antiresonance. The case of $\mathcal{N}$ coupled emitters will be presented in
Sec. A showing that the addressing of collective subradiant resonances can
lead to a superlinear $\propto\mathcal{N}^{4}$ scaling of
$\mathcal{C}_{\mathcal{N}}$.

###  B Input-output formalism for operators

The master equation formalism previously utilized is based on discarding the
state of the environment and focusing on the dynamics of the reduced density
operator of the much smaller system of interest. An alternative approach are
quantum Langevin equations (also called Heisenberg-Langevin equations)
Gardiner and Zoller (2004); Jacobs (2010) which follow the evolution at the
level of system operators instead of system states, i.e., its density
operator. In essence, the formalism consists in supplementing the Heisenberg
equations of motion with the proper dissipation and fluctuation terms.

Mapping the ME onto QLEs \- For open-system dynamics described by loss in
Lindblad form, a direct transformation between the master equation and the QLE
formulation exists that reads Gardiner and Zoller (2004)

$\displaystyle\dot{\mathcal{O}}=\frac{i}{\hbar}[\mathcal{H},\mathcal{O}]$
$\displaystyle-[\mathcal{O},c^{\dagger}]\left\\{\gamma_{\text{c}}c+\sqrt{2\gamma_{\text{c}}}c_{\text{in}}\right\\}$
(63)
$\displaystyle+\left\\{\gamma_{\text{c}}c^{\dagger}+\sqrt{2\gamma_{\text{c}}}c_{\text{in}}^{\dagger}\right\\}[\mathcal{O},c].$

The equation above for a generic system operator $\mathcal{O}$ is applied to
each individual Lindblad collapse operator $c$ acting at rate
$\gamma_{\text{c}}$ and with associated input noise $c_{\text{in}}$. The input
noise is a stochastic operator which is zero averaged
$\braket{c_{\text{in}}}=0$, delta-correlated in time
$\braket{c_{\text{in}}(t)c_{\text{in}}^{\dagger}(t^{\prime})}=\delta(t-t^{\prime})$
(all other correlations vanish) and satisfies the following commutation
relation
$[c_{\text{in}}(t),c_{\text{in}}^{\dagger}(t^{\prime})]=\delta(t-t^{\prime})$.
For linear evolution, such an equation can be formally integrated and
expectation values of correlations containing any number of operators can be
obtained.
We will proceed by first exemplifying the derivation and solution to a QLE
applied to a lossy bosonic field corresponding to a decaying and externally
driven optical cavity mode. Then we introduce the input-output relations
allowing the readout of intracavity light-matter interactions via the cavity
output field. Finally, we introduce a set of coupled linear QLEs for an
ensemble of coupled, collectively decaying quantum emitters within the same
optical cavity mode and show how quantum correlations can be computed. We will
then utilize QLEs in Sec. A for the description of the optical response of
subradiant arrays in optical cavities and in Sec. B for the description of
hybrid cavities. An extension of the QLEs to include electron-vibron coupling
will be introduced in Sec. VI to analytically tackle molecule-photon dynamics.

QLEs and input-output relations. \- For a driven, lossy, single-sided optical
cavity, Eq. (63) yields the following QLE

$\dot{a}=-(\kappa+i\omega_{c})a+\eta
e^{-i\omega_{\ell}t}+\sqrt{2\kappa}a_{\text{in}}.$ (64)

Notice that an expectation value of the equation above simply reproduces Eq.
(59) which has been alternatively derived directly from the master equation.
The more general Langevin equation also allows for an exact solution
obtainable by formal integration. We separate this solution into a transient
and steady state part $a(t)=a_{\text{tr}}(t)+a_{\text{ss}}(t)$. While the
transient part $a_{\text{tr}}(t)$ decays away on a timescale defined by the
cavity loss $\kappa^{-1}$, the steady state part $a_{\text{ss}}(t)$ dominates
in the long time (steady state) limit

$\displaystyle a_{\text{tr}}(t)$ $\displaystyle=a(0)e^{-\kappa
t}e^{-i\omega_{c}t},$ (65a) $\displaystyle a_{\text{ss}}(t)$
$\displaystyle=\int_{0}^{t}dt^{\prime}(\eta
e^{-i\omega_{\ell}t^{\prime}}+\sqrt{2\kappa}a_{\text{in}}(t^{\prime}))e^{-(\kappa+i\omega_{c})(t-t^{\prime})}.$
(65b)

It is important to notice that the inclusion of the noise terms in such
equations is crucial: while the commutator
$[a_{\text{tr}}(t),a_{\text{tr}}^{\dagger}(t)]=e^{-2\kappa t}$ decays away,
the commutator $[a_{\text{ss}}(t),a_{\text{ss}}^{\dagger}(t)]=1-e^{-2\kappa
t}$ maintains the standard bosonic commutation relation
$[a(t),a^{\dagger}(t)]=1$ at every time during evolution.
In addition to the equation of motion for the cavity field amplitude
operators, the full description of the optical cavity includes the input-
output relation

$\displaystyle a_{\text{out}}+a_{\text{in}}=\sqrt{\kappa}a(t),$ (66)

which allows for the derivation of the properties of the cavity output mode.
For example, the average shows that the cavity transmission can be simply
written as: $t=\braket{a_{\text{out}}}/\eta=\sqrt{\kappa}\braket{a}/\eta$ (as
already used and consistent with the result of Eq. (59) in Sec. A).

QLEs for cooperative cavity QED \- We now assume a cavity mode coupled to
$\mathcal{N}$ interacting emitters, as a system described by the Hamiltonian
$\mathcal{H}_{\text{TC}}+\mathcal{H}_{dd}+\mathcal{H}_{\ell}$ and loss rates
incorporated in the terms $\mathcal{L}_{\kappa}[\rho]+\mathcal{L}_{e}[\rho]$
(as introduced in Sec. A). In a first step, we notice that the decay terms are
not in Lindblad form, thus we perform a diagonalizing transformation to bring
it into the form of $\mathcal{N}$ independent decay channels as described by
Eq. (12) (such that Eq. (63) can then be directly applied). We also only
consider the low excitation limit $\braket{\sigma_{z}^{j}}\approx-1$. Notice
that this approximation discards the nonlinear response of the dipoles, which
could, for example, lead to a collective Kerr nonlinear effect Plankensteiner
_et al._ (2019). One then obtains a set of coupled linear differential
equations

$\displaystyle\dot{a}$
$\displaystyle=-\left(\kappa-i\Delta_{c}\right)a-i\sum_{j}g_{j}\sigma_{j}+\sqrt{2\kappa}a_{\text{in}}+\eta,$
(67a) $\displaystyle\dot{\sigma}_{j}$
$\displaystyle=i\Delta\sigma_{j}-ig_{j}a-\sum_{j^{\prime}}\mathcal{M}_{jj^{\prime}}\sigma_{j^{\prime}}+\sqrt{2\gamma}\sigma_{\text{in}}^{j},$
(67b)

with the emitter and cavity detunings defined as
$\Delta=\omega_{\ell}-\omega_{0}$ and $\Delta_{c}=\omega_{\ell}-\omega_{c}$
respectively. The input noise terms $\sigma_{\text{in}}^{j}$ (stemming from
the coupling of the emitters to the electromagnetic modes outside the solid-
state angle covered by the cavity field) are considered here as zero-averaged
and delta-correlated in time
$\braket{\sigma_{\text{in}}^{j}(t)\sigma_{\text{in}}^{j^{\prime},\dagger}(t^{\prime})}=(\gamma_{jj^{\prime}}/\gamma)\delta(t-t^{\prime})$
as a consequence of the low excitation assumption (for the more complex case
of high excitation we refer to Refs. Gardiner and Zoller (2004);
Plankensteiner _et al._ (2019)). We can write the equations above in a compact
form in terms of the matrix
$\textbf{M}(\Delta)=-i\Delta\mathbbm{1}+i\bm{\Omega}+\bm{\Gamma}$ and the
coupling vector $\textbf{G}=(g_{1},g_{2},...,g_{\mathcal{N}})^{\top}$.
Additionally defining the vectors
$\bm{\sigma}=(\sigma_{1},\ldots,\sigma_{\mathcal{N}})^{\top}$ and
$\bm{\sigma}_{\text{in}}=(\sigma_{\text{in}}^{1},\ldots,\sigma_{\text{in}}^{\mathcal{N}})^{\top}$,
Eqs. (67) express as

$\displaystyle\dot{a}$
$\displaystyle=-(\kappa-i\Delta_{c})a-i\mathbf{G}^{\top}\bm{\sigma}+\sqrt{2\kappa}a_{\text{in}}+\eta,$
(68a) $\displaystyle\dot{\bm{\sigma}}$
$\displaystyle=-\mathbf{M}(\Delta)\bm{\sigma}-i\mathbf{G}a+\sqrt{2\gamma}\bm{\sigma}_{\text{in}}.$
(68b)

This is the starting point for the derivations in Sec. A where we compute the
classical and quantum response of a weakly excited transverse 1D or 2D array
of coupled quantum emitters to a driven cavity mode.

Quantum correlations \- From the equations of motion for operators Eqs. (68),
one can go a step further and fully analyze the quantum properties (second-
order correlation function for light $g^{(2)}$, bipartite entanglement,
squeezing properties, etc.) of both the photon and matter counterparts. To
this end, let us first consider only fluctuation operators only by expansion
around expectation values $a=\alpha+\delta a$,
$\sigma_{j}=\beta_{j}+\delta\sigma_{j}$. We can then cast Eqs. (68) in
convenient vector form as
$\dot{\mathbf{v}}=\mathbf{A}\mathbf{v}+\mathbf{N}\mathbf{v}_{\text{in}}$ for
the vectors $\mathbf{v}=(\delta a,\delta
a^{\dagger},\delta\bm{\sigma},\delta\bm{\sigma}^{\dagger})^{\top}$,
$\mathbf{v}_{\text{in}}=(a_{\text{in}},a_{\text{in}}^{\dagger},\bm{\sigma}_{\text{in}},\bm{\sigma}_{\text{in}}^{\dagger})^{\top}$
with the drift matrix expressed in compact form as (for details also see Ref.
Plankensteiner _et al._ (2019))

$\displaystyle\mathbf{A}$
$\displaystyle=\begin{pmatrix}-\left(\kappa-i\Delta_{c}\right)&0&-i\textbf{G}^{\top}&\textbf{0}^{\top}\\\
0&-(\kappa+i\Delta_{c})&\textbf{0}^{\top}&i\textbf{G}^{\top}\\\
-i\textbf{G}&\textbf{0}&-\textbf{M}(\Delta)&\underline{\textbf{0}}\\\
\textbf{0}&i\textbf{G}&\underline{\textbf{0}}&-\textbf{M}^{*}(\Delta)\end{pmatrix}.$
(69)

We have defined the vector 0 containing $\mathcal{N}$ zeros and
$\underline{\textbf{0}}$ is a $\mathcal{N}\times\mathcal{N}$ matrix with only
zeros. The matrix multiplying the input noise operators is

$\displaystyle\mathbf{N}=\begin{pmatrix}\sqrt{2\kappa}&0&\textbf{0}^{\top}&\textbf{0}^{\top}\\\
0&\sqrt{2\kappa}&\textbf{0}^{\top}&\textbf{0}^{\top}\\\
\textbf{0}&\textbf{0}&\sqrt{2\gamma}\mathbbm{1}&\underline{\textbf{0}}\\\
\textbf{0}&\textbf{0}&\underline{\textbf{0}}&\sqrt{2\gamma}\mathbbm{1}\end{pmatrix}.$
(70)

Formally integrating this system of linearly coupled equations gives the
solution

$\displaystyle\mathbf{v}(t)=e^{\mathbf{A}t}\mathbf{v}(0)+\int_{0}^{t}dt^{\prime}e^{\mathbf{A}(t-t^{\prime})}\mathbf{N}\mathbf{v}_{\text{in}}(t^{\prime}).$
(71)

If all eigenvalues of the drift matrix are negative, the system is stable and
will go towards a steady state (i.e., the transient term $e^{\mathbf{A}t}$
vanishes). From the steady state solution only, one can then define a
correlation matrix $\mathbf{V}=\braket{\mathbf{v}(t)\mathbf{v}^{\top}(t)}$
which can be expressed as

$\displaystyle\mathbf{V}=\int_{0}^{t}dt^{\prime}e^{\mathbf{A}(t-t^{\prime})}\mathbf{D}e^{\mathbf{A}^{\top}(t-t^{\prime})},$
(72)

where we used that
$\braket{\mathbf{v}_{\text{in}}(t^{\prime})\mathbf{v}_{\text{in}}^{\top}(t^{\prime\prime})}=\mathbf{C}\delta(t^{\prime}-t^{\prime\prime})$
and defined the diffusion matrix as
$\mathbf{D}=\mathbf{N}\mathbf{C}\mathbf{N}^{\top}$. The matrix containing the
correlations can be readily obtained as

$\displaystyle\mathbf{C}=\begin{pmatrix}0&1&\textbf{0}^{\top}&\textbf{0}^{\top}\\\
0&0&\textbf{0}^{\top}&\textbf{0}^{\top}\\\
\textbf{0}&\textbf{0}&\underline{\textbf{0}}&\mathbf{\Gamma}/\gamma\\\
\textbf{0}&\textbf{0}&\underline{\textbf{0}}&\underline{\textbf{0}}\end{pmatrix}.$
(73)

Figure 9: Second-order optical correlation function $g^{(2)}(0)$ of a driven
emitter-cavity system as a function of the cavity-emitter detuning
$\omega_{c}-\omega_{0}$. Parameters: $\gamma=0.1\kappa$, $\eta=5\cdot
10^{-3}\kappa$, $g=0.3\kappa$.

For the covariance matrix, one can then derive the Lyapunov equation using
integration by parts

$\displaystyle\mathbf{A}\mathbf{V}$
$\displaystyle+\mathbf{V}\mathbf{A}^{\top}=\int_{0}^{t}dt^{\prime}\mathbf{A}e^{\mathbf{A}(t-t^{\prime})}\mathbf{D}e^{\mathbf{A}^{\top}(t-t^{\prime})}+\mathbf{V}\mathbf{A}^{\top}=$
$\displaystyle=-e^{\mathbf{A}(t-t^{\prime})}\mathbf{D}e^{\mathbf{A}^{\top}(t-t^{\prime})}\Big{|}_{0}^{t}-\mathbf{V}\mathbf{A}^{\top}+\mathbf{V}\mathbf{A}^{\top}=$
$\displaystyle=-\mathbf{D}.$ (74)

Alternatively, one can also perform a Fourier analysis to turn the system of
differential equations to an algebraic set of coupled equations
$i\omega\mathbf{v}(\omega)=\mathbf{A}\mathbf{v}(\omega)+\mathbf{N}\mathbf{v}_{\text{in}}(\omega)$,
which allows to relate the intracavity operators to the input noise via

$\displaystyle\mathbf{v}(\omega)=\left(i\omega\mathbbm{1}-\mathbf{A}\right)^{-1}\mathbf{N}\mathbf{v}_{\text{in}}(\omega).$
(75)

Furthermore, one can relate the input fields directly to the output fields by
using the fact that in the time domain
$\mathbf{v}_{\text{out}}(t)=\mathbf{N}^{\top}\mathbf{v}(t)-\mathbf{v}_{\text{in}}(t)$
and therefore one can express
$\mathbf{v}_{\text{out}}(\omega)=\mathbf{F}(\omega)\mathbf{v}_{\text{in}}(\omega),$
where
$\mathbf{F}(\omega)=\left[\mathbf{N}^{\top}(i\omega\mathbbm{1}-\mathbf{A})^{-1}\mathbf{N}-\mathbbm{1}\right]$.
In the frequency domain, the $\delta$-correlations are preserved for the
output fields
$\braket{\mathbf{v}_{\text{out}}(\omega)\mathbf{v}_{\text{out}}^{\top}(\omega^{\prime})}=\bm{\mathcal{S}}_{\text{out}}(\omega)\delta(\omega+\omega^{\prime})$,
where the two-operator correlations are completely encoded in the spectrum
matrix given by

$\displaystyle\bm{\mathcal{S}}_{\text{out}}(\omega)=\mathbf{F}(\omega)\mathbf{C}\mathbf{F}^{\top}(-\omega).$
(76)

While this suffices to characterize the two-operator correlations of the
system, also higher-order correlations (e.g. four-operator correlations) can
be of interest. We note that all higher order correlations can be expressed as
a sum over products of two-point correlations via the Isserlis’ theorem. An
example of such a correlation is the the $g^{(2)}$-function already considered
in Sec. C, characterizing the photon statistics of the system. It is defined
in steady state as

$\displaystyle
g^{(2)}(\tau)=\frac{\braket{a^{\dagger}(t)a^{\dagger}(t+\tau)a(t+\tau)a(t)}_{\text{ss}}}{\braket{a^{\dagger}(t)a(t)}_{\text{ss}}^{2}}.$
(77)

Nonclassical light sources show sub-Poissonian statistics $g^{(2)}(0)<1$
(antibunching), implying that it becomes unlikely for photons to be detected
in pairs. In Fig. 9 this is illustrated for the coupled cavity-emitter system
where the mechanism originates from the anharmonicity of the Jaynes-Cummings
ladder, allowing only single photons inside the cavity.

###  C Cavity QED with disordered ensembles

In the most general case of Eqs. (61) the frequencies of the emitters are not
identical. This is often encountered in experiments, as the coupling of
electronic systems to embedding matrices can lead to strong inhomogeneous
broadening, such as it is the case for quantum dots or molecular systems. To
quantify frequency disorder, we assume that the frequencies $\omega_{j}$ are
distributed around $\omega_{0}$ according to the distribution function
$p(\delta)$ normalized to unity $\int_{-\infty}^{\infty}p(\delta)d\delta=1$.
In particular we choose a Gaussian distribution of frequencies
$p(\delta)=(1/\sqrt{2\pi w^{2}})e^{-\delta^{2}/(2w^{2})}$. We write the
transition frequencies as $\omega_{j}=\omega_{0}+\delta_{j}$ with vanishing
classical average $\braket{\delta_{j}}_{\text{cl}}=0$ and variance
$\braket{\delta_{j}^{2}}_{\text{cl}}=w^{2}$. For simplicity, we restrict the
discussion to the case of identical couplings $g_{j}=g$ (for all $j$) while
more generally we refer the reader to Ref. Sommer _et al._ (2021).
For $w=0$ we note that the cavity couples only to a symmetric superposition
$\hat{B}=\textstyle\sum_{j}\sigma_{j}/\sqrt{\mathcal{N}}$, i.e., a bright
state, with a collective coupling strength
$g_{\mathcal{N}}=\sqrt{\mathcal{N}}g$. The other $\mathcal{N}-1$ combinations
define dark states which are obtainable by a Gram-Schmidt algorithm that leads
to all vectors orthogonal to the bright state and to each other. However, a
simple choice of coefficients is indicated by a discrete Fourier transform
$\hat{D}_{k}=1/\sqrt{\mathcal{N}}\textstyle\sum^{\mathcal{\mathcal{N}}}_{j=1}e^{-i2\pi
jk/\mathcal{N}}\sigma_{j}$. We index the dark-state manifold for
$k=1,\dots,\mathcal{N}-1$ and note that for $k=\mathcal{N}$ we have
$\hat{D}_{\mathcal{N}}=\hat{B}$. The equations of motion for all variables
$\mathcal{D}_{k}=\braket{\hat{D}_{k}}$ and $\alpha=\braket{a}$ become (in a
frame rotating at $\omega_{0}$)

$\displaystyle\dot{\mathcal{D}}_{k}$
$\displaystyle=-\gamma\mathcal{D}_{k}-i\sum^{\mathcal{N}}_{k^{\prime}=1}\Delta_{kk^{\prime}}\mathcal{D}_{k^{\prime}}-ig_{\mathcal{N}}\alpha\delta_{k\mathcal{N}},$
(78a) $\displaystyle\dot{\alpha}$ $\displaystyle=-(\kappa-i\Delta_{c})\alpha-
ig^{*}_{\mathcal{N}}\mathcal{D}_{\mathcal{N}}+\eta,$ (78b)

where the couplings between collective states are defined as
$\Delta_{kk^{\prime}}=1/\mathcal{N}\textstyle\sum^{\mathcal{N}}_{j=1}\delta_{j}e^{-i2\pi
j(k-k^{\prime})/\mathcal{N}}$. For $k=\mathcal{N}$, the equations above
indicate that the bright state couples to the cavity mode with the standard
$g_{\mathcal{N}}$ rate, while also being coupled to all dark states.

Figure 10: In the case of disorder of the electronic transition frequencies,
the dark modes are coupled to the bright mode and to each other while the
bright mode couples to the cavity mode with rate $g_{\mathcal{N}}$.
Elimination of the dark, cavity uncoupled modes then leads to an effective two
hybridized mode problem (cavity mode and bright collective state) with
renormalized loss rates.

Non-Markovian regime \- We proceed (see Fig. 10) to eliminate the dark-state
manifold in an exact way without making a Markovian approximation (which would
imply that the dark state reservoir has no memory and therefore it would allow
to set all derivatives of $\mathcal{D}_{k}$ to zero). Instead, we formally
integrate the equations for $\mathcal{D}_{k}$ to obtain

$\dot{\mathcal{B}}(t)=-(\gamma+i\bar{\delta})\mathcal{B}(t)-\int_{-\infty}^{\infty}dt^{\prime}f(t-t^{\prime})\mathcal{B}(t^{\prime})-ig_{N}\alpha,$
(79)

where the memory kernel generally describes a non-Markovian loss process. In
the mesoscopic limit $\mathcal{N}\to\infty$, one finds

$f(t-t^{\prime})\approx\Theta(t-t^{\prime})w^{2}e^{-i(\bar{\delta}-i\gamma)(t-t^{\prime})}\operatorname{sinc}(2w(t-t^{\prime})).$
(80)

This now allows to identify the Markovian regime discussed in the next
paragraph, where the condition that the disorder $w\gg\gamma$ implies that the
kernel $f(t-t^{\prime})$ tends towards a delta function.

Markovian limit \- In the Markovian limit, the elimination of the dark-state
manifold is straightforward as all the derivatives $\dot{\mathcal{D}}_{k}$ for
$k\neq\mathcal{N}$ can be set to zero. This means
$\mathcal{D}_{k}=-\sum_{k^{\prime}}({\mathcal{A}^{-1}})_{kk^{\prime}}\Delta_{k^{\prime}\mathcal{N}}\mathcal{B}$
with the matrix

$\mathcal{A}=\left(\begin{array}[]{cccc}(\bar{\delta}-i\gamma)&\Delta_{12}&\ldots&\Delta_{1(\mathcal{N}-1)}\\\
\Delta_{21}&(\bar{\delta}-i\gamma)&\ldots&\Delta_{2(\mathcal{N}-1)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\Delta_{(\mathcal{N}-1)1}&\Delta_{(\mathcal{N}-1)2}&\ldots&(\bar{\delta}-i\gamma)\\\
\end{array}\right),$ (81)

and
$\bar{\delta}=\Delta_{kk}=\textstyle\sum^{\mathcal{N}}_{j=1}\delta_{j}/\mathcal{N}$.
The Markovian kernel is explicitly
$f(t-t^{\prime})=(i\delta_{\text{dark}}+\gamma_{\text{dark}})\delta(t-t^{\prime})$,
containing an effective frequency shift $\delta_{\text{dark}}$ and a loss rate
$\gamma_{\text{dark}}$ derived from
$\delta_{\text{dark}}+i\gamma_{\text{dark}}=\sum_{k,k^{\prime}=1}^{\mathcal{N}-1}\Delta_{\mathcal{N}k}(\mathcal{A}^{-1})_{kk^{\prime}}\Delta_{k^{\prime}\mathcal{N}}$.
In the mesoscopic limit, one can further simplify the expression of the loss
rate and derive scaling laws where $\gamma_{\text{dark}}=w^{2}/\gamma$ for
$\gamma\gg w$ and $\gamma_{\text{dark}}=\pi w/4$ for $w\gg\gamma$. Moreover,
both $\bar{\delta}$ and $\delta_{\text{dark}}$ vanish in this regime. The
dependence of the VRS on disorder can then be obtained by diagonalizing the
dynamics in the reduced cavity-bright state subspace

$\displaystyle\text{VRS}\approx\Im\left\\{2\sqrt{\left(\gamma+\gamma_{\text{dark}}-\kappa\right)^{2}/4-g_{\mathcal{N}}^{2}}\right\\}.$
(82)

The expression above shows that a large degree of disorder (on the order of
$\kappa$) can lead to a strong reduction of the VRS and consequently pull the
system out of the collective strong coupling regime. In Ref. Sommer _et al._
(2021) a more detailed derivation is provided that shows that particles which
have a large detuning with respect to the cavity frequency, or a very lossy
behavior are effectively pulled out of the macroscopic superposition
participating in the strong coupling condition. This in turn leads to the
degradation of the VRS described by the equation above.

##  V Applications in cavity QED

We now discuss a few applications of cavity QED with coupled quantum-emitter
ensembles. For 1D and 2D arrays of emitters, we now have systems exhibiting
cooperative surface resonances, as already analyzed for free space
applications in Sec. III, but with an extra degree of freedom introduced by a
cavity-confined field mode. In the limit of low reflectivity, such systems
have been shown to act as quick frequency switchers Plankensteiner _et al._
(2017) and to give rise to enhanced optical nonlinearities Plankensteiner _et
al._ (2019). In the high reflectivity limit, such as theoretically predicted
in Refs. Bettles _et al._ (2016); Shahmoon _et al._ (2017) and experimentally
proven in Ref. Rui _et al._ (2020), we introduce the necessary steps, as laid
out in Ref. Cernotik _et al._ (2019), to derive the correct input-output
formalism for optical resonators made up of subwavelength arrays as end-
mirrors.
In the weak reflectivity limit, the array reacts as a strongly dispersive
optical element within a short frequency window but it does not, at the same
time, considerably change the spatial profile of the cavity mode (owing to its
weak reflectivity). The treatment is then perturbative, within the standard
Tavis-Cummings formalism and nontrivial effects occur such as enhanced optical
nonlinearities and the reach of an enhanced collective Purcell effect with a
cooperativity scaling up to $\propto\mathcal{N}^{4}$ (for a 1D arrangement).
In the opposite case, where the emitter array acts as a near-unity
reflectivity end-mirror, a hybrid cavity design emerges where asymmetric
transmission profiles can be achieved, potentially much narrower than those
obtained with frequency-independent mirrors of comparable reflectivity. For
such designs, the standard input-output theory and the master equation for the
photon mode losses (from Sec. A) loses validity, as the tunneling rate of
photons strongly depends on their frequency. We present a general roadmap to
derive the correct input-output relations for optical cavities comprised of
either one or two of such end-mirrors. We remark that this is not limited to
subradiant emitter arrays but also extendable to patterned subwavelength
gratings or photonic crystal structures Miroshnichenko _et al._ (2010); Chang-
Hasnain and Yang (2012); Zhou _et al._ (2014); Limonov _et al._ (2017) and
semiconducting monolayers Zeytinoğlu _et al._ (2017); Back _et al._ (2018);
Scuri _et al._ (2018).
Finally, for ensembles of cavity-embedded quantum emitters, externally and
incoherently pumped, we provide a minimal theory of lasing to illustrate how
the laser threshold condition can be easily derived using the formalism
introduced in Sec. A and by going beyond the linear regime. In particular, we
analyze the physics of superradiant lasers, where coherence is being stored in
the gain medium instead of in the cavity photon field. Finally, we connect the
dynamics of such a laser to the Dicke superradiance model introduced in Sec.
D.

###  A Antiresonance spectroscopy with 1D or 2D arrays

Let us first consider a weakly reflecting array placed transversely to the
axis of a standard optical cavity. In the presence of near-field couplings,
the system undergoes dynamics describable at the level of averages by Eqs.
(61) and at the fully quantum level by Eqs. (67). We will assume a vector G
containing all the cavity-emitter couplings and later particularize to
symmetric coupling (where all couplings are equal to $g$) and
antisymmetrically-phased coupling (with alternating $(-1)^{j}g$ couplings).
Assuming a single cavity drive, one can then easily deduce the cavity
transmission as

$\displaystyle
t=\displaystyle\frac{\kappa}{i\Delta_{c}+\kappa+\textbf{G}^{\top}\textbf{G}/\left[i\Omega_{\text{eff}}(\Delta)+\gamma_{\text{eff}}(\Delta)\right]},$
(83)

where the effective $\Delta$-dependent collective energy shifts and linewidths
are derived from the matrix
$\textbf{M}(\Delta)=-i\Delta\mathbbm{1}+i\bm{\Omega}+\bm{\Gamma}$ as real and
imaginary parts

$\displaystyle\gamma_{\text{eff}}(\Delta)+i\Omega_{\text{eff}}(\Delta)=\frac{\textbf{G}^{\top}\textbf{G}}{\textbf{G}^{\top}\textbf{M}^{-1}(\Delta)\textbf{G}}.$
(84)

These expressions are the equivalents of the decay rates and shifts for free
space arrays given in Sec. C where the addressing of collective resonances is
now controlled via the choice of cavity couplings instead of the external
drive.

Figure 11: (a) Effective cooperativity $C_{\text{eff}}$ of a cavity-embedded
chain of emitters as a function of $\mathcal{N}$ for the superradiant (red,
symmetric addressing $g_{i}=g$) and subradiant (blue, antisymmetric addressing
$g_{i}=(-1)^{i}g$) case. Parameters: $a=0.1\lambda$, $\gamma=\kappa/20$,
$g=\kappa/100$. The bottom half shows the cavity transmission functions for
(b) antisymmetric and (c) symmetric excitation of four equally spaced emitters
for $a=0.08\lambda$, $g=\kappa/10$. The gray dashed curves shows the result
for independent emitters. In all plots, the dipole moments are chosen
perpendicular to the chain.

Enhanced cooperativity \- Notice that one can now proceed by introducing a
modified $\mathcal{N}$ emitter effective cooperativity