Extending the applicability of Redfield theories into highly non-Markovian regimes

Useful approximate methods for the description of quantum dynamics and relaxation can often circumvent the large computational expense of numerically exact approaches while maintaining quantitative accuracy in certain regions of parameter space. The general applicability of such methods, however, is often limited and their domain of validity difficult to assess. The most widely used approximate approaches fall into two broad and general classes. The first class of methods involves techniques that employ uncontrolled approximations to yield dynamics which are non-perturbative in the various couplings (e.g., intra-system or system-bath) that characterize the problem. The second class of methods is systematically perturbative in a well-defined coupling parameter and is free from further classical or semiclassical approximations, at least for simply defined models such as the spin-boson model.

One of the most celebrated perturbative techniques is a lowest-order treatment of the system-bath coupling, known traditionally as Redfield theory.^1–3 As we will show later, the relevant dimensionless parameter characterizing the accuracy of Redfield theory is $η=max 2 λ / β ω c 2 , 2 λ / π ω c$⁠, where λ is the nuclear reorganization energy, β = 1/k_BT is the inverse temperature, and ω_c is a characteristic bath frequency (we henceforth work in units with ħ = 1). Redfield theory becomes unreliable when η ≳ 1. We emphasize that η is controlled by multiple bath parameters and, in particular, low-frequency degrees of freedom (small ω_c) limit the range of accessible reorganization energies. Indeed, violations of this condition explain the failures of Redfield theory found by Ishizaki and Fleming⁴ for certain models of excitation energy transfer, which appear to be characterized by low-frequency protein baths.

While the very lowest frequency degrees of freedom are thus most problematic for Redfield theory to handle (even in its non-Markovian forms), it is often the case that nuclear modes of such frequencies are effectively frozen on the time scale of relevance for the system’s dynamics. In this regard, the key function of such modes is simply to provide static energetic disorder for the more rapidly evolving degrees of freedom. This suggests a methodology whereby the very low-frequency phonons are approximated as static (and treated non-perturbatively as a source of static disorder), while the remaining portion of the bath is treated dynamically within Redfield theory. Here, we develop this “Redfield theory with frozen modes” (Redfield-FM) method and show that it greatly extends the applicability of Redfield theory into highly non-Markovian dynamical regimes at essentially no change in computational cost.

The outline of this paper is as follows. In Sec. II, we introduce the theoretical background for the Redfield equations and the derivation and general properties of the Redfield-FM extension. In addition, we also introduce the spin-boson Hamiltonian as the model system on which we test the methods developed in this paper. Section III A presents the computational details in the implementation of the Redfield-FM method, while Sec. III B presents illustrative results for the method. Sec. III C introduces the generalization of the Redfield-FM approach to a multi-site system and presents representative results. In Sec. III D, we relax the mode-freezing approximation via the derivation and implementation of a dynamical hybrid method (hybrid-Redfield) that combines Redfield dynamics for the high-frequency part of the bath coupled to the electronic system and Ehrenfest dynamics for the low-frequency modes. In Sec. IV, we conclude.

First, we briefly describe the model system we use to test the approximations developed in Secs. II C and III. This allows us to define notation and parameters that will be used in our numerical comparisons. We focus on the well-known spin-boson model, which consists of a two-level system coupled linearly to a harmonic bath. This model has been extensively used to investigate a wide variety of relaxation, charge, and energy transport processes in condensed phase systems.⁵ 

The total Hamiltonian is divided into system, bath, and interaction components, H = H_sys + H_bath + V. The system Hamiltonian takes the form

where σ_i, i = {x, y, z}, are the Pauli matrices, 2ε is the energy difference, and Δ is the coupling between the two electronic sites, which is here assumed to be static. The bath Hamiltonian consists of an infinite set of harmonic oscillators,

Finally, the system-bath interaction couples the electronic states linearly to coordinates of the bath oscillators,

Physically, the system-bath coupling acts as a (quantum) fluctuating field that shifts the origin of the bath harmonic oscillators by a magnitude that depends on the system’s electronic state and the strength of the coupling.

The spectral density, which completely determines the coupling between the bath and the system, is taken to be Ohmic with a Lorentzian cutoff (Debye form),

The cutoff frequency, ω_c, characterizes how quickly the bath relaxes toward equilibrium, while the reorganization energy, $λ= π − 1 ∫ 0 ∞ dωJ ( ω ) /ω$⁠, characterizes the energy dissipated by the environment after a Franck–Condon transition between electronic states. It is important to note that the methods studied here are neither limited to the spin-boson model nor to the Debye form for the spectral density.

Because of its simplicity in the time domain, we employ the time-local (i.e., time-convolutionless) form of the generalized Redfield equations. A full derivation of these equations is contained in  Appendix A. Here, our aim is to highlight important but often overlooked aspects pertaining to the applicability of the Redfield approach. For the spin-boson model, the time-local version of the Redfield theory takes the following form:

where all operators except the reduced density matrix (RDM) are evolved in the interaction picture, O(t) = e^{−i(H_sys+H_bath) t}Oe^{i(H_sys+H_bath) t}, and the free bath correlation function is given by

