# Subradiant entanglement in plasmonic nanocavities

Kalun Bedingfield  School of Physics and Astronomy, University of Birmingham,
Edgbaston, Birmingham, B15 2TT, United Kingdom  Benjamin Yuen  School of
Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham, B15
2TT, United Kingdom  Angela Demetriadou<EMAIL_ADDRESS>School of
Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham, B15
2TT, United Kingdom

###### Abstract

Plasmonic nanocavities are known for their extreme field enhancement and sub-
wavelength light confinement in gaps of just a few nanometers. Pairing this
with the ability to host quantum emitters, they form highly promising
platforms for controlling and ‘engineering’ quantum states at room
temperature. Here, we show how sub-radiant entangled states emerge between two
or more quantum emitters within plasmonic nanocavities. We develop a
theoretical description that directly links quantum variables to
experimentally measurable quantities, such as the extinction cross-section. We
show that the lossy nature of plasmonic nanocavities aids the emergence of
sub-radiant entangled states between quantum emitters, persisting for $\sim
100$ times longer than the plasmonic excitation. This work paves the way
towards designing and ‘engineering’ quantum entangled states in ambient
conditions with plasmonic nanocavities, for potential applications as rapid
quantum memories, quantum communications and quantum sensors.

Plasmonic nanocavities are formed by two or more nano-metallic structures with
a nanometer sized gap between them [1, 2, 3]. They have the ability to capture
incident electromagnetic fields and confine their energy at extremely sub-
wavelength volumes, whilst at the same time efficiently radiate energy to the
far-field [4, 5, 6, 7, 8]. A characteristic example is the nanoparticle-on-
mirror (NPoM) geometry, which consists of a spherical nanoparticle assembled
on a flat gold surface, creating a robust sub-nanometer gap between them
(Figure 1A). It has drawn significant attention in recent years due to its
extremely high field confinement and ability to reveal light-matter
interactions in ambient conditions. For example, single molecule light-matter
strong coupling at room temperature was experimentally realised in a NPoM
system [5], giving rise to coherent cyclic energy exchange between the plasmon
(light) and quantum emitter(s) (matter). Since then, these extreme plasmonic
nanocavities have led to a plethora of other novel realisations [9, 10], such
as imaging of a single molecule [6, 7], the interrogation of specific chemical
bonds within a single molecule, the formation of ‘pico-cavities’ [8] and
ultrafast single photon emission at room temperature [11].

