On detailed balance in nonadiabatic dynamics: From spin spheres to equilibrium ellipsoids

Amati, Graziano; Runeson, Johan E.; Richardson, Jeremy O.

doi:10.1063/5.0137828

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).¹ 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.² 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),³ ring-polymer molecular dynamics (RPMD),⁴ and path-integral Liouville dynamics.⁵ These have in many ways replaced previous methods that lacked this property, such as the linearized semiclassical-initial value representation (LSC-IVR).⁶ 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,⁹ 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)¹² obeys detailed balance in the limit of small diabatic coupling¹³ but not in general. Although simple rate descriptions, such as secular Redfield theory, do obey detailed balance by construction,¹⁴ 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.²⁹ 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,³⁰ 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- $\frac{1}{2}$ with radius $\frac{\sqrt{3}}{2}$ ⁠, 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_BT, 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_BT, reduce to adiabatic dynamics on the lower state. In the opposite limit of states close in energy compared to k_BT, 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

\hat{H} = \frac{{\hat{p}}^{2}}{2 m} + \hat{V} (\hat{x}), \hat{V} (x) = (\begin{matrix} V_{1} (x) & Δ^{*} (x) \\ Δ (x) & V_{2} (x) \end{matrix}),

(1)

where V₁ and V₂ are potential-energy surfaces and Δ is the coupling between the two levels. Typically, $\hat{x}$ and $\hat{p}$ 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.³²

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

⟨ \hat{A} ⟩ = T r [\hat{ρ} \hat{A}] = \frac{1}{Z} T r [e^{- β \hat{H}} \hat{A}],

(2)

where β = 1/k_BT, the partition function is $Z = T r [e^{- β \hat{H}}]$ ⁠, 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,⁴⁰ and scattering functions,⁴¹ can be expressed in terms of equilibrium time-correlation functions. The standard correlation function between two operators $\hat{A}$ and $\hat{B}$ is defined as

\begin{align} \begin{aligned} C_{A B} (t) & = \frac{1}{Z} T r [e^{- β \hat{H}} \hat{A} (0) \hat{B} (t)], \\ \hat{B} (t) & = e^{i \hat{H} t} \hat{B} (0) e^{- i \hat{H} t}, \end{aligned} \end{align}

(3)

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⁴² is not the standard one above, but the Kubo-transformed correlation function, defined by

K_{A B} (t) = \frac{1}{β Z} \int_{0}^{β} d λ T r [e^{- (β - λ) \hat{H}} \hat{A} (0) e^{- λ \hat{H}} \hat{B} (t)] .

(4)

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.⁴ 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.⁴³ 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_BT 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 $\hat{x}$ and $\hat{p}$ 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⁴⁴

\hat{ρ} (x, p) = \frac{1}{Z_{q c}} e^{- β \hat{H} (x, p)},

(5)

and the partition function

Z_{q c} = {T r}_{q c} [e^{- β \hat{H} (x, p)}],

(6)

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,

{T r}_{q c} [\hat{f}] = \int d x d p {T r}_{q} [\hat{f} (x, p)] .

(7)

Based on this prescription, expectation values are defined as

{⟨ \hat{A} ⟩}_{q c} = {T r}_{q c} [\hat{ρ} \hat{A}] = \frac{1}{Z_{q c}} \int d x d p {T r}_{q} [e^{- β \hat{H} (x, p)} \hat{A} (x, p)] .

(8)

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

K_{A B}^{(q c)} (0) = \frac{1}{Z_{q c}} \int d x d p \frac{1}{β} \int_{0}^{β} d λ {T r}_{q} [e^{- (β - λ) \hat{H}} \hat{A} e^{- λ \hat{H}} \hat{B}],

