Generalized spin mapping for quantum-classical dynamics

We recently derived a spin-mapping approach for treating the nonadiabatic dynamics of a two-level system in a classical environment [J. Chem. Phys. 151, 044119 (2019)] based on the well-known quantum equivalence between a two-level system and a spin-1/2 particle. In the present paper, we generalize this method to describe the dynamics of $N$-level systems. This is done via a mapping to a classical phase space that preserves the $SU(N)$-symmetry of the original quantum problem. The theory reproduces the standard Meyer--Miller--Stock--Thoss Hamiltonian without invoking an extended phase space, and we thus avoid leakage from the physical subspace. In contrast with the standard derivation of this Hamiltonian, the generalized spin mapping leads to an $N$-dependent value of the zero-point energy parameter that is uniquely determined by the Casimir invariant of the $N$-level system. Based on this mapping, we derive a simple way to approximate correlation functions in complex nonadiabatic molecular systems via classical trajectories, and present benchmark calculations on the seven-state Fenna--Matthews--Olson complex. The results are significantly more accurate than conventional Ehrenfest dynamics, at a comparable computational cost, and can compete in accuracy with other state-of-the-art mapping approaches.


I. INTRODUCTION
The full quantum dynamics of complex systems is in general far too complicated to be simulated in practice. Instead it is often necessary to separate the problem into a (smaller) subsystem that is treated quantummechanically and an environment that can be approximated by classical dynamics. In chemistry, the typical example is to treat a molecular system as a subsystem of N electronic levels coupled to an environment of classical nuclear modes. If the coupling between the electronic and nuclear motion cannot be neglected, methods based on the standard Born-Oppenheimer approximation are not applicable. Instead new methods are needed to describe such nonadiabatic processes, which are important for the study of solar cells, vision, and photosynthesis, among others. 1 One way to make large-scale simulations of these phenomena possible is to approximate the nuclear motion by an ensemble of independent trajectories that propagate under classical equations of motion. Among the simplest trajectory-based methods are Ehrenfest dynamics, in which the nuclei move on a mean-field potential defined by the instantaneous electronic populations, while the electronic variables follow exact subsystem dynamics according to the instantaneous nuclear configuration. This method has a number of known severe drawbacks, 2 but is still popular due to its simplicity and low computational cost. Other options of comparable cost include surface hopping 3 and mapping-based techniques. 4 In particular, the Meyer-Miller-Stock-Thoss (MMST) mapping 5,6 has recently regained attention. [7][8][9][10][11] As a generalization of the Schwinger bosonization to N -level systems, its basic principle is to represent the N electronic a) Electronic mail: johan.runeson@phys.chem.ethz.ch b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch states by N coupled harmonic oscillators that share a single excitation. This mapping is formally exact and has inspired a number of methods for calculating correlation functions, such as the linearized semiclassical initialvalue representation (LSC-IVR), 12 the Poisson-bracket mapping equation (PBME), 13,14 the symmetrical quasiclassical windowing approach (SQC), 7,15 partially linearized density matrix dynamics (PLDM), 16,17 and the forward-backward trajectory solution (FBTS) 18,19 of the quantum-classical Liouville equation. 20 These quasiclassical approaches all use a classical description of the nuclear dynamics, while preserving the exact quantum dynamics of an isolated subsystem.
Even though the MMST mapping is formally exact, its descendant methods are not, due to the quasiclassical approximation. In particular the classical dynamics may bring the system out of the singly-excited subspace. 14 One way to improve upon this is to introduce additional projectors. In principle one could do this at every time step, but in practice this is usually done only at the start and/or end of the simulation. 21 Another problem is that the zero-point energy of the fictitious harmonic oscillators is not respected by the classical dynamics. Historically it has been observed that this leakage can be mitigated by reducing the zero-point energy from 1 to a parametric value γ. 22,23 In the more recently introduced symmetrical quasiclassical windowing approach (SQC), 7 γ is determined via a window function, which is in turn freely chosen. In the case of two-level systems, there is a natural choice of γ that originates from the mapping of a spin vector, which was first proposed by Cotton and Miller 15 and was derived in our previous paper (paper I) 24 by mapping the two-level system to a spin- 1 2 instead of two harmonic oscillators. In the present paper we show that this spin mapping can be generalized to multiple levels. Its dynamics turns out to be equivalent to that of the MMST Hamiltonian, but with a new zero-point energy parameter γ, for which we derive a closed formula as a function of the number of levels.
The search for such a theory follows the intuition of Meyer and Miller, who originally considered the well-known equivalence between a two-level system and a spin- 1 2 system as an alternative derivation of their method. 25 This however turned out to be difficult to extend to many levels. Since their generalization no longer reduced to give the correct dynamics for an isolated subsystem, they abandoned this path in favour of the harmonic-oscillator mapping, which since then has inspired the rich field of mapping-based methods mentioned above. More recently, Cotton and Miller returned to the idea of a spin mapping by representing the twolevel problem in terms of two spins, in the hope of finding a more natural mapping than to harmonic oscillators. 26 Unfortunately, this approach did not reduce to the correct dynamics for isolated subsystems either, which has lead some authors to believe that spin is not a good classical analogue for a quantum system. 9 In the present paper we demonstrate how a spin mapping can indeed be generalized to multi-level systems, in a way that gives identical results to the Schrödinger equation for an isolated subsystem.
The main practical difference between our spin mapping and the MMST mapping lies in the definition of the phase-space distribution. While the 2N -dimensional phase space of MMST is unbounded, the spin-mapping phase space is confined to a sphere with 2N − 2 degrees of freedom. In this way, this phase space conserves the symmetries of the original quantum problem. The phase-space construction used in spin mapping was originally proposed by Stratonovich,27 and is now known as the Stratonovich-Weyl (SW) representation, which has found various applications in quantum optics. 28 It is a generalization of Weyl's correspondence rule 29 and the classical phase-space theories by Wigner and Moyal. 30,31 Early works of SW-representations for spin were made by Agarwal, Várilly and Gracia-Bondía. 32,33 These rely on the properties of the SU (2) Lie group, the fundamental symmetry of particles with spin. Brif and Mann have presented a construction for general Lie groups 34 and later Klimov and de Guise 35 as well as Tilma and Nemoto 36 for the case of SU (N ), which is the symmetry group of N -level systems. This has recently been used in the study of qudits (qubits generalized to multiple states). [37][38][39][40] In this paper we apply the Stratonovich-Weyl formalism to describe nonadiabatic dynamics in N -level molecular systems (but the resulting method is applicable for any quantum-classical problem). This leads to a straightforward generalization of our results for the two-level system. 24 The Stratonovich-Weyl representations could be formulated in spherical variables of coherent states, but like in the two-level case there is also a natural description in Cartesian variables, which leads to the same form of the Hamiltonian as in the MMST mapping, but with a more natural phase space that does not require projections and cannot suffer from unphysical leakage. In particular we derive a previously unknown closed for-mula for the zero-point energy parameter γ in terms of N .
The generalized spin mapping is not just a useful methodology in itself, but may also give insights about the standard MMST mapping. Recently it was found that the accuracy of MMST-based methods like LSC-IVR and PBME can be significantly improved by separating all observables into a linear combination of the identity operator and a traceless operator. 11,41 While it is not so obvious from the harmonic-oscillator picture why this would be a more natural choice, it is clear from the construction of the spin mapping that the identity must be treated separately. Therefore, the key to understand the success of traceless MMST might lie in the generalized spin mapping.
In Sec. II we use the generalized spin mapping to approximate correlation functions in a manner similar to classical Wigner dynamics. In Sec. III we apply the method to the seven-state Fenna-Matthews-Olson complex, which is a benchmark problem relevant for studies of light harvesting. The results can compete in accuracy with other state-of-the-art methods in the mapping community, and are far superior to conventional Ehrenfest dynamics with comparable cost.

