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

We present a new partially linearized mapping-based approach for approximating real-time quantum correlation functions in condensed-phase nonadiabatic systems, called spin-PLDM. Within a classical trajectory picture, partially linearized methods treat the electronic dynamics along forward and backward paths separately by explicitly evolving two sets of mapping variables. Unlike other previously derived partially linearized methods based on the Meyer-Miller-Stock-Thoss mapping, spin-PLDM uses the Stratonovich-Weyl transform to describe the electronic dynamics for each path within the spin-mapping space; this automatically restricts the Cartesian mapping variables to lie on a hypersphere and means that the classical equations of motion can no longer propagate the mapping variables out of the physical subspace. The presence of a rigorously derived zero-point energy parameter also distinguishes spin-PLDM from other previously proposed partially linearized approaches. These new features appear to give the method superior accuracy for computing dynamical observables of interest, when compared with other methods within the same class.


I. INTRODUCTION
The coupled dynamics of electrons and nuclei in molecular condensed-phase systems remains a challenging problem for computer simulation. [1][2][3] Nonadiabatic transitions can be induced when the electronic adiabatic states of a system become close in energy at some nuclear geometry, which results in a breakdown of the Born-Oppenheimer approximation.
One approach is to use numerically exact wavefunction-based methods to compute the dynamical observables of interest. The advantage of such methods is that the obtained results are systematically improvable on increasing the size of the basis; they can also in principle describe nuclear quantum effects, such as tunnelling and zero-point energy. Of particular promise are methods, like the time-evolving block decimation (TEBD) technique, 4 which can even be applied to fairly large systems. 5 However such methods often have restrictions on the type of systems that they can treat. For example, TEBD can only be efficiently applied to spatially linear systems with at most nearest-neighbour couplings.
Due to the continuous nature of the nuclear subspace, the dynamics of such degrees of freedom are ideally suited to be performed using classical trajectories. The advantage of using such methods is that they are numerically cheap and can be easily applied to large, condensedphase problems. Of these methods, the most popular are Tully's fewest-switches surface hopping, 6 which is heuristically motivated and hence cannot be rigorously derived and Ehrenfest dynamics, 7,8 which neglects all quantum entanglement between the electronic and nuclear subsystems. Recently, there has been a renaissance in mappingbased techniques, [10][11][12][13][14][15][16][17][18][19][20][21][22][23] which describe the dynamics of the electronic subsystem as well as the nuclei using classical trajectories within a continuous mapping space. Such methods have been shown to have superior accuracy compared to Ehrenfest dynamics, when calculating dynamical observables for a wide range of model systems. While these methods still treat the nuclear degrees of freedom classically and hence neglect nuclear quantum effects, these could be perhaps reintroduced by using ringpolymer based formulations. [24][25][26][27] Although none of these approaches will be able to describe quantum nuclear coherence effects, 28 such effects are typically unimportant in condensed-phase systems. 29 Historically, the vast majority of mapping-based approaches have used the so called Meyer-Miller-Stock-Thoss (MMST) mapping. 30,31 In this scheme, the electronic subsystem is mapped onto a set of singly-excited harmonic oscillators. The nuclear and electronic dynamics are then described as an average over many classical trajectories. While trajectory-based methods using this mapping have been able to qualitatively reproduce the correct dynamics in a range of model systems, such methods may also exhibit significant problems. One of the most important of these is zero-point energy leakage, where zero-point energy can unphysically flow between the mapping harmonic oscillators as a result of the trajectories leaving the physical subspace during the classical dynamics. 32,33 In other words, the trajectories can evolve into areas of phase space in the mapped system which do not correspond to a valid state in the electronic system. To alleviate this problem, related approaches have been suggested, such as reducing the zero-point energy parameter in the underlying theory, 34 symmetric windowing 13,35 of trajectories, an 'identity trick' [21][22][23] and using alternative classical mapping models. 16,36 Recently a new form of mapping, called spinmapping, 37,38 has been suggested. While this approach leads to exactly the same equations of motion for the tra-arXiv:2007.05047v1 [physics.chem-ph] 9 Jul 2020 jectories as MMST mapping, the Cartesian mapping variables are constrained to a hypersphere which is isomorphic with the phase space of the actual electronic subsystem. This means that unlike standard MMST mapping, spin-mapping does not need projection operators onto the physical subspace and trajectories no longer suffer from zero-point energy leakage. Also, when converted to Cartesian mapping variables, the spin-mapping Hamiltonian also has a different zero-point energy parameter to that of MMST mapping.
When derived from a path-integral formalism, previous mapping-based techniques generally fall into one of two categories. Fully linearized methods 23,29,39-41 result from performing a linearization approximation to the difference between the forward and backward paths for both the electronic and nuclear degrees a freedom; a semiclassical approximation that is expected to be valid in the classical limit. In contrast, partially linearized methods result from only performing a linearization approximation for the nuclear paths and then treating the dynamics of the forward and backward electronic paths separately. Examples of partially linearized methods using MMST mapping are the partially linearized density matrix (PLDM) approach 11,12,[42][43][44][45][46][47][48][49][50][51] and the forwardbackward trajectory solution (FBTS). [17][18][19] While spinmapping has already successfully improved the accuracy of so called fully linearized methods, 37,38 it has yet to be applied to partially linearized methods. Often partially linearized methods are better than their fully linearized counterparts, so we may optimistically expect a spin-mapping version of PLDM to be the best method of all.
In this paper, after a review of the standard PLDM based on MMST mapping, we use spin-mapping to derive a new partially linearized approach, which we call spin-PLDM. The method is derived generally for systems containing any number of electronic states by employing the framework introduced in Ref. 38. We show that such a method is able to more accurately reproduce the correlation functions for the spin-boson model than standard PLDM. To aid the reader, a schematic illustrating the relationship of spin-PLDM to other mapping-based classical-trajectory techniques is given in Fig. 1. In Paper II, 52 the spin-PLDM method is tested further and analyzed extensively. In addition, in this following paper we will compare spin-PLDM with a fully linearized spin method and show that spin-PLDM has greater accuracy than other mapping-based trajectory methods.

