Multi-level description of the vibronic dynamics of open quantum systems

Recent experimental observations of oscillatory features in the time-resolved two-dimensional electronic spectra of natural photosynthetic light-harvesting complexes have stimulated significant theoretical investigation of the origin of such features. Subsequent realization that these oscillations could arise from excitation of vibrational degrees of freedom (DOFs) that are coupled to electronic transitions probed in these experiments suggests that these coupled modes could play a role in framing the ultrafast resonant excited state vibronic excitation energy transfer pathways. In fact, recent experimental and theoretical investigations of model light-harvesting complexes, such as phycobiliproteins found in cryptophyte algae, have highlighted the existence of strongly coupled intramolecular vibrational modes with frequencies matching those observed in spectroscopic signals.^1–5 Motivated by these observations, this manuscript details the development of a novel coherent state density matrix (CSDM) propagation method for computing the vibronic dynamics based on the semiclassical Herman-Kluk propagator,^6,7 which is complementary to the previously developed partially linearized density matrix (PLDM) propagation scheme that is applicable to high-dimensional dissipative open quantum system dynamics.^8–10 

In previous studies,^11,12 a discrete vibronic system basis was employed to retain vibrational quantum effects. However, the unfavorable scaling with system basis set size restricts this treatment to very few vibrational modes and potentially unrealistic truncation of the number of basis states. The implementation outlined in this manuscript avoids this need for very large numbers of system basis states while retaining the full vibrational state space for all modes. It does this by using a highly efficient moving basis set of coherent states for describing the quantum vibrational states.

Modeling the excited state dynamics of molecular aggregates, such as photosynthetic light-harvesting complexes, requires accurate parameterization of the vibronic Hamiltonian. One commonly used Hamiltonian model for this purpose is the Frenkel-exciton Hamiltonian, a tight-binding model that is expressed in a localized diabatic basis. In the single excitation manifold and within the bilinear harmonic bath approximation, this Hamiltonian is expressed as

where $({Q^i},{P^i})$ are the phase-space operators of the bath oscillators, M_k are the bath oscillator masses, and $hα,β({Q^i})$ are the bath position-dependent electronic Hamiltonian matrix elements, which, for the model studied here, are assumed to take to form

where $ϵ̃α=ϵα+∑k=1Noscck(α)2/2MkΩk2$ is the excitation energy of molecule α at the Franck-Condon point, ϵ_α is the minimum energy of the potential in excited state α, $ck(α)$ is the local gradient on the excited state potential evaluated at the Franck-Condon point that defines the coupling between the excited electronic state of molecule α and the kth bath oscillator, and Δ_α,β is the electronic coupling between excited states of molecules α and β which, for simplicity, is assumed to be independent of bath configuration. Here, the state |α⟩ represents a single exciton manifold electronic basis function in which molecule α is in its first excited state, i.e., $ϕ1(α)$⁠, while all other molecules are in their electronic ground state, $ϕ0(n)$⁠, i.e., $|α=|ϕ0(1),…,ϕ1(α),ϕ0(α+1),…,ϕ0(Nst)$⁠.

In general, the spectral density for a specific electronic state that defines the coupling-weighted frequency profile of the environment, $Jα(ω)=π∑k=1Noscωkλk(α)δ(ω−ωk)=π2∑k=1Noscck(α)2ωkδ(ω−ωk)$⁠, can be decomposed into short-range (intramolecular) and long-range (intermolecular) components.^5,13–15 The long-range intermolecular component is often well-represented by a broad low-frequency continuum of weakly coupled modes, while the short-range intramolecular component generally consists of a discrete set of more strongly coupled modes.

Efficiently and accurately computing the ultrafast nonadiabatic dynamics of a system described by the above Hamiltonian represents an important challenge in the theory of open quantum systems. While there has been significant effort in the development of algorithms that treat problems like this exactly, e.g., the hierarchical equations of motion (HEOM) approach, which is exact for restricted spectral density forms,^16–18 or the more general grid based multiconfiguration time dependent Hartree (MCTDH) approach,¹⁹ or path integral based methods,²⁰ that scale exponentially with the number of environmental modes, there are also a host of approximate, perturbative methods (e.g., Förster and Redfield based theories) that can be reliable in various regimes of system-environment coupling, or strength of coupling between quantum subsystem states,^21–23 or that make different dynamical approximations (e.g., secular or Markovian) about relative relaxation time scales, or duration of memory effects in these complex systems.²⁴ Another class of popular approximate approaches that attempt to overcome the exponential scaling of quantum dynamics with the number of DOFs is the mixed quantum-classical techniques based on surface-hopping or Ehrenfest^25,26 ideas that attempt to partition the DOFs into different sets that are described with different quantum and classical levels of treatment and devise some approximate, often heuristic, means of allowing exchange of energy between these subsystems. A more rigorous approach for treating these types of systems is offered by the mixed-quantum classical Liouville formulation.^27–31 

