Improved population operators for multi-state nonadiabatic dynamics with the mixed quantum-classical mapping approach

The mapping approach addresses the mismatch between the continuous nuclear phase space and discrete electronic states by creating an extended, fully continuous phase space using a set of harmonic oscillators to encode the populations and coherences of the electronic states. Existing quasiclassical dynamics methods based on mapping, such as the linearised semiclassical initial value representation (LSC-IVR) and Poisson bracket mapping equation (PBME) approaches, have been shown to fail in predicting the correct relaxation of electronic-state populations following an initial excitation. Here we generalise our recently published modification to the standard quasiclassical approximation for simulating quantum correlation functions. We show that the electronic-state population operator in any system can be exactly rewritten as a sum of a traceless operator and the identity operator. We show that by treating the latter at a quantum level instead of using the mapping approach, the accuracy of traditional quasiclassical dynamics methods can be drastically improved, without changes to their underlying equations of motion. We demonstrate this approach for the seven-state Frenkel-Exciton model of the Fenna-Matthews-Olson light harvesting complex, showing that our modification significantly improves the accuracy of traditional mapping approaches when compared to numerically exact quantum results.


Introduction
Simulating nonadiabatic effects in quantum dynamics continues to pose a considerable challenge in theoretical chemistry and physics, especially in the condensed phase. Arising when the energies of two or more electronic states approach each other, resulting in the breakdown of the Born-Oppenheimer approximation, these effects have been found to have a profound impact on a wide range of systems spanning physics, chemistry and biology. [1][2][3] The development of simulation methods for nonadiabatic effects has thus continued to be the focus of considerable research efforts. Methods relying on an explicit expansion and propagation of the wavefunction, oen on a grid, have yielded highly accurate results. [4][5][6] However, many models inspired by condensedphase systems still prove too computationally expensive to treat with these methods, due to their unfavourable exponential scaling with system size. Despite recent efforts to overcome this scaling hurdle, 7,8 many systems from the elds of chemistry and biology, especially those in a condensed-phase environment, are simply too large to treat using a wavefunction-based approach. Mixed quantumclassical methods,  though inherently more approximate, are oen the only choice when seeking to simulate nonadiabatic dynamics in the condensed phase. As many of these scale linearly with system size, they can readily be applied to large and complex realistic systems, yielding highly valuable insights at reasonable computational costs.
The representation typically chosen for a nonadiabatic process consists of a continuous nuclear phase space and a set of discrete electronic states. The resulting Hamiltonian is given bŷ wherep j and m j are the momentum operator and mass of nuclear degree of freedom (DoF) j respectively, andx is a vector of length F consisting of the position operators for each nuclear DoF. U(x) is the state-independent potential, and the state-dependent potential is given bŷ where S is the number of electronic states. The diagonal elements ofV (x) are the diabatic potential energy surfaces, while its off-diagonal elements are the couplings between the electronic states. Everything that follows does not rely on a particular choice of potential, i.e. we are not limited to simple harmonic models. Note that we will use reduced units throughout, such that ħ ¼ 1.
The mismatch between the continuous nuclear phase space and the discrete space of electronic states constitutes a recurring challenge in mixed quantumclassical dynamics. The mapping approach solves this problem by projecting the electronic degrees of freedom into a space of singly-excited harmonic oscillators (SEOs). 12,30,31 In the space of the SEO wavefunctions, the representation of state |ni is given by where X and its conjugate P are vectors of length S, corresponding to the position and momenta of the SEOs. ðX n X m þ P n P m À d nm ÞV nm ðxÞ: In addition to the relative simplicity of the mapping approach, the extended phase space grows linearly with the number of electronic states. Furthermore, given the favourable scaling of classical trajectories with respect to the nuclear DoFs, a number of mixed quantum-classical dynamics approaches, aimed specically at large, realistic systems in the condensed phase, have been developed based on this formalism. [10][11][12]14,24,26 Note that in this work we will use the term quasiclassical to refer to mixed quantum-classical approaches which employ a single set of mapping variables per electronic state as well as a single set of positions and momenta for each nuclear degree of freedom.
Quasiclassical methods yield accurate results for most observables at short times. In the long time limit however, they are well known to degrade in accuracy, especially for the relaxation to thermal equilibrium following an initial electronic excitation. Attempts to address this shortcoming with the use of master equations have shown considerable promise. 32,33 Other approaches to improve quasiclassical dynamics have led to the development of related dynamics approaches. For instance, the symmetrical quasiclassical windowing method uses a windowing function to "bin" the electronic populations, ensuring that they have integer values at the beginning and end of each trajectory. This approach has been applied to the benchmark we study below, achieving accuracy comparable to that reported here. 34-38 A number of methods which depart from the equations of motion underlying quasiclassical dynamics, but remain close to its overall motivation, have also shown considerable promise in treating multi-state systems. Prominent examples, which have been very successfully applied to the benchmark studied here, include the forward-backward trajectory solution 19,20 (FBTS) and the partially linearised density matrix (PLDM) method. 16,17 In a recent publication, we have however shown that a simple modication, with a similar motivation as that underlying the use of master equations, can drastically improve the performance of quasiclassical methods, without changing the equations of motion. 39 We split the population operator into two parts, one of which is the identity. 18 We can then use our understanding of the exact behaviour of this operator to drastically improve traditional quasiclassical methods. The resulting approach has the benet of retaining all the advantages of these methods, as the underlying equations of motion are unchanged.
Here we extend this approach, which was originally presented for only two electronic states, to an arbitrarily large electronic space. We use the fact that the electronic population, like any Hermitian operator, can be expanded exactly into the identity and a purely traceless component. 18 Given that the behaviour of its quantum operator is well understood, we treat the identity exactly, resulting in a simpler phase-space representation of the population operator, involving only the traceless part. Population dynamics calculated using these modied operators are of drastically higher quality than those obtained from the traditional quasiclassical denition.
We apply this general formulation to the challenging benchmark model for the Fenna-Matthews-Olson (FMO) light harvesting complex. [40][41][42][43][44] Our results are signicantly more accurate than those obtained using the standard operator denitions and in excellent agreement with numerically exact quantum dynamics methods.

