Fewest-switches surface hopping (FSSH) has emerged as one of the leading methods for modeling the quantum dynamics of molecular systems. While its original formulation was limited to adiabatic populations, the growing interest in the application of FSSH to coherent phenomena prompts the question of how one should construct a complete density matrix based on FSSH trajectories. A straightforward solution is to define adiabatic coherences based on wavefunction coefficients. In this paper, we demonstrate that inconsistencies introduced in the density matrix through such treatment may lead to a violation of positivity. We furthermore show that a recently proposed coherent generalization of FSSH results in density matrices that satisfy positivity while yielding improved accuracy throughout much (but not all) of parameter space.

## I. INTRODUCTION

Many dynamical phenomena in chemistry and physics are effectively described by partitioning the problem of interest into a system and its environment while limiting an explicit quantum treatment to the system in the form of a reduced density matrix. Given that appropriate approximations are made, this enables the accurate computational modeling of phenomena that are otherwise intractable. The multitude of approximations that can be taken has led to a wide variety of quantum dynamical methods that have been proposed over the years. When invoking approximations, in addition to the sheer accuracy of the method, it is of particular importance that the physical properties of the (reduced) density matrix are conserved. A well-known example of a violation of such physical properties is the breaking of positivity of the density matrix under the Markov approximation.^{1–3} Both the Nakajima–Zwanzig and Redfield master equations violate positivity in principle,^{4} with manifestations of this being particularly well documented for the latter.^{5–7} Within Redfield theory, positivity can be preserved by additionally invoking the secular approximation,^{8–10} albeit at the cost of additional approximations.

In 1990, Tully introduced fewest-switches surface hopping (FSSH),^{11} which has grown into one of the most widely used quantum dynamical methods.^{12–15} FSSH is a mixed quantum–classical technique that represents the quantum system by means of a swarm of “active” surfaces that are allowed to switch between adiabats (instantaneous eigenstates of the quantum Hamiltonian) as a result of interactions with a classical environment. Intended for the modeling of (incoherent) scattering phenomena,^{11,16} FSSH in its original form exclusively specified adiabatic populations, i.e., diagonal elements of the density matrix. However, a growing interest in the application of FSSH to coherent phenomena^{17–28} prompted the question of how to construct the *entire* density matrix based on FSSH. Arguably the most straightforward and common approach to constructing coherent (off-diagonal) density matrix elements is by invoking the wavefunction coefficients that are propagated through the time-dependent Schrödinger equation and which govern the hopping probabilities.^{17,29} Although it was immediately recognized that this approach (henceforth simply referred to as FSSH) yields inconsistencies between populations and coherences, only recently has a systematic survey of the spin-boson model shown that it introduces inaccuracies in some cases, although no unphysical behavior was observed.^{30}

The purpose of this paper is twofold. First, it demonstrates that the aforementioned inconsistencies between populations and coherences may yield unphysical density matrices. This is exemplified in Fig. 1, where a violation of positivity is observed in the local (diabatic) basis for a trimer (a three level system) coupled to a harmonic environment. These observations underscore the importance of recent efforts aimed at finding consistent FSSH-like formulations of the entire density matrix,^{30–35} leading up to the second purpose of this paper of singling out a recently proposed generalization of FSSH, called coherent fewest-switches surface hopping (C-FSSH).^{30} C-FSSH uniquely provides a consistent formulation of the entire density matrix based on the existing (and practically unaltered) FSSH method by invoking pairs of active surfaces. Meanwhile, the classical coordinates are represented by Gaussian wavepackets, allowing zero-point fluctuations to be represented without invoking a Wigner distribution. We demonstrate that not only does C-FSSH resolve the positivity violations of FSSH, but it also produces comparatively superior results throughout much (but not all) of parameter space for a trimeric system.

The remainder of this paper is organized as follows: In Sec. II, we provide a summary of the physical properties of density matrices. In Sec. III A, we describe some general concepts of mixed quantum–classical dynamics, while in Secs. III B and III C, we introduce method details specific to FSSH and C-FSSH, respectively. In Sec. IV A, we introduce the model system used in this study, after which we discuss the results in Sec. IV B. Finally, we present our conclusions and outlook in Sec. V.

## II. PHYSICAL PROPERTIES OF DENSITY MATRICES