II. THEORY
In this paper, we study nonadiabatic dynamics using the specific example of the familiar electronic-nuclear coupled system. Naturally, the method we derive in this paper could also be applied to any other quantumclassical system. The Hamiltonian for such a system can always be written in the following form: where V 0 (x) is a purely nuclear potential andV (x) has an F × F matrix representation in the basis of electronic states, |1 , . . . , |F ; these are uniquely partitioned so that tr[V (x)] = 0. We define F as the dimensionality of the electronic subsystem and tr[·] as the partial trace over the electronic degrees of freedom. Both x and p are vectors of dimension f , which have been mass weighted such that all degrees of freedom have the same mass m. We set = 1 throughout. Many dynamical quantities of interest can be written in the form of a real-time quantum correlation function: 53 whereÂ is an electronic operator effectively initializing the system at t = 0, andρ b is the initial density matrix for the nuclear bath. These are known as factorized initial conditions. In addition,B is an electronic operator that measures the system at time t. Although we only consider the case whereÂ andB are purely electronic operators in this paper, the theory could be easily extended for more general real-time correlation functions. Inserting complete sets of nuclear position eigenstates (|x 0 , |x 0 and |x N ) and electronic basis states (|λ and |µ ) results in the following expression for the real-time 2. An illustration of the discretization of the forward (blue) and backward (red) paths used to obtain a pathintegral representation of the CAB(t) correlation function.
The variables x k and p k correspond to the position and momentum basis states inserted at each time-step of the realtime propagation respectively. correlation function: One of the main challenges in evaluating this expression for the real-time quantum correlation functions is obtaining an efficient expression for the forward and backward real-time propagators. By inserting a complete set of position and momentum eigenstates at each time-step, matrix elements of the real-time propagators can be written in a path-integral formalism: 11 where N is the number of time-steps. The discretization scheme used to obtained this expression is illustrated by the blue forward path in Fig. 2. The nuclear action is given by: with = t/N . This expression was obtained using the following splitting of the real-time propagator for a single time-step: and inserting complete sets of position and momentum states. Hence, Eq. (4) is exact in the N → ∞ limit, when → 0. Additionally, T [µ,λ] is the electronic transition amplitude, which as a function of the nuclear path, {x 1 , . . . , x N }, is given by: Similar expressions for the backwards path can be obtained by using the backwards propagator, e iĤt , in Eq. (4) and following the same steps outlined above. The equivalent expressions are marked with a prime. When derived from a path-integral representation, all mapping-based classical-trajectory approaches first implement a linearization approximation 29 for the nuclear degrees of freedom to avoid sampling complex phases. The linearization approach is also well-known for single-surface systems where it leads to classical Wigner dynamics. 29,54,55 Hence such a linearization approximation amounts to treating the nuclear degrees of freedom classically; such techniques will therefore not be able to describe nuclear quantum effects, such as delocalization, interference and tunnelling through potential energy barriers. First the nuclear degrees of freedom are transformed to sum and difference coordinates: where (x k , p k ) and (x k , p k ) are the nuclear phase-space variables at time-step k for the forward and backward paths respectively. In addition, we note that x N = x N as shown in Fig. 2, because there is no nuclear operator applied at time t within the correlation function C AB (t).
Using these new coordinates, the real-time correlation function is given exactly in the N → ∞ limit by: where the initial nuclear density matrix is now described by its Wigner transform: Defined in this way, the Wigner transform of the initial nuclear density is normalized such that dx 0 dp 0 ρ W b (x 0 ,p 0 ) = 1. The linearization approximation involves only retaining terms in the nuclear and electronic action up to first order in the difference coordinates, ∆x k . Transforming variables and using a Taylor expansion of the potential V 0 up to first order in ∆x k (i.e., V 0 (x k ) ≈ V 0 (x k ) + 1 2 ∇V 0 (x k )∆x k ), the difference in the nuclear action between the two paths can be approx-imated as: 11,54 where is the state-independent nuclear force and ∇ is the gradient, a vector of derivatives with respect to nuclear positions.
The electronic action is included within the definition of the electronic transition amplitude, given by Eq. (7). Transforming variables in a similar way and using a Taylor expansion of the potentialV (x k ) up to first order in ∆x k (i.e., e −iV (x) ≈ e −i∇V (x) ∆x/2 e −iV (x) ) the electronic transition amplitude for the forward and backward paths can be approximated in an analogous way to the nuclear action: where there is no term which depends on ∆x N , because the boundary conditions of the forward and backward nuclear paths satisfy x N = x N , as illustrated in Fig. 2. This quantity contains two types of propagator: the electronic propagator, e ±iV (x k ) , which describes the electronic dynamics; and the coupling propagator, e −i∇V (x k ) ∆x k /2 , which due to its dependence on ∆x k adds an electronic contribution to the nuclear force, as we will show. With these approximate expressions for the nuclear action between the two paths [Eq. (11)] and the transition amplitudes for the forward and backward paths [Eq. (13)], the real-time correlation function would be an exact solution of the quantum-classical Liouville equation (QCLE). [56][57][58] However, it is known that it is not possible to implement such an expression in an independent trajectory simulation. 56 The easiest way to further simplify the expression for the real-time correlation function is by introducing a mapping variable representation for the electronic transition amplitudes. In general this cannot be performed exactly within a classical-trajectory picture and hence different mapping-based classical-trajectory techniques differ by how they approximate the electronic transition amplitudes. We first consider the electronic transition amplitudes within the MMST mapping representation and then how these quantities are approximated within the standard partially linearized density matrix (PLDM) approach. As we now have only a single set of nuclear coordinates, we setx → x andp → p throughout the rest of the paper.
A. The standard partially linearized density matrix approach The standard partially linearized density matrix approach (PLDM) 11 represents the electronic transition amplitudes using the Meyer-Miller-Stock-Thoss mapping. 30,31 Within this MMST mapping, the electronic basis states, |λ , are described by the single phonon excitation subspace of a set of F harmonic oscillators, where F is the size of the electronic system: 31 In the case of operators, the mapping can be accomplished by replacing an electronic annihilation operator, a λ , with the harmonic oscillator annihilation operator for oscillator number λ: Harmonic oscillator coherent states, 59 |Z = |Z 1 , Z 2 , · · · , Z F , form a convenient basis with which to evaluate the electronic transition amplitudes given by Eq. (13). The harmonic oscillator coherent states are defined as eigenstates of the harmonic oscillator annihilation operator: where Z λ = X λ + iP λ is a complex number, which encodes the average nuclear position and momentum of degree of freedom λ. They are normalized such that First, these harmonic oscillator coherent states can be used to represent the coupling propagator, which appears in the definition of the electronic transition amplitudes: 60 where the (scalar) mapping representation of the traceless potential operator is The result given by Eq. (17) can be confirmed by performing a Taylor expansion of the exponential on both sides of the equation. As the error in this expression is of order 2 , the substitution is valid when the N → ∞ limit is taken. Second, inserting Eq. (17) into Eq. (13a) at each timestep along the path results in a harmonic oscillator coherent state appearing to the right of every electronic propagator (except for the e −iV (x1) operator, which we post-multiply byÎ = (2π) −F dZ 0 |Z 0 Z 0 |). The effect of these electronic propagators can be accounted for by simply evolving the harmonic oscillator mapping variables at time-step k along the path, Z k , as follows: 60 where F ( )} are the mapping variables evolved in time by according to: Eq. (20) is an exact solution of Eq. (19), for any size timestep . This is necessary condition for a mapping-based classical-trajectory technique to be exact for a purely electronic system. The coupling propagator is not treated in the same way as the electronic propagator in PLDM, because we require the ∆x k dependence of this quantity to be purely a complex exponential, so that it can be easily incorporated with the nuclear action in Eq. (11) (17) and (19) lead to an expression for the electronic transition amplitudes in terms of the harmonic oscillator mapping variables. However, the expression contains an integral over mapping variables at each time-step, k. To generate a deterministic trajectory picture for the electronic dynamics, the overlap of harmonic oscillator coherent states at different time-steps are approximated as: 60 which is not an equality as coherent states are actually overcomplete. 59 Additionally, the Dirac delta function of a complex argument, Z, is defined as: ). This approximation in Eq. (21) reduces the expression for the electronic transition amplitudes to: 11 where we now set Z 0 = Z and Z 0 = Z . The electronic action for the forward and backward paths are defined as: where t k = k is the time at time-step k. Additionally, |Z(t) is the time-evolved harmonic oscillator coherent state, which is defined as: whereÛ (t) is is a time-ordered propagator for the electronic states according to the time-dependent potential obtained by moving along a given path x(t k ) = x k : To complete the linearization approximation of the nuclear path, outlined in Section II, the integrals over ∆x k and ∆p k must be performed. Using the approximate expressions for the nuclear and electronic action, given by Eqs. (11) and (23), these integrals can be performed analytically to give: 11 The result is a product of Dirac delta functions which constrain the dynamics of the nuclear variables to follow Newton's equation of motion. Here, is the electronic dependent nuclear force, which depends on the electronic mapping variables of the two paths. This means that the approximate expression for the real-time correlation function becomes: where we now set x 0 = x and p 0 = p. Writing theÂ andB operators in terms of the electronic basis states |λ means that their associated matrix elements in the harmonic oscillator mapping space can then be evaluated using Eq. (14) and the well-known harmonic oscillator wavefunctions to give: 11,60 C AB (t) ≈ dx dp dZ dZ In this expression, the harmonic oscillator mapping variables Z and Z are initially sampled from the following Gaussian distribution: and the electronic operators within the real-time quantum correlation function are now described in terms of these mapping variables as: The dynamics of the mapping coherent states, given by Eq. (24), and the dynamics of the nuclear phase-space variables, defined by the Dirac delta functions in Eq. (26), can be described by evolving the classical nuclear and electronic degrees of freedom under the following equations of motion: 11,60 The discretized PLDM equations of motion are illustrated schematically in Fig. 3, for time-steps of size . The electronic mapping variables for the forward and backward paths (Z and Z respectively) are not directly coupled, but instead couple via the nuclear degrees of freedom, x, through the mapping Hamiltonian, V m . This coupling within the equations of motion means that the back action of the nuclear environment on the electronic subsystem is included within this technique. The equations of motion for the standard PLDM approach previously derived by Huo and Coker 11 differ slightly from the expression in Eq. (32). This is because within the standard PLDM approach, the potential is not separated into a purely nuclear potential, V 0 (x), and an electronic-dependent potential,V (x). This results in the whole nuclear force explicitly depending on the mapping variables Z and Z . The forward-backward trajectory solution (FBTS), 17,18 which does separate the nuclear potential into V 0 (x) andV (x) components, has identical equations of motion as Eq. (32) and also has the same expression for the real-time correlation function as Eq. (29). FBTS is however derived as an approximate solution to the QCLE 56 and not from a path-integral representation of the real-time correlation function. Hence the derivation presented in this section uses elements from the derivation of both the standard PLDM and FBTS approaches. Within the spin-mapping approach, the separation of these two nuclear potential components arises naturally from the underlying theory and hence it is no longer an arbitrary choice that one must make during the derivation of mapping-based methods.
As the QCLE is derived by taking a partial Wigner transform of the exact time-dependent Schrödinger equation, it describes the exact dynamics of electron-nuclear systems in the classical nuclear limit. The derivation of FBTS involves applying the same approximation as standard PLDM, given by Eq. (21), to an otherwise exact mapping variable solution of the QCLE. This means that in principle a more accurate solution of the QCLE dynamics could be found by improving upon the approximation given by Eq. (21). To do so, we derive PLDM within the spin-mapping space.

