Recent experiments have observed that the chemical and photophysical properties of molecules can be modified inside an optical Fabry–Pérot microcavity under collective vibrational strong coupling (VSC) conditions, and such modification is currently not well understood by theory. In an effort to understand the origin of such cavity-induced phenomena, some recent studies have focused on the effect of the cavity environment on the nonlinear optical response of the molecular subsystem. Here, we use a recently proposed protocol for classical cavity molecular dynamics simulations to numerically investigate the linear and the nonlinear response of liquid carbon dioxide under such VSC conditions following an optical pulse excitation. We find that applying a strong pulse of excitation to the lower hybrid light–matter state, i.e., the lower polariton (LP), can lead to an overall molecular nonlinear absorption that is enhanced by up to two orders of magnitude relative to the excitation outside the cavity. This polariton-enhanced multiphoton absorption also causes an ultrashort LP lifetime (0.2 ps) under strong illumination. Unlike usual polariton relaxation processes—whereby polaritonic energy transfers directly to the manifold of singly excited vibrational dark states—under the present mechanism, the LP transfers energy directly to the manifold of higher vibrationally excited dark states; these highly excited dark states subsequently relax to the manifold of singly excited states with a lifetime of tens of ps. Because the present mechanism is generic in nature, we expect these numerical predictions to be experimentally observed in different molecular systems and in cavities with different volumes.

## I. INTRODUCTION

Molecular polaritons, the hybrid quasi-particles stemming from strong light–matter interactions, open a new avenue to control molecular properties.^{1} Particular attention has recently been paid to the vibrational strong coupling (VSC) regime, where the interaction between a molecular ensemble and a cavity mode of a micrometer-size Fabry–Pérot cavity results in an observed Rabi splitting of order ∼100 cm^{−1}.^{2–4} Vibrational polaritons have been implicated in the observed modification of various molecular properties in the electronic ground state, including (i) changes in ground-state chemical reaction rates,^{5,6} reaction pathways’ selectivities,^{7} and even chemical equilibria^{8} without external pumping; (ii) modification of optical nonlinearities;^{9,10} and (iii) enhanced intermolecular vibrational energy transfer (VET) rates.^{11} Significantly, these modifications appear to be collective phenomena, originating from the interaction of cavity mode(s) with a large number of molecules.

From a theoretical perspective, collective optical response is easily understandable. In particular, vibrational Rabi splitting is easily modeled by mapping molecular vibrations and cavity modes onto harmonic oscillators.^{12} By contrast, cavity-induced chemical phenomena must involve a nuanced balance of collective and individual effects, and many questions remain as to the exact origin of the observed cavity-induced chemical effects.^{13–17} Indeed, the so-named “chemical catalysis” under VSC cannot be explained by a simple transition-state theory of chemical reactions because the potential of mean force along a reaction pathway is unchanged inside the cavity.^{14–16} Currently, descriptions of polariton relaxation,^{18–20} polariton-enhanced optical nonlinearity,^{21} and intermolecular VET rates^{11} under VSC are largely limited to either phenomenological master equations (which rely on the parameter fitting from experiments) or single-molecule analytical models that cannot address the collective aspects of the observed phenomena.

For a different perspective, we have recently proposed a numerical scheme to study cavity effects—cavity molecular dynamics (CavMD) simulations,^{22} in which classical dynamics is used to propagate coupled photon–nucleus dynamics for realistic molecular models. Compared with recently proposed theoretical methods in polariton chemistry that mainly focus on electronic structure calculations^{23–28} for a molecular electronic transition strongly coupled to a cavity mode (i.e., under electronic strong coupling), our approach serves as an affordable classical tool for describing VSC for a large ensemble of realistic molecules with full atomic resolution. Our approach^{22} captures the asymmetric Rabi splitting^{6} under VSC and even vibrational ultrastrong coupling when the O—H stretch mode of liquid water is strongly coupled to a cavity mode.

An important result from CavMD simulations of liquid water is that the static equilibrium properties of H_{2}O are completely unchanged inside vs outside the cavity, and dynamic response functions (evaluated in the linear response regime) of individual molecules also show very little or no effect. Thus, one must conclude that any significant VSC modification, e.g., VSC “catalysis” and/or the acceleration of an intermolecular VET rate, cannot be explained from an equilibrium or linear-response point of view. This conclusion suggests that cavity effects that apply to individual molecules may reflect cavity modification of the nonequilibrium dynamics of the molecular subsystem.

In this paper, we will present the results of such nonequilibrium CavMD simulations, looking at how a molecular system behaves following the application of an external pulse that excites a vibrational polariton; our goal is to explore if and how cavity effects modify molecular nonequilibrium properties. We will show that nonequilibrium CavMD simulations can recover experimental observations such as polariton relaxation to vibrational dark modes^{18} and a delay in the population gain of the singly excited states of vibrational dark modes after pumping the lower polariton (LP);^{19} moreover, such simulations also predict an intriguing process whereby a LP can enhance the molecular nonlinear absorption of light by up to two orders of magnitude, leading to very large molecular populations of highly excited vibrational states.

The model system studied is an ensemble of carbon dioxide molecules, where the C=O asymmetric stretch mode of liquid CO_{2} is nearly resonant with a cavity mode and forms lower and upper polaritons (UPs) under VSC; see Fig. 1 for the simulation setup. This model system resembles experimental systems studied recently with weak intermolecular interactions [such as W(CO)_{6}].^{18,19,29}

## II. METHODS

Here, we outline the theoretical considerations that underlie our CavMD simulations. A detailed account is given in Ref. 22. Under the Born–Oppenheimer approximation, the full-quantum photon–electron–nucleus Hamiltonian is projected onto the electronic ground state, which leads to a quantum Hamiltonian for the coupled photon–nucleus system only,

Here, $\u0124MG$ denotes the conventional ground-state molecular Hamiltonian,