Before evaluating density matrices within the FSSH and C-FSSH formalisms, we begin by iterating the physical properties a valid (reduced) density matrix should satisfy. First, it should be noted that the density matrix of a pure state in some basis is given by $\rho nm=cn*cm$, with *c*_{n} as the wavefunction coefficient in the basis at hand (for which basis states are labeled by *n*). For a mixed state, the density matrix instead takes the form $\rho nm=\u27e8cn\u2020cm\u27e9$, where ⟨⋯⟩ denotes an average over pure states. The physical properties emanating from these principles are as follows:

Hermitian property: The density matrix should be Hermitian,

*ρ*^{†}=*ρ*. It should be noted that the Hermitian property survives a unitary transformation of*ρ*, meaning it carries over from one basis to another.Normalization: Populations represent probabilities, and therefore their sum should add up to one. As a result, the trace of the density matrix should amount to unity, tr(

*ρ*) = 1. Like the Hermitian property, normalization survives a unitary transformation of*ρ*.Positivity: Positivity of a density matrix can be evaluated in various forms, as outlined in the following.

- (3a)
Positive populations: As populations represent probabilities, the values they assume should be strictly positive,

*ρ*_{nn}≥ 0. However, this positivity criterion does not necessarily survive a unitary transformation, meaning that positivity in one basis does not imply positivity in another basis. This is the result of a mixing of populations and coherences by the basis transformation. Hence, positivity in*all possible bases*requires these elements of*ρ*to be consistent, leading up to the next two properties. - (3b)Cauchy–Schwartz (CS) inequality: The CS inequality puts bounds on the off-diagonal elements of the density matrix in terms of the diagonal elements, asThe CS inequality is saturated for pure states, i.e.,(1)$\rho nn\rho mm\u2265\rho nm\rho mn.$
*ρ*_{nn}*ρ*_{mm}=*ρ*_{nm}*ρ*_{mn}, which is the result of the diagonal and off-diagonal elements being interrelated through $\rho nm=cn*cm$. For mixed states, on the other hand, it takes the form*ρ*_{nn}*ρ*_{mm}≥*ρ*_{nm}*ρ*_{mn}, since taking the average ⟨⋯⟩ may yield destructive interference for off-diagonal elements but not for diagonal elements. When a 2 × 2 density matrix satisfies the CS inequality, positivity is ensured in*any*basis, rendering the CS inequality a more complete requirement than positivity in a given basis (Property 3a). However, density matrices of size 3 × 3 and larger involve multiple CS inequalities (invoking pairs of basis states) that become basis dependent. The importance of having a basis-independent metric of positivity leads to the next and last property. - (3c)Positive eigenvalues: Eigenvalues of a matrix represent the extrema of diagonal elements as a function of all possible unitary rotation angles. Hence, the extrema for the density matrix can be found by introducing the density operator,and solving the eigenvalue equation $\rho \u0302|\alpha \u3009=r\alpha |\alpha \u3009$. Positivity of the density matrix then implies that(2)$\rho \u0302=\u2211m,n\rho m,n|m\u3009\u3008n|,$
*r*_{α}≥ 0 for all*α*. This forms an unambiguous and basis-independent definition of positivity, regardless of the size of the density matrix.

- (3a)

## III. THEORY AND METHODS

Having summarized the physical properties a density matrix should satisfy, we now proceed to evaluate density matrices within the FSSH and C-FSSH formalisms, after first introducing the general concept of mixed quantum–classical dynamics.

### A. Mixed quantum–classical dynamics

Mixed quantum–classical dynamics approximates the environment as classical, yielding a total Hamiltonian given by

where ** q** and

**are the classical position and momentum vectors, respectively. In Eq. (3), we differentiate contributions from the quantum system (q) and the classical environment (c), as well as their mutual interaction (q–c), the latter two depending parametrically on the classical coordinates. From here onward, this dependence will be omitted in order to simplify the notation.**

*p*A mixed quantum–classical simulation amounts to propagating a swarm of classical trajectories, while quantities of interest are evaluated as a trajectory average. (For simplicity, we will use the same notation for a trajectory average as for an average over pure states, i.e., ⟨⋯⟩, even though these two averages are not strictly interchangeable.) Within each trajectory, the classical coordinates evolve according to Hamilton’s equations using the Hamiltonian $H\u0302q\u2212c(q)+Hc(p,q)$. Appropriately taking the expectation value of $H\u0302q\u2212c(q)$ in order to arrive at a meaningful classical potential energy is a source of ambiguity, leading to the many variants by means of which mixed quantum–classical dynamics can be implemented.

