Investigating the Collective Nature of Cavity Modified Chemical Kinetics under Vibrational Strong Coupling

In this paper we develop quantum dynamical methods capable of treating the dynamics of chemically reacting systems in an optical cavity in the vibrationally strong-coupling (VSC) limit at finite temperatures and in the presence of a dissipative solvent in both the few and many molecule limits. In the context of two simple models we demonstrate how reactivity in the {\em collective} VSC regime does not exhibit altered rate behavior in equilibrium, but may exhibit resonant cavity modification of reactivity when the system is explicitly out of equilibrium. Our results suggest experimental protocols that may be used to modify reactivity in the collective regime and point to features not included in the models studied which demand further scrutiny.

Cavity-modified chemical reactivity in the VSC regime has two distinguishing characteristics.The first is that the chemical rate is only strongly modified (enhanced or suppressed) when the photon frequency is close to some characteristic molecular vibration frequency.This is marked by a sharp peak (or dip) in the cavity-modified chemical rate constant as a function of the cavity photon frequency, with the width of the rate profile matching the width of the IR absorption profile [1][2][3][4][5][6][7][8].The second is that such cavity modifications operate in the collective regime, where a macroscopic number of molecular vibrations are collectively coupled to the cavity radiation.Consequently, the light-matter coupling between each individual molecules and the cavity radiation is effectively minuscule.Thus, a crucial question arises: Can collective light-matter coupling, which couples cavity radiation and molecules in a delocalized fashion, lead to a modification of chemical reactivity which operates locally?
Our previous theoretical work [10], corroborated by recent work [21,23] demonstrates that a sharp resonant modification of cavity-modified chemical rate stems * drr2103@columbia.edufrom the quantum dynamical interplay between the cavity photon mode and molecular vibrations.However, these works operate at the single molecule level, where an individual molecule is assumed to strongly couple to the cavity radiation mode.This is achieved by artificially scaling the single molecular coupling by √ N (where N is the number of molecules coupled to cavity radiation in an experiment) such that the Rabi splitting observed in the single molecule-cavity setup is similar to that of the experiments.While such large single molecular coupling may be achieved in plasmonic cavities [26], this situation is not representative of the present experiments showing cavity-modified ground state chemistry [1][2][3][4][5][6][7][8] that operate in the collective regime.
In this work, with fully quantum dynamical simulations, we investigate the role of collective cavity coupling in the cavity-modified chemical kinetics of two different model systems under various initial conditions.Specifically, we explore two initial conditions that correspond to two significantly different physical scenarios.One is when the initial condition is uncorrelated, where the N molecules are thermalized in the absence of the cavity and at time t = 0 molecule-cavity interactions are introduced.The other is when the initial condition is correlated, where the N molecules and the cavity mode are thermalized in the presence of a dissipative environment, which at zero temperature corresponds to the polaritonic ground state.
Our results indicate that under uncorrelated initial conditions, the non-equilibrium decay of state population can be significantly modified when collectively and resonantly coupling molecules to the cavity mode, which corroborates previous works that operate in the classical regime [27] or those that operate at zero temperature but treat light-matter interactions quantum mechanically [19,29] 1 .We find that under such circumstances, 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,28] nonequilibrium relaxation dynamics is modified, pointing to a possible route to achieving mode-selective chemistry inside optical cavities.On the other hand, under equilibrium (correlated) initial conditions, we find that chemical reaction rates associated with a barrier crossing process, can be substantially modified only in the few-molecule limit, where the other molecules provide an effective source of cavity dissipation, leading to an enhancement of chemical reaction rate in the energy diffusion-limited regime.This modification is observed when the cavity mode is resonantly coupled to molecular vibrations.However, we find that this effect is negligible in the meanfield N → ∞ limit.At the same time, we find that chemical reactions that occur via direct nuclear tunnelling can be modified resonantly where it is found that the coherent oscillations are resonantly damped when coupled to the cavity.Overall, these results point to the possibility of modifying chemical dynamics by coupling molecular vibrations under non-equilibrium conditions while narrowing down the number of factors that might govern the behavior of cavity-modified ground state chemistry.
This paper is organized as follows: In Section II we describe the model systems considered in this study and the quantum dynamics methods that we used in this work.In Section III we present and discuss our numerical results.Finally, in Section IV we document the conclusions of this work and provide avenues for future investigations.

