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.

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 equilibria8 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 rates11 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 calculations23–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 approach22 captures the asymmetric Rabi splitting6 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 H2O 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 modes18 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 CO2 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

FIG. 1.

Sketch of the simulation setup where a large collection of CO2 molecules is confined between a pair of metallic mirrors; see Sec. III for details. As shown in the left cartoon, the main finding of this study is that strongly exciting the LP in a cavity can greatly enhance the multiphoton nonlinear absorption of molecules relative to that outside the cavity. This enhancement is largest when twice the LP frequency approximately matches the vibrational 0 → 2 transition. After one directly excites the LP, and indirectly excites a localized vibrational state with v = 2, there is subsequently a gradual transfer of population to the first vibrationally excited state with a timescale that is much slower than the LP lifetime.

FIG. 1.

Sketch of the simulation setup where a large collection of CO2 molecules is confined between a pair of metallic mirrors; see Sec. III for details. As shown in the left cartoon, the main finding of this study is that strongly exciting the LP in a cavity can greatly enhance the multiphoton nonlinear absorption of molecules relative to that outside the cavity. This enhancement is largest when twice the LP frequency approximately matches the vibrational 0 → 2 transition. After one directly excites the LP, and indirectly excites a localized vibrational state with v = 2, there is subsequently a gradual transfer of population to the first vibrationally excited state with a timescale that is much slower than the LP lifetime.

Close modal

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,

ĤQEDG=ĤMG+ĤFG.
(1a)

Here, ĤMG denotes the conventional ground-state molecular Hamiltonian,

ĤMG=n=1NjnP^nj22Mnj+V^g(n)({R^nj})+n=1Nl>nV^inter(nl),
(1b)

where P^nj, R^nj, and Mnj denote the momentum operator, position operator, and mass for the jth 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. ĤFG denotes the field-related Hamiltonian,14,30,31

ĤFG=k,λp̃^k,λ22mk,λ+12mk,λωk,λ2q̃^k,λ+n=1Nd^ng,λωk,λΩϵ0mk,λ2,
(1c)

where p̃^k,λ, q̃^k,λ, ωk,λ, and mk,λ 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,λmk,λq̃^k,λ, where q^k,λ is the standard mass-reduced photonic position operator, the final dynamics of the physical q^k,λ 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,λ denotes the electronic ground-state dipole operator for the molecule n projected along the direction of ξλ. In Eq. (1c), we can define

εk,λmk,λωk,λ2/Ωϵ0
(1d)

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,λ. 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,

  MnjR̈nj=Fnj(0)+Fnjcav+Fnjext(t),
(2a)
mk,λq̈k,λ=mk,λωk,λ2qk,λε̃k,λn=1Nsubdng,λ.
(2b)

Here, the subscript nj denotes the jth atom in the molecule n; Fnj(0)=Vg(n)/RnjlnVinter(nl)/Rnj denotes the molecular part of the force on each nucleus;

Fnjcav=k,λε̃k,λqk,λ+ε̃k,λ2mk,λωk,λ2l=1Nsubdlg,λdng,λRnj

denotes the cavity force on each nucleus; and Fnjext(t)=QnjEext(t) denotes an external driving force acting on each nucleus with partial charge Qnj under the pumping of a time-dependent electric field Eext(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 qk,λ=q̃k,λ/Ncell, where Ncell denotes the number of periodic simulation cells for molecules, and also an effective coupling strength

ε̃k,λ=Ncellεk,λ,
(2c)

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 Nsub, the total number of molecules is N = NsubNcell. Note that when CavMD is used to reproduce experimental observations such as polariton relaxation, given the value of Nsub, ε̃k,λNcell can be chosen by fitting the experimentally observed Rabi splitting ΩN; since ΩNN, ε̃k,λNcellΩ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 Nsub) to make sure that any observed dynamics are not an artifact of the simulation; after all, the effective coupling strength for molecules ε̃k,λ 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 Nsub, all while keeping the Rabi splitting ΩN fixed and adjusting ε̃k,λ accordingly. For such calculations, the asymptotic results when Nsub 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 Ncell = 1 (so that ε̃k,λ=εk,λ 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).

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 

n(ω)α(ω)=πβω22ϵ0Vc12π+dteiωti=x,yμS(0)eiμS(t)ei,
(3)

while the latter is defined by

