# Liquid water under vibrational strong coupling: An extended cavity Born–Oppenheimer molecular dynamics study

Jessica Bowles, Jaime de la Fuente Diez, Damien Laage, Johannes Richardi, Rodolphe Vuilleumier, Riccardo Spezia

► To cite this version:

Jessica Bowles, Jaime de la Fuente Diez, Damien Laage, Johannes Richardi, Rodolphe Vuilleumier, et al. Liquid water under vibrational strong coupling: An extended cavity Born–Oppenheimer molecular dynamics study. The Journal of Chemical Physics, 2025, 163 (4), pp.044104. 10.1063/5.0274246. hal-05193437

HAL Id: hal-05193437

https://hal.science/hal-05193437v1

Submitted on 30 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Liquid water under vibrational strong coupling: an extended cavity
Born-Oppenheimer molecular dynamics study

Jessica Bowles,¹ Jaime De La Fuente Diez,² Damien Laage,² Johannes Richardi,¹ Rodolphe Vuilleumier,² and Riccardo Spezia¹

¹Laboratoire de Chimie Théorique, Sorbonne Université, Paris, France.

²Laboratory CPCV, Department of Chemistry, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, Paris, France

(*Electronic mail: riccardo.spezia@sorbonne-universite.fr.)

A computational study of liquid water when the system is coupled with a (model) Fabry-Perot cavity is reported. At this end, the Cavity Born-Oppenheimer Molecular Dynamics approach proposed recently (Li et al., Proc. Nat. Acad. Sci. USA, 2020, 117, 18324–18331) is employed and different properties of water under vibrational strong coupling (VSC) are investigated. Different cavity frequencies are considered, corresponding to different modes in the IR spectrum of liquid water: high frequency (corresponding to O–H stretching modes), medium frequency (corresponding to water molecule bending) and low frequencies (corresponding to librational modes). Simulations were done both using classical and quantum nuclear dynamics, this last via Ring Polymer Molecular Dynamics. Similar effects of the cavity are obtained in both cases. Namely, whereas the infrared spectrum is observed to be split for all cavity frequencies, no effects on structural properties are detected. In addition, transport and dynamical properties, including the diffusion coefficient, molecular reorientation and hydrogen bond (HB) jump exchange times, show no effect due to cavity coupling when an extended statistical analysis is performed.

## I. INTRODUCTION

In recent years, there has been a growing interest in vibrational strong coupling (VSC) and in particular on its impact on chemical reactivity.¹⁻³ VSC is obtained by placing a molecular system inside an optical cavity such as a Fabri-Pérot cavity that is resonant with a selected vibrational mode of the matter, thereby generating a hybrid light-matter state. When the exchange of (virtual) photons is faster than their dissipation, the strong coupling regime is achieved.

Early approaches focused on the coupling of optical cavity modes with electronic states.⁴ However, subsequent research demonstrated the feasibility of coupling molecular vibrations with cavity modes in the infrared (IR) range, a phenomenon referred to as vibrational strong coupling (VSC).⁵,⁶

When the frequency of the cavity ($\omega_c$) is at resonance with the vibrational frequency ($\omega_m$) of the molecular system then a Rabi splitting ($\hbar\Omega_R$) is obtained in the infrared (IR) spectrum. It has been shown that VSC can also influence chemical reactivity.¹,⁷,⁸ Specifically, it has been observed that kinetics are altered, as highlighted in different studies,⁹⁻¹⁹ and more recently, some equilibrium properties have been found to undergo modifications. These include changes in conductivity,²⁰,²¹ and the equilibrium constant in the formation of charge-transfer complexes.²² On the other hand, dielectric constant is not affected by VSC, as recently deeply discussed by Michon and Simpkins.²³

While the Rabi splitting in IR spectrum is now well understood,²⁴⁻²⁸ it is much less clear how these other properties are modified in VSC regime. For this reason, the last years have been characterized by an important theoretical activity in the field.²⁹⁻³⁶ The general theoretical description of matter under VSC is typically based on the Pauli–Fierz (PF) quantum electrodynamics (QED) Hamiltonian:³⁷

$$ \hat{H} = \hat{H}_M + \hat{H}_F \quad (1) $$

where $\hat{H}_M$ is the Hamiltonian of the molecular system and $\hat{H}_F$ is the interaction with the vacuum field, which can be written as:²⁹,³⁵,³⁸⁻⁴²

$$ \hat{H}_F = \sum_{k,\lambda} \frac{1}{2} \hat{p}_{k,\lambda}^2 + \frac{1}{2} \omega_{k,\lambda}^2 \left( \hat{q}_{k,\lambda} + \sum_{s=1}^{N} \frac{\hat{\mu}_s \cdot \hat{\xi}_\lambda}{\sqrt{V\epsilon_0}} \right)^2, \quad (2) $$

where $k$ denotes the (virtual) photon with given frequency, $\omega_k$, position, $\hat{q}_{k,\lambda}$, and momentum, $\hat{p}_{k,\lambda}$, operators, $\lambda$ the two polarization directions orthogonal to the cavity walls (so typically $x$ and $y$), $V$ the cavity volume, $\epsilon_0$ the vacuum permittivity, $\hat{\mu}_s$ the dipole operator for the molecular system, $\hat{\xi}$ the polarization direction of the (virtual) photons and $N$ the number of modes interacting with the cavity field.

Various theoretical studies using this Hamiltonian have been performed in recent years. A first kind of approach focused on adapting transition-state theory (TST) for the presence of the cavity and thus obtained rate constants with an effective friction, in analogy with Kramers and Grote-Hynes theories.⁴³,⁴⁴ They generally considered only one mode for the matter coupled to the cavity with effective coupling parameters. For example, Li et al.¹⁸ applied this theory to explain changes in the reaction rate under strong coupling. Also Lindoy et al.³¹ used an analytical rate which shows that chemical kinetics can be modified (enhanced or suppressed) due to a modified friction under VSC. A quantum description via TST of cavity-catalysed adiabatic chemical reaction was then carried out by Yang and Cao,³² providing an explanation for the maximum modification when the cavity frequency is near the vibrational resonance frequency and also for the collective and selective aspects.

In a similar vein, quantum dynamics were performed using different models, but generally considering only one matter
mode coupled with the cavity. Morse dissociation or double
well models were studied in this way, showing an effect on
the rate.<sup>45,46</sup> Also quasi-classical trajectories were employed
to reveal the modification of water dimer dissociation due to
VSC.<sup>47</sup> Notably, using the open quantum system formalism,
Reichmann and co-workers showed that it is possible to obtain
the highest enhancement when the cavity is in resonance with
the reactant mode.<sup>45</sup>

Another approach, developed in Rubio’s group, is based
on the application of the coupling to the full nuclear-electron
Hamiltonian.33 This method was used, for example, to ex-
amine whether the collective coupling implies local mod-
ifications of chemical properties38 and applied on a HD+
molecule.34 Furthermore, quantum dynamics have been em-
ployed to describe a model of a chemical reaction in the con-
densed phase within the cavity.40,41,45 However, these meth-
ods based on quantum dynamics only allow the description of
model systems and most of the time, the matter is described
with only one explicit mode.

To study how VSC modifies the properties of a condensed
phase system, an approach based on molecular dynamics sim-
ulations was developed by Subotnik and co-workers,³⁰ called
Cavity Molecular Dynamics (CavMD). In this approach, the
atoms of the system interact with classical force fields and
the cavity acts on the dipole of the molecules using a simple
point charge representation of the molecular charge distribu-
tion (as in classical molecular dynamics). This approach was
first used to study liquid water³⁰ and then extended to other
systems. For example, it is possible to study the nonlinear re-
sponse of carbon dioxide under VSC,⁴⁸ revealing that a strong
pulse can lead to an enhanced overall molecular nonlinear ab-
sorption that causes an ultrashort lifetime (0.2 ps) of the lower
hybrid light-matter state. Furthermore, it has been employed
to show that intermolecular vibrational energy transfer is fa-
cilitated under VSC.⁴⁹ The same authors also numerically il-
lustrate that IR pumping a liquid confined in an optical cav-
ity under vibrational strong coupling can selectively excite a
small concentration of solute molecules.⁵⁰ Notably, this ap-
proach can be extended to include phenomenologically the
cavity loss.⁵¹ The CavMD approach has also been applied,
using a QM/MM approach, to study Fe(CO)₅ in n-dodecane⁵²
predicting a preferred vibrational energy transfer to the equa-
torial CO vibrations rather than the axial CO vibrations when
exciting the lower polariton.

Calculations on liquid water done by Subotnik and co-
workers where done setting the cavity at resonance with the
O–H vibrational mode and they do not observe any relevant
effect on properties like water structure, diffusion coefficient
or dielectric constant. However water has different vibrational
modes and our main interest is to see if a cavity in reso-
nance with these modes can be at the origin of the modifi-
cation of some properties. At this end we have considered
a simple molecular dynamics simulation of liquid water us-
ing the approach proposed by Li et al.,³⁰ and we investigated
different properties as the cavity is tuned at different frequen-
cies. Notably, based on the works of Laage, Hynes and co-
workers, it is known that many dynamical processes in liq-
id water, including, e.g., translational diffusion and reorien-

