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

The coupled dynamics of electrons and nuclei in molecular condensed-phase systems remains a challenging problem for computer simulation.^1–4 Nonadiabatic transitions can be induced when the electronic adiabatic states of a system become close in energy at some nuclear geometry, which results in a breakdown of the Born–Oppenheimer approximation.

One approach is to use numerically exact wavefunction-based methods to compute the dynamical observables of interest. The advantage of such methods is that the obtained results are systematically improvable on increasing the size of the basis; they can also in principle describe nuclear quantum effects, such as tunneling and zero-point energy. Of particular promise are methods such as the time-evolving block decimation (TEBD) technique,⁵ which can even be applied to fairly large systems.⁶ However, such methods often have restrictions on the type of systems that they can treat. For example, TEBD can only be efficiently applied to spatially linear systems with, at most, nearest-neighbor couplings.

Due to the continuous nature of the nuclear subspace, the dynamics of such degrees of freedom are ideally suited to be performed using classical trajectories. The advantage of using such methods is that they are numerically cheap and can be easily applied to large, condensed-phase problems. Of these methods, the most popular are Tully’s fewest-switches surface hopping,⁷ which is heuristically motivated and hence cannot be rigorously derived,⁸ and Ehrenfest dynamics,^9,10 which neglects all entanglement between the electronic and nuclear subsystems.¹¹ 

Recently, there has been a renaissance in mapping-based techniques,^12–25 which describe the dynamics of the electronic subsystem as well as the nuclei using classical trajectories within a continuous mapping space. Such methods have been shown to have superior accuracy compared to Ehrenfest dynamics when calculating dynamical observables for a wide range of model systems. While these methods still treat the nuclear degrees of freedom classically and hence neglect nuclear quantum effects, these could be perhaps reintroduced by using ring-polymer based formulations.^26–29 Although none of these approaches will be able to describe quantum nuclear coherence effects,³⁰ such effects are typically unimportant in condensed-phase systems.⁴ 

Historically, the vast majority of mapping-based approaches have used the so-called Meyer–Miller–Stock–Thoss (MMST) mapping.^31,32 In this scheme, the electronic subsystem is mapped onto a set of singly excited harmonic oscillators. The nuclear and electronic dynamics are then described as an average over many classical trajectories. While trajectory-based methods using this mapping have been able to qualitatively reproduce the correct dynamics in a range of model systems, such methods may also exhibit significant problems. One of the most important of these is zero-point energy leakage, where zero-point energy can unphysically flow between the mapping harmonic oscillators as a result of the trajectories leaving the physical subspace during the classical dynamics.^33,34 In other words, the trajectories can evolve into areas of phase space in the mapped system that do not correspond to a valid state in the electronic system. To alleviate this problem, related approaches have been suggested, such as reducing the zero-point energy parameter in the underlying theory,³⁵ symmetric windowing^{12,15,36–45} of trajectories, an MMST identity-trick,^23–25 and using alternative classical mapping models.^18,46

Recently, a different form of mapping based on classical spin dynamics^31,38 has been reformulated in terms of Stratonovich–Weyl kernels.^47,48 While this spin-mapping approach leads to exactly the same equations of motion for the trajectories as MMST mapping,³¹ the Cartesian mapping variables are constrained to a hypersphere, which is isomorphic with the phase space of the actual electronic subsystem. This means that unlike standard MMST mapping, spin-mapping does not need projection operators onto the physical subspace and trajectories no longer suffer from zero-point energy leakage. In addition, when converted to Cartesian mapping variables, the spin-mapping Hamiltonian also has a different zero-point energy parameter to that of MMST mapping. Different zero-point energy parameters within MMST mapping have also been previously investigated by Müller and Stock³⁵ and have also previously been justified using analogies to spin by Cotton and Miller.^36,41

When derived from a path-integral formalism, previous mapping-based techniques generally fall into one of two categories. Fully linearized methods^4,25,49–51 result from performing a linearization approximation to the difference between the forward and backward paths for both the electronic and nuclear degrees of freedom; a semiclassical approximation that is expected to be valid in the classical limit. While this linearization approximation of the electronic paths is for many forms of mappings exact when applied to a purely electronic system, it constitutes an additional approximation when there is electron–nuclear coupling. In order to try to improve upon this, partially linearized methods⁵² result from only performing a linearization approximation for the nuclear paths and then treating the dynamics of the forward and backward electronic paths separately. Examples of partially linearized methods using MMST mapping are the partially linearized density matrix (PLDM) approach^{13,14,53–60} and the forward–backward trajectory solution (FBTS).^19–21 While spin-mapping has already successfully improved the accuracy of so-called fully linearized methods,^47,48 it has yet to be applied to partially linearized methods. Often partially linearized methods are better than their fully linearized counterparts, so we may optimistically expect a spin-mapping version of PLDM to be the best method of all. To aid the reader, a schematic illustrating the relationship between spin-PLDM and other mapping-based classical-trajectory techniques is given in Fig. 1.

