Simple Quantum Dynamics with Thermalization.

In this paper, we introduce two simple quantum dynamics methods. One is based on the popular surface-hopping method, and the other is based on rescaling of the propagation on the bath ground-state potential surface. The first method is special, as it avoids specific feedback from the simulated quantum system to the bath and can be applied for precalculated classical trajectories. It is based on the equipartition theorem to determine if hops between different potential energy surfaces are allowed. By comparing with the formally exact Hierarchical Equations Of Motion approach for four model systems we find that the method generally approximates the quantum dynamics toward thermal equilibrium very well. The second method is based on rescaling of the nonadiabatic coupling and also neglect the effect of the state of the quantum system on the bath. By the nature of the approximations, they cannot reproduce the effect of bath relaxation following excitation. However, the methods are both computationally more tractable than the conventional fewest switches surface hopping, and we foresee that the methods will be powerful for simulations of quantum dynamics in systems with complex bath dynamics, where the system-bath coupling is not too strong compared to the thermal energy.

Introduction

Dissipation in quantum systems is a phenomenon playing a role in many physics and chemistry problems. An accurate description of quantum dynamics is thus crucial for understanding many physical and chemical phenomena. Such phenomena include photosynthesis, vibrational energy redistribution, and photoisomerizations. A popular approach to this problem is rooted in coupling a system described by a quantum mechanical approach with a bath described at a more approximate level of the theory.[1,2] Since Ehrenfest proposed the first such approach[3] numerous alternatives have been proposed. The Hierarchical Equations Of Motion (HEOM) theory[4−6] and the QUAsiadiabatic Propagator Path Integral (QUAPI) method[7] stand out as two formally exact methods. These methods, while accurate, are computationally prohibitively demanding as the size of the quantum system and bath increase. A popular alternative is the Surface Hopping family of approximations[8] with the so-called Fewest Switches Surface Hopping (FSSH) as proposed by Tully[9,10] as the most common approximation. Numerous other approximate methods for quantum dynamics in open systems have been developed and tested.[1,11,12,12−37] The key motivation for this paper is to develop simple approximate quantum dynamics algorithms that allow treating sizable and complex systems, for which the formally exact methods are computationally too expensive or too complicated to formulate.

In the following, a brief overview of the most popular of the currently existing methods for calculating quantum dynamics is given. Currently, one of the most popular quantum dynamics approaches is the HEOM. It was developed for different types of bath coordinates including overdamped[5,6] and underdamped[37] Brownian oscillator modes. Combinations of essentially arbitrary Gaussian baths can be modeled. However, for every added bath mode the computational cost of the HEOM increases. Significant work has been made to make this method computationally more tractable including approximations made in the hierarchical expansion used to treat the bath[38] and development of efficient GPU codes.[39] The key limitations of this method are, thus, that it scales very unfavorably with systems size and that it is limited to Gaussian bath dynamics. The bath is treated in a stochastic manner, and the bath density of states must be preparametrized.

For weakly coupled systems perturbative approaches for quantum dynamics as that developed by Redfield[40] and modifications of this[41] have become popular. A perturbative approach developed by Sumi is based on a generalization of Förster energy transfer to multichromophoric systems.[42−45] These approaches are computationally rather cheap but limited by their perturbative nature and not applicable to degenerate systems, where the energy gap between the donor and acceptor become small compared to the coupling in which the perturbation is made. These methods further build on the assumption of Gaussian dynamics and a stochastic description of the bath.

The Surface Hopping approach proposed by Tully[9,10] is simply a postulated solution to solve quantum dynamics even recent attempts have been at least partially successful in a formal derivation.[46] In essence, the system–bath interaction is treated in such a way that the bath at all times only feels a force from a specific state of the system, and through carefully chosen stochastic “hops” between these states trajectories are obtained, where correlations between bath and system are preserved. From our perspective, one of the powerful points of Surface Hopping is that it allows the calculation of quantum dynamics using bath trajectories. It is therefore not limited to Gaussian dynamics, and it can describe systems with very anharmonic bath modes and systems, where the notion of a bath density of states is not valid. This also allows a quite simple analysis of the correlation between bath and system dynamics.[47] Still, because of the nature of the method, quite a few fundamental questions on the application and details of how to apply the method exist, and many adapted versions of the method have emerged.[8,48−54] The method has still proven useful for studying many chemical problems, and it has, for example, been implemented in the Newton-X software package.[55] A particular drawback of the method is that it requires the simultaneous solution of the explicit quantum and classical dynamics, which for large systems can be rather time-consuming.

