# Strong coupling in chiral cavities:
nonperturbative framework for enantiomer discrimination

Rosario R. Riso1 Laura Grazioli2 Enrico Ronca3 Tommaso Giovannini4 Henrik
Koch1,4<EMAIL_ADDRESS>1Department of Chemistry, Norwegian University of
Science and Technology, 7491 Trondheim, Norway  2Institut für Physikalische
Chemie, Universität Mainz, D-55099 Mainz, Germany  3Istituto per i Processi
Chimico Fisici del CNR (IPCF-CNR), Via G. Moruzzi, 1, 56124, Pisa, Italy
4Scuola Normale Superiore, Piazza dei Cavalieri 7, 56126 Pisa, Italy

###### Abstract

Development of efficient techniques to distinguish mirror images of chiral
molecules (enantiomers) is very important in both chemistry and physics.
Enantiomers share most molecular properties except for instance different
absorption of circularly polarized light. Enantiomer purification is therefore
a challenging task that requires specialized equipment. Strong coupling
between quantized fields and matter (e.g. in optical cavities) is a promising
technique to modify molecular properties in a non-invasive way. The modulation
of molecular features is achieved by changing the field characteristics. In
this work, we investigate whether strong coupling to circularly polarized
electromagnetic fields is a viable way to discriminate chiral molecules. To
this end, we develop a nonperturbative framework to calculate the behaviour of
molecules in chiral cavities. We show that in this setting the enantiomers
have different energies – that is, one being more stable than the other. The
field-induced energy differences are also shown to give rise to
enantiospecific signatures in rotational spectra.

††preprint: APS/123-QED

##  I Introduction

The study of strongly coupled light-matter systems is becoming a well
established area in both physics and chemistry[1, 2, 3, 4, 5, 6]. Through the
interaction with the quantized field, it is indeed possible to modify both
macroscopic material features, for example conductivity or phase transitions
[7, 8], or to influence microscopic properties like the chemistry of single
molecules [9, 10, 11]. The easiest way to reach the strong coupling regime is
through optical cavities, composed of two highly reflective mirrors placed in
front of each other [12]. Because of the mirrors, the electromagnetic fields
inside the cavity must fulfill specific boundary conditions. The modulation of
the boundary conditions allows for a fine tuning of the main field properties.
For example, choosing a certain geometry of the device, specific sets of
frequencies can be selected. On the other hand, an opportune engineering of
the mirrors can lead to changes in the shape of the field or selection of its
polarization [13, 14, 15]. The confinement of photons in a small quantization
volume leads to an increase in the light-matter interaction [16, 17] with the
consequent formation of hybrid states called polaritons [18]. Polaritons
represent a very effective way to modulate matter properties in a non-invasive
way. In fact, the characteristics of the mixed electron-photon states can be
manipulated by tuning the field properties [19, 20]. In recent years,
experimental developments have dramatically improved our control of the cavity
field, opening the way to a wide range of new applications [14, 16, 13, 21,
22, 23, 24, 25, 26, 27, 28, 29]. Theoretical modeling of strong coupling
phenomena is urgently needed to advance our intuitive understanding of what
happens inside cavities and to assist the experiment-design phase. Many
remarkable efforts have already been presented in the literature, with the
introduction of both ab initio electronic structure methods [30, 31, 32, 33,
34, 35, 22] and molecular dynamics techniques [36, 37] for strongly coupled
systems. Nonetheless, a comprehensive framework for arbitrary field shapes is
still not available and requires further theoretical developments.

One of the most intriguing perspectives in molecular polaritonics is the
possibility of enhancing the spectroscopic techniques through the interaction
with quantized fields. Pioneering works have indeed shown that upon coupling
with the quantized fields it is possible to increase the spectral resolution
amplifying the signal intensities [38, 39], even at the level of single
molecule imaging [40, 41, 42]. In this regard, a particularly interesting
perspective is to use strong coupling to discriminate among chiral molecules
(systems that are non-superimposable with their mirror image). The two
different configurations of a chiral molecule, called enantiomers, share most
physical properties and can be distinguished when they interact with
circularly polarized light [43, 44]. In particular, chiral molecules interact
and absorb left-handed circular polarization (LHCP) and right-handed circular
polarization (RHCP) light differently [45, 46]. This suggests that in chiral
cavities [47, 14], devices where the electromagnetic field has a fixed
circular polarization, it might be possible to discriminate among chiral
molecules and even create energy differences between two enantiomers.

In this work, we present a theoretical framework to study the strong coupling
regime in chiral cavities. We first discuss the formal properties that a
cavity quantum electrodynamical (cQED) model must obey. Afterwards, we develop
a novel ab initio method that enables the study of such systems. We
demonstrate that, even when no real photons are present in the cavity,
coupling to a circularly polarized field induces energy differences between
the enantiomers. Finally, we show that the energy differences produce
enantiospecific signatures in molecular rotational spectra.

##  II Chiral cavities and the single-mode approximation

The interaction between photons and matter will be described using the minimal
coupling Hamiltonian [48, 49]

$\begin{split}H=&\frac{1}{2}\sum_{i}\left(\mathbf{p}_{i}-\mathbf{A}(\mathbf{r}_{i})\right)^{2}+\sum_{I}\frac{1}{2M_{I}}\left(\mathbf{p}_{I}+Z_{I}\mathbf{A}(\mathbf{R}_{I})\right)^{2}\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|R_{I}-r_{i}\right|}\\\
+&\frac{1}{2}\int\left(\mathbf{E}^{2}(\mathbf{r})+c^{2}\mathbf{B}^{2}(\mathbf{r})\right)d^{3}r,\end{split}$
(1)

where $i$ and $j$ label electrons while $I$ and $J$ label nuclei with charges
$Z_{I}$ and $Z_{J}$. The vector potential, the electric and the magnetic
fields are denoted by $\mathbf{A}(\mathbf{r})$, $\mathbf{E}(\mathbf{r})$ and
$\mathbf{B}(\mathbf{r})$, respectively. Since in the strong coupling regime
the electromagnetic field is a critical component of the system, it must be
treated on the same footing as the electrons by means of quantum
electrodynamics (QED) [46]. The optical properties of a cavity are encoded in
the vector potential $\mathbf{A}(\mathbf{r})$ (and consequently in
$\mathbf{E}(\mathbf{r})$ and $\mathbf{B}(\mathbf{r})$). Specifically, in a
chiral cavity only one of the two possible circular polarizations of the field
is allowed. The second quantized vector potential can therefore be written as

$\mathbf{A}(\mathbf{r})=\sum_{k}\frac{\lambda}{\sqrt{2\omega_{k}}}\left(\bm{\epsilon}_{k\pm}b_{k}e^{i\mathbf{k}\mathbf{r}}+\bm{\epsilon}^{*}_{k\pm}b^{\dagger}_{k}e^{-i\mathbf{k}\mathbf{r}}\right),$
(2)

where $k$ spans all the possible wave vectors and $\bm{\epsilon}_{k\pm}$
denotes the field polarization. In Eq. (2), $\omega_{k}$ is the frequency of
the field while $\lambda$ quantifies the strength of the light-matter coupling
and is related to the quantization volume $\lambda\propto\frac{1}{\sqrt{V}}$ .
In this paper we will assume, without loss of generality, that $\mathbf{k}$ is
aligned along the z-axis and that the field polarization is

$\bm{\epsilon}_{k\pm}=\bm{\epsilon}_{\pm}=\frac{1}{\sqrt{2}}\begin{pmatrix}1\\\
\pm i\\\ 0\end{pmatrix}.$ (3)

In the Born-Oppenheimer approximation, Eq. (1) takes the form

$\begin{split}H=&\frac{1}{2}\sum_{i}\left(\mathbf{p}_{i}-\sum_{k}\frac{\lambda}{\sqrt{2\omega_{k}}}\left(\bm{\epsilon}_{\pm}b_{k}e^{i\mathbf{k}\mathbf{r}_{i}}+\bm{\epsilon}^{*}_{\pm}b^{\dagger}_{k}e^{-i\mathbf{k}\mathbf{r}_{i}}\right)\right)^{2}\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|R_{I}-r_{i}\right|}\\\
+&\sum_{k}\omega_{k}\left(b^{\dagger}_{k}b_{k}+\frac{1}{2}\right),\end{split}$
(4)

where we have used that the nuclear mass $M_{I}>c$ in atomic units and the
nuclei are kept fixed. Since the full Hamiltonian in Eq. (4) involves an
infinite number of modes, it is computationally unfeasible to determine the
eigenvalues and eigenfunctions. However, we expect that restricting Eq. (4) to
include only one mode will be sufficient to obtain a qualitatively correct
description of the field effects. This restriction is known as the single-mode
approximation and is widely used together with the dipole approximation
($e^{\pm i\mathbf{k}\mathbf{r}}=1$) [50, 51, 52]. In the single-mode
approximation, the Hamiltonian in Eq. (4) becomes

$\begin{split}H_{k}=&\frac{1}{2}\sum_{i}\left(\mathbf{p}_{i}-\frac{\lambda}{\sqrt{2\omega_{k}}}\left(\bm{\epsilon}_{\pm}b_{k}e^{i\mathbf{k}\mathbf{r}_{i}}+\bm{\epsilon}^{*}_{\pm}b^{\dagger}_{k}e^{-i\mathbf{k}\mathbf{r}_{i}}\right)\right)^{2}\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|R_{I}-r_{i}\right|}\\\
+&\omega_{k}\left(b^{\dagger}_{k}b_{k}+\frac{1}{2}\right).\end{split}$ (5)