Instead of treating different DOFs with different approximate dynamical methods, there are more controlled approaches that manipulate the representation of the Hamiltonian and attempt to treat all DOFs consistently at some semiclassical level, for example, the linearized semiclassical initial value representation (LSC-IVR)^13,32,33 or the truncated Wigner approximation (TWA).³⁴ The approach we outline below is a hybrid of these methods that builds on semiclassical ideas, which become exact for harmonic environmental modes; it employs the mapping Hamiltonian description, providing a convenient way to implement the electronically nonadiabatic dynamics semiclassically, a coherent state description to treat high frequency vibrational modes, and a linearized approximation to derive classical-like dynamics for low frequency environmental modes giving a highly flexible approach that can be tuned to optimize accuracy of description and algorithmic efficiency. Often, treatment of high-dimensional open quantum system models necessitates restrictive approximations to the dynamical description of the vibrational bath that, when applied to more strongly coupled higher-frequency modes, can result in inaccurate representation of time evolution. Below, we will demonstrate in simple open quantum system models where accurate benchmarks studies can be performed that the highly flexible approach developed here provides a useful bridge to modeling multichromophore exciton transport systems that include the important influence of many higher frequency chromophore vibrations as well as low frequency environmental modes.

The Meyer-Miller-Stock-Thoss (MMST) mapping model provides an exact representation of discrete basis states in terms of harmonic oscillator phase-space operators. This is particularly convenient for semiclassical approximations to the time evolution of a vibronic system since it provides an equivalent dynamical footing for the required classical propagation of both electronic and vibrational DOFs. Within the mapping Hamiltonian formalism, a given basis state of the system, |α⟩, is mapped onto the excitation of a harmonic oscillator DOF associated with this state. We write the general transfer operator in Eq. (1) using the following mapping $|αβ|→âα†|00|âβ=|01,…,1α,0α+1,…,0Nst01,…,1β,0β+1,…,0Nst|$⁠, where $âα†=12ℏ(q^α−ip^α)$ (with unit mass and frequency for the mapping oscillator variables) is the harmonic oscillator raising operator representing quantum subsystem state |α⟩ and $|0=|01,…,0Nst$ is the ground state of the set of harmonic mapping oscillators. The equivalent mapping Hamiltonian in Eq. (1) can be expressed in these mapping oscillator phase space operators as^35,36

The physical space in which each mapping oscillator moves must be restricted to only the ground, 0_α, and first excited, 1_α, states of these oscillators to reliably describe the changing dynamical occupation of the electronic states these mapping oscillators represent. It can be shown³¹ that when the mixed quantum classical Liouville equation, for example, is cast in the mapping formulation, the mapping projector and the mixed quantum-classical Liouville operator commute so formally the dynamics governed by this equation must stay in the physical subspace. If approximations are made to implement this equation, however, they may give rise to an evolution in which the mapping variables leave the physical subspace and this approximate implementation may result in an unphysical description of the nonadiabatic dynamics. Here, however, we implement the mapping variable dynamics semiclassically within the coherent state representation, which provides an exact prescription for quadratic variables like the mapping oscillator DOFs appearing in the mapping Hamiltonian. The environmental variables will be treated with varying levels of approximation; however, with our implementation, the mapping variable dynamics can never leave the physical subspace by construction. To see this, we note that using Hamilton’s equations derived from Eq. (3), it is easily shown that $ξ=∑α(qα2+pα2)$ is a constant of motion for the semiclassical mapping Hamiltonian dynamics, independent of the dynamics of the bath subsystem variables. In previous work^8,12,38 developing our partial linearized path integral methods, we showed that the nonadiabatic evolution of the full density matrix can be approximated by propagating semiclassical trajectories of the mapping variables forward and backward in time representing the dynamics of the amplitudes of the bra and ket states moving under the influence of the mean of the forward and backward paths of the environmental DOF to linear order in the difference environmental path. This partial linearization in the environmental variables leads to expressions analogous to Eq. (16), discussed below, in which the contribution from each sampled trajectory is weighted by initial ground state mapping oscillator distributions [by sampling the $G(q→,p→)$ function defined below Eq. (16)] and explicitly weighting by the product of coherent state transformed first Hermite polynomials appearing in the factor $Ts(nf,n0)(qnft,pnft,qn00,pn00)$ defined under Eq. (12), thus incorporating the projection onto the relevant mapping oscillator excited states, e.g., $1nft$⁠. This sampling and weighting approach guarantees that our partial linearized implementation of the mapping variable dynamics never drifts from the physical subspace. Thus, a feature common to all the partial linearized density matrix propagation approaches that we have developed based on the mapping Hamiltonian formulation is that they must reliably remain in the physical projection space. The only error we make is in our treatment of the nuclear dynamics, which in PLDM involves including only terms to linear order in the difference paths of all nuclear DOFs, while in the CSDM approach, we go beyond this and employ a coherent state description for the forward and backward propagators for the high frequency harmonic modes and, as we outline in our multilevel implementation, we can incorporate low frequency anharmonic modes through a PLDM description, which can be combined with the new CSDM propagation approach to give a highly accurate hybrid propagation scheme.