In this paper, after a brief review of the standard PLDM method based on MMST mapping, we use spin-mapping to derive a new partially linearized approach, which we call spin-PLDM. The method is derived generally for systems containing any number of electronic states by employing the successful framework introduced in Ref. 48. We show that such a method is able to more accurately reproduce dynamical observables associated with both the spin-boson model and Fenna-Matthews–Olsen (FMO) complex than other commonly used mapping-based approaches. In Paper II,⁶¹ the spin-PLDM method is analyzed extensively to explain the reasons for its improved behavior over other methods within the same class.

In this paper, we study nonadiabatic dynamics using the specific example of the familiar electron–nuclear coupled system. Naturally, the method we derive in this paper could also be applied to any other quantum-classical system. The Hamiltonian for such a system can always be written in the following form:

where $V0(x^)$ is a state-independent potential function of nuclear configuration, x, whereas $V^(x^)$ has an F × F matrix representation in the basis of electronic states, |1⟩, …, |F⟩; these are uniquely partitioned so that $tr[V^(x^)]=0$⁠. We define F as the dimensionality of the electronic subsystem and tr[·] as the partial trace over the electronic degrees of freedom. Both x and p are vectors of dimension f, which have been mass weighted such that all degrees of freedom have the same mass m. We set ℏ = 1 throughout.

Many dynamical quantities of interest can be written in the form of a real-time quantum correlation function,⁶² 

where $Â$ is an electronic operator effectively initializing the system at t = 0 and $ρ^b$ is the initial density matrix for the nuclear bath. These are known as factorized initial conditions. In addition, $B^$ is an electronic operator that measures the system at time t. Although we only consider the case where $Â$ and $B^$ are purely electronic operators in this paper, the theory could be easily extended for more general real-time correlation functions. In addition, one can in principle go beyond factorized initial conditions using the more general density matrix $ρ^=∑jρ^b(j)Â(j)$⁠.

Inserting complete sets of nuclear position eigenstates (|x₀⟩, |$x0′$⟩, and |x_N⟩) and electronic basis states (|λ⟩ and |μ⟩) results in the following expression for the real-time correlation function:

where the unprimed coordinates are used for the forward propagation path, the primed coordinates are used for the backward propagation path, and x_N = $xN′$ because there is no nuclear operator applied at time t within the correlation function C_AB(t). One of the main challenges in evaluating this expression is obtaining an efficient expression for the forward and backward real-time propagators. When derived from a path-integral representation consisting of N time steps, all mapping-based classical-trajectory approaches first implement a linearization approximation⁴ for the nuclear degrees of freedom to avoid sampling complex phases. The linearization approach is also well known for single-surface systems where it leads to classical Wigner dynamics.^4,63,64 Hence, such a linearization approximation amounts to treating the dynamics of the nuclear degrees of freedom classically.

To implement the linearization approach,^49,65 the nuclear degrees of freedom must first be transformed to sum and difference coordinates, $x¯k=12(xk+xk′)$ and Δx_k = x_k − $xk′$⁠, respectively, where x_k (⁠$xk′$⁠) is the nuclear coordinate at the path-integral time step k, associated with the forward (backward) propagation path. Then, the approximation involves only retaining terms in the nuclear and electronic action up to first order in the difference coordinates, Δx_k. This results in the following approximate path-integral expression for the real-time correlation function:

where the initial nuclear density matrix is now described by its Wigner transform,

Defined in this way, the Wigner transform of the initial nuclear density is normalized such that $∫dx¯0dp¯0ρb(x¯0,p¯0)=1$⁠. Additionally, the linearization approximation to the difference of the nuclear action along the two paths is given by^13,49,63

where ϵ = t/N is the path-integral time step,

is the state-independent nuclear force, and ∇ is the gradient, a vector of derivatives with respect to nuclear positions.

The electronic action is included within the definition of the electronic transition amplitude, T_[μ,λ]. Retaining terms in the electronic action up to first order in the difference coordinates, Δx_k, the electronic transition amplitude for the forward and backward paths are approximated by

These quantities contain two types of propagators: the electronic propagator, $e±iV^(x¯k)ϵ$⁠, which describes the electronic dynamics, and the coupling propagator, $e−i∇V^(x¯k)ϵΔxk/2$⁠, which due to its dependence on Δx_k adds an electronic contribution to the nuclear force, as will be shown later.

With these approximate expressions for the nuclear action between the two paths [Eq. (6)] and the transition amplitudes for the forward and backward paths [Eq. (8)], the real-time correlation function would be an exact solution of the quantum-classical Liouville equation (QCLE).^66–68 However, it is known that it is not possible to implement such an expression in an independent trajectory simulation.⁶⁶ The easiest way to further simplify the expression for the real-time correlation function is by introducing a mapping variable representation for the electronic transition amplitudes. In general, this cannot be performed exactly within a classical-trajectory picture, and hence, different mapping-based classical-trajectory techniques differ by how they approximate the electronic transition amplitudes. We first consider the electronic transition amplitudes within the MMST mapping representation and then how these quantities are approximated within the standard partially linearized density matrix (PLDM) approach. We now set $x¯→x$ and $p¯→p$ throughout the rest of this paper.

