Spin-mapping approach for nonadiabatic molecular dynamics

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.

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.

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. 1 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, 2 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. [3][4][5][6] 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. 7,8 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) a) Electronic mail: johan.runeson@phys.chem.ethz.ch b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch mapping, 9,10 which maps the N electronic levels to an Ndimensional singly-excited harmonic oscillator, inspired by Schwinger's theory of angular momentum. 11 There are many other mappings as well, 12,13 including representations based on spin-coherent states. [14][15][16][17][18] 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. 19 One can derive a quasiclassical expression for correlation functions in the MMST basis by linearizing the semiclassical propagator (a method referred to as linearized semiclassicalinitial value representation, LSC-IVR). 20,21 The same dynamics can also be derived by neglecting a cubic term in the quantum-classical Liouville equation (QCLE) 22 when written in the MMST basis, an approach called the Poisson-bracket mapping equation (PBME). 23,24 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. 10,[25][26][27] This is however only feasible for small systems because these approaches get exponentially more difficult as the system size or simulation time increases. 2,22 A possible compromise is provided by partially linearized methods that are somewhere between linearized and semiclassical methods in terms of both accuracy and efficiency. [28][29][30] 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. 27,31,32 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. 33 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. 34,35 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. 36,37 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. 38 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, 39 but this approach is unable to describe Rabi oscillations correctly in semiclassical simulations. 16 Based on Moyal's phase-space theory of quantum mechanics, 40 Stratonovich formulated a set of general properties to be fulfilled by mappings between quantum and classical mechanics. 41 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. 42-45. In this paper, we propose a mapping of the two-level system based on the Stratonovich-Weyl transform of a spin- 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. 34 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 iŝ wherex andp are position and momentum operators of the nuclear modes with associated mass m. U and V are the state-independent and state-dependent potentials. All methods in this paper will employ the classicalpath approximation, in which one replaces the nuclear operators by classical phase-space variables,x,p → 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 = 1.