B. Spin coherent states
The main idea behind spin-mapping is the representation of the electronic states of a system by a set of spins, described by classical vectors with fixed radii. This can be achieved in several ways. The most analogous form of spin-mapping to that of MMST mapping is to represent a F electronic state system by F spin-1 2 particles. 61 Like MMST mapping, this spin-mapping is mathematically given by Eq. (14), except each degree of freedom can now only contain 0 or 1 excitation, corresponding to the two quantum states of each spin- 1 2 particle. While the size of the unphysical subspace available to spin-mapping is clearly smaller than that of the MMST mapping space, it is still larger than necessary, and more importantly it also has the unfavourable feature that the classical dynamics are not exact for the pure electronic subsystem. 36 This probably explains why this particular form of spinmapping has been found to give rise to less accurate dynamical observables compared to MMST mapping for many systems. 61 A seemingly more successful version of spin-mapping involves parameterizing the probability amplitudes associated with the electronic basis in terms of spin angles. Now the associated classical dynamics, which take the form of spins precessing around an effective magnetic field, exactly describe the pure electronic subsystem dynamics and the spin-mapping space now contains no unphysical regions. While this form of spinmapping has been around since the advent of mappingbased techniques, 62 the underlying equations have recently been reformulated using a Stratonovich-Weyl approach, which leads to a more rigorous way of determining the optimum spin radius (or equivalently the zeropoint energy parameter) 37 as well as a generalization to systems with more than two levels. 38 This reformulation has already been shown to produce more accurate observables compared to linearized MMST techniques. Hence, it is this approach that we will use to derive a spin-mapping version of PLDM.
In deriving PLDM within this spin-mapping space, we first restrict ourselves to systems containing two electronic states and generalize to multi-state problems in Sec. III. The most natural generalization of standard PLDM to the spin-mapping space is found by using spin coherent states. This can be easily achieved by performing the same steps as in Sec. II A, but using spin coherent states instead of harmonic-oscillator coherent states. Although it turns out that first this leads to a method which is not particularly accurate, the accuracy can be significantly improved by generalising the spin sphere using Stratonovich-Weyl kernels, as was similarly found for fully linearized spin-mapping methods. 37,38 This extra step will be performed in Sec. II C.
Within the two-level electronic state space, spin coherent states with similar properties to the harmonicoscillator coherent states can be defined as follows: 63 where 0 ≤ ϕ < 2π and 0 ≤ θ ≤ π are the angles of the spherical polar representation of a three-dimensional normalized spin vector, u, with components: u x = u|σ x |u = sin θ cos ϕ, (34a) u y = u|σ y |u = sin θ sin ϕ, (34b) u z = u|σ z |u = cos θ. (34c) In this expression,σ j correspond to the 2 × 2 Pauli spin matrices. This spin vector, u, along with its corresponding polar angles are given in Fig. 4. A minor problem with Eq. (33) is that this particular ansatz does not include an overall phase for the electronic state. It turns out that this phase information will be required for our spin-mapping version of PLDM, so as to be able to correctly describe the electronic transition amplitudes given by Eq. (13). One way of including an overall phase in the electronic state is by the alternative representation of the coherent state in terms of the probability amplitudes, c λ : where Ω = {c 1 , c 2 } encodes the complex values of these probability amplitudes, which are normalized such that |c 1 | 2 + |c 2 | 2 = 1. It was found that when generalising spin-mapping to F -level electronic systems, the underlying equations were also much simpler when the coherent states were parameterized in terms of their probability amplitudes, c λ . 38 Written in this way, the integral over the spin sphere can equally well be written as an integral over Ω. In other words, 1 2π 2π 0 dϕ π 0 dθ sin θ · · · = dΩ · · · , where the integration element is defined as: and dc λ = dRe[c λ ]dIm[c λ ]. In this expression, the factor of 2 appears such that the spin coherent state integrals satisfy dΩ = F . Physically, the Dirac delta function in Eq. (36) guarantees that the coherent state |Ω is correctly normalized. First these spin coherent states can be used to represent the coupling propagator, which appears in the definition of the electronic transition amplitudes: where the mapping representation (scalar) of the traceless potential operator (matrix) is now: Eq. (37), as for the harmonic oscillator coherent states, is exact up to first order in , which can again be confirmed by comparing the Taylor expansions of both sides and performing the integrals. The required integrals are easily performed by transforming to polar coordinates, with c 1 = cos θ 2 e −iϕ/2 and c 2 = sin θ 2 e iϕ/2 , as defined in Eq. (33). This transformation is valid here, because the integrand on the right hand side of Eq. (37) does not depend on the overall phase of the probability amplitudes, c λ . As the error in Eq. (37) is of order 2 , it is hence valid when the N → ∞ limit is taken. It will become clear why we have labelled the mapping representation of the traceless potential operator with the 'P' subscript in Sec. II C when we will also give an interpretation of the factor of 3 that appears in Eq. (38).
Second, like in the derivation of standard PLDM, inserting Eq. (37) into Eq. (13a) results in a spin coherent state appearing to the right of every electronic propagator (except for the e −iV (x1) operator, which we postmultiply byÎ = dΩ 0 |Ω 0 Ω 0 |). The effect of these evolution operators can be accounted for by simply evolving the mapping variables, Ω, as follows: where Ω k ( ) is defined by probability amplitudes, c λ ( ), at time , evolved according to Eq. (20) with c λ appearing instead of Z λ . The time-evolved electronic state, |Ω k ( ) , now contains an overall phase due to the time propagation and hence cannot be described by the spin vector, u, introduced in Eq. (33). Evolving the spin-mapping variables, Ω k ( ), in this way exactly solves Eq. (39), which means that PLDM in this spin-mapping space will also be exact for a purely electronic system. For the backward electronic transition amplitude, T [λ ,µ ] , inserting Eq. (37) into Eq. (13b) results in a spin coherent state appearing to the left of every electronic propagator, e iV (x k ) . Hence the complex transpose of Eq. (39) is used in this case. Therefore, Eqs. (37) and (39) lead to an expression for the electronic transition amplitudes in terms of the Ω mapping variables. As before, the expression contains an integral over mapping variables at each time-step, k. To generate a trajectory picture for the electronic dynamics, the overlap of spin coherent states at different time-steps are approximated as: which, like for harmonic oscillator coherent states, is not an equality as coherent states are also overcomplete. 63 As before, the Dirac delta function of a complex argument, Ω, is defined as: δ(Ω) = δ(c 1 )δ(c 2 ) and δ(c λ ) = δ(Re[c λ ])δ(Im[c λ ]). This approximation in Eq. (40) reduces the expression for the electronic transition amplitudes to: where we now set Ω 0 = Ω and Ω 0 = Ω . The electronic action for the forward and backward paths are defined as: Additionally, |Ω(t) is the time-evolved spin coherent state, which is defined by Eq. (24) with Ω instead of Z. These transition amplitudes are equivalent to the initial-value representation of the spin coherent propagator, which has been derived previously. [64][65][66] To complete the linearization approximation of the nuclear path, outlined in Sec. II, the integrals over ∆x k and ∆p k must be performed. Using the approximate expressions for the nuclear and electronic action, given by Eqs. (11) and (42), these integrals can be performed analytically using Eq. (26). The result is a product of Dirac delta functions which constrain the dynamics of the nuclear variables to follow Newton's equation of motion.
, Ω (t k ), x k ) is the nuclear force at time-step k and is the electronic dependent nuclear force, which depends on the electronic mapping variables of the two paths. This means that the approximate expression for the real-time correlation function when using spin coherent states becomes: The evolution of the classical nuclear and electronic degrees of freedom when using spin coherent states at first sight looks different from that for standard PLDM, given by Eq. (32). However, the equations of motion become identical if the spin coherent state equations of motion are written in terms of Cartesian mapping variables, using the transformation: c λ → Z λ / √ 6. This can be seen by comparing the two mapping expressions for the traceless potential operator, given by Eqs. (18) and (38). The only difference between the two methods, therefore, is the integration space, with the initial mapping variables now being confined to a spin sphere when PLDM is derived using spin coherent states. The spin sphere on which the dynamics are performed corresponds to the 'P' sphere, which has previously been used for fully linearized spinmapping in Ref. 37. This was found to give inaccurate results within the fully linearized theory and our preliminary tests showed that the partially linearized form also performs badly. A similar correlation function to Eq. (44) was also derived in Ref. 67 (except without applying the linearization approximation to the nuclear paths) where they also found the inaccuracies to be quite large. However, the fully linearized version comes in three flavours corresponding to Q, W and P-spheres, of which the symmetric W-sphere was found to be the most reliable for a wide range of benchmark systems. 38 We hence now consider how to generalize this partially linearized method to other spin spheres.

