Non-equilibrium reactive flux: A unified framework for slow and fast reaction kinetics

Reactive processes, such as charge transfer and barrier crossing, often occur on time scales that are much longer than those associated with the ro-vibrational dynamics of the reactants. Rather than following the slow transformation of reactants to products, it is advantageous in such cases to evaluate the reaction rate constant using the flux formalism.^1–15 The latter is based on equilibrium correlation functions that involve the reactive flux, and classical, quantum mechanical, as well as semiclassical, formulations are available. The main advantage of the flux correlation function formalism is that the need to follow the dynamics only up to the “plateau time,” i.e., the time required for initial transients to settle and the reactant population to enter its slow, exponential decay.

There are, of course, many cases of fast reactions characterized by low potential barriers, where there is no clear separation of time scales. In such situations, no plateau regime can be identified, thus the flux correlation function formalism cannot be used to obtain the reaction rate. In some cases of ultrafast reactions, transient dynamics persists almost until the reactant population attains its equilibrium value, causing strongly nonexponential kinetics for the duration of the reaction. These situations require full simulation of the population dynamics. Since one does not know a priori, one would typically attempt to infer the rate by evaluating the flux correlation function, and in the event this function does not appear to plateau, one would abandon this approach and proceed to simulate the evolution of the reactant population.

Further, fully quantum mechanical calculations of condensed phase reactions are impractical, and one has to resort to approximate methods. Quantum-classical approaches are particularly attractive because classical mechanics usually captures the dynamics of the quantum system’s environment with satisfactory accuracy, while offering linear scaling. Several intuitive and efficient quantum-classical approximations are available for simulating state populations. However, equilibrium correlation functions require the evaluation of the Boltzmann operator for all interacting degrees of freedom. Evaluation of this operator in a mixed quantum-classical representation presents a major challenge, making quantum-classical methods not directly suitable to this task.

In this paper, we propose a non-equilibrium, factorized reactant density formulation of the reactive flux, which (at least in the case of electron transfer reactions) addresses both of the above issues. By removing the need for evaluating the Boltzmann density of the total Hamiltonian, this formulation is easily amenable to quantum-classical treatments. In addition, the factorized reactant density formulation yields directly the early time transient dynamics and (when relevant) the rate coefficient associated with long-time exponential decay, allowing characterization of slow or fast reaction kinetics from a single calculation. Factorized forms of the Boltzmann operator have been employed in previous rate calculations¹⁶ as an approximation to the near-equilibrium density.¹⁷ We emphasize that the factorized reactant density formulation of the reactive flux presented in Sec. II is not an approximation and that (if evaluated exactly) it yields the same reaction rate value obtainable through the conventional equilibrium reactive flux formalism.

Section II describes the reactive flux with a factorized initial density and its relation to population dynamics. The characteristics of the non-equilibrium flux are illustrated through model calculations in Sec. III. Some concluding remarks, along with an outlook, are given in Sec. IV.

We denote the reactant and product states collectively as $R$ and $P$⁠, respectively. Miller has shown^6,7 that the (forward) rate constant for reactive processes in gas phase bimolecular collisions is given by the following expression:

where β = 1/k_BT is the inverse temperature, Z_R is the partition function of the reactants, $ĥR=RR$ is an operator that projects onto reactants, and $F^$ is the symmetrized flux operator,

The dividing surface that separates reactants from products is perpendicular to the reaction coordinate, which defines the “system” with Hamiltonian $Ĥsys$⁠, while the remaining degrees of freedom constitute the “environment,” defined by $Ĥenv=Ĥ−Ĥsys−V^int$⁠.

For typical barrier crossing processes in the condensed phase, where the reactant-product complex is characterized by a double well potential (or by two bound diabatic surfaces, as in electron transfer reactions), Eq. (2.1) must be modified to

