Tests for, origins of, and corrections to non-Gaussian statistics. The dipole-flip model Free

Understanding the chemical dynamics of liquids is critical to the study of chemical reactions, spectroscopy, energy transfer, and catalysis. Linear response theories,^1–5 in which perturbations are described by fluctuations in the equilibrium state, are the dominant method for describing nonequilibrium dynamics. In many linear-response approaches,^1–4 the perturbation is assumed to be weak—on the order of thermal energies (k_BT)—so that higher-order correlations in the equilibrium fluctuations are neglected. Such linear response approximations have been remarkably successful in predicting nonequilibrium dynamical properties for not only weak perturbations, but also for perturbations much larger than those on the order of the thermal energy. It has been appreciated for some time,^6,7 and even shown directly,^5,8 that the latter is a consequence of systems exhibiting Gaussian statistics—which predicts not just local but global, linear response. Thus, Gaussian statistics approximations underlie the linear response approaches that assume small perturbations and represent a more broadly applicable framework for thinking about linear response.

A key application of linear-response approximations has been in describing time-dependent fluorescence (TDF) experiments that are used to study solvation dynamics. In these experiments, a dye molecule with a charge-transfer transition is photoexcited, and the subsequent return to the ground state via fluorescence is followed as a function of time after excitation using femtosecond spectroscopy. Specifically, the fluorescence energy, $ΔE(t)$⁠, or Stokes shift at time t after excitation reports on the rearrangements of the solvent induced by the change in the dye-molecule charge distribution upon excitation. The result is often presented as the normalized time-dependent Stokes shift (TDSS)

where the $⟨⋯⟩ne$ indicates a nonequilibrium average and $⟨ΔE(∞)⟩ne$ is the relaxed fluorescence energy. This experiment can be modeled with nonequilibrium molecular dynamics (MDs) simulations in which the equilibrium ground-state dye molecule is excited, i.e., the charge distribution is changed at t = 0, and the resulting time-dependent fluorescence energy, $ΔE(t)=ΔEe(t)−ΔEg(t)$⁠, is monitored. Here, the subscripts indicate the excited (e) and ground (g) states and $ΔE(t)$ thus represents the difference in the dye-solvent interactions between the two electronic states (i.e., solute charge distributions).

Linear response approximations relate the nonequilibrium response to the fluctuations observed at equilibrium. There are two approaches to obtain linear response approximations which give S(t) in terms of only equilibrium simulations, i.e., in terms of a correlation function of the equilibrium fluctuations of the energy gap $ΔE$⁠. In the first approach, the dynamic linear-response approximation⁹ gives the TDSS in terms of a ground-state correlation function, $S(t)≈Cg(t)=⟨δΔE(t)δΔE(0)⟩g/⟨δΔE2⟩g$⁠, where the angled brackets denote an equilibrium average in the ground (g) state and $δΔE(t)=ΔE(t)−⟨ΔE⟩g$ is the fluctuation in the energy gap at time t from its equilibrium value. This approximation has been used to successfully describe the solvation dynamics of many systems and, indeed, most simulations of TDF measurements reported in the literature have been carried out using this method. The second approach, the static linear-response approximation, gives the TDSS in terms of the excited-state fluctuations, $S(t)≈Ce(t)=⟨δΔE(t)δΔE(0)⟩e/⟨δΔE2⟩e$⁠. This is obtained (vide infra) by assuming that fluctuations of the energy gap at equilibrium are small. Much like the dynamic linear-response approximation, the static linear-response approximation can accurately describe solvation dynamics in many systems. Interestingly, Maroncelli and Fleming found that it was generally superior to the dynamic linear-response approximation in tests against nonequilibrium simulations,¹⁰ a result that was subsequently explained more directly by Laird and Thompson.^5,8 However, in a number of notable systems, neither linear-response approximation is accurate.^11–15 The reasons for these breakdowns are not yet fully understood in terms of the linear-response approximation, making it difficult to know if an approximation is justified absent a full nonequilibrium treatment. It is these issues that are the subject of this paper.

Laird and Thompson have previously shown that the static linear-response approximation, i.e., $S(t)≈Ce(t)$⁠, actually represents the assumption that the energy gap, $ΔE$⁠, obeys Gaussian statistics.^5,8 Where the static linear response approximation ignores higher-order correlation functions (truncating the perturbation after the linear term), the Gaussian statistics assumption factorizes them into lower-order correlation functions. The Gaussian statistics assumption therefore does not rely upon the (usually unjustified) approximation that the energy gap fluctuations are small^1–4,14 and thus it represents a global approximation that is still applicable for systems driven far from equilibrium. We note that the connection between linear response and Gaussian statistics was widely recognized before this explicit relationship was derived (see, e.g., Refs. 6 and 7). For example, Marcus theory invokes linear response to illustrate the mechanism by which electron transfer reactions occur, that is, that solvent fluctuations on a harmonic free energy surface can adequately describe nonequilibrium dynamics.¹⁶ Throughout the rest of this paper, we will refer to the approximation $S(t)≈Ce(t)$ as the Gaussian statistics, rather than static linear-response, approximation except when noting the differences between the two approaches.

In this paper, we address some of the outstanding fundamental questions about nonequilibrium chemical dynamics and the Gaussian statistics approximation. First, can non-Gaussian statistics be predicted by relying only on equilibrium simulations? Second, what is the origin of the non-Gaussian statistics? Third, can the Gaussian statistics approximation be corrected to describe nonlinear response?