tation, are governed by hydrogen-bond exchanges.<sup>53-57</sup> The
rate-limiting step of these hydrogen-bond exchanges has been
shown to be the combination of the initial hydrogen-bond
elongation and of the new hydrogen-bond contraction. This
motion is a low-frequency motion and in previous theoreti-
cal studies on VSC on liquid water, the cavity frequency was
set only on high frequency values.<sup>30,58</sup> Here, we have tuned
the cavity frequency in different spectral regions correspond-
ing not only to O–H stretching but also to the bending and
to the vibrational modes. Furthermore, Sun and Rubeiro<sup>59</sup>
have pointed out that coupling the cavity with modes at lower
frequency (~500-600 cm<sup>-1</sup>) has a higher impact on an S<sub>N2</sub>
reaction. We have thus analyzed if different properties are
modified under VSC. In particular, dynamical properties, re-
lated to hydrogen bond dynamics were investigated as well as
transport properties. We are motivated by experimental stud-
ies showing that diffusion coefficient of water can be modified
under VSC.<sup>20,21</sup> Since also low vibrational frequencies can be
related to water diffusion (e.g. through the aforementioned
H-bond exchange mechanism), we have extended the study of
liquid water under VSC using the CavMD approach to other
bands present in the water IR spectrum.

In the first part, we restricted ourselves to classical evolu-
tion of particles, in order to strengthen the importance of a
large statistical sampling to minimize the uncertainties on the
calculated properties. It was shown that using Path-Integral
based MD, the IR splitting changes slightly, compared to clas-
sical simulations, as well as other properties like the H-bond
jump time.<sup>42,56,58</sup> Thus in a final part ring-polymer molecular
dynamics (RPMD) was employed to study nuclear quantum
effects on the properties.<sup>60</sup>

The remainder of the article is as follows: first we will de-
scribe the basic principles of the methods and the simulations
performed. Then we will show that polaritonic states are ob-
tained for all the cavity frequencies, and then we will report
results on structural and thermodynamic properties then trans-
port and dynamical properties. A discussion on the statistical
sampling needed in such simulations is presented, followed by
the study of possible nuclear quantum effects using RPMD.
The last section provides conclusions and outlook for future
studies.

II. METHODS

A. Molecular Dynamics Simulations

We simulated a system composed by a box of 216 water
molecules using the flexible q-TIP4P/F force field.⁶¹ The box
size was set at 18.64 Å to reproduce the experimental water
density at 300 K, corresponding to what used by Li et al.³⁰
Standard periodic boundary conditions with an Ewald sum-
mation for the intermolecular Coulombic interactions were
used. In the classical nuclear dynamics simulations, for each
cavity frequency (and for the simulation without the cavity),
we first equilibrated the system in the canonical ensemble at
300 K, using a Langevin thermostat with a friction caracteris-
tic time τ = 100 fs, and a time step of 1 fs, using the velocity
Verlet scheme.⁶² We then selected 30 initial conditions to run NVE simulations of 100 ps with a propagation time step of 0.5 fs to obtain all the properties except the dielectric constant where 20 NVE trajectories of 1 ns were performed and the IR and VDOS where 4 trajectories of 20 ps were carried out. In the RPMD simulations, the equilibration trajectories employed a white noise Langevin thermostat on the beads and a global velocity rescaling thermostat on the centroid mode, with $\tau = 100$ fs.⁶³,⁶⁴ For the dielectric constant, we then extracted 50 initial conditions and ran T-RPMD trajectories of 300 ps each, in which only the internal modes are thermostatied,⁶³ using the positions and velocities extracted from the equilibrated trajectories as initial conditions. For the dynamical properties, 30 NVE trajectories of 100 ps, also started from the equilibrated trajectories, were exploited when we used 8 beads and just two trajectories of 100 ps for 32 beads. All RPMD simulations employed a time step of 0.25 fs.

The description of the cavity followed the methodology suggested by Li et al.³⁰ using the CavMD software,⁶⁵ a modified version of i-PI.⁶⁶ This software is coupled here with the LAMMPS package to use classical force fields.⁶⁷ CavMD models the system through two equations of motions, one corresponding to the atoms and one corresponding to the virtual photons of the cavity. The equations of motion used are obtained after reducing all of the operators in Eq. 2 to the corresponding classical observables:

$$M_j \ddot{\mathbf{R}}_j = \mathbf{F}_j^{(0)} - \sum_{k,\lambda} \left( \varepsilon_{k,\lambda} \tilde{q}_{k,\lambda} + \frac{\varepsilon_{k,\lambda}^2}{m_{k,\lambda}} \frac{\omega_{k,\lambda}^2}{\sum_{s=1}^N \vec{\mu}_s \cdot \vec{\xi}_\lambda} \right) \frac{\partial \vec{\mu}_j \cdot \vec{\xi}_\lambda}{\partial \mathbf{R}_j}, \quad (3)$$

$$m_{k,\lambda} \ddot{\tilde{q}}_{k,\lambda} = -m_{k,\lambda} \omega_{k,\lambda}^2 \tilde{q}_{k,\lambda} - \varepsilon_{k,\lambda} \sum_{s=1}^N \vec{\mu}_s \cdot \vec{\xi}_\lambda, \quad (4)$$

$$\text{with the coupling constant: } \varepsilon_{k,\lambda} = \sqrt{m_{k,\lambda} \omega_{k,\lambda}^2 / V \varepsilon_0}. \quad (5)$$

In Eq. 3 and 4, an auxiliary mass $m_{k,\lambda}$ is introduced to align with standard MD simulations where particle mass information is required. Consequently, the positions representing the cavity are modified as follows: $q_{k,\lambda} = \tilde{q}_{k,\lambda} \cdot \sqrt{m_{k,\lambda}}$.

The positions and velocities of the atoms obtained from the dynamics are used to determine the properties of interest. Various values of the cavity frequency are chosen that correspond to frequencies where a signal is present in the IR spectrum, such as 3550, 1650, 600 and 400 cm⁻¹. For each cavity frequency, a different value of the coupling constant is chosen, to have the same interaction strength with the cavity for each frequency value, $\varepsilon_{k,\lambda} \omega_{k,\lambda} = \text{cte}$. For the cavity frequency at 3550 cm⁻¹ the same coupling constant is chosen as in Li et al.³⁰: 4 · 10⁻⁴ a.u. then the other values are adapted linearly from this following Equation 5: 1.9 · 10⁻⁴, 6.8 · 10⁻⁵, 4.5 · 10⁻⁵ a.u. respectively.

As pointed out by Li et al.³⁰, for lower cavity modes overtones of the cavity should also be included. We have then run additional trajectories setting the cavity fundamental frequency at 600 cm⁻¹, and including additional cavity overtone modes at 1200, 1800, 2400, 3000 and 3600 cm⁻¹ using the

multi-cavity method described by Li et al.³⁰ and for two values of coupling strength of 2.0 · 10⁻⁴ and 6.8 · 10⁻⁵ a.u. .

We should point out that while the water model used does not reproduce accurately all the desired liquid water properties, it has the great advantage of being widely used and is one of the simplest flexible water models. It was used by Li et al.³⁰ in their first study of water under VSC and thus it is a natural choice for an initial exploration of the impact of VSC on dynamical properties and with different cavity frequencies. We emphasize that our focus will be not much on the absolute values of the different properties, but rather on how (and if) they evolve under VSC.

## B. Calculations of Properties

From MD simulations we directly obtained different properties such as the infrared spectrum, the vibrational density of states (VDOS), the radial distribution function (RDF), the dielectric constant, the diffusion coefficient, the H-bond reorientation time and the H-bond jump times. Here we list briefly how they are obtained, while in the original articles more details are reported.

1. The infrared spectrum is determined as the Fourier transform of the auto-correlation function of the total electric dipole moment of the molecular system and the VDOS from the auto-correlation function of the velocities. As in the model, and shown in the original paper by Li et al.,³⁰ the direction of the cavity lies along the z-axis, so the IR spectra along x- and y-axis are those who should show the formation of a polaritonic state.

2. The dielectric constant, $\varepsilon$, is calculated from the usual equation for tin foil periodic boundary conditions with Ewald sums (dielectric constant of surrounding is $\infty$):⁶⁸⁻⁷⁰

$$\varepsilon = 1 + 3y_D g_K \quad (6)$$

where the Kirkwood factor $g_K = (\langle M^2 \rangle - \langle M \rangle^2)/(N\mu^2)$ and the dipole parameter $y_D = \rho\mu^2/9k_B T \varepsilon_0$. N is the number of molecules, $\mu$ is the absolute value of the molecular dipole averaged for each molecule and $M = \sum_i \vec{\mu}_i$, $\rho$ is the number density and $\varepsilon_0$ is the vacuum permittivity. Due to the well known slow convergence of the dielectric constant, we calculated that on longer trajectories (20 ns).

3. The diffusion coefficient⁷¹ is obtained through the slope of the mean square displacement (MSD) of the oxygen atoms. Here we limited ourselves to the diffusion coefficient corresponding to a given water box length, since we are interested into the impact of VSC keeping the same simulation conditions. In fact, it is well known that hydrodynamics suffers finite size effects and to recover the diffusion coefficient one should extrapolate the infinite box value.⁷²
4. To determine the reorientation time, the angle between the orientation of an OH bond at time $t_0$ and $t_1$ is examined. The order 2 Legendre polynomial of the cross correlation between these two vectors $P_2[v(t_0) \cdot v(t_1)]$ is determined. Then the natural logarithm of the Legendre polynomial is fitted over a time interval of 10 ps to obtain the slope as $ln(P_2) = -t/\tau_{reor}$.