In concert with the classical dynamics, the quantum wavefunction is propagated through the time-dependent Schrödinger equation (setting ℏ = 1),

At any given time, the quantum wavefunction can be expanded in the adiabatic basis as

with the adiabatic basis states following from the instantaneous eigensolution of the time-independent Schrödinger equation,

Here, ϵ_{k} is the associated (adiabatic) energy.

While it may serve to intuitively assess quantities of interest, the adiabatic basis strictly does not allow for a trajectory average to be taken, as the basis itself is trajectory dependent. Rather, a diabatic (trajectory independent) basis needs to be adopted, such as the local basis, within which the adiabatic states are expanded as

where the unitary transformation *U*_{k,n} itself is trajectory dependent. Here and henceforth, we use *k* and *l* to denote adiabatic states and *n* and *m* to denote local basis states. Moreover, as much as possible, we will use the local basis to assess whether a given density matrix satisfies the physical properties summarized in Sec. II.

### B. Fewest-switches surface hopping (FSSH)

We proceed with a quick summary of FSSH, referring to earlier literature for a more elaborate explanation.^{16} FSSH is an example of a broader class of surface hopping methods, wherein the classical potential energy is derived from an expectation value of the quantum–classical Hamiltonian with respect to a single adiabatic surface, i.e., $\u27e8a|H\u0302q\u2212c|a\u27e9$. FSSH allows for stochastic switches between surfaces to happen at any time, with the probability of switching from surface *k* to *l* being governed by the probability

Here, *c*_{k} is the wavefunction expansion coefficient from Eq. (5), Δ*t* is the time integration step, and *d*_{k,l} is the nonadiabatic coupling vector defined as

In order to limit the number of switches, a uniform random number *ξ* ∈ (0, 1) is drawn at each time step, and the switch from *k* to *l* is allowed to occur only when

Moreover, in order to conserve (quantum plus classical) energy, the change in quantum energy upon a switch, *ɛ*_{k} − *ɛ*_{l}, needs to be absorbed or donated in the form of an adjustment of the classical kinetic energy, *p*^{2}/2. Accordingly, the classical momentum is rescaled in the direction of the nonadiabatic coupling vector. Whenever there is insufficient classical kinetic energy available, the switching event is abandoned.

As mentioned in the introduction, FSSH was originally developed with the intention to study scattering phenomena. In such cases, the active surface can be initiated in some adiabatic state as *a* = *i,* and scattering probabilities can be assessed by taking a trajectory average of *a* at a later time. As pointed out in Sec. III A, such a trajectory average within the adiabatic basis is not strictly allowed, yet it provides a good approximation in the asymptotic regimes of scattering phenomena where the adiabatic basis becomes near-diabatic. A more rigorous approach to quantifying scattering outcomes, however, is to consider the density matrix in the local basis,^{29,36,37}

with the trajectory-dependent adiabatic contributions given by^{29,36,37}

In order to assess whether *ρ*_{n,m} satisfies the physical properties summarized in Sec. II, it is worth pointing out that the adiabatic contributions given by Eq. (12) each trivially satisfy these properties, which carries over to *ρ*_{n,m}.

Where Eqs. (11) and (12) fall short is in their inability to account for a sharing of the initial quantum population among different surfaces. This can be straightforwardly incorporated by taking $\rho n,m=\u27e8\u2211i\rho n,m(i)\u27e9$, with $\rho n,m(i)=\u2211k,lUk,n(i)\rho k,l(i)Ul,m(i)*$ and with

Here, *i* refers to a “branch” of the trajectory that was initiated at surface *i*, meaning that *a*^{(i)} = *i* initially. In addition to the active surface, this branch comes with its own set of classical coordinates *p*^{(i)} and *q*^{(i)}, as well as wavefunction coefficients $ck(i)$, each initiated using the same (branch-independent) value. Moreover, each branch is weighted by $\rho 0(i)$ in Eq. (13). If $\u2211i\rho 0(i)=1$ is enforced, it follows that ∑_{n}*ρ*_{n,n} = 1 and normalization is satisfied. Likewise, *ρ*_{n,m} can be shown to satisfy the other physical properties summarized in Sec. II.