A. 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.Fig. 1a (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 Fig. 1a-b.
The Hamiltonian describing the set of molecular vibrations coupled to cavity radiation is written as where ĤM is the matter Hamiltonian, ĤRM is the radiation-matter coupling Hamiltonian and ĤC is the Hamiltonian for the cavity radiation mode.In the singlemode approximation of a perfect lossless optical cavity, the cavity Hamiltonian becomes The presence of cavity loss can be modelled with a simple Caldeira-Leggett Hamiltonian [30], where ĤC is written as with qc = ( b † c + bc )/( √ 2ω c ), and with C k and ωk sampled from a spectral density function 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 μi = n i μi , 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, Fig. 1 illustrates two possible strategies to analyze the effect of increasing the number of molecules present in the simulations.The second strategy, illustrated in Fig. 1b, is the choice that has been employed in many recent theoretical works (e.g.[27,28]).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 kept constant.
In the second approach, we vary the number of particles and keep the quantisation volume constant while increasing the number of particles present in the volume.As a result, the average particle density per unit of crosssectional area is constant.In this case, we write where η c is independent of N , as has been done in recent work [18,27,28].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 (11) with ĤM describing the bare molecular Hamiltonian.On the other hand, the Hamiltonian for the constant V approach is written as We note that the inter-molecular 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.

Model I
This model is a multi-molecule generalisation of the cavity VSC system considered in our previous work [10].Each molecular system is described by a symmetric double well potential with the reaction coordinate coupled to a dissipative bath for which well-defined chemical rates can be obtained.In this model, each molecule consists of a single reaction coordinate R that is bilinearly coupled to an infinite set of harmonic modes, representing solvent degrees of freedom.We consider the few-molecule case of N = 1, 2, 3, and 4 molecules, as well as the thermodynamic limit of N → ∞.Here, for simplicity, we will only consider the case where all N molecules are aligned with the cavity polarization direction, that is n i = e. 2he molecular Hamiltonian is given by where the reaction coordinate Hamiltonian is given by with V ( R) describing a simple quartic potential where ω b and E b are the barrier frequency and height, respectively.Equivalently, the reaction coordinate Hamiltonian ĤR can be represented using its vibrational states, where {|v i ⟩} are the vibrational eigenstates of the molecular Hamiltonian ( ĤR |v i ⟩ = E i |v i ⟩).In the second line, we have introduced localized states , with an energy Ē0 = 1  2 (E 0 + E 1 ) and a coupling ∆ = 1  2 (E 1 − E 0 ).These states are the localized ground states of the left and the right wells (blue and red wave functions in Fig. 3a), respectively.
The solvent Hamiltonian is given by where the frequencies and coupling constants are determined by the spectral density of the bath, here taken to be of the Debye form ).We note that the quantum dynamical approach used here can be extended to simulate a more realistic molecular system by obtaining the spectral densities describing the bi-linear system-bath couplings and beyond (such as square-linear couplings) and the vibrational levels in Eq. 16 using ab initio approaches [38][39][40].Finally, in this model, the molecular dipole is taken as μ = R.