Alternatively, the quantum dynamics can be propagated on the bath potential energy surface corresponding to the ground state,[56] equivalent to Ehrenfest dynamics without quantum feedback.[3] This type of approach has recently been implemented to model linear absorption and two-dimensional spectroscopy with quite some success[57−62] and will here be denoted the Numerical Integration of the Schrödinger Equation (NISE) approach. An important advantage of the NISE method is that it can be applied to systems with arbitrary bath dynamics, including baths with multiple underdamped and non-Gaussian[63−66] modes without any additional computational cost of the quantum propagation. As long as a trajectory of a time-dependent Hamiltonian can be constructed the method can be applied. Furthermore, non-Condon effects[67] including effects of simple rotations can be trivially included in spectral calculations.[68] The method scales very favorably with system size,[60] and it is applicable to simulate spectroscopic observables of systems consisting of hundreds of coupled chromophores[69−73] and treat more than 1 × 105 coupled quantum states needed to simulate two-dimensional infrared and sum-frequency generation spectra of proteins.[60,74] The obvious drawbacks are that as in this method the quantum system feels the classical bath the propagation of the classical bath does not depend on the state of the quantum system. The quantum dynamics thermalize to an infinite temperature and bath relaxation, for example, causing spectral Stokes shifts are absent.[75]

The goal of the present paper is to develop computationally tractable methods for approximate quantum propagation at finite temperatures. The methods should allow the use of arbitrary bath dynamics, be applicable to large systems, and have approximately correct thermalization. For the first approach, the basic idea that will be elaborated on in the Methods Section is that from the equipartition theorem we know how much thermal energy is available on average in each bath coordinate. We will use this knowledge to determine when a hop is allowed in a Surface Hopping algorithm. We will neglect any feedback of the quantum system on the classical bath allowing the use of precalculated classical trajectories. A key assumption for this approximation to apply is that the dissipation of energy within the classical system is faster than the frequency of attempted surface hops in the Surface Hopping algorithm. This approximation will potentially solve the thermalization problem but obviously not the bath relaxation problem, for example, leading to Stokes shifts. We will denote the new Surface Hopping approach the Simple Equipartition Surface Hopping (SESH) approach. While the number of Surface Hopping approximations is very large[8] to the best of our knowledge this approximation has not been presented before. We further define an adapted NISE approach using simple rescaling of the nonadiabatic couplings. This approach will be denoted the thermalized numerical integration of the Schrödinger equation approach (TNISE). We foresee that such methods may be particularly useful for calculating two-dimensional infrared[76] (2DIR) and two-dimensional UV/visible[77] (2DUVvis) spectra for large systems, which explicitly depend and report on the quantum dynamics of the excited states. Here we will test the new models against the exact HEOM method for coupled two-level systems, which has become a standard test for approximate quantum propagation methods.[5,12,15,21,54]

The outline of the remainder of this paper is as follows. In the Methods Section, we will briefly discuss the basic coupled two-level models used for testing the new approximations, the HEOM and NISE implementations used for comparison will be summarized, and the new SESH and TNISE methods will be described in detail. In the Results Section, we will compare results for these methods for four different situations mimicking typical problems, where the new methods may be applicable. These are two problems for vibrational dynamics and two problems for electronic excitations. Finally, the conclusions will be presented along with a discussion of the potential applications of the new methods.

Methods

