The shape of the electric dipole function determines the sub-picosecond dynamics of anharmonic vibrational polaritons

Recent experimental demonstrations of strong and ultrastrong light–matter interaction with molecules and molecular materials in infrared cavities^1–18 have stimulated intense theoretical efforts for understanding the microscopic properties of hybrid photon-vibration states from a quantum mechanical perspective.¹⁹ Motivated by pioneering measurements in liquid-phase Fabry–Perot cavities,^5,17,20 theoretical studies have suggested several potential mechanisms that enable the modification of chemical reactivity in the ground electronic state under conditions of vibrational strong coupling.^21–23 Another theoretical focus is the study of linear and nonlinear spectroscopic signals of infrared cavities under strong coupling.^19,24–27

Vibrational polaritons are the hybrid light–matter states that emerge in infrared cavities under strong coupling.²⁸ Several models with varying degrees of complexity have been used to study these systems. In one of the earliest approaches,²⁹ molecular vibrations were treated as two-level systems with an energy gap given by the fundamental vibration frequency. The vibrational qubits were simultaneously coupled to a quantized harmonic oscillator describing a single-mode cavity field. This approach corresponds to the Tavis–Cummings model of cavity quantum electrodynamics (QED),³⁰ which implies the rotating-wave approximation for light–matter coupling. In Ref. 26, the authors went beyond the two-level approximation and also treated intramolecular vibrations within an electronic state as quantum harmonic oscillators. Light–matter coupling was also extended to include the counter-rotating and self-energy terms that are commonly taken into account under conditions of ultrastrong coupling,^31,32 broadly defined as the regime in which the light–matter interaction energy is comparable with the vibrational and cavity frequencies.

The anharmonicity of molecular vibrations was first taken into account in Ref. 25, in an effort to model the pump–probe spectrum of strongly coupled infrared cavities.¹⁶ Here, the authors supplemented the quantum harmonic oscillator picture with a weak anharmonic energy correction. The latter is introduced to capture the spectral anharmonicity of real vibrations (e.g., energy levels are not equally spaced). Later, in Ref. 27, this perturbative approach was further extended to include the effect of electrostatic anharmonicity (e.g., the dipole moment function is not symmetric relative to displacements from equilibrium). In both works, light–matter coupling was treated within the rotating-wave approximation.

Perturbative anharmonicity models can describe the frequencies and transition strengths of the fundamental (0 → 1) and overtone vibrational transitions (0 → ν) to excited states with a low vibrational quantum number ν.³³ Building on the intuition gained from the study of vibrational ladder climbing in strong laser fields,³⁴ one expects that under conditions of strong coupling in a resonant cavity, vibrational anharmonicity models should also properly describe light-induced couplings among all the vibrational bound states in the spectrum.

The first consistent approach that takes into account the entire bound state spectrum and electric dipole function of an electronic potential energy curve was developed in Ref. 19. The so-called multi-level quantum Rabi (MLQR) model describes an individual anharmonic vibration in an infrared cavity. It can be derived from a multipolar formulation of light–matter interaction³⁵ by projecting the composite system Hamiltonian into a complete vibrational energy basis. Given the electronic potential curve and electric dipole function, the MLQR model can be used to understand both material and photonic properties of vibrational polaritons both in the strong and ultrastrong coupling regimes.

In Ref. 19, only molecules without an electric dipole moment at equilibrium were discussed. Moreover, only off-diagonal dipole matrix elements in the vibrational energy basis (transition dipoles) were taken into account in the light–matter interaction process. We now significantly expand the MLQR model to include the contribution of diagonal elements of the dipole matrix to the vibrational basis (permanent dipoles) and compare the resulting polariton physics of molecular vibrations that are polar at equilibrium (e.g., CO) with those that are non-polar (NP) at equilbrium (e.g., CO₂). We show that the extended MLQR model is equivalent to a numerically exact^36,37 formulation of the problem in the coordinate representation. This equivalence has been widely ignored in the vibrational polariton literature. We focus on the evolution of vibrational polariton wavepackets and the corresponding dynamics of simple material and photonic observables of experimental interest.

We perform a systematic comparison between polar and non-polar molecules in infrared cavities, correlating the entire shape of the electric dipole function along the nuclear coordinate with the resulting vibrational polariton dynamics. Among other results, we predict that for a specific class of polar molecules, the light–matter system can evolve from a diabatic product state with a definite number of vibrational and cavity excitations (possibly vacuum) into a polariton wavepacket with a mean intracavity photon number that could be reliably measured using the current detector technology.

We show that not only transition dipole moments between vibrational levels are relevant for light–matter interaction inside infrared cavities—as in free-space IR absorption spectroscopy—but also vibrationally averaged permanent dipole moments. Although permanent dipole moments do not change the vibrational state of a molecule, they displace the cavity field from the vacuum state into a polariton wavepacket with a finite number of photons over femtosecond time scales. The displaced cavity field can then strongly drive transitions between higher vibrational levels. This is analogous to laser-induced vibrational ladder climbing,³⁴ but the number of photons involved is many orders of magnitude smaller.

In this article, we discuss the theoretical foundations of our model (Sec. II) and the details of our methods (Sec. III). We then describe the results obtained for the static and dynamical properties of vibrational polaritons that emerge under the various material and photonic conditions considered (Sec. IV). We finally conclude with a discussion of the fundamental physical principles that support our computational results and propose a quantum state preparation method that could be used to test our predictions with infrared nanoantennas (Sec. V).

Following Refs. 36–38, we model an individual molecular vibration coupled to a single quantized electromagnetic mode in the electric dipole approximation of light–matter interaction. In the coordinate representation, the system Hamiltonian reads (in atomic units)