Although this Hamiltonian seems to be theoretically consistent with Eq. (4),
the symmetry properties are not the same. In particular, while a cavity
rotation of $180$ degrees around an axis perpendicular to $\mathbf{k}$ (
$W(\pi)$ ) leaves Eq. (4) unchanged, the single-mode Hamiltonian is
transformed into

$\begin{split}\tilde{H}_{k}=&\;W^{-1}(\pi)H_{k}W(\pi)\\\
=&\;\frac{1}{2}\sum_{i}\left(\mathbf{p}_{i}-\frac{\lambda}{\sqrt{2\omega_{k}}}\left(\bm{\epsilon}_{\mp}b_{-k}e^{-ik\mathbf{r}_{i}}+\bm{\epsilon}^{*}_{\mp}b^{\dagger}_{-k}e^{i\mathbf{k}\mathbf{r}_{i}}\right)\right)^{2}\\\
+&\sum_{i,j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|r_{i}-R_{I}\right|}\\\
+&\omega_{k}\left(b^{\dagger}_{-k}b_{-k}+\frac{1}{2}\right).\end{split}$ (6)

The Hamiltonians in Eqs. (5) and (6) are different because the bosonic
operators are multiplied with different electronic operators. This
inconsistency suggests that the single-mode approximation cannot be used for
chiral cavities. When rotating the cavity by $180$ degrees the light-matter
operators transform as

$W^{-1}(\pi)\left[\bm{\epsilon}_{\pm}b_{k}e^{i\mathbf{k}\mathbf{r}_{i}}\right]W(\pi)\longrightarrow\bm{\epsilon}_{\mp}b_{-k}e^{-i\mathbf{k}\mathbf{r}_{i}},$
(7)

where the wave vector has been reversed and the polarization vector has
changed sign. We point out that since the transformation in Eq. (7) does not
involve reflections, the circular polarization of the field is unchanged. In
order to be invariant under $W(\pi)$, the vector potential must include
contributions from both $b_{k}$ and $b_{-k}$, that is

$\mathbf{A}(\mathbf{r})=\frac{\lambda}{\sqrt{2\omega_{k}}}\left(\bm{\epsilon}_{\pm}(b_{k}+b^{\dagger}_{-k})e^{i\mathbf{k}\mathbf{r}}+\bm{\epsilon}^{*}_{\pm}(b^{\dagger}_{k}+b_{-k})e^{-i\mathbf{k}\mathbf{r}}\right).$
(8)

Therefore, at least two different modes of the electromagnetic field are
needed to characterize a chiral cavity.

Figure 1: In a chiral cavity (left side of the picture) the left and right
moving photons are different even in dipole approximation (the two arrows
rotate in a different way) and both the modes must therefore be included in
the Hamiltonian. In a linearly polarized cavity, instead, once the dipole
approximation has been adopted the contribution from the two photons is
exactly the same. Therefore a one mode picture is enough to describe the
system.

Specifically, such modes should have the same frequency but opposite wave
vector. Using the definition in Eq. (8) the Hamiltonian for a molecular system
in a chiral cavity (Eq. (4)) is

$\begin{split}H=&\sum_{i}\frac{\mathbf{p}^{2}_{i}}{2}+\frac{N_{e}\lambda^{2}}{2\omega}\left(b_{k}+b^{\dagger}_{-k}\right)\left(b_{-k}+b^{\dagger}_{k}\right)\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|r_{i}-R_{I}\right|}\\\
+&\frac{\lambda}{\sqrt{2\omega}}\sum_{i}\left(\mathbf{p}_{i}\cdot\bm{\epsilon}_{\pm}\right)(b_{k}+b^{\dagger}_{-k})e^{i\mathbf{k}\mathbf{r}_{i}}\\\
+&\frac{\lambda}{\sqrt{2\omega}}\sum_{i}(\mathbf{p}_{i}\cdot\bm{\epsilon^{*}}_{\pm})(b^{\dagger}_{k}+b_{-k})e^{-i\mathbf{k}\mathbf{r}_{i}}\\\
+&\omega\left(b^{\dagger}_{k}b_{k}+b^{\dagger}_{-k}b_{-k}+1\right),\end{split}$
(9)

where $N_{e}$ is the number of electrons in the system and and $\omega$ is
used to denote $\omega_{k}$.
So far, we have shown that two modes are necessary to describe the chiral
field consistently. A similar line of arguments shows that the two-mode
treatment is also required when the dipole approximation is applied in Eq.
(4). However, the single-mode Hamiltonian has been extensively used with great
success for a linearly polarized field (real $\bm{\epsilon})$ in the dipole
approximation [53, 54, 55]. In this case, the two-mode Hamiltonian equals

$\begin{split}H=&\sum_{i}\frac{\mathbf{p}^{2}_{i}}{2}+\frac{N_{e}\lambda^{2}}{4\omega}\left(b_{k}+b^{\dagger}_{k}+b_{-k}+b^{\dagger}_{-k}\right)^{2}\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|r_{i}-R_{I}\right|}\\\
+&\frac{\lambda}{\sqrt{2\omega}}\sum_{i}(\mathbf{p}_{i}\cdot\bm{\epsilon})(b_{k}+b^{\dagger}_{-k})\\\
+&\frac{\lambda}{\sqrt{2\omega}}\sum_{i}\left(\mathbf{p}_{i}\cdot\bm{\epsilon}\right)(b^{\dagger}_{k}+b_{-k})\\\
+&\omega\left(b^{\dagger}_{k}b_{k}+b^{\dagger}_{-k}b_{-k}+1\right)\end{split}$
(10)

that can be rewritten as

$\begin{split}H=&\sum_{i}\frac{\mathbf{p}^{2}_{i}}{2}+\frac{N_{e}\lambda^{2}}{2\omega}\left(b+b^{\dagger}\right)^{2}+\frac{\lambda}{\sqrt{\omega}}\sum_{i}\left(\mathbf{p}_{i}\cdot\bm{\epsilon}\right)(b+b^{\dagger})\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{j,I}\frac{Z_{I}}{\left|r_{j}-R_{I}\right|}\\\
+&\omega\left(b^{\dagger}b+\tilde{b}^{\dagger}\tilde{b}+1\right),\end{split}$
(11)

where we have introduced the two symmetry-adapted boson operators

$\begin{pmatrix}b^{\dagger}\\\ \tilde{b}^{\dagger}\\\
\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}\;\;1&\;\;1\\\ -1&\;\;1\\\
\end{pmatrix}\begin{pmatrix}b^{\dagger}_{k}\\\ b^{\dagger}_{-k}\\\
\end{pmatrix}\\\ $ (12)

such that $W^{-1}(\pi)b^{\dagger}W(\pi)=b^{\dagger}$ and
$W^{-1}(\pi)\tilde{b}^{\dagger}W(\pi)=-\tilde{b}^{\dagger}$. Equation (11) is
the standard single-mode Hamiltonian with the coupling $\lambda$ scaled by a
factor $\sqrt{2}$, plus an additional $\omega\tilde{b}^{\dagger}\tilde{b}$
term. Therefore, as expected, for Fabry-Pérot cavities (where the field
polarization is real) the single-mode approximation (11) can be used if the
dipole approximation has been adopted. A pictorial representation of the
results discussed in this section is reported in Fig. 1. In particular, we
highlight that, when the field is circularly polarized, the photons moving to
the right are different from those moving to the left even in the dipole
approximation. Therefore, both modes need to be included explicitly in our
description of the vector potential. In the case of a linearly polarized
field, instead, the left and right moving photons can only be distinguished if
we go beyond the dipole approximation.

##  III Formal properties of the chiral cavity Hamiltonian

Figure 2: The L enantiomer of proline in a cavity with a LHCP field has the
same energy as the R enantiomer in a cavity with a RHCP. However, no symmetry
in the Hamiltonian ensures that the L enantiomer has the same energy in the
LHCP cavity as in a RHCP polarized cavity.

From the properties of the Hamiltonian in Eq. (9), we can infer the
characteristic features of the chiral cavity. Upon reflection of the entire
system (molecule plus cavity), the molecule transforms into its mirror image
while the cavity polarization is inverted ( LHCP$\leftrightarrow$RHCP). The
energy of the system remains the same because there is no parity-violating
interaction in the Hamiltonian [56]. However, while a non-chiral molecule is
indistinguishable from its mirror image, performing the reflection of a chiral
system will produce a different enantiomer, non-superimposable with the
original molecule. Therefore non-chiral molecules have the same energy inside
two differently circularly polarized cavities. On the other hand, a chiral
molecule in a chiral cavity has the same energy as the other enantiomer in a
cavity with the opposite circular polarization. Most importantly, there is no
symmetry of the Hamiltonian that requires the two enantiomers to have the same
energy in the same chiral cavity. Strong coupling between circularly polarized
light and chiral molecules might therefore be a viable way to discriminate
enantiomers, as shown in Fig. 2.
The discriminating properties are lost when the dipole approximation is
adopted, even if the field polarization remains chiral. This can be verified
by setting $e^{i\mathbf{k}\mathbf{r}}$ equal to one in Eq. (9). In this case,
the Hamiltonians for the LHCP and RHCP cavities have the same eigenvalues, as
one Hamiltonian can be transformed into the other through a relabeling of the
wave vectors ($\mathbf{k}\leftrightarrow-\mathbf{k}$). Thus, the exponential
$e^{i\mathbf{k}\mathbf{r}}$ in the field parameterization plays a critical
role for the discriminating power of the cavity.
We will employ the Hamiltonian in Eq. (9) to study a chiral molecule in a
chiral cavity. The eigenvalue problem of this Hamiltonian becomes easier if we
remove the quadratic term in the field. This can be accomplished carrying out
a Bogoljubov transformation [57] that introduces two new bosonic operators