Model II
This model describes a two-dimensional hydrogentransfer model constructed for thioacetylacetone, developed by Doslic et al. [41] which was recently explored in the context of cavity-modified chemical reactivity in a recent work [19].Explicitly, the molecular Hamiltonian for this two-dimensional model with a reaction coordinate q (describing H-transfer) and a spectator mode Q is written as where m q = 1914.028a.u. and m Q = 8622.241a.u.are the corresponding masses.Further, V (q, Q) is the interaction potential written as where ω Q = 0.0009728 a.u. is the spectator mode frequency and λ s (q) = a q q2 + b q q3 with a q = 0.794 a.u. and b q = −0.2688a.u.Here V 0 (q) is the reaction path potential and is written as, where V j (q) = 1 2 m j ω 2 j (q − q0 j ) 2 + ∆ j with j ∈ {OH, SH} and K(q) = k c exp −(q − q 0 ) 2 with k c = 0.15582 a.u. and q 0 = 0.2872 a.u.The remaining parameters are tabulated in Table I.Following Ref. [19], we consider a restricted one-dimensional reaction path Hamiltonian in which we restrict dynamics to the 1D path obtained by minimising V (q, Q) with respect to Q at each value of q.Following Ref. [19], the dipole moment operator for this one-dimensional model is taken as μ(q) = µ 0 +µ 1 (q− q 0 ) with µ 0 = 1.68 a.u and µ 1 = −0.129a.u.The reaction coordinate degrees of freedom are discretised using a Colbert-Miller DVR [42], containing 120 grid points, over the range q ∈ [−1.5, 2.1]a 0 .The cavity mode is discretised using a harmonic oscillator number basis containing up to 5000 states depending on the number of molecules, N , coupled to the cavity.

Hierarchical Equations of Motion Approach
We simulate exact quantum dynamics of the molecular subsystems, each of which are coupled to a dissipative bath, using the Hierarchical Equations of Motion (HEOM) approach.This well-established open-quantum system dynamics method provides an exact description of the dynamics of a quantum system that is linearly coupled to a set of N harmonic baths [43][44][45].The details of this method can be found in Ref. [43][44][45].In this work we use a customized version of the HEOM approach which is documented in the appendix of our recent paper [10].We use this approach to simulate the quantum dynamics in Model I when considering cases with N = 1 and 2. For N = 3 and 4, direct solutions of the HEOM is unfeasible, and in this regime we made use of a Multi-Layer Multiconfiguration Time-Dependent Hartree (ML-MCTDH) based solver for HEOM [46][47][48].Finally, for N → ∞ we use the mean-field approach described below.
Mean-Field Quantum Dynamics.With increase in N , it becomes impractical to perform a direct HEOM calculation of the dynamics of all molecules and their dissipative baths.Instead, we develop a mean-field approach by using the fact that the coupling between each molecule and the cavity mode tends to zero as 1/ √ N .We follow Ref. [49] in developing this approach for cavity modified chemical dynamics in the limit N → ∞, where this approach is expected to be valid [49][50][51].Numerically, we have set N = 10000 for obtaining the dynamics.We have verified that the results presented here are converged with respect to N .
Within this mean-field treatment we assume that the total density operator of the system comprising N molecules and a single cavity mode can be written in a factorised form for all times Here ρc (t) and ρi (t) are the time-dependent density matrices for the cavity radiation mode and the ith molecule, respectively.The resultant mean-field equations of motion are We observe that for the models of the type considered here, in the mean-field limit, we have a set of N + 1 coupled equations of motion for each of the individual density operators.Next, for simplicity we ignore any static disorder in the molecule dipole operators and assume that we start in the initial thermal mean-field state.With these assumptions, we have ρi (t) = ρM (t) for all values of i.We thus reduce this set of N + 1 coupled equations to a set of two coupled equations: where We approximate the dynamics generated by these ODEs by linearising the functions f c (t) and µ M (t) at each time point.Upon doings so, we see that the evolution at each time point is solely the quantum dynamics under a time-dependent Hamiltonian.Noting that the time-dependent terms only act on the system degrees of freedom in the system-bath model of the molecular Hamiltonian, we apply the HEOM approach to obtain the equations of motion for evolving each of the density operators.
In Appendix A, we explore the convergence of the mean-field approximation upon increasing the number of molecules, N , coupled to the cavity.These results show that the deviation between the mean-field and description of model I and allow for exploration of the accuracy of the mean-field treatment for large but finite N .
Evaluation of the Forward Reaction Rate.In order to evaluate the forward reaction rate, we have assumed that the system is initially in the reactant region with an initial density operator ρ(0) = ρR .The time-dependent reactant and product populations may be written as where ĥ is the side operator that projects onto the reactant states, and is only a function of the reaction coordinate position operator, R, in our work.If first-order kinetics provides a valid description of the reaction process, then in the long-time limit the reactant and product populations will evolve according to the kinetic equations [52][53][54] ṖR (t) = −κP R (t) + κ ′ P P (t), where κ and κ ′ are the forward and backward rate constants, respectively (and are related by κ⟨P R ⟩ = κ ′ ⟨P P ⟩, where ⟨P R ⟩ and ⟨P P ⟩ are the equilibrium reactant and product populations, which can be obtained from the steady state solution of the HEOM [55]).Rearranging the expression for the forward rate constant, we have [52,54,56] where the limit t → ∞ indicates that the kinetic description of the reaction process is only valid after some initial transient process.
In the evaluation of the rate constants, we have considered an initial condition in which the system of N reaction coordinates and a cavity mode are initially uncorrelated from their dissipative baths.Here we have taken an initial density operator of the form where Z R = Tr e −β ĤS /2 (1 − ĥ)e −β ĤS /2 , which allows for the direct application of the HEOM approach.The short-time transient dynamics depends on the choice of the initial reactant density operator.For the noninteracting molecule-solvent initial condition, a shorttime (∼ 200 fs) transient slippage in the population of the reactant state is observed.However, we have found that the long-time plateau value of Eq. 28 is independent of the choice of initial density operator for the models considered in this work.

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 employ an efficient projector splitting integrator [57][58][59][60][61][62][63] that avoids many of the issues that arise in the presence direct product state initial wavefunction conditions [63].In all calculations we used a single-site with subspace expansion integration scheme similar to those proposed in reference [64] for the standard ML-MCTDH algorithm to allow for growth of the bond-dimension (or number of singleparticle 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 Fig. 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.
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.
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], where ψ gsi is the ground state of the molecule Hamiltonian, Ĥ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 thermalised 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 optimisation algorithm presented in Ref. [65], however, with the use of subspace expansion allowing for adaptive control of bond dimension throughout the optimisation.All simulations presented here were performed using the tree tensor network library ttns lib [66].
We report the single-molecule side -total side correlation function 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 ⟩ = N j=1 Θ j |ψ⟩, and evaluating the matrix elements at each point in time.

