Using the Haken–Strobl–Reineker Model to Determine the Temperature Dependence of the Diffusion Coefficient

Stochastic quantum Liouville equations (SQLE) are widely used to model energy and charge dynamics in molecular systems. The Haken–Strobl–Reineker (HSR) SQLE is a particular paradigm in which the dynamical noise that destroys quantum coherences arises from a white noise (i.e., constant-frequency) spectrum. A system subject to the HSR SQLE thus evolves to its “high-temperature” limit, whereby all the eigenstates are equally populated. This result would seem to imply that the predictions of the HSR model, e.g., the temperature dependence of the diffusion coefficient, have no validity for temperatures lower than the particle bandwidth. The purpose of this paper is to show that this assumption is incorrect for translationally invariant systems. In particular, provided that the diffusion coefficient is determined via the mean-squared-displacement, considerations about detailed-balance are irrelevant. Consequently, the high-temperature prediction for the temperature dependence of the diffusion coefficient may be extrapolated to lower temperatures, provided that the bath remains classical. Thus, for diagonal dynamical disorder the long-time diffusion coefficient, D∞(T) = c1/T, while for both diagonal and off-diagonal disorder, D∞(T) = c1/T + c2T, where c2 ≪ c1. An appendix discusses an alternative interpretation from the HSR model of the “quantum to classical” dynamics transition, whereby the dynamics is described as stochastically punctuated coherent motion.


Introduction
Coherent exciton dynamics in static, ordered molecular systems was described by Merrifield 1 in 1958.Assuming an exciton created at time t = 0 on a monomer, say n = 0, he showed that the subsequent wavefunction is Ψ n (t) = J n (2βt), where J n is the nth order Bessel function of the firstkind and β is the intermonomer exciton transfer integral.The wavefunction (illustrated in Fig. 2) spreads ballistically with a constant speed and a mean-squared-displacement (MSD) increasing quadratically with time.
As Merrifield observed, 1 however, dynamics on a static, ordered system is an idealization.
Various physical processes, e.g., exciton-phonon coupling, and static and dynamic disorder destroy the coherent motion, eventually causing incoherent (or diffusive) motion where the MSD increases linearly with time.[4][5][6] The purpose of this paper is to expand on one particular rich area of investigation, namely the role of thermally-induced noise in destroying quantum coherences.A notable paradigm in this subject is the so-called HSR stochastic quantum Liouville equation (SQLE), developed and investigated by Haken, Strobl and Reineker. 2,7This equation was developed from the underlying time-dependent Schrödinger equation (TDSE) assuming that the dynamical fluctuations obey a white-noise spectrum, i.e., a constant power spectrum (or an Ohmic spectral function).Many important results have been derived from this model. 2,8,9In particular, it describes the 'quantum to classical' transition, in which exciton dynamics exhibits a crossover from ballistic to diffusive as a result of the noise destroying the coherent motion.
As stated, the HSR model assumes a white-noise spectrum, which implies that quantum tran-sitions can occur between any pair of system energy eigenstates.This is turn implies that a system subject to the HSR SQLE evolves to its 'high-temperature' limit, whereby all the eigenstates are equally populated.This result would seem to suggest that the predictions of the HSR model, e.g., the temperature dependence of the diffusion coefficient, have no validity for temperatures lower than the particle bandwidth.The purpose of this paper is to show that this assumption is incorrect for translationally invariant systems.In particular, provided that the diffusion coefficient is determined via the mean-squared-displacement, considerations about detailed-balance are irrelevant.
Consequently, the high-temperature prediction for the temperature dependence of the diffusion coefficient may be extrapolated to lower temperatures, provided that the bath remains classical.
This key result will be proved in Section 2. Since the HSR QSLE predictions for the diffusion coefficient in translationally invariant systems are valid for temperatures lower than the particle bandwidth, the predictions of the underlying stochastic TDSE are also equally valid.We use this realization to reinterpret the 'quantum to classical' transition as stochastically punctuated coherent motion.This is described in Appendix A.
Unlike the use of the MSD, Appendix B shows that the velocity autocorrelation function cannot be used to extrapolate the HSR predictions to temperatures lower than the bandwidth.Finally, Appendix C derives an expression for the temperature-dependence of the dephasing rate and Appendix D contains some details of the computational techniques.