We are motivated, in part, by the work of Ladanyi and co-workers, who studied solvation dynamics of model diatomic dye molecules in water and methanol to investigate linear response. They studied two models for the electronic transition. In the first model, excitations induced a dipole moment in the excited state from a nonpolar, ground-state diatomic solute.^11,12 In the second model, excitation induced a “dipole flip” in which the ground-state dipolar molecule charges switch in the excited state, as depicted in Fig. 1.^12,13 In both models, linear-response approximations failed to predict the nonequilibrium dynamics: the nonequilibrium simulations showed faster relaxation times than were computed from equilibrium correlation functions. These linear-response breakdowns were attributed to changes in hydrogen-bonding between the solvent and solute after the excitation. The linear-response approximation accurately described the inertial relaxation at short times but showed significant deviations at longer times associated with molecular reorientations. The latter, dipole-flip, model is particularly interesting in that the ground- and excited-state equilibrium dynamics are the same so that the dynamic linear-response and Gaussian statistics approximations are equivalent. Not only does this model provide a convenient test for the approximations, but because of the breakdown despite the symmetry, it also hints at the way in which linear response fails, as will be discussed throughout this paper.

The remainder of this paper is organized as follows. The linear-response and Gaussian statistics approximations are briefly outlined in Sec. II. The details of the solvation dynamics of the dipole-flip model in various solvents are given in Sec. III. The results of the simulations are then presented in Secs. IV–VI in the context of testing diagnostics for non-Gaussian statistics, uncovering the origins of the non-Gaussian behavior, and evaluating corrections to non-Gaussian dynamics. The present results are discussed in the context of other examples of nonlinear response in the literature in Sec. VII. Finally, some concluding remarks are offered in Sec. VIII.

In this section, we briefly outline the static linear response and Gaussian statistics approximations for the case of a change in solute state, e.g., solvation dynamics. Detailed derivations, including that for dynamic linear response, have been provided elsewhere.^5,8 We consider the change in some property A as a function of time after excitation from the ground, g, state to the excited, e, state. This can be expressed as

where H_g (H_e) is the ground-state (excited-state) Hamiltonian and $L^e$ is the excited-state Liouville operator, $L^e={He,⋅}$⁠, with ${F,G}=∑j[(∂F/∂pj)(∂G/∂qj)−(∂F/∂qj)(∂G/∂pj)]$ the Poisson bracket. Note that the time-dependence of the property A is determined by excited-state dynamics initiated from equilibrium ground-state configurations. The exact (classical) result for $⟨A(t)⟩ne$ can thus be obtained from non-equilibrium molecular dynamics in which ground-state equilibrium simulations are used to generate configurations from which excited-state dynamics are initiated. Then, $⟨A(t)⟩ne$ is calculated as an average of A over the non-equilibrium, excited-state trajectories as a function of time, t, after excitation.

Linear-response approximations, however, yield expressions for $⟨A(t)⟩ne$ in terms of only equilibrium dynamics, resulting in significantly less computational effort. In these approximations, the ground-state Hamiltonian is expressed in terms of the excited-state Hamiltonian and the fluorescence energy, $Hg=He−ΔE$⁠. Then, Eq. (2) can be rearranged as

where A(t) is propagated under the excited-state Hamiltonian (i.e., $A(t)=e−L^etA(0)$⁠). By dividing the numerator and the denominator by the excited-state partition function the dynamic nonequilibrium average of A is

where the angled brackets indicate an equilibrium average and the subscript e signifies that the average is in the excited state. Defining the average energy gap $ΔE¯≡⟨ΔE⟩e$ and property $A¯≡⟨A⟩e$⁠, it is straightforward to show that

where $δΔE(t)=ΔE(t)−ΔE¯$ and $δA(t)=A(t)−A¯$⁠.

So far no approximation has been invoked, but it is at this point that the derivations of the static linear-response and Gaussian statistics approximations diverge.

The static linear-response approximation^1–4 is obtained by expanding the exponential in Eq. (5) in a Taylor series and truncating at the linear term

This approximation can be introduced in both the numerator and denominator and is valid only when the equilibrium fluctuations in $ΔE$ are small relative to k_BT (⁠$βδΔE≪1$⁠). This gives the dynamic Stokes shift as

where the last equality results because $⟨δA(t)⟩e=0$ and $⟨δΔE⟩e=0$⁠.

The mathematical expression in Eq. (7) can, however, be obtained without invoking any approximations regarding the size of fluctuations. If, instead of truncating at the linear term the full Taylor series is included for the exponential, one has

Then, if $δA(t)$ and $δΔE(t)$ are Gaussian random variables, all of the higher-order correlation functions can be factorized into only two-point correlation functions, e.g., $⟨δA(ti)δΔE(tj)⟩e$⁠, according to Wick’s theorem.^17,18 Thus, the correlation function $⟨δA(t)δΔE(0)2n⟩e$ can be factorized into terms of the form $⟨δA(t)δΔE(0)⟩e⟨δΔE(0)⟩e$ and $⟨δΔE(0)2⟩e⟨δA(t)⟩e$⁠. However, because $δA$ and $δΔE$ both have zero mean, then $⟨δA(t)δΔE(0)2n⟩e=0$ for all n when the two variables obey Gaussian statistics. Applying similar factorizations to $⟨δA(t)δΔE(0)2n+1⟩e$ yields^5,8

If $δΔE$ is a Gaussian variable then the denominator in Eq. (4) (which is the characteristic function evaluated at $−iβ$⁠) can be similarly evaluated to give

which gives the result for the nonequilibrium change in the property A within the Gaussian statistics approximation as

This is the same expression as the static linear-response approximation to the Stokes shift in Eq. (7). It is important, however, to note that the physical assumptions used here are quite different from those used in static linear response. In particular, there is no requirement that $βδΔE≪1$⁠. As has been shown before^8,19 and will be demonstrated below, this assumption of the static linear-response approach is not satisfied and thus it is more accurate to consider the expression in Eqs. (7) and (11) as the Gaussian statistics approximation.

