# Molecules in Real Cavities with Quantum Electroynamical Density Functional
Theory

Mark Kamper Svendsen<EMAIL_ADDRESS>Max Planck Institute
for the Structure and Dynamics of Matter and Center for Free-Electron Laser
Science & Department of Physics, Luruper Chaussee 149, 22761 Hamburg, Germany
CAMD, Department of Physics, Technical University of Denmark, 2800 Kgs.
Lyngby, Denmark  Center for Computational Quantum Physics, Flatiron Institute,
New York, New York 10010, USA  Kristian Sommer Thygesen  CAMD, Department of
Physics, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark  Angel
Rubio  Max Planck Institute for the Structure and Dynamics of Matter and
Center for Free-Electron Laser Science & Department of Physics, Luruper
Chaussee 149, 22761 Hamburg, Germany  Center for Computational Quantum
Physics, Flatiron Institute, New York, New York 10010, USA  Nano-Bio
Spectroscopy Group and European Theoretical Spectroscopy Facility (ETSF),
Universidad del País Vasco (UPV/EHU), Av. Tolosa 72, 20018 San Sebastian,
Spain  Johannes Flick  Center for Computational Quantum Physics, Flatiron
Institute, New York, New York 10010, USA  Department of Physics, City College
of New York, New York, New York 10031, USA  Department of Physics, The
Graduate Center, City University of New York, New York, New York 10016, USA

###### Abstract

Rapid experimental progress in realizing strongly coupled light-matter systems
in complex electromagnetic environments necessitates the development of
theoretical methods capable of treating light and matter from first
principles. A popular such method is quantum electrodynamical density
functional theory (QEDFT) which is a generalization of density functional
theory to situations where the electronic system is coupled to quantized light
modes. While this method provides a powerful description of the electronic
system and the quantized modes of light, it has so far been unable to deal
correctly with absorbing and dispersing electromagnetic media in practice. In
addition, the cavity field strength parameters have not been linked to the
real electromagnetic environment in which the matter is embedded meaning that
these are effectively free parameters. In this paper, we discuss how
macroscopic QED (MQED) can be invoked to correctly parameterize QEDFT for
realistic optical cavity setups. To exemplify this approach, we consider the
example of a absorbing spherical cavity and study the impact of different
parameters of both the environment and the electronic system on the transition
from weak-to-strong coupling. As a result of our work, the coupling parameters
in general, lossy environments can be now expressed in terms of the
_classical_ Dyadic Green’s Function. Because the Dyadic Green’s Function is
completely determined by the electromagnetic environment and the boundary
conditions, it thus removes the light-matter coupling strengths as free
parameters. As part of this work, we also provide an easy to use tool that can
calculate the cavity coupling strengths for simple cavity setups.

††preprint: APS/123-QED††preprint: APS/123-QED

##  I Introduction

Recent years have seen rapid developments in the experimental realization of
novel setups, where light and matter are strongly coupled. Examples include
plasmonic nanocavities [1], optical cavities using four-wave mixing schemes
[2], metasurface systems [3], self-assembled Casimir microcavities [4], deep-
strong coupling in plasmonic nanoparticles [5], and many more [6, 7]. Due to
the inherent complexity of these general strongly coupled light-matter
systems, theoretical methods usually apply various simplifications to keep
simulations tractable. One first-principles method which has gained popularity
is the recently introduced quantum-electrodynamical density-functional
theory(QEDFT) [8]. QEDFT is a generalization of density-functional theory
(DFT) [9] for electronic systems interacting strongly with quantized modes of
the electromagnetic field. The method represents a good compromise between
accuracy and computational cost, and it has been successfully applied to
desribe both the ground-state [10, 11], and excited-states [12, 13] of single
(few) molecules strongly coupled to quantized modes of light, as well as for
applications in polaritonic chemistry [14].

The existence of the QEDFT formulation can be proven under very general
conditions [8] and it in principle allows for the treatment of coupled light-
matter systems with many electrons and many photonic modes under very general
conditions. However, currently, most practical implementations of QEDFT are
based on the dipole approximation for the light-matter coupling, and a
discrete mode expansion of the electromagnetic field. The latter implies that
the material or molecular system of interest is embedded in a lossless
electromagnetic medium. There have been previous studies at applying the QEDFT
formulation to lossy optical cavities to describe e.g. photon losses through
cavity mirrors [12, 15, 16]. Here, different models of the optical cavity were
used, but no general connection between the cavity field parameters and the
optical environment for absorbing and dispersing magnetoelectric bodies has
yet been established. As a result, the electron-photon coupling parameters,
while in principle connected to the physical quantity of the vacuum electric
field at the center of charge of the system, are then in practice often
treated as free parameters. This highlights the other current practical
limitation of QEDFT: Treating the cavity coupling parameters as essentially
free parameters makes quantitative calculations and direct comparison with
experiments hard.

In the presence of absorbing and dispersing magnetoelectric bodies in the
electromagnetic environments, it is not possible to correctly quantize the
electromagnetic field via discrete modes. Instead, the field quantization can
be realized as a continuous expansion using the _classical_ dyadic Green’s
function (DGF) in the framework of macroscopic quantum electrodynamics (MQED)
[17]. It has previously been shown how standard DFT can be combined with MQED
to provide a quantitative, first principles description of the quantum light-
matter interactions for real cavity setups in the weak and intermediate
coupling regimes beyond the dipole approximation [18]. However, the previous
work relied on a wave function ansatz which only considered a subset of the
electronic structure. This is expected to become problematic in the ultra- and
deep strong coupling regimes.

In this work, we overcome the shortcomings of the existing methodologies by
establishing a direct connection between the frameworks of MQED and QEDFT that
can be used to study the interaction of molecular systems with arbitrary
quantized electromagnetic fields while correctly considering absorbing and
dispersing electromagnetic environments. Specifically, we investigate how
QEDFT can be formulated in terms of the field expansion from MQED within the
dipole approximation for the light matter coupling. This in turn allows us to
express the cavity coupling parameters in terms of the _classical_ DGF, which
can be related to a real cavity setup. The result is that we arrive at a
formulation of QEDFT appropriate in lossy electromagnetic environments which
simultaneously relates the coupling parameters to a real cavity. This removes
the electron-photon coupling strengths as free parameters. We will exemplify
this scheme on a lossy spherical microcavity setup and study the impact of
different parameters of both the environment and the electronic system on the
transition from weak-to-strong coupling.

##  II Emitter-centered formulation of Macroscopic Quantum Electrodynamics

Figure 1: (a) Different radii of the spherical microcavity that host a
resonance that aligns with the $\Pi\rightarrow\Pi^{*}$ of the benzene
molecule. The total light-matter coupling strength of the modes is also shown
in the legend, and it can be clearly observed that the cavity-coupling
strength grows with decreasing cavity radius. The mode structure is shown with
a sampling density of $10$ points/meV. (b) The linear absorption spectrum of
benzene as a function of cavity radius showing a clear radius dependent
Purcell enhancement.

Macroscopic Quantum Electrodynamics (MQED) is a framework for quantizing the
electromagnetic field in the presence of arbitrary absorbing or dispersing
environments [17, 19, 18]. The central object in MQED is the classical Green’s
function that solves the Helmholtz equation for a point source, the so-called
dyadic Green’s function (DGF)[20],

$\displaystyle\left[\bm{\nabla}\times\kappa(\bm{r},\omega)\bm{\nabla}\times-\frac{\omega^{2}}{c^{2}}\epsilon(\bm{r},\omega)\right]\bm{G}(\bm{r},\bm{r^{\prime}},\omega)=\delta(\bm{r}-\bm{r^{\prime}}),$
(1)

where $\kappa(\bm{r},\omega)=\mu^{-1}(\bm{r},\omega)$ and
$\epsilon(\bm{r},\omega)$ are the spatially dependent dielectric function and
inverse magnetic permability respectively. The DGF is of central importance to
the quantized theory of electromagnetic fields in lossy environments because
it simultaneously carries the information about the electromagnetic boundary
conditions, and serves as a projector from the coupled light-matter system
onto the electromagnetic degrees of freedom [17].

For a spatially local magnetoelectric medium in the nonrelativistic limit, the
MQED expansion of the electric field in the Power-Zienau-Woolley (PZW) frame
[21, 22] (multipolar gauge) can be written as [17, 19],

$\displaystyle\bm{\hat{E}}(\bm{r})$ $\displaystyle=\int
d\omega\bm{\hat{E}}(\bm{r},\omega)+\mathrm{h.c.},$ (2)
$\displaystyle\bm{\hat{E}}(\bm{r},\omega)$
$\displaystyle=\sum_{\lambda=e,m}\int
d^{3}r^{\prime}\bm{G}_{\lambda}(\bm{r},\bm{r}^{\prime},\omega)\cdot\bm{\hat{f}}_{\lambda}(\bm{r}^{\prime},\omega).$
(3)

Here $\bm{\hat{f}}_{\lambda}(\bm{r}^{\prime},\omega)$ are the spatially
resolved polaritonic field operators of MQED which fulfill the commutation
relations of the quantum harmonic oscillator, and
$\bm{G}_{e}(\bm{r},\bm{r}^{\prime},\omega)$ and
$\bm{G}_{m}(\bm{r},\bm{r}^{\prime},\omega)$ are the electric and magnetic
components of the DGF respectively,

$\displaystyle\bm{G}_{e}(\bm{r},\bm{r}^{\prime},\omega)=i\frac{\omega^{2}}{c^{2}}\sqrt{\frac{\hbar}{\pi\epsilon_{0}}\mathrm{Im}\epsilon(\bm{r},\omega)}\bm{G}(\bm{r},\bm{r}^{\prime},\omega),$
(4)
$\displaystyle\bm{G}_{m}(\bm{r},\bm{r}^{\prime},\omega)=-i\frac{\omega}{c}\sqrt{\frac{\hbar}{\pi\epsilon_{0}}\mathrm{Im}\mu(\bm{r},\omega)}\bm{G}(\bm{r},\bm{r}^{\prime},\omega)\times\bm{\nabla^{\prime}}.$
(5)

In the following we neglect magnetic interactions and consider the coupling
between light and matter within the dipole approximation. Therefore, if we
consider a set of emitters $i$ with positions (centers of charge)
$\bm{r}_{i}$, the interaction only samples the electromagnetic field at these
positions. In this sense, the full electric field $\hat{\bm{E}}(\bm{r})$ in
Eq. 2 contains more information than strictly necessary to describe the light-
matter interaction completely. As discussed in Refs [19, 23, 24], it is
therefore possible to arrive at a significantly more compact expression by
alternatively expanding the electric field in terms of a set of $\mathcal{N}$
explicitly orthogonalized bright modes at each frequency,

$\displaystyle\bm{\hat{E}}(\bm{r})=\sum_{i=1}^{\mathcal{N}}\int_{0}^{\infty}d\omega\bm{E}_{i}(\bm{r},\omega)\hat{C}_{i}(\omega)+\mathrm{H.c.},$
(6)
$\displaystyle\bm{E}_{i}(\bm{r},\omega)=\frac{\hbar\omega^{2}}{\pi\epsilon_{0}c^{2}}\sum_{j=1}^{\mathcal{N}}V^{*}_{ij}(\omega)\frac{\mathrm{Im}\bm{G}(\bm{r},\bm{r}_{j},\omega)\cdot\hat{n}_{j}}{G_{j}(\omega)}.$
(7)

Here $\hat{C}^{(\dagger)}_{i}(\omega)$ destroys (creates) a photon in the
$i$th bright mode. $\bm{E}_{i}(\bm{r},\omega)$ describes the spatial mode
function of the electric field associated with mode $i$. The normalisation
factor $G_{j}(\omega)$ is the square root of the dipole spectral density[25],

$\displaystyle
G_{j}(\omega)=\left(\frac{\hbar\omega^{2}}{\pi\epsilon_{0}c^{2}}\bm{n}_{j}\cdot\mathrm{Im}\bm{G}(\bm{r}_{j},\bm{r}_{j},\omega)\cdot\bm{n}_{j}\right)^{1/2}.$
(8)

Finally, the matrix $V_{ij}(\omega)$ is the transformation matrix which obeys
$\bm{V}(\omega)\bm{S}(\omega)\bm{V}^{\dagger}(\omega)=\bm{I}$ where