where $P^nj$, $R^nj$, and *M*_{nj} denote the momentum operator, position operator, and mass for the *j*th nucleus in the molecule *n*, $V^g(n)$ denotes the intramolecular potential for the molecule *n*, and $V^inter(nl)$ denotes the intermolecular interactions between molecules *n* and *l*. $\u0124FG$ denotes the field-related Hamiltonian,^{14,30,31}

where $p\u0303^k,\lambda $, $q\u0303^k,\lambda $, *ω*_{k,λ}, and *m*_{k,λ} denote the momentum operator, position operator, frequency, and auxiliary mass for each cavity photon mode with the wave vector **k** and polarization direction *ξ*_{λ}. Note that the use of the auxiliary mass does not alter any dynamics and is necessary only because most MD packages require such mass; since the full light–matter coupling term in Eq. (1c) is proportional to $q^k,\lambda \u2261mk,\lambda \u2009q\u0303^k,\lambda $, where $q^k,\lambda $ is the standard mass-reduced photonic position operator, the final dynamics of the physical $q^k,\lambda $ operator and any molecular operator are not influenced by the magnitude of the auxiliary mass.^{32} Ω denotes the volume of the microcavity, *ϵ*_{0} denotes the vacuum permittivity, and $d^ng,\lambda $ denotes the electronic ground-state dipole operator for the molecule *n* projected along the direction of *ξ*_{λ}. In Eq. (1c), we can define

to characterize the coupling strength between each cavity photon and individual molecules. Note that Eq. (1c) assumes the long-wave approximation, i.e., the molecular ensemble is assumed to be much smaller than the wavelength of the cavity mode. Equation (1c) is exact when the cavity volume is large, e.g., in a microcavity where collective VSC is studied; otherwise, an additional self-dipole fluctuation term will also emerge due to the quantum nature of electrons.^{14} Compared with most model Hamiltonians such as the Tavis–Cummings Hamiltonian, the most different feature in Eq. (1c) is the inclusion of the self-dipole term, i.e., the term that is quadratic in $d^ng,\lambda $. We emphasize that including this self-dipole term is critical in CavMD simulations because this term preserves gauge invariance, contributes to the asymmetry of Rabi splitting,^{22} and, most importantly, influences the long-time molecular dynamics and the reliability of CavMD results. For example, for thermal equilibrium simulations, including the self-dipole term guarantees that static molecular properties will be unchanged under VSC,^{22} while neglecting such a term will cause numeric artifacts (such as changed static molecular properties). See also Refs. 27 and 33 that discuss the importance of the self-dipole term in other contexts.

After reducing all of the operators in Eq. (1) to classical observables and also applying periodic boundary conditions for the molecules, we arrive at equations of motion for the coupled photon–nucleus system,

Here, the subscript *nj* denotes the *j*th atom in the molecule *n*; $Fnj(0)=\u2212\u2202Vg(n)/\u2202Rnj\u2212\u2211l\u2260n\u2202Vinter(nl)/\u2202Rnj$ denotes the molecular part of the force on each nucleus;

denotes the cavity force on each nucleus; and $Fnjext(t)=\u2212QnjEext(t)$ denotes an external driving force acting on each nucleus with partial charge *Q*_{nj} under the pumping of a time-dependent electric field **E**_{ext}(*t*). Note that this $Fnjext(t)$ was not introduced in Ref. 22 and is included here to represent the optical pulse excitation that leads to the molecular nonequilibrium response.^{34} In order to apply periodic boundary conditions, in Eq. (2), we have redefined $q\u2248k,\lambda =q\u0303k,\lambda /Ncell$, where *N*_{cell} denotes the number of periodic simulation cells for molecules, and also an effective coupling strength

where *ε*_{k,λ}, the true coupling strength between each cavity photon and individual molecules, is defined in Eq. (1d). By invoking periodic boundary conditions as above, CavMD simulations can yield the same Rabi splitting when calculating the consequence of photons interacting with molecules in a simulation cell as found in the original system. With the number of molecules in a single simulation cell denoted as *N*_{sub}, the total number of molecules is *N* = *N*_{sub}*N*_{cell}. Note that when CavMD is used to reproduce experimental observations such as polariton relaxation, given the value of *N*_{sub}, $\epsilon \u0303k,\lambda \u221dNcell$ can be chosen by fitting the experimentally observed Rabi splitting Ω_{N}; since $\Omega N\u221dN$, $\epsilon \u0303k,\lambda \u221dNcell\u221d\Omega N/Nsub$. Also note that if CavMD is used to simulate VSC phenomena in Fabry–Pérot microcavities with *N* molecules using parameters designed to recapitulate an experimental Rabi splitting Ω_{N}, it is necessary to check the dependence on periodic boundary conditions (or the choice of *N*_{sub}) to make sure that any observed dynamics are not an artifact of the simulation; after all, the effective coupling strength for molecules $\epsilon \u0303k,\lambda $ has been amplified by a factor of $Ncell$ relative to the true coupling strength *ε*_{k,λ}. In Sec. IV D, we will study how our CavMD results depend on *N*_{sub}, all while keeping the Rabi splitting Ω_{N} fixed and adjusting $\epsilon \u0303k,\lambda $ accordingly. For such calculations, the asymptotic results when *N*_{sub} approaches a macroscopic number should correspond to Fabry–Pérot microcavities. For the same calculations, it may also be helpful to imagine a physical experiment with *N*_{cell} = 1 (so that $\epsilon \u0303k,\lambda =\epsilon k,\lambda $ becomes the true coupling strength); here, one can interpret our CavMD dynamics as reliably reporting on cavities with different effective volumes (and, therefore, different effective molecular numbers *N*).

### A. Molecular spectroscopy

Below, we will calculate two different spectroscopic response functions: the global infrared (IR) absorption spectrum and its “local” correspondence that is the spectrum obtained if the molecules respond to light individually. The former is calculated by Fourier transforming the dipole auto-correlation function,^{35–38}

while the latter is defined by