C. The Stratonovich-Weyl transform
The Stratonovich-Weyl transform 68 can be thought of as a discrete version of the Wigner transform within the spin-mapping space. Within this formalism, operators are represented by their Stratonovich-Weyl functions, which for two-level systems are given as follows in terms of the three-dimensional spin vector, u: whereŵ s (u) is the Stratonovich-Weyl kernel, which is defined as:ŵ Here, R s defines the spin sphere radius. Examples of different possible spin spheres for two-level systems are given in Table I. Like for the Wigner transform, the trace of a product of two operators can be exactly written as an integral over the product of their corresponding Stratonovich-Weyl functions: tr[ÂB] = du A s (u)Bs(u). Notice that the two operators,Â and B, are represented in 'dual' spin spheres s ands respectively. The 'dual' of each spin sphere, s, is given in Table I for a two electronic state system.
As in the previous section, we will want to represent the Stratonovich-Weyl kernel in terms of the probability amplitudes of the coherent states, Ω = {c 1 , c 2 }, so that the overall phase of the electronic state can be correctly described within the electronic transition amplitudes. This can be achieved by first writing the Stratonovich-Weyl kernel in terms of the spin coherent state, given by Eq. (33) and then re-expressing this coherent state in terms of these probability amplitudes, as given by Eq. (35). This results in the following expression: where γ s = (R 2 s −2)/2 is the zero-point energy parameter, whose values for different spin spheres is given in Table   I for a two electronic state system. From Table I, we see that R 2 Q = 2 and γ Q = 0, so thatŵ Q (Ω) = |Ω Ω|. Hence the Stratonovich-Weyl kernel can be thought of a generalization of the coherent state outer product used in the previous section. The previous section effectively uses Q-type Stratonovich-Weyl kernels, which must be paired with its dual P-type Stratonovich-Weyl function, which hence results in P-sphere dynamics for the classical trajectories.
The Stratonovich-Weyl kernel can also be used to represent electronic operators; namely: where V s (Ω, x) is the Stratonovich-Weyl function of the traceless electronic HamiltonianV (x), which from Eqs. (45) and (47) can be shown to be Hence the coupling propagator, which appears in the definition of the electronic transition amplitudes, can also therefore be obtained in terms of the Stratonovich-Weyl kernel as: In general, this expression is again exact up to first order in , which can be confirmed by comparing a Taylor expansion of both sides of Eq. (50), while using the properties of the Stratonovich-Weyl kernel given by Eqs. (48a) and (48b). However, for the case s =s = W, the same Taylor expansion shows that the value of the O( 2 ) term is also correctly obtained for this spin sphere. This is perhaps an indication that there is an advantage in using the W-sphere for the PLDM method within the spinmapping space, even though Eq. (50) becomes exact for all spin spheres when the N → ∞ limit is taken. As for the previous derivations, inserting Eq. (50) into Eq. (13a) results in a Stratonovich-Weyl kernel appearing to the right of every electronic propagator (except for the e −iV (x1) operator, which we post-multiply bŷ I = dΩ 0ŵs (Ω 0 )). We then define the time-evolved Stratonovich-Weyl kernel,ŵ s (Ω k , ), as follows: Unlike for spin coherent states, the time-evolved Stratonovich-Weyl kernel cannot be obtained purely by evolving the spin coherent probability amplitudes, c λ . This is because the Stratonovich-Weyl kernel contains a zero-point energy parameter, γ s , which must also be evolved forward in time. Hence, the time-evolved Stratonovich-Weyl kernel is given by: where Ω k ( ) are defined by the spin coherent state probability amplitudes, c λ ( ), at time , evolved according to Eq. (20) with c λ instead of Z λ . Additionally, e −iV (x k+1 ) is a 2 × 2 matrix in the electronic-state basis, which can easily be computed numerically. For the backward (Ω ) electronic transition amplitude, inserting Eq. (50) into Eq. (13b) results in a Stratonovich-Weyl kernel appearing to the left of every electronic propagator, e iV (x k ) . Hence the complex transpose of Eq. (51) is used in this case. Therefore, Eqs. (50) and (51) lead to an expression for the electronic transition amplitudes in terms of the Ω mapping variables. However, the expression contains an integral over mapping variables at each time-step, k. Following the standard PLDM/FBTS procedure, the overlap of Stratonovich-Weyl kernels at different time-steps are approximated as follows, to generate a trajectory picture for the electronic dynamics: (53) This approximation in Eq. (53) is exact if ∆x k = 0 (i.e., in the absence of electron-nuclear coupling) and is an identical approximation to Eq. (40) whens = Q, such thatŵ Q (Ω k , ) = |Ω k ( ) Ω k |. Employing the approximation given by Eq. (53) reduces the expression for the electronic transition amplitudes to: where we now set Ω 0 = Ω and Ω 0 = Ω . The electronic action for the forward and backward paths are defined as: where t k = k is the time at time-step k. Additionally, ws(Ω, t) is the time-evolved Stratonovich-Weyl kernel, which is defined as: where |Ω(t) is the time-evolved spin coherent state, which is defined by Eq. (24) with Z → Ω andÛ (t) is the time-ordered propagator, given by Eq. (25). The time-ordered propagator is a matrix product involving 2 × 2 matrices, e −iV (x k ) , which is hence easily calculated for each trajectory. The back action does not appear directly in the time-ordered propagator, but is treated by the evolution of the Cartesian mapping variables.
To complete the linearization approximation of the nuclear path, outlined in Sec. II, the integrals over ∆x k and ∆p k must be performed. Using the approximate expressions for the nuclear and electronic action, given by Eqs. (11) and (55), these integrals can be performed analytically from Eq. (26). The result is a product of Dirac delta functions which constrain the dynamics of the nuclear variables to follow Newton's equation of motion.
Here, F k = F 0 (x k ) + F e (Ω(t k ), Ω (t k ), x k ) is the nuclear force at time-step k and is the electronic dependent nuclear force, which depends on the electronic mapping variables of the two paths. This means that the approximate expression for the real-time correlation function when using Stratonovich-Weyl kernels becomes: For two-level systems, the trace in Eq. (58) only contains two terms and is hence easy to perform explicitly. For the generalized theory introduced in Sec. III, where the dimension of the electronic subsystem is F > 2, this trace contains F terms and is also easy to perform for the majority of systems of interest. As for PLDM using spin coherent states, the equations of motion become identical to those for standard PLDM if the equations of motion for the spin-mapping variables are written in terms of Cartesian mapping variables, using the transformation: c λ → Z λ /R s . One difference between the two methods is that the timeevolved Stratonovich-Weyl kernels, which are given by Eq. (51) and enter into the real-time correlation function given in Eq. (58), contains a zero-point energy parameter. This results in the correlation functions which contain the identity operator being treated differently within the underlying theory. The other difference is the integration space, with the initial mapping variables now being confined to different size spin spheres when PLDM is derived using Stratonovich-Weyl kernels. The different spin spheres involve different approximations, given by Eq. (53), in order to form a classical-trajectory picture for the dynamics. Therefore we wish to discover which spin sphere has the least severe approximation and hence most accurately reproduces the dynamics of the QCLE. For the fully linearized spin methods, 37,38 it was found that the W-sphere consistently gave the most accurate results. This is perhaps due to the symmetry that s =s = W and is also consistent with previous work that found that using the zero-point energy parameter for the W-sphere is in some sense optimum for MMST classical trajectories. 34,69,70 Preliminary tests performed by us have also shown that PLDM is more accurate on the W-sphere than on the Q-and P-spheres. Hence from now on, we will only consider PLDM on the W-sphere, which we will refer to as spin-PLDM. In the next section, we outline the underlying spin-PLDM algorithm to simulate dynamics in a F -level electronic subsystem.