Consider a molecular system with N electronic states and the general diabatic Hamiltonian
wherex andp are vectors of position and momentum operators of the nuclear modes with associated mass m, andV (x) is a Hermitian potential-energy matrix of shape N × N . We use the diabatic representation in this paper since it leads to the simplest formulation, but working in the adiabatic representation would also be possible. 42 Like in other trajectory-based methods, we will treat the nuclear variables classically (that is, replacex,p → x, p) but keep the quantum-mechanical evolution of the electronic operators. To handle the coupling between the two in a consistent fashion, we will map the electronic (subsystem) operators to a phase-space representation in which all variables are treated on the same footing. The mapping procedure will be similar to the spin mapping for two levels in paper I. 24 In each section we will therefore first remind the reader of the two-level case, before generalizing to N levels.
Throughout this paper we set = 1.

A. Generalization of the spin matrices
First we discuss the spin matrix decomposition of twolevel Hamiltonians, before we generalize to N levels.
Consider the Hamiltonian in Eq. (1) with a general (diabatic) potential matrix: It is well known that the Hamiltonian, or any other twolevel Hermitian operator, can be decomposed into a basis of spin operators and the identity: are the Pauli matrices multiplied by 1 2 . The explicit relations between the quantities in Eqs. (1) and (3) are: Without loss of generality we chooseV (x) to be real, so that H 2 = 0. Let us point out three important properties of the spin operators. First, they are traceless (i.e. tr[Ŝ i ] = 0) in contrast toÎ that has tr[Î] = 2, where lowercase tr denotes a trace over the subsystem degrees of freedom). As a consequence, the trace ofV (x) will only appear in H 0 , while H i≥1 only depends on the traceless part ofV (x). Note that this appears naturally and is not artificially imposed on the mapping, as is sometimes necessary for other mappings. 11,14,15 Second, the spin matrices are orthogonal: Other normalizations of the spin matrices are possible, but we shall keep the factor of 1/2 to maintain the connection to a spin system. The third relevant property of the spin matrices is that the sum of their squares is proportional to the identity: .
The reader probably recognizes the square-root of the proportionality constant, √ 3/2 = S(S + 1), as the magnitude of a classical spin vector for a spin S = 1/2. This observation will be important in the treatment of the N -level system.
Let us now generalize to an N -level potential. A general Hermitian (N × N )-matrix has N 2 independent elements, or N 2 − 1 for traceless matrices. Therefore the basis expansion can be written on the form whereŜ i are now (N × N )-matrices (also called the generators of the su(N ) Lie algebra). The matricesŜ i are necessarily traceless, and we keep the same normalization as in the two-level case, such that Eq. (5) is still fulfilled. Finally, it is well-known in the literature that the sum of the squares of the basis matrices is which is called the (quadratic) Casimir operator of su(N ). We include a short proof in Appendix A for completeness. The Casimir operator is invariant to unitary basis transformations, and therefore not dependent on the particular choice of decomposition in Eq. (7). This simple expression will be the key to defining the zeropoint energy parameter, which ultimately leads to significant improvements upon the MMST-mapping results.
There are many possible ways to choose the basis matrices, 43,44 but the theory of this paper will not depend on this choice. As an example for N = 3, a direct generalization of the Pauli matrices are the Gell-Mann matrices (which have been used in the SU (3)-symmetric theory of quarks 45 ): and the reader can easily confirm that iŜ . Note that the first three contain the two-level basis matrices as blocks padded with zeros. This construction can be generalized to higher N , and the details are given in Appendix B. In practice, we shall see that it is not necessary for the results of this paper to carry out the expansion in Eq. (7) at all, so that the basis does not have to be known explicitly.
One might ask what the basis matrices have to do with spins once N > 2. Already in the 1970s, Meyer and Miller proposed a mapping of the N -level system to a higher spin S = 1 2 (N − 1). 25 While their spin-matrix decomposition is the same as ours for N = 2, it is different for all N > 2, since in their construction not all S i are traceless (such that there is no Casimir invariant as in Eq. (8)). Their basis matrices are therefore not generators of su(N ) and cannot be used to derive the results of this paper. Nonetheless, this can be easily fixed such that it is possible to construct the su(N ) generators from a spin picture. For the interested reader we show in Appendix C how to obtain basis matrices for N = 3 by describing a spin-1 system as two interacting spin-1 2 particles in a triplet configuration.
An even more important difference between the approach introduced in this paper and that of Ref. 25 is the phase-space representation used to convert the spin matrices to classical variables. Meyer and Miller mapped each matrix to the same two variables for any number of levels, which again is equivalent to our work for the N = 2 case, 46 but not for N > 2. Although a two-variable phase-space is appropriate for the SU (2) symmetry of a single spin-1 particle, the true symmetry group of the three-level system is SU (3), which we represent by four phase-space variables, as explained in Sec. II B.
Thoss and Stock have also investigated a spin-1 2 mapping of two-level systems, 47 and derived a semiclassical initial-value representation of its corresponding propagator. Like our approach, they also use spin coherent states, and their dynamics is exact for an isolated subsystem. However, they did not generalize their method to more than two levels. In this paper we pursue a quasiclassical approach to such a generalization by using the Stratonovich-Weyl representation of the N -level problem.