$\displaystyle
S_{ij}(\omega)=\frac{\hbar\omega^{2}}{\pi\epsilon_{0}c^{2}}\frac{\bm{n}_{i}\cdot\mathrm{Im}\bm{G}(\bm{r}_{i},\bm{r}_{j},\omega)\cdot\bm{n}_{j}}{G_{i}(\omega)G_{j}(\omega)},$
(9)

and it is a result of the mode-orthogonalization inherent to the emitter-
centered representation [19]. The number of emitter centered modes per
frequency, $\mathcal{N}$, is equal to the number of emitter positions times
the number of dipole orientations considered. In this work, we consider a
single emitter position $\bm{r}_{0}$ and the full three-dimensional space of
dipole orientations. This results in three emitter-centered modes per
frequency.

We want to express the total Hamiltonian of the coupled light-matter system in
a similar form as used previously in Refs. [26, 12, 15]. Therefore we write
(see Supplementary note A),

$\displaystyle\mathcal{H}=$ $\displaystyle\mathcal{H}_{\mathrm{Mat}}+$
$\displaystyle\frac{1}{2}\sum_{i=1}^{\mathcal{N}}\int_{0}^{\infty}d\omega\left\\{\hat{p}_{i}(\omega)^{2}+\omega^{2}\left[\hat{q}_{i}(\omega)+\frac{\bm{\lambda}_{i}(\omega)}{\omega}\cdot{\hat{\bm{R}}}\right]^{2}\right\\}.$
(10)

The light-matter interaction now contains the dipole moment operator for
$N_{e}$ electrons with position $\hat{\bm{r}}_{i}$,
${\hat{\bm{R}}=\sum_{i=1}^{N_{e}}\hat{\bm{r}}_{i}}$, which gives rise to an
explicit electron-photon interaction and the dipole-self energy term.
$\mathcal{H}_{\mathrm{Mat}}$ is the standard matter Hamiltonian that describes
the electronic system. The cavity field strength $\bm{\lambda}_{i}(\omega)$ is
given by
$\bm{\lambda}_{i}(\omega)=-e\left(\frac{2}{\hbar\omega}\right)^{1/2}\bm{E}_{i}(\bm{r}_{0},\omega)$,
where $\bm{r}_{0}$ denotes the center of charge. We have further introduced
the photon field quantities $\hat{p}_{i}$ and $\hat{q}_{i}$ that are connected
to the magnetic- and electric field in their corresponding mode and are given
explicitly by
$\hat{q}_{i}(\omega)=\left(\frac{\hbar}{2\omega}\right)^{1/2}(\hat{C}_{i}(\omega)+\hat{C}^{\dagger}_{i}(\omega))$
and
$\hat{p}_{i}(\omega)=\left(\frac{\hbar\omega}{2}\right)^{1/2}(\hat{C}_{i}(\omega)-\hat{C}^{\dagger}_{i}(\omega))$.
In terms of these new quantities, the electric field expansion at the center
of charge $\bm{r}_{0}$ reads,

$\displaystyle\bm{\hat{E}}(\bm{r}_{0})=-\sum_{i=1}^{\mathcal{N}}\int_{0}^{\infty}{\omega}\bm{\lambda}_{i}(\omega)\hat{q}_{i}(\omega)d\omega.$
(11)

The interaction of the electronic system with the electric field within the
dipole approximation can thus be expressed as the interaction of the
electronic system and a continuous set of quantum harmonic oscillator modes
via the dipole moment of the electronic system. Note that this formulation
addresses both the problem of how to formulate the length gauge Hamiltonian in
the presence of losses, and allows for the explicit calculation of the cavity
field strengths in terms of the boundary conditions set by the cavity via the
DGF of the electromagnetic field.

##  III Quantum electrodynamic density functional theory

To calculate the excited-state properties of the coupled electron-photon
problem defined by the Hamiltonian in Eq. II, we employ the linear response
formulation of QEDFT, which has been described in detail in Refs. [12, 27,
13]. Importantly, this method considers the _full_ electronic structure of the
matter system and goes beyond the rotating wave approximation. By combining it
with the MQED quantization of the electromagnetic field we thus go beyond the
method presented in Ref. [18]. Specifically, we employ the generalized Casida
formulation of linear response QEDFT to calculate the oscillator strengths of
the many-body excitations of the coupled light-matter (polaritonic) system.
The generalized Casida equation reads,

$\displaystyle\begin{bmatrix}U&V\\\
V^{T}&\omega_{\alpha}\end{bmatrix}\begin{bmatrix}\bm{F}\\\
\bm{P}\end{bmatrix}=\Omega^{2}_{S}\begin{bmatrix}\bm{F}\\\
\bm{P}\end{bmatrix}.$ (12)

Here $\omega_{\alpha}$ is a diagonal matrix with the frequency of photon modes
in the diagonal. $U$ accounts for the coupling between the electrons.
Introducing the pair orbital index $S=(ia)$, corresponding to the pair orbital
$\Phi_{S}(\bm{r})=\phi_{i}(\bm{r})\phi_{a}^{*}(\bm{r})$ with energy
$\epsilon_{S}=\epsilon_{a}-\epsilon_{i}$, $U$ can be expanded as,

$\displaystyle
U_{SS^{\prime}}=\epsilon_{S}^{2}\delta_{SS^{\prime}}+2\epsilon_{S}^{1/2}K_{SS^{\prime}}(\Omega_{S})\epsilon_{S^{\prime}}^{1/2},$
(13) $\displaystyle K_{SS^{\prime}}(\Omega_{S})=\int\int
d\bm{r}d\bm{r^{\prime}}\Phi_{S}(\bm{r})f^{n}_{\mathrm{Mxc}}(\bm{r},\bm{r^{\prime}},\Omega_{S})\Phi_{S^{\prime}}(\bm{r^{\prime}}).$
(14)

If $N_{\mathrm{pair}}$ pair-orbitals are included in the calculation, $U$ is
therefore an $N_{\mathrm{pair}}\times N_{\mathrm{pair}}$ matrix. $V_{\alpha
S}$ is the matrix accounting for the coupling between the electrons and the
photon modes,

$\displaystyle V_{\alpha S}=2\sqrt{\epsilon_{S}M_{\alpha
S}(\Omega_{S})N_{\alpha S}\omega_{\alpha}},$ (15) $\displaystyle
M_{\alpha,S}(\Omega_{S})=\int
d\bm{r}\Phi_{S}(\bm{r})f_{\mathrm{Mxc}}^{q_{\alpha}}(\bm{r},\Omega_{S}),$ (16)
$\displaystyle N_{\alpha,S}=\frac{1}{2\omega_{\alpha}^{2}}\int
d\bm{r}\Phi_{S}(\bm{r})g_{\mathrm{M}}^{q_{\alpha}}(\bm{r}).$ (17)

The size of these matrices will be $N_{\mathrm{photon}}\times
N_{\mathrm{pair}}$. The $\bm{F}$ part of the eigenvector in Eq. 12 can be
interpreted as the matter part and the $\bm{P}$ can be interpreted as the
photon part. The square norm of the two gives the electronic- and photonic
fraction of the excitation respectively. As discussed in Ref. [12] it is
possible to calculate the oscillator strengths of the coupled system in terms
of these eigenvectors. Finally, note that when neglecting the quantized light,
only the top left block in Eq. 12 survives and one recovers the standard
Casida formulation of linear response TDDFT for finite systems [28]. We note
that these equations include the exchange correlation kernels $f_{Mxc}^{n}$,
$f_{Mxc}^{q}$ and $g_{\mathrm{M}}^{q_{\alpha}}$ which have to be approximated
in practice. So far the available exchange-correlation functionals for QEDFT
in the time-domain are still limited [10, 11], but recent developments based
on the photon-free formulation of QEDFT [29] or the QEDFT fluctuation-
dissipation theorem [30] are promising. In this work, we neglect any exchange-
correlation contribution in the kernels and apply the mean-field photonic
random phase approximation. We refer to Ref. [12] for a thorough discussion of
this approximation.

The excitation spectrum is then characterized by the strength function [12],

$\displaystyle S(\omega)=\sum_{I}f_{I}\delta(\omega-\omega_{I}),$ (18)

where $f_{I}$ is the oscillator strengths, and $I$ runs over the many-body
excitations of the coupled light-matter (polaritonic) system. We emphasize
that in this work we do not apply any broadening to the spectra and that the
width of the peaks is determined completely by the electromagnetic
environment. This highlights a further advantage of using the QEDFT-MQED
combination, namely that this provides a natural description of the transition
linewidths as they relate to decay induced by the electromagnetic environment.
This is relevant beyond cavity QED settings and paves the way towards TDDFT
without artificial linewidths.

##  IV Spherical microcavity

Figure 2: (a) The linear absorption spectrum of benzene in the gold cavity
from Fig. 1a with an inner radius $R=16$nm as a function of the Drude
dampening in the metallic mirror region. (b) The linear absorption spectrum in
a cavity with with an inner radius $R=16$nm and Drude dampening
$\gamma=\gamma_{\mathrm{gold}}/4$ as a function of the number of benzene
molecules.

To exemplify the developed approach, we consider the spherically layered
microcavity also considered in Ref. [31]. As shown in the inset of Fig. 1a,
the spherical cavity consists of three concentric spherical layers, each
characterized by a frequency dependent dielectric function,
$\epsilon_{n}(\omega)$. In general, the source and field points ($\bm{r}$ and
$\bm{r}^{\prime}$, respectively) can be in either the same layer or different
layers. Consequently, the DGF for the reflected field is a nine-component
object, where each of the components is a three-dimensional dyad. Labelling
the possible combinations of source and field points by two extra indices
$m,n=1,2,3$, the full DGF of the system therefore takes the following form,

$\displaystyle\bm{G}_{mn}(\bm{r},\bm{r^{\prime}},\omega)=\bm{G}^{\mathrm{vac}}(\bm{r},\bm{r}^{\prime},\omega)\delta_{nm}+\bm{G}^{\mathrm{ref}}_{mn}(\bm{r},\bm{r}^{\prime},\omega),$
(19)

where $\bm{G}^{\mathrm{vac}}$ and $\bm{G}^{\mathrm{ref}}$ are the components
of the DGF describing the free dipole field and the reflections of the dipole
field in the surroundings, respectively. In combination, the two fully
describe the electromagnetic environment of the cavity. Due to the spherical
symmetry of the problem, the DGF of the system is most efficiently expanded
onto vector spherical harmonics [20]. The nine different components of the
reflection contribution can be calculated by explicitly invoking the
electromagnetic boundary conditions at the material interfaces as appropriate.
In this work we consider the situation where the emitter is placed in the
inner region of the cavity, and it is therefore only necessary to consider the
$n=m=3$ component of the DGF.

For a general emitter position inside the cavity, it is necessary to carefully
converge the number of vector spherical harmonics used in the calculation of
the DGF. However, if the emitter is placed exactly in the center of the
cavity, the situation simplifies significantly. In this case, only the lowest
order magnetic mode contributes and the reflection DGF can be written as [31],

$\displaystyle\bm{G}^{\mathrm{ref}}_{33}(\bm{r},\bm{r^{\prime}},\omega)|_{r,r^{\prime}\rightarrow
0}=\frac{i\omega}{6\pi c}r^{N}_{n=1}(\omega)\bm{\mathbbm{1}}.$ (20)

The imaginary part of the spatial trace of the vacuum DGF is given by
$\mathrm{Im}\bm{G}^{\mathrm{vac}}(\bm{r},\bm{r},\omega)=\frac{\omega}{6\pi
c}\bm{\mathbbm{1}}$. The imaginary part of the full DGF, evaluated in the
center of the cavity $\bm{0}$, thus becomes,

$\displaystyle\mathrm{Im}\bm{G}(\bm{0},\bm{0},\omega)=\frac{\omega}{6\pi
c}\left[1+\mathrm{Re}\left(r^{N}_{n=1}(\omega)\right)\right]\bm{\mathbbm{1}}.$
(21)

Notice that
$\bm{n}_{i}\cdot\mathrm{Im}\bm{G}(\bm{0},\bm{0},\omega)\cdot\bm{n}_{j}\propto\delta_{ij}$
which means that the mode orthogonalization is trivial in this case, and the
cavity field strengths can be derived directly,

$\displaystyle\bm{\lambda}_{i}(\omega)=e\left(\frac{\omega^{2}}{3\pi^{2}\epsilon_{0}c^{3}}\left[1+\mathrm{Re}\left(r^{N}_{n=1}(\omega)\right)\right]\right)^{1/2}\hat{n}_{i}.$
(22)