The terms in the first bracket describe the vibrational motion of interest, characterized by a potential energy curve (PEC) V(q) along a normal mode coordinate q with reduced mass μ. In this work, the PEC is assumed to have a single equilibrium configuration at q = q_e and a well-defined dissociation energy in free space. This is the typical behavior of stretching modes in diatomic and polyatomic molecules.³³ 

The second bracket in Eq. (1) describes the energy of a single cavity field mode with frequency ω_c and quadrature operator $x^$⁠. The last term corresponds to light–matter interaction in a multipolar form,³⁵ truncated to the electric dipole contribution. $E0$ is the square-root amplitude of the vacuum fluctuations at ω_c, and $d^(q)$ is the electric dipole function along the vibrational coordinate. Quadrupole contributions to light–matter coupling are ignored under the assumption that the electric field gradient is weak at the location of the molecular dipole. This simplified description is not expected to describe all cavity systems.³⁹ 

The evolution of an arbitrary intracavity light–matter state $ψ(t)$ is obtained from $H^$ in Eq. (1) by solving $i(d/dt)ψ(t)=H^ψ(t)$ with a unitary propagator. As described in more details below, we solve the Schrödinger equation by representing the system wavefunction and the Hamiltonian both in the coordinate and Hilbert space representations. In the coordinate representation, we discretize both $H^$ and $ψ(t)$ along both nuclear and cavity coordinates (q, x) and propagate an initial state using the multi-configurational time-dependent Hartree (MCTDH) method.^40,41 In the Hilbert space representation, time evolution is carried out by projecting Eq. (1) into a set of energy eigenstates of the nuclear potential V(q) to give

where $Π^vib=∑ν=0νmaxνν$ is a projector operator into the vibrational energy basis $ν$⁠, with ν being the vibrational quantum number. If the energy basis set is complete, i.e., ν_max = ∞, then we have $H^′≡H^$⁠, and evolving a system state in the coordinate or Hilbert space representation must give exactly equivalent results. In practice, the projection operator $Π^vib$ can only reliably be truncated up to a cutoff energy by setting ν_max finite. We demonstrate below that despite this practical limitation, it is possible to find a value of ν_max that gives numerically equivalent results for the propagation of vibrational polariton wavepackets in the coordinate and the Hilbert space representations.

We are interested in the individual roles of diagonal and off-diagonal elements of the electric dipole matrix $ν′d(q)ν$ in the dynamics of vibrational polaritons. In order to have controllable access to this information, we partition the vibrationally projected Hamiltonian $H′^$ in Eq. (2) as $H′^=H1^+H2^$⁠, where

and

with summations over ν carried out up to ν_max.

$H1^$ in Eq. (3) corresponds to the multi-level quantum Rabi model from Ref. 19. In this model, light–matter coupling parameters are determined only by transition dipole moments through the state-dependent Rabi frequencies,

for ν ≠ ν′. Cavity field variables are described in terms of cavity annihilation operator â. For a two-level system, i.e., ν_max = 1, Eq. (3) reduces to the quantum Rabi model for a qubit in ultrastrong coupling.^42–44 

Noting that the intracavity electric field $E^$ is proportional to $E0(â†+â)$⁠,⁴⁵ the term $H^2$ in Eq. (4) can be interpreted as the contribution to d ⋅E from the permanent dipole moment of each vibrational level ν. The corresponding diagonal coupling strength is given by

The dependence on the vibrational quantum numbers of the diagonal (g_ν) and off-diagonal (g_ν′ν) coupling parameters is directly related to the potential energy curve V(q). Consider a potential that has even parity relative to the equilibrium mode length q_e, i.e., V(q) is invariant under the transformation q → −q and q_e = −q_e. This is the case for the harmonic potential $V(q)=ω0(q−qe)2/2$⁠, whose vibrational eigenstates are also eigenstates of parity.

Selection rules for electric dipole matrix elements can be derived by expanding the dipole function d(q) near q_e up to second order as

where d_e is the electric dipole moment at the equilibrium configuration q_e, c₁ is proportional to the slope of the dipole function at equilibrium, and c₂ to its curvature. The absolute magnitudes |c₁| and |c₂| can be inferred from the strengths of fundamental and overtone absorption peaks in stationary infrared spectroscopy.⁴⁶ Since the absorption line strengths are proportional to the square of the transition dipole moments, the signs of the expansion parameters in Eq. (7) are, therefore, not accessible in linear spectroscopy. For harmonic oscillators, the dipole expansion in Eq. (7) gives the selection rules Δν = ±1, ±2 for transition dipole moments. Diagonal elements (permanent dipoles) are only weakly dependent on the vibrational quantum number ν through the quadratic term in the expansion and are, thus, primarily given by $⟨ν|d^(q)|ν⟩≈de$ for small ν.

For realistic anharmonic molecular vibrations, the nuclear potential V(q) in general is not invariant under the transformation q → −q and q_e → −q_e. Anharmonic vibrational eigenstates $ν$⁠, thus, do not have a well-defined parity. This changes the structure of the dipole matrix elements $⟨ν′|d^(q)|ν⟩$ in comparison with the case of harmonic vibrations. In general, for anharmonic vibrations, there are no selection rules for transition dipole moments, and permanent dipole moments have a stronger dependence on the vibrational quantum number in comparison with harmonic vibrations. This occurs because the contribution of the linear term in Eq. (7) is not parity-forbidden, giving $⟨ν|d^(q)|ν⟩≈de+c1⟨ν|(q−qe)|ν⟩$ for small ν.

In this work, we consider the intracavity dynamics of an individual anharmonic vibration described by the Morse potential,⁴⁷ 