We will consider a simple system of coupled two-level systems described by the HamiltonianHere, the wave functions for the two sites are |1⟩ and |2⟩, while the difference in the energy between the two sites is ΔE. As the dynamics of the quantum system does not depend on the absolute energy of the sites the average energy of site 1 was set to 0. The site energies fluctuate depending on the bath degrees of freedom, x1, and x2, with a coupling constant, σ. The two sites are coupled with coupling constant J, which for simplicity is assumed constant. The model and the involved parameters are also illustrated in Figure .

Figure 1

Dimer of two-level systems is illustrated with the energy gap ΔE between the excitations, the coupling between the two sites J, and the bath-induced energy fluctuation magnitude σ.

The bath degrees of freedom were treated as an overdamped Brownian oscillator characterized by a Gaussian distribution of the bath coordinates x1 and x2. The connected frequency time-correlation function is C(t) = ⟨σx(t)σx(0)⟩= σ2 exp(−t/τ)δ = 2kBTλ exp(−tγ)δ. Here i and j enumerate the sites. The memory can be characterized by either of the complementary parameters τ or γ, and the fluctuation magnitude can be characterized by either σ or temperature and reorganization energy λ. The correlation function is illustrated in Figure . For the surface hopping and NISE simulation, trajectories of the fluctuating bath coordinates were constructed using a simple stochastic approach outlined in ref (68) ensuring the correct overdamped Brownian oscillator behavior. The Brownian oscillator dynamics is formally a solution to the equation of motion:where the magnitude of the fluctuations of the white-noise random force term f(t) determines the magnitude of the fluctuations in the coordinate x.[78,79] While we will use this model in the present paper to allow comparison with the HEOM method, it is important to realize that presented new methods do not depend on how the fluctuating frequencies are generated. The bath dynamics can thus be represented by virtually any stochastic or statistical model or may arise from frequency models[80−86] connected with molecular dynamics simulations.

For the HEOM implementation, we followed ref (5) and validated the implementation by comparing with the data presented in that paper. That implementation neglected fast fluctuation effects resulting in contributions from Matsubara modes at low temperature and is only valid when ℏ/τ< kBT. The implementation of the NISE procedure followed ref (17), while the Surface Hopping implementation followed ref (87) with the specific changes outlined below.

In the SESH method, a primary (ϕP(t)) and an auxiliary (ϕA(t)) wave function are used. Initially, the two are identical, and they are chosen to be identical to one of the instantaneous eigenstates of the Hamiltonian. The time-dependent Schrödinger equation determines the dynamics of the primary wave function. Therefore, it fulfills the equationTo solve this, the primary wave function is propagated in the diabatic site basis avoiding diverging nonadiabatic couplings, which may arise if propagating in the instantaneous adiabatic basis. For this propagation short time steps Δt are used assuming that the Hamiltonian can then be considered constant resulting in the wave function update:Here defining the time-evolution operator U(t + Δt,t). The auxiliary wave function is propagated by first assuming an adiabatic change. As a first step, the auxiliary wave function at time t + Δt is chosen as the instantaneous eigenstate ϕE(t + Δt) of the Hamiltonian at time t + Δt maximizing the overlap with the auxiliary wave function at time t. Then the possibility of a surface hop is considered following the standard FSSH protocol. A hopping probability for transferring to another instantaneous eigenstate ϕE(t + Δt) is given byIf the probability is negative it is nullified. Then a random number η between zero and one is generated, and a hop is executed to the state j ifThis is all according to the FSSH protocol. One of the key differences with the new approximation arises when a hop is accepted according to this criterion. In the FSSH protocol, a check will be performed if sufficient kinetic energy is available in the bath mode along the nonadiabatic coupling vector, and in the case of a jump the energy difference connected with the jump will be exchanged with the bath. In this approach, we deliberately want to use a pregenerated bath trajectory. Alternatively, we use the assumption that the kinetic energy available in any given bath mode is given by the equipartition theorem. We further assume that the energy dissipation in the bath is faster than the time between hops in the system. The kinetic energy available for hops can, thus, be considered to be independent of any previous successful hops. A second random number, ξ, is, thus, generated by taking the absolute value of a Gaussian distribution with standard deviation of kBT mimicking a random kinetic energy. If the kinetic energy required for moving from the instantaneous eigenstate corresponding to the present auxiliary state to the state ϕE(t + Δt) is smaller than ξ the hop is accepted, and the auxiliary wave function for the time t + Δt is changed to ϕE(t + Δt). In contrast to the FSSH algorithm nothing is done to the bath trajectory, neither in the case of an accepted nor in the case of a refused hop.[49] This is obviously an approximation that will fail in the cases of numerous successive upward hops that all require energy from the same bath mode. The other part of the present approximation is the neglect of the force exerted by the quantum system on the classical bath coordinates x at every time step as given by F(t) = ⟨ϕA(t)|dH/dx|ϕA(t)⟩. The consequence of neglecting this force is the loss of the inclusion of a Stokes shift induced by the bath. In contrast to the FSSH surface hopping the SESH algorithm outlined above does not ensure conservation of energy of the whole system.