In practice, to use the QEDFT Casida formalism, it is necessary to employ a
dense, discretized sampling of the continuous frequency expressions for the
coupling strengths. This procedure is discussed in Supplementary note C.
Furthermore, the details of the QEDFT calculations with these modes are given
in Supplementary note D.

###  IV.1 Drude metal shell

We now consider the case where the inner and outer regions consist of vacuum.
For the middle region we consider a simple, but realistic model of a metallic
mirror, namely a Drude metal with dielectric function,

$\displaystyle\epsilon_{2}(\omega)=1-\frac{\omega_{p}^{2}}{\omega^{2}+i\gamma\omega},$
(23)

where $\omega_{p}$ is the metal plasma frequency and $\gamma$ is the Drude
dampening rate. As a concrete example of a metal, we use the Drude parameters
for gold taken from Ref. [32]. Within the Drude mode, the plasma frequency of
gold is around 8.5 eV 111It should be noted that the Drude model is not a good
model for the dielectric properties of gold due to the neglect of d-band
absorption. Additionally, it should be mentioned that the Drude model will
fail to account for size dependent resonance shifts and broadening resulting
from the nonlocal response and surface enhanced Landau dampening in small
metallic structures [50, 51].. This results in the dielectric function shown
in Supplementary Figure. 1b. Below the plasma frequency, the real part of the
dielectric function will be negative and the material surface will
consequently be highly reflective. Above the plasma frequency, the real part
of the dielectric function becomes positive and the material will loose its
metallic characteristics resulting in a significant loss of surface
reflectivity.

In Supplementary Figure 1c the mode structure of the cavity setup is shown for
two different cavity radii, $R=140$nm and $R=450$nm. It is clearly observed
that the number of modes, as well as their spectral position, is directly
linked to the radius of the microcavity. Furthermore, we clearly observe that
above the gold plasma frequency, the mirrors lose their reflectivity which
results in the loss of the sharp mode structure which is replaced by a
continuum. This highlights that the formalism we present is able to directly
link the cavity field strengths to the real cavity setup made of real
materials. Supplementary Figure 1d zooms in on the mode around 7.1 eV in the
cavity with $R=140$nm and shows the effect of changing the Drude dampening
parameter. We clearly observe that the width of the cavity mode increases with
increasing dampening in the metal which highlights the connection between the
width of the cavity modes and the losses in the gold.

This example highlights how the use of the emitter centered representation of
MQED allows us to directly and uniquely relate the light-matter coupling
strength to a real electromagnetic environment and connect it to the QEDFT
formalism. While this is a relatively simple example, the approach is general
and works analogously for an arbitrary electromagnetic environment provided
that the DGF can be determined.

###  IV.2 Adding an emitter to the cavity

We next add a benzene molecule to the cavity. Benzene is chosen mainly because
of it prevalence as a test system in the existing TDDFT and QEDFT literature
on strong coupling [12, 15, 34, 11], but we emphasize that the method can
treat arbitrary finite electronic systems. We focus on finding a cavity
configurations with a mode resonant with the $\Pi\rightarrow\Pi^{*}$
transition of the benzene molecule. This transition is spectrally well
isolated from other transitions. It thus represents an effective two level
system for coupling strengths in the weak coupling and strong coupling regime,
while care has to be taken in the ultrastrong coupling regime. The first step
is to determine the spectral position of this transition. Using the Casida
linear response QEDFT framework without photons, we find that the transitions
occurs at 6.808 eV (182 nm) in free space, consistent with previous TDDFT
calculations for benzene [12, 15].

As shown in Fig. 1a, it is possible to find different radii of the gold
microcavity for which there is a cavity mode resonant with the benzene
$\Pi\rightarrow\Pi^{*}$ electronic transition. As expected, the cavity field
strength increases as the cavity is made smaller. All but the smallest cavity
are optical cavities in the sense that the characteristic dimension of the
cavity, the radius, is larger than half the wavelength of the transition. For
the smallest cavity of radius 16 nm, a significant increase in the coupling
strength relative to the other sizes is observed. This happens exactly because
this cavity is sub-wavelength sized and therefore significant near-field
coupling to the surface plasmon mode of the gold starts to occur.
Fig. 1b shows the linear absorption spectra of the coupled emitter-cavity
system calculated for the different cavity radii using the linear response
QEDFT method. We emphasize that all linewidths in the figure are true
linewidths in the sense that they are not related to any broadening parameters
in the QEDFT calculation and only reflect the density of states in the optical
environment. A radius-dependent Purcell enhancement with decreasing radius is
clearly observed, reflecting the reduction in radiative lifetime resulting
from the altered optical environment. However, it is not possible to achieve
strong coupling with a single benzene molecule using the gold-shell cavity. We
attribute this to the fact that as the coupling strength gets larger with
decreasing radius, the optical losses also increase, resulting in a broader
cavity resonance. We note in passing that the Purcell enhancement for the
smallest $R=16$ nm cavity is around 3500 which means that the local field
enhancement at the center of the spherical microcavity is comparable to what
is found in experiments with plasmonic microcavities [35].

The reason that it was impossible to reach strong coupling with benzene in the
gold cavity was the losses of the cavity mode. In an attempt to reach the SC
regime, we therefore next seek to reduce the losses in the cavity. As already
discussed above, the width of the cavity resonance is reduced for smaller
Drude dampening parameters, $\gamma$. For this reason, Fig. 2a shows the
absorption spectrum for the case with the true gold dampening, as well as 25%,
10% and 5% of the dampening respectively. We mention in passing that in
practice one could imagine realizing these lower losses for example using
metals specifically engineered to show weaker losses [36]. Furthermore, one
could also imagine exploring different cavity setups potentially leveraging
the lower losses in dielectric nano-optical setups [37]. Considering the wide
range of available materials this design space becomes enormous [38, 39]. As
shown in the inset of Fig. 2a, we find that reducing the losses results in a
narrower cavity mode without a significant reduction in the overall cavity
field strength. At around 25% of the true dampening, we begin to observe clear
indications of the two polariton peaks in the linear absorption spectrum.
Further reducing the dampening we see clear strong coupling with a Rabi
splitting of around $2g\sim 10\mathrm{\,eV}$. This emphasizes the importance
of the optical losses in reaching the strong coupling regime, and further
highlights the significant strength of our method that we are able to study
the effect of different cavity parameters from first principles via our
combination of QEDFT and MQED.

Another way to engineer the coupling strength is by changing the emitter.
There are two ways to do this, either by changing the number of emitters or by
changing the emitter itself. To investigate the first option, we take the
cavity with 25% of the true gold losses, and compute the absorption spectrum
for one, two and three benzene molecules, all of which we place in the center
of the cavity. As shown in Fig. 2b, we see a clear evolution of the Rabi
splitting with the number of benzene molecules indicating the onset of
collective strong coupling [40, 41, 42, 43]. To perform this analysis, one
needs to solve a coupled many electron, many photon problem and it again
highlights the strength of the method.

Figure 3: (a) The Casida spectrum without photon modes for different aromatic
compounds. The spectra here are shown with an artificial broadening of 0.1361
eV because the photon-free calculations fail to naturally describe the
linewidth of the transitions. (b) Transition dipole moment of the
$\Pi\rightarrow\Pi^{*}$-transition as a number of aromatic rings. (c) Mode
structure of the smallest cavities which host resonances aligned with the
$\Pi\rightarrow\Pi^{*}$-transition of the different molecules. The mode
structure is shown with a sampling density of $10$ points/meV. (d) Effective
light-matter coupling strength as a function of the number of benzene rings in
the aromatic molecule.

To investigate the latter option, we consider longer aromatic compounds with
$N$ aromatic rings; naphthalene $(N=2)$, anthracene $(N=3)$, tetracene $(N=4)$
and pentacene $(N=5)$. We first use the Casida method without photons
(standard TDDFT) to characterize the spectral properties of the aromatic
molecules. Fig. 3(a,b) shows respectively the spectrum and transition dipole
moment of the $\Pi\rightarrow\Pi^{*}$-transition (HOMO-LUMO) as a function of
the number of aromatic rings $N$. Note that in Fig. 3(a) the spectra are shown
with the Octopus default artificial broadening of 0.1361 eV. This artificial
broadening is necessary because unlike the combined MQED-QEDFT method, the
standard photon-free TDDFT formulation fails to describe the linewidth of the
transitions. We observe that with increasing number of aromatic rings $N$, the
transition energy redshifts and the transition dipole moment increases
linearly (Fig. 3(b)). The linearly increasing transition dipole moment would
suggest that the light-matter coupling strength can be monotonically increased
by simply using a larger aromatic molecule. However, because the transition
energy also redshifts with increasing molecule length, the cavity has to be
re-optimized to be resonant with the transition for each molecule, as shown in
Fig. 3 (c). Specifically, focusing on modifications of the $R=16$ nm cavity,
we find that this re-optimization of the cavity means that the cavity radius
needs to be increased. This increase in radius leads to a reduced field
concentration via an increased effective mode volume. This behavior highlights
the important point that the light-matter coupling strength is a joint
property of both the electronic system and the electromagnetic environment. A
proper treatment of both is therefore essential for quantitative predictions.

We can characterize the intricate interplay between the electromagnetic
environment and the electronic structure by looking at an effective coupling
strength for the cavity modes which we define as [15],

$\displaystyle g_{\mathrm{eff}}\propto\sqrt{\omega_{c}}\lambda_{c}|\bm{d}|.$
(24)

Here the total cavity field strength parameter is defined as the coupling
strength averaged across the cavity peak $p$,
$\lambda_{c}=\sqrt{\int_{p}d\omega|\bm{\lambda}_{d}(\omega)|^{2}}$, and
$\omega_{c}$ is the center frequency of the cavity mode. The $d$ subscript
indicates that we take the field strength parameter for dipole orientation
$\bm{d}/|\bm{d}|$. $g_{\mathrm{eff}}$ would thus be the true light-matter
coupling if the total spectral weight was concentrated in a single mode. As
shown in Fig. 3d, we observe that the increase in transition dipole moment is
counteracted by the reduced field concentration for the larger molecules,
effectively resulting in a weaker light-matter coupling strength. This is in
stark contrast to the intuitive argument based solely on the increased dipole
moment of the longer molecules. It should be noted that $g_{\mathrm{eff}}$ is
not a perfect measure of the light-matter coupling strength and it only gives
a rough idea of the cavity re-optimization’s effect. This is because the
widths of the modes are not taken into account which, as we have seen in Fig.
2a, is very important for the nature of the light-matter coupling.

We mention in passing that with the presented framework it would be possible
to perform further engineering of the electromagnetic environment to increase
the coupling for the larger aromatic compounds. The application of the
framework to general electromagnetic environments is discussed further below.

##  V Comment on Fabry Perot cavities

A common example of a cavity in the literature for both theory and experiments
is the layered Fabry Perot cavity (FPC). As discussed in Supplementary note E,
the FPC is a layered system and consequently its DGF is expanded in terms of
in-plane plane waves, augmented by a function accounting for the reflection at
the interfaces between the layers [20]. Because the FPC only constrains the
electromagnetic modes in one direction, it retains significant dispersion of
the modes in plane. Consequently, the concentration of electromagnetic density
of states is significantly less efficient than in the case of e.g. the
spherical microcavity. As discussed in Supplementary note E, this means that
the resulting coupling strength is weaker and the FPC will therefore generally
not be suited for single- or few emitter strong coupling [44]. For this
reason, we do not perform explicit QEDFT calculations for this cavity setup.
However, the FPC can still be suited for collective strong coupling [45] and
coupling to extended systems [46] where the extended modes of the
electromagnetic environment can be sampled more effectively.

##  VI Tool for cavity field parameters in simple cavities

As a part of this work, we are making the code to generate the cavity field
strengths available for everyone to use as part of the new PhotonPilot tool.
This tool currently allows the user to calculate cavity field strengths for
the spherical and layered cavity setups, and we plan to expand its
capabilities in the future.

##  VII General electromagnetic environments

We emphasize that the method we have presented here is general and applicable
to any electromagnetic environment as long as the DGF can be determined.
However, in general electromagnetic environments with lower symmetry it is not
possible to write down an analytical expression for the DGF. In such cases,
the DGF must be constructed numerically from e.g. a mode expansion based on
finite element simulations [20, 25]. We note in passing that in the general
setting the Helmholtz equation is not a Hermitian operator. Special care is
therefore needed when constructing the spectral representation of the DGF from
the modes. One solution is to use the biorthonormal construction discussed in
Ref. [47]. We envision that the method presented in this paper will eventually
be integrated with the existing Maxwell solver in the Octopus code[48, 49].
Such an integration would allow for the treatment of general electromagnetic
environments completely within Octopus.

