Spin-mapping approach for nonadiabatic molecular dynamics

TL;DR

This paper introduces a spin-mapping approach for nonadiabatic molecular dynamics that improves accuracy over traditional methods by using a quantum-mechanically exact spin representation, avoiding extra projection operators, and showing significant performance gains.

Contribution

The authors develop a spin-based mapping method for nonadiabatic dynamics that enhances accuracy and simplifies implementation compared to existing MMST-based approaches.

Findings

Significant improvement over Ehrenfest dynamics.

Better results than linearized semiclassical MMST mapping.

No additional computational complexity required.

Abstract

We propose a trajectory-based method for simulating nonadiabatic dynamics in molecular systems with two coupled electronic states. Employing a quantum-mechanically exact mapping of the two-level problem to a spin-1/2 coherent state, we construct a classical phase space of a spin vector constrained to a spherical surface with a radius consistent with the quantum magnitude of the spin. In contrast with the singly-excited harmonic oscillator basis used in Meyer-Miller-Stock-Thoss (MMST) mapping, the theory requires no additional projection operators onto the space of physical states. When treated under a quasiclassical approximation, we show that the resulting dynamics is equivalent to that generated by the MMST Hamiltonian. What differs is the value of the zero-point energy parameter as well as the initial distribution and the measurement operators. For various spin-boson models the results of our method are seen to be a significant improvement compared to both standard Ehrenfest dynamics and linearized semiclassical MMST mapping, without adding any computational complexity.

Tables

Table 1. Table 1: Spin radii, r s subscript 𝑟 𝑠 r_{s} , and dual symbols, s ¯ ¯ 𝑠 {\bar{s}} , of the three phase-space functions.

$s$	Q	W	P
$\overline{s}$	P	W	Q
$r_s$	1/2	$\sqrt{3}/2$	3/2

Table 2. Table 2: Definition of the terms used in Eq. ( 42 ) to compute general correlation functions C A ​ B ​ ( t ) subscript 𝐶 𝐴 𝐵 𝑡 C_{AB}(t) for various methods. We use the short-hand notation R 2 = X 1 2 + P 1 2 + X 2 2 + P 2 2 superscript 𝑅 2 superscript subscript 𝑋 1 2 superscript subscript 𝑃 1 2 superscript subscript 𝑋 2 2 superscript subscript 𝑃 2 2 R^{2}=X_{1}^{2}+P_{1}^{2}+X_{2}^{2}+P_{2}^{2} , ϕ = 16 ​ exp ⁡ ( − R 2 ) italic-ϕ 16 superscript 𝑅 2 \phi=16\exp(-R^{2}) and 𝒩 = [ 2 ​ π 2 ​ ( γ + 1 ) ] − 1 𝒩 superscript delimited-[] 2 superscript 𝜋 2 𝛾 1 1 \mathcal{N}=[2\pi^{2}(\gamma+1)]^{-1} .

Method	${\rho }_{\text{el}}(𝖷,𝖯)$	$\gamma$	$A$	$B$
Q	$𝒩\delta \left(R^2-2\left(\gamma +1\right)\right)$	0	$A_{\text{SW}}$	$B_{\overline{\text{SW}}}$
W	$𝒩\delta \left(R^2-2\left(\gamma +1\right)\right)$	$\sqrt{3}-1$	$A_{\text{SW}}$	$B_{\overline{\text{SW}}}$
P	$𝒩\delta \left(R^2-2\left(\gamma +1\right)\right)$	2	$A_{\text{SW}}$	$B_{\overline{\text{SW}}}$
PBME	$\varphi /\left(8{\pi }^2\right)$	1	${\tilde{A}}_{\text{SEO}}$	$B_{\text{wig}}$
LSC-IVR	${\varphi }^2/\left(8{\pi }^2\right)$	1	${\tilde{A}}_{\text{SEO}}$	${\tilde{B}}_{\text{SEO}}$

Table 3. Table 3: Definition of the terms used in Eq. ( 42 ) to compute correlation functions C A ​ B ​ ( t ) subscript 𝐶 𝐴 𝐵 𝑡 C_{AB}(t) , where the initial state is A ^ = | n ⟩ ​ ⟨ n | ^ 𝐴 ket 𝑛 bra 𝑛 \hat{A}=|n\rangle\langle n| , for various methods employing focused initial conditions. The electronic distribution in each case is given by ρ el ​ ( 𝖷 , 𝖯 ) = 𝒩 ​ δ ​ ( R 2 − 2 ​ ( γ + 1 ) ) subscript 𝜌 el 𝖷 𝖯 𝒩 𝛿 superscript 𝑅 2 2 𝛾 1 \rho_{\text{el}}(\mathsf{X},\mathsf{P})=\mathcal{N}\delta(R^{2}-2(\gamma+1)) and A = δ ​ ( X n 2 + P n 2 − ( γ + 2 ) ) 𝐴 𝛿 superscript subscript 𝑋 𝑛 2 superscript subscript 𝑃 𝑛 2 𝛾 2 A=\delta(X_{n}^{2}+P_{n}^{2}-(\gamma+2)) , with 𝒩 = π − 2 𝒩 superscript 𝜋 2 \mathcal{N}=\pi^{-2} .

Method	$\gamma$	$B$
Ehrenfest	0	$B_{\text{SW}}=B_{\text{MFT}}$
MMST-focused	$1$	$B_{\text{SW}}=B_{\text{wig}}$
W-focused	$\sqrt{3}-1$	$B_{\text{SW}}$

Table 4. Table 4: Parameter settings of the spin-boson models, referring to the panels in Figures  3 – 5 .

Panel	$\epsilon /\Delta$	$\xi$	$\beta \Delta$	${\omega }_{\text{c}}/\Delta$
(a)	0	$0.09$	$0.1$	2.5
(b)	0	$0.09$	$5$	2.5
(c)	1	$0.1$	$5$	2.5
(d)	1	$0.1$	$0.25$	1.0
(e)	0	$2$	$1$	1.0
(f)	5	$4$	$0.1$	2.0

Equations

$$\hat{\mathcal{H}}$$

$$\hat{V}(x)$$

$$|n⟩⟨m|\mapsto {\hat{a}}_n^{†}{\hat{a}}_m.$$

$$|n⟩⟨n|$$

$$|n⟩⟨m|+|m⟩⟨n|$$

$$\mathrm{i}(|n⟩⟨m|-|m⟩⟨n|)$$

$$H=\frac{p^2}{2m}+U+\sum _{nm}{V}_{nm}\frac{1}{2}(X_n{X}_m+P_n{P}_m-{\delta }_{nm}),$$

$$H=\frac{p^2}{2m}+U+\sum _{nm}{V}_{nm}\frac{1}{2}(X_n{X}_m+P_n{P}_m-\gamma {\delta }_{nm})$$