A. MMST mapping and Ehrenfest dynamics
The so-called Meyer-Miller-Stock-Thoss (MMST) mapping 9,10 maps the electronic basis states to singlyexcited multidimensional harmonic oscillator states and hence electronic operators to ladder operators A transition from state |m to |n 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â n = 1 √ 2 (X n + iP n ), so that where we used the commutation relation [X n ,P m ] = iδ nm . Up until this point, the mapping is exact. In the MMST method one now defines a classical limit by removing the hats from theX n ,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, (5) in which V nm = n|V (x)|m , contains a term nm V nm 1 2 δ nm = 1 2 trV (x) that can be interpreted as an electronic zero-point energy (ZPE). 33 In this expression we have assumed that ∆ 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. 2 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. 33 They introduced a parameter γ into the mapped Hamiltonian which creates a family of methods where γ = 1 corresponds to standard MMST mapping and γ = 0 to Ehrenfest dynamics, as we will see below. Sometimes γ is treated as a parameter that is tuned according to some criterion, for example so that average adiabatic populations stay positive. 33 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Ẋ From Eq. (7c), we see that each potential energy surface, V n , contributes to the nuclear force proportionally to 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,cl + J 2,cl = 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. 31,32 It is possible to eliminate γ from the equations of motion by choosing a traceless form ofV , so that γ 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 andV (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, 7,8 in which one replaces the electronic operators by their expectation values:Â → ψ|Â|ψ .
Let us write an arbitrary electronic state |ψ in terms of diabatic basis states as where c n = 1 √ 2 (X n + iP n ) are the complex coefficients. Note that if |ψ 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 is equivalent to the MMST Hamiltonian in Eq. (6) with γ = 0. 2 The electronic equations of motion (7a) are equivalent to the solution of the time-dependent Schrödinger equation for the subsysteṁ which is also true for any value of γ. Because γ = 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 , that is |c 1 | 2 = 1 2 (X 2 1 + P 2 1 ) = 1 and |c 2 | 2 = 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, 9 which leads to the MMST Hamiltonian with γ = 1 in Eq. (5). This is not the only possible solution though; any γ = 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-1 2 system.
B. 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:Ĥ

z}) and the Pauli spin matrices arê
Such a decomposition was used already in the derivation of the spin-matrix method by Meyer and Miller. 46 The explicit relations between the quantities in Eqs. (1) and (13) are Note that in this way, it does not matter how the initial split between U andV is chosen. For simplicity we again assume that ∆ is real, so that H y = 0. The Hamiltonian in equation (13) is identical to that of a spin-1 2 particle in a magnetic field −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-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 as 14 where u denotes a unit vector defined by the spherical coordinates θ, ϕ and the states are normalized such that u|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Ŝ i in this state are where Thus S i (u) are the Cartesian coordinates of a sphere with radius r = 1/2, see Figure 1. One can think of them as The expectation values of the spin operators form the Cartesian coordinates of a sphere with radius 1/2. Their equations of motion describe precession of the spin vector around the external magnetic field, which is in one-to-one correspondence with the two-level potential-energy matrix. The "north" and "south" poles indicate points where all population is concentrated in one of the basis states of the subsystem.
orthogonal functions analogous to p-orbitals. Likewise the populations are and it is clear that J 1 (u) + J 2 (u) = 1 is always fulfilled. Note that on the "poles" of the sphere along the S z -axis (θ = 0, π), 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 Replacing all operators with their expectation values S = (S x (u), S y (u), S z (u)) gives 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). 47 The dynamics preserves the magnitude |S|, i.e. the radius of the sphere. It is the x-component of the field, H x = 2∆, 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. 46. 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.
C. 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Â associate this procedure with the mappinĝ which in literature is known as the Q-function. 43 The Qfunction allows us evaluate a quantum-mechanical trace over the subsystem, tr[Â], as an integral over a classical phase space: where the integration measure is du = 1 2π dϕdθ sin θ (which is normalized to be consistent with tr[Î] = 2). However, in order to calculate correlation functions, we need a phase-space representation for products of operators of the form tr [ÂB]. The Q-function is not enough for this purpose since because of the uncertainty property u|ÂB|u = u|Â|u u|B|u . The only case when equality holds is whenÂ orB is (proportional to) the identity operator 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. 48 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: Next, introduce the P-function 49 which together with the Q-function allows us to represent tr[ÂB] as Interestingly, one can also define a self-dual W-function which fulfils We will refer to the formulas (27) and (29) as the traciality property of the mapping. It is easy to show (using the special caseB =Î) 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 phasespace 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. 41,44 It has been applied to various problems in quantum optics 43 but to our knowledge not yet in the context of nonadiabatic molecular dynamics. To summarize, let us write the Stratonovich-Weyl kernelŝ w s (u) (s ∈ {Q, P, W}) collectively aŝ where we introduced the spin radius r s , defined in Table I for each s.
As simple examples, we have and all other operators can be built as linear combinations of these. In other words, the Q-function of the quantum spin vectorŜ defines the same sphere with radius r Q = 1/2 that we saw in Figure 1, whereas the corresponding W-and P-functions define larger spheres of radii r W = √ 3/2 and r P = 3/2, respectively. In quantum mechanics, the squared magnitude of the spin isŜ Its eigenvalues are equal to r 2 W , i.e. r 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 Q r 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 The spin operators in the W-representation form the Cartesian components of a sphere with radius √ 3/2. The system is entirely in |1 at the northern arctic circle defined by cos θc = 1/ √ 3, instead of at the north pole as in Figure 1. For 0 < θ < θc we have J1 > 1 and J2 < 0, and vice versa for π − θc < θ < π. The ZPE parameter γ = √ 3−1 measures the size of the regions with negative populations. The corresponding P-sphere looks similar but with a larger radius 3/2 and γ = 2, meaning that the red circles lie closer to the equator (cos θc = 1/3).
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,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,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 = p 2 2m + U + V 1 J 1,s + V 2 J 2,s + 2r s ∆ sin θ cos ϕ.
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. 22 The resulting equations of motion arė which can be shown to conserve H s (u). Importantly the subsystem equations are exact and equal for all s, and preserve the magnitude |u|. Thus the only dependence on s is in the nuclear equations of motion. Firstly, the contribution from the electronic coupling, ∆, scales by r s . Secondly, the region of phase space corresponding to inverted potential have different size. In this aspect, H 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.