where t_plateau is the “plateau time,” when the initial transients have died out and the reactant population has just entered the slow, exponential decay. Further, Eq. (2.3) is valid under the assumption that all non-reactive intra-well processes occur on a time scale much shorter than the time scale for completion of the reaction. This is often the case, as typical potential barriers separating reactants and products are significantly larger than the thermal energy. Under these conditions, the plateau time occurs relatively early, such that k_ft_plateau$≪$ 1. On a longer time scale, the function inside the trace of Eq. (2.3) is not constant but decays exponentially.

Rearranging the trace and exploiting the Hermitian character of the flux operator, Eq. (2.3) can also be written in the form

where

Assuming that the coupling V_int between system and environment is diagonal in position, the flux operator acts within the space of the system. Then Eq. (2.5) may be rewritten in the form

where $ρ^eqred(t)$ is the reduced density operator of the quantum system,

and the initial condition of the density is

According to Eq. (2.6), the reaction rate is given by the (negative of) the expectation value of the flux operator in the plateau regime, with the initial density given by Eq. (2.8).

Invoking Onsager’s ideas,¹⁸ we argue that (after the initial transients have settled) the decay of the expectation value of the flux should approach its thermodynamic limit with the same rate, regardless of the initial condition. This invariance has been used to show the equivalence of several commonly used rate expressions that involve different symmetrizations of the Boltzmann-transformed flux.¹⁷ In this work, we consider replacing Eq. (2.8) by a simpler, factorized initial condition, which is physically meaningful and easier to evaluate. In particular, a factorized initial condition would allow a fully quantum mechanical treatment of the reaction coordinate and a classical (or quasiclassical) treatment of the its environment, by means of trajectories sampled from a phase space distribution. The particular form we choose is designed to mimic the early state of the reactive process of interest. We choose a factorized initial condition that describes the isolated reactants,

where $Ĥenvinit$ is the Hamiltonian for the degrees of freedom of the environment in (exact or approximate) equilibrium with the system in the reactant state and $Zenvinit$ is its partition function. In the particular case of an electron transfer reaction, $Ĥenvinit$ is the Hamiltonian for the solvent equilibrated with respect to the quantum state describing the electron donor, thus the factorized density of Eq. (2.9) describes the reactant state exactly in this case. We define the expectation value of the flux subject to this new, non-equilibrium initial density,

It is easy to see that

where

is the population of the reactants with the initial value equal to 1, according to Eq. (2.9). Thus, $Fnon-eq(t)$ is equal to the time derivative of the reactant population at all times. Once the transients die out, the reactive process enters the exponential decay regime, where the reactant population decays according to the form

where $texp$ indicates the onset of the exponential regime in the population dynamics and $kf$ and $kb$ are the forward and backward rates. The time derivative of the reactant population is

The time t_exp may be longer than the plateau time, but (as long as there is a separation of time scales) it is much shorter than (k_f + k_b)⁻¹, thus Eq. (2.14) will plateau. At the onset of the plateau regime, the reactant population has not yet changed much from its initial value, thus ρ_RR(t_plateau) $≃1$⁠. It follows that

where the last equality is a consequence of the detailed balance condition. Thus, if there is a separation of time scales, the expectation value of the flux with the non-equilibrium initial condition corresponding to the reactant density will plateau to the (negative of) the forward reaction rate constant. Thus, in line with the discussion of Ref. 17, the equilibrium and non-equilibrium flux expressions have the same plateau value (and subsequent decay), although they are expected to differ at earlier times. The ability to use a factorized thermal density without affecting the accuracy of the computed rate value is particularly useful in simulations employing quantum-classical methods.

Most importantly, the non-equilibrium flux expression with a factorized reactant density is particularly useful in the case of electron transfer reactions, where Eq. (2.9) provides an exact description of the reactive species at the onset of the chemical process. Even when the lack of a clear separation of time scales causes a late onset of the exponential regime, i.e., the absence of a plateau and/or nonexponential kinetics, the non-equilibrium flux can still be employed to infer the evolving reactant population,