$\begin{split}&\begin{pmatrix}b_{k}\\\ b^{\dagger}_{-k}\\\
\end{pmatrix}=\begin{pmatrix}\cosh{\theta}&-\sinh{\theta}\\\
-\sinh{\theta}&\cosh{\theta}\\\ \end{pmatrix}\begin{pmatrix}\alpha\\\
\beta^{\dagger}\\\ \end{pmatrix}\\\ &\begin{pmatrix}\alpha\\\
\beta^{\dagger}\\\ \end{pmatrix}=\begin{pmatrix}\cosh{\theta}&\sinh{\theta}\\\
\sinh{\theta}&\cosh{\theta}\\\ \end{pmatrix}\begin{pmatrix}b_{k}\\\
b^{\dagger}_{-k}\\\ \end{pmatrix},\end{split}$ (13)

that still satisfy the bosonic commutation relationships

$\begin{split}[\alpha,\alpha^{\dagger}]=1\hskip
28.45274pt[\alpha,\beta]=0\hskip
28.45274pt[\alpha,\beta^{\dagger}]=0.\end{split}$ (14)

If $\tanh{2\theta}=\frac{N_{e}\lambda^{2}}{N_{e}\lambda^{2}+2\omega^{2}}$ in
Eq. (13), the Hamiltonian can be rewritten as

$\begin{split}H=&\sum_{i}\frac{\mathbf{p}^{2}_{i}}{2}+\tilde{\omega}\left(\alpha^{\dagger}\alpha+\beta^{\dagger}\beta+1\right)\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{jI}\frac{Z_{I}}{\left|r_{j}-R_{I}\right|}\\\
+&\frac{\lambda}{\sqrt{2\tilde{\omega}}}\sum_{i}\left(\mathbf{p}_{i}\cdot\bm{\epsilon}_{\pm}\right)(\alpha+\beta^{\dagger})e^{i\mathbf{k}\mathbf{r}_{i}}\\\
+&\frac{\lambda}{\sqrt{2\tilde{\omega}}}\sum_{i}(\mathbf{p}_{i}\cdot\bm{\epsilon}^{*}_{\pm})(\alpha^{\dagger}+\beta)e^{-i\mathbf{k}\mathbf{r}_{i}},\end{split}$
(15)

where $\tilde{\omega}=\sqrt{\omega^{2}+N_{e}\lambda^{2}}$ and the
$\mathbf{A}^{2}$ in Eq. (9) term has been reabsorbed [57]. Note that, since
the transformation in Eq. (13) can also be expressed as a unitary
transformation with the operator

$U_{k}=\textrm{exp}\left[\theta_{k}\left(b^{\dagger}_{k}b^{\dagger}_{-k}-b_{k}b_{-k}\right)\right],$
(16)

the Hamiltonians in Eqs. (15) and (9) have the same eigenvalues. The form in
Eq. (15) is particularly interesting because the zero-point energy term
explicitly depends on the number of electrons. This result is relevant for two
main reasons:

  1. (i)

When the cavity frequency is approaching zero, numerical methods based on the
minimal coupling approach in Eq. (9) become unstable because the light-matter
coupling term diverges as $\omega^{-\frac{1}{2}}$ while the quadratic term in
the field diverges as $\omega^{-1}$. The Hamiltonian in Eq. (15) does not
diverge for $\tilde{\omega}=0$;

  2. (ii)

The Hamiltonian in Eq. (15) is explicitly non size-extensive due to the
contributions from $\tilde{\omega}$, which is non linear in the number of
electrons. We note that size-extensivity is recovered in the limit where
$\lambda=0$.

In Appendix A we show that
$\tanh{2\theta}=\frac{N_{e}\lambda^{2}}{N_{e}\lambda^{2}+2\omega^{2}}$
minimizes the zero-point energy and that the non size-extensive effects
obtained from the zero-point energy contribution in Eq. (15) have the correct
dependence on the number of electrons when $\lambda/\omega<<1$. The
transformation in Eq. (16) is also useful when more modes are included in the
minimal coupling Hamiltonian

$\displaystyle H=$
$\displaystyle\sum_{i}\frac{\mathbf{p}^{2}_{i}}{2}+\sum_{k>0}\omega_{k}\left(b^{\dagger}_{k}b_{k}+b^{\dagger}_{-k}b_{-k}+1\right)$
$\displaystyle+$
$\displaystyle\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|R_{I}-r_{i}\right|}$
$\displaystyle+$
$\displaystyle\lambda\sum_{i}\sum_{k>0}\left(\mathbf{p}_{i}\cdot\bm{\epsilon}_{\pm}\right)e^{i\mathbf{k}\mathbf{r}_{i}}\frac{(b_{k}+b^{\dagger}_{-k})}{\sqrt{2\omega_{k}}}$
$\displaystyle+$
$\displaystyle\lambda\sum_{i}\sum_{k>0}(\mathbf{p}_{i}\cdot\bm{\epsilon^{*}}_{\pm})e^{-i\mathbf{k}\mathbf{r}_{i}}\frac{(b^{\dagger}_{k}+b_{-k})}{\sqrt{2\omega_{k}}}$
(17) $\displaystyle+$
$\displaystyle\sum_{k>0}\frac{N_{e}\lambda^{2}}{2\omega_{k}}\left(b_{k}+b^{\dagger}_{-k}\right)\left(b_{-k}+b^{\dagger}_{k}\right)$
$\displaystyle+$ $\displaystyle\sum_{ik\neq
k^{\prime}>0}\frac{N_{e}\lambda^{2}e^{i(\mathbf{k}-\mathbf{k}^{\prime})\mathbf{r}_{i}}}{2\sqrt{\omega_{k}\omega_{k^{\prime}}}}\left(b_{k}+b^{\dagger}_{-k}\right)\left(b_{-k^{\prime}}+b^{\dagger}_{k^{\prime}}\right).$

In this case the squared term couples different field modes. Due to the
presence of electronic operators, the quadratic contributions in the field
cannot be fully reabsorbed by using a unitary bosonic transformation similar
to Eq. (16). However, a product of those rotations leads to

$U=\prod_{k>0}\textrm{exp}\left[\theta_{k}\left(b^{\dagger}_{k}b^{\dagger}_{-k}-b_{k}b_{-k}\right)\right],$
(18)

which can be used to reabsorb the purely photonic terms of Eq. (17) leading to

$\begin{split}H=&\sum_{i}\frac{\mathbf{p}^{2}_{i}}{2}+\sum_{k}\tilde{\omega}_{k}\left(\alpha^{\dagger}_{k}\alpha_{k}+\beta^{\dagger}_{k}\beta_{k}+1\right)\\\
+&\sum_{i>j}\frac{1}{\left|r_{i}-r_{j}\right|}+\sum_{I>J}\frac{Z_{I}Z_{J}}{\left|R_{I}-R_{J}\right|}-\sum_{i,I}\frac{Z_{I}}{\left|R_{I}-r_{i}\right|}\\\
+&\lambda\sum_{i}\sum_{k}\left(\mathbf{p}_{i}\cdot\bm{\epsilon}_{\pm}\right)e^{i\mathbf{k}\mathbf{r}_{i}}\frac{(\alpha_{k}+\beta^{\dagger}_{k})}{\sqrt{2\tilde{\omega}_{k}}}\\\
+&\lambda\sum_{i}\sum_{k}(\mathbf{p}_{i}\cdot\bm{\epsilon^{*}}_{\pm})e^{-i\mathbf{k}\mathbf{r}_{i}}\frac{(\alpha^{\dagger}_{k}+\beta_{k})}{\sqrt{2\tilde{\omega}_{k}}}\\\
+&\sum_{i\mathbf{k}\neq\mathbf{k}^{\prime}}\frac{N_{e}\lambda^{2}e^{i(\mathbf{k}-\mathbf{k}^{\prime})\mathbf{r}_{i}}}{2\sqrt{\tilde{\omega}_{k}\tilde{\omega}_{k^{\prime}}}}\left(\alpha_{k}+\beta^{\dagger}_{k}\right)\left(\beta_{k^{\prime}}+\alpha^{\dagger}_{k^{\prime}}\right).\end{split}$
(19)

Once more, the frequencies have been redefined as
$\tilde{\omega}_{k}=\sqrt{\omega^{2}_{k}+N_{e}\lambda^{2}}$.

##  IV QED coupled cluster for the minimal coupling Hamiltonian

We now present a QED coupled cluster (QED-CC) framework for the minimal
coupling Hamiltonian [58, 31, 35, 59, 60]. This approach is referred to as
minimal coupling QED-CCSD (MC-QED-CCSD). The starting point is the electron-
photon wave function