The standard partially linearized density matrix (PLDM)¹³ approach represents the electronic transition amplitudes using the Meyer–Miller–Stock–Thoss mapping.^31,32 Within this MMST mapping, the electronic basis states, |λ⟩, are described by the single phonon excitation subspace of a set of F harmonic oscillators, where F is the size of the electronic system.³² Within the standard PLDM/FBTS derivation, the electronic transition amplitudes are approximated within this MMST mapping space as follows. First, the coupling propagator is evaluated using a harmonic oscillator coherent state basis,¹⁹ 

where $Zk={Z1(k),Z2(k),…,ZF(k)}$ are the Cartesian mapping variables for the electronic subsystem and $Vm(Z,x)$ is the (scalar) mapping representation of the traceless potential operator,

Equation (9) is exact when the N → ∞ limit of the path-integral expression is taken. Second, the action of the electronic propagator on a harmonic oscillator coherent state is exactly described by evolving the Cartesian mapping variables under the following equations of motion:

Together, Eqs. (9) and (11) lead to an expression for the electronic transition amplitudes [Eq. (8)] that contain an integral over Cartesian mapping variables at each time step. In order to describe the electronic transition amplitudes approximately using continuous trajectories, both standard PLDM and FBTS implement the following approximation involving the overlap of harmonic oscillator coherent states:¹⁹ 

The standard PLDM/FBTS correlation function is then obtained by inserting the resulting approximate expressions for the electronic transition amplitudes into Eq. (4) and then performing the integrals over the nuclear difference coordinates, Δx_k and Δp_k, analytically,^13,19

In this expression, $ϕ(Z)$ is the normalized Gaussian distribution used to sample the Cartesian mapping variables,

and the electronic operators within the real-time quantum correlation function are now described in standard PLDM in terms of these mapping variables as

The dynamics of the mapping coherent states and the dynamics of the nuclear phase-space variables can be described by evolving the classical nuclear and electronic degrees of freedom under the following equations of motion:^13,19

where F₀(x) is the state-independent nuclear force, given by Eq. (7), and $Fe(Z,Z′,x)$ is the state-dependent nuclear force,

The discretized PLDM equations of motion are illustrated schematically in Fig. 2 for time steps of size ϵ. The electronic mapping variables for the forward and backward paths (⁠$Z$ and $Z′$⁠, respectively) are not directly coupled but instead couple via the nuclear degrees of freedom, x, through $Fe(Z,Z′,x)$⁠. This coupling within the equations of motion means that the back action of the electronic subsystem on the nuclear environment is included, at least within the approximations of the method.

The equations of motion for the standard PLDM approach previously derived by Huo and Coker¹³ differ slightly from the expression in Eq. (16). This is because within the standard PLDM approach, the potential is not separated into a purely nuclear potential, V₀(x), and an electronic-dependent potential, $V^(x)$⁠. This results in the whole nuclear force explicitly depending on the mapping variables $Z$ and $Z′$⁠. The forward–backward trajectory solution (FBTS),^19,20 which does separate the nuclear potential into V₀(x) and $V^(x)$ components, has identical equations of motion as Eq. (16) and also has the same expression for the real-time correlation function as Eq. (13). FBTS is however derived as an approximate solution to the QCLE⁶⁶ and not from a path-integral representation of the real-time correlation function. Within the spin-mapping approach (introduced in Sec. II B), the separation of these two nuclear potential components arises naturally from the underlying theory, and hence, it is no longer an arbitrary choice that one must make during the derivation of mapping-based methods.⁶⁹ 

As the QCLE is derived by taking a partial Wigner transform of the exact time-dependent Schrödinger equation, it describes the exact dynamics of electron–nuclear systems in the classical nuclear limit. In order to get a single trajectory picture of the electronic transition amplitudes, the derivation of both standard PLDM and FBTS applies the single approximation given by Eq. (12) for the overlap of harmonic oscillator coherent states at time steps k and k + 1 to an otherwise exact mapping variable solution of the QCLE. Hence, in principle, a more accurate solution of the QCLE dynamics could be found by improving upon the approximation given by Eq. (12). To do so, we derive spin-PLDM within the spin-mapping space.

The main idea behind spin-mapping is the representation of the electronic states of a system by a set of spins, described by classical vectors with fixed radii. This can be achieved in several ways, leading to different properties of the resulting dynamics. The most analogous form of spin-mapping to that of MMST mapping is to represent an F electronic state system by F spin-$12$ particles.³⁸ Like MMST mapping, this spin-mapping is mathematically described as a set of F independent degrees of freedom, except that each degree of freedom can now only contain 0 or 1 excitation, corresponding to the two quantum states of each spin-$12$ particle. While the size of the unphysical subspace available to this spin-mapping is clearly smaller than that of the MMST mapping space, it is still larger than necessary, and more importantly, it also has the unfavorable feature that the classical dynamics are not exact for the pure electronic subsystem.⁴⁶ This probably explains why this particular form of spin-mapping has been found to give rise to less accurate dynamical observables compared to MMST mapping for many systems.³⁸ 