So far, light-matter interactions of quantum emitters (QEs) in nanoplasmonic
devices have predominantly been studied experimentally, where far-field
measurements (i.e. scattering and extinction cross-sections) are usually
obtained for the combined system [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. The
quantum dynamics of the system are usually obtained with theoretical
descriptions, and contain vital information that can guide us on engineering
quantum states at room temperature [14, 15], and allow for more insightful
conclusions to be drawn from experiments. However, such theoretical quantum
descriptions face significant challenges in describing realistic systems that
correlate directly with experimental measurements. Plasmonic devices often
support multiple, spectrally close and broadband modes that cannot be easily
decomposed due to the dispersive nature of gold, with some recent attempts
trying to address this problem [16, 17, 18, 19, 20, 21, 22]. Furthermore, the
highly radiative and dissipative nature of nanoplasmonic devices creates a
major challenge in quantizing plasmonic modes, since energy is not conserved
during electromagnetic oscillations. Nevertheless, several methods have been
proposed so far that try to link the classical plasmonic behaviour to a
quantized description, with each method usually addressing one of these
problems [23, 24, 25, 26, 27, 28, 29, 30, 31, 32]. Their common denominator is
their tendency to focus on quantum scattering and spontaneous emission (i.e.
Purcell factor), so in general do not offer the quantum dynamics of the system
and cannot be extended to multiple QEs. Other theoretical methods that provide
the quantum dynamics of QEs in plasmonic devices [28], have so far been
limited to rapidly decaying systems, which is hindering our attempts to design
and build quantum devices with nanoplasmonics operating in ambient conditions.

In this paper, we propose a theoretical description of an open quantum system
consisted of multiple QEs placed in a plasmonic nanocavity that directly links
quantum variables to experimental measurable quantities, and we reveal semi-
persistent sub-radiant entanglement between two or more QEs. This sub-radiant
entanglement lives $\sim 100$ times longer than the lifetime of Rabi
oscillations and actually emerges because of the lossy nature of plasmons. It
is independent of the QEs’ locations in the nanocavity, as long as they are
strongly coupled to the plasmonic mode, and their decay is only dependent on
the decoherences of QEs alone. We use a Quasi-Normal mode (QNM) analysis
method [19, 33] to decompose and quantize the nanocavity’s plasmon modes to
form a cavity Quantum Electrodynamics (cQED) description of the open quantum
system, which is coupled to external reservoirs to represent the losses and
decoherences, and obtain the system’s full range of quantum electrodynamics.
Furthermore, we derive observable quantities from our quantum model that are
directly linked to experimental measurements and ensure these are consistent
in the classical limit. Hence, this work paves the way towards engineering
quantum states [34] and building quantum plasmonic devices operating in
ambient condition.

Figure 1:  (a) NPoM geometry for a $40$nm radius spherical gold nanoparticle
assembled $d=1$nm above an infinite flat gold mirror, separated by a spacer of
refractive index $n=2.5$. (b) NPoM schematic showing the coupling of energy in
and out of the system with rates $\kappa_{in}$ and $\kappa_{out}$,
respectively, and shows two emitters placed at different locations. (c) The
normalised mode profile of the $(1,0)$ mode on a plane through the centre of
the cavity, with the white and red dots respectfully indicating the centre of
the cavity and the dependence of the coupling strength on the emitter
positions. (d) The comparison of the classical absorption, scattering and
extinction cross-sections with the quantum extinction cross-section for a
single resonance model.

##  1 Quantum description of plasmonic nanocavities

The open and lossy nature of plasmonic environments creates a non-Hermitian
system and places significant theoretical impediments for using common cQED
descriptions that were developed for high-finesse cavities. Recently, several
quantum descriptions have been proposed that aim to overcome the non-Hermitian
nature of the problem [23, 24, 29]. Even though each method is a significant
breakthrough, they often focus on quantum scattering processes and spontaneous
emission, by presenting the Purcell factor modifications due to the QE
interaction with the plasmons [29, 30, 31, 32], with no explicit quantum
dynamics provided, or focus on rapidly decaying systems [29, 35, 36]. An
exception is the Weisskopf-Wigner theory that was recently employed to obtain
the quantum dynamics [23, 25, 24, 26, 27, 28], but one cannot easily
generalize this method to describe multiple QEs or entangled states.
Alternatively, one can use the semi-classical Maxwell-Bloch description [5,
37, 38, 39], but it only describes the problem within the first-quantization,
and therefore cannot describe spontaneous emission, entanglement or any non-
classical light that might emerge from nanoplasmonic devices. Hence, to
describe correlations or potential entanglement between multiple emitters
within a plasmonic nanocavity, one needs to introduce new theoretical
descriptions.

Here, we introduce an open quantum system formalism that allows us to obtain
the quantum dynamics of multiple QEs placed in a plasmonic nanocavity, driven
by an external source and access experimentally measurable quantities, such as
the extinction cross-section. We decompose the classical photonic response of
the plasmonic system into quasi-normal modes (QNMs), using the QNMEig
methodology [19, 33]. This method employs an auxiliary-field eigenvalue
approach to account for the dispersive nature of the metals in a linear
manner. The eigenvectors provide the near-field of each plasmon mode with
complex eigenfrequencies [40]
$\tilde{\omega}_{cav}=\omega_{cav}-i\kappa_{out}/2$, where $\omega_{cav}$ is
the cavity resonant frequency and $\kappa_{out}$ the total decay rate
(radiative and Ohmic losses) of each plasmonic mode. We then build a cavity
Quantum Electrodynamics (cQED) description with input and output fields, so
unlike previous methods we can relate to experimental measurements. Here, we
consider the Nano-Particle on Mirror (NPoM) plasmonic nanocavity as shown in
Figure 1(a)-(b). To ease the interpretation of the system’s quantum dynamics
and the correlations between multiple QEs residing within the plasmonic
nanocavity, we design a gold NPoM nanocavity for which it is safe to assume
that only one plasmon mode couples to the QEs. Our method though can be easily
extended to include many plasmonic modes.

We divide the Hamiltonian $\mathcal{H}=\mathcal{H}_{sys}+\mathcal{H}_{res}$
into two components: the system ($\mathcal{H}_{sys}$) and reservoir
($\mathcal{H}_{res}$) Hamiltonians, which we discuss separately. The system
Hamiltonian describes the plasmonic modes obtained from classical QNM
calculations [19] ($\mathcal{H}_{cav}$), the multiple QEs as two-level systems
($\mathcal{H}_{QEs}$), the interaction between the plasmon modes and the QEs
($\mathcal{H}_{int}$) and a driving external field coupled to the
cavity($\mathcal{H}_{ext}$), such as a pump, that here represents an external
illumination of the system. The system Hamiltonian is given by:

$\mathcal{H}_{sys}=\underbrace{\Delta_{cav}a^{\dagger}a}_{\mathcal{H}_{cav}}+\underbrace{\sum_{j}^{N}\frac{1}{2}\Delta_{d}\sigma^{(j)}_{z}}_{\mathcal{H}_{QEs}}+\underbrace{\sum_{j}^{N}g(\mathbf{r}^{(j)})\left(a^{\dagger}\sigma^{(j)}_{-}+a\sigma^{(j)}_{+}\right)}_{\mathcal{H}_{int}}+\underbrace{i\sqrt{\kappa_{in}}\left(\alpha^{*}a-\alpha
a^{\dagger}\right)}_{\mathcal{H}_{ext}}$ (1)

where $N$ is the number of QEs placed in the nanocavity,
$\Delta_{cav}=\omega_{cav}-\omega_{p}$ and $\Delta_{d}=\omega_{QE}-\omega_{p}$
are the detunings of the cavity resonance ($\omega_{cav}$) and each QE
resonant frequency ($\omega_{QE}$) from the external field’s frequency
($\omega_{p}$). The $\left\\{a,a^{\dagger}\right\\}$ are the bosonic operators
that create and annihilate the plasmon mode, and $\sigma_{z}^{(j)}$ is the
Pauli $z$-operator describing the $j$-th two-level QE (see Supp. Info. for
more information). The interaction of the plasmonic mode with the $j$-th QE at
position $\mathbf{r}^{(j)}$ is given by
$\left(a^{\dagger}\sigma^{(j)}_{-}+a\sigma^{(j)}_{+}\right)$ with the coupling
strength measured by:

$g(\mathbf{r}^{(j)})=\sqrt{\frac{\hbar\omega_{cav}}{2\varepsilon_{0}\varepsilon
V}}\ \mathbf{d}^{(j)}\cdot\mathbf{u}(\mathbf{r}^{(j)})$ (2)

where $\mathbf{u}(\mathbf{r}^{(j)})$ is the normalized electric field vector
of the QNM mode, which makes the coupling strength $g(\mathbf{r}^{(j)})$
dependent on the $j$-th QE’s location $\mathbf{r}^{(j)}$ within the plasmonic
nanocavity. The mode volume is obtained from the QNM calculations as
$V=\frac{\int
d\mathbf{r}|\mathbf{u}(\mathbf{r})|^{2}}{\text{Max}[|\mathbf{u}(\mathbf{r})|^{2}]}$
(see Supp. Info.), $\varepsilon$ is the electric permittivity of the material
hosting the QEs and $\mathbf{d}^{(j)}$ is the dipole moment of the $j$-th QE.
The last term ($i\sqrt{\kappa_{in}}\left(\alpha^{*}a-\alpha
a^{\dagger}\right)$ ) of equation (1) describes an external coherent
monochromatic source driving the system, where $\sqrt{\kappa_{in}}$ is the
rate that energy couples into the system and $\alpha$ is the amplitude of the
coherent state, defined by the source’s photon flux $c|\alpha|^{2}$ (see Supp.
Info. for more discussion).

While the system Hamiltonian describes the excitation of the system and the
interaction between QEs and the plasmon mode, it contains no description of
losses. Nanoplasmonic devices hosting QEs in general have two loss channels:
(i) the plasmon modes loose energy by both radiating to the far-field and via
Ohmic losses, and (ii) the decoherences of QEs due to the (ro)-vibrational
energy states of the molecule. We represent the plasmonic loss and QE
decoherences with two separate reservoirs as:

$\displaystyle\mathcal{H}_{res}=$
$\displaystyle\sum^{N}_{j}\int^{\infty}_{0}d\omega_{vib}\
\sigma^{(j)}_{z}\kappa_{vib}(\omega_{vib})\left[b_{vib}^{\dagger}e^{i\omega_{vib}t}+b_{vib}e^{-i\omega_{vib}t}\right]$
(3)
$\displaystyle-i\sum_{\lambda}\int^{\infty}_{0}d\omega\kappa_{res}(\omega)\left[b_{\lambda}(\omega)a^{\dagger}e^{i\Delta(\omega)t}-b^{\dagger}_{\lambda}(\omega)ae^{-i\Delta(\omega)t}\right]$

where the first term describes the coupling of the QEs to a phonon bath
reservoir, with phonon frequencies $\omega_{vib}$ and strength
$\kappa_{vib}(\omega_{vib})$, and $\left\\{b_{vib},b_{vib}^{\dagger}\right\\}$
the phonon reservoir operators. The second term of equation (3) accounts for
the coupling of the plasmonic mode to a background plasmon reservoir, with
photon frequencies $\omega$ and coupling strength $\kappa_{res}(\omega)$, to
account for the energy lost both radiatively and via Ohmic losses. Note that
the detuning $\Delta(\omega)=\omega-\omega_{p}$ is the detuning of frequency
$\omega$ with the incident external field frequency ($\omega_{p}$), and
$\left\\{b_{\lambda}(\omega),b_{\lambda}^{\dagger}(\omega)\right\\}$ are the
plasmon reservoir creation/annihilation operators.

The time evolution of the open quantum system with the plasmon loss and QE
decoherences is described by the Lindblad master equation as:

$\displaystyle\frac{d}{dt}\rho=$
$\displaystyle-i\left[\mathcal{H}_{sys},\rho\right]+2\kappa_{out}a\rho
a^{\dagger}-\kappa_{out}\left(a^{\dagger}a\rho+\rho a^{\dagger}a\right)$
$\displaystyle+\sum^{N}_{j}\left[\kappa_{vib}\sigma^{(j)}_{z}\rho\sigma^{(j)}_{z}-\frac{\kappa_{vib}}{2}\left(\sigma^{(j)2}_{z}\rho+\rho\sigma^{(j)2}_{z}\right)\right]$
(4)

where the density operator $\rho$ describes the statistically mixed state of
the system. The expectation value of any operator of the system can be
obtained via $\langle O\rangle=\text{Tr}[O\rho]$, for example the expected
plasmon excitation number $\left\langle a^{\dagger}a\right\rangle$.

###  1.1 Quantum Extinction Cross-Section

One can use the above cQED description to obtain classical quantities that are
commonly measured experimentally, such as the extinction cross-section. First,
to represent an external plane wave incident on the system, we use an external
driving field that continuously illuminates the system. This external field
originates by singling out an input field mode from the plasmon reservoir and,
under the mean-field approximation, replacing the field operator
$b_{\lambda}(\omega)$ with a complex number $\alpha$ that directly relates to
the photon flux of the field (see Supp. Info.). As the reservoir describes the
rest of the universe, we can assume that the reservoir remains unchanged due
to the quantum statistics of the external field mode not significantly
affecting the statistics of the reservoir. The rates with which energy is
coupled into ($\kappa_{in}$) and lost ($\kappa_{out}$) from the system is
schematically shown in Figure 1b and can be used to find the quantum
extinction cross-section. To compare our results with the classical extinction
cross-section of the NPoM, we remove the QEs and their interaction with the
plasmon mode from the Hamiltonian, solve equation (4) analytically in steady-
state by setting $d\rho/dt=0$, and take the solutions for coherent states.
This gives us the quantum extinction cross-section, which we then compare to
the classical extinction cross-section, to obtain the correct input coupling
parameter $\kappa_{in}$.

Classically, the extinction cross-section is defined as the total power lost
by the system ($P_{out}$) relative to the incident intensity ($S$) of the
plane wave. When considering their quantum analogues in steady-state, they
respectively become: $P_{out}=2\kappa_{out}\hbar\omega\langle
a^{\dagger}a\rangle$ and $S=c_{0}\hbar\omega|\alpha|^{2}$. The $\langle
a^{\dagger}a\rangle=\frac{\kappa_{in}|\alpha|^{2}}{\kappa_{out}^{2}+\Delta_{cav}^{2}}$,
and is proportional to $\langle a\rangle$ and $\langle a^{\dagger}\rangle$
(see Supp. Info.), which leads to the quantum extinction cross-section:

$\displaystyle\langle\sigma_{ext}^{q}\rangle_{noQEs}=\frac{2\kappa_{in}\kappa_{out}}{c_{0}}\frac{1}{\kappa_{out}^{2}+\Delta_{cav}^{2}}\
,$ (5)

where $c_{0}$ is the speed of light constant. We define the in-coupling
coefficient as: $\kappa_{in}=\kappa_{out}c_{0}\sigma_{ext}/2$, which ensures
that the quantum, $\langle\sigma_{ext}^{q}\rangle_{noQEs}$ and classical,
$\sigma_{ext}^{c}$, extinction cross-sections are equal
$\langle\sigma_{ext}^{q}\rangle_{noQEs}=\sigma_{ext}^{c}$ when the external
illumination is on resonance with the QNM mode (i.e. $\omega_{cav}=\omega_{p}$
and $\Delta_{cav}=0$).

Now, we apply this theoretical description to a specific NPoM nanocavity.
Throughout this article, we consider a NPoM geometry (see figure 1a-b) made of
a spherical, gold nanoparticle of radius $40\,nm$ assembled $1\,nm$ above a
flat gold mirror, separated by a molecular monolayer with refractive index
$n=2.5$, which is consistent with recent experimental set-ups [41, 42]. This
plasmonic nanocavity has the first order $(1,0)$ quasi-normal mode resonant at
$\tilde{\omega}_{cav}=(338.9-i9.8)THz$ (i.e. $\lambda_{cav}\sim 900\,nm$). The
normalised electric field eigenvector $\mathbf{u}(\mathbf{r})$ for the $(1,0)$
mode is shown on a slice through the centre of the cavity in Figure 1c, and
exhibits a large central extrema, with the mode volume being
$\text{Re[}V\text{]}=60.01$nm3.

In Figure 1d, we plot the classical extinction ($\sigma_{ext}$), scattering
($\sigma_{scat}$) and absorption ($\sigma_{abs}$) cross-sections that show the
lowest order mode to be spectrally isolated, and the $\kappa_{out}$ obtained
from the QNM calculations corresponding to the full-width half-maximum of
$\sigma_{ext}$. Therefore we can safely assume that only one plasmonic mode
interacts with the QEs for this specific plasmonic nanocavity. The quantum
extinction cross-section is in good agreement with its classical counterpart,
and since we considered only one plasmonic mode for the cQED description, it
only shows one peak. Note that the values of the quantum and classical
$\sigma_{ext}$ match exactly on resonance with the cavity (i.e. the $(1,0)$
QNM mode: $\lambda_{(1,0)}\sim 900\,$nm). The spectral shifts seen in Figure
1D, are due to the broadband contributions of higher order modes for far-field
quantities such as the classical $\sigma_{ext}$, and have been previously
reported and explained [43, 44]. With all the physical parameters of the
system obtained, our open quantum model now describes the interaction between
the NPoM geometry and multiple QEs placed at any position within the plasmonic
nanocavity. This driven formalism is also very flexible, allowing us to
initialize the system in any particular state.

##  2 Two-Emitter Quantum Dynamics

We first consider a single two-level QE placed at the centre of the cavity
($x=0$nm). Note that all QEs considered in this paper are Cy5 molecules with a
dipole moment $d=10.1$D and excited state decoherence rate of
$\kappa_{vib}=25$meV [42]. With the plasmon state initially excited and no
external illumination present, energy is passed between the plasmon and
emitter for several cycles (i.e. Rabi oscillations) before being dissipated
(either radiatively or via the non-radiative processes outlined above), and
the system eventually returns to its ground state (see Supp. Info.). These
Rabi oscillations reveal that $g^{2}>\kappa_{out}\kappa_{vib}$ and the system
is therefore comfortably within the strong coupling regime, with the QE
linewidth for its excited state ($\kappa_{vib}$) adding a ‘blurring’ effect on
the oscillations (see Supp. Info.).

Now we place a second QE (Cy5 molecule) in the nanocavity and for simplicity
initially assume that both QEs are at the centre of the nanocavity to ensure
that both have the same coupling strength with the plasmon mode. Both QEs show
identical Rabi oscillations, as expected, similar to the case of one emitter
at the nanocavity centre, but with the Rabi frequency increased by a factor of
$\sqrt{N}=\sqrt{2}$, as expected. However, the increased Rabi frequency is
limited by the total decay rate of the system, and instead follows
$\Omega_{Rabi}^{2}=4g^{2}-\kappa_{out}^{2}$ [45]. Although the system is
further into the strong coupling regime than the single emitter case, it is
worth noting that the population of the two emitter system decays faster; this
is due to the additional factor of $\kappa_{vib}$ supplied by the second
emitter, but also to the increased Rabi frequency (i.e. the interactions and
energy exchanges happen on shorter timescales).

Figure 2: The excited emitter state population dynamics of the NPoM geometry
interacting with two emitters, for symmetric (A-C) and asymmetric (D-G)
emitter arrangements—initialised in the excited field state. Schematics for
the symmetric (A) and asymmetric (D) emitter placements, along the $x$-axis on
a plane through the centre of the cavity. Population dynamics as a function of
emitter positioning for the: (B) degenerate symmetric emitters; (E) central
emitter and (F) displaced emitter for the asymmetric arrangement. Collective
coupling strength as a function of emitter positioning for the symmetric (C)
and asymmetric (G) arrangements.

However, in realistic systems, the two Cy5 molecules cannot be placed at the
same position within the nanocavity, and their coupling strength is position
dependent, according to the QNM mode’s eigenvector shown in Figure 1c. Due to
the cylindrical symmetry of the cavity and of the $(1,0)$ mode, it makes sense
to consider QE positions only along the $x$-axis—as indicated by the white
dashed line in Figures 2a. Figure 2a diagrammatically shows a symmetric
arrangement for the two QEs away from the centre of the cavity, which means
both QEs continue to experience the same coupling strength, producing
identical excited state populations for all positions $|x|$. These are shown
in Figure 2b for $|x|$ up to $10$nm away from the nanocavity centre. As the
emitters move away from the cavity centre, they experience a rapid reduction
of their coupling strength, and for $|x|>6$nm the system transitions into the
weak coupling regime. Another way of quantifying this is through the ratio of
an effective coupling strength $g_{E}$ relative to the coupling of a single
emitter at the centre of the cavity, $g_{0}$, which is shown in Figure 2c. Due
to the dependence of the coupling strength on the mode profile, the effective
coupling strength ($g_{E}$) changes with the position of the emitters as:
$g_{E}/g_{0}\propto\sqrt{u(x_{1})^{2}+u(x_{2})^{2}}$, which becomes
$\sqrt{2}u(x)$ for the symmetric placement of the QEs (since
$u(x_{1})=u(x_{2})$), and at $x_{1}=x_{2}=0$ the effective coupling strength
is enhanced by $\sqrt{2}$.

If instead the two QEs are placed asymmetrically within the plasmonic
nanocavity, then each QE is coupled to the plasmon mode with a different
coupling strength. We take the case where one QE always remains at the centre
of the cavity, while the second QE takes various positions along the
$x$-axis—as shown in Figure 2d. The population dynamics of the central emitter
and displaced emitter now differ, and are respectively shown in Figure 2e and
f. Nevertheless, we find that the Rabi oscillations are still governed by the
effective coupling strength
$g_{E}/g_{0}\propto\sqrt{u(x_{1})^{2}+u(x_{2})^{2}}$ (see Figure 2g). This now
decreases more slowly as $\sqrt{1+u(x_{2})^{2}}$ compared to the symmetric
case, and reduces to unity once the second QE is only weakly coupled to the
plasmon, which occurs at $x_{2}\sim 10\,nm$. Also, the central QE’s Rabi
oscillations slow down as the second emitter is placed further out from the
nanocavity, but its excited population increases as the system returns to a
single emitter interaction, which occurs when the displaced QE is at
$x_{2}\sim 10\,nm$.

###  2.1 Hybrid states

Experimental measurements of systems with QEs in plasmonic nanocavities
usually record the hybrid states of the system in the far-field. Usually an
external field illuminates the combined system, with the scattered light (or
extinction energy) been collected and analysed to deduce the hybrid states and
the interaction between the plasmon and QEs. As we said earlier, the power
lost from the cavity is: $P=2\kappa_{out}\hbar\omega\langle
a^{\dagger}a\rangle$, which gives the rate of photon lost. Due to reciprocity,
this energy loss in steady state is equivalent to the plasmon excitation in
the system by the external illumination, which is described by:
$P=-\sqrt{\kappa_{in}}(\alpha\langle a^{\dagger}\rangle+\alpha^{*}\langle
a\rangle)$. Normalizing this energy loss with the incident energy
$S=c_{0}\hbar\omega|\alpha|^{2}$ leads to the quantum extincion cross-section
for the combined system of QEs in a plasmonic nanocavity:
$\langle\sigma_{ext}^{q}\rangle=-\sqrt{\kappa_{in}}(\alpha\langle
a^{\dagger}\rangle+\alpha^{*}\langle a\rangle)/c_{0}|\alpha|^{2}$, which is in
terms of the field operators. Without loss of generality, $\alpha$ can be
taken to be real in our case, which reduces the quantum extinction cross-
section to:

$\displaystyle\langle\sigma^{q}_{ext}\rangle=-2\sqrt{\kappa_{in}}\frac{Re\\{\langle
a\rangle\\}}{c_{0}\alpha}$ (6)

where $\langle a\rangle$ describes the rate at which energy is removed from
the excited plasmon state (i.e. extinction energy) and can be found by first
solving the master equation for the density matrix in steady state, and then
taking the trace with the field operator as: $\langle
a\rangle=\text{Tr}[a\rho]$.

Figure 3: (A) The $\langle a\rangle$ of the driven single emitter NPoM system
as a function of both cavity
($\Delta_{cav}=\frac{\omega_{cav}-\omega_{em}}{\omega_{em}}$) and external
source ($\Delta_{p}=\frac{\omega_{p}-\omega_{em}}{\omega_{em}}$) detunings.
The black and white dashed lines respectively show $\Delta_{cav}=\Delta_{p}$
and $\Delta_{p}=0$, highlighting the two resonances which cross in a weakly
coupled regime but form an avoided crossing here. Rabi splitting in terms of
$\langle a\rangle$ as a function of both plasmon detuning and coupling
strength, for (B) one and (C) two emitters at the centre of the cavity. Grey
dashed lines trace the peaks of these split resonances, and the solid grey
lines follow $\omega_{\pm}$ [20]. In (B) and (C), $g_{0}=1$ and
$g^{\prime}=\sqrt{2}$ respectively mark the collective coupling strengths in
these systems.

For this externally driven regime, the system which is initially in the vacuum
state, is excited at a rate $\sqrt{\kappa_{in}}\alpha$—such that the average
photon number in the system is $N_{ph}=10^{-6}$. In Figure 3a, we plot the
$\langle a\rangle$ (which is proportional to $\langle\sigma_{ext}^{q}\rangle$)
for a single emitter at the centre of the nanocavity, with both the external
source ($\Delta_{p}$) and cavity ($\Delta_{cav}$) resonances detuned relative
to the QE resonance. This exhibits the signature avoided crossing to form two
hybrid states, which is characteristic of strong coupling [46, 47]. The black
and white dashed lines respectively show $\Delta_{cav}=\Delta_{p}$ and
$\Delta_{p}=0$. Many experimental analyses of strongly coupled systems measure
these hybrid states, with other theoretical descriptions also predicting
them[38, 26, 20]. Actually, it was recently shown [20] that one can consider a
complex coupling strength, which accounts for the dispersive nature of metals
and the complex mode volume, to obtain the complex frequency of the strong
coupling hybrid states as [20, 48, 49, 50]

$\displaystyle\tilde{\omega}_{\pm}=\frac{\tilde{\omega}_{e}+\tilde{\omega}_{cav}}{2}\pm\sqrt{g^{2}[1-iR]+\left(\frac{\tilde{\omega}_{e}-\tilde{\omega}_{cav}}{2}\right)^{2}}$
(7)

where $R=\text{Im[}V\text{]/Re[}V\text{]}$ is the ratio of the imaginary and
real components of the mode volume [18], and the imaginary components of
$\tilde{\omega}_{e}$ and $\tilde{\omega}_{cav}$ are respectively
$\kappa_{vib}$ and $\kappa_{out}$. For systems with a large quality factor,
$Q$, the negligibly small $\text{Im[}V\text{]}$ makes it apparent that the
standard description may be sufficient. For low Q systems though, such as our
plasmonic nanocavity, the complex nature of the mode volume has the ability to
significantly alter the observed dynamics. The solutions for
$\tilde{\omega}_{\pm}$ are plotted with grey lines in Figure 3a, demonstrating
their strong agreement with our non-Hermitian cQED description for light-
matter interactions.

Now, we obtain the strong coupling hybrid states for various coupling
strengths with respect to the external field detuning from the QE $\Delta_{p}$
(see Figure 3b). The Rabi splitting arises due to the coupling strength
exceeding the loss rate of the system, leading to the coherent exchange of
energy between the plasmon and the QE [51]. For an empty cavity ($g=0$), a
singular resonance is observed which is attributed to the absorption and re-
emission of the energy collected by the plasmonic mode alone. Increasing the
coupling between the QE and the plasmon, one quickly sees their hybridisation
and emergence of the characteristic Rabi splitting, even for small coupling
strengths. For small coupling strengths the strong coupling condition
$g^{2}>\kappa_{out}\kappa_{vib}$ is still valid, despite the QE population
dynamics showing a relatively small number of Rabi oscillations that dissipate
very quickly. However, to acquire a more intuitive picture of the occurring
interaction, it is more appropriate to separate this condition into two: one
for the plasmon $g>\kappa_{out}$ and one for the QE $g>\kappa_{vib}$. For
example, a single emitter at the centre of our NPoM cavity (i.e. at $g_{0}$)
has: $g_{0}^{2}/\kappa_{out}\kappa_{vib}=4.64$, placing the system well into
the strong coupling regime. With $g_{0}/\kappa_{vib}=9.74$, which allows the
Rabi splitting to occur and $g_{0}/\kappa_{out}=0.48$ that leads to a large
radiative emission and rapidly decaying Rabi oscillations. To lose the Rabi
splitting and the hybrid states for the overall system, the coupling strength
needs to reduce enough for the ratio $g/\kappa_{vib}<1$, which is in agreement
with the results shown in Figure 3(b). This is also in agreement with the Rabi
splitting frequencies obtained from equation (7) [20] (plotted with grey
lines), yielding the same Rabi splitting for low coupling strengths. For
multiple QEs in the cavity, the coupling strength scales with $\sqrt{N}$ which
causes the plasmon and QEs to hybridise even more intensely, increasing the
Rabi splitting. Figure 3(c) shows the Rabi splitting for varying coupling
strengths of two emitters in a plasmonic nanocavity, with the $\omega_{\pm}$
solutions of equation (7) also in agreement (plotted with grey lines).

##  3 Sub-radiant Entanglement

Figure 4: The population dynamics for the NPoM geometry coupling to two
emitters—initialised in one of the excited emitter states. S and A
respectively describe the super- and sub-radiant superposition states. The
dynamics for the emitters placed at the centre of the cavity—with (B) and
without (A) the QE linewidth ($\kappa_{vib}$) included in the model. The
symmetric (C) and asymmetric (D) placement of the two emitters, at a distance
of $x=|4|$nm.

We now take advantage of the flexibility of our cQED description to consider
no external illumination driving the system, by setting the input field
$\alpha$ to zero, such that the system is free to be initialised in any
particular state. This allows the free decay dynamics of specific light-matter
interactions to be probed. We specifically choose the initialisation of one QE
in the excited state and the other in the ground state. To better demonstrate
the unique dynamics presented by this preparation, Figure 4A ignores the
inherent molecular linewidths of the emitters ($\kappa_{vib}=0$), and assumes
both emitters are placed at the centre of the cavity. The primary energy
exchange and Rabi oscillations are now observed between the two emitters—with
twice as much energy passed to the initially non-excited emitter than to the
plasmonic mode. Due to the asymmetric initialisation, the emitters no longer
interact with the field coherently: they instead oscillate out of phase with
each other as energy is passed between them. The key observation here is that
the emitters exchange energy until they reach an equal, non-zero steady state
value. This long lived state originates from the symmetric and anti-symmetric
superpositions formed from the states of the two emitters:

$\displaystyle|S\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|0,e,g\rangle+|0,g,e\rangle\right)$ (8)
$\displaystyle|A\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|0,e,g\rangle-|0,g,e\rangle\right)\ ,$
(9)

