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

In the previous paper [J. R. Mannouch and J. O. Richardson, J. Chem. Phys. xxx, xxxxx (xxxx)] we derived a new partially linearized mapping-based classical-trajectory technique, called spin-PLDM. This method describes the dynamics associated with the forward and backward electronic path integrals, using a Stratonovich-Weyl approach within the spin-mapping space. While this is the first example of a partially linearized method derived using such a spin-mapping approach, fully linearized spin mapping has already been established as a method capable of accurately reproducing dynamical observables for a range of nonadiabatic model systems. Alongside numerical simulations for a set of two- and multi-state benchmark systems, we here present a thorough comparison of the terms that arise in the underlying expressions for the real-time quantum correlation functions for spin-PLDM and fully linearized spin mapping, in order to ascertain the relative accuracy of the two methods in computing dynamical observables. In particular, we show that spin-PLDM contains an additional term within the definition of its real-time correlation function, which diminishes many of the known errors that are ubiquitous for fully linearized approaches. We additionally implement focused initial conditions for the spin-PLDM method, which enables us to reduce the number of classical trajectories that are needed in order to reach convergence of dynamical quantities, with seemingly little difference to the accuracy of the result.


I. INTRODUCTION
Trajectory simulations offer a computationally cheap, as well as a physically motivated approach for calculating dynamical quantum-mechanical observables in condensed-phase systems. 1 This explains the current popularity of using such methods, like Ehrenfest dynamics 2,3 and Tully's fewest-switches surface hopping. 4 Certain methods within this class can be rigorously derived from a path-integral formulation of real-time quantum correlation functions. Upon linearizing the difference between the forward and backward nuclear paths, the resulting expression for such quantities is ideally suited for evaluation with classical trajectories. In order to simulate nonadiabatic processes, a quantum-classical approach is needed, by which the electronic dynamics can also be described within the same classical-trajectory picture.
One such category of quantum-classical approaches that fit into this formalism are 'fully linearized' methods because the same linearization approximation used for the nuclear paths is also applied to the electronic paths. Examples of such methods are the linearized semiclassical initial-value representation (LSC-IVR) [5][6][7][8] and the Poisson-bracket mapping equation (PBME) 9,10 approaches, which make use of the Meyer-Miller-Stock-Thoss (MMST) mapping 11,12 to describe the electronic degrees of freedom, as well as fully linearized spinmapping, 13,14 which uses a Stratonovich-Weyl approach a) Electronic mail: jonathan.mannouch@phys.chem.ethz.ch b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch for describing the electronic dynamics within the spinmapping space. Other forms of mapping also exist, from which fully-linearized methods can be derived. 15,16 However, the use of spin mapping to describe the electronic degrees of freedom appears to result in a more accurate quantum-classical method than MMST mapping; for example fully linearized spin-mapping almost always outperforms LSC-IVR and PBME. 13,14 Another category of quantum-classical approaches also exists, referred to as 'partially linearized' methods, so called because the electronic forward and backward paths, unlike the nuclear paths, are explicitly described by separate dynamical variables. An exact solution of the quantum-classical Liouville equation (QCLE) in terms of independent trajectories does not exist 17,18 and hence such methods necessarily provide only an approximate solution for the full dynamics of electrons coupled to classical nuclei. Most notable of these are the partially linearized density matrix (PLDM) [19][20][21][22][23][24][25][26][27][28][29][30] and the forward-backward trajectory solution (FBTS) [31][32][33][34] approaches, where the electronic dynamics of the two paths are described using MMST mapping. 11,12 In Paper I, 35 we developed a new partially linearized method, called spin-PLDM, derived using the Stratonovich-Weyl transform to describe the electronic dynamics within the spinmapping space. Based on results presented in that work and this, spin-PLDM appears to reproduce real-time correlation functions with a greater accuracy than other partially linearized methods.
Fully linearized and partially linearized quantumclassical methods are derived using quite different approximations to the underlying path integrals. It is therefore not a priori obvious which one of these classes of mapping-based classical-trajectory techniques will consistently produce the most accurate results, although re-sults in the literature suggest that PLDM and FBTS typically outperform PBME and LSC-IVR. [36][37][38][39] Hence one focus of the current paper is to compare the terms entering the real-time correlation functions for both spin-PLDM and fully linearized spin mapping, as well as using numerical simulations for a wide range of commonly used model systems, in order to benchmark the accuracy of computing dynamical observables with both methods. Numerical results from other mapping-based classical-trajectory techniques, such as standard PLDM and PBME, are also included for comparison. In particular, we show in the proceeding analysis that spin-PLDM has a term with no analogue in the linearized method, which is contained within the definition of the real-time correlation functions of traceless operators. This term diminishes many of the known errors that are ubiquitous for fully linearized methods.
One factor that often limits the efficiency and applicability of such mapping-based classical-trajectory methods is the issue of sampling. Often a large number of trajectories are required in order to converge results for dynamical quantities of interest, in particular for systems which contain a large number of degrees of freedom. In previous work, focused initial conditions have been developed for mapping-based techniques to alleviate this problem, by restricting the sampling space of initial electronic mapping variables to correspond to populating a single initial electronic state. Such focused conditions have previously been successfully implemented for fully linearized MMST mapping, 9,[40][41][42][43] fully linearized spinmapping 13,14 and partially linearized MMST mapping techniques, 23,[44][45][46] and typically require an order of magnitude fewer trajectories in order to reach the same level of convergence, without significantly reducing the accuracy of the result. We show in this paper that similar focused initial conditions can be rigorously derived and also easily implemented for spin-PLDM. Importantly, we show that the accuracy of the results is practically unaffected by this modification. We hence believe that spin-PLDM is one of the most promising methods for accurately and efficiently calculating dynamical quantities within condensed-phase nonadiabatic systems.