We further develop an approximate thermalizing version of the NISE approach. The starting point for this is the equation for propagating the wave function in the adiabatic basis:[58]where c is the wave function coefficient for each adiabatic basis function (instantaneous eigenstate), ϵ is the eigenvalue corresponding to that adiabatic basis function, and S is the nonadiabatic coupling, which is defined S ≡ ⟨ϕ̃(t)|ϕ̃̇(t)⟩. In this formulation of the nonadiabatic dynamics on the ground-state potential energy surface, the meaning of the two terms is quite simple. The first term accounts for the phase acquired with time for each individual instantaneous eigenstate, while the nonadiabatic coupling accounts for population transfer between these states. We know that the predicted equilibrium populations will be equal according to the high-temperature limit in this case, and to correct for this incorrect thermalization we redefine the nonadiabatic coupling. For simplicity, we define the propagation in matrix form:Here ψ̃(t) is the wave function in the adiabatic basis, the matrix Ũ(t + Δt,t) is the adiabatic propagation, and S(t + Δt,t) is the nonadiabatic propagation, which can be formulated in terms of the eigenvector coefficient matrices C at times t + Δt and t. We now redefine the nonadiabatic matrix with a temperature correction:In this way transfer requiring energy from the bath is slowed, and transfer dumping energy to the bath is enhanced. The motivation of the scaling is that downward transfer should be faster than upward transfer with a ratio determined by the Boltzmann factor to reach equilibrium. The scaling factor is motivated by making the scaling factor mathematically equivalent for upward and downward transfer. The factor four in the exponent arise as the scaling is for the wave function propagation and as the populations are proportional to the square of the wave function the scaling of the population transfer will include a factor 2. This is obviously an ad hoc choice of scaling similar in spirit to previous methods.[15]N is chosen to normalize the problem such that N2 + ∑|S(t + Δt,t)|2 = 1. Care must be made that if the order of the eigenstates is swapped between two successive steps the new eigenstates are matched to maximize the overlap. Despite the normalization performed by rescaling the diagonal elements of the nonadiabatic coupling matrix the matrix is not unitary anymore, and the wave functions are renormalized after every propagation time step. When the temperature is infinite the correction term will become unity, and the normal NISE approach will be recovered.

In summary, we expect that the new approximations may work well when the reorganization energy is small compared to the temperature (λ ≪ kBT), where it may approximate the thermal relaxation well. In this case, the Stokes shift caused by the coupling with the bath is significantly low compared to the witdh of the frequency distributions of the site energies and the effect on the population transfer, thus, negligible. In contrast the simpler NISE method neglects both the thermalization and the Stokes shift. Still, this method has proven quite successful for predicting vibrational energy redistribution, for example, in the OH-stretch manifold of water.[88,89]

Results