which respectively describe the super- and sub-radiant states, as shown in
Figure 4A with dot-dashed and dotted lines. Although the super-radiant state
rapidly decays to zero, the sub-radiant state reaches a steady state solution.
This occurs due to the asymmetry of the state, with the interactions of the
$|1,g,g\rangle$ state with the components of $|A\rangle$ ($|0,e,g\rangle$ and
$|0,g,e\rangle$) perfectly cancelling out. In fact, these states emerge
because of the high plasmon losses, causing the super-radiant state to decay
very rapidly, with the plasmon providing a conduit to form the entangled
state. To determine the impact of QE decoherences on these sub-radiant
entangled states, we introduce the decoherence of Cy5 molecules
($\kappa_{vib}=25$meV) and plot the same quantum dynamics in Figure 4B. These
states are formed, but now decay non-radiatively via $\kappa_{vib}$ and
eventually the system returns to its ground state. However, they live for
$\sim 100$ times longer than the plasmonic excitation, creating semi-
persistent sub-radiant entangled states in lossy plasmonic nanocavities.

For a more realistic arrangement of the emitters, Figure 4C and D respectively
refer to the symmetric and asymmetric placement of the emitters at $x=|4|$nm.
The equal coupling strength for the symmetric QE arrangement leads to energy
exchange between the two QEs of the same period that are perfectly out of
phase, to eventually form long lived entangled states. In fact, the long-lived
sub-radiant state has the same population as with both QEs at the centre of
the nanocavity (Supp. Info.). The only difference between the two cases
(Figure 4B and C) is the number of oscillations between the two emitters
before the system reaches the sub-radiant state. In contrast, the asymmetric
QE placement leads to unequal coupling strengths with the plasmon, and
different rates for the energy exchange between them, which distorts the sub-
radiant entangled states. This is an interesting conclusion to draw, that the
coherence between emitters and the formation of the sub-radiant entangled
states is mainly dependent on their relative coupling strength with each
other, and less on the coupling strength value that the nanocavity provides.