A seemingly more successful version of spin-mapping involves parameterizing the probability amplitudes associated with the electronic basis in terms of spin angles. Now, the associated classical dynamics, which take the form of spins precessing around an effective magnetic field, exactly describe the pure electronic subsystem dynamics, and the spin-mapping space now contains no unphysical regions. While this form of spin-mapping has been around since the advent of mapping-based techniques,^31,70 the underlying equations have recently been reformulated using a Stratonovich–Weyl approach, which leads to a more rigorous way of determining the optimum spin radius (or equivalently the zero-point energy parameter),⁴⁷ different expressions for the measurement operators relative to MMST mapping, as well as a generalization to systems with more than two levels.⁴⁸ This reformulation has already been shown to produce more accurate observables compared to MMST mapping when applied to fully linearized techniques. The multi-state generalization of this spin mapping was formulated in terms of Cartesian mapping variables instead of spin angles, and hence, it is this approach that we will use to derive a spin-mapping version of PLDM.

Within this reformulation of spin-mapping, it was also shown that the theory can be equally implemented on many different spin spheres. For the fully linearized spin methods,^47,48 it was found that the W-sphere consistently gave the most accurate results. This is perhaps because the Stratonovich–Weyl W-functions are self-dual⁴⁷ and is also consistent with previous work that found that using the zero-point energy parameter for the W-sphere is in some sense optimum for MMST classical trajectories.^35,36,41,71 Preliminary tests performed by us have also shown that PLDM is systematically more accurate on the W-sphere than on the Q- and P-spheres. Hence, we will only consider PLDM on the W-sphere, which we will refer to as spin-PLDM.

The Stratonovich–Weyl transform⁷² can be thought of as a discrete version of the Wigner transform within the spin-mapping space. Within this formalism, operators are represented by their Stratonovich–Weyl W-functions, which are given as follows in terms of the Cartesian mapping variables, $Z$⁠,

where $ŵW(Z)$ is the Stratonovich–Weyl W-kernel, which is defined as

and γ_W is the zero-point energy parameter for the W-sphere,⁴⁸ 

Like for the Wigner transform, the trace of a product of two operators can be exactly written as an integral over the product of their corresponding Stratonovich–Weyl W-functions: $tr[ÂB^]=∫dZρW(Z)AW(Z)BW(Z)$⁠. Note that only the expression for the zero-point energy parameter given by Eq. (20) results in this property of the Stratonovich–Weyl W-functions being satisfied. Within this relation, $AW(Z)$ and $BW(Z)$ are the W-functions of their associated electronic operator, given by Eq. (18), and $ρW(Z)$ is a distribution that confines the mapping variables to lie on the surface of the W-sphere of radius R_W,

where the W-sphere radius, R_W, is given generally for a F-level electronic system as^48,73

Additionally, a factor of F is included in this definition of the Cartesian mapping variable distribution [Eq. (21)] so that $∫dZρW(Z)=tr[I^]=F$⁠.

The Stratonovich–Weyl W-kernel can also be used to represent electronic operators, namely,

where $Vm(Z,x)$ is the MMST mapping representation of the traceless electronic Hamiltonian $V^(x)$ [Eq. (10)], which can be shown to be identical to the Stratonovich–Weyl W-function for the same operator from Eqs. (18) and (19). Hence, the coupling propagator, which appears in the definition of the electronic transition amplitudes, can also therefore be obtained in terms of the Stratonovich–Weyl W-kernel as

This result can be confirmed by comparing a Taylor expansion of the exponential on both sides of Eq. (24) and using the properties of the Stratonovich–Weyl kernel given by Eq. (23). As the error in this expression is of order ϵ², the Stratonovich–Weyl representation of the coupling propagator is hence valid when the N → ∞ limit is taken. This step is analogous to representing the coupling operator in terms of harmonic oscillator coherent states within the derivation of standard PLDM/FBTS [Eq. (9)].

Inserting Eq. (24) into Eq. (8a) results in a Stratonovich–Weyl kernel appearing to the right of every electronic propagator [except for the $e−iV^(x1)ϵ$ operator, which we post-multiply by $I^=∫dZ0ρW(Z0)ŵW(Z0)$]. We then define the time-evolved Stratonovich–Weyl kernel, $ŵW(Zk,ϵ)$⁠, as follows:

Unlike for the harmonic oscillator coherent states used in standard PLDM, the time-evolved Stratonovich–Weyl kernel cannot be obtained purely by evolving the Cartesian mapping variables, $Zk$⁠. This is because the Stratonovich–Weyl kernel contains a zero-point energy parameter, γ_W, multiplied by an identity operator, which must also be evolved forward in time. Hence, the time-evolved Stratonovich–Weyl kernel is given by

where $Zk(ϵ)={Z1(k)(ϵ),Z2(k)(ϵ),…,ZF(k)(ϵ)}$ are the Cartesian mapping variables at time ϵ, evolved according to Eq. (16). These equations of motion for the Cartesian mapping variables exactly describe the pure electronic sub-system dynamics for any size time step ϵ. Additionally, $e−iV^(xk+1)ϵ$ is an F × F matrix in the electronic-state basis, which can easily be computed numerically. For the backward (⁠$Z′$⁠) electronic transition amplitude, inserting Eq. (24) into Eq. (8b) results in a Stratonovich–Weyl kernel appearing to the left of every electronic propagator, $eiV^(xk)ϵ$⁠. Hence, the Hermitian conjugate of Eq. (25) is used in this case.