The adiabatic contribution to the density matrix embodied by Eq. (13) is still diagonal and, as such, completely neglects adiabatic coherences. Such coherences can be incorporated by supplementing Eq. (13) with off-diagonal elements based on the wavefunction coefficients, as^{17,18,21,22,24,29,38}

(This implementation was referred to as “Method 3” in Ref. 29.) It can be easily checked that this supplement does not contribute to ∑_{n}*ρ*_{n,n} so that normalization is satisfied as long as $\u2211i\rho 0(i)=1$. Moreover, if the wavefunction coefficients are initiated such that $\rho 0(i)=|ci|2$, it can be shown that initially $\rho n,m=cn*cm$ as a result of which the physical properties outlined in Sec. II are automatically satisfied. However, as time evolves, the diagonal and off-diagonal elements start to behave markedly different, being constructed out of active surfaces and wavefunction coefficients, respectively. As a result of such inconsistencies, positivity may become violated, which is most intuitively understood as a consequence of the inconsistencies breaking the CS inequality (Property 3b). As a result, negative populations may arise in some basis, as found in Fig. 1.

### C. Coherent fewest-switches surface hopping (C-FSSH)

If one propagates multiple branches for each trajectory as outlined in Sec. III B, one can *in principle* construct off-diagonal elements of the density matrix based on pairs of branches. This is the core idea behind C-FSSH. A flowchart of this approach is provided in the supplementary material, and further details can be found in Ref. 30. Within C-FSSH, the density matrix is constructed as $\rho n,m=\u27e8\u2211i,j\rho n,m(i,j)\u27e9$. Here, $\rho n,m(i,j)=\u2211k,lUk,n(i,j)\rho k,l(i,j)Ul,m(i,j)*$ with

where (*i*, *j*) labels the pair of branches, the contribution of which is weighted by $\rho 0(i,j)$. When this weight is derived from the initial wavefunction coefficients as $\rho 0(i,j)=ci*cj$, it can again be shown that initially $\rho n,m=cn*cm$. Interestingly, the branch-dependent wavefunction coefficients do not need to be uniformly initiated (as $ck(i)=ck$), and it has been found that $ck(i)=\delta k,i$ yields improved accuracy.^{30} In addition to the wavefunction coefficients and the active surface, initiated as *a*^{(i)} = *i*, each branch features a phase factor that is propagated as $\varphi \u0307k(i)=\epsilon k(i)$. Here, $\epsilon k(i)$ is the *k*th eigenvalue of the quantum Hamiltonian $H\u0302q+H\u0302q\u2212c(i)$, which depends parametrically on the branch-dependent classical coordinates *q*^{(i)} through $H\u0302q\u2212c(i)$. One can also define a Hamiltonian for each branch pair, $H\u0302q+H\u0302q\u2212c(i,j)$, from which the unitary transformation $Uk,n(i,j)$ is derived. This Hamiltonian depends parametrically on branch pair-dependent classical coordinates, which can simply be taken to be *q*^{(i,j)} = (*q*^{(i)} + *q*^{(j)})/2.

In Eq. (15), *F*^{(i,j)} represents a classical overlap factor. Accordingly, the classical coordinates are interpreted as describing the centers of wavepackets with a finite spread in position and momentum space. This allows incorporation of the Heisenberg uncertainty principle within the classical environment and, consequently, zero-point fluctuations. When constructing the reduced density matrix of the quantum system by tracing out the classical coordinates, the classical wavepackets contribute an overlap factor of the form^{39}

with

Here, *σ*_{q} represents the width of the classical wavepackets. Upon introducing C-FSSH,^{30} this width was taken to be a constant, consistent with the concept of frozen Gaussians introduced by Heller.^{40} However, whereas Heller’s approach shaped the Gaussian after the original wavepacket, *σ*_{q} was effectively treated as an arbitrary parameter in C-FSSH, following previous examples in the literature.^{41,42} In the following, we demonstrate that this parameter can be eliminated by maximizing the overlap, thus minimizing the argument in Eq. (17) and yielding

The resulting parameter-free treatment of the classical overlap is consistently adopted in the present paper, and based on the favorable results it generates (*vide infra*), we conclude that this treatment forms an improvement over the original C-FSSH method.