By going to the interaction picture with respect to H_sys + H_bath (to eliminate the free-evolution) and formally integrating equation of motion (5), one finds

Examination of the function $∫ 0 t d τ 1 ∫ 0 τ 1 d τ 2 C ( τ 2 )$ reveals the natural dimensionless parameter that determines the limit of validity of Redfield theory. In general, even in the non-Markovian case, we expect that the ordering of terms in expansion (7) is governed by the function $∫ 0 t d τ 1 ∫ 0 τ 1 d τ 2 C ( τ 2 ) =ηg ( t )$⁠, where η is a dimensionless constant and g(t) is a function expressed in terms of a scaled, dimensionless time variable. In the high-temperature limit (βω_c ≪ 1), where C(t) ≈ (2λ/β) e^−ω_ct, it is easy to show

In the low-temperature limit (βω_c ≫ 1), we assume that the low-frequency behavior of the spectral function dominates, so we choose J(ω) = (2λω/ω_c) e^{−2ω/πω_c} as an approximation to the Debye form in Eq. (4) that exactly matches the value of λ and its low-frequency asymptotic behavior. Using this spectral density, one can show

Thus, for a Debye spectral density, Redfield theory will be reliable as long as $max 2 λ / β ω c 2 , 2 λ / π ω c$ is not significantly larger than unity.⁶ It should be noted that recent work purported to be in the Redfield limit actually violates the above condition.⁷ As long as the relevant energy scales in the system Hamiltonian are not too large, we expect the above to hold. In cases where the system’s bare energy difference ε is the largest energy scale in the problem, the dynamics will be mediated by multi-phonon processes which are a challenge for lowest-order Redfield-like theories. However, in this limit, the problem acquires an increasing amount of “pure-dephasing” character, for which the time-local version of Redfield theory provides an exact multi-phonon resummation.

As discussed in the Introduction, low-frequency bath modes ω_k which lead to a violation of the validity of the Redfield theory frequently evolve so slowly as to be effectively static on the electronic time scale. Here, we develop the “Redfield theory with frozen modes” (Redfield-FM) method, based on the physically appealing notion of dividing modes into a low-frequency portion (treated as static disorder), and a high-frequency bath (treated by time-local Redfield theory). The separation of modes used here mirrors that utilized in previous work on a Förster-like dynamical hybrid approach.⁸ The approach presented here is not in any way limited to time-local Redfield theory (nor even to any specific flavor of Redfield theory). However, due to pathologies associated with a strictly Markovian Redfield theory, we suggest certain adjustments to the partitioning algorithm, as discussed in  Appendix B.

First, it is advantageous to assume that the total density matrix is multiplicatively separable into weakly interacting parts, i.e., ρ_tot(t) ≈ ρ_slow(0) ρ_sys+fast(t), where ρ_slow(0) is the density matrix of the frozen low-frequency “slow modes” and ρ_sys+fast(t) is the density matrix for the system and high-frequency “fast modes.” As in Ref. 8, a splitting function, S(ω), divides the spectral density into two components, J(ω) = J_slow(ω) + J_fast(ω), where

and

Here, we take the same form of the splitting function as that suggested in Ref. 8, namely,

which, by virtue of its smoothness, avoids problems associated with long-time oscillatory tails in the bath correlation function.⁹ While the above splitting induces no errors if the dynamics are treated exactly, ω^∗ clearly serves as a free parameter that allows one to tune the optimum percentage of frozen bath modes, and hence the accuracy of the results if the dynamics are treated within our approximate method. The utility of the present method is greatly enhanced if a physical a priori prescription for choosing ω^∗ based only on the parameters of the initial Hamiltonian can be put forth. In this work, we choose ω^∗ = max[ω_c, ω_R/4], where $ω R =2 ε 2 + Δ 2$ is the system Rabi frequency. This choice for ω^∗ is simple, yields non-trivial improvements over standard Redfield theory, and may easily be generalized to multiple electronic states. Physically, this choice partitions the bath into modes that evolve more slowly than the system (to be treated as frozen) and modes that evolve faster than the system (to be treated via Redfield theory). However, it should be noted that this choice is not always optimal. Future work will be devoted to the goal of arriving at an optimal choice of ω^∗. Other choices for ω^∗ that fit within the general physical guidelines discussed above will be discussed in Sec. III.

With a prescription for choosing ω^∗ in hand, it is possible to separate the fast and slow portions of the Hamiltonian, starting with the interaction, $V fast = σ z ∑ k ∈ fast c k Q k$ and $V slow = σ z ∑ k ∈ slow c k Q k$⁠. Regrouping terms, it is evident that freezing the slow part, V_slow, yields a classical reorganization energy that renormalizes the bias for every realization of the bath’s initial conditions. The modified total Hamiltonian now takes the following form:

where the classical reorganization energy is defined as $λ c l ( 0 ) = ∑ k ∈ slow c k Q k ( 0 )$ and the set of Q_k(0) is sampled from a bath distribution function after the discretization of J_slow(ω). Physically, each realization of the frozen bath degrees of freedom constitutes a local, rigid environment that modifies the site energies for the system Hamiltonian. The time-local Redfield dynamics, under the Hamiltonian in Eq. (13) with an interaction given by $V fast = σ z ∑ k ∈ fast c k Q k$⁠, are subsequently ensemble averaged with respect to the slow frozen modes. Thus, there are two important differences for the Redfield equations used in each realization of the frozen modes: (i) the bias is given by $ε ̃ ≡ε+ λ c l ( 0 )$ and (ii) the bath correlation function given by Eq. (6) is modified, with J(ω) replaced by J_fast(ω).