Therefore, Eqs. (24) and (25) lead to an expression for the electronic transition amplitudes in terms of the Cartesian mapping variables. However, the expression contains an integral over mapping variables at each time step, k. Following the standard PLDM/FBTS procedure, in analogy to Eq. (12), the overlap of Stratonovich–Weyl kernels at different time steps is approximated as follows to generate a trajectory picture for the electronic dynamics:

where the Dirac delta function of a complex argument, $Z$⁠, is defined as $δ(Z)=δ(Re[Z])δ(Im[Z])$⁠. Employing the approximation given by Eq. (27) reduces the expression for the electronic transition amplitudes to

which is exact for a purely electronic system and where we now set $Z0=Z$ and $Z0′=Z′$⁠. The electronic action for the forward and backward paths is defined as

where t_k = kϵ is the time at time step k. Additionally, $ŵW(Z,t)$ is the time-evolved Stratonovich–Weyl kernel, which is defined as

where $Û(t)$ is a time-ordered propagator for the electronic states according to the time-dependent potential obtained by moving along a given path x(t_k) = x_k,

Assuming that F is not too large, $Û(t)$ is easy to evaluate as it is just the dynamics of the bare electronic space according to a time-dependent Hamiltonian. The back action does not appear directly in the time-ordered propagator but is treated by the evolution of the Cartesian mapping variables.

To complete the linearization approximation of the nuclear path, outlined in Sec. II, the integrals over Δx_k and Δp_k must be performed. Using the approximate expressions for the nuclear and electronic action, given by Eqs. (6) and (29), these integrals can be performed analytically in a similar fashion to the standard PLDM approach¹³ to give

The result is a product of Dirac delta functions that constrain the dynamics of the nuclear variables to follow Newton’s equation of motion. The equations of motion for the classical variables within spin-PLDM are identical to those of standard PLDM [Eq. (16)], although the correlation functions will of course still give different results due to the change in the initial distribution of the Cartesian mapping variables.

This means that the approximate expression for the spin-PLDM real-time correlation function becomes

where we now set x₀ = x and p₀ = p. The trace in Eq. (33) only contains F terms and is hence easy to perform explicitly for the majority of systems of interest. This expression for the spin-PLDM correlation function looks quite different compared to that of other mapping-based approaches, particularly compared to fully linearized methods. In Paper II,⁶¹ we rewrite the correlation function given by Eq. (33) into a form that allows direct comparison with other mapping-based methods. While this reformulation is not practical for use in numerical simulations, it shows that the spin-PLDM correlation function contains additional terms that are absent within standard PLDM and fully linearized methods, which are shown numerically to lead to improved accuracy. Therefore, this analysis also justifies why such a partially linearized approach can be seen as advantageous compared to a fully linearized one.

The spin-PLDM algorithm can be easily implemented as follows. First, the integrals contained in the correlation function given by Eq. (33) can be evaluated numerically by sampling the initial classical coordinates using Monte Carlo. This means sampling the nuclear phase-space variables from the Wigner distribution of the nuclear initial density [ρ_b(x, p) given by Eq. (5)] and sampling the Cartesian mapping variables for the forward and backward paths (⁠$Z$ and $Z′$⁠) independently from uniform hyperspheres of radius R_W. For each instance of the sampling, the classical coordinates (x, p, $Z$⁠, and $Z′$⁠) can then be evolved in time using the equations of motion given by Eq. (16), and additionally, the time-ordered propagator, $Û(t)$⁠, can be obtained by using Eq. (31), where ϵ is the time step. Finally, the W-kernel for both the forward and backward paths [$ŵW(Z,t)$ and $ŵW(Z′,t)$] can then be constructed using Eq. (30), and the contribution to the real-time correlation function for each sample can be obtained by explicitly evaluating the matrix multiplications and trace in the electronic basis. Such an algorithm can also be implemented with only a few minor changes into a preexisting standard PLDM code.

The main differences between spin-PLDM and standard PLDM [whose correlation function is given by Eq. (13)] are as follows. In spin-PLDM, the initial mapping variables are constrained to a hypersphere $|Z|2=RW2$⁠, which guarantees that the electronic state being represented is correctly normalized. Additionally, spin-PLDM has a zero-point energy parameter, which means that it treats correlation functions containing the identity operator slightly differently from those of traceless operators. These two changes are expected to lead to improved results in a similar way in which fully linearized spin-mapping^47,48 demonstrates improved accuracy over fully linearized MMST mapping-based approaches.

In this section, we consider a set of challenging systems far from the Born–Oppenheimer limit. In order to improve the convergence of the spin-PLDM method with respect to the number of trajectories, we have obtained the following results using so-called focused initial conditions to sample the Cartesian mapping variables at t = 0. The details concerning the implementation of these focused initial conditions for spin-PLDM can be found in Paper II⁶¹ along with numerical comparisons of the two approaches. In contrast to MMST mapping-based methods, we have found that the spin-PLDM results obtained using focused initial conditions are essentially indistinguishable from those obtained using the original sampling of the Cartesian mapping variables, as outlined in Sec. II B, although the latter require more trajectories. For the MMST mapping-based techniques, only the most accurate form of the methods are used, so as to act as a fair comparison with spin-mapping; hence, we present results only for the non-focused variants of these methods.

