Spin-Mapping Methods for Simulating Ultrafast Nonadiabatic Dynamics

: Many chemical reactions exhibit nonadiabatic effects as a consequence of coupling between electronic states and/or interaction with light. While a fully quantum description of nonadiabatic reactions is unfeasible for most realistic molecules, a more computationally tractable approach is to combine a classical description of the nuclei with a quantum description of the electronic states. Combining the formalisms of quantum and classical dynamics is however a difficult problem for which standard methods (such as Ehrenfest dynamics and surface hopping) may be insufficient. In this article, we review a new trajectory-based approach developed in our group that is able to describe nonadiabatic dynamics with a higher accuracy than previous approaches but for a similar level of computational effort. This method treats the electronic states with a phase-space representation for discrete-level systems, which in the two-level case is analogous to a spin- . We point out the key features of the method and demonstrate its use in a variety of applications, including ultrafast transfer through conical intersections, damped coherent excitation under coupling to a strong light field, and nonlinear spectroscopy of light-harvesting complexes.


Introduction
Molecules live in a quantum world and, in principle, knowledge of the quantum-mechanical wavefunction, , contains sufficient information to determine the outcome of almost any experiment (at least in the low-energy physics domain). However, the complexity of solving the full time-dependent Schrödinger equation, , limits state-of-the-art calculations to very small molecules. [1] It is therefore common in computational chemistry to completely neglect quantum effects and employ classical molecular dynamics ( ), which enables simulations of thousands (or even millions) of atoms. A particularly challenging problem is the simulation of nonadiabatic chemical reactions, including (but not limited to) ultrafast dynamics initiated by light. In contrast to adiabatic processes, which evolve on a single Born-Oppenheimer potential-energy surface, it is not possible to model nonadiabatic processes directly using classical molecular dynamics, due to the fact that multiple potential-energy surfaces are coupled together. Thus, in order to describe experiments involving complex molecular systems, it is necessary to develop methods that bridge the quantum and classical worlds by including the most important quantum effects but retaining as much as possible of the simplicity, computational efficiency and interpretability afforded by classical mechanics. The present article summarizes our recent developments together with a few representative applications.
In many cases, we wish to use a classical description for the nuclei and a quantized treatment of the electronic-state population amplitudes. The mixed quantum-classical Hamiltonian may be written as where is a state-independent potential, and is a matrix of potential-energy surfaces and diabatic couplings. Note that such mixed quantum-classical models are not just limited to describing electrons and nuclei but are applicable also in many other cases, for example atoms interacting with photons or the decoherence of qubits in a noisy environment. In itself, mixing quantum and classical dynamics is not a new idea and one can trace similar concepts back to the foundations of quantum mechanics. The complication is that the equations of motion generated by this mixed Hamiltonian are not uniquely determined.
The simplest example of a mixed quantum-classical approach is provided by a mean-field approach called Ehrenfest dynamics, in which one simply propagates the electronic wavefunction, , using the Schrödinger equation and the nuclear positions and momenta, (x, p), using Hamilton's equations of motion with the expectation value of the force: Although this approach can correctly describe the electronic and nuclear motion separately, it neglects most of the correlations that should exist between them. [2] It is therefore widely known to lead to incorrect predictions for nonadiabatic dynamics as well as to severely violate detailed balance. [3] Much effort has been invested in developing methods that improve upon Ehrenfest dynamics. The most popular today is surface hopping, [4] which has become the standard tool for simulating ultrafast processes encountered in photochemistry. Unfortunately, no rigorous derivation has been given and the algorithm is known to suffer from a number of serious problems, such as a lack of decoherence. [5] Although many different "decoherence corrections" have been suggested, it is clear that this is an ad hoc fix which may help in some cases but is not guaranteed to work in general.
classical magnitude = . This choice guarantees the classical expectation values of products of operators to be equal to their quantum traces, [22,23] where × sin d d is an integral over the sphere. This provides a phasespace formulation for two-state systems, which is known as a Stratonovich-Weyl representation and is analogous to the Wigner representation of continuous systems. The relation does not hold when using the Q-sphere but is unique for the W-sphere. A comparison between the two spheres is shown in Fig. 1. Note that the W-sphere (right panel) is purely in state |1〉 on the "northern polar" circle (or in state |2〉 on the "southern" circle) rather than at the poles. These circles are positioned at an angle c = arccos ≈ 54.7°which is notably identical to the "magic angle" used in NMR spectroscopy.
Our spin-mapping dynamics are initialized by sampling x and p from a Wigner distribution and uniformly from the W-sphere (in contrast, Ehrenfest dynamics would start from a single point). For each set of initial conditions, we run a trajectory according to Eq. (4), measure an observable B( ) at time t and average the results over the ensemble. In the case that the zero-time observable is a population, one can alternatively employ so-called "focused" initial conditions by sampling the initial from the corresponding polar circle rather than the entire sphere.
As an example of how the different spin constructions impact the dynamics, we have calculated the population dynamics in a model of a two-level system linearly coupled to a harmonic bath (called a spin-boson model), for which the fully quantum result can be computed using the quasiadiabatic path-integral method. [24] As shown in Fig. 2, the spin-mapping approach is found to give more accurate results than the corresponding calculations using Ehrenfest dynamics or MMST mapping (which correspond to dynamics on a smaller/larger sphere, respectively). All these three approximations employ the same form of the equations of motion and have roughly the same level of computational effort, yet the results are of different accuracy. Further results for a range of spinboson models lead to similar conclusions and are given in ref. [19].
Formally, the classical-like dynamics presented above is known as a linearization of the fully quantum propagation. [25] One can go one step further and apply this linearized approxima- An alternative strategy is to map the quantum system onto a problem that has a well-defined classical analogue, in order to define a consistent classical limit where the quantum and classical degrees of freedom are treated on the same footing. One such mapping, proposed by Meyer and Miller and formalized by Stock and Thoss (MMST), [6,7] expresses the electronic operators in a second-quantization picture. The resulting creation and annihilation operators are written in terms of the positions and momenta of fictitious harmonic oscillators, which are then treated classically in the same way as nuclear modes. This approach has been used to systematically derive approximations that improve upon Ehrenfest dynamics. [8] However, these approximations introduce new inconsistencies associated with the fact that only the singly-excited oscillator states map back to the electronic states, whereas all other excitations are unphysical. [9] Although certain tricks have been able to overcome or at least mitigate some of these problems, [10][11][12][13][14][15] it is clear that mapping to harmonic oscillators cannot be the final answer.

