Detailed balance in mixed quantum–classical mapping approaches

The violation of detailed balance poses a serious problem for the majority of current quasiclassical methods for simulating nonadiabatic dynamics. In order to analyze the severity of the problem, we predict the long-time limits of the electronic populations according to various quasiclassical mapping approaches, by applying arguments from classical ergodic theory. Our analysis conﬁrms that regions of the mapping space that correspond to negative populations, which most mapping approaches introduce in order to go beyond the Ehrenfest approximation, pose the most serious issue for reproducing the correct thermalization behaviour. This is because inverted potentials, which arise from negative electronic populations entering into the nuclear force, can result in trajectories unphysically accelerating oﬀ to inﬁnity. The recently developed mapping approach to surface hopping (MASH) provides a simple way of avoiding inverted potentials, while retaining an accurate description of the dynamics. We prove that MASH, unlike any other quasiclassical approach, is guaranteed to describe the exact thermalization behaviour of all quantum–classical systems, conﬁrming it as one of the most promising methods for simulating nonadiabatic dynamics in real condensed-phase systems.


I. INTRODUCTION
Quantum nonadiabatic effects play a crucial role in the description of many relevant physical processes, including radiationless decay, 1 energy and charge transfer 2 and coherence control in quantum gates. 3Unfortunately the computational cost of brute-force quantum calculations scales exponentially with the number of system degrees of freedom, making it impossible to simulate the vast majority of realistic systems in this way.An increased interest in the development of quasiclassical methods has therefore arisen, aimed at mapping nonadiabatic dynamics onto effective quasiclassical models that can be simulated as efficiently as classical systems.However, improvements in efficiency are often achieved at the expense of accuracy.Although rigorous quantum-classical dynamics can be defined in terms of coupled trajectories, 4,5 we wish to reduce the numerical complexity by evolving trajectories independently.These quasiclassical methods are often accurate at short times, but can suffer from significant long-time errors. 6In particular, if detailed balance is violated, the final populations will not reflect the correct thermal distribution.
Ehrenfest dynamics is the prototypical quasiclassical nonadiabatic theory that propagates the nuclear dynamics on a mean-field potential energy surface corresponding to a dynamical superposition of electronic states. 7,80][11][12][13][14][15] This is in part due to the fact that the a) These authors contributed equally; Present Address: Institute of Physics, University of Freiburg, Freiburg, Germany b) These authors contributed equally; Present Address: Max Planck Institute for the Structure and Dynamics of Matter, Hamburg, Germany c) Electronic mail: jeremy.richardson@phys.chem.ethz.chmethod drastically violates detailed balance in thermal equilibrium. 6,16n order to improve upon the accuracy of Ehrenfest dynamics, a wide array of mean-field mapping approaches have been developed.The main idea is to construct a quasiclassical continuous phase space for the electrons, so that the electronic and nuclear degrees of freedom can be treated on an equal footing.8][19][20][21][22][23] Further improvements in the accuracy of quasiclassical approaches can be achieved by treating the identity operator exactly within the mapping, by representing it by the number one. 11,24,25A similar improvement is obtained more naturally by departing from the MMST formalism and instead describing the electronic degrees of freedom in terms of a set of spherical spin coordinates, known as spin mapping. 12,26While all these developments often lead to improved dynamics over Ehrenfest theory, they all clearly still violate detailed balance in general.
One important difference between mapping-based approaches and Ehrenfest dynamics is that mapping-space regions corresponding to negative electronic populations are present in the former.Although mapping approaches are typically more accurate than Ehrenfest simulations, these problematic regions can pose a serious issue for the dynamics, as the resulting inverted potentials may give rise to an unphysical nuclear force that can cause trajectories to accelerate off to infinity. 27n addition to these problems, there is a conceptual difficulty with using mapping approaches in strongly asymmetric systems in which one initializes the simulation in an excited state and expects the dynamics to thermalize to the ground adiabatic state.Here, the short-time dynamics must be described by a multi-state nonadiabatic system, whereas the long-time dynamics is an effective one-state problem.A method that attempts to tackle this issue head-on is the ellipsoid-mapping approach, 28 which replaces the spherical spin-mapping phase space with an anisotropic ellipsoid geometry.The shape and orientation of the ellipsoid is dynamically adjusted so as to best represent the effective structure of the local electronic Hamiltonian.This is achieved in such a way that the dynamics rigorously obey detailed balance by construction.The approach is however only applicable for computing thermal correlation functions and not able treat the thermalization of systems initialized out of equilibrium, as we wish to do in this paper.
Another approach designed to go beyond meanfield methods is the symmetric quasiclassical (SQC) approach. 29Although the dynamics of SQC are identical to those of the mean-field methods introduced previously, the electronic populations are instead measured by windowing the electronic phase space in such a way that the obtained values are guaranteed to lie in the physical range between zero and one.An advantage of the method compared to other quasiclassical theories is that SQC is known to obey detailed balance in the limit of zero electron-nuclear coupling. 30However, although it is often an improvement over the original MMST methods, SQC may still lead to inaccurate long-time dynamics in the more general coupled regime, 31 and is typically of a similar level of accuracy to the spin-mapping methods. 12,26his is because although SQC avoids measuring negative populations, its trajectories still suffer from their contribution to the nuclear force.
Fewest switches surface hopping (FSSH) 32 ensures that the nuclear force is always physical by propagating the nuclei on a single adiabatic surface at any given time.In order to describe nonadiabatic transitions, the active surface is changed stochastically with a probability that mimics the dynamics of the underlying electronic wavefunction.After a transition, the nuclear momenta are rescaled in the direction of the nonadiabatic coupling vector, consistent with earlier semiclassical scattering theories. 33,34Remarkably, FSSH has been shown to approximately fulfill detailed balance, even in cases when the electronic and nuclear degrees of freedom are strongly coupled (although rigorously so only in the limit of strong nonadiabatic coupling). 35Although much is known about the approximations underlying surface hopping (see in particular the work of Subotnik and Kapral), 36,37 the main problem of FSSH is that the approach still lacks a complete formal derivation.This means that there is still some disagreement over the correct way of performing certain aspects of the method that could not be derived from first principles, such as the treatment of frustrated hops.Additionally, the stochastic nature of FSSH means that the electronic wavefunction and the current active surface can become inconsistent during the dynamics, leading to the so-called 'inconsistency error' that is known to significantly degrade the accuracy of FSSH observables. 38t seems that the ultimate quasiclassical approach for accurately describing both the short-and long-time dynamics in nonadiabatic systems is one that combines the best features of SQC and FSSH.The mapping approach to surface hopping (MASH) is a newly developed quasiclassical method that achieves just that. 38MASH consistently windows both the observables and the nuclear force, eliminating the possibility of obtaining negative populations in either.This results in MASH having deterministic dynamics, for which the active surface and the electronic degrees of freedom remain consistent throughout the entire time evolution.The approach is also completely derivable from first principles, leading to a unique prescription for frustrated hops that ensures that the exact short-time dynamics is correctly reproduced.Finally, MASH has been tested on a range of commonly-used condensed-phase model systems, where it appears to describe the correct long-time thermalization behaviour.
In this paper, we utilize classical ergodic theory to predict the long-time thermalization behaviour of a wide range of quasiclassical approaches, including mean-field mapping methods, SQC and MASH.By comparing our predictions to those of the expected quantum-classical outcome, we provide a simple and rigorous procedure for benchmarking the long-time dynamics.We apply our theory to spin-boson models in challenging parameter regimes and also to an anharmonic model, which together illustrate why negative populations pose a serious problem for reproducing the correct thermalization behaviour for the majority of quasiclassical approaches.Despite this, we show that MASH is guaranteed to exactly reproduce the correct thermal populations in the long-time limit for all quantum-classical systems, confirming its potential for being one of the best methods for accurately simulating nonadiabatic dynamics in real condensed-phase systems.

