Investigating the collective nature of cavity-modified chemical kinetics under vibrational strong coupling

2.1 Model system

In this work, we consider a minimal model of a set of N molecules collectively coupling to a single radiation mode supported by a cavity. Figure 1(a) (see left panel) shows an example of a setup where N molecules in a simulation cell (red shaded area) are placed inside a cavity with width L. There are two ways to analyze the quantum dynamics in the thermodynamic limit: A → ∞ and N → ∞, and these two approaches are illustrated in Figure 1(a) and (b).

Figure 1:

A graphical representation of two possible ways for increasing the number of molecules included in the cavity QED system. (a) We consider a system at a fixed density ρ and simply increase the volume of the simulation cell (thereby increasing the number of particles included but also increasing the volume). (b) Treating a fixed simulation volume but varying the number of particles present in it (i.e., changing the density of particles in the cavity by altering N). This is the approach taken in references [27], [29].

The Hamiltonian describing the set of molecular vibrations coupled to cavity radiation is written as

where H ̂ M is the matter Hamiltonian, H ̂ R M is the radiation-matter coupling Hamiltonian, and H ̂ C is the Hamiltonian for the cavity radiation mode. In the single-mode approximation of a perfect lossless optical cavity, the cavity Hamiltonian becomes

The presence of cavity loss can be modeled with a simple Caldeira–Leggett Hamiltonian [30], where H ̂ C is written as

with q ̂ c = b ̂ c † + b ̂ c / ( 2 ω c ) , and with C k and ω ̃ k sampled from a spectral density function J loss ( ω ) = π 2 ∑ k C k 2 ω ̃ k 2 δ ( ω ̃ k − ω ) that is taken to be of the Debye form, approximating the result for a simple 1D model cavity [31].

In the long-wavelength and single cavity mode limit, the light–matter interaction Hamiltonian can be written as

where e and ω _c are the polarization direction and frequency of the radiation mode, respectively,

with V the quantization volume, and we take the dipole moment operator of molecule i as

with n _i a normalized vector that specifies the orientation of molecule i. Here, η _c quantifies the light–matter interaction strength and may depend on N depending on how we choose to increase the number of molecules present in the system. Note that the light–matter interactions described in Eq. (5) ignores the spatial dependence of the radiation field (long-wavelength approximation), which may break down when considering a large number of molecules filling the entirety of the optical cavity.

As mentioned above, Figure 1 illustrates two possible strategies to analyze the effect of increasing the number of molecules present in the simulations. The second strategy, illustrated in Figure 1(b), is the choice that has been employed in many recent theoretical works (e.g., [27], [29]). Within the context of the first approach, we consider a simulation box with cross-sectional area (perpendicular to the cavity direction) A that contains N molecules and use the volume of this box as the quantization volume present in χ. Doing so, we find that

where we have used that 2L = λ and that λω _c = 2πc and have absorbed all of the constants into λ. Introducing an average particle density per unit of cross-sectional area (which is a constant in the experimental setup) ρ = N/A, we find that η _c

and so the final light–matter coupling Hamiltonian in this case can be written

where g = λ ρ . This choice gives rise to a constant Rabi splitting as the number of molecules is increased.

In the second approach, we vary the number of particles and keep the quantization volume constant while increasing the number of particles present in the volume. As a result, the average particle density per unit of cross-sectional area is constant. In this case, we write

where η _c is independent of N, as has been done in recent work [18], [27], [29]. This choice of Hamiltonian gives rise to a Rabi splitting that increases as N .

In the absence of direct dipole–dipole interactions, the total Hamiltonians in each of these two cases are provided below. The Hamiltonian for the constant ρ scenario is written as

with H ̂ M describing the bare molecular Hamiltonian. On the other hand, the Hamiltonian for the constant V approach is written as

We note that the intermolecular dipole self-energy term is known to cancel with the direct Coulomb interaction term when considering all radiation modes in the light matter Hamiltonian beyond the long-wavelength approximation [32]. However, recent work [33,34] suggest that when performing a radiation mode truncation, as is done here (considering a single radiation mode), this term should be explicitly included within the dipole gauge Hamiltonian [35,36]. We also find that this term does not contribute to resonantly modifying the chemical reaction rate.