##  4 Multi-Emitter systems

Figure 5: The excited emitter state population dynamics of the NPoM geometry
interacting with eight emitters, for ring (A-B) and encircled emitter (C-E)
arrangements, initialised in the excited field state. Schematics for the ring
(A) and encircled emitter (C) placements, at a radius $r$ on a plane through
the centre of the cavity. Population dynamics as a function of emitter
positioning for the: (B) degenerate ring of emitters; (D) central emitter and
(E) degenerate surrounding emitters of the encircled emitter arrangement.

Of course one can have more than just two QEs in a plasmonic nanocavity, and
with DNA origami [41, 42] can be placed with a spatial precision of $\pm 1$nm.
Here we choose eight QEs ($N=8$), but keep the discussion generalised as it is
applicable for any reasonable number of emitters. We place the emitters in two
different spatial configurations. We initially consider a symmetric QE ring
configuration of radius $r$, shown diagrammatically in Figure 5A, where all
N-QEs experience the same coupling strength due to the cylindrical symmetry of
the plasmon mode. The quantum dynamics of this system can be described using
only three states: the vacuum state $|0,\vec{g}\rangle$, the excited field
state $|1,\vec{g}\rangle$, and an $N$-fold degenerate state with any one of
the $N$ emitters excited $|0,\phi\rangle$. Here
$|\vec{g}\rangle=\prod^{N}_{i}\otimes|g_{i}\rangle$ is a state in which all of
the emitters are in their ground states, and
$|\phi_{j}\rangle=\sigma_{+}^{(j)}|\overrightarrow{g}\rangle$ is the state
where emitter $j$ is excited whilst the rest remaining in their ground states.
The state $|0,\phi_{j}\rangle$ is shown in Figure 5B for this symmetric ring
configuration when the system is initialized with the plasmon in the excited
state; this exhibits a $\sqrt{N}$ enhancement to the effective coupling
strength, and therefore faster Rabi oscillations. However, each emitter has a
relatively low peak population due to the single plasmon quanta being shared
across all QEs. As before, the coupling strength and Rabi frequency decrease
as the emitters move away from the centre of the cavity; this continues until
the collective coupling is so weak that the Rabi oscillations are lost.