As before, in order to assess whether the density matrix satisfies the physical properties summarized in Sec. II, one can examine the adiabatic contributions, which are now trajectory *and branch-pair* dependent. First, it can easily be verified that Eq. (15) by construction satisfies Property 3b, i.e., the CS inequality,

which is a result of the preserved consistency between diagonal and off-diagonal elements. However, it can also easily be verified that Eq. (15) does not strictly satisfy normalization, since a switching resulting in *a*^{(i)} = *j* within branch (*i*, *j*) yields an excess contribution to the trace of the density matrix. Although this can be remedied by renormalizing the density matrix, for C-FSSH, it was instead suggested to decouple each off-diagonal density matrix elements from *all other* elements, such that for *k* ≠ *l* only the branch pair with *i* = *k* and *j* = *l* contributes,

This leaves diagonal density matrix elements of the form

upon which the conservation of the trace is restored. As a downside of this decoupling, it can no longer be trivially shown that the CS inequality is automatically satisfied, as Eq. (19) can no longer be formulated. What can be shown, however, is that at instances where all branches happen to share the same classical coordinates, *q*^{(i,j)} = ** q**, meaning there is a well-defined, branch-pair-independent adiabatic basis, we have

Although rather restrictive, this example suggests that the CS inequality would have to be satisfied at least approximately in the more general case. For a less ambiguous demonstration, however, numerical calculations are necessary, as presented in Sec. IV. This also applies to Property 3c, for which no analytical relationship can generally be given.

## IV. APPLICATION TO A TRIMER

In Sec. III C, we contributed analytical arguments of why unphysical density matrices arise in FSSH but not in C-FSSH. In the present section, we substantiate these arguments with numerical results for a trimer system. We should point out that no pronounced unphysical behaviors were observed in dimeric systems (not shown here, but for a comprehensive characterization, we refer to Ref. 30). It therefore seems that the combination of multiple pairwise interactions between quantum levels is necessary for such behaviors to become prevalent. In addition to being the simplest system to incorporate multiple pairwise interactions, a trimer is an important model for donor–bridge–acceptor systems widely used for energy and spin transfer processes.^{43–45}

### A. Model

The trimer model considered here is effectively a three-site tight-binding model with open boundary conditions, such that interactions are limited to nearest-neighbors. Within the local basis, the Hamiltonian is given by

where *V* is the interaction strength and *E*_{n} is the diagonal energy associated with local basis state *n*.

The classical coordinates are represented by harmonic oscillators, as detailed in Ref. 30. Bilinear coupling between these coordinates and the quantum system is incorporated through the quantum–classical interaction Hamiltonian

with *K*_{n} ≡ ∑_{α}*g*_{n,α}*q*_{α}, and where *g*_{n,α} is the coupling constant associated with quantum level *n* and classical coordinate *q*_{α}. We assume that each coordinate couples exclusively to a single level *n*, so that the different *K*_{n} are fully uncorrelated. The coupling constants are determined based on a Drude–Lorentz spectral density (also known as the overdamped Brownian oscillator model^{46,47}), taken to be identical for each *n*,

where *ω*_{α} is the frequency associated with classical coordinate *q*_{α}, and where Ω and *λ* are the spectral density characteristic frequency and reorganization energy, respectively. It should be stressed, however, that mixed quantum–classical methods are generally capable of treating arbitrary spectral densities as well as anharmonic effects.

For a given temperature *T*, the initial classical coordinates *p*_{0} and *q*_{0} can be sampled from a Boltzmann distribution