##  VIII Conclusion

In this paper, we have presented a methodology combining macroscopic quantum
electrodynamics with quantum-electrodynamical density-functional theory which
provides a fully ab-initio description of coupled quantum light-matter
systems. To exemplify this approach, we have consided a benzene molecule
strongly coupled to a metallic spherical cavity and investigated the impact of
the both cavity radius and cavity loss on the nature of the light-matter
coupling. We have further investigated the effect of adding more molecules and
exchanging benzene with larger aromatic molecules. Together, these results
highlight the intricate interplay between the electronic structure of the
emitter and the environment in determining the nature of the light-matter
coupling. Our work therefore illustrates the importance of having a proper
description of both the electronic system and the electromagnetic environment
for a proper description of quantum light-matter interactions. This work sets
out the direction for more quantitative calculations in the future and also
opens the possibility for the proper treatment of real experimental setups. We
emphasize that the connection between the optical environment and the DGF is
not limited to setups similar to cavities, but instead provides a general way
to determine the electromagnetic spectral density of an arbitrary environment.
In addition to the QED setup, our method therefore also provides a way to
perform time-dependent density-functional theory in a lossy optical
environment and removes the need for artificial spectral broadening.

Finally, we have provided an easy to use tool that everyone can use to
generate cavity parameters for simple cavities such as the spherical
microcavity or a layered cavity.

##  IX Acknowledgements:

All calculations were performed using the computational facilities of the
Flatiron Institute. The Flatiron Institute is a division of the Simons
Foundation. We acknowledge support from the Max Planck-New York City Center
for Non-Equilibrium Quantum Phenomena.

K.S.T. acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program grant no. 773122
(LIMA). K.S.T. is a Villum Investigator supported by VILLUM FONDEN (grant no.
37789)

## References

  * Benz  _et al._ [2016] F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg, Single-molecule optomechanics in picocavities, Science 354, 726 (2016), https://www.science.org/doi/pdf/10.1126/science.aah5243 .
  * Pscherer  _et al._ [2021] A. Pscherer, M. Meierhofer, D. Wang, H. Kelkar, D. Martín-Cano, T. Utikal, S. Götzinger, and V. Sandoghdar, Single-molecule vacuum rabi splitting: Four-wave mixing and optical switching at the single-photon level, Phys. Rev. Lett. 127, 133603 (2021).
  * Do  _et al._ [2022] T. T. H. Do, M. Nonahal, C. Li, V. Valuckas, A. I. Kuznetsov, H. S. Nguyen, I. Aharonovich, and S. T. Ha, Room-temperature strong coupling in a single photon emitter-dielectric metasurface system (2022).
  * Munkhbat  _et al._ [2021] B. Munkhbat, A. Canales, B. Küçüköz, D. G. Baranov, and T. O. Shegai, Tunable self-assembled casimir microcavities and polaritons, Nature 597, 214 (2021).
  * Mueller  _et al._ [2020] N. S. Mueller, Y. Okamura, B. G. M. Vieira, S. Juergensen, H. Lange, E. B. Barros, F. Schulz, and S. Reich, Deep strong light–matter coupling in plasmonic nanoparticle crystals, Nature 583, 780 (2020).
  * Basov  _et al._ [2020] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio, Polariton panorama, Nanophotonics 10, 549 (2020).
  * Simpkins  _et al._ [2023] B. S. Simpkins, A. D. Dunkelberger, and I. Vurgaftman, Control, modulation, and analytical descriptions of vibrational strong coupling, Chemical Reviews (2023).
  * Ruggenthaler  _et al._ [2014] M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, Quantum-electrodynamical density-functional theory: Bridging quantum optics and electronic-structure theory, Physical Review A 90, 012508 (2014).
  * Kohn [1999] W. Kohn, Nobel lecture: Electronic structure of matter-wave functions and density functionals, Rev. Mod. Phys. 71, 1253 (1999).
  * Pellegrini  _et al._ [2015] C. Pellegrini, J. Flick, I. V. Tokatly, H. Appel, and A. Rubio, Optimized effective potential for quantum electrodynamical time-dependent density functional theory, Physical Review Letters 115, 093001 (2015).
  * Flick  _et al._ [2018a] J. Flick, C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, Ab initio optimized effective potentials for real molecules in optical cavities: Photon contributions to the molecular ground state, ACS photonics 5, 992 (2018a).
  * Flick  _et al._ [2019] J. Flick, D. M. Welakuh, M. Ruggenthaler, H. Appel, and A. Rubio, Light–matter response in nonrelativistic quantum electrodynamics, ACS photonics 6, 2757 (2019).
  * Yang  _et al._ [2021] J. Yang, Q. Ou, Z. Pei, H. Wang, B. Weng, Z. Shuai, K. Mullen, and Y. Shao, Quantum-electrodynamical time-dependent density functional theory within gaussian atomic basis, The Journal of Chemical Physics 155, 064107 (2021).
  * Schäfer  _et al._ [2022] C. Schäfer, J. Flick, E. Ronca, P. Narang, and A. Rubio, Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity, Nature Communications 13, 10.1038/s41467-022-35363-6 (2022).
  * Wang  _et al._ [2021] D. S. Wang, T. Neuman, J. Flick, and P. Narang, Light–matter interaction of a molecule in a dissipative cavity from first principles, The Journal of Chemical Physics 154, 104109 (2021).
  * Kudlis  _et al._ [2021] A. Kudlis, I. Iorsh, and I. Tokatly, Dissipation and spontaneous emission in quantum electrodynamical density functional theory based on optimized effective potential: A proof of concept study, arXiv preprint arXiv:2111.04523 (2021).
  * Scheel and Buhmann [2008] S. Scheel and S. Y. Buhmann, Macroscopic qed-concepts and applications, acta physica slovaca , 675 (2008).
  * Svendsen  _et al._ [2021] M. K. Svendsen, Y. Kurman, P. Schmidt, F. Koppens, I. Kaminer, and K. S. Thygesen, Combining density functional theory with macroscopic qed for quantum light-matter interactions in 2d materials, Nature communications 12, 1 (2021).
  * Feist  _et al._ [2021] J. Feist, A. I. Fernández-Domínguez, and F. J. García-Vidal, Macroscopic qed for quantum nanophotonics: emitter-centered modes as a minimal basis for multiemitter problems, Nanophotonics 10, 477 (2021).
  * Chew [1999] W. C. Chew,  _Waves and fields in inhomogenous media_ , Vol. 16 (John Wiley & Sons, 1999).
  * Woolley [1980] R. G. Woolley, Gauge invariant wave mechanics and the power-zienau-woolley transformation, J. Phys. A: Math. Gen 13, 2795 (1980).
  * Babiker and Loudon [1983] M. Babiker and R. Loudon, Derivation of the power-zienau-woolley hamiltonian in quantum electrodynamics by gauge transformation, Proc. R. Soc. London, A 385, 439 (1983).
  * Buhmann and Welsch [2008] S. Y. Buhmann and D. G. Welsch, Casimir-Polder forces on excited atoms in the strong atom-field coupling regime, Physical Review A - Atomic, Molecular, and Optical Physics 77, 10.1103/PhysRevA.77.012110 (2008).
  * Hümmer  _et al._ [2013] T. Hümmer, F. J. García-Vidal, L. Martín-Moreno, and D. Zueco, Weak and strong coupling regimes in plasmonic QED, Physical Review B - Condensed Matter and Materials Physics 87, 10.1103/PhysRevB.87.115419 (2013), arXiv:1209.1724 .
  * Novotny and Hecht [2012] L. Novotny and B. Hecht,  _Principles of nano-optics_ (Cambridge university press, 2012).
  * Tokatly [2013] I. V. Tokatly, Time-dependent density functional theory for many-electron systems interacting with cavity photons, Physical review letters 110, 233001 (2013).
  * Welakuh  _et al._ [2022] D. M. Welakuh, J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Frequency-dependent sternheimer linear-response formalism for strongly coupled light–matter systems, Journal of Chemical Theory and Computation 18, 4354 (2022), pMID: 35675628, https://doi.org/10.1021/acs.jctc.2c00076 .
  * Casida [1995] M. E. Casida, Time-dependent density functional response theory for molecules, in  _Recent Advances In Density Functional Methods: (Part I)_ (World Scientific, 1995) pp. 155–192.
  * Schäfer  _et al._ [2021] C. Schäfer, F. Buchholz, M. Penz, M. Ruggenthaler, and A. Rubio, Making ab initio qed functional (s): Nonperturbative and photon-free effective frameworks for strong light–matter coupling, Proceedings of the National Academy of Sciences 118, e2110464118 (2021).
  * Flick [2022] J. Flick, Simple exchange-correlation energy functionals for strongly coupled light-matter systems based on the fluctuation-dissipation theorem, Phys. Rev. Lett. 129, 143201 (2022).
  * Dung  _et al._ [2000] H. T. Dung, L. Knöll, and D.-G. Welsch, Spontaneous decay in the presence of dispersing and absorbing bodies: General theory and application to a spherical cavity, Physical Review A 62, 053804 (2000).
  * Blaber  _et al._ [2009] M. G. Blaber, M. D. Arnold, and M. J. Ford, Search for the ideal plasmonic nanoshell: the effects of surface scattering and alternatives to gold and silver, The Journal of Physical Chemistry C 113, 3041 (2009).
  * Note [1] It should be noted that the Drude model is not a good model for the dielectric properties of gold due to the neglect of d-band absorption. Additionally, it should be mentioned that the Drude model will fail to account for size dependent resonance shifts and broadening resulting from the nonlocal response and surface enhanced Landau dampening in small metallic structures [50, 51].
  * Rossi  _et al._ [2019] T. P. Rossi, T. Shegai, P. Erhart, and T. J. Antosiewicz, Strong plasmon-molecule coupling at the nanoscale revealed by first-principles modeling, Nature communications 10, 3336 (2019).
  * Akselrod  _et al._ [2014] G. M. Akselrod, C. Argyropoulos, T. B. Hoang, C. Ciracì, C. Fang, J. Huang, D. R. Smith, and M. H. Mikkelsen, Probing the mechanisms of large purcell enhancement in plasmonic nanoantennas, Nature Photonics 8, 835 (2014).
  * Gjerding  _et al._ [2017] M. N. Gjerding, M. Pandey, and K. S. Thygesen, Band structure engineered layered metals for low-loss plasmonics, Nature communications 8, 15133 (2017).
  * Kuznetsov  _et al._ [2016] A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S. Kivshar, and B. Luk’yanchuk, Optically resonant dielectric nanostructures, Science 354, aag2472 (2016).
  * Baranov  _et al._ [2017] D. G. Baranov, D. A. Zuev, S. I. Lepeshov, O. V. Kotov, A. E. Krasnok, A. B. Evlyukhin, and B. N. Chichkov, All-dielectric nanophotonics: the quest for better materials and fabrication techniques, Optica 4, 814 (2017).
  * Svendsen  _et al._ [2022] M. K. Svendsen, H. Sugimoto, A. Assadillayev, D. Shima, M. Fujii, K. S. Thygesen, and S. Raza, Computational discovery and experimental demonstration of boron phosphide ultraviolet nanoresonators, Advanced Optical Materials 10, 2200422 (2022).
  * Ebbesen [2016] T. W. Ebbesen, Hybrid light–matter states in a molecular and material science perspective, Accounts of chemical research 49, 2403 (2016).
  * Genet  _et al._ [2021] C. Genet, J. Faist, and T. W. Ebbesen, Inducing new material properties with hybrid light–matter states, Physics Today 74, 42 (2021).
  * Sidler  _et al._ [2022] D. Sidler, M. Ruggenthaler, C. Schäfer, E. Ronca, and A. Rubio, A perspective on ab initio modeling of polaritonic chemistry: The role of non-equilibrium effects and quantum collectivity, The Journal of Chemical Physics 156, 230901 (2022).
  * Ruggenthaler  _et al._ [2022] M. Ruggenthaler, D. Sidler, and A. Rubio, Understanding polaritonic chemistry from ab initio quantum electrodynamics, arXiv preprint arXiv:2211.04241 (2022).
  * Dutra and Knight [1996] S. Dutra and P. Knight, Spontaneous emission in a planar fabry-pérot microcavity, Physical Review A 53, 3587 (1996).
  * Thomas  _et al._ [2019] A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, T. C. J. George, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, Tilting a ground-state reactivity landscape by vibrational strong coupling, Science 364, 615 (2019).
  * Gu  _et al._ [2021] J. Gu, V. Walther, L. Waldecker, D. Rhodes, A. Raja, J. C. Hone, T. F. Heinz, S. Kéna-Cohen, T. Pohl, and V. M. Menon, Enhanced nonlinear interaction of polaritons via excitonic rydberg states in monolayer WSe2, Nature Communications 12, 10.1038/s41467-021-22537-x (2021).
  * Chen  _et al._ [2010] Y. Chen, T. R. Nielsen, N. Gregersen, P. Lodahl, and J. Mørk, Finite-element modeling of spontaneous emission of a quantum emitter at nanoscale proximity to plasmonic waveguides, Physical Review B 81, 125431 (2010).
  * Tancogne-Dejean  _et al._ [2020] N. Tancogne-Dejean, M. J. Oliveira, X. Andrade, H. Appel, C. H. Borca, G. Le Breton, F. Buchholz, A. Castro, S. Corni, A. A. Correa, _et al._ , Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems, The Journal of chemical physics 152, 124119 (2020).
  * Jestädt  _et al._ [2019] R. Jestädt, M. Ruggenthaler, M. J. Oliveira, A. Rubio, and H. Appel, Light-matter interactions within the ehrenfest–maxwell–pauli–kohn–sham framework: fundamentals, implementation, and nano-optical applications, Advances in Physics 68, 225 (2019).
  * Mortensen  _et al._ [2014] N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, and S. I. Bozhevolnyi, A generalized non-local optical response theory for plasmonic nanostructures, Nature communications 5, 3809 (2014).
  * Svendsen  _et al._ [2020] M. Svendsen, C. Wolff, A.-P. Jauho, N. A. Mortensen, and C. Tserkezis, Role of diffusive surface scattering in nonlocal plasmonics, Journal of Physics: Condensed Matter 32, 395702 (2020).
  * Mukamel [1999] S. Mukamel,  _Principles of nonlinear optical spectroscopy_ , 6 (Oxford University Press on Demand, 1999).
  * Flick  _et al._ [2018b] J. Flick, N. Rivera, and P. Narang, Strong light-matter coupling in quantum chemistry and quantum photonics, Nanophotonics 7, 1479 (2018b).
  * Medina  _et al._ [2021] I. Medina, F. J. García-Vidal, A. I. Fernández-Domínguez, and J. Feist, Few-mode field quantization of arbitrary electromagnetic spectral densities, Physical Review Letters 126, 093601 (2021).
  * Sánchez-Barquilla  _et al._ [2022] M. Sánchez-Barquilla, F. J. García-Vidal, A. I. Fernández-Domínguez, and J. Feist, Few-mode field quantization for multiple emitters, Nanophotonics 11, 4363 (2022).

