# Ab-Initio Vibro-Polaritonic Spectra in Strongly Coupled Cavity-Molecule
Systems

Thomas Schnappinger<EMAIL_ADDRESS>Department of Physics,
Stockholm University, AlbaNova University Center, SE-106 91 Stockholm, Sweden
Markus Kowalewski<EMAIL_ADDRESS>Department of Physics,
Stockholm University, AlbaNova University Center, SE-106 91 Stockholm, Sweden

###### Abstract

Recent experiments have revealed the profound effect of strong light-matter
interactions in optical cavities on the electronic ground state of molecular
systems. This phenomenon, known as vibrational strong coupling (VSC), can
modify reaction rates and induce the formation of molecular vibrational
polaritons, hybrid states involving both photon modes and vibrational modes of
molecules. We present an ab-initio methodology, based on the cavity Born-
Oppenheimer Hartree-Fock ansatz, which is specifically powerful for ensembles
of molecules, to calculate vibro-polaritonic IR spectra. This method allows a
comprehensive analysis of these hybrid states. Our semi-classical approach,
validated against full quantum simulations, reproduces key features of the
vibro-polaritonic spectra. The underlying analytic gradients also pave the way
for optimizing cavity-coupled molecular systems and performing semi-classical
dynamics simulations.

##  1 Introduction

Strong light-matter interactions between a molecular system and the
electromagnetic field of an optical cavity offer new possibilities to modify
chemical reactivity and selectivity, as demonstrated in recent experiments 1,
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12. A particularly intriguing but not yet
fully understood situation occurs when a cavity mode is strongly coupled to
molecular vibrations, called vibrational-strong coupling (VSC) 13, 14, 15, 16,
17, 18, 19. For VSC the rate constants of the ground state reactions can be
modified even without external driving, i.e. without explicitly adding photons
into the cavity. Probably one of the most striking features observed in such
experiments is the change of the vibrational spectra due to the formation of
molecular vibrational polaritons, hybrid states that involve both cavity modes
and vibrational modes of molecules. These vibrational polaritons are
characterized by the appearance of a lower polariton (LP) transition and a
upper polariton (UP) transition separated by the Rabi splitting frequency
$\Omega_{R}$ originating from a single vibrational peak in the spectrum of the
coupled molecular-cavity system. Both the presence of vibrational polaritons
and the change in chemical reactivity are characterized by a sharp resonance
20, 21, 22, which can be achieved if one of the cavity frequencies
$\omega_{c}$ is in resonance with a vibrational mode in the reaction mixture.

Despite the large number of theoretical studies in the literature, the
understanding of the underlying microscopic and macroscopic physical
mechanisms, especially with respect to the effects of VSC, is still incomplete
and under discussion 23, 15, 16, 17, 18, 19, 24, 25. To better understand the
formation of molecular vibrational polaritons and their role in polaritonic
chemistry, a reliable theoretical description is needed that can accurately
treat also the in practice most relevant case of collective strong coupling.
One way to achieve such a description is to simulate molecules in the VSC
regime in the time domain, capturing the dynamics of the system to analyze
optical spectra and chemical reactivity 26, 27, 28, 29. These explicit time
domain calculations have their advantages in simulating complex and anharmonic
dynamics, but suffer from high computational cost, especially when the
dynamics is calculated in a ab-initio framework. Following the idea of Bonini
and Flick 30, the use of the generalized force-constant matrix, often called
the mass-weighted Hessian matrix $\bm{H}^{M}$ in quantum chemistry, offers an
alternative way to obtain vibro-polaritonic IR spectra without simulating the
temporal evolution of the system. Within the well-known harmonic
approximation, the eigenvectors and eigenvalues of $\bm{H}^{M}$ give the
vibro-polaritonic normal modes of the correlated matter-photon system and the
frequencies of the vibro-polaritons.

In our recent work 19 we have introduced a Hartree-Fock ansatz in the cavity
Born-Oppenheimer approximation (CBOA) 31, 32, 26, 33, capable of describing
the electronic ground state of single molecules as as well as of an ensemble
of molecules coupled to an optical cavity. Within the framework of this cavity
Born-Oppenheimer Hartree-Fock (CBO-HF) ansatz, we now derive analytic
expressions for the first derivatives of the energy with respect to the
nuclear and photonic degrees of freedom. These analytic gradients are then
used to construct the mass-weighted Hessian matrix $\bm{H}^{M}$ via finite
differences and subsequently ab-initio vibro-polaritonic IR spectra within the
harmonic approximation are calculated. After introducing the CBO-HF formalism,
the first part of this work describes how the analytic expressions for the
first derivatives of the energy are obtained and used to calculate vibro-
polaritonic IR spectra in the CBO-HF framework. Next, the vibro-polaritonic IR
spectra in the harmonic approximation are carefully compared against a full
quantum mechanical approach for the rather anharmonic case of one and two
diatomic hydrogen fluoride HF molecules. Even in this rather simple molecular
example, we observe that the interaction between an optical cavity and the
molecule(s) leads to a detuning and a change in the optical length of the
cavity for both the harmonic approximation and the full quantum treatment, in
agreement with our recent findings 18. Since both the self-consistent field
(SCF) treatment and the correct description of the full dipole self-energy
(DSE) contribution are crucial to capture relevant aspects for the description
of strongly coupled molecules and their chemical properties 34, 17, 18, 19,
35, we analyze their direct influence on the vibro-polaritonic IR spectra of
small ensembles of HF molecules. In the last part, we leave the case of
diatomic molecules and discuss the vibro-polaritonic IR spectra for ammonia
(NH3). Here, the interaction with the confined light mode of an optical cavity
influences the whole vibrational spectrum not only by creating vibrational
polaritons but also by allowing symmetry breaking.

##  2 Theory

To describe the interaction of a confined light mode with atoms, molecules,
and condensed matter systems, the non-relativistic Pauli–Fierz Hamiltonian 36,
37, 38, 22, 35 in the length gauge, assuming the dipole approximation is an
accurate starting point. If the quantized cavity modes are coupled via their
characteristic frequencies to vibrational degrees of freedom of molecules, the
situation is described as VSC, for which the CBOA 31, 32, 26, 33 is a well
suited theoretical approach. Within the CBOA, the cavity modes are grouped
with the nuclei in a generalized Born-Huang expansion 39, 35, and then one can
subsequently solve the quantum problem of the electrons and then of the nuclei
and photons. Here we first focus on the electronic degrees of freedom that
have been shown to be strongly affected even for cavity modes that are tuned
to the vibrational degrees of freedom 18, 19. In the following bold symbols
denote vectors, and atomic units ($\hbar=4\pi;\varepsilon_{0}=m_{e}=1$) are
used throughout the paper, unless otherwise noted. For a single mode cavity
the electronic CBOA Hamiltonian takes the form

$\begin{split}\hat{H}_{CBO}&=\hat{H}_{el}+\frac{1}{2}\omega_{c}^{2}q_{c}^{2}-\omega_{c}q_{c}\left(\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\right)+\frac{1}{2}\left(\bm{\lambda}_{c}\cdot\bm{\hat{\mu}}\right)^{2}\,,\end{split}$
(1)

where

$\begin{split}\quad\bm{\hat{\mu}}=\bm{\hat{\mu}}_{el}+\bm{\mu}_{Nuc}=-\sum_{i=1}^{N_{el}}\bm{\hat{r}}_{i}+\sum_{A=1}^{N_{Nuc}}Z_{A}\bm{R}_{A}\,,\end{split}$
(2)

represents the molecular dipole operator, which is defined by the operators of
the $N_{el}$ electron coordinates $\bm{\hat{r}}$, the classic coordinates
$\bm{R}$ of the $N_{Nuc}$ nuclei and the nuclear charges $Z_{A}$.
$\hat{H}_{el}$ is the Hamiltonian for the field-free many-electron system, and
the second term defines the harmonic potential introduced by the photon
displacement field, with the photon displacement coordinate $q_{c}$ and
$\omega_{c}$ being the frequency of the cavity mode. The third term of Eq. (1)
describes the dipole coupling between the molecular system and the photon
displacement field, which is characterized by the coupling strength
$\bm{\lambda}_{c}$. The last term is the DSE operator 40, 34, 17, which is an
energy contribution that describes the self-polarization of the molecule-
cavity system. The coupling parameter $\bm{\lambda}_{c}$ for a cavity with an
effective mode volume $V_{c}$ is defined as follows:

$\bm{\lambda}_{c}=\bm{e}\lambda_{c}=\bm{e}\sqrt{\frac{4\pi}{V_{c}}}\,.$ (3)

The unit vector $\bm{e}$ denotes the polarization axis of the cavity mode.

###  2.1 CBO-Hartree-Fock

In our recent work 19 we have introduced the CBO-HF approach representing a
formulation of the well-known Hartree-Fock ansatz in the context of the CBOA.
The resulting energy expectation value $\bigl{\langle}E_{CBO}\bigr{\rangle}$
consists of four energy contributions:

$\bigl{\langle}E_{CBO}\bigr{\rangle}=E_{el}+E_{lin}+E_{dse}+E_{dis}\quad\text{
with}\quad E_{dis}=\frac{1}{2}\omega_{c}^{2}q_{c}^{2}\,.$ (4)

The first term $E_{el}$ contains all Hartree-Fock energy components of the
many-electron system 41, 42 and $E_{lin}$ describes the linear dipole coupling
between the photon displacement field, the electrons, and the nuclei. The
energy contribution $E_{dse}$ is due to the DSE operator in the Hamiltonian
and $E_{dis}$ is the energy resulting from the photon displacement field 34.
Following the standard procedure the transformation of
$\bigl{\langle}E_{CBO}\bigr{\rangle}$ in the basis of atomic orbitals results
in

$\bigl{\langle}E_{CBO}\bigr{\rangle}=\sum_{\alpha,\beta}^{M}D_{\alpha,\beta}F_{\alpha,\beta}+\tilde{E}_{n,c}=\sum_{\alpha,\beta}^{M}D_{\alpha,\beta}\bigl{\langle}\alpha\big{|}\tilde{h}\big{|}\beta\bigr{\rangle}+\frac{1}{2}\sum_{\alpha,\beta,\gamma,\delta}^{M}D_{\alpha,\beta}D_{\gamma,\delta}\bigl{\langle}\alpha\beta\big{|}\tilde{g}\big{|}\gamma\delta\bigr{\rangle}+\tilde{E}_{n,c}$
(5)

where $\alpha,\beta,\gamma,\delta$ denote atomic orbitals, $F_{\alpha,\beta}$
the CBO-HF Fock matrix and $D_{\alpha,\beta}$ the corresponding density matrix
elements. Due to the classical nature of the nuclei and the photon
displacement field in the CBOA, when we determine the electronic ground state,
all their energy contributions are scalar and summed up in $\tilde{E}_{n,c}$:

$\tilde{E}_{nuc,c}=V_{nn}+E_{dis}+E^{(nuc)}_{lin}+E^{(nuc)}_{dse}.$ (6)

The modified one-electron integrals
$\bigl{\langle}\alpha\big{|}\tilde{h}\big{|}\beta\bigr{\rangle}$ and two-
electron integrals
$\bigl{\langle}\alpha\beta\big{|}\tilde{g}\big{|}\gamma\delta\bigr{\rangle}$
used to build the CBO-HF Fock matrix elements $F_{\alpha,\beta}$ consist of
the standard one-electron integral
$\bigl{\langle}\alpha\big{|}h\big{|}\beta\bigr{\rangle}$ and two-electron
integrals $\bigl{\langle}\alpha\beta\big{|}g\big{|}\gamma\delta\bigr{\rangle}$
as well as terms describing the linear cavity-electron interaction and the
electronic DSE contributions:

$\begin{split}\bigl{\langle}\alpha\big{|}\tilde{h}\big{|}\beta\bigr{\rangle}&=\bigl{\langle}\alpha\big{|}\hat{h}\big{|}\beta\bigr{\rangle}-\omega_{c}q_{c}\bigl{\langle}\alpha\big{|}\bm{\lambda}_{c}\cdot\bm{\hat{r}}\big{|}\beta\bigr{\rangle}-\bm{\lambda}_{c}\cdot\bm{\mu}_{Nuc}\bigl{\langle}\alpha\big{|}\bm{\lambda}_{c}\cdot\bm{\hat{r}}\big{|}\beta\bigr{\rangle}+\frac{1}{2}\bigl{\langle}\alpha\big{|}\left(\bm{\lambda}_{c}\cdot\bm{\hat{r}}\right)^{2}\big{|}\beta\bigr{\rangle}\\\
\bigl{\langle}\alpha\beta\big{|}\tilde{g}\big{|}\gamma\delta\bigr{\rangle}&=\bigl{\langle}\alpha\beta\big{|}g\big{|}\gamma\delta\bigr{\rangle}+\bigl{\langle}\alpha\big{|}\bm{\lambda}_{c}\cdot\bm{\hat{r}}\big{|}\beta\bigr{\rangle}\bigl{\langle}\gamma\big{|}\bm{\lambda}_{c}\cdot\bm{\hat{r}}\big{|}\delta\bigr{\rangle}-\bigl{\langle}\alpha\big{|}\bm{\lambda}_{c}\cdot\bm{\hat{r}}\big{|}\delta\bigr{\rangle}\bigl{\langle}\gamma\big{|}\bm{\lambda}_{c}\cdot\bm{\hat{r}}\big{|}\beta\bigr{\rangle}\end{split}$
(7)

###  2.2 CBO-Hartree-Fock Gradients

The first derivative of the energy $\bigl{\langle}E_{CBO}\bigr{\rangle}$ with
respect to a nuclear or photon displacement coordinate $\zeta_{i}$ is

$\frac{\partial}{\partial\zeta_{i}}\bigl{\langle}E_{CBO}\bigr{\rangle}=\sum_{\alpha,\beta}^{M}D_{\alpha,\beta}\frac{\partial}{\partial\zeta_{i}}F_{\alpha,\beta}+\sum_{\alpha,\beta}^{M}F_{\alpha,\beta}\frac{\partial}{\partial\zeta_{i}}D_{\alpha,\beta}+\frac{\partial}{\partial\zeta_{i}}\tilde{E}_{n,c}.$
(8)

The first term in Eq. (8) is the Hellmann-Feynman term, the second term is the
wavefunction derivative, or Pulay term, and the last part is the derivative of
all scalar energy contributions. As the CBO-HF wave function is variational
optimized, the explicit calculation of the density derivatives can be avoided
43 and the Pulay term can be written in terms of orbital energies
$\epsilon_{i}$, orbital coefficients $c_{i}$ and overlap integral derivatives:

$\begin{split}\sum_{\alpha,\beta}^{M}F_{\alpha,\beta}\frac{\partial
D_{\alpha,\beta}}{\partial\zeta_{i}}&=-\sum_{\alpha,\beta}^{M}W_{\alpha,\beta}\frac{\partial
S_{\alpha,\beta}}{\partial\zeta_{i}}\\\
&=\sum_{\alpha,\beta}^{M}\sum_{i=1}^{N_{oc}}\epsilon_{i}c_{i,\alpha}^{*}c_{i,\beta}\left(\bigl{\langle}\frac{\partial\alpha}{\zeta_{i}}\big{|}\beta\bigr{\rangle}+\bigl{\langle}\alpha\big{|}\frac{\partial\beta}{\zeta_{i}}\bigr{\rangle}\right)\end{split}$
(9)

In the case of $\zeta_{i}$ being a photon displacement coordinate, Eq. (9)
becomes zero, since the atomic orbitals $\alpha$, $\beta$ are independent of
$q_{c}$. A detailed derivation of all ”new” terms in Eq. (8) for
$\zeta_{I}=R_{i}$ and $\zeta_{I}=q_{c}$ can be found in section S1 of the
supporting informations. All derivatives for
$\bigl{\langle}\hat{h}\bigl{\rangle}$, $\bigl{\langle}\hat{g}\bigl{\rangle}$,
and $V_{nn}$ with respect to nuclear coordinates are well known and can be
found elsewhere in the literature, for example 42. The resulting analytic CBO-
HF gradient $\bm{g}_{CBO}$ is a $\left(3N_{A}+1\right)$ vector with $N_{A}$
being the number of atoms in the molecule:

$\bm{g}_{CBO}=\begin{pmatrix}\frac{\partial}{\partial
R_{1}}&\frac{\partial}{\partial R_{2}}&\ldots&\frac{\partial}{\partial
q_{c}}\end{pmatrix}^{T}\bigl{\langle}E_{CBO}\bigr{\rangle}$ (10)

###  2.3 CBO-Hartree-Fock Frequencies

In order to calculate frequencies in the harmonic approximation for the CBO-HF
ansatz we need to calculate the Hessian matrix, $\bm{H}$, with dimensions
$\left(3N_{A}+1\right)\text{x}\left(3N_{A}+1\right)$ for the case of a single
photon mode. Each element $\bm{H}_{ij}$ is the second derivative of
$\bigl{\langle}E_{CBO}\bigr{\rangle}$ with respect to nuclear coordinates
respectively the photon displacement coordinate 30 that is

$\bm{H}_{ij}=\bm{H}_{ji}=\frac{\partial^{2}}{\partial\zeta_{i}\partial\zeta_{j}}\bigl{\langle}E_{CBO}\bigr{\rangle}.$
(11)

Using the analytic CBO-HF gradient $\bm{g}_{CBO}$ the second derivatives can
be determined using finite differences:

$\bm{H}_{ij}=\frac{\partial}{\partial\zeta_{i}}\bm{g}_{cbo}^{j}\approx\frac{\bm{g}_{cbo}^{j}\left(\zeta_{i}+\Delta\right)-\bm{g}_{cbo}^{j}\left(\zeta_{i}-\Delta\right)}{2\Delta}$
(12)

The harmonic frequencies and corresponding normal modes can be obtained by
transforming the Hessian matrix $\bm{H}$ from the $\zeta_{i}$ coordinates into
their mass-weighted version $\bm{H}^{M}$ and subsequently solving the
eigenvalue problem:

$\bm{A}^{\dagger}\bm{H}^{M}\bm{A}=\bm{\mathcal{L}}$ (13)

Here $\bm{\mathcal{L}}$ is the diagonal matrix of eigenvalues $\ell_{i}$ and
$\bm{A}$ is the matrix that diagonalizes the mass-weighted Hessian
$\bm{H}^{M}$. $\bm{A}$ is formed by the juxtaposition of its eigenvectors
$\bm{a}_{i}$ which define the normal modes. From the $\left(3N_{A}+1\right)$
eigenvalues, only $\left(3N_{A}+1-6\right)$, or $\left(3N_{A}+1-5\right)$ if
the molecule is linear, correspond to the harmonic frequencies of vibrational
motions, whereas the others correspond to translations and rotations. The
harmonic frequencies, in , are given by:

${}^{H}\nu_{i}=\sqrt{\frac{\ell_{i}}{4\pi^{2}c^{2}}}$ (14)

The corresponding intensities ${}^{H}\mathcal{I}_{i}$ in the harmonic
approximation are calculated as the projection of the dipole moment gradient
on the normal mode vectors $\bm{a}_{i}$:

${}^{H}\mathcal{I}_{i}=\left(\bm{\nabla}\bigl{\langle}\bm{\hat{\mu}}\bigr{\rangle}\cdot\bm{a}_{i}\right)^{2}$
(15)

The necessary dipole moment gradient
$\bm{\nabla}\bigl{\langle}\bm{\hat{\mu}}\bigr{\rangle}$ is calculated via
finite differences. To validate the quality of the harmonic approximation and
to go beyond the semi-classical treatment we determined the vibartional
spectra of the cavity coupled systems by then also determining the quantum
mechanical wavefunction of the nuclei and photons in the Born-Huang expansion.
Therefore, we not only have access to the fundamental transitions but also to
overtones. The anharmonic frequency ${}^{A}\nu_{j}$ for a given fundamental
transition is the difference of the energies of the eigenstate $j$ and the
nuclear-photonic ground state. The corresponding intensity
${}^{A}\mathcal{I}_{j}$ are calculated with the nuclear-photonic
eigenfunctions $\chi_{i}$ of the coupled molecular-cavity system:

${}^{A}\mathcal{I}_{0j}=\bigl{\langle}\chi_{0}\big{|}\bm{\hat{\mu}}\big{|}\chi_{j}\bigr{\rangle}$
(16)

##  3 Computational Details

The analytic gradients for the CBO-HF ansatz have been implemented in the
Psi4NumPy environment 44, which is an extension of the PSI4 45 electronic
structure package. All calculations were performed using the aug-cc-pVDZ basis
set 46 and the geometry of the isolated single HF and NH3 molecule have been
optimized at the Hartree-Fock level of theory. Consistent with our previous
studies 18, 19, the optimized structure of a single HF molecule is used to
define small ensembles of up to four molecules to study the collective effects
on the infrared spectra of molecules in cavities. These ensembles are
constructed by placing replicas of the single HF molecule separated by
$800\text{\,}\mathrm{\text{Å}}$ inside the cavity. Two different orientations
of the molecular HF ensembles are studied: The all-parallel configuration, in
which all HF molecules are aligned parallel and the antiparallel
configuration, which is characterized by a pairwise antiparallel orientation
of the $N_{mol}$ HF molecules both with respect to the cavity polarization
axis of the cavity. For NH3 only the single molecule case is studied. As
discussed in the supporting information of Ref. 19, the aligned orientation of
HF is not the most energetically favorable when coupled to a cavity, but a
transition state along the rotation with respect to the polarization axis.
However, the two corresponding minima lead to a situation where the molecule
is not coupled to the cavity field. For the case of NH3 the reoptimization in
the presence of the cavity is mainly characterized by a rotation in the
laboratory frame, resulting in an orientation with reduced coupling. Since the
effects on the internal coordinates of HF and NH3 are negligible for the
coupling strengths studied, we do not re-optimize the geometries of the
molecular systems in the cavity and align the molecules for maximum coupling
to the cavity. In all CBO-HF calculations performed in this work, we consider
a single mode, lossless cavity. We keep the collective coupling strength
$\bm{\lambda}_{c}$ constant by applying a scaling factor of
${1}/{\sqrt{N_{mol}}}$ to obtain a fixed Rabi splitting for different ensemble
sizes and treat $\lambda_{0}$ as a tunable coupling parameter:

$\bm{\lambda}_{c}=\frac{\lambda_{0}}{\sqrt{N_{mol}}}\bm{e}$ (17)

Here $\lambda_{0}$ is equivalent to $\lambda_{c}$ in Eq. (3) in the single
molecule case. As a result, we increase the mode volume $V_{c}$ of the cavity,
but by including more molecules, we keep the average density of molecules
$N_{mol}/V_{c}$ fixed. We use an artificially increased coupling strength
$\lambda_{0}$ in the range of 0.004 to 0.04, which corresponds to effective
mode volumes, see Eq. (3), in the single-molecule case as large as
$125.27\text{\,}{\mathrm{nm}}^{3}$ (for $\lambda_{0}=0.004$) or as small as
$1.25\text{\,}{\mathrm{nm}}^{3}$ (for $\lambda_{0}=0.04$).

To validate the harmonic approximation, the vibrational spectra of a single HF
molecule and a pair of HF molecules coupled to a cavity are calculated in a
fully quantum mechanical way. The necessary nuclear-photonic eigenfunctions
$\chi_{i}$ of the coupled molecular-cavity system are obtained by using the
imaginary time propagation method 47 with the Arnoldi propagation scheme 48
both implemented in our in-house quantum dynamics code (QDng). The used CBO-HF
cavity potential energy surfaces (cPESs) are represented by a two-dimensional
or three-dimensional numerical grid with 128 points along the bond lengths $R$
and 64 points along the photon displacement coordinate $q_{c}$. For all
vibrational spectra shown in this work the underling signals are broadened by
a Lorentzian function with the width of $10\text{\,}{\mathrm{cm}}^{-1}$. All
calculations were performed in a reproducible environment using the Nix
package manager together with NixOS-QChem 49 (commit f5dad404) and Nixpkgs
(nixpkgs, 22.11, commit 594ef126).

##  4 Results and Discussion

###  4.1 Validation of the Harmonic Approximation

Let us start with the validation of the harmonic approximation for vibro-
polaritonic IR spectra of a single HF molecule coupled to an optical cavity.
The spectra of the bare HF molecule calculated in the harmonic approximation
and fully quantum mechanical (anharmonic) are shown in Fig. S1 in the
supporting informations. The fundamental vibrational transition in the
harmonic approximation has a frequency of
${}^{H}\nu_{1}=$4467\text{\,}{\mathrm{cm}}^{-1}$$. The fully quantum
mechanical treatment leads to an expected strong red-shifted of
$186\text{\,}{\mathrm{cm}}^{-1}$ for the fundamental vibrational transition
(${}^{A}\nu_{1}=$4281\text{\,}{\mathrm{cm}}^{-1}$$). The vibro-polaritonic
infrared spectra for different coupling strengths calculated in the harmonic
approximation and full-quantum (anharmonic) as well as the trends in the Rabi
splitting $\Omega_{R}$ are shown in Fig. 1.

Figure 1: a) Vibro-polaritonic IR spectra of a single HF molecule calculated
in the harmonic approximation (redish) and in the full-qunatum/anharmonic
setup (bluish), individually normalized for each coupling strength. Black
dashed-dotted lines indicate the bare molecular frequencies of the harmonic
(${}^{H}\nu_{1}=$4467\text{\,}{\mathrm{cm}}^{-1}$$) and anharmonic
(${}^{A}\nu_{1}=$4281\text{\,}{\mathrm{cm}}^{-1}$$) fundamental transitions.
The cavity frequency $\omega_{c}$ is resonant with the corresponding
fundamental transition, and the coupling strength $\lambda_{c}$ is increased
from $0.009\text{\,}\mathrm{\text{au}}$ to $0.039\text{\,}\mathrm{\text{au}}$
(from lightest to darkest color). The Rabi splitting $\Omega_{R}$ (b) and its
asymmetry $\Delta\Omega_{R}=\omega_{c}-0.5\left(\nu^{LP}+\nu^{UP}\right)$ (c)
as a function of $\lambda_{c}$.