where D_e is the potential depth that defines the dissociation energy and a is a parameter that contributes to the anharmonicity of the cavity-free vibrational spectrum. The nuclear Schrödinger equation with a Morse potential can be solved analytically in terms of associated Laguerre polynomials.^46,47 By comparing the exact expression for the vibrational energies E_ν with the Dunham expansion,³³ the vibrational spectrum up to second order in ν can be written as

where $ω0≡[ℏ]a2De/μ$ is the fundamental vibration frequency in the harmonic approximation and χ_e ≡ [ℏ²]a²/2μ is the spectral anharmonicity parameter.

Having defined a model for the nuclear potential and a structure for the light–matter coupling Hamiltonian both in the coordinate space [Eq. (1)] and Hilbert space [Eqs. (3) and (4)], we compute the dynamics of selected intracavity light–matter wavepackets. Vibrational wavepackets can in general be written as

where $βνnc(t)$ are time-dependent wavepacket coefficients in the diabatic product basis $νnc$⁠, where $nc$ is an eigenstate of the photon number operator â^†â with eigenvalue n_c (Fock state). The coordinate-space analog of Eq. (10) can also be defined.

In order to gain physical intuition about the fate of the nuclear motion and the cavity field under conditions of strong and ultrastrong light–matter coupling, we focus on the short-time dynamics of the mean mode length $⟨ψ(t)|q^|ψ(t)⟩$ and the mean intracavity photon number ⟨ψ(t)|â^†â|ψ(t)⟩ as functions of the light–matter coupling strength and the shape of the dipole function d(q). We devote special attention to the latter, as one would expect that for strong enough light–matter interaction, the spectral observables should depend not only on the magnitudes of the dipole expansion parameters in Eq. (7)—as is the case in cavity-free infrared absorption spectroscopy—but also on their signs.

Radiative and non-radiative dissipative processes are not taken into account in this work. Therefore, our results can only accurately describe the unitary dynamics of vibrational polariton wavepackets over sub-picosecond time scales, before dissipative effects become relevant.¹³ 

We consider a model anharmonic oscillator described by a Morse potential [Eq. (8)] with parameters in atomic units given by D_e = 0.23 a.u., q_e = 4.0 a.u., and α = 1.4465 a.u. The reduced mass of the vibrational mode is μ = 1.43764 amu. The same Morse potential V(q) is used for the system Hamiltonian represented in coordinate space and Hilbert space. For our chosen parameters, the potential has 24 vibrational bound states, and the fundamental vibration period is 2π/ω₁₀ = 8.27 fs, where ω₁₀ is the frequency of the 0 → 1 transition. The first five Morse levels and eigenfunctions are shown in Fig. 1.

We numerically compute the vibrational energies ω_ν and eigenstates $ν$ of the potential using a discrete variable representation (DVR) method with a uniform grid and Fourier basis functions.⁴⁸ We use up to N_q = 721 grid points over the interval 2.5 < q < 20.5 (a.u.) along the nuclear coordinate. We can construct a quasi-complete nuclear basis $ν$ with up to ν_max ∼ 700, including states above the dissociation threshold. However, in most static and dynamical calculations considered, converged results can be obtained with ν_max ∼ 20–80, depending on the coupling strength and molecular species.

For the definition of the electric dipole function along the nuclear coordinate, we follow Ref. 19 and consider the universal model,

where the set of parameters (d₀, c₀, q₀, σ) can be chosen such that d(q) can equally well describe the qualitative behavior of IR-active molecular species that are polar or non-polar at equilibrium. We use this to compare three main types of molecular species according to the form of their dipole function d(q): non-polar molecules for which d(q_e) = 0 (e.g., iron pentacarbonyl), polar molecules with |d(q_e)| > 0 and $(d/dq)[d(q)]|qe>0$ (e.g., hydrogen fluoride), and finally polar molecules with |d(q_e)| > 0 and $(d/dq)[d(q)]|qe<0$ (e.g., sodium nitroprusside). In what follows, we, respectively, denote these cases as non-polar (NP), polar-right (PR), and polar-left (PL). The set of model function parameters used in this work is given in Table I (in atomic units). In the last row, we also include the value of the transition dipole moment $d10=1d(q)0$⁠.

In Fig. 2, we plot the dipole functions d(q), permanent dipole moments d_ν and transition dipole moments d_ν0 for the first 20 bound states of each of the three types of molecular species parameterized in Table I. The squares of the transition dipoles d_ν0 are proportional to the oscillator strength of the vibrational transition 0 → ν in linear infrared absorption. The three model functions qualitatively reproduce the typical behavior in infrared absorption spectroscopy, with a strong fundamental peak (0 → 1) and weaker overtones (0 → 2, 3, …). The permanent dipole moments also behave as expected, with their magnitude decreasing rapidly with ν near the dissociation threshold.

In all our calculations, the cavity frequency ω_c is set to be on exact resonance with the fundamental vibrational frequency ω₁₀. The vacuum field amplitude $E0$ is considered as a tunable parameter, simulating the fact that in real cavities, the magnitude of the light–matter coupling strength can be tuned by changing the intracavity molecular density, for fixed cavity geometry and material composition.⁴ For consistency between our coordinate-space and Hilbert space calculations, throughout this work, we parameterize the light–matter coupling strength by the dimensionless quantity

This parameter corresponds to the light–matter coupling ratio g/ω_c used in the ultrastrong coupling literature,^31,32 if we define $d10E0≡g$⁠.

The dimensionality of the cavity Hilbert space is set to ensure convergence of static and dynamical calculations. For the values of λ_g considered in this work, converged results were obtained by including Fock states $nc$ with up to n_c ∼ 80 photons. Convergence of the calculations in the coordinate-space representation is discussed below.