##  X Supplementary Information

##  Supplementary note A Length-gauge Pauli Fierz Hamiltonian in terms of the
emitter centered modes

Our starting point is the PZW gauge Hamiltonian describing an electronic
system coupled to the MQED field within the dipole approximation [17].
Neglecting magnetic interactions, the Hamiltonian reads,

$\displaystyle\mathcal{H}=$
$\displaystyle\sum_{\alpha}\frac{\bm{\hat{p}}_{\alpha}^{2}}{2m_{\alpha}}+\frac{1}{2\epsilon_{0}}\int
d\bm{r}\bm{\hat{P}}(\bm{r})^{2}$
$\displaystyle+\sum_{\lambda}\int_{0}^{\infty}d\omega\hbar\omega\int
d\bm{r}\bm{\hat{f}}_{\lambda}^{\dagger}(\bm{r},\omega)\cdot\bm{\hat{f}}_{\lambda}(\bm{r},\omega)-\bm{\hat{d}}\cdot\bm{\hat{E}}(\bm{r}_{0}).$
(25)

To connect this to the Hamiltonian used in the Octopus code QEDFT
implementation[12],

$\displaystyle\hat{H}_{0}=$
$\displaystyle\hat{H}_{e}+\frac{1}{2}\sum^{N}_{\alpha=1}\left(\hat{p}_{\alpha}^{2}+\omega_{\alpha}^{2}\left[\hat{q}_{\alpha}+\frac{\bm{\lambda}_{\alpha}}{\omega_{\alpha}}\cdot\hat{\bm{R}}\right]^{2}\right),$
(26)

we follow Ref. [15] in defining the cavity field strength,
$\bm{\lambda}_{i}(\omega)=-\left(\frac{2}{\hbar\omega}\right)^{1/2}e\bm{E}_{i}(\bm{r}_{0},\omega)$,
and rewrite the emitter-centered representation of the MQED field in Eq. 6
from the manuscript as,

$\displaystyle\bm{\hat{E}}(\bm{r}_{0})=-\sum_{i}\int_{0}^{\infty}d\omega\left(\frac{\hbar\omega}{2e^{2}}\right)^{1/2}\bm{\lambda}_{i}(\omega)\hat{C}_{i}(\omega)+\mathrm{h.c.}.$
(27)

Next we want to change from the ladder operators to the canonical operators,
$\hat{q}_{i}(\omega)=\left(\frac{\hbar}{2\omega}\right)^{1/2}(\hat{C}_{i}(\omega)+\hat{C}^{\dagger}_{i}(\omega))$
and
$\hat{p}_{i}(\omega)=\left(\frac{\hbar\omega}{2}\right)^{1/2}(\hat{C}_{i}(\omega)-\hat{C}^{\dagger}_{i}(\omega))$,
again to be consistent with the Octopus implementation of QEDFT. This is
possible as long as the orthogonalization scheme used in the construction of
$\hat{C}_{i}(\omega)$ results in real valued $\bm{V}(\omega)$ matrices. With
this change, the electric field operator becomes,

$\displaystyle\bm{\hat{E}}(\bm{r}_{0})=-\sum_{i}\int_{0}^{\infty}d\omega\frac{\omega}{e}\bm{\lambda}_{i}(\omega)\hat{q}_{i}(\omega).$
(28)

The idea is then to use Eq. 28 to put Eq. A in a form similar to the
Hamiltonian in Eq. 26. The intuitive generalization of the Hamiltonian to the
continuous case suggests that,

$\displaystyle\mathcal{H}=$
$\displaystyle\mathcal{H}^{\mathrm{C}}_{\mathrm{Mat}}+$
$\displaystyle\frac{1}{2}\sum_{i}\int_{0}^{\infty}d\omega\left\\{\hat{p}_{i}(\omega)^{2}+\omega^{2}\left[\hat{q}_{i}(\omega)+\frac{\bm{\lambda}_{i}(\omega)}{\omega}\cdot\bm{\hat{R}}\right]^{2}\right\\}.$
(29)

Here $\hat{\bm{R}}=\sum_{i}\hat{\bm{r}}_{i}$ is the center of mass position of
the charges making up the emitter which is related to the dipole moment as
$\hat{\bm{d}}=-e\hat{\bm{R}}$. The question is then under what conditions this
holds. Inserting the field expansion into the emitter-field coupling term in
the Hamiltonian results in,

$\displaystyle\hat{H}_{\mathrm{int}}=-\hat{\bm{d}}\cdot\bm{\hat{E}}(\bm{r}_{0})=\sum_{i}\int_{0}^{\infty}d\omega\omega\left[\bm{\lambda}_{i}(\omega)\cdot\bm{\hat{R}}\right]\hat{q}_{i}(\omega).$
(30)

which is the straightforward generalization of the lossless case.

The other term in which the field parameters appear is the dipole self energy
term. This term should generally be derived in terms of the transverse
projection of the full polarization term[52],

$\displaystyle\bm{\hat{P}}_{\bot}(\bm{r})=\frac{2}{3}\bm{\hat{P}}(\bm{r})+\int
d\bm{r^{\prime}}\,\bm{\hat{T}}_{\bot}(\bm{r},\bm{r^{\prime}})\cdot\bm{\hat{P}}(\bm{r^{\prime}}).$
(31)

However, the direct application of this projection is computationally
cumbersome. In the lossless case, the Helmholtz equation is Hermitian and its
transverse solutions form a orthogonal and complete basis for the transverse
space. Consequently, the transverse polarization can be expanded in terms of
these modes and the projection can be simplified significantly. This procedure
leads to the standard form of the dipole self-energy term,

$\displaystyle\frac{1}{2\epsilon_{0}}\int
d\bm{r}\left[\bm{\hat{P}}_{\bot}(\bm{r})\right]^{2}=\frac{1}{2}\sum_{\alpha}\left[\bm{\lambda}_{\alpha}\cdot\hat{\bm{R}}\right]^{2}.$
(32)

The fact that the transverse modes can be used to span the transverse
projector indicates that any truncation of the photon Hilbert space needs to
be performed carefully as it indirectly results in a truncation of the
transverse basis and therefore also the transverse polarization. The two
should therefore generally always be truncated consistently[53].

The expansion of the transverse polarization field in the presence of losses
is more complicated. If the cavity losses are not true absorptive losses, but
simply a result of e.g. finite reflectivity of the cavity mirrors, the modes
will globally obey,

$\displaystyle\int
d\bm{r}\bm{E}_{i}(\bm{r},\omega_{1})\cdot\bm{E}_{j}(\bm{r},\omega_{2})=\delta_{ij}\delta(\omega_{1}-\omega_{2}),$
(33)

and thus remain orthogonal and complete. In this case the transverse
polarization can be straightforwardly generalized,

$\displaystyle\frac{1}{2\epsilon_{0}}\int
d\bm{r}\left[\bm{\hat{P}}_{\bot}(\bm{r})\right]^{2}=\frac{1}{2}\sum_{i}\int_{0}^{\infty}d\omega\left[\bm{\lambda}_{i}(\omega)\cdot\bm{R}\right]^{2}.$
(34)

In the presence of true absorption losses it is not possible to define the
transverse projector in terms of the field modes because these no longer form
a complete, orthonormal set. We speculate that the transverse projection
should instead by defined in terms of a biorthonormal construction as
discussed in Ref. [47]. From the point of the emitter it should not be
possible to distinguish these two cases because it only samples the field at
is position. We thus speculate that, at least for a single emitter position,
the true form of the dipole self energy term in the presence of losses will
functionally be no different than Eq. 34. For multiple emitter positions it is
unlikely that the dipole self-energy term can be expanded directly in terms of
the mode functions, and we will explore this further in future research.

##  Supplementary note B Dyadic Green’s Function of the spherical microcavity

We consider the spherically layered microcavity from Ref. [31]. To use the
emitter centered framework, the DGF must be derived. In general, the source
and field points ($\bm{r}$ and $\bm{r^{\prime}}$ respectively) can be in
either the same layer or different layers. Consequently, the DGF for the
reflected field is a 9 component object, where each component is a 3D dyad.
Labelling the possible combinations of source and field points by two extra
indices $m,n=1,2,3$, the full DGF of the system can be written as,

$\displaystyle\bm{G}_{mn}(\bm{r},\bm{r^{\prime}},\omega)=\bm{G}^{\mathrm{vac}}(\bm{r},\bm{r^{\prime}},\omega)\delta_{nm}+\bm{G}^{\mathrm{ref}}_{mn}(\bm{r},\bm{r^{\prime}},\omega).$
(35)

Due to the spherical symmetry of the problem, the DGF of the system is most
efficiently expanded onto vector spherical harmonics. The 9 different
components of the reflection contribution can be worked out by invoking the
electromagnetic boundary conditions at the material interfaces. In this work
we consider the situation where the emitter is placed in the inner region of
the cavity, and it is therefore only necessary to consider the $n=m=3$
component of the DGF.

For a general emitter position inside the cavity, it is necessary to carefully
converge the number of vector spherical harmonics used in the calculation of
the DGF. However, if the emitter is placed exactly in the center of the cavity
things simplify significantly. In this case, only the lowest order transverse
magnetic mode contributes and the reflection DGF can be written as[31],

$\displaystyle\bm{G}^{\mathrm{ref}}_{33}(\bm{r},\bm{r^{\prime}},\omega)|_{r,r^{\prime}\rightarrow
0}=\frac{i\omega}{6\pi c}r^{N}_{n=1}(\omega)\bm{\mathbbm{1}}.$ (36)

where $r^{N}_{n=1}(\omega)\bm{\mathbbm{1}}$ is the reflection coefficient for
the lowest order transverse magnetic mode. If we assume that region 2 is very
thick we can write[31],