Spin Mapping
As an alternative to the MMST mapping, a number of authors have considered mapping the system to spin degrees of freedom. [16][17][18] However, these methods originally suffered from various problems and were therefore abandoned in favour of the harmonic-oscillator mapping mentioned above. In recent work, we have shown how these problems can be overcome to obtain a practical and accurate "spin-mapping" method. Here, we only briefly introduce this theory; more details can be found in refs. [19,20].

Outline of the Theory
For a two-level system, it is well-known that the system and its potential are equivalent to a spin-in a (fictitious) magnetic field. The system state is then represented by a vector on the Bloch sphere. In comparison to a set of harmonic oscillators, the spinhas the advantage that its Hilbert space has the same size as the original problem (two states) whereas the oscillators introduce an infinite number of unwanted states. It is common to represent the state of a two-level system by the vector , whose components are phase-space representations of the Pauli spin operators. The corresponding classical representation of the Hamiltonian is leading to the equations of motion In other words, all degrees of freedom are coupled such that the spin precesses around the fictitious magnetic field V(x), which changes as the nuclei evolve on a potential defined by the instantaneous direction of the spin.
Different methods can be distinguished by their construction of . If one chooses to be a vector on the usual Bloch sphere, then one recovers Ehrenfest dynamics. Here the system is purely in state |1〉 on the north pole, in state |2〉 on the south pole, or in a superposition state everywhere else. In this way, any wavefunction of a two-level system can be represented by a spin coherent state. In phase-space terminology, this is analogous to the Q symbol, also known as a Husimi distribution. [21] This is however not the only possible choice: in particular, it is beneficial to pick the magnitude of the spin to be a factor √3 larger than the Q-sphere, where the numerical factor appears because a spin s = 1/2 has the Fig. 1: Comparison of the "spin spheres" used in Ehrenfest dynamics (left) and spin mapping (right). The Ehrenfest sphere can be thought of as the usual Bloch sphere of radius 1, where the populations |1〉 〈1| and |2〉 〈2| correspond to the north and the south pole, respectively. The spin-mapping sphere (right) is a factor √3 larger, and populations correspond to "polar" circles.