$\displaystyle\ket{\psi}=$
$\displaystyle\;\textrm{exp}\left((S_{1\alpha}+\gamma_{\alpha})\alpha^{\dagger}+(S_{1\beta}+\gamma_{\beta})\beta^{\dagger}\right)$
$\displaystyle\times$
$\displaystyle\;\textrm{exp}\left(T_{1}+T_{2}\right)\ket{\textrm{HF}}\otimes\ket{0,0},$
(20)

where $\ket{\textrm{HF}}$ is the no cavity Hartree-Fock Slater determinant
while $\ket{0,0}$ denotes the photonic vacuum for the $\alpha$ and $\beta$
photons introduced in Eq. (13). The operators in the exponentials of Eq. (20)
are electron ($T_{1}$ and $T_{2}$) and electron-photon (
$S^{\alpha}_{1}\alpha^{\dagger}$ and $S^{\beta}_{1}\beta^{\dagger}$)
excitation operators defined explicitly as

$\begin{split}T_{1}=&\sum_{ai}t^{a}_{i}E_{ai}\hskip
42.67912ptT_{2}=\frac{1}{2}\sum_{abij}t^{ab}_{ij}E_{ai}E_{bj}\\\
S_{1\alpha}\alpha^{\dagger}=&\sum_{ai}s^{a}_{i\alpha}E_{ai}\alpha^{\dagger}\hskip
19.91684ptS_{1\beta}\beta^{\dagger}=\sum_{ai}s^{a}_{i\beta}E_{ai}\beta^{\dagger}.\end{split}$
(21)

The amplitude parameters $t^{a}_{i}$,$t^{ab}_{ij}$, $s^{a}_{i\alpha}$,
$s^{a}_{i\beta}$ as well as $\gamma_{\alpha}$ and $\gamma_{\beta}$ are
determined in the ground state calculation. In Eq. 21, we use second
quantization for the electrons [61] with

$E_{pq}=\sum_{\sigma}a^{\dagger}_{p\sigma}a_{q\sigma},$ (22)

where $a^{\dagger}_{p\sigma}$ creates and $a_{p\sigma}$ annihilates an
electron in orbital $p$ with spin $\sigma$. The indices $i,j$ and $a,b$ label
occupied and virtual orbitals, respectively. In the limit where all
excitations are included in the wave function in Eq.20, the coupled cluster
expansion is exact and gives the same result as QED full configuration
interaction [24, 31, 32]. The optimal values for the coupled cluster
parameters are obtained by solving the projection equations [61]

$\begin{split}\Omega_{\mu,n,m}=&\bra{\mu,n,m}\textrm{exp}\left(-T\right)H\;\textrm{exp}\left(T\right)\ket{\textrm{HF},0,0}\\\
=&0,\end{split}$ (23)

with

$\begin{split}&\ket{\textrm{HF},0,0}=\ket{\textrm{HF}}\otimes\ket{0,0}\\\\[3.0pt]
&\ket{\mu,n,m}=\ket{\mu}\otimes\ket{m,n}\\\\[2.0pt]
&T=T_{1}+T_{2}+(S_{1\alpha}+\gamma_{\alpha})\alpha^{\dagger}+(S_{1\beta}+\gamma_{\beta})\beta^{\dagger},\end{split}$
(24)

where $\mu$ labels an electronic excitation while $n$ and $m$ are photonic
excitations in $\alpha$ and $\beta$. The ground state energy is equal to

$E=\bra{\textrm{HF},0,0}\textrm{exp}\left(-T\right)H\;\textrm{exp}\left(T\right)\ket{\textrm{HF},0,0}.$
(25)

Equation (25) can be rewritten in terms of the cluster amplitudes and the two
electron integrals $g_{pqrs}$ [61] as

$\begin{split}E=&\;E_{\textrm{HF}}+\sum_{aibj}(t^{ab}_{ij}+t^{a}_{i}t^{b}_{j})(2g_{iajb}-g_{ibja})\\\
+&\frac{\lambda}{\sqrt{2\tilde{\omega}}}\sum_{ai}\left[(\mathbf{p}\cdot\bm{\epsilon})e^{i\mathbf{k}\mathbf{r}}\right]_{ia}(s^{a}_{i\alpha}+\gamma_{\alpha}t^{a}_{i})\\\
+&\frac{\lambda}{\sqrt{2\tilde{\omega}}}\sum_{ai}\left[(\mathbf{p}\cdot\bm{\epsilon}^{*})e^{-i\mathbf{k}\mathbf{r}}\right]_{ia}(s^{a}_{i\beta}+\gamma_{\beta}t^{a}_{i}),\end{split}$
(26)

where

$E_{\textrm{HF}}=\bra{\textrm{HF}}H_{e}\ket{\textrm{HF}}.$ (27)

The electronic Hamiltonian $H_{e}$ is given by the standard expression

$H_{e}=\sum_{pq}h_{pq}E_{pq}+\frac{1}{2}\sum_{pqrs}g_{pqrs}(E_{pq}E_{rs}-\delta_{qr}E_{ps}),$
(28)

and $h_{pq}$ is the one electron integral [61].
The coupled cluster method can give complex energies for molecules in external
magnetic fields [62, 63, 64]. Since this case is implicitly contained in our
framework, we must ensure that the energy in Eq. (26) is real. This condition
is satisfied if $S_{1\alpha}$ and $S_{1\beta}$ (as well as $\gamma_{\alpha}$
and $\gamma_{\beta}$ ) fulfill the relations

$s^{a}_{i\alpha}=-s^{a*}_{i\beta}\hskip
34.14322pt\gamma_{\alpha}=-\gamma^{*}_{\beta}$ (29)

for real $t^{a}_{i}$ and $t^{ab}_{ij}$. These conditions are not accidental as
the same relationship holds between the integrals
$[(\mathbf{p}\cdot\bm{\epsilon}^{*})e^{-i\mathbf{k}\mathbf{r}}]_{pq}$ and
$[(\mathbf{p}\cdot\bm{\epsilon})e^{i\mathbf{k}\mathbf{r}}]_{pq}$ when the
orbitals $p$ and $q$ are real. Therefore, the conditions in Eq.(29) must arise
from a symmetry in the Hamiltonian.
In the following, we demonstrate that Eq. (29) holds in the exact limit

$\displaystyle\ket{\psi}=$
$\displaystyle\;\textrm{exp}(\sum_{mn}(S_{m\alpha,n\beta}+\gamma_{m\alpha,n\beta})\alpha^{\dagger\;m}\beta^{\dagger\;n})$
(30) $\displaystyle\times$
$\displaystyle\;\textrm{exp}(T)\ket{\textrm{HF},0,0},$ (31)

with $T$ and $S_{m\alpha,n\beta}$ containing all the electronic excitations.
The exact eigenfunction $\ket{\psi}$ satisfies

$H\ket{\psi}=E\ket{\psi}$ (32)

with real $E$. Applying a unitary transformation ($V$)

$\displaystyle V=$
$\displaystyle\textrm{exp}\left(i\pi(\alpha^{\dagger}\beta+\beta^{\dagger}\alpha)/2\right)\textrm{exp}\left(i\pi(\alpha^{\dagger}\alpha+\beta^{\dagger}\beta)/2\right)$
$\displaystyle\hskip 56.9055ptV^{\dagger}\alpha V=-\beta$ (33)
$\displaystyle\hskip 56.9055ptV^{\dagger}\beta V=-\alpha.$

on both sides of Eq. (32) we obtain

$V^{\dagger}HV\;V^{\dagger}\ket{\psi}=EV^{\dagger}\ket{\psi}.$ (34)

We can now use Eq. (33) to exchange $\alpha$ with $-\beta$ and viceversa in
both $H$ and the wave function. Complex conjugation of the transformed
Hamiltonian in Eq.(34) gives

$\begin{split}(V^{\dagger}HV)^{*}=&H_{e}+\tilde{\omega}\left(\alpha^{\dagger}\alpha+\beta^{\dagger}\beta+1\right)\\\
+&\frac{\lambda}{\sqrt{2\tilde{\omega}}}\sum_{pq}\left[\left(\mathbf{p}_{i}\cdot\bm{\epsilon}\right)e^{i\mathbf{k}\mathbf{r}_{i}}\right]_{pq}E_{pq}(\alpha+\beta^{\dagger})\\\
+&\frac{\lambda}{\sqrt{2\tilde{\omega}}}\sum_{pq}\left[(\mathbf{p}_{i}\cdot\bm{\epsilon}^{*})e^{-i\mathbf{k}\mathbf{r}_{i}}\right]_{pq}E_{pq}((\alpha^{\dagger}+\beta),\end{split}$
(35)

where we assume the orbitals $p,q$ are real. On the other hand, the
transformed wave function is equal to

$\begin{split}(V^{\dagger}\ket{\psi})^{*}=&\;\textrm{exp}(\sum_{mn}(S^{*}_{m\alpha,n\beta}+\gamma^{*}_{m\alpha,n\beta})\alpha^{\dagger\;n}\beta^{\dagger\;m}(-1)^{m+n})\\\
\times&\;\textrm{exp}(T^{*})\ket{\textrm{HF},0,0}.\end{split}$ (36)

Since the Hamiltonians in Eq. (35) and in Eq. (15) are identical, the
eigenfunctions in Eqs. (31) and (36) must also be equal to each other. This
implies that the amplitudes follow the relations

$\begin{split}S_{m,n}=&(-1)^{m+n}S^{*}_{n,m}\\\
\gamma_{m,n}=&(-1)^{m+n}\gamma^{*}_{n,m}\\\
\textrm{T}=&\textrm{T}^{*},\end{split}$ (37)

which proves the assertion.

##  V Results

Figure 3: Identification of the absolute configuration of a chiral molecule.
We first assign a priority number to each atom attached to the chiral center
based on the atomic number (from largest to smallest). The structure is then
rotated such that the lightest element (H in this case) is pointing backwards.
Finally, an arrow is drawn from the highest to the lowest priority element. If
the arrow rotates clockwise, we refer to the enantiomer as R, otherwise we
refer to the enantiomer as S. This procedure provides a unique identification
of the molecular structure [65] Figure 4:  a) Frequency dispersion of the
discriminating power in a LHCP cavity at coupling $\lambda=0.005$. The
discriminating power reaches its maximum value around $\omega=4$ eV. b)
Frequency dispersion of the discriminating power for different values of the
coupling strength. The size of the discriminating power increases with the
coupling but the qualitative shape of the dispersion is the same for all the
analyzed cases. The maximum in $\Delta(\textit{R}-\textit{S})$ is shifted at
higher frequencies as the coupling increases. Figure 5: Collective effects on
the cavity discriminating power for R and S for 1-fluoroethylamine. Figure 6:
Solvent induced enhancement of the chiral discriminating effect. The solvent
arranges in a chiral configuration around the chiral solute and the field
interacts with a larger chiral system.