II. FULLY LINEARIZED SPIN MAPPING
Fully linearized spin mapping 13,14 describes the electronic dynamics by evolving a single set of classical variables that are constrained to the surface of a hypersphere, where the possible radii of the hypersphere are determined by the Stratonovich-Weyl Transform. Within this method, the accuracy of any obtained dynamical quantities, such as real-time correlation functions, depends strongly on the spin-sphere radius that is used. In previous work, it has been found that the so-called W-sphere consistently produces the most accurate results for a wide range of model systems. 13,14 It is simplest to analyze the method and perform the calculations using the Cartesian representation for the spin variables. In terms of Cartesian mapping variables, this fully linearized W-sphere method, which we will refer to as spin-LSC, has the following expression for the realtime correlation function: where the t = 0 value is implied for any quantity for which time, t, is not explicitly stated. In this expression, Z = {Z 1 , Z 2 , · · · , Z F } are the Cartesian mapping variables for an F -level system, where Z λ = X λ + iP λ are complex numbers associated with electronic state λ. These mapping variables represent the electronic degrees of freedom of the system. Along with the nuclear phasespace variables x and p, these quantities are initially sampled in the definition of the real-time correlation function from the spin-LSC average, defined as: where dZ = λ dX λ dP λ and ρ W b (x, p) is the Wigner transform of the initial nuclear density matrix, normalized such that dx dp ρ W b (x, p) = 1. The factor of F ensures the correct normalization such that 1 spin-LSC = tr[Î] = F , where tr[· · · ] is the partial trace over the electronic degrees of freedom. Within the definition of the spin-LSC average, the initial Cartesian mapping variables, Z, are constrained to a hypersphere of radius R W : Because the Cartesian mapping variables are sampled uniformly from the hypersphere in Eq. (2), we refer to this as full-sphere initial conditions for the electronic degrees of freedom within the spin-LSC technique. For the spin-LSC correlation function, given in Eq. (1), the electronic operatorsÂ andB are represented by their Stratonovich-Weyl W-functions. The spin-LSC method can also be applied to real-time correlation functions whereÂ andB contain nuclear operators, although we will not consider such correlation functions in this paper. The W-function of operatorÂ is: which is defined in terms of the Stratonovich-Weyl kernel,ŵ W (Z), whose matrix elements are The Stratonovich-Weyl kernel contains a zero-point energy parameter, γ W , which for the W-spin sphere is given by: Inserting Eq. (5) into Eq. (4) leads to an equivalent expression for the W-function of operatorÂ: The spin-LSC correlation function [Eq. (1)] approximates theB operator at time t by evolving in time the Cartesian mapping variables, Z, for the electronic degrees of freedom. These variables, along with the nuclear phase-space variables, are propagated under the following equations of motion: which despite the different derivation are completely equivalent to those of the MMST fully linearized mapping approaches. 13 Within these equations of motion, the nuclear potential has been partitioned into two components: V 0 (x) +V (x). The nuclear potential term V 0 (x) corresponds to the electronic independent nuclear potential, whereasV (x) is the electron-nuclear coupling and is defined such that tr[V (x)] = 0 for all values of x. The electronic independent and dependent nuclear forces, F 0 (x) and F e (Z, x) are given by: where ∇ is the gradient, a vector of derivatives with respect to the nuclear positions. The spin-LSC correlation function given by Eq. (1) has several nice properties. First, although it is in general an approximate theory, the spin-LSC correlation function is always exact at t = 0. This arises from the following property of the W-functions: 14,47 Second, because the electronic dynamics are exactly described by the equation of motion in Eq. (8), the spin-LSC correlation function is also exact in the absence of electron-nuclear coupling. Finally, because |Z| 2 = R 2 W is a constant, spin-LSC correlation functions withB =Î (and hence I W = 1) are also exact: All of these properties are also satisfied by MMST mapping based approaches. However, it has been observed from numerical simulations that spin-LSC generally gives rise to more accurate correlation functions than fully linearized MMST mapping based approaches, in particular forÂ =Î, where C AB (t) does not tend to zero in the long-time limit. For these identity-containing correlation functions, it can be shown that within spin-LSC they have the following simple form: which means that the identity operator is treated exactly within the spin-LSC method. Whereas this appears naturally in fully linearized spin-mapping methods, a modification of MMST mapping was required to obtain a similar result. [48][49][50] As with all mapping-based classical-trajectory techniques, the efficiency of the method is often limited by the number of trajectories one requires to converge the results. Such convergence has been improved previously by using focused initial conditions, which we will now introduce.

A. Focused sampling
Focused sampling has been previously implemented for spin-LSC in Refs. 13 and 14. In this subsection, we present spin-LSC focused initial conditions in a way that makes clear the connection with focused sampling used in spin-PLDM, which will be introduced in Sec. III A. However, for the case of real-time correlation functions with an initial population operator, our approach becomes identical to these previously implemented focused initial conditions.
The spin-LSC method described above makes use of the properties of the Stratonovich-Weyl kernel. In particular, the method relies on the fact that the trace of two operators can be written as an integral over their corresponding W-functions, as given by Eq. (10). However, this property can also be satisfied when using alternative sampling distributions of the Cartesian mapping variables, such as focused conditions: A proof showing that these focused conditions do satisfy the W-functions property in Eq. (10), for any electronic operatorsÂ andB, is given in Appendix A. These focused conditions are similar in form to the full-sphere sampling of the initial Cartesian mapping variables, given by Eq. (2), but with the following differences. Instead of simply constraining the mapping variables to lie on a hypersphere defined by |Z| 2 = R 2 W , the mapping variables in the focused initial conditions are instead now constrained by ρ (λ) foc (Z), the focused sampling distribution for electronic state λ: This distribution, ρ (λ) foc (Z), still constrains the mapping variables onto this hypersphere, but also further constrains the mapping variables so that they entirely occupy electronic state |λ : Additionally, the factor of F within the definition of the full-sphere initial sampling in Eq. (2) is now incorporated within the focused variant into the sum over λ, a complete set of orthonormal electronic states for which the Cartesian mapping variables are focused onto. From the definition of the focused sampling in Eq. (13), the focused spin-LSC correlation function then becomes: The definition of this real-time correlation function is defined so that the focusing, given by Eq. (13) is only performed at t = 0. This is in contrast to the symmetrical quasi-classical (SQC) windowing method, where windowing functions are used to essentially focus trajectories at time t as well. [41][42][43] For real-time correlation functions with an initial population operator (i.e.,Â = |n n|), the focused initial sampling defined in Eq. (2) can be simplified as follows. Using the expression for the Wfunction of a population operator under focused initial conditions, given by Eq. (15), the expression for this focused spin-LSC correlation function then becomes: which is the identical expression for the spin-LSC focused conditions used in Refs. 13 and 14. The fact that only one term in the sum over λ in Eq. (2) contributes to the spin-LSC real-time correlation function whenÂ = |n n| means that focused initial conditions are particularly simple in this case. To treat correlations with off-diagonal (coherence) operators, one must include all of the terms in this sum over λ in Eq. (2). However, focused initial conditions can also be implemented more efficiently in this case by choosing a basis which diagonalizes the initial operatorÂ, as was suggested in Refs. 13 and 51. The constraints imposed within focused initial conditions can easily be implemented using the following parameterization of the Cartesian mapping variables, where r µ=λ = √ 2 + γ W , r µ =λ = √ γ W and φ µ is uniformly sampled from 0 to 2π. In order to illustrate these focused conditions further, we consider the case of a twolevel system, in which the electronic state can be described in terms of the expectation values of the Pauli spin matrices. From Eq. (18), these expectation values for focused initial conditions are: Hence each Cartesian mapping variable is focused onto one of two 'polar circles', as illustrated in Fig. 1, and (φ 2 − φ 1 ) can be thought of as the azimuthal angle around this circle, whereas (φ 1 + φ 2 ) is an unimportant cyclic variable. In addition, the upper arctic circle satisfies [σ z ] W (Z) = 1 and the lower antarctic circle satisfies The focused initial conditions are basis-dependent and hence there is a choice of which complete set of states to focus onto. Focusing onto the electronic states, which diagonalize the nuclear force operator, guarantees that the nuclear dynamics are correct in the absence of offdiagonal diabatic couplings. This is therefore advantageous for systems far from the Born-Oppenheimer limit. Additionally, focusing onto the adiabatic electronic states would be advantageous for systems close to the Born-Oppenheimer limit. Note also that the dynamics can be carried out equivalently in this representation. 52 The model systems that we consider in this paper are all far from the Born-Oppenheimer limit, as this is the regime that is traditionally most difficult for mappingbased classical-trajectory techniques to describe correctly. Hence, from now on we will only use diabatic focusing.
For the spin-LSC method, the difference in the converged results between full-sphere and focused sampling was hardly noticeable in our tests, but the latter required an order of magnitude fewer trajectories. This is consistent with the results presented in Refs. 13 and 14, where the difference between computed correlations functions using full-sphere and focused sampling within the spin-LSC method (referred to in these papers as the Wmethod) are again hardly noticeable. Note that this is by no means true of other mapping approaches; for example spin-mapping on the Q-sphere appears to be fairly accurate when using full-sphere sampling, but with focusing reduces to Ehrenfest dynamics, which is known to be inaccurate. From now on, the numerical results presented in this paper will only be for the focused variant of spin-LSC. However, the proceeding analysis and conclusions that we present are equally valid for both full-sphere and focused variants.