Now, instead consider the emitters placed in an encircled emitter
configuration, shown in Figure 5C, with one emitter at the centre of the
nanocavity and the other $N-1$ forming a ring of radius $r$. The central and
surrounding QEs experience different coupling strentghs and therefore have
different quantum dynamics, where $|\phi_{1}\rangle$ of the central emitter is
no longer degenerate with $|\phi_{j>1}\rangle$ of the ring emitters, shown in
Figure 5D and E respectively. Analogously to the asymmetric displacement of
the two emitter system, the surrounding emitters in this arrangement
experience a reducing Rabi frequency and population as $r$ increases; however,
they do so at a slower rate due to the maintained strong coupling of the
central, encircled emitter. This central emitter also experiences a reduced
Rabi frequency, but its population rises with $r$ as it approaches the single
emitter case.

We now explore the emergence of the sub-radiant entangled states for this
multi-emitter system. We initialize the system with one QE in its excited
state and show in Figure 6A that the sub-radiant states are formed, when all
N-emitters are at the centre of the nanocavity and there is no molecular
decoherence ($\kappa_{vib}=0$). Note that Figure 6 shows the population of the
state $|0,\phi_{1}\rangle$ for the initially excited QE (blue) together with
the total population of the other $(N-1)$ QEs (magenta). Figure 6B
incorporates the decoherence for Cy5 molecules of $\kappa_{vib}=25$meV, and
shows the emergence of the sub-radiant entangled states between the
N-emitters, that now decay more rapidly because of the higher total QE
decoherence. If the QEs are placed on a symmetric ring at $r=4$nm, the sub-
radiant entangled states are formed, but there are fewer oscillations between
the QEs since the coupling strength is much weaker. However, the coherent
dynamics of the QEs helps maintain their long-lived nature, with the sub-
radiant state population matching the same state’s population for all QEs at
the nanocavity centre. For the asymmetric encircled emitter configuration, the
system is initialized with the central emitter in an excited state, which
maintains the quantum dynamic degeneracy of the $(N-1)$ QEs forming the
circle. However, the varying coupling strength between the central and ring
QEs distorts the sub-radiant entangled states.

Figure 6: The population dynamics for the NPoM geometry coupling to eight
emitters—initialised in one of the excited emitter states (blue) and compared
to the collective excited state of the other emitters (magenta). The dynamics
for the emitters placed at the centre of the cavity—with (B) and without (A)
the QE linewidth ($\kappa_{vib}$) included in the model. The ring (C) and
encircled emitter (D) arrangements, at a radius $r=4$nm on a plane through the
centre of the cavity.

Hence, multiple QEs within a plasmonic nanocavity give rise to semi-persistent
sub-radiant states at room temperature, if the system is initialized
appropriately. These sub-radiant state form within just $\sim 100$fsec and
persist for $\sim 100$ times longer than the plasmonic excitation. They are
also formed regardless of the coupling strength value, as long as the QEs are
strongly coupled to the plasmon mode. The high loss of the plasmonic
nanocavity actually enables the emergence of the sub-radiant states, since it
causes the super-radiant state to decay very rapidly. Their decay depends only
on the molecular decoherence, and one can choose molecules, quantum dots or
other QEs with more favourable decoherence values. Such systems pave the way
towards engineering quantum states at room temperature with plasmonic
nanocavities, for quantum memory and quantum communication applications,
avoiding complex and cumbersome experimental set-ups.

##  5 Conclusion

We introduce a theoretical description suitable for multiple QEs in plasmonic
devices that allows to obtain both the quantum dynamics of the system and
experimentally measured quantities, such as the extinction cross-section.
Using this method, we demonstrate that multiple emitters in a plasmonic
nanocavity form semi-persistent sub-radiant states, if the system is
initialized appropriately. The sub-radiant states are independent of the high
plasmonic losses and the coupling strength, as long as the QEs are strongly
coupled to the plasmon, and survive for $\sim 100$ times longer than the
plasmonic excitation. This work paves the way towards building quantum
memories and quantum communication systems operating at room temperature.

## Acknowledgements

AD gratefully acknowledges support from the Royal Society University Research
Fellowship URF\R1\180097, Royal Society Research Fellows Enhancement Award RGF
\EA\181038, Royal Society Research grants RGS \R1\211093 and funding from
EPSRC for the CDT in Topological Design EP/S02297X/1.

## Supporting Information

The Supporting Information provides technical details on: (1) Mode volume; (2)
Hamiltonian for the driven cavity QED model; (3) Determination of the in-
coupling constant; (4) Quantum extinction cross-section; (5) Single emitter
dynamics; (6) Two emitter dynamics and (7) Multi-emitter dynamics.