To account for the classical frozen modes, the nonequilibrium population dynamics takes the following form:

where ρ_slow(P, Q, 0) could be either the classical distribution function or the Wigner transform of the equilibrium density operator of the slow bath degrees of freedom, and ρ_sys+fast(t) is the reduced density matrix of the system and the fast bath degrees of freedom. Observables such as Tr_sys+fast[ρ_sys+fast(t) σ_z] may then be calculated via the Redfield equations.

To understand the relaxation processes with mode freezing for finite ω^∗, it is first useful to investigate the effect of the approximation at its most extreme, namely, the adiabatic limit, where all bath modes are assumed to be static (ω^∗ → ∞). In this limit, we are effectively in the Born–Oppenheimer regime, where excitations in the electronic subspace move along the potential energy surface determined by the frozen reservoir. Analytical evaluation of the trace within Eq. (14) leads to the following expression for the nonequilibrium population dynamics:

where $ξ=2 ε ̃ 2 + Δ 2$ is the Rabi frequency of the modified system Hamiltonian. The integration over the ensemble of equilibrium configurations of the bath reduces to averaging over different values of λ^cl(0) that are consistent with the bath distribution function. For some realizations of the bath, Eq. (15) recovers the Rabi oscillations characteristic of the isolated system if λ^cl(0) = 0. Conversely, when λ^cl(0)≠0, the population starts from 1 at t = 0 and oscillates between $( ε ̃ 2 + Δ 2 ) / ( ε ̃ 2 + Δ 2 ) =1$ and $( ε ̃ 2 − Δ 2 ) / ( ε ̃ 2 + Δ 2 )$⁠. Taking for simplicity ε = 0, one notes that the lower bound of the population oscillations increases with increasing λ^cl(0), approaching 1 as λ^cl(0) → ∞. This limit corresponds to an infinitely rigid bath that completely localizes the excitation on its initial site.

Averaging over different realizations of the slow modes decreases the amplitude of oscillations in the population dynamics due to the decoherence between the functions with distinct oscillation frequencies. In general, Redfield theory has difficulty describing non-Markovian, multi-step relaxation dynamics. However, when $λ slow = π − 1 ∫ 0 ∞ dω J slow ( ω ) /ω$ is sufficiently large, the full dynamics produced by Redfield theory with frozen modes at finite ω^∗ will include both a slow, perhaps oscillatory component as well as a more rapid decay induced by the high-frequency modes in J_fast(ω). These qualitative considerations suggest this approach may correct certain deficiencies of conventional Redfield-like approaches. In Sec. III, we test the approach quantitatively.

An additional concern regarding any approximate dynamical theory is whether it leads to the appropriate long-time limit. While the asymptotic behavior in the standard (time-local and time-nonlocal) Redfield equations can only be studied numerically, it is simple to show that Markovian Redfield theory within the secular approximation, where the populations and coherences (in the system’s eigenbasis) are assumed to evolve independently of each other, obeys detailed balance with respect to the isolated system.¹⁰ In this case, one may write the long-time limit of the RDM as $lim t → ∞ ρ sys ( t ) = e − β H sys / Tr sys [ e − β H sys ]$⁠, which is only correct in the weak system-bath coupling limit. Making the analogous approximations for the Redfield-FM approach, calculation of the long-time limit of system observables requires tracing over the appropriate operator and averaging over the static disorder imposed by the arrested modes. Hence, the equilibrium value of the population difference in the diabatic basis takes the form

where $ξ ̃ = ε ̃ ( Q ) 2 + Δ 2$ is the (positive) eigen-energy of the isolated system for each realization of the environmental disorder. In the case of an unbiased system (ε = 0), the integrand in the above expression is odd with respect to Q, leading to 〈σ_z(t → ∞)〉 = 0, agreeing with the result obtained via the Markovian secular Redfield approach. The two methods yield different results for biased cases, ε≠0. While the current discussion of detailed balance is strictly applicable only to Markovian Redfield-FM under the secular approximation, the result in Eq. (16) should serve as an approximate guide to the long-time limit of the populations in the Redfield-FM scheme presented here.¹¹ In this regard, it should be noted that the long-time behavior exhibited by the Redfield-FM approach appears to be in better agreement with exact results than is standard Redfield theory.

Finally, we remark that while the idea of dynamically arrested modes has been used in various contexts leading to much notable work, our approach exploits this idea in a novel way. Previous implementations of this idea have ranged from a partial freezing of the modes, as in the case of electron transfer in proteins,¹² to complete arrest of the environment, leading to the Gaussian disorder model for studying charge carrier mobility in disordered organic semiconductors.¹³ As with the two examples cited above, other uses of the mode-freezing approximation (to the best of the authors’ knowledge) have been concerned with the modification of incoherent, classical transfer rates of excitations among sites. In contrast, we employ the partial arrest of the bath in conjunction with a fully quantum mechanical treatment to obtain real-time dynamics that account for both coherent and incoherent motion of the excitations.