The new method was tested for different parameters typical for realistic systems presented in Table . For each system, we calculated the time-dependent probability that a population created in the highest excited state was still in that state at a later time. For the HEOM the system Hamiltonian was used to determine the eigenstates used for this analysis, while in the SESH and NISE methods the instantaneous eigenstates were used for this analysis. For the HEOM and NISE further, the probability of starting at the site with the highest site energy and still being there at a later time was calculated as well. As the SESH formulated here only applies to the transfer between instantaneous eigenstates the site transfer could not be obtained. For the present simulations, the time step was set to 10 fs both when solving the HEOM, NISE, and SESH. For TNISE a smaller time step of 2 fs was employed to guarantee the correct identification of the eigenstates in successive time steps. For the hierarchy, a depth of five was found to be sufficient. For the NISE and SESH methods, 10 000 disorder realizations of the stochastic bath were used. All simulations were performed at T = 300 K unless explicitly stated.

Table 1

Parameters Selecteda for Two-Level Pair Systems Designed to Mimic Realistic Systems

system	J (cm–1)	ΔE (cm–1)	σ (cm–1)	τ (fs)
water in MeCN[68]	–43	0	60	50
amide I/amide II[92]	36	70	25	50
FMO[105]	–106	140	150	140
LH2[123]	47	300	100	100

References are given to papers on which these parameters are based. The temperature was set to T = 300 K corresponding to a thermal energy of 208.5 cm–1.

We first tested the SESH for parameters similar to those found for the OH-stretch vibrations of water in acetonitrile.[68] In this system, the two local OH-stretch vibrations are on average identical, and the coupling leads to a low-frequency symmetric mode and a high-frequency asymmetric mode (Δ < |J| < σ < kBT). The population dynamics of an initially excited asymmetric mode is followed in Figure a. On the one hand, for the NISE, TNISE, and SESH methods, the population of the instantaneous eigenstate is provided, while for the HEOM the population of the average eigenstate is provided, as this method is intrinsically treating the whole ensemble and does not allow looking at the instantaneous eigenstates. The SESH on the other hand at least as implemented here only allows excitation of the instantaneous eigenstates. Despite this difference, the TNISE, SESH, and HEOM populations are very similar, and they decay to the same final thermalized population, while the NISE method exhibits slower dynamics and decays to the high-temperature approximation 50% population. At short times the TNISE and SESH decay a bit faster than the HEOM. This is here attributed to the difference between the instantaneous and average eigenstates. In the site representation presented for the NISE, TNISE, and HEOM in Figure b the quantum dynamics is essentially identical. The small differences observed here between the NISE and HEOM methods are well in line with the fact that FTIR and 2DIR spectra for this system match the experimental observations for this system when modeled with the NISE approach.[63] This is not so surprising, considering that both the reorganization energy and the energy gap between the eigenstates are smaller than the thermal energy.

Figure 2

(a) The population of the highest (instantaneous) eigenstate following excitation of that state for the waterlike parameters according to the HEOM (red), NISE (blue), SESH (green), and TNISE (dashed black) approaches. (b) The population of one of the two sites following initial excitation of that site according to the HEOM (red), NISE (blue), and TNISE (dashed black) approaches.

We proceeded to model a system with an amide I mode coupled with an amide II mode. This mimics the most important relaxation channel from the amide I band of proteins[76,90] and, in particular, of the N-methylacetamide molecule[35,76,91−98] used for parametrizing and benchmarking protein models. The two amide modes are separated by ∼70 cm–1, which is still smaller than the thermal energy at 300 K (σ < |J| < ΔE < kBT). The population of the amide I eigenstate following initial excitation of this mode is presented in Figure a, and the results from the three methods are quite similar. Again the NISE population approach is at 50% occupation at long times as expected, while the HEOM goes to 40% reflecting the higher energy of the amide I mode. The initial decay of the SESH is faster than for the other methods, as was the case for the water example. The population in the long time limit is slightly higher than in the HEOM approach but considerably better than the NISE prediction. For the TNISE approach, the transfer is essentially identical to that observed in the HEOM approach. The population dynamics following the excitation of the isolated amide I vibration is shown in Figure b. Coherent vibrations between the two amide modes in the “site” basis are observed for both HEOM, TNISE, and NISE. In this case, the site basis corresponds essentially to the CO-stretch for amide I and CN-stretch for amide II.[91] In contrast to water, the predictions from NISE and HEOM start to deviate after ∼500 fs, where the HEOM begin to decay toward the correct thermal equilibrium. The TNISE, in this case, follows the HEOM. Again the relatively small deviation between the different methods supports the success of the NISE method in predicting FTIR and 2DIR spectra of coupled amide I-amide II systems;[99] however, it is also clear that the new SESH and TNISE methods would provide results in better agreement with the exact HEOM approach.