Frontiers in ultraFast spectroscopy and dynamics
Here, 〈·〉 0 denotes the average with respect to the initial distribution 0 , whereas 〈·〉eq employs the thermal equilibrium distribution eq = e − H /Z according to the mapping Hamiltonian [Eq. (3)], where = 1/k B T is the inverse temperature and Z the partition function. Equation (6) relates dynamical properties to statistical averages. It implies that the accuracy of the long-time predictions of an ensemble of mapping trajectories depends on how accurately the method approximates the correct thermal distribution. To investigate this issue, we consider as a simple example a quantum system consisting of two levels separated by an energy gap 2 in thermal equilibrium. The correct population of the excited state according to the quantum Boltzmann distribution should be e −2 /(1+e −2 ), where = denotes the only relevant dimensionless parameter in this problem.
In Fig. 3 we compare the average population of the excited state for this system versus from three mapping methods, namely spin mapping, Ehrenfest and focused MMST mapping, the key difference being the radius of the sphere used in each case. Formally, it can be shown that spin mapping is able to predict the correct population up to second order in powers of , while the other two methods are only correct to zeroth order. [30] This provides a mathematical explanation for why spin mapping typically gives accurate predictions of population relaxation in condensed-phase systems, at least at high temperatures, and allows us to determine the parameter regime in which it is expected to be a reliable approximation.
For increasing values of we notice a dramatic break-down of the predictions of both spin mapping and MMST. As the relative energy difference between the two levels increases, the system is more likely to occupy the two "polar" regions of the electronic phase space. Those regions correspond to states where the population of the excited level acquires a negative value. For small values of the existence of these regions does not cause a problem for spin mapping (in fact they are the reason for its success) and even if some individual trajectories enter the polar region, the ensemble average remains accurate. For larger values of , however, the strong magnetic field distorts the ensemble average causing the populations to become negative. In the case of Ehrenfest these regions do not exist and hence the population of this method is guaranteed to remain positive for arbitrary values of . However, although the absolute error of Ehrenfest is small in the large-regime, the asymptotic behaviour is nonetheless completely wrong, leading to large relative errors.
Luckily, in many cases, the most interesting regime for nonadiabatic dynamics is where is not too large. This is because tion only to the nuclei and not to the electrons, called partial linearization, which is more accurate but typically requires an order of magnitude more trajectories to converge. [26,27] While the fully linearized dynamics can predict population dynamics reasonably well, the partially linearized dynamics are especially well-suited for the computation of linear and nonlinear spectra, as demonstrated in Section 3.4 below. This is a particular advantage over surfacehopping methods, which are not naturally applicable to simulate dynamics of the coherences created in optical spectroscopy.
How can this approach be generalized to more than two levels? As a first attempt, one may expect to look for a mapping to higher spin quantum numbers, such as a spin-1 for three levels. However, there is an important distinction between these two systems, in that the spin-1 in a magnetic field always has equally spaced energy levels, which is not sufficiently general to describe arbitrary threelevel systems. From a group-theoretical perspective, all spins have symmetry SU(2), while an N-level system has symmetry SU(N). For this reason, higher spins do not provide the most natural extension of the theory outlined above. Instead, one can construct a phase-space representation of an N-level system by replacing the factor √3 by . This generalized spin-mapping approach guarantees that Eq. (5) remains true for any number of levels. [20]