Model of Exciton Dynamics in Linear Molecular Systems
We formulate the problem in terms of Frenkel exciton dynamics in one-dimensional molecular systems, e.g., J-aggregates or conjugated polymers.However, the analysis applies equally to triplet excitons and charges.

The total Hamiltonian is
where ĤS , ĤSB and ĤB are the system, system-bath and bath Hamiltonians, respectively.ĤB is defined in Appendix C (eqn (47)).
The system Hamiltonian is defined as The ket |m⟩ represents an exciton on monomer m, denoting a 'site', where N is the number of sites.
α and β are the onsite potential and nearest-neighbor exciton transfer integral, respectively.
For a system with translational invariance, the eigenstates of ĤS are the Bloch states with eigenvalues E a = α + 2β cos k a , where k a = 2πa/N is the wavevector.The quantum numbers that label the eigenstates satisfy 1 ≤ a ≤ N. The particle bandwidth in one-dimension is 4|β |.
For a classical, harmonic bath, γ × ∝ σ 2 × ∝ k B T .More specifically, as shown in Appendix C, for linear system-bath coupling where E r × is the reorganization energy arising from the system-bath coupling and ω c is a highfrequency cut-off for the spectral function.

Determining the Temperature-Dependent Diffusion Coefficient
The one-dimensional thermal diffusion coefficient as a function of time is defined as where ⟨• • • ⟩ indicates a thermal average, x is the operator for the particle position, and ρ(t, T ) is the system's reduced density operator.In the long-time limit the asymptotic diffusion coefficient is where is the equilibrium density operator.
Evaluating the trace over the eigenstate basis of ĤS , Eqn (9) becomes where and p a (T ) is the Boltzmann factor.
The mean-squared-displacement of a particle prepared in an arbitrary state |ψ⟩ at t = 0 is defined as where ℓ is the intermonomer separation and ρ mm (t) are the diagonal elements of the system density matrix in the site basis, {|m⟩}.As described in Section 2.3, the density matrix, ρ mn (t), evolves according to an appropriate quantum Liouville equation with the initial condition ρ mn (0) = ⟨m|ψ(0)⟩⟨ψ(0)|n⟩.
Transforming to the eigenstate basis via a transformation matrix, S, Eqn (13) becomes where ρab is the density matrix in the eigenstate basis and S ma = ⟨m|a⟩.For a system with translational invariance, the transformation matrix elements are the Bloch factors Splitting the double sum in Eqn ( 14) over a and b into the separate sums of a = b and a ̸ = b, and using Eqn (15) we obtain where we have also used ∑ a ρaa = 1.The significance of this result is that it shows that for a translationally invariant system the diffusion coefficient of a particle prepared in the state |ψ⟩ depends on the evolution of eigenstate coherences and not directly on the evolution of eigenstate populations.Moreover, as will be shown in Section 2.3, for a translationally invariant system, eigenstate coherences are decoupled from eigenstate populations, and thus for such systems the diffusion coefficient is completely independent of eigenstate populations.This means that when evaluating Eqn (11) to determine ⟨D ∞ (T )⟩ enforcing detailed balance -or ensuring that eigenstate populations satisfy their thermal values -is unnecessary.Thus, ⟨D ∞ (T )⟩ only depends on temperature parametrically via the temperature-dependence of the dephasing rate, Eqn (7).
Finally, as shown in Section 2.3, the asymptotic diffusion coefficient is independent of the initial state.Denoting this asymptotic value by D ∞ (T ), setting D a = D ∞ (T ), ∀a, and using ∑ a p a = 1, Eqn (11) becomes ⟨D ∞ (T )⟩ = D ∞ (T ).Thus, our task now is to determine D ∞ (T ) for a given dephasing rate.This is achieved via the Haken-Strobl-Reineker model as described in the next section.