B. Stratonovich-Weyl representations
Again, we will start with the two-level case that was previously presented in paper I. 24 As is commonly done in textbooks, one can think of the diabatic states |1 and |2 as the eigenstates of a (fictitious) spin-1 2 degree of freedom. In order to map these to a phase-space, we introduce the spin-1 2 coherent states 48 where u denotes a unit vector with spherical coordinates (θ, ϕ) and the states are normalized such that u|u = 1.
The expectation values of the spin operators in this state have the simple form It is then clear that { u|Ŝ i |u } 3 i=1 are the Cartesian coordinates of a sphere with radius 1/2. An even more important observation is that u|Ŝ i |u are orthogonal functions on the sphere: where we have defined the integration measure as du = 1 2π sin θ dθ dϕ. We refer to u|Ŝ i |u as the Q-representation (or Q-function) of the operatorŜ i .
Likewise the Qrepresentation of a general operator is defined as It is easy to show that this is equivalently written as whereŵ Q is the Stratonovich-Weyl kernel of the Qrepresentation. 28 The Q-representation is analogous to the Husimi representation in the nuclear variables. 49 It is dual to the P-representation, which is analogous to the Glauber-Sudarshan representation in the nuclear variables. 49 What is meant by 'dual' is that any quantum-mechanical trace of a product of operators can be expressed as an integral over a product of Q-and P-symbols as (15) Most importantly, there is also a W-representation, The W-representation is analogous to the Wigner representation in the nuclear variables.
In particular [Î] W (u) = 1 and [Ŝ i ] W (u) = √ 3 u|Ŝ i |u , so that for the general operator in Eq. (13) we have The introduction of the Stratonovich-Weyl representations Q, P and W for the spin operators constitutes the major difference between our spin-mapping approach of paper I 24 and the spin-mapping models of Refs. 9, 25, 26, and 47. We will employ the Stratonovich-Weyl representation in a similar way for the N -level system.
For the N -level case, the coherent state will have a more complicated form than in Eq. (10), but for the following treatment it is not necessary to write out its explicit expression. We shall denote the generalized coherent state by |Ω , parametrized by 2N − 2 angles, and leave the details of its construction to Appendix D. (One can in the following always replace |Ω by |u to recover the results of the two-level case.) For now all we need to know is that Ω|Ŝ i |Ω are orthogonal functions such that which is proved in Appendix E. It turns out that the Stratonovich-Weyl kernels are remarkably simple to generalize for N levels, as has been shown by Tilma and Nemoto for general SU (N )symmetric coherent states. 36 With the choice of normalization in Eq. (5), the kernels arê giving the SW-representations [Â] s (Ω) = tr[Âŵ s (Ω)] for s ∈ {Q,P,W}. The readers can easily convince themselves that with this construction, traces of products still obey Eqs. (15) and (16) (but with |Ω instead of |u ), as a consequence of Eq. (18).
In the two-level case in paper I, 24 we interpreted the W-functions of the spin operators as the components of a classical spin vector with magnitude √ 3/2. Let us define a generalized magnitude as the square-root of where the sum is worked out in Appendix E and relies on the fact that theŜ i operators are traceless. Thus the Casimir invariant in Eq. (8) is the generalization of the (squared) spin magnitude, and it is preserved by the mapping.
Many authors have described Ω|Ŝ i |Ω as the components of a generalized Bloch vector in N 2 − 1 dimensions. 44,50,51 The impact of the W-representation would then be to scale this vector to the length N 2 −1 2N . However, note that not all of its N 2 − 1 components can be independent for N > 2, since |Ω depends only on 2N − 2 spherical variables. Consequently, only a subset of points on such an (N 2 − 1)-sphere correspond to physical states for N > 2. [52][53][54] We shall therefore not pursue that picture in this paper.
A more natural picture would be to think of the 2N −2 spherical variables {θ n , ϕ n } N −1 n=1 of |Ω as the orientations of N − 1 spin- 1 2 vectors (their interpretation as spins is explained in Appendix C). Fig. 1 shows an example for N = 3. In the two-level case, we saw in paper I 24 that the orientations corresponding to single basis states were found at the poles in the Q-representation, but at "polar circles" with fixed θ in the W-representation. In Appendix D it is worked out that the W-representation of the basis states can be represented by polar circles also for N > 2, but on different latitudes from the N = 2 case. This is in contrast to the Q-representation, where these circles would be replaced by points at the respective poles.
While this picture is instructive as a generalization to Fig. 2 in paper I, 24 it does not help us find the equations of motion of the system. To describe the dynamics in a simple way, we shall therefore from now on switch to a Cartesian representation. This will also reveal the link between the generalized spin mapping and the MMST mapping.