$${\dot{X}}_n=\sum _m{V}_{nm}{P}_m,{\dot{P}}_n=-\sum _m{V}_{nm}{X}_m,$$

$$\dot{x}$$

$$\dot{p}$$

$$J_{n,\text{cl}}=\frac{1}{2}(X_n^2+P_n^2-\gamma )$$

$$\hat{A}\mapsto ⟨\psi |\hat{A}|\psi ⟩.$$

$$|\psi ⟩=c_1|1⟩+c_2|2⟩,$$

$$\hat{H}\mapsto \frac{p^2}{2m}+U(x)+\sum _{nm}{V}_{nm}{c}_n^{\ast }{c}_m$$

$${\dot{c}}_n=-\mathrm{i}\sum _m{V}_{nm}{c}_m,$$

$$\hat{\mathcal{H}}$$

$$=H_0\hat{\mathcal{I}}+\mathbf{H}\cdot \hat{\mathbf{S}}$$

$${\hat{\sigma }}_x=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\hat{\sigma }}_y=\left(\begin{array}{cc}0& -\mathrm{i}\\ \mathrm{i}& 0\end{array}\right),{\hat{\sigma }}_z=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$$H_0$$

$$H_x$$

$$H_y$$

$$H_z$$

$$|\mathbf{u}⟩=\cos \frac{\theta }{2}{\mathrm{e}}^{-\mathrm{i}\phi /2}|1⟩+\sin \frac{\theta }{2}{\mathrm{e}}^{\mathrm{i}\phi /2}|2⟩,$$

$$S_i(\mathbf{u})\equiv ⟨\mathbf{u}|{\hat{S}}_i|\mathbf{u}⟩=\frac{1}{2}{u}_i,$$

$$u_x$$

$$u_y$$

$$u_z$$

$$J_1(\mathbf{u})$$

$$J_2(\mathbf{u})$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{S}}=\mathbf{H}\times \hat{\mathbf{S}}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{S}=\mathbf{H}\times \mathbf{S},$$

$$\hat{A}\mapsto A_{\text{Q}}(\mathbf{u})\equiv ⟨\mathbf{u}|\hat{A}|\mathbf{u}⟩,$$

$$\mathrm{tr}[\hat{A}]=\int \mathrm{d}\mathbf{u}{A}_{\text{Q}}(\mathbf{u})=\frac{1}{2\pi }{\int }_0^{2\pi }\mathrm{d}\phi {\int }_0^{\pi }\mathrm{d}\theta \sin \theta A_{\text{Q}}(\mathbf{u}),$$

$$\mathrm{tr}[\hat{A}\hat{B}]\ne \int \mathrm{d}\mathbf{u}{A}_{\text{Q}}(\mathbf{u})B_{\text{Q}}(\mathbf{u}),$$

$$A_{\text{Q}}(\mathbf{u})\equiv \mathrm{tr}[\hat{A}{\hat{w}}_{\text{Q}}(\mathbf{u})],{\hat{w}}_{\text{Q}}(\mathbf{u})=\frac{1}{2}\hat{\mathcal{I}}+\frac{1}{2}\mathbf{u}\cdot \hat{\mathbf{\sigma }},$$

Full text

I Introduction

The simulation of electronically nonadiabatic dynamics is important for a wide range of applications across physics, chemistry and biology, including the study of solar cells, vision, photosynthesis, radiation damage and many more.Tully (2012) What characterizes these nonadiabatic processes is that the electronic and nuclear motion are coupled, which causes the Born-Oppenheimer approximation to break down. The computational task is particularly challenging in condensed-phase systems, where the number of nuclear degrees of freedom is typically too large for an exact quantum solution to be feasible. Approximate methods are therefore necessary to reach a reasonable balance between accuracy and computational cost.

This has led to the ongoing development of mixed quantum-classical dynamics,Stock and Thoss (2005) in which the nuclear degrees of freedom (called the environment) are treated by classical molecular dynamics, and the evolution of the electronic states (called the subsystem) is treated quantum mechanically. It is, however, not trivial to treat the coupling between the two in a rigorous manner. Popular methods such as Tully’s fewest-switches surface hopping are computationally cheap but only heuristically motivated.Tully (1990); Subotnik, Ouyang, and Landry (2013); Martens (2019); Kapral (2015) Another simple approach is to neglect all dynamical correlations between the classical environment and the quantum subsystem, which is called the Ehrenfest or mean-field trajectory method.Ehrenfest (1927); McLachlan (1964)

To go beyond Ehrenfest dynamics, a number of methods are based on mapping the quantum-mechanical system to an equivalent problem, for which one can find a well-defined classical limit. In this way, the dynamics of the quantum and classical degrees of freedom can be treated on an equal footing. Among the most commonly used examples is the Meyer-Miller-Stock-Thoss (MMST) mapping,Meyer and Miller (1979a); Stock and Thoss (1997) which maps the $N$ electronic levels to an $N$-dimensional singly-excited harmonic oscillator, inspired by Schwinger’s theory of angular momentum.Schwinger (1965) There are many other mappings as well,Garbaczewski (1978); Blaizot and Marshalek (1978) including representations based on spin-coherent states.Radcliffe (1971); Kuratsuji and Suzuki (1980); Thoss and Stock (1999); Lucke, Mak, and Stockburger (1999); Song and Van Voorhis (2006) Each of these mappings is exact on a quantum-mechanical level, but can lead to different levels of approximation when used in a classical trajectory-based simulation.

Among the simplest trajectory-based approaches are quasiclassical methods defined in a mapping basis, which are based on an ensemble of uncoupled trajectories with no phase-dependence on the paths. These can describe electronic coherence effects and have a favourable scaling with system size, at the expense of excluding interference effects in the nuclear dynamics.Miller (2012) One can derive a quasiclassical expression for correlation functions in the MMST basis by linearizing the semiclassical propagator (a method referred to as linearized semiclassical-initial value representation, LSC-IVR).Miller (2001); Sun, Wang, and Miller (1998); *Wang1999mapping The same dynamics can also be derived by neglecting a cubic term in the quantum-classical Liouville equation (QCLE)Kapral (2006) when written in the MMST basis, an approach called the Poisson-bracket mapping equation (PBME).Kim, Nassimi, and Kapral (2008); Kelly et al. (2012) If one goes beyond quasiclassical methods by using the full QCLE or various semiclassical methods, which weight trajectories by phases, the dynamics can also include nuclear interference effects, which is necessary to predict the correct final momentum distributions in scattering theory. Herman and Kluk (1984); Stock and Thoss (1997); Sun and Miller (1997); Bonella and Coker (2001) This is however only feasible for small systems because these approaches get exponentially more difficult as the system size or simulation time increases.Stock and Thoss (2005); Kapral (2006) A possible compromise is provided by partially linearized methods that are somewhere between linearized and semiclassical methods in terms of both accuracy and efficiency. Huo and Coker (2011); Hsieh and Kapral (2012, 2013) Luckily, interference effects become less important in large condensed-phase systems due to fast decoherence. Our objective in this paper is therefore to find a new quasiclassical method with the same computational cost as the linearized methods mentioned above, but with an improved accuracy.