$\displaystyle
r^{N}_{n=1}(\omega)=\frac{\left[i+\rho(n(\omega)+1)-i\rho^{2}n(\omega)-\frac{\rho^{3}n(\omega)^{2}}{n(\omega)+1}\right]e^{i\rho}}{\mathrm{sin}\rho-\rho(\mathrm{cos}\rho+in(\omega)\mathrm{sin}\rho)+i\rho^{2}n(\omega)\mathrm{cos}\rho-\rho^{3}(\mathrm{cos}\rho-
in(\omega)\mathrm{sin}\rho)\frac{n(\omega)^{2}}{n(\omega)^{2}-1}},$ (37)

where $\rho=R\omega/c$ is the standard size parameter from Mie theory[25].

The contribution from the vacuum DGF is
$\mathrm{Im}\bm{G}^{\mathrm{vac}}(\bm{r},\bm{r},\omega)=\frac{\omega}{6\pi
c}\bm{\mathbbm{1}}$. The imaginary part of the full DGF, evaluated in the
center of the cavity, thus becomes,

$\displaystyle\mathrm{Im}\bm{G}(\bm{0},\bm{0},\omega)=\frac{\omega}{6\pi
c}\left[1+\mathrm{Re}r^{N}_{n=1}(\omega)\right]\bm{\mathbbm{1}}.$ (38)

Denoting the unit vector in the $i$’th direction $\bm{n}_{i}$ it can be
noticed that,
$\bm{n}_{i}\cdot\mathrm{Im}\bm{G}(\bm{0},\bm{0},\omega)\cdot\bm{n}_{j}\propto\delta_{ij}$
which means that the mode orthogonalization needed in the emitter centered
formulation is trivial in this case. The cavity field parameters can therefore
be derived directly,

$\displaystyle\bm{\lambda}_{i}(\omega)=-e\left(\frac{\omega^{2}}{3\pi^{2}\epsilon_{0}c^{3}}\left[1+\mathrm{Re}r^{N}_{n=1}(\omega)\right]\right)^{1/2}\bm{n}_{i}.$
(39)

Notice that the spherical symmetry of the cavity means that at each frequency
there are three orthogonal modes with identical field strengths, differing
only by their spatial orientation $i$. Any dipole orientation therefore
couples identically to the cavity modes if the emitter is placed in the center
of the cavity.

Supplementary Figure 1: Spherical microcavity setup: (a) Illustration of the
three layer spherical microcavity setup. (b) The Drude model dielectric
function of gold. (c) The modes of the spherical microcavity plotted for two
different radii, $R=140$nm in blue and $R=450$nm in orange. The mode structure
is shown with a sampling density of $10$ points/meV. (d) The impact from the
Drude dampening parameter on the cavity resonances, illustrated by focusing on
the mode around 7.1 eV in the cavity with $R=140$nm.

##  Supplementary note C Discretization of the cavity coupling parameters

In practice, we need to discretize the continuous expressions for the cavity
field strengths, $\bm{\lambda}(\omega)$. The idea is then to sample densely
enough to mimic a true continuum of modes. This means that the sampling
density should be converged in the simulations. If we assume that we have
uniformly spaced modes we can perform the simple discretization,

$\displaystyle\int_{0}^{\infty}d\omega\rightarrow\sum_{k}\Delta\omega.$ (40)

In the continuous formulation, the operators $\hat{C}_{i}(\omega)$ have units
of $\mathrm{s}^{1/2}$ and the field, $\bm{E}(\bm{r},\omega)$ has units of
$\mathrm{V}\mathrm{s}^{1/2}/\mathrm{m}$. It therefore makes sense to define
the a new discretized field and creating/annihilation operators as,

$\displaystyle\bm{E}_{k}(\omega_{k})\equiv\sqrt{\Delta\omega}\bm{E}_{i}(\bm{0},\omega_{k}),$
(41)
$\displaystyle\hat{a}^{(\dagger)}_{k}\equiv\sqrt{\Delta\omega}\hat{C}^{(\dagger)}_{i}(\omega).$
(42)

This also means that we should write,

$\displaystyle\bm{\lambda}_{k}=-e\left(\frac{\Delta\omega\omega_{k}^{2}}{3\pi^{2}\epsilon_{0}c^{3}}\left[1+\mathrm{Re}\,r^{N}_{m=1}(\omega_{k})\right]\right)^{1/2}\hat{n}_{i}.$
(43)

##  Supplementary note D Computational details for the QEDFT calculations

We use the publicly available real space psuedo-potential DFT code
OCTOPUS[48]. Molecular geometries are optimized with the LDA exchange-
correlation functional, on a real space grid consisting of spheres of radius
6Å around each atom with a grid spacing of 0.16Å.

We also calculate the ground state using the LDA exchange-correlation
functional on a real space grid consisting of spheres of radius 6 Å around
each atom with a grid spacing of 0.08 Å. We perform a fixed-density
calculation on top of the ground state calculation to determine the excited
states needed for the Casida calculation for which we use 500 excited states.
These parameters ensure that the relevant transitions are converged to within
$1$ meV/atom.

We employ uniform sampling of the electromagnetic environment and converge
each calculation separately. In general, we find good convergence when
sampling each cavity peak with around 1000 photon modes. While the current
work employs uniform sampling, we note that more efficient methods for
sampling the electromagnetic environment has recently been introduced in the
context of master equation approaches [54, 55]. In the future, it would
therefore be interesting to investigate whether similar approaches can be
applied to photon mode sampling within the framework presented in the current
work.

##  Supplementary note E The Fabry Perot cavity

In this supplementary note we briefly discuss the widely used Fabry Perot
cavity (FPC). As shown in Supplementary Figure 2, the stereotypical FPC is a
stratified system. Consequently, its DGF can be conveniently represented in
the angular spectrum representation. Considering just the cavity region, the
DGF can be written as[20],

$\displaystyle\bm{G}(\bm{r},\bm{r}^{\prime},\omega)=\frac{i}{8\pi^{2}}\int\frac{d\bm{q}}{k_{z}q^{2}}\left[\bm{M}(\bm{q},\bm{r},\bm{r^{\prime}})+\bm{N}(\bm{q},\bm{r},\bm{r^{\prime}})\right].$
(44)

Here $\bm{q}$ is the in-plane wave vector, $k_{z}=\sqrt{k^{2}-q^{2}}$ is the
out of plane wave vector where $k$ is the wave number in the medium and
$q=|\bm{q}|$. $\bm{M}$ and $\bm{N}$ are the tensorial contributions from the
TE and TM polarized modes respectively. These tensors can be found by solving
the scattering from a point source in a layered medium as discussed in detail
in e.g. Refs. [20, 25]. Importantly, if we consider a single emitter position
in the FPC, we only need information about the $\bm{r}=\bm{r^{\prime}}$ part
of the DGF to describe the coupling to the electromagnetic environment.
Conveniently, at $\bm{r}=\bm{r^{\prime}}$ both $\bm{M}$ and $\bm{N}$ are
diagonal matrices when expressed in cylindrical coordinates with the $z$-axis
perpendicular to the mirrors,

$\displaystyle k_{x}=q\,\mathrm{cos}\,\phi,\hskip
28.45274ptk_{y}=q\,\mathrm{sin}\,\phi,\hskip 28.45274ptk_{z}=k_{z}.$ (45)

This is also the natural choice because the in-plane wave number
$q=\sqrt{k_{x}^{2}+k_{y}^{2}}$ is the natural variable in the stratified FPC.
The integration over $\phi$ can be performed analytically which results in,

$\displaystyle\bm{G}^{\mathrm{TE}}(\bm{r},\bm{r},\omega)=\int
dqq\int_{0}^{2\pi}d\phi\,\bm{M}(\bm{q},\bm{r},\bm{r})=\frac{i}{8\pi}\int
dq\frac{q}{k_{z}}\begin{bmatrix}1&0&0\\\ 0&1&0\\\
0&0&0\end{bmatrix}R_{\parallel}^{\mathrm{TE}}(z),$ (46)

and,

$\displaystyle\bm{G}^{\mathrm{TM}}(\bm{r},\bm{r},\omega)=\int
dqq\int_{0}^{2\pi}d\phi\,\bm{N}(\bm{q},\bm{r},\bm{r})=\frac{i}{8\pi k^{2}}\int
dq\begin{bmatrix}qk_{z}R_{\parallel}^{\mathrm{TM}}(z)&0&0\\\
0&qk_{z}R_{\parallel}^{\mathrm{TM}}(z)&0\\\
0&0&2q^{3}/k_{z}R_{\bot}^{\mathrm{TM}}(z).\end{bmatrix}$ (47)

Taking the emitter to be placed at $z=0$ and denoting the distance between the
emitter and the top and bottom mirrors as $t$ and $b$ respectively, the
reflection functions can be written as,

$\displaystyle
R_{\bot}^{\mathrm{TM}}(z,q)=\frac{(1+r^{\mathrm{TM}}_{cb}e^{2ik_{z}b})(1+r^{\mathrm{TM}}_{ct}e^{2ik_{z}t})}{1-r^{\mathrm{TM}}_{cb}r^{\mathrm{TM}}_{ct}e^{2ik_{z}d}}$
(48) $\displaystyle
R_{\parallel}^{\mathrm{TE}}(z,q)=\frac{(1+r^{\mathrm{TE}}_{cb}e^{2ik_{z}b})(1+r^{\mathrm{TE}}_{ct}e^{2ik_{z}t})}{1-r^{\mathrm{TE}}_{cb}r^{\mathrm{TE}}_{ct}e^{2ik_{z}d}}$
(49) $\displaystyle
R_{\parallel}^{\mathrm{TM}}(z,q)=\frac{(1-r^{\mathrm{TM}}_{cb}e^{2ik_{z}b})(1-r^{\mathrm{TM}}_{ct}e^{2ik_{z}t})}{1-r^{\mathrm{TM}}_{cb}r^{\mathrm{TM}}_{ct}e^{2ik_{z}d}}$
(50)

The subscripts $cb$ and $ct$ refers to the interfaces between the cavity
region and the bottom- and top mirror respectively. We use the following
definition of the Fresnel coefficients for light incident on the interface
between two regions $a$ and $b$ from region $a$,

$\displaystyle
r_{ab}^{\mathrm{TE}}=\frac{k^{a}_{z}-k^{b}_{z}}{k^{a}_{z}+k^{b}_{z}},$ (51)
$\displaystyle
r_{ab}^{\mathrm{TM}}=\frac{\epsilon_{b}k^{a}_{z}-\epsilon_{a}k^{b}_{z}}{\epsilon_{b}k^{a}_{z}+\epsilon_{a}k^{b}_{z}}.$
(52)

Note that these Fresnel coefficients could be replaced with the generalized
Fresnel coefficients if more complicated mirrors such as e.g. distributed
Bragg reflectors were considered.

Supplementary Figure 2: The Fabry Perot cavity: a) Paradigmatic sketch of a
Fabry Perot Cavity. b) Purcell enhancement as a function of $d/\lambda_{0}$
for both vertical and horizontal dipole orientations when considering
idealized mirrors with a frequency independent reflectivity $R=0.95$. b)
Purcell enhancement as a function of $d/\lambda_{0}$ for both vertical and
horizontal dipole orientations when considering metal mirrors with dielectric
functions described by the Drude model in Eq. 23.

###  E.1 Light-matter interaction in the Fabry-Perot cavity

Because the DGF is diagonal at $\bm{r}=\bm{r^{\prime}}$, there is no coupling
between the horizontal and vertical dipole orientations by feedback from the
electromagnetic environment. Consequently, the local density of states (LDOS)
for the two dipole orientations,

$\displaystyle\rho(\bm{r},\omega)=\frac{6\omega}{\pi
c^{2}}\hat{n}_{D}\cdot\mathrm{Im}\left[\bm{G}(\bm{r},\bm{r},\omega)\right]\cdot\hat{n}_{D}.$
(53)

for a given dipole orientation $\hat{n}_{D}$ becomes a direct proxy for the
light-matter coupling strength. It can therefore be used evaluate the
potential of FPCs for realizing single/few emitter strong coupling. One
concrete measure of the FPCs alterations of the electromagnetic environment is
its modifications of the free space density of states. This is measured
directly by the Purcell enhancement,

$\displaystyle
P_{F}(\bm{r},\omega)=\frac{\rho(\bm{r},\omega)}{\rho_{0}(\omega)},$ (54)