Theory
Here we extend our previous work 39 by presenting a general formalism, which can be applied to any system of multiple electronic states.

Quasiclassical population operators
In the mapping formalism, the operator |nihn|, which measures the population of electronic state |ni can be written as where X n and P n are the mapping variables associated with state |ni. In the quasiclassical approximation, the Wigner transform is used to dene a phasespace representation for the operators of interest. The Wigner transform of a general operator,Ô, is given by When considering the population operator, there are two representations one can choose to Wigner transform, corresponding to either the le-or right-hand side of eqn (5). Using the le-hand side is equivalent to including a projection on the SEO subspace, which yields an expression in terms of the harmonic oscillator wavefunctions as in eqn (3). 18 The resulting phase-space representations, identied as A w n and A SEO n respectively, are Note that crucially, A w n s A SEO n Â f. Each of these phase-space representations is derived via the Wigner transform of a formally exact mapping form of the |nihn| operator. Therefore, a clear choice of which to use when calculating observables is not obvious a priori. A more detailed discussion of the possible combinations of phase-space representations and electronic initial conditions can be found in our recent work. 39 An observable commonly computed using quasiclassical methods is the population of a given electronic state, |ni, given that the system was initially in a pure state, |mi. In quantum mechanics this is dened by P n)m ðtÞ ¼ Tr where r b is a density matrix which denes the initial state of the nuclei, normalised such that the trace over nuclear DoFs only is Tr b [r b ] ¼ 1.