The main difference between the MMST and Ehrenfest formulations is that the former adds a term to the Hamiltonian which can be interpreted as a zero-point energy (ZPE) of the mapping degrees of freedom. Although this term is rigorously derived, an unfortunate consequence is that populations can become negative in individual trajectories, which means that the system evolves on inverted potentials.Bonella and Coker (2001); Coronado, Xing, and Miller (2001); Bellonzi, Jain, and Subotnik (2016) A fundamental problem is that in classical mechanics, this ZPE can flow unphysically between different degrees of freedom, referred to as the ZPE-leakage problem. Müller and Stock have reported that this problem can be mitigated by reducing the ZPE, treating it as an adjustable parameter.Müller and Stock (1999) Recently, Cotton and Miller have also used a ZPE-reduced MMST Hamiltonian in their symmetrical quasiclassical windowing approach, finding a particular value to be near optimal for most cases, while still treating it as a parameter that can be tuned for each system or window geometry.Cotton and Miller (2013, 2016)

What MMST and Ehrenfest have in common is that they replace the electronic levels by a quadratic Hamiltonian, whose classical equations of motion are equivalent to the time-dependent Schrödinger equation of the subsystem. This idea dates back to Dirac and was further investigated by Strocchi.Dirac (1927); Strocchi (1966) However, the quantum phase space of an $N$-level system has a fundamental difference from the classical phase space of a set of $N$ harmonic oscillators. In quantum mechanics the phase space is described by $N-1$ complex variables (with $2N-2$ real components), whereas the classical phase space uses $2N$ free real variables.Zhang and Feng (1995) The quantum system therefore has one degree of freedom fewer than the classical system. One way to eliminate this degree of freedom from the MMST mapping is via the Holstein-Primakoff transformation,Holstein and Primakoff (1940) but this approach is unable to describe Rabi oscillations correctly in semiclassical simulations.Thoss and Stock (1999) Based on Moyal’s phase-space theory of quantum mechanics,Moyal (1949) Stratonovich formulated a set of general properties to be fulfilled by mappings between quantum and classical mechanics.Stratonovich (1957) This approach is now known as the Stratonovich-Weyl transform and can be thought of as a discrete version of the Wigner transform. For an overview of its properties and applications in quantum optics, see Refs. Várilly and Gracia-Bondía, 1989; Klimov and Chumakov, 2009; Brif and Mann, 1999; Gadella, 1995.

In this paper, we propose a mapping of the two-level system based on the Stratonovich-Weyl transform of a spin-$\tfrac{1}{2}$ system. This leads to a family of new quasiclassical methods for approximating electronic correlation functions. We show that they are connected with previous approaches by formulating them in terms of an MMST Hamiltonian with reduced ZPE. One of the methods yields the same value of the ZPE parameter as was successfully used by Cotton and Miller.Cotton and Miller (2013) We finally investigate the accuracy of our method compared to previous approaches for spin-boson models.

II Theory

Before we present our spin-mapping approach, we briefly revisit previous mappings that are commonly used in quasiclassical calculations. We shall limit the discussion to two states only, although both MMST mapping and Ehrenfest are easily extendable to any finite number of electronic states.

The Hamiltonian of a molecular system with two electronic states in the diabatic representation is

[TABLE]

where $\hat{x}$ and $\hat{p}$ are position and momentum operators of the nuclear modes with associated mass $m$. $U$ and $\hat{V}$ are the state-independent and state-dependent potentials. All methods in this paper will employ the classical-path approximation, in which one replaces the nuclear operators by classical phase-space variables, $\hat{x},\hat{p}\mapsto x,p$. What distinguishes the methods is their treatment of the electronic operators and in particular the coupling between the electronic states and the nuclei. Throughout this paper we set $\hbar=1$.

II.1 MMST mapping and Ehrenfest dynamics

The so-called Meyer-Miller-Stock-Thoss (MMST) mappingMeyer and Miller (1979a); Stock and Thoss (1997) maps the electronic basis states to singly-excited multidimensional harmonic oscillator states and hence electronic operators to ladder operators

[TABLE]

A transition from state $|m\rangle$ to $|n\rangle$ thus corresponds to moving the single excitation from oscillator $m$ to oscillator $n$. The ladder operators can be written in terms of electronic position and momentum variables as $\hat{a}_{n}=\frac{1}{\sqrt{2}}(\hat{X}_{n}+\mathrm{i}\hat{P}_{n})$, so that

[TABLE]

where we used the commutation relation $[\hat{X}_{n},\hat{P}_{m}]=\mathrm{i}\delta_{nm}$. Up until this point, the mapping is exact. In the MMST method one now defines a classical limit by removing the hats from the $\hat{X}_{n},\hat{P}_{n}$ operators and treating them as classical canonical variables, which is analogous to the classical-path approximation of the nuclear variables. In other words, electronic and nuclear degrees of freedom are treated on the same footing. The mapped Hamiltonian,

[TABLE]

in which $V_{nm}=\langle n|\hat{V}(x)|m\rangle$, contains a term $\sum_{nm}V_{nm}\frac{1}{2}\delta_{nm}=\frac{1}{2}\mathrm{tr}\,\hat{V}(x)$ that can be interpreted as an electronic zero-point energy (ZPE).Müller and Stock (1999) In this expression we have assumed that $\Delta$ is chosen to be real.

A problem with the MMST method is that in a classical simulation, zero-point energy can flow unrestrictedly between different states, which should not be allowed in a quantum-mechanical simulation.Stock and Thoss (2005) This is referred to as the ZPE-leakage problem. Stock and Müller have found that the ZPE-leakage problem can be mitigated by reducing the overall ZPE.Müller and Stock (1999) They introduced a parameter $\gamma$ into the mapped Hamiltonian

[TABLE]

which creates a family of methods where $\gamma=1$ corresponds to standard MMST mapping and $\gamma=0$ to Ehrenfest dynamics, as we will see below. Sometimes $\gamma$ is treated as a parameter that is tuned according to some criterion, for example so that average adiabatic populations stay positive.Müller and Stock (1999) In this paper we shall derive a new type of mapping and find that it leads to a Hamiltonian of the form in Eq. (6), but a different phase space.

The equations of motion generated by this Hamiltonian are

[TABLE]

From Eq. (7c), we see that each potential energy surface, $V_{n}$, contributes to the nuclear force proportionally to