If the transient dynamics survive long enough for the reactant population to drop substantially from its initial value, Eq. (2.16) captures this early nonexponential dynamics faithfully, while the rate constant obtained from Eq. (2.15) (along with the backward rate, which is available through the detailed balance condition) can be used to infer the subsequent population decay through an exponential function, i.e.,

Given the intimate connection between the non-equilibrium flux and the population dynamics, an obvious question is whether one could obtain the same information by computing the population and evaluating its time derivative numerically. The main problem with that approach is that it can be unstable, as numerical derivatives are very sensitive to statistical noise. Since most simulation methods applicable to molecular/condensed phase systems employ Monte Carlo sampling,¹⁹ it is preferable to differentiate the population expression analytically. By contrast, the integration procedure required inferring the population from the non-equilibrium flux is stable; in fact integration tends to wipe out the effects of random noise that may be present in the derivative function.

As mentioned earlier, one of the difficulties associated with numerical evaluation of the equilibrium flux expression is often the need to evaluate the Boltzmann operator in the appropriate representation. For example, if a mixed quantum-classical approximation to the time evolution is adopted, one needs to obtain the Wigner quantized phase space distribution²⁰ for the particles comprising the system’s environment, in equilibrium with the reaction coordinate which should be treated quantum mechanically. In the case of molecular systems with a well-defined harmonic zeroth order Hamiltonian, we have shown that the Wigner function can be obtained approximately from adiabatically switched trajectories;²¹ however, accounting for the interaction between classical and quantum degrees of freedom presents a challenge. On the other hand, the use of a non-equilibrium flux expression with a factorized reactant density is easily amenable to a variety of numerical treatments and is ideally suited to quantum-classical methods. Even when the factorized form is only an approximation of the true reactant density, the procedure yields the exact reaction rate. Additionally, when the density factorization is exact, as in the case of electron transfer, integration of the non-equilibrium flux yields the early population dynamics, which is crucial for a correct characterization of the reaction prior to the onset of exponential kinetics and when a plateau cannot be identified.

We illustrate the ideas presented in Sec. II with several numerical examples on model symmetric two-level systems (TLSs) coupled to harmonic dissipative environments. The TLS Hamiltonian is given by

while the harmonic bath and system-bath interaction have the usual form,

The bath is described by an Ohmic spectral density,²² 

where ξ is the Kondo parameter and ω_c is the “cutoff” frequency, which corresponds to the maximum of Eq. (3.3). The equilibrium flux was evaluated using the quasi-adiabatic propagator path integral²³ (QuAPI) methodology of complex-time flux correlation functions.^24,25 The non-equilibrium flux with a factorized initial condition was obtained by propagating the reduced density matrix using the iterative decomposition of the QuAPI expression.^26–28 The factorized initial condition is an exact description of the reactant species in the case of a TLS.

The parameters of the first example are chosen from earlier work,²⁹ where the TLS coupling corresponds to a tunneling splitting $2ℏΩ$ = 0.001 05 cm⁻¹ and the bath cutoff frequency has the value 500 cm⁻¹. These parameters are characteristics of many proton transfer or isomerization reactions, where a relatively high potential barrier leads to a small tunnelling splitting, while the vibrational frequencies of the environment are much higher. Thus, there is a clear separation of time scales and one expects a well-defined flux plateau.

The equilibrium and non-equilibrium flux functions are compared in Fig. 1. It is seen that both functions plateau within a fairly short time, of the order of the characteristic time scale $ωc−1$ of the degrees of freedom comprising the environment. While the equilibrium flux appears to plateau somewhat earlier, as expected (given the equilibrium initial condition), the flux with a factorized reactant density does not take much longer to reach its plateau. The short-time behavior of the non-equilibrium flux differs somewhat, reflecting the true transient dynamics of the reactant population, and the difference between the two forms increases with system-bath coupling strength. The rate constants obtained from the plateau values of both flux functions agree with the results of Ref. 29.