Traceless projection operators
There are two differences between the phase-space representations given in eqn (7a) and (7b): the factor of f(X,P), which is only present in A SEO n , and the differing constant terms, which are related to zero-point energy (ZPE) of the mapping DoFs. 45 The origin of the latter is that both the projected and unprojected forms of |nihn| have a non-zero trace. We propose a form of the quantum population operator in which the trace is shied to the identity operator, which in turn is treated exactly using quantum mechanics. 39 The result is a phase-space representation of the quantum population operator which is traceless.
There is a unique expansion of the population operator |nihn|, such that: |mihm| is the identity operator, andQ n is, by design, traceless, Note that in a two-level system, this operator is the Pauli spin matrix, i.e.Q 1 ¼ ŝ z , such that |1ih1|¼(Î+ ŝ z )/2, which was used in our previous work. 39 Substituting this denition for the quantum population operator into eqn (9) and expanding yields where we have used Tr[r bQm ] ¼ 0 and Tr[r b Î] ¼ S. The nal two terms in this expression are quantum correlation functions which can be approximated by well-known quasiclassical dynamics methods. Following the standard quasiclassical procedure, in order to calculate the value of the population operator given in eqn (12), we Wigner transform the operators in these two constituent correlation functions. The phase-space representation of the traceless operatorQ n is If we had performed the Wigner transform on the projected operator, the phase-space representation would simply be Q n (X,P)f(X,P). Note that neither expression contains constant terms which play the role of ZPE-parameters.
It would be possible to arrive at a phase-space representation of the identity operator via similar Wigner transforms. However, we suggested in our previous work 39 that we can instead use our understanding of its behaviour in quantum mechanics, which is to leave its operand unchanged. We therefore simply avoid directly computing the identity altogether.
Starting from the exact expression forP n)m ðtÞ in eqn (12), we thus arrive at our nal quasiclassical expression for the population of electronic state |ni, assuming the system was initially in state |mi, The constituent correlation functions C IQ n and C Q m Q n are given by where formed density matrix of the nuclear DOFs. In practice the values of these correlation functions are averaged over an ensemble of trajectories, with initial conditions for the mapping variables being drawn from either f(X,P) or f 2 (X,P), depending on whether the projected forms of one or both ofQ m andQ n were Wigner transformed. This corresponds to a ¼ 1 and a ¼ 2 respectively. Note that we can include the factors of f(X,P) at time zero, because this function is constant over the course of any trajectory evolving under the Hamiltonian H; given in eqn (4). Also the two constituent correlation functions can be calculated for all values of m and n in a single simulation. Just as in traditional quasiclassical methods, 11,14 the values of these constituent correlation functions, and therefore P n)m ðtÞ; are exact in the limit of t ¼ 0.

Traditional quasiclassical dynamics methods
The traditional quasiclassical approach does not involve treating the identity quantum mechanically as we have done above. There are two standard approaches which differ in whether both population operators are projected onto the subspace, or just one. These methods were derived in different ways 11,14 and are called the Poisson bracket mapping equation 14,18 (PBME) and the linearised semiclassical initial value representation 10,11 (LSC-IVR) methods. LSC-IVR commonly involves projecting both operators prior to Wigner transforming them, i.e. using |mihm|1A SEO m ðX; PÞ and |nihn|1A SEO n ðX; PÞ: The Wigner transform of each operator yields, as per eqn (7b), a factor of f(X,P). Initial conditions for the mapping variables are therefore sampled from f 2 (X,P). In PBME on the other hand, traditionally only the operator for the initial population is Wigner transformed in its projected form. The operators are therefore |mihm|1A SEO m ðX; PÞ and |nihn|1A w n ðX; PÞ: Only the transform of |mihm| yields a factor of f(X,P). Consequently, electronic initial conditions are sampled from f(X,P). Using these denitions, the electronic population can be calculated from We note that the differences between these two methods do not actually stem from the derivations, but are mere convention. It would in principle be possible to derive a PBME method using two projections. However for convenience, we will use eqn (16) as the denition of PBME and LSC-IVR throughout this work. Finally, it is important to note that at least one of the operators has to be Wigner transformed in its projected form in order to ensure that the dynamics are initialised to the physical subspace.
While both LSC-IVR and PBME, as well as other mixed quantum-classical methods, have been applied to challenging systems with considerable success, their failure to accurately reproduce population dynamics in the long time limit has been well documented. 18,39,46,47 As mentioned above, a number of modications to quasiclassical methods which aim to address this issue have been proposed. [34][35][36][37][38][39]48 In practice, both eqn (16) and (15) are evaluated by averaging over an ensemble of trajectories, propagated with Hamilton's equations of motion dened by H: Initial conditions for each trajectory are sampled from r w b (x,p) for the nuclei and f a (X,P) for the mapping variables.