In this section we apply the MC-QED-CCSD methodology on a set of chiral
molecules interacting with a handedness-preserving cavity. Specifically, we
perform a detailed quantitative analysis of the cavity-induced discriminating
effects and we show that rotational spectra present enantioselective
signatures in a chiral cavity. The calculations have been performed using a
development version of the eT program [66] using the cc-pVDZ basis set [67,
68]. The molecular structures have been optimized with the ORCA software
package [69] at the DFT-B3LYP/def2-SVP level111The geometries can be found in
the repository https://doi.org/10.5281/zenodo.7035887  [71] . The basis set
effects on the cavity calculations are discussed in Appendix B. To uniquely
appoint the enantiomers of a chiral molecule we use the absolute configuration
notation [65] explained in Fig. 3.

###  V.1 Cavity induced chiral discriminating power

Strong coupling with circularly polarized electromagnetic fields creates
energy differences between enantiomeric pairs. Since this is a field-induced
effect, the cavity properties strongly influence the magnitude of the
discriminating power. In Fig. 4, we show the frequency dispersion of the
energy difference between R and S enantiomers of proline,
$\Delta(\textit{R}-\textit{S})$, in a LHCP cavity. The sign of
$\Delta(\textit{R}-\textit{S})$ remains the same for every frequency of the
field and for every coupling value. The function changes sign when the
circular polarization is reversed. The magnitude of the discrimination is
strongly affected by both $\omega$ and $\lambda$ and shows the same
qualitative behaviour for all coupling values. Specifically,
$\Delta(\textit{R}-\textit{S})$ is equal to zero for $\omega\longrightarrow
0$, then reaches a minimum (or a maximum) for an intermediate frequency and
approaches zero again for $\omega\longrightarrow\infty$. The curve shape can
be explained using the previously developed theory. The discriminating power
of the cavity approaches zero at small $\omega$ because the exponential
$e^{ikr}\approx 1$ (see Section III). On the other hand,
$\Delta(\textit{R}-\textit{S})$ also approaches zero for large cavity
frequencies because the photonic component of the wave function is small for
highly energetic fields. The location of the minimum moves to higher
frequencies when the coupling increases as $\lambda$ enhances the population
of higher-energy photonic states.
The chiral discriminating powers reported in Fig.4 are very small for
realistic values of the light matter coupling ($\lambda\leq 0.05$ a.u.). In
particular, the field-induced energy differences in Fig.4 are significantly
smaller than the average thermal energy at room temperature. However, it is
well-known that field-induced modifications can be amplified by many orders of
magnitudes through cooperation mechanisms between multiple molecules and field
modes [72, 73, 74, 75, 19, 76]. This class of phenomena is commonly called
collective effects [77, 75]. In Fig.5 we investigate how the discriminating
power of the cavity is influenced by the number of chiral centers.
Specifically, we perform calculations on a set of 1-fluoroethylamine (both R
and S) separated by 200 Å either along the wave vector direction, Fig.5a, or
perpendicularly to $\mathbf{k}$, Fig.5b. The coupling is fixed to
$\lambda=0.05$ a.u. while the frequency of the cavity is equal to $13.6$ eV.
In both directions the discriminating power of the cavity is enhanced as the
number of chiral centers increases. The dispersions in Fig.5 are slower than
linear and the effect is dominated by the square root of the number of
enantiomers in the cavity. This trend is in line with the behaviour shown by
other field-induced effects, i.e. the Rabi splittings. We note, however, that
the enantiomeric discrimination discussed here is a non-resonant effect. The
fitting functions in Fig.5 have been chosen considering that the main energy
contribution from the light-matter interaction term scales as $\sqrt{N}$.
Moreover, for $N_{e}\lambda^{2}\gg\omega^{2}$ we have:

$\frac{\lambda}{\sqrt{2\sqrt{\omega^{2}+N_{e}\lambda^{2}}}}\approx\frac{1}{\sqrt{2\lambda\sqrt{N_{e}}}}\left(1-\frac{\omega^{2}}{2N_{e}}\right),$
(38)

which justifies the second contribution in the fit. The collective effects are
slightly stronger when molecules are displaced along $\mathbf{k}$.
Specifically, we note that the coefficient for the square root term in Fig.5a
is twice as large as in Fig.5b. Although the effect becomes infinitely large
as the number of molecules increases, the cavity-induced differences are not
extensive. If the chiral molecule is dissolved in a liquid, the solvent too
plays a significant role in enhancing the field-discriminating power. Indeed,
when the solute is chiral the solvent arranges in a chiral structure itself at
least in the first few solvation shells. Under these conditions the field
interacts with a significantly larger chiral system increasing the
discriminating effect. For example, in Fig.6 we show that when 10 water
molecules solvate 1-fluoroethylamine, the cavity effect almost doubles.
Finally, a molecule in chiral cavity interacts with more than one mode at a
time. Since the sign of the discriminating power is independent of the mode
frequency, the interaction with multiple cavity modes enhances the overall
cavity effects. All the factors discussed above will increase the field-
induced energy differences in chiral cavities, potentially to the kJ/mol
range.

###  V.2 Rotational spectra in chiral cavities

Figure 7: Orientational effects for R and S enantiomers of 1-fluoroethylamine
upon rotation around an arbitrary axis shown above of an angle $\theta$.
Figure 8: S enantiomer of the alanine molecule. Figure 9: Pictorial
representation of the Euler angles Figure 10: a) Rotational surface for S
alanine in a LCHP cavity. The surface is not flat, instead some configurations
are stabilized by the field. b) Difference in the rotational surfaces of R and
S alanine. We observe that the difference has a non-fixed sign. Figure 11:
Rotational spectra of R and S alanine inside an LCHP cavity.

Inside an optical cavity, the molecular energy is highly dependent on the
orientation because the quantization direction for the cavity field,
$\mathbf{k}$, naturally introduces a spacial anisotropy [24]. For example, in
Fig. 7 we show how the energy of 1-fluoroethylamine in a LHCP field changes
upon rotation of the molecule around an arbitrary axis. Since different
molecular orientations have different energies, the rotational spectrum is
modified by the cavity. We point out that modifications in the rotational
spectrum are to be expected even in the dipole approximation. However,
different enantiomers would still exhibit the same rotational levels. The
orientational effects reported in Fig. 7 for the R and S enantiomers are not
exactly the same, suggesting that the rotational surfaces might differ for the
two mirror images. These variations should induce enantioselective signatures
in rotational spectra. In this section we compute the cavity-induced
orientational effects on the two enantiomers of alanine (see Fig.8) and the
relative rotational spectra. All the possible molecular orientations are
obtained varying the three Euler angles $\phi$, $\theta$ and $\chi$ describing
a rotation around $z$, a rotation around $x$ and a rotation on around $z$
again, respectively (see Fig.9). The energy is invariant with respect to
$\chi$ for symmetry reasons. The Hamiltonian used to describe the nuclear
motion is

$\begin{split}H=&\sum_{I}\frac{p^{2}_{I}}{2M_{I}}+V(R),\end{split}$ (39)

where $V(R)$ is the potential energy surface obtained solving the electron-
photon Schrödinger equation. The interactions between the nuclei and the field
in Eq. (39) are neglected because the nuclear motion is much slower than the
electronic one. Once the motion of the center of mass has been removed, the
kinetic contribution in Eq. (39) can be split into a vibrational $H_{vib}$ and
a rotational part $H_{rot}$ (the mixed rovibrational contribution will be
neglected in this analysis). The Hamiltonian in Eq. (39) is therefore
rewritten as:

$\begin{split}&H=H_{rot}+H_{vib}\\\
H_{rot}=&\frac{J^{2}_{\xi}}{2I_{\xi}}+\frac{J^{2}_{\eta}}{2I_{\eta}}+\frac{J^{2}_{\zeta}}{2I_{\zeta}}+V_{rot}(\theta,\phi)\end{split}$
(40)