n(ω)αlocal(ω)=πβω22ϵ0Vc12π+dteiωt1Nsubn=1Nsubμn(0)μn(t).
(4)

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), ei 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) = ∑njQnjRnj(t). Note that, in our force field calculations, the partial charges (Qnj, 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) = ∑jQnjRnj(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 Nsub = 216 is taken during the simulation, the bright state contribution to the local IR spectrum in Eq. (4) is negligibly small, equal to 1/Nsub = 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, ÂB^(t)ÂB^(t). 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.

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., SiO2, 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 mk,λ = 1 a.u. (atomic units) though, as mentioned above, this mass does not affect any dynamics. Under periodic boundary conditions, we simulate 216 CO2 molecules in a cubic cell with cell length 24.292 Å (45.905 a.u.); the density of the liquid CO2 is 1.101 g/cm3. Unless stated otherwise, by default, we set the cavity mode frequency as 2320 cm−1 with an effective coupling strength ε̃=2×104 a.u. Compared with VSC experiments for which usually N = 109–1011 molecules are involved, our choice of the effective coupling strength ε̃Ncell should correspond to the involvement of Ncell = 107–109 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:

Eext(t)=E0cosωt+ϕex,
(5)

where the phase ϕ ∈ [0, 2π) is set as random. This pulse is turned on at tstart = 0.1 ps and is turned off at tend = 0.6 ps. Below, we will show results obtained after weak pumping, i.e., E0 = 3.084 × 106 V/m (6 × 10−4 a.u.) and the input pulse fluence F=12ϵ0cE02(tendtstart)=6.32 mJ/cm2, and also after strong pumping, i.e., E0 = 3.084 × 107 V/m (6 × 10−3 a.u.) and F = 632 mJ/cm2. 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 CO2 (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:

VCO(r)=Drαr2Δr2αr3Δr3+712αr4Δr4,
(6)

where Δr = rreq. Equation (6) is a fourth-order Taylor expansion of a Morse potential VM(r)=Dr[1exp(αrΔr)]2. The parameters are taken as follows: req = 1.162 Å (2.196 a.u.), Dr = 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 Dr 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. (104 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 

In Fig. 2(a), we plot the IR spectrum obtained from Eq. (3) for different values of the effective coupling strength ε̃; the value of ε̃ (in a.u.) is labeled on each lineshape, and the case of molecular system outside the cavity corresponds to ε̃=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 ε̃ scales as ε̃N with the total number of molecules N; see Sec. II for details. Therefore, increasing ε̃ 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 ε̃. The inset plots the Rabi frequency, or the frequency difference between the UP and LP peaks, as a function of ε̃N. 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.

FIG. 2.

Rabi splitting from equilibrium simulations. (a) IR spectrum as a function of the increased (from bottom to top) effective coupling strength ε̃ when the cavity mode frequency (2320 cm−1) is nearly resonant with the C=O asymmetric stretch (2327 cm−1). Inset: Linear relationship between Rabi splitting (ΩN) and ε̃N. (b) Avoided crossing of the IR spectrum as a function of the cavity mode frequency given ε̃=2×104 a.u. The dashed white line denotes the asymmetric C=O stretch outside a cavity, the dashed green line denotes the cavity mode, and the intensities of IR spectra are plotted on a logarithmic scale; see the colorbar. Note that the IR spectra here are calculated from equilibrium trajectories by evaluating the dipole auto-correlation function in Eq. (3); see Sec. III for other simulation details.

FIG. 2.

Rabi splitting from equilibrium simulations. (a) IR spectrum as a function of the increased (from bottom to top) effective coupling strength ε̃ when the cavity mode frequency (2320 cm−1) is nearly resonant with the C=O asymmetric stretch (2327 cm−1). Inset: Linear relationship between Rabi splitting (ΩN) and ε̃N. (b) Avoided crossing of the IR spectrum as a function of the cavity mode frequency given ε̃=2×104 a.u. The dashed white line denotes the asymmetric C=O stretch outside a cavity, the dashed green line denotes the cavity mode, and the intensities of IR spectra are plotted on a logarithmic scale; see the colorbar. Note that the IR spectra here are calculated from equilibrium trajectories by evaluating the dipole auto-correlation function in Eq. (3); see Sec. III for other simulation details.

Close modal

Figure 2(b) plots Rabi splitting as a function of the cavity mode frequency for the condition ε̃=2×104 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.

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 (λ=x,ymk,λωk,λqk,λ2/2+pk,λ2/2mk,λ) 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 ε̃=2×104 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) = E0 cos(ωt + ϕ)ex (with fluence F = 6.32 mJ/cm2), 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/cm2), 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.

FIG. 3.

Time-resolved dynamics for the cavity photon energy after resonantly exciting the UP (magenta; peaked at 2428 cm−1) or LP (cyan; peaked at 2241 cm−1) with a (a) weak or (b) strong pulse E(t) = E0 cos(ωt + ϕ)ex of 0.5 ps duration. The pulse fluence is F = 6.32 mJ/cm2 or F = 632 mJ/cm2, respectively. The cavity mode frequency is ωc = 2320 cm−1, and the effective coupling strength is ε̃=2×104 a.u.; see Sec. III for other details. The polariton lifetimes (τLP and τUP labeled at each subplot) are obtained by an exponential fit of the decaying energy of the cavity photon after the pulse.

FIG. 3.

Time-resolved dynamics for the cavity photon energy after resonantly exciting the UP (magenta; peaked at 2428 cm−1) or LP (cyan; peaked at 2241 cm−1) with a (a) weak or (b) strong pulse E(t) = E0 cos(ωt + ϕ)ex of 0.5 ps duration. The pulse fluence is F = 6.32 mJ/cm2 or F = 632 mJ/cm2, respectively. The cavity mode frequency is ωc = 2320 cm−1, and the effective coupling strength is ε̃=2×104 a.u.; see Sec. III for other details. The polariton lifetimes (τLP and τUP labeled at each subplot) are obtained by an exponential fit of the decaying energy of the cavity photon after the pulse.

Close modal

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.

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 Ti is calculated by evaluating Eq. (3) in a time window [Ti, Ti + Δ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 cm1

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/cm2). 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 CO2 force field). In particular, for the anharmonic potential that we use to simulate CO2 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.