III. SPIN-PLDM
At this point we generalize to multi-state problems and also present the method in a way in which it can be most easily implemented. Propagating the generalized spin-PLDM equations of motion is made more convenient by transforming back to Cartesian mapping variables, but retaining the integration space and dynamics derived in Sec. II C. The Cartesian mapping variables, Z λ , are related to the coherent state probability amplitudes, c λ , by: In this expression, the W-sphere radius, R W , is given generally for a F -level electronic system as: 38,71 To express the spin-PLDM correlation function in terms of the Cartesian mapping variables, we first use Eq. (59) to rewrite the time-evolved W-kernel. The matrix elements of this kernel in terms of the Cartesian mapping variables are then: whereÛ (t) is the time-ordered propagator, given by Eq. (25). Assuming that F is not too large,Û (t) is easy to evaluate as it is just the dynamics of a F -level space according to a time-dependent Hamiltonian. Additionally, γ W is the zero-point energy parameter for the Wsphere: 38 Defining the Cartesian mapping variables in this way results in spin-PLDM equations of motion which are identical to those of standard PLDM [Eq. (32)], although the correlation functions will of course still give different results due to the change in the initial distribution. Now that the time-evolved W-kernel has been generated in terms of Cartesian mapping variables, the spin-PLDM correlation function can then be obtained as: where N is a normalization constant for the Cartesian mapping variables, given by: The factor of F 2 in Eq. (63) comes from the fact that the two spin coherent state integrals must each satisfy dΩ = F for multi-state problems. The spin-PLDM algorithm can be easily implemented as follows. First, the integrals contained in the correlation function given by Eq. (63) can be obtained numerically by sampling the initial classical coordinates using Monte Carlo. This means sampling the nuclear phasespace variables from the Wigner distribution of the nuclear initial density [Eq. (10)] and sampling the Cartesian mapping variables for the forward and backward paths (Z and Z ) independently from uniform hyperspheres of radius R W . For each instance of the sampling, the classical coordinates can then be evolved in time using the equations of motion given by Eq. (32) and additionally the time-ordered propagator can be obtained by using Eq. (25), where is the time-step. Eq. (25) involves successive matrix multiplications of F × F matrices, which can be easily performed explicitly as long as F is not too large. Finally, the W-kernel for both the forward and backward paths can then be constructed using Eq. (61) and the contribution to the real-time correlation function for each sample can be obtained by explicitly evaluating the matrix multiplications and trace in the electronic basis. Such an algorithm can also be implemented with only a few minor changes into a preexisting standard PLDM code.
The main differences between spin-PLDM and standard PLDM (whose correlation function is given by Eq. (29)) are as follows. In spin-PLDM, the initial mapping variables are thus constrained to a hypersphere |Z| 2 = R 2 W , which guarantees that the electronic state being represented is correctly normalized. Additionally, spin-PLDM has a zero-point energy parameter, which means that it treats correlation functions containing the identity operator differently from those of traceless operators. It has been shown for MMST mapping that an 'identity trick', 21-23 which treats the identity operator differently within the underlying theory, significantly increases the accuracy of correlation functions withÂ =Î. Fully linearized spin-mapping 37,38 also treats the correlation functions containing the identity operator differently through a zero-point energy parameter.