[TABLE]

which we call the classical population of state $n$. The initial $\{X_{n},P_{n}\}$ variables are typically sampled from either a Wigner distribution or a Dirac delta distribution, depending on the model one wants to simulate. In the first case individual trajectories can have a total population different from one, $J_{1,\text{cl}}+J_{2,\text{cl}}\neq 1$. In both cases populations can even be negative. Negative populations means that the system evolves on inverted potentials, which can lead to problems in particular for systems with steep potentials.Coronado, Xing, and Miller (2001); Bellonzi, Jain, and Subotnik (2016) It is possible to eliminate $\gamma$ from the equations of motion by choosing a traceless form of $\hat{V}$, so that $\gamma$ only affects the initial distribution and the observables, but does not appear directly in the equations of motion. In other words the MMST mapping is not unique, but depends on a particular choice of splitting between $U$ and $\hat{V}$ (unless the initial conditions are chosen to fix the total population to be one).

A widely used method that gives similar dynamics is the Ehrenfest or mean field trajectory method,Ehrenfest (1927); McLachlan (1964) in which one replaces the electronic operators by their expectation values:

[TABLE]

Let us write an arbitrary electronic state $|\psi\rangle$ in terms of diabatic basis states as

[TABLE]

where $c_{n}=\frac{1}{\sqrt{2}}(X_{n}+\mathrm{i}P_{n})$ are the complex coefficients. Note that if $|\psi\rangle$ is normalized and has a fixed global phase, only two of the four real variables $X_{1},X_{2},P_{1},P_{2}$ are independent (this fact will be used in the next section). Using these coefficients, it is clear that the Ehrenfest Hamiltonian

[TABLE]

is equivalent to the MMST Hamiltonian in Eq. (6) with $\gamma=0$.Stock and Thoss (2005) The electronic equations of motion (7a) are equivalent to the solution of the time-dependent Schrödinger equation for the subsystem

[TABLE]

which is also true for any value of $\gamma$. Because $\gamma=0$, Ehrenfest dynamics does not suffer from the inverted potential problem. Instead, one finds the one-trajectory problem: if the system is initially in state $|1\rangle$, that is $|c_{1}|^{2}=\frac{1}{2}(X_{1}^{2}+P_{1}^{2})=1$ and $|c_{2}|^{2}=\frac{1}{2}(X_{2}^{2}+P_{2}^{2})=0$, then for each initial nuclear position and momentum, the dynamics will consist of a single unique trajectory, which on its own cannot capture the correct quantum dynamics. Historically Meyer and Miller solved this problem by introducing a Langer correction,Meyer and Miller (1979a) which leads to the MMST Hamiltonian with $\gamma=1$ in Eq. (5). This is not the only possible solution though; any $\gamma\neq 0$ solves the one-trajectory problem.

Throughout this paper, we will repeatedly refer back to the MMST mapping and Ehrenfest dynamics. Before we present our new mapping, let us look at how the Ehrenfest method can be rephrased in terms of a spin-$\tfrac{1}{2}$ system.

II.2 Equivalence of a two-level system and a spin-1/2 particle in a magnetic field

It is well known that the Hamiltonian (and in general any Hermitian operator of the electronic states) can be decomposed into a basis of spin operators and the identity:

[TABLE]

where $\hat{S}_{i}=\frac{1}{2}\hat{\sigma}_{i}$ ($i\in\{x,y,z\}$) and the Pauli spin matrices are

[TABLE]

Such a decomposition was used already in the derivation of the spin-matrix method by Meyer and Miller.Meyer and Miller (1979b) The explicit relations between the quantities in Eqs. (1) and (13) are

[TABLE]

Note that in this way, it does not matter how the initial split between $U$ and $\hat{V}$ is chosen. For simplicity we again assume that $\Delta$ is real, so that $H_{y}=0$.

The Hamiltonian in equation (13) is identical to that of a spin-$\tfrac{1}{2}$ particle in a magnetic field $-\bm{H}$ (if we choose its gyromagnetic ratio to be one). We emphasize that this spin is only a theoretical tool and has no relation to any real spins of the physical system. Since the Hamiltonians in Eqs. (1) and (13) are in one-to-one correspondence, any classical phase space of spin-$\tfrac{1}{2}$ particles will also be a phase space of the general two-level system.

In order to construct a classical phase space, we use spin coherent states defined asRadcliffe (1971)

[TABLE]

where $\bm{u}$ denotes a unit vector defined by the spherical coordinates $\theta,\varphi$ and the states are normalized such that $\braket{\bm{u}}{\bm{u}}=1$. Note that this defines an arbitrary state using two real variables, in contrast with the MMST mapping that is defined with four real variables. The expectation values of the spin operators $\hat{S}_{i}$ in this state are

[TABLE]

where

[TABLE]

Thus $S_{i}(\bm{u})$ are the Cartesian coordinates of a sphere with radius $r=1/2$, see Figure 1. One can think of them as orthogonal functions analogous to p-orbitals. Likewise the populations are

[TABLE]

and it is clear that $J_{1}(\bm{u})+J_{2}(\bm{u})=1$ is always fulfilled. Note that on the “poles” of the sphere along the $S_{z}$-axis ($\theta=0,\pi$), all population is concentrated in one of the electronic basis states, as shown in Figure 1.

Let us now briefly discuss the dynamics of the spin. It is known from quantum mechanics that the vector of spin operators follows the (Heisenberg) equation of motion

[TABLE]

Replacing all operators with their expectation values $\bm{S}=(S_{x}(\bm{u}),S_{y}(\bm{u}),S_{z}(\bm{u}))$ gives

[TABLE]

which describes precession of the spin vector around the magnetic field, see Figure 1. These equations of motion are also equivalent to the time-dependent Schrödinger equation of the subsystem in Eq. (12).111This can be shown by setting $2S_{x}=\text{Re}[c_{1}^{*}c_{2}]$, $2S_{y}=\text{Im}[c_{1}^{*}c_{2}]$, $2S_{z}=\tfrac{1}{2}(|c_{1}|^{2}-|c_{2}|^{2})$. The dynamics preserves the magnitude $|\bm{S}|$, i.e. the radius of the sphere. It is the $x$-component of the field, $H_{x}=2\Delta$, that drives population transfer between the two states, whereas the $z$-component, $H_{z}=V_{1}-V_{2}$, preserves all populations and only contributes with a relative phase.

What we have presented in this section is an alternative formulation of the subsystem dynamics in the Ehrenfest method, which also appears in Ref. Meyer and Miller, 1979b. In the next section, we will generalize the above procedure to allow for a phase-space representation of correlation functions, leading to an accurate quasiclassical methodology.

II.3 Stratonovich-Weyl transform of the spin-1/2 system