III. SPIN-PLDM
Within the spin-LSC method, the Stratonovich-Weyl kernels are used to describe the observable electronic operators,Â andB, that appear within the real-time correlation function. The electronic dynamics is thus described by a single set of Cartesian electronic mapping 1. An illustration of the two 'polar circles' (red) on the W-spin sphere for a two-level system, from which the Cartesian mapping variables are sampled when using focused initial conditions. The upper arctic circle, having a latitude of cos θc = 1/ √ 3, corresponds to the electronic system solely occupying state |1 (i.e., [σz]W(Z) = 1). The region above this 'arctic circle' corresponds to the system having a negative population associated with state |2 (i.e., [σz]W(Z) > 1).
In contrast, the spin-PLDM method, derived in Paper I, 35 is a partially linearized approach, where the Stratonovich-Weyl kernels are used to describe the forward and backward real-time propagators. This leads to an expression for the real-time correlation function which contains two sets of Cartesian mapping variables, Z and Z : The initial sampling of the Cartesian mapping variables, Z, along with the nuclear phase-space variables x and p, is defined by the spin-PLDM average, given by: Like before, the Wigner transform of the initial nuclear density matrix is defined so that dx dp ρ W b (x, p) = 1. Because each set of Cartesian mapping variables in Eq. (21) is uniformly sampled from the surface of a hypersphere, we refer to this as full-sphere initial conditions for the electronic degrees of freedom within the spin-PLDM technique, like for spin-LSC.
In the spin-PLDM correlation function, the timeevolved W-kernel,ŵ W (Z, t), can be obtained by applying the time-ordered propagator,Û (t), to the left of the Wkernel, given by Eq. (5): where: In addition, t k = kt/N is the time at each time-step of the propagation and N is the number of time-steps. As for all mapping-based classical-trajectory methods, we wish to approach the N → ∞ limit when performing numerical calculations. The Cartesian mapping variables, Z and Z , are evolved in spin-PLDM under the same equations of motion as Eq. (8), except that the expression for the electronic dependent nuclear force is now given by: Hence the Cartesian mapping variables for the forward and backward paths, Z and Z respectively, are coupled in the equations of motion via this nuclear force term. The spin-PLDM correlation function given by Eq. (20) has similar desirable properties to spin-LSC, as well as to MMST based mapping approaches. First, the spin-PLDM correlation function is always exact at t = 0. This arises from the following property of the W-kernel: where the averages in this equation can equally correspond to those for spin-LSC or spin-PLDM, defined in Eqs. (2) and (21) respectively. We have also used that Z and Z are uncorrelated at t = 0. Additionally, as for both spin-LSC and MMST approaches, the spin-PLDM correlation function is also exact in the absence of electron-nuclear coupling. Finally, the correlation function for all these approaches is also exact whenB =Î, which can be shown for spin-PLDM by using the following result: In Paper I 35 it was shown that for an Ohmic spin-boson model, spin-PLDM is also able to accurately reproduce real-time correlation functions, even away from the limits for which the method is formally exact. One potential disadvantage of spin-PLDM is that it typically requires more trajectories than spin-LSC to converge results, because the initial sampling now contains integrals over two sets of Cartesian mapping variables. 53 However, as stated before, such convergence issues can be alleviated using focused sampling of the initial Cartesian mapping variables.

A. Focused sampling
For spin-LSC, the key property of any allowed initial mapping variable sampling must be that it preserves the properties of the Stratonovich-Weyl kernel in Eq. (10). For spin-PLDM, we require that the focused sampling can correctly describe electronic quantum operators in terms of the W-kernel. In other words: (27) This is also satisfied by the focused initial conditions defined by Eq. (13), as can be proved using the same arguments that are presented in Appendix A. In spin-PLDM, two sets of Cartesian mapping variables are used, in order to represent both the forward and backward electronic paths. Hence in focused spin-PLDM, both sets of mapping variables must be initially sampled independently using these focused conditions. This leads to the following expression for the correlation function: where the ensemble average now amounts to using focused conditions for both sets of Cartesian mapping variables: In a similar fashion to how the focused conditions were defined for spin-LSC, ρ (λ) foc (Z) additionally constrains the Cartesian mapping variables so that they entirely occupy electronic state |λ . Additionally, the factor of F 2 within the definition of the spin-PLDM full-sphere sampling, given in Eq. (21), is incorporated within the focused variant by two independent sums over the complete set of electronic states (|λ and |λ ).
In contrast to the focused spin-LSC method, the focused spin-PLDM method contains terms where the initial Cartesian mapping variables for the forward and backward propagator paths are focused onto different electronic states, |λ and |λ . Because fully linearized methods only consider the average of these two paths, such configurations cannot be represented in spin-LSC. The focused conditions previously implemented for standard PLDM and FBTS 46 also focus both the forward and backward paths onto the same initial electronic state, which can lead to poor results for the long-time dynamics or for systems with relatively strong electron-nuclear coupling. 23,46 Because the Stratonovich-Weyl kernels are used to represent the forward and backward propagators in spin-PLDM, rather than the observable operators, none of the terms in the {λ, λ } sum in Eq. (29) are identically zero for the correlation function given in Eq. (28), even whenÂ = |n n|. This means that although focused spin-PLDM is a much more efficient method than performing spin-PLDM with full-sphere sampling, the focused conditions for spin-PLDM are nevertheless less efficient than for spin-LSC. 54 As for the spin-LSC focused conditions, the sampling of focused spin-PLDM is now basis-dependent. However, using either a diabatic or adiabatic basis seems the obvious choice, depending on whether the system is close to or far away from the Born-Oppenheimer limit. As stated before, we will only use diabatic focusing in this paper.