IV. RESULTS
To test this newly derived spin-PLDM method, we calculate real-time correlation functions for an Ohmic spin-boson model 72 previously considered in Refs. 73 and 21. The parameters used for this spin-boson model correspond to a non-trivial intermediate regime between strongly incoherent decay and coherent oscilla-tions. As in Ref. 21, we calculate in particular C Iσz (t) and C σzσz (t).By taking linear combinations of these two functions, one can recover the more usual measures of state-to-state population transfer. These results for standard and spin-PLDM are given in Fig. 5, along with numerically exact results obtained using the quasiadiabatic path-integral (QUAPI) technique. 74 The results for the mapping-based methods were performed using f = 100 nuclear bath modes and 10 6 trajectories. In Paper II 52 we will present a more complete set of results including multi-state systems with further analysis.
Considering first the identity-containing correlation function, C Iσz (t), the spin-PLDM method produces significantly more accurate results than standard PLDM, which has the incorrect asymptote at long times. For fully linearized methods, spin-mapping was also found to significantly improve the accuracy of identitycontaining correlation functions relative to linearized MMST mapping. 37 In addition, spin-PLDM also computes correlation functions of traceless operators, such as C σzσz (t), more accurately than standard PLDM. In particular, Fig. 5 shows that spin-PLDM corrects for the 'overdamped' oscillations observed in correlation functions of traceless operators calculated using standard PLDM, perhaps due to better preserving of the physical subspace. A more detailed analysis on how these differences effect the calculated correlation functions with spin-PLDM will be presented in the following paper. 52

V. CONCLUSIONS
In this paper we have derived a new partially linearized mapping approach based on classical trajectories, which uses elements of previously derived partially linearized techniques (such as standard PLDM and FBTS) and applies them to the spin-mapping space using the Stratonovich-Weyl transform. We show that this method, called spin-PLDM, exhibits improved accuracy over standard PLDM, when calculating correlation functions. In particular, spin-PLDM corrects for the large error observed in the long-time limit of identity-containing correlation functions calculated using standard PLDM and also appears to solve the issue of overdamped oscillations observed in many dynamical quantities calculated using standard PLDM. This suggests that spin-PLDM is in some sense closer to an exact solution of the QCLE than the standard PLDM method. In Paper II, 52 the spin-PLDM method is tested and analysed further against a wide range of model systems and is compared with the fully linearized spin-mapping method. We also introduce focused initial conditions for spin-PLDM, so as to increase computational efficiency.