In this work, we investigate two model molecular systems. In the first, Model I, we consider a molecular system described by a double-well potential energy surface with the reaction coordinate coupled to a set of dissipative solvent degrees of freedom. In the second, Model II, we consider a proton transfer reaction where the reaction coordinate is coupled only to a spectator mode. Below, we describe each of these models.

2.2 Methods

2.2.2 Multi-Layer Multiconfiguration Time-Dependent Hartree approach

For simulating the exact quantum dynamics of Model II, which is a closed quantum system, we use the Multi-Layer Multiconfiguration Time-Dependent Hartree (ML-MCTDH) approach. ML-MCTDH uses a tensor network based ansatz for the total wavefunction and here we employ an efficient projector splitting integrator [56–62] that avoids many of the issues that arise in the presence direct product state initial wavefunction conditions [62]. In all calculations, we used a single-site with subspace expansion integration scheme similar to those proposed in reference [63] for the standard ML-MCTDH algorithm to allow for growth of the bond-dimension (or number of single-particle functions) throughout the dynamics. In all calculations, we consider a ML-MCTDH tree consisting of a balanced binary tree for representing the N molecule system that is connected to the cavity mode at its root. An example of the topology of the ML-MCTDH wavefunction for Model II with N = 16 is show in Figure 2. Here, we find converged results allowing a maximum bond dimension (or number of single particle functions) in the tensor network to be 32. The method described below can be generalized to finite temperature, but for simplicity, we work at zero temperature.

Figure 2:

Multilayer tree used for the N = 16 molecule ML-MCTDH calculations for Model II. Bond dimension (number of single particle functions) is shown on the edges connecting circular nodes, and the primitive Hilbert space dimensions are shown for bonds connecting circular and square nodes. Here, C denotes the cavity degree of freedom and as M _i denotes the ith molecules degree of freedom.

For Model I, we have made use of an ML-MCTDH representation of the HEOM for all simulations with N > 2. Here, we used the same balanced binary tree structure for the ML-MCTDH tree, with each leaf mode of the molecule representing the full space of auxillary density operators (ADOs) for that molecule and its environment. In this case, we find converged results allowing a maximum bond dimension in the tensor network of 48.

2.2.2.1 Initial state preparation

For this model (Model II) molecular system, we explore two types of initial conditions. The first is an uncorrelated initial condition in which the cavity is prepared in the vacuum state in the dipole gauge, as has been employed in recent works on collective VSC [19],

(30) ψ u c ( 0 ) = 0 ⊗ | ψ g s , i ⟩ ,

where ψ g s i is the ground state of the molecule Hamiltonian, H ̂ M , for molecule i, and 0 corresponds to the cavity vacuum state. Such an initial condition is relevant to the case in which the N molecules are thermalized in the absence of the light–matter coupling, and at t = 0, the light–matter coupling is introduced.

Second, we consider a correlated ground state of the total N molecules and the cavity radiation mode, e.g., the polaritonic ground state. This choice for initial condition relates to the case in which interactions between the N molecules and the cavity is introduced at t → −∞ and the total system is allowed to evolve under the composite Hamiltonian to reach the ground state at t = 0 (e.g., the thermalized case).

We obtain the ground state wavefunction using a single-site tree tensor network state optimization algorithm presented in Ref. [64], however, with the use of subspace expansion allowing for adaptive control of bond dimension throughout the optimization. All simulations presented here were performed using the tree tensor network library ttns_lib [65].

We report the single-molecule side – total side correlation function

(31) ⟨ Θ i ( t ) Θ ( 0 ) ⟩ = ψ e i H ̂ t Θ i e − i H ̂ t ⊗ j = 1 N Θ j ψ ,

where Θ _i is the side operator projecting onto the reactants for molecule i, and ψ the initial wave functions described above. Here, we evaluate this quantity by independently evolving two initial wavefunctions, ψ 1 = ψ and ψ 2 = ⊗ j = 1 N Θ j ψ and evaluating the matrix elements