IV. COMPARISON BETWEEN SPIN-PLDM AND SPIN-LSC
In spin-LSC, the W-kernels are used to represent the observable operators,Â andB, within the definition of the real-time quantum correlation function, while in spin-PLDM, the W-kernels are used to represent the timeordered propagators for the forward and backward paths. Hence because the spin-LSC and spin-PLDM methods use the W-kernels to represent different objects within the real-time correlation function, the associated underlying expressions for both techniques, given by Eqs. (16) and (28), look quite different. In order to compare the methods to each other, we can re-express the spin-PLDM correlation function in terms of the mapping-variable representations of the underlying operatorsÂ andB: We choose this definition of the spin-PLDM operator representation so that it becomes identical to the standard PLDM operator representation, B m (Z, Z ), when γ W = 0 and also identical to the spin-LSC operator representation [Eq. (7)] when Z = Z .
Here, we only consider correlation functions that satisfy tr[B] = 0, because any operator can be decomposed into a traceless part and the identity and the realtime correlation function C AI (t) is time-independent. To make the comparison, the identity-containing correlation functions (withÂ =Î) and the correlation functions of traceless operators (with tr[Â] = 0) are considered separately. The discussion in this section can be equally applied to either the full-sphere or focused variants of the methods.
A. Identity-containing correlation functions From Eq. (22), the time evolution of the spin-PLDM correlation function in general does not just depend on the time-evolved mapping variables, but also on the timeordered electronic propagator,Û (t). However, when A =Î, such identity-containing correlation functions can be exactly written in terms of the evolved mapping variables: as derived in Appendix B. In this expression, B W (Z) is the W-function for an electronic operatorB, given by Eq. (7) and B W (Z, Z ) is the spin-PLDM expression for the same operator, given by Eq. (30). The term I W (Z, Z ) = 1 2 (Z * · Z − R 2 W ) + 1, which is also defined by Eq. (30), acts as an 'overlap factor', between the Cartesian mapping variables for the forward and backward electronic paths(Z and Z respectively).
There are two noteworthy points of Eq. (31). First, if we set the zero-point energy parameter, γ W , to zero, the expression for the identity-containing correlation function depends entirely on the standard PLDM expressions for the underlying electronic operators, I m (Z, Z ) and B m (Z (t), Z(t)). This illustrates that spin-PLDM differs from standard PLDM not just through the initial sampling of the Cartesian mapping variables, which are constrained to the hypersphere, but also through the existence of a zero-point energy parameter, which results in additional terms being present in the expression for the identity-containing correlation function.
Second, in the case that Z = Z , so that the method has only one set of Cartesian mapping variables as for fully-linearized techniques, the form of Eq. (31) becomes similar to the expression for the identity-containing correlation function within the spin-LSC method, given by Eq. (12). This follows from: . This result is perhaps unsurprising, because Eq. (12) can be derived as an 'identity trick' 48,49 for spin-PLDM, as shown in Appendix B 1. Hence unlike standard PLDM, both spin-LSC and spin-PLDM treat the identity-containing correlation functions on a similar footing.

B. Correlation functions of traceless operators
In general, the expression for the correlation functions of traceless operators within spin-PLDM cannot be easily simplified, although it is still simple to evaluate numerically. This is because the evolution of the zero-point energy parameter within Eq. (22) depends on the timeordered electronic propagator, which cannot be expressed directly in terms of the time-evolved Cartesian mapping variables. However, the correlation functions of traceless operators within spin-PLDM can be written as follows: as shown in Appendix C. In this expression ∆C AB (t) is a term, which contains the contribution to the correlation function from the propagation of the zero-point energy parameter. This separation is defined so that both ∆C AB (0) = 0 and in the absence of electron-nuclear coupling that ∆C AB (t) = 0 . While a convenient and general expression for ∆C AB (t) is hard to obtain, such an 2. An illustration of the ζ λ mapping variables for an F = 2 system, which lie on the same hypersphere as Z λ , but are also orthogonal to it. The corresponding |ζ state is the inversion of |Z through the origin.
expression is easier to derive if just two-level systems are considered.
For two-level systems, any operator can be written in terms of a spin coherent state and a state orthogonal to it. We therefore define a state, |ζ , whose Cartesian mapping variables lie on the same hypersphere as |Z , but are also orthogonal to it: It can be shown that |ζ corresponds to the inversion of |Z through the origin. In other words: This relationship between the Z λ and ζ λ mapping variables is shown pictorially in Fig. 2. Using these two sets of mapping variables, the term ∆C AB (t) can be rewritten as: where the ζ λ mapping variables are in the same way orthogonal to Z λ . The derivation of this expression is outlined in Appendix C. Such a term has never been included within a mapping-based classical-trajectory technique before. For example, Eq. (35) is zero when γ W = 0 and hence such a term is not present within standard PLDM. In addition, the first term on the right hand side of Eq. (32) has a similar form to the full spin-LSC correlation function [Eq. (1)] when Z = Z and therefore spin-LSC also has no analogue to ∆C AB (t). This additional ∆C AB (t) term therefore constitutes the main difference between the spin-LSC and spin-PLDM correlation functions, in addition to the fact that both methods also contain different numbers of Cartesian mapping variables.
The presence of ζ and ζ Cartesian mapping variables in Eq. (35) confirms that this additional term contains effects arising physically from propagating the electronic subsystem on the opposite side of the spin-sphere to that of the coherent state from which the nuclear force is calculated. Because such a term has been rigorously derived within a Stratonovich-Weyl approach to spin-mapping, this suggests that spin-mapping based techniques which only describe the electronic dynamics at a single point of the spin-sphere are in some sense deficient; electronic dynamics at the 'antipode' is therefore necessary in order to give a fuller description of the exact dynamics within coupled electron-nuclear systems. In quantum mechanics, the nuclear degrees of freedom can induce instantaneous excitations of the electronic subsystem and hence the various dynamics of the 'antipode' do correspond to physical allowed processes within the real system. Hence spin-PLDM is able describe physical processes which are completely neglected in other mapping-based classicaltrajectory techniques.
For systems with an arbitrary number of electronic states, the term ∆C AB (t) will involve propagation of a set of F − 1 orthogonal states to Z and hence this term physically corresponds to propagating all possible instantaneous excitations of the underlying coherent state. In the next section, the properties and effects of this additional term will be investigated numerically.

