Communication: Fully coherent quantum state hopping

In this paper, we describe a new and fully coherent stochastic surface hopping method for simulating mixed quantum-classical systems. We illustrate the approach on the simple but unforgiving problem of quantum evolution of a two-state quantum system in the limit of unperturbed pure state dynamics and for dissipative evolution in the presence of both stationary and nonstationary random environments. We formulate our approach in the Liouville representation and describe the density matrix elements by ensembles of trajectories. Population dynamics are represented by stochastic surface hops for trajectories representing diagonal density matrix elements. These are combined with an unconventional coherent stochastic hopping algorithm for trajectories representing o ff -diagonal quantum coherences. The latter generalizes the binary (0,1) “probability” of a trajectory to be associated with a given state to allow integers that can be negative or greater than unity in magnitude. Unlike existing surface hopping methods, the dynamics of the ensembles are fully entangled, correctly capturing the coherent and nonlocal structure of quantum mechanics. C LLC. initial

Probably, the most popular method for combining classical dynamics with quantum transitions is trajectory surface hopping, [1][2][3][4][5][6]21 and, in particular, variants of Tully's fewest switches surface hopping (FSSH). 2 Together, these methods are well-known for their mistreatment of coherence. The origin of the problem is that each trajectory carries with it its own complete quantum state, which is fully coherent. This is fundamentally incorrect, as quantum effects in a trajectory representation emerge as the inability to separate the unified whole of a quantum state into independent ensemble members. Quantum mechanics entangles classical trajectory evolution. 22,23 From a trajectory ensemble perspective, quantum coherence results from the relationships between ensemble members and is not a property of any single trajectory. This perspective has been emphasized in our previous work. [12][13][14][15][16]22,23 Recently, a great deal of attention has been brought to correct the overcoherence problem while retaining the general surface hopping approach (see, e.g., Ref. 21).
In this paper, we present a fully coherent stochastic quantum state hopping method for simulating quantum dynamics. We consider the simple but unforgiving two-state system evolving coherently in a pure state without external perturbations and then add random noise to the ensemble evolution to demonstrate the method's ability to accurately capture decoherence effects in both stationary and nonstationary environments. a) Email: cmartens@uci.edu We formulate our approach in the Liouville representation and describe the quantum evolution by a density matrix. 24 The Hamiltonian of the system is given bŷ where the diagonal elements satisfy E 2 − E 1 = ω. The offdiagonal coupling is taken to be a real constant V . The state of the system is described by the density matrix where ρ 11 (t) and ρ 22 (t) are the time-dependent populations of states |1⟩ and |2⟩, respectively. The coupling V induces transfer of population between the states. The off-diagonal elements ρ 12 (t) and ρ 21 (t) = ρ * 12 (t) are the coherences. The density operatorρ(t) undergoes pure state evolution under the quantum Liouville equation, This can be broken down into coupled ordinary differential equations for the density matrix elements, We represent the time evolution of the density matrix in terms of ensembles of stochastic realizations, which we refer to as trajectories. First, consider the evolution of populations, which are represented collectively by an ensemble of N trajectories. Each trajectory carries a discrete random variable σ j (t), j = 1, 2, . . . , N which can take on values of either 0 or 1. If σ j (t) = 1 at time t, then the jth trajectory occupies state 1, while if σ j (t) = 0, the trajectory occupies state 2. The density matrix elements can then be written as sums over the trajectory ensemble, where N 1 (t) and N 2 (t) are the number of trajectories associated with states 1 and 2, respectively, and N 1 (t) + N 2 (t) = N. This conservation of trajectories is honored by the equations of motion, Eqs. (4) and (5), asρ 11 +ρ 22 = 0. The dynamics of the σ j (t) are represented by a stochastic hopping algorithm. In terms of N 1 (t) = N ρ 11 (t) we have, for a time step ∆t, This gives the number of ensemble members that will join or leave the state 1 portion of the ensemble during the time step ∆t.
The direction of change depends on the sign of the right side of the equation. In practice, ∆t is chosen so that the magnitude of ∆N 1 during the time step is much less than N. The probability that a given member of the ensemble will hop is determined as follows. If the jth trajectory is already a member of the state 1 part of the ensemble with σ j (t) = 1 and ifρ 11 is positive, nothing happens to the status of the trajectory. Ifρ 11 is negative, the probability that the trajectory will undergo a transition to state 2 is given by The probability P(2 → 1) for hops from state 2 to state 1 is computed in an analogous manner. The numerical implementation is achieved by generating a uniform random number ξ between 0 and 1 for each candidate population ensemble member at each time t and executing the hop if ξ ≤ P. The coherence ρ 12 (t) is also represented by an ensemble of N trajectories. This aspect of the method is novel and differs from previous surface hopping approaches. As we will see, it correctly captures the nature of quantum coherent dynamics.
It is convenient to write the coherence in amplitude-phase form: ρ 12 (t) = R(t)e iΦ(t) . In terms of the trajectory ensemble, we have This is a complex number that results from summing N complex contributions from each trajectory. The nature of the contributing trajectory amplitudes r j and phases φ j will be addressed below. We consider the increment of ρ 12 during a time interval ∆t: ∆ρ 12 =ρ 12 ∆t. In terms of increments of the amplitudes and phases of the ensemble, The real and imaginary parts give the relations From the equations of motion, Eqs. (4)-(6), we have Equating Eqs. (13) to (15) and (14) to (16) gives the conditions for determining the increments ∆r j and ∆φ j . This is a highly under-determined set of equations. A great deal of freedom is thus present for creating an algorithm for consistently updating their values. This is reminiscent of the gaugelike freedom of quantum trajectory methods in general that we discussed previously. 25 To establish an algorithm, we impose constraints, but note that a different choice of constraints will generate a different dynamics for the trajectory ensemble that nonetheless will be an equivalent representation of the quantum state evolution.
With that constraint imposed, conditions for ∆r j can be derived by equating Eqs. (13) with (15) and (14) with (16). After some algebra, these can be reduced to A solution for ∆r j is then given by where We now implement Eq. (19) with a stochastic algorithm. We interpret ∆r j as a generalized probability of change during the time interval ∆t. In particular, we generate a random number ξ between 0 and 1 and compare the magnitude of ∆r j with this number. If |∆r j | ≤ ξ, the trajectory undergoes a coherence transition. This is accomplished by changing the value of r j by an integer amount. If ∆r j > 0, then r j → r j + 1. For a negative ∆r j , we let r j → r j − 1. Unlike the population case, the "probability" that a trajectory is in a coherent state is not necessarily a positive number nor is its magnitude bounded by unity. Rather, the coherence state of a given trajectory is an integer: r j = . . . ,  1, 0, 1, 2, 3, . . .. This is reminiscent of the "generalized probability" status of the Wigner function, 26,27 where its magnitude can be larger than allowed by a classical probability density, while its value can be negative. Restricting the sign or magnitude of r j in the numerical implementation leads to poor results in practice.
We test the method by treating the coherent evolution of the two-state system with Hamiltonian given by Eq. (1). The dynamics of this system consists of Rabi oscillations induced by the coupling. We consider the initial density matrix with ρ 11 (0) = 1, ρ 22 (0) = ρ 12 (0) = ρ 21 (0) = 0. In Fig. 1(a), we show results for the degenerate case ω = 0, V = 1, while the dynamics for ω = 5, V = 1 case are given in Fig. 1(b). The population ρ 11 (t) is shown in blue, while the imaginary part of the coherence Im ρ 12 is shown in red. These are compared with the exact solution (solid curves). An ensemble of N = 10 000 is used in all cases, which gives effectively quantitative agreement with exact results. (We note that ensembles as small as N = 100 give reasonable and stable, if noisy, results.) In Fig. 2, we show results for the two-state system perturbed by a stationary random environment. This is modeled by adding to ω a fluctuating term δω(t), characterized by Gaussian white noise with a standard deviation η. This corresponds in the modified Liouville equation governing the quantum system as damping of the coherences, adding a term −γ ρ 12 (t) to the equation for ρ 12 (t) and a corresponding term for ρ 21 (t), where γ = η 2 /2. In Fig. 2(a), the simulated evolution of ρ 11 (t) and Im ρ 12 (t) is compared with exact results for the case ω = 0, V = 1, and η = 0.5, while Fig. 2(b) gives results for the more strongly damped η = 1.0 case. Both systems show decaying coherent oscillations which relax to an incoherent mixture with equal populations in each state.
For both the undamped and damped systems shown above, the coherent surface hopping method gives nearly quantitative agreement with exact results. These are very simple systems that nonetheless provide rigorous tests for our stochastic method for modeling populations and coherence, as coherent quantum effects dominate the dynamics.
We now apply our approach to the nontrivial and largely unexplored problem of quantum dynamics coupled to a nonequilibrium environment. We consider a nonstationary stochastic model introduced by us previously to describe the decay of initial coherence coupled to a local nonequilibrium bath, for which analytic results could be obtained. 28,29 The model is based on a nonequilibrium local phonon oscillator, which modulates the difference frequency ω by a term where Ω is the local phonon frequency and t j corresponds to relative time origin for which the oscillator is excited by a quantum transition. The oscillator phase θ(t) takes on a definite initial value θ(0) = θ o , which is related to the mechanical properties of the system and bath and which leads to classical mechanical coherence of the ultrafast bath response. As time progresses, the phase θ(t) undergoes diffusion from the initial value θ o with diffusion constant D. The interplay of quantum coherence and the classical mechanical coherence of the bath oscillator can influence the quantum dynamics, and the behavior of the system can depend significantly on the control parameter θ o . 28,29 A number of possible realizations of this nonequilibrium model in the context of our stochastic method suggest FIG. 2. Evolution of ρ 11 (t) and Re ρ 12 (t), as in Fig. 1. (a) η = 0.5 and (b) η = 1.0 (right). For both systems, ω = 0, V = 1, and N = 10 000. themselves. In the present paper, our focus is on the quantum simulation method itself, and so, we adopt one simple prescription and leave a detailed analysis of alternatives for future work.
In our approach, each trajectory (r j (t), φ j (t)) of the coherence ensemble interacts with its own environmental ensemble member. We include the effect of fluctuations of the instantaneous difference frequency δω(t) resulting from both stationary Gaussian noise with strength η, as above, and the nonstationary local phonon.
The first hop that creates coherence induces a nonequilibrium bath response and so adds to the fluctuating difference frequency a term δω(t) = δω o cos[Ω(t − t j ) + θ(t − t j )], where t j is the time of the transition of r j . Subsequent transitions of r j (t) to other nonzero values impulsively add to the response of the evolving local phonon. In general, this will generate a series of trigonometric terms, each with its own relative time zero t (k) j and stochastic phase history θ (k) (t − t (k) j ). We represent this approximately as a single effective oscillatory term, which simplifies the method while retaining a nonstationary bath with control parameter θ o .
Immediately after a transition of r j (t), the term δω(t) is updated to become Fig. 3(a), we show the dynamics of the state 1 population ρ 11 , while the evolution of the linear entropy S L = Tr (1 −ρ 2 ) = 2(ρ 11 ρ 22 − | ρ 12 | 2 ) is given in Fig. 3(b). We investigate the dependence of the dynamics on the initial phase θ o for the case V = 0.5, ω = 1, Ω = 1, δω o = 0.05, D = 0.1, and η = 0.5. The trajectory ensemble size is again N = 10 000. Three initial phase values are shown: θ o = 0 (blue), θ o = 1.0 (red), and θ o = π (green), along with the results for which the oscillator is absent, with only the stationary Gaussian component present (black dashed). The pure dephasing of the coupled system ultimately degrades the Rabi oscillations and leads to equilibration of the populations to ρ 11 = 0.5. The linear entropy S L measures the purity of the state: a pure state is characterized by S L = 0, while a completely incoherent and equilibrated system decays to S L = 0.5.
A pronounced dependence of the population dynamics on the initial oscillator phase θ o is revealed in Fig. 3(a). In particular, the θ o = 0 case exhibits a significantly slower decay to equilibrium than the purely damped system. The θ o = 1.0 shows a small positive effect, while θ o = π leads to an increase in the population transfer at intermediate times, depressing ρ 11 below the pure dephasing case.
In Fig. 3(b), we show the corresponding time dependence of the linear entropy S L . The damped system in the absence of the local nonstationary oscillator exhibits a time-dependent growth of S L from its initial value of S L = 0, indicating a pure state, to the equilibrium value of S L = 0.5. The nonequilibrium system with θ o = 0 shows a strong suppression of the rate of growth of the linear entropy, especially at intermediate to long times. The effect is less pronounced for θ o = 1.0. In contrast, the system with θ o = π exhibits an enhanced rate of increase of S L .
Persistent nonequilibrium bath coherence is created by and subsequently fully entangled with the quantum evolution. The interplay between quantum system and classical bath has a strong influence on coupled dynamics. As discussed in Refs. 28 and 29, the initial phase θ o of local environmental oscillator acts as a control parameter for this bath-mediated "coherent control" process. We are unaware of other numerical methods that can simulate fully entangled nonequilibrium system-bath dynamics, which will be explored in more detail in future publications.
In the present method, population and coherence hops are determined stochastically from a global ensemble level representation. The evolution of a given trajectory is "entangled" with that of all the other trajectories. This is consistent with the underlying exact quantum equations of motion, where the functions representing the state of the system-the individual density matrix elements-are coupled directly to each other. This entanglement and ensemble level nature of the evolution does not single out a given trajectory's history as observable. Each member has a definite history, but different stochastic realizations would give different histories without changing anything measurable or observable. Quantum coherence in the present approach is represented at the ensemble level through interrelationship between contributions of each trajectory that determines the coherence of the state. Energy conservation is also not obeyed by individual trajectories, but only holds at the ensemble level. Random perturbations that affect each trajectory in a distinct manner lead to decoherence by eroding systematic relations between different ensemble members. These perturbations, which single out members of the ensemble by the distinctions of each stochastic realization, are essentially acts of measurement, and so modifyand ultimately destroy-the fragile coherences that characterize quantum dynamics. In contrast with the conventional surface hopping methodology, this fundamental characteristic of quantum mechanics is correctly incorporated in the present approach.