Communication: Partial linearized density matrix dynamics for dissipative, non-adiabatic quantum evolution

An approach for treating dissipative, non-adiabatic quantum dynamics in general model systems at finite temperature based on linearizing the density matrix evolution in the forward-backward path difference for the environment degrees of freedom is presented. We demonstrate that the approach can capture both short time coherent quantum dynamics and long time thermal equilibration in an application to excitation energy transfer in a model photosynthetic light harvesting complex. Results are also presented for some nonadiabatic scattering models which indicate that, even though the method is based on a "mean trajectory" like scheme, it can accurately capture electronic population branching through multiple avoided crossing regions and that the approach offers a robust and reliable way to treat quantum dynamical phenomena in a wide range of condensed phase applications.

The mixed quantum-classical strategy for modeling systems in which quantum coherent dynamics and electronically non-adiabatic transitions play important roles describes the nuclear degrees of freedom (DOF) classically or semiclassically, while treating the electronic DOF quantum mechanically with an evolution operator parameterized by trajectories of the nuclear DOF. 1 Though many schemes for implementing this idea have been developed, 2 fundamental questions arise about the accuracy of dynamics methods that treat electronic and nuclear DOF on different dynamical footing. 3 To overcome this difficulty, the mapping Hamiltonian idea that exactly maps discrete quantum states onto continuous coordinates was proposed 4 and enables consistent treatment for all DOF. This idea replaces the evolution of the electronic subsystem by the dynamics of a system of fictitious mapping harmonic oscillators. With this approach, e.g., the quantum amplitude transfer operator transforms as |β λ| →â † βâ λ , wherê a λ = 1 √ 2¯(q λ − ip λ ), and a general electronic Hamiltonian in the diabatic representation,ĥ el = β,λ |β β|ĥ el (R)|λ λ|, can be rewritten asĥ map (R) = 1 2¯ β h β,β (R)(q 2 β +p 2 β −¯) + 1 2¯ λ =β h βλ (R)(q βqλ +p βpλ ), where (P ,R) and (p,q) represent the nuclear and mapping oscillator phase space DOF, respectively. Direct implementation of this mapping Hamiltonian by stationary phase approximation and classical trajectories encounters a fatal problem when q 2 β + p 2 β <¯, since some of the classical DOF can evolve on an inverted potential surface, proportional to −h β, β (R). 5 Moreover, the population, ρ ββ =â † , is not guaranteed to have a positive expected value when the approach is implemented approximately, e.g., with the linearized semiclassical initial value representation (LSC-IVR) applied to multi-state systems. 5, 6 LSC-IVR linearizes in the difference a) Author to whom correspondence should be addressed. Electronic mail: coker@bu.edu. between forward and backward paths for both the mapping and nuclear DOF. To overcome these problems we use the coherent state representation of the mapping DOF, 7 and linearize only in the nuclear DOF; allowing different mapping variable paths for the forward and backward propagators, i.e., partial linearization. The quantity of interest is the evolution of the density matrix involving forward and backward propagation: Here, the total Hamiltonian isĤ =P 2 /2M + h map (R,p,q), and n t label the basis states at time t. The propagator matrix elements in discrete path integral form are where is a time step, the nuclear kinetic ac- ], and T [n t ,n 0 ] = n t |e − ī ĥ map (R N−1 ) ....e − ī ĥ map (R 0 ) |n 0 is the nuclear path dependent quantum transition amplitude.
In the coherent state representation (with coherent state width parameter, γ = 1/2 and setting units so that¯=1), the transition amplitude can be expressed as 7 T [n t ,n 0 ] = dq 0 dp 0 1 4 (q n t + ip n t )(q n 0 − ip n 0 )c t e iS 1 (t) The term2 β h ββ (R) in the action, S 1 , that gives rise to the problem of inverting the potential can be eliminated as this term is cancelled exactly by the pre-factor, c t , 7 leaving S cl Here we use the idea of partial linearization 9 in the nuclear DOF that involves transforming the forward and backward nuclear path variables, R and R , to mean and difference variables:R = (R + R )/2 and Z = (R − R ), respectively (with similar definitions for the nuclear momenta,P and Y). The nuclear kinetic action difference becomes The central approximation with this approach involves truncating the phase difference in the combined transition amplitude terms to linear order in Z, based on the assumption that for short times, forward and backward nuclear paths will remain close to each other. This may appear to be a restrictive approximation that will only be valid for very short times for high dimensional problems, but such linearization approximations have been shown to be reliable even when forward and backward paths differ significantly in some degrees of freedom. 10 With the linearization approximation, the key term is written as Using this to expand the action difference in Z, and the fact that With this expression the first term in Eq. (6) cancels the boundary terms in T [n t ,n 0 ] given in Eq. (3) with similarity for the backward path transition amplitude, T [n 0 ,n t ] .
Combining forward and backward phase factors ) and performing the integrals over Z 0 ....Z N − 1 gives our approximation for ρ(t): Here, β 0 ) are the initial distributions for the forward and backward mapping variables that satisfyq n t = ∂h cl map (R t )/∂p n t andṗ n t = −∂h cl map (R t )/∂q n t , and the nuclear trajectories are determined by a "mean field" like force resulting from the different forward and backward mapping paths: The mean nuclear DOF initial distribution is the partial Wigner transform: (ρ) We use factorized initial conditions, ρ 0 = ρ eq bath (R)ρ sys , though the non-separable case can be treated. 11 Numerical implementation of Eq. (7) involves sampling initial nuclear DOF from (ρ) n 0 ,n 0 W (R 0 ,P 1 ) and mapping variables from the Gaussian functions. The product of δ-functions in Eq. (7) gives a time-stepping prescription for evolving the mean nuclear DOF with the force in Eq. (8). We refer to the approach presented above as partial linearized density matrix (PLDM) propagation.
The PLDM propagation scheme is a "mean trajectory" approach, but according to Eq. (8) it is different from "Ehrenfest" or LSC-IVR dynamics where the classical DOF force depends only on one set of mapping variables as 8,12 The PLDM force in Eq. (8) is governed by both forward and backward mapping DOF and the relationship between these different semiclassical schemes will be detailed in a forthcoming publication. 13 By following a different sequence of canceling and linearizing in the phase difference, a "surface hopping" version of the present theory was derived known as iterative linearized density matrix (ILDM) propagation. It has been tested in various quantum dynamical problems. 14, 15 While ILDM uses a single surface, linearized short time approximation that is iterated to treat longer times, its convergence with large numbers of iterations can be problematic. The LANDmap approach 7 uses this same linearized approximation but without iteration. This improves statistical convergence (for examples below typically 10 8 trajectories can be reduced to 10 4 ) but the linearized propagator underlying these methods is generally only reliable for short times and a balance between many iterations and statistical convergence must be considered. The mean trajectory linearized approximation underlying the PLDM scheme developed here, however, is generally accurate for much longer times as demonstrated below. The PLDM approach thus offers a significant improvement in statistical convergence (requiring only 10 4 -10 5 trajectories) while preserving high accuracy, even at long times.
The example applications presented below demonstrate that the PLDM propagation mean trajectory approach provides a generally accurate representation of electronic branching in model scattering situations though such methods usually do not give reliable nuclear information. 8 We also present results which show that in a dissipative, non-Markovian open quantum system model the PLDM propagation method reliably captures both the coherent short time evolution and the approach to thermal equilibrium at long times. We compare with benchmark results obtained with the hierarchical coupled reduced master equation (HCRME) (Refs. 16 and 17) approach that can be converged to give exact results for restricted forms of the model system-bath coupling and spectral density. The PLDM approach can, in principle, be applied to arbitrary model system bath interactions. Thus, the method offers a very favorable tradeoff between accuracy, computation efficiency, and statistical convergence for applications to general high dimensional, multi-state, quantum dissipative systems. On the other hand the HCRME approach assumes bilinear system-bath coupling, a harmonic bath, and has very large demands on memory. 18 In Fig. 1 we present results for two nonadiabatic multistate scattering models that have been used in benchmark studies of other approaches. Results from PLDM propagation are compared with those of other approximate schemes and exact results. These studies demonstrate that PLDM propagation provides reliable population branching results for nonadiabatic dynamical problems over a range of models and parameters.
Next the PLDM propagation is applied to a dissipative non-Markovian model of excitation energy transfer in the Fenna-Matthews-Olsen (FMO) light harvesting complex. This model system Hamiltonian is detailed in Ref. 16. Figure 2 compares the evolution of initially occupied site 1 population computed using various approximate methods. , HCRME original, 16 re-scaled, 17 LANDmap, 15 LSC-IVR, 6 PBME, 20 and ILDM propagation. 14  PLDM propagation gives results in excellent agreement with results from the HCRME and ILDM approaches 16,17 that, in principle, are numerically exact.
The new approach accurately captures both the short time coherent dynamics and, as demonstrated in Fig. 3, the method is capable of generating stable converged results at very long times, as the system and bath exchange energy and approach thermal equilibrium. The Boltzmann equilibrium populations of the site basis states, ∼ e −β α , are also displayed in the figure. We see that by 10 ps the state populations are accurately approaching their equilibrium values. The other approximate implementations of mapping, 6, 15, 20 on the other hand, fail to give the correct thermal equilibrium limit at long times (e.g. LSC-IVR (Ref. 6) even gives negative populations). Similar equilibration problems have been reported with the stochastic Liouville equation approach, 21 the original "Ehrenfest dynamics" scheme 22 and the Haken-Strobl-Reineker approximation 23 for such complex quantum dissipative systems. A detailed analysis of precisely why the PLDM approach is capable of describing thermal equilibration will be presented in a forthcoming full publication on the approach. 13 The results displayed in the different panels of Fig. 3 explore the effects of starting in different initial electronic states of the quantum network and demonstrate that the PLDM approach reproduces the correct equilibration from different initial states.
This new approach is generally applicable for arbitrary system-bath interactions (beyond the bi-linear form assumed here) and can be readily applied to treat non-Markovian spectral densities. 15,21 The new approach could also be interfaced with a general quantum mechanics/molecular mechanics (QM/MM) description of a complex dissipative quantum system. The PLDM propagation approach can be efficiently applied to treat both the short time transient coherent quantum dynamics that is addressable with ultrafast non-linear spec-troscopy and reliably capture the long time behavior associated with approach to thermal equilibrium.