Quantum entanglement from classical trajectories

A long-standing challenge in mixed quantum-classical trajectory simulations is the treatment of entanglement between the classical and quantal degrees of freedom. We present a novel approach which describes the emergence of entangled states entirely in terms of independent and deterministic Ehrenfest-like classical trajectories. For a two-level quantum system in a classical environment, this is derived by mapping the quantum system onto a path-integral representation of a spin-1/2. We demonstrate that the method correctly accounts for coherence and decoherence and thus reproduces the splitting of a wavepacket in a nonadiabatic scattering problem. This discovery opens up a new class of simulations as an alternative to stochastic surface-hopping, coupled-trajectory or semiclassical approaches.

Introduction.-Manyimportant phenomena across physics and chemistry are best described by a small quantum system and a large classical environment, for example light-matter interaction, chemical reactions, and qubits.As it is intractable to treat the entire problem with quantum mechanics, it is necessary to simulate the coupled quantum-classical dynamics directly [1].Deriving an approach which is both computationally efficient and accurate is, however, a highly non-trivial task.The simplest method based on classical trajectories is Ehrenfest dynamics, also known as mean-field theory (MFT).While this approach is computationally efficient, it completely neglects quantum-classical entanglement, such as the branching of a nuclear wavepacket in a nonadiabatic scattering problem [2].
Over the last decades, a considerable effort has been invested in the development of more accurate trajectorybased methods.A popular approach, especially in simulations of photochemistry, is Tully's fewest switches surface hopping (FSSH) [3][4][5], whose trajectories take stochastic jumps to simulate wavepacket branching.Although its original form is known to suffer from overcoherence, there have been many suggestions to introduce decoherence corrections [6,7] with little consensus that there is a definitive solution.Another known way to include entanglement is to use coupled trajectories, either on top of Ehrenfest [8] or surface hopping [9] or through methodologies such as the exact-factorization framework [10][11][12], ab-initio multiple spawning [13], the quantumclassical Liouville equation [14], or Bohmian dynamics [15].A third possibility is to use interference between path histories and weight Ehrenfest-like trajectories (obtained from a mapping scheme [16][17][18] which has a close relation to the Stratonovich-Weyl representation used in the present paper [19,20]) by phases and prefactors derived from a semiclassical propagator based on a realtime path integral [21][22][23][24].At first sight, decoherence and entanglement appear to be inherently quantum phenomena which cannot be described with a fully classical simulation [25].However, in this Letter, we introduce a new approach that, in contrast to the three approaches described above, can capture these effects based on independent and deterministic classical trajectories.
Our theory is based on the Stratonovich-Weyl (SW) phase-space representation of the quantum system, which is a Wigner representation of discrete spaces [26,27].For simplicity, we consider only the two-level case, which employs the well-known isomorphism to a spin S = 1 2 system, representing the spin by a classical vector of length S(S + 1).We propose to extend this approach to a path integral of spin vectors, where the centroid of the spin path determines the dynamics and the initial configuration specifies the weight of each trajectory.This weight, which is preserved along the trajectory, contains the information necessary for quantum-classical entanglement.
Method.-First, consider an isolated two-level quantum system with density matrix ρ.A convenient classical analogue for this system is given by the Stratonovich-Weyl W-representation [28], which expresses the expectation value of an operator Â as an integral, The classical functions are defined as ρ(s) = tr[ρ ŵ(s)] (and likewise for A(s)) where ŵ(s) = 1 2 Î + s • σ is the SW kernel, Î is the 2 × 2 identity matrix, σ = [σ x , σy , σz ] are the Pauli matrices, and s is a vector with magnitude For the integration measure we use the shorthand notation d 2 s = 1 2π dϕ dθ sin θ, where ϕ and θ are the spherical coordinates of s.Since each Cartesian component, s j , is the SW representation of the spin operator Ŝj = 1 2 σj , one can think of s as a classical spin vector with the familiar quantum magnitude S(S + 1) of a spin S = 1 2 (where = 1 throughout).Next, consider the time evolution of the density matrix.As is well known, the dynamics of a two-level system is equivalent to that of a spin- 1  2 in an effective magnetic field H, where the Hamiltonian is Ĥ = H 0 Î + H • Ŝ.Using this decomposition, it is straightforward to write where ρ 0 = 1 2 is fixed by the normalization.When the Liouville-von-Neumann equation, d dt ρ = i[ρ, Ĥ], is converted to its phase-space equivalent, it follows that the standard precession formula for the classical spin vector, ṡ = s × H, generates the correct quantum dynamics.When coupled to a general classical environment (described by coordinates x, mass m and conjugate momenta p), the total Hamiltonian, The corresponding equations of motion are [19] in addition to the spin dynamics as before.While these equations of motion are equivalent to those of Ehrenfest dynamics [29], the SW treatment differs in the initial distribution: while standard Ehrenfest starts from a unique vector s of length 1 2 (as in the Bloch-sphere picture), the SW approach averages over all initial spin directions in Eq. ( 1) and uses the magnitude √ 3 2 .We have recently found that the latter, called the linearized spin-mapping method, leads to a better prediction of population dynamics [19,20,30].Other mapping approaches have also found an effective spin magnitude of √ 3 2 to be optimal [31,32], and averaging over initial directions to be beneficial [33], even with the Ehrenfest spin length [34].However, one important drawback is present in both Ehrenfest and linearized spin mapping, namely that the dynamical quantization is lost.This has the unfortunate consequence that after a scattering event, the trajectories evolve on a weighted average of the two product potential energy surfaces, in contrast with the correct entangled state which splits into parts on one or the other surface [2].We will now show that such quantization can be systematically reintroduced by representing the system by a path integral of spins.In contrast to standard spin coherent-state path integrals, we do not require paths to be continuous in the N → ∞ limit and therefore do not have to deal with the difficulties that arise when restricting to such paths [35][36][37].
By construction, the SW representation has an inversion formula ρ = d 2 s ρ(s) ŵ(s), (5) with the particular example of the identity, Î = d 2 s ŵ(s).By applying Eq. ( 5) to both operators in tr[ρ Â] and inserting resolutions of the identity, we can generalize Eq. (1) to a path integral of N spins, (6) where we symmetrized over the indices l and m and used Eq. ( 1) for terms with l = m.Due to the linearity of the SW representation, it follows that 1 N l ρ(s l ) = ρ(s) (and similar for A), where we introduced the centroid s = 1 N l s l .The expression looks like a classical phase-space average with a weight function g Note that if we had used |s| = 1 2 , the weight function would reduce to that of standard spin coherent-state path integrals [38].However, we find that our choice |s| = converges quicker in N .
A practical consideration is that g N ({s k }) is a complicated complex-valued function that varies rapidly for high N .However, since the observables depend only on s and not on the relative geometry, it is possible to rigorously integrate all degrees of freedom other than the centroid.Explicitly, we define Note that the weight function g N ({s k }) is invariant under global rotations of the spin vectors [39], so that G N is spherically symmetric, G N (s) = G N (s), where s = |s|.Equation ( 6) thus simplifies to Since the centroid of N > 1 points on a sphere can reach any point inside the sphere, the integration domain of s is the ball |s| ≤ , which recovers Eq. ( 1).The resulting universal function G N (s) has several important properties: (1) it depends only on the centroid magnitude not on its direction, (2) it is real-valued, (3) it is independent of the Hamiltonian and the initial conditions.In other words, even though its computation becomes exponentially hard with increasing N , it only has to be computed once for a given N , hence the name universal.We have evaluated G N (s) numerically up to N = 16 using Monte Carlo [40]. Figure 1(a-b) shows that the universal function consists of a few positive and negative domains, but the number of nodes seems to remain small for high N .The simulation will thus include trajectories with both positive and negative weights but this does not lead to a severe sign problem [40].
Next, we consider the distribution of the spin components, Ŝj .Quantum-mechanically these are expected to be quantized with the eigenvalues ± 1  2 , but the integrand in Eq. ( 1) is smeared over all spin directions.However, as N increases, the centroid distribution of Eq. ( 8) becomes peaked around sj = ± 1 2 for all j ∈ {x, y, z}, with heights that are consistent with the components of ρ, as shown in Figure 1(c-d).In other words, the path-integral weight function G N (s) reintroduces the quantization to the system that is necessary for quantum-classical entanglement.
Finally, consider time-dependent expectation values.Using a similar argument as in Eq. ( 2), one can describe the dynamics in the spin path-integral representation by a homogeneous precession of all spins, ṡk = s k × H. Consequently, the centroid evolves in the same way, s = s × H and we do not need to keep track of the individual spin vectors.Since G N (s) is invariant under global rotations, its value is preserved by the dynamics, which has the important implication that Eq. ( 8) is valid for all times.
For an isolated system, this gives exact time-dependent expectation values for any value of N .For the coupled quantum-classical problem, we propose the approximation (9) where the phase-space version of the density operator ρ(x, p, s) involves a Wigner transform of the environment in addition to the SW transform of the quantum system (and likewise for A).This equation is the main result in this Letter and will be referred to as the spin pathintegral method.It is exact at t = 0 and in the limit of an isolated system for all N .The N = 1 case uses the same dynamics and spin distribution as the linearized spinmapping method [19] and in this more general formula, the accuracy is expected to increase with N due to the quantization of the spin vectors.
Results.-We have applied the spin path-integral method to Tully's seminal scattering problems [3], which are well-known benchmark models but also proxies for realistic chemical reactions [41].The results are compared against calculations using numerically exact quantum mechanics as well as Ehrenfest dynamics and surface hopping.Simulation details are included in the Supplemental Material (SM) [40].
First, we consider the single avoided crossing (model I) with diabatic surfaces and coupling shown in the inset of Fig. 2.An initial Gaussian wavepacket enters from the left on the lower surface, ψ 1 (x, t = 0) ⊗ |1 , with enough kinetic energy that both product channels are open.Due to nonadiabatic coupling, the wavepacket splits into two separate wavepackets on the two surfaces and emerges as an entangled state, , on the right.By quantum-classical entanglement, we mean that if the nuclei are in a certain region of phase space, then we know with certainty the quantum state of the electrons and vice versa.In Fig. 2 we show the distribution of the final momentum for two different initial energies.As is well known, Ehrenfest dynamics (MFT) is unable to capture the branching of the wavepacket, whereas surface hopping (FSSH) provides a reasonably accurate description for this model.The N = 1 simulation predicts an envelope that covers the full range of momenta allowed by energy conservation, but lacks the two-peak structure.However, by increasing N , we find that the distribution smoothly splits into two parts and thus recovers the quantum-mechanical entanglement.
We emphasize that the dynamics consists of independent and deterministic trajectories on a weighted average of the two states, similar to both Ehrenfest dynamics and the linearized spin-mapping method, and the key difference lies in the weighting of the trajectories.Because the weights may be positive or negative, some of these cancel out in such a way that the ensemble branches when it emerges on uncoupled surfaces.This cancellation is reminiscent of more involved semiclassical methods such as Miller's forward-backward propagator, which is also known to capture wavepacket splitting in the present model [21][22][23].However, these approaches are inherently semiclassical, not classical, and include nuclear-coherence effects (to some level of approximation) via phases and prefactors that depend sensitively on the trajectory his- x FIG. 2. Probability distribution of the nuclear momentum after an avoided crossing (model I).The inset shows the diabatic potentials (solid lines) and coupling (dashed line).A wavepacket enters from the left on the lower surface with a narrow distribution of kinetic energies at roughly 1.5 (left panels) or 5 (right panels) times the asymptotic energy difference.Ehrenfest (MFT) gives a single peak around the average momentum, while linearized spin mapping (N = 1) envelopes the exact wavepacket distribution.For higher N , the spin pathintegral method correctly reproduces the wavepacket branching.
tories and make sampling difficult.The results of the simpler spin path-integral method demonstrate that only electronic coherence is necessary to recover the correct result.Although the trajectories also carry a sign, this depends in a relatively simple manner on a single degree of freedom, is fixed by the initial sampling, and is preserved by the dynamics.For Tully's dual avoided crossing (model II) we reach the same conclusions, and in the SM we show that the scattering probabilities are in good agreement with exact wavepacket calculations for a wide range of initial momenta [40].
Next, consider the more challenging extended coupling model (model III) shown in the inset of Fig. 3.As before, an initial wavepacket enters on the lower surface from the left but now the total energy is low enough for the upper channel to be closed on the right.During the collision, it thus splits into a transmitted part on the lower surface and a part on the upper surface which reflects and passes through the interaction region a second time.Surface hopping is well known to fail dramatically for systems with recrossing, because the electronic amplitudes picked up during the first crossing are inconsistent with the active surfaces [4].This 'overcoherence' problem arises because the assumption of a unique trajectory for each electronic density matrix is not valid [42] and is related to neglecting quantum-classical entanglement.Ehrenfest and various linearized mapping The inset shows the adiabatic surfaces (solid) and nonadiabatic coupling (dashed).A wavepacket enters from the left and is partly reflected.Only the spin path-integral method with N = 4 is able to correctly describe the second crossing of the interaction region (results for N = 8, not shown, overlay with N = 4).
approaches have also been unable to describe this model [43].
To quantify the overcoherence problems of quantumclassical simulations, we have calculated the time evolution of the impurity which is a measure of entanglement and is related to the decoherence indicator studied in Ref. 11 with coupledtrajectory simulations.Here, ρ nm denotes elements of the reduced density matrix in the adiabatic representation and the results are shown in Fig. 3. Ehrenfest completely misses the second crossing at about 100 fs (since its trajectories do not reflect), and although some surface-hopping trajectories do reflect, FSSH is unable to correctly describe the entanglement in this system.The spin path-integral method on the other hand reproduces the correct result for this system.
Another well-known consequence of the overcoherence problems in surface hopping are erroneous oscillations [44] in the scattering probabilities as shown in Fig. 4. For the spin path-integral method, we observe that the calculated scattering probabilities converge towards the correct values with increasing N (although reproducing the step as the upper channel opens appears to be difficult).Note that we did not need to add 'decoherence corrections' for each trajectory (as is commonly done to fix surface hopping), but nevertheless do not observe the problems of overcoherence for the ensemble as a whole.
Finally, we note that unlike surface hopping, the results of the present method (like Ehrenfest and other mapping approaches [21,45]) are not dependent on whether the adiabatic or diabatic representation is used.
Conclusions.-In this Letter we have showed that features of quantum-classical entanglement, such as wavepacket branching and impurity measurements, can indeed be captured by an ensemble of independent and deterministic classical trajectories.This discovery opens While surface hopping (FSSH) suffers from erroneous oscillations, the spin path-integral method appears to be converging smoothly with increasing N towards the exact scattering probabilities.
up for a new class of mixed quantum-classical methods, as an alternative to surface hopping, coupled-trajectory or semiclassical simulations.It also extends the applicability of mapping approaches, which have been successful for predicting electronic coherences but so far have struggled to describe the nuclear dynamics of scattering problems.The presented method relies on positive and negative trajectory weights whose sign cancellation does not become more difficult for larger systems or longer simulation time.We therefore expect it to be applicable to complex molecular systems and condensed-phase problems.
Here we have limited the treatment to two-level systems, but a multi-level extension already exists for linearized spin mapping [20] and the spin path-integral extension is straightforward (although there is no guarantee that G N will depend on only a scalar variable).Since the SW formalism can be applied to any symmetry group [46], a similar treatment could be made also in systems with different symmetries.
Particularly interesting is the case where ρ is a thermal density matrix.Since the weights are preserved by the dynamics, we expect this to be useful for equilibrium dynamics as the quantum Boltzmann distribution will automatically be conserved.The details are left to a forthcoming paper.
Quantum entanglement from classical trajectories: Supplemental Material In the comma-separated-values file that is part of this Supplementary Material, we provide values of 4πs 2 G N (s) on a grid of s for N = 2, 4, 6, 8, 12, and 16.We obtained these results with very long Monte Carlo runs to sample g N ({s k }) and histogrammed the centroid length s.To use the results in a spin path-integral simulation, we recommend interpolating it with a spline fit.If one samples s uniformly from [0, √ 2 ] and the direction s/s uniformly from the unit sphere, then the correct weight is given by the interpolated function.