## References

  * [1] A. Emboras, J. Niegemann, P. Ma, J. J. Baumberg, and D. R. Smith.  Atomic scale plasmonic switch.  Nano Letters, 16(4):709–714, 2016.
  * [2] D. O. Sigle, J. Mertens, L. O. Herrmann, J. J. Baumberg, and D. R. Smith.  Monitoring morphological changes in 2d monolayer semiconductors using atom-thick plasmonic nanocavities.  ACS Nano, 9(1):825–830, 2015.
  * [3] F. Benz, C. Tserkezis, L. O. Herrmann, J. J. Baumberg, and D. R. Smith.  Nanooptics of molecular-shunted plasmonic nanojunctions.  Nano Letters, 15(3):669–674, 2015.
  * [4] T. B. Hoang, G. M. Akselrod, C. Ciraci, J. Huang, D. R. Smith, and M. H. Mikkelsen.  Ultrafast spontaneous emission source using plasmonic nanoantennas.  Nature Communications, 6(1):7788, 2015.
  * [5] R. Chikkaraddy, B. De Nijs, F. Benz, A. Demetriadou, L. Marrucci, J. J. Baumberg, and D. R. Smith.  Single-molecule strong coupling at room temperature in plasmonic nanocavities.  Nature, 535(7614):127–130, 2016.
  * [6] Y. Zhang, Q. Meng, L. Zhang, J. Li, J. J. Baumberg, and D. R. Smith.  Sub-nanometre control of the coherent interaction between a single molecule and a plasmonic nanocavity.  Nature Communications, 8(1):15225, 2017.
  * [7] B. Yang, G. Chen, A. Ghafoor, X. Zhu, J. Zhang, J. J. Baumberg, and D. R. Smith.  Sub-nanometre resolution in single-molecule photoluminescence imaging.  Nature Photonics, 14(9):693–699, 2020.
  * [8] F. Benz, M. Schmidt, A. Dreismann, T. B. Hoang, J. J. Baumberg, and D. R. Smith.  Single-molecule optomechanics in “picocavities”.  Science, 354(6314):726–729, 2016.
  * [9] G. Zengin, M. Wersåll, S. Nilsson, T. J. Antosiewicz, M. Käll, and T. Shegai.  Realizing strong light-matter interactions between single-nanoparticle plasmons and molecular excitons at ambient conditions.  Phys. Rev. Lett., 114(15):157401, 2015.
  * [10] J. Huang, A.J. Traverso, G. Yang, and M.H. Mikkelsen.  Real-time tunable strong coupling: From individual nanocavities to metasurfaces.  ACS Photonics, 6(6):838–843, 2019.
  * [11] T. B. Hoang, G. M. Akselrod, and M. H. Mikkelsen.  Ultrafast room-temperature single photon emission from quantum dots coupled to plasmonic nanocavities.  Nano Letters, 16(2):270–275, 2016.
  * [12] J. J. Baumberg, J. Aizpurua, M. H. Mikkelsen, and D. R. Smith.  Extreme nanophotonics from ultrathin metallic gaps.  Nature Materials, 18(7):668–678, 2019.
  * [13] Hiroyuki Kishida and M.H. Mikkelsen.  Ultrafast lifetime and bright emission from graphene quantum dots using plasmonic nanogap cavities.  Nano Letters, 22(3):904–910, 2022.
  * [14] J. Fregoni, F.J. Garcia-Vidal, and J. Feist.  Theoretical challenges in polaritonic chemistry.  ACS Photonics, 9(4):1096–1107, 2022.
  * [15] M. Sanchez-Barquilla, A.I. Fernandez-Dominguez, J. Feist, and F.J. Garcia-Vidal.  A theoretical perspective on molecular polaritonics.  ACS Photonics, 9(6):1830–1841, 2022.
  * [16] Rong-Chun Ge and Stephen Hughes.  Design of an efficient single photon source from a metallic nanorod dimer: a quasi-normal mode finite-difference time-domain approach.  Opt. Lett., 39(14):4235–4238, Jul 2014.
  * [17] David A. Powell.  Resonant dynamics of arbitrarily shaped meta-atoms.  Phys. Rev. B, 90:075108, 2014.
  * [18] P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin.  Light interaction with photonic and plasmonic resonances.  Nanophotonics, 12(10):1700113, 2018.
  * [19] Wei Yan, Rémi Faggiani, and Philippe Lalanne.  Rigorous modal analysis of plasmonic nanoresonators.  Physical Review B, 97(20):205422, 2018.
  * [20] T. Wu, M. Gurioli, and P. Lalanne.  Nanoscale light confinement: the q‘s and v‘s.  ACS Photonics, 8(6):1522–1538, 2021.
  * [21] Christophe Sauvan, Tong Wu, Rachid Zarouf, Egor A. Muljarov, and Philippe Lalanne.  Normalization, orthogonality, and completeness of quasinormal modes of open systems: the case of electromagnetism.  Optics Express, 30(5):6846–6885, 2022.
  * [22] P. Fauche, S.K. Kosionis, and P Lalanne.  Collective scattering in hybrid nanostructures with many atomic oscillators coupled to an electromagnetic resonance.  Physical Review B, 95:195418, 2017.
  * [23] Rui-Qi Li, D. Hernangomez-Perez, F.J. Garcia-Vidal, and A.I. Fernandez-Dominguez.  Transformation optics approach to plasmon-exciton strong coupling in nanocavities.  Physical Review Letters, 117:107401, 2016.
  * [24] A. Cuartero-Gonzalez and A.I. Fernandez-Dominguez.  Light-forbidden transitions in plasmon-emitter interactions beyond the weak coupling regime.  ACS Photonics, 5(8):3415–3420, 2018.
  * [25] Rui-Qi Li, F.J. Garcia-Vidal, and A.I. Fernandez-Dominguez.  Plasmon-exciton coupling in symmetry-broken nanocavities.  ACS Photonics, 5(1):177–185, 2018.
  * [26] A. Cuartero-Gonzalez and A.I. Fernandez-Dominguez.  Dipolar and quadrupolar excitons coupled to a nanoparticle-on-mirror cavity.  Physical Review B, 101:035403, 2020.
  * [27] A. Cuartero-Gonzalez, A. Manjavacas, and A.I. Fernandez-DomÃnguez.  Distortion of the local density of states in a plasmonic cavity by a quantum emitter.  New Journal of Physics, 23:073011, 2021.
  * [28] I. Medina, F.J. Garcia-Vidal, A.I. Fernandez-Dominguez, and J. Feist.  Few-mode field quantization of arbitrary electromagnetic spectral densities.  Physical Review Letters, 126:093601, 2021.
  * [29] S. Franke, S. Hughes, M.K. Dezfouli, P.T. Kristensen, K. Busch, A. Knorr, and M. Richter.  Quantization of quasinormal modes for open cavities and plasmonic cavity quantum electrodynamics.  Physical Review Letters, 122:213901, 2019.
  * [30] J. Ren, S. Franke, A. Knorr, M. Richter, and S. Hughes.  Near-field to far-field transformations of optical quasinormal modes and efficient calculation of quantized quasinormal modes for open cavities and plasmonic resonators.  Physical Review B, 101:205402, 2020.
  * [31] M.K. Schmidt, R. Esteban, A. Gonzalez-Tudela, G. Giedke, and J. Aizpurua.  Quantum mechanical description of raman scattering from molecules in plasmonic cavities.  ACS Nano, 10(6):6291–6298, 2016.
  * [32] T. Neuman, R. Esteban, G. Giedke, M.K. Schmidt, and J. Aizpurua.  Quantum description of surface-enhanced resonant raman scattering within a hybrid-optomechanical model.  Physical Review A, 100:043422, 2019.
  * [33] Philippe Lalanne.  Light-in-complex-nanostructures/MAN: Versions 7.1 of QNMEig and QNMPole, 2020.
  * [34] F. Verstraete, M.M. Wolf, and J. Ignacio Cirac.  Quantum computation and quantum-state engineering driven by dissipation.  Nature Physics, 5:633–636, 2009.
  * [35] S. Franke, M. Richter, J. Ren, A. Knorr, and S. Hughes.  Quantized quasinormal-mode description of nonlinear cavity-qed effects from coupled resonators with a fano-like resonance.  Physical Review Research, 2:033456, 2020.
  * [36] C. Carlson, R. Salzwedel, M. Selig, A. Knorr, and S. Hughes.  Strong coupling regime and hybrid quasinormal modes from a single plasmonic resonator coupled to a transition metal dichalcogenide monolayer.  Physical Review B, 104:125424, 2021.
  * [37] S. Wuestner, A. Pusch, K.L. Tsakmakidis, J.M. Hamm, and O. Hess.  Gain and plasmon dynamics in active negative-index metamaterials.  Philosophical Transcations of the Royal Society A, 369(1950):3525–3550, 2011.
  * [38] A. Demetriadou, J.M. Hamm, Y. Luo, J.B. Pendry, J.J. Baumberg, and O. Hess.  Spatiotemporal dynamics and control of strong coupling in plasmonic nanocavities.  ACS Photonics, 4(10):2410–2418, 2017.
  * [39] N. Kongsuwan, A. Demetriadou, R. Chikkaraddy, L. Marrucci, J. J. Baumberg, and D. R. Smith.  Suppressed quenching and strong-coupling of purcell-enhanced single-molecule emission in plasmonic nanocavities.  ACS Photonics, 5(2):186–191, 2018.
  * [40] P. T. Kristensen, R. C. Ge, and S. Hughes.  Normalization of quasinormal modes in leaky optical cavities and plasmonic resonators.  Physical Review A, 92(5):053810, 2015.
  * [41] S. Simoncelli, E.M. Roller, P. Urban, R. Schreiber, A.J. Turberfield, T. Leidl, and T. Lohmuller.  Quantitative single-molecule surface-enhanced raman scattering by optothermal tuning of dna origami-assembled plasmonic nanoantennas.  ACS Nano, 10(11):9809–9815, 2016.
  * [42] R. Chikkaraddy, V. Turek, N. Kongsuwan, F. Benz, C. Carnegie, T. van de Goor, B. de Nijs, A. Demetriadou, O. Hess, and J. J. Baumberg.  Mapping nanoscale hotspots with single-molecule emitters assembled into plasmonic nanocavities using dna origami.  Nano Letters, 18(2):405–411, 2018.
  * [43] Anna Lombardi, Angela Demetriadou, Lee Weller, Patrick Andrae, Felix Benz, Rohit Chikkaraddy, Javier Aizpurua, and Jeremy J. Baumberg.  Anomalous Spectral Shift of Near- and Far-Field Plasmonic Resonances in Nanogaps.  ACS Photonics, 3(3):471–477, 2016.
  * [44] P. Alonso-Gonzalez, P. Albella, F. Neubrech, C. Huck, J. Chen, F. Golmar, F. Casanova, L. E. Hueso, A. Pucci, J. Aizpurua, and R. Hillenbrand.  Experimental verification of the spectral shift between near- and far-field peak intensities of plasmonic infrared nanoantennas.  Phys. Rev. Lett., 110(20):203902, 2013.
  * [45] R. Loudon.  The Quantum Theory of Light.  Oxford University Press, 1973.
  * [46] T. G. Eck.  Level crossings and anti-crossings.  Physica, 33(1):157–162, 1967.
  * [47] P. Törmä and W. L. Barnes.  Strong coupling between surface plasmon polaritons and emitters: a review.  Journal of Physics: Condensed Matter, 26(11):113201, 2014.
  * [48] L. C. Andreani, G. Panzarini, and J. M. Gérard.  Strong-coupling regime for quantum boxes in pillar microcavities: Theory.  Phys. Rev. B, 60(15):13276–13279, 1999.
  * [49] K. Santhosh, O. Bitton, L. Chuntonov, and G. Haran.  Vacuum rabi splitting in a plasmonic cavity at the single quantum emitter limit.  Nature Communications, 7(1):11823, 2016.
  * [50] K.-D. Park, M. A. May, H. Leng, J. Wang, J. A. Kropp, T. Gougousi, M. Pelton, and M. B. Raschke.  Tip-enhanced strong coupling spectroscopy, imaging, and control of a single quantum emitter.  Science Advances, 5(11):eaay5931, 2019.
  * [51] Z. Liu, J. Li, Z. Liu, W. Li, J. Li, C. Gu, and Z. Li.  Fano resonance rabi splitting of surface plasmons.  Scientific Reports, 7(1):13013, 2017.


0 Research please specify All authors are with Ohio State University. Hamoud
Alhazmi. E-mail<EMAIL_ADDRESS>Marcelo Morales. E-mail:
<EMAIL_ADDRESS>Jiachen Jiang. E-mail<EMAIL_ADDRESS>Jinxin Zhou.
E-mail<EMAIL_ADDRESS>Jian Chen. E-mail<EMAIL_ADDRESS>
# Analysis of the Causes of Car Accidents in the United States of America in
2023: Gauge People Understanding of Data Visualisation