FIG. 4.

Time-resolved spectra after exciting the LP with a strong pulse (F = 632 mJ/cm2). Three cases are compared: (a) and (d) exciting the LP (2241 cm−1) inside a 2320 cm−1 cavity; (b) and (e) exciting the LP (2167 cm−1) inside a 2200 cm−1 cavity; (c) and (f) off-resonant excitation at 2241 cm−1 outside the cavity. The corresponding incident exciting frequency is also labeled by a red arrow in each subplot, and (c) and (f) have been multiplied by a factor of four for better visualization. Here, the IR spectra [top panel, evaluated with Eq. (3)] reflect information about the molecular collective bright state, while the local IR spectra [bottom panel, evaluated with Eq. (4)] reflect mostly information about the molecular dark modes. At every time snapshot Ti, the IR or local IR spectrum is calculated by averaging over the time period [Ti, Ti + ΔT] with Eq. (3) or (4), where ΔT = 5 ps. See Sec. III for other simulation details. To better distinguish between linear and nonlinear absorption, the region of the linear absorption is labeled with blue horizontal lines (i.e., from 2220 cm−1 to 2360 cm−1) and the region of the nonlinear absorption is labeled with cyan horizontal lines (i.e., from 2150 cm−1 to 2220 cm−1). Note that, inside the cavity, when exciting the LP, the nonlinear absorption can be greatly enhanced than that outside the cavity; see Figs. 4(d) and 4(f) for comparison. After the pulse, the system temperature is increased from 300 K to 505 K, 366 K, 331 K (from left to right), respectively.

FIG. 4.

Time-resolved spectra after exciting the LP with a strong pulse (F = 632 mJ/cm2). Three cases are compared: (a) and (d) exciting the LP (2241 cm−1) inside a 2320 cm−1 cavity; (b) and (e) exciting the LP (2167 cm−1) inside a 2200 cm−1 cavity; (c) and (f) off-resonant excitation at 2241 cm−1 outside the cavity. The corresponding incident exciting frequency is also labeled by a red arrow in each subplot, and (c) and (f) have been multiplied by a factor of four for better visualization. Here, the IR spectra [top panel, evaluated with Eq. (3)] reflect information about the molecular collective bright state, while the local IR spectra [bottom panel, evaluated with Eq. (4)] reflect mostly information about the molecular dark modes. At every time snapshot Ti, the IR or local IR spectrum is calculated by averaging over the time period [Ti, Ti + ΔT] with Eq. (3) or (4), where ΔT = 5 ps. See Sec. III for other simulation details. To better distinguish between linear and nonlinear absorption, the region of the linear absorption is labeled with blue horizontal lines (i.e., from 2220 cm−1 to 2360 cm−1) and the region of the nonlinear absorption is labeled with cyan horizontal lines (i.e., from 2150 cm−1 to 2220 cm−1). Note that, inside the cavity, when exciting the LP, the nonlinear absorption can be greatly enhanced than that outside the cavity; see Figs. 4(d) and 4(f) for comparison. After the pulse, the system temperature is increased from 300 K to 505 K, 366 K, 331 K (from left to right), respectively.