(9)

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⁴⁵ (although it is exact for the spin–boson model). One strategy to define a quantum–classical limit for $\hat{B} (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_BT, 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,

C_{A B}^{b} (t) = T r [{\hat{ρ}}_{b} \hat{A} (0) \hat{B} (t)],

(10)

where ${\hat{ρ}}_{b}$ 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.³¹

As is well known, a two-level system is isomorphic to a spin- $\frac{1}{2}$ in a magnetic field since the Hamiltonian is of the form

\hat{H} = H_{0} {\hat{σ}}_{0} + H_{1} {\hat{σ}}_{1} + H_{2} {\hat{σ}}_{2} + H_{3} {\hat{σ}}_{3} \equiv H_{μ} {\hat{σ}}_{μ},

(11)

where ${\hat{σ}}_{μ}$ are the Pauli matrices (including ${\hat{σ}}_{0}$ 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 $H_{0} = \frac{p^{2}}{2 m} + \frac{1}{2} (V_{1} + V_{2})$ ⁠, H₁ = Re Δ, H₂ = Im Δ, $H_{3} = \frac{1}{2} (V_{1} - V_{2})$ ⁠, and the corresponding “magnetic field” is the vector H = (H₁, H₂, H₃).⁴⁶ The operators $\hat{A}$ and $\hat{B}$ 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,⁴⁷ which has recently gained attention^48,49 after having been largely overlooked for many decades.⁵⁰ 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

\{\begin{cases} {\hat{σ}}_{0} \mapsto σ_{0}^{W} = 1, \\ {\hat{σ}}_{i} \mapsto σ_{i}^{W} = g_{W} u_{i} & i = 1, 2, 3, \end{cases})

(12)

where u = (u₁, u₂, u₃) = (sin θ cos φ, sin θ sin φ, cos θ) and $g_{W} = \sqrt{3}$ ⁠. Here, θ and φ are the usual Bloch-sphere angles. Then, $H_{W} = H_{μ} σ_{μ}^{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,

{T r}_{q} [{\hat{σ}}_{μ} {\hat{σ}}_{ν}] = 2 δ_{μ ν} = \int d u σ_{μ}^{W} σ_{ν}^{W},

(13)

where $d u = \frac{1}{2 π} d φ d θ \sin θ$ denotes integration over the unit sphere (normalized such that $\int d u = {T r}_{q} [{\hat{σ}}_{0}] = 2$ ⁠). Note that the factor g_W also appears in the familiar expression for the magnitude $\frac{\sqrt{3}}{2} = \sqrt{S (S + 1)}$ of a spin $S = \frac{1}{2}$ ⁠, 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

C_{A B}^{b} (t) \approx \int d x d p d u ρ_{b} (x, p) A_{W} (x, p, u) B_{W} (x_{t}, p_{t}, u_{t}),

(14)

where the time evolution is given by

\dot{x} = p / m,

(15a)

\dot{p} = - \frac{\partial}{\partial x} H_{W},

(15b)

\dot{u} = 2 H \times u,

(15c)

and H_W = H₀(x) + g_WH(x) ⋅ u. Equation (14) is exact in the case of an uncoupled subsystem. In general, however, this independent-trajectory treatment is only exact at 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^{- β \hat{H}}$ in Eq. (3) or (4), one might (naively) attempt to use the form $\frac{1}{Z_{W}} \int d x d p d u e^{- β H_{W}} A_{W} B_{W} (t)$ ⁠, where the symbol “W” refers to the SW representation of subsystem operators [Eq. (12)] and $Z_{W} = \int d x d p d u e^{- β H_{W}}$ ⁠. 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^{- β H_{W}} \neq {[e^{- β \hat{H}}]}_{W} = \sum_{n} \frac{1}{n!} {(- β)}^{n} {[{\hat{H}}^{n}]}_{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 $\frac{1}{β Z_{q c}} \int_{0}^{β} d λ e^{- (β - λ) \hat{H}} \hat{A} e^{- λ \hat{H}}$ in place of ${\hat{ρ}}_{b} \hat{A}$ 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^{- β H_{W}}$ [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_BT, 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.

FIG. 1.

View large Download slide

Equilibrium population of the high-energy level within different quasiclassical mapping formalisms as a function of the dimensionless parameter ζ = βH₃. Spin mapping is closer to the quantum result than Ehrenfest or MMST for high T (low ζ) but predicts negative populations for low T (large ζ). The reason is that for large ζ, the spin distribution localizes to a single direction (purple dots on the spheres) antiparallel to the field (top right corner), which is different from the correct ground state (red circles). In this paper, we generalize spin mapping to reproduce the quantum populations for any ζ. Figure adapted from Ref. 8.

Since it is always possible to find a basis in which H₁ = H₂ = 0 (and because the theory is basis-independent), it suffices to consider the case $\hat{H} = H_{3} {\hat{σ}}_{3}$ ⁠. Then, ${⟨ {\hat{σ}}_{1} ⟩}_{q} = {⟨ {\hat{σ}}_{2} ⟩}_{q} = 0$ ⁠, while the remaining expectation value is

{⟨ {\hat{σ}}_{3} ⟩}_{q} = \frac{{T r}_{q} [e^{- β H_{3} {\hat{σ}}_{3}} {\hat{σ}}_{3}]}{{T r}_{q} [e^{- β H_{3} {\hat{σ}}_{3}}]} = - \tanh ζ,

(16)

where ζ ≡ βH₃. This quantum result is to be compared with the corresponding classical expectation value that arises from the mapping ${\hat{σ}}_{3} \mapsto σ_{3} = g \cos θ$ ⁠, which is

{⟨ σ_{3} ⟩}_{c} = \frac{\int d u e^{- ζ g \cos θ} g \cos θ}{\int d u e^{- ζ g \cos θ}} = - g L (g ζ),

(17)

where $L (x) \equiv \coth x - \frac{1}{x}$ is the Langevin function (a limiting case of the Brillouin function, which appears frequently in the paramagnetic theory of spins⁵⁵). For $g = g_{W} = \sqrt{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 (ζ^{3})$ ⁠, a result unique to $g_{W} = \sqrt{3}$ ⁠, which is thus the optimal scaling in this regard.⁵⁶ For comparison, Ehrenfest (corresponding to a sphere with scaling g = 1) and the original MMST mapping (g = 2) both have an error of $O (ζ)$ ⁠. 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 $g_{W} = \sqrt{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,²⁹ 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.¹³ 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.⁵⁷

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:⁵⁰

Linearity: $\hat{A} \mapsto A (u)$ is a one-to-one linear map.
Reality:
$[{\hat{A}}^{†}] (u) = {[A (u)]}^{*} .$
(18a)
Normalization:
${T r}_{q} [\hat{A}] = \int d u A (u) .$
(18b)
Tracing:
${T r}_{q} [\hat{A} \hat{B}] = \int d u A (u) B (u) .$
(18c)
Covariance:
$[\hat{U} (G^{- 1}) \hat{A} \hat{U} (G)] (u) = A (G \cdot u),$
(18d)
where G ∈ SU(2) is a linear rotation operator and $\hat{U} (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 $\hat{B}$ is the identity). In other words, it connects the trace product of quantum operators, ${(\hat{A}, \hat{B})}_{q} \equiv {T r}_{q} [\hat{A} \hat{B}]$ ⁠, with the “classical” inner product of functions on a sphere, (A,B)_c ≡ ∫du A(u)B(u). 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

{⟨ \hat{A}, \hat{B} ⟩}_{q} \equiv \int_{0}^{1} d λ {T r}_{q} [{\hat{ρ}}^{1 - λ} \hat{A} {\hat{ρ}}^{λ} \hat{B}],

(19)

where $\hat{ρ} = \frac{1}{Z_{q}} e^{- β H_{i} {\hat{σ}}_{i}}$ and $Z_{q} = {T r}_{q} [e^{- β H_{i} {\hat{σ}}_{i}}]$ ⁠. The classical analog of this inner product is

{⟨ A, B ⟩}_{c} \equiv \int d u ρ (u) A (u) B (u),

(20)

where $ρ (u) = \frac{1}{Z_{c}} e^{- β H_{i} σ_{i} (u)}$ and $Z_{c} = \int d u e^{- β H_{i} σ_{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:
$\int d u ρ (u) A (u) = {T r}_{q} [\hat{ρ} \hat{A}] .$
(21a)
4′.
Preservation of inner products:
$\int d u ρ (u) A (u) B (u) = \int_{0}^{1} d λ {T r}_{q} [{\hat{ρ}}^{1 - λ} \hat{A} {\hat{ρ}}^{λ} \hat{B}] .$
(21b)

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 $\hat{B}$ is the identity. We leave property 5 unchanged, since it is the key to defining dynamics in phase space.⁵⁸ The two modified properties are not obeyed by the standard mapping [Eq. (12)]. The problem of finding a new mapping ${\hat{σ}}_{μ} \mapsto σ_{μ} (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

F_{q} (a) = \ln {T r}_{q} [e^{- (β H_{i} - a_{i}) {\hat{σ}}_{i}}],

(22a)

F_{c} (a) = \ln \int d u e^{- (β H_{i} - a_{i}) σ_{i} (u)} .

(22b)

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

{((\frac{\partial}{\partial a_{i}} F_{q})|}_{a \to 0} = {((\frac{\partial}{\partial a_{i}} F_{c})|}_{a \to 0},

(23a)

{((\frac{\partial^{2}}{\partial a_{i} \partial a_{j}} F_{q})|}_{a \to 0} = {((\frac{\partial^{2}}{\partial a_{i} \partial a_{j}} F_{c})|}_{a \to 0} .

(23b)

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 $\hat{A} \mapsto A (u)$ that incorporates the properties described in Sec. II C. Since one can decompose any Hermitian operator as $\hat{A} = A_{μ} {\hat{σ}}_{μ}$ ⁠, 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,⁵⁹ which says that any trace-preserving quantum operation of the Bloch vector can be expressed as an affine map

u \mapsto g u + c,

(24)

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 u, 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

{\hat{σ}}_{i} \mapsto σ_{i} = g_{i j} u_{j} + c_{j},

(25)

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, ${\hat{σ}}_{0} \mapsto 1$ ⁠, 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 $\tilde{H}$ in the mapping Hamiltonian,

\hat{H} \mapsto H = H_{i} σ_{i} + H_{0} + \tilde{H}

(26a)

= H_{i} g_{i j} u_{j} + H_{i} c_{i} + H_{0} + \tilde{H} .

(26b)

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).⁶⁰

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

Z_{q} \overset{!}{=} Z_{c},

(27a)

〈 {{\hat{σ}}_{i} 〉}_{q} \overset{!}{=} {〈 σ_{i} 〉}_{c},

(27b)

〈 {\hat{σ}}_{i}, {{\hat{σ}}_{j} 〉}_{q} \overset{!}{=} 〈 σ_{i}, {σ_{j} 〉}_{c},

(27c)

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

FIG. 2.

View large Download slide

Visualization of the ellipsoid for an isolated two-level system. The ellipsoid is centered at c and $σ = g u + c$ is a vector on the surface (where u is a vector on the unit sphere). It is easiest to construct the ellipsoid in the principal-axis basis defined by the direction of the magnetic field, H.

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 $\hat{H}$ ⁠.) In this basis, the only non-zero expectation values are ${⟨ {\hat{σ}}_{3} ⟩}_{q}$ and ${⟨ {\hat{σ}}_{i}, {\hat{σ}}_{j} ⟩}_{q}$ for i = j. To match the system of equations above, the same needs to be true for ${⟨ σ_{i} ⟩}_{c}$ and ${⟨ σ_{i}, σ_{j} ⟩}_{c}$ ⁠. 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 $g_{11}^{p} = g_{22}^{p}$ ⁠. These considerations lead to the simpler form

c^{p} = (\begin{matrix} 0 \\ 0 \\ c_{3}^{p} \end{matrix}), g^{p} = (\begin{matrix} g_{11}^{p} & 0 & 0 \\ 0 & g_{11}^{p} & 0 \\ 0 & 0 & g_{33}^{p} \end{matrix}),

(28)

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 (⁠ $\tilde{H}$ ⁠, $c_{3}^{p}$ ⁠, $g_{11}^{p}$ and $g_{33}^{p}$ ⁠). 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 $ζ = β H_{3}^{p}$ ⁠, 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, $c_{3}^{p}$ and $\tilde{H}$ vanish while $g_{33}^{p}$ and $g_{11}^{p}$ approach $\sqrt{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 $\tilde{H}$ is vanishingly small compared to $H_{3}^{p}$ ⁠, 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., $g_{33}^{p} < g_{11}^{p} = g_{22}^{p}$ ⁠), which changes smoothly along the transition between two-level and effectively one-level systems. The parameter $\tilde{H}$ does not influence the shape of the ellipsoid but can be thought of as a shift of the energy.

FIG. 3.

View large Download slide

Ellipsoid-mapping parameters in the principal-axis basis for ζ > 0. In the limit of high temperature or low energy separation (left side), the ellipsoid reduces back to the W-sphere, and in the opposite limit (right side), it deforms into a point that represents the adiabatic ground state. The insets depict the gradual transition between these two limits. For ζ < 0, the signs of $c_{3}^{p}$ and $\tilde{H} / H_{3}^{p}$ are reversed, i.e., $\tilde{H}$ remains positive.

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

\dot{σ} = 2 H \times σ

(29a)

or equivalently in terms of u as

\dot{u} = 2 H \times u .

(29b)

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

〈 {\hat{σ}}_{μ} (0), {\hat{σ}}_{ν} {(t) 〉}_{q} = 〈 σ_{μ} (0), σ_{ν} {(t) 〉}_{c} .

(30)

In particular, the Boltzmann distribution is conserved, ${⟨ σ_{i} (t) ⟩}_{c} = {⟨ σ_{i} ⟩}_{c} = {⟨ {\hat{σ}}_{i} ⟩}_{q} = {⟨ {\hat{σ}}_{i} (t) ⟩}_{q}$ ⁠.

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

Z_{q c} \overset{!}{=} Z_{c c},

(31a)

〈 {{\hat{σ}}_{i} 〉}_{q c} \overset{!}{=} {〈 σ_{i} 〉}_{c c},

(31b)

〈 {\hat{σ}}_{i}, {{\hat{σ}}_{j} 〉}_{q c} \overset{!}{=} 〈 σ_{i}, {σ_{j} 〉}_{c c}

(31c)

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

{⟨ f ⟩}_{c c} = \frac{1}{Z_{c c}} \int d x d p \int d u e^{- β H (x, p, u)} f (x, p, u),

(32)

where

Z_{c c} = \int d x d p \int d u e^{- β H (x, p, u)} .

(33)

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 $\tilde{H}$ 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 $\tilde{H} (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|² = 1.
Obey phase-space incompressibility¹
$\nabla_{u} \cdot \dot{u} + \frac{\partial}{\partial x} \dot{x} + \frac{\partial}{\partial p} \dot{p} = 0$
(34)
such that Liouville’s theorem is valid.
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
$\frac{1}{Z_{q c}} {T r}_{q} [e^{- β \hat{H} (x, p)}] .$
(35)

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 H leaves the surface of the ellipsoid, whereas precession of u about H (which does ensure |u|² = 1) does not preserve the mapping Hamiltonian

\hat{H} (x, p) \mapsto H_{0} (x, p) + H_{i} (x) σ_{i} (u) + \tilde{H} .

(36)

In order to preserve energy as well as remain on the phase space, we instead propose to let u precess around the direction of $g H$ ⁠. 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

\dot{u} = \frac{2}{g_{eff}} (g H) \times u,

(37)

where $g_{eff} = \frac{| g H |}{| 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,⁶¹ 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 $\nabla_{u} \cdot \dot{u} = 0$ ⁠). Note that the term $\tilde{H}$ 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 ${⟨ σ_{μ}, σ_{ν} ⟩}_{c c} = {⟨ {\hat{σ}}_{μ}, {\hat{σ}}_{ν} ⟩}_{q c}$ ⁠, there is not enough flexibility to ensure that the equilibrium distribution over the nuclear phase space, $\frac{1}{Z_{c c}} \int d u e^{- β 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), $\tilde{H} (x)$ are functions of nuclear configuration. Effectively, this means that each x defines a two-level system with a local potential matrix $\hat{V} (x)$ ⁠. To continue, define the local density $\hat{ρ} (x) = \frac{1}{Z_{q} (x)} e^{- β \hat{V} (x)}$ ⁠, where $Z_{q} (x) = {T r}_{q} [e^{- β \hat{V} (x)}]$ and the trace is taken only over the electronic subsystem (leaving x unaffected). Next, define the local quantities

{〈 {\hat{σ}}_{i} 〉}_{q} (x) = {T r}_{q} [\hat{ρ} (x) {\hat{σ}}_{i}],

(38a)

〈 {\hat{σ}}_{i}, {{\hat{σ}}_{j} 〉}_{q} (x) = \int_{0}^{1} d λ {T r}_{q} [{\hat{ρ}}^{1 - λ} (x) {\hat{σ}}_{i} {\hat{ρ}}^{λ} (x) {\hat{σ}}_{j}] .

(38b)

Analogously, define the phase-space density $ρ (x, u) = \frac{1}{Z_{c} (x)} e^{- β V (x, u)}$ ⁠, where Z_c(x) = ∫due^−βV(x,u) and the local potential is

V (x, u) = V_{0} (x) + H_{i} (x) σ_{i} (x, u) + \tilde{H} (x),

(39)

with $V_{0} = \frac{1}{2} (V_{1} (x) + V_{2} (x))$ ⁠. Finally, we define

{〈 σ_{i} 〉}_{c} (x) = \int d u ρ (x, u) σ_{i} (x, u),

(40a)

〈 σ_{i}, {σ_{j} 〉}_{c} (x) = \int d u ρ (x, u) σ_{i} (x, u) σ_{j} (x, u) .

(40b)

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

Z_{q} (x) \overset{!}{=} Z_{c} (x),

(41a)

{〈 {\hat{σ}}_{i} 〉}_{q} (x) \overset{!}{=} 〈 {σ_{i} 〉}_{c} (x),

(41b)

〈 {\hat{σ}}_{i}, {{\hat{σ}}_{j} 〉}_{q} (x) \overset{!}{=} {〈 σ}_{i}, {σ_{j} 〉}_{c} (x) .

(41c)

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:

\dot{x} = \frac{p}{m},

(42a)

\dot{p} = - \frac{\partial}{\partial x} V (x, u) .

(42b)

\dot{u} = 2 H \times u,

(42c)

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_BT, 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

\dot{σ} = g \dot{u} + \dot{g} u + \dot{c} = 2 H \times σ + \dot{x} (\frac{\partial g}{\partial x} g^{- 1} (σ - c) + \frac{\partial c}{\partial x}) .

(43)

In practice, we sample and evolve the u variables and only convert to σ when calculating observables. The rotations of u were carried out in the standard way using quaternions.

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,

\hat{H} = {\hat{H}}_{s} + H_{b} {\hat{σ}}_{0} + {\hat{H}}_{s b},

(44)

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

{\hat{H}}_{s} = Δ {\hat{σ}}_{1} + ε {\hat{σ}}_{3},

(45a)

H_{b} (x, p) = \frac{1}{2} \sum_{α = 1}^{F} (p_{α}^{2} + m_{α} ω_{α}^{2} x_{α}^{2}),

(45b)

{\hat{H}}_{s b} (x) = {\hat{σ}}_{3} \sum_{α = 1}^{F} c_{α} x_{α} .

(45c)

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

J_{F} (ω) = \frac{π}{2} \sum_{α = 1}^{F} \frac{c_{α}^{2}}{m_{α} ω_{α}} δ (ω - ω_{α}),

(46)

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

J (ω) = η \frac{ω ω_{c}}{ω^{2} + ω_{c}^{2}},

(47)

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.⁷⁷

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.⁷⁸ 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

y = \sum_{α = 1}^{F} c_{α} x_{α} .

(48)

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₃₃(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₃₃(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

\lim_{t \to + \infty} K_{A B}^{W} (t) = 〈 {\hat{A} 〉}_{q c} 〈 {B 〉}_{eq}^{W},

(49)

where

{⟨ B ⟩}_{eq}^{W} \equiv \frac{\int d x d p d u e^{- β H_{W}} B}{\int d x d p d u e^{- β H_{W}}} .

(50)

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 ${⟨ \hat{A} ⟩}_{q c} {⟨ \hat{B} ⟩}_{q c}$ (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

FIG. 4.

View large Download slide

Thermal autocorrelation function of ${\hat{σ}}_{3}$ using the global and local ellipsoid (green and red, respectively), spin-LSC (yellow), and fully quantum HEOM results (black). Triangles indicate the theoretical predictions of the long-time limit of spin-LSC from Eq. (49), while the dashed–dotted light grey lines show the correct quantum–classical long-time limits. Inset: expanded results for ɛ = 4 and η = 1, including the generalization of spin mapping to an optimized sphere (light blue). The lower panel of the inset shows the time evolution of ⟨σ₃(t)⟩, which is expected to be constant provided detailed balance is obeyed. Shaded regions for the global and local ellipsoid indicate the 95% confidence interval for the statistical average. Note that time-averaging procedure [Eq. (51a)] was not used in this case, since it is not formally valid for the spherical methods, which do not obey detailed balance. In the case of the ellipsoid methods, it is valid, but we choose not to use it here in order to provide numerical evidence of this assertion.

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 $K_{33}^{g} (t)$ and ${⟨ σ_{3} (t) ⟩}_{c c}^{g}$ 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₃₃(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 ${⟨ σ_{3} (t) ⟩}_{c c}$ (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 ${⟨ σ_{1} (t) ⟩}_{c c}$ [and hence also of K₃₁(t) and K₁₃(t), while K₂₁(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¹

〈 A (0) B {(t) 〉}_{c c} = \frac{1}{τ_{\max}} \int_{0}^{τ_{\max}} d τ 〈 A (τ) B {(t + τ) 〉}_{c c} .

(51a)

For these models, we found that the order of 10³ trajectories of length 2τ_max = 10³ 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,⁶³ 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,⁴⁰

P_{d} k_{d \to a} = P_{a} k_{a \to d} .

(52)

For the sake of the present discussion, we identify the reactants with a donor (d) and the products with an acceptor (a). Then $P_{d} = {⟨| d 〉〈 d |⟩}_{q c}$ 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

| d 〉 〈 d | = \frac{1}{2} ({\hat{σ}}_{0} + {\hat{σ}}_{3}), | a 〉 〈 a | = \frac{1}{2} ({\hat{σ}}_{0} - {\hat{σ}}_{3}) .

(53)

These definitions lead to the thermal side–side correlation functions⁸⁴

K_{d a} (t) = \frac{1}{4} (K_{00} (t) - K_{03} (t) + K_{30} (t) - K_{33} (t)),

(54a)

K_{a d} (t) = \frac{1}{4} (K_{00} (t) + K_{03} (t) - K_{30} (t) - K_{33} (t)) .

(54b)

Based on these, the rate constant for the transition between |d⟩ and |a⟩ is defined by⁴⁰

k_{d \to a}^{(s)} = \frac{1}{P_{d}} \lim_{t \to t_{p}} {\dot{K}}_{d a} (t)

(55)

and vice versa for $k_{a \to d}^{(s)}$ (where “s” stands for “side”). Here, t_p is the time required for the derivative ${\dot{K}}_{d a} (t)$ to relax to a plateau after an initial transient on a shorter timescale.

Due to time-translational invariance, quantum correlation functions obey K₀₃(t) = K₃₀(t) = ⟨σ₃⟩, 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 $K_{30} (t) = {⟨ σ_{3} ⟩}_{q c}$ ⁠, K₀₃(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 $K_{03} (t) = K_{30} (t) = {⟨ σ_{3} ⟩}_{q c}$ 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 $\frac{d}{d t} {\hat{σ}}_{3} = 2 Δ {\hat{σ}}_{2}$ ⁠, one can write the rate constant as an integral of the electronic (“e”) flux–flux correlation function,⁸⁴

k_{d \to a}^{(e)} = \frac{Δ^{2}}{P_{d}} \int_{0}^{t_{p}} d t K_{22} (t) .

(56)

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⁴⁰

k_{d \to a}^{(n)} = \frac{1}{P_{d}^{(n)}} \lim_{t \to t_{p}} K_{fs}^{(n)} (t),

(57)

where

K_{fs}^{(n)} (t) = {⟨ \dot{y} δ (y - y^{‡}) θ (y_{t} - y^{‡}) ⟩}_{q c}

(58)

is the nuclear (“n”) flux–side correlation function and $P_{d}^{(n)} = {⟨ θ (y^{‡} - y) ⟩}_{q c}$ ⁠. 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 $k_{d \to a}^{(n)} = k_{d \to a}^{(s)}$ (because the rate is formally independent of the choice of dividing surface),⁸⁵ $k_{d \to a}^{(e)}$ differs from the others. As a practical example, we study a symmetric spin–boson model (ɛ = 0) with a Brownian oscillator spectral density, $J (ω) = \frac{Λ}{2} \frac{γ Ω^{2} ω}{{(ω^{2} - Ω^{2})}^{2} + γ^{2} ω^{2}}$ ⁠, where Λ = 60 and βγ = βΩ = 0.5.⁸⁶ For this model, we use y as defined in Eq. (48) and y^‡ = 0. The results for k⁽ⁿ⁾ and k^(e) are shown as a function of Δ in Fig. 5. We have verified numerically that k^(s) tends towards k⁽ⁿ⁾ across the whole range of parameters studied, but these results require a longer time to plateau and are not shown here.

FIG. 5.

View large Download slide

Reaction rates as a function of βΔ for a spin–boson model with Λ = 60 and ɛ = 0, comparing two different measurements using the local ellipsoid to numerically exact HEOM results from Ref. 86. If the rates are calculated from the plateau value of the nuclear flux–side correlation function (orange squares) one obtains a trend akin to mean-field dynamics (valid only for large Δ). Conversely, if the rates are calculated from the integral of the electronic flux–flux correlation function, one obtains correct behavior, including in the Marcus theory limit for small Δ (dashed line).

In the adiabatic limit (large Δ), k^(e) is found to agree with k⁽ⁿ⁾ = 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)],

\dot{σ} = g \dot{u} + \dot{g} u + \dot{c} .

(59)

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 ${\dot{σ}}_{3} \neq 2 Δ σ_{2}$ ⁠. Hence, the time derivative ${\dot{K}}_{da}$ is not equal to the electronic flux–side correlation function (and ${\ddot{K}}_{da}$ 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 ${\dot{K}}_{da}$ 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⁸⁷ 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.⁸⁸ 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 $- \frac{1}{β} \log {T r}_{q} [e^{- β \hat{V} (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,⁷⁶ 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.¹⁶ 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¹⁹ 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.⁹⁴ 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,⁸⁸ we found that these properties are not sufficient to give the correct relaxation timescales in the golden-rule limit. Specifically, we observed that ${\dot{σ}}_{3} \neq 2 Δ σ_{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,

\frac{d \hat{B}}{d t} = 2 ϵ_{ijk} H_{i} B_{j} {\hat{σ}}_{k} + \frac{p}{m} \frac{\partial \hat{B}}{\partial x} - \frac{1}{2} {[\frac{\partial \hat{H}}{\partial x}, \frac{\partial \hat{B}}{\partial p}]}_{+},

(60)

which in the ellipsoid approach takes the form

\frac{d B}{d t} = 2 ϵ_{ijk} H_{i} B_{j} (g_{k l} u_{l} + c_{k}) + \frac{p}{m} \frac{\partial B}{\partial x} - \frac{\partial H_{μ}}{\partial x} \frac{\partial B_{μ}}{\partial p} .

(61)

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

\frac{d B}{d t} \approx {\dot{u}}_{i} \frac{\partial B}{\partial u_{i}} + \dot{x} \frac{\partial B}{\partial x} + \dot{p} \frac{\partial B}{\partial p} .

(62)

There are many ways in which to perform this approximation. Our hope is that in the future an approximation can be found that fulfills $\dot{σ} = 2 H \times σ$ 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|² = 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.⁵⁴ The previous proof that GQME cannot improve the results of methods that obey detailed balance⁹⁵ does not hold in this case because it would still correct for the issue that ${\dot{σ}}_{3} \neq 2 Δ σ_{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.³² 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 $\tilde{H}$ 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₂ = 0).

Rotating to the principal-axis basis corresponds to diagonalizing the Hamiltonian, $\hat{H} = \hat{U} {\hat{H}}^{p} {\hat{U}}^{†}$ ⁠, where

\hat{U} = (\begin{matrix} \cos ϑ & - \sin ϑ \\ \sin ϑ & \cos ϑ \end{matrix}) .

(A1)

The diagonal elements of ${\hat{H}}^{p}$ 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, $H_{3}^{p} = \frac{1}{2} ({[{\hat{H}}^{p}]}_{11} - {[{\hat{H}}^{p}]}_{22})$ is positive. However, we shall present the following solution so that it is valid for both positive and negative $H_{3}^{p}$ ⁠.

We define the rotation matrix $R$ such that ${\hat{U}}^{†} {\hat{σ}}_{i} \hat{U} = R_{i j} {\hat{σ}}_{j}$ ⁠. Hence, $H_{i}^{p} = H_{j} R_{j i}$ and $R_{i j} = \frac{1}{2} {T r}_{q} [{\hat{U}}^{†} {\hat{σ}}_{i} \hat{U} {\hat{σ}}_{j}]$ ⁠, or more explicitly

R = (\begin{matrix} \cos 2 ϑ & 0 & - \sin 2 ϑ \\ 0 & 1 & 0 \\ \sin 2 ϑ & 0 & \cos 2 ϑ \end{matrix}) .

(A2)

Given that scalar quantities are basis-independent, $H^{⊺} g u = {(H^{p})}^{⊺} g^{p} u^{p}$ and c · u = c^p · u^p, which is fulfilled by the transformation rules $u = R^{⊺} u^{p}$ ⁠, $g = R^{⊺} g^{p} R$ ⁠, $c = R^{⊺} c^{p}$ ⁠, and $\tilde{H} = {\tilde{H}}^{p}$ ⁠.

In the principal-axis basis, the only non-zero expectation value is ${⟨ {\hat{σ}}_{3} ⟩}_{q} = - \tanh ζ$ ⁠. The corresponding phase-space average is

{⟨ σ_{3} ⟩}_{c} = \frac{\int d u e^{- β H_{3}^{p} g_{33}^{p} u_{3}} (g_{33}^{p} u_{3} + c_{3}^{p})}{\int d u e^{- β H_{3}^{p} g_{33}^{p} u_{3}}} = - g_{33}^{p} L (g_{33}^{p} ζ) + c_{3}^{p},

(A3)

where $L (x) = \coth x - \frac{1}{x}$ is the Langevin function. Equating the expectation values leads to an expression for $c_{3}^{p}$ in terms of $g_{33}^{p}$ ⁠,

c_{3}^{p} = g_{33}^{p} L (g_{33}^{p} ζ) - \tanh ζ .

(A4)

Note that reversing the sign of ζ (i.e., flipping the direction of the principal z axis) also reverses the sign of $c_{3}^{p}$ (assuming that $g_{33}^{p} > 0$ ⁠).

Next, we match the zero-time correlations. Recall that there is a circular symmetry about the field, so that $g_{11}^{p} = g_{22}^{p}$ ⁠. By evaluating the traces and integrals analytically, the equation ${⟨ {\hat{σ}}_{1}, {\hat{σ}}_{1} ⟩}_{q} = {⟨ σ_{1}, σ_{1} ⟩}_{c}$ reduces to

\frac{\tanh ζ}{ζ} = {(g_{11}^{p})}^{2} \frac{g_{33}^{p} ζ \coth (g_{33}^{p} ζ) - 1}{{(g_{33}^{p} ζ)}^{2}},

(A5)

which means that $g_{11}^{p}$ can also be expressed in terms of $g_{33}^{p}$ ⁠. The only remaining unknown is therefore $g_{33}^{p}$ ⁠, which we obtain by solving ${⟨ {\hat{σ}}_{3}, {\hat{σ}}_{3} ⟩}_{q} = {⟨ σ_{3}, σ_{3} ⟩}_{c}$ ⁠, or

{s e c h}^{2} ζ = \frac{1}{ζ^{2}} - {(g_{33}^{p})}^{2} {c s c h}^{2} (g_{33}^{p} ζ) .

(A6)

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 ${\tilde{H}}^{p} = \tilde{H}$ ⁠. This is done by matching the partition functions

Z_{q} = 2 e^{- β H_{0}} \cosh ζ,

(A7a)

Z_{c} = 2 e^{- c_{3}^{p} ζ - β (H_{0} + \tilde{H})} \frac{\sinh (g_{33}^{p} ζ)}{g_{33}^{p} ζ},

(A7b)

so that

\frac{\tilde{H}}{H_{3}^{p}} = - c_{3}^{p} - \frac{1}{ζ} \log \frac{g_{33}^{p} ζ \cosh ζ}{\sinh (g_{33}^{p} ζ)} .

(A8)

The right-hand side is positive for ζ > 0 and negative for ζ < 0. Because ζ and $H_{3}^{p}$ have the same sign, it follows that $\tilde{H}$ 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

H_{b} = \frac{1}{2} \sum_{α = 1}^{F - 1} [{\tilde{p}}_{α}^{2} + {\tilde{ω}}_{α}^{2} {({\tilde{x}}_{α} - \frac{{\tilde{c}}_{α} y}{{\tilde{ω}}_{α}^{2}})}^{2}],

(B1a)

{\hat{H}}_{s b} = κ y {\hat{σ}}_{3},

(B1b)

where

κ = {(\sum_{α = 1}^{F} c_{α}^{2})}^{1 / 2}, Ω^{2} = κ^{2} {(\sum_{α = 1}^{F} \frac{c_{α}^{2}}{ω_{α}^{2}})}^{- 1} .

(B2)

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 ${\tilde{c}}_{α}$ and frequencies ${\tilde{ω}}_{α}$ can be calculated from the relation between the spectral densities of the primary and secondary bath.⁶² The quantum–classical partition function in the new variables becomes

\begin{align} Z_{q c} & = \int d \tilde{x} d \tilde{p} d y d p_{y} \prod_{α = 1}^{F - 1} \exp \{- \frac{β}{2} [{\tilde{p}}_{α}^{2} + {\tilde{ω}}_{α}^{2} {({\tilde{x}}_{α} - \frac{{\tilde{c}}_{α} y}{{\tilde{ω}}_{α}^{2}})}^{2}]\} \\ \times \exp \{- β (\frac{p_{y}^{2}}{2} + Ω^{2} \frac{y^{2}}{2 κ^{2}})\} {T r}_{q} [e^{- β ({\hat{H}}_{s} + κ y {\hat{σ}}_{3})}], \end{align}

(B3)

and analogous expressions can be written for ${⟨ {\hat{σ}}_{i} ⟩}_{q c}$ and ${⟨ {\hat{σ}}_{i}, {\hat{σ}}_{j} ⟩}_{q c}$ ⁠. The reduced partition function

{\tilde{Z}}_{q c} = \int d y e^{- β Ω^{2} \frac{{(y - ε)}^{2}}{2 κ^{2}}} {T r}_{q} [e^{- β (Δ {\hat{σ}}_{1} + y {\hat{σ}}_{3})}]

(B4)

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

\hat{ρ} (y) = e^{- β Ω^{2} \frac{{(y - ε)}^{2}}{2 κ^{2}}} e^{- β (Δ {\hat{σ}}_{1} + y {\hat{σ}}_{3})},

(B5)

we can rewrite exact quantum–classical averages as

{\tilde{Z}}_{q c} = \int d y {T r}_{q} [\hat{ρ} (y)],

(B6a)

{〈 {\hat{σ}}_{i} 〉}_{q c} = \frac{1}{{\tilde{Z}}_{q c}} \int d y {T r}_{q} [\hat{ρ} (y) {\hat{σ}}_{i}],

(B6b)

〈 {\hat{σ}}_{i}, {{\hat{σ}}_{j} 〉}_{q c} = \frac{1}{{\tilde{Z}}_{q c}} \int d y \int_{0}^{1} d λ {T r}_{q} [\hat{ρ} {(y)}^{1 - λ} {\hat{σ}}_{i} \hat{ρ} {(y)}^{λ} {\hat{σ}}_{j}] .

(B6c)

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 $\hat{A}$ and $\hat{B}$ can be represented as an integral over a spherical phase space

{T r}_{q} [\hat{A} \hat{B}] = \int d u A_{W} (u) B_{W} (u),

(C1a)

where

A_{W} (u) = {T r}_{q} [\hat{A} {\hat{w}}_{W} (u)],

(C1b)

B_{W} (u) = {T r}_{q} [\hat{B} {\hat{w}}_{W} (u)],

(C1c)

{\hat{w}}_{W} (u) = \frac{1}{2} (\hat{I} + g_{W} {\hat{σ}}_{i} u_{i}),

(C1d)

and $g_{W} = \sqrt{3}$ in standard spin-LSC.³¹ Based on the mapping prescription in Eq. (C1), we define thermal expectation values and correlation functions within spin-LSC by

{⟨ σ_{j} (t) ⟩}^{W} = \frac{1}{Z_{q c}} \int d x d p e^{- β H_{b}} \int d u T_{0}^{W} (x, u) σ_{j}^{W} (u_{t}),

(C2a)

K_{i j}^{W} (t) = \frac{1}{Z_{q c}} \int d x d p e^{- β H_{b}} \int d u T_{i}^{W} (x, u) σ_{j}^{W} (u_{t}),

(C2b)

where

T_{μ}^{W} (x, u) = {T r}_{q} [{\hat{T}}_{μ} (x) {\hat{w}}_{W} (u)],

(C3a)

{\hat{T}}_{μ} (x) = \frac{1}{β} \int_{0}^{β} d λ e^{- (β - λ) ({\hat{H}}_{s} + {\hat{H}}_{s b} (x))} {\hat{σ}}_{μ} e^{- λ ({\hat{H}}_{s} + {\hat{H}}_{s b} (x))},

(C3b)

and $σ_{i}^{W} (u)$ has been defined in Eq. (12). At t = 0, ${⟨ σ_{j} (0) ⟩}^{W}$ and $K_{i j}^{W} (0)$ are equal by construction to their quantum–classical correspondents ${⟨ {\hat{σ}}_{i} ⟩}_{q c}$ and ${⟨ {\hat{σ}}_{i}, {\hat{σ}}_{j} ⟩}_{q c}$ ⁠, respectively. When the variables (x, p, u) are evolved under the Hamiltonian $H_{W} (x, p, u) = {T r}_{q} [\hat{H} {\hat{w}}_{W} (u)]$ ⁠, the time evolution of $K_{i j}^{W} (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

ρ_{b} = \frac{e^{- β H_{b}}}{Z_{b}}, Z_{b} = \int d x d p e^{- β H_{b}}

(C4)

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

\begin{align} \lim_{t \to + \infty} {⟨σ_{j} (t)⟩}^{W} & = \lim_{t \to + \infty} {⟨T_{0}^{W} σ_{j}^{W} (t)⟩}_{0} \\ = {⟨{T_{0}}^{W}⟩}_{0} {⟨σ_{j}⟩}_{eq}^{W} = {⟨σ_{j}⟩}_{eq}^{W}, \end{align}

(C5a)

\begin{align} \lim_{t \to + \infty} K_{i j}^{W} (t) & = \lim_{t \to + \infty} {⟨{T_{i}}^{W} {σ_{j}}^{W} (t)⟩}_{0} \\ = {⟨{T_{i}}^{W}⟩}_{0} {⟨σ_{j}⟩}_{eq}^{W} = {⟨{\hat{σ}}_{i}⟩}_{q c} {⟨σ_{j}⟩}_{eq}^{W}, \end{align}

(C5b)

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

{⟨ A ⟩}_{0} = \frac{1}{Z_{q c}} \int d x d p e^{- β H_{b}} \int d u A (x, u),

(C6)

with Z_qc as defined in Eq. (6), and with respect to the equilibrium distribution

{⟨ A ⟩}_{eq}^{W} = \frac{1}{Z_{eq}^{W}} \int d x d p \int d u e^{- β H_{W} (x, p, u)} A_{W} (x, u),

(C7a)

Z_{eq}^{W} = \int d x d p \int d u e^{- β H_{W} (x, p, u)} .

(C7b)

Also, we used that

\begin{align} {⟨{T_{μ}}^{W}⟩}_{0} & = \frac{1}{Z_{q c}} \int d x d p e^{- β H_{b}} \int d u {T_{μ}}^{W} (u, x) \\ = \frac{1}{Z_{q c}} \int d x d p e^{- β H_{b}} {T r}_{q} [{\hat{T}}_{μ} (x)] = {⟨{\hat{σ}}_{μ}⟩}_{q c} . \end{align}

(C8)

In general, ${⟨ σ_{j} ⟩}_{eq}^{W} \neq {⟨ {\hat{σ}}_{j} ⟩}_{q c}$ ⁠, 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

{T r}_{q} [\hat{A} \hat{B}] = \int d u A_{\bar{s}} (u) B_{s} (u),

(C9a)

where

A_{\bar{s}} (u) = {T r}_{q} [\hat{A} {\hat{w}}_{\bar{s}} (u)],

(C9b)

B_{s} (u) = {T r}_{q} [\hat{B} {\hat{w}}_{s} (u)],

(C9c)

{\hat{w}}_{s} (u) = \frac{1}{2} (\hat{I} + g_{s} {\hat{σ}}_{i} u_{i}),

(C9d)

and where ${\hat{w}}_{\bar{s}}$ is defined in the same way but with $g_{\bar{s}} = 3 / g_{s}$ ⁠. 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 ${⟨ {\hat{σ}}_{j} (t) ⟩}_{q c}$ 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

{⟨σ_{j}⟩}_{eq}^{(s)} \overset{!}{=} {⟨{\hat{σ}}_{j}⟩}_{q c},

(C10)

where ${⟨ σ_{j} ⟩}_{eq}^{(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 $σ_{3} (z) = g_{s} z$ ⁠. In the specific case of j = 3 [with the aim of studying K₃₃(t)], we find

{⟨ σ_{3} ⟩}_{eq}^{(s)} = \frac{\int_{- 1}^{1} d z g_{s} z f (z)}{\int_{- 1}^{1} d z f (z)},

(C11)

f (z) = I_{0} (- β g_{s} Δ \sqrt{1 - z^{2}}) e^{β g_{s} z (Λ g_{s} z / 4 - ε)},

(C12)

where we introduced the modified Bessel function of the first kind,⁹⁷ I₀(x), and the reorganization energy $Λ = \sum_{α = 1}^{F} 2 \frac{c_{α}^{2}}{m_{α} ω_{α}^{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⁹⁸ calculated with the open-source pyrho package.⁹⁹ 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.¹⁰⁰ 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

{\hat{ρ}}_{0} = \frac{1}{2} {\hat{σ}}_{0} \otimes {\hat{ρ}}_{b},

(D1)

where ${\hat{ρ}}_{b}$ 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,

〈 {\hat{σ}}_{j} {(t) 〉}_{q}^{b} = {〈 {\hat{σ}}_{j} 〉}_{q c}, \forall t \geq t_{eq},

(D2)

where ${⟨ \cdot ⟩}_{q}^{b}$ 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^{- i \hat{H} t} {\hat{ρ}}_{0} e^{i \hat{H} t} = \frac{1}{Z} e^{- β \hat{H}}$ ⁠.

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⁴

{\tilde{K}}_{i j} (ω) = {\tilde{C}}_{i j} (ω) \frac{1 - e^{- β ω}}{β ω},

(D3)

where $\tilde{f} (ω) = \int_{- \infty}^{+ \infty} d t e^{- i ω t} f (t)$ ⁠. Finally, Eq. (D3) was transformed back to time domain to get K_ij(t).

REFERENCES

1.

M. E.

Tuckerman

,

Statistical Mechanics: Theory and Molecular Simulation

(

Oxford University Press

,

2010

).

Google Scholar

2.

J. S.

Bader

and

B. J.

Berne

, “

Quantum and classical relaxation rates from classical simulations

,”

J. Chem. Phys.

100

,

8359

–

8366

(

1994

).

https://doi.org/10.1063/1.466780

Google Scholar

Crossref

3.

J.

Cao

and

G. A.

Voth

, “

A new perspective on quantum time correlation functions

,”

J. Chem. Phys.

99

,

10070

–

10073

(

1993

).

https://doi.org/10.1063/1.465512

Google Scholar

Crossref

4.

I. R.

Craig

and

D. E.

Manolopoulos

, “

Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics

.”

J. Chem. Phys.

121

,

3368

–

3373

(

2004

).

https://doi.org/10.1063/1.1777575

Google Scholar

Crossref

PubMed

5.

J.

Liu

, “

Path integral Liouville dynamics for thermal equilibrium systems

,”

J. Chem. Phys.

140

,

224107

(

2014

).

https://doi.org/10.1063/1.4881518

Google Scholar

Crossref

PubMed

6.

W. H.

Miller

, “

The semiclassical initial value representation: A potentially practical way for adding quantum effects to classical molecular dynamics simulations

,”

J. Phys. Chem. A

105

,

2942

–

2955

(

2001

).

https://doi.org/10.1021/jp003712k

Google Scholar

Crossref

7.

P. V.

Parandekar

and

J. C.

Tully

, “

Detailed balance in Ehrenfest mixed quantum-classical dynamics

,”

J. Chem. Theory Comput.

2

,

229

–

235

(

2006

).

https://doi.org/10.1021/ct050213k

Google Scholar

Crossref

PubMed

8.

J. E.

Runeson

,

J. R.

Mannouch

,

G.

Amati

,

M. R.

Fiechter

, and

J. O.

Richardson

, “

Spin-mapping methods for simulating ultrafast nonadiabatic dynamics

,”

Chimia

76

,

582

(

2022

).

https://doi.org/10.2533/chimia.2022.582

Google Scholar

Crossref

9.

J. C.

Tully

, “

Molecular dynamics with electronic transitions

,”

J. Chem. Phys.

93

,

1061

–

1071

(

1990

).

https://doi.org/10.1063/1.459170

Google Scholar

Crossref

10.

P. V.

Parandekar

and

J. C.

Tully

, “

Mixed quantum-classical equilibrium

,”

J. Chem. Phys.

122

,

094102

(

2005

).

https://doi.org/10.1063/1.1856460

Google Scholar

Crossref

PubMed

11.

J. R.

Schmidt

,

P. V.

Parandekar

, and

J. C.

Tully

, “

Mixed quantum-classical equilibrium: Surface hopping

,”

J. Chem. Phys.

129

,

044104

(

2008

).

https://doi.org/10.1063/1.2955564

Google Scholar

Crossref

PubMed

12.

W. H.

Miller

and

S. J.

Cotton

, “

Classical molecular dynamics simulation of electronically non-adiabatic processes

,”

Faraday Discuss.

195

,

9

–

30

(

2016

).

https://doi.org/10.1039/c6fd00181e

Google Scholar

Crossref

PubMed

13.

W. H.

Miller

and

S. J.

Cotton

, “

Communication: Note on detailed balance in symmetrical quasi-classical models for electronically non-adiabatic dynamics

,”

J. Chem. Phys.

142

,

131103

(

2015

).

https://doi.org/10.1063/1.4916945

Google Scholar

Crossref

PubMed

14.

V.

May

and

O.

Kühn

,

Charge and Energy Transfer Dynamics in Molecular Systems

, 3rd ed. (

Wiley VCH

,

Weinheim

,

2011

).

Google Scholar

Crossref

15.

J. O.

Richardson

and

M.

Thoss

, “

Communication: Nonadiabatic ring-polymer molecular dynamics

,”

J. Chem. Phys.

139

,

031102

(

2013

).

https://doi.org/10.1063/1.4816124

Google Scholar

Crossref

PubMed

16.

N.

Ananth

, “

Mapping variable ring polymer molecular dynamics: A path-integral based method for nonadiabatic processes

,”

J. Chem. Phys.

139

,

124102

(

2013

).

https://doi.org/10.1063/1.4821590

Google Scholar

Crossref

PubMed

17.

J. O.

Richardson

,

P.

Meyer

,

M.-O.

Pleinert

, and

M.

Thoss

, “

An analysis of nonadiabatic ring-polymer molecular dynamics and its application to vibronic spectra

,”

Chem. Phys.

482

,

124

–

134

(

2017

).

https://doi.org/10.1016/j.chemphys.2016.09.036

Google Scholar

Crossref

18.

S. N.

Chowdhury

and

P.

Huo

, “

Coherent state mapping ring polymer molecular dynamics for non-adiabatic quantum propagations

,”

J. Chem. Phys.

147

,

214109

(

2017

).

https://doi.org/10.1063/1.4995616

Google Scholar

Crossref

PubMed

19.

X.

Tao

,

P.

Shushkov

, and

T. F.

Miller

III, “

Path-integral isomorphic Hamiltonian for including nuclear quantum effects in non-adiabatic dynamics

,”

J. Chem. Phys.

148

,

102327

(

2018

).

https://doi.org/10.1063/1.5005544

Google Scholar

Crossref

PubMed

20.

H.-D.

Meyer

and

W. H.

Miller

, “

A classical analog for electronic degrees of freedom in nonadiabatic collision processes

,”

J. Chem. Phys.

70

,

3214

–

3223

(

1979

).

https://doi.org/10.1063/1.437910

Google Scholar

Crossref

21.

G.

Stock

and

M.

Thoss

, “

Semiclassical description of nonadiabatic quantum dynamics

,”

Phys. Rev. Lett.

78

,

578

–

581

(

1997

).

https://doi.org/10.1103/physrevlett.78.578

Google Scholar

Crossref

22.

W. H.

Miller

, “

Electronically nonadiabatic dynamics via semiclassical initial value methods

,”

J. Phys. Chem. A

113

,

1405

–

1415

(

2009

).

https://doi.org/10.1021/jp809907p

Google Scholar

Crossref

PubMed

23.

A.

Kelly

,

R.

van Zon

,

J.

Schofield

, and

R.

Kapral

, “

Mapping quantum-classical Liouville equation: Projectors and trajectories

,”

J. Chem. Phys.

136

,

084101

(

2012

).

https://doi.org/10.1063/1.3685420

Google Scholar

Crossref

PubMed

24.

C.-Y.

Hsieh

and

R.

Kapral

, “

Nonadiabatic dynamics in open quantum-classical systems: Forward-backward trajectory solution

,”

J. Chem. Phys.

137

,

22A507

(

2012

).

https://doi.org/10.1063/1.4736841

Google Scholar

Crossref

PubMed

25.

P.

Huo

and

D. F.

Coker

, “

Consistent schemes for non-adiabatic dynamics derived from partial linearized density matrix propagation

,”

J. Chem. Phys.

137

,

22A535

(

2012

).

https://doi.org/10.1063/1.4748316

Google Scholar

Crossref

PubMed

26.

S. J.

Cotton

and

W. H.

Miller

, “

Symmetrical windowing for quantum states in quasi-classical trajectory simulations: Application to electronically non-adiabatic processes

,”

J. Chem. Phys.

139

,

234112

(

2013

).

https://doi.org/10.1063/1.4845235

Google Scholar

Crossref

PubMed

27.

J.

Liu

, “

A unified theoretical framework for mapping models for the multi-state Hamiltonian

,”

J. Chem. Phys.

145

,

204105

(

2016

).

https://doi.org/10.1063/1.4967815

Google Scholar

Crossref

PubMed

28.

M. A. C.

Saller

,

A.

Kelly

, and

J. O.

Richardson

, “

On the identity of the identity operator in nonadiabatic linearized semiclassical dynamics

,”

J. Chem. Phys.

150

,

071101

(

2019

).

https://doi.org/10.1063/1.5082596

Google Scholar

Crossref

PubMed

29.

G.

Stock

and

U.

Müller

, “

Flow of zero-point energy and exploration of phase space in classical simulations of quantum relaxation dynamics

,”

J. Chem. Phys.

111

,

65

–

76

(

1999

).

https://doi.org/10.1063/1.479254

Google Scholar

Crossref

30.

X.

He

,

B.

Wu

,

Z.

Gong

, and

J.

Liu

, “

Commutator matrix in phase space mapping models for nonadiabatic quantum dynamics

,”

J. Phys. Chem. A

125

,

6845

–

6863

(

2021

).

https://doi.org/10.1021/acs.jpca.1c04429

Google Scholar

Crossref

PubMed

31.

J. E.

Runeson

and

J. O.

Richardson

, “

Spin-mapping approach for nonadiabatic molecular dynamics

,”

J. Chem. Phys.

151

,

044119

(

2019

).

https://doi.org/10.1063/1.5100506

Google Scholar

Crossref

PubMed

32.

J. E.

Runeson

and

J. O.

Richardson

, “

Generalized spin mapping for quantum-classical dynamics

,”

J. Chem. Phys.

152

,

084110

(

2020

).

https://doi.org/10.1063/1.5143412

Google Scholar

Crossref

PubMed

33.

J. R.

Mannouch

and

J. O.

Richardson

, “

A partially linearized spin-mapping approach for nonadiabatic dynamics. I. Derivation of the theory

,”

J. Chem. Phys.

153

,

194109

(

2020

).

https://doi.org/10.1063/5.0031168

Google Scholar

Crossref

PubMed

34.

J. R.

Mannouch

and

J. O.

Richardson

, “

A partially linearized spin-mapping approach for nonadiabatic dynamics. II. Analysis and comparison with related approaches

,”

J. Chem. Phys.

153

,

194110

(

2020

).

https://doi.org/10.1063/5.0031173

Google Scholar

Crossref

PubMed

35.

X.

He

and

J.

Liu

, “

A new perspective for nonadiabatic dynamics with phase space mapping models

,”

J. Chem. Phys.

151

,

024105

(

2019

).

https://doi.org/10.1063/1.5108736

Google Scholar

Crossref

PubMed

36.

O.

Bramley

,

C.

Symonds

, and

D. V.

Shalashilin

, “

Quantum system-bath dynamics with quantum superposition sampling and coupled generalized coherent states

,”

J. Chem. Phys.

151

,

064103

(

2019

).

https://doi.org/10.1063/1.5100145

Google Scholar

Crossref

37.

R.

Kubo

,

M.

Toda

, and

N.

Hashitsume

,

Statistical Physics II

(

Springer

,

New York

,

1978

).

Google Scholar

38.

A.

Nitzan

,

Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems

(

Oxford University Press

,

Oxford

,

2006

).

Google Scholar

Crossref

39.

S.

Mukamel

,

Principles of Nonlinear Optical Spectroscopy

(

Oxford University Press

,

1995

).

Google Scholar

40.

D.

Chandler

,

Introduction to Modern Statistical Mechanics

(

Oxford University Press

,

New York

,

1987

).

Google Scholar

41.

G.

Amati

,

H.

Meyer

, and

T.

Schilling

, “

Memory effects in the Fermi–Pasta–Ulam model

,”

J. Stat. Phys.

174

,

219

–

257

(

2019

).

https://doi.org/10.1007/s10955-018-2207-6

Google Scholar

Crossref

42.

R.

Zwanzig

,

Nonequilibrium Statistical Mechanics

(

Oxford University Press

,

New York

,

2001

).

Google Scholar

43.

S. C.

Althorpe

, “

Path-integral approximations to quantum dynamics

,”

Eur. Phys. J. B

94

,

1

–

17

(

2021

).

https://doi.org/10.1140/epjb/s10051-021-00155-2

Google Scholar

Crossref

44.

F.

Mauri

,

R.

Car

, and

E.

Tosatti

, “

Canonical statistical averages of coupled quantum-classical systems

,”

Europhys. Lett.

24

,

431

–

436

(

1993

).

https://doi.org/10.1209/0295-5075/24/6/001

Google Scholar

Crossref

45.

S.

Nielsen

,

R.

Kapral

, and

G.

Ciccotti

, “

Statistical mechanics of quantum-classical systems

,”

J. Chem. Phys.

115

,

5805

–

5815

(

2001

).

https://doi.org/10.1063/1.1400129

Google Scholar

Crossref

46.

Note that some of the present expressions differ by a factor 2 from previous work which used a basis of ${\hat{S}}_{i} = \frac{1}{2} {\hat{σ}}_{i}$ ⁠; the form presented here is equivalent but leads to slightly neater equations.

47.

R. L.

Stratonovich

, “

On distributions in representation space

,”

Sov. Phys. JETP

4

,

891

–

898

(

1957

); available at http://jetp.ras.ru/cgi-bin/e/index/e/4/6/p891?a=list.

Google Scholar

48.

T.

Tilma

,

M. J.

Everitt

,

J. H.

Samson

,

W. J.

Munro

, and

K.

Nemoto

, “

Wigner functions for arbitrary quantum systems

,”

Phys. Rev. Lett.

117

,

180401

(

2016

).

https://doi.org/10.1103/physrevlett.117.180401

Google Scholar

Crossref

PubMed

49.

R. P.

Rundle

,

T.

Tilma

,

J. H.

Samson

,

V. M.

Dwyer

,

R. F.

Bishop

, and

M. J.

Everitt

, “

General approach to quantum mechanics as a statistical theory

,”

Phys. Rev. A

99

,

012115

(

2019

).

https://doi.org/10.1103/physreva.99.012115

Google Scholar

Crossref

50.

C.

Brif

and

A.

Mann

, “

Phase-space formulation of quantum mechanics and quantum-state reconstruction for physical systems with Lie-group symmetries

,”

Phys. Rev. A

59

,

971

(

1999

).

https://doi.org/10.1103/physreva.59.971

Google Scholar

Crossref

51.

J. E.

Runeson

and

J. O.

Richardson

, “

Quantum entanglement from classical trajectories

,”

Phys. Rev. Lett.

127

,

250403

(

2021

).

https://doi.org/10.1103/physrevlett.127.250403

Google Scholar

Crossref

PubMed

52.

J. R.

Mannouch

and

J. O.

Richardson

, “

A partially linearized spin-mapping approach for simulating nonlinear optical spectra

,”

J. Chem. Phys.

156

,

024108

(

2022

).

https://doi.org/10.1063/5.0077744

Google Scholar

Crossref

PubMed

53.

J. E.

Runeson

,

J. E.

Lawrence

,

J. R.

Mannouch

, and

J. O.

Richardson

, “

Explaining the efficiency of photosynthesis: Quantum uncertainty or classical vibrations?

,”

J. Phys. Chem. Lett.

13

,

3392

–

3399

(

2022

).

https://doi.org/10.1021/acs.jpclett.2c00538

Google Scholar

Crossref

PubMed

54.

G.

Amati

,

M. A.

Saller

,

A.

Kelly

, and

J. O.

Richardson

, “

Quasiclassical approaches to the generalized quantum master equation

,”

J. Chem. Phys.

157

,

234103

(

2022

).

https://doi.org/10.1063/5.0124028

Google Scholar

Crossref

PubMed

55.

C.

Kittel

and

P.

McEuen

,

Introduction to Solid State Physics

(

Wiley

,

New York

,

1996

), 8th edition.

Google Scholar

56.

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

57.

S. J.

Cotton

and

W. H.

Miller

, “

Trajectory-adjusted electronic zero point energy in classical Meyer-Miller vibronic dynamics: Symmetrical quasiclassical application to photodissociation

,”

J. Chem. Phys.

150

,

194110

(

2019

).

https://doi.org/10.1063/1.5094458

Google Scholar

Crossref

PubMed

58.

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.

59.

M. A.

Nielsen

and

I.

Chuang

,

Quantum Computation and Quantum Information

(

Cambridge University Press

,

2002

).

Google Scholar

60.

A.

Abragam

and

B.

Bleaney

,

Electron Paramagnetic Resonance of Transition Ions

, International Series of Monographs on Physics (

Clarendon Press

,

Oxford

,

1970

).

Google Scholar

61.

A. G.

Maryasov

and

M. K.

Bowman

, “

Spin dynamics of paramagnetic centers with anisotropic g tensor and spin of 1/2

,”

J. Magn. Reson.

221

,

69

–

75

(

2012

).

https://doi.org/10.1016/j.jmr.2012.05.011

Google Scholar

Crossref

PubMed

62.

A.

Garg

,

J. N.

Onuchic

, and

V.

Ambegaokar

, “

Effect of friction on electron transfer in biomolecules

,”

J. Chem. Phys.

83

,

4491

–

4503

(

1985

).

https://doi.org/10.1063/1.449017

Google Scholar

Crossref

63.

C. H.

Mak

and

D.

Chandler

, “

Coherent-incoherent transition and relaxation in condensed-phase tunneling systems

,”

Phys. Rev. A

44

,

2352

–

2369

(

1991

).

https://doi.org/10.1103/physreva.44.2352

Google Scholar

Crossref

PubMed

64.

N.

Makri

, “

Numerical path integral techniques for long time dynamics of quantum dissipative systems

,”

J. Math. Phys.

36

,

2430

–

2457

(

1995

).

https://doi.org/10.1063/1.531046

Google Scholar

Crossref

65.

M.

Thoss

,

H.

Wang

, and

W. H.

Miller

, “

Self-consistent hybrid approach for complex systems: Application to the spin-boson model with Debye spectral density

,”

J. Chem. Phys.

115

,

2991

(

2001

).

https://doi.org/10.1063/1.1385562

Google Scholar

Crossref

66.

A. O.

Caldeira

and

A. J.

Leggett

, “

Path integral approach to quantum Brownian motion

,”

Physica A

121

,

587

–

616

(

1983

).

https://doi.org/10.1016/0378-4371(83)90013-4

Google Scholar

Crossref

67.

D.

Xu

and

K.

Schulten

, “

Coupling of protein motion to electron transfer in a photosynthetic reaction center: Investigating the low temperature behavior in the framework of the spin—boson model

,”

Chem. Phys.

182

,

91

–

117

(

1994

).

https://doi.org/10.1016/0301-0104(94)00016-6

Google Scholar

Crossref

68.

M.

Merkli

,

G. P.

Berman

, and

R.

Sayre

, “

Electron transfer reactions: Generalized spin-boson approach

,”

J. Math. Chem.

51

,

890

–

913

(

2013

).

https://doi.org/10.1007/s10910-012-0124-5

Google Scholar

Crossref

69.

L.

Magazzù

,

P.

Forn-Díaz

,

R.

Belyansky

,

J.-L.

Orgiazzi

,

M. A.

Yurtalan

,

M. R.

Otto

,

A.

Lupascu

,

C. M.

Wilson

, and

M.

Grifoni

, “

Probing the strongly driven spin-boson model in a superconducting quantum circuit

,”

Nat. Commun.

9

,

1403

(

2018

).

https://doi.org/10.1038/s41467-018-03626-w

Google Scholar

Crossref

PubMed

70.

J. T.

Stockburger

, “

Simulating spin-boson dynamics with stochastic Liouville–von Neumann equations

,”

Chem. Phys.

296

,

159

–

169

(

2004

).

https://doi.org/10.1016/j.chemphys.2003.09.014

Google Scholar

Crossref

71.

A.

Mitra

and

A. J.

Millis

, “

Spin dynamics and violation of the fluctuation dissipation theorem in a nonequilibrium ohmic spin-boson model

,”

Phys. Rev. B

72

,

121102

(

2005

).

https://doi.org/10.1103/physrevb.72.121102

Google Scholar

Crossref

72.

N.

Boudjada

and

D.

Segal

, “

From dissipative dynamics to studies of heat transfer at the nanoscale: Analysis of the spin-boson model

,”

J. Phys. Chem. A

118

,

11323

–

11336

(

2014

).

https://doi.org/10.1021/jp5091685

Google Scholar

Crossref

PubMed

73.

S.

Kundu

and

N.

Makri

, “

Time evolution of bath properties in spin-boson dynamics

,”

J. Phys. Chem. B

125

,

8137

–

8151

(

2021

).

https://doi.org/10.1021/acs.jpcb.1c03861

Google Scholar

Crossref

PubMed

74.

H.

Wang

,

D. E.

Skinner

, and

M.

Thoss

, “

Calculation of reactive flux correlation functions for systems in a condensed phase environment: A multilayer multiconfiguration time-dependent Hartree approach

,”

J. Chem. Phys.

125

,

174502

(

2006

).

https://doi.org/10.1063/1.2363195

Google Scholar

Crossref

PubMed

75.

A.

Montoya-Castillo

and

D. R.

Reichman

, “

Approximate but accurate quantum dynamics from the Mori formalism. II. Equilibrium time correlation functions

,”

J. Chem. Phys.

146

,

084110

(

2017

).

https://doi.org/10.1063/1.4975388

Google Scholar

Crossref

PubMed

76.

A.

Montoya-Castillo

and

D. R.

Reichman

, “

Path integral approach to the Wigner representation of canonical density operators for discrete systems coupled to harmonic baths

,”

J. Chem. Phys.

146

,

024107

(

2017

).

https://doi.org/10.1063/1.4973646

Google Scholar

Crossref

PubMed

77.

I. R.

Craig

,

M.

Thoss

, and

H.

Wang

, “

Proton transfer reactions in model condensed-phase environments: Accurate quantum dynamics using the multilayer multiconfiguration time-dependent Hartree approach

,”

J. Chem. Phys.

127

,

144503

(

2007

).

https://doi.org/10.1063/1.2772265

Google Scholar

Crossref

PubMed

78.

Y.

Tanimura

and

R.

Kubo

, “

Time evolution of a quantum system in contact with a nearly Gaussian-Markoffian noise bath

,”

J. Phys. Soc. Jpn.

58

,

101

–

114

(

1989

).

https://doi.org/10.1143/jpsj.58.101

Google Scholar

Crossref

79.

A.

Bastida

,

C.

Cruz

,

J.

Zúñiga

,

A.

Requena

, and

B.

Miguel

, “

The Ehrenfest method with quantum corrections to simulate the relaxation of molecules in solution: Equilibrium and dynamics

,”

J. Chem. Phys.

126

,

014503

(

2007

).

https://doi.org/10.1063/1.2404676

Google Scholar

Crossref

PubMed

80.

N.

Bellonzi

,

A.

Jain

, and

J. E.

Subotnik

, “

An assessment of mean-field mixed semiclassical approaches: Equilibrium populations and algorithm stability

,”

J. Chem. Phys.

144

,

154110

(

2016

).

https://doi.org/10.1063/1.4946810

Google Scholar

Crossref

PubMed

81.

E. A.

Coronado

,

J.

Xing

, and

W. H.

Miller

, “

Ultrafast non-adiabatic dynamics of systems with multiple surface crossings: A test of the Meyer–Miller Hamiltonian with semiclassical initial value representation methods

,”

Chem. Phys. Lett.

349

,

521

–

529

(

2001

).

https://doi.org/10.1016/s0009-2614(01)01242-8

Google Scholar

Crossref

82.

J.

Zheng

,

Y.

Xie

,

S.-s.

Jiang

,

Y.-z.

Long

,

X.

Ning

, and

Z.-g.

Lan

, “

Ultrafast electron transfer with symmetrical quasi-classical dynamics based on mapping Hamiltonian and quantum dynamics based on ML-MCTDH

,”

Chin. J. Chem. Phys.

30

,

800

–

810

(

2017

).

https://doi.org/10.1063/1674-0068/30/cjcp1711210

Google Scholar

Crossref

83.

G.

Tao

and

N.

Shen

, “

Mapping state space to quasiclassical trajectory dynamics in coherence-controlled nonadiabatic simulations for condensed phase problems

,”

J. Phys. Chem. A

121

,

1734

–

1747

(

2017

).

https://doi.org/10.1021/acs.jpca.6b10936

Google Scholar

Crossref

PubMed

84.

W. H.

Miller

,

S. D.

Schwartz

, and

J. W.

Tromp

, “

Quantum mechanical rate constants for bimolecular reactions

,”

J. Chem. Phys.

79

,

4889

–

4898

(

1983

).

https://doi.org/10.1063/1.445581

Google Scholar

Crossref

85.

W. H.

Miller

, “

Beyond transition-state theory: A rigorous quantum theory of chemical reaction rates

,”

Acc. Chem. Res.

26

,

174

(

1993

).

https://doi.org/10.1021/ar00028a007

Google Scholar

Crossref

86.

J. E.

Lawrence

,

T.

Fletcher

,

L. P.

Lindoy

, and

D. E.

Manolopoulos

, “

On the calculation of quantum mechanical electron transfer rates

,”

J. Chem. Phys.

151

,

114119

(

2019

).

https://doi.org/10.1063/1.5116800

Google Scholar

Crossref

PubMed

87.

J. E.

Subotnik

,

A.

Jain

,

B.

Landry

,

A.

Petit

,

W.

Ouyang

, and

N.

Bellonzi

, “

Understanding the surface hopping view of electronic transitions and decoherence

,”

Annu. Rev. Phys. Chem.

67

,

387

–

417

(

2016

).

https://doi.org/10.1146/annurev-physchem-040215-112245

Google Scholar

Crossref

PubMed

88.

S. C.

Althorpe

,

N.

Ananth

,

G.

Angulo

,

R. D.

Astumian

,

V.

Beniwal

,

J.

Blumberger

,

P. G.

Bolhuis

,

B.

Ensing

,

D. R.

Glowacki

,

S.

Habershon

et al, “

Non-adiabatic reactions: General discussion

,”

Faraday Discuss.

195

,

311

–

344

(

2016

).

https://doi.org/10.1039/c6fd90078j

Google Scholar

Crossref

PubMed

89.

T. J. H.

Hele

, “

An electronically non-adiabatic generalization of ring polymer molecular dynamics

,”

Master’s thesis

,

Oxford University

,

2011

; arXiv:1308.3950.

Google Scholar

90.

M. A. C.

Saller

,

J. E.

Runeson

, and

J. O.

Richardson

, “

Path-integral approaches to non-adiabatic dynamics

,” in

Quantum Chemistry and Dynamics of Excited States

, edited by

L.

González

and

R.

Lindh

(

Wiley

,

2020

), Chap. 20, pp.

629

–

653

.

https://doi.org/10.1002/9781119417774.ch20

Google Scholar

Crossref

91.

D.

Bossion

,

S. N.

Chowdhury

, and

P.

Huo

, “

Non-adiabatic ring polymer molecular dynamics with spin mapping variables

,”

J. Chem. Phys.

154

,

184106

(

2021

).

https://doi.org/10.1063/5.0051456

Google Scholar

Crossref

PubMed

92.

X.

Tao

,

P.

Shushkov

, and

T. F.

Miller III

, “

Simple flux-side formulation of state-resolved thermal reaction rates for ring-polymer surface hopping

,”

J. Phys. Chem. A

123

,

3013

–

3020

(

2019

).

https://doi.org/10.1021/acs.jpca.9b00877

Google Scholar

Crossref

PubMed

93.

X.

Tao

,

P.

Shushkov

, and

T. F.

Miller III

, “

Microcanonical rates from ring-polymer molecular dynamics: Direct-shooting, stationary-phase, and maximum-entropy approaches

,”

J. Chem. Phys.

152

,

124117

(

2020

).

https://doi.org/10.1063/1.5144307

Google Scholar

Crossref

PubMed

94.

J. E.

Lawrence

and

D. E.

Manolopoulos

, “

An analysis of isomorphic RPMD in the golden rule limit

,”

J. Chem. Phys.

151

,

244109

(

2019

).

https://doi.org/10.1063/1.5138913

Google Scholar

Crossref

PubMed

95.

A.

Kelly

,

A.

Montoya-Castillo

,

L.

Wang

, and

T. E.

Markland

, “

Generalized quantum master equations in and out of equilibrium: When can one win?

,”

J. Chem. Phys.

144

,

184105

(

2016

).

https://doi.org/10.1063/1.4948612

Google Scholar

Crossref

PubMed

96.

H.

Wang

and

M.

Thoss

, “

A multilayer multiconfiguration time-dependent Hartree simulation of the reaction-coordinate spin-boson model employing an interaction picture

,”

J. Chem. Phys.

146

,

124112

(

2017

).

https://doi.org/10.1063/1.4978901

Google Scholar

Crossref

PubMed

97.

F.

Bowman

,

Introduction to Bessel Functions

(

Courier Corporation

,

New York

,

2012

).

Google Scholar

98.

Y.

Tanimura

, “

Numerically ‘exact’ approach to open quantum dynamics: The hierarchical equations of motion (HEOM)

,”

J. Chem. Phys.

153

,

020901

(

2020

).

https://doi.org/10.1063/5.0011599

Google Scholar

Crossref

PubMed

99.

T. C.

Berkelbach

, “

A python package for reduced density matrix techniques

,” https://github.com/berkelbach-group/pyrho, 2020.

100.

Y.

Tanimura

, “

Reduced hierarchical equations of motion in real and imaginary time: Correlated initial states and thermodynamic quantities

,”

J. Chem. Phys.

141

,

044114

(

2014

).

https://doi.org/10.1063/1.4890441

Google Scholar

Crossref

PubMed

2023

Author(s)

On detailed balance in nonadiabatic dynamics: From spin spheres to equilibrium ellipsoids

I. INTRODUCTION

II. THEORY

A. Summary of (spherical) spin mapping

1. Problems of spin mapping in equilibrium

B. Spherical mapping with an optimized radius

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

D. Ellipsoid mapping for isolated subsystems

E. Ellipsoid mapping for mixed quantum–classical systems

1. Global construction

2. Local construction

III. MODEL

IV. RESULTS

A. Detailed balance in rate theory

V. OUTLOOK TO NONADIABATIC PATH-INTEGRAL METHODS

VI. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGEMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: PRINCIPAL-AXIS BASIS

APPENDIX B: CALCULATION OF QUANTUM–CLASSICAL THERMAL AVERAGES

APPENDIX C: SPHERICAL SPIN MAPPING FOR THERMAL CORRELATION FUNCTION

1. Spin-LSC

2. Optimized sphere

APPENDIX D: BENCHMARK KUBO-TRANSFORMED THERMAL CORRELATION FUNCTIONS

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

On detailed balance in nonadiabatic dynamics: From spin spheres to equilibrium ellipsoids

I. INTRODUCTION

II. THEORY

A. Summary of (spherical) spin mapping

1. Problems of spin mapping in equilibrium

B. Spherical mapping with an optimized radius

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

D. Ellipsoid mapping for isolated subsystems

E. Ellipsoid mapping for mixed quantum–classical systems

1. Global construction

2. Local construction

III. MODEL

IV. RESULTS

A. Detailed balance in rate theory

V. OUTLOOK TO NONADIABATIC PATH-INTEGRAL METHODS

VI. CONCLUSIONS

SUPPLEMENTARY MATERIAL

ACKNOWLEDGEMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: PRINCIPAL-AXIS BASIS

APPENDIX B: CALCULATION OF QUANTUM–CLASSICAL THERMAL AVERAGES

APPENDIX C: SPHERICAL SPIN MAPPING FOR THERMAL CORRELATION FUNCTION

1. Spin-LSC

2. Optimized sphere

APPENDIX D: BENCHMARK KUBO-TRANSFORMED THERMAL CORRELATION FUNCTIONS

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only