Figure 3

(a) The population of the amide I (instantaneous) eigenstate following excitation of that state for the amide-like parameters according to the HEOM (red), NISE (blue), SESH (green), and TNISE (dashed black) approaches. (b) The population of the isolated amide I vibration following initial excitation of that site according to the HEOM (red), NISE (blue), and TNISE (dashed black) approaches.

The energy differences between electronic states are larger than those of vibrational states, and more significant deviations are expected in such cases. The next system to examine is, thus, a dimer that mimics the states in the much-studied Fenna–Matthews–Olson (FMO) complex.[82,83,100−110] In this case, we used a dimer with an energy gap of 140 cm–1 and a frequency fluctuation magnitude of 150 cm–1. As we further consider a coupling between the sites of −106 cm–1 the energy gap between the average eigenstates is 254 cm–1 and, thus, larger than the thermal energy (|J| < ΔE < σ < kBT). Following the initial excitation of the highest eigenstate in the population of that state decay as shown in Figure a, again the population predicted by the NISE method approached the 50% high-temperature equilibrium population at long times. The SESH and HEOM both end at populations of ∼22%, and the predicted transfer is approximately twice as fast as for NISE, as the backward transfer from the lowest- to the highest-energy eigenstate is heavily suppressed. The TNISE in this case only decays to 30% and does not recover the full thermalization. At short times, the HEOM exhibit a fast nonexponential decay feature (T < 100 fs). This is attributed to the use of the average eigenstate population instead of the instantaneous eigenstate population in this case. For the site population following the initial excitation of the highest-frequency site (see Figure b) initial coherent oscillations are observed at early times for both NISE, TNISE, and HEOM, and again the equilibrium populations for the two methods deviate due to the high-temperature approximation employed in the NISE method. The TNISE relax to a final population between the NISE and HEOM results. The present results suggest that using the SESH for modeling quantum dynamics in FMO-like systems will be a significant advantage as compared to using the NISE approximation. The recovered results are very close to those found with the formally exact HEOM method. The TNISE performs worse even though it does present an improvement compared to NISE.

Figure 4

(a) The population of the highest-energy FMO-like (instantaneous) eigenstate following excitation of that state for the amide-like parameters according to the HEOM (red), NISE (blue), SESH (green), and TNISE (dashed black) approaches. (b) The population of the highest-energy FMO chromophore following initial excitation of that site according to the HEOM (red), NISE (blue), and TNISE (dashed black) approaches.