In the previous section we replaced spin operators by their expectation values. Let us for a general operator $\hat{A}$ associate this procedure with the mapping

[TABLE]

which in literature is known as the Q-function.Klimov and Chumakov (2009) The Q-function allows us evaluate a quantum-mechanical trace over the subsystem, $\mathrm{tr}[\hat{A}]$, as an integral over a classical phase space:

[TABLE]

where the integration measure is $\mathrm{d}\bm{u}=\tfrac{1}{2\pi}\mathrm{d}\varphi\mathrm{d}\theta\sin\theta$ (which is normalized to be consistent with $\mathrm{tr}[\hat{\mathcal{I}}]=2$). However, in order to calculate correlation functions, we need a phase-space representation for products of operators of the form $\mathrm{tr}[\hat{A}\hat{B}]$. The Q-function is not enough for this purpose since

[TABLE]

because of the uncertainty property $\langle\bm{u}|\hat{A}\hat{B}|\bm{u}\rangle\neq\langle\bm{u}|\hat{A}|\bm{u}\rangle\langle\bm{u}|\hat{B}|\bm{u}\rangle$. The only case when equality holds is when $\hat{A}$ or $\hat{B}$ is (proportional to) the identity operator $\hat{\mathcal{I}}$. This indicates that we need to treat the spin operators differently from the identity. Note that the idea of such a separation was recently investigated in the MMST representation.Saller, Kelly, and Richardson (2019) Now we will show how the uncertainty property of the spin operators can be included directly in the construction of the phase space. First, note that Eq. (22) can written equivalently as the result of the quantum-mechanical trace:

[TABLE]

Next, introduce the P-function222The Q- and P-functions are sometimes called Berezin’s covariant and contravariant symbols Berezin (1975)

[TABLE]

which together with the Q-function allows us to represent $\mathrm{tr}[\hat{A}\hat{B}]$ as

[TABLE]

Interestingly, one can also define a self-dual W-function

[TABLE]

which fulfils

[TABLE]

We will refer to the formulas (27) and (29) as the traciality property of the mapping. It is easy to show (using the special case $\hat{B}=\hat{\mathcal{I}}$) that the P- and W-functions also fulfil equations of the form of Eq. (23).

It is natural to think of the mapping to the phase-space representation in Eq. (28), and the pair of Eqs. (25) and (26), as a discrete version of the Wigner transform. In literature this is known as the Stratonovich-Weyl transform.Stratonovich (1957); Brif and Mann (1999) It has been applied to various problems in quantum opticsKlimov and Chumakov (2009) but to our knowledge not yet in the context of nonadiabatic molecular dynamics. To summarize, let us write the Stratonovich-Weyl kernels $\hat{w}_{s}(\bm{u})$ ($s\in\{\text{Q},\text{P},\text{W}\}$) collectively as

[TABLE]

where we introduced the spin radius $r_{s}$, defined in Table 1 for each $s$. As simple examples, we have

[TABLE]

and all other operators can be built as linear combinations of these. In other words, the Q-function of the quantum spin vector $\hat{\bm{S}}$ defines the same sphere with radius $r_{\text{Q}}=1/2$ that we saw in Figure 1, whereas the corresponding W- and P-functions define larger spheres of radii $r_{\text{W}}=\sqrt{3}/2$ and $r_{\text{P}}=3/2$, respectively.

In quantum mechanics, the squared magnitude of the spin is

[TABLE]

Its eigenvalues are equal to $r_{\text{W}}^{2}$, i.e. $r_{\text{W}}$ is the quantum magnitude of the spin vector. Thus the W-sphere in Figure 2 can be interpreted as a classical phase space of spin. Also the Q- and P-functions give the correct quantum magnitude of spin as long as they are used in pairs ($r_{\text{Q}}r_{\text{P}}=3/4$) as in Eq. (27).

Let us now take a step back and see what consequences this has for the general two-level system. Populations are in the spin-mapping space represented as

[TABLE]

and $J_{1,s}+J_{2,s}=1$ for all $s$. As a result, the points on the phase-space sphere where all population is concentrated to one eigenstate are no longer on the poles but are identified with “polar circles” defined by $\cos\theta_{\text{c}}=\pm 1/(2r_{s})$, as indicated for the W-function in Figure 2. In the region “north” of the “arctic” circle one finds $J_{1,s}>1$ and $J_{2,s}<0$, and vice versa to the “south” of the “antarctic” circle. In the Q-case, $J_{n,\mathrm{Q}}=1$ only on the poles as in Figure 1, so that the are no points with negative populations, and in the P-case the situation is similar to in Figure 2, but the circles with $J_{n,\text{P}}=1$ are closer to the equator, like tropical circles.

This has consequences for the dynamics in the three spin-mapping representations. The Hamiltonian in equation (13) maps to

[TABLE]

In Appendix A we derive the equations of motion using a spin-mapping representation of the QCLE (within an uncoupled-trajectory approximation), which allows us describe the dynamics of both the subsystem and the environment.Kapral (2006) The resulting equations of motion are

[TABLE]

which can be shown to conserve $H_{s}(\bm{u})$. Importantly the subsystem equations are exact and equal for all $s$, and preserve the magnitude $|\bm{u}|$. Thus the only dependence on $s$ is in the nuclear equations of motion. Firstly, the contribution from the electronic coupling, $\Delta$, scales by $r_{s}$. Secondly, the region of phase space corresponding to inverted potential have different size. In this aspect, $H_{\text{Q}}$ is interesting because it never gives inverted potentials (it is the same Hamiltonian as is used in Ehrenfest dynamics). The W-representation on the other hand has the attractive property of being self-dual.

II.4 Comparison with previous methods

Let us now connect spin mapping to the previous methods introduced in section II.1. We introduce the coordinate transformation

[TABLE]

where $\bm{\tilde{\sigma}}=(\tilde{\sigma}_{x},\tilde{\sigma}_{y},\tilde{\sigma}_{z})$ is the MMST representation of the Pauli matrices. Using $|\bm{u}|^{2}=1$ one can show that

[TABLE]

i.e. the new variables are also constrained to a sphere (a 3-dimensional hypersphere embedded in $\mathbb{R}^{4}$) but with a different radius $R=\sqrt{4r_{s}}$. The transformation used together with Eq. (37) brings the Hamiltonian in Eqs. (34) to the form

[TABLE]