5. To determine the jump time, the cross-correlation function of the probability for a water OH group to donate a stable hydrogen-bond to an initial acceptor I at time 0, $n_I(0)$, and to form a hydrogen-bond to a new stable partner F at time t, $n_F(t)$, is used:⁵⁴

$$C_{I,F}(t) = \langle n_I(0)n_F(t) \rangle. \qquad (7)$$

Here, $n_{I,F}(t) = 1$ when the OH group under consideration donates a stable H-bond to I, F at time t and $n_{I,F}(t) = 0$ otherwise. The H-bond presence is determined from strict geometric criteria ($R_{OO} < 2.95$ Å, $R_{HO} < 2.05$ Å, and $\theta_{HOO} < 20^{\circ}$)⁵⁴. This correlation function gives the probability to form a stable final H-bonded state at time t when the system was engaged in a stable H-bond with a different acceptor in the initial state at time 0. Absorbing boundary conditions in the product state are applied to consider the first acceptor exchange only. Thus the decay time of the complementary probability, $1-C_{I,F}(t) = \exp(-t/\tau_{jump})$, provides the jump time $\tau_{jump}$ corresponding to the time for a water molecule to switch stable H-bond acceptors.⁵⁴

6. Finally, it was shown that the free energy barrier corresponding to an hydrogen bond exchange can be estimated from the oxygen-oxygen radial distribution function.⁵⁶ This free energy barrier is the sum of two terms, which correspond respectively to the elongation of the distance between the oxygen atom of the jumping H-bond donor water and the initial acceptor molecule, and to the compression of the distance between the oxygen atom of the same H-bond donor molecule and the new acceptor molecule. Each term can be determined from the potential of mean force along the oxygen-oxygen distance, $W(r)$, which is related to the oxygen-oxygen radial distribution function (RDF), $g(r)$:

$$W(r) = -k_B T \ln[g(r)] \qquad (8)$$

where $g(r)$ corresponds to the O-O RDF. We denote $r_{max1}$, $r_{max2}$ and $r_{min}$ the distances corresponding to the first two maxima and the first minimum respectively of the RDF.

In the reactant state the average oxygen-oxygen distance corresponds to the first maximum of the $g(r)$ ($g(r_{max1})$, a minimum in the $W(r)$) and, before the jump, the oxygen atom of the final state will be at an average separation corresponding to the second maximum in the $g(r)$ ($g(r_{max2})$, a second minimum in the $W(r)$). It is then possible to estimate the free energy barrier of

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Scheme of vibrational strong coupling (b) IR spectra with cavity frequencies of 400, 600, 1650, 3550 cm⁻¹ and excluding the cavity (c) VDOS spectra with the same cavity frequencies. The black vertical lines represent the cavity frequency and the positions and velocities were collected over 80 ps.</figcaption></figure>

the H-bond exchange transition state from $g(r)$ calculated from simulations as:

$$\Delta G^{\ddagger} = -k_B T \left[ \ln \left( \frac{g(r_{min})}{g(r_{max1})} \right) + \ln \left( \frac{g(r_{min})}{g(r_{max2})} \right) \right]. \qquad (9)$$

Thus, we estimated the free energy barrier for the hydrogen-bond switch reaction simply from the RDF calculated from the simulations as a function of the cavity frequency.

In all the article, the error bars are reported as a confidence interval of 95%.

### III. RESULTS AND DISCUSSION

#### A. Polariton Formation

As previously discussed, the CavMD approach is able to mimic the formation of a hybride state in which the vibrational modes of liquid water are mixed with the resonant cavity mode (see a pictorial explanation in Figure 1a). The IR spectrum of liquid water, obtained from 4 NVE trajectories with the q-TIP4P/F model was calculated from the Fourier
transform of the dipole-dipole correlation function on the (x,y) plane. In the absence of cavity the IR spectrum exhibits the well-known bands (as shown in the bottom of Figure 1b) corresponding to the librations (around 600 cm⁻¹), the bending mode (at 1650 cm⁻¹) and OH stretching (at 3550 cm⁻¹). We chose to turn the cavity to each of these modes successively. Since the librational band is quite broad, we have also set the cavity to a lower value, 400 cm⁻¹. Results are summarized in Figure 1b.

When the cavity frequency is at resonance with a water vibrational mode, a splitting can be seen in the resulting IR spectrum, corresponding to the formation of two polaritonic modes. The two new peaks correspond to the P⁺ and P⁻ polaritonic modes that arise from the hybrid state as shown in the IR spectra reported in Figure 1(b). On the other hand, for the VDOS spectra, shown in Figure 1(c), no splitting occurs, since no selection rules apply and so also the modes which are not hybridized with matter, the dark states (DS), are present. These states are numerically dominant over the bright states as depicted in Figure 1(a) and explained numerically in details by de la Fuente et al.⁷³

In Figure 2 (panels A and B), we show in details the IR spectra obtained under resonance conditions: in particular we show both the spectrum calculating from dipoles in the (x,y) plane, corresponding to the plane on which the field is coupled to the matter, and in the z-axis. In panel A, we show the spectra obtained when the cavity is set to 600 cm⁻¹, while in panel B when it is set to 1650 cm⁻¹. Interestingly, we notice that for the signals in the light polarization plane, where the polaritons are formed, the other bands, which are out-of-resonance with the cavity, are also shifted. When the cavity is in resonance with the librational motion, we observe that the bending is blue shifted by about 80 cm⁻¹ (the maximum being now at about 1730 cm⁻¹) and the stretching is also blue shifted by about 100 cm⁻¹ (the maximum is now at about 3650 cm⁻¹). The same phenomenon is observed when the cavity is resonant with the bending: now the broad librational broad band is red-shifted, while the stretching is still blue shifted. This phenomenon is probably due to the fact that the cavity is off-resonance with the other modes, while it remains close enough to slightly shift them. These spectral features, namely the asymmetry of the polaritons and the shifts of the other bands, are obtained also by employing a simple harmonic oscillator model proposed by Lieberherr et al.⁴² in which the cavity-free IR spectrum was given as unperturbed input. Results are shown in the same Figure 2, panels A and B.

To understand these particular features of polaritonic bands, we have also employed a model Pauli-Fierz Hamiltonian where the matter is represented by three modes corresponding to the frequencies of liquid water. To represent the broadening of liquid phase, a distribution of frequencies is considered for each mode, reflecting the distribution width obtained from cavity-free simulations. The cavity is then added at a frequency of 1650 cm⁻¹. With this Hamiltonian, Path Integral Monte Carlo simulations allow to obtain IR spectrum and polaritons in very good agreement with exact quantum Hamiltonian, as shown recently by some of us.⁷³ More details

are reported in Appendix A.

In panel C of Figure 2, we report results obtained with this 3-mode model with the cavity set at 1650 cm⁻¹ (note that the non exactly unimodal IR bands of the uncoupled system is due to the numerical sampling of 100 vibrational modes). Notably, we observe that the high frequency peak is blue shifted while the low frequency is red shifted. This confirms that the Born-Oppenheimer cavity-MD approach is able to correctly reproduce spectroscopic properties for systems under VSC, even concerning detailed features such as the shift of modes not directly at resonance with the cavity.

We also obtain that the lower polariton (LP) is much less intense than the upper polariton (UP), as in molecular simulations. From PIMC simulations, we can obtain the coefficients of the cavity and of the matter modes in the coupled system which reflect the light-matter hybridization and can be used to get deeper details. In particular, we have analyzed how the cavity coefficients are distributed as a function of frequency. In the simple case of a "isolated" polariton, the cavity contribution is 50% on the LP and 50% on the UP.⁷³ Here, however, we observe that the cavity has a non negligible contribution also away from the polaritonic band. Notably we divided the spectrum in four regions: (I) $\omega < 1300$ cm⁻¹; (II) $1300 < \omega < 1650$ cm⁻¹; (III) $1650 < \omega < 2500$ cm⁻¹ and (IV) $\omega > 2500$ cm⁻¹. Regions II and III correspond to LP and UP. From the cavity coefficients obtained from PIMC simulations we got that while the UP (region III) has almost 50% of the cavity, the LP (region II) has about 35% of the cavity and region I has about 15%. This explains the decrease in intensity of the LP.

Furthermore, from PIMC simulations it is possible to evaluate the degree of delocalization of the different bands by calculating the Inverse Participation Ratio (IPR):

$$IPR_i = \sum_{j=1}^{N} |c_j^{(i)}|^4 \qquad (10)$$

where $c_j^{(i)}$ are the coefficients for each matter mode j of the final band i and N are the number of matter modes. In our previous study where only one vibrational band was considered, we obtained that the polaritonic bands have very low values, corresponding to high delocalization (note that by construction the maximum delocalization is 1/N, so 0.01 in the present model). Here we obtain that the upper polariton has a $IPR_{UP} = 0.012$, while the lower polariton has $IPR_{LP} = 0.16$. This means that the upper polariton is a "normal" polariton (i.e. as obtained when studying a single isolated vibrational mode), such that its delocalization is close to the maximal one. On the other hand, the lower polariton is more localized, with a value similar to what we obtained for a dark state.⁷³ It is still bright (but much less intense than UP) because the cavity coefficient is not zero even if lower than for UP because it is spread over the LP and lower frequencies (the librational band in the present case).
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 2. IR spectra under VSC as obtained from signal in the (x,y) plane (black solid line) and z-axis (red line) for two cavity frequencies: 600 cm⁻¹ (panel A) and 1650 cm⁻¹ (panel B). In Panel B we added also the spectrum obtained from the harmonic oscillator model developed by Lieberherr et al.⁴² (blue line). In Panel C we report the results obtained from a simple 3-mode model with the cavity set at 1650 cm⁻¹: red line the IR without coupling with the cavity and in black the coupled spectrum. In all panels we report as vertical dashed lines the cavity frequency.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 3. (a) O-O Radial distribution function (RDF) and (b) potential of mean force (PMF). The values without the cavity and for various cavity frequencies are plotted but are indistinguishable since they are all superimposed.</figcaption></figure>