Equilibrium and nonequilibrium simulations were carried out for a dipolar solute in bulk solutions of methanol, water, and acetonitrile. Each system consisted of 256 total molecules: 1 solute molecule and 255 water, methanol, or acetonitrile molecules. The solute was modeled as a diatomic molecule with each atom having equal and opposite charges (Fig. 1). The molecules interact through site-site Lennard-Jones and Coulombic potentials; the force-field parameters are given in Table I. The solvent molecules are flexible with the water, methanol, and acetonitrile potentials similar to the similar to the extended simple point charge (SPC/E),²⁰ all-atom optimized potentials for liquid simulations (OPLS-AA),²¹ and Nikitin and Lyubartsev²² models, respectively; the intramolecular potential uses harmonic bond and angle bending potentials using the parameters given in Table II. Long-range interactions were computed using an Ewald sum with a tolerance of 10⁻⁴ and a potential cutoff radius of 12.0 Å. The solute bond length was held fixed using the SHAKE algorithm.²³ 

Two types of MD simulations were carried out: (1) equilibrium trajectories with the solute in the excited state and (2) nonequilibrium trajectories where the solute is promoted to the excited state from initial configurations obtained from an equilibrium ground state simulation and the dynamics are subsequently followed. Note that the dipole-flip solute has the same ground- and excited-state dynamics due to the symmetry of the charges. For each system, five equilibrium trajectories were sampled. Each simulation was propagated for a time between 400 and 600 ps in the NVE ensemble in which velocity rescaling was used to maintain the temperature. Following the NVE period, the system was propagated using the NVT ensemble at a temperature of 298.15 K with a Verlet integrator and a time step of 1 fs. The Nosé-Hoover thermostat was used to maintain temperature with a time constant of 100 fs. An equilibration time of 40–60 ps was followed by a production run of 5 ns for the solute in the ground state. The potential energy was then computed for the solute in the excited state in the configurations of the ground-state trajectory. Data for the simulations were recorded every 100 fs.

The nonequilibrium MD simulations were performed from the configurations produced by the ground state equilibrium trajectory. The solute at a particular configuration was promoted to the excited state and the dynamics in the excited state were followed for 50 ps. The ground-state dynamics were then computed along this excited-state trajectory to give the fluorescence energy, $ΔE(t)$⁠. This process was repeated for a total of 1000 nonequilibrium trajectories, taken from the ground-state equilibrium configurations every 5 ps. Data for the nonequilibrium simulations were recorded every 10 fs. Error bars were computed using block averaging. For the equilibrium simulations, five blocks were used where each block was taken to be an individual equilibrium simulation. Ten blocks were used for the nonequilibrium MD simulations, where each block was composed of 100 individual trajectories.

We start by addressing an important practical issue involved in applying the Gaussian statistics approximation. Namely, how can one be confident that the results are accurate, i.e., the system exhibits Gaussian statistics, without carrying out a full nonequilibrium simulation?

We first show, as has been done before,^6,12,13,24 that the dipole-flip model exhibits a breakdown of the Gaussian statistics approximation for water and methanol solvents. It is first worthwhile, however, to note some general properties of the normalized time-dependent Stokes shift, S(t) for these two systems. The TDSS results obtained from nonequilibrium MD simulations are shown in Fig. 2. For the solute immersed in either solvent, S(t) is well fit by a tri-exponential function. In water, the shortest time scale is 20 fs and represents the inertial motion of the nearby solvent molecules. The intermediate time scale is $∼0.35$ ps and represents the oscillating librational motion of water due to its hydrogen bonds. The longest time scale, which is associated with molecular rearrangements, is $∼1.7$ ps. Methanol has an inertial time scale of $∼65$ fs, a librational time scale of $∼0.64$ ps, and a longest time scale of $∼5$ ps. These relaxation times resulting from the simulations conducted in this study are consistent with those found previously.¹³ 

The difference in relaxation at short times can be attributed to key differences in the solvation dynamics between the two solvents. First, the three moments of inertia for water are all small in comparison to methanol which has one small and two large moments of inertia.^13,24 The larger moments of inertia describe the slower response time of methanol in the inertial regime. Second, and more relevant to the issue of the Gaussian behavior, is the number of solvent-solute hydrogen bonds. Water, because of its smaller size, donates more hydrogen bonds to the negatively charged solute site and thus exhibits a larger amplitude for the librational component of the solvation dynamics compared to methanol. The differences in the longer time scales, associated with molecular rearrangements, are consistent with the faster OH reorientation for water compared to methanol.²⁵ 

For the dipole-flip model considered in this paper, the Gaussian statistics approximation has been implemented for the fluorescence energy, i.e., $A=ΔE$⁠, to approximate the normalized TDSS response function in Eq. (1). This gives

where the subscript “e” has been suppressed on the thermal average due to the symmetry of the model. The full nonequilibrium S(t) and the Gaussian statistics approximation, C(t), for water and methanol are shown in Fig. 3. The data show that Gaussian statistics is an imperfect approximation for the nonequilibrium solvation dynamics of the dipole-flip model in water and a very poor one in methanol. For each system, the relaxation predicted by Gaussian statistics is slower and does not agree (within error bars) with the nonequilibrium dynamics until $t≈0.5$ ps for water and $t>8$ ps for methanol.

Clearly, the error of the Gaussian statistics approximation is greater for methanol than it is for water. However, the larger question is how one could evaluate the approximation in either case without the direct comparison to the exact, nonequilibrium result, as presented in Fig. 3. As we now discuss, the Gaussian statistics approximation outlined in Sec. II provides an answer. Namely, Wick’s theorem states that if $ΔE$ is a Gaussian random variable then

if m is odd, and, if m is even