In the following, we compare the numerically exact population dynamics reported by Thoss et al.¹⁴ for the spin-boson model with a Debye spectral density and the initial condition Γ(0) = |1〉 〈1|exp(−βH_bath)/Z_bath with the results obtained from the Redfield-FM method. Subsequently, we examine the effect of relaxing the mode-freezing approximation by treating the low-frequency modes dynamically via the Ehrenfest method. We call this latter approach the hybrid-Redfield method, in analogy with the previously developed hybrid-NIBA method.^8,9

To treat the frozen portion of the spectral density, J_slow(ω), we have discretized the bath into f = 300 modes with frequencies and couplings given by^9,15

and

Initial conditions for the reservoir of frozen modes were sampled from a Wigner distribution. Sampling from this distribution becomes particularly important at very low temperatures, where quantum effects become significant. However, for most cases, sampling from a Boltzmann distribution is sufficient since the modes being sampled are always of low-frequency. For convergence, up to 10⁴ trajectories have been run for the results presented.

As mentioned in Sec. II B, the validity of Redfield theory is limited to the small η regime. Fig. 2(a) shows the results for a slow bath (ω_c = 0.25) with small reorganization energy (λ = 0.25) at low temperature (β = 5.0); here and in the following, all energies are in units of Δ. In spite of the slow bath, the dimensionless applicability parameter is only slightly larger than unity (η = 1.6) suggesting that Redfield theory should be reasonably accurate, in agreement with the numerical results. The Redfield-FM method provides an even better estimate of the dynamics, almost quantitatively correcting the already accurate Redfield dynamics.

Fig. 2(b) considers a biased system (ε = 1.0) with the same parameters, except at much higher temperature (β = 0.5), yielding an applicability parameter which is significantly larger than unity (η = 16). Here, it is evident that the Redfield dynamics relaxes far too quickly, suppressing the coherence and missing the slower relaxation process revealed by the exact dynamics. The improvement afforded by the Redfield-FM method compared to standard Redfield theory is clear. In particular, the Redfield-FM approach accurately reproduces the short- to intermediate-time dynamics, the frequency of the oscillations, and the initial rate of decoherence. The terminal decay rate is slightly underestimated due to the mode-freezing approximation, as discussed in Section II B. Yet despite these shortcomings, the improvement derived from a simple scheme like the Redfield-FM approach with the numerical complexity of the original Redfield theory is noteworthy.

For cases exemplified by Fig. 2(c), serious problems such as the violation of the positivity of the RDM dynamics can occur within standard (non-secular) Redfield theory. Fig. 2(c) corresponds to a slow bath (ω_c = 0.25) and a large reorganization energy (λ = 5.0) again at high temperature (β = 0.5), for which the applicability parameter is very large (η = 320). Despite the evident failure of the Redfield equations to even maintain positivity, the Redfield-FM method is able to correct the positivity issue and almost quantitatively reproduce the two-step relaxation process in the exact dynamics up to intermediate times.

It is possible to understand the surprising success presented in Fig. 2(c) in the context of the analysis of Section II B. Using the definitions given in that section, the effective parameters for the Redfield equation are λ_fast = 2.5, ω_c = 0.5, yielding η = 40. Although η ≫ 1, the reduction by an order of magnitude from the initial value, η = 320, is sizeable, and likely responsible for solving the positivity problem evident in the bare Redfield dynamics. The reproduction of the two-step relaxation process is a direct result of the trapping effect that arises from freezing a large portion of the low-frequency bath in the presence of large coupling. This example indicates that the trapping effect can partially reproduce slow relaxation dynamics associated with strong system-bath interactions.

Figure 2(d) shows the regime of intermediate bath speeds (ω_c = 1), large reorganization energy (λ = 2.5), and intermediate temperature (β = 1). In contrast to Fig. 1(c), the Redfield-FM method is not capable of significantly improving the Redfield dynamics in this regime, missing the two-step relaxation process visible in the exact dynamics. In light of the previous case, it is evident that the slowing down of the RDM dynamics can be caused by freezing a large portion of the strongly coupled modes, an effect which is absent in this case. In this case, λ_fast = 2.1 and λ_slow = 0.4, which indicates that most of the reorganization energy is included already in the high-frequency portion of the bath. In such cases, the Redfield-FM method yields results that are similar to bare Redfield theory.