If $\{X_{n},P_{n}\}$ are identified as canonical variables of the subsystem, we get equations of motion equivalent to Eq (35), and $R$ is constant.Therefore this Hamiltonian is equivalent to the ZPE-reduced MMST Hamiltonian in equation (6) with $\gamma=2r_{s}-1$. In particular $H_{\text{W}}$ has $\gamma=\sqrt{3}-1\approx 0.732$, which is the value Cotton and Miller gave a heuristic argument for in Ref. Cotton and Miller, 2013, also based on an analogy with spin (note that they define $\gamma$ as one half times our definition). It is also close to values that Müller and Stock found optimal to reproduce the correct level density of various models,Müller and Stock (1999) and Golosov and Reichman have showed that it gives correct short-time dynamics up to seventh order in time for some observables, compared to fifth order with standard MMST mapping.Golosov and Reichman (2001) Note however that in both the last two references, the mapping Hamiltonian is defined as $\sum_{nm}H_{nm}\tfrac{1}{2}(X_{n}X_{m}+P_{n}P_{m}-\gamma\delta_{nm})$ where $H_{nm}=(\tfrac{p^{2}}{2m}+U)\delta_{nm}+V_{nm}$. This is different from our approach, which treats terms proportional to the identity differently than other operators.

The split between $U$ and $\hat{V}$ is therefore unambiguous in our approach. This is clear from Eq. (34) that only multiplies the traceless part of $\hat{V}$ by $r_{s}$. In this way, when $\hat{V}$ goes to $\hat{\mathcal{I}}$, the mapped Hamiltonian smoothly reduces to the single-state Hamiltonian. In LSC-IVR however, different choices of the split between $U$ and $\hat{V}$ generally lead to different results since the total population $J_{1,\text{cl}}+J_{2,\text{cl}}$ is allowed to deviate from one. A further difference is that Wigner initial conditions in MMST mapping corresponds to a Gaussian distribution in four dimensions with unbounded variables. In our case, the phase-space variables are bounded by the constraint to a sphere. This leads us to identify the most important difference: since the entire spin-mapping space is isomorphic to the physical space, there is no need to introduce projections onto a physical subspace as in the MMST mapping.

Compared to another recent spin mapping by Cotton and Miller,Cotton and Miller (2015) our method maps two levels to one spin-$\tfrac{1}{2}$ particle, whereas Cotton and Miller map to two spin-$\tfrac{1}{2}$ particles. Thus their mapping space has, just like in standard MMST mapping, more degrees of freedom than the underlying quantum system.

A link between the MMST mapping and spin coherent states was found already by Stock and Thoss,Thoss and Stock (1999) who related their work to the spin-coherent state propagator by Suzuki.Suzuki (1983) They expressed a semiclassical propagator in the MMST representation as an integral over the coordinates of the spin coherent state333Stock and Thoss used a complex coherent state variable $\mu$ which is effectively a stereographic projection of our $\bm{u}$ vector onto a complex plane. and what is in our notation the radius $R$. They then performed the approximation of fixing $R$ in order to relate it to Suzuki’s propagator. In our work we directly derive the spin coherent state representation, which naturally leads to a fixed $R$, and additionally opens up the path towards defining correlation functions via the Q-, P- and W-functions, as we describe in the next section.

II.5 Correlation functions from the various methods

We will now use our phase-space representation from the previous sections to approximate correlation functions of the type

[TABLE]

where $\mathrm{Tr}$ denotes a full quantum-mechanical trace, $\hat{\rho}$ is a density matrix normalized to give $\mathrm{Tr}[\hat{\rho}]=1$, and $\hat{B}(t)=\mathrm{e}^{\mathrm{i}\hat{H}t}\hat{B}\mathrm{e}^{-\mathrm{i}\hat{H}t}$. In this paper we will only consider the case when $\hat{A},\hat{B}$ are electronic operators and the density matrix is state-independent, $\hat{\rho}=\hat{\rho}_{\text{nuc}}\otimes\hat{\rho}_{\text{el}}$, where $\hat{\rho}_{\text{el}}=\frac{1}{2}\hat{\mathcal{I}}$ (the normalization ensures that $\mathrm{tr}[\hat{\rho}_{\text{el}}]=1$). The initial projection onto one of the electronic states is defined by $\hat{A}$.

We will use a Wigner distribution of the $F$ nuclear degrees of freedom,

[TABLE]

and a quasiclassical propagator analogous to classical Wigner dynamics:

[TABLE]

where the dual indices ${\bar{s}}$ are specified in Table 1. The trajectory $\bm{u}(t)$ is generated by $H_{s}$, with the same index as the operator at time zero. (We also considered the opposite index choice $A_{\bar{s}}B_{s}$ but found that it gives less accurate results at $t>0$.)

We integrate over $\frac{\mathrm{d}\bm{u}}{2}=\frac{\mathrm{d}\varphi}{4\pi}\mathrm{d}\theta\sin\theta$ using a Monte Carlo scheme, drawing samples uniformly from the unit sphere $|\bm{u}|^{2}=1$. Using Eq. (37), this can be done directly in the $\{X_{n},P_{n}\}$ variables by sampling from a 3-sphere with radius $R=\sqrt{4r_{s}}=\sqrt{2(\gamma+1)}$. An easy way to do this in practice is to sample $\{X_{n},P_{n}\}$ independently from a standard normal distribution and multiply them with a common factor to fulfil Eq. (37).

Eq. (41) gives us three new methods classified by $s$: we will call them the Q-, P- and W-methods. To evaluate their performance, we will compare them to two standard quasiclassical methods in the MMST-basis. These differ in how the system is projected onto the physical subspace:Kelly et al. (2012) the first, LSC-IVR, has projection operators both at initial and final times,Miller (2001); Sun, Wang, and Miller (1998) and the second, PBME, has a projection only at initial times.Kim, Nassimi, and Kapral (2008) As a general formula we write

[TABLE]

where $\mathsf{X}=(X_{1},X_{2})$, $\mathsf{P}=(P_{1},P_{2})$. The explicit expressions for $\rho_{\text{el}}(\mathsf{X},\mathsf{P})$ are given in Table 2 together with the appropriate observables in the various methods. For a general operator $\hat{A}=A_{0}\hat{\mathcal{I}}+\tfrac{1}{2}\bm{A}\cdot\hat{\bm{\sigma}}$, the MMST mapping representations used in the standard formulations of LSC-IVR and PBME are defined as

[TABLE]

where $A_{\text{SEO}}$ is used whenever one projects the harmonic oscillators to the physical subspace, and $\phi=16\exp(-R^{2})$. The special case of a population operator, $\hat{A}=\ket{n}\!\bra{n}=(\hat{\mathcal{I}}+\hat{\sigma}_{z})/2$ gives

[TABLE]

Similarly we can also write the Stratonovich-Weyl observables in terms of $\{\mathsf{X},\mathsf{P}\}$ using Eq. (36), which gives

[TABLE]

Note that $A_{\text{SW}}$ gives a standard expression for a population operator, $A_{\text{SW}}=\tfrac{1}{2}(X_{n}^{2}+P_{n}^{2}-\gamma)=J_{n,\text{cl}}$, whereas $B_{\overline{SW}}$ has not previously been used in quasiclassical mapping dynamics.