where n = m/2 − 1 (or m = 2n + 2). Thus, the quality of the Gaussian statistics approximations depends on the accuracy of these factorizations and can be directly evaluated by the comparison of C_m(t) and $Cmfac(t)$ for the higher-order (⁠$m>2$⁠) correlation functions.

The first of these higher-order correlation functions and their factorizations based on Wick’s theorem are displayed in Fig. 4, where the non-normalized Gaussian statistics approximation, $C2(t)=β⟨δΔE(t)δΔE(0)⟩e$⁠, and the first term that is neglected according to Eq. (13), $C3(t)=β2⟨δΔE(t)δΔE(0)2⟩e/2$⁠, are shown. For both water and methanol solvents, the latter, higher-order correlation function is nearly, but not fully, zero within error bars, which is expected for Gaussian statistics. For water it is nonzero up to $t≈0.4$ ps, while for methanol the results suggest that C₃(t) may be non-zero for much longer times, but the (substantial) error bars overlap zero for all times. The next highest, odd-order correlation functions, $C5(t)=β2⟨δΔE(t)δΔE(0)4⟩e/24$ (not shown), are essentially zero for all times for water but non-zero and positive for methanol with a rapid (⁠$<30$ fs) initial drop followed by an $∼1.5$ ps decay.

The even-order time correlation functions are all non-zero and it is the factorization in Eq. (14) that can be checked for deviations from Gaussian statistics. The first example is shown in Fig. 5 where C₄(t) and its factorized form, $C4fac(t)$⁠, are compared for both water and methanol solvents. Similarly, the sixth-order correlation functions, C₆(t), are compared to the factorized forms, $C6fac(t)$⁠, in Fig. 6. Note that the magnitude of these correlation functions increases with the order in direct contradiction of the approximation invoked in static linear response. Interestingly, when water is the solvent, the factorization of these even-order time correlation functions is an excellent approximation, with no discernable difference with the full time correlation functions. In contrast, the results for the methanol case show significant disagreement between the full and factorized forms. For both C₄(t) and C₆(t), the corresponding factorized correlation function is larger than the full version. The difference is larger for $C6fac(t)$ than $C4fac(t)$⁠, but in both cases the factorized form decays slightly faster than the full correlation function such that they are in better agreement at longer times. These results are consistent with the comparison of the nonequilibrium response, S(t), to the Gaussian statistics approximation, C(t), in Fig. 2, where the latter is significantly larger than the former.

Taken together, these comparisons of the full, higher-order time correlation functions to the factorized forms assumed in the Gaussian statistics approximation indicate a strong breakdown of the approximation for methanol and a quite weak one for water. This is generally consistent with the comparison of the Gaussian statistics result with the nonequilibrium response in Fig. 2, where the methanol case shows significant disagreement while for water the differences are smaller and more quantitative.

It is interesting to compare these results to an analysis of the equilibrium energy gap distribution, $P(δΔE)$⁠. This is done in Fig. 7, where the distributions for water and methanol are shown. For both solvents, the results are interesting. In water, the energy gap distribution is well described by a Gaussian, with deviations only observed in the very wings of the distribution. In methanol, however, the statistics of the energy gap are non-Gaussian exhibiting a shoulder on the low energy side of the distribution. These results make prescribing a cause of non-Gaussian statistics difficult: in some cases, like methanol, the breakdown is expressly found in the energy gap distribution, but in others, like water, the breakdown is subtle, or even “hidden,” and not clearly evident in the distribution.²⁶ This was previously observed for the dipole-flip model in water by Geissler and Chandler who observed a similar Gaussian behavior but showed, by analysis of the nonequilibrium energy gap distributions after excitation, that it was nonstationary, i.e., the distribution changed with time.⁶ 

The previous results show how breakdowns of the Gaussian statistics approximation can be predicted without resorting to full nonequilibrium simulations. However, it is also important to understand the origin of the breakdown. For the dipole-flip model, it has been proposed that the breakdown of the Gaussian statistics approximation is caused by changes in the solute-solvent hydrogen-bonding (H-bonding) network upon excitation.^11,13 Here, we investigate this hypothesis in two ways. First, we consider an aprotic solvent that cannot hydrogen bond with the dipole-flip solute. Second, we examine how the Gaussian statistics approximation allows a partitioning between different contributions to the solvation energy and thus an isolation of the non-Gaussian motions responsible for the breakdown.

Choosing a coordinate to partition, as will be shown in Secs. V A and V B, for the dipole-flip model is fairly intuitive, and previous studies have made the identification simple. For other systems, choosing an adequate partition may not be as straightforward and may require further analysis of equilibrium simulations. These coordinates must be sufficiently small in number such that the central limit theorem does not apply but their contributions must be sufficiently large that their effect is important in the overall nonequilibrium response.^14,27,28

If H-bonding rearrangements are the origin of the failure of the Gaussian statistics approximation, there should be no such breakdown for solvation in an aprotic solvent. We have tested this idea by comparing the nonequilibrium TDSS result, S(t), for the dipole-flip model in acetonitrile to the Gaussian statistics result, C(t), as shown in Fig. 8. The data show that the Gaussian statistics approximation is in excellent agreement with the full, nonequilibrium MD result. This supports the implication of H-bonding dynamics in the non-Gaussian behavior.

This can be examined further by comparison of the higher-order correlation functions and their Wick’s factorizations. These are shown for the dipole-flip solute in acetonitrile in Fig. 9. The results clearly show that for this non-H-bonding solvent, the Gaussian statistics approximation is an excellent description of the nonequilibrium solvation dynamics, with only a small deviation from zero of the lowest odd-order correlation function, C₃(t), at short times. The data also reinforce the diagnostic utility of the comparison of the full and factorized higher-order correlation functions for evaluating the Gaussian statistics approximation from purely equilibrium simulations. Finally, they illustrate why the static linear-response approximation is not the appropriate interpretation of Eq. (12). Even in this case where $C(t)≃S(t)$⁠, each higher-order correlation function is larger than all the lower order ones and is thus not remotely negligible, while, in contrast, the factorization assumed in the Gaussian statistics approximation is highly accurate.