Close modal

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 learned47 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 (Ti = 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.

FIG. 5.

The corresponding density distribution of the C=O bond potential energy per molecule at time Ti = 1 ps for Fig. 4. Three cases are compared: exciting the LP inside a (a) 2320 cm−1 or (b) 2200 cm−1 cavity and (c) exciting 2241 cm−1 [the LP frequency for (a)] outside the cavity. All parameters are the same as in Fig. 4. The C=O bond potential energy is calculated according to Eq. (6) and is shown in units of ℏω0 (ω0 = 2327 cm−1). Note that exciting the LP in a 2320 cm−1 cavity induces a meaningful fraction of vibrationally highly excited molecules, which corresponds to the nonlinear local IR peak ∼2170 cm−1 in Fig. 4(d).

FIG. 5.

The corresponding density distribution of the C=O bond potential energy per molecule at time Ti = 1 ps for Fig. 4. Three cases are compared: exciting the LP inside a (a) 2320 cm−1 or (b) 2200 cm−1 cavity and (c) exciting 2241 cm−1 [the LP frequency for (a)] outside the cavity. All parameters are the same as in Fig. 4. The C=O bond potential energy is calculated according to Eq. (6) and is shown in units of ℏω0 (ω0 = 2327 cm−1). Note that exciting the LP in a 2320 cm−1 cavity induces a meaningful fraction of vibrationally highly excited molecules, which corresponds to the nonlinear local IR peak ∼2170 cm−1 in Fig. 4(d).

Close modal

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 cavity48 and quantum dots near surface plasmons,49,50 while its possibility under collective VSC has not been extensively studied19,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/cm2), 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 CO2 molecules in a thermal distribution) and 2170 cm−1 (the frequency of the highly excited CO2 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 Ti = 1 ps. This early-time nonlinear signal (INL) 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), INL 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.

FIG. 6.

Integrated nonlinear absorption (INL) at early times plotted as a function of the fluence of the exciting pulse of frequency 2241 cm−1 corresponding to the LP frequency in Figs. 4(a) and 4(d). INL is calculated by integrating the local IR spectra over frequency for the nonlinear absorption [i.e., the cyan region in Fig. 4(d)] at time Ti = 1 ps. All other parameters are the same as in Figs. 3(a) and 3(b). Note that compared with the nonlinear absorption outside the cavity (black squares) under the same pulse, the nonlinear signal inside the cavity (blue circles) can be enhanced by up to two orders of magnitude by polariton-enhanced multiphoton absorption.

FIG. 6.

Integrated nonlinear absorption (INL) at early times plotted as a function of the fluence of the exciting pulse of frequency 2241 cm−1 corresponding to the LP frequency in Figs. 4(a) and 4(d). INL is calculated by integrating the local IR spectra over frequency for the nonlinear absorption [i.e., the cyan region in Fig. 4(d)] at time Ti = 1 ps. All other parameters are the same as in Figs. 3(a) and 3(b). Note that compared with the nonlinear absorption outside the cavity (black squares) under the same pulse, the nonlinear signal inside the cavity (blue circles) can be enhanced by up to two orders of magnitude by polariton-enhanced multiphoton absorption.

Close modal

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 ε̃=Ncellε is artificially amplified proportionally to the number of periodic simulation cells (Ncell) 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 (Nsub) and study the polariton lifetime dependence on Nsub. Note that fixing Rabi splitting (i.e., fixing N = NsubNcell) while increasing Nsub will lead to a decrease in ε̃ (ε̃Ncell=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/Nsub 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/cm2] or strong [Fig. 7(c); F = 632 mJ/cm2] 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/Nsub, suggesting that the reported simulation results are not sensitive to renormalization applied to the light–matter interaction (ε̃) 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.

FIG. 7.