II. THEORY
In the present work, we focus on the dynamics of a quantum subsystem consisting of N states, k⟩, coupled to a classical environment, in this case a heat bath.The general form of the Hamiltonian is where (x, p) are the multidimensional classical phasespace variables of the environment.U (x) and V (x) are state-independent and state-dependent potentials respectively; the latter is a traceless N × N matrix in the diabatic basis, { k⟩} N k=1 .Let us remark that it is valid to treat the bath classically provided that the energy associated with the thermal fluctuations of the environmental modes is large compared to their zero-point energies, although not necessarily large compared to the energy scales of the quantum subsystem.In the framework of atomistic systems in the condensed phase, the subsystem and the environment typically refer to electronic and nuclear degrees of freedom respectively, although this particular situation is not strictly required in the following.
The aim of quasiclassical approaches is to accurately predict the evolution of quantum correlation functions, with a computational cost comparable to classical simulations.Approximate quasiclassical dynamics associated with Eq. ( 1) can be derived by first taking the partial Wigner transform of the quantum Hamiltonian and then taking the classical limit of the bath. 39,40The quantum subsystem is mapped onto a continuous phase-space representation, via Cartesian mapping variables, X, P ∈ R N . 9,41This leads to a representation of the full system in terms of the phase-space points, Γ = {X, P, x, p}.The Hamiltonian [Eq.( 1)] is then mapped onto a phase-space function, H(Γ) = p 2 2m + U (x) + V(x, X, P ), from which the dynamics are obtained.
Quasiclassical approaches derived in this way result in equations of motion of the form Here we have assumed that the matrix elements ⟨k V (x) k ′ ⟩ are real.The first two lines are equivalent to the Schrödinger equation for the real and imaginary parts of the electronic wavefunction.The last two lines are Newton's equations of motion for the bath modes.The expressions for the state-dependent potential, V(x, X, P ) and nuclear force, F j (x, X, P ), depend on the specific method considered.Apart from the total energy, H(Γ), these equations of motion additionally conserve the norm of the Cartesian mapping variables, r = 1 2 ∑ N k=1 (X 2 k +P 2 k ).From these dynamics, quantum correlation functions, Tr[ρ 0 Â B(t)], are then approximated by where the observable operators are mapped onto phasespace functions [i.e., Â(x, p) ↦ A(Γ)], the specifics of which are method dependent.Note that in some approaches, the quasiclassical representations of the initial and time-evolved operators can differ, even if quantum mechanically the operators are identical (due to different projections of electronic operators 11,22 or due to the inclusion of weighting factors into correlation functions 38 ).
A well known limitation of the majority of quasiclassical methods is their inability to obey detailed balance and recover the correct long-time populations. 6,16,28,42This issue results in dynamical correlation functions relaxing to incorrect long-time limits.As we will see, the way in which the correlation function operators and the statedependent force are represented by phase-space functions can lead to significant differences in the thermalization properties of different methods.