The Fenna-Matthews-Olson Hamiltonian
The Fenna-Matthews-Olson complex is a pigment protein biomolecule found in green sulfur bacteria adapted for low-light environments. It consists of three identical trimers, each containing seven bacteriochlorophyll (BChl) pigments supported by a protein backbone. In photosynthesis the task of FMO is to transport the excitation gained from absorbing sunlight to the reaction centre where it is converted into electrochemical energy. [40][41][42]44,49 The Frenkel-exciton model for the energy transfer in FMO is a challenging benchmark for quantum dynamics methods. In comparison to the Spin-Boson systems studied in our previous work, 39 the FMO Hamiltonian presents a different kind of challenge to quasiclassical dynamics methods: the electronic subsystem is comprised of more electronic states and the system-bath coupling is different. We note that the key challenge resulting from a larger electronic state space is the possibility of reaction chains involving more than two states. As a result this benchmark and the FMO system in general has been extensively studied using a considerable number of approaches. 16,20,[36][37][38]48,[50][51][52][53][54][55][56][57][58][59][60] In addition, though computationally challenging, numerically exact results are available, e.g. from hierarchical equations of motion (HEOM). [42][43][44]61,62 In the seven-site model (S ¼ 7), the full FMO Hamiltonian is given bŷ whereĤ s is the electronic sub-system Hamiltonian, given bŷ where 3 n is the energy of BChl site |ni and D nm is the electronic coupling between sites |ni and |mi. The values of site energies and couplings used in the matrix representation ofĤ s are given bŷ with all energies being in units of cm À1 . The protein environment around every BChl site is modelled by a bath of harmonic oscillators. The system-bath Hamiltonian,Ĥ sb , denes the coupling between the electronic sub-system and these baths. It is given byĤ where c j (n) is the vibronic coupling coefficient between site |ni and bath mode j. B is the number of modes per bath, such that B ¼ S Â F. The position coordinate of bath mode j of the nth bath is x j (n) . Finally, the Hamiltonian for the baths,Ĥ b , is given byĤ where p j (n) and u j (n) are the momentum coordinate and frequency of bath mode j associated with site |ni. The choice of masses does not affect results, so one can effectively set m j ¼ 1. Note that, following previous work, 38,52,61,62 each BChl site is coupled to an identical bath, which in turn is uncoupled from all other baths. The coupling between sites is thus contained purely inĤ s . The frequencies, u j (n) , and coupling coefficients, c j (n) , which are identical for each bath, are drawn from a spectral density of the Debye form, given by where u c is the characteristic frequency of the bath, related inversely to the phonon relaxation time, u c À1 ¼ s c , and we use l ¼ 35 cm À1 throughout, following previous work. 38,52,61,62 We discretize this function using a scheme known to reproduce exact reorganisation energies. 63,64 We dene the initial nuclear density matrixr b ¼ e ÀbĤb /Z b , where the partition function Z b is dened such that the trace over bath modes only is Tr b [r b ] ¼ 1. Nuclear positions and momenta were sampled from the thermal Wigner distribution of the uncorrelated bath, given, for any bath, by

Simulation parameters
In order to test our alternative denition of the quasiclassical population in eqn (15), we investigated three parameter regimes of the FMO Hamiltonian, which have been studied extensively, including using the numerically exact HEOM approach. [42][43][44]61,62 All our simulations used a timestep of dt ¼ 1 fs, which was found to be numerically converged. The results presented here are averaged over an ensemble of 10 6 trajectories in order to demonstrate the converged performance of our approach. We found however that using as few as 10 3 trajectories was enough to qualitatively capture all signicant features of the population dynamics and already exhibits a clear improvement over a fully converged traditional quasiclassical result. We note that in all our simulations, we used the traceless form of theV (x) matrix to propagate our trajectories, absorbing the remainder into U(x). Finally, while the results presented here were obtained using the integrator presented in ref. 23, we have also tested the integrators discussed in ref. 65 and found our approach to be insensitive to the choice of integration scheme.

