Trajectory-based methods that propagate classical nuclei on multiple quantum electronic states are often used to simulate nonadiabatic processes in the condensed phase. A long-standing problem of these methods is their lack of detailed balance, meaning that they do not conserve the equilibrium distribution. In this article, we investigate ideas for restoring detailed balance in mixed quantum–classical systems by tailoring the previously proposed spin-mapping approach to thermal equilibrium. We find that adapting the spin magnitude can recover the correct long-time populations but is insufficient to conserve the full equilibrium distribution. The latter can however be achieved by a more flexible mapping of the spin onto an ellipsoid, which is constructed to fulfill detailed balance for arbitrary potentials. This ellipsoid approach solves the problem of negative populations that has plagued previous mapping approaches and can therefore be applied also to strongly asymmetric and anharmonic systems. Because it conserves the thermal distribution, the method can also exploit efficient sampling schemes used in standard molecular dynamics, which drastically reduces the number of trajectories needed for convergence. The dynamics does however still have mean-field character, as is observed most clearly by evaluating reaction rates in the golden-rule limit. This implies that although the ellipsoid mapping provides a rigorous framework, further work is required to find an accurate classical-trajectory approximation that captures more properties of the true quantum dynamics.

## I. INTRODUCTION

Detailed balance is a crucial concept in both classical and quantum statistical mechanics. It underpins the validity of many important microscopic relations, such as the fluctuation–dissipation theorem, which have lead to computationally powerful algorithms, especially for classical molecular dynamics (MD).^{1} Classical MD provides an internally consistent treatment of detailed balance because Hamiltonian evolution conserves the classical Boltzmann distribution. In a similar way, the exact quantum-dynamical propagator commutes with the quantum Boltzmann operator, which proves that detailed balance should be a property of all quantum systems in thermal equilibrium.

However, when attempting to include quantum effects within classical-trajectory methods, detailed balance turns out to be a much more elusive property.^{2} Within the adiabatic Born–Oppenheimer limit, it is now possible to respect quantum detailed balance (including zero-point energy and tunneling) through path-integral methods, such as centroid molecular dynamics (CMD),^{3} ring-polymer molecular dynamics (RPMD),^{4} and path-integral Liouville dynamics.^{5} These have in many ways replaced previous methods that lacked this property, such as the linearized semiclassical-initial value representation (LSC-IVR).^{6} The situation is, however, quite different when it comes to *mixed quantum–classical* systems. The standard example of such a system is a molecule with classical nuclei and multiple quantum electronic states, but the classical part could as well be any low-frequency environment (not necessarily harmonic), and the quantum part could represent excitons, polaritons, high-frequency vibrations, and so on. Despite significant effort, simulating nonadiabatic dynamics of such systems while obeying detailed balance remains a crucial unsolved problem of statistical mechanics.

After decades of development, a variety of trajectory-based methods have been proposed to simulate such systems. However, to our knowledge, none of these methods is guaranteed to preserve the correct equilibrium distribution without breaking other reasonable limits (such as recovering Rabi oscillations for an isolated quantum system), even when the classical-nuclear approximation is valid. For example, the Ehrenfest approach (also known as mean-field theory) is known to severely violate detailed balance and can therefore not be used to describe relaxation to thermal equilibrium.^{7,8} Another prominent example, fewest-switches surface hopping,^{9} is known to recover detailed balance only in certain limits (small adiabatic splitting or large nonadiabatic coupling).^{10,11} Similarly, the symmetric quasiclassical windowing approach (SQC)^{12} obeys detailed balance in the limit of small diabatic coupling^{13} but not in general. Although simple rate descriptions, such as secular Redfield theory, do obey detailed balance by construction,^{14} these methods are only valid in certain parameter regimes (e.g., weak system–bath coupling and Markovian dynamics).

The aim of this article is to investigate ways of constructing trajectory-based nonadiabatic dynamics that strictly obey detailed balance. Note that the term *detailed balance* is used to encompass several different concepts in the literature. In this paper, we say that a quantum–classical method obeys detailed balance if it initializes equilibrium systems in the correct quantum–classical Boltzmann distribution and preserves this distribution over time, so that the dynamics are microscopically reversible and time-translationally invariant. Apart from being necessary for internal consistency, detailed balance also enables a powerful arsenal of sampling tools developed for classical MD, which are currently not formally applicable to the conventional nonadiabatic techniques mentioned above. Apart from this property, the quantum–classical dynamics should fulfill a few other relevant features in order to provide a realistic description. In the limit of zero nuclear–electronic coupling, they should recover the correct result for an isolated quantum system (known as Rabi oscillations in the two-level case). Finally, the dynamics should reduce to classical MD on the ground-state Born–Oppenheimer potential in the adiabatic limit.

In this paper, we investigate a potential solution to this problem. Here, we limit ourselves to situations where the nuclei can be treated classically, although we note that because no known nonadiabatic extension to RPMD obeys all the required prescriptions above,^{15–19} the present development is also likely to be relevant to tackle quantum nuclei with the path-integral formalism in future work.

Our treatment is based on mapping the quantum subsystem onto a classical counterpart in order to treat both quantum and classical degrees of freedom on an equal footing. This idea has a long history, originally through the Meyer–Miller–Stock–Thoss (MMST) mapping,^{20,21} which uses a set of harmonic oscillators as the mapping space. This mapping has lead to a variety of classical-trajectory methods^{12,22–28} that evolve the nuclei with a mean-field force, reminiscent of the Ehrenfest approach, but start from a different initial distribution and typically lead to higher accuracy. However, these methods still break detailed balance, and in certain situations, the populations (weights of nuclear forces) may even become negative, which means that the nuclei effectively move on inverted potentials and can lead to unphysical predictions. A possible solution proposed by Müller and Stock in the late 1990s is to modify a parameter that can be thought of as the zero-point energy (ZPE) of the oscillators, so as to decrease the likelihood of negative populations.^{29} They proposed a criterion for the optimal value of this parameter based on the long-time limit of the dynamics, or in other words, an attempt to recover detailed balance. Further generalizations of this parameter have recently been suggested,^{30} but without addressing the issue of detailed balance.

In this paper, we do not employ mapping to harmonic oscillators but instead use a phase space for *N*-level systems known as the Stratonovich–Weyl (SW) representation, which generalizes the concept of a Wigner representation to systems with symmetry SU(*N*). Several of the methods originally developed for the MMST mapping have now been adapted to the SW framework and typically lead to improved accuracy.^{31–34} For two-level systems, the SW phase space is closely connected to the classical-vector model of a spin-$12$ with radius $32$, which is why the corresponding methods are referred to as “spin mapping.” The spin-mapping equations of motion have the same form as in Ehrenfest and MMST mapping but start from a different (“spherical”) initial distribution and use a value of the zero-point energy parameter that is uniquely determined by the spin magnitude. (One can also construct a spherical distribution within the framework of MMST mapping by constraining the total population,^{35,36} although the special value of the ZPE parameter is less apparent in this approach.) Incidentally, the value used in spin mapping provides the optimal high-temperature approximation to the long-time populations, as will be demonstrated in this paper. At low temperature, however, spin mapping is known to suffer from the same inverted-potential problem as other mapping methods. One way to understand this problem is that when the upper state is very high in energy compared to *k*_{B}*T*, the system should effectively be treated as a single-level system (the ground state) rather than a two-level system, and in such cases, an SU(2) mapping is no longer appropriate.

In this paper, we show that the problem of detailed balance, as well as the inverted-potential problem, can be solved (at least in the two-level case) by generalizing the SW representation to a form that is more appropriate for equilibrium dynamics. This leads to a new phase-space framework that takes into account which states are thermally accessible and preserves detailed balance by construction. In the new theory, one can think of the “spin” as evolving not on a sphere, but on an ellipsoid whose center and shape can adapt along the trajectories. The ellipsoid dynamics preserve Rabi oscillations and, in the limit that the states are well-separated compared to *k*_{B}*T*, reduce to adiabatic dynamics on the lower state. In the opposite limit of states close in energy compared to *k*_{B}*T*, the theory reduces back to the original (spherical) spin mapping. Away from these two limits, the theory provides a gradual transition between two- and one-level systems.

After deriving the theory in Sec. II, we assess the accuracy of the ellipsoid mapping in Sec. IV by computing thermal correlation functions for a spin–boson model, particularly in the strongly asymmetric regime where the standard mapping methods break down. Such correlation functions are relevant to calculate a variety of experimental properties, from spectra to reaction rates, through their connection to linear-response theory. The present theoretical framework can therefore be applied to the study of a variety of physical processes and phenomena, such as exciton transfer, vibronic spectroscopy, molecular junctions, and heat transport.^{37,38}