Hamoud Alhazmi  Marcelo Morales  Jiachen Jiang  Jinxin Zhou  and Jian Chen

###### Abstract

This paper presents a comprehensive examination of interactive data
visualization tools and their efficacy in the context of United States car
accident data for the year 2023. We developed interactive heatmaps,
histograms, and pie charts to enhance the understanding of accident severity
distribution over time and location. Our research included the creation and
distribution of an online survey, consisting of nine questions designed to
test participants’ comprehension of the presented data. Fifteen respondents
were recruited to complete the survey, with the intent of assessing the
effectiveness of both static and interactive versions of each visualization
tool. The results indicated that participants using interactive heatmaps
showed a greater understanding of the data, as compared to those using
histograms and pie charts. In contrast, no notable difference in comprehension
was observed between users of static and interactive histograms. Unexpectedly,
static pie charts were found to be slightly more effective than their
interactive counterparts. These findings suggest that while interactive
visualizations can be powerful, their utility may vary depending on the type
and complexity of the data presented. Future research is recommended to
explore the influence of socioeconomic factors on the understanding of car
accident data, potentially leading to more tailored and effective
visualization strategies. This could provide deeper insights into the patterns
and causes of car accidents, facilitating better-informed decision-making for
stakeholders. Visit our website to explore our interactive plots and engage
directly with the data for a more comprehensive understanding of our findings.

  * •

URL: https://cse-5544-datavis-project.vercel.app/

###### keywords:

Visual Analytics, Traffic Accident Dataset, US-Accidents

Introduction and Motivation

Visualizing complex data remains crucial for several sectors, including public
safety, urban planning, and transportation management. Among the many types of
data, one that remains very important is the distribution of the severity of
the accident at different times of the day and in various geographic areas.
This understanding can offer great help to the decision process and therefore
in the implementation of interventions and safety measures. Many static and
interactive plots are used equally in the context of data visualization, each
with its advantages and disadvantages. This research examines these different
visualization techniques and compares them to determine which method best
provides the accident distribution information.

There have been recent advances in interactive visualization technology, of
which much has improved in terms of how a user interacts with data toward
dynamically exploring complex information with understanding. An interactive
plot allows the user to directly interact with the data. This involves
dragging views, drilling down and filtering out information. These actions
enable the user to gain deep insight from the data. This prompts a few
critical inquiries. Firstly, do these interactive features continue to hold
value? More essentially, does the interactivity truly augment our
comprehension of the distribution of accident severity, particularly when
juxtaposed with a conventional static plot?

This research is motivated by the need to:

  * •

Evaluate the Effectiveness of Interactive Visualization: By directly comparing
interactive and static visualizations, we aim to assess which method better
supports comprehension of complex datasets.

  * •

Identify Optimal Plot Types: This study aims to determine the most effective
form of interactive plots-HeatMap, Histogram, and PieChart-to communicate data
on the distribution of accident severity. The selection of the type of plot
can substantially affect the efficiency of information transmission and
comprehension.

  * •

Guide Future Visualization Tools: The results of this research will offer
crucial insights that can steer the creation of upcoming visualization
instruments, guaranteeing their ease of use and efficacy in conveying vital
data points.

By providing a comprehensive view of how visualization techniques can be
leveraged, this study offers practical insights into improving public safety
through improved data interpretation and decision processes. These findings
can be directly applied to the development of user-friendly visualization
tools.

##  1 Research Questions

The primary objective of this study is to examine the effectiveness of
different visualization techniques in helping people understand the
distribution of the severity of the accident over time and locations. To
achieve this, we have formulated the following specific research questions:

  * •

Is an interactive plot more effective than a static plot in helping people
understand the distribution of accident severity over time and locations?

  * •

Among interactive graphs such as HeatMap, Histogram, and PieChart, which is
the most effective to aid in understanding the distribution of the severity of
the accident?

The first question aims to determine whether the dynamic display of
interactive plots would yield better data understanding than that of static
plots with a fixed format. The second question will establish which type of
interactive visualization, HeatMap, Histogram, or PieChart, could best be used
to display complex spatial and temporal data. It should be used to understand
how the structure and presentation styles of the different plots facilitate
user comprehension and interaction.

##  2 Preliminary Hypotheses

Based on the research questions outlined, we propose the following hypotheses
to guide our experimental investigation:

  * •

Interactive plots will be more effective than static plots in helping people
understand the distribution of accident severity over time and locations.

  * •

Among the types of interactive plots, HeatMap will prove to be the most
effective in helping people understand the distribution of severity of
accidents.

Interactive plots allow users to use functionalities to engage deeply in the
data, e.g. zooming, filtering, selecting specific parameters, etc. This
capability improves comprehension and information retention since the facility
enables learners to learn at their own pace and style.

By nature, HeatMap offers an intuitive geographical context and is thus very
suitable for visualizing spatial data. It will be hypothesized that the
familiar format of maps, combined with interactive elements, offers faster and
more precise information on regional and temporal trends of accident severity.

Figure 1: A visualization of car accidents by regions to compare the severity
per region.

##  3 Related Work

In recent years, there has been a notable increase in the incidence of fatal
car accidents in all states. In the United States, there were 42,795 fatal
vehicle crashes in 2022 compared to 36,835 fatal car crashes in 2018, marking
an increase of more than 16% [5, 6]. The dataset we intend to utilize,
consisting of more than 2.25 million traffic accident reports collected for
the contiguous United States over a three-year period, represents the first
nationwide dataset of this magnitude.

Related work focuses on the prediction of real-time traffic accidents [4, 3].
However, our approach diverges as we aim to present the temporal distribution
of car accidents by showcasing crashes by time of day. Furthermore, we intend
to present various visual representations to effectively tackle the issue.
This includes utilizing a histogram to depict the severity of the crash. We
also plan to use a histogram showing the distribution of accidents over time
and severity levels. In addition, a HeatMap will be used to indicate accident
locations. Lastly, a gradient map will be used to illustrate the timing of the
most severe accidents, among others.

This project addresses the relationship between crash statistics in the United
States of America and how people perceive car accidents.

##  4 Evaluation Consideration

There are numerous visualization and analysis techniques that can be
considered to explore the US-Accidents dataset. However, we consider the
utilization of environmental factors, temporal analysis, comparative analysis,
and spatial analysis. The aforementioned analysis will provide us with a
complete picture of the problem.

Additionally, our intention is to present data on car accidents in each
contiguous state, followed by their classification according to regions, to
compare the severity within each region as illustrated in Figure 1. On the
other hand, to improve the performance, we think of creating multiple
interactive plots and combining them using Tableau. We will explore evaluation
methods similar to [2, 7, 8, 1].

##  5 Experimental Design

The experiment design aims to systematically compare the effectiveness of
interactive and static visualizations to help to understand the distribution
of the severity of the accident over time and locations. The design is
structured around two main components: figure design and survey design.

###  5.1 Hypothesis 1: Interactive vs. Static Maps

To effectively test the first hypothesis, that interactive maps are more
effective than static plots in helping people understand the distribution of
the severity of accidents over time and locations—the experiment design will
include a specific focus on comparing interactive and static maps. Here is a
step-by-step breakdown of how this will be implemented:

The interactive map will allow users to select different time segments
(morning, afternoon, evening, night) and regions (West, Midwest, South,
Northeast). This can be done through dropdown menus. Interactive elements such
as zoom, hover details (showing accident severity data upon hovering over a
region), and dynamic filtering (to isolate specific data points based on user
selection) will be included. Color gradients will be used to represent varying
levels of accident severity, with tooltips providing quantitative details.

For each of the timeslices, we will create a sequence of static maps to
distribute the severity of the accidents over the four areas of interest. Each
map will use an identical color gradient to the interactive one, only for
consistency. However, none of these maps will have interactive capabilities.
The maps will be presented in sequence or grid format, allowing the
participant to view differences between time and region of the forecast
without interaction with a view.

###  5.2 Hypothesis 2: Different Types of Plots

To test the second hypothesis -that graphs are the most effective type of
interactive graph to help people understand the distribution of accident
severity - our experimental design will compare interactive HeatMap,
Histogram, and PieChart. The detailed experimental setup is below for
evaluating the effectiveness of each visualization type.

HeatMap will provide geographical context with zoom, filter, and hover-over
details (such as specific accident severity data). Users can select different
times and regions to view specific data, enhancing spatial understanding.

The histogram will be designed to show the totals of the severity of the
accident for different regions and times. Users can click on specific bars to
drill down into more detailed temporal or regional data. Users can also use
the filter to select the range of total severity numbers of different states.

PieChart will depict accident severity data as a function of time of day, with
different pie segments representing different states. Clicking on a segment
can reveal more detailed data, and users can toggle between different layers
of information.

###  5.3 Survey Design

The study used a methodological approach involving 15 participants, segmented
into three groups of five. Each group was tasked with evaluating a different
data visualization format: one for heatmaps, another for histograms, and the
third for pie charts. This approach facilitated a targeted comparison of user
interaction and comprehension between the various visualization methods. After
interacting with the maps, participants will answer a series of survey
questions that assess their comprehension and ability to use the information
presented effectively.

These will include both direct questions about specific data points and more
general questions about trends and patterns. Tasks include identifying the
time of day or region with the highest severity of the accident or comparing
severity between different regions or times.