We are interested in the dynamics of expectation values of the form ⟨ψ(t)|Ô|ψ(t)⟩, where Ô is any molecular or photonic observable of experimental interest such as the photon number operator (â^†â) and the mode distance operator (⁠$q^$⁠). The system state $ψ(t)$ is obtained by propagating numerically an initial light–matter wavepacket $ψ0$ with a unitary propagator Û(t), i.e., $ψ(t)=Û(t)ψ0$⁠. The propagator and the wavefunction can be accurately represented in the coordinate-space representation using the MCTDH method and also in Hilbert space by projecting the total Hamiltonian into a vibrational and photonic product basis. Since MCTDH results are numerically exact, they serve as a benchmark. However, it is not clear how to distinguish the roles of vibrationally averaged permanent dipole moments from transition dipoles in the light–matter dynamics using the coordinate-space picture. This is straightforward to do in the Hilbert space representation.

State evolution in Hilbert space is carried out by projecting the time-evolution operator $Û(t)=exp[−iH^′t]$ into a truncated energy basis of the system Hamiltonian $H′^$⁠. For calculations that only take into account the contribution of transition dipole moments to the light–matter interaction, we set $H^2=0$ [Eq. (4)]. Energy eigenstates $ϵj$ satisfy $H^′ϵj=ϵjϵj$⁠, where j labels discrete and quasi-discrete polariton energy levels. Physically, energy eigenstates correspond to anharmonic vibrational polaritons.¹⁹ 

We project an arbitrary initial light–matter state into the polariton basis to read $ψ0=∑j⟨ϵj|ψ0⟩ϵj$⁠. In this basis, polariton wavepackets undergo a trivial phase evolution of the form

where the summation over the discrete label j includes the absolute ground state (j = 0) and runs up to a spectral cutoff (j = j_max). The latter is chosen such that observables ⟨Ô(t)⟩ are simultaneously converged over the entire evolution time. In this approach, convergence of the polariton spectrum is essential for computing accurate wavepacket dynamics.

In the coordinate-space representation, the system Hamiltonian in Eq. (1) is assumed to describe a two-dimensional potential energy surface (2D-PES) with respect to nuclear and photonic coordinates, together with the corresponding kinetic energy terms.^36,38,49,50 The light–matter interaction term proportional to x × d(q) is regarded as non-adiabatic coupling between the anharmonic nuclear potential V(q) and the harmonic cavity potential $V(x)≡ωc2x2/2$⁠.

In this 2D-PES representation, the evolution of an arbitrary light–matter wavepacket ψ₀(x, q) is accurately computed using the multi-configurational time-dependent Hartree method (MCTDH),^40,51 as implemented in Ref. 41. By projecting a wavepacket ψ(x, q, t) into the polariton eigenbasis, static properties such as the polariton spectrum (ϵ_j) can also be obtained.

Static and dynamic calculations are carried out using DVR with a sine primitive basis for the molecular coordinate q, in a grid of N_q = 721 points in the range 2.5 ≤ q ≤ 20.5 a.u. For the photonic coordinate x, we use harmonic oscillator primitive basis functions on a dimensionless DVR grid within the range −90 ≤ x ≤ 90. The number of photonic grid points N_x is chosen such that all the Fock states needed for convergence can be properly described.

Technical details of the MCTDH method are given in Ref. 40. Briefly, the time-dependent Schrödinger equation is solved by introducing the ansatz,

which is an expansion of the system state into a time-varying diabatic product basis composed of the function $ϕjq(q,t)$ in the nuclear coordinate and $ϕjx(x,t)$ in the photonic coordinate, each labeled by the integers j_q and j_x, respectively. We set the number of basis functions in the expansion to n_q = n_x = 30. The equations of motion for the coefficients $Ajq,jx(t)$ and the basis functions on the corresponding 2D product grid are obtained using an open-source MCTDH package⁴¹ and then solved using a numerical integrator. The evolution of system observables ⟨ψ(x, q, t)|Ô(x, q)|ψ(x, q, t)⟩ is obtained by numerically evaluating the corresponding integrals on the (x, q)-grid.

For static calculations (e.g., polariton spectrum), we use a routine in the MCTDH package that diagonalizes the Hamiltonian matrix on the 2D grid using the Lanczos algorithm.⁵² The dimensionality of the Hamiltonian matrix is N_qN_x, where N_q and N_x are the number of grid points in the nuclear and photonic coordinates, respectively. In the Lanczos method, the eigenvalue problem is transformed into an iterative eigendecomposition. After the number of iterations is set, a defined number of eigenvalues M are computed, with M < N_qN_x.

Our first goal is to compare the static properties of intracavity vibrational polaritons that emerge under strong light–matter coupling for molecular species that are either polar or non-polar at equilibrium (see the definition in Table I). This comparison is best carried out in the Hilbert space representation because we can controllably neglect the contribution of permanent dipole moments to the light–matter interaction, by setting $H2^=0$ in Eq. (4).

In Fig. 3(a), we plot the polariton spectrum as a function of the dimensionless light–matter coupling strength λ_g, for a non-polar molecule treated in two alternative ways: (i) both transition and permanent dipole moments in the vibrational eigenbasis are taken into account (solid curves); (ii) only transition dipoles are considered (dashed curves). For a given value of λ_g, energies are shown relative to the energy of the absolute ground state (E_GS). The ground level in general exhibits a red shift relative to its energy at λ_g = 0 (not shown). For λ_g ≲ 0.1, the lowest excited polariton manifold exhibits the usual polariton doublet. The second excited manifold has a well-defined triplet structure. At larger coupling strengths λ_g > 0.1, the spectrum develops into several real and avoided crossings between excited polariton levels.¹⁹ If we consider only the contribution of transition dipole moments to the Hamiltonian, the low-energy excitation spectrum remains mostly unaltered in comparison with the spectrum of the full Hamiltonian for a wide range of coupling strengths. However, for higher coupling strengths λ_g ≳ 0.3, ignoring the permanent dipole moments results in a qualitatively different polariton spectrum.