III. RESULTS AND DISCUSSION
Few Molecule Limit.We first discuss how a cav-  1a).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.
ity 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 Fig. 3, we present the frequency dependent reaction rate, normalised 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 reaciton rate (blue solid line in Fig. 3c) 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 thermalalization, 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 Eqn. 3.
In Fig. 3c-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 Fig. 3c by keeping the per molecule coupling to the cavity mode a constant, which would correspond to keeping V fixed as illustrated in Fig. 1b.In Fig. 3c the cavity enhances the chemical reaction rate by ∼25% at N = 1 (Rabi-splitting of ≈ 26.58 cm −1 ).At N = 2, 3 and 4 ( with the Rabi splittings of 37.541, 45.95 and 53.0 cm −1 ), 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, it saturates for a small number of molecules coupled to the cavity.Fig. 3d 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 Fig. 3d.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.

k (𝜔
Fig. 3e-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 −1 .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 Fig. 3e).As before, whether an increase in N will lead to an increase in chemical reactivity also depends on the cavity lifetime.
Fig. 3f 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. Thermodynamic Limit.We explore the N → ∞ using the mean-field HEOM approach explained in Sec.II B 1. As hinted at from the trends observed in Fig. 3, we expect no cavity modification in the N → ∞ case.Numerically, we have set N = 10000 (N → ∞) and have ensured that these results are converged with respect to N .This expectation is confirmed in Fig. 4a, 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 selfenergy terms (see Eqn. 11) whose inclusion within the light-matter Hamiltonian has been a subject of ongoing debate [34,[67][68][69][70][71].
Fig. 4b 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. [71].We observe that this effect increases with the increase in the light-matter coupling strength η c .Note that, in comparison to Ref. [71], here we explicitly simulate the dynamics of the solvent degrees of freedom.
In Fig. 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 lightmatter coupling which is introduced at t = 0.Such an initial condition is naturally non-equilibrium and shows different short time dynamics when compared to the dynamics starting with an equilibrated cavity field.Fig. 5 presents the time-dependent reactant population dynamics when starting from a non-equilibrium initial condition in the thermodynamics limit.The population dynamics presented in Fig. 5a 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 offresonant.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.
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.
Fig. 5b 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 non-equilibrium initial condition used here to the equilibrium, thermalized density matrix.While the effects shown in Fig. 5b  they do suggest the tantalizing possibility of modifying chemical reactivity under non-equilibrium scenarios, a topic we return to below and in the conclusions.
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 Sec.II A 2).For the following results presented, we define a dimensionless lightmatter coupling strength of η = 0.05, related to the lightmatter 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 (Eqn.11).This corresponds to a Rabi-Splitting of ≈ 12.59 cm −1 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, ĤM , and explore the dependence of the dynamics on N for uncorrelated and correlated (here at T = 0 K) initial conditions.
FIG. 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 non-equilibrium (uncorrelated) initial condition at various N .Light-matter introductions were included using the constant ρ light-matter coupling Hamiltonian.
In Fig. 6, we compare the time-dependence of the sideside 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 Eqn. 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.Fig. 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 stabilises, 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 non-equilibrium initial condition.
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 [72][73][74][75][76][77][78][79].Here, in Fig. 6 we show that the non-equilibrium 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 non-equilibrium 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 non-equilibrium state close to the correlated ground state, thereby reducing the non-equilibrium 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.Fig. 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 Fig. 7a where the amplitude of the crossing dynamics is largely unmodified for the photon frequency ω c = 60 cm −1 (blue dashed line) or 180 cm −1 (blue solid line).However, when the photon frequency is close to the vibrational transition of the molecular sub-system, at ω c = 129 cm −1 we observe a decay of the amplitude of the side-side correlator.
In Fig. 7b we confirm that this decay persists in the collective limit by comparing the N = 1 and N = 32 cases.Fig. 7b shows that the overall dynamics are identical between the two scenarios.
To analyze the resonant behavior of the cavity modi- !(cm -1 ) Time (fs) Decay Rate (34) where a, b, κ are the fitting parameters, θ 0 (t) is the normalised side-side correlation function without the cavity and θ(t; ω c ) is the side-side correlation function obtained in the presence of a cavity mode.
Fig. 7c 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 Fig. 6 and Fig. 7 also shed light on reports on the presence or absence of collective cavity modification of chemical reactivity in recent work [15,27,80].Specifically, a number of studies employing classical trajectory-based simulations [27,80], 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 , in contrast, in Ref. [15,80] 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 demonstrate 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 Fig. 8, we present the time-dependence of the dipolegauge cavity number operator nc = b † c bc (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 essesntially 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 con- sidered 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.

IV. SUMMARY AND CONCLUSIONS
In this paper, we use exact quantum dynamics methods to look into the collective nature of cavity-modified chemical reactivity in the vibrational strong coupling regime.In particular, we investigate how a cavity can modify chemical reactivity when coupling a set of identical molecules to cavity radiation mode in several simplified models.Specifically, we consider two model molecular systems, which we call Model I and Model II.Model I describes a set of molecules embedded in a dissipative environment which are also coupled to a cavity radiation mode.Each molecule is described by a barrier crossing dynamics that is described with a double-well potential.Model II describes a proton transfer reaction with a reaction coordinate coupled to a spectator mode.
We simulate the cavity-modified chemical dynamics in the few-molecule limit with our customized version of the HEOM approach (for N = 1, 2) and through the use of a ML-MCTDH based HEOM solver (for N > 2).From the point of view of one reactive molecule, the other molecules can be seen as a dissipative bath connected to a cavity mode.The addition of several molecules renders dissipation stronger overall, which in turn can lead to an enhancement of chemical reactivity, depending on the molecule-solvent couplings and cavity lifetime.
To simulate the quantum dynamics of this system as N → ∞, we develop a mean-field HEOM approach for addressing the thermodynamic limit.In comparison to recent work that allows for simulating polaritonic quantum dynamics for zero temperatures [29], our mean-field HEOM approach simulates dissipative dynamics at finite temperatures in the presence of solvent degrees of freedom.Using this approach, we find the dissipative effect originating from the other molecules becomes negligible when N → ∞.These results agree with a previous work that treated all degrees of freedom classically [15], as well as the results of Ref. [29].
In this work, we also explore the importance of the choice of initial condition on the short to intermediatetime dynamics of a set of N molecules interacting with a cavity mode.We find that whether or not coupling to a cavity modifies the dynamics of an individual molecule in the limit that N → ∞, depends explicitly on the initial condition.For systems prepared in uncorrelated states, deviations between the inside cavity and out-of-cavity dynamics are observed to persist for large N , however, for correlated initial conditions no such deviations are observed.
Overall, these results indicate that a non-equilibrium effects can lead to both collective and resonant modifications of chemical reactivity when coupling to cavity radiation.This finding should prompt the consideration of explicitly out of equilibrium techniques such as Floquet or pulsed laser methods, to induce non-equilibrium steady-states where altered long-time reactivity might be observed.The mean-field HEOM approach developed in this work is particularly well-suited for investigations.In this regard, we note a recent work [81] which shows that VSC can alter chemistry via a nonequilibrium preparation of initial state.On the other hand, in the context of interpreting experiments demonstrating collective VSC modification of adiabatic chemical reaction rates, our findings indicate that terms that are not incorporated in the models studied here might have relevance for the experimentally observed phenomena.Below, we list a few approximations whose role in cavity-modified chemistry should be investigated in the future.
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,27,36,80,82,83]. 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 [84][85][86].In the future, we will explore setups where an ensemble of cavity modes is coupled to an ensemble of molecules.
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 [76,85,87], which may potentially impact the modelling 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.
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 inter-molecular 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 [76].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 [88][89][90] to attempt to build in this spatial dependence.
Overall, our results point to the possibility of modifying chemical dynamics by coupling molecular vibrations under non-equilibrium 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 in recent experiments await further investigation.Convergence of the Mean-field cavity modified chemical dynamics in Model I under non-equilibrium (uncorrelated) initial conditions at T = 300K.Time-dependent reactant population obtained for various values of N using the numerically exact HEOM approach (with a ML-MCTDH based solver) compared to the mean-field dynamics (black solid line) obtained in the limit that N → ∞ for the constant ρ light-matter coupling Hamiltonian.The cavity frequency is taken to be ωc = 1185 cm −1 .the mean-field population dynamics to the exact HEOM dynamics for this simplified variant of model I for various values of N .As was seen in Fig. 5, a significant deviation, from the out-of-cavity population dynamics, is observed in the mean-field limit.For finite N , the exact HEOM calculations exhibit larger deviations than are observed in the N → ∞ mean-field limit, however, these deviations are observed to decrease with increasing N .These results suggest that the resonant enhancement effect observed in the few molecule regime, arises due to beyond mean-field correlations in the dynamics, that decay in the N → ∞ limit.