Thermalization Properties
The results in Fig. 2 indicate that spin mapping is more accurate than both Ehrenfest and MMST for the prediction of timedependent populations. In this section, we investigate the reasons behind these observations and show that the radius used for the spin-mapping sphere is optimal for predicting long-time thermalization properties for an important class of nonadiabatic systems.
The mapping approaches discussed in the present work approximate quantum nonadiabatic dynamics with the evolution of Hamilton's equations of motion in an extended phase space. This feature allows us to use concepts and results of classical ergodic theory. In particular, we expect that the long-time limit of any arbitrary initial distribution will be of the canonical form, provided energy can be identified as the only constant of motionthat is, if the dynamical system is ergodic. [28] We additionally assume that the mapping dynamics is strong mixing [29] on the surfaces of constant energy, meaning that the time-correlation functions between any two observables A and B will relax to the limit sue is not further investigated. An advantage of the spin-mapping method is that we can obtain more accurate results without these ad hoc corrections.

Population Dynamics in Light-harvesting Complexes
Over the last decades, there has been an intense discussion in the literature about whether quantum effects are necessary in order for biological systems to function. In our opinion, one of the reasons for why this has become such a controversial question is because there has not yet been a clear-cut way to provide an answer. For instance, one might find a large discrepancy between an Ehrenfest simulation and the experimental result and conclude that it arises from the neglect of nuclear quantum effects. However, an alternative explanation would simply be that the well-known problems with detailed balance are the cause of the discrepancy. It is therefore necessary to use more accurate mixed quantum-classical approaches like spin mapping to answer these questions.
One of the prototypical examples in the study of photosynthesis is the Fenna-Matthews-Olson (FMO) light-harvesting complex found in green sulfur bacteria. A simple excitonic model of this complex consists of seven sites coupled to a protein environment of hundreds of vibrational modes. Fig. 5 demonstrates for 1 the problem would be better and more simply described by a single adiabatic ground state, as the excited state is so much higher in energy. Nonetheless, it would still be useful to develop a method which is uniformly accurate across the whole range. Further improvements of the thermal predictions of spin-mapping methods are currently under investigation in our research group.

Ultrafast Dynamics in the Benzene Cation
A challenging problem for nonadiabatic simulation methods is to describe transfer through a conical intersection. As an example of such a process, we investigate the nonradiative decay of the benzene cation from its second excited state to the ground state . Fig. 4 compares the population dynamics using a few different trajectory-based methods for a linear vibronic-coupling model with three states and five modes, for which exact quantum benchmarks are available. Again, it is clear that the dynamics on the spin-mapping (W) sphere is more accurate than using Ehrenfest dynamics or MMST mapping. Also shown for comparison is surface hopping, which is found to deviate from the correct behaviour already after the first period of population oscillations. It is possible that the surface hopping results reported here [8] could be improved by employing a decoherence correction, but this is- Fig. 4: Population dynamics in the benzene cation following an initial excitation to . The dynamics was propagated on the three-state five-mode model by Köppel. [31] Figure adapted from ref. [30]; surface-hopping results were taken from ref. [8]. Both the spin-and MMST-mapping results used focused initial conditions. shows the spatial configuration of the seven sites included in the model. The simulations were carried out using "focused" initial conditions as described in ref. [20]. The fully quantum HEOM results are from ref. [32]. that spin-mapping dynamics together with a completely classical description of the nuclear vibrations is sufficient to obtain the fully quantum result, which in this special case can be computed using the hierarchical equations of motion (HEOM). [33] In a world that obeyed Ehrenfest dynamics, photosynthesis could not function efficiently as there would be no preferential flow of energy through the complex from site 1 to sites 3 and 4. However, in this case, spin mapping provides a better approximation to reality and can be used to understand the mechanism in the simplest way. This conclusively demonstrates that this key step of photosynthesis does not require quantum-mechanical effects of the vibrations in order to function [34] and thus refutes a commonly-held belief to the contrary. [35]

