A new approximate coherent state path integral approach, which enables accurate and efficient dynamical treatment of model Hamiltonians that incorporate excited electronic states of multiple chromophores that are coupled to discrete high frequency harmonic vibrational modes, is presented. The approach is based on the mapping Hamiltonian formalism for the electronic states together with semiclassical coherent state expressions for the forward and backward propagators describing the quantum bath modes. The density matrix dynamics is propagated in the full coherent state basis for the electronic mapping and discrete vibrational mode oscillators using ensembles of weighted trajectories. An effective scheme for projecting the ensemble onto selected vibronic basis states is presented enabling the evolution of the reduced system density matrix to be monitored as well as exploring the importance of selected vibronic relaxation pathways in the multichromophore system dynamics. The approach is demonstrated for simple model Hamiltonians, and we show how this coherent state density matrix propagation approach for high frequency discrete harmonic vibrational modes can be combined with partial linearized density matrix propagation to treat an additional continuum bath of low frequency environmental modes that could, in principle, include anharmonicity.

## I. INTRODUCTION

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.

## II. THEORY

### A. Exciton models and the mapping Hamiltonian formulation

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\alpha ,\beta ({Q^i})$ are the bath position-dependent electronic Hamiltonian matrix elements, which, for the model studied here, are assumed to take to form

where $\u03f5\u0303\alpha =\u03f5\alpha +\u2211k=1Noscck(\alpha )2/2Mk\Omega k2$ is the excitation energy of molecule *α* at the Franck-Condon point, *ϵ*_{α} is the minimum energy of the potential in excited state *α*, $ck(\alpha )$ 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 *k*th 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., $\varphi 1(\alpha )$, while all other molecules are in their electronic ground state, $\varphi 0(n)$, i.e., $|\alpha \u2009=|\varphi 0(1),\u2026,\varphi 1(\alpha ),\varphi 0(\alpha +1),\u2026,\varphi 0(Nst)\u2009$.

In general, the spectral density for a specific electronic state that defines the coupling-weighted frequency profile of the environment, $J\alpha (\omega )=\pi \u2211k=1Nosc\omega k\lambda k(\alpha )\delta (\omega \u2212\omega k)=\pi 2\u2211k=1Noscck(\alpha )2\omega k\delta (\omega \u2212\omega 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,^{19} or path integral based methods,^{20} 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.^{24} 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).^{34} 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 $|\alpha \u2009\u2009\beta |\u2192\xe2\alpha \u2020|0\u2009\u20090|\xe2\beta =|01,\u2026,1\alpha ,0\alpha +1,\u2026,0Nst\u2009\u200901,\u2026,1\beta ,0\beta +1,\u2026,0Nst|$, where $\xe2\alpha \u2020=12\u210f(q^\alpha \u2212ip^\alpha )$ (with unit mass and frequency for the mapping oscillator variables) is the harmonic oscillator raising operator representing quantum subsystem state |*α*⟩ and $|0\u2009=|01,\u2026,0Nst\u2009$ 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^{31} 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 $\xi =\u2211\alpha (q\alpha 2+p\alpha 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\u2192,p\u2192)$ 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.