In Eqs. (3) and (4), *α*(*ω*) or *α*_{local}(*ω*) denotes the absorption coefficient, *n*(*ω*) denotes the refractive index, *V* is the volume of the system (i.e., the simulation cell), **e**_{i} denotes the unit vector along the direction *i* = *x*, *y*, and *μ*_{S}(*t*) denotes the *total* dipole moment of the molecules at time *t*, which is *μ*_{S}(*t*) = ∑_{nj}*Q*_{nj}**R**_{nj}(*t*). Note that, in our force field calculations, the partial charges (*Q*_{nj}, the true nuclear charge for the nucleus *nj* plus the electronic shielding effect) of nuclei are fixed; see Appendix C for the values of partial charges. *μ*_{n}(*t*) = ∑_{j}*Q*_{nj}**R**_{nj}(*t*) denotes the dipole moment for the molecule *n*. In Eq. (3), the summation over *x* and *y* is a summation over the two possible polarizations of the relevant cavity mode of the *z*-oriented cavity (see Fig. 1 for the simulation setup). In Eq. (4), the inner product implies that a summation over three dimensions has been applied, though summing over only two dimensions *x*, *y* yields the same local IR spectra up to a prefactor.

It is important to note that *α*(*ω*) of Eq. (3) represents the molecular IR absorption (for a molecular sample much smaller than the wavelength of light). Its local correspondence, *α*_{local}(*ω*) of Eq. (4), corresponds to the absorption of a fictitious molecular system in which each molecule responds to the field individually. Thus, Eq. (3) describes the collective behavior of the molecular dipole system, while Eq. (4) provides information about the dynamics of individual molecules. Just as the collective bright and dark modes of the molecular ensemble can be expressed as a linear combination of individual molecular modes, the individual response can also be expressed as a linear combination of the collective modes—note, though, that the latter is predominately dominated by the dark modes. For instance, if *N*_{sub} = 216 is taken during the simulation, the bright state contribution to the local IR spectrum in Eq. (4) is negligibly small, equal to 1/*N*_{sub} = 1/216. Therefore, one can roughly interpret the local IR spectrum as reporting dark-mode dynamics.

Last, a few words are appropriate as far as interpreting the response functions above. Quantum mechanically, one calculates absorption and emission spectra differently; even though one must propagate the same Hamiltonian in both cases, absorption and emission can be differentiated by their respective initial conditions and the fact that quantum correlation functions are not time-reversible, $\xc2B^(t)\u2260\xc2B^(\u2212t)$. Classically, correlation functions (and their spectra) report on the energy present in a given mode, which cannot be naturally dissected into absorption or emission components, and this inability can lead to some confusion when analyzing MD simulations and looking to connect with experimental spectra (which reflect true quantum mechanical dynamics). Nevertheless, by following the vibrational energy evolving in time during a given trajectory, we will be able to semiclassically rationalize how the cavity mode, the bright mode, and the dark modes relax and exchange energy.

## III. SIMULATION DETAILS

A sketch of the cavity structure used in our simulation is plotted in Fig. 1. The cavity is placed along the *z*-axis. In addition to the (assumed perfect) mirrors that define the Fabry–Pérot cavity, we assume that additional parallel inert layers (e.g., SiO_{2}, which was used in Ref. 5) confine the molecular system in a region small enough to guarantee the validity of the long wave approximation and also far enough from the mirrors so that image charge effects can be disregarded. Both assumptions are made in writing the model Hamiltonian (1). Only one cavity mode with frequency near the C=O asymmetric stretch is considered. This cavity mode contains two polarization directions along the *x*-direction and *y*-direction. We set the auxiliary mass for the cavity mode as *m*_{k,λ} = 1 a.u. (atomic units) though, as mentioned above, this mass does not affect any dynamics. Under periodic boundary conditions, we simulate 216 CO_{2} molecules in a cubic cell with cell length 24.292 Å (45.905 a.u.); the density of the liquid CO_{2} is 1.101 g/cm^{3}. Unless stated otherwise, by default, we set the cavity mode frequency as 2320 cm^{−1} with an effective coupling strength $\epsilon \u0303=2\xd710\u22124$ a.u. Compared with VSC experiments for which usually *N* = 10^{9}–10^{11} molecules are involved, our choice of the effective coupling strength $\epsilon \u0303\u221dNcell$ should correspond to the involvement of *N*_{cell} = 10^{7}–10^{9} periodic simulation cells.

When calculating equilibrium properties, we perform simulations as follows. At 300 K, we first run the simulation for 150 ps to guarantee thermal equilibrium under an NVT (constant particle number, volume, and temperature) ensemble where a Langevin thermostat with a lifetime (i.e., inverse friction) of 100 fs is applied to the momenta of all particles (nuclei + photons). The resulting equilibrium configurations are used as starting points for 40 consecutive NVE (constant particle number, volume, and energy) trajectories of length 20 ps. At the beginning of each trajectory, the velocities are resampled by a Maxwell–Boltzmann distribution under 300 K. The intermolecular Coulombic interactions are calculated by an Ewald summation. The simulation step is set as 0.5 fs, and we store the snapshots of trajectories every 2 fs.

When performing nonequilibrium simulations under an external pulse, we start each simulation with an equilibrium geometry, which is chosen from the starting configurations of the above 40 NVE trajectories. Each nonequilibrium trajectory is run for 100 ps under an NVE ensemble, and the physical properties are calculated by averaging over these 40 nonequilibrium trajectories. Note that the use of NVE trajectories when calculating equilibrium or nonequilibrium physical properties implies the assumption that cavity losses are small on the timescale of simulations and can be ignored. This situation is usually valid when considering Fabry–Pérot microcavities where the cavity loss lifetime usually takes ∼5 ps, while polariton relaxation to vibrational dark modes usually occurs on a timescale <5 ps; see Fig. 3 and also experiments.^{19} Of course, for cavities with larger losses or when the polariton relaxation process occurs on a longer timescale, the polariton relaxation might be accelerated and the resulting linear and nonlinear vibrational dark-mode energy absorption (in Fig. 4) would then also be weakened.

The form of the external pulse is taken as follows:

where the phase *ϕ* ∈ [0, 2*π*) is set as random. This pulse is turned on at *t*_{start} = 0.1 ps and is turned off at *t*_{end} = 0.6 ps. Below, we will show results obtained after weak pumping, i.e., *E*_{0} = 3.084 × 10^{6} V/m (6 × 10^{−4} a.u.) and the input pulse fluence $F=12\u03f50cE02(tend\u2212tstart)=6.32$ mJ/cm^{2}, and also after strong pumping, i.e., *E*_{0} = 3.084 × 10^{7} V/m (6 × 10^{−3} a.u.) and *F* = 632 mJ/cm^{2}. The choice of an *x*-polarized pumping pulse implies that molecules with the dipole component in the *x*-direction can be excited, and because the cavity sits along the *z*-direction, the *x*-polarized pulse [in Eq. (5)] can excite the polariton associated with the cavity mode polarized along the same *x*-direction. Throughout this paper, we will refer to exciting the polaritons as exciting the molecular ensemble with such an *x*-polarized pulse. Similarly, polaritons can also be excited with a *y*-polarized pulse.

For the force field of CO_{2} (see Appendix C for details), we largely follow Ref. 39. The only difference is that while Ref. 39 uses a harmonic potential for the C=O bond, we change this harmonic potential to the following anharmonic form:

where Δ*r* = *r* − *r*_{eq}. Equation (6) is a fourth-order Taylor expansion of a Morse potential $VM(r)=Dr[1\u2212exp(\u2212\alpha r\Delta r)]2$. The parameters are taken as follows: *r*_{eq} = 1.162 Å (2.196 a.u.), *D*_{r} = 127.13 kcal mol^{−1} (0.2026 a.u.), and *α*_{r} = 2.819 Å^{–1} (1.492 a.u.) are chosen to fit the harmonic potential used in Ref. 39 in the harmonic limit and the value of *D*_{r} takes the bond dissociation energy of O=CO at room temperature.^{40} A comparison of Eq. (6), the Morse potential, and the harmonic limit is plotted in Fig. 8. When the C=O bond energy is smaller than ∼0.05 a.u. (10^{4} cm^{−1}), Eq. (6) agrees with the Morse potential very well. However, the Morse potential is not used for the present simulation as we wish to avoid any potential molecular dissociation. At the same time, note that the use of an anharmonic potential (rather than a harmonic one) is critical for the present paper because, as will be shown in Sec. IV C, polariton-enhanced molecular nonlinear absorption will rely on a non-uniform distribution of molecular energy level spacing. This reliance also highlights the fact that, in order to capture (at least some) nontrivial VSC phenomena, it is crucial to use realistic molecules instead of conventional harmonic models.

As far as the technical details are concerned, the initial configuration is prepared with PACKMOL,^{41} the CavMD scheme is implemented by modifying the I-PI package,^{42} and the nuclear forces are evaluated by calling LAMMPS.^{43} A toolkit including the source code, input and post-processing scripts, and the corresponding tutorials is available in Github.^{44}

## IV. RESULTS AND DISCUSSION

### A. Rabi splitting and avoided crossing

In Fig. 2(a), we plot the IR spectrum obtained from Eq. (3) for different values of the effective coupling strength $\epsilon \u0303$; the value of $\epsilon \u0303$ (in a.u.) is labeled on each lineshape, and the case of molecular system outside the cavity corresponds to $\epsilon \u0303=0$. The Rabi splitting seen in Fig. 2(a) confirms the existence of a strong coupling between the C=O asymmetric stretch and the cavity mode at frequency 2320 cm^{−1}. Note that, with the periodic boundary conditions applied during CavMD simulations, the effective coupling strength $\epsilon \u0303$ scales as $\epsilon \u0303\u221dN$ with the total number of molecules *N*; see Sec. II for details. Therefore, increasing $\epsilon \u0303$ provides a simple way to study the Rabi splitting as a function of the total number of molecules *N*. As shown in Fig. 2(a), unlike the lineshape outside a cavity (bottom), inside a cavity, a pair of the UP and LP is formed under VSC and the Rabi splitting increases with $\epsilon \u0303$. The inset plots the Rabi frequency, or the frequency difference between the UP and LP peaks, as a function of $\epsilon \u0303\u221dN$. Here, a linear relationship is observed, which agrees with both analytical models of coupled harmonic oscillators and many experiments. Note that the asymmetry in the positions and amplitudes of the UP and LP seen in Fig. 2(a) is discussed in detail in Ref. 22.

Figure 2(b) plots Rabi splitting as a function of the cavity mode frequency for the condition $\epsilon \u0303=2\xd710\u22124$ a.u. When the cavity mode frequency is highly negatively detuned—i.e., when the cavity frequency is much smaller than the bare C=O asymmetric stretch at 2327 cm^{−1} [white dashed line (which is very close to the experimental value, e.g., 2333 cm^{−1} in Ref. 45)]—the UP is close to the bare C=O asymmetric stretch and its character is dominated by this molecular vibrational mode; at the same time, the LP is close to the cavity mode frequency (green dashed line) and its character is dominated by the cavity mode. By contrast, when the cavity mode frequency is highly positively detuned, the LP (UP) is mostly contributed by the molecular vibration (cavity mode). The avoided crossing seen between these limits expresses the Rabi splitting that measures the collective coupling strength.

### B. Polariton relaxation and ultrashort LP lifetime

Next, we use nonequilibrium CavMD simulations to explore the polariton relaxation following pulse excitation of the UP or the LP [see Eq. (5)]. Because polaritons are hybrid light–matter quasi-particles, their relaxation can be captured by monitoring either the photonic or the matter side. The photonic energy ($\u2211\lambda =x,ymk,\lambda \omega k,\lambda q\u2248k,\lambda 2/2+p\u2248k,\lambda 2/2mk,\lambda $) is simpler to calculate and is used in this calculation.^{46}

Figure 3(a) shows the time-resolved photonic energy, in a system where the cavity mode frequency is set to 2320 cm^{−1} and effective coupling strength $\epsilon \u0303=2\xd710\u22124$ a.u., after resonantly exciting the UP (magenta; peaked at 2428 cm^{−1}) or LP (cyan; peaked at 2241 cm^{−1}) with a weak incoming pulse **E**(*t*) = *E*_{0} cos(*ωt* + *ϕ*)**e**_{x} (with fluence *F* = 6.32 mJ/cm^{2}), where the yellow-shadowed region denotes the 0.5 ps time window during which the pulse is applied. The fast oscillations (with a period ∼0.2 ps) of the cavity photon signals appear to reflect oscillations between the cavity photons and the vibrational bright mode, as here the Rabi splitting is 187 cm^{–1} = 0.18 ps. By fitting the energy decay observed after the pulse to an exponential function, the polariton lifetime can be captured: the UP lifetime is *τ*_{UP} = 2.5 ps and the LP lifetime is *τ*_{LP} = 1.0 ps. By contrast, under a strong incoming pulse (*F* = 632 mJ/cm^{2}), the same plot in Fig. 3(b) shows that while the UP lifetime is largely unchanged (*τ*_{UP} = 3.3 ps), the LP shows an ultrafast decay with a very short lifetime, *τ*_{LP} = 0.2 ps. Note that this residual decay follows the 0.5 ps pumping pulse, implying that much of the LP relaxation has already taken place during the pumping stage. Under weak illumination [Fig. 3(a)], the cavity photon population can reach 0.01 a.u. (≈*ℏω*_{0} = 2327 cm^{−1}), implying that the system receives about one quantum of energy and stays on the singly excited manifold. By contrast, under strong illumination [Fig. 3(b)], the cavity photon population signals can reach 1.0 a.u. (≈100*ℏω*_{0}), meaning that the system is very highly excited.

Because we have assumed no cavity loss and we find that, outside a cavity, vibrational relaxation of an individual molecule to the ground state takes much longer than several ps, the fast (<5 ps) polariton relaxation observed in Fig. 3 must reflect energy transfer to the vibrational dark modes of the asymmetric C=O stretch. Therefore, the fact that the LP relaxation is faster than that of the UP implies either a stronger interaction between the LP and individual molecular motions or the existence of a decay channel for the LP that is not open for the UP. Below, we provide evidence in support of the latter scenario. This new decay channel, which we will call polariton-enhanced molecular nonlinear absorption, can exist under both weak and strong illumination of the LP (see Fig. 6 for details). Due to this new decay channel, in both Figs. 3(a) and 3(b), the LP signal exhibits a shortened height and lifetime compared to the UP signal. A more detailed study of vibrational polariton relaxation combined with analytical theory and CavMD simulations will be given elsewhere.

### C. Polariton-enhanced nonlinear absorption

In order to demonstrate that polariton-enhanced molecular nonlinear absorption is the origin of the ultrashort LP lifetime in Fig. 3(b), we next study the molecular response to the pulse excitation as expressed by its transient IR spectrum. Here, the transient IR spectrum about time *T*_{i} is calculated by evaluating Eq. (3) in a time window [*T*_{i}, *T*_{i} + Δ*T*], where Δ*T* = 5 ps. Note that, in this classical calculation, the spectrum in Eq. (3) yields information about the molecular frequency distribution, and the corresponding transient spectrum corresponds to this distribution in the excited molecular ensemble.

#### 1. Cavity mode at 2320 cm^{−1}

For the same parameters as in Fig. 3(b), Fig. 4(a) shows the time-resolved IR spectra after exciting the LP at 2241 cm^{−1} (the red arrow) with a strong pulse (*F* = 632 mJ/cm^{2}). Two observations can be made: (a) there is an ultrafast relaxation of the LP signal that disappears almost immediately after the exciting pulse and (b) a distribution of lower frequencies, peaked at ∼2170 cm^{−1}, emerges. The existence of such lower frequencies indicates the appearance of molecules with higher vibrational energies (see Appendix B where we establish an explicit semiclassical relationship between the vibrational frequency and vibrational quanta for our anharmonic CO_{2} force field). In particular, for the anharmonic potential that we use to simulate CO_{2} dynamics, the frequency ∼2200 cm^{−1} roughly corresponds to the motion of a classical anharmonic oscillator whose amplitude is determined by having two quanta of vibrational energy. Quantum mechanically, this frequency corresponds to the energy of the 1 → 2 vibrational transition. This suggests that the additional peak ∼2170 cm^{−1} is a signal of nonlinear multiphoton absorption.

For the same transient state that yields Fig. 4(a), Fig. 4(d) plots the corresponding time-resolved local IR spectra obtained from Eq. (4) that reflect information about the individual molecular modes (which are mostly composed of vibrational dark modes) rather than a collective bright state. During and after excitation of the LP at 2241 cm^{−1} (red arrow), a large fraction of the energy is transferred to the higher vibrational excited states of the C=O asymmetric stretch, leading to a broad distribution of frequencies and showing a peak at ∼2170 cm^{−1} (see the frequency region in the cyan rectangle); we also find a peak around 2327 cm^{−1} (see the frequency region in the blue rectangle) that corresponds to the absorption of individual molecules, where the asymmetric C=O stretches are in a thermal distribution. As time increases, the excess energy in the aforementioned higher excited states gradually relaxes to the latter thermal distribution of vibrational states. This interpretation can be validated by noting that the fitted lifetimes for the energy decay (39 ps) and gain (34 ps) processes are consistent. In other words, due to the polariton-enhanced molecular nonlinear absorption mechanism, the population gain of the C=O asymmetric stretches in a thermal distribution (with predominantly 0 or 1 vibrational quanta) is significantly delayed compared with the ultrafast LP lifetime.

Interestingly, a recent experiment performed by Xiang *et al.*^{19} reported that, for W(CO)_{6}, after an excitation of the LP, there is a significant delay in the gain of population in the vibrationally first excited state of the dark-state manifold [which corresponds to the lower vibrational excited states in the blue rectangle of Fig. 4(d)]. In Ref. 19, the authors proposed that the underlying mechanism should be the direct transition from the LP to second or higher vibrational excited states and the subsequent relaxation to the first excited state. In fact, after the submission of this article, via private communication, we have learned^{47} of experimental data that would appear to suggest that, for W(CO)_{6}, pumping the LP leads to excitation of a vibrationally second excited state, followed by a delay and then population of a first vibrationally excited state (in agreement with our numerical results above). Note that their results have a weaker nonlinear signal presumably because the excitation of the LP was not very strong in their work.

Now, although the experiments in Ref. 19 would appear to validate the CavMD results in Fig. 4(d), a keen reader may still be very surprised that a classical simulation can produce a two-peak feature in Fig. 4(d). While, within a quantum model, there are several different vibrational transition energies for an anharmonic oscillator, it is widely known that a classical response function should predict only a single vibrational peak for a single anharmonic oscillator driven by an external field.

In order to understand the origin of the two-peak feature in Fig. 4(d), we directly study the probability density for the C=O bond potential energy. Figure 5(a) plots this probability distribution (in the logarithmic scale) at an early time just after the pulse (*T*_{i} = 1 ps)—under the same conditions as in Fig. 4(d), Fig. 5(a) clearly shows that the C=O vibrational energy probability distribution has not only a strong peak in the low (thermal) energy regime but also a large shoulder for states above 2*ℏω*_{0} (*ω*_{0} = 2327 cm^{−1}). As shown in Fig. 9, since a higher vibrational energy corresponds to a lower vibrational frequency within the anharmonic C=O force field, the strong peak and the large shoulder in Fig. 5(a) will lead to two peaks in Fig. 4(d) if we consider an ensemble average of all molecules.

#### 2. Excitation outside the cavity

The data above suggest that, for some experiments, exciting the system at the LP frequency can facilitate large nonlinear absorption of energy. To confirm that the LP is in fact facilitating such an effect, in Figs. 4(c) and 4(f), we plot the corresponding IR and local IR spectra when the molecules are excited outside the cavity with the same strong pulse (centered at 2241 cm^{−1}) as in Figs. 4(a) and 4(d). A very weak nonlinear absorption peaked at ∼2200 cm^{−1} can be detected; see the dashed cyan region. As mentioned above, this ∼2200 cm^{−1} weak peak may correspond to a small number of molecules with two or more quanta of vibrational energy. Note that this interpretation agrees with the corresponding C=O bond potential energy distribution in Fig. 5(c), where we find a very small fraction of molecules have C=O bond potential energies near 2*ℏω*_{0}, and the large population of highly excited molecules seen in Fig. 4(d) is not reproduced here. Comparing the results inside vs outside the cavity, we conclude that the LP in Figs. 4(a) and 4(d) can greatly enhance the multiphoton nonlinear absorption of molecules. Note that polariton-enhanced multiphoton absorption under strong illumination has been shown experimentally for various setups, such as organic excitons in a Fabry–Pérot cavity^{48} and quantum dots near surface plasmons,^{49,50} while its possibility under collective VSC has not been extensively studied^{19,21} especially under strong illumination. Note also that the reader should not be hesitant in deriving such a conclusion based on the fact that our calculations are entirely classical; classical simulations have long been known to capture molecular multiphoton nonlinear absorption outside a cavity.^{51}

#### 3. Cavity mode at 2200 cm^{−1}

Finally, let us consider a cavity with a non-resonant frequency mode; such a consideration will lead to a new understanding of the conditions necessary for a polariton to facilitate multiphoton nonlinear absorption. Consider the case when the cavity mode frequency is slightly off resonance with respect to the molecular ground-state frequency. Figures 4(b) and 4(e) show results of similar simulations for a system with cavity mode frequency 2200 cm^{−1}. Now, when the LP (peaked at 2167 cm^{−1}) is resonantly excited with a strong pulse (*F* = 632 mJ/cm^{2}), as shown in Figs. 4(b) and 4(e), the large nonlinear absorption at initial times does not appear; see also Fig. 5(b) for the corresponding C=O bond potential energy distribution.

Here, readers might be curious why the LP cannot facilitate molecular nonlinear absorption even when the LP frequency 2167 cm^{−1} is very close to the frequency of higher vibrational excited states. We hypothesize that this fact can be best understood within a quantum picture. In a quantum picture, when twice the LP frequency matches the 0 → 2 vibrational transition, one can expect an enhancement of molecular nonlinear absorption, so that the LP can function as a “virtual state” for exciting dark modes. By contrast, if the LP frequency matches the 1 → 2 transition (similar to the present case), the total effect of exciting the LP should be small, since vibrational dark modes obey a thermal distribution and the $v$ = 1 population is small at room temperature. Overall, as a rule of thumb, according to Fig. 5, it would appear that molecular nonlinear absorption is facilitated when the LP frequency sits in between 2327 cm^{−1} (the frequency of the CO_{2} molecules in a thermal distribution) and 2170 cm^{−1} (the frequency of the highly excited CO_{2} molecules).

In the present case, the LP frequency does not match the frequency needed to move the molecule above its first excited state (i.e., the virtual state near 2241 cm^{−1} as mentioned above), so the LP cannot function as a resonant gateway for nonlinear molecular absorption. Consequently, the LP lifetime becomes much longer: An exponential fit of the LP intensity in the time-resolved IR spectra [Fig. 4(b)] gives *τ*_{LP} = 8.5 ps, while a consistent result of *τ*_{LP} = 7.5 ps can be obtained by an exponential fit of the photonic energy (which is calculated in the same way as in Fig. 3). In this relatively slow relaxation process, energy is directly transferred from the LP to the lower excited states of individual molecules and there is no delay in the population gain of the vibrational singly excited manifold.

#### 4. Pulse intensity dependence of nonlinear absorption

The cavity effect in enhancing nonlinear absorption can be further seen by studying the nonlinear dependence on the pulse fluence in and outside the cavity. Figure 6 shows the local IR spectrum from Fig. 4(d), integrated over the frequency region that corresponds to the nonlinear absorption inside a 2320 cm^{−1} cavity [the same as in Figs. 4(a) and 4(d)] and outside a cavity, as a function of the pulse fluence. In order to make a quantitative analysis, we define the relevant integration region as the area colored cyan in Fig. 4(d) and evaluate this integral at the early time *T*_{i} = 1 ps. This early-time nonlinear signal (*I*_{NL}) provides a direct means to quantify the magnitude of the nonlinearity at the beginning stages of absorption. When the incoming pulse is weak (the green-shadowed region), the nonlinearity is weakly enhanced inside the cavity (blue circles) relative to the molecular response outside the cavity (black squares).^{52} By contrast, when the incoming pulse is strong (the red-shadowed region), *I*_{NL} increases at both in and outside the cavity, and under the same pulse fluence, the nonlinearity can be enhanced by up to two orders of magnitude inside the cavity in comparison with the free space case. This provides a direct demonstration of the cavity role in promoting and maintaining multiphoton nonlinear absorption inside an optical cavity.

### D. Effect of periodic boundary conditions

Before ending our article, we emphasize that there is a huge numerical gap between realistic VSC experiments and our CavMD simulations. In experiments, a macroscopic number of molecules form VSC in Fabry–Pérot cavities and the coupling to a cavity mode for each molecule (*ε*) is very small. By contrast, in CavMD simulations, hundreds of molecules are explicitly simulated and the light–matter interaction per molecule $\epsilon \u0303=Ncell\epsilon $ is artificially amplified proportionally to the number of periodic simulation cells (*N*_{cell}) due to the invoking of periodic boundary conditions (see Ref. 22 for details). In order to validate the intriguing simulation results above, we have performed additional simulations to investigate the effect of the imposed periodic boundary conditions on the simulation results. When macroscopic observables (such as Rabi splitting and molecular density) are kept the same, we increase the number of molecules in the simulation cell (*N*_{sub}) and study the polariton lifetime dependence on *N*_{sub}. Note that fixing Rabi splitting (i.e., fixing *N* = *N*_{sub}*N*_{cell}) while increasing *N*_{sub} will lead to a decrease in $\epsilon \u0303$ ($\epsilon \u0303\u221dNcell=N/Nsub$) and, therefore, the reduction of any artifacts due to the use of periodic boundary conditions.

Figure 7 shows that the polariton lifetimes are unchanged against 1/*N*_{sub} under the same conditions as in Fig. 3: the left panel plots the LP lifetime under the weak [Fig. 7(a); *F* = 6.32 mJ/cm^{2}] or strong [Fig. 7(c); *F* = 632 mJ/cm^{2}] pulse as in Fig. 3; the right panel plots the UP lifetime accordingly. For all situations [including the ultrashort LP lifetime in Fig. 7(c)], the results show no dependence on 1/*N*_{sub}, suggesting that the reported simulation results are not sensitive to renormalization applied to the light–matter interaction ($\epsilon \u0303$) and implying that, as long as the observed Rabi splitting is fixed, these predictions may hold in cavities with different volumes, ranging from plasmonic cavities to Fabry–Pérot cavities, the latter of which are usually used for VSC experiments.

Of course, because small-volume cavities usually have large cavity losses (while we have ignored cavity losses in our simulations), energy transfer between the polariton and the dark modes may not be as strong for such experiments as in the present results. In Fabry–Pérot cavities, however, because the cavity loss lifetime (∼5 ps) is longer than the polaritonic relaxation lifetimes in Fig. 7, the presented results should not be meaningfully altered even when the cavity loss is considered.

## V. CONCLUSION

To summarize our observations, by studying how molecules respond both individually and collectively after a polariton has been weakly excited, we have found that polariton relaxation to molecular dark modes usually occurs on a timescale of several ps, in agreement with previous experiments.^{18,19} However, when a strong pulse is applied to the LP in a suitable cavity and the LP energy can support transitions to higher molecular states, the LP lifetime can become ultrashort (0.2 ps) and one can find individual molecules in very highly excited vibrational states. This so-called multiphoton nonlinear molecular absorption of light can be as large as two orders of magnitude relative to the excitation outside the cavity and arises only in concert with the LP (not the UP) because, for a realistic vibration, the 0 → 1 transition always has a larger frequency than the 1 → 2 transition so that the 0 → 2 transition can approximately match twice the LP frequency.

Given that this LP relaxation behavior is robust within our simulation, and especially to the choice against of periodic boundary conditions used in our CavMD simulations, we have every reason to believe that our microscopic simulations presented here will have real macroscopic experimental consequences in cavities with different volumes, and we expect such intriguing ultrashort LP lifetime and LP-enhanced molecular nonlinear absorption will soon be experimentally verified in usual VSC setups such as Fabry–Pérot cavities (where a macroscopic number of molecules form VSC). As highlighted in Sec. IV C, recent experiments by Xiang *et al.* have observed a delayed population gain in the singly excited dark state manifold following an excitation of the LP,^{19} which would appear to be a strong endorsement of the current CavMD simulations (while also raising the possibility of even more dramatic findings if the LP is strongly illuminated).

In the present study, we have investigated VSC using an exclusively classical approach. Although a quantum approach would be ideal for studying light–matter interactions, any brute-force approach is not feasible because VSC phenomena involve a large number of molecules. Moreover, given the large number of vibrationally highly excited states, and the complex interactions between bright and dark modes as caused by short-range intermolecular interactions, no easy simplification of the problem seems feasible. Importantly, the effects discussed in the present study are all classical in nature. Given the excellent agreement between our classical CavMD simulations and VSC relaxation experiments, classical CavMD simulations would appear to be a promising tool to investigate VSC-related phenomena, at least qualitatively. This statement is consistent with previous findings that classical theories can qualitatively capture many intriguing light–matter interaction processes that have classical analogs.^{53,54} Obviously, purely quantum effects such as entanglement will need a quantum treatment. Looking forward, the open question remains as to if and how the present CavMD simulations relate to the key outstanding and unexplained VSC phenomenon: chemical catalysis.

## ACKNOWLEDGMENTS

This material is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award No. DE-SC0019397 (J.E.S.) and U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division (A.N.). This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

## DATA AVAILABILITY

The data that support the findings of this study (including the source code, input and post-processing scripts, and the corresponding tutorials) are openly available in Github at https://github.com/TaoELi/cavity-md-ipi.^{44}

### APPENDIX A: POTENTIAL OF C=O BOND

The potential of the C=O bond in our force field is plotted in Fig. 8.

### APPENDIX B: CORRESPONDENCE BETWEEN VIBRATIONAL FREQUENCY AND ENERGY

In simulations, we have found the emergence of a low-frequency absorption after exciting the LP. The appearance of low frequencies comes from the use of an anharmonic C=O potential as described in Fig. 8. Due to the intrinsically anharmonic nature of molecular vibrations, a higher quantum vibrational transition (e.g., *n* → *n* + 1, where *n* is large) has a lower frequency than the fundamental 0 → 1 transition. Classically speaking, an anharmonic oscillator also demonstrates a lower, red-shifted frequency when this oscillator contains a higher energy (or classical action); see Ref. 55 for details. In order to quantify the vibrational state for the corresponding frequencies, we perform a simple simulation for a bare CO_{2} in the gas phase after an external pulse peaked at 2327 cm^{−1} with duration 0.1 ps; see Eq. (5). Because this pulse has a shorter duration than that used in the main text (0.5 ps), it contains a wide spectrum band and can excite the molecule to a wide range of vibrational excited states by increasing the pulse fluence. The simulation is performed for 20 ps in an NVE ensemble, and the initial configuration of the molecule is set at the global minimum of potential energy surface with no initial velocity. Under different pulse amplitudes, the C=O asymmetric stretch can oscillate with different vibrational energies.

Figure 9 plots the corresponding C=O asymmetric peak frequency by evaluating the dipole auto-correlation function in Eq. (3) as a function of the vibrational energy of the molecule, which is predominately contributed by the vibrational energy of the C=O asymmetric stretch mode. While, at thermal equilibrium, the C=O asymmetric stretch peaks at *ω*_{0} = 2342 cm^{−1}, Fig. 9 shows that the peak frequency exhibits a negative relationship with the vibrational potential energy. For example, the frequency near 2200 cm^{−1} corresponds to roughly two vibrational quanta. According to the correspondence between frequency and vibrational energy, we classify the frequency range (2220–2360) cm^{−1} as lower excited states (or linear absorption; see the blue region) and the frequency range (2150–2220) cm^{−1} as higher excited states (or nonlinear absorption; the cyan region).

### APPENDIX C: CO_{2} FORCE FIELD

For completeness, below we will provide the details of the CO_{2} force field. This force field largely resembles that in Ref. 39 except for the use of an anharmonic C=O bond potential in Eq. (6).

In this CO_{2} force field, the pairwise intermolecular potential is characterized by the Lennard-Jones potential ($VLJ(nl)$) plus the Coulombic interactions between atoms ($VCoul(nl)$),

where the superscripts *n*, *l* denote the indices of different molecules. The form of $VLJ(nl)$ is

where *R*_{ij} (*i* ∈ *n* and *j* ∈ *l*) denotes the distance between atoms in molecules *n* and *l*. Here, *ϵ*_{CC} = 0.0559 kcal mol^{−1} (8.9126 × 10^{−5} a.u.), *σ*_{CC} = 2.800 Å (5.291 a.u.), *ϵ*_{OO} = 0.1597 kcal mol^{−1} (2.5454 × 10^{−4} a.u.), *σ*_{OO} = 3.028 Å (5.722 a.u.), $\u03f5CO=\u03f5CC\u03f5OO$, and *σ*_{CO} = (*σ*_{CC} + *σ*_{OO})/2. The form of $VCoul(nl)$ is

Here, *Q*_{C} = 0.6512 |*e*| and *Q*_{O} = −0.3256 |*e*| (where *e* denotes the charge of the electron).

The intramolecular interaction is characterized by

Here, the form of *V*_{CO} is defined in Eq. (6). *R*_{n1} and *R*_{n2} denote the lengths of two C=O bonds, and *θ*_{n} and *θ*_{eq} denote the O—C—O angle and the equilibrium angle. Here, *k*_{θ} = 108.0 kJ mol^{−1} rad^{−2} (0.0861 a.u.) and *θ*_{eq} = 180°.

Given the CO_{2} force field defined above, one can easily calculate the cavity-free force $Fnj(0)$ as a function of the nuclear configurations by standard molecular dynamics packages. The molecular dipole moment projected along the direction *ξ*_{λ} = **e**_{x}, **e**_{y} (*d*_{ng,λ}) is given by

and the derivative *∂d*_{ng,λ}/*∂***R**_{nj} is straightforward. In calculating the IR spectrum, the total dipole moment *μ*_{S} projected along the direction *ξ*_{λ} is given by $\mu S\u22c5\xi \lambda =\u2211n=1Nsubdng,\lambda $.

## REFERENCES

As a practical matter, even though the raw value of $q\u0303^k,\lambda $ is different from the raw value of $q^k,\lambda $, the spectrum and the energy of the photon can be calculated with either $q^k,\lambda $ or $q\u0303^k,\lambda $ and the same results can be obtained.

Although, in principle, an external pulse can also directly excite cavity photons, this feature is not included in Eq. (2). Including such an excitation would not qualitatively change the simulation results in this manuscript since pumping either the molecular or the photonic part equivalently pumps polaritons. More importantly, in order to include this feature, one would need to introduce a new phenomenological term, namely, the effective transition dipole of cavity photons, the magnitude of which varies for different cavities. Therefore, introducing this term would bring an additional manipulatable parameter and would hinder the universality of our simulation results: in Sec. IV D, we will show that the presented results are universal for cavities with different volumes.

_{2}

Note that consistent polaritonic relaxation dynamics can also be captured by evaluating the square norm of the molecular total dipole moment ($|\mu S(t)|2\u2212|\mu S(t)|2$).

Note that this weak enhancement of nonlinearity might be responsible for the relatively short LP lifetime as found in Fig. 3(a).