These results implicate the H-bonds between the solute and solvent molecules in the breakdown of the Gaussian statistics approximation in water and methanol. To explicitly demonstrate this we take a different tack. Specifically, we have partitioned the trajectories into a hydrogen-bonded partition—which is presumably non-Gaussian—in which only those solvent molecules that are H-bonded to the solute and the solute itself are included and a non-H-bonded partition—which is presumably Gaussian—in which only those solvent molecules that are not H-bonded to the solute and the solute itself are included. Solute-solvent H-bonds were identified using geometric criteria for both methanol and water. Specifically the criteria were based on the oxygen-site distance $rOα≤3.5$ Å, the hydrogen-site distance $rHα≤2.45$ Å, and the hydrogen-oxygen-site angle $∠HOα≤30°$⁠, where $α=A$ or B of the solute. Transiently broken H-bonds were also included in the H-bonding partition.

One can easily show that the Gaussian statistics approximation for the nonequilibrium response of a property A can be decomposed into various contributions to A by an additive partitioning. Specifically, we can write the fluorescence energy as $ΔE(t)=ΔEHB(t)+ΔENHB(t)$⁠, where $ΔEHB$ is the contribution to the energy gap due to solvent molecules H-bonded to the solute and $ΔENHB$ the contribution from all other, non-H-bonded, molecules. The Gaussian statistics approximation is distributive such that

where the Gaussian statistics approximation gives

and similarly $⟨ΔENHB(t)⟩ne−⟨ΔENHB(∞)⟩ne≃CNHB(t)$⁠; note that C_HB(t) and C_NHB(t), unlike C(t), are not normalized. The correlation functions C_HB(t) and C_NHB(t) can be further divided into components using $ΔE(0)=ΔEHB(0)+ΔENHB(0)$ in Eq. (15). However, as is clear from the above derivation, this is not necessary for the present purposes because C_HB(t) (C_NHB(t)) itself represents the contribution of the H-bonded (non-H-bonded) molecules to the unnormalized TDSS. The same diagnostics applied in Sec. IV can be independently applied to C_HB(t) and C_NHB(t) by comparison of higher-order correlation functions to the corresponding factorized forms. In this way, the Gaussian statistics approximation can be independently tested in different partitions of the relevant coordinates in the same way as was done for the total fluorescence energy.

Note that the energy gaps, $ΔEHB$ and $ΔENHB$⁠, involve only the Coulombic contributions to the energy difference of the solute-solvent interactions between the excited and ground states. The Lennard-Jones interactions are the same in the ground and excited states and thus cancel in the energy gap. Although there is a contribution to the energy gap from the long-range electrostatics due to the solute interaction with itself in different periodic images in the Ewald sum, we have not included it here.

The even-order time correlation functions for the H-bonded and non-H-bonded partitions are compared with their factorized forms in Figs. 10 and 11. In Sec. IV we found that for water, the factorized forms were in excellent agreement with the full correlation functions (see Figs. 5(a) and 6(a)). Thus it is not surprising that the factorization is accurate for both the fourth-order and sixth-order partitioned correlation functions shown in Figs. 10(a) and 11(a), respectively. On the other hand, the same is not true for the odd-order correlation function which is zero in the Gaussian statistics approximation but was shown to be non-zero at short times for the dipole-flip model in water (see Fig. 4(a)). The results for the corresponding partitioned correlation functions are presented in Fig. 12 and show that $C3,NHB(t)≃0$ for all times and the deviations from Gaussian statistics are isolated in the H-bonding contributions represented in C_3,HB(t). This indicates that the non-Gaussian character of the dipole-flip model in water is present primarily in an asymmetry of the energy gap distribution as shown in Fig. 7(a) and is associated solely with the waters H-bonded to the solute.

The results are different for the methanol solvent, which can be seen from the comparison of the fourth- and sixth-order correlation functions to their factorized forms in Figs. 10(b) and 11(b) for the H-bonding and the non-H-bonding partitions. As in the case of water, the non-H-bonding partitions obey Gaussian statistics based on the close agreement between the correlation functions and their factorizations. However, in methanol, the H-bonding partition shows a significant deviation from the Gaussian statistics approximation in both the fourth- and sixth-order correlation functions. A similar behavior is observed for the partitioned third-order correlation functions shown in Fig. 12(b), for which the non-H-bonding contribution is effectively zero at all times (and thus consistent with Gaussian statistics) while the H-bonding correlation function is non-zero (albeit with significant error bars). Thus, it is clear that the non-Gaussian behavior is associated with the H-bonding of solvent molecules to the dipole-flip solute.

This analysis represents a general approach to identifying the origins of non-Gaussian statistics from only equilibrium simulations. It is interesting to note that in the case of water, the non-Gaussian behavior only arises from the odd-order correlation functions that are related to the asymmetry of the energy gap distribution. In contrast, for methanol, the non-Gaussian behavior is observed at all orders (even and odd), indicative of the stronger failure of the Gaussian statistics approximation that is also reflected in the energy gap distribution. This is presumably related to the smaller number of methanol molecules that contribute to the dipole-flip solute solvation through H-bonding, compared to the case of water, that is, the central limit theorem does not apply for the H-bonding solvent molecules in methanol but nearly does in water. It is important to note that the partitioning used here was relatively easy to find and motivated by the previous work on the dipole-flip model.^6,13 Indeed, we cannot be sure that the partitioning has completely identified the origin of the non-Gaussian statistics in that the H-bonding partition we have defined may still contain degrees-of-freedom that are Gaussian. The approach outlined here does not provide a prescription for how to determine the non-Gaussian coordinates but rather gives a straightforward method for evaluating hypotheses for such coordinates.