The resonant coupling of the cavity mode to the fundamental transition leads
to the expected formation of a LP transition and UP transition in both the
harmonic approximation and full-anharmonic treatment, see Fig. 1 a). The same
anharmonic shift as in the cavity-free spectra is observed in the vibro-
polaritonic spectra due to the anharmonicity along the nuclear coordinate
independent of the used coupling strength $\lambda_{c}$. In addition to the
red shift, the harmonic and full-anharmonic spectra differ mainly in their
intensity patterns. In the harmonic approximation, the UP transition is more
intense and the difference between LP and UP is smaller compared to the full-
quantum case. For the latter, the LP transition is more intense. This
discrepancy is probably due to the fact that the intensities in the harmonic
approximation are only calculated as the first derivative of the dipole
moment. The increase in $\lambda_{c}$ leads to a linear increase of the Rabi
splitting $\Omega_{R}$ between LP and UP which is identical for both the
harmonic and full-anharmonic spectra (Fig. 1 b)). Consistent with our previous
work 18, we observe an asymmetry
$\Delta\Omega_{R}=\omega_{c}-0.5\left(\nu^{LP}+\nu^{UP}\right)$ in the Rabi
splitting (Fig. 1 c)), where the LP is stronger red-shifted than the UP with
respect to $\omega_{c}$. This asymmetry is observed for the harmonic
approximation and the full-anharmonic treatment and in both cases
$\Delta\Omega_{R}$ increases quadratically with $\lambda_{c}$. However, the
asymmetry is more pronounced in the full-anharmonic case. This observed
asymmetry of $\Omega_{R}$ can also be understood as an detuning and change in
the optical length of the cavity when interacting with the molecule. Similar
results are found for two HF molecules in the all-parallel and antiparallel
configuration coupled to an optical cavity. The corresponding vibro-
polaritonic spectra for both case are shown in Fig. S2 and Fig. S3 of the
supporting informations. Due to the rescaling of $\bm{\lambda}_{c}$ (see Eq.
(17)) the spectra for the two molecule case look very similar to the single
molecule case and both configurations are nearly indistinguishable in their
spectra.

To further extend the validation of the harmonic approximation, we next
compare the effect of detuning the cavity frequency with respect to the
fundamental transition of the single HF molecule on the vibro-polaritonic IR
spectra. Besides the spectral properties, the harmonic approximation gives
access to normal modes $\bm{a}_{i}$. In the standard Born-Oppenheimer
approximation (BOA), these normal mode vectors describe the displacement of
classical nuclei associated with the corresponding vibrational transition. In
the CBOA for the case of single photon mode these vectors have an additional
term $a_{c}$ describing the change in the classical photon displacement field
$q_{c}$. The value of $|a_{c}|^{2}$ for a given normal mode is a measure of
how strongly the corresponding vibrational transition interacts with the
photon field. For an uncoupled light-matter system a pure molecular transition
is characterized by an $|a_{c}|^{2}$ value of zero, while the pure photon mode
has a value of one. Note that due to the length gauge $q_{c}$ and the
corresponding value $a_{c}$ are no longer a pure photonic quantity if light
and matter are coupled 40, 34. Nevertheless, $|a_{c}|^{2}$ can still be used
as a probe to identify how photonic the corresponding vibrational transition
is. The obtained information is comparable with the coefficients in Hopfield
models 50. In Fig. 2 the Rabi splitting $\Omega_{R}$ in harmonic
approximation, the magnitude $|a_{c}|^{2}$ and the cavity-induced energy
changes $E_{lin}$, $E_{dse}$ and $E_{dis}$ are plotted as function of the
cavity frequency $\omega_{c}$ for a fixed coupling strength.

Figure 2: a) Difference between the Rabi splitting $\Omega_{R}$ and
$\Delta=\omega_{c}-^{H}\nu_{1}$ as a function of the cavity frequency
$\omega_{c}$ for a single HF molecule using the CBO-HF harmonic approximation.
Here $\Delta$ describes the difference between the cavity frequency
$\omega_{c}$ and the fundamental bare molecular transition ${}^{H}\nu_{1}$. b)
The normal mode value $|a_{c}|^{2}$ describing the change in $q_{c}$ for the
LP and UP vibrational transition as function of $\omega_{c}$. c) Energy
contributions due to the cavity interaction as function of the cavity
frequency $\omega_{c}$ for the CBO-HF ansatz. Black dashed-dotted lines
indicate the frequencies of the harmonic fundamental transition
(${}^{H}\nu_{1}=$4467\text{\,}{\mathrm{cm}}^{-1}$$). A constant coupling
strength $\lambda_{c}$ of $0.039\text{\,}\mathrm{\text{au}}$ is used.

For better visualization, the difference between $\Omega_{R}$ and
$\Delta=\omega_{c}-^{H}\nu_{1}$, which describes the detuning of the cavity
mode, is shown in Fig. 3 a). For an uncoupled molecular cavity system, this
difference would be zero and independent of $\omega_{c}$. Note that
$\Omega_{R}-\Delta$ is only a measure of the total size of $\Omega_{R}$, not
of the asymmetry $\Delta\Omega_{R}$, which is discussed in Fig. 1. As
expected, the largest Rabi splitting $\Omega_{R}$ is obtained for $\omega_{c}$
being resonant with the fundamental bare molecular transition
(${}^{H}\nu_{1}$). Interestingly, the difference between $\Omega_{R}$ and the
detuning $\Delta$ is not symmetric with respect to ${H}\nu_{1}$, and tends to
a finite non-zero value even for large detunings. These results clearly show
that the observation of a large Rabi splitting is a rather sharp resonance
effect. A similar resonance is visible for the normal mode values
$|a_{c}|^{2}$ of the vibrational transitions of the LP state and UP state, as
shown in Fig. 3 b). For $\omega_{c}$ equal to ${}^{H}\nu_{1}$ the
$|a_{c}|^{2}$ values of LP and UP are close to $0.5$, indicating that both
vibrational transitions have a hybrid light-matter character. In contrast to
the prominent resonance features observed in Fig. 3 a) and b), the cavity-
induced energy modifications $E_{lin}$, $E_{dse}$ and $E_{dis}$ shown in Fig.
3 c) are constant for all values of $\omega_{c}$. This result indicates that
the energy changes induced by the cavity do not depend on any resonance
condition, consistent with recent theoretical findings 51, 52, where the
coupling to a non-resonant cavity alters the ground state potential energy
surface (PES). Therefore, the experimentally observed change in chemical
reactivity near resonance conditions cannot be explained by simple shifts of
the underlying potential energy surfaces alone.

The direct comparison of the Rabi splitting and the light-matter character of
the LP and UP transitions between the harmonic and full anharmonic simulations
is shown in Fig. 3. Since we treat the coupled cavity-molecular system in full
quantum fashion, we have access to the nuclear-photonic eigenfunctions of the
LP state ($\chi_{LP}$) and the UP state ($\chi_{UP}$). Figs. S5 and S6 showing
examples of uncoupled and coupled nuclear-photonic eigenfunctions can be found
in section S3 of the supporting information. The coupled eigenfunctions can be
expanded in terms of the uncoupled pure matter eigenfunction
($\chi_{M}^{(0)}$) and the photonic eigenfunction ($\chi_{C}^{(0)}$):

$\begin{split}\chi_{LP}&=a_{M}^{LP}\chi_{M}^{(0)}+a_{c}^{LP}\chi_{c}^{(0)}\\\
\chi_{UP}&=a_{M}^{UP}\chi_{M}^{(0)}+a_{c}^{UP}\chi_{c}^{(0)}\end{split}$ (18)

The absolute square of obtained expansion coefficients $|a_{c}^{LP}|^{2}$ and
$|a_{c}^{UP}|^{2}$ can be used as a measure of the light-matter character
similar to their classical pendant normal mode values $|a_{c}|^{2}$. For a
better comparison, the results of the harmonic simulations are red-shifted by
$186.0\text{\,}{\mathrm{cm}}^{-1}$ and the black dashed dotted line in Fig. 3
indicates the anharmonic fundamental transition.

Figure 3: a) Difference between the Rabi splitting $\Omega_{R}$ and
$\Delta=\omega_{c}-^{H}\nu_{1}$ as a function of the cavity frequency
$\omega_{c}$ for a single HF molecule in the full anharmonic treatment (blue)
and the harmonic approximation (red). Here $\Delta$ describes the difference
between the cavity frequency $\omega_{c}$ and the fundamental bare molecular
transition ${}^{H}\nu_{1}$. b) Change of the cavity normal mode magnitude
$|a_{c}|^{2}$ (red) for LP and UP and the values $|a_{c}^{LP}|^{2}$ and
$|a_{c}^{UP}|^{2}$ of the one-photon cavity vacuum state (blue). The harmonic
results are shifted to match the anharmonic fundamental transition
(${}^{A}\nu_{1}=$4281\text{\,}{\mathrm{cm}}^{-1}$$), indicated by the black
dashed dotted line. A constant coupling strength $\lambda_{c}$ of
$0.039\text{\,}\mathrm{\text{au}}$ is used.