We now proceed to look at a more complex bacterial light-harvesting system, LH2. This system is characterized by two separate bands denoted the B800 and the B850 bands.[111,112] The B800 band is understood to consist of eight or nine weakly coupled bacteriochlorophyll molecules, while the B850 band consists of 16 or 18 strongly coupled bacteriochlorophyll molecules.[42,113−126] The number of molecules depend on the bacterial species.[127−129] Here, we again simplify this to a representative two-level dimer system, with one site representing a site in the B800 band and a site representing a site in the B850 band. The gap is set to 300 cm–1, which is larger than the thermal energy. At the same time, the frequency fluctuations are set to 100 cm–1, which is also quite large. The relatively small coupling of 47 cm–1, in this case, means that the eigenstates are essentially localized on the sites in contrast to the FMO-like system described above, where the coupling was comparable to the energy gap (|J| σ < ΔE < kBT). In this case the quantum dynamics at short times (T < 1 ps) is essentially identical for the HEOM and SESH approaches (see Figure a). The longtime behavior is, however, different, with the HEOM predicting a lower equilibrium population. This was also seen for the amide I–amide II system, albeit to a much smaller extent. The explanation is that for these systems the reorganization energy (λ = σ2/2kBT) is non-negligible compared to the thermal energy, and the final excitation is localized on one of the sites, which results in a Stokes shift of this low-energy site in the exact HEOM approach. In the SESH this effect is not included. For the water model and the FMO model the excitation on the low-energy eigenstate is much more delocalized, and the Stokes shift present is affecting both bath coordinates more or less equally, and then the quantum dynamics is less affected, as the change of the energy gap in the Hamiltonian is not affected by the lowering of the energy of both sites. Still, here the SESH provide a considerable improvement compared to the NISE approach, and we can easily understand why the applied approximation breaks down at longer times. For TNISE the decay is between that of the NISE and the HEOM as was the case for the FMO parameters. In the real LH2 complex the dynamics is actually much faster than seen here due to the presence of mixed B800/B850 states,[123] when all chromophores are included. It is also known that the excitation in the B850 band is very delocalized, and the effect of the Stokes shift is smaller in reality. We finally examine the quantum dynamics of the site populations predicted by the NISE, TNISE, and HEOM approaches as presented in Figure b. The results are similar for the first 100 fs, after which the decay in the HEOM become significantly faster than for both the TNISE and NISE approach. This can be understood as the time-scale for the bath coordinates is 100 fs, and following this time the bath has time to respond to the excitation of the quantum system. The initial effect is that the energy of the initially excited high-energy state is slightly lowered due to the bath reorganization, and at the same time the back-transfer of excitation from the low-energy state is suppressed, as there is not sufficient energy in the bath to support such transfer.

Figure 5

(a) The population of the highest-energy LH2, B800 bandlike (instantaneous) eigenstate following excitation of that state for the amide-like parameters according to the HEOM (red), NISE (blue), SESH (green), and TNISE (dashed black) approaches. (b) The population of the highest-energy LH2 B800-like chromophore following initial excitation of that site according to the HEOM (red), NISE (blue), and TNISE (dashed black) approaches.

The temperature dependence was further examined for the amide I–amide II parameters. The instantaneous eigenstate populations calculated for four different temperatures are presented in Figure . Both the SESH and the TNISE methods provide good results down to 150 K, but not at 70 K, where the thermal energy is smaller than the gap between the eigenstates. We conclude that the new methods can be safely applied when the temperature is not too much smaller than the width of the band. In the present approaches, the bath is assumed to behave in a classical manner, and it may thus be possible to improve the agreement using quantum correction factors[130,131] to correct for the quantization of the bath energy at low temperature.

Figure 6

Population of the amide I (instantaneous) eigenstate following excitation of that state for the amide-like parameters according to the HEOM (red), NISE (blue), SESH (green), and TNISE (dashed black) approaches at 75 K (a), 150 K (b), 300 K (c), and 600 K (d).

The TNISE method scaling the elements of the nonadiabatic coupling matrix is quite similar to the scaling method of Bastida et al.,[15] for which an improved scaling version was presented by Kleinekathöfer.[21] In the methods by Bastida and Kleinekathöfer a symmetrization of the nonadiabatic couplings was employed to account for the nonhermitian structure of the scaled matrix, while here a renormalization of the wave function was used for the TNISE method. In Figure the present methods are further compared with the improved scaling method (following the implementation as described in ref (21)) revealing that the TNISE results are at least for this case in better agreement with the exact HEOM result. For the site population, the decoherence of the improved scaling method is clearly too fast, while the correct thermal equilibrium is reached. The slower convergence of the SESH method due to its stochastic nature is evident.

Figure 7

(a) The population of the amide I (instantaneous) eigenstate following excitation of that state for the amide-like parameters according to the HEOM (red), NISE (blue), SESH (green), TNISE (dashed black), and improved scaling (magenta) approaches at 150 K. (b) The population of the amide I site following excitation of that site. The same color scheme is used as in (a), but there are no data for SESH.