### B. Semiclassical description of vibronic dynamics

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*)^{2} + *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\Omega 2\u210f(Q+iM\Omega 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\Omega 2Q2+\Gamma (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^{37}

Therefore, the quantities of interest are found to be $\u2202acl(Qt,Pt)\u2202(Re(acl\u2020(Q0,P0)))=exp[\u2212i\Omega t]$ and $\u2202acl(Qt,Pt)\u2202(Im(acl\u2020(Q0,P0)))=\u2212i\u2061exp[\u2212i\Omega 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[\u2212i2\u2211k=1Nosc\Omega kt]$. A similar result has previously been derived for the phase space mapping operators of the electronic subsystem, where $det[12Ctqq]12=exp[\u2212i2\u2211\alpha \u222b0td\tau h\alpha ,\alpha ({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 $\u27e8q,p|0\u27e9\gamma =m\omega 2\u210f=\pi \u210fm\omega 14\u2061exp[\u2212m\omega 4\u210fq2\u221214\u210fm\omega p2+i2\u210fpq]$ 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 $\gamma i=mi\omega i2\u210f$ 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\u0303t[q\u2192,p\u2192,Q\u2192,P\u2192]=\u222b0td\tau \u2009p\u2192\u22c5q\u0307\u2192+P\u2192\u22c5Q\u0307\u2192\u2212Hmap\u221212\u2211\alpha h\alpha ,\alpha ({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\u2009)(qn00\u2212ipn00\u2009)$ and $Tb(\nu \u2192f,\nu \u21920)(Q\u2192,P\u2192)=\u220fk=1Nosc(1\nu 0(k)!acl\u2020(Qk0,Pk0))\nu 0(k)1\nu f(k)!acl(Qkt,Pkt)\nu f(k)\u2009\u2009$, 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\u2192t\u22c5q\u2192t\u2212p\u21920\u22c5q\u21920=\u222b0td\tau \u2009dd\tau (p\u2192\u22c5q\u2192)=\u222b0td\tau \u2009(p\u0307\u2192\u22c5q\u2192+p\u2192\u22c5q\u0307\u2192)$ 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

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

where $G(q\u2192,p\u2192)\u2009\u2009\u2009\u2009=\u2009\u2009\u2009\u2009exp[\u221212\u2211\alpha =1Nst(q\alpha 2\u2009\u2009\u2009\u2009+\u2009\u2009\u2009\u2009p\alpha 2)]$ and $F(Q\u2192,P\u2192)=exp[\u2212(\u2211i=1NoscMi\Omega i2Qi2+12Mi\Omega 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 $\rho ^0=\rho ^s\u2297e\u2212\beta \u0124bZ$, where $\u0124b=\u2211k=1NoscP^k22Mk+12Mk\Omega k2Q^k2$, $\beta =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 $\Phi t[Q\u2192,P\u2192,Q\u2192\u2032,P\u2192\u2032,q\u2192,p\u2192,q\u2192\u2032,p\u2192\u2032]=\Theta t[Q\u2192,P\u2192,p\u2192,q\u2192]+\Theta t*[Q\u2192,P\u2192,p\u2192\u2032,q\u2192\u2032]\u2212\Theta t*[Q\u2192\u2032,P\u2192\u2032,p\u2192\u2032,q\u2192\u2032]\u2212\Theta t[Q\u2192\u2032,P\u2192\u2032,p\u2192,q\u2192]$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

The first two lines above can be sampled directly via Cholesky transformation, and we include the complex phase factor as a weight.

## III. RESULTS

### A. Simple vibronic model demonstration

The results presented in Fig. 1 explore the implementation of the CSDM propagation scheme, outlined above, for a simple vibronic model Hamiltonian consisting of two excited electronic states represented as spin up, ↑, and spin down, ↓, quantum states coupled in different ways to a discrete harmonic vibrational DOF, *Q*. Thus, the equilibrium position of the spin down state is shifted by *c*/(*M*Ω^{2}) relative to the spin up state, which has its equilibrium position at *Q* = 0. A schematic of the potentials and the lowest few diabatic vibronic states, together with details of the model parameters, are presented in Fig. 1(a).

Figure 1(b) compares population dynamics computed using CSDM propagation with results obtained using the PLDM approach in which the dynamics of the discrete vibrational DOF is treated assuming the density matrix propagation can be linearized in the difference between the forward and backward paths of the vibrational DOFs. Also shown in the figure are the results of the truncated Wigner approximation (TWA), in which all DOFs are linearized in difference paths. For the model parameters studied here, the coupling between the electronic states and the coupling between the electronic subsystem and the vibrational DOF are comparable, and it is found that photosynthetic light harvesting systems often function in this so-called intermediate coupling regime where there is no clear small parameter to justify applying perturbation theory. Under these conditions, we see that the classical bath approaches (PLDM and TWA) are reliable only at very short times. Also presented in the figure are results obtained by propagating the mapping Hamiltonian dynamics in the limited vibronic basis including just the 5 lowest-energy vibronic states [see states indicated schematically in Fig. 1(a)]. While the mapping variabledynamics corresponding to the restricted vibronic basis is exactly captured by the Herman-Kluk propagator due to the quadratic dependence of the mapping Hamiltonian, the tensor product Hilbert space grows exponentially with the number of vibrational modes included in this description. Thus, the exponentially scaling size of the basis needed to describe the vibronic wavefunction quickly becomes problematic with an increasing number of vibrational modes. The CSDM approach, however, avoids this issue while retaining linear scaling with respect to the number of vibrational modes at the cost of making a semiclassical approximation to the vibronic propagator. The classical bath approaches also achieve linear scaling but offer a less accurate description of high-frequency vibrational modes than the CSDM approach, as shown in Fig. 1(b).

Figure 1(c) explores a parameter regime where the classical bath approaches are valid, and the approximation made to the Herman-Kluk prefactor breaks down. In such parameter regimes, a classical description is valid, suggesting a partitioning scheme for a multilevel description of vibronic dynamics.

### B. Multilevel description of vibronic dynamics

Since the system mapping variable transition amplitude used in the current analysis is the same as that used in the PLDM algorithm, the combination of these two methods for a multilevel description of vibronic dynamics follows straightforwardly. We consider the time-evolution of the composite vibronic system, consisting of a discrete set of high frequency harmonic modes as well as a continuum of potentially anharmonic modes, as governed by

where $V0({R^j})$ is the bare bath potential and $h\alpha ,\beta ({Q^i},{R^j})=\u03f5\u0303\alpha \u2212\u2211k=1NQck(\alpha )Q^k+V\alpha ,\beta ({R^j})\delta \alpha ,\beta +\Delta ({R^j})\alpha ,\beta (1\u2212\delta \alpha ,\beta )$ describes the system-bath coupling, with no restriction to the dependence of the potential on the bath DOFs, ${R^j}$. Following an identical procedure described in Ref. 8, the approximate time-evolution of the reduced density matrix for the multilevel vibronic system can be expressed as

where $\rho Wn0n0\u2032(R\xaf0,P\xafR1)$ is the Wigner transform of the initial continuum bath density, *k* labels a discrete time index, $Fk=\u221212\u2207R\xafk(H\u0303({Qi},{R\xafj},q\u2192,p\u2192)+H\u0303({Q\u2032i},{R\xafj},q\u2032\u2192,p\u2032\u2192))$ is the instantaneous bath force with $H\u0303=Hmap+12\u2211\alpha =1Nsth\alpha ,\alpha ({Qi},{Rj})$, where all operators in *H*_{map} in Eq. (21) are replaced by their corresponding classical variables. We leave the vector symbol off of the $(R\xaf,P\xaf)$ DOFs for notational simplicity.

Figure 2 explores the behavior of these different open system quantum dynamics methods for a more challenging, realistic model involving the same two level electronic subsystem, but now, bilinearly coupled strongly to *N*_{d} discrete harmonic vibrational modes, and *N*_{c} = 100 more weakly coupled continuum environmental modes sampled from a Drude-Lorentz Ohmic spectral density. The parameters employed in the model are summarized in the caption. Figure 2(a) shows that for *N*_{d} = 1 weak to moderately coupled discrete mode, with frequency comparable to the difference in electronic state energies, the PLDM dynamics spin state populations agree reasonably well with the results obtained from the CSDM propagation approach. When a second discrete vibrational mode with a higher frequency is included [Fig. 2(b)], we see that for the strongest electronic-vibrational coupling, there are significant deviations of the PLDM results from the CSDM approach. We have also presented results obtained without including the linearized continuum of low-frequency modes, treating the two discrete modes with the CSDM approach (solid line) or in an expanded vibronic basis (dashed line). To achieve convergence, the vibronic basis, in this case, was required, including the 10 lowest energy vibronic diabatic states. For larger numbers of discrete vibrational DOFs, this representation becomes prohibitively expensive, compared to the CSDM scaling, as noted above. Since the Herman-Kluk propagator is not norm-conserving for the potential considered here, all results have been renormalized by the sum of populations.

## IV. CONCLUSIONS

While the CSDM propagation method outlined here is restricted to the description of harmonic oscillators, which often provide a reliable approximation to the strongly coupled, higher frequency intramolecular vibrational DOFs that this method was developed to treat, the approximate PLDM propagation scheme is not restricted to such potentials. While PLDM struggles to treat discrete high frequency vibrational modes, it has been previously shown to accurately describe thermalization effects from coupling to a continuum of lower frequency, possibly strongly anharmonic vibrational modes.^{8}

In contrast to standard implementations of the Herman-Kluk propagator, our approximation to the stability prefactor allows us to avoid propagating monodromy matrix elements. Because of this, the time taken to compute a single trajectory scales linearly with the number of harmonic DOFs represented by coherent states. This represents a significant improvement with respect to the exponential scaling of the vibronic basis state representation, where converged results often require the inclusion of several discrete vibrational energy states for each mode included in the representation, as opposed to the highly efficient moving basis set of coherent states for describing the quantum vibrational states within the CSDM description. Moreover, the trajectory weight functional defined in Eq. (15) is solely dependent on a system-bath couplinglike term that is linear in bath position and momentum variables. Since this term is proportional to the system-bath coupling strength, we anticipate numerical convergence issues in the limit of very strong coupling. However, we have not encountered this issue in the calculations presented here, which align well with parameter regimes relevant for the study of photosynthetic complexes.^{4,5,13–15,41} It is also worth noting that the imaginary part of the exponent for the trajectory weighting functional, $Im(\Theta t[Q\u2192,P\u2192,p\u2192,q\u2192])$, which determines an overall real weight, is inversely proportional to the frequency of a bath oscillator. In the high frequency limit, this corresponds to a trajectory weight with unit magnitude. However, if used to treat low-frequency modes, which do not align well with the approximation made to the stability factor in Eq. (10), we expect this weighting magnitude to become rather volatile with respect to bath momentum variables. In such cases, we could imagine employing iterative procedures^{9,42} or the modified Filinov filtering technique to suppress the sampling of initial bath phase space points that lead to ill-behaved trajectory weights at long times, while retaining an acceptable description of phase interference effects, to reduce the number of initial sampling points required to reach convergence.^{43,44} For low frequency modes, however, a linearized, classical-like description is often reliable, so the hybrid approach combining the CSDM propagation method developed here for higher frequency discrete modes with the PLDM approach for lower frequency continuum bath modes as demonstrated in the example model considered above should be a useful approach.

In upcoming work, we plan to explore the efficiency and applicability of the composite CSDM-PLDM propagation scheme for describing excitation energy transfer dynamics in realistic first-principle vibronic models of light-harvesting complexes,^{5} including the computation of spectroscopic signals. Moreover, the convenient trajectory weighting functions defined by the transition amplitudes [detailed under Eq. (12)] of the bath also allow for a straightforward analysis of the impact of pulse-shaping for optimal quantum control and state preparation.

## ACKNOWLEDGMENTS

We gratefully acknowledge support for this research from the National Science Foundation under Grant No. CHE-1665367. J.P. acknowledges support from the Molecular Sciences Software Institute under Grant No. ACI-1547580.