V. RESULTS
In this section, we consider a set of challenging systems far away from the Born-Oppenheimer limit. Although all the spin-mapping results presented in this section are obtained using focused sampling, we have found that they hardly differ from the same results obtained using fullsphere sampling. For the MMST mapping-based techniques, only the most accurate form of the methods are used, so as to act as a fair comparison with spin-mapping; hence we present results only for the non-focused variants of these methods.

A. The spin-boson model
To test the accuracy of our focused spin-PLDM method, we have applied the method to a series of spinboson models. We will first consider the same spin-boson model that was studied in Paper I 35 and then we will consider other parameter regimes, in order to try and cover a large part of the parameter space, including both symmetric and asymmetric systems, low and high temperature limits and strong and weak system-bath coupling. Such spin-boson models are commonly used to benchmark new methods, as numerically exact results are available for comparison. The model essentially consists of two electronic states coupled to an initially thermalized harmonic bath, whose Hamiltonian is given by: 55 Here, ε is the energy bias and ∆ is the constant diabatic coupling. The bath contains f nuclear modes, each with frequency ω j and electron-nuclear coupling coefficient c j . The mass, m j , has been incorporated into the definition of the nuclear position,x j , and momentum,p j , operators. As the nuclear bath is harmonic, the dynamics of the spin-boson model is exactly described by the quantum-classical Liouville equation (QCLE). 56,57 Hence any errors arising from calculating correlation functions using mapping-based classical-trajectory techniques arise solely from approximations to the QCLE and are not inherently due to the classical treatment of the nuclear degrees of freedom. The spectral density, J bath (ω), determines the distribution of nuclear frequencies within the bath. One of the most commonly used spectral densities is the Ohmic bath form: where ω c is the characteristic frequency and ξ is the Kondo parameter. In order to perform numerical simulations, the continuous bath must first be discretized into a finite number of modes. The spectral density can then be represented as follows: We have used the discretization scheme employed in Ref. 58. Additionally, we consider solely correlation functions in which the uncoupled bath is initially thermalized. This means that for mapping-based classical-trajectory techniques, the initial nuclear coordinates are sampled from the following Wigner distribution: where α j = tanh( 1 2 βω j ) and β = 1/k B T . For comparison, numerically exact results for the real-time quantum correlation functions have been obtained using the quasiadiabatic path-integral (QUAPI) technique. 59 We begin by comparing the C Iσz (t) and C σzσz (t) correlation functions obtained with our focused spin-PLDM method, to those for focused spin-LSC. All calculations are performed with f = 100 nuclear degrees of freedom. The differences in the underlying expressions for these correlation functions has already been discussed in Sec. IV. To illustrate the differences numerically, we will first discuss the Ohmic spin-boson model studied previously in Refs. 37, 48 and Paper I. 35 This model is both asymmetric and at a relatively low temperature, which makes it one of the more challenging spin-boson models to study. The parameters used for this spin-boson model, given by the first row in Table I, also correspond to an intermediate regime between strongly incoherent decay and coherent oscillations.
First, comparing the focused spin-PLDM results in Fig. 3 (orange) with the spin-PLDM results in Fig. 1 of Paper 1 (orange) shows that using focused sampling in spin-PLDM does not significantly change the obtained results. As there is no disadvantage in using focused initial conditions, we will only present focused spin-PLDM results in this paper. Note that Refs. 13 and 14 also show that focusing appears to barely change the results of spin-LSC (as long as the W-sphere is used). In Fig. 3, we compare these correlation functions for the focused spin-PLDM and focused spin-LSC methods. For the identitycontaining C Iσz (t) correlation function, both methods produce similarly accurate results. This is not surprising, as the expression for the C Iσz (t) correlation functions for both methods, given by Eqs. (12) and (31), have a similar form. Additionally, the focused spin-LSC expression (given by Eq. (12)) can be derived from an 'identity trick' applied to focused spin-PLDM, as shown in Appendix B 1, which suggests that both methods treat the identity-containing correlation functions on an equal footing. Therefore Fig. 3 shows that having two electronic mapping variables in focused spin-PLDM only results in a small improvement to the C Iσz (t) correlation function for this parameter regime.
For the C σzσz (t) correlation function, the solid green line shows the focused spin-PLDM result, given by Eq. (32), without including the term ∆C σzσz (t). Hence the only difference between this result and the C σzσz (t) correlation function for the focused spin-LSC method is that the focused spin-PLDM expression contains two electronic mapping variables.
As for the identitycontaining correlation functions, having two sets of mapping variables only leads to a small improvement in the C σzσz (t) correlation function for this model. The main improvement of focused spin-PLDM over focused spin-LSC is hence the term ∆C σzσz (t), which corrects for the overdamped coherences observed in the correlation functions of traceless operators calculated using focused spin-LSC.  To further test the spin-PLDM method, we apply the method to a selection of other spin-boson problems. The parameters for the spin-boson models we consider are given by rows (a)-(f) in Table I. For these calculations, we have again used f = 100 nuclear degrees of freedom, except for the strong coupling systems (e) and (f), where f = 400 nuclear degrees of freedom were needed for convergence. Such a set of spinboson models cover a wide range of different physical regimes, from symmetric to asymmetric, weak to strong system-bath coupling and low to high temperature limits. Hence these models offer a comprehensive test for our new focused spin-PLDM method. Additionally, all of these spin-boson models lie in a challenging regime for mean-field, mapping-based classical-trajectory methods to describe correctly, which is far from the Born-Oppenheimer limit. This set of spin-boson models have already been used to test other mapping-based classicaltrajectory methods. 13,50,60,61 While Ehrenfest and standard linearized mapping methods are known to perform relatively well for the symmetric models, these methods fail to correctly predict the long-time populations in the asymmetric models. As in Ref. 50, we also make a point of computing the time-dependent coherences, as well as populations. In this paper, Eq. (39) initializes the nuclear coordinates from the Boltzmann distribution of nuclear potential V 0 (x). This is subtly different from Refs. 60 and 61, where the nuclei are sampled from the Boltzmann distribution for potential V 0 (x) + 1|V (x)|1 . Nonetheless, these different bath initial conditions do not appear to lead to significant differences in the results. Fig. 4 shows the calculated dynamical expectation values of all Pauli spin matrices for a range of mappingbased techniques. All results correspond to the system initially occupying the higher energy diabatic state, |1 . The Poisson-bracket mapping equation (PBME) approach is the fully linearized version of standard PLDM 62 and is described in Ref. 63. While we do not show results for LSC-IVR in this paper, it has been observed that both PBME and LSC-IVR produce real-time correlation functions to a similar degree of accuracy. 48 Table I. The dashed black lines give the numerically exact results, obtained using QUAPI.
the weak system-bath coupling spin-boson models at high temperature (given by (a) and (c)). We would expect to describe the dynamics of these models correctly, as the classical approximation that is employed in deriving mapping-based classical-trajectory techniques should be valid at high temperatures. We find that only PBME is unable to accurately describe the dynamics within these spin-boson models. The error arises because PBME does not correctly calculate the identity-containing correlation functions. This problem however can be alleviated by using an 'identity trick' recently derived for PBME 48,49 or by focusing trajectories at time t using SQC. [41][42][43] The fact that the focused spin-PLDM and focused spin-LSC results are similar for these models illustrates that the term in focused spin-PLDM, ∆C σzσz (t), can be neglected at high temperatures. The second set of mapping variables within focused spin-PLDM also offers no improvement in the accuracy of the result for this parameter regime. Additionally, because the standard PLDM approach also correctly describes dynamical properties of the system in this parameter regime, the results are insensitive to the distribution from which the initial Cartesian mapping variables are sampled from. Next, we consider models with weak system-bath coupling at low temperatures (given by (b) and (d)). For these models, spin-PLDM appears to consistently produce the most accurate results. As discussed above, focused spin-PLDM outperforms focused spin-LSC for these systems due to the term, ∆C σzσz (t), which corrects for the overdamped coherences found in the focused spin-LSC correlation functions of traceless operators. Additionally, focused spin-PLDM also outperforms standard PLDM in this parameter regime, because restraining the Cartesian mapping variables to a hypersphere fixes the large errors observed in the identity-containing correlation functions. This different initial condition for the mapping variables also corrects for the overdamped co-herences found in the standard PLDM correlation functions of traceless operators.
Finally, systems (e) and (f) are the most challenging spin-boson models that we consider in this paper, due to the strong system-bath coupling. In fact, the coupling in system (f) is so strong that we struggled to converge the numerically exact QUAPI results. Hence we only present the QUAPI results for times for which we are certain that they are numerically exact. The main conclusion from these figures is that spin-PLDM is the best method at consistently obtaining the short time dynamics correctly and also minimizes the error of the dynamics in the long-time limit. The figures show that spin-PLDM provides a larger correction to spin-LSC for the strong-coupling spin-boson models. This improvement even holds in the high temperature limit, where spin-LSC was already quite good (i.e., model (f), where ∆C σzσz (t) ≈ 0). This suggests that the two sets of electronic mapping variables in focused spin-PLDM are particularly important for describing the dynamics in systems with relatively strong system-bath coupling.