The TNISE, the improved scaling method, and the SESH all require transformations to the instantaneous adiabatic basis. This may cause problems when the instantaneous eigenstates between neighboring times steps need to be matched. This problem also exists for the original surface-hopping method, and small time steps may be required, in particular, for systems involving conical intersections.[8] One could consider improving the behavior, in particular, of the TNISE method by adapting the temperature dependence of the scaling factor. In the SESH one could likewise try to use the quantum thermal energy average. This is, however, not so trivial, as the kinetic energy of a quantum harmonic oscillator depends on the frequency of the oscillator and not just on the temperature. More explicit knowledge on the degrees of freedom is, thus, needed than for the present approximations, where the bath fluctuations are effectively treated through one effective bath coordinate. For anharmonic bath coordinates, the assumptions are even more complex.

The two new methods for calculating quantum dynamics with thermalization were found to predict the dynamics quite well when the system parameters are not too large compared to the thermal energy. As such the approaches may work very well, in particular, for modeling thermal relaxation in vibrational modes or when the thermal effect is just a secondary effect and not the primary point of interest. The SESH method has two important drawbacks. First, because of its stochastic nature it converges quite slowly, and a large number of realizations are needed. Second, the method, at least as developed here, only consider the dynamics in terms of instantaneous eigenstates. Coherent superpositions may be possible to consider combining the present approach with recently developed surface-hopping approaches that include coherent superpositions.[46,54] The TNISE method performed very well for the vibrational systems and required the use of smaller time steps to guarantee identification of the right eigenstates in the electronic systems, where the performance was also generally worse than for the SESH method. Still the TNISE approach is very simple in implementation, and it converges faster than the SESH. This will, in particular, be important if the methods are combined with response function calculations to model coherent multidimensional spectroscopic signals (see refs (63, 76, 87, 90, 98, 105, 123, and 132−137)) and when a large number of classical coordinates are involved.

The new methods have some similarity to the classical path approximation (CPA),[11] where one first integrates the time-dependent Schrödinger equation for a predetermined classical trajectory and then use Tully’s hopping procedure to estimate populations to make sure one reaches detailed balance. This method is also only applicable for small reorganization energies. The SESH should be easy to combine with a recently developed generalization for treating coherences in the fewest-switches surface-hopping framework.[54]

Conclusions

In summary, we developed a new surface-hopping-based approach for simulating quantum dynamics that accounts for thermalization and a rescaled NISE approach, both without explicit feedback to the bath trajectory. The new approaches were compared for a number of coupled two-level systems mimicking the expected conditions in typical vibrational and electronic systems. Clearly, they both considerably improve the thermalization dynamics as compared to the simple NISE approach, while the bath relaxation is neglected resulting in deviations from the exact HEOM results when the reorganization energy is large. Here, we only considered very simple model systems, which can also easily be modeled with the HEOM approach; the new approaches can both easily be extended to systems with general Hamiltonian fluctuations, including fluctuations of the couplings, correlated couplings between sites, and non-Gaussian fluctuations. They can further be extended to complex systems including hundreds of coupled sites[69,138] or systems with conical intersections,[137] as it is the case for the simpler but not correctly thermalizing NISE approach.[60] The new SESH method was found to provide the best overall accuracy especially when the energy gaps are large. However, this method does converge slower than the new TNISE method, which will thus be the method of choice when the energy gaps are smaller than the thermal energy. One other important limitation of the SESH method is that it treats the system as always being in an eigenstate, which complicates implementation for simulating spectroscopic signals.[87]

The new approaches can be implemented to model pump–probe, two-dimensional spectroscopies,[76,77] and time-resolved fluorescence[139,140] spectroscopies that report on the quantum dynamics. The methods may further prove useful in the development of new methods for modeling spectroscopies involving both vibrational and electronic degrees of freedom.[141,142] With the development of broad-band laser sources,[143,144] experimental spectra spanning very broad frequency regions have become available, and the new method has the potential to allow modeling such spectra. As the Achilles heel of the methods is that the system–bath interaction should not be too large compared to the thermal energy the most straightforward solution to this problem is to include all strongly coupled modes explicitly in the quantum Hamiltonian. Obviously, there is a limit for how many degrees of freedom can be included; however, one could imagine combining the present approaches with the multiconfiguration time-dependent Hartree[20] method that allows treating very large quantum systems.