All methods presented so far treat initial conditions via a weighting procedure. We will contrast these to methods that use focused initial conditions, which are used in Ehrenfest dynamics and sometimes also with the MMST mapping. By focused we mean that to start in for example $\hat{A}=|1\rangle\langle 1|$, an Ehrenfest calculation would use the initial conditions $|c_{1}|^{2}=1$ and $|c_{2}|^{2}=0$, which in the MMST representation is

[TABLE]

where $\mathcal{N}=\pi^{-2}$ (although the normalization is not explicitly needed when evaluating with a simple Monte Carlo procedure). In a similar way, we define a focused initial condition for the W-representation by sampling $\bm{u}$ from the northern polar circle in Figure 2, which corresponds to Eq. (47) with $\gamma=\sqrt{3}-1\approx 0.732$. A focused Q-distribution is identical to Ehrenfest dynamics. Unlike the approaches in Table 2, we cannot derive the focused methods in a rigorous way from the Stratonovich-Weyl formalism but instead follow the same line of reasoning as Refs. Bonella and Coker, 2003; Müller and Stock, 1999.

Note that when using focused initial conditions, the Stratonovich-Weyl operators reduce to the usual Ehrenfest (MFT) and MMST operators. For $\gamma=0$

[TABLE]

and similarly for $\gamma=1$ one has $1=(R^{2}-2)/2$, so that $B_{\text{SW}}=B_{\text{wig}}$.

III Results and discussion

To test the accuracy of our method, we have applied it on a spin-boson model,Garg, Onuchic, and Ambegaokar (1985) which is a standard benchmark problem consisting of two electronic levels coupled to a harmonic bath of nuclear modes. Note however that all quasiclassical mapping approaches can be applied to more general problems with anharmonic potentials. The potentials defined in Eq. (1) are of the form

[TABLE]

Here $\varepsilon$ is the energy bias and $\Delta$ the constant diabatic coupling between the two electronic states. The bath consists of $F$ nuclear modes, each with frequency $\omega_{j}$ and vibronic coupling coefficient $c_{j}$. In this form the matrix $\hat{V}$ is already traceless. We choose units such that $\hbar=1$.

The so-called spectral density of the bath, $J_{\text{bath}}(\omega)$, determines the distribution of nuclear frequencies. Among the most common forms of this function is the Ohmic bath

[TABLE]

where $\omega_{\text{c}}$ is called the characteristic frequency and $\xi$ the Kondo parameter, which determines the strength of the friction. Numerically one uses a discretization

[TABLE]

with a finite $F$. In this work we used the discretization scheme employed in Ref. Craig and Manolopoulos, 2005.

To sample the initial nuclear variables $\{x_{j},\;p_{j}\}$, we used the thermal Wigner distribution

[TABLE]

with $\alpha_{j}=\tanh\left(\frac{1}{2}\beta\omega_{j}\right)$ and $\beta=1/k_{\text{B}}T$.

To facilitate comparison with other recent papers, we have used the same parameter settings as in Ref. Cotton and Miller, 2016, which are given in Table 4. Note that we define our initial nuclear distribution to be centred around the minimum of $U$, which for these systems gives practically identical results to when starting around $U+V_{1}$ as in Ref. Cotton and Miller, 2016 (except for subtle changes in the strong-bath case). In our calculations we used $F=100$ nuclear modes and a timestep of $0.01\Delta^{-1}$. To ensure very tight convergence, we used $10^{5}$ trajectories for the methods with focused initial conditions and $10^{6}$ trajectories for the others. However, like other mapping approaches, $10^{3}-10^{4}$ is typically enough to observe the correct physical behaviour with only a little numerical noise.

The results for the linearized methods specified in Table 2 are compared in Figure 3 with numerically exact path-integral calculations (QUAPI).Makarov and Makri (1994) One sees that the Q-method captures both correct oscillations and correct long-time limits in all systems, except for the strong bath systems (e-f) where nuclear quantum effects are expected to be important. Also the W-method yields fairly accurate results with correct long-time limits, but slightly underestimates the oscillation amplitudes. This too quick decoherence is even more pronounced in the P-method. We conclude that a minimal $r_{s}$ in the Hamiltonian is optimal, at least for the problems considered here. The traditional MMST-based methods with double (LSC-IVR) and single (PBME) projection perform well in some cases but can be completely wrong in others, especially for final populations in asymmetric problems.

It should be pointed out that a recently published trickSaller, Kelly, and Richardson (2019) to treat the identity operator in LSC-IVR or PBME gives essentially the same results as our Q-method. There is however no direct connection between the two, since Ref. Saller, Kelly, and Richardson, 2019 infers no restriction on the radius $R$. Our results suggest that the separation of the identity is what leads to correct long-time limits, whereas the choice of $\gamma$ controls the decoherences. It should also be noted that other methods like SQC, FBTS and PLDM, which go beyond the simplest quasiclassical approximation, also perform well for the spin-boson model.Cotton and Miller (2016); Hsieh and Kapral (2013); Huo and Coker (2012)

When computing correlation functions of off-diagonal operators like $C_{\sigma_{x}\sigma_{x}}$ (see Fig. 4), both Q and LSC-IVR are essentially exact. This is a general observation for all correlation functions of traceless observables. Note that our method is able to compute all electronic correlation functions from the same set of trajectories, rather than different ones for each initial operator $\hat{A}$.

For the methods with focused initial conditions in Table 3 and Figure 5, the focused W-method is clearly superior to both Ehrenfest and standard MMST ($\gamma=1)$. It uses $\gamma=\sqrt{3}-1\approx 0.732$, which is close to the values that Müller and Stock found optimal in their simulations using action-angle (focused) initial conditions.Müller and Stock (1999) Both Ehrenfest ($\gamma=0$, which is the same as a focused Q-method) and standard MMST $(\gamma=1)$ with focused initial conditions are much less accurate for asymmetric problems. We also tried a focused P-method $(\gamma=2)$, but this deviates badly from the exact results and is not shown in the plots. This is expected from the trend seen in the asymmetric problems that a larger $\gamma$ leads to a lowering of the curve $C_{1\sigma_{z}}(t)$.

Our focused W-method results are of similar quality as with Cotton and Miller’s symmetric quasi-classical (SQC) windowing approach, which interestingly uses the same value of $\gamma$ in their method involving square windows.Cotton and Miller (2013, 2016) Compared to SQC, however, our method does not require any choice of window functions, and no trajectories are discarded from the averaging procedure. Also note that while the spin mapping method of Ref. Cotton and Miller, 2015 (which maps to two spins instead of one) has been found to give worse results than SQC defined in the MMST-basis, our spin mapping is an improvement compared to standard focused MMST. The W-function also appears to be more accurate than FBTS or PLDM with focused initial conditions.Hsieh and Kapral (2013); Huo and Coker (2012)