where $\xi$, $\eta$ and $\zeta$ are the principal axes of inertia of the
molecule, treated as a rigid body [78]. A detailed discussion of how to
compute the rotational energy levels from Eq. (40) is reported in Appendix C.
The rotational potential energy surface, $V_{rot}(\theta,\phi)$, for the S
enantiomer of alanine in a LHCP cavity is shown in Fig. 10a. Differently from
the no-cavity case, the surface is not flat and presents two maxima and two
minima, with some configurations being stabilized and others being
destabilized by the cavity. The difference in the orientational effects for
the two enantiomers ($\Delta(\textit{R}-\textit{S})$) is plotted in Fig. 10b.
In particular, the results shown in Fig. 10 are obtained summing up the cavity
orientational effects obtained from the first 20 modes of a chiral cavity with
fundamental excitation at $2.7$ eV. The coupling has been set to
$\lambda=0.05$ a.u. We note that in Fig. 10b, $\Delta(\textit{R}-\textit{S})$
does not have a constant sign implying that the surfaces are not just shifted
by a constant. Specifically, the potential of the S enantiomer has higher
maxima and lower minima when compared to the surface obtained for the R
enantiomer. The rotational spectra of the two mirror images are shown in Fig.
11. Due to the orientational effects in the cavity, the peaks are shifted to
higher energies than standard rotational spectroscopy. We highlight that R and
S spectra are not identical. While the intensities are mostly unchanged, the
peak positions are slightly different with the largest modifications observed
on the signals around 55, 70 and 100 cm-1. The shifts are not always in the
same direction and we observe that the R enantiomer spectrum has a lower
excitation around 70 cm-1 while the S enantiomer spectrum shows lower
transitions at 55 and 100 cm-1. Differences in the peak positions are on the
order of 0.5 cm-1, which is still large enough to be detected experimentally.
However, the shifts can be enhanced by solvent and collective effects as
discussed in the previous section. Together with the data shown before, these
results confirm that chiral cavities can be used to discriminate enantiomers
and that they induce enantiospecific shifts in their rotational spectra.

##  VI Conclusion

In this work, we present the first ab initio framework to model strong
coupling between chiral molecules and circularly polarized light. We show that
the simplest theoretical approach to properly describe these systems requires
the inclusion of two modes to ensure the correct field symmetry. Therefore, we
investigate the formal properties of the chiral cavity Hamiltonian and argue
that, using circularly polarized electromagnetic fields, it is indeed possible
to discriminate between the two mirror images of a chiral molecules. To
perform numerical simulations on realistic systems we develop a complex QED
coupled cluster approach and show that the energy is real if the appropriate
form of the cluster operator is employed. Our results demonstrate that while
non-chiral molecules have the same energy in cavities with different light
polarization, chiral cavities create energy differences between enantiomers.
Moreover, the sign of the chiral discrimination does not depend on the
frequency of the field or the coupling strength. Instead, it is only affected
by the circular polarization of the field. The dependence of the
discrimination on the number of chiral molecules has also been investigated.
In particular, we observe an increase in the discriminating power with the
number of chiral systems in the cavity. The solvent also enhances the
enantioselective effect by creating a chiral solvation shell around the chiral
solute. Finally, we demonstrate that enantiomers do not have the same
rotational spectrum in chiral cavities. Specifically, the circularly polarized
field induces enantiospecific shifts in the peak positions. We believe our
results pave the way for further investigations of strong coupling in
enantiomeric separation. Moreover, they suggest that interesting phenomena
should be observable upon placing molecules in chiral cavities.

### Acknowledgements

We acknowledge Tor S. Haugland for insightful discussions. R.R.R and H.K.
acknowledge funding from the Research Council of Norway through FRINATEK
Project No. 275506. We acknowledge computing resources through UNINETT
Sigma2—the National Infrastructure for High Performance Computing and Data
Storage in Norway, through Project No. NN2962k. This work has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 Research and Innovation Programme (grant agreement No. 101020016).

### Appendix A: Zero-point energy and non size-extensive effects

The $\theta$ value in the Bogoljubov transformation, Eq. (16), has been chosen
to remove the quadratic contribution in Eq.(9). The same $\theta$ value is
obtained when minimizing the zero-point energy. Indeed, after the Bogoljubov
transformation, the zero-point energy becomes equal to

$\mathscr{E}(\theta,N_{e})=\frac{N_{e}\lambda^{2}}{2\omega}\left(\cosh\theta-\sinh\theta\right)^{2}+2\omega\sinh^{2}\theta+\frac{\omega}{2},$
(41)

which has a minimum for

$\tanh 2\theta=\frac{\lambda^{2}N_{e}}{\lambda^{2}N_{e}+2\omega^{2}}.$ (42)

Using Eq.(42), the zero-point energy becomes equal to:

$\mathscr{E}_{0}(N_{e})=\frac{1}{2}\sqrt{\omega^{2}+N_{e}\lambda^{2}}=\frac{\omega}{2}\sqrt{1+\frac{N_{e}\lambda^{2}}{\omega^{2}}},$
(43)

which is non size-extensive due to the square root dependence on the number of
molecules. The non size-extensize contribution to the zero-point energy equals

$E_{nse}=\mathscr{E}_{0}(N_{e})-N_{e}\mathscr{E}_{0}(1)+\frac{(N_{e}-1)\omega}{2},$
(44)

where the last term has been added to avoid overcounting of the cavity zero-
point energy. When $\frac{N_{e}\lambda^{2}}{\omega^{2}}<<1$, Eq. (44) can be
expanded in Taylor series leading to

$E_{nse}\approx-\frac{N_{e}(N_{e}-1)\lambda^{4}}{8\omega^{3}}$ (45)

which shows the same $N_{e}(N_{e}-1)$ behaviour already reported in Ref.[24].

### Appendix B: Basis set effects in QED calculations

Figure 12: Dispersion of the QED-Effects for a system of water molecules
displaced of 10 Å from each other. The coupling is equal to 0.1 a.u. while the
frequency is equal to $13.6$ eV.

Usually, the minimal coupling Hamiltonian is not employed to describe the
interaction between light and matter, as it involves the complex momentum
operator. Instead, it is customary to perform calculations using the unitarily
transformed version of the minimal coupling Hamiltonian, called the PZW
Hamiltonian. When the dipole approximation is adopted, the PZW Hamiltonian is
equal to

$\begin{split}H_{\textrm{length\;gauge}}=&H_{e}+\lambda\sqrt{\frac{\omega}{2}}\sum_{i}(\mathbf{d}_{i}\cdot\bm{\epsilon})(b+b^{\dagger})\\\
+&\frac{\lambda^{2}}{2}\sum_{ij}(\mathbf{d}_{i}\cdot\bm{\epsilon})(\mathbf{d}_{j}\cdot\bm{\epsilon})+\omega
b^{\dagger}b,\end{split}$ (46)

where $H_{e}$ is the standard electronic Hamiltonian and $\mathbf{d}_{i}$ is
the molecular dipole operator. This Hamiltonian is referred to as length gauge
form. Extension of Eq. (46) beyond the dipole approximation is a relatively
easy task if a multipole expansion is applied. However, this approach breaks
the origin invariance of the problem. This is the reason why we have chosen a
modified minimal coupling Hamiltonian to describe strongly coupled systems in
chiral cavities. The two approaches (multipole expansion and minimal coupling)
are completely equivalent in theory, but the basis truncation in real
calculations makes them differ.

In this section we therefore compare the QED-effects obtained using the
minimal coupling Hamiltonian in dipole approximation Eq. (10), also called the
velocity-gauge Hamiltonian, and the length-gauge Hamiltonian in Eq.(46). In
Fig.12 we show the dispersion of the QED-effects for an increasing number of
water molecules spaced by 10 Å along the $\bm{\epsilon}$ direction. As
expected, the results obtained using the different gauges are not the same. In
particular, the velocity gauge consistently predicts larger field
contributions than length gauge. However, the agreement improves as the basis
set is enlarged. While the QED-effects computed using velocity gauge decrease
when the basis is enlarged, the effects computed using length gauge increase.
From Fig.12 we infer that both gauges converge to the complete basis limit,
one from above, and the other from below. Because of the non-variational
character of the coupled cluster approach, gauge invariance is not necessarily
reached in the complete basis limit unless the full set of excitations are
included [79]. However, as shown in previous works [80], reasonable agreement
between the two gauges can be reached with large basis sets and the two
frameworks describe the same qualitative effects regardless of the basis size.

### Appendix C: Calculation of the rotational spectrum

The molecular rotational spectra shown in this paper are obtained treating the
molecule as an asymmetrical top, following the theory in Ref. [78]. In
particular, the rotational Hamiltonian, $H_{\textrm{rot}}$ is equal to

$\begin{split}H_{rot}=&\frac{1}{4}\left(\mathbf{J}^{2}-J^{2}_{\zeta}\right)\left(\frac{1}{I_{\xi}}+\frac{1}{I_{\eta}}\right)+\frac{J^{2}_{\zeta}}{2I_{\zeta}}\\\
+&\frac{1}{8}\left(J^{2}_{+}+J^{2}_{-}\right)\left(\frac{1}{I_{\xi}}-\frac{1}{I_{\eta}}\right)+V(\theta,\phi),\end{split}$
(47)