In Figs. 3(b) and 3(c), we show the mean bond distance $⟨q^⟩$ and mean intracavity photon number ⟨â^†â⟩, respectively, for the polariton ground state. The model predicts that at moderate values of λ_g ∼ 0.1, the ground state behaves qualitatively different from the cavity-free diabatic ground state $ν=0nc=0$⁠.

Figure 3(b) shows that the mean bond distance can be significantly higher in the polariton ground state (solid line) than the bond distance of the bare ν = 0 vibrational level in free space. This type of mode lengthening effect is not captured when only transition dipole moments are taken into account in light–matter interaction. By ignoring the permanent dipole moments of the vibrational states, the ground state is seen to experience bond shortening instead (dashed line). This is consistent with the results obtained in Ref. 19, which neglect state-dependent permanent dipoles.

Figure 3(c) shows that the polariton ground state is composed of Fock states with n_c ≥ 1 photons, even for coupling strengths as low as λ_g ∼ 0.01. Photon buildup in the ground state is a conventional signature of the ultrastrong coupling regime,^31,32 for two-level quantum emitters (qubits). This unexpected behavior of molecular vibrations is largely insensitive to the inclusion or neglect of permanent dipole moments in the system Hamiltonian, for non-polar species. However, we show below how this relative insensitivity does not hold for vibrational modes that are polar at equilibrium.

In Fig. 4, we show the polariton spectrum and ground state properties as functions of the coupling strength λ_g, for nuclear modes that are polar at equilibrium. For concreteness, we consider the polar-right electric dipole function in Table I. In general, the results are more sensitive to the presence or absence of permanent dipole moments in the Hamiltonian than for the case of non-polar molecules. For instance, Fig. 4(a) shows that the exact polariton energies in the first and second excited manifolds (solid lines) differ by a few percent from the energies obtained by neglecting permanent dipole moments (dashed lines). The difference is more pronounced for λ_g > 0.1.

Figure 4(b) shows that the mode lengthening effect predicted earlier for non-polar molecules (Fig. 3) is stronger in the polariton ground state of polar molecular species (solid line). Again, ignoring permanent dipole moments gives qualitatively different results (dashed line).

Figure 4(c) shows that for polar vibrations, the cavity field can build up a significant amount of photons in the polariton ground state, reaching up to ⟨â^†â⟩ ∼ 10 for λ_g ∼ 0.4 (solid line). We find this prediction to be sensitive to the presence or absence of permanent dipole moments in the Hamiltonian. Considering only the contribution of transition dipole moments to the light–matter interaction gives photon numbers that are consistently lower by about two orders of magnitude relative to the full dipole matrix results, over the entire range of coupling strengths (dashed line). This should be compared with Fig. 3(c), where only small differences are found.

Let us now consider the short-time unitary dynamics of the molecular and photon observables $⟨q^⟩$ and ⟨â^†â⟩, for a molecule-cavity system initially prepared in a general state of the form

where $βνnc$ are complex wavepacket coefficients in the diabatic basis $νnc$⁠. Transforming Eq. (15) into the coordinate-space representation preserves the values of $βνnc$⁠. This initial state can be propagated either in Hilbert space or in coordinate space as described in Sec. III C.

We first consider an initial condition that describes a system in which an individual molecular vibration in its ground level (ν = 0) is embedded at t = 0 into a cavity that was previously prepared in a coherent state (e.g., via laser pumping). The molecular vibration is assumed to be non-polar at equilibrium (see Table I), and the cavity initially has |α|² = 2 photons on average, where α is a real coherent state amplitude. We set λ_g = 0.125.

Physical intuition suggests that the ground state molecule should interact with the finite cavity field and become vibrationally excited by absorbing cavity photons. This intuitive picture is, indeed, reproduced by converged MCTDH calculations in Fig. 5. Over the first 150 fs, the molecule is seen to increase its mode length from the bare value at ν = 0 to the value for ν = 1, by absorbing a single cavity photon on average.

The behavior over larger time scales (∼10² fs) replicates at much shorter time scales on the order of the vibrational period (2π/ω₁₀ = 8.27 fs). In the first 40 fs, nuclear oscillations resulting in an overall mode lengthening [Fig. 5(a)] are accompanied by a stepwise decrease in the cavity photon number [Fig. 5(b)]. The partial recurrence of the photon number toward its initial value in turn results in the overall shortening of the vibrational mode. This qualitative behavior tends to repeat in cycles, but the amplitudes of the oscillations in the nuclear and photonic observables become reduced at later times.

Figure 5 also shows that the evolution of polariton observables using the Hilbert space approach gives results that are indistinguishable from those obtained using the MCTDH method, as long as the entire dipole matrix is used, i.e., both transition dipoles and permanent dipoles contribute to light–matter coupling. Small deviations from the exact MCTDH evolution are found when we neglect the permanent dipole moments [Fig. 5(b)].

We now consider the polariton dynamics of polar molecules, specifically polar-right species (see Table I). We again consider the vibration initially in ν = 0 and the cavity field in a coherent state with a real amplitude $α=2$⁠. Figure 6 shows the resulting evolution of the mean mode length and the intracavity photon number. The coupling strength λ_g is the same as that in Fig. 5. The results suggest that the intuitive physical picture of a vibrational mode overall being excited by absorbing a cavity photon is not universal. For polar species, MCTDH calculations predict an overall increase in both the mean mode length and cavity photon number over a few hundred femtoseconds, as the result of light–matter coupling (solid lines). Quantitatively equivalent results are obtained by computing the dynamics in Hilbert space with both diagonal and off-diagonal dipole matrix elements taken into account (crosses). If we ignore the contribution of permanent dipoles to the Hamiltonian, we obtain qualitatively different results (dashed lines). Specifically, both the mean photon number and mode length are consistently underestimated if permanent dipoles are neglected.