The second model employs parameters encountered in some ultrafast condensed phase electron transfer reactions. Here the parameters are ω_c = 2.5Ω, $ℏωcβ=0.5$⁠, and ξ = 1.2. Figure 2 shows the non-equilibrium flux as a function of time, along with its time integral, which is seen to be in excellent agreement with the reactant population obtained through direct propagation of the reduced density matrix. The proximity of TLS and bath time scales leads to rapid population decay. Figure 2(b) shows that the flux reaches a large negative value around ω_ct ≈ 1. This region corresponds to an inflection point of the population, which shortly thereafter settles into an exponential decay with a constant rate, as evidenced from the logarithmic plot shown in Fig. 2(c). (Note that the deviation of the integrated flux from the directly computed population observed on the logarithmic plot is a consequence of the numerical error in the trapezoid rule integration of the flux.) However, since the duration of the transient dynamics is not negligible compared to the population decay time, the flux enters its exponential decay regime on the same time scale and thus does not display a plateau.

As a third example, we discuss the quantum-classical path integral^30,31,33 (QCPI) all-atom simulation of the ferrocene-ferrocenium self-exchange electron transfer reaction in liquid hexane.³² This is an ultrafast reaction that completes in a few picoseconds. In that case, the non-exponential transient dynamics persist until the reaction is nearly complete, leading to much faster decay than the rate constant associated with the long-time exponential regime would seem to predict.

One of the flux formulations of reaction dynamics expresses the rate as the cross correlation function of a projector and the flux operator, or (equivalently) as the expectation value of the reactive flux with respect to a symmetrically “Boltzmannized” operator. These expressions involve the Boltzmann operator with respect to all degrees of freedom, whose computation often presents challenges, in particular if mixed quantum-classical methods are to be employed. To remove this obstacle, one would like to replace the Boltzmann operator by a factorized form, such that the quantum system and its environment may be treated at different levels of approximation. However, such a modification naturally leads to a different function.

Since a change of initial conditions does not affect the rate of approach to equilibrium, one may replace the Boltzmann operator by a non-equilibrium form without affecting the plateau value of the flux. As long as the real time dynamics is evaluated in a numerically exact manner, the non-equilibrium flux expression yields the exact value of the rate constant, irrespective of any approximations made to the initial thermal density. This feature offers considerable flexibility, allowing convenient approximations of the Boltzmann density without sacrificing accuracy in the computed value of the reaction rate. This flexibility is particularly valuable in simulations employing quantum-classical methods.

In particular, if the reactant density involves a single quantum state, as in the case of electron transfer reactions, the factorized initial condition describes the true initial state of the reactant species and the flux expression gives precisely the time derivative of the reactant population. Thus, in addition to its convenient factorized initial density, this flux expression yields the short-time transient dynamics and (through its plateau value) the rate of exponential decay, allowing characterization of slow or ultrafast reaction in a unified framework. The examples presented in Sec. III illustrate different types of ultrafast dynamics, where determination of the rate through the equilibrium flux may not be possible (because the flux does not plateau, even though the population follows exponential decay for the most part) or where such a rate may not even be meaningful (because the non-exponential population decays within a time that is much shorter than the inverse of the long-time rate).

As expected, the non-equilibrium flux exhibits more pronounced transients and may settle to its plateau value somewhat slower compared to the equilibrium flux. However, the examples presented in Sec. III show that the plateau time is still reached quite rapidly. This is to be expected even in the case of very slow processes, since the plateau time of the non-equilibrium flux coincides with the onset of exponential population decay, which occurs very early (on the time scale of intra-well dynamics) compared to the time for completion of the reaction. Thus, the use of the non-equilibrium flux expression does not require simulation of the dynamics for very long times.

The non-equilibrium formulation of the reactive flux will facilitate simulations of fast and slow reactions in solution or biological systems. A particularly attractive possibility is its evaluation using the rigorous QCPI methodology,^30,31 which treats the interaction between the quantum system and its classical environment in full detail and without approximation. Applications to charge transfer reactions are in progress.

This material is based on work supported by the National Science Foundation under Award No. CHE 13-62826.