To test the accuracy of our spin-PLDM method, we have applied the method to a series of spin-boson models. The parameters associated with these spin-boson models are chosen in order to cover a large part of the parameter space, including both symmetric and asymmetric systems, low and high temperature limits, and strong and weak system–bath coupling. Such spin-boson models are commonly used to benchmark new methods as numerically exact results are available for comparison. The model essentially consists of two electronic states coupled to an initially thermalized harmonic bath, whose Hamiltonian is given by⁷⁴ 

Here, ε is the energy bias and Δ is the constant diabatic coupling. The bath contains f nuclear modes, each with frequency ω_j and electron–nuclear coupling coefficient c_j. The mass, m_j, has been incorporated into the definition of the nuclear position, $x^j$⁠, and momentum, $p^j$⁠, operators. As the nuclear bath is harmonic, the dynamics of the spin-boson model is exactly described by the quantum-classical Liouville equation (QCLE).^68,75 Hence, any errors arising from calculating correlation functions using mapping-based classical-trajectory techniques arise solely from approximations to the QCLE and are not inherently due to the classical treatment of the nuclear degrees of freedom. The spectral density, J_bath(ω), determines the distribution of nuclear frequencies within the bath. One of the most commonly used spectral densities is the Ohmic bath form

where ω_c is the characteristic frequency and ξ is the Kondo parameter. In order to perform numerical simulations, the continuous bath must first be discretized into a finite number of modes. The spectral density can then be represented as follows:

We have used the discretization scheme employed in Ref. 76. Additionally, we consider solely correlation functions in which the uncoupled bath is initially thermalized. This means that for mapping-based classical-trajectory techniques, the initial nuclear coordinates are sampled from the following Wigner distribution:

where $αj=tanh(12βωj)$ and β is the inverse temperature. For comparison, numerically exact results for the real-time quantum correlation functions have been obtained using the quasiadiabatic path-integral (QUAPI) technique.⁷⁷ 

The parameters for the spin-boson models we consider are given by rows (a)–(f) in Table I. For these calculations, we have used f = 100 nuclear degrees of freedom, except for the strong coupling systems (e) and (f), where f = 400 nuclear degrees of freedom were needed for convergence. Such a set of spin-boson models covers a wide range of different physical regimes, from symmetric to asymmetric, weak to strong system–bath coupling, and low to high temperature limits. Hence, these models offer a comprehensive test for our new spin-PLDM method. Additionally, all of these spin-boson models lie in a challenging regime for classical-trajectory methods to describe correctly as they are far from the Born–Oppenheimer limit and have already been used to test other methods, which allows us to directly compare the accuracy of various methods.^25,41,47,78 While Ehrenfest and standard linearized mapping methods are known to perform relatively well for the symmetric models, these methods fail to correctly predict the long-time populations in the asymmetric models due to a failure to obey detailed balance. As in Ref. 25, we also make a point of computing the time-dependent coherences, as well as populations. In this paper, Eq. (37) initializes the nuclear coordinates from the Boltzmann distribution of nuclear potential V₀(x). This is subtly different from Refs. 41 and 78, where the nuclei are sampled from the Boltzmann distribution for potential $V0(x)+⟨1|V^(x)|1⟩$⁠. Nonetheless, these different bath initial conditions do not appear to lead to significant differences in the results.

Figure 3 shows the calculated dynamical expectation values of all Pauli spin matrices for a range of mapping-based techniques. All results correspond to the system initially occupying the higher energy diabatic state, |1⟩. The spin linearized semiclassical (spin-LSC) approach is the fully linearized version of spin-PLDM and corresponds to the method described in Ref. 48 for the W-sphere. Additionally, the Poisson-bracket mapping equation (PBME) approach is the fully linearized version of standard PLDM¹⁹ and is described in Ref. 51. While we do not show results for LSC-IVR in this paper, it has been observed that both PBME and LSC-IVR produce real-time correlation functions to a similar degree of accuracy.^23–25 Consider first the weak system–bath coupling spin-boson models at high temperatures [given by (a) and (c)]. We would expect to describe the dynamics of these models correctly as the classical approximation that is employed in deriving mapping-based classical-trajectory techniques should be valid at high temperatures. We find that only PBME is unable to accurately describe the dynamics within these spin-boson models. The error arises because PBME does not correctly calculate the identity-containing correlation functions. Note that this is not a fundamental problem of linearized MMST dynamics as it is known that the error can be almost completely removed by using an “identity-trick” recently derived for PBME^23–25 or by windowing trajectories at time t using the symmetrical quasi-classical (SQC) approach.^12,36,44 Therefore, for this parameter regime, spin-PLDM does not offer any advantage over spin-LSC or standard PLDM, which are essentially numerically exact already for these systems.