We now address the dependence of the dynamics on the choice of ω^∗. Eschewing the simple criteria for choosing ω^∗ presented above, one may ask how closely the Redfield-FM dynamics can be made to agree with exact dynamics when ω^∗ is allowed to vary. To address this question, we include two extreme cases in Fig. 3. First, Fig. 3(a), which corresponds to the same parameters as those of Fig. 2(c), shows that optimization of ω^∗ can result in quantitative agreement between the Redfield-FM result and the exact dynamics. Such agreement may be understood as the result of fortuitous cooperation between strongly dissipative Redfield dynamics that damp the frozen mode-generated oscillations and the trapping effect from the mode-freezing approximation that prevents immediate relaxation to the equilibrium population. Conversely, Fig. 3(b), which corresponds to the parameters in 2(d), is an example of when perfect agreement is impossible. Clearly, attempts at optimizing ω^∗ result in better agreement of the two-step relaxation process at the cost of long-lived oscillations, a direct result of including a large fraction of modes into the slow part of the bath. In freezing a sufficiently large portion of the reservoir to reproduce the trapping effect, λ_fast is reduced to the point where the Redfield dynamics are no longer sufficiently dissipative to damp the frozen mode-generated oscillations. In addition, $λ slow = π − 1 ∫ 0 ∞ dω J slow ( ω ) /ω$ is also not large enough to ensure that the oscillations dephase sufficiently rapidly. Overall, it is clear that although it may be possible to optimize the results, the simple initial criteria presented represent a robust approach to frozen mode dynamics that essentially always yields results that are as good or better than bare Redfield dynamics without a significantly increased computational cost.

While the spin-boson model serves as a simple yet nontrivial model to test of the performance of new approximate dynamical theories, the usefulness of a method also largely lies in its applicability to multi-site systems where numerically exact methods become difficult to perform. Generalization of the Redfield-FM approach to multi-site systems is straightforward. To illustrate this, we focus on the theoretically¹⁶ and experimentally¹⁷ well characterized example of the FMO complex.

Within the standard FMO model, the total Hamiltonian takes the form H = H_sys + H_bath + V, where

with parameters E_i and J_ij taken from Ref. 18. The bath Hamiltonian consists of independent reservoirs for each chromophore,

The system-bath coupling is bilinear, same as in the spin-boson model,

where c_i,k is the coupling constant between the ith site and the kth mode in the local bath. All local baths are assumed to be equivalent, characterized by Debye spectral densities, with λ = 35 cm⁻¹.

To define the splitting frequencies that generalize our previous choice, one first defines a set of “Rabi” frequencies associated with each pair of sites, {ω_R,ij}, where $ω R , i j =2 ( E i − E j ) 2 + J i j 2$⁠. The splitting frequency for each bath is then taken, as before, as $ω i ∗ =max [ ω c , i , max [ ω R , i j ] / 4 ]$⁠. As we will demonstrate, despite the fact that this definition is essentially unaltered from that used to treat the two-site cases, the results are in remarkable agreement with exact calculations for more complex multi-site examples.

Although Redfield theory has been criticized for its inability to recover the correct dynamics in prototypical electronic energy transfer systems where intrasystem and system-bath couplings are comparable,⁴ its performance for the FMO model is surprisingly good as long as temperature is low and the bath relaxation time scale is short. This is in harmony with the discussion contained in Sec. II B. In Fig. 4, we consider one such favorable case corresponding to a fast bath with correlation time τ_c = 50 fs (⁠$ω c = τ c − 1$⁠) and low temperature T = 70 K. As is clear from the figure, the Redfield equations recover the dynamics quantitatively, including the correct oscillation frequency, amplitude, and long-time limit for all but sites 6 and 7. Fig. 5 corresponds to the case of high temperature (T = 300 K) and sluggish baths with correlation time τ_c = 166 fs. Here, the parameter η is significantly larger, rendering this a much more difficult parameter regime for the Redfield approach, which overestimates the decoherence and leads to incorrect behavior in the long-time limit for all but one of the site populations. The Redfield-FM approach corrects the dynamics quantitatively, alleviating the difficulties associated with a sluggish bath.

The examples presented here clearly illustrate the ingredients that determine the success or failure of standard Redfield theory. In contrast to the pervasive claim that Redfield theory should fail when the reorganization energy exceeds the intersite couplings, we have shown in Fig. 4 that as long as the temperature is relatively low and the bath responds rapidly, Redfield theory can be quantitative even in the “intermediate” coupling regime. On the other hand, when the parameter η becomes large, in particular in cases when the bath is sluggish, the Redfield-FM approach quantitatively corrects standard Redfield theory. The fact that this is the case for multi-site examples with no adjustment to the criteria for choosing ω^∗ lends credence to the robustness of the Redfield-FM methodology.

On first inspection, the mode-freezing approximation appears extreme. To thoroughly assess its effect, we develop a dynamical hybrid method in which we evolve the previously frozen low-frequency modes in J_slow(ω) via classical Ehrenfest dynamics. The derivation and the implementation details of this approach may be found in  Appendix C.

This hybrid-Redfield method is similar in spirit to the successful hybrid-NIBA developed and implemented in Ref. 9. Evolution of the low-frequency modes using Ehrenfest dynamics in such hybrid approaches only requires two assumptions: (i) that Γ(t) ≈ ρ_slow(t) ρ_sys+fast(t) and (ii) that the motion of the low-frequency modes is well-captured by classical mechanics. For such a factorization to be valid, the reorganization energy due to the low-frequency bath needs to be small, i.e., $λ slow = π − 1 ∫ 0 ∞ dω J slow ( ω ) /ω<Δ$⁠. The applicability of classical dynamics relies on the low energies of the reservoir modes and sufficiently high temperatures that help suppress quantum effects.^9,19 However, even when the Ehrenfest approximation is valid, problems may arise. Most prominent among these is that the final populations approach those of the infinite-temperature limit.¹⁹ 