B. The Fenna-Mathews-Olsen complex
To demonstrate that the spin-PLDM method can be easily applied to multistate systems, we have applied the method to a F = 7 model for the Fenna-Matthews-Olsen (FMO) complex, which is a commonly studied photosynthetic pigment-protein complex found in green sulfur bacteria. 65 This is a challenging benchmark problem that has already been used to test a wide range of mapping-based techniques. 14,19,20,22,26,28,43,46,49,[66][67][68][69][70][71][72][73] However, as for the spin-boson model, the dynamics of this FMO model is exactly described by the QCLE, because the nuclear bath is harmonic, which means that nuclear quantum effects do not need to be included in order to describe the quantum dynamics of this model exactly. 33,56,57 Numerically exact results have been computed for this model from the hierarchical equation of motion (HEOM) approach. [74][75][76][77][78][79] The Hamiltonian of this model is: 74Ĥ where the electronic system Hamiltonian is given, in units of cm −1 , as: 80 (41) Each diabatic electronic state is coupled to its own independent harmonic bath that consists of f nuclear modes with frequencies ω j . The mass-weighted Hamiltonian for all seven independent baths is therefore given by: where each bath has the same set of frequencies, ω j . The model then contains a linear coupling term, which connects each harmonic bath to its corresponding diabatic state as follows: where c j are the exciton-nuclear coupling coefficients. Such a Hamiltonian can be partitioned into a stateindependent nuclear potential, V 0 (x), and a traceless electronic operator,V (x), as follows: where tr[· · · ] signifies the partial trace over the electronic degrees of freedom. The spectral density, J bath (ω), determines the distribution of nuclear frequencies within each of the baths and their couplings. This model employs the Debye spectral density: where ω c is the characteristic frequency of the bath (with τ c = ω −1 c its corresponding time scale) and Λ is the reorganization energy. We use Λ = 35 cm −1 , to be consistent with previous work, and only consider the situation of a 'fast bath' (i.e., τ c = 50 fs), which is the most difficult previously studied parameter regime. In order to perform numerical simulations, the continuous bath must be first discretized into a finite number of modes, as shown in Eq. (38). In this paper for the Debye spectral density, we have used the discretization scheme employed in Ref. 81 with f = 60 bath modes per state. We consider the dynamics after the initial excitation of a single site, with the baths in thermal equilibrium before the excitation. This means that for mapping-based classicaltrajectory techniques, the initial nuclear coordinates associated with each excitonic site are sampled from the same Wigner distribution given by Eq. (39).
To test the accuracy of our newly developed focused spin-PLDM method, we calculate the time-dependence of the exciton dynamics for this FMO model. Of particular interest is the picosecond exciton population transfer to the lowest energy site/chromophore (n = 3), from which the excitons can be harvested at the reaction centre. We consider both the dynamics at high and low temperature (300 K and 77 K respectively) and also with the initial exciton residing on different sites/chromophores of the complex (i.e., site 1 and 6). Fig. 5 shows the population dynamics of the FMO at high temperature (300 K). As stated before, this is the regime in which the approximations involved in deriving mapping-based classicaltrajectory techniques are expected to be valid. From  Fig. 5, it can be seen that focused spin-LSC, standard PLDM and focused spin-PLDM are each able to describe the short time dynamics very well (i.e., up to 1 ps), irrespective of whether the exciton is initialized on site 1 (the top row of figures) or site 6 (the bottom row of figures). Even for this short-time behaviour, the accuracy of the focused spin-PLDM approach however appears superior, as the calculated dynamics from this technique are essentially indistinguishable from the exact HEOM results (solid lines) for this high temperature parameter regime. The accuracy of focused spin-PLDM continues to be good, even when the long-time dynamics (up to 10 ps) are considered. The focused spin-LSC method also gives an excellent description of the long time limits of the exciton populations. However, standard PLDM is unable to describe the long time dynamics correctly for this high temperature regime and exhibits a large deviation from the exact result, in particular for the site 3 population. This illustrates a clear advantage in using spin-mapping instead of the standard MMST approach. Constraining the mapping variables to a hypersphere ensures that these Cartesian mapping variables remain in the physical subspace during the dynamics and this appears to minimize the errors produced in the long-time limit and also allows the method to better approximate the correct Boltzmann distribution.
The more challenging regime corresponds to studying the dynamics of the model at a low temperature (77 K). Fig. 6 compares the calculated population dynamics for the same mapping-based classical-trajectory techniques.
In the short-time dynamics, focused spin-LSC now shows significant errors compared to the exact HEOM results. Additionally some of the calculated populations with focused spin-LSC now become negative, even at quite short times (see for example the n = 6 population, when starting in site 1), which is unphysical. While the standard PLDM results appear fairly accurate at short times, the method suffers from overdamped coherences, similar to what was observed when applying the method to the spin-boson model. Focused spin-PLDM produces almost exact results for this short-time dynamics, which again illustrates its apparent superiority. In the long-time limit (i.e., up to 10 ps), the standard PLDM approach again leads to large errors in the populations. In contrast, focused spin-LSC appears to produce less severe errors in the populations in the long-time limit. In particular, the n = 3 population is well reproduced by the method. This is however probably partly fortuitous, due to the method's inability to describe the short-time dynamics of this population correctly. The errors within focused spin-PLDM in the long-time limit also appear not too severe and the method also suffers less severely with negative populations compared to focused spin-LSC. Standard PLDM, in contrast, always calculates positive populations, because the theory has no zero-point energy (γ) parameter. Although the difference is not as great as for the 300 K case, nonetheless the absolute errors in the populations are larger than those for spin-PLDM, which again shows itself to be the most accurate of the three methods.