Figure 6(a) shows that the mean mode length behaves qualitatively similar to the case of non-polar molecules, i.e., the molecule experiences an overall mode lengthening in comparison with its initial ν = 0 configuration. However, for a polar molecule, the mode is lengthened by a greater extent at equal coupling strength. In Fig. 6(a), the mode length reaches a quasi-steady value that is comparable with the length of a bare Morse vibrational state with ν = 3 vibrational quanta.

Figure 6(b) shows that over short time scales on the order of the bare vibrational period (8.27 fs), the system exhibits alternating patterns of mode lengthening at the expense of photon absorption and mode shortening accompanied by re-emission of photons into the cavity field. These cycles were also found to occur for non-polar species (Fig. 5). However, for polar vibrations, we predict a qualitatively new physical behavior: The light–matter system rapidly develops a sizable number of cavity photons over the first few vibrational periods. In less than 5 fs, the cavity field amplitude rapidly evolves to a state with about ⟨â^†â⟩ ≈ 6.3 photons. Subsequent evolution of the photon number occurs in cycles of decreasing amplitude over the first hundred femtoseconds.

Coherent states are known to have classical field statistics.⁵³ Let us now consider a qualitatively different scenario in which the cavity field is initially prepared in a state with no classical analog. We are interested in the differences that can be expected in the dynamics of vibrational polaritons that evolve from a field state with quantum statistics (e.g., Fock state), relative to the dynamics of a cavity initially in a coherent state (Figs. 5 and 6). If differences are found, we want to understand the role of bond polarity in distinguishing these two scenarios.

We first consider a light–matter system in which a polar molecular vibration is prepared in its ν = 0 level and then suddenly embedded into a cavity with no photons (vacuum). The initial state in the diabatic basis is simply

In Fig. 7(a), we show the evolution of the mean intracavity photon number starting from this initial state, for a polar-right molecular vibration (see Table I). The coupling parameter is λ_g = 0.25. Numerically exact MCTDH results show that starting from the diabatic vacuum at t = 0, the cavity field develops a significant amplitude over the first few femtoseconds, with the photon number reaching up to ⟨â^†â⟩ ≈ 7.3 in 5 fs (solid line). The polariton wavepacket then evolves into a state with a photon number that undergoes moderate fluctuations about a relatively large number (e.g., ⟨â^†â⟩ ≈ 4.6 in 300 fs). Calculations performed with the Hilbert space approach taking into account both transition and permanent dipole moments (crosses) give results that are equivalent to those obtained with MCTDH. However, when we ignore the contribution of permanent dipole moments to the light–matter interaction (dashed line), the results consistently underestimate the intracavity photon number by at least two orders of magnitude.

For comparison purposes, in Fig. 7(b), we show the evolution of photon number for a molecule in an initial vibrational wavepacket placed in the diabatic cavity vacuum, i.e., $ψ0=(∑νcνν)⊗nc=0$⁠. The vibrational wavepacket has a mean vibrational energy of Ē_vib = 2.126ω₁₀ and is initialized at equilibrium (⟨q⟩₀ = q_e) with a positive mean velocity. The overall photon number evolution is found to be similar to Fig. 7(a), which suggests that cavity photon generation is more related to the structure of the light–matter coupling Hamiltonian than the initial degree of molecular excitation.

We perform simulations for non-polar molecules, also initialized in the diabatic ground state $ν=0nc=0$ (not shown). For the same light–matter coupling parameters as those in Fig. 7, we find that for this class of molecule, the mean mode length does not significantly vary in comparison with free space, and the average number of cavity photons that can be produced remains at least two orders of magnitude smaller than the case of polar molecules, even after several hundred femtoseconds.

We have established that the sub-picosecond dynamics of vibrational polaritons formed with polar vibrations is qualitatively different from the evolution of polaritons formed by vibrations that are non-polar at equilibrium. We now discuss what appears to be an unexpected feature of vibrational polaritons for polar molecules: not only the value of the dipole moment at equilibrium is relevant to determine the dynamics of vibrational polaritons but also equally important is the sign of the slope of the electric dipole function.

We support this conjecture by comparing the dynamics of polariton observables for polar-right and polar-left vibrational modes (see definitions in Table I). In our Morse oscillator model, the transition dipole moments (d_ν′ν) of polar-left and polar-right vibrations have the same qualitative structure as a function of the vibrational quantum number ν (see Fig. 2), but their permanent dipole moments (d_ν) behave qualitatively different. d_ν decreases monotonically with ν for polar-left vibrations. For polar-right species, d_ν increases with ν for low quantum numbers, then decays to small values for higher vibrational levels. In other words, the slope Δd_ν/Δν has opposite signs at low ν, for our two types of polar vibrations.

In Fig. 8, we correlate the sign of Δd_ν/Δν for polar-left and polar-right species, with the evolution of the mean mode length and intracavity photon number under strong coupling. The polariton dynamics is computed with the MCTDH method. We consider the same initial condition as that in Fig. 6, i.e., molecules are initially in their ground vibrational level and the cavity is in a coherent state with two photons on average.

Figure 8(a) shows that starting from the same initial condition, polar-right vibrations experience an overall increase in the mode length over the first few hundred femtoseconds and polar-left species undergo an effective compression of their mode length over the same time interval. For the coupling strength λ_g = 0.25, the mode length of a polar-right molecule increases to a value comparable with a ν = 8 vibrational eigenstate. On the contrary, a polar-left vibrational mode can only transiently exceed the initial mode length of the ν = 0 vibrational level.

