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

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, 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 for a comprehensive analysis of these hybrid states. Our semiclassical approach, validated against full quantum simulations, reproduces key features of the vibro-polaritonic spectra. The underlying analytic gradients also allow for optimization of cavity-coupled molecular systems and performing semiclassical dynamics simulations.


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 Ω 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 ω 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. [15][16][17][18][19][23][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 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 H M give the vibro-polaritonic normal modes of the correlated matterphoton 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), 26,[31][32][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 H M via finite differences and subsequently abinitio 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 vibropolaritonic 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, [17][18][19]34,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 (NH 3 ). 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.

Theory
To describe the interaction of a confined light mode with atoms, molecules, and condensed matter systems, the non-relativistic Pauli-Fierz Hamiltonian 22,[35][36][37][38] 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 26,31-33 is a well suited theoretical approach. Within the CBOA, the cavity modes are grouped with the nuclei in a generalized Born-Huang expansion, 35,39 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 (ℏ = 4π; ε 0 = m e = 1) are used throughout the paper, unless otherwise noted. For a single mode cavity the electronic CBOA Hamiltonian takes the form whereμ represents the molecular dipole operator, which is defined by the operators of the N el electron coordinatesr, the classic coordinates R of the N N uc nuclei and the nuclear charges Z A .Ĥ 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 pho-ton displacement coordinate q c and ω 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 λ c . The last term is the DSE operator, 17,34,40 which is an energy contribution that describes the self-polarization of the molecule-cavity system. The coupling parameter λ c for a cavity with an effective mode volume V c is defined as follows: The unit vector e denotes the polarization axis of the cavity mode.

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 E CBO consists of four energy contributions: 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 E CBO in the basis of atomic orbitals results in where α, β, γ, δ denote atomic orbitals, F α,β the CBO-HF Fock matrix and D α,β 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Ẽ n,c : The modified one-electron integrals α h β and two-electron integrals αβ g γδ used to build the CBO-HF Fock matrix elements F α,β consist of the standard one-electron integral α h β and two-electron integrals αβ g γδ as well as terms describing the linear cavityelectron interaction and the electronic DSE contributions:

CBO-Hartree-Fock Gradients
The first derivative of the energy E CBO with respect to a nuclear or photon displacement 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 ϵ i , orbital coefficients c i and overlap integral derivatives: In the case of ζ i being a photon displacement coordinate, Eq. (9) becomes zero, since the atomic orbitals α, β are independent of q c . A detailed derivation of all "new" terms in Eq. (8) for ζ I = R i and ζ I = q c can be found in section S1 of the supporting informations.
All derivatives for ĥ , ĝ , 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 g CBO is a (3N A + 1) vector with N A being the number of atoms in the molecule:

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, H, with dimensions (3N A + 1) x (3N A + 1) for the case of a single photon mode. Each element H ij is the second derivative of E CBO with respect to nuclear coordinates respectively the photon displacement coordinate 30 that is Using the analytic CBO-HF gradient g CBO the second derivatives can be determined using finite differences: The harmonic frequencies and corresponding normal modes can be obtained by transforming the Hessian matrix H from the ζ i coordinates into their mass-weighted version H M and subsequently solving the eigenvalue problem: Here L is the diagonal matrix of eigenvalues ℓ i and A is the matrix that diagonalizes the mass-weighted Hessian H M . A is formed by the juxtaposition of its eigenvectors a i which define the normal modes. From the (3N A + 1) eigenvalues, only (3N A + 1 − 6), or 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 cm −1 , are given by: The corresponding intensities H I i in the harmonic approximation are calculated as the projection of the dipole moment gradient on the normal mode vectors a i : The necessary dipole moment gradient ∇ μ 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 ν 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 I j are calculated with the nuclear-photonic eigenfunctions χ i of the coupled molecular-cavity system:

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 Here λ 0 is equivalent to λ 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 λ 0 in the range of 0.004 to 0.04, which corresponds to effective mode volumes, see Eq.  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 λ 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 λ c leads to a linear increase of the Rabi splitting Ω 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 ∆Ω R = ω c − 0.5 ν LP + ν U P in the Rabi splitting ( Fig. 1 c)), where the LP is stronger red-shifted than the UP with respect to ω c . This asymmetry is observed for the harmonic approximation and the full-anharmonic treatment and in both cases ∆Ω R increases quadratically with λ c . However, the asymmetry is more pronounced in the full-anharmonic case. This observed asymmetry of Ω R can also be 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 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. 34,40 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 Ω 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 ω c for a fixed coupling strength.
For better visualization, the difference between Ω R and ∆ = ω c − H ν 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 ω c . Note that Ω R − ∆ is only a measure of the total size of Ω R , not of the asymmetry ∆Ω R , which is discussed in Fig. 1. As expected, the largest Rabi splitting Ω R is obtained for ω c being resonant with the fundamental bare molecular transition ( H ν 1 ). Interestingly, the difference between Ω R and the detuning ∆ is not symmetric with respect to Hν 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). 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 (χ LP ) and the UP state (χ U P ).
The absolute square of obtained expansion coefficients |a LP c | 2 and |a U P c | 2 can be used as a 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 Ω R and detuning ∆ (Fig. 3 a)) are close to resonance nearly identical and for large detunings qualitatively similar. The asymmetry ∆Ω R , as well as the asymmetry of the difference between Ω R and ∆ 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 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 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.

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 ap-proximation 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 α,β , and orbital energies ϵ 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 Ω R are shown in Fig. 4.
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 Ω 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 Ω R is reproduced quite well. However, its asymmetry ∆Ω R is smaller and characterized by blue-shift of LP and UP with respect to ω 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-  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  Comparing the all-parallel spectra with and without DSE, the differences are rather small. The observed Rabi splitting Ω 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 ∆Ω R is also independent of N M ol and constant (Fig. 6 c) blue line). The observed changes are more drastic in the case