Model systems
Model I was defined by the diabatic potential with A = 0.01, B = 1.6, C = 0.005, D = 1.This is a slightly modified version of Tully's first model problem [1] in the form used by Miller and coworkers [2].
Model II is Tully's dual avoided crossing model defined by the diabatic potential with A = 0.1, B = 0.28, C = 0.015, D = 0.06, E 0 = 0.05.Model III was defined by the diabatic potential with A = 6 × 10 −4 , B = 0.1, C = 0.9.This is identical to the third of Tully's model problems [1] although the state labels have been swapped compared to the original paper so that 1 is always the index of the lowest state.
In each case, the mass was m = 2000 and all quantities are reported in atomic units.

Initial conditions
The initial condition in all models was a product state ρ = ρn ⊗ ρe , where ρn = |ψ ψ| and ρe = |1 1|.For the simulations shown in Figs.2-3 of the main text, the nuclear wavefunction was a Gaussian wavepacket of the form S4) where x 0 = −15 is well to the left of the interaction region.The numerically-exact wavepacket simulations were thus initiated with |ψ ⊗ |1 .
For MFT and FSSH, the initial values of x and p were sampled from Wigner distribution of the nuclear part, together with a method-specific treatment of the electronic part ρe = |1 1|.For MFT, this corresponds to all trajectories starting with a spin at the north pole of a sphere of radius 1 2 .Note that the spin path-integral dynamics are initialized from all parts of the spin sphere, in contrast with Ehrenfest and FSSH.Therefore, in order for the total energy distribution to be consistent with the wavepacket calculation, the initial conditions need to take account of the sampled value of V (s).It is necessary to do this because the same set of spin path-integral trajectories can in principle be used to study the alternative scattering problem with the same total energy distribution initialized on the upper state.The simplest algorithm for obtaining the appropriate energy distribution is to sample a momentum p from the distribution Eq.(S5) to obtain E = (p ) 2 2m + V 11 .Then the corrected momentum p is defined by the solution of E = p 2 2m + V (s), and this is used to initialize the trajectory dynamics.Although this modification appears to change the initial momentum distribution, it can be seen that in the N → ∞ limit the resulting distribution is formally correct (see Figure S1) due to the fact that V (s) → V 11 on account of the peaking of the distribution of s.Simple extensions of this algorithm can be constructed to treat multidimensional systems by sampling a distribution of E and the direction of the momentum vector.
In all spin path-integral dynamics calculations, we used the single-variable universal function G N (s) rather than the more complicated g N ({s k }).The results show that the linearized spin-mapping method (N = 1) is already very accurate for this problem except in the low-energy regime and that the remaining error is significantly decreased when using the path-integral spin-mapping for N > 1.The conclusions from this model are thus in line with those discussed in the main text.