The Herman-Kluk semiclassical propagator provides an exact quantum dynamical treatment of systems with quadratic potentials, like the Hamiltonian in Eq. (3). As such, it can be used as an exact propagator for system or bath DOFs within the MMST mapping model. However, as soon as the coupled dynamics of the system and bath begins, the excited vibronic wavefunction can no longer be factorized into electronic and vibrational components, and the Herman-Kluk propagator is approximate for the composite wavefunction. The Herman-Kluk propagator is expressed as^6,7

where S_t is the action evaluated along a classical trajectory connecting the coherent state phase space boundary values in a time, t, and |q, p⟩ is a coherent state defined in the coordinate representation as ⟨x∣q, p⟩ = exp[−γ(q − x)² + ip(x − q)/ℏ]. The factor, C_t, provides information about the stability of a classical trajectory’s endpoint in phase-space with respect to a change in its initial phase-space point and is defined as

where we have chosen the coherent state width matrix, γ, to be diagonal. Computation of C_t can be troublesome due to its potentially highly oscillatory complex nature. However, for the composite vibronic system and for the specific choice of γ_i = m_iω_i/2ℏ (for a general harmonic oscillator, either environment or mapping, with mass m_i and frequency ω_i), the stability factor takes the form (omitting oscillator labels)

Here, for example, the diagonal block corresponding to the bath DOFs is defined as

where $acl(Q,P)=MΩ2ℏ(Q+iMΩP)$ is the classical analogue of the harmonic oscillator lowering operator. The time evolution of this quantity for uncoupled classical linearly driven harmonic oscillators [with $V(Q)=12MΩ2Q2+Γ(t)Q$] can be expressed as

where Γ(t) is proportional to the mapping variable portion of the system-bath coupling term in Eq. (3). Imposing the boundary condition at t = 0, we find

In the high-frequency bath limit, the endpoint in phase space of the semiclassical trajectory used to evolve the coherent state basis functions for these discrete bath oscillator DOFs will be weakly dependent on the initial mapping variable phase space point so the stability factor can be approximated as³⁷ 

Therefore, the quantities of interest are found to be $∂acl(Qt,Pt)∂(Re(acl†(Q0,P0)))=exp[−iΩt]$ and $∂acl(Qt,Pt)∂(Im(acl†(Q0,P0)))=−i⁡exp[−iΩt]$⁠. Using these results in Eq. (7), together with the fact that det[exp[A]] = exp[Tr A], the stability factor is seen to be $det[12CtQQ]12=exp[−i2∑k=1NoscΩkt]$⁠. A similar result has previously been derived for the phase space mapping operators of the electronic subsystem, where $det[12Ctqq]12=exp[−i2∑α∫0tdτhα,α({Q^i})]$⁠.^37,38 Thus, while approximate, we avoid the necessity of propagating the stability factor for the computing the time-evolution of the vibronic density matrix.

It can be shown that, with our choice of γ_i = m_iω_i/2ℏ, the projection of an arbitrary harmonic oscillator eigenstate into the coherent state representation can be computed as

where $⟨q,p|0⟩γ=mω2ℏ=πℏmω14⁡exp[−mω4ℏq2−14ℏmωp2+i2ℏpq]$ is the coherent state representation of the ground harmonic oscillator eigenstate. This provides a useful recursive definition for higher excited harmonic oscillator eigenstates in this representation. For the remainder of this manuscript, we choose the coherent state width parameter to be $γi=miωi2ℏ$ for all harmonic oscillators, i, without loss of generality.