## II. THEORY

We start by considering a general two-level nonadiabatic system described by the Hamiltonian in the diabatic representation

where *V*_{1} and *V*_{2} are potential-energy surfaces and Δ is the coupling between the two levels. Typically, $x\u0302$ and $p\u0302$ would represent the nuclear degrees of freedom in an electronic–nuclear problem, but the same model could be used to represent any two-state quantum subsystem in an environment. Although the present treatment is limited to two levels, we expect that it can be generalized to multiple levels in the future, similarly to the standard spin-mapping approach.^{32}

When the system is in thermal equilibrium, its expectation values are given by

where *β* = 1/*k*_{B}*T*, the partition function is $Z=Tr[e\u2212\beta H\u0302]$, and the quantum-mechanical trace is taken over both subsystem and environmental degrees of freedom. Many experimentally measurable properties of the system, such as spectra,^{2,39} rates,^{40} and scattering functions,^{41} can be expressed in terms of equilibrium time-correlation functions. The standard correlation function between two operators $A\u0302$ and $B\u0302$ is defined as

where, throughout this paper, we use units for which *ℏ* = 1. In general, thermal correlation functions depend on the ordering of the operators inside the trace. The ordering that is most closely related to linear-response theory and experimentally measurable quantities^{42} is not the standard one above, but the Kubo-transformed correlation function, defined by

Both *C*_{AB} and *K*_{AB} contain the same information and either of them can be computed from the other through a simple relation between their Fourier transforms.^{4} However, *K*_{AB} is more closely related to classical correlation functions as they are both real and have the same symmetry under time-reversal, in contrast to *C*_{AB}. Adiabatic path-integral methods, such as CMD and RPMD, have therefore focused on this quantity when including nuclear quantum effects. Even though their dynamics are fictitious, what makes these approximations to quantum dynamics particularly appealing is that they reproduce quantum statistics not only at *t* = 0, but also at later times. One therefore says that they conserve the Boltzmann distribution.^{43} This is a direct consequence of obeying detailed balance.

Following the same spirit, we will in this paper present a theory to approximate *K*_{AB} for nonadiabatic problems. For now, we limit the discussion to classical nuclei, which is a reasonable assumption if *k*_{B}*T* is large compared to the zero-point energies of the nuclear modes (but not necessarily compared to the electronic energy scales). The assumption of classical nuclei means that we replace $x\u0302$ and $p\u0302$ by classical phase-space variables *x* and *p*, while the electronic dynamics are still treated quantum-mechanically.

To define equilibrium properties within this mixed quantum–classical framework, we introduce the quantum–classical density matrix^{44}

and the partition function

where the quantum–classical trace is to be understood as a classical phase-space integral over the nuclear (environmental) part and a quantum trace over the electronic (subsystem) part,

Based on this prescription, expectation values are defined as

Likewise, we define the quantum–classical limit of the Kubo-transformed correlation function at zero time,

where the dependence on the nuclear variables has been suppressed for brevity.

The expressions in Eqs. (8) and (9) only report on statistical properties and it is not obvious how to generalize them to time-evolved quantities without breaking time-translational invariance. For example, even the quantum–classical Liouville equation (QCLE), which has been used to derive many classical-trajectory methods, is known to break this property in general^{45} (although it is exact for the spin–boson model). One strategy to define a quantum–classical limit for $B\u0302(t)$ is to express not only the environment but also the subsystem in a classical phase-space picture, so that nuclear and electronic degrees of freedom are treated on the same footing. This is the philosophy of quasiclassical methods, such as the MMST and spin mappings. In Secs. II A and II B, we summarize the main features of spin mapping and demonstrate that it tends to the wrong thermal distribution when the energy gap is large compared to *k*_{B}*T*, which makes it prone to the same inverted-potential problems as the MMST mapping. To overcome these problems, we will generalize the SW representation to the equilibrium case in Sec. II C and use it to develop a new ellipsoid mapping method in Secs. II D and II E.

### A. Summary of (spherical) spin mapping

Most previous work on mixed quantum–classical dynamics has focused on calculating nonequilibrium correlation functions with factorized initial conditions,

where $\rho \u0302b$ is a purely nuclear density. This density is typically represented either as a Wigner function or as a classical Boltzmann distribution for the bath Hamiltonian (i.e., not the total Hamiltonian). The heart of the classical mapping technique is to define an analogous representation for the electronic operators. In this section, we summarize the main points of one such method, called spin mapping.^{31}

As is well known, a two-level system is isomorphic to a spin-$12$ in a magnetic field since the Hamiltonian is of the form

where $\sigma \u0302\mu $ are the Pauli matrices (including $\sigma \u03020$ as the 2 × 2 identity) and repeated Greek indices imply summation from 0 to 3, whereas Latin indices will be used for 1 to 3. The explicit relations between Eqs. (1) and (11) are $H0=p22m+12(V1+V2)$, *H*_{1} = Re Δ, *H*_{2} = Im Δ, $H3=12(V1\u2212V2)$, and the corresponding “magnetic field” is the vector ** H** = (

*H*

_{1},

*H*

_{2},

*H*

_{3}).

^{46}The operators $A\u0302$ and $B\u0302$ can be expanded in a similar way.

A convenient phase-space construction for the spin degree of freedom is provided by the so-called Stratonovich–Weyl representation,^{47} which has recently gained attention^{48,49} after having been largely overlooked for many decades.^{50} The SW representation can be thought of as the finite-level version of the more widely known Wigner representation of continuous degrees of freedom. In the SW representation, the Pauli operators are replaced by classical functions on the Bloch sphere according to the prescription

where ** u** = (

*u*

_{1},

*u*

_{2},

*u*

_{3}) = (sin

*θ*cos

*φ*, sin

*θ*sin

*φ*, cos

*θ*) and $gW=3$. Here,

*θ*and

*φ*are the usual Bloch-sphere angles. Then, $HW=H\mu \sigma \mu W$ and similarly for other operators. (Note that functions of the Pauli operators, such as the exponential, need to be expanded into a linear combination before this rule can be applied.) The additional factor

*g*

_{W}guarantees that traces of products of operators are equal to the corresponding classical integrals,

where $du=12\pi d\phi d\theta sin\u2061\theta $ denotes integration over the unit sphere (normalized such that $\u222bdu=Trq[\sigma \u03020]=2$). Note that the factor *g*_{W} also appears in the familiar expression for the magnitude $32=S(S+1)$ of a spin $S=12$, hence the name “spin mapping.”

Because this formalism ensures that quantum traces are equal to classical integrals, Eq. (10) can be approximated by a classical correlation function

where the time evolution is given by

and *H*_{W} = *H*_{0}(*x*) + *g*_{W}** H**(

*x*) ⋅

**. Equation (14) is exact in the case of an uncoupled subsystem. In general, however, this independent-trajectory treatment is only exact at**

*u**t*= 0 as it neglects higher-order terms in the Moyal series analogously to the LSC-IVR approximation (classical Wigner dynamics). It is therefore referred to as “spin-LSC.”

In recent work, the SW representation for the electronic states (and its generalization to more than two states) has been employed to simulate various types of nonadiabatic dynamics initialized in nonequilibrium factorized states.^{31–34,51–54} However, the theory is not yet in a suitable form for equilibrium correlation functions, as we will demonstrate in the following.

#### 1. Problems of spin mapping in equilibrium

So far in Sec. II A, we have focused on the case when the zero-time observable is a product state between the subsystem and the environment. To represent the thermal correlation function based on an initial density $e\u2212\beta H\u0302$ in Eq. (3) or (4), one might (naively) attempt to use the form $1ZW\u222bdxdpdue\u2212\beta HWAWBW(t)$, where the symbol “W” refers to the SW representation of subsystem operators [Eq. (12)] and $ZW=\u222bdxdpdue\u2212\beta HW$. However, this does not recover the correct initial distribution because, like Wigner transforms, SW can only represent traces of two spin operators correctly but not higher-order combinations, i.e., $e\u2212\beta HW\u2260[e\u2212\beta H\u0302]W=\u2211n1n!(\u2212\beta )n[H\u0302n]W$. Hence, this does not recover *K*_{AB}(*t*) even for *t* = 0.

To get the correct zero-time correlations, one could, in principle, use the Kubo-transformed operator $1\beta Zqc\u222b0\beta d\lambda e\u2212(\beta \u2212\lambda )H\u0302A\u0302e\u2212\lambda H\u0302$ in place of $\rho \u0302bA\u0302$ in Eq. (10). (This expression is easily evaluated since all operators have 2 × 2 matrix representations within the classical-nuclear treatment). More details on this “Kubo-transformed” formulation of spin mapping can be found in Appendix C 1. However, even though the initial values of thermal correlations are now corrected , the ensemble dynamics under *H*_{W} still do not preserve thermal expectation values. In the long-time limit (assuming ergodicity), the distribution would actually tend towards $e\u2212\beta HW$ [see Eq. (C5b)], which we have already concluded is incorrect. As a consequence, spin-mapping time-evolved expectation values ⟨*B*(*t*)⟩^{W} [as defined in Eq. (C2a)] are not constant in time, in complete contradiction to the concept of equilibrium.

A further issue is that spin mapping in its current form is known to suffer from negative populations, which violates the rule that the density matrix should be positive definite. On the level of individual trajectories, negative populations cause the nuclei to effectively evolve on inverted potentials, which can lead to unphysical dynamics at least for steep potentials. This problem is also present for quasiclassical methods based on MMST mapping (but notably not in the Ehrenfest approach).

### B. Spherical mapping with an optimized radius

In this section, we demonstrate that populations predicted from spin mapping become negative when the energy separation of the two levels is large compared to *k*_{B}*T*, and subsequently discuss possible solutions. According to ergodic theory, it is assumed that all nontrivial systems (with a non-zero coupling to a thermal environment) will relax to an equilibrium distribution regardless of their initial (nonequilibrium) state. It is thus necessary to determine the relevant expectation values of this equilibrium distribution. Our point is most simply illustrated by considering an isolated subsystem, so that we only have to consider traces over the electronic degrees of freedom.

Since it is always possible to find a basis in which *H*_{1} = *H*_{2} = 0 (and because the theory is basis-independent), it suffices to consider the case $H\u0302=H3\sigma \u03023$. Then, $\u27e8\sigma \u03021\u27e9q=\u27e8\sigma \u03022\u27e9q=0$, while the remaining expectation value is

where *ζ* ≡ *β**H*_{3}. This quantum result is to be compared with the corresponding classical expectation value that arises from the mapping $\sigma \u03023\u21a6\sigma 3=gcos\theta $, which is

where $L(x)\u2261cothx\u22121x$ is the Langevin function (a limiting case of the Brillouin function, which appears frequently in the paramagnetic theory of spins^{55}). For $g=gW=3$, the two expressions agree closely for small *ζ*, as shown in Fig. 1, meaning that the classical phase-space average is valid if temperature is large or the energy separation between the two levels is small. A simple Taylor expansion shows that the error is of the order $O(\zeta 3)$, a result unique to $gW=3$, which is thus the optimal scaling in this regard.^{56} For comparison, Ehrenfest (corresponding to a sphere with scaling *g* = 1) and the original MMST mapping (*g* = 2) both have an error of $O(\zeta )$. For larger *ζ*, however, the classical spin analogy breaks down when the high-energy level becomes energetically inaccessible. Although it may seem that Ehrenfest improves again in this limit, one should keep in mind that the infinite-*ζ* case corresponds to an adiabatic system, where the direct use of Born–Oppenheimer MD would be preferable.

To improve the spin-mapping approach at larger values of *ζ*, a simple solution would be to replace the factor $gW=3$ with a *ζ*-dependent scaling factor *g*_{ζ}, defined as the unique (numerical) solution to tanh *ζ* = *g*_{ζ}*L*(*g*_{ζ}*ζ*). This procedure can be extended to also include the environment in the optimization of *g*, as described in Appendix C 2. A similar idea has been considered already by Müller and Stock in the context of the MMST mapping,^{29} where the so-called zero-point energy parameter *γ* = *g* − 1 was optimized to obtain the correct long-time populations. However, in both mapping formalisms, this value is dependent on the choice of basis and only fixes one of the three expectation values (except in the special case of an isolated system). Furthermore, this solution is not sufficient to preserve the equilibrium distribution (i.e., it does not give time-independent expectation values), as we will demonstrate in Sec. IV.

It is also worth noting that by windowing the mapping variables, the SQC method predicts the exact long-time populations in this particular example.^{13} Nonetheless, there are number of reasons why we do not employ the SQC approach in this work. First, it has not been generalized to calculate thermal correlation functions. Second, the results depend on the choice of the basis. Finally, even when it predicts the correct (positive) populations on average, SQC can still suffer from the inverted-potential problem for individual trajectories.^{57}

### C. Desired properties of a phase-space representation for thermal equilibrium

In order to modify spin mapping such that it is more suitable to treat equilibrium correlation functions (as opposed to the nonequilibrium case of Sec. II A), we start by analyzing the defining properties of the self-dual or “W” Stratonovich–Weyl representation, which is the mathematical framework underlying the current spin mapping. These are (for the subsystem degrees of freedom) as follows:^{50}

Linearity: $A\u0302\u21a6A(u)$ is a one-to-one linear map.

- Reality:(18a)$[A\u0302\u2020](u)=[A(u)]*.$
- Normalization:(18b)$Trq[A\u0302]=\u222bduA(u).$
- Tracing:(18c)$Trq[A\u0302B\u0302]=\u222bduA(u)B(u).$
- Covariance:where(18d)$[U\u0302(G\u22121)A\u0302U\u0302(G)](u)=A(G\u22c5u),$
*G*∈ SU(2) is a linear rotation operator and $U\u0302(G)$ its unitary representation.

Of these, property 4 is the one that enables an exact phase-space expression for the correlation functions [Eq. (10)] at zero time (and 3 is just a special case where $B\u0302$ is the identity). In other words, it connects the trace product of quantum operators, $(A\u0302,B\u0302)q\u2261Trq[A\u0302B\u0302]$, with the “classical” inner product of functions on a sphere, (*A*,*B*)_{c} ≡ *∫*d*u**A*(** u**)

*B*(

**). However, there is flexibility in the definition of the inner product and these are not necessarily the best choice for the thermal case. In fact, an inner product that is more closely related to the canonical (Kubo-transformed) correlation in Eq. (4) is**

*u*where $\rho \u0302=1Zqe\u2212\beta Hi\sigma \u0302i$ and $Zq=Trq[e\u2212\beta Hi\sigma \u0302i]$. The classical analog of this inner product is

where $\rho (u)=1Zce\u2212\beta Hi\sigma i(u)$ and $Zc=\u222bdue\u2212\beta Hi\sigma i(u)$.

To make the Stratonovich–Weyl representation more suitable for this equilibrium problem, we propose to modify properties 3 and 4 to

- 3′.Preservation of averages:(21a)$\u222bdu\rho (u)A(u)=Trq[\rho \u0302A\u0302].$
- 4′.Preservation of inner products:(21b)$\u222bdu\rho (u)A(u)B(u)=\u222b01d\lambda Trq[\rho \u03021\u2212\lambda A\u0302\rho \u0302\lambda B\u0302].$

Note that in the limit *β* → 0, these reduce back to the original SW properties 3 and 4. Again, 3′ is just a special case of 4′ when $B\u0302$ is the identity. We leave property 5 unchanged, since it is the key to defining dynamics in phase space.^{58} The two modified properties are not obeyed by the standard mapping [Eq. (12)]. The problem of finding a new mapping $\sigma \u0302\mu \u21a6\sigma \mu (u)$ to replace Eq. (12) in order to fulfill these generalized Stratonovich–Weyl properties is the topic of Sec. II D (for the case of an isolated subsystem) and Sec. II E (for a mixed quantum–classical system).

As a final remark, an alternative way to express the generalized conditions is in terms of the *generating functions*

Then, properties 3′ and 4′ are equivalent to matching the first two cumulants,

Finally, we will require *F*_{q}(**0**) = *F*_{c}(**0**) in order to match the partition functions.

### D. Ellipsoid mapping for isolated subsystems

We are now ready to construct a mapping $A\u0302\u21a6A(u)$ that incorporates the properties described in Sec. II C. Since one can decompose any Hermitian operator as $A\u0302=A\mu \sigma \u0302\mu $, it suffices to define the mapping for each basis operator. As was made clear in Sec. II B, a simple scaling as in Eq. (12) is not enough. To proceed, it is helpful to consult a central theorem of quantum information theory,^{59} which says that any trace-preserving quantum operation of the Bloch vector can be expressed as an *affine map*

where $g$ is a real 3 × 3 matrix (with elements *g*_{ij}) and ** c** a translation vector. Furthermore, the matrix $g$ can be decomposed into a real orthogonal matrix (pure rotation) and a real symmetric matrix (deformation). Since thermal averages are computed from integrals over

**, they are unaffected by pure rotations. Therefore, we can take the matrix to be symmetric. Based on these considerations, we propose a mapping of the form**

*u*where *g*_{ij} = *g*_{ji}. Geometrically, one can think of this mapping as a transformation of the unit sphere into an ellipsoid, where the elements *c*_{i} specify its center and *g*_{ij} its shape.

For completeness, we map the identity to one, $\sigma \u03020\u21a61$, as usual (so that the theory is invariant to a global shift of all energy levels). In order for the mapping to preserve the partition function (equivalent to the zeroth cumulant of the generating function), we additionally introduce a scalar energy parameter $H\u0303$ in the mapping Hamiltonian,

If ** u** is thought of as an effective spin direction, $g$ is analogous to the anisotropy tensor used to describe effective spins in electron paramagnetic resonance (EPR).

^{60}

In total, we have introduced ten independent parameters through the quantities $g$, ** c**, and $H\u0303$. These are defined to fulfill the requirements,

which constitute a system of ten nonlinear equations (because the inner products are symmetric).

This problem is most easily solved in a principal-axis basis where the *z*-axis is aligned with the magnetic field, as indicated in Fig. 2. (Rotating to this basis corresponds to diagonalizing $H\u0302$.) In this basis, the only non-zero expectation values are $\u27e8\sigma \u03023\u27e9q$ and $\u27e8\sigma \u0302i,\sigma \u0302j\u27e9q$ for *i* = *j*. To match the system of equations above, the same needs to be true for $\u27e8\sigma i\u27e9c$ and $\u27e8\sigma i,\sigma j\u27e9c$. This is achieved if $g$ is chosen to be diagonal in the same basis because then the Boltzmann factor does not depend on *φ* and the integral over this variable is non-zero only for *i* = *j*. Furthermore, the rotational symmetry around the *z*-axis also implies that $g11p=g22p$. These considerations lead to the simpler form

where “p” refers to the principal-axis basis. The transformation step to and from the original basis is given in Appendix A. There, we also show that in the principal-axis basis, the system of Eq. (27a) reduces to a single nonlinear equation that is easily solved numerically and uniquely determines all the remaining parameters ($H\u0303$, $c3p$, $g11p$ and $g33p$). The solution is then converted back to the original basis to perform the simulation.

To analyze the solution, Fig. 3 shows all non-zero elements as a function of $\zeta =\beta H3p$, which is the only relevant parameter for an isolated subsystem. The insets depict the ellipsoid for three different temperatures or field strengths (i.e., separation of the eigenenergies). First, one may note that for high *T* or small energy separation, $c3p$ and $H\u0303$ vanish while $g33p$ and $g11p$ approach $3$, meaning that we recover the original “W-sphere” spin mapping (Sec. II A). In the opposite limit (low *T* or large energy separation), the ellipsoid collapses to a point at distance 1 from the origin in the opposite direction from the magnetic field, and since $H\u0303$ is vanishingly small compared to $H3p$, this puts the spin vector into the adiabatic ground state. Between these limits, we observe that the spin vectors ** σ** lie on an oblate ellipsoid (i.e., $g33p<g11p=g22p$), which changes smoothly along the transition between two-level and effectively one-level systems. The parameter $H\u0303$ does not influence the shape of the ellipsoid but can be thought of as a shift of the energy.

Before we couple the two-level system to the nuclear degrees of freedom, we make a brief comment about the spin dynamics. At this point, one may define the equations of motion either in terms of ** σ** as

or equivalently in terms of ** u** as

These equations describe the usual precession of the classical spin vector around the magnetic field. The circular motions generated by both equations are equivalent because ** H** is always aligned with one of the principal axes of the ellipsoid, around which there is a circular symmetry. Just like the spherical spin mapping, the dynamics are exact in the case of an isolated subsystem such that

In particular, the Boltzmann distribution is conserved, $\u27e8\sigma i(t)\u27e9c=\u27e8\sigma i\u27e9c=\u27e8\sigma \u0302i\u27e9q=\u27e8\sigma \u0302i(t)\u27e9q$.

### E. Ellipsoid mapping for mixed quantum–classical systems

We are now ready to include coupling to nuclear degrees of freedom (or in general to any type of classical environment). Our objective is to define a phase-space density for which

and define dynamics that preserve equilibrium expectation values and reduce to the exact result in the isolated case. By the notation “cc,” we mean “fully classical” expectation values

where

In order to make the fully classical and quantum–classical expressions agree at zero time, we need to further generalize the construction of the ellipsoid beyond that of Sec. II D, to take the additional nuclear degrees of freedom into account. To this end, we consider two different approaches. First, we attempt a *global* construction in which $g$, ** c** and $H\u0303$ are independent of

*x*. This is straightforward to implement for simple harmonic models where the nuclei can be integrated out analytically, but would be unpractical for more general potentials. Second, we consider a

*local*construction with position-dependent $g(x)$,

**(**

*c**x*), and $H\u0303(x)$, which makes use of the solution from Sec. II D for each configuration

*x*. This approach can be applied relatively easily to general potentials and uses an ellipsoid that changes shape and position along the trajectories.

Before we describe these two constructions in detail together with their equations of motion, let us make a list of properties that we wish the dynamics to obey:

Conserve energy, Eq. (26).

Remain on the space |

|*u*^{2}= 1.- Obey phase-space incompressibility
^{1}such that Liouville’s theorem is valid.(34)$\u2207u\u22c5u\u0307+\u2202\u2202xx\u0307+\u2202\u2202pp\u0307=0$ Recover the exact Rabi oscillations in the limit of zero electron–nuclear coupling.

Reduce to evolution on the adiabatic ground state in the limit of large

*β*||.*H*- Visit the phase-space point (
*x*,*p*) with probability(35)$1ZqcTrq[e\u2212\beta H\u0302(x,p)].$

Together, these properties guarantee preservation of the equilibrium distribution and ensure physical dynamics in important limiting cases. In particular, property 2 is equivalent to the requirement that ** σ** remains on the ellipsoid surface. Along with Eq. (31), the classical equilibrium distribution is constructed to reproduce the correct quantum–classical results. Property 6 additionally ensures that the expectation value of nuclear operators are also correct.

#### 1. Global construction

First, we attempt a global solution where $g$ is independent of the nuclear variables. For this purpose, the solution in Sec. II D is no longer valid, but instead we propose to optimize $g$ numerically as to solve Eq. (31) (a nonlinear multidimensional root problem of ten unknowns). The additional nuclear integrals make this problem considerably harder to solve than the isolated case, but on the other hand, it only has to be solved once at the start of the simulation. For the special case of a spin–boson model, most of the bath degrees of freedom can be integrated out analytically using the reaction-coordinate representation of the Hamiltonian (see Sec. III and Appendix B).

Compared to the isolated case, an important difference is that the local magnetic field is no longer aligned with any of the principal axes of the ellipsoid. This means that Eqs. (29a) and (29b) are no longer equivalent, and one needs to construct the electronic equations of motion with more care. In fact, neither of these two options fulfills the requirements that we set up above: precession of ** σ** about

**leaves the surface of the ellipsoid, whereas precession of**

*H***about**

*u***(which does ensure |**

*H***|**

*u*^{2}= 1) does not preserve the mapping Hamiltonian

In order to preserve energy as well as remain on the phase space, we instead propose to let ** u** precess around the direction of $gH$. Then, requirements 1 and 2 are simultaneously fulfilled. In addition, we note that the correct Rabi oscillations are obtained in the absence of electron–nuclear coupling if we pick the precession frequency such that

where $geff=|gH||H|$. One may think of *g*_{eff} as an effective radius and the equations of motion correspond to Hamiltonian dynamics with conjugate variables (*φ*, *g*_{eff} cos *θ*). The direction of precession is the same as for “physical” spins with an anisotropic $g$-tensor in EPR,^{61} but the precession frequencies differ by a factor *g*_{eff} due to the different physical significance of this quantity.

Finally, it is clear that, along with the nuclear equations of motion in Eqs. (15a) and (15b), the dynamics preserve the Hamiltonian and fulfill point 3 (using $\u2207u\u22c5u\u0307=0$). Note that the term $H\u0303$ does not influence the dynamics because it is independent of *x* and may therefore be omitted within the global construction.

However, although the parameters of the global ellipsoid can be chosen to fulfill $\u27e8\sigma \mu ,\sigma \nu \u27e9cc=\u27e8\sigma \u0302\mu ,\sigma \u0302\nu \u27e9qc$, there is not enough flexibility to ensure that the equilibrium distribution over the nuclear phase space, $1Zcc\u222bdue\u2212\beta H(x,p,u)$, will recover the correct quantum–classical result according to property 6. Thus, it only fulfills points 1–5 above but not point 6. This means that equilibrium averages will be correct only for electronic operators, but that in general nuclear expectation values will be incorrect.

#### 2. Local construction

To overcome the drawbacks of the global construction and fulfill all points 1–6, we consider an alternative local construction where the mapping parameters $g(x)$, ** c**(

*x*), $H\u0303(x)$ are functions of nuclear configuration. Effectively, this means that each

*x*defines a two-level system with a local potential matrix $V\u0302(x)$. To continue, define the local density $\rho \u0302(x)=1Zq(x)e\u2212\beta V\u0302(x)$, where $Zq(x)=Trq[e\u2212\beta V\u0302(x)]$ and the trace is taken only over the electronic subsystem (leaving

*x*unaffected). Next, define the local quantities

Analogously, define the phase-space density $\rho (x,u)=1Zc(x)e\u2212\beta V(x,u)$, where *Z*_{c}(*x*) = *∫*d** u**e

^{−βV(x,u)}and the local potential is

with $V0=12(V1(x)+V2(x))$. Finally, we define

Then, the solution from Sec. II D can be used to ensure that

In contrast to the global construction, the local ellipsoid will always be oriented with one principal axis along the magnetic field. Then, the spin equation of motion, Eq. (37), reduces to the simpler Eq. (29b). For the nuclear dynamics, the force needs to be modified compared to Eq. (15b) to take into account that *σ*_{i} depend on *x*. Unfortunately, there is not a unique way to determine the equations of motion. Here, we suggest the following:

By construction, these equations guarantee that not only electronic but also nuclear expectation values are correct and conserved, ensure incompressibility of phase space, recover Rabi oscillations for an uncoupled system, and reduce to single-surface dynamics on the lower adiabat in the large-|** H**| limit. Thus, this proposal obeys all six of the properties listed above. Along a trajectory, the ellipsoid deforms and shifts according to the local potential-energy matrix. In regions of space where the excited state is high in energy compared to

*k*

_{B}

*T*, the ellipsoid contracts so that the nuclei feel the ground-state force. When the states are close in energy, the ellipsoid becomes spherical and allows population transfer through precession of the spin vector. Due to the reshaping of the ellipsoid, Eq. (29a) is

*not*equivalent to Eq. (29b), but instead

**has the more complicated equations of motion**

*σ*In practice, we sample and evolve the ** u** variables and only convert to

**when calculating observables. The rotations of**

*σ***were carried out in the standard way using quaternions.**

*u*## III. MODEL

In the last decades, the spin–boson model proved to be a valuable tool to investigate decoherence and population transfer in nonadiabatic dynamics.^{62–65} Despite its simplicity, this model has been shown to accurately describe a large class of relevant processes, from electron transfer to current flux in superconducting circuits.^{66–69} The out-of-equilibrium dynamics of the spin–boson model has been studied extensively in the literature.^{70–73} However, fewer studies have tackled the problem of thermal equilibrium correlation functions^{15,63,74–76} as we do in the present work.

The Hamiltonian of the model can be split into three terms,

describing the uncoupled electronic subsystem, a harmonic bath of *F* modes and the system–bath coupling interaction, respectively. These terms are defined by

Here, *ɛ* denotes half the energy bias between the two electronic levels and Δ is the coupling constant between them. The frequencies *ω*_{α} and the coupling constants *c*_{α} are related by the nuclear spectral density

and we set the masses of all nuclear modes to *m*_{α} = 1. Equation (46) is constructed as a discretization of a Debye spectrum

where *η* denotes the strength of the system–bath coupling (half the reorganization energy) and *ω*_{c} the characteristic frequency. We observed satisfactory convergence of the correlation functions studied in this work with *F* = 200 bath modes using a logarithmic discretization scheme.^{77}

To justify a classical treatment of the nuclear statistics and dynamics, we consider here the regime of high temperature *β* = 0.2 and small nuclear frequencies *ω*_{c} = 1. As we will discuss in later sections, the agreement between full quantum and quantum–classical results validates such choice of parameters. The electronic coupling constant is fixed to Δ = −1. Numerically exact benchmark correlation functions are calculated using the hierarchical equation of motion (HEOM) method.^{78} The standard thermal correlations obtained by HEOM are then converted to Kubo-form, as explained in Appendix D.

For the spin–boson model it is particularly convenient to calculate the exact quantum–classical thermal averages defined in Eqs. (8) and (9). Given that the electronic–nuclear coupling in Eq. (45c) is a linear function of the configurations, it is possible to integrate out analytically all nuclear modes in these phase-space averages, except for a one-dimensional reaction coordinate

We refer to Appendix B for details on these calculations.

## IV. RESULTS

In this section, we test the global and local ellipsoid methods as well as the optimized sphere and compare their results with the original spin-LSC method (adapted for Kubo-transformed thermal correlation functions). For the spin–boson model described in Sec. III, we consider a wide range of values of the electronic–nuclear coupling constant *η* and of the energy bias *ɛ*.

Of all correlation functions *K*_{ij}(*t*), the hardest to reproduce with mapping techniques tends to be *K*_{33}(*t*), especially for large *ɛ*. This is due to the fact that it relaxes to a nontrivial positive limit, which can be poorly captured as a consequence of the negative-population problem introduced in the discussion of Fig. 1. The time evolution of *K*_{33}(*t*) is shown in Fig. 4 as computed with the various methods for comparison with the HEOM benchmark. First, conventional spin-LSC (yellow lines) is found to predict incorrect long-time dynamics, especially for large *ɛ*. The details of how to apply spin-LSC to an equilibrium problem are given in Appendix C 1. There, we show that the long-time limit of this method is

where

The values predicted by Eq. (49) are marked in Fig. 4 by yellow triangles, which agree with the long-time limits of the yellow lines obtained from trajectory simulations. They are, however, *not* equal to the correct long-time limit $\u27e8A\u0302\u27e9qc\u27e8B\u0302\u27e9qc$ (shown as grey dashed–dotted lines, calculated as described in Appendix B). The error is observed to become more severe for higher values of *ɛ* and *η*, and for the most extreme system (*ɛ* = 8, *η* = 4), the long-time limit of spin-LSC is even found to be larger than 1. This corresponds to an unphysical situation where the higher-energy state has a negative thermal population. We remark that the failure to relax to the correct limit (especially for strongly biased systems) is a problem not just for spin-LSC, but is present also in several other mapping techniques including Ehrenfest mean-field methods^{75,79} and approaches based on the MMST mapping,^{28,80,81} including SQC.^{13,82,83}

Next, we consider optimizing the spin radius in spin-LSC with respect to the long-time population (see Appendix C 2 for details). This idea is closely connected to the approach by Müller and Stock of optimizing the zero-point energy parameter in the MMST mapping. The analogous strategy in the spin-mapping framework is to choose the sphere radius *g* such that $K33g(t)$ and $\u27e8\sigma 3(t)\u27e9ccg$ relax to the correct long-time limits. We include results for this approach in the inset on the top right corner of Fig. 4, for a representative system with *ɛ* = 4 and *η* = 1. The long-time limit of *K*_{33}(*t*) is correct by construction. However, a caveat of the method is that the system is formally sampled from an out-of-equilibrium initial distribution (relative to the Hamiltonian that generates the dynamics). This implies that the time-dependent average $\u27e8\sigma 3(t)\u27e9cc$ (shown with light blue line in the lower panel of the inset) is not preserved by the dynamics, which breaks the desired time-translational invariance. The optimized-sphere approach is also limited to fixing a single expectation value and will, in general, lead to an incorrect long-time limit of $\u27e8\sigma 1(t)\u27e9cc$ [and hence also of *K*_{31}(*t*) and *K*_{13}(*t*), while *K*_{21}(*t*) is anyway guaranteed to relax to the correct zero long-time limit by symmetry in this case].

In contrast, the ellipsoid mapping approach recovers all the correct long-time limits by construction, as shown with green and red lines for the global and local constructions, respectively. Because these methods preserve their respective distributions and obey Liouville’s theorem, it is possible to obtain converged results with a low number of trajectories by using the time-averaging procedure^{1}

For these models, we found that the order of 10^{3} trajectories of length 2*τ*_{max} = 10^{3} was sufficient for convergence. Note, however, that while the ellipsoid methods relax by construction to the correct long-time limits, the intermediate dynamics are found to be less accurate than with spin-LSC or the optimized sphere for large *η* and *ɛ*. In the models considered here, the global construction is found to overestimate and the local construction to underestimate the timescale of population transfer.

Finally, we remark that all methods are able to capture the transition from coherent to incoherent relaxation, which is seen to occur when increasing the electronic–nuclear coupling constant *η*. This behavior was first studied using real-time path-integral Monte Carlo,^{63} and it is interesting to note that it can be well described by simpler classical-trajectory methods.

### A. Detailed balance in rate theory

Having identified the local ellipsoid as a method that fulfills all the required properties in Sec. II E, we next discuss to what extent it is useful for calculating nonadiabatic rates. In particular, we show that the ellipsoid approach formally obeys detailed balance as typically defined in a rate-theory framework,^{40}

For the sake of the present discussion, we identify the reactants with a donor (d) and the products with an acceptor (a). Then $Pd=|d\u3009\u3008d|qc$ is the equilibrium probability of finding the system in the donor, and similarly *P*_{a} = 1 − *P*_{d}. These probabilities are multiplied by the rate constants *k*_{d→a} and *k*_{a→d} to obtain the forward and backward reaction rates, which according to Eq. (52) must be equal in equilibrium.

First, we consider identifying the donor |d⟩ and acceptor |a⟩ with the two electronic states such that

These definitions lead to the thermal side–side correlation functions^{84}

Based on these, the rate constant for the transition between |d⟩ and |a⟩ is defined by^{40}

and vice versa for $ka\u2192d(s)$ (where “s” stands for “side”). Here, *t*_{p} is the time required for the derivative $K\u0307da(t)$ to relax to a plateau after an initial transient on a shorter timescale.

Due to time-translational invariance, quantum correlation functions obey *K*_{03}(*t*) = *K*_{30}(*t*) = ⟨*σ*_{3}⟩, such that *K*_{ad}(*t*) = *K*_{da}(*t*) and hence detailed balance [Eq. (52)] is obeyed. Spin-LSC and the optimized-sphere method do not fulfill this property: although $K30(t)=\u27e8\sigma 3\u27e9qc$, *K*_{03}(*t*) is not correct at intermediate times in these methods (as demonstrated in the lower panel of the inset of Fig. 4), and thus within these approaches *K*_{ad}(*t*) ≠ *K*_{da}(*t*), meaning that the forward and backward reaction rates differ and Eq. (52) does not hold. Ellipsoid mapping, on the other hand, being time-translationally invariant by construction, fulfills $K03(t)=K30(t)=\u27e8\sigma 3\u27e9qc$ and *K*_{ad}(*t*) = *K*_{da}(*t*); hence, it obeys detailed balance [Eq. (52)].

Note, however, that there are alternative ways to compute the rate in practice. Using $ddt\sigma \u03023=2\Delta \sigma \u03022$, one can write the rate constant as an integral of the electronic (“e”) flux–flux correlation function,^{84}

Another definition of the rate follows from using the nuclear side operator *θ*(*y* − *y*^{‡}), where *y* denotes a nuclear reaction coordinate and *y*^{‡} a dividing surface between donor and acceptor. This leads to the expression^{40}

where

is the nuclear (“n”) flux–side correlation function and $Pd(n)=\u27e8\theta (y\u2021\u2212y)\u27e9qc$. Both of these alternative definitions trivially obey the detailed-balance condition [Eq. (52)] regardless of the underlying dynamics.

One would like the rates calculated from the three approaches to be identical, but in the following, we show that they are not. In particular, although $kd\u2192a(n)=kd\u2192a(s)$ (because the rate is formally independent of the choice of dividing surface),^{85} $kd\u2192a(e)$ differs from the others. As a practical example, we study a symmetric spin–boson model (*ɛ* = 0) with a Brownian oscillator spectral density, $J(\omega )=\Lambda 2\gamma \Omega 2\omega (\omega 2\u2212\Omega 2)2+\gamma 2\omega 2$, where Λ = 60 and *βγ* = *β*Ω = 0.5.^{86} For this model, we use *y* as defined in Eq. (48) and *y*^{‡} = 0. The results for *k*^{(n)} and *k*^{(e)} are shown as a function of Δ in Fig. 5. We have verified numerically that *k*^{(s)} tends towards *k*^{(n)} across the whole range of parameters studied, but these results require a longer time to plateau and are not shown here.

In the adiabatic limit (large Δ), *k*^{(e)} is found to agree with *k*^{(n)} = *k*^{(s)} and with the exact benchmark. In this regime, the ellipsoid shrinks to a “pancake”-like shape with all population on the lower adiabat leading to Born–Oppenheimer dynamics. In the nonadiabatic limit (small Δ), the electronic measure of the flux operator leads to rate constants in good agreement with Marcus theory, whereas the nuclear-flux approach fails to reproduce the correct result. Although it appears promising that the rate based on the electronic flux is accurate in this case, it is clear that the underlying nuclear dynamics have the wrong physical behavior.

Next, we present the reason for why there is a difference between the formulation based on the electronic fluxes (*k*^{(e)}) and that based on the electronic side operators (*k*^{(s)}). This can be understood by analyzing the time derivative of ** σ** [Eq. (43)],

The first term describes rotation under a fixed ellipsoid shape, while the last two terms describe reshaping and translation of the ellipsoid. Due to the presence of these terms, the time derivative of the electronic side operator is not equal to the electronic flux operator, since $\sigma \u03073\u22602\Delta \sigma 2$. Hence, the time derivative $K\u0307da$ is not equal to the electronic flux–side correlation function (and $K\u0308da$ is not equal to the electronic flux–flux correlation function). The rate Eq. (57) obtained from the nuclear flux, however, is equivalent to Eq. (55) obtained from $K\u0307da$ as the rate is formally independent of the choice of dividing surface within the extended mapping phase space. This explains why the rate based on the electronic flux is different from that based on the nuclear flux.

Finally, we consider why the rate based on the nuclear flux approach appears to be less accurate than that based on the electronic flux. The nuclear flux reports directly on the dynamics of the nuclei, whereas the electronic flux is relatively insensitive to them. It is, therefore, clear that the nuclei are not behaving as expected. In the small-Δ limit, away from the diabatic crossing, the ellipsoid quickly reduces to the lower adiabat due to the large reorganization energy and in this way forces too much population transfer to occur. This is a consequence of our choice of equations of motion, in which the ellipsoid parameters are determined solely by *x* and are independent of the current value of ** u**. As was stated in Sec. II E 2, the equations of motions used here are not unique but could possibly be adapted to correct this issue in the future. We also remark that other nonadiabatic trajectory methods such as surface hopping do not generally recover Marcus theory without additional decoherence corrections

^{87}and that this may provide inspiration for an alternative solution.

## V. OUTLOOK TO NONADIABATIC PATH-INTEGRAL METHODS

In this paper, we have investigated equilibrium dynamics of a mixed quantum–classical system in which the nuclei are treated with classical statistics. This neglects nuclear zero-point energy and tunneling, which may become important for large values of *β**ω*_{c}. A range of methods have been proposed to include such effects in nonadiabatic dynamics by means of imaginary-time path-integral techniques. However, it has turned out to be difficult to make such methods obey detailed balance (according to the full quantum Boltzmann distribution), without breaking other important limits such as Rabi oscillations.^{88} We can compare these to the methods developed in the present paper by considering the high-temperature (single-bead) limit, where nuclear quantum effects are unimportant.

Perhaps the simplest way to formulate nonadiabatic path-integral dynamics is using a path integral over the nuclear degrees of freedom but a discrete sum over the electronic degrees of freedom. This approach is sometimes referred to as mean-field RPMD.^{89,90} In its single-bead limit, the method would evolve nuclei on the mean-field potential defined by $\u22121\beta logTrq[e\u2212\beta V\u0302(x)]$. Given that it does not provide any real-time information of electronic operators, it cannot capture the rate in the nonadiabatic limit (similar to those obtained with the nuclear flux operator in Fig. 5). A more elaborate mean-field method has been developed by Montoya-Castillo and Reichman,^{76} based on a path-integral simulation of the Wigner representation of the quantum Boltzmann distribution. Their approach propagates the trajectories with Ehrenfest dynamics, which does not preserve the equilibrium distribution or guarantee that correlation functions will relax to the correct long-time limits.

Another type of method represents the electronic as well as the nuclear degrees of freedom with path integrals using the MMST or spin mapping.^{15,17,18,91} Since these approaches use a different Hamiltonian for dynamics and for statistics, they also do not preserve the equilibrium distribution even in the high-temperature limit. Although one formulation of RPMD derived in Ref. 16 does preserve the distribution, it does so at the expense of breaking ergodicity. In particular, the distribution is not positive definite and trajectories are unable to cross between positive and negative regions of the phase space. Importantly, the same method fails to reproduce the correct Rabi oscillations in the limit of an isolated subsystem.^{16} The distribution of ellipsoid mapping, on the other hand, is positive definite and conserved by the dynamics. In addition, the ellipsoid method is capable of reproducing the correct Rabi oscillations in the limit of zero system–bath coupling.

Finally, an approach known as isomorphic-RPMD^{19} converts the path-integral problem into an isomorphic nonadiabatic system in which classical Boltzmann sampling over the effective Hamiltonian yields the correct statistics. The approach makes use of standard methods such as surface hopping^{92,93} to solve the dynamics. This method would only obey detailed balance provided that the underlying classical nonadiabatic dynamics does, which in general is not the case for surface hopping. In principle, one might attempt to combine the isomorphic-RPMD approach with ellipsoid mapping in order to fulfill detailed balance according to the full quantum Boltzmann distribution. However, even if the ellipsoid method could be improved to give consistent rate constants, isomorphic-RPMD cannot capture tunneling in the golden-rule limit regardless of the choice of the underlying dynamics.^{94} In conclusion, the search for a nonadiabatic version of RPMD is still an open question.

## VI. CONCLUSIONS

In this article, we have investigated the problem of detailed balance in mixed quantum–classical systems and explored new ways to calculate thermal equilibrium correlation functions. The direct extension of spin-LSC to thermal correlation functions as well as its optimized-sphere generalization are both found to break detailed balance in general. The global ellipsoid preserves electronic expectation values but does not sample the nuclear phase space correctly. These issues are solved by the local ellipsoid, which rigorously conserves the quantum–classical equilibrium distribution and relaxes to the correct long-time limits, while predicting the correct Rabi oscillations in the limit of zero electron–nuclear coupling. These properties are not, in general, simultaneously fulfilled by any other classical-trajectory method that we are aware of. By smoothly reducing *g*_{ij} for energetically separated states, the local ellipsoid prevents trajectories from running off on inverted potentials and is therefore expected to avoid the difficulties for anharmonic systems experienced by previous mappings.

Although the local ellipsoid is a solution to the long sought after method that simultaneously obeys detailed balance as well as recovers Rabi oscillations,^{88} we found that these properties are not sufficient to give the correct relaxation timescales in the golden-rule limit. Specifically, we observed that $\sigma \u03073\u22602\Delta \sigma 2$, meaning that the time derivative of the electronic side operator is different from the electronic flux operator. We noted in Sec. II E 2 that the equations of motion were not unique, and it may be possible to adapt them to fulfill further properties than the ones listed in Sec. II E. To outline how the dynamics may be improved, we point out that the ellipsoid approach can (like other mappings) be made formally exact by taking the phase-space representation of the full propagator. In the classical-nuclear limit, a reasonable starting point is the QCLE,

which in the ellipsoid approach takes the form

To allow a solution in terms of independent trajectories, this needs to be approximated as

There are many ways in which to perform this approximation. Our hope is that in the future an approximation can be found that fulfills $\sigma \u0307=2H\xd7\sigma $ such that the time derivative of the electronic side operator becomes the electronic flux operator. This may require constraints or other specialized techniques in order to not break the condition |** u**|

^{2}= 1. In any case, the ellipsoid approach presented here constitutes a rigorous framework on which further improvements can be based.

An alternative development could be to process the current ellipsoid dynamics with the generalized quantum master equation (GQME), as has recently been done to improve the accuracy of nonequilibrium spin-mapping methods.^{54} The previous proof that GQME cannot improve the results of methods that obey detailed balance^{95} does not hold in this case because it would still correct for the issue that $\sigma \u03073\u22602\Delta \sigma 2$ in the current dynamics.

Although in the present work we focused on a two-state model for simplicity, we expect that the approach could be extended to more than two levels similarly to the original (spherical) spin mapping.^{32} Finally, we restricted our analysis to purely classical nuclear dynamics neglecting nuclear zero-point energy and tunneling, which is justified at high temperatures and low nuclear frequencies. It may be possible to include such effects in the future by combining it with path-integral methods, such as RPMD.

## SUPPLEMENTARY MATERIAL

See the supplementary material for numerical data points of the results shown in Figs. 3 and 4.

## ACKNOWLEDGEMENTS

This project has received funding from European Union’s Horizon 2020 under MCSA Grant No. 801459, FP-RESOMUS.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

G.A. and J.E.R. contributed equally to this work.

**Graziano Amati:** Formal analysis and data curation (lead of HEOM, thermal spin-LSC, optimized sphere in Appendix C 2, thermal averages), writing – original draft (equal), writing – review & editing (equal). **Johan E. Runeson:** Formal analysis and data curation (lead of optimized sphere in Sec. II B, SW generalization, ellipsoids, rates), writing – original draft (equal), writing – review & editing (equal). **Jeremy O. Richardson:** Formal analysis (supporting), writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

### APPENDIX A: PRINCIPAL-AXIS BASIS

In this section, we show how to obtain the ellipsoid parameters $g$, ** c**, and $H\u0303$ in the isolated subsystem and in the local construction for mixed quantum–classical systems. For simplicity, we assume that the Hamiltonian is chosen to be real (i.e.,

*H*

_{2}= 0).

Rotating to the principal-axis basis corresponds to diagonalizing the Hamiltonian, $H\u0302=U\u0302H\u0302pU\u0302\u2020$, where

The diagonal elements of $H\u0302p$ are the adiabatic energies, sorted in decreasing order if the *z*^{p}-axis is chosen to be aligned with the field, as in Fig. 2. In this orientation, $H3p=12([H\u0302p]11\u2212[H\u0302p]22)$ is positive. However, we shall present the following solution so that it is valid for both positive and negative $H3p$.

We define the rotation matrix $R$ such that $U\u0302\u2020\sigma \u0302iU\u0302=Rij\sigma \u0302j$. Hence, $Hip=HjRji$ and $Rij=12Trq[U\u0302\u2020\sigma \u0302iU\u0302\sigma \u0302j]$, or more explicitly

Given that scalar quantities are basis-independent, $H\u22bagu=(Hp)\u22bagpup$ and ** c** ·

**=**

*u*

*c*^{p}·

*u*^{p}, which is fulfilled by the transformation rules $u=R\u22baup$, $g=R\u22bagpR$, $c=R\u22bacp$, and $H\u0303=H\u0303p$.

In the principal-axis basis, the only non-zero expectation value is $\u27e8\sigma \u03023\u27e9q=\u2212tanh\zeta $. The corresponding phase-space average is

where $L(x)=cothx\u22121x$ is the Langevin function. Equating the expectation values leads to an expression for $c3p$ in terms of $g33p$,

Note that reversing the sign of *ζ* (i.e., flipping the direction of the principal *z* axis) also reverses the sign of $c3p$ (assuming that $g33p>0$).

Next, we match the zero-time correlations. Recall that there is a circular symmetry about the field, so that $g11p=g22p$. By evaluating the traces and integrals analytically, the equation $\u27e8\sigma \u03021,\sigma \u03021\u27e9q=\u27e8\sigma 1,\sigma 1\u27e9c$ reduces to

which means that $g11p$ can also be expressed in terms of $g33p$. The only remaining unknown is therefore $g33p$, which we obtain by solving $\u27e8\sigma \u03023,\sigma \u03023\u27e9q=\u27e8\sigma 3,\sigma 3\u27e9c$, or

This is a single transcendental equation which we solve numerically with a one-dimensional root search. It has a unique positive solution that is shown in Fig. 3.

What remains is to determine $H\u0303p=H\u0303$. This is done by matching the partition functions

so that

The right-hand side is positive for *ζ* > 0 and negative for *ζ* < 0. Because *ζ* and $H3p$ have the same sign, it follows that $H\u0303$ is always positive.

### APPENDIX B: CALCULATION OF QUANTUM–CLASSICAL THERMAL AVERAGES

Here, we discuss how to efficiently calculate the quantum–classical thermal averages on the left-hand side of Eqs. (31b) and (31c), for the specific case of the spin–boson model.

To reduce the dimensionality of the phase-space averages, we note from Eq. (45c) that the system and the bath are coupled only via a scalar nuclear reaction coordinate *y* as defined in Eq. (48). All other *F* − 1 nuclear degrees of freedom can be identified as a secondary bath, which can be analytically integrated out in the calculation of static correlation functions.

The bath and system–bath coupling terms in the spin–boson Hamiltonian are rewritten as a function of *y* as^{65,96}

where

Here, *p*_{y} is the momentum conjugate to *y*, while all variables related to the secondary bath are denoted with a tilde. The coupling constants $c\u0303\alpha $ and frequencies $\omega \u0303\alpha $ can be calculated from the relation between the spectral densities of the primary and secondary bath.^{62} The quantum–classical partition function in the new variables becomes

and analogous expressions can be written for $\u27e8\sigma \u0302i\u27e9qc$ and $\u27e8\sigma \u0302i,\sigma \u0302j\u27e9qc$. The reduced partition function

is obtained after integrating out in Eq. (B3) all nuclear variables apart from the reaction coordinate *y*. An unimportant prefactor from the integration of the secondary bath is omitted in Eq. (B4), as it cancels out in the normalization of phase-space averages of the electronic operators we are interested in. In addition, by defining

we can rewrite exact quantum–classical averages as

The trace in the above expressions is easily evaluated and the integral over *λ* can be performed analytically. Finally, the numerical solution of the integral over *y* can be calculated via one-dimensional quadrature.

### APPENDIX C: SPHERICAL SPIN MAPPING FOR THERMAL CORRELATION FUNCTION

Here, we discuss two spherical mapping approaches suited to the study of thermal correlation functions. A comparison between these two methods and ellipsoid mapping is discussed in Sec. IV.

#### 1. Spin-LSC

We introduce here a formulation of spin-LSC^{31,32} suited to the study of thermal correlation functions. The approach is based on the notion [already introduced in Eq. (13)] that the electronic trace of the product of two electronic operators $A\u0302$ and $B\u0302$ can be represented as an integral over a spherical phase space

where

and $gW=3$ in standard spin-LSC.^{31} Based on the mapping prescription in Eq. (C1), we define thermal expectation values and correlation functions within spin-LSC by

where

and $\sigma iW(u)$ has been defined in Eq. (12). At *t* = 0, $\u27e8\sigma j(0)\u27e9W$ and $KijW(0)$ are equal by construction to their quantum–classical correspondents $\u27e8\sigma \u0302i\u27e9qc$ and $\u27e8\sigma \u0302i,\sigma \u0302j\u27e9qc$, respectively. When the variables (*x*, *p*, ** u**) are evolved under the Hamiltonian $HW(x,p,u)=Trq[H\u0302w\u0302W(u)]$, the time evolution of $KijW(t)$ is an approximation to the exact dynamics of

*K*

_{ij}(

*t*) analogous to LSC-IVR.

To evaluate Eq. (C2), we sample (*x*, *p*) from the uncoupled bath distribution

and ** u** uniformly from the unit sphere. Since this initial distribution is effectively out of equilibrium [like Eq. (10)], spin-LSC will not conserve expectation values in time or lead to the correct long-time limits. If we assume that correlation functions decorrelate at long times, the expressions in Eq. (C2) will instead approach

where we defined expectation values with respect to the initial distribution,

with *Z*_{qc} as defined in Eq. (6), and with respect to the equilibrium distribution

Also, we used that

In general, $\u27e8\sigma j\u27e9eqW\u2260\u27e8\sigma \u0302j\u27e9qc$, and thus Eq. (C5) leads to the wrong long-time limits.

#### 2. Optimized sphere

The mapping prescriptions in Eq. (C1) can be generalized to *any* radius *g*_{s} according to

where

and where $w\u0302s\u0304$ is defined in the same way but with $gs\u0304=3/gs$. We can use the additional free variable to correct one of the equilibrium expectation values of spin-LSC. In this way, one can make the mapping approximations of both *K*_{ij}(*t*) and $\u27e8\sigma \u0302j(t)\u27e9qc$ relax to the correct long-time limits for a chosen value of *j* (and all *i* = 1, 2, 3). Here, the time dependence refers to evolution under *H*_{s}. In particular, we look for a numerical solution of *g*_{s} to the nonlinear equation

where $\u27e8\sigma j\u27e9eq(s)$ is defined as in Eq. (C2a), but with a radius *g*_{s}.

In the case of the spin–boson model, Eq. (C10) can be simplified by integrating out all degrees of freedom except for an electronic coordinate *z* = cos *θ*, defined such that $\sigma 3(z)=gsz$. In the specific case of *j* = 3 [with the aim of studying *K*_{33}(*t*)], we find

where we introduced the modified Bessel function of the first kind,^{97} *I*_{0}(*x*), and the reorganization energy $\Lambda =\u2211\alpha =1F2c\alpha 2m\alpha \omega \alpha 2$. The *z* integrals were evaluated numerically.

### APPENDIX D: BENCHMARK KUBO-TRANSFORMED THERMAL CORRELATION FUNCTIONS

To assess the accuracy of our quantum–classical calculations, we used a numerically exact solution of the HEOM^{98} calculated with the open-source pyrho package.^{99} Even though HEOM treats both the electronic and the nuclear degrees of freedom quantum-mechanically, the method is relevant for our quantum–classical framework by restricting to the high temperature regime. The solution of HEOM needs to converge over two parameters: the number of Matsubara modes, *K*, and the truncation level, *L*. Given that we restrict our analysis to high temperatures, it sufficed for us to fix *K* = 0. This regime corresponds to an exponentially fast decorrelation of the autocorrelation function of the nuclear reaction coordinate, as expected for a classical bath.^{100} For all the systems whose HEOM solution is shown in Fig. 4, the results converged with a truncation level between *L* = 13 and *L* = 19. For two of the most “extreme” systems (marked with the label “no HEOM” in the figure), the solution of the HEOM failed to converge in the parameter *L* before calculations became too demanding; hence, we could not include benchmark results for those two cases. This issue stresses the importance of the development of efficient classical-trajectory techniques that can tackle systems inaccessible by numerically exact methods.

To calculate Kubo-transformed correlation functions with HEOM, we first prepared the system in an arbitrary out-of-equilibrium density

where $\rho \u0302b$ denotes the quantum thermal distribution for the uncoupled bath. We then propagated the dynamics for an equilibration time *t*_{eq} ≫ 1, chosen such that each expectation value plateaus to a constant,

where $\u27e8\u22c5\u27e9qb$ denotes an average over the initial distribution Eq. (D1), and the right-hand side corresponds to the expected quantum–classical equilibrium average as defined in Eq. (8). The condition Eq. (D2) can be seen as a consistency test on whether the assumption of classical nuclei is valid for a given system. This means that after the time *t*_{eq}, the initial density has thermalized to $e\u2212iH\u0302t\rho \u03020eiH\u0302t=1Ze\u2212\beta H\u0302$.

The output of the HEOM simulations is the standard thermal correlation function *C*_{ij}(*t*) in Eq. (3), which was then Kubo-transformed using the Fourier-space identity^{4}

where $f\u0303(\omega )=\u222b\u2212\u221e+\u221edte\u2212i\omega tf(t)$. Finally, Eq. (D3) was transformed back to time domain to get *K*_{ij}(*t*).

## REFERENCES

Note that some of the present expressions differ by a factor 2 from previous work which used a basis of $S\u0302i=12\sigma \u0302i$; the form presented here is equivalent but leads to slightly neater equations.

For *N* levels, the optimal choice is the scaling $g=N+1$, but for three or more levels the error is of order $O(\beta 2)$ instead of $O(\beta 3)$. The two-level case is special because there the second-order contribution is identically zero due to the properties of the Pauli matrices.

Note that by setting *Û* = e^{−iĤt}, where *Ĥ* is the Hamiltonian of the subsystem, property 5 says that quantum time evolution should correspond to a rotation in the phase space of ** u**.