A key aim of investigating the origins of failures of the Gaussian statistics approximation is to ultimately develop approaches for correcting for non-Gaussian effects and thereby obviate the need for performing nonequilibrium MD simulations. A possible approach for this is offered by the deviations from Gaussian statistics manifested in the difference between the higher-order time correlation functions and their factorized forms.

As shown in Sec. II B, the Gaussian approximation is obtained by replacing the higher-order correlation functions with their factorized forms. Corrections to this approximation are most straightforwardly examined by considering the even- and odd-order correlation functions separately. The latter have no contribution in the Gaussian approximation and thus the corrections for odd-order terms can be included simply by addition of any non-zero higher-order correlation functions. Thus, the lowest-order correction to the nonequilibrium energy gap would be

where the subscript “3,corr” indicates that this correlation function represents a correction up to the third-order. This corrected result has been applied to both the water and methanol solvent cases and the results are shown (dashed violet lines) in Fig. 13, plotted in normalized form. In the case of water, the effect of the correction is to predict a significantly more rapid decay in the energy gap compared to the Gaussian statistics approximation alone followed by a nearly constant but non-zero value. Despite the fact that the correction is in better quantitative agreement with the full, nonequilibrium TDSS at most times compared to the Gaussian statistics result, it more poorly predicts the shape of S(t). In contrast, for the methanol solvent, the C_3,corr(t) corrected result is in significantly better qualitative and quantitative agreement with the nonequilibrium TDSS compared to the Gaussian statistics approximation. In particular, it exhibits a larger amplitude for the initial, inertial component as well as a faster decay at longer times. We note that the higher-order correlation functions are obtained with essentially no additional effort above that required for the Gaussian statistics approximation. However, a key issue with this correction is the large statistical errors associated with the C₃(t) correlation function it is based on as illustrated in Fig. 4.

It is less clear how to correct for deviations from Gaussian statistics in the even-order correlation functions. One approach is to note that the Gaussian approximation assumes

Then an approximation to the Gaussian statistics result for a given n can be obtained as

Clearly, for n = 0 this gives the Gaussian statistics approximation, C(t), but for $n>0$ it yields a correction to this C(t) based on the time-dependence of the higher-order correlation function. Then a reasonable scheme for correcting the Gaussian statistics approximation is to average the estimates of this correlation function from the different orders, e.g.,

This result is shown for methanol solvent in Fig. 13(b) and is in agreement with the full TDSS result for all times (within error bars). It represents a slight improvement over the C_3,corr(t) correction, particularly for times less than $∼3$ ps. However, it is clear that the corrections for the even-order terms are less significant than the first odd-order correction represented in C_3,corr(t). For water, the higher-order, even correlation functions are in good agreement with the Gaussian approximations so that C_3–6,corr(t) = C_3,corr(t) within our calculations.

Finally, we note that an alternative approach to correct the Gaussian statistics approximation is to add the cumulants, i.e., the differences between the higher-order correlation functions and their factorized forms, to the Gaussian statistics approximation result.²⁹ If, however, the resulting series of higher-order terms does not converge (as has been the case in previous studies^8,19), then this construction will fail at least theoretically if not also practically. Indeed, this is what we have observed for both the water and methanol solvents for the dipole-flip solute.

The dynamic linear response and Gaussian statistics approximations are used frequently, particularly in applications to solvation dynamics. Because of the general success of these approaches, the cases where they fail have drawn significant theoretical and experimental interest.

The dipole-flip model, in particular, has been studied by others following the early work by Ladanyi and co-workers establishing it as a key example of nonlinear response. Geissler and Chandler computed the nonequilibrium distribution of the excited- and ground-state energy gaps for the dipole-flip model in water and found that the statistics of the energy gap are Gaussian but are nonstationary.⁶ These nonstationary Gaussian statistics were attributed to the change in the H-bond network due to the sensitivity of dielectric response on boundary conditions. As the solvent molecules rearrange themselves upon solute excitation, boundary conditions of the solvent region change, influencing the solvent response. These results also indicate that the non-Gaussian statistics result from some alternative, non-Gaussian degree(s) of freedom, which cause the overall energy gap to be non-Gaussian and result in nonlinear response. The present work is consistent with these results, particularly in explicitly isolating the role of the H-bonding network, through the projection of the energy gap fluctuations (Sec. V), as the few degrees of freedom responsible for the non-Gaussian behavior.

Breakdowns in Gaussian statistics for solvation dynamics have also been reported in a number of other systems. Schwartz and co-workers used MD to study a single atom capable of carrying a charge dissolved in water.^27,28 They found that solute size changes (e.g., induced by ionization) cause a breakdown in linear response. In particular, both the dynamic linear-response and Gaussian statistics approximations give incorrect results in the inertial regime of relaxation. This is associated with the size change of the solute, which results in a restructuring of the H-bond network due to inertial motion of proximal OH groups.

Systems of the solvated electron serve as another class of non-Gaussian solvation dynamics. Rossky and co-workers have shown that the Gaussian statistics approximation fails to describe the solvation dynamics of an excess electron in methanol.^15,30 They conducted mixed quantum-classical MD (MQC-MD) simulations and computed the normalized TDSS from nonequilibrium simulations as well as linear response correlation functions from equilibrium simulations in two states: the initial state (the neat solvent) and the final state (the solvated electron system). Both linear response correlation functions differed from the TDSS, which was, like the previous example, attributed to large structural changes in the H-bond network.¹⁵ Upon excitation, mechanical forces from localization of the electron push H-bonds away, so at short times (<500 fs) the number of H-bonds around the electron decreases.