where *β* = 1/*T* and the Boltzmann constant and classical masses are taken to be unity. For FSSH, this distribution can be replaced by the Wigner distribution,

in order to account for zero-point fluctuations that may become particularly significant at low temperatures. For C-FSSH, however, such fluctuations are already represented by the classical wavepackets, such that no Wigner sampling is necessary (or desirable).

### B. Results

For a systematic survey of the trimer system, we refer to the supplementary material, where results from C-FSSH are compared against those from FSSH under Boltzmann sampling [FSSH (B)] and Wigner sampling [FSSH (W)] (through Eqs. (26) and (27), respectively) as well as numerically-exact results from the hierarchical equations of motion^{48} (as implemented in the Parallel Hierarchy Integrator^{49}). Here, a high-dimensional sweep over parameter space is performed including variations in *V*, *λ*, Ω, and *T*. Moreover, three different trimer systems are considered, characterized by the diagonal energies *E*_{n}, the first being a homogeneous trimer with *E*_{1} = *E*_{2} = *E*_{3}, the second being a biased trimer with *E*_{1} − *E*_{2} = *E*_{2} − *E*_{3} = 0.5, and the third being a donor–bridge–acceptor system with *E*_{2} − *E*_{1} = *E*_{1} − *E*_{3} = 1.0, where the reference unit of energy is taken to be 208.5 cm^{−1} (which is the thermal quantum at 300 K). In all cases, the quantum system is initiated at the first level, *n* = 1. In the following, we will be discussing select cases of interest.

The result shown in Fig. 1 was obtained using FSSH (W) for a homogeneous trimer with *V* = 1.0, *λ* = 0.005, Ω = 1.0, and *T* = 0.1, which corresponds to the adiabatic regime at low temperature. It should be pointed out that the low temperature regime, where Ω > *T*, is notoriously demanding for mixed quantum–classical methods, although a previous survey for the spin-boson model has shown C-FSSH to reach surprising accuracy in this regime.^{30} Regardless, it is of interest to assess whether violations of positivity also occur in the high temperature regime where instead Ω < *T* and mixed quantum–classical dynamics is generally expected to behave well. To this end, we show in Fig. 2 the complete set of local populations for a biased trimer with *V* = 1.0, *λ* = 0.005, Ω = 0.1, and *T* = 1.0. Results are shown for FSSH (B), FSSH (W), as well as C-FSSH, alongside numerically exact results.

Interestingly, in Fig. 2, C-FSSH is shown to yield improved accuracy compared to FSSH (B) and FSSH (W), but none of the methods predict negative populations in the local basis. However, as pointed out in Sec. II, this does not necessarily mean that positivity is preserved in *any* basis. To assess this, one should instead consider Property 3c, that is, the eigenvalues of the density matrix, which are depicted in Fig. 3. From this, violations of positivity are observed for both FSSH (B) and FSSH (W), as eigenvalues assume negative values.

Whereas Property 3c is the most robust means to assess positivity violations, a more intuitive assessment is provided by Property 3b, i.e., the CS inequality, as it directly invokes the density matrix elements in some given basis. To this end, we plot in Fig. 4 the quantity *ρ*_{nn}*ρ*_{mm} − *ρ*_{nm}*ρ*_{mn} for each pair of states. Once this quantity becomes negative, we have a violation of the CS inequality and thus a violation of positivity. From Fig. 4, FSSH (B) and FSSH (W) are seen to violate the CS inequality, particularly for *n* = 2 and *m* = 3, although the degree at which this violation occurs is very small.

An even more intuitive analysis of positivity is provided by the CS inequality in the adiabatic basis, as this is the basis in which analytical expressions for the involved density matrix elements are given; cf. Eqs. (14) and (15). As noted in Sec. III A, however, such analysis is only meaningful when, at a given time, the classical coordinates are taken to be the same for every branch and every trajectory, so that a well-defined adiabatic basis exists. We have performed an additional simulation, continuously reinforcing this condition by initiating all classical coordinates identically and by neglecting the quantum contribution to the classical potential energy as well as the adjustment of the classical kinetic energy upon a switch of the active surface. Under these constraints, C-FSSH rigorously satisfies the CS inequality, cf. Eq. (22), which is borne out in Fig. 5, where the CS inequality is assessed for the same parameters as in Fig. 4. Here, FSSH (B) and FSSH (W) are seen to significantly violate the CS inequality as the product of wavefunction coefficients statistically exceeds the average contributions due to active surfaces. Overall, the CS inequalities behave markedly different in the adiabatic basis as compared to the local basis, underscoring the basis dependence of this property. It should be pointed out that the entanglement between the system and environment prohibits the extraction of numerically exact results within the adiabatic basis, as a result of which no numerically exact results are shown in Fig. 5.

Finally, we revisit the low-temperature regime, and present results for a homogeneous trimer with *V* = 1.0, *λ* = 0.025, Ω = 1.0, and *T* = 0.1. Shown in Figs. 6 and 7 are the local populations and eigenvalues of the density matrix, respectively. Significant violations of Property 3c are seen for FSSH (B) and FSSH (W), with eigenvalues being strongly negative. For FSSH (B), this leads to pronounced negative values of *ρ*_{11} and *ρ*_{33}. Interestingly, FSSH (W) leads to an overall improvement in accuracy while enforcing positive values for *ρ*_{11} and *ρ*_{33}, but yields negative values for *ρ*_{22} instead. C-FSSH is once again seen to satisfy positivity throughout, while yielding superior accuracy.

C-FSSH consistently satisfies positivity throughout all of the parameter space covered in the systematic survey presented in the supplementary material, underscoring the robustness of this method in reinforcing the physical properties of the density matrix. FSSH (B) and FSSH (W), on the other hand, are seen to frequently violate positivity in the form of negative density matrix eigenvalues. At high temperatures, this is particularly prevalent for small reorganization energies, whereas at low temperatures violations are observed more consistently. It should be noted that imposing physical properties on the density matrix does not necessarily imply an improvement in accuracy. As reported before,^{30} C-FSSH is oftentimes seen to yield radically improved accuracy compared to FSSH (B and W), but not always, which is underscored by the survey presented in the supplementary material. In particular, FSSH oftentimes performs better in the Marcus regime, where the interaction strength *V* is weak.

## V. CONCLUSIONS AND OUTLOOK

With the growing interest in the application of FSSH to evaluate the quantum dynamics of an entire reduced density matrix, this method is increasingly prompted to satisfy the physical properties that alternative quantum dynamical methods have been subject to. In this paper, we have demonstrated a violation of positivity for density matrices obtained within FSSH following the prevalent implementation^{17,29} where adiabatic coherences are constructed based on wavefunction coefficients rather than the active surfaces that determine populations. We have furthermore shown that C-FSSH,^{30} which invokes a density matrix constructed entirely out of active surfaces, does not suffer from such positivity violations. While a formal proof of positivity within C-FSSH is complicated by the secular approximation, taken to conserve the trace of the density matrix, our numerical survey for a trimeric system revealed this property to be consistently satisfied.

In many cases, the reinforcement of positivity by C-FSSH yields an overall improvement in accuracy over FSSH when compared to numerically exact results. This renders C-FSSH an attractive formalism, especially considering that it comes at virtually the same cost as FSSH. Moreover, as shown in the present paper, it does not require the introduction of additional parameters. However, there are instances where FSSH provides better accuracy, especially when interactions between quantum levels are weak, reorganization energies are large, and quantum levels are degenerate. In such cases, inaccuracies for C-FSSH are predominantly manifested as excessive coherence decay, likely due to an overestimation of hopping rates that cause the adiabatic off-diagonal density matrix elements to become zero [cf. Eq. (20)]. Interestingly, overestimations of hopping rates within this very parameter regime have previously been shown to be resolvable by augmented FSSH (A-FSSH),^{50–52} which was introduced with the aim to address overcoherence known to occur for FSSH. Correcting for overcoherence involves a damping of off-diagonal matrix elements within the adiabatic basis, which should mitigate violations of the CS inequality within this basis and could, therefore, lead to an overall improvement in physical behavior in addition to enhanced accuracy. It would therefore be of interest to combine A-FSSH with the coherent generalization embodied by C-FSSH in order to take advantage of the beneficial properties of both approaches.

It would furthermore be of interest to assess positivity for other methods providing consistent formulations of the entire density matrix based on surface-hopping approaches.^{30–35} Assuring physical and consistent density matrices based on surface hopping is of particular importance for applications to spectral simulations as well as phenomena that lie markedly outside the Marcus regime. As such, the present work may also be of interest to the recently proposed reciprocal-space formulation^{53} of FSSH^{54} for the modeling of bandlike phenomena.

## SUPPLEMENTARY MATERIAL

See the supplementary material for a flowchart of the C-FSSH algorithm and a systematic survey of a trimeric system.

## ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation under Grant No. 2145433.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**A. S. Bondarenko**: Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (lead); Writing – review & editing (equal). **R. Tempelaar**: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

## REFERENCES

We note that an absolute value was taken in the previous work introducing C-FSSH^{30} based on empirical grounds, although it can be argued that this was not physically motivated. Such has been omitted here, bearing in mind that the final form of the overlap factor is real and positive by construction.