We have no proof that our approach will always be better than (or as good as) these alternative methods. However, since it avoids the problem of leaving the physical subspace, and otherwise relies on similar quasiclassical approximations, it can be expected to lead to an improvement. The spin-mapping theory presented in this paper is limited to only two electronic levels, but since the Stratonovich-Weyl transform exists also for higher dimensional spaces, we believe it should be possible to extend to more states. Like all quasiclassical mapping approaches, it is expected to give the wrong asymptotic limit for trajectories emerging from a nonadiabatic crossing, which will generally evolve on a superposition of electronic states instead of on a single state.Miller (2001) It will also suffer from the usual classical-Wigner problem of ZPE-leakage between nuclear modes. For studies of ultrafast dynamics this is typically an acceptable approximation.

Even though the Ehrenfest method performs poorly for these systems, its accuracy can be improved by combining it with the generalized quantum master equation.Kelly, Brackbill, and Markland (2015); *Kelly2016master In this way, observables are calculated via a memory kernel that generally decays much faster than the observables themselves, so that one can make the most out of the short-time dynamics. It would be interesting to try the same trick for our method, to compute long-time dynamics from short simulations and (potentially) further improve the accuracy.

All together, we find that the most accurate method for these spin-boson systems is the Q-method without focusing. This only fails to quantitatively reproduce the exact result for the strong bath, where nuclear quantum effects are expected to be important. Further tests will be carried out to discover if this behaviour carries over to other systems.

IV Conclusions and future prospect

We have presented a family of three new mappings for two-level molecules, based on the Stratonovich-Weyl transform of spin-$\tfrac{1}{2}$ systems. Just like the MMST mapping, these are exact on a quantum-mechanical level. In a quasiclassical treatment, both the spin and MMST mappings give exact results for a bare two-level system, but the two approaches have different accuracy when there is a coupling to an environment.

By using the traciality property of the Stratonovich-Weyl transform [Eqs. (27) and (29)] we can approximate electronic correlation functions in a manner similar to classical Wigner approaches. The Q-method seems particularly promising because it gives accurate results and avoids the possibility of inverted potentials. The self-dual W-representation also appears to be very useful, in particular in implementations using focused initial conditions.

Unlike most other quasiclassical methods, our scheme only requires two free real variables to represent the two-level system, so that it has the same number of (electronic) degrees of freedom as the quantum problem. The Hamiltonian is unique, without ambiguity in the splitting between the state-dependent and state-independent potentials. Via a coordinate transformation we show that the same dynamics can be generated with the usual quadratic MMST Hamiltonian, leading to numerically stable linear equations of motion. Our approach is therefore closely related to the MMST mapping but with a couple of important differences. Firstly, the total population is always fixed to one, because the $\{\mathsf{X},\mathsf{P}\}$ variables are constrained to a three-dimensional sphere. Secondly, it uses a different value of the ZPE-parameter, $\gamma$, which scales the operators in a way that is unexpected from an MMST point of view. Because of this, it is not enough to simply fix the total population in an MMST simulation $(\gamma=1)$, since one would still need to rescale the traceless part of the $\hat{B}$ operator by $3/4$.

A fundamental advantage compared to MMST is that the system cannot leave the physical subspace, so there is no need for additional projections. We believe that this fact will be particularly useful for the development of a nonadiabatic version of ring-polymer molecular dynamics (RPMD),Ananth (2013); Chowdhury and Huo (2017); Richardson and Thoss (2013); Richardson et al. (2017) which would add nuclear quantum effects to the simulations, as well as to allow sampling from an exact thermal distribution.

Acknowledgements.

We are grateful for Michael Thoss for valuable comments on the manuscript. J.E.R. is supported by the Hans H. Günthard scholarship, and both authors acknowledge support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network.

Appendix A Derivation of the equations of motion from the QCLE

The quantum-classical Liouville equationKapral and Ciccotti (1999) is an equation of motion for mixed quantum-classical dynamics that, in its full implementation, is exact for spin-boson models.Kapral (2015) It has previously been used to analyse both surface hopping methodsKapral and Ciccotti (1999); Horenko et al. (2002); Kelly and Markland (2013) and methods based on the MMST representation.Kim, Nassimi, and Kapral (2008); Nassimi, Bonella, and Kapral (2010); Kelly et al. (2012) After performing a partial Wigner transform of the nuclear coordinates, it states that the time evolution of an operator $\hat{B}=\hat{B}(x,p)$ follows the equation

[TABLE]

where $[\cdot,\cdot]$ indicates the commutator and $\{\hat{A},\hat{B}\}=\frac{\partial\hat{A}}{\partial x}\frac{\partial\hat{B}}{\partial p}-\frac{\partial\hat{A}}{\partial p}\frac{\partial\hat{B}}{\partial x}$ is the Poisson bracket in the nuclear phase space. We shall transform this to the Stratonovich-Weyl representation by using its definition

[TABLE]

for a general decomposition $\hat{B}=B_{0}\hat{\mathcal{I}}+\tfrac{1}{2}\bm{B}\cdot\hat{\bm{\sigma}}$, together with the inversion formula

[TABLE]

The left-hand side of Eq. (55) is directly written as

[TABLE]

For the right-hand side we first need the SW-transform of a product. Repeated use of Eqs. (56)–(57) gives

[TABLE]

Using the result

[TABLE]

together with $r_{s}r_{{\bar{s}}}=3/4$, we can (after some algebra) express the commutator as

[TABLE]

where we used the integrals

[TABLE]

Similarly for the Poisson bracket we get

[TABLE]

We will omit the terms on the last line in order to find a solution of uncoupled trajectories. Such an approximation is used also to derive PBME,Kelly et al. (2012) but compared to the MMST-basis used there, the Stratonovich-Weyl representation has the advantage that it requires no projectors. It is therefore not possible for the system to leave the physical subspace, even without the last line in Eq. (63). A detailed analysis shows that the omitted term integrates to

[TABLE]

which interestingly vanishes for $r_{s}=\sqrt{3}/2=r_{\text{W}}$, indicating a potential advantage of the symmetric W-function. It is nevertheless an approximation to omit this term in an uncoupled trajectory-solution, since it does not vanish for all higher derivatives.

Combining all the terms of Eq. (55) together, we get

[TABLE]

The easiest solution to this equation is obtained by setting the integrands on each side to be equal. Using the chain rule to expand the total derivative of $B_{s}(\bm{u},x,p)$ as

[TABLE]

and comparing the expressions term by term, we obtain the equations of motion (35). Therefore an ensemble of these uncoupled trajectories can be used to obtain quasiclassical dynamics as described in the main text.