expressed in terms of angular momentum operators $J_{+}$ and $J_{-}$. The
angular momentum component M along the fixed $z$ axis is a good quantum number
for the stationary states of Eq. (47). The eigenvalues and eigenvectors of
Eq.(47) are obtained solving the eigenvalue problem

$\sum_{Jk^{\prime}}(\bra{JMk}H_{rot}\ket{J^{\prime}Mk^{\prime}}-E\delta_{kk^{\prime}})c_{J^{\prime}k^{\prime}}=0,$
(48)

where $\\{\ket{JMk}\\}$ are the eigenfunctions of the symmetric top
problem[78]

$D^{j}_{Mk}(\chi,\theta,\phi)=\textrm{exp}(iM\chi)d^{j}_{Mn}\textrm{exp}(ik\phi).$
(49)

While the matrix elements for the $J$ operators are obtained using the angular
momentum properties, the contribution arising from the cavity-induced
potential is more difficult to account for. To include the contribution from
the cavity potential we perform a discrete Fourier transform in $\theta$ and
$\phi$

$V(\theta,\phi)\approx\sum_{ln}V_{ln}e^{il\theta}e^{in\phi},$ (50)

where we have only retained the $V_{ln}$ such that
$\absolutevalue{V_{ln}}>0.001$. The diagonalization is performed in
Mathematica for quantum number M = 0, 1, 2 and states with $J\leq 4$. The
intensities in the spectrum have been computed using the transition dipole

$f_{eg}=\frac{2}{3}(E_{e}-E_{g})\sum_{n=x,y,z}\absolutevalue{\bra{\psi_{e}}d_{n}\ket{\psi_{g}}}^{2},$
(51)

where $\psi_{g}$ and $\psi_{e}$ are the ground and excited state wave
functions, respectively, with energies $E_{g}$ and $E_{e}$.