C. Dynamics in Cartesian variables
An alternative to using spherical variables is to write the coherent states in terms of complex coefficients {c n }: Given the constraint n |c n | 2 = 1 and an arbitrary choice of global phase, the {c n } have 2N − 2 real degrees of freedom. Starting as usual with the two-level coherent state, |u = c 1 |1 + c 2 |2 , the orthogonal functions in Eq. (11) take the form Let us insert these into the W-representation of an arbitrary operator in Eq. (17): We now introduce the Cartesian variables X n and P n via 3 1/4 c n ≡ 1 √ 2 (X n + iP n ). Because of n |c n | 2 = 1, these are constrained to a sphere with squared radius The general W-representation becomes for which a particularly important case is the two-level Hamiltonian (which we choose to be real), With the use of Eq. (24), this can also be written where γ = √ 3 − 1 is called the zero-point energy parameter. (Note that some authors use an alternative convention used for γ, which is half the value of ours.) This is similar to the Hamiltonian derived by Meyer, Miller, Stock and Thoss, 5,6 the difference being that γ = 1 in their formulation. This value emerged from a Langer correction in Meyer and Miller's formulation, and from the commutation relations of harmonic-oscillator operators in the formulation by Stock and Thoss. Stock and Müller have however observed that decreasing the value of γ often gives more accurate results, 22,23 and suggested choosing γ ≈ 1 2 as a general rule of thumb. 4,55 The particular value γ = √ 3−1 ≈ 0.732 has previously been proposed by Cotton and Miller 15 and appears naturally from the Stratonovich-Weyl formalism for a two-level system mapped to a spin-1 2 . 24 We will now derive the value of γ from the Stratonovich-Weyl representation for a general N -level system.
Again we focus our attention on the W-representation and simply state the Q-and P-versions at the end of the section. By construction, it is clear that Given the kernel in Eq. (19c), the W-representation ofŜ i is (29) Let us therefore introduce X n and P n through from which follows that X n and P n are constrained to a hypersphere with squared radius Note that when N = 2, this result reduces to Eq. (24). To find the value of γ in the Hamiltonian, consider its decomposition in basis matrices according to Eq. (7) with By inserting Eq. (30), we finally recover the MMST-form of the Hamiltonian (which is again chosen to be real): This is always smaller than the value γ = 1 that Stock and Thoss obtained from the commutation relations [X n ,P n ] = 1 of the harmonic oscillator operatorsX n ,P n of the extended mapping space. Note that we have not invoked such an extended space, but instead derived Eq. (34) purely based on the commutation relations of the original problem that enter through the Casimir invariant. The new value of γ as well as the constraint in Eq. (31) are both minor modifications to any code that uses MMST mapping or Ehrenfest dynamics, but lead to significant improvements in accuracy, as we shall show in Sec. III. The derivation holds under the assumption that the symmetry group of the subsystem is SU (N ). This means for example that we assume that no level is decoupled from all the others, since that would lead to a reduction of the symmetry group and change the value of γ. A treatment in terms of time-dependent symmetry groups would be much more involved 56 and is outside the scope of this paper. In the case of N = 1 there is only one basis matrix (the identity of shape 1×1) so that H W = H 0 and γ does not appear, which recovers the standard singlesurface Born-Oppenheimer Hamiltonian.
The same analysis can be done for the Q-and the Pfunctions, and the results for R 2 and γ are summarized in Table I. Note that γ is independent of N for the Qand the P-representations, but decreases with N in the W-representation, and that all quantities coincide with our previous results 24 for N = 2.
What we have presented is a major step forward from the previously suggested spin-mapping approaches, 9,25,26 which did not recover the MMST Hamiltonian and did not reduce to the exact quantum dynamics for an isolated subsystem. One could even say that the generalized spin mapping is a more natural derivation of the MMST Hamiltonian, as it requires no extended phase space and directly gives a γ closer to those found optimal in numerical simulations. In the rest of this section we will give two further reasons for this point of view.
The equations of motion that correspond to the Hamiltonian in Eq. (33) look identical to those of the MMSTmapping: but in contrast to the MMST mapping, our X n and P n are constrained to a hypersphere with squared radius R 2 (see Table I) and so the ensemble dynamics will be subtly different. Each trajectory exactly preserves R 2 , so that there is no leakage out of the mapping space, thereby solving a problem of the original MMST mapping. Like in the MMST dynamics, the classical equations of motion for the X n and P n variables exactly correspond to the electronic Schrödinger equation for the uncoupled system.
In literature it is common to separateV (x) into a stateindependent and a state-dependent part, and the particular choice of splitting can influence the results of some methods (such as LSC-IVR and SQC). Typically authors suggest to separate the traced and traceless parts 11,14,18 although other choices have also been used. 57,58 In our approach, the weights of the potential-energy surfaces always sum up to one: which means that it is independent of such a splitting.