Figure 8(b) in turn shows that the lengthening of the vibrational mode for a polar-right species is accompanied by a significant increase in the intracavity photon number. For the parameters considered, the cavity evolves from the vacuum state into a polariton wavepacket with ⟨â^†â⟩ ≈ 7.0 in about 300 fs. In contrast, for polar-left species, the cavity on average gives away about one photon into the material system.

Until now, we have discussed the dynamics of idealized Morse vibrations under strong light–matter coupling. Although Morse oscillators are good approximations to the stretching modes of many molecular species, we have yet to show that the light–matter physics predicted above is also relevant for real molecules.

In order to prove this, let us first consider the polariton dynamics of an individual hydrogen fluoride (HF) molecule in an infrared cavity resonant with its fundamental frequency (⁠$ν̃10=3990cm−1$⁠, 2π/ω₁₀ = 8.36 fs). In Fig. 9(a), we plot the ab initio potential energy curve V(q) and electric dipole function d(q) for the ground electronic state (X¹Σ⁺). These were obtained using the electronic structure package MOLPRO. In both cases, a complete active space (CAS) calculation of the lowest electronic state has been performed using a multi-configurational self-consistent field (MCSCF) method, then followed by the multireference configuration interaction (MRCI) method, using Dunning’s correlation consistent basis set with diffuse functions aug-cc-pVQZ. According to our nomenclature, HF belongs to the polar-right class of molecules.

In Fig. 9, we show the evolution of the mean mode length $⟨q^⟩$ and the intracavity photon number ⟨â^†â⟩, for a vibrational polariton wavepacket that evolves from an HF molecule initially prepared in its ground vibrational level (ν = 0) with a cavity in a vacuum Fock state (n_c = 0). The results show that the molecule experiences an overall increase in its mode length of about 10% in 300 fs, accompanied by a significant buildup of intracavity photons. This behavior is consistent with the Morse model calculations in Fig. 8.

As another example, let us consider the dynamics of an individual carbon dioxide molecule (CO₂) in an infrared cavity resonant with the asymmetric stretching mode (⁠$ν̃10=3860.7cm−1$⁠, 2π/ω₁₀ = 8.64 fs). This vibrational mode is not polar at equilibrium but acquires a dipole moment away from it. Therefore, the molecule belongs to the non-polar class. In Fig. 10(a), we plot the ab initio ground state potential energy curve (¹A′ state) and electric dipole function near equilibrium along the asymmetric stretching mode. These were obtained using the same ab initio method used for HF (MCSCF/aug-cc-pVQZ).

In Fig. 10, we show the evolution of the polariton mode length and photon number, over the first few hundred femtoseconds. The polariton wavepacket is also assumed to evolve from a vibrationless molecule in a vacuum Fock state. Unlike the cases considered in previous sections, the mode length in CO₂ is found to remain invariant at its bare equilibrium value throughout the evolution of the polariton state, which we attribute to the negligible anharmonicity of the potential energy curve. In contrast, the intracavity photon number does vary significantly from its initial value zero, rapidly reaching up to ⟨â^†â⟩ ∼ 10⁻³ after a few vibrational periods, from where it undergoes quasi-stationary oscillations.

We have shown that the resonant interaction of an individual molecular vibration with a quantized cavity field can have very different physical observable consequences depending on the dipolar properties of the molecular electron density. This is not completely unexpected if we recall that the oscillator strength S_ν′ν of the infrared absorption band for a vibrational transition ν → ν′ is proportional to the square of the slope of the electric dipole function d_e(q) at the equilibrium configuration.³³ However, our results show that under conditions of strong and ultrastrong light–matter coupling, not only the slope magnitude but also the entire shape of the electric dipole function are important to understand the static and dynamical properties of vibrational polaritons. In other words, two molecules that have nominally identical infrared absorption bands in free space may experience qualitatively different polariton dynamics inside an infrared cavity.

Consider polar-right and non-polar molecules, following the nomenclature in Fig. 2. These two species essentially have the same dipole transition strength for the fundamental IR absorption peak |d₁₀|². Therefore, if they also have the same fundamental frequency ω₁₀, they would be indistinguishable in IR spectroscopy. However, their evolution inside a strongly coupled resonant cavity would be qualitatively different. For instance, our results show that when a polar-right vibration in its ground vibrational level outside the cavity is suddenly placed inside a cavity in total darkness (no photons), the system evolves into a hybrid light–matter wavepacket that behaves as if the molecular mode lengthens over sub-picosecond time scales, while the cavity spontaneously generates photons at the same time interval (see Fig. 7).

We also unexpectedly find that polar molecules can either experience bond lengthening or bond shortening, depending on the shape of their electric dipole function (see Fig. 8). For interaction strengths at the conventional onset of ultrastrong coupling,^31,32 polar-right molecules can increase its length by up to 10% from its equilibrium value (see Fig. 9 for hydrogen fluoride), while developing an intracavity field with up to about 10 photons on average in a few hundred femtoseconds.

Mode lengthening in a strongly coupled infrared cavity may in turn result in different chemical reactivities of polar molecules, in comparison with free space. Our theoretical and numerical analyses may, thus, provide a consistent basis for the development of a reaction rate theory for vibrational polaritons that can offer a microscopic understanding of the observed chemical reaction enhancements under vibrational strong coupling,⁸ which remains to be an outstanding goal in the field.