Analysis of the impurity in Model III
Figure S3 shows separately the quantities |ρ 12 |(t) = (Re ρ 12 (t)) 2 + (Im ρ 12 (t)) 2 and ρ 11 (t) • ρ 22 (t) that are used in computing the impurity in the main text.(All quantities refer to the reduced density matrix in the adiabatic representation.)The off-diagonal element is a direct measure of coherence, and the second crossing gives rise to a recoherence that is correctly described by spin path-integral method, but not by Ehrenfest or surface hopping.These decoherence indicators are related to that used in Ref. 12 to test the quantum-classical coupled-trajectory solutions of the exact factorization formalism.

Convergence details
The number of trajectories used in Model I was 10 6 in order to provide fully-converged high-resolution histograms, but the overall trend can be observed already with 10 3 trajectories (i.e.just tens of trajectories per histogram bin), see Fig. S4.Converged results for model II were computed from 10 5 trajectories.All spin-dynamics calculations for model III used 10 6 trajectories, although the important trends can already be clearly seen with just 10 4 , see Fig. S5.
The convergence rate is mostly determined by the average sign of f (s) [Eq.(S6)], which is reported in Table I.It is clear that it is far more efficient to use the  single-variable universal function G N (s) instead of the path-integral weight g N ({s k }), especially as N increases.These results are universal for any two-level quantumclassical system given a pure initial state in the quantum system, and as stated in the main text, independent of both the trajectory length and the number of environmental degrees of freedom.The values imply that typically a simulation with N = 16 will require only one order of magnitude more trajectories than N = 1 for the same level of convergence.Fig. S6 shows the contributions of positively and negatively weighted trajectories separately, and it is clear that the two-peak structure emerges from cancellation of trajectories with intermediate momenta.