Next, we consider the more difficult models with weak system–bath coupling at low temperatures [given by (b) and (d)]. For these models, spin-PLDM appears to consistently produce the most accurate results compared to other mapping-based approaches. First, spin-PLDM outperforms spin-LSC for these systems by correcting the overdamped coherences found in the spin-LSC correlation functions of traceless operators. Additionally, spin-PLDM also outperforms standard PLDM in this parameter regime because restraining the Cartesian mapping variables to a hypersphere reduces the large errors observed in the identity-containing correlation functions. This in turn results in a reduction in the error associated with the long-time limit of any real-time correlation function because all correlation functions containing traceless operators decay to zero as t → ∞. The spin-PLDM initial conditions for the mapping variables also appear to correct for the overdamped coherences found in the standard PLDM correlation functions of traceless operators. While both SQC and the MMST identity-trick can to some extent correct for the errors in the long-time limits of the PBME correlation functions, they typically still predict overdamped coherences similar to spin-LSC.^25,41

Finally, systems (e) and (f) are the most challenging spin-boson models that we consider in this paper due to the strong system–bath coupling. In fact, the coupling in system (f) is so strong that we struggled to converge the numerically exact QUAPI results. Hence, we only present the QUAPI results for times for which we are certain that they are numerically exact. The main conclusion from these figures is that spin-PLDM is the best method at obtaining the short time dynamics correctly and therefore is also expected to be the most trustworthy method for predicting the long-time limit of dynamical quantities. We assume that any method that is less accurate than spin-PLDM at short times but appears more accurate in the long-time limit in isolated cases must be due to a lucky cancellation of errors. The figures show that spin-PLDM provides a larger correction to spin-LSC for the strong-coupling spin-boson models. This improvement even holds for the high-temperature model (f), where spin-LSC was already quite good. This suggests that the explicit treatment of the forward and backward electronic mapping variables in spin-PLDM is particularly important for describing the dynamics in systems with relatively strong system–bath coupling.

To demonstrate that the spin-PLDM method can be easily applied to multistate systems, we have applied the method to an F = 7 model for the FMO complex, which is a commonly studied photosynthetic pigment–protein complex found in green sulfur bacteria.⁷⁹ This is a challenging benchmark problem that has already been used to test a wide range of mapping-based techniques.^{13,15,20,24,40,44,48,53,54,56,58,80–85} However, as for the spin-boson model, the dynamics of this FMO model is exactly described by the QCLE because the nuclear bath is harmonic, which means that nuclear quantum effects do not need to be included in order to describe the quantum dynamics of this model exactly.^3,68,75 Numerically exact results have been computed for this model from the hierarchical equation of motion (HEOM) approach.^86–91 

The Hamiltonian of this model is⁸⁶ 

where the electronic system Hamiltonian is given, in units of cm⁻¹, as⁹² 

Each diabatic electronic state is coupled to its own independent harmonic bath that consists of f nuclear modes with frequencies ω_j. The mass-weighted Hamiltonian for all seven independent baths is therefore given by

where each bath has the same set of frequencies, ω_j. The model then contains a linear coupling term, which connects each harmonic bath to its corresponding diabatic state as follows:

where c_j are the exciton–nuclear coupling coefficients.

The spectral density, J_bath(ω), determines the distribution of nuclear frequencies within each of the baths and their couplings. This model employs the Debye spectral density

where ω_c is the characteristic frequency of the bath (with $τc=ωc−1$ being its corresponding time scale) and Λ is the reorganization energy. We use Λ = 35 cm⁻¹, to be consistent with previous work, and only consider the situation of a “fast bath” (i.e., τ_c = 50 fs), which is the most difficult previously studied parameter regime. In order to perform numerical simulations, the continuous bath must be first discretized into a finite number of modes, as shown in Eq. (36). We have used the discretization scheme employed in Ref. 93 with f = 60 bath modes per state. We consider the dynamics after the initial excitation of a single site, with the baths in thermal equilibrium before the excitation. This means that the initial nuclear coordinates associated with each excitonic site are sampled from the same Wigner distribution given by Eq. (37).

To test the accuracy of our newly developed spin-PLDM method, we calculate the time-dependence of the exciton dynamics for this FMO model. Of particular interest is the picosecond exciton population transfer to the lowest energy site/chromophore (n = 3), from which the excitons can be harvested at the reaction center. We consider both the dynamics at high and low temperatures (300 K and 77 K, respectively) and also that with the initial exciton residing on different sites/chromophores of the complex (i.e., sites 1 and 6). Figure 4 shows the population dynamics of the FMO at high temperature (300 K). As stated before, this is the regime in which the approximations involved in deriving mapping-based classical-trajectory techniques are expected to be valid. From Fig. 4, it can be seen that spin-LSC, standard PLDM, and spin-PLDM are each able to describe the short time dynamics very well (i.e., up to 1 ps), irrespective of whether the exciton is initialized on site 1 (the top row of figures) or site 6 (the bottom row of figures). Even for this short-time behavior, the accuracy of the spin-PLDM approach, however, appears superior as the calculated dynamics from this technique are essentially indistinguishable from the exact HEOM results (solid lines) for this high-temperature parameter regime. The accuracy of spin-PLDM continues to be good, even when the long-time dynamics (up to 10 ps) are considered. The spin-LSC method also gives an excellent description of the long-time limits of the exciton populations. However, standard PLDM is unable to describe the long-time dynamics correctly for this high-temperature regime and exhibits a large deviation from the exact result, in particular, for the site 3 population. This illustrates a clear advantage in using spin-mapping instead of the standard MMST approach. Constraining the mapping variables to a hypersphere ensures that these Cartesian mapping variables remain in the physical subspace during the dynamics, and this appears to minimize the errors produced in the long-time limit and also allows the method to better approximate the correct Boltzmann distribution. While both the SQC⁴⁴ and MMST identity-trick²⁴ results for this model are clearly superior to those for PBME,²⁴ both these methods show greater errors in the dynamical populations compared to spin-PLDM, even when just considering the evolution up to 1 ps.