Polariton lifetime against the number of molecules in the simulation cell (1/Nsub). (a) and (b) Lifetimes of the (a) LP (cyan circles) and (b) UP (magenta triangles) under weak excitation (F = 6.32 mJ/cm2). (c) and (d) Lifetimes of the (c) LP and (d) UP under strong excitation (F = 632 mJ/cm2). Dashed lines denote the linear fitting of each dataset. For parameters, when Nsub is increased from 216 to 1200, molecular density and the observed Rabi splitting are fixed the same by increasing the simulation cell size and decreasing the effective coupling strength ε̃ accordingly. All other simulation details are the same as in Fig. 3.

FIG. 7.

Polariton lifetime against the number of molecules in the simulation cell (1/Nsub). (a) and (b) Lifetimes of the (a) LP (cyan circles) and (b) UP (magenta triangles) under weak excitation (F = 6.32 mJ/cm2). (c) and (d) Lifetimes of the (c) LP and (d) UP under strong excitation (F = 632 mJ/cm2). Dashed lines denote the linear fitting of each dataset. For parameters, when Nsub is increased from 216 to 1200, molecular density and the observed Rabi splitting are fixed the same by increasing the simulation cell size and decreasing the effective coupling strength ε̃ accordingly. All other simulation details are the same as in Fig. 3.

Close modal

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.

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.

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.

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 

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

FIG. 8.

The C=O bond potential in our force field. The fourth-order anharmonic potential in Eq. (6) (red solid) is compared with the corresponding Morse potential (black dashed-dotted) and the harmonic approximation (green dashed).

FIG. 8.

The C=O bond potential in our force field. The fourth-order anharmonic potential in Eq. (6) (red solid) is compared with the corresponding Morse potential (black dashed-dotted) and the harmonic approximation (green dashed).

Close modal

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., nn + 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 CO2 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).

FIG. 9.

The C=O asymmetric peak frequency as a function of the potential energy of a single CO2 molecule in the gas phase. The potential energy is plotted in units of ℏω0, where ω0 = 2342 cm−1 denotes the peak frequency in thermal equilibrium at 300 K. See the text for simulation details.

FIG. 9.

The C=O asymmetric peak frequency as a function of the potential energy of a single CO2 molecule in the gas phase. The potential energy is plotted in units of ℏω0, where ω0 = 2342 cm−1 denotes the peak frequency in thermal equilibrium at 300 K. See the text for simulation details.

Close modal

For completeness, below we will provide the details of the CO2 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 CO2 force field, the pairwise intermolecular potential is characterized by the Lennard-Jones potential (VLJ(nl)) plus the Coulombic interactions between atoms (VCoul(nl)),

Vinter(nl)=VLJ(nl)+VCoul(nl),
(C1)

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

VLJ(nl)=injl4ϵijσijRij12σijRij6,
(C2)

where Rij (in and jl) 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.), ϵCO=ϵCCϵOO, and σCO = (σCC + σOO)/2. The form of VCoul(nl) is

VCoul(nl)=injlQiQj4πϵ0Rij.
(C3)

Here, QC = 0.6512 |e| and QO = −0.3256 |e| (where e denotes the charge of the electron).

The intramolecular interaction is characterized by

Vg(n)=VCO(Rn1)+VCO(Rn2)+12kθθnθeq2.
(C4)

Here, the form of VCO is defined in Eq. (6). Rn1 and Rn2 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 CO2 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 ξλ = ex, ey (dng,λ) is given by

dng,λ=QORnO1+RnO2+QCRnCξλ,
(C5)

and the derivative ∂dng,λ/Rnj is straightforward. In calculating the IR spectrum, the total dipole moment μS projected along the direction ξλ is given by μSξλ=n=1Nsubdng,λ.