(32) ⟨ Θ i ( t ) Θ ( 0 ) ⟩ = ψ 1 ( t ) Θ i ψ 2 ( t ) ,

at each point in time.

3.1 Few molecule limit

We first discuss how a cavity modifies chemical reactivity in the few molecule limit. Via exact quantum dynamics simulations, we find that additional molecules coupled to the cavity radiation mode provide additional dissipation, leading to the enhancement of chemical reactivity for solvent interactions in the energy diffusion-limited regime.

In Figure 3, we present the frequency-dependent reaction rate, normalized by the out-of-cavity rate, associated with the individual molecules for N = 1, 2, 3, and 4 molecules in the cavity. As we demonstrated in our recent work [10], in the single molecule limit, cavity coupling leads to a resonant enhancement of the chemical reaction rate (blue solid line in Figure 3(c)) due to the additional dissipation originating from the bath that describes cavity loss. This is because the overall chemical reaction rate is limited by the rate of thermal relaxation when the molecule-solvent coupling is of a magnitude such that the reaction falls within the energy diffusion-limited regime. Additional sources of dissipation, such as the dissipation from coupling to a lossy cavity mode, naturally increase the rate of thermalization, leading to an enhancement of the reaction rate [10]. It is worth noting that the rate increases with the increase in the cavity loss rate (or a decrease in the cavity lifetime) achieved by increasing the coupling between a cavity mode and its bath in Eq. (3).

Figure 3:

Cavity-modified chemical dynamics in the few molecule limit of Model I at T = 300 K. (a) Double-well potential describing a model molecular system with the vibrational eigenstates. (b) Schematic illustration of multiple molecular systems coupled to a lossy cavity radiation mode. (c) Cavity photon frequency-dependent normalized chemical reaction rate constant when coupling N = 1, 2, 3, and 4 molecules to a radiation mode using the constant V approach (see Figure 1(b)). (d) Cavity-modified chemical rate constant as a function cavity lifetime τ _c at various N. (e)–(f) Same as (c) and (d) but with a constant ρ approach (see Figure 1(a)). For N = 1 and 2 molecules, the results were obtained using HEOM calculations. For N = 3 and 4 molecules, the HEOM/ML-MCTDH approach was used.

In Figure 3(c)–(f), we show that the cavity modification of chemical reaction rate is further increased as the number of molecules in the cavity is increased. We explore this in Figure 3(c) by keeping the per molecule coupling to the cavity mode a constant, which would correspond to keeping V fixed as illustrated in Figure 1(b). In Figure 3(c), the cavity enhances the chemical reaction rate by ∼ 25 % at N = 1 (Rabi splitting of ≈ 26.58  cm⁻¹). At N = 2, 3, and 4 (with the Rabi splittings of 37.541, 45.95, and 53.0 cm⁻¹), the cavity enhancement rises to ∼ 55 % . This interesting effect can be explained in terms of enhancement of the overall dissipation. From the perspective of one reactive molecule, additional molecules coupling to cavity can be viewed as additional dissipative bath degrees of freedom that are coupled to the cavity mode, thereby increasing cavity loss. However, the extent of this effect quickly saturates and as a result cavity modification of chemical reaction rate marginally changes when going from N = 3 to N = 4. This reveals that even though there is an N dependence of the cavity-modified chemical rate, this effect begin to saturate even for a small number of molecules coupled to the cavity. Whether this enhancement could be significant for a larger N remains an open question.

Figure 3(d) shows that the extent to which an additional molecule coupling to cavity can enhance chemical reactivity depends on the cavity loss rate. Here, we observe that increasing the number of molecules N reduces the sensitivity of the peak rate modification to the cavity loss rate. Intuitively, if the cavity mode coupling strength to its dissipative bath is already very large (i.e., very lossy), the impact of additional sources of dissipation would be comparatively negligible. Similarly, for a nearly perfect lossless cavity, the coupling to an additional molecule can lead to a significant increase in dissipation, and hence a substantial modification of the reaction rate for the reactive molecule. This can be seen in the cavity-modified reaction rate as a function of cavity lifetime, as shown in Figure 3(d). For a small cavity lifetime, τ _c < 400 fs, the cavity-modified reaction rate is nearly the same for N = 1, 2, 3, and 4. The presence of additional molecule starts to significantly alter chemical reactivity for τ _c > 400 fs.