In contrast to the Hamiltonian derived under the mode-freezing approximation in Eq. (13), the modified Hamiltonian that needs to be treated via the Redfield equation in the hybrid-Redfield method is time-dependent,

where the disorder due to the low-frequency bath is no longer static as it is in the Redfield-FM method, but rather dynamic, namely, $λ c l ( t ) = ∑ k c k Q k ( t )$⁠.

Since the system part of this Hamiltonian is nondiagonal and time-dependent, evolution with respect to the system Hamiltonian requires diagonalization at every time step, significantly increasing the computational cost associated with the method proposed here. The need to evolve the low-frequency bath also adds to the computational cost of the approach. Importantly, under the mode-freezing approximation, we circumvent these costly requirements. This means that aside from the trivial cost of parallelization for the ensemble averaging over the slow bath, the Redfield-FM method scales as gracefully with system size as the original Redfield equation.

For completeness, we remark that the nonequilibrium population dynamics under the hybrid-Redfield approximation now takes the form

Fig. 4 shows two sets of parameters for which the hybrid-Redfield scheme yields results that illustrate the issues at play in comparing the hybrid-Redfield approach to the Redfield-FM method. Extensive testing of the hybrid method suggests that an approximately optimal form for the splitting frequency can be taken as

Physically, this form encodes the interplay between the Redfield and Ehrenfest methods, favoring a larger portion of the modes to be treated classically with increasing Rabi frequency ω_R, which is a measure of how rapidly the electronic system evolves. Furthermore, this form of ω^∗ ensures that in the limit of small λ and large ω_c, the hybrid method correctly reproduces the more appropriate Redfield dynamics, whereas in the limit of large λ and small ω_c, it reproduces the Ehrenfest results. It is expected that generally nontrivial results may be obtained from this method for cases where ω^∗ ∼ ω_c, as is the case for the choice of ω^∗ used in the Redfield-FM approach.

Fig. 6(a) corresponds to the parameters in Fig. 2(c) and illustrates that by means of the suggested form for ω^∗, hybrid-Redfield automatically tunes itself to yield nearly optimal results achievable from the two parent methods. This example, for which $ω h y ∗ ≫ ω c$⁠, illustrates that the hybrid-Redfield method trivially reproduces the Ehrenfest result when it is appropriate. It is noteworthy that the Redfield-FM method obtains similar agreement at a much lower computational cost without evolving the reservoir modes, indicating that dynamic treatment of these modes may not be generally necessary. Indeed, it is rather remarkable that the Redfield-FM approach basically recapitulates the Ehrenfest results even though no Ehrenfest dynamics are used.

Fig. 6(b) shows the analogue of Fig. 2(d), where the Redfield-FM method fails to correct the Redfield dynamics. In contrast, the hybrid-Redfield results are in very good agreement with the exact dynamics. Indeed, the hybrid method is able to qualitatively and almost quantitatively reproduce the shape of the two-step relaxation process evident in the exact dynamics, an effect that both Ehrenfest and Redfield dynamics independently miss.

As the above considerations indicate, there are cases where the dynamical hybrid-Redfield method can provide a substantial improvement over the Redfield-FM method, albeit at a much higher computational cost. In most regions of parameter space we have studied, however, we find that hybrid-Redfield theory offers little accuracy gain over the Redfield-FM approach. Thus, the benefits of the hybrid-Redfield approach compared to the Redfield-FM method do not justify its use when accuracy and cost are factored together.

In this work, we have presented a new scheme for simulating dynamics in quantum dissipative systems. Our approach, which we call the Redfield-FM method, recognizes that standard Redfield theory becomes inaccurate for slow bath degrees of freedom. By partitioning the bath into high- and low-frequency components, we propose solving the Redfield equations for the high-frequency partition in the statically disordered field of the low-frequency components. Such an approach may greatly increase the accuracy of Redfield theory in highly non-Markovian regimes at essentially the same computational cost. In addition, we find that this simple approach can fundamentally cure positivity problems associated with standard non-secular Redfield theory. The straightforward generalization of the Redfield-FM method to multi-site models and its application to the FMO complex have proven efficient and highly successful, underlining the method’s capacity to tackle problems of significant complexity. We have further discussed a scheme (the hybrid-Redfield approach) whereby the previously frozen degrees of freedom are instead evolved with classical Ehrenfest dynamics. While this method can improve upon the dynamics as described by the Redfield-FM approach, the increase in accuracy is incremental and comes at a significantly larger computational cost. Overall, while the Redfield-FM method does not cure all of the ills of Redfield theory, it does provide a simple and efficient framework for improving its accuracy and range of validity, especially for sluggish bath degrees of freedom such as those implicated in biological energy transfer.

D.R.R. acknowledges support from NSF No. CHE-1464802. A.M.C. thanks Hsing-Ta Chen for useful discussions. T.C.B. was supported by the Princeton Center for Theoretical Science.

Here, for completeness, we review the derivation of the Redfield equations. For a more detailed discussion of Redfield theory, we refer the reader to Refs. 20 and 21.