1.
F.
Herrera
and
J.
Owrutsky
, “
Molecular polaritons for controlling chemistry with quantum optics
,”
J. Chem. Phys.
152
,
100902
(
2020
).
2.
A.
Shalabney
,
J.
George
,
J.
Hutchison
,
G.
Pupillo
,
C.
Genet
, and
T. W.
Ebbesen
, “
Coherent coupling of molecular resonators with a microcavity mode
,”
Nat. Commun.
6
,
5981
(
2015
).
3.
J.
George
,
A.
Shalabney
,
J. A.
Hutchison
,
C.
Genet
, and
T. W.
Ebbesen
, “
Liquid-phase vibrational strong coupling
,”
J. Phys. Chem. Lett.
6
,
1027
1031
(
2015
).
4.
J.
George
,
T.
Chervy
,
A.
Shalabney
,
E.
Devaux
,
H.
Hiura
,
C.
Genet
, and
T. W.
Ebbesen
, “
Multiple Rabi splittings under ultrastrong vibrational coupling
,”
Phys. Rev. Lett.
117
,
153601
(
2016
).
5.
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
, “
Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field
,”
Angew. Chem., Int. Ed.
55
,
11462
11466
(
2016
).
6.
R. M. A.
Vergauwe
,
A.
Thomas
,
K.
Nagarajan
,
A.
Shalabney
,
J.
George
,
T.
Chervy
,
M.
Seidel
,
E.
Devaux
,
V.
Torbeev
, and
T. W.
Ebbesen
, “
Modification of enzyme activity by vibrational strong coupling of water
,”
Angew. Chem., Int. Ed.
58
,
15324
15328
(
2019
).
7.
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
, “
Tilting a ground-state reactivity landscape by vibrational strong coupling
,”
Science
363
,
615
619
(
2019
).
8.
Y.
Pang
,
A.
Thomas
,
K.
Nagarajan
,
R. M. A.
Vergauwe
,
K.
Joseph
,
B.
Patrahau
,
K.
Wang
,
C.
Genet
, and
T. W.
Ebbesen
, “
On the role of symmetry in vibrational strong coupling: The case of charge-transfer complexation
,”
Angew. Chem., Int. Ed.
59
,
10436
10440
(
2020
).
9.
B.
Xiang
,
R. F.
Ribeiro
,
Y.
Li
,
A. D.
Dunkelberger
,
B. B.
Simpkins
,
J.
Yuen-Zhou
, and
W.
Xiong
, “
Manipulating optical nonlinearities of molecular polaritons by delocalization
,”
Sci. Adv.
5
,
eaax5196
(
2019
); arXiv:1901.05526.
10.
R. F.
Ribeiro
,
A. D.
Dunkelberger
,
B.
Xiang
,
W.
Xiong
,
B. S.
Simpkins
,
J. C.
Owrutsky
, and
J.
Yuen-Zhou
, “
Theory for nonlinear spectroscopy of vibrational polaritons
,”
J. Phys. Chem. Lett.
9
,
3766
3771
(
2018
).
11.
B.
Xiang
,
R. F.
Ribeiro
,
M.
Du
,
L.
Chen
,
Z.
Yang
,
J.
Wang
,
J.
Yuen-Zhou
, and
W.
Xiong
, “
Intermolecular vibrational energy transfer enabled by microcavity strong light–matter coupling
,”
Science
368
,
665
667
(
2020
).
12.
S.
Rudin
and
T. L.
Reinecke
, “
Oscillator model for vacuum Rabi splitting in microcavities
,”
Phys. Rev. B
59
,
10227
10233
(
1999
).
13.
J.
Galego
,
C.
Climent
,
F. J.
Garcia-Vidal
, and
J.
Feist
, “
Cavity Casimir-Polder forces and their effects in ground-state chemical reactivity
,”
Phys. Rev. X
9
,
021057
(
2019
).
14.
T. E.
Li
,
A.
Nitzan
, and
J. E.
Subotnik
, “
On the origin of ground-state vacuum-field catalysis: Equilibrium consideration
,”
J. Chem. Phys.
152
,
234107
(
2020
); arXiv:2002.09977.
15.
J. A.
Campos-Gonzalez-Angulo
and
J.
Yuen-Zhou
, “
Polaritonic normal modes in transition state theory
,”
J. Chem. Phys.
152
,
161101
(
2020
).
16.
V. P.
Zhdanov
, “
Vacuum field in a cavity, light-mediated vibrational coupling, and chemical reactivity
,”
Chem. Phys.
535
,
110767
(
2020
).
17.
G. D.
Scholes
,
C. A.
DelPo
, and
B.
Kudisch
, “
Entropy reorders polariton states
,”
J. Phys. Chem. Lett.
11
,
6389
6395
(
2020
).
18.
B.
Xiang
,
R. F.
Ribeiro
,
A. D.
Dunkelberger
,
J.
Wang
,
Y.
Li
,
B. S.
Simpkins
,
J. C.
Owrutsky
,
J.
Yuen-Zhou
, and
W.
Xiong
, “
Two-dimensional infrared spectroscopy of vibrational polaritons
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
4845
4850
(
2018
).
19.
B.
Xiang
,
R. F.
Ribeiro
,
L.
Chen
,
J.
Wang
,
M.
Du
,
J.
Yuen-Zhou
, and
W.
Xiong
, “
State-selective polariton to dark state relaxation dynamics
,”
J. Phys. Chem. A
123
,
5918
5927
(
2019
).
20.
A. B.
Grafton
,
A. D.
Dunkelberger
,
B. S.
Simpkins
,
J. F.
Triana
,
F. J.
Hernandez
,
F.
Herrera
, and
J.
Owrutsky
, “
Excited-state vibration-polariton transitions and dynamics in nitroprusside
,”
Nat. Commun.
12
,
214
(
2021
).
21.
R. F.
Ribeiro
,
J. A.
Campos-Gonzalez-Angulo
,
N. C.
Giebink
,
W.
Xiong
, and
J.
Yuen-Zhou
, “
Enhanced optical nonlinearities under strong light-matter coupling
,” arXiv:2006.08519 (
2020
).
22.
T. E.
Li
,
J. E.
Subotnik
, and
A.
Nitzan
, “
Cavity molecular dynamics simulations of liquid water under vibrational ultrastrong coupling
,”
Proc. Natl. Acad. Sci. U. S. A.
117
,
18324
18331
(
2020
); arXiv:2004.04888.
23.
F.
Herrera
and
F. C.
Spano
, “
Cavity-controlled chemistry in molecular ensembles
,”
Phys. Rev. Lett.
116
,
238301
(
2016
).
24.
J.
Flick
,
M.
Ruggenthaler
,
H.
Appel
, and
A.
Rubio
, “
Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
3026
3034
(
2017
).
25.
H. L.
Luk
,
J.
Feist
,
J. J.
Toppari
, and
G.
Groenhof
, “
Multiscale molecular dynamics simulations of polaritonic chemistry
,”
J. Chem. Theory Comput.
13
,
4324
4335
(
2017
).
26.
G.
Groenhof
,
C.
Climent
,
J.
Feist
,
D.
Morozov
, and
J. J.
Toppari
, “
Tracking polariton relaxation with multiscale molecular dynamics simulations
,”
J. Phys. Chem. Lett.
10
,
5476
5483
(
2019
).
27.
C.
Schäfer
,
M.
Ruggenthaler
,
V.
Rokaj
, and
A.
Rubio
, “
Relevance of the quadratic diamagnetic and self-polarization terms in cavity quantum electrodynamics
,”
ACS Photonics
7
,
975
990
(
2020
).
28.
A.
Mandal
,
S.
Montillo Vega
, and
P.
Huo
, “
Polarized Fock states and the dynamical Casimir effect in molecular cavity quantum electrodynamics
,”
J. Phys. Chem. Lett.
11
,
9215
9223
(
2020
).
29.
A. D.
Dunkelberger
,
R. B.
Davidson
,
W.
Ahn
,
B. S.
Simpkins
, and
J. C.
Owrutsky
, “
Ultrafast transmission modulation and recovery via vibrational strong coupling
,”
J. Phys. Chem. A
122
,
965
971
(
2018
).
30.
C.
Cohen-Tannoudji
,
J.
Dupont-Roc
, and
G.
Grynberg
,
Photons and Atoms: Introduction to Quantum Electrodynamics
(
Wiley
,
New York
,
1997
), pp.
280
295
.
31.
T. S.
Haugland
,
E.
Ronca
,
E. F.
Kjønstad
,
A.
Rubio
, and
H.
Koch
, “
Coupled cluster theory for molecular polaritons: Changing ground and excited states
,”
Phys. Rev. X
10
,
041043
(
2020
).
32.

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