## B. Structural and Thermodynamic Properties

In the following section, various structural and thermodynamic properties such as the oxygen-oxygen radial distribution function (RDF), the potential of mean force (PMF), the free energy barrier corresponding to H-bond exchange $\Delta G^{\ddagger}$ and the dielectric constant are examined.

The water structure is typically probed by the oxygen-oxygen radial distribution function which is reported in Figure 3(a). Here we report also the RDFs obtained when the cavity is resonant with different modes obtaining almost no effect, similarly to what was reported by Li et al.³⁰ for the cavity set only at 3550 cm⁻¹. Note that in the figure the curves overlap almost perfectly, such that they cannot be distinguished.

From the RDF, it is possible to estimate the free energy barrier for H-bond exchange, as summarized previously and discussed in details by Wilkins et al.⁵⁶ The barriers as a function of cavity frequency are obtained from the potential of mean force reported in Figure 3(b). Also the PMFs are not affected by the presence of a resonant cavity. Unsurprisingly, since there is no large change when including the cavity on both RDF and PMF, there is also no effect on the $\Delta G^{\ddagger}$. Values are summarized in Table I. This means that the thermodynamics of H-bond exchange is not modified when including a cavity mode with the CavMD formalism.

We should also notice that our values are in agreement with the free energy barrier obtained previously by Wilkins et al. (1.22 kcal/mol)⁵⁶ with the same water model but with no cavity.

To continue the study of water within the cavity, we look at the dielectric constant, $\epsilon$, obtained as function of cavity frequency. In this case, the total simulation time is increased from 3 ns to 20 ns due the well-known slow convergence of this property. Results are summarized in Table I which also includes values reported in the literature using the same water model. Also in this case, we cannot see any significant effect on the cavity. The results found in the literature with the
same water model (without the cavity or with the cavity set at 3550 cm⁻¹) also correspond to our findings.

The fact that no effect is observed for structural and thermodynamic properties, while searched extensively, should be expected for classical particles under the Pauli-Fierz Hamiltonian, as a consequence of classical statistical mechanics. A concise derivation is given by Li et al.³⁰ and in the Appendix B we provide a more extended derivation.

### C. Transport and Dynamical Properties

We now dive into transport and dynamical properties of liquid water that can be extracted from molecular dynamics simulations, and examine the diffusion coefficient, the water reorientation time and the H-bond jump time. As discussed previously, to obtain transport and dynamical properties we averaged, for each value of cavity frequency ($\omega_c$) over 30 NVE trajectories generated from a long NVT trajectory. The cavity is tuned to the same frequencies as for the structural and thermodynamic properties.

We first consider the diffusion coefficients, shown in Figure 4(a) and Table I as a function of the cavity frequency. The value without the cavity, and the corresponding error bars, are shown as horizontal lines to help guide the eye when comparing to with the cavity. In panel (b) of the same Figure 4, the water reorientation time is represented and corresponding values are reported Table I. Finally, in panel (c), we report the jump time obtained from the same simulations (again the corresponding values are listed in Table I).

As one can see from the Figure and in details from the specific values, there is no measurable effect from the cavity for the three properties considered. It is important to notice (and we will discuss more this aspect later) that the error bars are crucial quantities to be carefully considered (and here we report the 95 % confidence interval). One could notice that for $\omega_c = 1650$ cm⁻¹ $\tau_{reor}$ and $\tau_{jump}$ are slightly slower than when the cavity is absent. Even if the average value is outside the error bars of the results obtained without the cavity, the difference is very small and the interval of confidence overlap. As we will show in next subsection, a more detailed statistical investigation confirms that the deviation is statistically irrelevant.

Furthermore, even when increasing the coupling strength, we do not see any variation of such properties. In Figure 5, we report results obtained when tuning the cavity at 400 and 3550 cm⁻¹ at different coupling strengths: variations are inside the uncertainty interval and they do not show any clear trend which could be a sign of some possible effect.

Finally, we have tested if the multi-cavity description could have an impact on such properties, with a fundamental frequency of the cavity at 600 cm⁻¹. Notably, for the 6.8 · 10⁵ a.u coupling strength we obtained $D = 0.195 \pm 0.030$ Å ps⁻¹, $\tau_{reor} = 3.26 \pm 0.23$ ps and $\tau_{jump} = 3.97 \pm 0.24$ ps, while for the coupling strength of 2.0 · 10⁴ a.u., we got $D = 0.209 \pm 0.018$ Å ps⁻¹, $\tau_{reor} = 3.04 \pm 0.24$ ps and $\tau_{jump} = 3.67 \pm 0.26$ ps. These values are, within the uncertainty, equivalent to what obtained without the cavity or with

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 4. (a) Diffusion coefficients, (b) OH reorientation times (c) H-bond jump characteristic times obtained from CavMD simulations as a function of cavity frequency. The error bars correspond to 95% of interval of confidence. The value outside cavity and corresponding error bars are reported as horizontal lines.</figcaption></figure>

a single cavity mode.

### D. Statistical Discussion

As pointed out previously, the evaluation of statistical uncertainty is a key aspect to decipher if a property can be
<table><caption>TABLE I. For different values of cavity frequency, the HB exchange free energy barrier ($\Delta G^{\ddagger}$), the dielectric constant, the diffusion coefficient ($D$), the reorientation time ($\tau_{reor}$) and the jump time ($\tau_{jump}$) are given as obtained from CavMD simulations with classical nuclear dynamics. The reported uncertainty is the 95% confidence interval.</caption><thead><tr><th>ω<sub>c</sub><br>(cm<sup>-1</sup>)</th><th>Coupling<br>(a.u.)</th><th>ΔG<sup>‡</sup><br>(kcal/mol)</th><th>ε</th><th>D<br>(10<sup>-3</sup> Å<sup>2</sup>·ps<sup>-1</sup>)</th><th>τ<sub>reor</sub><br>(ps)</th><th>τ<sub>jump</sub><br>(ps)</th></tr></thead><tbody><tr><td>No cavity</td><td>0.000</td><td>1.22 ± 0.03</td><td>56 ± 2<br>52.7 ± 6.5<sup>a</sup><br>54.5<sup>b</sup></td><td>195 ± 6</td><td>3.28 ± 0.08</td><td>3.95 ± 0.07</td></tr><tr><td>400</td><td>4.5·10<sup>-5</sup></td><td>1.21 ± 0.02</td><td>55 ± 2</td><td>196 ± 5</td><td>3.25 ± 0.07</td><td>3.92 ± 0.06</td></tr><tr><td>600</td><td>6.8·10<sup>-5</sup></td><td>1.21 ± 0.03</td><td>57 ± 2</td><td>197 ± 5</td><td>3.22 ± 0.07</td><td>3.93 ± 0.06</td></tr><tr><td>1650</td><td>1.9·10<sup>-4</sup></td><td>1.20 ± 0.02</td><td>55 ± 2</td><td>196 ± 5</td><td>3.19 ± 0.07</td><td>3.86 ± 0.05</td></tr><tr><td>3550</td><td>4.0·10<sup>-4</sup></td><td>1.20 ± 0.02</td><td>56 ± 2<br>54.5<sup>b</sup></td><td>194 ± 5</td><td>3.22 ± 0.05</td><td>3.88 ± 0.05</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>FIG. 5. Diffusion coefficient (D), water reorientation time ($\tau_{reor}$) and water jump time ($\tau_{jump}$) as a function of coupling strength with the cavity for two values of cavity frequency: 400 cm$^{-1}$ (black) and 3550 cm$^{-1}$ (red). Error bars correspond to 95% of interval of confidence.</figcaption></figure>

considered modified or not. It is also important to evaluate how many trajectories are necessary to have a reasonable convergence. At this end, we have employed the bootstrapping method, increasing the number of trajectories up to 60 in evaluating transport and dynamical properties (D, $\tau_{reor}$ and $\tau_{jump}$). Using bootstrapping, a dataset of size $m$ with $m \in [2, 60]$ is randomly selected from the extended dataset, and the mean of this sample of size $m$ is then calculated. For each $m$, this random selection is repeated 1000 times from which we can determine an error reported here as the 95% confidence interval. As discussed previously, one could argue that for $\omega_c = 1650$ cm$^{-1}$ there can be some effect of the cavity: we have thus considered this cavity frequency in the bootstrapping analysis.

In Figure 6, we report the results as a function of trajectories, where the mean and 95% confidence interval are reported

as a function of trajectories considered. Note that for each value a random re-sampling is performed, such that values are not exactly the same as in Table I, where only the (first) full set of 30 trajectories is considered. The results from simulations without the cavity are compared to simulations with a cavity frequency of 1650 cm$^{-1}$.

As shown in Figure 6, the cavity has no effect even increasing the number of trajectories to 60. Notably also the size of the error bars does not change dramatically from 30 to 60. This is important, since it will allow us to use 30 trajectories to fully converge results and to estimate a reasonable uncertainty when using less trajectories (which can be necessary when the computational cost increases).