Bedard-Hearn et al. used MQC-MD to study the hydrated electron to show how their previous findings on Gaussian statistics breakdowns from solute size changes apply to the system.³¹ Using Steele theory,³² they found that the Gaussian statistics breakdown for size changes in the hydrated electron from excitation occurs from significant deviations in the solvent librations that interfere with the electron p-like wavefunction. These studies on solute size changes are possible explanations for the breakdown of Marcus theory (i.e., linear response) for reactions of cation transition metals with the hydrated electron.³³ Specifically, Kanjana et al. studied the kinetics of hydrated electron reactions with a number of dication transition metals using UV-Vis absorption spectroscopy, which showed deviations from Marcus theory.³³ That Marcus theory fails for these reactions is perplexing; the reactions are Arrhenius, are not diffusion limited, and are thermally activated, all of which point to the validity of Marcus theory for this system.³⁴ 

Laird and Thompson studied TDF and dye molecule diffusion in nanoconfined solvents.⁸ The results showed that the Gaussian statistics approximation agrees reasonably with nonequilibrium MD results for the TDSS of a nanoconfined dye molecule (whereas dynamic linear response does not) but more poorly for the dye molecule position. Additionally, in the case of nanoconfined solvation dynamics, the assumption invoked by static linear response (neglect of higher-order terms) is not justified. In fact, the higher-order correlation functions not only do not decrease but also increase by orders of magnitude with each term in the series. The energy gap was found to obey Gaussian statistics for a fixed dye molecule position but with position-dependent properties leading to the non-Gaussian behavior. This is consistent with the nonstationary Gaussian statistics observed by Geissler and Chandler.⁶ 

Li and co-workers have studied the validity of linear response for TDF of Staphylococcus Nuclease (SNase).^19,35 They performed equilibrium and nonequilibrium simulations and found that the Gaussian statistics approximation proved more favorable than dynamic linear response based on comparison with nonequilibrium simulations.³⁵ Furthermore, they showed that the static linear-response approximation was not valid for the system they studied by computing higher-order correlation functions, which also increased with increasing order.¹⁹ 

These examples should not be used to infer that breakdowns of the Gaussian statistics approximation are limited to systems where H-bonding is present. There are, indeed, many systems displaying nonlinear response in aprotic solvents. One example is the solvation of neutral sodium atoms in tetrahydrofuran (THF). Schwartz and co-workers have investigated the solvation of a neutral Na atom in liquid THF using two approaches: one in which a sodium atom was produced by removing an electron from Na⁻ and the other in which an electron was added to Na⁺.³⁶ Were the system to obey Gaussian statistics, the dynamics of both experiments should exhibit the same time scales because they occur on different sides of the same harmonic free energy surface. Bragg et al. however, observed that the solvation dynamics of the sodium atom prepared from Na⁻ were much slower than that of the sodium atom prepared from Na⁺.³⁷ Bedard-Hearn et al. used MD to study sodium deionizing in THF; however, instead of an explicit breakdown in linear response resulting from comparison of C(t) and S(t), they found that the breakdowns were hidden in the mechanism of relaxation for each property, that is, S(t) and C(t) were in good agreement but the components for each differed.²⁶ Using Steele theory, they analyzed the contributions to C(t) and S(t) from the Coulomb and Lennard-Jones partitions of the energy and found that the relaxation from each contribution differed between the nonequilibrium value and linear-response approximation.

Another example is the relaxation of rotational energy. Moskun et al. studied the energy dissipation in highly rotationally excited CN dissolved in liquid argon.³⁸ The rotational excitation in the gas phase would be observed for an extended period compared to the liquid phase in which, as linear response would predict, the rotational energy would rapidly dissipate. They observed, however, that the dynamics went against the intuition of both phases: initially a linear-response-like dissipation was observed followed by a sudden switch to gas phase-like dynamics. This switch corresponds to a molecular event that cannot be predicted by linear response theory. In fact, the rotations of CN produce a distinct cavity in the liquid at short times, which extends the coherence of the rotations. This was confirmed by Tao and Stratt in which they derived a linear-response approximation to the normalized response function and computed both linear-response correlation functions and nonequilibrium MD simulations.⁷ Similar findings were reported for simulations of OH rotors.³⁹ An analogous result was observed in a combined experimental and computational study of vibrational relaxation of HCN formed as the product of a chemical reaction.⁴⁰ 

In each of the aforementioned examples of a breakdown in linear response, i.e., non-Gaussian statistics, specific molecular rearrangements of the solvent are reported that are strongly coupled to the solvation dynamics. It is from these molecular events—cavitation of a solvated rotor, H-bond restructuring, solute cavity expansion or contraction, etc.—that a breakdown of the Gaussian statistics approximation can be attributed. A key commonality is that each involves only a small number of degrees-of-freedom preventing Gaussian behavior via the central limit theorem,¹ as has been previously suggested by Schwartz and co-workers.^27,28 This is supported in the case of the dipole-flip model by the more nonlinear response observed in methanol compared to water (Fig. 3) and observed by Fonseca and Ladanyi for a smaller solute compared to a larger one in methanol.¹² That is, the response is less linear when fewer solvent molecules are involved in the molecular rearrangements, which for the dipole-flip model are associated with H-bonding to the solute and influenced by both the solvent and solute sizes. Similarly, the nonlinear response observed for the solvated electron^15,30 and solute size changes^27,28 is associated with rearrangements of the liquid structure in the first solvation shell, i.e., involving a small number of molecules.