33.
J.
Feist
,
A. I.
Fernández-Domínguez
, and
F. J.
García-Vidal
, “
Macroscopic QED for quantum nanophotonics: Emitter-centered modes as a minimal basis for multiemitter problems
,”
Nanophotonics
10
,
477
489
(
2020
); arXiv:2008.02106.
34.

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.

35.
D. A.
McQuarrie
,
Statistical Mechanics
(
HarperCollins Publishers
,
New York
,
1976
).
36.
M.-P.
Gaigeot
and
M.
Sprik
, “
Ab initio molecular dynamics computation of the infrared spectrum of aqueous uracil
,”
J. Phys. Chem. B
107
,
10344
10358
(
2003
).
37.
S.
Habershon
,
G. S.
Fanourgakis
, and
D. E.
Manolopoulos
, “
Comparison of path integral molecular dynamics methods for the infrared absorption spectrum of liquid water
,”
J. Chem. Phys.
129
,
074501
(
2008
).
38.
A.
Nitzan
,
Chemical Dynamics in Condensed Phases: Relaxation, Transfer and Reactions in Condensed Molecular Systems
(
Oxford University Press
,
New York
,
2006
).
39.
R. T.
Cygan
,
V. N.
Romanov
, and
E. M.
Myshakin
, “
Molecular simulation of carbon dioxide capture by montmorillonite using an accurate and flexible force field
,”
J. Phys. Chem. C
116
,
13079
13091
(
2012
).
40.
B.
de Darwent
,
Bond Dissociation Energies in Simple Molecules
(
U.S. National Bureau of Standards
,
1970
).
41.
L.
Martínez
,
R.
Andrade
,
E. G.
Birgin
, and
J. M.
Martínez
, “
PACKMOL: A package for building initial configurations for molecular dynamics simulations
,”
J. Comput. Chem.
30
,
2157
2164
(
2009
).
42.
V.
Kapil
,
M.
Rossi
,
O.
Marsalek
,
R.
Petraglia
,
Y.
Litman
,
T.
Spura
,
B.
Cheng
,
A.
Cuzzocrea
,
R. H.
Meißner
,
D. M.
Wilkins
,
B. A.
Helfrecht
,
P.
Juda
,
S. P.
Bienvenue
,
W.
Fang
,
J.
Kessler
,
I.
Poltavsky
,
S.
Vandenbrande
,
J.
Wieme
,
C.
Corminboeuf
,
T. D.
Kühne
,
D. E.
Manolopoulos
,
T. E.
Markland
,
J. O.
Richardson
,
A.
Tkatchenko
,
G. A.
Tribello
,
V.
Van Speybroeck
, and
M.
Ceriotti
, “
i-PI 2.0: A universal force engine for advanced molecular simulations
,”
Comput. Phys. Commun.
236
,
214
223
(
2019
).
43.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
44.
T. E.
Li
, “
Cavity molecular dynamics simulations tool sets
,” https://github.com/TaoELi/cavity-md-ipi,
2020
.
45.
T.
Seki
,
J.-D.
Grunwaldt
, and
A.
Baiker
, “
In situ attenuated total reflection infrared spectroscopy of imidazolium-based room-temperature ionic liquids under ‘supercritical’ CO2
,”
J. Phys. Chem. B
113
,
114
122
(
2009
).
46.

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