D. Comparison with previous methods
Let us now connect spin mapping to the previous methods introduced in section II A. We introduce the coordinate transformation 2r s u x = X 1 X 2 + P 1 P 2 ≡σ x (X, P) (36a) 2r s u y = X 1 P 2 − X 2 P 1 ≡σ y (X, P) (36b) whereσ = (σ x ,σ y ,σ z ) is the MMST representation of the Pauli matrices. Using |u| 2 = 1 one can show that i.e. the new variables are also constrained to a sphere (a 3-dimensional hypersphere embedded in R 4 ) but with a different radius R = √ 4r s . The transformation used together with Eq. (37) brings the Hamiltonian in Eqs. (34) to the form V n 1 2 (X 2 n + P 2 n − 2r s + 1) + ∆(X 1 X 2 + P 1 P 2 ). (38) 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 γ = 2r s − 1. In particular H W has γ = √ 3 − 1 ≈ 0.732, which is the value Cotton and Miller gave a heuristic argument for in Ref. 34, also based on an analogy with spin (note that they define γ 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, 33 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. 51 Note however that in both the last two references, the mapping Hamiltonian is defined as nm H nm 1 2 (X n X m + P n P m − γδ nm ) where H nm = ( p 2 2m + U )δ nm + V nm . This is different from our approach, which treats terms proportional to the identity differently than other operators.
The split between U andV is therefore unambiguous in our approach. This is clear from Eq. (34) that only multiplies the traceless part ofV by r s . In this way, whenV goes toÎ, the mapped Hamiltonian smoothly reduces to the single-state Hamiltonian. In LSC-IVR however, different choices of the split between U andV generally lead to different results since the total population J 1,cl + J 2,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, 52 our method maps two levels to one spin-1 2 particle, whereas Cotton and Miller map to two spin-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, 16 who related their work to the spin-coherent state propagator by Suzuki. 53 They expressed a semiclassical propagator in the MMST representation as an integral over the coordinates of the spin coherent state 54 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.

E. 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 where Tr denotes a full quantum-mechanical trace,ρ is a density matrix normalized to give Tr[ρ] = 1, and B(t) = e iĤtB e −iĤt . In this paper we will only consider the case whenÂ,B are electronic operators and the density matrix is state-independent,ρ =ρ nuc ⊗ρ el , wherê ρ el = 1 2Î (the normalization ensures that tr[ρ el ] = 1). The initial projection onto one of the electronic states is defined byÂ.
We will use a Wigner distribution of the F nuclear degrees of freedom, and a quasiclassical propagator analogous to classical Wigner dynamics: C AB (t) ≈ dx dp du 2 ρ nuc (x, p)A s (u)Bs(u(t)), (41) where the dual indicess are specified in Table I. The trajectory u(t) is generated by H s , with the same index as the operator at time zero. (We also considered the opposite index choice AsB s but found that it gives less accurate results at t > 0.) We integrate over du 2 = dϕ 4π dθ sin θ using a Monte Carlo scheme, drawing samples uniformly from the unit sphere |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 = √ 4r s = 2(γ + 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: 24 the first, LSC-IVR, has projection operators both at initial and final times, 20,21 and the second, PBME, has a projection only at initial times. 23 As a general formula we write C AB (t) ≈ dx dp dX dP ρ nuc (x, p)ρ el (X, P) × A(X, P)B(X(t), P(t)), (42) where X = (X 1 , X 2 ), P = (P 1 , P 2 ). The explicit expressions for ρ el (X, P) are given in Table II together with the appropriate observables in the various methods. For a general operatorÂ = A 0Î + 1 2 A ·σ, the MMST mapping representations used in the standard formulations of LSC-IVR and PBME are defined as where A SEO is used whenever one projects the harmonic oscillators to the physical subspace, and φ = 16 exp(−R 2 ). The special case of a population operator, A = |n n| = (Î +σ z )/2 gives Similarly we can also write the Stratonovich-Weyl observables in terms of {X, P} using Eq. (36), which gives Note that A SW gives a standard expression for a population operator, A SW = 1 2 (X 2 n + P 2 n − γ) = J n,cl , whereas B SW has not previously been used in quasiclassical mapping dynamics.  (42) to compute general correlation functions CAB(t) for various methods. We use the short-hand notation R 2 = X 2 1 + P 2 1 + X 2 2 + P 2 2 , φ = 16 exp(−R 2 ) and N = [2π 2 (γ + 1)] −1 . Method TABLE III. Definition of the terms used in Eq. (42) to compute correlation functions CAB(t), where the initial state iŝ A = |n n|, for various methods employing focused initial conditions. The electronic distribution in each case is given by ρ el (X, P) = N δ(R 2 − 2(γ + 1)) and A = δ(X 2 n + P 2 n − (γ + 2)), 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Â = |1 1|, an Ehrenfest calculation would use the initial conditions |c 1 | 2 = 1 and |c 2 | 2 = 0, which in the MMST representation is ρ el A = N δ(R 2 − 2(γ + 1)) δ(X 2 1 + P 2 1 − (γ + 2)), (47) where N = π −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 u from the northern polar circle in Figure 2, which corresponds to Eq. (47) with γ = √ 3 − 1 ≈ 0.732. A focused Q-distribution is identical to Ehrenfest dynamics. Unlike the approaches in Table II, 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. 33 and 55.
Note that when using focused initial conditions, the Stratonovich-Weyl operators reduce to the usual Ehrenfest (MFT) and MMST operators. For γ = 0 and similarly for γ = 1 one has 1 = (R 2 − 2)/2, so that B SW = B wig .

III. RESULTS AND DISCUSSION
To test the accuracy of our method, we have applied it on a spin-boson model, 56 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 Here ε is the energy bias and ∆ the constant diabatic coupling between the two electronic states. The bath consists of F nuclear modes, each with frequency ω j and vibronic coupling coefficient c j . In this form the matrix V is already traceless. We choose units such that = 1. The so-called spectral density of the bath, J bath (ω), determines the distribution of nuclear frequencies. Among the most common forms of this function is the Ohmic bath where ω c is called the characteristic frequency and ξ the Kondo parameter, which determines the strength of the friction. Numerically one uses a discretization with a finite F . In this work we used the discretization scheme employed in Ref. 57.
To sample the initial nuclear variables {x j , p j }, we used the thermal Wigner distribution with α j = tanh 1 2 βω j and β = 1/k B T . To facilitate comparison with other recent papers, we have used the same parameter settings as in Ref. 35, which are given in Table IV. 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. 35 (except for subtle changes in the strong-bath case). In our calculations we used F = 100 nuclear modes and a timestep of 0.01∆ −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 II are compared in Figure 3 with numerically exact path-integral calculations (QUAPI). 58 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 trick 48 to treat the identity operator in LSC-IVR or PBME gives essentially the same results as our Qmethod. There is however no direct connection between the two, since Ref. 48 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 γ 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. 30,35,59 When computing correlation functions of off-diagonal operators like C σxσ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Â.
For the methods with focused initial conditions in Table III and Figure 5, the focused W-method is clearly superior to both Ehrenfest and standard MMST (γ = 1). It uses γ = √ 3 − 1 ≈ 0.732, which is close to the values that Müller and Stock found optimal in their simulations using action-angle (focused) initial conditions. 33 Both Ehrenfest (γ = 0, which is the same as a focused Q-method) and standard MMST (γ = 1) with focused initial conditions are much less accurate for asymmetric problems. We also tried a focused P-method (γ = 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 γ leads to a lowering of the curve C 1σ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 γ in their method involving square windows. 34,35 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. 52 (which maps to two spins instead of one) has been found to give worse results than SQC defined in the MMSTbasis, 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. 30, 59 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. 20 It will also suffer from the usual classical-Wigner problem of ZPEleakage 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. 60 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  Table II and the model parameters are given in Table IV. 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-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 selfdual 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 twolevel 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   Figure 3, but comparing the focused methods defined in Table III. couple of important differences. Firstly, the total population is always fixed to one, because the {X, P} variables are constrained to a three-dimensional sphere. Secondly, it uses a different value of the ZPE-parameter, γ, 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 (γ = 1), since one would still need to rescale the traceless part of theB 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), 61-64 which would add nuclear quantum effects to the simulations, as well as to allow sampling from an exact thermal distribution.
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 (u, x, p) as 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.