VI. CONCLUSIONS
In this paper we have tested the spin-PLDM method derived in Paper I 35 on a range of model systems in various different regimes. Comparison has also been made with other mapping-based classical-trajectory techniques.
Spin-PLDM appears to be able to improve upon previously derived fully linearized approaches such as spin-LSC and also upon the standard PLDM approach based on MMST mapping. In particular, we have shown that spin-PLDM contains a term which is absent in spin-LSC, and which appears to solve the problem of 'overdamped coherences' observed in the correlation functions of traceless operators calculated using spin-LSC. A focused sampling variant of spin-PLDM is also introduced, which reduces the number of trajectories needed to converge results, while still producing dynamical observables to the same level of accuracy. For the FMO model, 10 7 trajectories were needed to fully converge the focused spin-PLDM results compared to 10 6 trajectories for the other methods; however we also observed that only 10 4 focused spin-PLDM trajectories were needed to qualitatively reproduce the main features of the dynamics. Even though focused spin-PLDM generally requires an order of magnitude more trajectories than other mapping-based classical-trajectory methods, the superior accuracy of the method makes the added computational expense worthwhile, while more efficient sampling schemes could easily be implemented in the future in order to reduce the number of required trajectories. 82 In general, our results show that generalized spin-mapping methods outperform their MMST analogues and partially linearized methods appear superior to fully-linearized ones; hence spin-PLDM is the best mapping method of them all. As every model used in this paper could in principle (although not in practice) be exactly solved using the quantum-classical Liouville equation (QCLE), the fact that spin-PLDM is able to consistently produce such accurate results suggests that the theory is in some sense closer to an exact solution of the QCLE. The superior accuracy of spin-PLDM is evident not just for calculating electronic population dynamics, but also for coherences.
In particular, spin-PLDM seems to reproduce the rel-atively short time properties of real-time quantum correlation functions extremely accurately compared to other mapping-based classical-trajectory techniques. Because of this there are a number of applications which may therefore be well suited for spin-PLDM. For instance, memory kernels within the generalized quantum master equation formalism normally decay relatively rapidly and hence it is expected that such quantities can be accurately computed using spin-PLDM. Other classicaltrajectory based techniques have already been successively used to calculate memory kernels, in order to accurately obtain the long-time dynamics of real-time correlation functions. 51,[83][84][85] Dipole-dipole correlation functions and other optical response functions, from which linear and non-linear spectra can be obtained, also often have short coherence times and hence perhaps can also be accurately computed accurately using spin-PLDM. Several questions still remain unanswered. First, is there some framework by which the accuracy of mappingbased techniques can be compared, so that some definitive consensus can be made on which techniques are the most reliable and accurate? Second, how close can mapping-based classical-trajectory techniques get to an exact solution of the QCLE? In other words, what would be the ultimate mapping-based technique and in what limits would it be able to exactly describe the dynamics of systems?

ACKNOWLEDGMENTS
The authors would like to acknowledge the support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network. We also thank Johan Runeson for helpful discussions and for his comments made on the original manuscript.

DATA AVAILABILITY
The data that supports the findings of this study are available within the article.

Appendix A: Focused initial conditions
In this appendix, we prove that the focused initial conditions introduced in Sections II A and III A satisfy the Stratonovich-Weyl property given by Eq. (10) for any offdiagonal (coherence) operators. While this is unimportant for spin-LSC when calculating correlation functions involving an initial population operator, this is nevertheless necessary for spin-PLDM, because the W-kernels are used to represent the propagators rather than the observ-able operators. In particular, we wish to show that: spin-LSC = 1, (A1) where the states |m and |n are chosen to be a complete set of orthonormal basis states for the electronic system (i.e, m|n = δ mn ). Additionally, the W-functions associated with these electronic coherence operators, [|m n|] W (Z), can be obtained from Eq. (7) as: Within the focused conditions, the constraints on the Cartesian mapping variables imposed by the projection function, ρ (λ) foc (Z) can be easily implemented using the expression for Z µ = X µ +iP µ given in Eq. (18). Inserting this into the definition of the phase space average for focused initial conditions, given by Eq. (13), gives rise to the following result: The first term on the right hand side of this expression comes from the two terms within the sum over λ in Eq. (13), where the Cartesian mapping variables are focused onto states m and n. The second term on the right hand side of Eq. (A3) then comes for the remaining F − 2 terms within this sum over λ. The result in Eq. (A3) can be further simplified by using the expression for the zero-point energy parameter, γ W , given in Eq. (6): Notice that the value of the t = 0 focused real-time correlation function involving off-diagonal (coherence) operators depends on the spin-sphere radius, R W . For the special case of the W-sphere, with the corresponding value for R 2 W given by Eq. (3), the value of this correlation function becomes: which satisfies the Stratonovich-Weyl property given by Eq. (A1). Hence off-diagonal (coherence) operators can be correctly described by these focused initial conditions on the W-sphere.
Appendix B: Spin-PLDM identity-containing correlation functions Within this paper, we show that the expression for the spin-PLDM real-time correlation function can be simplified in the case whenÂ =Î. In order to do this, the following expression for the W-kernel can be used: which is simply an analogous expression to Eq. (22). Inserting this expression for the W-kernel into the expression for the spin-PLDM identity-containing correlation functions [Eq. (20), withÂ =Î] results in: Additionally, because the trace is invariant to cyclic permutations, the last term on the right hand side of Eq. (B2) is zero, using again that tr[B] = 0. To simplify the third and fourth term on the right hand side of Eq. (B2), we make use of the following property involving the W-kernel: which follows from Eq. (B1). Inserting Eq. (B3) into the third and fourth term on the right hand side of Eq. (B2) shows that the time-dependence of C IB (t) within spin-PLDM can be described entirely by evolving the Cartesian mapping variables, Z(t) and Z (t). Finally, using the definition of the W-function of operatorB, given by Eq. (4), results in the expression for C IB (t) given by Eq. (31).

