# Spectroscopic properties under vibrational strong coupling in disordered matter from path-integral Monte Carlo simulations

Jaime de la Fuente Diez, Riccardo Spezia, Rodolphe Vuilleumier

► To cite this version:

Jaime de la Fuente Diez, Riccardo Spezia, Rodolphe Vuilleumier. Spectroscopic properties under vibrational strong coupling in disordered matter from path-integral Monte Carlo simulations. The Journal of Chemical Physics, 2024, 161 (18), pp.184114. 10.1063/5.0226938. hal-05193397

HAL Id: hal-05193397

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

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
Spectroscopic Properties under Vibrational Strong Coupling in Disordered Matter
from Path-Integral Monte Carlo Simulations.

Jaime de la Fuente Diez,¹ Riccardo Spezia,², a) and Rodolphe Vuilleumier¹, b)

¹) PASTEUR, Département de chimie, École normale supérieure,
PSL University, Sorbonne Université, CNRS, 75005 Paris,
France

²) Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS,
4 Place Jussieu, 75005 Paris (France)

(Dated: 30 July 2025)

Vibrational strong coupling (VSC), the strong coupling between a Fabry-Perrot cavity and molecular vibrations at mid-infrared frequencies, has received important attention in last years due to its capacity of modifying both vibrational spectra and chemical reactivity. VSC is a collective effect and in this work we introduce Path Integral Monte Carlo (PIMC) simulations that not only take into account the quantum character of the molecular vibrations and of the optical resonance of the cavity, but also reproduces this collective behaviour by considering multiple replicas of the molecular system. Moreover, we show that it is possible to extract from the PIMC simulations the decomposition of the hybrid optical and molecular states in terms of the uncoupled modes. On a model system of an ensemble of disordered Morse oscillators coupled to a single cavity through the Pauli-Fierz Hamiltonian, PIMC can retrieve known features obtained from analytical modes such as the Tavis-Cummings model, and obtain a very close agreement with exact diagonalization for small number of Morse oscillators. We also find that notwithstanding the anhamonic character of the Morse oscillators, the collective mode coupled to the cavity behaves as a harmonic oscillator, following the quantum central limit theorem.

Keywords: Vibrational Strong Coupling — Path Integral Theory — Disorder — Spectroscopy — Light-matter Hybridization

<footnote>a) Electronic mail: riccardo.spezia@sorbonne-universite.fr</footnote>

<footnote>b) Electronic mail: rodolphe.vuilleumier@ens.psl.eu</footnote>
# I. INTRODUCTION

Understanding how vibrational strong coupling (VSC) modifies the spectroscopic and reactive properties of matter has a quickly growing interest in physical chemistry.¹ This modifications are due to the coupling between the zero-point fluctuations of the electric field of an optical cavity and a molecular vibrational transition, thus generating hybrid light-matter states with particular properties such as enhanced energy transport independent from the distance between molecules.²⁻⁴

Experimentally, this hybrid light-matter state is clearly observable from vibrational spectroscopy: in fact, when the matter is in a Fabry-Perot cavity with the correct size, two vibro-polaritonic bands (called $P^+$ and $P^−$) appear in the transmittance. This is due to the resonance between the cavity and the matter vibrational modes and it was shown for different systems that it is possible to tune the cavity size to be resonant with vibrational molecular modes in a wide range of frequencies, from ~840 $cm^{-1}$ for a Si-C stretching mode⁵,⁶ to 3400 $cm^{-1}$ for coupling to an O-H stretching mode.⁷⁻⁹ Upon coupling, significant Rabi splittings of up to 20% of the mode energy have been reported, which are also larger than the full width at half maximum (FWHM) of the targeted vibrational bands, proving that the systems are in the strong or ultra strong coupling regime.¹⁰

Recently, it has been shown experimentally that chemical reactions can be affected by VSC: first studies, by Ebbesen and co-workers, have shown that the reaction rate constant decreases when the cavity is in resonance with some vibrational modes of the reactants.⁵⁻⁷ In few cases the reaction was accelerated by VSC.⁸ Ebbesen and co-workers have recently summarized the different experimental findings.¹¹ More recently, also the thermodynamics of a chemical reaction were modified by VSC.¹² Other groups have then observed a modification of chemical reactivity when the reaction is done under VSC conditions.¹³,¹⁴ Other experiments have also pointed out that ion conductivity can be modified under VSC.⁹

Theoretically, different works have been published in recent years, many of them motivated by the aforementioned experiments done by Ebbesen and co-workers on the modification of reaction rate constants under VSC.¹⁵⁻²³. In most cases, only one molecule was considered for the matter which is then coupled to the cavity, while the rest of the system is sometimes added as a bath: this approach was able to show how chemical kinetics is enhanced or suppressed as a function of
the intensity of the friction.²⁴ Furthermore, a quantum approach using transition state theory by Yang and Cao considered one explicit reactive mode in the cavity.²⁵ Similarly, quantum dynamics was done using different models, but typically considering only a single system coupled to the cavity.²²,²⁶,²⁷ It is also possible to directly work with the nuclear-electron Hamiltonian, as done by Rubio and co-workers.²⁸⁻³⁰

Another theoretical approach was proposed by Li et al. which allows to study how VSC modifies condensed phase systems.¹⁷,¹⁸,³¹ This is based on classical molecular dynamics where a model cavity is introduced interacting with the dipoles of the molecules present in the simulation box and using a point charge approximation for the electrostatic force field. In this way, it was possible to study liquid water¹⁸, carbon dioxide¹⁷ and also polaritonic relaxation dynamics.³¹,³² It was extended to include QM/MM force field³³ and Ring Polymer Molecular Dynamics to take into account nuclear quantum effects.²³,³⁴ In this approach, the different matter modes coupled with the cavity are explicitly considered and they also interact with each other. However, it is difficult to disentangle from such simulations the different contributions to the polaritonic bands and to the reaction rate.

Despite the fact that theoretical studies are very active and variate, there is not a clear picture when considering an assembly of molecules to understand and predict the different chemical reactivity, although VSC is a collective effect. In particular, different aspects, like the role and nature of dark states and the cooperativity between matter modes are still open issues.

In the present work, we introduce a Path Integral Monte Carlo (PIMC) simulation³⁵ of an optical cavity coupled to an ensemble of replicas of the molecular system. It allows to include the quantum character of the electromagnetic field in the cavity together with the quantum nature of molecular vibrations, while, by using replicas of the molecular systems, it accounts for collective effects among molecular states.

To obtain the infrared spectrum from PIMC, we have extended the method presented recently by Morresi et al.³⁶. It further allows to describe the modes of the system under VSC in terms of the modes of the non-coupled system and obtain detailed information on the light-matter hybridization. We apply this approach to the simulation of an ensemble of Morse oscillators coupled to a Fabry-Perrot cavity and study the dependence of the vibrational states on the number of explicit modes in the model and of the Rabi splitting. We also include disorder among the Morse
oscillators to mimic inhomogenous broadening.

Although simplified, this model is already very rich and different studies have examined the role of disorder. For example, Random Matrix Monte Carlo simulations have been performed by Scholes<sup>16</sup> on a system composed by molecular transitions exhibiting frequency disorder coupled to an imperfect cavity which also displays frequency disorder. The cavity disorder was shown not to break the effective mode approximation such that the matter can be considered coupled to a single effective cavity mode due to the highly coherent nature of the cavity vacuum fluctuations. Also the onset of the strong coupling regime happened after a certain critical number of molecules become coherently coupled to the effective cavity mode, effectively raising the Rabi Splitting beyond the width of the molecular disorder. Later, the disorder was considered in studying a quantum two-level system under VSC by Botzung et al.,<sup>37</sup> showing the role of resulting semi-localized dark states in transport properties. Recently, this approach was extended by Du and Yuen-Zhou to build a simple model of catalysis, starting from a Tavis-Cumming model for the VSC. Also in this case, such semi-localized dark states are responsible in enhancing the channels of vibrational cooling.<sup>38</sup> The Tavis-Cummings Hamiltonian was also studied through the Green function framework in the limit of infinite number of modes by Engelhardt et al.<sup>39</sup>.

The Tavis-Cummings Hamiltonian used in these studies is approximate, as it is restricted to single excitations, and can equally apply to coupled harmonic oscillators or to a harmonic cavity coupled to an ensemble of two-states systems, a very anharmonic situation. The harmonic problem can however be solved exactly for an assembly of modes, see Ref.<sup>25</sup>.

This problem is examined here without approximations and with anharmonic molecular vi-
brations, describing a typical OH bond, up to a 100 replicas. For comparison, we also present
results from the Tavis-Cummings model, both numerical at finite number of modes and from the
Green function in the limit of infinite number of modes. To evaluate the accuracy of the PIMC
calculations we also performed exact diagonalization of the full quantum Hamiltonian for a small
number of Morse oscillators coupled to the optical cavity. We finally examine the properties of the
collective dipole composed of the dipoles of the replicas of the molecular vibration.

The remainder of the article is organised as follows. In section II we recall the general theory,
and in particular: the Hamiltonians used to describe the vibrational strong coupling (section II A),
the Path Integral theory for the sake of fixing the notations (section II B), the extraction of vibra-
tional modes from Monte Carlo simulations (section II C). For direct comparison we introduce the Green function theory on the Tavis-Cummings Hamiltonian in section II D. We then provide computational details (section III) on the Monte Carlo approach we employed, together with the calculation of error bars, and the exact diagonalization of the anharmonic Pauli-Fierz Hamiltonian for further validation of the PIMC results. Section IV presents the results, firstly without disorder as a benchmark (sections IV A) and then including disorder in the matter by employing a distribution of initial matter mode frequencies (section IV B). In section IV C we discuss the resulting light-matter hybridization and finally we give some insights on the behavior for a larger number of matter modes (section IV D). At the end we give some general conclusions and perspectives for further studies (section V).

## II. THEORY

### A. Hamiltonian for vibrational strong coupling

To describe the coupling between the cavity and the matter degrees of freedom, we consider two possible Hamiltonians, one in which the coupling between radiation and matter is linear, akin to the Tavis-Jaynes-Cummings Hamiltonian,⁴⁰ and the Pauli-Fierz (PF) Hamiltonian²⁸,⁴¹,⁴² which can be derived from Quantum Electro Dynamics (QED) theory and gives rise to a dipole self-interaction. Also, the Tavis-Cummings Hamiltonian was considered because it can be numerically diagonalized and studied analytically³,³⁷,⁴³.