D. Correlation functions
We will now use the results of Secs. II B-II C to approximate the correlation function where capitalized Tr means a trace over both electronic and nuclear states. The trace over the electronic degrees of freedom can be written as integrals of Stratonovich-Weyl functions (see Eqs. (15) and (16)). Likewise, the trace over the nuclei can be expressed in terms of the Wigner distribution and we choose initial conditions such thatρ =ρ nuc ⊗Î.
(The initial electronic state is defined byÂ.) Then the correlation function can be exactly written as where s ∈ {Q, P, W},s is the dual of s, and · · · = dx dp dX dP · · · δ(X 2 + P 2 − R 2 s )ρ nuc (x, p) dx dp dX dP δ( .
The additional factor of N appears because we have defined Tr[ρ] = N but [ρ] s = 1, the subscript s on R s specifies which radius to use from Table I, and we used the shorthand notation X 2 ≡ N n=1 X 2 n and similarly for P 2 .
We now propose to approximate the correlation function by where the dynamics is driven by the Hamiltonian H s in the s-representation. This formula is the multi-level generalization of the quasiclassical spin-mapping method in paper I. 24 Note that for s = W, we are using classical Wigner dynamics in both the nuclear and the electronic degrees of freedom. In paper I, 24 we showed how the dynamics follows from approximating the time derivative of B with a Poisson bracket, similar to the approximation used in PBME. 14 Typically we will be interested in population transfer from a state n to a state m, i.e.Â = |n n| andB = |m m|. The corresponding Stratonovich-Weyl functions are the population observables where we use subscripts on γ and R 2 to distinguish between the s-and thes-symbols. The factor R 2 s /R 2 s appears when s = Q or P, because the X n and P n variables are sampled from a hypersphere with radius R s but measured on a sphere with radius Rs. In the symmetric case of s =s = W, Eqs. (41a) and (41b) reduce to the same expression. It is however no more difficult to calculate off-diagonal elements of the density matrix, for example: and thes symbols are again obtained by multiplying with R 2 s /R 2 s . We thus have the alternative of calculating the correlation function in a symmetric way, meaning (s,s)=(W,W), or in an asymmetric way, that is (s,s)=(Q,P) or (P,Q). In paper I 24 we saw that (Q,P) and (W,W) both gave accurate results for a wide range of spin-boson models, while (P,Q) was always less accurate. After running tests on further systems we have observed that (Q,P) is not always so reliable but that the symmetric definition (W,W) is the most robust. This is confirmed by the results we show in Sec. III and in the Supplementary Material.
The initial distribution of the mapping variables that follows from the Stratonovich-Weyl formalism is a uniform distribution over the sphere X 2 + P 2 = R 2 s . We will call this full-sphere initial conditions. Note that it gives the results of all n → m transitions in a single simulation. The sampling of the distribution δ(X 2 + P 2 − R 2 s ) is easy to implement in practice by drawing {X n , P n } N n=1 from a standard normal distribution and rescaling them with a common factor so that X 2 + P 2 = R 2 s . Previous authors in the mapping community have also used an approximation called focused initial conditions 13,19,22,23,[59][60][61] in which the initial distribution only includes points that directly correspond to coherent states that diagonalizeÂ, rather than by weighting as in Eq. (40). In the case ofÂ = |n n|, this means that 1 2 (X 2 n + P 2 n − γ) = 1 while 1 2 (X 2 k + P 2 k − γ) = 0 for all k = n, or equivalently with r k=n = √ 2 + γ, r k =n = √ γ and uniformly sampled φ k ∈ [0, 2π). This can be interpreted as sampling from the "polar circles" in Fig. 1. In the correlation function, A(X, P ) = 1 by construction, so that where ... foc = dx dp dX dP · · · ρ foc (X, P )ρ nuc (x, p) dx dp dX dP ρ foc (X, P )ρ nuc (x, p) , uses the focused distribution ρ foc (X, P ) = δ(X 2 n +P 2 n −γ −1) and trajectories are defined according to H s . Note that for focused methods the observable B s (X(t), P (t)) must be calculated with the same index s as the Hamiltonian, in contrast with methods that use full-sphere initial conditions. 62 It is also possible to define focused initial conditions when starting from off-diagonal elements of the density matrix, as explained in paper I. 24 This prescription is the same as that of Müller and Stock, 23 apart from the value of γ in Eq. (34), which we derived instead of treating it as a free parameter. Our γ decreases with N , but the convergence to zero is so slow that only with N ≥ 360 does it reach γ ≤ 0.1. It is interesting to note that the limit N → 0 gives the standard MMST-value γ = 1, and that the spin-mapping value is therefore somewhere between that of MMST and the Ehrenfest value γ = 0. A comparison of numerical values for γ is given in Table II. One might ask if the N → ∞ limit of focused W would be equivalent to the Ehrenfest method, which is the same as focused Q. Bearing in mind that R W diverges as N → ∞, while γ → 0, this is still an open question.
Finally note that if one applies focused initial conditions to a problem without inter-state couplings, the prescription reduces to Born-Oppenheimer dynamics on the initial state.

III. APPLICATION TO THE FENNA-MATTHEWS-OLSON MODEL
We have tested the theory of Sec. II on a seven-level model of the Fenna-Matthews-Olson (FMO) complex, which is a well-studied light-harvesting pigment-protein complex found in green sulphur bacteria. 65 Each diabatic state represents an exciton localized on one of the sites. This is a challenging benchmark problem for electronically nonadiabatic dynamics and allows our method to be compared with other mapping approaches, 16,41,64,66 as well as to numerically exact results obtained via the hierarchical equation of motion (HEOM) approach. [67][68][69][70][71]

A. Model description
The model Hamiltonian is of the subsystem-bath type: where the subsystem Hamiltonian in a diabatic basis is given in units of cm −1 as 72 (47) Each level is coupled to its own bath of F nuclear modes with unit mass and frequencies ω j , so that the total bath Hamiltonian iŝ The baths of the different sites have identical frequencies and are not directly coupled to each other. The systembath coupling is in turn with coupling coefficients ξ j . The frequencies are distributed according to a Debye spectral density where ω c is the characteristic frequency of the bath (τ c = ω −1 c is its corresponding timescale) and λ is the reorganization energy. In accordance with previous work, 41,64,68-70 we used λ = 35 cm −1 in all our simulations. We used a discretization of the bath with F = 60 modes per site (in total 420 modes), according to the discretization scheme in Ref. 73. The discretization also determines the coupling coefficients ξ j .
To initialize the nuclear bath, we sampled the Wigner distribution in Eq. (38), which is explicitly where α j = tanh βωj 2 . Each simulation was run with timestep 1 fs and 10 6 trajectories to guarantee convergence, although we point out that it was possible to observe the trend of each line already with 10 3 trajectories.

B. Population dynamics
We have tested our theory on the FMO model using the same parameters as in the original paper by Ishizaki and Fleming. 68 In each figure we will compare six methods, where upper panels correspond to full-sphere initial conditions and lower panels to focused initial conditions (these are defined in Sec. II D). Within each row γ increases from left to right, so that left panels display methods with γ = 0 (the lowest possible), middle panels the W-value of γ from Eq. (34), and right panels the large-γ case. In the upper right panel this is the P-value, while in the lower right panel we show the standard MMSTvalue γ = 1 for the focused method (the P-value would typically be worse). Note that the focused method with γ = 0 (lower left panels) is identical to Ehrenfest dynamics. Fig. 2 shows the results for low temperature (T = 77 K) and a fast bath (τ c = 50 fs), which is the hardest of the model problems since it has the strongest quantum effects. It is clear from each row that dynamics using γ derived in the W-representation is generally more accurate than the other cases. For W, the difference between using full-sphere or focused initial conditions is negligible. Typically the focused methods converge with an order of magnitude fewer trajectories than the full-sphere methods.
These observations become even clearer when looking at the long-time limit of the same model in Fig. 3. Again the W-value of γ derived in Eq. (34) is the most accurate for the final populations, while Ehrenfest and focused MMST are very unreliable. The W-methods may still predict unphysical negative populations, but the absolute error is still typically smaller than in the other methods (and note that such negative populations are also possible in other mapping approaches 12,13,16,18,41 ).
Finally we show the results for a higher temperature (T = 300 K) in Fig. 4. This problem is not as hard as the previous, in the sense that all methods have decent accuracy, but it is still clear that W is the most accurate. In the Supplementary Material we show the case of starting from state 6 instead of 1, as well as the case of a slow bath (τ c = 166 fs), and they give further weight to our conclusions.
It should be noted that the middle column (W) results shown here clearly outperform both PBME, 75 LSC-IVR 76 and Ehrenfest dynamics. The W results are of similar accuracy to other state-of-the-art mapping approaches such as SQC 64 and traceless MMST. 41 They are slightly less accurate than PLDM for this system, 16,74 but the improvement of PLDM compared to linearized MMST suggests that a (future) partially linearized version of the spin-mapping method might perform even better. As discussed in Sec. II C, the W approach has the advantages compared to MMST-based methods that it does not require choosing a window function, it has no leakage from the mapping space, and it is independent of the splitting of the potential-energy matrix.

C. Bipartite entanglement
We now turn to the problem of calculating off-diagonal elements of the density matrix. A quantity for which numerically exact benchmarks exist is the bipartite entanglement (sometimes called concurrence) between state n and m, defined as 2|ρ nm (t)|, after an initial excitation to one of the states. Here ρ nm (t) denotes elements of the reduced density matrix of the subsystem. More specifically, we compute the correlation functions C AS + and C AS − withÂ = |k k| being the initial state andŜ + = |n m|+|m n| for 2 Re[ρ nm ] andŜ − = i(|n m|−|m n|) for 2 Im[ρ nm ]. The time-dependent concurrence is then given by In Fig. 5 we show the concurrences that are largest in magnitude for an FMO model with τ c = 100 fs. As before, the middle panels are the most accurate, while both Ehrenfest and focused MMST deviate significantly from the benchmark for at least one of the concurrences. All of the methods begin to disagree with the benchmark after about 0.2 ps, but the error for W is smaller than what has previously been reported with PBME. 75 It is noteworthy how the qualitative shapes of all lines can be predicted by our quasiclassical method, which cannot be done with Redfield theory (since that requires λ to be much smaller than the electronic couplings). 77 In the Supplementary Material we also show the long-time limit of Fig. 5, the equivalent calculations starting from state 6, as well as the higher temperature case, which illustrate the same trends as have already been pointed out.
All together, the symmetric W approach that we propose in this paper is seen to be a promising method that is both simple to compute and resolves several drawbacks of traditional MMST methods.

IV. CONCLUSIONS
In this paper we have generalized the spin mapping of a two-level system in paper I 24 to N -level systems, which is a problem that had not been satisfactorily solved since is was first posed by the seminal works of Meyer and Miller in 1979. 25 The general idea is to make the classical phase space inherit the SU (N )-symmetry properties of the quantum system. Of particular significance is the Casimir invariant, which plays the role of a generalized (squared) spin magnitude. This quantity is independent of basis representation and controls the overall strength of the nuclear forces. In contrast with previous spin mapping attempts, 9,25,26 we have shown how the dynamics can be generated by a quadratic Hamiltonian of the same form as in the standard harmonic-oscillator mapping, but with a new formula for the zero-point energy parameter γ. Originally γ was included as a Langer correction, then justified through the commutation relations of a set of harmonic-oscillator operators. Now we recommend that this term is changed to become dependent on N , with values close to what was previously found optimal when it was treated as a free parameter. 23 One can therefore say that the generalized spin mapping is a more natural derivation of the MMST Hamiltonian than the original harmonic-oscillator mapping. We have shown that the spin mapping solves the problem of leakage from the physical space, so that there is no need for additional projectors. The present theory also does not assume any particular form of the Hamiltonian other than that the subsystem belongs to the symmetry group SU (N ), while MMST-based approaches often depend on howV (x) is split in a state-dependent and a state-independent part.
We have demonstrated that the resulting method can predict population dynamics in benchmark systems to similar accuracy as other state-of-the-art mapping approaches such as SQC and traceless MMST. In the future we expect that the accuracy can be extended to longer times by combining the dynamics with a general-ized quantum master equation, as has been successfully done for other mapping approaches. [78][79][80][81] Another natural extension would be to develop an FBTS or PLDM method based on spin mapping. Finally, we also believe that the spin mapping will be relevant in the search for a nonadiabatic extension to ring-polymer molecular dynamics. 58,[82][83][84][85]

SUPPLEMENTARY MATERIAL
Supplementary Material is available with additional results for the FMO model. Here we give a proof for the formula for the Casimir invariant that is known from many textbooks (for example p. 500 in Ref. 86). The generators of the Lie algebra su(N ) have commutation relations of the form where f ijk is totally antisymmetric and contains the structure constants of su(N ). Now define the quadratic Casimir operatorĈ 2 ≡ i . It is easy to show that it commutes with all generators: and the proportionality constant is easily found as Other bases can be constructed from linear combinations of these. However, as previously mentioned, it is not actually necessary to choose a particular basis in order to obtain the results of this paper.

Appendix C: Formulation in spin-1 matrices
It is well known that an N -level system can also be described in terms of operators of a spin S = 1 2 (N − 1) system. 87 Here we show the specific example of how a three-level system can be related to a spin-1 particle, in a slightly different way than Meyer and Miller in Ref. 25. One way to represent a spin-1 is via symmetrized product states of two spin-1 2 particles. A natural basis of this (triplet) space is {|↑↑ , 1 √ 2 (|↑↓ + |↓↑ ), |↓↓ }. In this basis, the total spin projection along each of the coordinate axes areŜ which we can take as the first three basis matrices. Their phase-space functions are analogous to p-orbitals. As the remaining five basis matrices, we take the following dorbital analogues: which are clearly linear combinations of the Gell-Mann matrices in Eq. (9), but with a physical meaning in terms of spin-1. These matrices are subtly different from those of Meyer and Miller in Ref. 25 which definedŜ 8 asŜ 2 z . This has a non-zero trace and therefore does not comply with our requirements.
In these variables the W-representation of the singlestate projectors are Note that their sum is one for all angles, meaning that the total population is identically one. To find when the system is entirely in state n, we solve the system of equations [|n n|] W = 1 and [|k k|] W = 0 for k = n.
The solutions define circles with fixed θ k for k ≤ n (and k ≤ N − 1), see Fig. 1. Lastly we point out that even though we visualize the spherical coordinates using multiple "spins", the phase space constructed from SU (N )-symmetric coherent states is different from when mapping multiple spins independently to a classical phase space (as in Ref. 93).

Appendix E: Auxiliary formulas
To prove the preservation of the Casimir invariant in Eq. (20), we need some further properties of the coherent states. It is known from the theory of harmonic functions that | Ω|Ω ′ | 2 can be expanded as 34 where τ ν are constants and Y ν (Ω) = 2 τν 1/2 Ω|Ŝ ν |Ω are generalized spherical harmonics that fulfil dΩ Y ν (Ω)Y ν ′ (Ω) = δ νν ′ . The index ν runs from 0 to N 2 − 1, whereŜ 0 = 1 √ 2NÎ . One can show that τ ν is invariant to transformations within an irreducible subspace, 34 which for our purposes means that are all equal (to, say, τ 1 ), and only τ 0 is different. 94 Let us insert Y ν (Ω) into Eq. (E1) and set Ω ′ = Ω: so that which is used in Eq. (20). Further, the normalization of Y i (Ω) gives Summation over all i and j leads to  Fig. 3 in the paper (T = 77 K, τc = 50 fs) but starting from state 6. Solid lines show numerically exact HEOM results. 5 The system is the same as in Fig. 6