47.
W.
Xiong
, private communication (
2020
).
48.
K.
Wang
,
M.
Seidel
,
K.
Nagarajan
,
T.
Chervy
,
C.
Genet
, and
T. W.
Ebbesen
, “
Large optical nonlinearity enhancement under electronic strong coupling
,” arXiv:2005.13325 (
2020
).
49.
M.
Fu
,
K.
Wang
,
H.
Long
,
G.
Yang
,
P.
Lu
,
F.
Hetsch
,
A. S.
Susha
, and
A. L.
Rogach
, “
Resonantly enhanced optical nonlinearity in hybrid semiconductor quantum dot–metal nanoparticle structures
,”
Appl. Phys. Lett.
100
,
063117
(
2012
).
50.
N.
Rivera
,
G.
Rosolen
,
J. D.
Joannopoulos
,
I.
Kaminer
, and
M.
Soljačić
, “
Making two-photon processes dominate one-photon processes using mid-IR phonon polaritons
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
13607
13612
(
2017
).
51.
S.
Scandolo
and
F.
Bassani
, “
Nonlinear sum rules: The three-level and the anharmonic-oscillator models
,”
Phys. Rev. B
45
,
13257
13261
(
1992
).
52.

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

53.
M.
Gross
and
S.
Haroche
, “
Superradiance: An essay on the theory of collective spontaneous emission
,”
Phys. Rep.
93
,
301
396
(
1982
).
54.
T. E.
Li
,
A.
Nitzan
,
M.
Sukharev
,
T.
Martinez
,
H.-T.
Chen
, and
J. E.
Subotnik
, “
Mixed quantum–classical electrodynamics: Understanding spontaneous decay and zero–point energy
,”
Phys. Rev. A
97
,
032105
(
2018
).
55.
M. E.
Goggin
and
P. W.
Milonni
, “
Driven Morse oscillator: Classical chaos, quantum theory, and photodissociation
,”
Phys. Rev. A
37
,
796
806
(
1988
).