A. Mean-field approaches
For mean-field approaches, the state-dependent potential is approximated by its expectation value with respect to the Cartesian mapping variables.This gives rise to dynamics that correspond to the nuclei being propagated on a weighted average of the potential energy surfaces Let us remark that Eqs. ( 2) are simply Hamilton's equations of motion given the mean-field form of the statedependent potential [Eq.(4a)].
The prototypical mean-field approach is Ehrenfest dynamics, for which both the nuclear force and the electronic observables are given by expectation values associated with the normalized time-evolved electronic wavefunction (here expressed in terms of mapping variables).Improved mean-field methods have been obtained by mapping the electronic subsystem onto a fictitious system containing continuous degrees of freedom, (X, P ), such that the quantum and nuclear subsystems can be treated on an equal footing.Whereas in Ehrenfest theory, X and P simply correspond to the real and imaginary parts of the electronic wavefunction, in the mapping approaches, they span a space in which representations of the wavefunction are constructed.While the mean-field expression for the state-dependent nuclear force remains the same [Eq.(4b)], the phase-space representations of the observable operators generally depend on the mapping space used.
In the case of the Meyer-Miller-Stock-Thoss (MMST) mapping, 17,18 electronic states are mapped onto the single-excitation subspace of a set of N harmonic oscillators.This leads to (at least) two possible representations of the electronic operators in terms of phase-space functions. 11,22In the Wigner approach, the phase-space functions are obtained from the Wigner transform of the associated operators in the mapped harmonic-oscillator system, analogous to the representation of the potential energy operator used in Eq. (4a).In the singly-excited oscillator (SEO) approach, 18 projection operators onto the physically relevant single-excitation subspace are added, which results in an additional factor of φ(r) = 16e −2r appearing from the Wigner transform.At least one factor of φ(r) is required for the integral over mapping variables to converge.These exponential factors are normally incorporated into the definition of the initial density, ρ 0 (Γ), as they define the distribution from which the Cartesian mapping variables are initially sampled.We will call this contribution ρ 0,s .
Both of these approaches lead to the representation of the identity operator, I s , being a function of the mapping variable norm, r.It was however previously observed that representing the identity operator with the number one (which can be thought of as the exact mapping for this operator) can lead to more accurate results. 11,24,25e refer to MMST methods that represent the identity operator in this way as 'unity approaches'.
Another variant of MMST mapping utilizes so-called focused initial conditions when the dynamics are initialized in an electronic population, 21 by initially sampling the mapping variables from the regions that satisfy A(Γ) = 1.These thus have delta functions for their initial r-distributions.All Ehrenfest trajectories have r = 1 due to the absence of electronic zero-point energy, whereas in focused MMST methods, which utilize the energy levels of a harmonic oscillator, the populated state contributes 3 2 and the N − 1 unpopulated states contribute 1 2 each, leading to r = 3 2 + 1 2 (N − 1) = 1 2 (N + 2).For spin mapping, 12,26 the electronic states are described by a set of spin variables, which results in the following advantageous features over MMST mapping.First, the associated mapping space is isomorphic with the original quantum subspace, such that the use of projection operators is no longer required and the representation of observable operators is therefore unique.Second, the correspondence between the identity and the number one arises naturally from the underlying theory and does not need to be imposed, as in the case of the unity approaches.The algorithm is practically identical to the MMST-focused method except that the hypersphere radius is r = √ N + 1. Spin mapping in its original formulation is a linearized quasiclassical method. 12,26Spin-PLDM is a partially linearized extension 13,43 developed within the framework of the partially linearized density matrix (PLDM) approach. 44,45The equations of motion are similar to those of the fully-linearized mapping methods except that two sets of mapping variables are employed to describe the forward and backward propagation separately.
Table I summarizes the different ways in which observable operators are represented by the methods studied in this work.We remark that double-SEO is also known as the linearized semiclassical initial value representation (LSC-IVR), 19,20 while single-SEO is often referred to as the Poisson bracket mapping equation (PBME). 22,46ost cases are easily obtained by rewriting the expressions from the original papers 11,12,26 in terms of r.The more obscure spin-PLDM expression for ρ 0,s given in Table I has first been derived in Ref. 47 by integrating out all other spin degrees of freedom apart from the centroid.Note that this particular expression is only valid for twostate systems.
One of the major problems of most mean-field approaches is the presence of negative populations, which can lead to both unphysical values for the long-time limit of correlation functions and unphysical dynamics corresponding to propagating the nuclei on inverted potentials.In the following, we introduce two approaches that are designed to alleviate the problem of negative populations, namely the symmetric quasiclassical theory (SQC) and the mapping approach to surface hopping (MASH).

B. Symmetric quasiclassical windowing (SQC)
The symmetric quasiclassical approach (SQC) attempts to improve upon the standard MMST mapping approaches by measuring the electronic populations with 'windows' that guarantee that population observables always physically lie between zero and one.These windows are usually defined in terms of the action-angle variables (n k ≥ 0 and q k respectively), which are related to the Cartesian mapping variables introduced in Sec.II as follows: 48 The most successful choice for the windowing functions are the so-called 'triangular windows', which for state k take the form 49 where h(⋅) denotes the Heaviside step function, the zeropoint energy parameter is set to γ = 1 3, and the angle variables, q k , are uniformly sampled independently from the interval [0, 2π).The success of this windowing scheme is due to the fact that the windows associated with different electronic states touch, allowing the approach to correctly describe population transfer even between weakly-coupled states. 50However, because the windowing functions given by Eq. ( 6) do not fill the entire mapping space, the sum of the populations is not a constant of motion.Therefore, all SQC observables must be renormalized to obtain a total population of one whenever they are measured.While this scheme guarantees that the population observables always return physical values, the approach still suffers from the effects of inverted potentials, because the dynamics evolve via the same Hamiltonian as for the mean-field methods.This means that while it has been shown that SQC correctly thermalizes in the limit of weak system-bath coupling, 30 it is expected that this is in general not true in other parameter regimes, particularly for strongly-coupled and anharmonic systems.One suggestion to alleviate the problem of inverted potentials is to introduce a state-dependent zero-point energy parameter, γ k , into the definition of the state-dependent nuclear force [Eq.(4b)], and adjust their values independently for each trajectory so that the contribution to the force from each electronic population is correctly reproduced. 51On TABLE I.A comparison between the different quasiclassical methods studied in this work.For each method, the initial mapping-variable distribution, ρ0,s and the identity operator representation, Is are given.Here, φ(r) = 16e −2r .Derivations of the radial distributions and representations of the identity for different mean-field approaches can be found in Refs.11 and 12. Additionally, the accuracy of the thermalization behaviour of each method in the high-temperature (β → 0), low-temperature (β → ∞) and weak electron-nuclear coupling (α → 0) limits is indicated, valid for any two-state system.We mark with ticks all entries that agree with the correct quantum-classical result, along with any erroneous multiplicative factors for the β → 0 and β → ∞ limits.For example, Ehrenfest predicts the long-time adiabatic population difference a factor of 3 too small when β → 0 but is correct when β → ∞.Further details on the representation of the operators in different methods are given in Sec.IV.

Method
r a Formulas and results for the case of N = 2 only b Mean-field MMST mapping method the one hand, this strategy can reduce the likelihood of trajectories unphysically accelerating off to infinity (at least at short times).On the other hand, this means that the values of γ k used for measuring the observables and performing the dynamics are inconsistent.