For polar and non-polar molecular vibrations, we predict an ultrafast dynamical buildup of intracavity photons during the evolution of a strongly coupled light–matter system with initially no vibrational or photonic excitations. This spontaneous generation of photons is a natural outcome of the wavepacket evolution in the polariton energy basis, since diabatic initial product states $νnc$ are not eigenstates of the light–matter Hamiltonian for any finite coupling strength (λ_g ≠ 0). In the adiabatic polariton basis, vibration–photon product states can, thus, be seen as wavepackets that have broad distribution polariton energies. As long as the wavepacket has contributions from polariton excited states, then the system can radiate intracavity photons into the far field through radiative decay.

The energy needed to produce cavity field excitations together with the vibrational excitation of a molecule comes from the light–matter interaction Hamiltonian itself. This is simpler to visualize in the Hilbert space representation [see Eqs. (3) and (4)]. The counter-rotating terms of the light–matter can directly couple diabatic states $νnc$ that can be interpreted as the simultaneous excitation of both the vibrational mode and the cavity field. If these counter-rotating couplings are significant, we can, thus, expect an overall increase in the photon number and mode length of the system.

Experimentally, it may be challenging to initialize a strongly coupled cavity–vibration system in a diabatic product state with a definite number of vibrational excitations and a definite photon number $νnc$⁠. We propose a suitable preparation scheme in Fig. 11. At early times (t ≪ t₀), an individual molecule is placed within a relevant plasmonic mode volume in total darkness (i.e., no laser driving). Although the vacuum field amplitude at the molecular–nanoparticle interface can be large, the molecular vibration does not exchange energy with the near-field vacuum because either the fundamental frequency $ωvib′$ is far detuned from the relevant plasmon frequency ω_cav (Δω ≡ ω_c − $ωvib′$ ≠ 0) or the dipole moment for the relevant vibrational transition is such that the light–matter coupling is weak (λ_g ≈ 0), or both conditions occur simultaneously (as in Fig. 11). Under these conditions, the molecule–nanoparticle system will simply thermalize with its environment and remain unaltered in the absence of additional external perturbations. In particular, the number of photons in the near field will not exceed the level imposed by the background radiation, which is negligibly small at infrared frequencies (⟨â^†â⟩ ≈ 0).

Strong light–matter coupling is suddenly activated over a time interval T_chem = t₁ − t₀ by chemically converting the adsorbed molecule into a species with a vibrational mode that either becomes resonant with the near-field vacuum (Δω ≈ 0) or the relevant vibrational transition dipole is such that light–matter coupling becomes strong (λ₀ ≳ 0.1), or both situations occur simultaneously (as in Fig. 11). UV–Vis photochemistry⁵⁴ or electron tunneling⁵⁵ can be used to activate the bond-forming reaction on the surface of the nanoparticle at t₀. After the reaction is complete (t ≥ t₁), the molecule–nanoparticle system is left in the strong coupling regime of light–matter interaction. Assuming that the strongly coupled bond is formed in the vibrational ground state, the system at t₁ is left in the diabatic product state $ν=0nc=0$⁠, from where it evolves into a polariton wavepacket that can eventually radiate a number of infrared photons into the far field over sub-picosecond time scales (⟨â^†â⟩ ≳ 1). The results in Fig. 7 show that even if the molecule is left in an excited vibrational state at t₁, a similar polariton evolution can be expected. In the absence of external infrared driving fields, the generation of near-field quantum light stops when the polariton wavepacket relaxes to the polariton ground state.

In the proposed scheme, the reaction time scale T_chem is very important. If this is comparable or much shorter than the Rabi oscillation period for the vibrational populations under strong light–matter coupling, the system evolves as described in the main text. However, if T_chem is much larger than the relevant Rabi period, then the light–matter system evolves adiabatically from an uncoupled product eigenstate $νnc$ into a polariton eigenstate, not a polariton wavepacket. For instance, the uncoupled bare vibrationless vacuum $ν=0nc=0$ would evolve adiabatically into the polariton ground state at t₁, from where no photons can be emitted. In a realistic experimental setup, the evolution would be neither sudden nor adiabatic. In this most general case, a diabatic product state $νnc$ prepared at t₀ evolves into an excited wavepacket in the polariton energy basis.

As long as the polariton wavepacket prepared at t₁ has a non-vanishing component in an excited polariton eigenstate, the state will relax by emitting photons into the far field. The frequencies of these emitted photons will match the Bohr frequencies between polariton levels, as energy conservation dictates. However, modeling the rates of vibrational polariton relaxation and dephasing under conditions of ultrastrong light–matter coupling is expected to be a challenge, as the standard quantum optical master equation is known to fail for two-level atoms in this regime.^56,57

In summary, we have developed the first detailed theoretical framework for understanding the sub-picosecond dynamics of anharmonic vibrational polaritons. Starting from a fundamental light–matter interaction model in the electric dipole approximation, we correlate the dynamics of both material and photonic observables of a strongly coupled cavity–vibration system with the underlying electrostatic properties of vibrations determined by the molecular electronic structure. Using numerically exact quantum dynamic methods, we show that the entire shape of the electric dipole function of anharmonic vibrational modes is relevant for understanding the femtosecond dynamics of vibrational polaritons in infrared cavities. Our single-molecule analysis may stimulate further developments on vibrational strong coupling with nanophotonics.^58–60 The results can also be extended in order to take into account many-body effects, as well as dissipative photonic and material processes. Such extensions would enable a more direct comparison with available experimental evidence in liquid-phase cavities.^8,13,18,61

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

J.F.T. is supported by ANID through Postdoctoral Fellowship Grant No. 3200565. F.H. is supported by ANID through grants REDES ETAPA INICIAL, Convocatoria 2017 No. REDI 170423, and FONDECYT Regular No. 1181743. All authors acknowledge the generous support of Iniciativa Científica Milenio (ICM) through the Millennium Institute for Research in Optics (MIRO).