The responses will be evaluated for correctness. The time each participant
takes to complete the tasks will be recorded. This will help determine whether
interactive features enhance quicker comprehension compared to static visuals.

The data collected will be analyzed using statistical methods to compare
performance metrics (accuracy and response time) between all groups.
Effectiveness will be measured based on the speed and accuracy of responses,
providing quantitative evidence on whether interactive maps lead to better or
faster understanding compared to static maps.

##  6 Implementation Method

In this section, we detail the specific methods and technologies employed to
support our hypothesis. The strength of our execution approach is crucial to
guaranteeing the dependability and reproducibility of our findings.

###  6.1 Interactive Plots

Interactive plots are particularly effective in multivariate data analysis,
where dynamic adjustment of variables can reveal hidden relationships and
influences between data points. By incorporating interactive graphs into our
analysis, we provided a user-friendly interface using Tableau that allowed
reviewers and readers to explore data dependencies and correlations at their
discretion.

The website provides an intuitive and user-friendly interface, allowing
visitors to effortlessly navigate through various interactive plots and gain
insights without needing extensive technical knowledge. With its clean layout
and responsive design, the website ensures that users of all devices can
access and interact with the plots seamlessly, enhancing the user experience
across different platforms.

Each interactive plot on the website is accompanied by concise, informative
captions and contextual help, guiding users through the data exploration
process without overwhelming them. There are three interactive plots on the
website: HeatMap, Histogram, and PieChart as shown in Figures Figure 2, Figure
3, and Figure 4. The accidents and their severity can be viewed by users by
navigating through them and choosing the region and time.

Figure 2: Interactive visualization of HeatMap provided in our website. Figure
3: Interactive visualization of Histogram provided in our website. Figure 4:
Interactive visualization of PieChart provided in our website.

Visit our website to explore our interactive plots and engage directly with
the data for a more comprehensive understanding of our findings.

  * •

URL: https://cse-5544-datavis-project.vercel.app/

  * •

The source code is available at: https://github.com/mm587517/cse-5544-datavis-
project

###  6.2 Survey Study

The surveys were conducted using Google Forms, chosen for its accessibility
and ease of use. Google Forms facilitated the rapid deployment of our survey,
allowing us to collect responses in real time and immediately begin analyzing
trends and patterns in the data. We developed the survey based on the
hypothesis. In particular, the survey includes nine questions, as depicted in
Figure 5.

To assess understanding, the survey included multiple choice questions that
tested participants’ knowledge of key concepts, allowing us to measure their
grasp of the information presented. In addition, the survey used scenario-
based questions to simulate real-life situations, helping us to assess how
well participants understood and could apply information in practical
contexts. Our research team crafted a series of distinct surveys, three in
total, each aimed at measuring the level of understanding on specific topics
(i.e., HeatMap, Histogram, and PieChart) among diverse participant groups.

Figure 5: Survey questions to gauge participants’ perceptions.

##  7 Results and Analysis

Three distinct surveys were designed to evaluate the effectiveness of
different data visualization tools: heatmaps, histograms, and pie charts. Each
survey was administered to a separate group of five participants, specifically
selected to assess one type of visualization tool.

Consequently, a group of five participants evaluated the heatmap, another
group assessed the histogram, and the final group analyzed the pie chart. This
methodological approach ensured a focused and direct comparison of user
interaction and comprehension across different visualization formats.

###  7.1 Static vs. Interactive Charts

####  7.1.1 HeatMap

According to our analysis in comparing the effectiveness of interactive versus
static heatmaps, survey respondents have been observed to demonstrate
significantly higher proficiency with interactive models. More than 75% of the
participants accurately answered questions when using interactive heatmaps,
indicating a substantial improvement in understanding and participation, as
shown in Figure 6.

This increased performance can be attributed to the dynamic nature of
interactive heatmaps, which allows users to manipulate the data visually,
explore different layers of information, and gain insights through active
learning. Such interactive features are notably absent in static heatmaps,
where data presentation is fixed, limiting the viewer’s ability to delve
deeper into the specifics of the data. This clear disparity underscores the
importance of interactive elements in data visualization tools for improved
comprehension and analytical capabilities.

####  7.1.2 Histogram

Interestingly, in Figure 6, when the survey respondents were presented with
static and interactive histograms to analyze and answer the questions, the
results did not show significant differences in their performance. This
unexpected outcome suggests that the addition of interactive features in
histograms may not universally improve user understanding or engagement,
unlike other forms of data visualization such as heatmaps.

The similarity in the results between static and interactive histograms could
indicate that the simplicity and straightforward presentation of histograms,
which typically display distributions and frequencies, may be sufficient for
effective data comprehension. This finding challenges the assumption that
interactivity is always beneficial in data visualization tools and highlights
the need to consider the specific data type and user context when designing
these tools.

####  7.1.3 PieChart

Surprisingly, the results of the survey in Figure 6 revealed that respondents
performed slightly better when using static pie charts compared to their
interactive counterparts. This counterintuitive finding suggests that the
simplicity and directness of static pie charts may facilitate faster
comprehension and decision-making. Static pie charts provide a clear and
unchanging representation of the proportions of data, which might reduce the
cognitive load and streamline the process of data interpretation.

In contrast, interactive features in pie charts, while potentially enriching
the user experience with additional data layers or insights upon interaction,
could introduce complexities that momentarily distract or confuse users. This
suggests that the benefits of interactivity in data visualization are not
universally applicable and may depend on the context and nature of the data
presented.

Figure 6: The distribution of responses of Questions 1-9 (HeatMap, Histogram,
and PieChart).

###  7.2 Optimizing Data Insights with Comparison Charts

As suggested in the preliminary hypothesis section, heat maps are hypothesized
to enhance the comprehension of the accident severity distribution more
effectively than histograms and pie charts, as illustrated in Figure 7.
Quantitatively, the median values extracted from the surveys underscore this
assertion: the median severity rating observed in the heatmap is 72.5,
significantly higher than the 36.5 recorded for the histogram, and
substantially surpassing the 21.5 noted in the pie chart. This data supports
the hypothesis that heatmaps provide a superior visualization mechanism to
detect complex data distributions.

Figure 7: Comparative Analysis of Visualization Tool Efficacy: A Boxplot of
Interactive Heatmap, Histogram, and Pie Chart Scores.

##  8 Conclusion and Future Work

In the current study, we developed a series of interactive visual
representations - namely, heat maps, histograms, and pie charts - to analyze
vehicular accidents across the United States for the year 2023. These
interactive visual tools have proven to be instrumental in improving our
understanding of the complexities surrounding the distribution of the severity
of the accident in temporal and spatial dimensions. To evaluate the
effectiveness of these visual tools, we constructed an online survey
consisting of nine meticulously crafted questions aimed at gauging the
participants’ grasp of the severity distribution data. Fifteen respondents
were selected to participate in this survey.

The findings of the survey are illuminating. Interactive heat maps yielded the
most favorable results, suggesting a higher efficacy in conveying complex data
compared to histograms and pie charts. Interestingly, when comparing the
histograms, no significant differences in performance were discerned between
their interactive and static forms. In contrast, static pie charts
unexpectedly surpassed their interactive variants in terms of user
comprehension.

Taking these insights into account, future research might benefit from
examining socioeconomic variables, such as income, age, marital status,
occupation, and educational background, to enrich the analysis. An in-depth
understanding of these factors may prove pivotal in unraveling the
multifaceted nature of the severity of traffic accidents and its determinants.

## References

  * [1] R. Johansson, D. Williams, A. Berglund, and P. Nugues.  Carsim: a system to visualize written road accident reports as animated 3d scenes.  In Proceedings of the 2nd Workshop on Text Meaning and Interpretation, pp. 57–64, 2004.
  * [2] C. A. Kahn.  National highway traffic safety administration (nhtsa) notes.  Ann. Emerg. Med., 66:669–669, 2015.
  * [3] S. Moosavi, M. H. Samavatian, S. Parthasarathy, and R. Ramnath.  A countrywide traffic accident dataset.  arXiv preprint arXiv:1906.05409, 2019.
  * [4] S. Moosavi, M. H. Samavatian, S. Parthasarathy, R. Teodorescu, and R. Ramnath.  Accident risk prediction based on heterogeneous sparse data: New dataset and insights.  In Proceedings of the 27th ACM SIGSPATIAL international conference on advances in geographic information systems, pp. 33–42, 2019.
  * [5] D. Murfet.  United states. (2023) national highway traffic safety administration nhtsa.  https://crashstats.nhtsa.dot.gov/Api/Public/ViewPublication/813428.
  * [6] D. Murfet.  United states. (2023) national highway traffic safety administration nhtsa.  https://crashstats.nhtsa.dot.gov/Api/Public/ViewPublication/813482.
  * [7] J. Setiawan et al.  Visualization of road accident data in jabodetabek area year 2014-2017.  In Journal of Physics: Conference Series, vol. 1196, p. 012050. IOP Publishing, 2019.
  * [8] N. Shaadan, M. A. Suhaimi, M. Hazmir, and E. Hamzah.  Road accidents analytics with data visualization: a case study in shah alam malaysia.  In Journal of Physics: Conference Series, vol. 1988, p. 012043. IOP Publishing, 2021.