Constituent correlation functions
Fig . 1 shows the constituent correlation functions, C IQ n (t) and C Q m Q n (t), calculated with electronic initial conditions having been sampled from both f(X,P) and f 2 (X,P). Considering the overall expression for the population given in eqn (14), the magnitudes of the constituent correlation functions are as one might expect. Notably the negative values observed for both C IQ n (t) and C Q m Q n (t) are not unphysical, as exact quantum mechanics would yield similar magnitudes for both correlation functions. We note that there is a noticeable difference between the correlation functions obtained from the two different initial distributions of the mapping variables we investigated. In order to assess their  comparative accuracy however they must be combined, using eqn (14), into a population and compared to exact results. Fig. 2 shows the population dynamics resulting from combining the constituent correlation functions using initial conditions sampled from f(X,P), i.e. the solid lines in Fig. 1. In addition to results obtained for an initial excitation of the |1i site, populations starting in site |6i are also shown. In both cases the traditional PBME populations, calculated as per eqn (16a), are shown for comparison, along with numerically exact HEOM results. [42][43][44]61,62 The comparison to PBME is a natural one here, as, in practice, the electronic initial conditions of PBME are also sampled from f(X,P), i.e. a ¼ 1. It is worth noting that this parameter regime of the FMO Hamiltonian is the most challenging for quasiclassical methods, due to the signicant impact of quantum effects at low temperatures. Our alternative denition of the population operators results in dynamics strikingly close in accuracy to the HEOM benchmark. We reproduce not only the correct ordering of states, even in the long time limit, but also capture all features present in the benchmark. We note that while we do observe negative populations for P 6)1 ðtÞ; their magnitude is almost negligible. In addition, our values for P 6)1 ðtÞ are within the same margin of error of the exact HEOM result as every other population. This is especially encouraging when comparing the accuracy of our approach to that achieved with traditional PBME. The latter, while capturing the population dynamics at short times rather well, completely fails to reproduce the long time behaviour. Notably the distribution and ordering of states beyond the short time limit degrades drastically with this approach. Considering that the low-temperature parameter regime of the FMO Hamiltonian poses a considerable challenge to quasiclassical methods, the accuracy of our results is highly encouraging. Fig. 3 shows populations for the same parameter regime of the FMO Hamiltonian as Fig. 2 however with electronic initial conditions now having been  sampled from f 2 (X,P). This corresponds to adding the dotted lines of Fig. 1 as per eqn (14). Also shown are standard LSC-IVR results, which again are a natural source of comparison as they use the same electronic initial conditions. Comparing to Fig. 2, LSC-IVR is seen to perform better than PBME for this particular parameter regime. Interestingly, the results achieved with our approach, using initial conditions sampled from f 2 (X,P), do not give a considerable improvement over LSC-IVR. We thus conclude that the most consistently accurate method is that shown in Fig. 2 using our improved population operators and sampling from f(X,P). This is the approach which we shall continue to employ for the remaining calculations. Fig. 4 and 5 show the population dynamics of the FMO Hamiltonian, calculated with our alternative denition of the population operators and, as in Fig. 2, with electronic initial conditions sampled from f(X,P), for two additional parameter regimes at T ¼ 300 K. Traditional PBME is again shown for comparison along with numerically exact HEOM benchmark results. [42][43][44]61,62 We note that the dynamics in the parameter regimes shown in these two gures are, owing to the higher temperature, less likely to be affected by nuclear quantum effects. Fig. 4 FMO site populations, with T ¼ 300 K, and a slower bath, s c ¼ 166 fs. Initial conditions for both constituent correlation functions sampled from f(X,P). Results are presented as in Fig. 2 and compared to numerical HEOM benchmark results. 44 Nevertheless it is clear that in the long-time limit, PBME diverges signicantly from the HEOM benchmark and yields an incorrect distribution of states.

Population dynamics
Using our alternative denition of the population operator again drastically improves the traditional quasiclassical result in both cases. Our approach in fact yields dynamics which now approach quantitative accuracy with respect to the exact HEOM benchmark. We furthermore note that the issue of small negative populations observed for our approach observed in Fig. 2 has now disappeared. This is not surprising, given that low temperature systems are well known to constitute a more considerable challenge for quasiclassical approaches. We recognise that quasiclassical approaches are well known not to capture nuclear effects, owing to the fact that the trajectories underlying them are driven by classical equations of motion. We note however that our alternative denition of the population operator is in fact not limited to quasiclassical methods, but may be applicable to other approaches based on mapping, which can capture nuclear quantum effects, such as nonadiabatic ring polymer molecular dynamics. [21][22][23][27][28][29]

Populations in the long time limit
One of the well known failings of quasiclassical dynamics methods is that they do not preserve detailed balance of the populations and are therefore inaccurate at long times. In order to investigate whether our alternative denition of the population operator can improve on this, we have carried out a longer simulation of the parameter regime shown in Fig. 2 and 3. We used the same timestep, dt ¼ 1 fs, as in the simulations above and averaged over the same number of trajectories (10 6 ). Fig. 6 shows the FMO site populations, following an initial excitation of either the |1i or the |6i site, up to 10 ps, calculated with our denition of the population operator. Electronic initial conditions for both were sampled from f(X,P).
Comparing our results to the long-time limit of the HEOM benchmark, 44 it is clear that our method does not rigorously preserve detailed balance. In addition, one of the populations is unphysically predicted to be slightly negative. We note however that as in Fig. 2 and 3, our negative result is within the same margin of error of the exact result as every other state population. It is well known that the long-time limits of the populations obtained from traditional quasiclassical approaches such as PBME and LSC-IVR can be much worse, predicting more strongly negative populations, and some greater than 1. 50 We are therefore encouraged that the equilibrium distribution predicted using our approach is qualitatively similar to the HEOM benchmark. We note furthermore that our approach predicts identical equilibrium distributions, whether site |1i or site |6i is initially excited. We have not computed the long-time dynamics of the other two parameter regimes we studied above. We do however expect that, owing to their higher temperature and the lower impact nuclear quantum effects are therefore likely to have on them, our approach would perform even better than in the T ¼ 77 K system. The fact that we do not observe negative populations in Fig. 4 and 5 further supports this hypothesis.
Overall we consider these population results, along with the others shown above, to be highly encouraging. They clearly demonstrate that our denition of