In good agreement with the previous discussion, the results obtained using the
harmonic approximation are very similar to the results of the full-anharmonic
simulations. After accounting for anharmonicity, the differences between
$\Omega_{R}$ and detuning $\Delta$ (Fig. 3 a)) are close to resonance nearly
identical and for large detunings qualitatively similar. The asymmetry
$\Delta\Omega_{R}$, as well as the asymmetry of the difference between
$\Omega_{R}$ and $\Delta$ with respect to the fundamental transition, is more
pronounced in the full-quantum simulation. Also, the $|a_{c}|^{2}$ values of
the LP and UP are beside the overall anharmonic-shift in very good agreement
between the harmonic approximation and the full anharmonic treatment. Overall,
we can conclude that for the coupling strengths studied, the semiclassical
harmonic approximation qualitatively reproduces the main features of the
vibro-polaritonic IR spectra very well, even for a very anharmonic molecule
like HF, except for the well-known general limitation of the harmonic
approximation 53, 54.

Of course, the normal mode vectors $\bm{a}_{i}$ obtained in the harmonic
approximation contain more information than only the $|a_{c}|^{2}$ value
discussed in this section. A detailed analysis of the whole $\bm{a}_{i}$
vectors and a comparison withe the nuclear-photonic eigenfunctions for the
case of one HF molecule and two HF molecules can be found in section S3 of the
supporting informations.

###  4.2 Influence of the SCF-Treatment and the Dipole Self-Energy on Vibro-
Polaritonic Spectra

As shown in our previous work 18, 19, both self-consistent treatment and the
consideration of the full DSE operator are crucial to capture relevant aspects
in the description of strongly coupled molecules and molecular ensembles.
Therefore, we want to determine the influence of both factors on the vibro-
polaritonic IR spectra for a single HF molecule and small ensembles with up to
four HF molecules. To analyze the influence of the SCF treatment, the vibro-
polaritonic IR spectra of a single HF molecule are calculated in the harmonic
approximation and full quantum (anharmonic), neglecting the SCF of the coupled
electronic photonic system. For the harmonic approximation, this is achieved
by using cavity-free dipole moments, density matrix elements
$D_{\alpha,\beta}$, and orbital energies $\epsilon_{i}$ obtained by a standard
Hartree-Fock calculation to determine the CBO-HF gradients and Hessians. In
the full-anharmonic treatment, we use cavity-free Hartree-Fock energies and
expectation values of dipole moments and DSE terms to construct the necessary
potential energy surfaces in a extended molecular Jaynes-Cummings model (EJCM)
55, 56, 57, 58. The vibro-polaritonic infrared spectra obtained for different
coupling strengths calculated in the harmonic approximation and the full
anharmonic as well as the trends in Rabi splitting $\Omega_{R}$ are shown in
Fig. 4.

Figure 4: Vibro-polaritonic IR spectra of a single HF molecule calculated in
the harmonic approximation (redish) and in the setup (bluish), individually
normalized for each coupling strength. In both cases the full SCF treatment is
neglected; for details, see text. Black dashed-dotted lines indicate the
frequencies of the harmonic fundamental transition
(${}^{H}\nu_{1}=$4467\text{\,}{\mathrm{cm}}^{-1}$$) and the anharmonic
fundamental transition (${}^{A}\nu_{1}=$4281\text{\,}{\mathrm{cm}}^{-1}$$).
The cavity frequency $\omega_{c}$ is resonant with the corresponding
fundamental transition in both cases, and the coupling strength $\lambda_{c}$
increases from $0.009\text{\,}\mathrm{\text{au}}$ to
$0.039\text{\,}\mathrm{\text{au}}$ (from light to dark color). b) Rabi
splitting $\Omega_{R}$ as a function of $\epsilon_{c}$. c) Asymmetry
$\Delta\Omega_{R}=\omega_{c}-0.5\left(\nu^{LP}+\nu^{UP}\right)$ of the Rabi
splitting.

Comparing the vibro-polaritonic IR spectra without SCF (Fig. 4 a) with those
shown in Fig. 1 a), both the harmonic and the full anharmonic spectra differ.
Especially for the harmonic spectra, SCF has a massive influence. Without a
fully converged CBO-HF the Rabi splitting $\Omega_{R}$ is larger, see Fig. 4
b), and nearly symmetric with respect to the cavity frequency, see Fig. 4 c).
Also, the intensity of the LP transition and the UP transition are almost
identical. For the full-anharmonic spectra in the EJCM model the discrepancies
seem smaller, but still clear differences can be recognized. The value of
$\Omega_{R}$ is reproduced quite well. However, its asymmetry
$\Delta\Omega_{R}$ is smaller and characterized by blue-shift of LP and UP
with respect to $\omega_{c}$. In addition, also the ratio of the intensities
is inverted.

Another important aspect in describing VSC for molecules is the DSE which
gives rise to a cavity-induced interaction between molecules in an ensemble
and is very sensitive to the molecular orientation 18, 19, 59. To quantify its
influence on spectral features, we calculated the vibro-polaritonic IR spectra
in the harmonic approximation for different numbers of HF molecules with and
without the DSE terms included in the underlying CBO-HF ansatz.In the
following, the case without the DSE terms is called liner CBO-HF. The results
for the all-parallel and the antiparallel configuration are shown in Fig. 5.
The observed Rabi splitting $\Omega_{R}$ and its asymmetry $\Delta\Omega_{R}$
for the all-parallel and the antiparallel configurations with and without the
DSE term are shown as a function of the number of molecules in Fig. 6

Figure 5:  Vibro-polaritonic IR spectra calculated in the harmonic
approximation for different numbers of HF molecules (color-coded) shown with
respect to the cavity frequency $\omega_{c}$. The cavity is resonant with the
harmonic fundamental transition
(${}^{H}\nu_{1}=$4467\text{\,}{\mathrm{cm}}^{-1}$$, black dashed-dotted lines)
and a rescaled coupling strength of $\lambda_{0}$ of
$0.057\text{\,}\mathrm{\text{au}}$ is used (see Eq. (17)). a) Full CBO-HF
simulation in the all-parallel configuration. b) Full CBO-HF simulation in the
antiparallel configuration. c) Linear CBO-HF simulation (without DSE terms) in
the all-parallel configuration. c)Linear CBO-HF simulation (without DSE terms)
in the antiparallel configuration. Figure 6:  Rabi splitting $\Omega_{R}$ and
asymmetry $\Delta\Omega_{R}=\omega_{c}-0.5\left(\nu^{LP}+\nu^{UP}\right)$ of
the Rabi splitting as a function of $N_{Mol}$. Both calculated in the harmonic
approximation applying full CBO-HF and linear CBO-HF (without DSE terms) for
the all-parallel configuration shown in a) and c) and in the antiparallel
orientation shown in b) and d). The cavity frequency $\omega_{c}$ is resonant
with the harmonic fundamental transition ($4467\text{\,}{\mathrm{cm}}^{-1}$)
and a rescaled coupling strength of $\lambda_{0}$ of
$0.057\text{\,}\mathrm{\text{au}}$ is used.

The full CBO-HF spectra of the all-parallel configuration (Fig. 5 a)) and the
antiparallel configuration (Fig. 5 b)) are indistinguishable. In line with our
recent work 18 increasing the number of molecules ($N_{Mol}$) in the cavity
while keeping the density constant (rescaling of the coupling strength, see
Eq. (17)) leads to an increasing asymmetry $\Delta\Omega_{R}$ in the
splitting, see Fig. 6 c) and d) red lines. At the same time, the Rabi
splitting itself decreases, see Fig. 6 a) and b) red lines, although the
applied rescaling should keep it constant. This may be an effect of the rather
large coupling strength applied. For weaker coupling, the applied rescaling of
the coupling strength leads to an almost constant value of $\Omega_{R}$ as the
number of molecules increases, see Ref. 18. If the DSE contribution is
neglected the two studied configurations become distinguishable, see Fig. 5 c)
and d). Comparing the all-parallel spectra with and without DSE, the
differences are rather small. The observed Rabi splitting $\Omega_{R}$ is
slightly larger in the case without DSE and independent of the number of
molecules, as shown in Fig. 6 a) blue line. Furthermore, the asymmetry
$\Delta\Omega_{R}$ is also independent of $N_{Mol}$ and constant (Fig. 6 c)
blue line). The observed changes are more drastic in the case of the
antiparallel configuration. In Fig. 5 d) an oscillatory shifting of the LP
signal and UP signal is observed while changing the number of molecules.
Consequently, this oscillation can be found in the Rabi splitting $\Omega_{R}$
and its asymmetry, see Fig. 6 b) and d) blue line. The oscillating pattern can
be attributed to the fact that odd or even number of molecules in the
antiparallel configuration create two different situations. For even $N_{Mol}$
the whole ensemble reduces to an effective antiparallel bimolecular case, and
the situation of odd $N_{Mol}$ is equivalent to the case of a single molecule.
Consistent with the literature 40, 34, 17, 18, 60, 19 these results clearly
show how important the DSE terms and the SCF treatment are to obtain a
physically meaningful result.