Figure 3(e) and (f) presents the cavity-modified chemical reaction rate in the presence of a small number of molecules where the total collective coupling is kept constant with a constant Rabi splitting of 26.58 cm⁻¹. This is achieved by rescaling the per-molecule coupling constant to the cavity mode by a factor of 1 / N . Thus, the overall cavity modification of the rate occurs via the interplay of two competing effects: the scaling down of individual molecular coupling to the cavity versus the additional dissipation due to the coupling of other molecules to the cavity mode. For the particular solvent couplings chosen here, the cavity modification of the rate is maximized at N = 2 when using τ _c = 1000 fs (see Figure 3(e)). As before, whether an increase in N will lead to an increase in chemical reactivity also depends on the cavity lifetime.

Figure 3(f) presents the cavity-modified chemical reaction rate as a function of the cavity lifetime with N = 1, 2, 3, and 4. We observe that while for larger cavity lifetimes (τ _c > 800 fs) additional molecules further enhance the chemical reaction rate, chemical reactivity is less enhanced for smaller cavity lifetimes (τ _c < 500 fs). This is because the cavity mode is already strongly connected to a dissipative bath and adding another molecule does not change this dissipation much. As a result, for smaller cavity lifetimes, as the per-molecule coupling is reduced, the overall enhancement decreases when increasing the number of molecules. Interestingly, for larger cavity lifetimes (τ _c > 800 fs), the additional dissipation when increasing the number of molecules enhances the chemical reaction rate despite the lower per-molecule light–matter coupling.

3.2 Thermodynamic limit

We explore the N → ∞ using the mean-field HEOM approach explained in Section 2.2.1. As hinted at from the trends observed in Figure 3, we expect no cavity modification in the N → ∞ case. Numerically, we have set N = 10,000 (N → ∞) and have ensured that these results are converged with respect to N.

This expectation is confirmed in Figure 4(a), which compares the cavity-modified chemical kinetics at N = 1 (blue solid line) and N = ∞. The coupling to the cavity radiation field leads to a resonant modification of chemical kinetics at N = 1 and at the same time, this modification vanishes for N = ∞ in the mean-field limit (red solid line). Note that here we properly include the dipole self-energy terms (see Eq. (11)) whose inclusion within the light–matter Hamiltonian has been a subject of ongoing debate [34], [66–70].

Figure 4:

Cavity-modified chemical reactivity in Model I in the thermodynamic limit at T = 300 K obtained using the constant ρ light–matter coupling Hamiltonian. (a) Cavity-modified (normalized) chemical rate constant as function of photon frequency at N = 1 (blue solid line) and at N = ∞. (b) Cavity-modified (normalized) chemical rate constant as function of photon frequency at N = ∞ but in the absence of the dipole self-energy term.

Figure 4(b) presents the cavity-modified chemical dynamics in the absence of the dipole self-energy (DSE) terms. In the absence of DSE terms, we do observe a large collective cavity modification of chemical kinetics for N → ∞. Specifically, we observe a suppression of the chemical kinetics, which is also cavity frequency dependent. However, we do not observe a resonant effect, where the cavity suppresses chemical reactivity strongly at a certain photon frequency. This observed effect can be explained in terms of the modification of the reaction barrier in the absence of the DSE terms, as explained in Ref. [70]. We observe that this effect increases with the increase in the light–matter coupling strength η _c . Note that, in comparison to Ref. [70], here we explicitly simulate the dynamics of the solvent degrees of freedom.

In Figure 4, the initial state of the system has been prepared in thermal equilibrium with the molecules placed in the reactant well. Next, we explore the cavity modification of chemical dynamics, where the cavity radiation field is in thermal equilibrium in the absence of the light–matter coupling, which is introduced at t = 0. Such an initial condition is naturally nonequilibrium and shows different short time dynamics when compared to the dynamics starting with an equilibrated cavity field.