Beyond Diatomic Molecules: Vibro-Polaritonic Spectra of NH 3
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 NH 3 molecule. The bare NH 3 molecule has six vibrational degrees of freedom: One symmetric bending mode (ν 1 ), two degenerate asymmetric bending modes (ν 2,3 ), one symmetric stretching mode (ν 4 ) and two degenerate asymmetric stretching modes (ν 5,6 ).
Since NH 3 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 NH 3 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 ω c is set to be resonant with the symmetric bending mode ν 1 , which has the lowest energy of 1103 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 NH 3 molecule are shown in Fig. 7 for the bending modes and in Fig. S14 of the supporting information for the stretching modes.
As expected, the splitting into a LP and a UP transition is observed for the symmetric bending mode ν 1 when coupled to a resonant optical cavity, see Fig. 7

b) and d). The Rabi
splitting Ω 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 (ν 2 and ν 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 λ c larger than 0.078 au the degeneracy of the two transitions is lifted. Figure 7: a) Low energy part of the vibronic IR spectra of a single NH 3 molecule calculated in the harmonic approximation. The low energy part of the vibro-polaritonic IR spectra of a single NH 3 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 ω c is resonant with the symmetric bending mode (1103 cm −1 ) and the cavity field strength λ c increases from 0.039 au to 0.155 au.
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 NH 3 , their magnitude is plotted for both polarization directions as a function of the coupling strength λ 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 λ 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. 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 ω c is resonant with the symmetric bending mode (1103 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 Ω R is larger and scales almost linear with ϵ 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 ν 2 and ν 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

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 vibropolaritonic IR spectra obtained for small ensembles of HF molecules with full quantummechanical 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 NH 3 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/fieldfree 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.  References 24 S1 Derivation of the CBO-Hartree-Fock gradients

Supporting Information Available
The energy expectation value E CBO in the basis of the atomic orbitals has the following form: dse . (S1) Wherex and X are the electronic and nuclear part of the projected molecular dipole moment: The first derivative of the energy E CBO with respect to a nuclear or photon displacement The first term in Eq. (S3) 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. The full Hellmann-Feynman term has the following structure: Since the cavity Born-Oppenheimer Hartree-Fock (CBO-HF) wave function is variational optimized, the explicit calculation of the density derivatives can be avoided 1 and the Pulay term can be written in terms of overlap integral derivatives: The remaining derivative is simpler since it dose not involve electron coordinates: In the following two cases will be discussed, the first where ζ i is a nuclear coordinate R i and the second where ζ i is a photon displacement coordinate q c . For the first case, the nuclear derivatives of ĥ , ĝ (in Eq. (S4)), V nn (in Eq. (S6)) and the overlap integral derivatives in (in Eq. (S5)) are identical to the standard Hartree-Fock gradient terms and can be found in the literature. 1-4 Therefore, we focus on the "new" terms introduced by the CBO-HF ansatz in Eq. (S4) and Eq. (S6). The nuclear derivative of the linear cavity-electron interaction term takes the following form: Here the central or Hellmann-Feynman term becomes equal zero sincex is independent of R i .
The nuclear derivatives of the pure electronic contributions to the dipole self-energy (DSE) follow a similar scheme. Whereby the one-electron contributions (Eq. (S8)) are simpler than the two electron contributions (Eq. (S9)).
The nuclear derivative for the mixed electron-nuclear DSE is Regarding Eq. (S6) there are only two "new" nuclear derivatives: For the second case ζ i being a photon displacement coordinate q c the Pulay term in Eq. (S5) is zero, since the atomic orbitals used are independent of q c . This is also the reason why for the Hellmann-Feynman part of the derivative (Eq. (S4)) only terms that explicitly depend on q c are of relevance. In combination with relevant parts of Eq. (S6) the following expressions for the derivative with respect to q c is found: This equivalent to the result obtained using the Hellmann-Feynman theorem. 5,6

S2 Validation of the Harmonic Approximation
The spectra of an individual HF molecule without coupling to an optical cavity calculated in the harmonic approximation and full quantum mechanical (anharmonic) are shown in Fig. S1.
The fundamental vibrational transition in the harmonic approximation has a frequency of to an optical cavity, 6 the vibro-polaritonic IR spectra are nearly indistinguishable for the two configurations studied. The only noticeable difference is the asymmetry of the Rabi splitting ∆Ω R = ω c − 0.5 ν LP + ν U P shown in Fig. S2 c) and Fig. S3 c). The value of ∆Ω R is slightly smaller for the antiparallel configuration. Fig. S4 shows the difference between Ω R and ∆ = ω c − H ν 1 in the harmonic approximation as well as the cavity-induced energy changes E lin , E dse and E dis as a function of the cavity frequency ω c , keeping λ c fixed for both configurations. Also for the bimolecular case, the largest Rabi splitting Ω R is obtained for ω c resonant with the fundamental molecular transition ( H ν 1 ). In agreement with the single molecule results, the difference between Ω R and the detuning ∆ is not symmetric with respect to H ν 1 and tends to a finite nonzero value even for large detunings, see Fig. S4 a).
As for the resonant spectra ( Fig. S2 and Fig. S3), the effect of cavity detuning is the same for both configurations, see gray line Fig. S4 a). The cavity-induced energy modifications E lin , E dse and E dis shown in Fig. S4 b) for the all-parallel configuration and in Fig. S4 c) for the antiparallel configuration are different for each configuration, but are constant for all values of ω c .    is the pure first excited molecular vibrational eigenstate with a nodal plan orthogonal to the R coordinate. The vibro-polaritonic normal modes give the classical pendants. The to-be LP normal mode (reduced normal mode shown in Fig. S5 d)) has only a contribution along the photon displacement coordinate q c , for visualization purposes shown on both atoms. The to-be UP normal mode (reduced normal mode shown in Fig. S5 f)) has only a contribution along the z coordinate axis in the laboratory frame, which describes the stretching mode of the HF molecule. Since we are discussing mass-weighted normal modes, the vector describing the stretching mode, shown in Fig. S5  Overall, both descriptions, the nuclear-photonic eigenfunctions and the vibro-polaritonic normal modes, contain similar information about the formation of hybrid-light matter states. Figure S5: a), c), and d) first three eigenfunctions for a single HF molecule described on two-dimensional CBO-HF surface. b) classic nuclear configuration, d), and e) both massweighted normal modes of a single HF molecule in harmonic approximation. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to zero. Figure S6: a), c), and d) first three eigenfunctions for a single HF molecule described on two-dimensional CBO-HF surface. b) classic nuclear configuration, d), and e) both both mass-weighted normal modes of a single HF molecule in harmonic approximation. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to 0.019 au. Figure S7: First excited state eigenfunction for two parallel HF molecules described on a three-dimensional CBO-HF surface. For a), b), and c, the eigenfunction is integrated over one coordinate each. d) classic nuclear configuration and first mass-weighted normal mode of two parallel HF molecules. The molecules are separated by a distance of 800Å. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to zero. Figure S8: Second excited state eigenfunction for two parallel HF molecules described on a three-dimensional CBO-HF surface. For a), b), and c, the eigenfunction is integrated over one coordinate each. d) classic nuclear configuration and second mass-weighted normal mode of two parallel HF molecules. The molecules are separated by a distance of 800Å. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to zero. Figure S9: Third excited state eigenfunction for two parallel HF molecules described on a three-dimensional CBO-HF surface. For a), b), and c, the eigenfunction is integrated over one coordinate each. d) classic nuclear configuration and third mass-weighted normal mode of two parallel HF molecules. The molecules are separated by a distance of 800Å. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to zero.
Figure S10: First excited state eigenfunction for two parallel HF molecules described on a three-dimensional CBO-HF surface. For a), b), and c, the eigenfunction is integrated over one coordinate each. d) classic nuclear configuration and first mass-weighted normal mode of two parallel HF molecules. The molecules are separated by a distance of 800Å. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to 0.019 au.
Figure S11: Second excited state eigenfunction for two parallel HF molecules described on a three-dimensional CBO-HF surface. For a), b), and c, the eigenfunction is integrated over one coordinate each. d) classic nuclear configuration and second mass-weighted normal mode of two parallel HF molecules. The molecules are separated by a distance of 800Å. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to 0.019 au.
Figure S12: Third excited state eigenfunction for two parallel HF molecules described on a three-dimensional CBO-HF surface. For a), b), and c, the eigenfunction is integrated over one coordinate each. d) classic nuclear configuration and third mass-weighted normal mode of two parallel HF molecules. The molecules are separated by a distance of 800Å. The cavity frequency ω c is resonant with the corresponding fundamental transition (harmonic 4467 cm −1 and anharmonic 4281 cm −1 ). The coupling strength λ c is set to 0.019 au.