## References

  * Flick  _et al._ [2018] J. Flick, N. Rivera, and P. Narang, Strong light-matter coupling in quantum chemistry and quantum photonics, Nanophotonics 7, 1479 (2018).
  * Li  _et al._ [2021a] T. E. Li, B. Cui, J. E. Subotnik, and A. Nitzan, Molecular polaritonics: Chemical dynamics under strong light–matter coupling, Annu. Rev. Phys. Chem. 73 (2021a).
  * Ebbesen [2016] T. W. Ebbesen, Hybrid light–matter states in a molecular and material science perspective, Accounts Chem. Res. 49, 2403 (2016).
  * Garcia-Vidal  _et al._ [2021] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Manipulating matter by strong coupling to vacuum fields, Science 373, eabd0336 (2021).
  * Hertzog  _et al._ [2019] M. Hertzog, M. Wang, J. Mony, and K. Borjesson, Strong light–matter interactions: a new direction within chemistry, Chem. Soc. Rev. 48, 937 (2019).
  * Bloch  _et al._ [2022] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Strongly correlated electron–photon systems, Nature , 1 (2022).
  * Sentef  _et al._ [2018] M. A. Sentef, M. Ruggenthaler, and A. Rubio, Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity, Sci. Adv. 4, eaau6969 (2018).
  * Ashida  _et al._ [2020] Y. Ashida, A. Imamoglu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, Quantum electrodynamic control of matter: Cavity-enhanced ferroelectric phase transition, Phys. Rev. X 10, 041027 (2020).
  * Mandal and Huo [2019] A. Mandal and P. Huo, Investigating new reactivities enabled by polariton photochemistry, J. Phys. Chem. Lett. 10, 5519 (2019).
  * Fregoni  _et al._ [2018] J. Fregoni, G. Granucci, E. Coccia, M. Persico, and S. Corni, Manipulating azobenzene photoisomerization through strong light–molecule coupling, Nat. Comm. 9, 1 (2018).
  * Fregoni  _et al._ [2020] J. Fregoni, S. Corni, M. Persico, and G. Granucci, Photochemistry in the strong coupling regime: A trajectory surface hopping scheme, J. Comput. Chem. 41, 2033 (2020).
  * Loudon [2000] R. Loudon,  _The quantum theory of light_ (OUP Oxford, 2000).
  * Viviescas and Hackenbroich [2003] C. Viviescas and G. Hackenbroich, Field quantization for open optical cavities, Phys. Rev. A 67, 013805 (2003).
  * Plum and Zheludev [2015] E. Plum and N. I. Zheludev, Chiral mirrors, Appl. Phys. Lett. 106, 221901 (2015).
  * Favero and Karrai [2009] I. Favero and K. Karrai, Optomechanics of deformable optical cavities, Nat. Photonics 3, 201 (2009).
  * Chikkaraddy  _et al._ [2016] R. Chikkaraddy, B. De Nijs, F. Benz, S. J. Barrow, O. A. Scherman, E. Rosta, A. Demetriadou, P. Fox, O. Hess, and J. J. Baumberg, Single-molecule strong coupling at room temperature in plasmonic nanocavities, Nature 535, 127 (2016).
  * Fitzgerald  _et al._ [2016] J. M. Fitzgerald, P. Narang, R. V. Craster, S. A. Maier, and V. Giannini, Quantum plasmonics, P. IEEE 104, 2307 (2016).
  * Flick  _et al._ [2017] J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (qed) chemistry, P. Natl. Acad. Sci. 114, 3026 (2017).
  * Sidler  _et al._ [2020] D. Sidler, C. Schafer, M. Ruggenthaler, and A. Rubio, Polaritonic chemistry: Collective strong coupling implies strong local modification of chemical properties, J. Phys. Chem. Lett. 12, 508 (2020).
  * Galego  _et al._ [2015] J. Galego, F. J. Garcia-Vidal, and J. Feist, Cavity-induced modifications of molecular structure in the strong-coupling regime, Phys. Rev. X 5, 041022 (2015).
  * Forbes and Andrews [2021] K. A. Forbes and D. L. Andrews, Orbital angular momentum of twisted light: chirality and optical activity, J. Phys. Photon. 3, 022007 (2021).
  * Galego  _et al._ [2019] J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, Cavity casimir-polder forces and their effects in ground-state chemical reactivity, Phys. Rev. X 9, 021057 (2019).
  * Hutchison  _et al._ [2012] J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, Modifying chemical landscapes by coupling to vacuum fields, Angew. Chem. Int. Edit. 51, 1592 (2012).
  * Haugland  _et al._ [2021] T. S. Haugland, C. Schafer, E. Ronca, A. Rubio, and H. Koch, Intermolecular interactions in optical cavities: An ab initio qed study, J. Chem. Phys. 154, 094113 (2021).
  * Riso  _et al._ [2022a] R. R. Riso, T. S. Haugland, E. Ronca, and H. Koch, On the characteristic features of ionization in qed environments, J. Chem. Phys. 156, 234103 (2022a).
  * Vergauwe  _et al._ [2016] R. M. Vergauwe, J. George, T. Chervy, J. A. Hutchison, A. Shalabney, V. Y. Torbeev, and T. W. Ebbesen, Quantum strong coupling with protein vibrational modes, J. Phys. Chem. Lett. 7, 4159 (2016).
  * Lather  _et al._ [2019] J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen, and J. George, Cavity catalysis by cooperative vibrational strong coupling of reactant and solvent molecules, Angew. Chem. 131, 10745 (2019).
  * Thomas  _et al._ [2019] A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, _et al._ , Tilting a ground-state reactivity landscape by vibrational strong coupling, Science 363, 615 (2019).
  * Thomas  _et al._ [2016] A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, _et al._ , Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field, Angew. Chem. 128, 11634 (2016).
  * Schafer  _et al._ [2021a] C. Schafer, 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, P. Natl. Acad. Sci. 118 (2021a).
  * Haugland  _et al._ [2020] T. S. Haugland, E. Ronca, E. F. Kjonstad, A. Rubio, and H. Koch, Coupled cluster theory for molecular polaritons: Changing ground and excited states, Phys. Rev. X 10, 041043 (2020).
  * Riso  _et al._ [2022b] R. R. Riso, T. S. Haugland, E. Ronca, and H. Koch, Molecular orbital theory in cavity qed environments, Nat. Comm. 13, 1 (2022b).
  * Mandal  _et al._ [2020] A. Mandal, S. Montillo Vega, and P. Huo, Polarized fock states and the dynamical casimir effect in molecular cavity quantum electrodynamics, J. Phys. Chem. Lett. 11, 9215 (2020).
  * Ashida  _et al._ [2021] Y. Ashida, A. Imamoglu, and E. Demler, Cavity quantum electrodynamics at arbitrary light-matter coupling strengths, Phys. Rev. Lett. 126, 153603 (2021).
  * DePrince III [2021] A. E. DePrince III, Cavity-modulated ionization potentials and electron affinities from quantum electrodynamics coupled-cluster theory, J. Chem. Phys. 154, 094112 (2021).
  * Flick and Narang [2020] J. Flick and P. Narang, Ab initio polaritonic potential-energy surfaces for excited-state nanophotonics and polaritonic chemistry, J. Chem. Phys. 153, 094116 (2020).
  * Schafer  _et al._ [2021b] C. Schafer, J. Flick, E. Ronca, P. Narang, and A. Rubio, Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity, arXiv preprint arXiv:2104.12429 (2021b).
  * Petryayeva and Krull [2011] E. Petryayeva and U. J. Krull, Localized surface plasmon resonance: Nanostructures, bioassays and biosensing—a review, Anal. Chim. Acta 706, 8 (2011).
  * Maccaferri  _et al._ [2021] N. Maccaferri, G. Barbillon, A. N. Koya, G. Lu, G. P. Acuna, and D. Garoli, Recent advances in plasmonic nanocavities for single-molecule spectroscopy, Nanoscale Adv. 3, 633 (2021).
  * Yang  _et al._ [2020] B. Yang, G. Chen, A. Ghafoor, Y. Zhang, Y. Zhang, Y. Zhang, Y. Luo, J. Yang, V. Sandoghdar, J. Aizpurua, _et al._ , Sub-nanometre resolution in single-molecule photoluminescence imaging, Nat. Photonics 14, 693 (2020).
  * Romanelli  _et al._ [2021] M. Romanelli, G. Dall’Osto, and S. Corni, Role of metal-nanostructure features on tip-enhanced photoluminescence of single molecules, J. Chem. Phys. 155, 214304 (2021).
  * Zhang  _et al._ [2013] R. Zhang, Y. Zhang, Z. Dong, S. Jiang, C. Zhang, L. Chen, L. Zhang, Y. Liao, J. Aizpurua, Y. e. Luo, _et al._ , Chemical mapping of a single molecule by plasmon-enhanced raman scattering, Nature 498, 82 (2013).
  * Schnoering  _et al._ [2021] G. Schnoering, S. Albert, A. Canaguier-Durand, and C. Genet, Chiral thermodynamics in tailored chiral optical environments, Phys. Rev. X 11, 041022 (2021).
  * Sato  _et al._ [2004] I. Sato, R. Sugie, Y. Matsueda, Y. Furumura, and K. Soai, Asymmetric synthesis utilizing circularly polarized light mediated by the photoequilibrium of chiral olefins in conjunction with asymmetric autocatalysis, Angew. Chem. 116, 4590 (2004).
  * Barron [2004] L. Barron,  _Molecular Light Scattering and Optical Activity_ (Cambridge University Press, 2004).
  * Craig and Thirunamachandran [1998] D. P. Craig and T. Thirunamachandran,  _Molecular quantum electrodynamics: an introduction to radiation-molecule interactions_ (Courier Corporation, 1998).
  * Hubener  _et al._ [2021] H. Hubener, U. De Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, Engineering quantum materials with chiral optical cavities, Nat. Mater. 20, 438 (2021).
  * Wang  _et al._ [2019] X. Wang, E. Ronca, and M. A. Sentef, Cavity quantum electrodynamical chern insulator: Towards light-induced quantized anomalous hall effect in graphene, Phys. Rev. B 99, 235156 (2019).
  * Schafer  _et al._ [2018] C. Schafer, M. Ruggenthaler, and A. Rubio, Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling, Phys. Rev. A 98, 043801 (2018).
  * Rokaj  _et al._ [2018] V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy, J. Phys. B-At. Mol. Opt. 51, 034005 (2018).
  * Di Stefano  _et al._ [2019] O. Di Stefano, A. Settineri, V. Macrì, L. Garziano, R. Stassi, S. Savasta, and F. Nori, Resolution of gauge ambiguities in ultrastrong-coupling cavity quantum electrodynamics, Nat. Phys. 15, 803 (2019).
  * Schafer  _et al._ [2020] C. Schafer, M. Ruggenthaler, V. Rokaj, and A. Rubio, Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics, ACS photonics 7, 975 (2020).
  * Li  _et al._ [2021b] X. Li, A. Mandal, and P. Huo, Theory of mode-selective chemistry through polaritonic vibrational strong coupling, J. Phys. Chem. Lett. 12, 6974 (2021b).
  * Fregoni  _et al._ [2021] J. Fregoni, T. S. Haugland, S. Pipolo, T. Giovannini, H. Koch, and S. Corni, Strong coupling between localized surface plasmons and molecules by coupled cluster theory, Nano Lett. 21, 6664 (2021).
  * Pannir-Sivajothi  _et al._ [2022] S. Pannir-Sivajothi, J. A. Campos-Gonzalez-Angulo, L. A. Martínez-Martínez, S. Sinha, and J. Yuen-Zhou, Driving chemical reactions with polariton condensates, Nat. Comm. 13, 1 (2022).
  * Appelquist  _et al._ [1986] T. Appelquist, M. J. Bowick, D. Karabali, and L. Wijewardhana, Spontaneous breaking of parity in (2+ 1)-dimensional qed, Phys. Rev. D 33, 3774 (1986).
  * Rokaj  _et al._ [2022] V. Rokaj, M. Ruggenthaler, F. G. Eich, and A. Rubio, Free electron gas in cavity quantum electrodynamics, Phys. Rev. Res. 4, 013012 (2022).
  * Liebenthal  _et al._ [2022] M. D. Liebenthal, N. Vu, and A. E. DePrince III, Equation-of-motion cavity quantum electrodynamics coupled-cluster theory for electron attachment, J. Chem. Phys. 156, 054105 (2022).
  * Pavosevic and Flick [2021] F. Pavosevic and J. Flick, Polaritonic unitary coupled cluster for quantum computations, J. Phys. Chem. Lett. 12, 9100 (2021).
  * Pavosevic  _et al._ [2022] F. Pavosevic, S. Hammes-Schiffer, A. Rubio, and J. Flick, Cavity-modulated proton transfer reactions, J. Am. Chem. Soc. 144, 4995 (2022).
  * Helgaker  _et al._ [2014] T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, 2014).
  * Stopkowicz  _et al._ [2015] S. Stopkowicz, J. Gauss, K. K. Lange, E. I. Tellgren, and T. Helgaker, Coupled-cluster theory for atoms and molecules in strong magnetic fields, J. Chem. Phys. 143, 074110 (2015).
  * Hampe and Stopkowicz [2017] F. Hampe and S. Stopkowicz, Equation-of-motion coupled-cluster methods for atoms and molecules in strong magnetic fields, J. Chem. Phys. 146, 154105 (2017).
  * Piecuch  _et al._ [1990] P. Piecuch, S. Zarrabian, J. Paldus, and J. Čížek, Coupled-cluster approaches with an approximate account of triexcitations and the optimized-inner-projection technique. ii. coupled-cluster results for cyclic-polyene model systems, Phys. Rev. B 42, 3351 (1990).
  * Clayden  _et al._ [2012] J. Clayden, N. Greeves, and S. Warren,  _Organic Chemistry_ (OUP Oxford, 2012).
  * Folkestad  _et al._ [2020] S. D. Folkestad, E. F. Kjønstad, R. H. Myhre, J. H. Andersen, A. Balbi, S. Coriani, T. Giovannini, L. Goletto, T. S. Haugland, A. Hutcheson, _et al._ , e t 1.0: An open source electronic structure program with emphasis on coupled cluster and multilevel methods, J. Chem. Phys. 152, 184103 (2020).
  * Pritchard  _et al._ [2019] B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibsom, and T. L. Windus, A new basis set exchange: An open, up-to-date resource for the molecular sciences community, J. Chem. Inf. Model. 59, 4814 (2019).
  * Dunning [1989] T. H. Dunning, Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen, J. Chem. Phys. 90, 1007 (1989).
  * Neese  _et al._ [2020] F. Neese, F. Wennmohs, U. Becker, and C. Riplinger, The orca quantum chemistry program package, J. Chem. Phys. 152, 224108 (2020).
  * Note [1] The geometries can be found in the repository https://doi.org/10.5281/zenodo.7035887.
  * Weigend  _et al._ [1998] F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, Ri-mp2: optimized auxiliary basis sets and demonstration of efficiency, Chemical physics letters 294, 143 (1998).
  * Sáez-Blázquez  _et al._ [2018] R. Sáez-Blázquez, J. Feist, F. García-Vidal, and A. Fernández-Domínguez, Photon statistics in collective strong coupling: Nanocavities and microcavities, Phys. Rev. A 98, 013839 (2018).
  * Ribeiro  _et al._ [2021] R. F. Ribeiro, J. A. Campos-Gonzalez-Angulo, N. C. Giebink, W. Xiong, and J. Yuen-Zhou, Enhanced optical nonlinearities under collective strong light-matter coupling, Phys. Rev. A 103, 063111 (2021).
  * Li  _et al._ [2021c] T. E. Li, A. Nitzan, and J. E. Subotnik, Collective vibrational strong coupling effects on molecular vibrational relaxation and energy transfer: Numerical insights via cavity molecular dynamics simulations, Angew. Chem. 133, 15661 (2021c).
  * Mandal  _et al._ [2022] A. Mandal, X. Li, and P. Huo, Theory of vibrational polariton chemistry in the collective coupling regime, J. Chem. Phys. 156, 014101 (2022).
  * 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, J. Chem. Phys. 156, 230901 (2022).
  * Herrera and Owrutsky [2020] F. Herrera and J. Owrutsky, Molecular polaritons for controlling chemistry with quantum optics, J. Chem. Phys. 152, 100902 (2020).
  * Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , Vol. 3 (Elsevier, 2013).
  * Pedersen and Koch [1997] T. B. Pedersen and H. Koch, Coupled cluster response functions revisited, J. Chem. Phys. 106, 8059 (1997).
  * Pedersen and Koch [1998] T. B. Pedersen and H. Koch, Gauge invariance of the coupled cluster oscillator strength, Chem. Phys. Lett. 293, 251 (1998).