Conclusion
We have outlined an extension of our previous work, 39 presenting an alternative denition of the electronic population operator for any system of multiple electronic states. We rely on the fact that any Hermitian operator can be split into two terms, one of which is the identity, the other being traceless. 18 We then use our understanding of the exact behaviour of the quantum identity instead of a quasiclassical treatment. The combination of this splitting and exact treatment of the identity results in a new form of the electronic population operator. Our approach retains the excellent scaling with respect to system size of traditional quasiclassical methods as well as their underlying equations of motion. Notably, as the constituent correlation functions into which our new operator is split can be calculated for all states in a single simulation, our approach is no more computationally expensive than the traditional methods it seeks to improve.
We have applied our approach to the challenging seven-state Frenkel-exciton model of the FMO light harvesting complex. [40][41][42][43][44] In addition to having been studied extensively with traditional quasiclassical methods, 50-52 the fact that numerically exact quantum results are available 42-44,61,62 makes this system an ideal benchmark for our modication of the traditional quasiclassical population operators.
Overall we nd that using our alternative denition of the electronic population operator drastically improves on the results obtainable with other quasiclassical methods. In addition, our results actually approach the exact quantum benchmark in accuracy for the three parameter regimes we study. Finally, we nd that rather encouragingly, our method reproduces the long-time distribution of the electronic states with far higher accuracy than existing quasiclassical approaches.
We recognise that there have been other efforts to x the well documented shortcomings of traditional quasiclassical dynamics methods. For instance, the symmetrical quasiclassical windowing approach uses "binning" to convert the population into an integer, using a windowing function applied symmetrically at the beginning and end of each classical trajectory. [34][35][36][37][38] This approach has been applied to the FMO Hamiltonian we study here with considerable success, yielding results comparable in accuracy to those presented here. [36][37][38] In recent work a post-processing method for the traditional LSC-IVR quasiclassical method has been proposed. 48 The dynamics resulting from the traditional approach are shied by a function which imposes the long-time Boltzmann distribution of the FMO subsystem Hamiltonian. We note that while the results obtained from this long-time correction do constitute an improvement over the traditional approach, they fail to address any inaccuracies at short to medium times. In addition, this approach relies on either having prior knowledge of the correct long-time distribution of states or approximating it.
Other mixed quantum-classical dynamics methods have also been applied to the FMO Hamiltonian with considerable success. In recent work both the forward-backward trajectory solution 19,20 and the partially linearised density matrix 16,17 approaches have yielded accurate results for the FMO systems studied here. We note however that these methods use a different set of equations of motion to propagate the classical trajectories from the quasiclassical approaches used here. This does however not disqualify them from also beneting from our alternative denition of the population operator, as the latter is independent of the equations of motion.
While we have shown that our alternative denition of the electronic population operator can drastically improve on the dynamics obtained from traditional quasiclassical methods, it cannot address all their shortcomings. Notably, it cannot capture nuclear quantum effects, which are fundamentally inaccessible to an approach relying on purely classical trajectories to calculate operators. Due to its simplicity and generality however, our approach could be combined with methods which can capture some nuclear quantum effects such as tunnelling. The nonadiabatic ring polymer molecular dynamics 21-23,27-29 method would seem to be a logical candidate for beneting from our denition of the population operator, given that it is also based on the mapping formalism. We anticipate that the low temperature regime of the FMO Hamiltonian studied here may particularly benet from such a combination, as nuclear quantum effects are likely to have a greater impact in this system.
Overall we have shown that using our alternative denition of the population operator results in a considerable improvement over traditional quasiclassical approaches. Our results in fact approach the accuracy of the numerically exact benchmark for the systems we have studied. This is highly encouraging, given that the traditional methods we compare to have previously been considered inadequate for the simulation of long-time nonadiabatic dynamics. Our modication does not involve changing the equations of motion underlying quasiclassical methods and in fact also scales identically. We therefore hope that this work will further the development of dynamics methods based on the quasiclassical approach, and the progress of mixed quantum-classical dynamics as a whole.

Conflicts of interest
There are no conicts to declare.