###  4.3 Beyond Diatomic Molecules: Vibro-Polaritonic Spectra of NH3

So far we have discussed only diatomic molecules, which have only one
vibrational degree of freedom, or, in the case of the ensembles, linear
combinations of the same vibrational degree. However, since we have validated
the harmonic approximation and the associated normal mode analysis, its semi-
classical nature also allows us to treat larger molecular systems. Therefore,
in the last part of this work we will study the vibro-polaritonic IR spectra
of a single NH3 molecule. The bare NH3 molecule has six vibrational degrees of
freedom: One symmetric bending mode ($\nu_{1}$), two degenerate asymmetric
bending modes ($\nu_{2,3}$), one symmetric stretching mode ($\nu_{4}$) and two
degenerate asymmetric stretching modes ($\nu_{5,6}$). Since NH3 is not a
linear molecule and the vibrational modes have different spatial orientation
(symmetry), the polarization axis of the single cavity mode plays a
significant role. In our simulations, the single NH3 molecule oriented with
respect to the center of nuclear charges and the lone pair orbital of nitrogen
is pointing along the $z$ axis of the laboratory frame. Since for this chosen
molecular orientation the $x$ axis and the $y$ axis are equivalent, only the
$y$ polarization and the $z$ polarization of the single cavity mode are
discussed in the following. The cavity frequency $\omega_{c}$ is set to be
resonant with the symmetric bending mode $\nu_{1}$, which has the lowest
energy of $1103\text{\,}{\mathrm{cm}}^{-1}$ and is the strongest transition in
the cavity-free case. The bare molecular vibronic IR spectrum and the vibro-
polaritonic IR spectra of a single NH3 molecule are shown in Fig. 7 for the
bending modes and in Fig. S14 of the supporting information for the stretching
modes.

Figure 7: a) Low energy part of the vibronic IR spectra of a single NH3
molecule calculated in the harmonic approximation. The low energy part of the
vibro-polaritonic IR spectra of a single NH3 molecule zoomed into the
symmetric mode (b) and d)) and the two asymmetric bending modes shown in c)
and e). The polarization axis of the cavity mode is the $y$ axis for b) and c)
and equal to the $z$ axis for d) and e). The cavity frequency $\omega_{c}$ is
resonant with the symmetric bending mode ($1103\text{\,}{\mathrm{cm}}^{-1}$)
and the cavity field strength $\lambda_{c}$ increases from
$0.039\text{\,}\mathrm{\text{au}}$ to $0.155\text{\,}\mathrm{\text{au}}$.

As expected, the splitting into a LP and a UP transition is observed for the
symmetric bending mode $\nu_{1}$ when coupled to a resonant optical cavity,
see Fig. 7 b) and d). The Rabi splitting $\Omega_{R}$ increases with
increasing coupling strength and is larger for $z$ polarization than for $x/y$
polarization due to the chosen molecular orientation. Interestingly, the
signal of the asymmetric bending modes ($\nu_{2}$ and $\nu_{3}$) is also
affected by the coupling to the cavity mode, as shown in Fig. 7 c) and e). For
the $x/y$ polarization the signal is blue-shifted and for a coupling strength
$\lambda_{c}$ larger than $0.078\text{\,}\mathrm{\text{au}}$ the degeneracy of
the two transitions is lifted. This shift and the splitting into two separate
signals is also observed for the $z$ polarization. However, it is smaller for
this polarization and occurs only at higher coupling strengths.

To better understand and characterize the two observed splittings in the
vibro-polaritonic IR spectra of NH3, their magnitude is plotted for both
polarization directions as a function of the coupling strength $\lambda_{c}$
in Fig. 8 a) and b). Furthermore, the corresponding values $|a_{c}|^{2}$ for
the associated normal modes for both splittings are plotted in Fig. 8 c) and
d) as a function of $\lambda_{c}$. These $|a_{c}|^{2}$ values describe the
change in the classical photon displacement field $q_{c}$ and probes how
photonic the corresponding vibrational transition is.

Figure 8: a) Rabi splitting $\Omega_{R}$ for the symmetric bending mode of a
single NH3 molecule as a function of $\lambda_{c}$. b) Splitting $\Delta\nu$
of the asymmetric bending modes due to cavity interaction. The normal mode
value $|a_{c}|^{2}$ describing the change in $q_{c}$ for the UP and LP modes
c) and for the two asymmetric bending modes d). The cavity frequency
$\omega_{c}$ is resonant with the symmetric bending mode
($1103\text{\,}{\mathrm{cm}}^{-1}$) and both $y/x$ (bluish) and $z$ (reddish)
polarization of the cavity mode are shown.

Clearly, the formation of a LP and a UP transition for the symmetric bending
mode is more efficient for the $z$ polarization since the observed Rabi
splitting $\Omega_{R}$ is larger and scales almost linear with $\epsilon_{c}$,
see Fig. 8 a). For the $y/x$ polarization the value is smaller but the scaling
is close to quadratic. The stronger effect observed for $z$ polarization is
understandable, since the symmetric bending mode mainly involves motion along
the $z$ axis in the laboratory frame and only very small changes along the
$y/x$ direction. A closer look at the $|a_{c}|^{2}$ values for the LP and UP
normal modes (Fig. 8 c) shows a similar behavior to the case of the HF
molecule. Once the cavity coupling is present, the $|a_{c}|^{2}$ values are
clearly non-zero for both LP and UP. Since for the $z$ polarization the
coupling between the cavity mode and the symmetric bending mode is more
efficient, the $|a_{c}|^{2}$ value of the UP mode is closer to that of the LP
mode and clearly larger compared to the $y/x$ case. In other words, the light
and the molecule are less hybridized for the $y/x$ polarization. The splitting
and the $|a_{c}|^{2}$ values for the two asymmetric bending modes shown in
Fig. 8 b) and d) behave differently. In the $z$ polarization, where a smaller
splitting is observed, both $|a_{c}|^{2}$ values are lower than $0.001$. Even
for the $y/x$ case, they do not become larger than $0.002$, which is still
orders of magnitude smaller than for the LP and the UP transition. This is
clear evidence that both asymmetric bending modes mostly are matter
transitions and no light-matter hybrid states are formed. The observed
splitting of $\nu_{2}$ and $\nu_{3}$ when coupled to a cavity can be explained
by a symmetry breaking by the confined photon mode. Its polarization axis in
$x$ or $y$ makes these two directions, which are identical in the case without
cavity, distinguishable. Note that in our current simulations we only consider
a single photon mode, and in a more general setup, there are many modes with
usually more than one polarization direction. Therefore, we assume that the
observed splitting due to symmetry breaking is overestimated in our
calculations. However, the fact that the interaction with an optical cavity
not only affects near-resonance transitions by forming LP and UP states, but
also modifies all signals by slightly changing the energetics and breaking
molecular symmetries is probably still important in more general cases and
already observed earlier in literature 30.

##  5 Summary and Conclusion

Based on the recently formulated cavity Born-Oppenheimer Hartree-Fock ansatz
19 we have introduced a wave function-based methodology to calculate the
vibro-polaritonic IR spectra in a ab-initio manner. By applying the well-known
harmonic approximation, we have access to the vibrational frequencies and
normal modes of systems when light and matter are strongly coupled. Using the
cavity Born-Oppenheimer approximation, the obtained normal modes combine both
nuclear and photonic degrees of freedom, allowing a detailed analysis of the
vibro-polariton states. The necessary second derivatives of the CBO-HF energy
are calculated via finite differences of the analytic first derivatives
(gradients) introduced in this work.

We demonstrate the capability of our framework by carefully comparing the
vibro-polaritonic IR spectra obtained for small ensembles of HF molecules with
full quantum-mechanical calculations. Overall, the semi-classical harmonic
approximation qualitatively reproduces the main features of the vibro-
polaritonic IR spectra very well. It even captures to some extent polarization
effects such as cavity detuning previously observed 18 when the full DSE
contribution is included and the coupled electronic-photonic problem is solved
self-consistently. The semi-classical nature of the harmonic approximation
allows efficient description of vibro-polaritonic IR spectra of large
molecular systems or small ensembles. As a first test case, we simulated the
vibro-polaritonic spectra of a single NH3 molecule. In addition to the
expected formation of vibro-polaritonic states, the coupling to an optical
cavity changes the whole vibrational spectrum. The presence of the confined
light mode leads to an energy shift of most signals and, depending on the
cavity polarization and the molecular orientation, the molecular symmetry is
reduced.