CFIG. 2 .
FIG. 2. Multilayer tree used for the N = 16 molecule ML-MCTDH calculations for Model II.Bond dimension (number of single particle functions) are 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 Mi denotes the ith molecules degree of freedom.

ΔFIG. 3 .
FIG.3.Cavity modified chemical dynamics in the few molecule limit of Model I at T = 300K.(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 Fig.1b).(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 Fig.1a).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.

FIG. 4 .
FIG. 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.

1 FIG. 5 .
FIG. 5. Cavity modified chemical dynamics in Model I under non-equilibrium (uncorrelated) initial conditions at T 300K.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 −1 the resonant frequency.(b) Reactant population at 1 ps as a function of photon frequency, note the size of the effect.

FIG. 7 .
FIG.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 −1 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.

FIG. 8 .
FIG.8.The time-dependent dipole-gauge cavity number operator expectation value ⟨nc⟩, 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.
FIG. 9.Convergence of the Mean-field cavity modified chemical dynamics in Model I under non-equilibrium (uncorrelated) initial conditions at T = 300K.Time-dependent reactant population obtained for various values of N using the numerically exact HEOM approach (with a ML-MCTDH based solver) compared to the mean-field dynamics (black solid line) obtained in the limit that N → ∞ for the constant ρ light-matter coupling Hamiltonian.The cavity frequency is taken to be ωc = 1185 cm −1 .

TABLE I .
Parameters for Model II (a.u.)