C. Mapping approach to surface hopping (MASH)
The mapping approach to surface hopping (MASH) 38 is a new nonadiabatic trajectory approach that offers the best of both worlds between FSSH and quasiclassical mapping dynamics.MASH can be thought of as going beyond SQC by windowing both the observables and the nuclear force, thus solving the problem of inverted potentials completely.Windowing forces in this way introduces hops similar to those of FSSH, although the inconsistency error of FSSH is not present in MASH.These appealing features make MASH one of the most promising approaches for performing accurate nonadiabatic dynamics within realistic ab initio simulations of molecules at a relatively low computational cost.
MASH was originally derived for systems involving two electronic states using the spin-mapping variables in the adiabatic basis, given by where ± refer to the upper and lower adiabatic states.Like FSSH, MASH works within the kinematic picture, where the nuclear momentum retains the intuitive meaning of mass times velocity, unlike the canonical momentum associated with the Hamiltonian in the adiabatic representation. 52Thus the MASH expression for the energy takes the same form as in Eq. ( 1), except that the state-dependent potential is now given in the adiabatic basis.In terms of these mapping variables, the MASH representation of the state-dependent potential and the nuclear force operator are given by 38 where sgn(⋅) returns the sign of its argument, the adiabatic potential energy surfaces are is the nonadiabatic coupling vector 53 and δ(⋅) is the Dirac delta function.From Eq. (8a) it can be seen that the adiabatic windows used in MASH correspond to the spin-hemispheres.These windows, like the triangular windows of SQC, touch at the equator, but unlike SQC windows, they fill the whole mapping space, which is necessary in MASH to guarantee that the force is well defined throughout the entire time-propagation.The last term on the right-hand side of Eq. (8b) ensures energy conservation by applying an impulse at the equator.As MASH is guaranteed to exactly reproduce the short-time dynamics of the quantum-classical Liouville equation, 4 this term therefore constitutes a unique and correct prescription for the momentum rescaling and frustrated hops, which is a feature that is lacking in almost all previous surface-hopping approaches.In practice the algorithm is simply a deterministic version of FSSH with momentum rescalings at each attempted hop.Finally, the initial mapping-variable distribution, ρ 0,s = W AB (S ad ), is determined so that the MASH windows exactly reproduce the Rabi oscillations of a bare electronic system. 38ecause of its use of the kinematic picture, the equations of motion associated with MASH are not obviously generated by a Hamiltonian.This means that H(Γ) with the MASH representation of the potential [Eq.(8a)] is technically a conserved energy rather than the generator of the quasiclassical dynamics.This raises the question of whether applying classical ergodic theory to predict the long-time limit of the MASH correlation functions is justified.Despite this, the method does possess many of the features associated with Hamiltonian dynamics, such as the conservation of phase-space volume, suggesting that the use of classical ergodic theory may still be applicable.In Sec.III B, we will show numerically that the long-time limits of MASH correlation functions do agree with the predictions of classical ergodic theory, which justifies our use of it for MASH.

D. Ergodic Theory
Quasiclassical approaches that are both deterministic and based on independent trajectories can be analyzed within the framework of classical ergodic theory.In this section, we expand this idea to derive a general expression for the long-time limits of correlation functions, Eq. ( 12).As previously discussed, the equations of motion [Eq.( 2)] conserve the norm of the Cartesian mapping variables, We cannot therefore simply apply classical ergodic theory directly to our problem.In fact, according to the so-called ergodic hierarchy, 54 a necessary requirement for the dynamics to exhibit mixing behavior on a given manifold of the phase space is the condition of ergodicity, i.e., only the Hamiltonian, H(Γ), is allowed to be conserved on that manifold.Thus, in order to isolate the additional conserved variable, r, we define a set of hyperspherical coordinates for the mapping variables, (X, P ) ↦ (r, Ω), where Ω denotes the solid angle.By introducing Γ = (Ω, x, p), we can rewrite Eq. (3) as To evaluate the long-time limit of Eq. ( 9), we define ρ (r) t,0 ( Γ′ Γ) as the conditional probability of occupying the state Γ′ at time t, given that the dynamics propagate on the hypersurface of constant r and that the system was initialized in state Γ.Equation (9b) is then expanded as We will assume that the dynamics fulfill the strong mixing condition 54,55 lim t→∞ ρ (r) where Z(r) = ∫ d Γ e −βH(r, Γ) .The equilibrium canonical distribution in Eq. ( 11) is the long-time distribution expected for a classical ergodic system with finite temporal correlations. 56The inverse temperature, β, in Eq. ( 11) is obtained using the assumption of equipartition of energy in thermal equilibrium.Equation (11) implies that as t → ∞, the probability of reaching any point Γ′ on a given manifold at constant r does not depend on the initial conditions.This is expected to be valid in most models that include a large number of bath degrees of freedom that couple to the relevant subsystem.With Eqs. ( 10) and ( 11), we obtain the following expression for the long-time limit of quasiclassical correlation functions: The correct long-time limit expected for the quantumclassical correlation function as ̵ h → 0 is however given by where tr s [⋅] denotes the partial trace with respect to the subsystem and p) ]. 57 This provides the benchmark result against which we will test the various quasiclassical predictions.While higher-order contributions in ̵ h to the quantum-classical Boltzmann operator do in principle exist, 58 we do not consider them here, as almost all quasiclassical approaches cannot even reproduce the dominant zeroth-order term correctly in all cases.Additionally we show in Appendix A that Eq. ( 13) is independent of the Hamiltonian representation, as long as the ̵ h → 0 limit is taken.The main difficulty for quasiclassical approaches in reproducing this long-time limit is found in the term given by Eq. (13c).This is because the majority of mappings implemented in quasiclassical approaches can only at most reproduce the correct trace relations for a product of two operators, ∫ dxdp tr s [ Ĥ(x, p) B(x, p)] = ∫ dΓ H(Γ)B(Γ), and will therefore not be able to correctly describe all the terms arising from the Taylor expansion of the Boltzmann operator in Eq. (13c).In order to better understand the relative accuracy of different quasiclassical approaches in reproducing the correct longtime limit of correlation functions, we apply this analysis to the specific case of two-state systems coupled to a heat bath.

III. APPLICATION TO TWO-LEVEL SYSTEMS
The arguments derived so far hold for a subsystem consisting of an arbitrary number of quantum levels.For the sake of simplicity, here and in the following we apply our analysis to two-level quantum systems, although the majority of our arguments apply equally well to multi-state problems.In this case, the state-dependent potential in Eq. ( 1) can be written as where ∆(x) and κ(x) denote diabatic Hamiltonian parameters and σi , for i = x, y, z, are the Pauli operators in the diabatic basis.These, together with the 2 × 2 identity, Îs , form a complete set of Hermitian operators for the two-level system.The σad i (x) operators are the Pauli matrices in the adiabatic basis, related to the diabatic σi operators by a linear transformation. 38Finally, V ad z (x) = ∆(x) 2 + κ(x) 2 denotes half the energy gap in the adiabatic basis.
We test our theoretical predictions on two nonadiabatic models in which the state-dependent potential depends on a one-dimensional reaction coordinate, x c .This convenient choice implies that the nuclear phase-space integrals in Eq. ( 12) become one-dimensional, simplifying the interpretation of our results.An f -dimensional secondary bath, with coordinates x 1 , . . ., x f , is introduced to provide friction on the reaction coordinate, such that the dynamics thermalize in the long-time limit.The secondary bath interacts with the subsystem only via the reaction coordinate and can therefore be easily integrated out of the long-time limit expressions. 28The two contributions to the potential in Eq. ( 1) are then expressed as where U RC (x c ) and VRC (x c ) depend on the specific model considered and m = 1 for all degrees of freedom, which corresponds to working in mass-weighted coordinates.We also choose a purely Ohmic spectral density for the secondary bath The use of an Ohmic spectral density means that the dynamical simulations can be easily performed by evolving the reaction coordinate using Langevin equations of motion, which implicitly describe the effects of the secondary bath. 59n the following, we focus on the long-time dynamics of B = σad z (x c ), which corresponds to the population difference between the two adiabatic states.Additionally, we consider a factorized initial condition, where the electronic subsystem is initialized in Â = 1 2 Îs and the phasespace variables associated with the reaction coordinate are sampled from an initial Gaussian distribution where Ω defines its physical width.We call this correlation function C Iz (t).Note that other correlation functions with observables orthogonal to σad z (x c ) [e.g., the coherences σad x (x c ) and σad y (x c )] have to zeroth-order in ̵ h a zero expectation value in the long-time limit by symmetry.This implies that if the long-time limit of C Iz (t) is captured correctly, so will all other correlation functions, as long as they are constructed from the appropriate linear combinations of the adiabatic observables.