The new analytical gradients used here to obtain vibro-polaritonic IR spectra
open many avenues for future exploration. The gradients and the numerical
Hessian can be used to optimize the molecular system coupled to an optical
cavity. In contrast to the cavity/field-free case, rotation with respect to
the laboratory frame is no longer trivial for optimizations in the presence of
a cavity mode with a defined polarization axis 19. Beside the possibility to
optimize molecules under VSC, the analytical gradients can be used to perform
ab-initio semi-classical dynamics simulations for single molecules or small
ensembles.

We thank Dominik Sidler and Michael Ruggenthaler for inspiring discussions and
helpful comments. M.K. acknowledges funding from the European Union’s Horizon
2020 research and innovation program under the Marie Skłodowska-Curie grant
agreement No. 860553. This project has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (grant agreement no. 852286).

See the supplementary material for the details of the derivation of the CBO-HF
energy gradients, additional figures for the discussion of HF and NH3 coupled
to an optical cavity and a detailed comparison of nuclear-photonic
eigenfunctions and vibro-polaritonic normal modes. All data underlying this
study are available from the corresponding author upon reasonable request.

## References

  * Hutchison et al. 2012 Hutchison, J. A.; Schwartz, T.; Genet, C.; Devaux, E.; Ebbesen, T. W. Modifying chemical landscapes by coupling to vacuum fields. _Angew. Chem. Int. Ed Engl._ 2012, _51_ , 1592–1596
  * George et al. 2016 George, J.; Chervy, T.; Shalabney, A.; Devaux, E.; Hiura, H.; Genet, C.; Ebbesen, T. W. Multiple Rabi Splittings under Ultrastrong Vibrational Coupling. _Phys. Rev. Lett._ 2016, _117_ , 153601
  * Thomas et al. 2016 Thomas, A.; George, J.; Shalabney, A.; Dryzhakov, M.; Varma, S. J.; Moran, J.; Chervy, T.; Zhong, X.; Devaux, E.; Genet, C.; Hutchison, J. A.; Ebbesen, T. W. Ground-State Chemical Reactivity under Vibrational Coupling to the Vacuum Electromagnetic Field. _Angew. Chem. Int. Ed Engl._ 2016, _55_ , 11462–11466
  * Thomas et al. 2019 Thomas, A.; Lethuillier-Karl, L.; Nagarajan, K.; Vergauwe, R. M. A.; George, J.; Chervy, T.; Shalabney, A.; Devaux, E.; Genet, C.; Moran, J.; Ebbesen, T. W. Tilting a ground-state reactivity landscape by vibrational strong coupling. _Science_ 2019, _363_ , 615–619
  * Hirai et al. 2020 Hirai, K.; Hutchison, J. A.; Uji-I, H. Recent Progress in Vibropolaritonic Chemistry. _Chempluschem_ 2020, _85_ , 1981–1988
  * Hirai et al. 2020 Hirai, K.; Takeda, R.; Hutchison, J. A.; Uji-I, H. Modulation of Prins Cyclization by Vibrational Strong Coupling. _Angew. Chem. Int. Ed Engl._ 2020, _59_ , 5332–5335
  * Kumar et al. 2023 Kumar, S.; Biswas, S.; Rashid, U.; Mony, K. S.; Vergauwe, R. M. A.; Kaliginedi, V.; Thomas, A. Extraordinary Electrical Conductance of Non-conducting Polymers Under Vibrational Strong Coupling. 2023,
  * Ahn et al. 2023 Ahn, W.; Triana, J. F.; Recabal, F.; Herrera, F.; Simpkins, B. S. Modification of ground-state chemical reactivity via light-matter coherence in infrared cavities. _Science_ 2023, _380_ , 1165–1168
  * Bhuyan et al. 2023 Bhuyan, R.; Mony, J.; Kotov, O.; Castellanos, G. W.; Gómez Rivas, J.; Shegai, T. O.; Börjesson, K. The Rise and Current Status of Polaritonic Photochemistry and Photophysics. _Chem. Rev._ 2023,
  * Zhong et al. 2023 Zhong, C.; Hou, S.; Zhao, X.; Bai, J.; Wang, Z.; Gao, F.; Guo, J.; Zhang, F. Driving DNA Origami Coassembling by Vibrational Strong Coupling in the Dark. _ACS Photonics_ 2023,
  * Gu et al. 2023 Gu, K.; Si, Q.; Li, N.; Gao, F.; Wang, L.; Zhang, F. Regulation of Recombinase Polymerase Amplification by Vibrational Strong Coupling of Water. _ACS Photonics_ 2023,
  * Weight et al. 2023 Weight, B. M.; Krauss, T. D.; Huo, P. Investigating Molecular Exciton Polaritons Using Ab Initio Cavity Quantum Electrodynamics. _J. Phys. Chem. Lett._ 2023, _14_ , 5901–5913
  * Nagarajan et al. 2021 Nagarajan, K.; Thomas, A.; Ebbesen, T. W. Chemistry under Vibrational Strong Coupling. _J. Am. Chem. Soc._ 2021, _143_ , 16877–16889
  * Dunkelberger et al. 2022 Dunkelberger, A. D.; Simpkins, B. S.; Vurgaftman, I.; Owrutsky, J. C. Vibration-Cavity Polariton Chemistry and Dynamics. _Annu. Rev. Phys. Chem._ 2022, _73_ , 429–451
  * Simpkins et al. 2023 Simpkins, B. S.; Dunkelberger, A. D.; Vurgaftman, I. Control, Modulation, and Analytical Descriptions of Vibrational Strong Coupling. _Chem. Rev._ 2023,
  * Campos-Gonzalez-Angulo et al. 2022 Campos-Gonzalez-Angulo, J. A.; Poh, Y. R.; Du, M.; Yuen-Zhou, J. Swinging between shine and shadow: Theoretical advances on thermally-activated vibropolaritonic chemistry (a perspective). 2022,
  * Sidler et al. 2022 Sidler, D.; Ruggenthaler, M.; Schäfer, C.; Ronca, E.; Rubio, A. A perspective on ab initio modeling of polaritonic chemistry: The role of non-equilibrium effects and quantum collectivity. _J. Chem. Phys._ 2022, _156_ , 230901
  * Sidler et al. 2023 Sidler, D.; Schnappinger, T.; Obzhirov, A.; Ruggenthaler, M.; Kowalewski, M.; Rubio, A. Unraveling a cavity induced molecular polarization mechanism from collective vibrational strong coupling. 2023,
  * Schnappinger et al. 2023 Schnappinger, T.; Sidler, D.; Ruggenthaler, M.; Rubio, A.; Kowalewski, M. Cavity Born-Oppenheimer Hartree-Fock Ansatz: Light-Matter Properties of Strongly Coupled Molecular Ensembles. _J. Phys. Chem. Lett._ 2023, _14_ , 8024–8033
  * Li et al. 2021 Li, X.; Mandal, A.; Huo, P. Cavity frequency-dependent theory for vibrational polariton chemistry. _Nat. Commun._ 2021, _12_ , 1315
  * Schäfer et al. 2022 Schäfer, C.; Flick, J.; Ronca, E.; Narang, P.; Rubio, A. Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity. _Nature Communications_ 2022, _13_ , 7817
  * Lindoy et al. 2023 Lindoy, L. P.; Mandal, A.; Reichman, D. R. Quantum dynamical effects of vibrational strong coupling in chemical reactivity. _Nature Communications_ 2023, _14_ , 2733
  * Schütz et al. 2020 Schütz, S.; Schachenmayer, J.; Hagenmüller, D.; Brennen, G. K.; Volz, T.; Sandoghdar, V.; Ebbesen, T. W.; Genes, C.; Pupillo, G. Ensemble-Induced Strong Light-Matter Coupling of a Single Quantum Emitter. _Phys. Rev. Lett._ 2020, _124_ , 113602
  * Davidsson and Kowalewski 2023 Davidsson, E.; Kowalewski, M. The role of dephasing for dark state coupling in a molecular Tavis-Cummings model. 2023,
  * Davidsson and Kowalewski 2020 Davidsson, E.; Kowalewski, M. Atom Assisted Photochemistry in Optical Cavities. _J. Phys. Chem. A_ 2020, _124_ , 4672–4677
  * Flick and Narang 2018 Flick, J.; Narang, P. Cavity-Correlated Electron-Nuclear Dynamics from First Principles. _Phys. Rev. Lett._ 2018, _121_ , 113002
  * Schäfer et al. 2022 Schäfer, C.; Flick, J.; Ronca, E.; Narang, P.; Rubio, A. Shining light on the microscopic resonant mechanism responsible for cavity-mediated chemical reactivity. _Nat. Commun._ 2022, _13_ , 7817
  * Yu 2022 Yu, Q. Multidimensional quantum calculation of the infrared spectra under polaritonic vibrational strong and ultrastrong coupling. 2022,
  * Mondal et al. 2023 Mondal, M. E.; Koessler, E. R.; Provazza, J.; Vamivakas, A. N.; Cundiff, S. T.; Krauss, T. D.; Huo, P. Quantum dynamics simulations of the 2D spectroscopy for exciton polaritons. _J. Chem. Phys._ 2023, _159_
  * Bonini and Flick 2022 Bonini, J.; Flick, J. Ab Initio Linear-Response Approach to Vibro-Polaritons in the Cavity Born-Oppenheimer Approximation. _J. Chem. Theory Comput._ 2022,
  * Flick et al. 2017 Flick, J.; Ruggenthaler, M.; Appel, H.; Rubio, A. Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry. _Proc. Natl. Acad. Sci. U.S.A._ 2017, _114_ , 3026–3034
  * Flick et al. 2017 Flick, J.; Appel, H.; Ruggenthaler, M.; Rubio, A. Cavity Born–Oppenheimer approximation for correlated electron–nuclear-photon systems. _J. Chem. Theory Comput._ 2017, _13_ , 1616–1625
  * Fischer and Saalfrank 2023 Fischer, E. W.; Saalfrank, P. Vibro-Polaritonic Chemistry beyond the Cavity Born-Oppenheimer Approximation and the Role of Electron-Photon Correlation. 2023,
  * Schäfer et al. 2020 Schäfer, C.; Ruggenthaler, M.; Rokaj, V.; Rubio, A. Relevance of the Quadratic Diamagnetic and Self-Polarization Terms in Cavity Quantum Electrodynamics. _ACS Photonics_ 2020, _7_ , 975–990
  * Ruggenthaler et al. 2023 Ruggenthaler, M.; Sidler, D.; Rubio, A. Understanding Polaritonic Chemistry from Ab Initio Quantum Electrodynamics. _Chem. Rev._ 2023,
  * Spohn 2004 Spohn, H. _Dynamics of charged particles and their radiation field_ ; Cambridge university press, 2004
  * Ruggenthaler et al. 2018 Ruggenthaler, M.; Tancogne-Dejean, N.; Flick, J.; Appel, H.; Rubio, A. From a quantum-electrodynamical light–matter description to novel spectroscopies. _Nature Reviews Chemistry_ 2018, _2_ , 1–16
  * Jestädt et al. 2019 Jestädt, R.; Ruggenthaler, M.; Oliveira, M. J.; Rubio, A.; Appel, H. Light-matter interactions within the Ehrenfest–Maxwell–Pauli–Kohn–Sham framework: fundamentals, implementation, and nano-optical applications. _Adv. Phys._ 2019, _68_ , 225–333
  * Schäfer et al. 2018 Schäfer, C.; Ruggenthaler, M.; Rubio, A. Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling. _Phys. Rev. A_ 2018, _98_ , 043801
  * Rokaj et al. 2018 Rokaj, V.; Welakuh, D. M.; Ruggenthaler, M.; Rubio, A. Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy. _J. Phys. B At. Mol. Opt. Phys._ 2018, _51_ , 034005
  * Szabo and Ostlund 1996 Szabo, A.; Ostlund, N. S. _Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory_ , 1st ed.; Dover Publications, Inc.: Mineola, 1996
  * Jensen 2017 Jensen, F. _Introduction to computational chemistry_ , 3rd ed.; John Wiley & Sons: Nashville, TN, 2017
  * Pulay 1969 Pulay, P. Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules. _Mol. Phys._ 1969, _17_ , 197–204
  * Smith et al. 2018 Smith, D. G. A.; Burns, L. A.; Sirianni, D. A.; Nascimento, D. R.; Kumar, A.; James, A. M.; Schriber, J. B.; Zhang, T.; Zhang, B.; Abbott, A. S.; Berquist, E. J.; Lechner, M. H.; Cunha, L. A.; Heide, A. G.; Waldrop, J. M.; Takeshita, T. Y.; Alenaizan, A.; Neuhauser, D.; King, R. A.; Simmonett, A. C.; Turney, J. M.; Schaefer, H. F.; Evangelista, F. A.; DePrince, A. E., 3rd; Crawford, T. D.; Patkowski, K.; Sherrill, C. D. Psi4NumPy: An Interactive Quantum Chemistry Programming Environment for Reference Implementations and Rapid Development. _J. Chem. Theory Comput._ 2018, _14_ , 3504–3511
  * Smith et al. 2020 Smith, D. G. A.; Burns, L. A.; Simmonett, A. C.; Parrish, R. M.; Schieber, M. C.; Galvelis, R.; Kraus, P.; Kruse, H.; Di Remigio, R.; Alenaizan, A.; James, A. M.; Lehtola, S.; Misiewicz, J. P.; Scheurer, M.; Shaw, R. A.; Schriber, J. B.; Xie, Y.; Glick, Z. L.; Sirianni, D. A.; O’Brien, J. S.; Waldrop, J. M.; Kumar, A.; Hohenstein, E. G.; Pritchard, B. P.; Brooks, B. R.; Schaefer, H. F., 3rd; Sokolov, A. Y.; Patkowski, K.; DePrince, A. E., 3rd; Bozkaya, U.; King, R. A.; Evangelista, F. A.; Turney, J. M.; Crawford, T. D.; Sherrill, C. D. Psi4 1.4: Open-source software for high-throughput quantum chemistry. _J. Chem. Phys._ 2020, _152_ , 184108
  * Kendall et al. 1992 Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the first‐row atoms revisited. Systematic basis sets and wave functions. _J. Chem. Phys._ 1992, _96_ , 6796–6806
  * Kosloff and Tal-Ezer 1986 Kosloff, R.; Tal-Ezer, H. A direct relaxation method for calculating eigenfunctions and eigenvalues of the schrödinger equation on a grid. _Chem. Phys. Lett._ 1986, _127_ , 223–230
  * Smyth et al. 1998 Smyth, E. S.; Parker, J. S.; Taylor, K. T. Numerical integration of the time-dependent Schrödinger equation for laser-driven helium. _Comput. Phys. Commun._ 1998, _114_ , 1–14
  * Kowalewski and Seeber 2022 Kowalewski, M.; Seeber, P. Sustainable packaging of quantum chemistry software with the Nix package manager. _Int. J. Quant. Chem._ 2022, _122_ , e26872
  * Hopfield 1958 Hopfield, J. J. Theory of the Contribution of Excitons to the Complex Dielectric Constant of Crystals. _Phys. Rev._ 1958, _112_ , 1555–1567
  * Pavosevic et al. 2023 Pavosevic, F.; Smith, R. L.; Rubio, A. Cavity Click Chemistry: Cavity-Catalyzed Azide-Alkyne Cycloaddition. 2023,
  * Pavošević et al. 2023 Pavošević, F.; Smith, R. L.; Rubio, A. Computational study on the catalytic control of endo/exo Diels-Alder reactions by cavity quantum vacuum fluctuations. _Nat. Commun._ 2023, _14_ , 2766
  * Bloino and Biczysko 2015 Bloino, J.; Biczysko, M. _Reference Module in Chemistry, Molecular Sciences and Chemical Engineering_ ; Elsevier, 2015
  * Raj et al. 2022 Raj, A.; Chao, Y.-B.; Witek, H. A. Testing the limitations of harmonic approximation in the determination of Raman intensities. _Mol. Phys._ 2022, _120_ , e2069613
  * Jaynes and Cummings 1963 Jaynes, E. T.; Cummings, F. W. Comparison of quantum and semiclassical radiation theories with application to the beam maser. _Proc. IEEE_ 1963, _51_ , 89–109
  * Kowalewski et al. 2016 Kowalewski, M.; Bennett, K.; Mukamel, S. Cavity Femtochemistry: Manipulating Nonadiabatic Dynamics at Avoided Crossings. _J. Phys. Chem. Lett._ 2016, _7_ , 2050–2054
  * Gudem and Kowalewski 2021 Gudem, M.; Kowalewski, M. Controlling the Photostability of Pyrrole with Optical Nanocavities. _J. Phys. Chem. A_ 2021, _125_ , 1142–1151
  * Couto and Kowalewski 2022 Couto, R. C.; Kowalewski, M. Suppressing non-radiative decay of photochromic organic molecular systems in the strong coupling regime. _Phys. Chem. Chem. Phys._ 2022,
  * Haugland et al. 2023 Haugland, T. S.; Philbin, J. P.; Ghosh, T. K.; Chen, M.; Koch, H.; Narang, P. Understanding the polaritonic ground state in cavity quantum electrodynamics. 2023,
  * Schnappinger and Kowalewski 2023 Schnappinger, T.; Kowalewski, M. Nonadiabatic Wave Packet Dynamics with Ab Initio Cavity-Born-Oppenheimer Potential Energy Surfaces. _J. Chem. Theory Comput._ 2023, _19_ , 460–471