Importantly, Tao and Stratt noted that changes in liquid structure were a common feature in many of the studies cited above and further suggested “…that the relative slowness of structural relaxation could be an important prerequisite for linear-response breakdowns.”⁷ The observation of non-Gaussian structural rearrangements occurring on slower time scales has been observed in the case of rotational relaxation, where it is associated with the formation of a solute cavity^7,38,39 as well as for vibrational relaxation which exhibits an analogous behavior.⁴⁰ In the case of solvation dynamics, a similar picture was posited by Turi et al. for the solvated electron in methanol, again associated with a change in solvent structure around the electron.¹⁵ For the dipole-flip model, the picture is more complicated. The linear response breakdown tends to appear at longer times, after the inertial dynamics, as was pointed out by Fonseca and Ladanyi,¹² who noted “…both fast and slow solvation dynamics have substantial contributions from structural rearrangements within the first solvation shell. After the perturbation is applied, the destruction of an existing local structure is rapid, but a build-up of a new one is much slower.” The present analysis shows deviations from Gaussian statistics at all time scales, which can be primarily observed in the odd-order correlation functions for both water and methanol (Fig. 4) but also in the even-order correlation functions for methanol (Figs. 5 and 6). However, it is also the case that the H-bonding rearrangements responsible for the nonlinear response in this model do occur on longer time scales. Thus, while the current results appear to suggest that structural relaxation need not be slow on the time scale of the response to lead to nonlinear response, this topic deserves additional examination for other systems and from other perspectives.

The static linear response and Gaussian statistics approximations have been outlined and investigated in the context of the dipole-flip model. Where the static linear response approximation assumes all higher-order correlation functions are negligible, the Gaussian statistics approximation factorizes them into lower-order correlation functions based on Wick’s theorem and, upon summation, gives the same mathematical expression as static linear response.

For a dipolar solute immersed in water or methanol, the approximation fails, confirming prior knowledge on the subject,^6,12,13 but it has also been shown that the physical approximation made in static linear response is unjustified: the higher-order correlation functions are not negligible and, in fact, grow by orders of magnitude. Thus, what is often referred to as the static linear response approximation is really the assumption of Gaussian statistics.

The Gaussian statistics approximation assumes that the relevant dynamical variable, e.g., the energy gap in solvation dynamics, is a Gaussian random variable. An accepted approach for evaluating this approximation is to compare it with full nonequilibrium MD simulations. This is rarely done, however, because the nonequilibrium simulations are significantly more computationally expensive than the equilibrium MD simulations required for the Gaussian statistics result. Thus, we have addressed the question can non-Gaussian statistics be predicted by relying only on equilibrium simulations? The answer is yes. Specifically, Gaussian statistics assumes that higher-order time correlation functions can be factorized into two-point time correlation functions; both the higher-order and factorized time correlation functions can be evaluated from the same equilibrium MD simulation. The dipole-flip model in water and methanol exhibits a breakdown of the Gaussian statistics approximation based on comparison with nonequilibrium MD results. This breakdown is also reflected in differences between the equilibrium higher-order and factorized time correlation functions, which thus offers a method for the evaluation of the Gaussian statistics approximation. This complements a direct examination of the distribution of the energy gap, which for simulations in methanol is clearly non-Gaussian, but for water is nearly Gaussian.

When there is a breakdown of the Gaussian statistics approximation it is natural to ask what is the origin of the non-Gaussian statistics? It has been proposed that the nonlinear response for the dipole-flip model in H-bonding solvents can be attributed to the solvent-solute H-bonds^6,12,13 which must rearrange upon the electronic excitation that changes the identity of the dipole-flip solute atom that can act as an H-bond acceptor. Here, we have shown that the Gaussian statistics approximation agrees with nonequilibrium dynamics for a non-H-bonding solvent, acetonitrile, indicating that the lack of H-bonding recovers agreement between the Gaussian statistics approximation and nonequilibrium dynamics.

An additional, more direct, test to the non-Gaussian nature of the solute-solvent hydrogen bonds was performed in which the equilibrium trajectories were partitioned into a solvent-solute H-bonding component and a bath (non-H-bonding) component. The same comparisons of higher-order correlation functions and their factorizations that, as discussed above, can be used to identify non-Gaussian statistics can then be applied to each partition separately. This analysis clearly shows that solvent-solute H-bonds are responsible for the breakdowns of the Gaussian statistics approximation for methanol and water. Moreover, it shows how the origins of the non-Gaussian behavior can be identified and isolated in general through such partitioning.

Ideally, one would like to find a way to improve upon the Gaussian statistics approximation when necessary. Thus, we asked can the Gaussian statistics approximation be corrected to describe nonlinear response? This is a more challenging problem than identifying the non-Gaussian behavior and its origins. Nevertheless, the results indicate that there is some promise in using the information about the non-Gaussian character of the dynamics embodied in the higher-order time correlation functions to improve upon the Gaussian approximation. The results for methanol, where the non-Gaussian behavior is more prominent, are encouraging, while those for water are less so. This is an area that is ripe for future exploration.

We have placed the present study of the dipole-flip model in various solvents in the context of other examples of breakdowns of the Gaussian statistics approximation throughout the literature. In each case discussed, one can reasonably attribute the origin of this breakdown to some specific solvent mechanism, involving a small number of relevant coordinates, that occurs in the nonequilibrium dynamics that does not occur at equilibrium. Though there have been no tests of these systems based on the assumptions invoked in the static linear response and Gaussian statistics approximations such as that presented here, they represent both model and real, experimentally accessible systems to which the present analysis could be applied to advantage.

Support for this work was provided by the Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy (Grant No. DE-FG02-05ER15708). The calculations presented were performed at the University of Kansas Advanced Compute Facility. W.H.T. thanks Professor Branka M. Ladanyi for helpful comments in the early stages of this work. The authors thank Professor Brian B. Laird for many useful discussions.