The same approach is then considered for the dielectric constant, where now we consider a maximum of 20 trajectories of 1 ns each. The comparison between with and without the cavity (again $\omega_c = 1650$ cm$^{-1}$) is shown in Figure 7. Here we show that 15 simulations is sufficient for convergence and further simulations would not change the results dramatically.

Finally, the error bars can be estimated for different number of trajectories and used for simulations where very few trajectories are used, which is typically the case when they are computationally expensive as is the case of RPMD with 32 beads. In the present case, we could use the error bars estimated from the bootstrapping procedure to obtain an interval of confidence when performing only two trajectories.

## E. Nuclear Quantum Effects

As we have shown, Newtonian Born-Oppenheimer MD simulations of liquid water do not show any remarkable effect of the cavity on both equilibrium, transport and dynamical properties. As we have discussed (and demonstrated in the Appendix B), the equilibrium thermodynamics properties should not modified via the interaction with a cavity when using the Pauli-Fierz Hamiltonian. This is not the case for dynamical properties, but our study shows that, at least for the system and the properties investigated, there is no effect of the cavity neither for them. However, this is not necessarily true

<footnote>a From Eltareb et al.⁷⁴</footnote>

<footnote>b From Li et al.⁵⁸</footnote>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 6. Mean and 95% interval of confidence of the (a) Diffusion coefficients, (b) OH reorientation times (c) H-bond jump characteristic times obtained for a rising number of CavMD simulations using bootstrapping to randomly sample the data set of 60 trajectories without the cavity and for $\omega_c=1650 \text{ cm}^{-1}$.</figcaption></figure>

when considering nuclei as quantum particles. A recent study has pointed out that the inclusion of the cavity can introduce a modification of the barrier due to a difference in zero-point energy.⁵⁹ Cavity Born-Oppenheimer molecular dynamics can include nuclear quantum effects (NQE) through ring-polymer MD (RPMD).⁵⁸ These simulations are much more expensive than classical MD, and this can be an issue given the required statistical sampling (see previous subsection). We have thus considered RPMD with a limited number of beads ($P = 8$) to investigate if NQE are relevant and tested qualitatively on few simulations with a larger number of beads ($P = 32$).

We first compare our RPMD simulations results with $P = 8$

<figure><img src="image_15.png" /><figcaption>FIG. 7. Mean and 95% interval of confidence of the dielectric constant obtained for a rising number of CavMD simulations using bootstrapping to randomly sample the data set of 20 trajectories for $\omega_c=1650 \text{ cm}^{-1}$.</figcaption></figure>

with previous RPMD simulations in which a larger number of beads was used. Results are summarized in Table II. Notably, Eltareb et al.⁷⁴ report values of dielectric constant and diffusion coefficient which are consistent with our results. In particular we have a very good agreement for the dielectric constant. The diffusion coefficient found in our simulations is slightly larger than those reported by them, and increasing from classical dynamics obtained in the present study. However, simulations reported by Eltareb et al. have a different number of water molecules (512) and the diffusion coefficient was corrected for the box finite box site effect. While this difference can be attributed to the different number of beads and statistics, a similar behaviour (i.e. the increasing of diffusion coefficient from classical to RPMD simulations) was obtained by Habershon et al.⁶¹.

We have then turned the cavity on at the $600 \text{ cm}^{-1}$ frequency, which corresponds to librational motion also when adding NQEs. Also in this case we do not observe any effect of the cavity on both dielectric constant and diffusion coefficient as reported in Table II. We notice that Li et al. have also found no effect of the cavity on dielectric constant when tuning the cavity at the O-H stretching frequency.⁵⁸

Moving to the water reorientation and jump times, Wilkins et al. determined a small acceleration when moving from classical to RPMD simulations.⁵⁶ Our simulations show the same behaviour. When then adding the cavity with RPMD simulations, these two dynamical quantities are also unchanged as obtained in Newtonian simulations. Note that the cavity is tuned to $600 \text{ cm}^{-1}$ as this lower frequency of the librational modes contributes towards the mechanism responsible for H-bond jumps through the angular motion of the O-H bond.⁵³,⁵⁴

Results with $P = 8$ are reported by averaging 50 trajectories of 300 ps for the dielectric constant and from 30 trajectories of 100 ps for the other properties. According to our statistical discussion, these simulations times are enough to be able to give a conclusion on the results. Simulations with $P = 32$ were done only on two trajectories of 100 ps in order to calculate diffusion coefficient, reorientation and jump times. They are estimated at $0.216 \pm 0.015 \text{ Å}^2 \cdot \text{ps}^{-1}$, $2.7 \pm 0.3$ ps and
<table><caption>TABLE II. Results obtained from RPMD simulations with P = 8 without the cavity and with the cavity set at 600 cm⁻¹ with a coupling strength of 6.8 · 10⁻⁵ au. Uncertainties are given as 95% of confidence interval. We also report results obtained by Li et al.⁵⁸ when setting the cavity resonant with O–H stretching and using P = 32.</caption><thead><tr><th></th><th>No Cavity</th><th>ω<sub>c</sub> = 600 cm<sup>-1</sup></th><th>ω<sub>c</sub> = 3400 cm<sup>-1</sup></th></tr></thead><tbody><tr><td>ε</td><td>55 ± 2<br>53 ± 7<sup>a</sup></td><td>58 ± 2</td><td></td></tr><tr><td>D (10<sup>-3</sup> Å<sup>2</sup>·ps<sup>-1</sup>)</td><td>59<sup>c</sup><br>246 ± 6<br>220<sup>a</sup><br>221<sup>b</sup> ± 1</td><td>247 ± 7</td><td>56.6<sup>c</sup></td></tr><tr><td>τ<sub>reor</sub> (ps)</td><td>2.54 ± 0.07</td><td>2.51 ± 0.07</td><td></td></tr><tr><td>τ<sub>jump</sub> (ps)</td><td>3.13 ± 0.04</td><td>3.13 ± 0.05</td><td></td></tr></tbody></table>

3.5 ± 0.3 ps (the uncertainties are estimated from the boot- strap analysis done in the classical trajectories). These re- sults show that increasing the number of beads do not change the fact that the cavity, within this simulation method and for these properties, has no effect them.

## IV. CONCLUSIONS

In the present paper, we have studied equilibrium and dynamical properties of liquid water under vibrational strong coupling by using the Born-Oppenheimer cavity MD method proposed by Li et al.³⁰,⁵⁸ in which the q-TIP-4P/F model is used. At this end, we have set the cavity frequency to the different modes present in the water IR spectrum, corresponding to the O–H stretching, H–O–H bending and librational modes.

Similarly to what observed by Li et al.³⁰,⁵⁸ we do not observe any effect of the cavity on any property investigated, also monitoring other cavity frequencies in resonance with bending or librational motion (previous works studied only the cavity in resonance with O–H stretching). We extended the analysis to other properties, with respect to what present in the literature, and namely the diffusion coefficient and dynamical properties of water motion (reorientational and jump times). Also these properties are not modified by the presence of the cavity at any frequency in resonance with different liquid water IR bands. The inclusion of nuclear quantum effects via RPMD approach does not change the picture.

We should notice, that the experimental possibility of a modification of the dielectric constant due to VSC had a controversy history in the recent literature: while a 2023 experimental paper reported an effect on it,⁷⁵ it was retracted very recently⁷⁶ and an explanation was given by Michon and Simpkins.²³ Our simulations are in agreement with those results and our statistical analysis shows that one has to be very careful in statistical convergence and uncertainty evaluation also in simulations. This is true not only for dielectric constant (which needs long trajectories to converge as it is well-known)

but also for other properties, for which one can see illusory effects without a correct statistical analysis.

Concerning the diffusion coefficient, experiments have shown that the diffusion of ions can be modified by VSC and in case of water this is likely due to a difference in proton conductivity.²⁰,²¹ The non reactive q-TIP4P/F model clearly cannot provide any support on it and our study shows that with this popular water model, the diffusion coefficient is not altered by VSC. Unfortunately there are no experimental studies on the H-bond reorientational times under VSC. H-bond jump times are not accessible experimentally, but were studied here since they hold all the properties of a chemical reaction even if no covalent bond is broken/made, but only H-bonds. Also on those properties related to H-bond dynamics, using the present model, we cannot highlight any effect of VSC at any cavity frequency.

As pointed out, we have extensively studied liquid water under VSC employing the q-TIP4P/F water model. This model has two limitations which could be at the origin of the absence of any relevant effect of VSC on dynamical properties: (i) it is a fixed charge model (the cavity field does not change the charge distribution of water molecules); (ii) the water model does not allow bond breaking/making. Nevertheless, the model is able to recover the spectroscopic features, not only in terms of formation of polaritonic bands in the IR signal but also the small blue and red shifts of other bands not in resonance with the cavity mode. As pointed out by Lieberherr et al.⁴², the coupling strength used to reproduce the typical Rabi splitting using the q-TIP4P/F water model is unphysically high. This provides additional support to the conclusions of our work: since an extremely large coupling strength does not lead to any significant impact on any of the investigated water properties, one can safely suggest that lower coupling values should not cause any meaningful change in the liquid water properties. In addition, we note that while the qTIP4P-F model already provides an excellent description of the properties of water, future extensions of this work will examine VSC effects with more sophisticated water models.