In a number or previous publications,^8,12,38 the Herman-Kluk propagator has been used to evolve the mapping variable dynamics and compute transition amplitudes. Combining this with the resulting transition amplitude for the bath oscillators, the time evolution of the combined system-bath density matrix can be expressed as (with ℏ = 1, as will be used throughout the remainder of this manuscript)

where $S̃t[q→,p→,Q→,P→]=∫0tdτp→⋅q̇→+P→⋅Q̇→−Hmap−12∑αhα,α({Qi})$ is the modified action including the mapping variable stability factor and the system and bath transition amplitudes are defined as $Ts(nf,n0)(qnft,pnft,qn00,pn00)=12(qnft+ipnft)(qn00−ipn00)$ and $Tb(ν→f,ν→0)(Q→,P→)=∏k=1Nosc(1ν0(k)!acl†(Qk0,Pk0))ν0(k)1νf(k)!acl(Qkt,Pkt)νf(k)$⁠, respectively. Here, the subscripts 0 and t label initial and final time-evolved values of the phase space variables.

Combining the complex phases from the ground harmonic state projections and the modified action in the first and second lines of Eq. (12), we find that the total phase for the forward propagator, e.g., is

where we have used the fact that $p→t⋅q→t−p→0⋅q→0=∫0tdτddτ(p→⋅q→)=∫0tdτ(ṗ→⋅q→+p→⋅q̇→)$ and inserted Hamilton’s equations for the time derivatives of phase space points according to the Hamiltonian in Eq. (3) for mapping and harmonic oscillator phase space variables. Thus, the complex phases associated with trajectory endpoints and the action are replaced solely with the system-bath coupling portion of the Hamiltonian. Moreover, the amplitude factors controlled by the Gaussian functions that arise from final projections of the ground state harmonic oscillator onto the coherent state basis can be rewritten utilizing Hamilton’s equations, resulting in

The results in Eqs. (13) and (14) can then be combined to define a complex weight functional as

which, for the bilinear coupling model, is linear in bath phase space variables. Combining this with the fact that, for the mapping variables, $ddt∑α=1Nst(qα2+pα2)=0$⁠, the time evolving density matrix can be expressed as

where $G(q→,p→)=exp[−12∑α=1Nst(qα2+pα2)]$ and $F(Q→,P→)=exp[−(∑i=1NoscMiΩi2Qi2+12MiΩiPi2)]$ are Gaussian phase space distributions for the system mapping variables and bath DOFs, respectively. The proportionality is used because of the absence of an overall normalization factor. The algorithm then consists of sampling the initial phase space variables from the appropriate initial distributions, running classical trajectories by integrating Hamilton’s equations for the Hamiltonian in Eq. (3), and computing the desired final projections by evaluating the transition amplitude functions and weighting each the trajectory by the exponential factor in Θ_t.

In what follows, we will consider initial conditions where the initial bath density is provided by the Boltzmann distribution for an isolated system of harmonic oscillators as $ρ^0=ρ^s⊗e−βĤbZ$⁠, where $Ĥb=∑k=1NoscP^k22Mk+12MkΩk2Q^k2$⁠, $β=1kBT$ is the inverse temperature, and Z is the partition function. Such an initial condition is appropriate for modeling photoexcitation of the system. Because of the convenient recursive definition of the coherent state representation of harmonic oscillator eigenstates given in Eq. (11), we can concisely represent the Boltzmann operator in terms of coherent states as

Often, we are only interested in observing properties of the system, while averaging over the states of the bath. Following a similar procedure used above for the Boltzmann distribution, the trace over a harmonic bath DOF can be computed according to

Combining these results and utilizing Hamilton’s equations for the time-evolution of the coherent state phase space variables, the reduced density matrix, σ(t), initialized from the Boltzmann distribution of the isolated bath, can be expressed as

where $Φt[Q→,P→,Q→′,P→′,q→,p→,q→′,p→′]=Θt[Q→,P→,p→,q→]+Θt*[Q→,P→,p→′,q→′]−Θt*[Q→′,P→′,p→′,q→′]−Θt[Q→′,P→′,p→,q→]$is the complex functional, the exponential of which weights each trajectory’s contribution to the reduced density matrix. The correlated initial bath distribution function is given by