A. Spin-boson model
The spin-boson model 60 is commonly used to describe charge transfer in the condensed phase.Despite its apparent simplicity, it provides a stringent test for quasiclassical methods, especially in strongly asymmetric cases due to the problem of negative populations.In the reaction-coordinate picture, the potentials associated with the spin-boson model are 61 where Ω is the frequency of the reaction coordinate (which we also choose to coincide with the width of the initial Gaussian distribution [Eq.( 17)]), ∆ is the coupling between the two diabatic states and α is the system-bath coupling strength.We employ reduced units where mass and ̵ h are 1.The long-time limits of the adiabatic population difference are shown in Fig. 1 as a function of the energy bias, ε, and are calculated with several different quasiclassical approaches.The other parameters are fixed to β = 0.3, ∆ = 1, α = 1, Ω = 1, so that the condition for classical nuclei, βΩ < 1, is satisfied. 62,63This parameter choice is also justified by calculations discussed in Ref. 28, where quasiclassical dynamics in the same parameter regime were found to agree well with numerically exact results with all degrees of freedom treated quantum mechanically. 64,65he validity of our theoretical formula, Eq. ( 12), for predicting the long-time limit of quasiclassical correlation functions has already been demonstrated by us for this model in another recent work, 15 where we found that its predictions matched well with dynamical simulations.
The general robustness of quasiclassical approaches for thermalizing correctly in any parameter regime of this model can be ascertained by considering the two parameter limits of ε → 0 and ε → ∞.Of these extremes, one would expect the ε → 0 limit to be the least challenging for quasiclassical approaches, as for our chosen parameter set of the model, the thermal energy is large compared to all other energy scales, such that the quasiclassical description of both electronic and nuclear degrees of freedom should be valid.It is therefore slightly surprising that a significant number of the approaches we tested fail in this regime, including the commonly used Ehrenfest approach.The ε → ∞ limit is even more challenging, as this leads to unphysical negative populations associated with the excited state in many quasiclassical approaches [i.e., C Iz < −1], as can be seen in Fig. 1.We do however note that the double SEO and single Wigner approaches, despite still having phase-space regions that correspond to negative populations, do thermalize to the correct physical value in the ε → ∞ limit.
The best approaches with regards to their thermalization behaviour for this model are single Wigner, SQC and MASH, which all thermalize reasonably accurately in both the ε → 0 and ε → ∞ limits.As a result, these approaches are observed to thermalize relatively well for all parameter regimes of the model, as illustrated in Fig. 1.In particular, MASH is seen to exactly reproduce our , for several values of the coupling parameter, α.The potentials have been vertically shifted such that the lowest energy minimum of the diabats is located at zero.In the case of α = 0, the two diabatic potentials are identical.
benchmark result across the whole range.

B. An anharmonic model
In order to fully test the effect of inverted potentials on the thermalization behaviour of mean-field approaches, an anharmonic model is needed for which at least one of the diabatic potentials can become unbounded.We choose the reaction-coordinate potentials appearing in Eq. ( 15) to be where 0 ≤ α ≤ 1 determines the strength of the electronnuclear coupling and we choose the other parameters to be ∆ = 1, Ω = 1, xc = 5 and β = 0.3.This model is constructed so that the electronic and nuclear subsystems are uncoupled at α = 0 and the model becomes identical to a previously used electronic predissociation model 66,67 at α = 1.On increasing α, we see in Fig. 2 that the red diabatic well at lower energy becomes progressively less bounded.Additionally, the energy gap between the lower and upper adiabatic states increases, which more significantly drives the system into the lower energy adiabat in the long-time limit.
If a trajectory has a large enough negative population in the purple diabat [Fig.2], so that the contribution from its potential becomes inverted, then the total nuclear force can become unbounded, resulting in the trajectory accelerating off to infinity.This problem can be understood in terms of the mean-field approximation to the FIG. 3. The long-time limits of the adiabatic population difference for the anharmonic model introduced in Sec.III B, as a function of the coupling parameter, α.The lines in the figure correspond to our theoretical predictions given by Eq. ( 12), while the stars are results from dynamical simulations.The simulations were performed by sampling the initial nuclear phase-space variables from Eq. ( 17) and then propagating the trajectories for t = 400 with η = 2Ω in Eq. ( 16).Results are only shown in regions where trajectories cannot become unbounded due to inverted potentials.
state-dependent potential [Eq.(4a)], the magnitude of which can become unphysically large through the multiplication of the Cartesian mapping variables with a norm, r, greater than one.For this model, this results in the following definition of the effective electron-nuclear coupling associated with these methods: α eff = αr.Given that that the red diabatic potential shown in Fig. 2 becomes unbounded when α ≥ 1, we see that this will first occur for the mean-field methods when αr max ≥ 1, where r max is the maximum allowed value of the Cartesian mapping variable norm for that method.We therefore note that the inverted potential problem is particularly problematic for MMST approaches, for which r max = ∞, so that the nuclear force can become unbounded for any non-zero value of α.
Results for the long-time limits of the adiabatic population difference for this model are shown in Fig. 3 for different quasiclassical methods.The lines in the picture correspond to the theoretical predictions from ergodic theory, while the stars are numerical results from dynamical simulations.Results for each method are only provided for the values of α for which the nuclear force cannot become unbounded because of the problem of inverted potentials.We first note that the long-time limits from the dynamical simulations agree well with our simple formula [Eq.( 12)], further confirming the validity of applying classical ergodic theory to predict the longtime limit of quasiclassical correlation functions.Importantly, the agreement confirms that the same formula can also be used to accurately predict the long-time limit of MASH correlation functions, even though MASH does not formally have Hamiltonian generated dynamics, as discussed previously in Sec.II C.
Of all the quasiclassical methods we have tested, MASH is the only one in complete agreement with the quantum-classical benchmark.Its success is largely due to its windowing scheme applied consistently to both the observables and the nuclear force, such that the method does not suffer from the problem of negative populations in either.Although the SQC approach does not deviate too strongly from the exact long-time limit in the parameter regime 0 ≤ α < 0.5 (in part due to the windowing of its observables), the problem of unbounded inverted potentials contributing to the nuclear force (which unlike for MASH are not windowed) means that SQC becomes unstable in the 0.5 ≤ α ≤ 1 regime (because r max = 2 for SQC, as shown in Appendix B).It would of course be possible to simply discount unbounded trajectories, but this would introduce an ad hoc modification to the method, which may affect some of its formal properties.In addition to MASH, Ehrenfest also does not suffer from the problem of inverted potentials (because r = 1) and can also be applied in all parameter regimes of the model.However, we see from Fig. 3 that the Ehrenfest long-time populations are highly inaccurate, as was also observed in Fig. 1.
In Sec.IV, we further analyze our long-time limit formula to better understand the deficiencies in the thermalization behaviour of certain methods, as well as the excellent thermalization properties of MASH.