The fixed cavity model could be improved in the future considering molecular polarizability, as a first step. In fact, recent theoretical studies by Rubio and co-workers have shown that VSC could modify the polarizability of molecules present in the cavity.⁷⁷ Unfortunately, their method is based on very expensive calculations which cannot be presently applied to liquid systems. Water models including polarizability can be employed as a first approximated description of fluctuating molecular electronic density (via induced dipoles, for example) and the Pauli-Fierz Hamiltonian can be adapted. Additionally, an interesting question is whether the inclusion of polarization, in classical MD, through models such as the fluctuating charges model, the Drude model or the induced dipoles model is sufficient or is the consideration of the electronic structure necessary. In any case, it will be important to both consider a number of molecules under VSC which is sufficiently large (at least 100)⁷³ and have a statistical convergence of the properties. Our research are now going in that direction.

<footnote>a From Eltareb et al.⁷⁴</footnote>

<footnote>b From Habershon et al.⁶¹</footnote>

<footnote>c From Li et al.⁵⁸</footnote>
# ACKNOWLEDGMENTS

The authors thank Dr. Tao Li for providing support in utilising the CavMD software. This work was supported by the French National Research Agency ANR (ANR-20-CE29-0008-02, MoMaVSC). We also thank IDRIS for generous allocation of computing time (project 103847).

# AUTHOR DECLARATIONS

## Conflict of Interest

The authors have no conflicts to disclose.

# DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

# Appendix A: 3-Mode Model Pauli-Fierz Hamiltonian and Spectra Calculations from Path Integral Monte Carlo

We have considered the Pauli-Fierz Hamiltonian, Equation 2, where now the matter potential is expressed via three harmonic oscillators centered at $\omega_{lib} = 650 \text{ cm}^{-1}$, $\omega_{bend} = 1650 \text{ cm}^{-1}$ and $\omega_{str} = 3500 \text{ cm}^{-1}$:

$$V_M = \frac{1}{2} m \sum_{i=1}^{N_{lib}} \sum_{j=1}^{N_{bend}} \sum_{k=1}^{N_{str}} (\omega_{lib,i}^2 q_{lib,i}^2 + \omega_{bend,j}^2 q_{bend,j}^2 + \omega_{str,k}^2 q_{str,k}^2) \quad (A1)$$

where for each mode we have sampled a Gaussian distribution such that the integral of the cavity-free signal obtained from simulations is obtained at best. Here we set $N_{lib} + N_{bend} + N_{str} = 100$. This finite distribution is at the origin of the small discrepancies on the cavity-free spectrum which in any case do not modify the main results. The coupling constant in the Pauli-Fierz Hamiltonian is empirically set to reproduce the Rabi splitting observed in the present CavMD simulations.

With this Hamiltonian, we set the cavity at $1650 \text{ cm}^{-1}$ and run Path Integral Monte Carlo (PIMC) simulations from which vibrational frequencies are obtained. Details are described in our recent work$^{73}$. Here we give a short overview.

Specifically, we run $10^8$ PIMC steps at 300 K using 64 beads and then we averaged over 50 replicas. From these simulations we calculated the displacement-displacement correlation function in the Path Integral formalism and solved the generalized eigenvalue problem following Morresi et al.$^{78}$ In this way we obtained the transition frequencies and the IR intensity is obtained from the matter components of the eigenvectors as detailed in Ref.$^{73}$

# Appendix B: Invariance of Structural and Thermodynamic Properties for Classical Particles under VSC

In this section, we will show how, as long as a classical canonical statistical mechanics are considered, the structural or thermodynamic properties of matter remain invariant under VSC, described by the Pauli-Fierz Hamiltonian in Eq.(2). A concise derivation is given by Li et al.$^{30}$ and here we provide an extended version of it.

The structural or thermodynamic value of an arbitrary property A, well defined at every point in phase space, can be computed by:

$$ \langle A \rangle = \int_{-\infty}^{\infty} dp_1 \dots \int_{-\infty}^{\infty} dq_1 \dots P(q_1, q_2, \dots, q_{3N}; p_1, p_2, \dots, p_{3N}) \\ A(q_1, q_2, \dots, q_{3N}; p_1, p_2, \dots, p_{3N}) \quad (B1) $$

Where $P(\{p_i, q_i\})$ is the probability of visiting the point in phase space given by the phase space vector $(q_1, q_2, \dots, q_{3N}; p_1, p_2, \dots, p_{3N})$. We will, for the sake of compactness, from now on abbreviate a phase space vector as $(\mathbf{q}, \mathbf{p})$. We will show that the phase space distribution function for matter is not affected by coupling to the cavity. To do so, consider that in the canonical ensemble, the solution to the Liouville equation that gives the phase space distribution function is:

$$ P(\mathbf{q}, \mathbf{p}) = \frac{e^{-\beta H}}{Z_T} \quad (B2) $$

Where $\beta = 1/k_B T$, and $Z_T$ is the partition function over all degrees of freedom, defined as $Z_T = \int_{-\infty}^{\infty} dp \int_{-\infty}^{\infty} dq e^{-\beta H(\mathbf{q}, \mathbf{p})}$. Let us now consider separately the degrees of freedom of the cavity, $\mathbf{q}_c$ and $\mathbf{p}_c$, which are the field coordinates and it conjugated momenta, respectively, for all the field modes and polarizations involved in the model. We also denote collectively $\mathbf{q}_M$ and $\mathbf{p}_M$ all the degrees of freedom relating to the matter part. If one is interested in an structural or thermodynamic property which is defined only on the matter part of the system, one must compute the marginal probability density function associated with the matter subsystem by integrating over all cavity degrees of freedom.

$$ P(\mathbf{q}_M, \mathbf{p}_M) = \frac{1}{Z_T} \int_{-\infty}^{\infty} d\mathbf{q}_c \int_{-\infty}^{\infty} d\mathbf{p}_c e^{-\beta H} \quad (B3) $$

In order to perform the integration, we note that the Hamiltonian in Equation (2) can be written, in terms of our collective notation for the modes, as:

$$ H = \sum_{k, \lambda} \left[ \frac{1}{2} p_{k, \lambda}^2 + \frac{1}{2} \omega_{k, \lambda}^2 (q_c + f(\mathbf{q}_M, \mathbf{p}_M))^2 \right] + H_M(\mathbf{q}_M, \mathbf{p}_M) \quad (B4) $$

Where we stress that the last term in the parenthesis is only a function of matter degrees of freedom, and $H_M$ is the Hamiltonian that describes the isolated matter subsystem. Now, going back to the integral in Eq(B3), as we restrict ourselves to
classical statistical mechanics, we can factorize the sum over modes and polarizations in a product of exponentials:

$$P(\mathbf{q}_M, \mathbf{p}_M) = \frac{e^{-\beta H_M(\mathbf{q}_M, \mathbf{p}_M)}}{Z_T} \int_{-\infty}^{\infty} dp_{1,1} e^{-\frac{\beta}{2} p_{1,1}^2} \\ \qquad \int_{-\infty}^{\infty} dp_{1,2} e^{-\frac{\beta}{2} p_{1,2}^2} \int_{-\infty}^{\infty} dp_{2,1} e^{-\frac{\beta}{2} p_{2,1}^2} \dots \\ \qquad \int_{-\infty}^{\infty} dq_{1,1} e^{-\frac{\beta}{2} \omega_{1,1} (q_{1,1} + f(\mathbf{q}_M, \mathbf{p}_M))^2} \\ \qquad \int_{-\infty}^{\infty} dq_{1,2} e^{-\frac{\beta}{2} \omega_{1,2} (q_{1,2} + f(\mathbf{q}_M, \mathbf{p}_M))^2} \\ \qquad \int_{-\infty}^{\infty} dq_{2,1} e^{-\frac{\beta}{2} \omega_{2,1} (q_{2,1} + f(\mathbf{q}_M, \mathbf{p}_M))^2} \dots \quad (B5)$$

We can here recognize that each integral involving a field momentum is a gaussian integral, that is solved by $\int_{-\infty}^{\infty} dx e^{-ax^2} = \sqrt{\frac{\pi}{a}}$. If we let M be the number of modes considered, and 2 polarizations are allowed for each mode:

$$P(\mathbf{q}_M, \mathbf{p}_M) = \left(\frac{2\pi}{\beta}\right)^M \frac{e^{-\beta H_M(\mathbf{q}_M, \mathbf{p}_M)}}{Z_T} \int dq_{1,1} e^{-\frac{\beta}{2} \omega_{1,1} (q_{1,1} + f(\mathbf{q}_M, \mathbf{p}_M))^2} \\ \qquad \int_{-\infty}^{\infty} dq_{1,2} e^{-\frac{\beta}{2} \omega_{1,2} (q_{1,2} + f(\mathbf{q}_M, \mathbf{p}_M))^2} \\ \qquad \int_{-\infty}^{\infty} dq_{2,1} e^{-\frac{\beta}{2} \omega_{2,1} (q_{2,1} + f(\mathbf{q}_M, \mathbf{p}_M))^2} \dots \quad (B6)$$

A second set of gaussian integrals can also be recognized, since the terms involving the exponential of the field coordinate involve a term dependent exclusively on the matter coordinates, which is a constant to the integral. A change of variable of the form $x = q + f(\mathbf{q}_M, \mathbf{p}_M)$ allows to integrate each cavity coordinate term:

$$P(\mathbf{q}_M, \mathbf{p}_M) = \left(\frac{2\pi}{\beta}\right)^{2M} \left(\frac{1}{\prod_{k,\lambda} \omega_{k,\lambda}^2}\right)^{1/2} \frac{e^{-\beta H_M(\mathbf{q}_M, \mathbf{p}_M)}}{Z_T} \\ = \frac{e^{-\beta H_M(\mathbf{q}_M, \mathbf{p}_M)}}{Z_M} \quad (B7)$$