Driven Damped Rabi Oscillations
The photo-induced reactions in benzene and FMO presented above were described with the assumption that the excitation pulse was instantaneous. However, this is not always a valid approximation, and in order to study the dynamics during as well as after the photo-excitation step, in strong fields as well as for arbitrary-shaped pulses, we need to explicitly include the coupling to a light field in our theory. To this end, we add the electric-dipole term -to the Hamiltonian, where is the dipole operator and E(t) the time-dependent electric field at the position of the molecule. [37] In this way, the spin-mapping equations of motion [Eq. (4)] obtain an explicit time dependence.
Note that we do not need to make any assumptions regarding the functional form of E(t), meaning we can simulate any pulse shape, strength or sequence of pulses we like. Also, our approach does not rely on the rotating-wave approximation (RWA ), so that we do not have to constrain ourselves to studying low light intensities and (near-)resonant frequencies. By taking only the electricdipole term as the interaction Hamiltonian, we have however effectively introduced the long-wavelength approximation. This is typically valid as long as the wavelength of the light is very large compared to the system dimensions. However, it would in principle be possible to include a multipole expansion [38] in order to treat special cases in which magnetic effects also play a role, such as in strong-field ionization. , where = 0.025, which is coupled to an Ohmic bath with Kondo parameter = 0.1, characteristic frequency c = 0.5 and inverse temperature = 1.0. The system is further coupled to a light field by , with the frequency = 2 on resonance, and = 1.0 > ω, such that the RWA is not valid. Also note that the characteristic timescale of the bath (~1/ c ) is comparable to that of the system (~1/ ), such that the dynamics is non-Markovian. We computed the numerically exact dynamics using HEOM with the open-source pyrho package, [36] both with and without the RWA.
As an example, we simulate the population dynamics of a twolevel system coupled to an environment of nuclear vibrational modes that is continuously driven by a resonant light field. This results in damped Rabi oscillations, as shown in Fig. 6. We have chosen the parameters such that we are in a regime where the commonly-used Markovian and RWA approximations break down. This means that if we do make the RWA (indicated with "RWA " in the figure), this will not yield a good estimation of the exact dynamics, and this is indeed what can be observed: the HEOM result within the RWA does not even qualitatively agree with the HEOM result without the RWA.
As mentioned, we do not need to make the RWA to be able to use spin mapping, and the spin-mapping result without the RWA can be seen to agree very well with the exact result.