In the following development, we utilize a projection operator technique to derive an equation of motion for the reduced density matrix (RDM) of the system, defined as ρ(t) = Tr_bath[Γ(t)], where Γ(t) = e^−iHtΓ(0) e^iHt and Γ(0) is the initial density matrix of the full system and bath. Moreover, we assume that the initial condition for the (total) density matrix contains no system-bath correlation, such that Γ(0) = ρ(0) ρ_bath, where ρ(0) is an arbitrary Hermitian system operator, ρ_bath = e^−βH_bath/Z, Z = Tr_bath[e^−βH_bath], and β = 1/k_BT is the thermal energy. Treatment of general initial conditions is also possible via the projection operator technique at the expense of the introduction of additional inhomogeneous terms in Eq. (5).^22–24 In the following, we ignore initial correlations, but note that their inclusion in the present framework is straightforward.

We start from the Liouville equation for the full density matrix in the interaction picture where the total Hamiltonian is divided into a zeroth order part and an interaction part, H = H₀ + H₁, such that

where Γ_I(t) = e^iH₀tΓ(t) e^−iH₀t and $L I ( t ) = [ e − i H 0 t H 1 e i H 0 t , … ]$⁠. To obtain the dynamics of the RDM, we define a projection operator of form $P≡ ρ bath Tr bath [ … ]$ with $Q≡1−P$⁠. We note that action of $P$ on the full density matrix followed by trace over the bath results in the RDM in the interaction picture, $ρ I ( t ) = Tr bath [ P Γ I ( t ) ]$⁠. Using these definitions, we obtain the following exact equations of motion:

Formal integration of Eq. (A2) yields

where $g ( t , τ ) = exp + [ − i ∫ τ t d s Q L I ( s ) ]$ and the time ordering (+) implies that time arguments increase from right to left. Substitution of this expression in Eq. (A1) results in the Nakajima–Zwanzig equation,^22,23 which is expressed in terms of the time convolution of a memory term with ρ_I(t) at earlier times as

where $K ( t − τ ) = Tr bath [ L I ( t ) g ( t , τ ) Q L I ( τ ) ρ bath ]$ is the (time-nonlocal) memory function.

If, instead, we use the formal solution of Eq. (A1) to evolve Γ_I(t) backwards in time to an earlier time τ, we obtain

where $G ( t , τ ) = exp − [ i ∫ τ t d s L I ( s ) ]$ and the time ordering (−) requires that time arguments increase from left to right. Replacing this expression in Eq. (A4) and solving for $Q Γ I ( t )$ yields

where

We note that a crucial requirement for the validity of this derivation is the existence of [1 − Σ(t)]⁻¹.

Substitution of Eq. (A7) into Eq. (A2) and subsequent trace over the bath degrees of freedom results in the following time-local equation of motion for the RDM:²⁴ 

where $R ( t ) =−i Tr bath [ L I ( t ) [ 1 − Σ ( t ) ] − 1 ρ bath ]$ is the (time-local) rate function.

The expression for the dynamical evolution in either the time-nonlocal (Eq. (A5)) or time-local (Eq. (A9)) form is exact but prohibitively difficult to evaluate without resorting to approximation schemes, such as truncated generalized cumulant expansions. Perturbative expansion to second order in the system-bath coupling (where H₁ = V from Eq. (3)) results in a non-Markovian generalization of the Redfield theory. Alternatively, one may derive both forms of the Redfield equations via resummations of differently time-ordered cumulants.^25,26,33 These derivations explicitly show that both forms of generalized Redfield theory account for non-Markovian behavior and have similar applicability requirements.^26–28 Specifically, since Redfield theory is tantamount to second-order perturbation theory in the system-bath coupling, truncation at low order is only accurate for η < 1, where $η=max [ 2 λ β ω c 2 , 2 λ π ω c ]$ is the validity parameter introduced in Sec. II B. Despite this restriction, the Redfield equations have been shown to perform surprisingly well, often beyond the small-λ and large-ω_c regimes.^29,30 Nevertheless, for inappropriate regions of parameter space, severe problems can arise, such as violation of positivity in the reduced density matrix.³ 

We wish to consider the performance of the Redfield-FM method for a strictly Markovian version of Redfield theory, i.e., with a rate tensor R = R(t → ∞). In this limit, the time integrals become Fourier–Laplace transforms, such that the Redfield tensor elements can be expressed algebraically in terms of the spectral density J(ω) evaluated at energy differences ħ ω_ij ≡ (E_i − E_j).^10,16 More specifically, we are interested in the dephasing terms of the Redfield tensor, which in general contain an elastic contribution

At low frequencies, the Bose–Einstein distribution, n_BE(ω) ∼ kT/ω, such that for a spectral density of the form J(ω) ∼ ω^s, we find

For all “super-Ohmic” spectral densities with s > 1, this elastic contribution to the dephasing rate vanishes. However, for an Ohmic spectral density with s = 1, there is a pure dephasing rate which vanishes only at T = 0. This contribution to the dephasing rate in the system’s eigenbasis can significantly affect both the population and coherence dynamics in the original basis of the problem.