Where $Z_M$ is the partition function of a system consisting exclusively on matter, $Z_M = \int d\mathbf{q}_M \int d\mathbf{p}_M e^{-\beta H_M(\mathbf{q}_M, \mathbf{p}_M)}$. All contributions to the cavity have been integrated, and all have yielded constant factors, that do not favor any particular matter configuration. The phase space distribution inside the matter manifold is controlled solely by $e^{-\beta H_M}$, as it would be in the case of a ordinary system without interaction to field variables. Therefore, this shows that within the Pauli-Fierz Hamiltonian description, VSC does not affect the structural or thermodynamic values of matter-only properties following classical statistical mechanics.

<footnote n="1">¹K. Nagarajan, A. Thomas, and T. W. Ebbesen, “Chemistry under vibrational strong coupling,” J. Am. Chem. Soc. **143**, 16877–16889 (2021).</footnote>

<footnote n="2">²K. Hirai, J. A. Hutchison, and H. Uji-i, “Molecular chemistry in cavity strong coupling,” Chem. Rev. **123**, 8099–8126 (2023).</footnote>

<footnote n="3">³R. Bhuyan, J. Mony, O. Kotov, G. W. Castellanos, J. Gómez Rivas, T. O. Shegai, and K. Börjesson, “The rise and current status of polaritonic photochemistry and photophysics,” Chem. Rev. **123**, 10877–10919 (2023).</footnote>

<footnote n="4">⁴J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, “Modifying chemical landscapes by coupling to vacuum fields,” Angew. Chem. Int. Ed. **51**, 1592–1596 (2012).</footnote>

<footnote n="5">⁵A. Shalabney, J. George, J. Hutchison, G. Pupillo, C. Genet, and T. W. Ebbesen, “Coherent coupling of molecular resonators with a microcavity mode,” Nature Comm. **6**, 5981 (2015).</footnote>

<footnote n="6">⁶B. Simpkins, K. P. Fears, W. J. Dressick, B. T. Spann, A. D. Dunkelberger, and J. C. Owrutsky, “Spanning strong to weak normal mode coupling between vibrational and fabry–pérot cavity modes through tuning of vibrational absorption strength,” ACS Photonics **2**, 1460–1467 (2015).</footnote>

<footnote n="7">⁷W. Ahn, J. F. Triana, F. Recabal, F. Herrera, and B. S. Simpkins, “Modification of ground-state chemical reactivity via light–matter coherence in infrared cavities,” Science **380**, 1165–1168 (2023).</footnote>

<footnote n="8">⁸K. Hirai, H. Ishikawa, T. Chervy, J. A. Hutchison, and H. Uji-i, “Selective crystallization via vibrational strong coupling,” Chem. Sci. **12**, 11986–11994 (2021).</footnote>

<footnote n="9">⁹A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, et al., “Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field,” Angew. Chem. Int. Ed. **128**, 11634–11638 (2016).</footnote>

<footnote n="10">¹⁰R. M. Vergauwe, A. Thomas, K. Nagarajan, A. Shalabney, J. George, T. Chervy, M. Seidel, E. Devaux, V. Torbeev, and T. W. Ebbesen, “Modification of enzyme activity by vibrational strong coupling of water,” Angew. Chem. Int. Ed. **58**, 15324–15328 (2019).</footnote>

<footnote n="11">¹¹J. Lather, P. Bhatt, A. Thomas, T. W. Ebbesen, and J. George, “Cavity catalysis by cooperative vibrational strong coupling of reactant and solvent molecules,” Angew. Chem. Int. Ed. **58**, 10635–10638 (2019).</footnote>

<footnote n="12">¹²K. Hirai, R. Takeda, J. A. Hutchison, and H. Uji-i, “Modulation of prins cyclization by vibrational strong coupling,” Angew. Chem. Int. Ed. **132**, 5370–5373 (2020).</footnote>

<footnote n="13">¹³J. Lather and J. George, “Improving enzyme catalytic efficiency by cooperative vibrational strong coupling of water,” J. Phys. Chem. Lett. **12**, 379–384 (2020).</footnote>

<footnote n="14">¹⁴H. Hiura, A. Shalabney, and J. George, “Vacuum-field catalysis: Accelerated reactions by vibrational ultra strong coupling,” ChemRxiv Preprint, 7234721 (2021).</footnote>

<footnote n="15">¹⁵J. Lather, A. N. Thabassum, J. Singh, and J. George, “Cavity catalysis: modifying linear free-energy relationship under cooperative vibrational strong coupling,” Chem. Sci. **13**, 195–202 (2022).</footnote>

<footnote n="16">¹⁶W. Ahn, J. F. Triana, F. Recabal, F. Herrera, and B. S. Simpkin, “Modification of ground-state chemical reactivity via light–matter coherence in infrared cavities,” Science **380**, 1165–1168 (2023).</footnote>

<footnote n="17">¹⁷J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, “Cavity casimir-polder forces and their effects in ground-state chemical reactivity,” Phys. Rev. X **9**, 021057 (2019).</footnote>

<footnote n="18">¹⁸X. Li, A. Mandal, and P. Huo, “Theory of mode-selective chemistry through polaritonic vibrational strong coupling,” J. Phys. Chem. Lett. **12**, 6974–6982 (2021).</footnote>

<footnote n="19">¹⁹C. Climent and J. Feist, “On the sn 2 reactions modified in vibrational strong coupling experiments: reaction mechanisms and vibrational mode assignments,” Phys. Chem. Chem. Phys. **22**, 23545–23552 (2020).</footnote>

<footnote n="20">²⁰T. Fukushima, S. Yoshimitsu, and K. Murakoshi, “Inherent promotion of ionic conductivity via collective vibrational strong coupling of water with the vacuum electromagnetic field,” J. Am. Chem. Soc. **144**, 12177–12183 (2022).</footnote>

<footnote n="21">²¹T. Fukushima, S. Yoshimitsu, and K. Murakoshi, “Unlimiting ionic conduction: manipulating hydration dynamics through vibrational strong coupling of water,” Chem. Sci. **14**, 1141 (2023).</footnote>

<footnote n="22">²²Y. Pang, A. Thomas, K. Nagarajan, R. M. A. Vergauwe, K. Joseph, B. Patrahau, K. Wang, C. Genet, and T. W. Ebbesen, “On the role of symmetry in vibrational strong coupling: The case of charge-transfer complexation,” Angew. Chem. Int. Ed. **59**, 10436–10440 (2020), https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.202002527.</footnote>

<footnote n="23">²³M. A. Michon and B. S. Simpkins, “Impact of cavity length non-uniformity on reaction rate extraction in strong coupling experiments,” J. Am. Chem. Soc. **146**, 30596–30606 (2024), pMID: 39466594, https://doi.org/10.1021/jacs.4c12269.</footnote>
²⁴T. W. Ebbesen, "Hybrid light-matter states in a molecular and material science perspective," *Acc. Chem. Res.* **49**, 2403–2412 (2016).

²⁵W. Ahn, I. Vurgaftman, A. Dunkelberger, J. Owrutsky, and B. Simpkins, "Vibrational strong coupling controlled by spatial distribution of molecules within the optical cavity," *ACS Photonics* **5**, 158–166 (2018).

²⁶F. Hernandez and F. Herrera, "Multi-level quantum rabi model for anharmonic vibrational polaritons," *J. Chem. Phys.* **151**, 144116 (2019).

²⁷G. D. Scholes, "Emergence of collective coherent states from strong-light coupling of disordered systems," *J. Phys. Chem. A* **125**, 6739–6750 (2021).

²⁸J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura, C. Genet, and T. Ebbesen, "Multiple rabi splittings under ultrastrong vibrational coupling," *Phys. Rev. Lett.* **117**, 153601 (2016).

²⁹A. Mandal, M. A. Taylor, B. M. Weight, E. R. Koessler, X. Li, and P. Huo, "Theoretical advances in polariton chemistry and molecular cavity quantum electrodynamics," *Chem. Rev.* **123**, 9786–9879 (2023).

³⁰T. E. Li, J. E. Subotnik, and A. Nitzan, "Cavity molecular dynamics simulations of liquid water under vibrational ultrastrong coupling," *Proc. Nat. Acad. Sci. USA* **117**, 18324–18331 (2020).

³¹L. P. Lindoy, A. Mandal, and D. R. Reichman, "Resonant cavity modification of ground-state chemical kinetics," *J. Phys. Chem. Lett.* **13**, 6580–6586 (2022).

³²P.-Y. Yang and J. Cao, "Quantum effects in chemical reactions under polaritonic vibrational strong coupling," *J. Phys. Chem. Lett.* **12**, 9531–9538 (2021).

³³C. Schäfer, M. Ruggenthaler, and A. Rubio, "Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling," *Phys. Rev. A* **98**, 043801 (2018).

³⁴D. Sidler, M. Ruggenthaler, and A. Rubio, "Numerically exact solution for a real polaritonic system under vibrational strong coupling in thermodynamic equilibrium: Loss of light-matter entanglement and enhanced fluctuations," *J. Chem. Theory Comput.* **19**, 8801–8814 (2023).

³⁵J. Sun and O. Vendrell, "Suppression and enhancement of thermal chemical rates in a cavity," *J. Phys. Chem. Lett.* **13**, 4441–4446 (2022).