The Haken-Strobl-Reineker Model
Haken and Strobl showed that ensemble averages of observables determined via the stochastic TDSE may be evaluated by a SQLE. 2,7For the case of diagonal noise within the site basis, the SQLE reads Rotating to the eigenstate basis of ĤS via Eqn (15), the SQLE becomes for the populations, and for the coherences, where From Eqn (18) and Eqn (19) we note the following: 1.In the HSR model k ab = k ba = 2γ α /N, thus guaranteeing equal eigenstate populations in the long-time limit.In principle, one could impose detailed balance on the rates, thus ensuring that the populations equilibrate to their thermal values.However, as now shown, this is unnecessary.
2. The equations of motion for the populations and coherences are decoupled.This proves, by virtue of Eqn (16), that ⟨D ∞ (T )⟩ is independent of temperature-dependent populations and only depends on T parametrically via γ α (T ).
3. Eqn ( 18) is a single N-coupled equation, whereas Eqn ( 19) are (N − 1) × N-coupled equations.Thus, there are N × N-coupled equations in total to solve.A numerical solution is described in Appendix D.

Results
Haken, Reineker and co-workers derived an expression for D(t) for an initial δ -function source. 2,10r diagonal noise this is Reineker 2 also derived the asymptotic value, D ∞ , for arbitrary initial conditions.For diagonal noise These analytical predictions were confirmed by solving both the HSR SQLE and the stochastic TDSE numerically (as described in Appendix D).Fig. 1 shows the numerically evaluated D(t) for different initial conditions, namely a δ -function, Gaussian and particle-in-a-box sources.The δfunction result satisfies Eqn (20), while for all other sources, having a larger initial mean-squaredsize and thus a smaller initial mean-squared-speed, the diffusion coefficient increases more slowly with time.However, Eqn ( 21) is satisfied for all initial conditions.
For translationally invariant systems, the evolution of the populations and coherences are also The diffusion coefficient at short times for the Gaussian wavepackets is superlinear, as also reported in ref. 11 decoupled in the presence of off-diagonal noise.Thus, for both diagonal and off-diagonal noise, the 'high-temperature' limit determined from the HSR SQLE can be extrapolated to temperatures lower than the bandwidth.Reineker 2 showed that for arbitrary initial conditions 12 where γ α and γ β are the diagonal and off-diagonal dephasing rates, respectively.
Substituting the temperature-dependence of the dephasing rates given by Eqn (7), we thus have the following prediction for the temperature-dependence of the diffusion coefficient in the presence of white noise: where and Equation ( 23) is valid for translationally invariant systems subject to white noise for all temperatures -including temperatures lower than the particle bandwidth, 4|β | -provided that the bath remains classical.It remains valid if uniform long-range couplings are included, although the coefficients c 1 and c 2 are altered.

Conclusions
This paper has shown that for translationally invariant systems the predictions of the HSR model for the thermal diffusion coefficient can be extrapolated to temperatures lower than the particle bandwidth, provided that the bath remains classical.The proof relies on the observation that for such systems the mean-squared-displacement is independent of eigenstate populations.Consequently, considerations about detailed balance are irrelevant, and thus the diffusion coefficient depends on temperature only parametrically via the dephasing rates.When diagonal disorder dominates, D(T ) ∼ T −1 .This 'high-temperature' limit for D(T ) is a common prediction in many theories of charge and energy transport. 13anslationally invariant systems subject to white noise are an idealization of more realistic systems, where static disorder, correlated (non-Markovian) noise and electron-phonon interactions causing polaron formation are all important processes that will modify the predictions of this paper.
For such systems, numerically solutions of the TDSE or QLE should explicitly ensure that detailed balance and stationarity are maintained during the system's evolution (e.g., the Redfield quantum chical equations for open systems, 19 and stochastic Liouville equation methods 20 ).The recently proposed MASH surface-hopping schemes 21,22 are also promising techniques for such simulations.
We conclude by noting that in the presence of static, diagonal disorder the scaling of D(T ) with temperature deviates from the HSR prediction. 8,9,23In particular, for γ ∼ T < β , D(T ) ∝ T .This increase of the diffusion coefficient with temperature at low temepatures, sometimes known as environment-assisted quantum transport, 8 is a consequence of the thermal fluctuations destroying the Anderson localization arising from coherent superposition of the particle wavefunction.