which is given as a rescaling of the LDOS by the free space density of states
of the electromagnetic field, $\rho_{0}(\omega)$. To explore this, we consider
two different FPCs, both of which consist of a central cavity region made of
vacuum, surrounded by mirrors. In one instance, we consider idealized mirrors
with a constant reflectivity of 0.95, and in the other we consider mirrors
made of a Drude metal. In the latter case, we consider frequencies below the
plasma frequency such that the mirrors retain their reflectivity.
Supplementary Figure 2 shows the Purcell enhancement as a function of mirror
distance, $d$, for a fixed emitter wavelength, $\lambda_{0}$, for idealized
mirrors (b) and for Drude mirrors (c). In both cases, the horizontal and
vertical dipole orientations are considered separately because, as mentioned
above, the cavity does not couple the two.

Starting with Supplementary Figure 2 b) for the horizontal dipole orientation,
we observe only relatively broad cavity resonances with modest Purcell
enhancements of maximum $\sim 2.5$. This reflects the relatively weak
concentration of the electromagnetic field in the FPC and emphasizes that the
case with idealized mirrors is ill-suited for single point emitter strong
coupling. We emphasize that the widths of these resonances are only weakly
linked to the finite mirror reflectivity, and instead reflect the fact that
the cavity features dispersion in plane. As such, these peaks would retain a
finite width even in the limit of perfect mirror reflectivity [44]. For
$d<\lambda_{0}/2$ we observe a near complete suppression of the coupling of
the horizontal dipole orientation to the electromagnetic environment. This
happens because the horizontal dipole will always couple to the $q=0$ mode.
This mode will not exist for $d<\lambda_{0}/2$ with perfect mirrors and be
strongly suppresses for near-perfect mirrors. This suppression can also be
understood as a result of the mirror charge effect which will result in
emission suppression for emitters with horizontal transition dipole moments
when these are placed close to idealized reflective surfaces [17]. For the
vertical dipole orientation, we generally observe relatively weak, near unity
Purcell enhancement, except when $d$ becomes smaller than $\sim\lambda_{0}/2$.
From this point, the Purcell enhancement diverges as $d\rightarrow 0$. This
divergence can be understood as a result of concentration of energy into the
mode propagating parallel to the mirrors when the mirror distance is decreased
[44]. Alternatively, it can also be more intuitively understood as a result of
the mirror charge effect which enhances emission for vertical dipole
orientations. For a single reflective interface, this enhancement would be a
simple factor of 2 [17]. However, in the double mirror case, the image charge
dipole in one mirror is enhanced by its mirror charge dipole in the second
mirror and vice versa to infinity. Therefore, the Purcell enhancement diverges
as $d\rightarrow 0$ instead of approaching a simple factor of 2. When we
consider the more realistic situation with mirrors described by the Drude
model in Figure 2 c) we observe more or less the same as was the case with
idealized mirrors. The only major difference is that we observe a divergence
for the horizontal dipole orientation at small mirror distances. This happens
due to near field coupling to the metal primarily representing quenching. We
note that this also happens to the vertical dipole orientation but that this
is less visible because of the mirror charge effect already causing a
divergence like behaviour for $d\rightarrow 0$ as was also the case with
idealized mirrors.

The above discussion shows that the FPC provides only modest modification of
the electromagnetic environment relative to the free space case with broad
resonances whose width arise from the in-place dispersion of the cavity. The
exception to this statement is for very small mirror spacing. However, in this
latter regime the electromagnetic environment does not show resonant
behaviour. For these reasons, we conclude that the FPC is ill-suited for
single- and few emitter vacuum strong coupling. However, the FPC can still be
suited for collective strong coupling [45] and coupling to extended systems
[46] where the extended modes of the electromagnetic environment can be
sampled more effectively.


# Effect of an uniaxial single-ion anisotropy on the quantum and thermal
entanglement of a mixed spin-(1/2,$S$) Heisenberg dimer

Hana Vargová11footnotemark: 1<EMAIL_ADDRESS>Jozef Strečka22footnotemark: 2
Natália Tomašovičová33footnotemark: 3 Institute of Experimental Physics,
Slovak Academy of Sciences, Watsonova 47, 040 01 Košice, Slovakia  Department
of Theoretical Physics and Astrophysics, Faculty of Science, P. J. Šafárik
University, Park Angelinum 9, 040 01 Košice, Slovakia

###### Abstract

Exact analytical diagonalization is used to study the bipartite entanglement
of the antiferromagnetic mixed spin-(1/2,$S$) Heisenberg dimer (MSHD) with the
help of negativity. Under the assumption of uniaxial single-ion anisotropy
affecting higher spin-$S$ ($S\\!>\\!1/2$) entities only, the ground-state
degeneracy $2S$ is partially lifted and the ground state is two-fold
degenerate with the total magnetization per dimer $\pm(S\\!-\\!1/2)$. It is
shown that the largest quantum entanglement is reached for the
antiferromagnetic ground state of MSHD with arbitrary half-odd-integer spins
$S$, regardless of the exchange and single-ion anisotropies. Contrary to this,
the degree of a quantum entanglement in MSHD with an integer spin $S$ for the
easy-plane single-ion anisotropy, exhibits an increasing tendency with an
obvious spin-$S$ driven crossing point. It is shown that the increasing spin
magnitude is a crucial driving mechanism for an enhancement of a threshold
temperature above which the thermal entanglement vanishes. The easy-plane
single-ion anisotropy together with an enlargement of the spin-$S$ magnitude
is other significant driving mechanism for an enhancement of the thermal
entanglement in MSHD.

###### keywords:

Heisenberg dimer , bipartite entanglement , exact diagonalization

††journal: Journal of Magnetism and Magnetic Materials

##  1 Introduction

One of the most fascinating feature of quantum mechanical systems is an
entanglement, referred to non-local correlations among the system
constituents. A discovery of its importance in a quantum information
processing [1] has initiated an intensive investigation of this unconventional
phenomenon in various quantum systems. From the application point of view
there naturally arises a question how to stabilize an entanglement at non-zero
(ideally room) temperature, because an entanglement is rapidly reduced through
thermal fluctuations. The Heisenberg dimer represents a profound theoretical
model for analysing the bipartite entanglement, because its mathematical
simplicity and quantum nature merge together. It was elucidated during last
several years, that the implementation of additional stimuli, like non-uniform
magnetic field, seems to be a relevant mechanism for an accrument of
entanglement stability [2, 3, 4, 5, 6]. Another mechanism for an enhancement
of entanglement lies in an anisotropy of a quantum system like the exchange
anisotropy [7, 3, 5, 8, 9, 10], Dzyaloshinskii-Moriya anisotropy [8, 4, 9] or
diversity of spin parity [8, 11, 12]. Although all aforementioned mechanisms
may enlarge the threshold temperature, they unfortunately simultaneously
reduce a strength of the entanglement. Inspired by the findings of Solano-
Carrillo et al. [13, 14], which imply that the uniaxial single-ion anisotropy
may enhance the degree of thermal entanglement, we will investigate the
simultaneous effect of an uniaxial single-ion anisotropy and the spin
magnitude on a bipartite entanglement of MSHD.

##  2 Model and Method

Let us consider MSHD consisting of two different spins defined through the
Hamiltonian

$\displaystyle\hat{\cal H}$
$\displaystyle=J\left\\{\Delta(\hat{\mu}^{x}\hat{S}^{x}\\!+\\!\hat{\mu}^{y}\hat{S}^{y})\\!+\\!\hat{\mu}^{z}\hat{S}^{z}\right\\}\\!+\\!D(\hat{S}^{z})^{2},$
(1)

where $\hat{\mu}^{\alpha}$ and $\hat{S}^{\alpha}$ ($\alpha\\!=\\!x,y,z$)
denote spatial components of spin-1/2 and spin-$S$ ($S\\!\geq\\!1$) operators,
respectively. The parameter $J$ is the Heisenberg exchange interaction,
$\Delta$ is the XXZ exchange anisotropy and $D$ is the uniaxial single-ion
anisotropy, which affects higher-spin $S\\!\geq\\!1$ only. The complete
spectrum of the eigenvalues and eigenvectors of the Hamiltonian (1) has been
derived in our previous paper [15], so we will list here for simplicity the
final formulas for eigenenergies and eigenvectors only

$\displaystyle\varepsilon_{\pm(S+\frac{1}{2})}\\!=\\!\frac{S}{2}(J\\!+\\!2DS)\rightarrow|\pm(S\\!+\\!1/2)\rangle\\!=\\!|1/2,S\rangle,$
(2)
$\displaystyle\varepsilon^{\mp}_{S^{z}_{t}}=-\frac{P_{S^{z}_{t}}}{4}\\!\mp\\!\frac{1}{4}\sqrt{R_{S^{z}_{t}}^{2}\\!+\\!Q_{S^{z}_{t}}}$
(3) $\displaystyle\hskip
28.45274pt\rightarrow\;|\left({S^{z}_{t}}\right)_{\mp}\rangle=c^{\mp}_{S^{z}_{t}}|1/2,S^{z}\rangle\\!\mp\\!c^{\pm}_{S^{z}_{t}}|\\!-1/2,S^{z}\\!+\\!1\rangle,$

where $S^{z}\\!=\\!-S,\dots,S$ and
$S^{z}_{t}\\!=\\!S^{z}\\!+\\!\mu^{z}\\!=\\!-S\\!-\\!1/2,\dots,S\\!+\\!1/2$.
The coefficients $P_{S^{z}_{t}}$, $R_{S^{z}_{t}}$, $Q_{S^{z}_{t}}$,
$c^{\pm}_{S^{z}_{t}}$, emergent in Eq. (3) are defined as follows

$\displaystyle
P_{S^{z}_{t}}=(J\\!-\\!2D)\\!-\\!D(2S^{z}_{t}\\!-\\!1)(2S^{z}_{t}\\!+\\!1),$
(4) $\displaystyle
Q_{S^{z}_{t}}=(J\Delta)^{2}[4S(S\\!+\\!1)\\!-\\!(2S^{z}_{t}\\!-\\!1)(2S^{z}_{t}\\!+\\!1)],$
(5) $\displaystyle
R_{S^{z}_{t}}=2(J\\!-\\!2D)S^{z}_{t},c_{S^{z}_{t}}^{\mp}\\!=\\!\frac{1}{\sqrt{2}}\\!\sqrt{\\!1\\!\mp\\!\frac{R_{S^{z}_{t}}}{\sqrt{R_{S^{z}_{t}}^{2}\\!+\\!Q_{S^{z}_{t}}}}\hskip
8.5359pt}.$ (6)

To study the bipartite entanglement we will use the concept of negativity [16,
17, 18], which is defined in terms of negative eigenvalues $\lambda_{i}$ of
the partially transposed density matrix $\hat{\rho}^{T_{1/2}}$, ${\cal
N}\\!=\\!\sum_{\lambda_{i}<0}|\lambda_{i}|$. The zero value of the negativity
corresponds to separable (disentangled) states, whereas its non-zero value
corresponds to inseparable (entangled) ones. For MSHD, the negativity of the
maximally entangled state equals to one-half (${\cal N}\\!=\\!1/2$). Detailed
derivations of the negativity of MSHD can readers find in Ref. [15].

##  3 Results and discussion

In the following, our discussion will be focused on the physically most
interesting antiferromagnetic case ($J\\!>\\!0$), for which a great variety of
ground states is expected. Fig. 1 illustrates the ground-state phase diagram
in the $D/J\\!-\\!\Delta$ plane commonly with the degree of entanglement
exemplified through the respective zero-temperature density plot of the
negativity. One immediately identifies that MSHDs with both half-odd-integer
constituents (right panels) remain maximally entangled if the non-degenerate
antiferromagnetic ground state $|(0)_{-}\rangle$ is favoured. A non-degenerate
character of a corresponding energy level $\varepsilon_{0}^{-}$ determines a
simple structure of the negativity, ${\cal N}\\!=\\!c_{0}^{-}c_{0}^{+}$, which
results in the negativity ${\cal N}\\!=\\!1/2$ completely independent of the
strength of both assumed anisotropies.

Figure 1:  Density plots of the negativity in the plane $D/J$ vs. $\Delta$ for
integer (left panels) and half-odd-integer (right panels) spin $S$, which
correspond to different ground states $|(\pm S^{z}_{t})_{-}\rangle$.