³⁶Q. Yu and S. Hammes-Schiffer, "Multidimensional quantum dynamical simulation of infrared spectra under polaritonic vibrational strong coupling," *J. Phys. Chem. Lett.* **13**, 11253–11261 (2022).

³⁷J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, "Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (qed) chemistry," *Proc. Nat. Acad. Sci. USA* **114**, 3026–3034 (2017).

³⁸D. Sidler, C. Schäfer, M. Ruggenthaler, and A. Rubio, "Polaritonic chemistry: Collective strong coupling implies strong local modification of chemical properties," *J. Phys. Chem. Lett.* **12**, 508–516 (2020).

³⁹D. Sidler, T. Schnappinger, A. Obzhirov, M. Ruggenthaler, M. Kowalewski, and A. Rubio, "Unraveling a cavity induced molecular polarization mechanism from collective vibrational strong coupling," *arXiv preprint* arXiv:2306.06004 (2023).

⁴⁰M. R. Fiechter, J. E. Runeson, J. F. Lawrence, and J. O. Richardson, "How quantum is the resonance behavior in vibrational polariton chemistry?" *J. Phys. Chem. Lett.* **14**, 8261–8267 (2023).

⁴¹W. Ying and P. Huo, "Resonance theory and quantum dynamics simulations of vibrational polariton chemistry," *J. Chem. Phys.* **159**, 084104 (2023).

⁴²A. Z. Lieberherr, S. T. Furniss, J. E. Lawrence, and D. E. Manolopoulos, "Vibrational strong coupling in liquid water from cavity molecular dynamics," *J. Chem. Phys.* **158** (2023).

⁴³H. A. Kramers, "Brownian motion in a field of force and the diffusion model of chemical reactions," *Physica* **7**, 284–304 (1940).

⁴⁴R. F. Grote and J. T. Hynes, "The stable states picture of chemical reactions. ii. rate constants for condensed and gas phase reaction models," *J. Chem. Phys.* **73**, 2715–2732 (1980).

⁴⁵L. P. Lindoy, A. Mandal, and D. R. Reichman, "Quantum dynamical effects of vibrational strong coupling in chemical reactivity," *Nature Comm.* **14**, 2733 (2023).

⁴⁶S. Mondal, D. S. Wang, and S. Keshavamurthy, "Dissociation dynamics of a diatomic molecule in an optical cavity," *J. Chem. Phys.* **157**, 244109 (2022).

⁴⁷Q. Yu and J. M. Bowman, "Manipulating hydrogen bond dissociation rates and mechanisms in water dimer through vibrational strong coupling," *Nature Comm.* **14**, 3527 (2023).

⁴⁸T. E. Li, A. Nitzan, and J. E. Subotnik, "Cavity molecular dynamics simulations of vibrational polariton-enhanced molecular nonlinear absorption," *J. Chem. Phys.* **154**, 094124 (2021).

⁴⁹T. E. Li, A. Nitzan, and J. E. Subotnik, "Collective vibrational strong coupling effects on molecular vibrational relaxation and energy transfer: Numerical insights via cavity molecular dynamics simulations," *Angew. Chem. Int. Ed.* **60**, 15533–15540 (2021).

⁵⁰T. E. Li, A. Nitzan, and J. E. Subotnik, "Energy-efficient pathway for selectively exciting solute molecules to high vibrational states via solvent vibration-polariton pumping," *Nature Comm.* **13**, 4203 (2022).

⁵¹T. E. Li, A. Nitzan, and J. E. Subotnik, "Polariton relaxation under vibrational strong coupling: Comparing cavity molecular dynamics simulations against fermi's golden rule rate," *J. Chem. Phys.* **156**, 134106 (2022).

⁵²T. E. Li and S. Hammes-Schiffer, "QM/MM modeling of vibrational polariton induced energy transfer and chemical dynamics," *J. Am. Chem. Soc.* **145**, 377–384 (2022).

⁵³D. Laage and J. T. Hynes, "A molecular jump mechanism of water reorientation," *Science* **311**, 832–835 (2006).

⁵⁴D. Laage and J. T. Hynes, "On the molecular mechanism of water reorientation," *J. Phys. Chem. B* **112**, 14230–14242 (2008).

⁵⁵D. Laage, G. Stirnemann, F. Sterpone, R. Rey, and J. T. Hynes, "Reorientation and allied dynamics in water and aqueous solutions," *Ann. Rev. Phys. Chem.* **62**, 395–416 (2011).

⁵⁶D. M. Wilkins, D. E. Manolopoulos, S. Pipolo, D. Laage, and J. T. Hynes, "Nuclear quantum effects in water reorientation and hydrogen-bond dynamics," *J. Phys. Chem. Lett.* **8**, 2602–2607 (2017).

⁵⁷A. Gomez, Z. A. Piskulich, W. H. Thompson, and D. Laage, "Water diffusion proceeds via a hydrogen-bond jump exchange mechanism," *J. Phys. Chem. Lett.* **13**, 4660–4666 (2022).

⁵⁸T. E. Li, A. Nitzan, S. Hammes-Schiffer, and J. E. Subotnik, "Quantum simulations of vibrational strong coupling via path integrals," *J. Phys. Chem. Lett.* **13**, 3890–3895 (2022).

⁵⁹K. Sun and R. F. Ribeiro, "Theoretical formulation of chemical equilibrium under vibrational strong coupling," *Nature Comm.* **15**, 2405 (2024).

⁶⁰I. R. Craig and D. E. Manolopoulos, "Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics," *J. Chem. Phys.* **121**, 3368–3373 (2004).

⁶¹S. Habershon, T. E. Markland, and D. E. Manolopoulos, "Competing quantum effects in the dynamics of a flexible water model," *J. Chem. Phys.* **131**, 024501 (2009).

⁶²L. Verlet, "Computer "experiments" on classical fluids. I. thermodynamical properties of lennard-jones molecules," *Phys. Rev.* **159**, 98 (1967).

⁶³M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, "Efficient stochastic thermostatting of path integral molecular dynamics," *J. Chem. Phys.* **133**, 124104 (2010).

⁶⁴M. Rossi, M. Ceriotti, and D. E. Manolopoulos, "How to remove the spurious resonances from ring polymer molecular dynamics," *J. Chem. Phys.* **140** (2014).

⁶⁵T. E. Li, "Cavity molecular dynamics simulations tool sets," (2020), https://github.com/TaoELi/cavity-md-ipi.

⁶⁶V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, et al., "i-pi 2.0: A universal force engine for advanced molecular simulations," *Comp. Phys. Comm.* **236**, 214–223 (2019).

⁶⁷A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, "LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales," *Comp. Phys. Comm.* **271**, 108171 (2022).

⁶⁸M. P. Allen and D. J. Tildesley, *Computer simulation of liquids* (Oxford university press, 2017).

⁶⁹J. Richardi, P. Fries, and C. Millot, "Fast hybrid methods for the simulation of dielectric constants of liquids," *J. Mol. Liq.* **117**, 3–16 (2005).

⁷⁰S. W. d. Leeuw, J. W. Perram, and E. Smith, "Simulation of electrostatic systems in periodic boundary conditions. i. lattice sums and dielectric constants," *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences* **373**, 27–56 (1980).

⁷¹W. Zhang and A. C. Van Duin, "Second-generation reaxff water force field: Improvements in the description of water density and oh-anion diffusion," *J. Phys. Chem. B* **121**, 6021–6032 (2017).

⁷²I.-C. Yeh and G. Hummer, "System-size dependence of diffusion coefficients," *J. Chem. Phys.* **154**, 094124 (2021).
cients and viscosities from molecular dynamics simulations with periodic boundary conditions,” J. Phys. Chem. B **108**, 15873–15879 (2004).

<footnote n="73">⁷³J. de la Fuente, R.Spezia, and R. Vuilleumier, “Spectroscopic properties under vibrational strong coupling in disorder matter from path-integral monte carlo simulations,” J. Chem. Phys. **161**, 184114 (2024).</footnote>

<footnote n="74">⁷⁴A. Eltareb, G. E. Lopez, and N. Giovambattista, “Nuclear quantum effects on the thermodynamic, structural, and dynamical properties of water,” Phys. Chem. Chem. Phys. **23**, 6914–6928 (2021).</footnote>

<footnote n="75">⁷⁵M. Piejko, B. Patrahau, K. Joseph, C. Muller, E. Devaux, T. W. Ebbesen, and J. Moran, “Solvent polarity under vibrational strong coupling,” J. Am. Chem. Soc. **145**, 13215–13222 (2023).</footnote>

<footnote n="76">⁷⁶M. Piejko, B. Patrahau, K. Joseph, C. Muller, E. Devaux, T. W. Ebbesen, and J. Moran, “Retraction of solvent polarity under vibrational strong coupling,” J. Am. Chem. Soc. **146**, 30716–30716 (2024), pMID: 39436019, https://doi.org/10.1021/jacs.4c12901.</footnote>

<footnote n="77">⁷⁷D. Sidler, T. Schnappinger, A. Obzhirov, M. Ruggenthaler, M. Kowalewski, and A. Rubio, “Unraveling a cavity-induced molecular polarization mechanism from collective vibrational strong coupling,” J. Phys. Chem. Lett. **15**, 5208–5214 (2024).</footnote>

<footnote n="78">⁷⁸T. Morresi, L. Paulatto, R. Vuilleumier, and M. Casula, “Probing an-harmonic phonons by quantum correlators: A path integral approach,” J. Chem. Phys. **154**, 224108 (2021).</footnote>