S4 Beyond Diatomic Molecules: Vibro-Polaritonic Spectra of NH 3
The bare molecular vibronic IR spectrum and the vibro-polaritonic IR spectra of a single NH 3 molecule are shown in Fig. S13 for the three stretching modes. The cavity frequency ω c is set to be resonant with the symmetric bending mode ν 1 (1103 cm −1 ). The vibro-polaritonic IR spectra are calculated for the y polarization (Fig. S13 b) and d)) and the z polarization ( Fig. S13 c) and d)) of the single cavity mode. Since the chosen cavity frequency w c is off-resonant with the stretching modes, the cavity interaction leads to only modifications of the signals by changing the energetics and breaking the molecular symmetries. These effects are stronger for the y polarization because of the molecular orientation chosen. Figure S13: a) High energy part of the vibronic IR spectra of a single NH 3 molecule calculated in the harmonic approximation. High energy part of the vibro-polaritonic IR spectra of a single NH 3 molecule zoomed into the symmetric stretching (b) and d)) and asymmetric stretching modes (c) and e)). The polarization axis of the cavity mode is the y axis for b) and c) and the z axis for d) and e). The cavity frequency ω c is resonant with the symmetric bending mode (1103 cm −1 ) and the cavity field strength λ c is increased from 0.039 au to 0.155 au.