Appendix A. The Quantum to Classical Transition
In free space the one-dimensional TDSE reads where m is the particle mass.As is well known, 15,24 a particle prepared in a Gaussian wavepacket at t = 0, i.e., undergoes a ballistic (or coherent) spread.Thus, its MSD, σ 2 x (t), increases as .
In contrast, the diffusion equation describes the evolution of the density of an ensemble of particles, ρ(x,t).In free space the one-dimensional diffusion equation reads Although this equation has precisely the same mathematical form as Eqn (26), the absence of imaginary i means that its physical predictions are quite different.In particular, an initially localized distribution of particles spreads diffusively, i.e., where the MSD increases as The HSR equation describes how white noise causes a quantum particle's coherent motion to become incoherent.As shown by Reineker, 2 assuming a δ -function source on a lattice, i.e., (as predicted by the TDSE on a lattice 1 ) where J n is the nth-order Bessel function of the first kind.
In contrast, for t ≫ 2γ, (as predicted by the diffusion equation) where D = σ 2 x (t)/2t = β 2 /γ α for diagonal noise.In general for a δ -function source (where henceforth in this Appendix, we set h = ℓ = 1) 2,25 In this interpretation, noise causes a crossover from coherent to incoherent transport.A different interpretation of the 'quantum to classical' transition, however, is afforded by the Lindblad formulation of the quantum Liouville equation.As we now show, in this interpretation the particle's dynamics can be viewed as stochastically punctuated coherent motion.
2. An ensemble of quantum trajectories reproduces the observables obtained via a QLE if each trajectory undergoes non-unitary evolution via an effective non-hermitian Hamiltonian, 26 defined by 3. The simulation of a quantum trajectory then proceeds as follows: (a) Given |Ψ(t)⟩ at a time t, compute where δt is the time step.(See Appendix D for details.)(b) Determine δ p, defined via ⟨Ψ trial |Ψ trial ⟩ = 1 − δ p.
(c) Then, (i) with a probability (1 − δ p) define In this case the wavefunction evolves according to Ĥeff .
(d) If a quantum jump occurs in 3(c)(ii), the site m is chosen with a probability Thus, the effect of the Lindblad jump operator in step is 3(c)(ii) to collapse the coherently evolving wavefunction onto the site m, after which it then again expands coherently until the next collapse.This process is illustated in Fig. 2: At time t = 0 the particle is created on site n 0 = 0.
Its wavefunction then evolves coherently until the first quantum jump at t = t 1 .For t ≤ t 1 , Ψ n (t) = J n (2βt), where J n is the nth order Bessel function. 1At t = t 1 the wavefunction collapses onto site n = n 1 , whereupon it again evolves coherently until a time t 2 when it again collapses onto site n = n 2 .
Figure 2: An illustration of the evolution of a particle wavefunction subject to the TDSE with stochastic wavefunction collapses caused by white noise.The particle probability, P n (t) = J 2 n (2βt).The cumulative effect of the stochastic wavefunction collapses is to cause particle diffusion, where D = β 2 /γ α , in exact agreement with the HSR model.γ α = 0.1β /h.Time is in units of h/β .As these collapses are stochastic in space and time, their cumulative effective is to cause particle diffusion.We can determine the diffusion coefficient as follows.For a diffusive process the MSD is defined as where N(t) = t/⟨τ⟩ is the number of jumps in a time t and ⟨τ⟩ = 1/Γ = 1/2γ α is the average time interval between jumps.⟨ℓ 2 ⟩ is the mean-squared jump size, determined by ⟨ℓ where v is the particle's coherent speed.For a δ -function source on a lattice, v = √ 2β .As determined by simulating the TDSE with Lindblad jump operators, for the dynamics described here, Hence, ⟨ℓ 2 ⟩ = β 2 /γ 2 α and therefore D = N(t)⟨ℓ 2 ⟩/2t = β 2 /γ α , thus -on average -rederiving the prediction of the HSR quantum Liouville equation.Ensemble averaging over many particle trajectories will give precisely the same observables (e.g., average particle density) as the stochastic quantum Liouville approach, but the interpretation of the dynamics is different, namely stochastically punctuated coherent dynamics.

Appendix B. Using the Velocity Autocorrelation Function
This is easily accomplished via the Short Iterative Lanczos Propagator method, 24 whereby a Krylov space of N vectors is generated via the Lanczos method.Diagonalizing the tridiagonal Hamiltonian, H Lanczos , within the Krylov space yields where S is the is the matrix whose columns are the eigenvectors of H Lanczos , the diagonal elements of D are its eigenvalues, and Ψ(t) is the first vector in the Krylov space.
The expectation value of the velocity autocorrelation function is where |Φ⟩ = v|Ψ⟩ and the velocity operator on a lattice is