IV. ANALYSIS
For the two-state models considered in this paper, the benchmark quantum-classical predictions for the longtime limits of the adiabatic populations are obtained from Eq. ( 13) using Â = 1 2 Îs and B = σad z (x).Inserting the expressions for the Hamiltonian [Eqs.( 1) and ( 14)] and explicitly performing both the quantum traces and the integrals over the nuclear momenta, we find that The phase-space average associated with the stateindependent nuclear potential is defined as and we have assumed that ρ0 (x, p) is normalized, ∫ dxdp tr s [ρ 0 (x, p)] = 1.This expression is most easily derived within the kinematic representation.We show in Appendix A that the same result can be obtained using the Hamiltonian in the adiabatic representation, in the ̵ h → 0 limit.
In the following we calculate the equivalent long-time limit expressions for several quasiclassical methods.In order to perform the electronic phase-space integrals necessary to compare our predictions with Eq. ( 20), we choose to work with the spin-mapping variables in the adiabatic basis [Eq.( 7)].

A. Mean-Field Approaches
For mean-field approaches, the long-time limit involves the following operator representations: A(r, Γ) = 1 2 I s (r), B(r, Γ) = rS ad z and ρ 0 (r, Γ) = ρ 0,s (r)ρ b (x, p).Inserting these into Eq.( 12) and performing some of the phasespace integrals, we find Comparing Eq. ( 20) and Eq. ( 22), we see that the expression for the long-time limit of the mean-field correlation function has the wrong functional form, such that there is no universal expression for ρ 0,s (r)I s (r) that guarantees that the corresponding mean-field approach would always thermalize to the quantum-classical benchmark.For a given temperature, a system-specific expression for these functions can however be determined so that Eq. ( 20) and Eq. ( 22) agree.This idea has recently been used to design the "ellipsoid mapping", a new mean-field approach for computing thermal correlation functions that rigorously obeys detailed balance. 28uch an approach however requires a static calculation to be carried out at thermal equilibrium for each system and temperature of interest, in order to calculate the correct associated expression for ρ 0,s (r)I s (r), before any dynamical simulations can be performed.
If we instead decide to only use universal functions for ρ 0,s (r)I s (r), then their form can at best be chosen to ensure correct thermalization in certain parameter limits.In the following, we consider the limiting cases of high temperature, low temperature and weak systembath coupling.For the high-temperature limit (β → 0), we show in Appendix C 1 that mean-field approaches are at best capable of reproducing the benchmark [Eq.(20)] up to first order in β.To achieve this for any two-level system, the condition must be satisfied.Mean-field methods that fulfill Eq. ( 23) are marked with a tick in the β → 0 column of Table I.For methods that do not satisfy this condition, the value associated with the left-hand side of Eq. ( 23) is given, which by comparing Eq. (C1a) and Eq. ( C1b) is seen to be the multiplicative error in the long-time population difference of the method when β → 0. Incidentally, these errors qualitatively explain the deviations of the long-time limit populations from the quantum-classical benchmark in the ε → 0 limit of Fig. 1.While the ε → 0 and the β → 0 limits of a spin-boson model do not always coincide, they do so approximately for the parameter regime that we study here.In the low-temperature limit (β → ∞), mean-field approaches can correctly reproduce the long-time limit of correlation functions to at best zeroth order in e −βV ad z (x) .As shown in Appendix C 2, they must then satisfy Mean-field methods that satisfy Eq. ( 24) are marked with a tick in the β → ∞ column of Table I.For methods that do not satisfy this condition, the value associated with the left-hand side of Eq. ( 24) is given, which is the multiplicative error in the long-time population difference of the method when β → ∞.For strongly asymmetric systems, relaxation in the long-time limit will predominately occur into the lowest potential-energy well, like in the low-temperature limit.This means that our results for β → ∞ can also be used to explain the thermalization behaviour of mean-field methods in the spin-boson model for ε → ∞, given in Fig. 1.We note however that the asymptotic approach to the low-temperature limit, which corresponds to the term that is first-order in e −βV ad z (x) , is always wrong for mean-field approaches, as is also observed numerically from our results in Fig. 1.
We show in Appendix C 3 (and also illustrate with crosses for all the mean-field methods in the α → 0 column of Table I) that it is not possible to design meanfield methods which consistently thermalize correctly in the weak-coupling limit.Despite this, spin mapping and spin-PLDM appear very accurate in the α → 0 limit of the anharmonic model [see Fig. 3].This is because the weak-coupling limit for this particular parameter set of the model also coincides with the β → 0 limit.These methods would not have been as accurate in this limit if we would have instead studied a lower-temperature regime.Let us remark that if the electron-nuclear coupling is identically zero, all methods discussed are able to correctly capture the Rabi oscillations of the isolated electronic subsystem.In this case ergodic theory does not apply, given that the dynamics do not thermalize to equilibrium at long times.Note that this condition is different from the α → 0 limit discussed in Table I and Appendix C 3, where the coupling is assumed to be infinitesimal but nonzero.