Figure 5 presents the time-dependent reactant population dynamics when starting from a nonequilibrium initial condition in the thermodynamics limit. The population dynamics presented in Figure 5(a) show a short-time (sub-picosecond) relaxation, which exhibits a photon frequency dependence. When the photon frequency is close to the molecular vibrational transition, the relaxation is faster compared to when the photon frequency is off-resonant. The suppression of the reactant population (or enhancement of chemical reactivity) observed here is due to fact that the cavity radiation is “hotter” than the molecular subsystem at the initial time and this excess heat results in enhanced reactivity at short times.

Figure 5:

Cavity-modified chemical dynamics in Model I under nonequilibrium (uncorrelated) initial conditions at T 300 K. Results were obtained with the Mean-Field HEOM approach for the constant ρ light–matter coupling Hamiltonian. (a) Time-dependent reactant population at various cavity photon frequencies ω _c with ω _c = 1185 cm⁻¹ as the resonant frequency. (b) Reactant population at 1 ps as a function of photon frequency, note the size of the effect.

At longer times, due to the presence of the solvent bath, the cavity radiation mode as well as the molecular system thermalize, and as a result this leads to chemical dynamics that is same as when starting from a thermalized initial condition. Consequently, the long-time chemical rate constant (compare the slopes of the three curves t > 1000 fs) shows no cavity frequency dependence and is identical to the results in the absence of the cavity.

Figure 5(b) presents the reactant population at t = 1 ps at various photon frequencies. The reactant population shows a sharp cavity resonance feature, signifying a collective and resonant cavity enhancement of chemical reactivity. However, the effect is exceedingly small. The exceedingly small magnitude of this effect reflects the closeness of the nonequilibrium initial condition used here to the equilibrium, thermalized density matrix. While the effects shown in Figure 5(b) are very tiny, they do suggest the tantalizing possibility of modifying chemical reactivity under nonequilibrium scenarios, a topic we return to below and in the conclusions.

3.3 Cavity-modified proton transfer reaction

Below we explore cavity-modified dynamics in a model molecular system describing a proton transfer reaction in thioacetylacetone (see details in Section 2.1.2). For the following results presented, we define a dimensionless light–matter coupling strength of η = 0.05, related to the light–matter interaction by

where μ 10 = ψ 0 μ ̂ ψ 1 = 0.042 a.u. [19]. In all calculations for Model II, we make use of the constant ρ form for the light–matter coupling Hamiltonian (Eq. (11)). This corresponds to a Rabi splitting of ≈ 12.59  cm⁻¹ at the resonant photon frequency.

Here, we aim to explore the importance of appropriate initial conditions for the N molecules and the cavity radiation mode when exploring collective VSC. We consider the dynamics of N molecules interacting with a cavity radiation mode that is resonant with the first vibrational transition of the molecular Hamiltonian, H ̂ M and explore the dependence of the dynamics on N for uncorrelated and correlated (here at T = 0 K) initial conditions.

In Figure 6, we compare the time dependence of the side-side correlation function defined in Eq. (31) calculated in the absence of a cavity to that with a cavity mode for various values of N = 1, 2, 4, 8, 16, and 32. Here, we have prepared an uncorrelated ground state as defined in Eq. (30) that is out of equilibrium (with the light–matter coupling introduced at t = 0) and consider a scenario where all the molecules are aligned along the cavity radiation polarization direction. We observe that the presence of a strongly coupled cavity significantly perturbs the coherent crossing dynamics of the molecular system. Figure 6 shows a significant enhancement of the crossing dynamics indicated by the much more rapid decay of the correlation function through the first 100 fs, as well as a more substantial transfer of population at later times. These results are consistent with those observed in reference [19] and the molecular dynamics simulations presented in reference [27]. Further, we observe that upon increasing the number of molecules, this dynamics stabilizes, with minimal deviations observed between the results obtained with N = 16 and 32, even up to times of 1000 fs. Given the convergence of the results with respect to N, we conclude that the observed modification of chemical dynamics occurs in the collective (or in the thermodynamic limit, N → ∞) regime when using a nonequilibrium initial condition.

Figure 6:

Cavity-modified proton transfer dynamics in a model thioacetylacetone system obtained using ML-MCTDH. The upper panel schematically illustrates how N molecules are coupled to a radiation mode. The bottom panel presents the time-dependent reactant population when starting with a nonequilibrium (uncorrelated) initial condition at various N. Light–matter introductions were included using the constant ρ light–matter coupling Hamiltonian.

A number of previous works have shown that static and dynamic disorder can play a crucial role in the cavity-modified dynamics of molecules and materials [71–78]. Here, in Figure 6, we show that the nonequilibrium effects observed when considering the uncorrelated initial conditions decrease significantly in the presence of angular disorder (see gray solid lines) where molecules are randomly oriented with respect to the cavity radiation polarization direction. We observe that the disordered N = 32 scenario shows negligible modification of the mean crossing dynamics when coupled to the cavity radiation, even when starting from a nonequilibrium initial condition. The reason for the disappearance of the cavity modification can be understood by considering the initial (expected) value ⟨ e ⋅ ∑ i μ ̂ i ⟩ = 0 , which brings the initial nonequilibrium state close to the correlated ground state, thereby reducing the nonequilibrium effect seen in the crossing dynamics.

Similar to the dynamics in the disordered scenario, the modification of reactivity becomes small when starting from a correlated initial condition corresponding to equilibrium in the reactant well. Interestingly, while the cavity modification of the crossing dynamics is much smaller in this scenario, we do observe a decay of the oscillation amplitude that exhibits a resonant feature.

Figure 7 presents the time-dependent reactant population at three different cavity photon frequencies when starting from a correlated initial condition. We observe that when the cavity photon frequency is in resonance with the vibrational transition, the amplitude of the crossing dynamics is modified. In Figure 7(a) where the amplitude of the crossing dynamics is largely unmodified for the photon frequency ω _c = 60 cm⁻¹ (blue dashed line) or 180 cm⁻¹ (blue solid line). However, when the photon frequency is close to the vibrational transition of the molecular subsystem, at ω _c = 129 cm⁻¹, we observe a decay of the amplitude of the side-side correlator.

Figure 7:

Cavity-modified proton transfer dynamics in a model thioacetylacetone system when starting from a correlated initial condition. (a) Reactant population dynamics at three different photon frequencies with ω _c = 129 cm⁻¹ as the resonant photon frequency. (b) Reactant population dynamics at various N compared to the no coupling scenario. (c) Reactant population decay as a function of the photon frequency obtained for N = 8. Results were obtained using ML-MCTDH and using the constant ρ light–matter coupling Hamiltonian.

In Figure 7(b), we confirm that this decay persists in the collective limit by comparing the N = 1 and N = 32 cases. Figure 7(b) shows that the overall dynamics are nearly identical between the two scenarios.

To analyze the resonant behavior of the cavity-modified dynamics, we compute a decay rate, κ, obtained by minimizing the function

where a, b, κ are the fitting parameters, θ ₀(t) is the normalized side-side correlation function without the cavity, and θ(t; ω _c ) is the side-side correlation function obtained in the presence of a cavity mode.

Figure 7(c) presents the decay rate κ as a function of the photon frequency ω _c . Overall, the decay rate clearly shows a resonant peak with ω _c close to the vibrational transition of the molecular system. This decay can be understood as a decoherence process due to the presence of the cavity radiation mode, which is acting as a bath degree of freedom. Note that this decay is transient, and at longer times, the populations will eventually return since this is still a small, closed quantum system. At the same time, we do expect this decay to persist when considering a lossy cavity mode.

While such cavity-modified dynamics are collective and resonant, whether or not such effects can be substantial when considering solvent interactions with the molecular system remains an open question, which we briefly return to before concluding. Further, we also note that the damping of the crossing dynamics does not necessarily indicate a modification of chemical reactivity, since the average product population (when averaging over the oscillations) remains the same for the time-range presented here.