1 2
FIG. 1.(a) Weight of the centroid spin magnitude s (arbitrary scaling).The negative regions grow in importance for increasing number of spins, N , but the number of nodes remains small.(b) Weight of s with positive contributions in blue and negative in red.(c-d) Distribution of the x, y (c) and z (d) components of the spin centroid.The distributions become peaked around the quantum-mechanical values ± 1 2 with increasing N , and the relative peak heights approach the corresponding expectation values of the density matrix, here plotted for ρ = 2 3 |1 1| + 1 3 |2 2|.

FIG. 3 .
FIG. 3. Impurity in the extended coupling model (model III).The inset shows the adiabatic surfaces (solid) and nonadiabatic coupling (dashed).A wavepacket enters from the left and is partly reflected.Only the spin path-integral method with N = 4 is able to correctly describe the second crossing of the interaction region (results for N = 8, not shown, overlay with N = 4).

FIG. 4 .
FIG. 4. Transmission and reflection probabilities onto the adiabatic states in the extended coupling model (model III).While surface hopping (FSSH) suffers from erroneous oscillations, the spin path-integral method appears to be converging smoothly with increasing N towards the exact scattering probabilities.
Johan E. Runeson * and Jeremy O. Richardson † Laboratory of Physical Chemistry, ETH Zürich, 8093 Zürich, Switzerland (Dated: November 5, 2021) FIG. S2.Probability of transmission on the upper state in the dual avoided crossing model (model II).Inset shows adiabatic surfaces (solid lines) and nonadiabatic coupling (dashed line, divided by 12).Our surface hopping (FSSH) results are in agreement with the equivalent calculations of Ref.[11].
FIG. S4.Final momentum distribution in Model I for an initial wavepacket centred around an energy of p 2 0 /2m = 0.03.Higher N generally requires more trajectories for convergence although the two-peak structure can already be observed with only 1000 trajectories.

2 0
FIG. S6.Final momentum distribution in Model I for an initial energy p 2 0 2m = 0.03.Cancellation of positively weighted (dashed line) and negatively weighted (dotted line) samples is what yields the two-peak structure (solid line).

TABLE I .
FIG. S5.Scattering probabilities inModel III (extended coupling model) with 10 4 trajectories for each initial momentum.Top row: average sign when sampling f (s) in Eq. (S6) for various N .Bottom row: the corresponding average sign if one instead integrated over N spin variables weighted by gN ({s k }).In our simulations we use GN (s), corresponding to the top row.These values are modelindependent and valid for any pure initial state.