The spin-PLDM 'identity trick'
An alternative expression for the identity-containing correlation function can also be obtained by usingÂ =Î right at the beginning of the derivation of the spin-PLDM correlation function. In Paper I, 35 it was shown that applying the linearization approximation to the nuclear degrees of freedom within the path-integral representation of the real-time correlation function, C AB (t), leads to the following expression: In the case whereÂ =Î, the summation over the λ and λ indices in Eq. (B4) can be performed explicitly. Hence, the contribution to the real-time quantum correlation function from the electronic transition amplitudes, where we have usedÎ = λ |λ λ|. This expression can then be approximated in terms of Cartesian mapping variables by insertingÎ = F N dZ 0ŵW (Z 0 )δ(|Z 0 | 2 −R 2 W ) in between the e ±V (x1) operators. Here, N is a normalization constant for the Cartesian mapping variables, given by N = dZ 0 δ(|Z 0 | 2 − R 2 W ). Performing the same steps as in Paper I 35 leads to the following expression: (B5) Now the time-dependence of the identity-containing correlation function is completely described by the timeevolved Cartesian mapping variables, because timeordered propagators are now positioned either side of the W-kernel, as in Eq. (B3). Additionally, the electronic action is defined as: where t k = k is the time at time-step k. Hence evaluating the identity operator explicitly before inserting the Cartesian mapping variables now leads to an expression which contains a single set of Cartesian mapping variables. Following the rest of the spin-PLDM derivation with Eq. (B5) leads to an expression for the identitycontaining correlation function which is identical to that for spin-LSC, given by Eq. (12). The fact that the form of the identity-containing correlation function within spin-LSC can be derived from the underlying spin-PLDM equations suggests that both spin-LSC and spin-PLDM will obtain the identity-containing correlation function with a similar accuracy. This is indeed in agreement with the spin-boson results presented in Section V.
Such an 'identity trick' can also be performed for fully linearized MMST mapping. 48 For example, the same analysis performed in deriving the spin-PLDM 'identity trick' presented here can be equally applied to the derivation of the identity-containing correlation function within LSC-IVR (except that an additional linearization approximation within the MMST mapping space must also be performed). This results in the 'single unity' method introduced in Ref. 48, which was found to offer a significant increase in accuracy when calculating identity-containing correlation functions with fully linearized MMST mapping based techniques. [48][49][50] Appendix C: Spin-PLDM correlation functions of traceless operators For the correlation functions of traceless operators within spin-PLDM, the time-dependence does not in general entirely arise from the time-evolution of the Cartesian mapping variables, Z(t) and Z (t). This is because the time-evolved W-kernel, given by Eq. (B1), contains a contribution from the time-ordered propagator, U (t), proportional to the zero-point energy parameter, γ W . However, in the absence of electron-nuclear coupling, the time-dependence of the correlation functions of traceless operators can be given solely in terms of the time-evolution of the Cartesian mapping variables. This is because PLDM, derived using a spin coherent state basis without a zero-point energy parameter, is also exact in this case: 35 C AB (t) ≈ dx dp dΩ dΩ ρ W b (x, p) Ω|Â|Ω Ω (t)|B|Ω(t) .
(C1) Eq. (C1) is also always exact at t = 0 for any value of the electron-nuclear coupling. In this expression, |Ω = λ c λ |λ is a spin coherent state, given in terms of the amplitudes c λ . In addition, dΩ is defined by Eq. (35) in Paper I. 35 Eq. (C1) can be rewritten in terms of Cartesian mapping variables by using the relation: c λ = Z λ /R W . This results in: When t = 0 or for systems with electron-nuclear coupling, Eq. (C2) is no longer exact and the spin-PLDM correlation function must instead be used to obtain accurate results. The difference between Eq. (C2) and the spin-PLDM correlation function given by Eq. (B1) can be incorporated into a term, ∆C AB (t), which generally must be calculated in order to accurately obtain correlation functions of traceless operators. This term is hence given by: which through the definition of the W-kernel in Eq. (B1) depends on the time-ordered propagator,Û (t).
For the case of F = 2, Eq. (C3) can be further simplified. First, the electronic identity operator can be expressed in terms of Cartesian mapping variables: The first term in this expression corresponds to the spincoherent state outer product, |Ω Ω|, as can be seen from the relation c λ = Z λ /R W . The second term contains new Cartesian mapping variables, ζ λ , which as defined in Eq. (33) to be orthogonal to Z λ . Defined in this way, the ζ λ mapping variables correspond to the electronic state on the opposite side of the spin-sphere to |Ω , as illustrated in Fig. 2. Hence this second term is just the outer product of an electronic state orthogonal to the coherent state, |Ω . Applying the time-evolved propagator, U (t), to the left of Eq. (C4) leads to: (C5) Inserting this expression for the time-evolved propagator into the expression for the time-evolved W-kernel, given by Eq. (B1), leads to a new expression for the W-kernel solely in terms of Cartesian mapping variables: where we have additionally used R 2 W = 2γ W + 2, which is valid when F = 2. Inserting Eq. (C6) into the definition of ∆C AB (t), given by Eq. (C3) and expanding the terms leads to the expression for this term given in Eq. (35). We have additionally used the expression for the spin-PLDM representation of the operatorB in terms of the Cartesian mapping variables, B W (Z, Z ), given by Eq. (30) and we have again used that tr[Â] = tr[B] = 0.