It is noteworthy that the degree of entanglement changes discontinuously at
the relevant phase boundary of the $|(0)_{-}\rangle$ ground state and then it
continuously decreases with decreasing $\Delta$ and/or $D/J$. For half-odd-
integer spins $S$ and the specific exchange anisotropy $\Delta$ one identifies
$S\\!-\\!1/2$ phase boundaries, which merge together for a fully isotropic
case ($\Delta\\!=\\!1$ and $D/J\\!=\\!0$) as a consequence of energy
equivalence between all $|(S^{z}_{t})_{-}\rangle$ ground states
($|S^{z}_{t}|\\!\leq\\!S\\!-\\!1/2$). On the other hand, the two-fold
degeneracy of each $|(S^{z}_{t})_{-}\rangle$ ground state
($|S^{z}_{t}|\\!\leq\\!S\\!-\\!1/2$) for MSHDs with integer spins $S$ (left
panels) determines the different behaviour of the negativity in the most
entangled ground state $|(\pm 1/2)_{-}\rangle$. Here
$\varepsilon^{-}_{1/2}\\!=\\!\varepsilon^{-}_{-1/2}$ and probability
amplitudes $c_{1/2}^{\pm}\\!=\\!c_{-1/2}^{\mp}$. Because
$S^{z}_{t}\\!=\\!1/2$, the probability amplitude
$c_{S^{z}_{t}-1}^{\pm}\\!=\\!c_{1/2}^{\mp}$ and then ${\cal
N}\\!=\\!c_{1/2}^{-}\left(c_{1/2}^{-}\\!-\\!\sqrt{(c_{1/2}^{-})^{2}\\!+\\!(c_{1/2}^{-}c_{1/2}^{+})^{2}}\right)/4$.
Using definitions (4)-(6) one obtains the relation

$\displaystyle\allowdisplaybreaks{\cal
N}\\!=\\!\frac{1}{4}\left(1\\!-\\!\frac{1\\!-\\!2D/J}{\sqrt{\alpha}}\right)\left(\sqrt{\frac{5\sqrt{\alpha}+3(1\\!-\\!2D/J)}{\sqrt{\alpha}-(1\\!-\\!2D/J)}}\\!-\\!1\right),$
(7)

where $\alpha\\!=\\!(1\\!-\\!2D/J)^{2}\\!+\\!4\Delta^{2}S(S\\!+\\!1)$.

Figure 2: The dependence of a quantum negativity as a function of $D/J$ for
different integer spin-$S$. Three exchange anisotropy limits are assumed. Open
circles visualised in panel $(b)$ denote the negativity at fully isotropic
case, $\Delta\\!=\\!1$ and $D/J\\!=\\!0$.

It is clear that in the region of $|(\pm 1/2)_{-}\rangle$ the degree of
negativity strongly depends on the value of $\Delta$, $D/J$ and spin
magnitude. Detailed numerical analysis (Fig. 2) exhibits existence of
invariant negativity point at $D/J\\!=\\!1/2$ below which the negativity of
dimer with a higher integer spin $S$ gradually decreases, whereas above
$D/J\\!=\\!1/2$ a pronounced enhancement of entanglement is observed. In Fig.
2 one may additionally identify a singularity at a fully isotropic case
($\Delta\\!=\\!1$, $D/J\\!=\\!0$), which originates from the $2S$-fold
degeneracy of $|(\pm S^{z}_{t})_{-}\rangle$ eigenstates. The negativity at
this point rapidly falls down with increasing spin value (open circles in Fig.
2$(b)$) in accordance with the identity (30) in Ref. [15]. Because the dimer
with $S\\!=\\!1$ can be arranged only in two-fold degenerate $|(\pm
1/2)_{-}\rangle$ ground state, the respective negativity in a fully isotropic
case is a smooth function with absence of an evident sharp minimum. Due to the
same reason the negativity of MSHD in the region of
$|(\pm(S\\!-\\!1/2))_{-}\rangle$ behaves differently than the negativity of
MSHD with higher spins $S\\!>\\!1$. Focusing on the negativity at
$|(\pm(S\\!-\\!1/2))_{-}\rangle$ phase, we confirm furthermore a previously
known fact derived for $D/J\\!=\\!0$ [8, 11, 12] that the increasing spin
value $S$ ($S\\!>\\!1$) of MSHD leads to the reduction of a degree of
entanglement. The explanation of a different behaviour for the negativity in
$|(\pm 1/2)_{-}\rangle$ and $|(\pm(S\\!-\\!1/2))_{-}\rangle$ phase can be
found in an analytical expression of respective negativity functions. For
$|(\pm(S\\!-\\!1/2))_{-}\rangle$ the negativity fulfils the relation ${\cal
N}\\!=\\!c_{(S-1/2)}^{-}c_{(S-1/2)}^{+}\\!=\\!(1/2)\\{(8S\Delta^{2})/(8S\Delta^{2}\\!+\\!(1\\!-\\!2D/J)^{2}(2S\\!-\\!1)^{2})\\}^{1/2}$,
which is evidently inversely proportional to the spin magnitude $S$. Thus, the
negativity decreases as the spin value $S$ enlarges.

Figure 3:  Thermal behaviour of negativity for different spin $S$ in the mixed
spin-(1/2,$S$) dimer for various $D/J$ anisotropy.

Analysing the influence of thermal fluctuations on the system entanglement we
detect a few interesting observations collected in Fig. 3. It should be
emphasized, that Fig. 3 involves $\Delta\\!=\\!1$ for a simplicity, but all
conclusions are qualitatively valid for an arbitrary $\Delta$. (i) The
increasing spin $S$ values of MSHD always enlarges the threshold temperature
above which the entanglement disappears; (ii) If the interplay between quantum
fluctuations favours the $|(\pm(S\\!-\\!1/2))_{-}\rangle$ ground state (Fig.
3$(a)$), the negativity gradually decreases upon increasing spin size $S$ with
exception of the specific case $S\\!=\\!1$; (iii) In the fully isotropic limit
(Fig. 3$(b)$), where all $|(S^{z}_{t})_{-}\rangle$ states have an identical
energy, the thermal negativity smoothly decreases upon temperature rise.
Similarly as in the previous case, the higher spin $S$ is responsible for a
higher thermal stability and lower degree of bipartite entanglement.
Nevertheless, the highest negativity in the fully isotropic case is
significantly smaller in comparison to its anisotropic counterpart due to a
degeneracy of the ground state discussed above; (iv) The reduction effect of
an increasing spin magnitude $S$ on the degree of thermal entanglement is
observed equally in case of $|(\pm 1/2)_{-}\rangle$ or $|(0)_{-}\rangle$
ground state if and only if $D/J\\!<\\!1/2$ (Fig. 3$(c)$). In accordance to
the zero-temperature analysis, the largest negativity of MSHDs with an
arbitrary half-odd-integer spin $S$ is ${\cal N}\\!=\\!1/2$, whereas the
negativity of MSHD with an integer spin $S$ saturates to a smaller value
depending on both anisotropies; (iv) The stability as well as the degree of
thermal entanglement of MSHD with $D/J\\!>\\!1/2$ can be enhanced by
increasing the spin size $S$. As illustrated in Fig. 3$(d)$ the thermal
entanglement can be additionally stabilized by half-odd-integer spins $S$, for
which the low-temperature asymptotic limit of the negativity achieves the
maximal possible value 1/2. The increasing degree of thermal entanglement
driven by an increasing spin $S$ magnitude is also present in MSHDs with
integer spin value $S$, however, the largest negativity does not exceed a half
of the golden ratio $(\\!\sqrt{5}\\!-\\!1)/4$.

Finally, let us look at the behaviour of a threshold temperature upon the
variation of $D/J$ (Fig. 4). For a better lucidity the half-odd-integer spins
$S$ are assumed only. As a result, the increasing spin value $S$ as well as
the increasing $\Delta$ enlarges the threshold temperature between the
entangled and separable state. It is evident that the single-ion anisotropy
$D/J$ has an additional enlarging effect on the threshold temperature with a
significant sharp minimum at a proximity of $D/J\\!=\\!0$. While the threshold
temperature for $D/J$ lying below the local minimum is a smooth continuous
function, the threshold temperature for $D/J$ above the local minimum involves
a few pronounced flections, whose number increases as the spin magnitude $S$
increases. It is supposed that the origin of such different behaviour relates
to the different ground-state spin arrangements, where the increasing
temperature affecting the starting (almost) antiferromagnetic structure can
temporarily favour one of other $|(S^{z}_{t})_{-}\rangle$ ones, until the
disentangled state is achieved.

Figure 4:  The behaviour of the threshold temperature of the model (1) for
different spin $S$ and three values of $\Delta$.

In conclusion, it was identified that the integer or half-odd-integer
character of the spin $S$ significantly determines the bipartite entanglement
both, at zero as well as non-zero temperatures. The half-odd-integer character
of both spins allows one to stabilize the maximally entangled
antiferromagnetic ground state, for which the negativity is completely kept
constant. On the other hand, the negativity of MSHDs with integer spins $S$ is
not higher than ${\cal N}\\!=\\!(\sqrt{5}\\!-\\!1)/4$ and it can be tuned by
magnetic anisotropies. Moreover, if $D/J\\!>\\!1/2$, the negativity can be
enhanced by increasing of the spin $S$ magnitudes though it does not reach the
maximal value ${\cal N}\\!=\\!1/2$. In addition, the negativity in the fully
isotropic MSHD ($\Delta\\!=\\!1$, $D/J\\!=\\!0$) is significantly smaller in
comparison to its anisotropic counterpart as a consequence of $2S$-fold
degeneracy. The increasing spin value $S$ rapidly reduces its maximal value
from ${\cal N}\\!=\\!1/3$ achieved for the particular case with $S\\!=\\!1$.
In contradiction to the knowledge for $D/J\\!=\\!0$, the increasing spin
magnitude $S$ in region of $D/J\\!>\\!1/2$ can be utilized to enhance not only
the thermal stability but also the degree of thermal entanglement of MSHD.
From the application point of view, MSHD with half-odd-integer spin $S$ is
more favourable due to its saturation to a maximally entangled state. At the
same time, the increasing spin $S$ additionally enlarges the threshold
temperature between the entangled and disentangled state, which enhancement
can be furthermore tuned by uniaxial single-ion anisotropy (of the easy-axis
as well as the easy-plan type).

This work was financially supported by the grant Nos. APVV-16-0186 and VEGA
1/0105/20.

## References

  * [1] M.A. Nielsen, I.L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, 2000.
  * [2] F. Wang, H. Jia, H. Zhang, X. Zhang, S. Chang, Sci. China Ser G. 52 (2009) 1919.
  * [3] K.T. Guo, M.Ch. Liang, H.Y. Xu, Ch.B. Zhu, J. Phys. A: Math. Theor. 43 (2010) 505301.
  * [4] S. Xu, X. Song, L. Ye, Quantum Inf. Process 13 (2014) 1013.
  * [5] K.T. Guo, S.H. Xiang, H.Y. Xu, X.H. Li, Quantum Inf. Process 13 (2014) 1511.
  * [6] Ch.B. Zhou, S.Y. Xiao, C. Zhang, G. Wu, Y.-Q. Ran, Physica B 477 (2015) 40.
  * [7] L. Zhou, H. S. Sonf, Y.Q. Guo, C. Li, Phys. Rev. A 68 (2003) 02430.
  * [8] S.S. Li, T.Q. Ren, X.M. Kong, K. Liu, Physica A 391 (2012) 35.
  * [9] Ch.B. Zhou, S.Y. Xiao, C. Zhang, G. Wu, Y.Q. Ran, Int. J. Theor. Phys. 55 (2016) 875.
  * [10] H. Čenčariková, J. Strečka, Phys. Rev. B 102 (2010) 184419.
  * [11] H. Huang, X. Wang, Z. Sun, G. Yang, Physica A 387 (2008) 2736.
  * [12] X. Hao and S. Zhu, Phys. Lett. A 366 (2007) 206.
  * [13] E. Solano-Carrilo, R. Franco, J. Silva-Valencia, Phys. Lett. A 375 (2011) 1032.
  * [14] E. Solano-Carrilo, R. Franco, J. Silva-Valencia, Phys. Lett. A 390 (2011) 2208.
  * [15] H. Vargová, J. Strečka, arXiv:2107.05628v1.
  * [16] A. Peres, Phys. Rev. Lett. 77 (1996) 1413.
  * [17] M. Horodecki, P. Horodecki, R. Horodecki, Phys. Lett. A 223 (1996) 1.
  * [18] G. Vidal, R.F. Werner, Phys. Rev. A 65 (2002) 032314.