We now return to the idea of a frozen modes variant of Markovian Redfield theory. Consider specifically an Ohmic spectral density with any non-zero splitting frequency ω^∗. After partitioning, the fast spectral density has the low-frequency behavior J_fast(ω) ∼ ω^s, with s > 1, which yields no elastic contribution to the dephasing rate. For this reason, a frozen modes version of Markovian Redfield theory does not reduce to the Redfield limit until the singular point ω^∗ = 0. Instead, as ω^∗ → 0, the result approaches a Redfield result which neglects the Ohmic pure dephasing rate. We emphasize that the time-dependent variants of Redfield theory are not significantly affected by this problem until very long times, and that all methods are only affected for strictly Ohmic spectral densities.

We propose a very simple solution to this pathological behavior in Markovian Redfield theory, by modifying the fast spectral density via

where W(ω, ϵ) is a rectangular window function centered at the origin with width ϵ, and ϵ should be chosen very small. In this way, the “fast” part of the bath will always produce a pure dephasing rate for arbitrary splitting frequency ω^∗. Thus, the Markovian Redfield-FM dynamics will smoothly interpolate towards the standard Markovian Redfield result as ω^∗ approaches zero. In Fig. 7, we compare the results of standard Markovian Redfield, Markovian Redfield-FM without this dephasing correction, and Markovian Redfield-FM with the correction. Results are presented for the model excitonic dimer discussed by Ishizaki and Fleming.⁴ Importantly, we find that this correction typically improves the results of the Markovian Redfield-FM variant, quite significantly in cases of strong system-bath coupling.

Here, we relax the mode-freezing approximation by deriving a fully hybrid method that separates the complete system into a slow part consisting of the low-frequency component of the bath, and a rapidly evolving part that includes both the electronic system and the high-frequency portion of the phonon bath. In this hybrid scheme, the slow part is treated quasi-classically, while the fast part is treated at the level of Redfield theory. Overall, the fast (slow) component of the system evolves in the mean field of the slow (fast) one.

Other hybrid approaches that combine classical and quantum dynamics include the self-consistent hybrid method of Wang and co-workers,^14,31 which yields numerically exact dynamics, and the approximate hybrid-NIBA approach of Refs. 8 and 9. In the former method, ω^∗, which is the energy scale that determines the splitting of the bath into slow and fast parts, is strictly a convergence parameter. In the latter, ω^∗ is an empirically determined adjustable parameter. As an approximate method, the hybrid-Redfield scheme derived here is akin to the hybrid-NIBA method. For a more detailed discussion of the hybrid RDM method, we refer the reader to Ref. 9.

As in the Redfield-FM approach, we make the approximation that Γ(t) ≈ ρ_slowρ_sys+fast(t), where ρ_sys+fast(t) is the density matrix for the system and fast bath degrees of freedom and ρ_slow(t) is the density matrix for the slow bath degrees of freedom. The system and fast bath modes obey the following effective Liouville equation:

where

and λ^cl is a dynamically fluctuating bias, $λ c l ( t ) = ∑ k ∈ slow c k Q k ( t )$⁠.

A classical treatment of the reservoir leads to the following equations of motion:

and

Employing the Ehrenfest approach demands that each part of the system evolves in the mean field of the other. For the quantum portion, the classical mean field consists of the time-dependent contribution to the bias, λ^cl(t). For the classical portion, the force term $− c k σ ̃ z ( t ) =− c k Tr sys+fast [ σ z ρ sys+fast ( t ) ]$ in the equations of motion embodies the mean-field “back-reaction.” This force term moves the classical oscillators from the ground state minima to the displaced minima associated with the excited state.

Using Eqs. (5) and (C2), the time-local Redfield equation takes the form

where $σ z ( t ) = U 0 † ( t ) σ z U 0 ( t )$⁠, $U 0 ( t ) = exp + [ − i ∫ 0 t d τ H sys ″ ( τ ) ]$⁠, and $H sys ″ ( t ) = [ ε + λ c l ( t ) ] σ z +Δ σ x + ∑ k ∈ slow [ P k 2 + ω k 2 Q k 2 ] /2$⁠. The bath correlation, as in the case of the Redfield-FM method, takes the following form:

To obtain the results shown in Fig. 3, trajectories corresponding to a set of initial conditions sampled from the Wigner distribution³² were calculated via a second-order Runge–Kutta scheme, using a step size of δt = 0.01Δ⁻¹. As required by the Runge–Kutta procedure, $σ ̃ z ( t )$ was kept constant during the evolution of the bath while λ^cl(t) was kept constant during the evolution of the system. Explicitly, over a half-time step, the equations for the bath become

and

where

In the hybrid-NIBA method of Ref. 9, the zeroth-order propagator necessary to evolve the perturbation in the interaction picture, $U 0 ( t ) =exp [ − i ∫ 0 ∞ d τ H sys ″ ( τ ) ]$⁠, was simple to calculate since $H sys ″$ was diagonal. In contrast, $H sys ″ ( t )$ for hybrid-Redfield contains off-diagonal elements. Within the Runge–Kutta scheme, this obstacle is easy to overcome, though it requires diagonalization of the time dependent $H sys ″ ( t )$ at every time step. Because numerical diagonalization at every time step is necessary for systems with more than two degrees of freedom, this can become computationally expensive for sufficiently large systems.