B. Symmetric quasiclassical windowing (SQC)
For SQC, the potential still retains the same representation as for the mean-field methods [i.e., σad z (x) ↦ rS ad z in Eq. ( 14)], but the observable operators are now represented by the triangular windows, given by Eq. (B4) in terms of the spin-mapping variables.Inserting these operator representations into Eq.( 12) and performing some of the phase-space integrals gives where N SQC (t) is the normalization factor, which is required to ensure that the electronic populations sum to one.Note that we have applied the ergodic theory outlined in Sec.II D separately to each term in the ratio of Eq. (25a).
A comparison with Eq. ( 20) clearly shows that Eq. ( 25) is not exact, again having the wrong functional form compared to the benchmark quantum-classical long-time limit.However, the advantage of SQC over the other mean-field approaches is that the long-time limit of its correlation function also reproduces the benchmark in the weak-coupling limit, 30 in addition to the high-and low-temperature limits (as indicated in Table I).More details on the thermalization of SQC in these limits is given in Appendix C. The fact that SQC can simultaneously describe the thermalization correctly in all of these limits explains why it is reasonably accurate for all values of ε for the spin-boson model in Fig. 1 and accurate in the α → 0 limit of the anharmonic model in Fig. 2. However for larger values of α, there is a clear deviation even before the unbounded inverted potentials appear at α ≥ 0.5.

C. Mapping approach to surface hopping (MASH)
In MASH, both the potential and observable operators are consistently described through the windowing scheme, σad z (x) ↦ sgn(S ad z ).Following a strategy similar to the one outlined above for the other methods, we find for the long-time population difference of MASH that  17) and propagated for t = 500 with the optimal damping parameter, η = 2Ω [Eq.( 16)].Additionally, the black dashed lines correspond to the expected results from classical ergodic theory, given by Eq. (26b).
where we have additionally used that the initial distribution for the spin-mapping variables for this correlation function is rρ 0,s = W PP (S ad ) = 2 S ad z . 38The marginal equilibrium distribution function for MASH, ρ MASH eq (S ad z ), corresponds to the long-time limit of the conditional probability, given by Eq. (11), with the bath degrees of freedom integrated out.
We first note that our expression for the marginal equilibrium distribution function for MASH agrees with the long-time distribution of S ad z , histogrammed from an ensemble of dynamical MASH trajectories.This is shown in Fig. 4 for the same spin-boson model that was introduced in Sec.III A, for several values of ε.This marginal distribution constitutes an exact discrete-valued representation of the quantum-classical Boltzmann factors, which demonstrates how MASH effectively quantizes the electronic states, enabling it to thermalize to the benchmark quantum-classical distribution in all cases.The ultimate proof of the long-time accuracy of the MASH dynamics is the fact that Eq. ( 20) and Eq.(26a) are seen to be identical once the integrals over the spin-mapping variables have been performed analytically.To the best of our knowledge, MASH is the only nonequilibrium quasiclassical approach that is guaranteed to thermalize correctly in all nonadiabatic systems in the ̵ h → 0 limit.This result may come as a surprise as it is known that MASH does not rigorously obey detailed balance, due to the non-conservation of its weighting factor by the dynamics. 38This raises the question of how MASH can still thermalize correctly.To describe the issue quanti-tatively, a measure of the microscopic reversibility error (MRE) was defined in Ref. 38.Based on the same ergodicity arguments used here, we show in Appendix D that in the long-time limit the MRE vanishes.This shows that the violation of detailed balance by MASH can only lead to errors at intermediate times, with no adverse effects on its long-time relaxation behaviour.