In general, we can write the total Hamiltonian as follows:

$$ \mathcal{H} = \mathcal{H}_m + \mathcal{H}_{VSC} \qquad (1) $$

where $\mathcal{H}_m$ is the Hamiltonian of the matter and $\mathcal{H}_{VSC}$ couples the matter with the cavity. Here we consider the matter as composed by a set of $N$ independent oscillators all coupled with the cavity:

$$ \mathcal{H}_m(q, p) = \sum_{i}^{N} \left( \frac{p_i^2}{2m} + V_i(q) \right) \qquad (2) $$
For the matter potential, $V_i$, we have considered a Morse function:

$$V_i^{\text{Morse}}(q_i) = D_e (1 - e^{-a_i q_i})^2 \quad (3)$$

where $D_e$ is constant and the $a_i$ are related to the frequencies $\omega_i$ at the bottom of the well by $a_i = \sqrt{\frac{m}{2D_e}}\omega_i$ and $q_i$ is the one dimensional coordinate associated with the i-th oscillator. This allows to define a frequency associated to each of the modes described by the Morse potential.

In both cases we consider that all the modes have the same mass, $m$. However, the vibrational frequencies of the matter are considered in two ways: (i) all the modes are set to be identical, or (ii) we model static disorder with a frequency distribution centered around a given value, $\omega_0$. In the latter case we have employed a Gaussian distribution with a given standard deviation, $\sigma$, which mimics the frequency broadening in condensed phase. We also considered an equally spaced uniform distribution, with $N$ frequencies centered on a given value in the interval $[\omega_0 - \Delta\omega/2, \omega_0 + \Delta\omega/2]$.

To describe the light-matter coupling, we used two representations. First, the linear Hamiltonian, which is the classical counterpart of the quantum Tavis-Cummings Hamiltonian:

$$\mathcal{H}_{\text{VSC,linear}} = \frac{p_c'^2}{2} + \frac{1}{2}\omega_c^2 q_c'^2 + \gamma' q_c' \sum_i^N \mu_i \quad (4)$$

where $p_c'$ and $q_c'$ are the momentum and position of the cavity mode in mass-weighted coordinates, $\omega_c$ is the cavity frequency, $\gamma'$ the coupling between the cavity and the matter (in mass-weighted coordinates) and $\mu_i$ the dipole of each matter mode. Note that this Hamiltonian can be unbound since the last term can assume negative values.41,42 This occurs when the coupling parameter becomes very strong, which is the case when the coupling constant is increased to compensate for a low number of molecules. In this situation, if the Rabi splitting is large enough, the lower polariton could go to a negative energy.

We have then considered another Hamiltonian for the VSC, the Pauli-Fierz Hamiltonian which derives from Quantum Electro Dynamics (QED).19,28,32 Its classical counterpart reads:
$$ \mathcal{H}_{VSC,PF} = \frac{p_c'^2}{2} + \frac{1}{2}\omega_c^2 \left( q_c' + \frac{\gamma}{\omega_c^2} \sum_{i}^{N} \mu_i \right)^2 \quad (5) $$

Here, as the Pauli-Fierz Hamiltonian includes a quadratic self interaction term, the energies of the eigenstates are bounded from below, meaning that the lower polariton will never be assigned a negative energy in spite of the value of the coupling constant $\gamma$. We have then introduced a mass for the cavity, $m_c$, in order to treat it in the same footing as the matter modes and be able to write the Hamiltonian in the Path Integral formalism. This amounts just to a re-scaling in both Hamiltonians by the cavity mass:

$$ \mathcal{H}_{VSC,linear} = \frac{p_c^2}{2m_c} + \frac{1}{2}m_c\omega_c^2 q_c^2 + \gamma q_c \sum_{i}^{N} \mu_i \quad (6) $$

$$ \mathcal{H}_{VSC,PF} = \frac{p_c^2}{2m_c} + \frac{1}{2}m_c\omega_c^2 \left( q_c + \frac{\gamma}{m_c\omega_c^2} \sum_{i}^{N} \mu_i \right)^2 \quad (7) $$

In both cases the dipoles of the matter mode are considered here simply as $\mu_i = q_i e_i$ where $e_i$ is the unit charge and $q_i$ the matter coordinate (which is a displacement from the equilibrium condition, here fixed to $q_i^{eq} = 0$).

Note that the Pauli-Fierz Hamiltonian includes a quadratic term, mediated by the coupling to the cavity $\gamma$, that considers a direct interaction between matter modes through their positions. This term keeps the Hamiltonian Gauge-invariant. If this term is dropped, however, we revert to linear Hamiltonian, where the interaction is linear in the coupling constant$^{28,40}$.

The coupling term, $\gamma$, is phenomenologically estimated. In fact, if we consider one cavity mode coupled to one molecular mode, described as a harmonic oscillator, and we rescale the cavity coordinates by the mass of the matter oscillator (whose value can be established arbitrarily), the Rabi splitting of the VSC-targeted mode is given by:

$$ \Omega = \frac{\gamma}{\sqrt{m_c m \omega_c}} \quad (8) $$
which, for the general case of N matter modes, becomes:

$$ \Omega_o = \frac{\gamma}{\sqrt{m_c m \omega_c}} \sqrt{N} \qquad (9) $$

This means that, if we choose an observed Rabi splitting $\Omega_o$, fixing the number of modes and the cavity frequency, we can derive the coupling. In other words, the coupling constant is determined to reproduce the experimental Rabi splitting of the vibrational transition under VSC for a given number of matter modes, N.

Note that, in principle, the coupling can be determined through the properties of the cavity. In fact, it reads (in mass-weighted coordinates) as follows:¹⁸

$$ \gamma = \sqrt{\frac{m_c}{\varepsilon_0 V_c}} \omega_c \qquad (10) $$

where $\varepsilon_0$ is the vacuum permittivity, $\omega_c$ is the resonant frequency of the cavity, and $V_c$ is the mode quantization volume inside of the cavity, that can be approximated as the geometric volume of the cavity. Note that then $\gamma\sqrt{N}$ is intensive as expected.

Therefore, in all our simulations, and similarly to what done in other theoretical works,¹⁶⁻¹⁸ the coupling constant will be obtained phenomenologically and it will result artificially large to compensate for the limitations in the number of matter modes explicitly modeled. This point will be further discussed later.

## B. Path Integral Monte Carlo

Path Integral theory is an exact theory developed by Feynman⁴⁴ in which the quantum partition function is written in the position basis as:

$$ \mathcal{Z} = \lim_{P \to \infty} \left( \frac{mP}{2\pi\beta\hbar^2} \right)^{P/2} \int dx_1 \dots dx_P e^{-\beta\phi(x_1 \dots x_P)} \qquad (11) $$
with

$$ \phi(x_1 \dots x_P) = \sum_{k=1}^{P} \left[ \frac{1}{2} m \omega_P^2 (x_k - x_{k+1})^2 + \frac{1}{P} U(x_k) \right] \quad (12) $$

where $P$ is the number of replicas (called often beads), $m$ the corresponding mass, $\beta = 1/k_bT$ (where $k_b$ is the Boltzmann constant and $T$ the temperature of the system), $\{x_k\}_{k=1,...,P}$ the positions of the beads and $\omega_P = \frac{\sqrt{P}}{\beta\hbar}$ and $x_1 = x_{P+1}$.

In real applications of the theory, the number of beads, $P$, is a finite number. For the present case, the resulting Pauli-Fierz Hamiltonian (combining equations 2 and 7) can be re-written using the Trotter factorization as

$$ \mathcal{H}_P = \sum_{i=1}^{N} \sum_{j=1}^{P} \left[ \frac{1}{2} m_i \omega_P^2 (q_i^{(j)} - q_i^{(j-1)})^2 \right]_{q_i^{(P+1)}=q_i^{(1)}} + \sum_{j=1}^{P} V(q_1^{(j)}, q_2^{(j)}, \dots, q_N^{(j)}) \\ + \sum_{j=1}^{P} \frac{m_c \omega_c^2}{2} \left( q_c^{(j)} + \frac{\gamma}{m_c \omega_c^2} \mu_T^{(j)} \right)^2 \quad (13) $$

where the index $j$ runs over the $P$ "beads" and the bead frequency, $\omega_P$ is $\omega_P^2 = \frac{P}{\beta^2\hbar^2}$. The quantum evolution of the system is mapped in the classical evolution of a set of correlated replicas of the system.

The action corresponding to this Hamiltonian is then:

$$ \mathcal{A} = \frac{1}{P} \sum_{i=1}^{N} \sum_{j=1}^{P} \left[ \frac{m_i P^2}{2\beta^2 \hbar^2} (q_i^{(j)} - q_i^{(j-1)})^2 \right]_{q_i^{(P+1)}=q_i^{(1)}} + \frac{1}{P} \sum_{j=1}^{P} V(q_1^{(j)}, q_2^{(j)}, \dots, q_N^{(j)}) \\ + \sum_{j=1}^{P} \frac{m_c \omega_c^2}{2P} \left( q_c^{(j)} + \frac{\gamma}{m_c \omega_c^2} \mu_T^{(j)} \right)^2 \quad (14) $$

This will be used to determine the acceptance probability of a given configuration in the Path Integral Monte Carlo (PIMC) simulations. Analogous expressions are obtained for the linear coupling Hamiltonian of Equation 6.

Details on the Monte Carlo procedure will be given in Section III.
### C. Vibrational states from Monte Carlo simulations

To obtain vibrational states and the corresponding spectra, we extended the method recently proposed by Morresi et al.³⁶ for Path Integral Molecular Dynamics to Path Integral Monte Carlo.

In the Path Integral theory, the force constant matrix is defined as:

$$ \langle\langle V_{i_1, i_2} \rangle\rangle = \left\langle\left\langle \frac{1}{P^2} \sum_{j_1, j_2=1}^{P} \frac{\partial^2 V}{\partial x_{i_1}^{(j_1)} \partial x_{i_2}^{(j_2)}} \right\rangle\right\rangle \quad (15) $$

where $i_1$ and $i_2$ indicate the matter degrees of freedom, $j_1$ and $j_2$ run over the $P$ beads and the double brackets indicate a first average over the members of the ring polymer, obtaining centroid properties, and the second average over all the ensemble.

We can now define the quantum force fluctuations by the Kubo-transformed correlation function:

$$ \tilde{c}_{F_{i_1}F_{i_2}} = \int d^f \mathbf{p} \int d^f \mathbf{x} \frac{e^{-\beta_p \mathcal{H}_p(\mathbf{x}, \mathbf{p})}}{\mathcal{Z}} F_{i_1}(\mathbf{x}) F_{i_2}(\mathbf{x}) = \langle\langle F_{i_1} F_{i_2} \rangle\rangle \quad (16) $$

where $\mathcal{Z}$ is the partition function and $\beta_P = \frac{1}{k_B T_P}$

In the work by Morresi et al.³⁶, it is shown that

$$ \langle\langle F_{i_1} F_{i_2} \rangle\rangle = \frac{1}{\beta} \langle\langle V_{i_1, i_2} \rangle\rangle \quad (17) $$

where now $\tilde{c}_{F_{i_1}F_{i_2}} = \langle\langle F_{i_1}F_{i_2} \rangle\rangle$. This means that we can obtain the Path Integral force constant matrix by computing the Kubo transformed force correlation function.

So, since the zero-time Kubo-transformed autocorrelation fulfills the equipartition theorem as well as for the effective ring-polymer systems, we can write a Path Integral generalized eigenvalue problem to evaluate the vibrational modes:
$$ \langle\langle F_{i_1} F_{i_2} \rangle\rangle = \omega_{FF,i_3}^2 \langle\langle p_{i_1} p_{i_2} \rangle\rangle Y_{i_2,i_3} \quad (18) $$

$$ = \frac{1}{\beta} \omega_{i_3}^2 m_{i_1} \delta_{i_1, i_2} Y_{i_2, i_3} \quad (19) $$

where *F* and *p* are averaged over the whole ring polymer and *Y* contains the normal modes as a linear combination of original modes.

In addition to the force-force correlation function, the displacement-displacement correlation function also allows to estimate the PI force constant matrix (see again Morresi et al. for more details³⁶):

$$ \langle\langle\delta x_{i_1} \delta x_{i_2}\rangle\rangle = \left\langle\left\langle\left(\frac{1}{P}\sum_{j_1=1}^{P} \delta x_{i_1}^{(j_1)}\right)\left(\frac{1}{P}\sum_{j_2=1}^{P} \delta x_{i_2}^{(j_2)}\right)\right\rangle\right\rangle \quad (20) $$

$$ \langle\langle\delta x_{i_1} \delta x_{i_2}\rangle\rangle^{-1} \approx \beta \langle\langle V_{i_1, i_2} \rangle\rangle^{-1} \quad (21) $$

where in this case the relation is only approximate in the limit of small displacements from equilibrium, leading to the generalized eigenvalue problem:

$$ [\langle\langle \delta \mathbf{x}^T \delta \mathbf{x} \rangle\rangle^{-1}]_{i_1, i_2} W_{i_2, i_3} = \omega_{\delta, i_3}^2 [\langle\langle \dot{\mathbf{x}}^T \rangle\rangle^{-1}]_{i_2, i_3} W_{i_2, i_3} \quad (22) $$

Now the normal mode patterns are contained inside the *W* matrix, as for *Y*. Using again the equipartition theorem for the velocities, the eigenvalue equation can be reformulated as:

$$ [\langle\langle \delta x_{i_1} \delta x_{i_2} \rangle\rangle]^{-1} W_{i_2, i_3} = \omega_{i_3}^2 m_{i_1} \beta W_{i_2, i_3} \quad (23) $$

Note that now the terms needed to solve eigenvalue problems of Equations (19) and (23) can be obtained also from Monte Carlo sampling, calculating the averages along the accepted configuration space. Despite the fact that both methods reflect two apparently equivalent ways to compute the mode energies of the system, they provide different results. The solution of Eq. (19) provides
the ground state energies from the system, while using Eq. (23) allows us to probe its lowest lying excitation and thus it will provide the peaks that are observed on a vibrational spectrum.³⁶ Therefore, only equation (23), the one that produces the transition energies, will be considered from now on.

### D. Green Function analysis of the polaritonic states.

We can obtain valuable insights from the analysis of a simplified case of the Pauli-Fierz Hamiltonian, which is the Tavis-Cummings model,⁴⁰ in which a set of two-level emitters are coupled to a 2-level radiation field. The only interaction between particles is the direct matter-cavity and terms between matter modes that are of second order in the coupling constant $\gamma$ in (7) are neglected. Note that in this case, the 2-level emitters are assimilated to be vibrational modes of the matter under strong coupling

$$ \hat{H}_{TC} = \hbar \sum_{i}^{N} \omega_i \hat{a}_i^\dagger \hat{a}_i + \hbar \omega_c \hat{a}_c^\dagger \hat{a}_c + g \sum_{i}^{N} (\hat{a}_c \hat{a}_i^\dagger + \hat{a}_i \hat{a}_c^\dagger) \quad (24) $$

where $g = \frac{\Omega_0}{2\sqrt{N}}$ is now the coupling constant, $\hat{a}_c, \hat{a}_c^\dagger$ denote the creation and annihilation operators for the cavity mode of frequency $\omega_c$ and $\hat{a}_i, \hat{a}_i^\dagger$ denote analogous operators for the molecular mode of frequency $\omega_i$.

This Hamiltonian conserves the excitation number $\sum_i^N \hat{a}_i^\dagger \hat{a}_i + \hat{a}_c^\dagger \hat{a}_c$. It can be thus diagonalized numerically within the single excitation subspace for a number of matter modes ($N$) which is in the range of what we can simulate via PIMC. This can provide us a direct comparison and benchmark for simulations results.

Moreover, quasi-analytical values for the polaritons and mode contributions to the polaritons can be obtained by analyzing the Green's function of the system.⁴³,⁴⁵ In fact, when the system is submitted to a unit perturbation at time $t_0$, it will evolve in time according to:

$$ \left( i\hbar \frac{\partial}{\partial t} - \hat{H}_{TC} \right) G(t, t_0) = \delta(t - t_0) \quad (25) $$
where $G(t, t_0)$ is the Green's function of the system under the perturbation at $t_0$ and $\delta$ is the delta function. The formal solution to (25) is:

$$ \tilde{G}(\omega) = \frac{I}{\hbar\omega - \hat{H}_{TC}} = \sum_n \frac{|n\rangle\langle n|}{\hbar\omega - \epsilon_n} \qquad (26) $$

where $\tilde{G}(\omega)$ is the Green's function in the frequency domain and $n$ runs over the entire spectrum of $\hat{H}_{TC}$, $\epsilon_n$ being the energies of the eigenstates of that Hamiltonian. Considering $\omega \in \mathbb{C}$, $\omega \to \omega + is$, and taking the limit where the imaginary part vanishes, we obtain:

$$ \lim_{s \to 0} (\tilde{G}(\omega + is)) = \sum_n \frac{\hbar\omega - \epsilon_n}{(\hbar\omega - \epsilon_n)^2} |n\rangle\langle n| - i\pi \sum_n |n\rangle\langle n| \delta(\hbar\omega - \epsilon_n) \qquad (27) $$

The imaginary part of the Green's function allows to access information about the composition of a given eigenstate of the Hamiltonian in the basis of uncoupled oscillators, where $|c\rangle$ and $|i\rangle$ are the cavity oscillator and the i-th matter mode respectively:

$$ \mathrm{Im}(\tilde{G}_{c,c}(\omega)) = -\pi \sum_n \delta(\hbar\omega - \epsilon_n) \langle c|n\rangle \langle n|c\rangle = -\pi \sum_n \delta(\hbar\omega - \epsilon_n) |C_{c,n}|^2 \qquad (28) $$

Similarly:

$$ \mathrm{Im}(\tilde{G}_{c,i}(\omega)) = -\pi \sum_n \delta(\hbar\omega - \epsilon_n) \langle c|n\rangle \langle n|i\rangle = -\pi \sum_n \delta(\hbar\omega - \epsilon_n) C_{c,n} C_{i,n}^* \qquad (29) $$

The pre-factors of the delta functions contain the information about the composition of each state of the coupled system. In order to obtain the expressions for $\tilde{G}_{c,c}$ and $\tilde{G}_{c,i}$, starting from (26), we can write:

$$ \begin{pmatrix} \hbar\omega - \hbar\omega_c & -g & -g & -g & \cdots \\ -g & \hbar\omega - \hbar\omega_1 & 0 & 0 & 0 \\ -g & 0 & \hbar\omega - \hbar\omega_2 & 0 & 0 \\ -g & 0 & 0 & \hbar\omega - \hbar\omega_3 & 0 \\ \vdots & 0 & 0 & 0 & \cdots \end{pmatrix} = \mathbb{I} \qquad (30) $$
Equation (30) allows to obtain expressions for the components of the Green function $\tilde{G}_{c,c}$ and $\tilde{G}_{c,i}$:

$$ (\hbar\omega - \hbar\omega_c) \tilde{G}_{c,c} - g\sum_i \tilde{G}_{c,i} = 1 \qquad (31-a) $$

$$ (\hbar\omega - \hbar\omega_c) \tilde{G}_{c,i} - g\sum_j \tilde{G}_{i,j} = 0 \qquad (31-b) $$

$$ (\hbar\omega - \hbar\omega_i) \tilde{G}_{c,i} - g\tilde{G}_{c,c} = 0 \qquad (31-c) $$

$$ (\hbar\omega - \hbar\omega_j) \tilde{G}_{i,j} - g\tilde{G}_{c,i} = \delta_{ij} \qquad (31-d) $$

Through the combination of (31-c) and (31-a), the components of the Green's function required at equations (28) and (29) can be found:

$$ \tilde{G}_{c,c} = \frac{1}{\hbar\omega - \hbar\omega_c - g^2 \int d\omega' \rho(\omega') \frac{1}{\hbar\omega - \hbar\omega'}} \qquad (32-a) $$

$$ \tilde{G}_{c,i} = \frac{g\tilde{G}_{c,c}}{\hbar\omega - \hbar\omega_i} \qquad (32-b) $$

where a continuous frequency distribution density $\rho(\omega)$ has been assumed for simplicity. The integral in (32-a) can be written as:

$$ \int d\omega' \rho(\omega') \frac{1}{\hbar\omega - \hbar\omega'} = \bar{\rho}(\omega) - i\pi\rho(\omega) \qquad (33) $$

where $\bar{\rho}$ is the principal part. As the polaritons are well outside the density distribution of the uncoupled matter, the imaginary part can be set to zero, and the real part can be approximated as $\bar{\rho}(\omega) \approx \frac{N}{\hbar\omega_{\pm} - \hbar\bar{\omega}}$, $\bar{\omega}$ being the average frequency of the uncoupled matter distribution. This approximation neglects the effect of frequency distributions. The case of weak coupling, when the line widths are bigger than the Rabi splitting, is excluded here. Therefore, equations (32) become:
$$ \tilde{G}_{c,c} = \frac{1}{\hbar\omega - \hbar\omega_c - g^2 \frac{N}{\hbar\omega_{\pm} - \hbar\bar{\omega}}} \quad (34-a) $$

$$ \tilde{G}_{c,i} = \frac{g \tilde{G}_{c,c}}{\hbar\omega - \hbar\omega_i} \quad (34-b) $$

Note that taking into account Eq. (26), we see that the energies for the eigenstates of the Hamiltonian will be the poles of $\tilde{G}$. If we assume the eigenstates to be far from the bulk of the frequency distribution, as in Eq. (33), then the poles correspond to the polaritons. The poles of (34-a) are:

$$ \hbar\omega_{\pm} = \frac{\hbar\omega_c + \hbar\bar{\omega}}{2} \pm \sqrt{\frac{(\hbar\omega_c - \hbar\bar{\omega})^2 + \Omega_o^2}{4}} \quad (35) $$

The subindex $\pm$ has been added to stress that the range of frequencies has now been limited to those corresponding to those well outside of the uncoupled frequency distribution, that is, to the polaritonic frequencies. It is evident that, as long as the polaritons are separate enough from the bulk of the distribution, the polariton frequencies will only depend on the central frequency of the distribution and on the relative frequency of the cavity.

Note that Eq. (35) proves that, concerning the polaritons, as long as the previous assumptions hold, the system is reducible to a simple one: a single cavity mode of frequency $\omega_c$ and a single matter mode of frequency $\bar{\omega}$, at the average of the distribution. Indeed, as the zeros of the denominator must be the polaritonic frequencies, we see that they have the same form as those coming from the Tavis-Cummings Hamiltonian of such a simple system:

$$ \hat{H}_{TC,1+1} = \begin{pmatrix} \hbar\omega_c & \hbar\Omega_o/2 \\ \hbar\Omega_o/2 & \hbar\bar{\omega} \end{pmatrix} \quad (36) $$

We can now finally obtain the polariton coefficients, by deriving expressions equivalent to Eq.(28) and Eq.(29). By again considering the limit in which the complex part of the frequency vanishes in Eq.(34-a) and taking the imaginary part, an expression isomorphic to Eq (28) is ob-
tained:

$$ Im(\tilde{G}_{c,c}(\omega)) = -\pi \left( \frac{\hbar\bar{\omega} - \hbar\omega_{+}}{\hbar\omega_{-} - \hbar\omega_{+}} \delta(\hbar\omega - \hbar\omega_{+}) + \frac{\hbar\bar{\omega} - \hbar\omega_{-}}{\hbar\omega_{+} - \hbar\omega_{-}} \delta(\hbar\omega - \hbar\omega_{-}) \right) \quad (37-a) $$

$$ |C_{c,\pm}|^2 = \frac{\hbar\bar{\omega} - \hbar\omega_{\pm}}{\hbar\omega_{\mp} - \hbar\omega_{\pm}} \quad (37-b) $$

$$ |C_{i,\pm}|^2 = \left( \frac{g}{\hbar\omega - \hbar\omega_i} \right)^2 |C_{c,\pm}|^2 \quad (37-c) $$

Thus, analytical results for the polaritons will be compared with numerical results in the following chapters. It is worth noting that equation (36) is also found by other works dealing with the exact solution of the dynamics of N harmonic oscillators under VSC²⁵, albeit it was obtained in the context of determining the Hessian matrix in the reactant well. This suggests that the Tavis-Cummings model is a good approximation to the N-harmonic oscillator problem under VSC in the single excitation assumption.

Alternative treatments³⁹ can describe the properties of the dark states, using disorder-averaged Green's functions to obtain an effective block-diagonal Hamiltonian for the system. As we have neglected all the features of the matter but its central energy, what we will recover here is exclusively the first block of this Hamiltonian, and thus are expected to be equivalent to the literature in the polariton-related properties.

### III. COMPUTATIONAL DETAILS

#### A. Monte Carlo Simulations

We used the Monte Carlo method to sample the Hamiltonians of Equations (6) and (7). Notably, we combined it with Path Integral theory, resulting in additional degrees of freedom in Monte Carlo sampling, namely the beads. At this aim, we have considered 2 types of Monte Carlo moves:

1. A 1 bead - 1 oscillator movement: a random update of the coordinate of a single oscillator,
in a single imaginary time point (a single bead) is performed, with random orientation and magnitude. This corresponds to applying a transformation of the form:

$$ \mathbf{x}_{i,j}^{k+1} = \mathbf{x}_{i,j}^{k} + d\mathbf{x}^{k} \qquad d\mathbf{x}_{i,j}^{k} \sim U(0, 0.2 \AA) \qquad (38) $$

where i stands for the oscillator index and j for the bead index

2. A P beads - 1 oscillator movement: the same random displacement, in magnitude and orientation, is applied to all the P beads corresponding to the same oscillator.

$$ \{\mathbf{x}_{i,j}^{k+1}\}_{j=1}^{P} = \{\mathbf{x}_{i,j}^{k} + d\mathbf{x}^{k}\}_{j=1}^{P} \qquad d\mathbf{x}_{i,j}^{k} \sim U(0, 0.2 \AA) \qquad (39) $$

where again i stands for the oscillator index and j for the bead index

The relative frequency of these Monte Carlo moves has been adjusted in order to obtain a total ratio of accepted/rejected configurations of 70%. This was achieved by applying Type 2 movements to generate 20% of the configurations.

The configurations proposed as a result of both trial movements have been evaluated through their action, as defined in Eq. (14) for the Pauli-Fierz Hamiltonian. Analogously, the action for the linearized Hamiltonian has been calculated by neglecting the quadratic term. The acceptance probability for a given configuration is:

$$ P(x_k) = \begin{cases} 1 & \mathcal{A}_k < \mathcal{A}_{k-1} \\ \exp(-(\mathcal{A}_k - \mathcal{A}_{k-1})) & \mathcal{A}_k \ge \mathcal{A}_{k-1} \end{cases} \qquad (40) $$

as in standard Monte Carlo algorithms.

Simulations were done using different parameters, reported in Table I. PIMC were run using 64 beads (P = 64) while classical MC is recovered using P = 1. For other parameters (like the number of matter modes, N) we used a range of values to better study how PIMC results evolve as a function of them.

In general, multiple PIMC runs are performed and the corresponding force-force and displacement-displacement correlation functions are calculated for each run, averaged and then used to solve
<table><caption>TABLE I. Summary of simulation parameters used in PIMC simulations.</caption><thead><tr><td>Parameter</td><td>Values</td><td>Description</td></tr></thead><tbody><tr><td>P</td><td>[1,64]</td><td>N of beads, for classical and quantum runs.</td></tr><tr><td>Ncyc</td><td>10<sup>8</sup></td><td>N of Monte Carlo trials</td></tr><tr><td>N</td><td>[1,225]</td><td>N of matter modes</td></tr><tr><td>mass</td><td>0.9411au</td><td>Reduced mass of oscillators.</td></tr><tr><td>T</td><td>300K</td><td>Temperature</td></tr><tr><td>D<sub>e</sub></td><td>500kJ/mol</td><td>Dissociation energy of Morse potential</td></tr><tr><td>Step</td><td>0.2Å</td><td>Max. size of a single MC step.</td></tr><tr><td>Estimator interval</td><td>10</td><td>Sample estimators every 10 configs</td></tr><tr><td>Start sample</td><td>0.3</td><td>First 30% of configs. to equilibrate</td></tr><tr><td>Coll. freq</td><td>0.2</td><td>Frequency of Type 2 MC movements</td></tr><tr><td>N<sub>runs</sub></td><td>[25,400]</td><td>Number of MC replicas averaged</td></tr></tbody></table>

the generalized eigenvalue problem of Eq (23). From these eigenvalues we obtain frequencies of the spectra. The vibrational density of states (VDOS) are calculated as a histogram of such frequencies, while infrared (IR) signals are obtained as a histogram of the frequencies, weighted by the sum of coefficients of the correspondings eigenvectors. Then, both VDOS and IR spectra have been convoluted with a Lorentzian function to have the typical vibrational lineshapes.

## B. Error calculation in PIMC

The estimation of the errors from PIMC simulations has been done through a bootstrapping method, as described in the following. Our typical PIMC simulation consists of 200 parallel PIMC runs, that differ only in their random seed. This produces 200 displacement-displacement correlation matrices used then in the eigenvalue equation (23). To improve the accuracy of our calculations, the matrices are averaged prior to diagonalization. To estimate the error we have estimated the distribution of each matrix elements, assuming normal and independent distributions. We then constructed a sample of matrices, which are then averaged in blocks of 200 and diagonalized. This generates a distribution for any of the observables from which we can evaluate the standard deviation.
### C. Diagonalization of Pauli-Fierz Hamiltonian

For testing the PIMC method proposed here, we have compared results between the direct diagonalization of the PF Hamiltonian Eq.(1)-(3) and Eq. (7) and the PIMC results, for small number of matter modes, up to 4 matter modes. At this end we have expressed the wavefunction on a finite basis, product of harmonic oscillator eigenstates, similarly to what done previously by us.⁴⁶. The Morse potential was applied using a discrete variable representation with a Gaussian quadrature. For the kinetic energy and harmonic terms of the PF Hamiltonian, Eq.(1), exact relations of Hermite polynomials were used.⁴⁶

The basis set size was 6 basis set per mode, including the cavity, so a total of $6^5 = 7776$ basis set elements. For the cavity, the harmonic oscillator parameters for setting the basis set were chosen equal to the free cavity parameters. For matter modes, the basis set parameters were optimized to minimize the isolated Morse oscillator ground state energy. These optimized parameters are a displacement of $x_0 = 0.02$ for the center of the basis set and an effective pulsation $\omega_{basis} = 3366 \text{ cm}^{-1}$. The resulting first transition for the isolated Morse potential is $\omega_{01} = 3262.4 \text{ cm}^{-1}$, while the theoretical value is $3261.4 \text{ cm}^{-1}$.

The intensity of each $0 \rightarrow n$ transition was computed from the transition dipole moment as

$$I_n = \frac{m}{2N\hbar^2} |\langle 0 | \sum_i \mu_i | n \rangle|^2 (\varepsilon_n - \varepsilon_0), \quad (41)$$

where N is the number of matter modes, m their mass, and $\varepsilon_n$ the eigenstate energies. The normalization $\frac{1}{2Nm\hbar^2}$, is introduced so that the intensities should satisfy the sum rule$^{47} \sum_n I_n = 1$, which was numerically verified with precision $10^{-5}$. It allows for a direct comparison with the PIMD results without any scaling parameter. The individual dipole moment operators were applied using the exact relation of Hermite polynomials.
# IV. RESULTS

## A. Properties without disorder

Before studying in detail how Vibrational Strong Coupling modifies the eigen-energies and states of the model Hamiltonians, we shall show how Monte Carlo is able to provide such information and what is the role of nuclear quantum effects included via Path Integral theory.

To this end we have considered a simple set of 4 Morse oscillators with $\omega_0 = 3404 \text{ cm}^{-1}$ and ran $10^8$ Monte Carlo steps, which are enough to have converged values. Using Equation (23) we can obtain the transition energies and corresponding effective states, and thus we can calculate the electric dipole and therefore the IR intensity. Results obtained from classical MC (cl-MC) and PIMC are shown in Figure 1.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Density of states (full lines) and corresponding IR spectra (dashed lines) obtained from cl-MC (in red) and PIMC (in blue) simulations using the Morse model uncoupled Hamiltonian.</figcaption></figure>

Notably, the cl-MC results show a red shift of 34 cm⁻¹ (the peak is found now at 3370 cm⁻¹) with respect to the generating frequency of 3404 cm⁻¹ while PIMC simulations provide a larger red shift of 109 cm⁻¹ (the peak is now at about 3295 cm⁻¹).
These frequencies correspond to vibrational transitions at finite temperature and we can compare them with the pure eigen-energies of a simple Morse oscillator:

$$E_n = \hbar\omega_0\left(n + \frac{1}{2}\right) - \frac{[\hbar\omega_0(n + 1/2)]^2}{4D_e} \quad (42)$$

where $\omega_0$ is the generating frequency, $D_e$ is the Morse parameter, $n$ is the vibrational quanutm number and $\hbar$ is the reduced Plank constant. Using our Morse parameters, we obtain a red shift of 17 cm⁻¹ for $n = 0$ and 156 cm⁻¹. The $0 \to 1$ transition is then red-shifted by 139 cm⁻¹ which is similar to what obtained in PIMC calculations (about 114 cm⁻¹). PIMC with the approach of Morresi et al.³⁶ recalled in section II C is thus able to capture the correct vibrational frequency.

We then run PIMC simulations using the Pauli-Fierz Hamiltonian (Equation (7)), adjusting $\gamma$ according to Equation (9), to obtain a Rabi Splitting of 750 cm⁻¹ which roughly corresponds to what is found in the case of resonance with O-H stretching in experiments.⁹ The cavity frequency, $\omega_c$, is set to 3290 cm⁻¹ which corresponds to the infrared band of the system according to PIMC. Density of states (DOS) and IR intensities are reported in Figure 2-A for PIMC together with results obtained from exact diagonalization of a system of 4 Morse oscillators using the Pauli-Fierz Hamiltonian. In PIMC simulations the two polaritons appear at 3679 cm⁻¹ (P⁺) and 2928 cm⁻¹ (P⁻), corresponding to a Rabi splitting of 751 cm⁻¹. These results are very close to the exact diagonalization results, in which the polaritons are located at 3653 cm⁻¹ (P⁺) and 2912 cm⁻¹, respectively.

In the same Figure 2 (panel B) we report DOS and IR intensities obtained using the same number of modes and $\gamma$, as in PIMC calculations, but then we performed the diagonalization of an equivalent Tavis-Cummings Hamiltonian (Equation (24))³⁷,⁴⁰,⁴³ For further comparison, the infrared spectrum obtained from the Green Function analysis (Section II D) of this Hamiltonian is presented in both panels. We notice that for this Tavis-Cummings Hamiltonian both lower and upper polaritons are shifted by exactly the same amount, 375 cm⁻¹, which is half of the Rabi splitting, while in the case of PIMC and exact diagonalization, an asymmetry is observed.

Furthermore, we obtained in both cases dark states at the frequency corresponding to the uncoupled system which are populated by $N-1$ modes. We will discuss in more details the contribution of matter and cavity (light) in the resulting peaks when we consider disorder. Here, given the
<figure><img src="image_3.png" /><figcaption>FIG. 2. Density of states (in full lines) and corresponding IR intensities (in dashed lines) as obtained from a system composed by 4 identical matter modes and 1 cavity mode with Ω₀ =750cm⁻¹. Panel A: PIMC simulations using Morse matter modes and Pauli-Fierz Hamiltonian (Equation (5)); Panel B: diagonalization of the Tavis-Cummings Hamiltonian (Equation (24)). In panel A, we report in orange the polaritonic bands from the diagonalization of the Pauli-Fierz Hamiltonian (Section III C). In both panels we also report the Green function results for the polaritons (in dashed green lines).</figcaption></figure>
simplicity of the model (all matter modes have the same frequency) the Tavis-Cummings Hamil-
tonian provides a dark state which is composed only by matter modes and polaritons which are
50% populated by the cavity and 50% by matter, these lasts with equal coefficients. However, we
note that in Figure 2-A the upper polariton is more intense than the lower polariton. The intensities
of the two polaritons are there too well reproduced by the PIMC calculations. This asymmetry in
intensity, together with the asymmetry in frequency mentioned above, will be investigated in more
details in section IV C.

Overall, we note that both VDOS and IR obtained from PIMC simulations are in excellent
agreement with what obtained by quantum Hamiltonians.

Now, we investigate how the vibrational transitions are modified when the cavity frequency is
detuned from exact resonance. Notably, we have displaced it by up to 1000 cm⁻¹ lower and higher
from the exact resonance using PIMC with 25 Morse matter modes. Results are shown in Figure 3
where we report both VDOS (panel A) and IR intensities (panel B).

We observe a behavior which is typical of an avoided crossing in which one of the two polari-
tons decays as a function of ω_c and the other gathers all the intensity until reverting to being akin
to the uncoupled vibrational band: bright and fully composed of matter. The rate at which the
polariton decays as the cavity is detuned is controlled by the parameter γ. Notably, the vibrational
bands are coupled on a long range of frequencies: we can observe an effect on the position almost
up to 1000 cm⁻¹ away from resonance, and in the IR intensity up to a about 250 cm⁻¹. Also in
this case, results obtained from PIMC simulations are in excellent agreement with what obtained
from analytical Green function theory.

Summarizing, the PIMC results in absence of disorder are in excellent agreement with theo-
retical analysis and exact Hamiltonian diagonalizaton. We can now consider the disorder in the
uncoupled matter modes and study how different properties will be (eventually) affected.

B. Polaritons with disorder in the matter

We now move to a model in which the matter modes are distributed around a central value in
order to mimic the inhomogeneous broadening on vibrational transitions occurring in condensed
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Density of states (DOS, panel A) and IR intensity (panel B) as obtained from PIMC simulations with 25 Morse matter modes while detuning the cavity frequency. The matter frequencies are reported as a red dashed line, the cavity frequencies as a red dashed line. In the same panels, we show the Green function results (in yellow squares).</figcaption></figure>
phase. We consider the O-H stretching inhomogeneous broadening using a distribution of fre-
quencies with a width at half maximum (FWHM) of 200 cm⁻¹,⁹ that we generate using a normal
distribution of initial matter frequencies (i.e. frequencies without the VSC coupling). In Figure 4-
A we show for later reference how the vibrational density of states (VDOS) and IR spectrum are
recovered from PIMC simulations when including disorder in the matter but for an uncoupled
Hamiltonian (γ = 0).

We now introduce VSC, where the cavity is in resonance with the maximum of the distribution
and we consider the a regime where Ω₀ > σ. In other words, the width of the initial distribution is
smaller than the Rabi splitting, but of the same order of magnitude.

VDOS and IR spectra under VSC are reported in Figure 4-B. Firstly, we notice that dark states
stay mostly dark also in case of disorder. Secondly the Rabi splitting is the same as in the corre-
sponding system without disorder (i.e. using the same coupling strength and with all the modes
with the same frequency as the central frequency in the disorder case). Notably, we have a Rabi
splitting of 755 cm⁻¹ with two polaritonic bands in the IR at 2921 and 3676 cm⁻¹ for P⁻ and P⁺,
respectively.

We now want to check if Eq. (9) also holds when disorder is introduced. At this end we run
simulations with different values of N = 4, 9, 16, 25, 36, 49, 64 and 81, while keeping the same
coupling constant γ which corresponds to the Rabi splitting of 750 cm⁻¹ for N = 1. In Figure 5
we plot the ratio between the N=1 Rabi splitting Ω, and the value obtained for each different N,
Ω₀. As clearly shown, we obtained the expected linear dependence. This strengthens the picture
that our Monte Carlo simulations are able to provide the correct properties of the Pauli-Fierz
Hamiltonian and that, if we limit to the simple Rabi splitting observation, the disorder (modeled
with a simple frequency distribution) has not a relevant effect on the resulting frequencies as
the same behaviour is obtained for simulations without disorder. As predicted by Eq. (35) our
numerical calculation shows that the Rabi splitting is independent of disorder. To understand this,
as in the case of disorder only a fraction of matter states are in resonance with the cavity, we will
study hybridization of light and matter modes. Nevertheless, if the disorder grows to the point that
it becomes comparable to the Rabi splitting, the homogeneous regime condition is broken, and
our Green function theory does no longer apply. In this situation, strong coupling is lost, and the
polaritonic bands will no longer be present in the spectrum³⁹.
<figure><img src="image_6.png" /><figcaption>FIG. 4. Density of states (blue solid lines) and IR intensities (red dashed lines) as obtained from PIMC calculations with 100 matter modes, P=64 path integral beads. Panel A: uncoupled system, panel B: coupled with a cavity with $\omega_c=3290 \text{ cm}^{-1}$ (grey vertical dashed line) with a target Rabi splitting of $\Omega_o=750 \text{ cm}^{-1}$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 5. Ratio between the collective Rabi Splitting ($\Omega_o$) and the predefined single molecule splitting ($\Omega$) as a function of the square root of the number of molecules ($N$) as obtained from PIMC simulations with the Pauli-Fierz Hamiltonian. In blue dashed line we show a linear fit of simulation data.</figcaption></figure>

### C. Light-matter hybridization

In Section IV B we have shown how PIMC is able to provide both Vibrational Density of States and in the IR intensities in the case of a model describing disordered matter under VSC. We now focus on the composition of these states after coupling in terms of the original modes before the coupling is applied. Notably, the effect of hybridization involves the suppression of the IR activity of the most part of states, without a large effect in their energies.37,43 The analysis of the effective eigenstates obtained as described in Section II C, as well as the diagonalization of the Tavis-Cummings Hamiltonian, allows to decompose the final states in terms of the original modes. We will focus on both the cavity contribution to bright and dark states, as well as on the delocalization of the coupled states over the matter modes.
# 1. Analysis of the Polaritons

We first look at the two spectral characteristics of the polaritons that we can access: their energies and their IR intensities as a function of the coupling, expressed as $\Omega_o$. It can be seen, from Figure 2A and Figure 4, that the spectra are asymmetric, with a more intense upper polariton. Also, the polaritonic splitting appears not to be exactly centered on the cavity frequency. On the contrary, the results reported from the Tavis-Cummings Hamiltonian diagonalization in Figure 2B, show that the intensities of both polaritons are equal and their frequencies centered on the cavity mode. This suggests that the asymmetry depends on the choice of Hamiltonian.

In order to describe this asymmetries, we consider two quantities: (i) the difference in intensity, $I_{P+} - I_{P-}$; (ii) the asymmetry in frequency, obtained from:

$$f = \frac{\omega_{P+} + \omega_{P-}}{2} - \omega_o \quad (43)$$

where $\omega_{P+}$ and $\omega_{P-}$ are the frequencies of the upper and lower polariton, respectively, and $\omega_o$ is the average frequency of the matter without the coupling.

Asymmetries in frequency as a function of Rabi splitting using the different Hamiltonians are reported in Figure 6-A, while asymmetries in intensity are reported in Figure 6-B. For both Pauli-Fierz and linear Hamiltonian we observe that the asymmetry grows with the Rabi Splitting but with substantial differences. Furthermore, from the diagonalization of the Tavis-Cummings Hamiltonian we do not observe any asymmetry (as already noticed from the spectrum reported in Figure 2B). We must notice that, when increasing the Rabi Splitting, we move from the VSC to the ultra strong vibrational coupling regime (V-USC)$^{48}$. This regime is characterized by the condition $\frac{\Omega_o}{2\omega_o} > 0.1$, which in the present case corresponds to $\Omega_o = 660 \text{ cm}^{-1}$. Interestingly, we do not see a striking difference for Rabi splitting values larger than this threshold (at least up to $1500 \text{ cm}^{-1}$). Our analysis done for $\Omega_o = 750 \text{ cm}^{-1}$ is slightly in the V-USC regime but results are qualitatively the same in the VSC regime. We now analyze in more details these differences by considering a simple 2-mode model, in which one collective matter mode couples to the mode of the cavity. We construct this model starting from the Pauli-Fierz and Tavis-Cummings Hamiltonians, transforming the basis of isolated molecular modes into another one in which a single collective mode has a dipole and thus can couple to the cavity while the rest remains dark.$^{39,49-52}$ In this new basis, the
<figure><img src="image_8.png" /><figcaption>FIG. 6. Panel A: frequency asymmetries (as defined in Eq.(43)) as a function of the collective Rabi splittings. Panel B: intensity asymmetry, as difference between upper and lower polaritons. Red dots: PIMC with Pauli-Fierz Hamiltonian simulations; blue dots: PIMC with linear Hamiltonian; red crosses: 2-mode model from equation (45); blue crosses: 2-mode model from equation (44); black diamonds: exact diagonalization of the Pauli Fierz Hamiltonian for 4 matter modes as detailed in Section III C. As dashed lines we report the results obtained from the diagonalization of Tavis-Cummings Hamiltonian. The vertical grey dotted line marks the transition between VSC and V-USC.</figcaption></figure>
spectroscopy of the system can be explained by coupling a single bright mode to the cavity. For the Tavis-Cummings Hamilmtonian we obtain the following 2 × 2 matrix:

$$ \hat{H}_{2mode} = \begin{pmatrix} \omega_c^2 & \omega_c \Omega_o \\ \Omega_o \omega_c & \omega_o^2 \end{pmatrix} \qquad (44) $$

where $\omega_o$ is the average frequency of matter.

In the case of the Pauli- Fierz Hamiltonian, the asymmetries are positive, that is, $|P+\rangle$ reaches higher energies and intensities than it would, would the spectrum be symmetric. This arises from the self interaction energy of the matter in the dipole gauge: it amounts to a harmonic contribution to the bare mode energies, a "blue shift". In the two coupled oscillators model, this translates into:

$$ \hat{H}_{2mode} = \begin{pmatrix} \omega_c^2 & \omega_c \Omega_o \\ \Omega_o \omega_c & \omega_o^2 + \Omega_o^2 \end{pmatrix} \qquad (45) $$

Therefore, the coupling modifies the bare energies of the modes, and introduces a positive detuning between the matter distribution and the cavity mode. This can be reproduced by simply adding a $\Omega_o^2$ term to the energy of the bright collective mode of matter. This asymmetry also has an effect on the intensities, since the most intense polariton will be the closest one in frequency to the matter. The $\hat{H}_{TC}$ Hamiltonian does not show any of this asymmetries, as is expected. The linear Hamiltonian also produces asymmetries with the coupling strength, but different: the lower polariton becomes overly low in energy, while the intensity is barely affected. This can be explained through the asymmetric character of the polaritons in a resonant 2-mode model. In such a model, the polaritonic frequencies can be written as:

$$ \omega_{P\pm} = \omega \left( 1 \pm \frac{\Omega_o}{\omega} \right)^{1/2} \qquad (46) $$

where $\omega$ refers to the frequency at which matter and cavity are resonant. It is evident in (46) that the eigen frequencies from the 2 mode naturally lead to an asymmetric splitting with respect to the coupling.
We can further verify that the asymmetries obtained from PIMC are in very good agreement with the exact diagonalization of a reduced number of modes under the Pauli-Fierz Hamiltonian, also reported in Figure 6.

We now turn to the study of the composition of the polaritons in terms of the bare molecular modes, as we report in Figure 7.

<figure><img src="image_9.png" /><figcaption>FIG. 7. Panel A: contributions to $|P-\rangle$ of matter modes ($N = 100$) as a function of distance from the polariton peak. Panel B: inverse of the coefficients of the matter mode as a function of distance from the polariton peak. Red dots: $PIMC - \hat{H}_{PF}$ simulations; blue dots: $PIMC - \hat{H}_{lin}$, dashed black line: Tavis-Cummings Hamiltonian diagonalized; solid green line: Green's function analytical theory.</figcaption></figure>

We obtain the contributions ($C_i^2$) of every matter mode hybridizing in the low polariton through the three Hamiltonians that we discuss, as well as through the analytic study the Green's function of $\hat{H}_{TC}$ as detailed in Section II D. Note that, since the systems studied using PIMC are composed by anharmonic vibrational modes, they will have a red shift with respect to their hamonic frequency $\omega_0$. Moreover, as we have seen in Figure 6, there will be asymmetries in the frequencies, resulting from the form of the Hamiltonian. In order to make different types of simulation comparable, we have represented in the horizontal axis of Figure 7 the difference in frequency between the lower polariton and the frequency of each mode, correcting the red shift of the Morse potentials and the "blue shift" of the Pauli-Fierz Hamiltonian:

$$ \omega_{eff} = \omega_{mode} - \Delta\omega_{anh} + \Delta\omega_{PF} \quad (47) $$

where $\omega_{eff}$ is the effective frequency of a given matter mode, $\Delta\omega_{anh}$ is the frequency red-shift
of the vibrational frequency induced by the Morse potential captured by PIMC, and $\Delta\omega_{PF}$ is the higher frequency shift due to the self interaction term in $\hat{H}_{PF}$, which has been already presented in equation (45). Note that, in the case of the linear Hamiltonian, $\Delta\omega_{PF} = 0$, and in the case of the Tavis-Cummings Hamiltonian and the green function analysis, $\Delta\omega_{anh} = \Delta\omega_{PF} = 0$.

In this way, we obtain a good overlap of the data from different simulation types, proving that the only difference between them lies precisely in anharmonicity and the self interaction from $\hat{H}_{PF}$ that renormalize the matter mode frequencies.

We observe that the matter mode contributing the most to the polariton is not the one in resonance to the cavity, but rather the one closest in energy to $|P-\rangle$. Moreover, by representing the inverse coefficients in panel B, we recover the relation between the mode coefficients in the bright states and the mode frequencies in Eq.(34-b). This observation provides an explanation to the fact that disorder does not play a significant role on the value of the Rabi splitting, as the mode contribution is not a function of the distance to the cavity frequency, but of the distance to the polariton frequency.

## 2. Dark states composition

While much attention has been devoted to the "bright states", i.e. the polaritons,$^{30,31,39}$ the role of the dark states in VSC has been explored to the same extent.$^{37}$ Many of the experimental attempts at VSC chemistry report a decrease in reaction rate, which can be explained by the "quenching" of the reactive vibration mode of the reagent, through some energy redistribution mechanism.$^{22,53,54}$ The dark states would be able to open additional de-excitation pathways if they are delocalized over enough vibrational modes of the system, that is, if multiple different modes participate on them. For this purpose, we calculate the Inverse Participation Ratio (IPR) to quantify delocalization on the dark states (Eq. (48)):

$$IPR_i = \sum_{j=1}^{N_{\text{matter}}} |c_j^{(i)}|^4 \qquad (48)$$

where $c_j^{(i)}$ is the coefficient of the mode $j$ to the final state $i$. Note that low values of the IPR indicate high delocalization (the lowest possible value is $1/N$, where $N$ is the number of modes in
the state) while IPR = 1 reveals a fully localized state.

In Figure 8 we report the IPR values as a function of the number of Monte Carlo runs and compare it with what is obtained from the numerical diagonalization of the Tavis-Cummings Hamiltonian. In particular we focus on the most delocalized dark state and on the average value. In both cases, we need a large number of MC runs to converge to what is obtained from Tavis-Cummings Hamiltonian diagonalization. This shows that the numerical noise intrinsic in MC procedure can be critical for evaluating IPR.

<figure><img src="image_10.png" /><figcaption>FIG. 8. Convergence with the number of independent Monte Carlo runs of the IPR at 2 key points: the average resonant dark state and the most delocalized dark states. In blue: minimum dark state IPR. In red: average dark state IPR at the cavity frequency. In solid lines: IPRs from PIMC with Pauli-Fierz simulations. In dashed lines the IPRs associated with the same distribution of bare modes, calculated by exact diagonalization of the Tavis-Cummings Hamiltonian.</figcaption></figure>

At convergence, we have calculated the IPR for each hybrid mode, while using the Gaussian distribution and also when using a uniform distribution of modes with a width of 200 cm⁻¹. Results obtained from PIMC and Tavis-Cummings diagonalization are reported in Figure 9. The polaritonic states are also reported for comparison, and, as expected, they have an almost maximal degree of delocalization, which corresponds to a situation in which all molecules are in resonance with the cavity. Even if most modes are not at perfect resonance, the polaritons are still highly delocalized. Notably, the dark states are also delocalized, in particular those closer to the center
of the distribution. Figure 9-A displays the IPR distribution that stems from a set of 100 equally spaced matter modes over a range of 200 cm⁻¹, next to the equivalent values from the diagonalization of the Tavis-Cummings Hamiltonian. It is clearly visible that the noise inherent to the Monte Carlo procedure generates artificially low IPR values. Nevertheless, through a parabolic fit, the average IPR at every energy can be recovered and it is in good agreement with the Tavis-Cummings Hamiltonian as well as with the results from previous work. The situation in Figure 9-B is fundamentally different, as the matter vibrations are now normally distributed, aiming to achieve a FWHM of 200 cm⁻¹, which makes this simulation comparable with the previous one and tries to imitate the OH stretching of water. In this panel, the dark states at the resonant frequency can reach IPR values as low as 0.1, which translate in 10 modes under hybridization in both the TC Hamiltonian and PIMC, 10% of the ensemble. This is due to the higher density of the Gaussian in the central region, and it also explains why the PIMC and the Tavis Cummings results are similar: as the normal distribution is already random, it is less affected by the stochastic behavior of PIMC. We also see that, while the minimum IPR is lower in the normal distribution case, the average value at each energy is close to its equivalent in the uniform distribution case. This highlights the fact that, while average IPRs have been previously reported³⁷,³⁸,⁴³ in the ranges corresponding to delocalization over mainly 3 modes, and they are indeed consistent among different types of distribution, this overlooks the fact that there are states that will be more extensively delocalized, and that they correspond to the resonant frequency.

The IPR analysis gives a global picture. To study more in details the delocalization of the dark states, we now consider how the coefficients are distributed for the most delocalized states (i.e. the frequencies which are closer to the cavity frequency, $\omega_c$). Results are shown in Figure 10. As we can notice, there are two modes which mostly contribute to this final state and they are those which initial frequencies ($\tilde{\omega}_{D_0}$) are closest to the frequency of the final state ($\omega_{mode}$). This behavior is independent on the model Hamiltonian and the way the coefficients are obtained and in agreement with what has been found by Pupillo et al.³⁷,⁴³ while analyzing the simpler Tavis-Cummings Hamiltonian⁴⁰.

As can be seen in Figure 10 the IPR values reported in Figure 9 correspond to delocalizations mainly over three or four modes closest in energy to the dark state. Similar dark state delocalizations have been found theoretically to increase electron transfer rates up to 50%, and show that the dark states can also have an impact in polaritonic chemistry that is often overlooked.³⁸ Finally,
<figure><img src="image_11.png" /><figcaption>FIG. 9. IPR values of the eigenstates of the model as a function of the displacement from the resonant frequency ($\omega_o$). Panel A: uniform matter distribution (200 cm$^{-1}$ width); Panel B: Gaussian distribution of matter modes (FWHM = 200 cm$^{-1}$). Red dots: Tavis-Cummings Hamiltonian diagonalization; blue dots: PIMC simulations with Pauli-Fierz Hamiltonian. The solid black line shows a quadratic fit of the dark states IPRs. The dashed horizontal grey line is set to the minimum (theoretical) IPR value corresponding to a fully delocalized state.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. Squares of mode coefficients contributing to the most delocalized dark state ($\omega_{D_0}$). Panel A: uniform matter distribution (200 cm$^{-1}$ width); Panel B: Gaussian distribution of matter modes (FWHM = 200 cm$^{-1}$). Red dots: Tavis-Cummings Hamiltonian exact diagonalization; blue dots: PIMC simulations with Pauli-Fierz Hamiltonian; dashed grey line: mode coefficients in a fully delocalized state.</figcaption></figure>

we should noticed that the frequency distribution has a certain impact on the IPR distribution, and in particular on the lenght of delocalization, as visible from comparing gaussian and uniform distributions in Figures 9 and 10.
## D. Large number of molecules limit

In this work, we have introduced the quantum simulation of a cavity coupled to a number of identical replicas of the molecular system. We have presented results with up to 100 replicas however, in a realistic setting, with a cavity of about 1 µL,⁵⁵ the number of molecules coupled to the cavity is $6.02 \times 10^{17}$ molecules for a concentration of 1 mol.L⁻¹. For a Rabi splitting $\Omega_0$ of a few 100's of wavenumber, say $\Omega_0 = 750$ cm⁻¹ as in Figure 4, the corresponding coupling constant $g_0$ for a single molecule is then $g_0 = \frac{\Omega_0}{\sqrt{N}} = 10^{-6}$ cm⁻¹. It appears then than even in such small cavities, one should consider the limit of a large number of molecules with $g_0 \to 0$ at $\Omega_0 = g_0\sqrt{N}$ constant.

In this limit, we can apply the quantum central limit theorem⁵⁶. It is then expected that the statistics of an extensive operator like the total dipole moment, properly normalized, i.e.:

$$ \mu_N = \frac{1}{\sqrt{N}} \sum_{i}^{N} \mu_i \qquad (49) $$

where $i$ runs over the $N$ matter modes and $\mu_i$ are the individual dipoles of each oscillator, becomes Gaussian and analogous to what is expected for a harmonic oscillator, even though each individual oscillator is not harmonic.

In Figure 11 we show the collective dipole distribution as obtained from 1 and 100 oscillators, as well as the distribution of the individual dipole for one single mode belonging to the ensemble of $N = 100$. The occurrence frequency is reported in logarithmic scale to better appreciate the parabolic behavior. We observe that when using only one matter mode, the dipole distribution is not symmetric, while it becomes symmetric and Gaussian when more modes are considered (here for simplicity we report only the case of $N = 100$, but similar results are obtained for larger $N$ values). The variance of the dipole distribution is constant as a function of $N$. Nevertheless, the distribution of the dipole of a single molecule remains asymmetric even in the $N = 100$ ensemble. In other terms, while the ensemble distribution becomes gaussian as the number of modes increases, each individual mode shows its own anharmonicity. Figure 11 and the central limit theorem suggest that at $N = 100$ convergence should already be attained, not only for the dipole distribution but also the polariton asymmetry, the IPR or the contribution of matter oscillators to the polaritons. However, it should be noted that the individual molecule coupling constant, $g_0$,
remains artificially large.

The quantum central limit theorem also justifies the use of perturbation theory, which we apply here by considering only the first excited state. Indeed, when the cavity is put inside a spectrometer for measuring transmittance, the laser power is such that many photons are present in the cavity, not just one. Usually, extensivity is more or less explicitly employed to justify that only a microscopic portion of the system is considered, in which on average the number of photons adsorbed is much lower than one (in the linear regime). Here, when we consider the whole of the cavity, this argument can not be used. However, when the quantum central limit theorem can be applied, the linear response theory becomes exact$^{57}$. This can be understood by the fact that matter inside the cavity behaves as a harmonic oscillator in this limit, hence the Gaussian distribution of the normalized dipole moment, and together with the harmonic structure of the electromagnetic field, the situation is analogous to coupled harmonic oscillators and the transition between the ground state and first excited state is sufficient to describe all transitions in the system.

In terms of extensivity and intensity, note that the Rabi splitting, being proportional to $\gamma\sqrt{N} \propto \sqrt{\frac{N}{V_S}}$ is intensive for a fixed concentration of the chromophores inside the cavity. Note again that while our results appear converged for $N = 100$ oscillators, the corresponding cavity volume $V_S$ is much smaller than in realistic settings, consistent with a large individual molecule coupling constant, $g_0$.

V. CONCLUSIONS

We have performed PIMC calculations of an ensemble of replicas of a molecular vibration, modeled by Morse oscillators, coupled to a Fabry-Perot cavity using the Pauli-Fierz Hamiltonian. The large number of replicas used, up to a hundred, allowed to simulate the effect of disorder among the vibrational modes, while PIMC does not make approximations on the system Hamiltonian. This number of replicas is much lower than the typical number of chromophores in a realistic settings, about $10^8$, yet we show that the collective behaviour is correctly captured, with a $\sqrt{N}$ dependence of the Rabi splitting. The vibrational spectrum was obtained from the PIMC calculations from the Kubo transformed auto-correlation of the mode displacements. It also provides a decomposition of the vibrational modes of the coupled system in terms of the original uncoupled
<figure><img src="image_13.png" /><figcaption>FIG. 11. Dipole moment histogram of matter in PIMC simulations under VSC. In blue crosses: dipole histogram of a simulation containing 100 matter modes; in red dots: dipole histogram of a simulation containing a single matter mode; in black crosses, the dipole histogram of one out of the 100 modes represented in blue crosses.</figcaption></figure>

modes, allowing for interpreting the features of the vibrational spectra.

These results were compared with semi-analytical and analytical results based on the Tavis-Cummings Hamiltonian, and with numerical diagonalization of the Pauli-Fierz Hamiltonian for four vibrational modes. The latter showed the great accuracy of the PIMC simulations that can reproduce quantitatively the excitation frequencies and intensities, including the asymmetry of the polaritonic states that arise with the Pauli-Fierz Hamiltonian and that was studied in details. In agreement with semi-analytical or analytical theories based on the Tavis-Cummings Hamiltonian it was found that disorder does not prevent the polaritonic states, in the range of disorder studied here. The composition of these polaritonic states in presence of disorder extracted from PIMC also agree with those approaches, as long as the vibrational frequencies of the uncoupled modes are properly renormalized. Furthermore, when considering the dark states, the mode decomposition applied to the PIMC simulation recovers that coupling to the cavity mode together with disorder induce a delocalization of these states over a finite number of molecules, typically three
molecules although with significant fluctuations, even for large number of molecular modes. Finally, notwithstanding the anharmonicity of each single mode, the collective mode that forms the polaritonic states behaves as a harmonic oscillator, with a Gaussian distribution of the total dipole moment, following the central limit theorem applied to this quantum system.

Altogether, augmenting path-integral simulations of matter coupled to a resonant mode of a Fabry-Perot cavity with an ensemble of replicas of the matter offers a versatile settings to simulate VSC. It can readily be applied to more realistic simulations, like the recent study VSC of liquid water,$^{23,34}$ employing Ring Polymer Molecular Dynamics (RPMD)$^{58}$. Mode decomposition applied to PIMC or RPMD simulations further provides a picture of mode delocalization induced by the coupling with the cavity. Including multiple resonant modes of the Fabry-Perot cavity (which are known to provide a different physics$^{59}$) can be readily done in this framework as well, by considering not one but several electromagnetic modes as harmonic oscillators.

## ACKNOWLEDGMENTS

This work was supported by the French National Research Agency ANR (ANR-20-CE29-0008). We are thankful to Damien Laage and Yong Chen for insightful discussions about this work.

## DATA AVAILABILITY STATEMENT

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

## REFERENCES

¹K. Hirai, J. A. Hutchinson, and H. Uji-i, Chem. Rev. (2023), 10.11021/acs.chemrev.2c00748.

²J. Feist and F. J. Garcia-Vidal, Phys. Rev. Lett. **114**, 196402 (2015).

³J. Schachenmayer, C. Genes, E. Tignone, and G. Pupillo, Phys. Rev. Lett. **114**, 196403 (2015).
⁴X. Zhong, T. Chervy, L. Zhang, A. Thomas, J. George, C. Genet, J. A. Hutchison, and T. W. Ebbesen, Angew. Chem. Int. Ed. **56**, 9034 (2017).

⁵A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, J. A. Hutchison, and T. W. Ebbesen, Angew. Chem. Int. Ed. **55**, 11462 (2016).

⁶A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, and T. W. Ebbesen, Science **363**, 615 (2019).

⁷R. M. A. Vergauwe, A. Thomas, K. Nagarajan, A. Shalabney, J. George, T. Chervy, M. Seidel, E. Devaux, V. Torbeev, and T. W. Ebbesen, Angew. Chem. Int. Ed. **58**, 15324 (2019).

⁸J. Lather and J. George, J. Phys. Chem. Lett. **12**, 379 (2021).

⁹T. Fukushima, S. Yoshimitsu, and K. Murakoshi, J. Am. Chem. Soc. **144**, 12177 (2022).

¹⁰J. George, T. Chervy, A. Shalabney, E. Devaux, H. Hiura, C. Genet, and T. W. Ebbesen, Phys. Rev. Lett. **117**, 153601 (2016).

¹¹K. Nagarajan, A. Thomas, and T. W. Ebbesen, J. Am. Chem. Soc. **143**, 16877 (2021).

¹²Y. Pang, A. Thomas, K. Nagarajan, R. M. A. Vergauwe, K. Joseph, B. Patrahau, K. Wang, C. Genet, and T. W. Ebbesen, Angew. Chem. Int. Ed. **59**, 10436 (2020).

¹³W. Ahn, J. F. Triana, F. Recabal, F. Herrera, and B. S. Simpkins, Science **380**, 1165 (2023).

¹⁴K. Hirai, H. Ishikawa, T. Chervy, J. A. Hutchison, and H. Uji-i, Chem. Sci. **12**, 11986 (2021).

¹⁵J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. X **9**, 021057 (2019).

¹⁶G. D. Scholes, J. Phys. Chem. A **125**, 6739 (2021).

¹⁷T. E. Li, A. Nitzan, and J. E. Subotnik, J. Chem. Phys. **154** (2021).

¹⁸T. E. Li, J. E. Subotnik, and A. Nitzan, Proc. Natl. Acad. Sci. USA **117**, 18324 (2020).

¹⁹X. Li, A. Mandal, and P. Huo, J. Phys. Chem. Lett. **12**, 6974 (2021).

²⁰J. Sun and O. Vendrell, J. Phys. Chem. Lett. **13**, 4441 (2022).

²¹C. Climent and J. Feist, Phys. Chem. Chem. Phys. **22**, 23545 (2020).

²²L. P. Lindoy, A. Mandal, and D. R. Reichman, Nature Comm. **14**, 2733 (2023).

²³A. Z. Lieberherr, S. T. E. Furniss, J. E. Lawrence, and D. E. Manolopoulos, J. Chem. Phys. **158**, 234106 (2023).

²⁴L. P. Lindoy, A. Mandal, and D. R. Reichman, J. Phys. Chem. Lett. **13**, 6580 (2022).

²⁵P.-Y. Yang and J. Cao, J. Phys. Chem. Lett. **12**, 9531 (2021).

²⁶S. Mondal, D. S. Wang, and S. Keshavamurthy, J. Chem. Phys. **157**, 244109 (2022).

²⁷Q. Yu and J. M. Bowman, Nature Comm. **14**, 3527 (2023).
<sup>28</sup>C. Schäfer, M. Ruggenthaler, and A. Rubio, Phys. Rev. A **98**, 043801 (2018).

<sup>29</sup>D. Sidler, M. Ruggenthaler, and A. Rubio, J. Chem. Theory Comput. **19**, 8801 (2023).

<sup>30</sup>D. Sidler, C. Schäfer, M. Ruggenthaler, and A. Rubio, J. Phys. Chem. Lett. **12**, 508 (2020).

<sup>31</sup>T. E. Li, A. Nitzan, and J. E. Subotnik, J. Chem. Phys. **156**, 134106 (2022).

<sup>32</sup>T. E. Li, A. Nitzan, and J. E. Subotnik, Angew. Chem. Int. Ed. **60**, 15533 (2021).

<sup>33</sup>T. E. Li and S. Hammes-Schiffer, J. Am. Chem. Soc. **145**, 377 (2022).

<sup>34</sup>T. E. Li, A. Nitzan, S. Hammes-Schiffer, and J. E. Subotnik, J. Phys. Chem. Lett. **13**, 3890 (2022).

<sup>35</sup>B. J. Berne, J. Stat. Phys. **43**, 911 (1986).

<sup>36</sup>T. Morresi, L. Paulatto, R. Vuilleumier, and M. Casula, J. Chem. Phys. **154**, 224108 (2021).

<sup>37</sup>T. Botzung, D. Hagenmüller, S. Schütz, J. Dubail, G. Pupillo, and J. Schachenmayer, Phys. Rev. B **102**, 144202 (2020).

<sup>38</sup>M. Du and J. Yuen-Zhou, Phys. Rev. Lett. **128**, 096001 (2022).

<sup>39</sup>G. Engelhardt and J. Cao, Phys. Rev. B **105**, 064205 (2022).

<sup>40</sup>M. Tavis and F. W. Cummings, Phys. Rev. **170**, 379 (1968).

<sup>41</sup>W. Liu, Int. J. Quantum Chem. **116**, 971 (2016).

<sup>42</sup>V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, J. Phys. B: At., Mol. and Opt. Phys. **51**, 034005 (2018).

<sup>43</sup>J. Dubail, T. Botzung, J. Schachenmayer, G. Pupillo, and D. Hagenmüller, Phys. Rev. A **105**, 023714 (2022).

<sup>44</sup>R. Feynman, Rev. Mod. Phys. **20**, 367 (1948).

<sup>45</sup>P. W. Anderson, Phys. Rev. **124**, 41 (1961).

<sup>46</sup>H. Bessone, R. Vuilleumier, and R. Spezia, J. Chem. Phys. **159**, 024112 (2023).

<sup>47</sup>S. Wang, Physical Review A **60**, 262 (1999).

<sup>48</sup>A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Nature Reviews Physics **1**, 19 (2019).

<sup>49</sup>T. E. Li, B. Cui, J. E. Subotnik, and A. Nitzan, Annual Review of Physical Chemistry **73**, 43 (2022).

<sup>50</sup>B. Cohn, S. Sufrin, A. Basu, and L. Chuntonov, The Journal of Physical Chemistry Letters **13**, 8369 (2022).

<sup>51</sup>B. M. Garraway, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences **369**, 1137 (2011).
⁵²C. Sommer, M. Reitz, F. Mineo, and C. Genes, Phys. Rev. Res. **3**, 033141 (2021).

⁵³D. S. Wang, T. Neuman, S. F. Yelin, and J. Flick, The Journal of Physical Chemistry Letters **13**, 3317 (2022).

⁵⁴S. Mondal and S. Keshavamurthy, J. Chem. Phys. **159**, 074106 (2023).

⁵⁵W. Wang, J. d. I. F. Diez, N. Delsuc, J. Peng, R. Spezia, R. Vuilleumier, and Y. Chen, Lab on a Chip **24**, 2497 (2024), publisher: Royal Society of Chemistry.

⁵⁶D. Goderis, A. Verbeure, and P. Vets, Journal of Statistical Physics **62**, 759 (1991).

⁵⁷D. Goderis, A. Verbeure, and P. Vets, Communications in Mathematical Physics **136**, 265 (1991).

⁵⁸I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. **121**, 3368 (2004).

⁵⁹G. Engelhardt and J. Cao, Phys. Rev. Lett. **130**, 213602 (2022).