Nonlinear Spectra
Many approaches exist for computing optical spectra involving electronic transitions between the ground and a single excited state. [39][40][41][42] A particularly convenient approach is the Wigneraveraged classical limit (WACL), which computes the associated real-time correlation functions by propagating nuclear trajectories on the time-independent arithmetic-mean surface of the ground and excited state. [43,44] However, calculating optical spectra for systems involving multiple coupled excited states is much more challenging, as the underlying correlation functions now involve nonadiabatic dynamics within the coupled excited-state manifold. To tackle this problem, we have recently developed a spin-mapping approach that reduces to WACL in the limit of uncoupled excited states, but otherwise incorporates the correct nonadiabatic spectral features away from this limit. [45] So far in this article we have considered electronic population dynamics, which we have shown can be well described by fully linearized approaches. This is because electronic population oper- ators are diagonal and so their time-evolution is well approximated by considering only the average electronic path associated with the forward and backward propagators, . To calculate optical spectra, correlation functions involving the electronic transition dipole operator must be computed. Because this operator is not diagonal (i.e. it connects the uncoupled electronic subspaces containing different numbers of excitations, commonly referred to as excitons), its time-evolution is best described using partially linearized approaches, [46,47] where the distinct subspace dynamics associated with the forward and backward electronic propagators can be treated separately. The partially linearized approach based on spin mapping is called spin-PLDM. [26,27] A key advantage of partially linearized approaches over their fully linearized counterparts (or alternative approaches such as surface hopping) is that they can also be used to compute multitime correlation functions required for nonlinear ultrafast spectroscopy, [45,48] such as pump-probe and 2D optical spectra.
The various multi-time correlation functions that contribute to these spectra are illustrated in Fig. 7 by their Feynman diagrams. [44] While the rephasing (RP) and non-rephasing (NR) signals can be separated in an experiment, the different Liouvillespace pathways associated with the ground-state bleach (GB), stimulated emission (SE) and excited-state absorption (ESA) cannot. Because these pathways are evaluated independently of each other in our approach, the final signal can be decomposed into sums and differences of various mechanistic contributions, offering further insight beyond what can be obtained directly from experiment. Fig. 8 shows the pump-probe spectra for FMO at T = 300 K for different t 2 delay times. The fact that a static approximation (in which the nuclei do not evolve during the entire simulation) is unable to correctly reproduce the qualitative features of the spectra illustrates the importance of including nonadiabatic effects in describing nonlinear ultrafast spectroscopy. In contrast,    Fig. 9: The SE contribution to the 2D optical spectra for the sevenstate FMO model at T = 300 K, calculated for different t 2 delay times. [49] Numerically exact results are computed using HEOM. [36] spin-PLDM accurately reproduces the pump-probe spectra. By separating the full signal into its various Liouville-space pathway components we can obtain a better understanding of the t 2 dependence of the signal. For example, the GB signal is independent of the delay time, while although the initial SE signal is identical to GB, in the large t 2 limit it corresponds to fluorescence from the lowest-energy singly-excited electronic states (in accordance with Kasha's rule). The ESA signal is similar but originates from absorption into the higher-energy unoccupied states. 2D optical spectra offer additional information compared to the pump-probe, such as the couplings between the various electronic states. Spin-PLDM can also be used to simulate these results (Fig. 9), where it reproduces all of the key features of the quantum-mechanical benchmark for the same FMO model. These positive results suggest that, given a realistic Hamiltonian, our approach can be used to interpret and better understand experimental spectroscopic results.

Outlook
In this article, we have shown that the spin-mapping framework can be used to obtain accurate predictions of nonadiabatic dynamics. The method is based on classical trajectories, which provide interpretations of reaction mechanisms and are computationally affordable also for large systems. We have applied the method to simulate ultrafast population dynamics, either initiated in an excited state or driven by a strong laser field, as well as nonlinear spectra.
All of the examples given above were found to be in good agreement with fully quantum benchmark calculations. In particular, the method can dramatically improve upon standard approaches based on mean-field, static, Markovian or rotatingwave approximations. In future work, we will apply the method to larger and more complicated systems including anharmonicity, for which it is not feasible to employ fully quantum calculations.
Finally, we point out that one can go beyond the classical approximations of spin-mapping dynamics by representing the electronic state by a path integral of spins rather than a single spin. [50] This path-integral version of the method is able to properly describe wavepacket splitting and other examples of entanglement between the system and the environment, and is free from the previously observed overcoherence problems in surface hopping. As well as being a useful method in its own right, the spin-mapping concept therefore also provides a foundation for the development of improved theories of nonadiabatic processes.

License and Terms
This is an Open Access article under the terms of the Creative Commons Attribution License CC BY 4.0. The material may not be used for commercial purposes. The license is subject to the CHIMIA terms and conditions: (https://chimia.ch/chimia/about).
The definitive version of this article is the electronic one that can be found at https://doi.org/10.2533/chimia.2022.582