V. CONCLUSIONS
In this paper, we presented an analysis of the long-time limits of quasiclassical approaches for simulating nonadiabatic dynamics.We exploited the Hamiltonian structure of the dynamics to make use of results from classical ergodic theory, rigorously accounting for the conservation of the norm of the mapping variables.This allowed us to assess and compare the accuracy of different methods in capturing the correct detailed balance at long times.Our theoretical framework can be applied to any quasiclassical method which is incompressible (i.e., it conserves phase-space volume) non-integrable (i.e., the dynamics cannot be solved analytically) and deterministic. 68sing our theory, we tested the thermalization behaviour of a wide array of quasiclassical approaches on both the harmonic and anharmonic models.Our analysis revealed that most of the commonly used quasiclassical approaches violate detailed balance and thus do not recover the correct quantum-classical thermal distribution.Errors are particularly large for the anharmonic model, where the problem of inverted potentials meant that trajectories from the majority of quasiclassical approaches were unphysically accelerated off to infinity.Among all of the methods considered, only MASH is guaranteed to predict the exact quantum-classical correlation functions in the long-time limit, to zeroth-order in ̵ h.This method is therefore capable of solving the long-standing issue of detailed balance (at least in the long-time limit), a problem which has plagued the field of quasiclassical dynamics to date.In addition, because MASH does not suffer from the problem of inverted potentials, it is ideally suited for performing ab-initio nonadiabatic simulations of molecules, where the strongly anharmonic adiabatic potentials would pose a significant challenge to the other quasiclassical approaches.
The main limitation of the original version of MASH, introduced in Ref. 38, is that it was derived only for two-state problems.Recently a multi-state generalization of MASH has been proposed by Runeson and Manolopoulos. 69In this version, they modified the observable operators to be more similar to the mean-field approaches, which made it easier for them to generalize the underlying theory to multi-state problems.Like the original MASH approach, it is guaranteed to thermalize correctly in the long-time limit.However, it is not guaranteed to exactly reproduce the short-time dynamics of the quantum-classical Liouville equation, which was the means by which the original MASH approach obtained unique and rigorous prescriptions for the momentum rescalings, treatment of frustrated hops and decoherence corrections.In addition, the measurement of observables is not performed using the same windowing procedure as for the nuclear force, meaning that the two can become inconsistent and the observables may measure negative populations.Whilst the accuracy of the numerical results look promising, 69 given that the rigorous nature and internal consistency of the original MASH approach was its main advantage over FSSH, we still think that there is more work to be done in order to obtain the ultimate multi-state generalization of MASH.
There may however still be a place for mean-field mapping methods.In contrast to typical molecules in the gas phase, many solid-state systems are characterized by a high degree of electronic coherence.Due to the dense bands of electronic states, the treatment of the nonadiabatic dynamics using mean-field approaches seems the most suitable.From our analysis, we conclude that single-Wigner and SQC would be the best mean-field approaches to use in this case.Given that the electronic structure of solids is fairly harmonic, we also expect that the problem of inverted potentials would not be as severe.
It is interesting to compare the conclusions from our analysis to those obtained in previous work.In Refs.11, 24, and 25, it was observed that the unity methods led to improved accuracy compared to other MMST approaches.Similar improvements were observed in Refs.12, 13, 26, 43, and 70 for spin-mapping methods, which also treat the identity operator exactly within the mapping.However, most previous studies were performed on harmonic models with reasonably small energy biases.From our analysis in this paper, we observe that for more extreme systems, the quasiclassical approaches that treat the identity operator exactly do not always perform the best, suggesting that other criteria may also be important to consider.For example, a recent analysis performed on the ab-initio exited-state dynamics in ethylene found that the best approaches in this case were those that had the least phase-space volume associated with negative populations. 713][74] This approach describes the non-unitary dissipative dynamics of the subsystem in terms of a non-Markovian equation of motion.Memory effects in the GQME are captured by kernels that decay on timescales that are in general much shorter than the typical electronic relaxation times.This means that it is often advantageous to calculate the kernels using quasiclassical methods, instead of obtaining the full dynamics directly.However it has recently been shown that using the GQME is not guaranteed to fix detailed balance 15 if the input quasiclassical method is not sufficiently accurate at short times.Given the remarkable long-time accuracy exhibited by MASH, it could be insightful to couple this method to the GQME for de-termining the importance of non-Markovian dissipative effects in the dynamics, as well as to reduce the cost of quasiclassical simulations through computing short-lived memory kernels.
Recently there has been increased interest in understanding the effect of strong light-matter coupling on the properties of matter.6][77] While the short-time dynamics is observed to be fairly accurate, the long-time dynamics is known to suffer as a result of unphysical zero-point energy leakage from these quantum modes. 78For classical light, the presence of an external driving field breaks the thermalization to the Boltzmann distribution.The dynamics converge at long times to nonequilibrium stationary states which can be described in certain parameter regimes by the Floquet-Gibbs distribution. 79,80In both cases, our analysis based on classical ergodic theory could be extended in order to analyze how well quasiclassical approaches describe the correct light-modified long-time properties.
where in the final line we take the ̵ h → 0 limit.Is it then seen that Eqs.(A2) and (A4) are identical in this limit.
we are already assuming in this paper that we are in the high-temperature limit with respect to the nuclear degrees of freedom (so that a classical treatment of the bath is valid), we clarify that by the high-temperature limit here we mean relative to the electronic energy scales.
Using the expressions Eqs. ( 20), ( 22) and ( 25), we obtain the following results for the first-order term in the corresponding high-temperature Taylor expansions ⟨V ad z (x)⟩ b dr r 3 ρ 0,s (r) I s (r), We thus see that SQC is able to correctly thermalize in the high-temperature limit, although we note that this is not the case for all mean-field methods.Comparing Eqs.(C1a) and (C1b) leads to the condition expressed in Eq. ( 23) that must be satisfied for mean-field methods to describe this limit correctly.While the ε → 0 and the β → 0 limits of a spin-boson model do not always coincide, they approximately do for the parameter regime that we study here.This means that these high-temperature limit formulas can be used to understand the thermalization behaviour of quasiclassical approaches in the ε → 0 limit of the spin-boson model discussed in Sec.III A.

The β → ∞ limit
Here we investigate the low-temperature (β → ∞) limit with respect to the electronic energy scales.We still assume that the temperature is high with respect to the nuclear degrees of freedom, such that it is valid to treat them classically.The zeroth-order terms of these expansions are C MF Iz (t → ∞) ∼ − dr r 2 ρ 0,s (r) I s (r), (C2b) We note that SQC also describes the thermalization in the low-temperature regime correctly, while mean-field methods do not unless they satisfy the condition Eq. ( 24).In the long-time limit, the system will predominately relax into the lowest potential-energy well in both the β → ∞ and the ε → ∞ limits of the spin-boson model.This means that our long-time limit formulas for β → ∞ can also be used to understand the thermalization behaviour of quasiclassical approaches in the ε → ∞ limit of the spin-boson model discussed in Section III A. We additionally remark that quasiclassical approaches that reproduce the correct β → ∞ limit are not in general guaranteed to reproduce the right asymptotic approach (apart from MASH), as observed in our numerical results [see in particular Ehrenfest and double SEO in Fig. 1].
3. The α → 0 limit The weak-coupling limit (α → 0) leads to a significant of our long-time limit formulas, because in this limit V ad z (x) = V (0) z becomes independent of the nuclear positions and such terms can be taken outside the phase-space averages, ⟨⋯⟩ b .This means that identical terms in the numerator and denominator of a fraction that previously belonged in different phase-space integrals now cancel.Performing these simplifications leads to the zeroth-order terms We find that SQC is also exact in the weak-coupling limit, as was also shown in Ref. 30. 82 This time, the mean-field condition depends on βV (0) z and so cannot be satisfied in general unless the method is explicitly dependent on this system-dependent variable (similar to Ref. 28).These results are summarized in the corresponding column of Table I. able to exactly reproduce the long-time limit of correlation functions despite not rigorously obeying detailed balance.

FIG. 1 .
FIG. 1.The long-time limits of the adiabatic population difference as a function of the energy bias, ε, for the spin-boson model introduced in Sec.III A. The other parameters are fixed to β = 0.3, ∆ = 1, α = 1, Ω = 1.Inset: the difference between the benchmark quantum-classical result [from Eq. (13)] and the quasiclassical predictions for the same methods shown in the main panel.

FIG. 2 .
FIG.2.The diabatic potentials associated with the anharmonic model,1 4[1 ± α]Ω 2 x 2 c + 1 2 [1 ∓ α]e −Ω(xc−xc), for several values of the coupling parameter, α.The potentials have been vertically shifted such that the lowest energy minimum of the diabats is located at zero.In the case of α = 0, the two diabatic potentials are identical.

FIG. 4 .
FIG.4.The marginal equilibrium distribution function for the mapping variable S ad z , calculated for the spin-boson model defined in Sec.III A, for several values of the energy bias ε.The colored lines are obtained from histogramming MASH trajectories initialized from Eq. (17) and propagated for t = 500 with the optimal damping parameter, η = 2Ω [Eq.(16)].Additionally, the black dashed lines correspond to the expected results from classical ergodic theory, given by Eq. (26b).