The more challenging regime corresponds to studying the dynamics of the model at a low temperature (77 K). Figure 5 compares the calculated population dynamics for the same mapping-based classical-trajectory techniques. In the short-time dynamics, spin-LSC now shows noticeable errors compared to the exact HEOM results. Additionally, some of the calculated populations with spin-LSC now become negative, even at quite short times (see, for example, the n = 6 population, when starting in site 1), which is unphysical. While the standard PLDM results appear fairly accurate at short times, the method suffers from overdamped coherences, similar to what was observed when applying the method to the spin-boson model. Spin-PLDM produces almost exact results for this short-time dynamics, which again illustrates its apparent increased accuracy, not only over spin-LSC and PLDM but also when compared with the MMST identity-trick²⁴ and SQC.⁴⁴ In the long-time limit (i.e., up to 10 ps), the standard PLDM approach again leads to large errors in the populations. In contrast, spin-LSC appears to produce less severe errors in the populations in the long-time limit. In particular, the n = 3 population is well reproduced by the method. This is however probably partly fortuitous due to the method’s inability to describe the short-time dynamics of this population correctly. The errors within spin-PLDM in the long-time limit also appear not too severe, and the method also suffers less severely with negative populations compared to spin-LSC. Standard PLDM, in contrast, always calculates positive populations because the theory has no zero-point energy (γ) parameter. Although the difference is not as great as for the 300 K case, nonetheless, the absolute errors in the populations are larger than those for spin-PLDM, which again shows itself to be the most accurate of the three methods.

In this paper, we have derived a new partially linearized mapping approach based on classical trajectories, which uses elements of previously derived partially linearized techniques (such as standard PLDM and FBTS) and applies them to the spin-mapping space using the Stratonovich–Weyl transform. We have then tested this method, called spin-PLDM, on a range of model systems in various different regimes. A comparison has also been made with other mapping-based classical-trajectory techniques.

We show that this method exhibits improved accuracy over standard PLDM, when calculating correlation functions. In particular, spin-PLDM corrects for the large error observed in the long-time limit of identity-containing correlation functions calculated using standard PLDM and also appears to solve the issue of overdamped oscillations observed in many dynamical quantities. This suggests that spin-PLDM is in some sense closer to an exact solution of the QCLE than the standard PLDM method. Spin-PLDM also appears to be able to improve upon previously derived fully linearized approaches, such as spin-LSC. In particular, we have observed that spin-PLDM appears to solve the problem of “overdamped coherences” observed in the correlation functions of traceless operators calculated using spin-LSC.

For the FMO model, 10⁷ trajectories were used to fully converge the focused spin-PLDM results compared to 10⁶ trajectories for the other methods; however, we also observed that only 10⁴ trajectories were needed to qualitatively reproduce the main features of the dynamics. Although we found that focused spin-PLDM generally requires one order of magnitude more trajectories than other mapping-based classical-trajectory methods, we expect that more efficient sampling schemes could easily be implemented in the future in order to reduce the number of required trajectories.⁹⁴ In general, our results show that generalized spin-mapping methods outperform their MMST analogs and partially linearized methods appear superior to fully linearized ones; hence, spin-PLDM is the best mapping method of them all. The superior accuracy of spin-PLDM is evident not just for calculating electronic population dynamics but also for coherences.

In particular, spin-PLDM seems to reproduce the relatively short time properties of real-time quantum correlation functions extremely accurately compared to other mapping-based classical-trajectory techniques. Because of this, there are a number of applications that may therefore be well suited for spin-PLDM. For instance, memory kernels within the generalized quantum master equation formalism normally decay relatively rapidly, and hence it is expected that such quantities can be accurately computed using spin-PLDM, as has been done with other methods.^83,95–98 Dipole–dipole correlation functions and other optical response functions from which linear and non-linear spectra can be obtained also often have short coherence times and, hence, perhaps can also be accurately computed using spin-PLDM.

In Paper II,⁶¹ the spin-PLDM method is extensively analyzed in order to understand the key differences between spin-PLDM and standard PLDM and between spin-PLDM and spin-LSC. We also outline the implementation of focused initial conditions for spin-PLDM, which were used to generate the results of this work, and introduce a jump spin-PLDM scheme, which, in principle, allows results to be systematically improved to those of QCLE.

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

The data that support the findings of this study are available within the article.