The results in Figures 6 and 7 also shed light on reports on the presence or absence of collective cavity modification of chemical reactivity in recent work [15], [27], [79]. Specifically, a number of studies employing classical trajectory-based simulations [27], [79], classical rate theories [15], and quantum dynamical simulations [19] have obtained conflicting results concerning vibrational strong coupling in the collective coupling regime. In Ref. [19], [27], cavity modification of chemical dynamics is observed that persists to large N, and in contrast, in Ref. [15], [79], no modification of the chemical reaction rate is observed in the large N limit. Our results illustrate that the collective effects observed in previous works arise from the choice of the initial condition in which the N molecules and cavity mode are initially uncorrelated. We further demonstrate that such collective effects vanish upon the use of an appropriate correlated initial condition consistent with the linear response definition of rate theory. While we have demonstrated this in the zero-temperature case, the underlying argument holds for finite temperatures as well.

Finally, to explore the origins of these deviations with respect to correlated and uncorrelated initial conditions, we consider the dynamics of the cavity mode itself. In Figure 8, we present the time-dependence of the dipole-gauge cavity number operator n ̂ c = b ̂ c † b ̂ c (note that this does not correspond to the number of photons in the cavity mode). For the case of the correlated (or equilibrated) initial conditions, ⟨ n ̂ c ⟩ increases linearly with the number of particles, and essentially no dynamics are observed in ⟨ n ̂ c ⟩ , over the timescales considered. In contrast, for the uncorrelated initial conditions with the cavity initially prepared in the vacuum state, significant coherent oscillations are observed, consistent with the dynamics of a displaced harmonic oscillator. For all times t considered here, we observe a deviation in the expectation value from that of the correlated initial state case, which scales with N. As a consequence, the mean-field value for the cavity–molecule interaction observed for the correlated initial state and the vacuum initial states differ by a factor of N . This factor exactly cancels the 1 N scaling of the cavity–molecule interaction strength and so in the N → ∞ limit, the cavity modifies chemical dynamics in the collective regime when using uncorrelated initial conditions.

Figure 8:

The time-dependent dipole-gauge cavity number operator expectation value ⟨ n ̂ c ⟩ , for the correlated (dashed) and uncorrelated (solid lines) initial conditions for various N. Results were obtained using ML-MCTDH and using the constant ρ light–matter coupling Hamiltonian.

4.1 Single mode approximation

Most theoretical works investigating the VSC modification of chemistry assumes the coupling of molecules to a single cavity radiation mode [10], [12], [13], [22], [27], [36], [79], [81], [82]. In reality, there are an infinite set of cavity radiation modes that also have a characteristic dispersion. Recent work shows that going beyond the single-mode approximation is necessary to capture various effects in light–matter hybrid systems [83–85]. In the future, we will explore setups where an ensemble of cavity modes is coupled to an ensemble of molecules.

4.2 Long-wavelength approximation

Within this approximation, the spatial variation of the radiation field is ignored. The spatial variation also plays a crucial role in the description of cavity-modified transport in materials and molecules [75], [84], [86], which may potentially impact the modeling of chemical reactivity. In the near future, we will explore how the spatial variation of the radiation field impacts the chemical reactivity in the VSC regime.

4.3 Interactions between molecules

Another approximation employed within our present simulation is that we ignore the interactions between the molecules, and they only interact via the dipole self-energy term. In previous work, it was found that cavities can modify chemical reactivities in the collective regime when coupling the rest of the molecules to a single reactive molecule [14]. In contrast, here we ignore the (spatially dependent) interactions between molecules. Future work will be devoted to exploring how intermolecular interactions may play a role beyond such extreme scenarios. It should be noted that the three aspects described above are necessary to enable and describe polariton transport in solids [75]. In addition, the spatial scale of interactions draws into question the validity of the mean-field limit, which explicitly assumes no spatial scale of interactions. One, however, may still employ cluster mean-field techniques [87]–[89] to attempt to build in this spatial dependence.

Overall, our results point to the possibility of modifying chemical dynamics by coupling molecular vibrations under nonequilibrium conditions and, at the same time, indicate the inability of simple models of collective VSC to exhibit modified reactivity in the equilibrium limit. The physical factors missing in our simple models that might allow for a more complete understanding of the cavity-modified ground state chemical reactivity observed in recent experiments await further investigation.