Efficient construction of generalized master equation memory kernels for multi-state systems from nonadiabatic quantum-classical dynamics

TL;DR

This paper introduces an efficient method for constructing generalized master equation memory kernels for multi-state quantum systems using quantum-classical dynamics, enabling accurate long-time dynamics with reduced computational effort.

Contribution

The authors develop an algorithm to accelerate the quantum-classical GQME approach by selectively sampling key kernel matrix elements, improving efficiency for multi-state systems.

Findings

The memory kernel decays within ~65 fs, capturing dynamics over tens of picoseconds.

The method accurately reproduces dynamics of FMO and LHCII complexes.

Quadratic scaling of trajectory number with system size improves computational feasibility.

Abstract

Methods derived from the generalized quantum master equation (GQME) framework have provided the basis for elucidating energy and charge transfer in systems ranging from molecular solids to photosynthetic complexes. Recently, the non-perturbative combination of the GQME with quantum-classical methods has resulted in approaches whose accuracy and efficiency exceed those of the original quantum-classical schemes while offering significant accuracy improvements over perturbative expansions of the GQME. Here we show that, while the non-Markovian memory kernel required to propagate the GQME scales quartically with the number of subsystem states, the number of trajectories required scales at most quadratically when using quantum-classical methods to construct the kernel. We then present an algorithm that allows further acceleration of the quantum-classical GQME by providing a way to selectively sample the kernel matrix elements that are most important to the process of interest. We demonstrate the utility of these advances by applying the combination of Ehrenfest mean field theory with the GQME (MF-GQME) to models of the Fenna-Matthews-Olson (FMO) complex and the light harvesting complex II (LHCII), with 7 and 14 states, respectively. This allows us to show that MF-GQME is able to accurately capture all the relevant dynamical time-scales in LHCII: the initial nonequilibrium population transfer on the femtosecond time-scale, the steady state-type trapping on the picosecond time-scale, and the long time population relaxation. Remarkably, all of these physical effects spanning tens of picoseconds can be encoded in a memory kernel that decays after only $\sim$65 fs.

Equations

$$\hat{H}={\hat{H}}_s+{\hat{H}}_b+{\hat{H}}_{sb},$$

$${\hat{H}}_{sb}=\sum _{j=1}^{N_{sb}}{\hat{S}}_j\otimes {\hat{\Gamma }}_j,$$

$${\hat{\rho }}_s(t)=\sum _{n=1}^{N_s^2}[{\rho }_s{]}_n(t){\hat{A}}_n,$$

$${\left[{\hat{\rho }}_s(t)\right]}_n=\mathrm{Tr}\left\{\hat{\rho }(0)e^{i\mathcal{L}t/\hslash }{\hat{A}}_n\right\}=\mathrm{Tr}\left\{\hat{\rho }(t){\hat{A}}_n\right\}.$$

$${\dot{\rho }}_s(t)={\rho }_s(t)\mathcal{X}-{\int }_0^{t}d\tau {\rho }_s(t-\tau )\mathcal{K}(\tau ),$$

$$\mathcal{K}(t)={\mathcal{K}}_1(t)-i{\int }_0^{t}d\tau {\mathcal{K}}_3(t-\tau )\mathcal{K}(\tau ).$$

$${\mathcal{K}}_3(t)$$

$${\mathcal{K}}_1(t)$$

$$\displaystyle\qquad\qquad\qquad\quad+c^{(j,+)}q_{1}^{(jk,-)}(t)c^{(k,-)}\Big{]}.$$

$$c_{mn}^{(j,\pm )}=\mathrm{T}{\mathrm{r}}_{\mathrm{s}}\left\{{\hat{A}}_m^{†}{\left[{\hat{S}}_j,{\hat{A}}_n\right]}_{\pm }\right\},$$

$$q_{3,mn}^{(j,\pm )}(t)$$

$$q_{1,mn}^{(jk,\pm )}(t)$$

$$⟨\left|\delta {\rho }_s\right|⟩=\frac{1}{t_{max}{N}_s^2}\sum _{n=1}^{N_s^2}{\int }_0^{t_{max}}d\tau \left|[\delta {\rho }_s(\tau ){]}_n\right|,$$

$$⟨\left|\delta {\dot{\rho }}_s\right|⟩=\frac{1}{t_{max}{N}_s^2}\sum _{n=1}^{N_s^2}{\int }_0^{t_{max}}d\tau \left|[\delta {\dot{\rho }}_s(\tau ){]}_n\right|,$$

$$\left|\delta {\dot{\rho }}_s(t)\right|\lesssim \frac{1}{t_{max}}\sum _m{\int }_0^{t_{max}}dt[ϵ(t){]}_m$$

$$[ϵ(t){]}_m=\frac{1}{\sqrt{N_m}}{\int }_0^{t}d{t}^{\prime }\sum _{\beta =1}^{N_s^2}\left|[{\rho }_s(t-t^{\prime }){]}_m\right|[{\sigma }_{K_1}(t^{\prime }){]}_{m\beta },$$

$$a_m={\int }_0^{t_{max}}[ϵ(t){]}_m,$$

$$w_m=\frac{a_m}{{\sum }_k{a}_k}.$$

$$N_m^{(j)}=w_m\times j\times \Delta N.$$

$${\hat{H}}_s=\sum _{k=1}^N{E}^{(k)}|k⟩⟨k|+\sum _{k\ne k^{\prime }}{\Delta }_{k{k}^{\prime }}|k⟩⟨k^{\prime }|,$$

$$\begin{array}{rl}{\hat{H}}_{sb}& =\sum _{k=1}^{N_s}\sum _{j=1}^{N_{osc}}{c}_{j,k}{\hat{R}}_{j,k}\otimes |k⟩⟨k|\\ & \equiv \sum _{k=1}^{N_s}{\hat{\gamma }}_k\otimes |k⟩⟨k|,\end{array}$$

$$\hat{H}_{b}=\frac{1}{2}\sum_{k=1}^{N_{s}}\sum_{j=1}^{N_{osc}}\Big{[}\hat{P}_{j,k}^{2}+\omega^{2}_{j,k}\hat{R}^{2}_{j,k}\Big{]},$$

$$J_k(\omega )=\frac{\pi }{2}\sum _{j=1}^{N_{osc}}\frac{c_{j,k}^2}{{\omega }_{j,k}}\delta (\omega -{\omega }_{j,k}).$$

$$J(\omega )=\frac{2\lambda {\omega }_c\omega }{{\omega }_c^2+{\omega }^2},$$

$$J(\omega )=\sum _{k=1}^{N_{peaks}}\frac{{\omega }_{c,k}{\lambda }_k\omega }{{\omega }_{c,k}^2+(\omega +{\Omega }_k{)}^2}+\frac{{\omega }_{c,k}{\lambda }_k\omega }{{\omega }_{c,k}^2+(\omega -{\Omega }_k{)}^2}$$

$${\omega }_j={\omega }_c\tan \left[\frac{\pi }{2N_{osc}}\left(j-\frac{1}{2}\right)\right].$$

$$\begin{split}j-\frac{1}{2}=\sum_{k=1}^{N_{peaks}}\frac{N_{osc}\lambda_{k}}{\pi\lambda_{tot}}\Bigg{[}&\arctan\left(\frac{\omega_{j}+\Omega_{k}}{\omega_{c,k}}\right)\\ &\ +\arctan\left(\frac{\omega_{j}-\Omega_{k}}{\omega_{c,k}}\right)\Bigg{]},\end{split}$$

$$c_j={\omega }_j\sqrt{\frac{2{\lambda }_{tot}}{N_{osc}}}.$$

$$\begin{split}[\mathcal{K}_{1}(t)]_{mn}&=-\frac{1}{\hbar^{2}}\mathrm{Tr}\Big{[}(\mathcal{L}_{sb}\hat{\rho}_{b}^{eq}\hat{A}_{m}^{\dagger})e^{i\mathcal{L}t/\hbar}(\mathcal{L}_{sb}\hat{A}_{n})\Big{]}\\ &=-\frac{1}{2\hbar^{2}}\sum_{j,k}^{N_{sb}}\Bigg{\{}\mathrm{Tr}\Big{[}[\hat{\Gamma}_{j},\hat{\rho}_{b}^{eq}]_{-}[\hat{S}_{j},\hat{A}_{m}^{\dagger}]_{+}e^{i\mathcal{L}t/\hbar}[\hat{S}_{k},\hat{A}_{n}]_{-}\hat{\Gamma}_{k}\Big{]}+\mathrm{Tr}\Big{[}[\hat{\Gamma}_{j},\hat{\rho}_{b}^{eq}]_{+}[\hat{S}_{j},\hat{A}_{m}^{\dagger}]_{-}e^{i\mathcal{L}t/\hbar}[\hat{S}_{k},\hat{A}_{n}]_{-}\hat{\Gamma}_{k}\Big{]}\Bigg{\}}\\ &=\frac{1}{2\hbar^{2}}\sum_{j,k}^{N_{sb}}\sum_{r,s}^{N_{s}^{2}}\Bigg{\{}c^{(j,+)}_{m,r}\mathrm{Tr}\Big{[}[\hat{\rho}_{b}^{eq},\hat{\Gamma}_{j}]_{-}\hat{A}_{r}^{\dagger}e^{i\mathcal{L}t/\hbar}\hat{A}_{s}\hat{\Gamma}_{k}\Big{]}+c_{mr}^{(j,-)}\mathrm{Tr}\Big{[}[\hat{\Gamma}_{j},\hat{\rho}_{b}^{eq}]_{+}\hat{A}_{r}^{\dagger}e^{i\mathcal{L}t/\hbar}\hat{A}_{s}\hat{\Gamma}_{k}\Big{]}\Bigg{\}}\mathbf{c}_{sn}^{(k,-)}\\ &=\frac{1}{\hbar^{2}}\sum_{j,k}^{N_{sb}}\sum_{r,s}^{N_{s}^{2}}\Big{\{}c_{mr}^{(j,+)}[q_{1}^{(jk,-)}]_{rs}(t)-c_{mr}^{(j,-)}[q_{1}^{(jk,+)}]_{rs}(t)\Big{\}}c_{sn}^{(k,-)}.\end{split}$$

$$\begin{split}[\mathcal{K}_{3}(t)]_{mn}&=-\frac{1}{\hbar}\mathrm{Tr}\Big{[}(\mathcal{L}_{sb}\hat{\rho}_{b}^{eq}\hat{A}_{m}^{\dagger})e^{i\mathcal{L}t/\hbar}\hat{A}_{n}\Big{]}\\ &=\frac{1}{\hbar}\sum_{j}^{N_{sb}}\sum_{r}^{N_{s}^{2}}\Big{\{}c_{mr}^{(j,+)}[q_{3}^{(j,-)}]_{rn}(t)+c_{mr}^{(j,-)}[q_{3}^{(j,+)}]_{rn}(t)\Big{\}},\end{split}$$

$$c^{(j,\pm)}_{mn}=\mathrm{Tr}_{s}\Big{[}\hat{A}_{m}^{\dagger}[\hat{S}_{j},\hat{A}_{n}]_{\pm}\Big{]}.$$

$$[{\hat{S}}_j,{\hat{A}}_n{]}_{\pm }$$

Full text

Non adiabatic dynamics, quantum dynamics, light harvesting complex, electronic energy transfer

I Introduction

Theoretical modelling of processes ranging from energy transfer in photosynthetic light harvesting complexes and organic electronics to proton coupled electron transfer, require the development of methods that allow for the accurate and efficient calculation of the dynamics of open quantum systems, i.e., a subsystem of interest coupled to an external environment. Generalized quantum master equations (GQMEs) provide a formally exact way of rewriting the equation of motion of a set of degrees of freedom of particular interest via the use of the appropriate projection operator Nakajima (1958); Zwanzig (1960); Berne and Pecora (1976). For example, in open quantum systems, the projection operator is frequently chosen to encompass all the states of the electronic subsystem, although this choice need not align with simple partitions at the level of the Hamiltonian, i.e., the projector could contain only a subset of electronic states or even include some nuclear degrees of freedom. The effect of the remaining degrees of freedom on the dynamics of those of interest is encoded through the memory kernel.Nakajima (1958); Zwanzig (1960); Mori (1965) However, evaluating the memory kernel requires overcoming the complications that arise from the presence of the projected propagator.

A common approach to avoid the projected propagator is to treat the memory kernel using second order perturbation theory. This leads to methods including Markovian and non-Markovian Redfield theory and its variants Bloch (1957); Redfield (1965); Chang and Skinner (1993); Zhang et al. (1998), the noninteracting blip approximation (NIBA) Leggett et al. (1987); Dekker (1987), Förster theory Förster (1960), and the polaronic quantum master equation Jang et al. (2008); Jang (2011); Nazir (2009); McCutcheon and Nazir (2011), which have been used to treat a wide variety of multi-state electronic energy transfer systems. In these approaches, the memory kernel depends only on correlation functions of the isolated system and bath, and hence can be straightforwardly generated even for systems with many subsystem states.Breuer and Petruccione (2007); Nitzan (2006) However, these approaches are limited in their applicability to systems where the perturbative parameter remains small.

More recently, it has been shown that the projected propagator can be exactly removed by employing the Dyson identity, leading to a representation of the full (non-perturbative, non-Markovian) memory kernel in terms of unprojected correlation functions of the system and bath.Shi and Geva (2003) This allows one to gain accurate results over a much wider range of regimes than perturbative master equations, albeit at the expense of requiring cross-correlation functions of system and bath operators to generate the memory kernel. Evaluating these cross-correlation functions using trajectory-based quantum-classical approaches with the GQME formalism has been shown to give considerable increases in both accuracy and efficiency compared to using the same quantum-classical method alone.Shi and Geva (2004); Kelly and Markland (2013); Kelly, Brackbill, and Markland (2015); Pfalzgraff, Kelly, and Markland (2015); Montoya-Castillo and Reichman (2016); Kelly et al. (2016); Montoya-Castillo and Reichman (2017)

However, while these studies showed the efficacy of such a combination in systems containing two subsystem states, for a system with $N_{s}$ subsystem states, the reduced density matrix (RDM) has $N_{s}^{2}$ elements and hence the size of the memory kernel needed to propagate it grows as $N_{s}^{4}$. One might therefore imagine that combining quantum-classical methods with the GQME would be limited to problems with small numbers of subsystem states. Here, we demonstrate that for quantum-classical methods that treat the subsystem using wavefunctions, e.g., Ehrenfest mean field theory (MFT) McLachlan (1964); Tully (1998); Stock (1995) and fewest switches surface hopping (FSSH) Tully and Preston (1971); Tully (1990); Hammes-Schiffer (1996), that the full memory kernel can be obtained using only $\mathcal{O}(N_{s}^{2})$ trajectories.

We then show how, in cases where one is interested in a particular (or a small set of) subsystem initial conditions, one can further improve the efficiency by focusing resources into accurately capturing the matrix elements of the memory kernel that most strongly determine the evolution of a given initial condition. We achieve this by noting that, while the fully converged memory kernel contains the information necessary propagate the GQME from any subsystem initial condition contained in the projector, often many of these may not be of interest. Here we introduce an algorithm that exploits this observation to identify and selectively converge the kernel elements that are most important to capture the relaxation from a given subsystem initial condition. This allows one to minimize the computational effort expended on the elements that only exert a minor influence on the process of interest.

We then illustrate the scalability of the quantum-classical GQME to multi-state systems by applying the mean field GQME (MF-GQME) method Kelly, Brackbill, and Markland (2015) to treat the 7-state Fenna-Matthews-Olsen (FMO) complex Adolphs and Renger (2006); Ishizaki and Fleming (2009a) and the 14-state light harvesting complex II (LHCII).Novoderezhkin et al. (2005); van Grondelle and Novoderezhkin (2006); Novoderezhkin, Marin, and van Grondelle (2011) We show that MF-GQME gives quantitative accuracy in regimes where direct MFT and modified Redfield theory fail, while remaining less expensive than a direct application of MFT even in these larger systems. This development thus opens the door to the application of quantum-classical GQME approaches to multi-state systems coupled to atomistic environments that lie beyond linear system-bath coupling and harmonic environments.Pfalzgraff, Kelly, and Markland (2015)

II Theory

To show how the GQME formalism can be efficiently used with trajectory-based quantum dynamics methods in multi-state systems, here we give the general expression for the memory kernel for a system with arbitrary system-bath coupling. By analyzing the form of this general expression, we demonstrate that only $\mathcal{O}(N_{s}^{2})$ trajectories are required to calculate the memory kernel in the MF-GQME method in a system with $N_{s}$ subsystem states.

To begin, we consider a Hamiltonian of the general form,

[TABLE]

where $\hat{H}_{s}$ is the isolated subsystem Hamiltonian, $\hat{H}_{b}$ the isolated bath Hamiltonian, and $\hat{H}_{sb}$ contains all the terms that couple the system and the bath. For any $\hat{H}_{sb}$, one can write it such that it is in the form of a sum of direct products of system and bath operators,

[TABLE]

where $\hat{S}_{j}$ is a pure subsystem operator, $\hat{\Gamma}_{j}$ is a pure bath operator, and $N_{sb}$ is the number of terms in the sum.

Here we consider initial conditions for the full system density operator $\hat{\rho}(0)$ that are of the spectroscopic form, $\hat{\rho}(0)=\hat{\rho}^{eq}_{b}\otimes\hat{\rho}_{s}(0)$, where $\hat{\rho}_{s}(0)$ encodes the initial state of the subsystem and $\hat{\rho_{b}}^{eq}$ is the equilibrium bath density operator, $\hat{\rho}^{eq}_{b}=e^{-\beta\hat{H}_{b}}/\mathrm{Tr}_{b}\left\{e^{-\beta\hat{H}_{b}}\right\}$. While this choice of initial condition simplifies the expressions for the memory kernel, we emphasize that the scaling arguments presented here do not depend on this choice and are straightforward to generalize to any factorizable initial condition.

The quantity of interest is the RDM, $\hat{\rho}_{s}(t)=\mathrm{Tr}_{b}\{\hat{\rho}(t)\}$, from which any subsystem observable, such as the population relaxation or electronic spectra, can be generated. In a given basis, one can write the RDM as

[TABLE]

where $\hat{A}_{n}\in\left\{\ket{\alpha}\bra{\alpha^{\prime}}\right\}$, which is a set of $N_{s}^{2}$ subsystem operators that span the subsystem Liouville space, and $[\rho_{s}]_{n}$ denotes a particular matrix element of the RDM. These matrix elements can be expressed as quantum mechanical expectation values,

[TABLE]

In the Mori-Nakajima-Zwanzig formalism, the dynamics of the RDM can be obtained by integrating the GQME Zwanzig (1960); Nakajima (1958); Mori (1965); Berne and Pecora (1976); Breuer and Petruccione (2007),

[TABLE]

where $\rho_{s}(t)$ is a row vector with $N_{s}^{2}$ time-dependent elements and $\mathcal{X}$ and $\mathcal{K}(t)$ are matrices of size $N_{s}^{2}\times N_{s}^{2}$, and $\dot{\rho}_{s}(t)=\frac{d}{dt}\rho_{s}(t)$. The static matrix $\mathcal{X}_{mn}$ = $\frac{i}{\hbar}\mathrm{Tr_{s}}\left\{\hat{A}_{m}^{\dagger}\mathcal{L}_{s}\hat{A}_{n}\right\}$ ensures that the GQME exactly recovers the free subsystem evolution when the memory kernel, $\mathcal{K}(t)$, is zero, i.e. when the subsystem and bath do not interact. The memory kernel, which encodes the effect of the bath on the subsystem dynamics, can be constructed using correlation functions in the full system and bath Hilbert space. Using the Dyson decomposition of the propagatorEvans and Morriss (2007) one obtains,Shi and Geva (2003); Montoya-Castillo and Reichman (2016)

[TABLE]

Here $\mathcal{K}_{3}$ and $\mathcal{K}_{1}$ are auxiliary kernels, which can be obtained using one’s choice of dynamical method, and the subscripts 1 and 3 are used to be consistent with earlier work.Shi and Geva (2003, 2004); Kelly and Markland (2013); Kelly, Brackbill, and Markland (2015); Pfalzgraff, Kelly, and Markland (2015); Montoya-Castillo and Reichman (2016) For any system-bath coupling, $\mathcal{K}_{3}$ and $\mathcal{K}_{1}$ can be written as linear combinations of correlation functions of system and bath operators of the form (see Appendix A),

[TABLE]

Here, like the memory kernels, $c^{(j,\pm)}$, $q_{3}^{(j,\pm)}(t)$, and $q_{1}^{(jk,\pm)}(t)$ are matrices of size $N_{s}^{2}\times N_{s}^{2}$. The transformation matrices $c^{(j,\pm)}$ are time-independent and simple to evaluate analytically, with their elements given by,

[TABLE]

where $[\cdot,\cdot]_{-}$ denotes the commutator and $[\cdot,\cdot]_{+}$ the anti-commutator. The matrix elements of the correlation functions $q_{3}^{(j,\pm)}(t)$ and $q_{1}^{(jk,\pm)}(t)$ are (see Appendix A)

[TABLE]

Equations (7) and (8) provide a formal basis for the calculation of the memory kernel, in which one obtains the auxiliary kernels $\mathcal{K}_{3}$(t) and $\mathcal{K}_{1}$(t) by evaluating $c^{(j,\pm)}$ exactly and calculating $q_{3}^{(j,\pm)}(t)$ and $q_{1}^{(jk,\pm)}(t)$ using one’s choice of exact or approximate dynamical method. These quantities are then used to to obtain the auxiliary kernels $\mathcal{K}_{1}$ and $\mathcal{K}_{3}$ using Eqs.(7)–(8), and the full memory kernel is obtained using Eq. (6) (see Appendix B) In Sec. II.1, we use these expressions to analyze the scaling with the number of subsystem states.

II.1 Scaling of the MF-GQME method with number of subsystem states

In the MF-GQME method, one approximates the correlation functions defined in Eqs. (10) and (11) using MFT (see Appendix C). The memory kernel can then be used to solve the GQME for the dynamics of the RDM, given by Eq. (5). In practice, the computational cost of generating the full memory kernel using Eqs. (6)–(8) and integrating the GQME is negligible compared to the cost of generating the MFT trajectories necessary to calculate $\mathcal{K}_{3}$ and $\mathcal{K}_{1}$. The computational cost is therefore determined by the number of MFT trajectories required to evaluate $q_{3}^{(j,\pm)}(t)$ and $q_{1}^{(jk,\pm)}(t)$, and how this cost scales with the number of subsystem states.

The number of correlation functions that need to be calculated to evaluate $q_{3}^{(j,\pm)}(t)$ and $q_{1}^{(jk,\pm)}(t)$ using Eqs. (10) and (11) can at first appear formidable, as the indices $m$ and $n$ in those equations each run over $N_{s}^{2}$ elements. This might lead one to erroneously conclude that the MF-GQME method scales as $\mathcal{O}(N_{s}^{4})$, which would severely limit its applicability to treat problems with many subsystem states. However, closer inspection of these expressions in the context of MFT dynamics reveals that one only needs to sample the distinct electronic initial conditions, labelled by the index $m$, to calculate all the necessary correlation functions. This is because the MFT trajectories generated from a given electronic initial condition, $\hat{A}_{m}^{\dagger}$, provide the time-dependent observables necessary to calculate all of the correlation functions corresponding to that initial condition. Consequently, since only $N_{s}^{2}$ electronic initial conditions need to be sampled, only $\mathcal{O}(N_{s}^{2})$ quantum-classical trajectories are required. In practice, for MFT the number of initial conditions can be shown to be $\frac{1}{2}N_{s}(N_{s}+1)$ by following the procedure outlined in Appendix D, which exploits the fact that the correlation functions corresponding to a given off-diagonal initial condition and its Hermitian conjugate can be obtained from the same set of trajectories.

In general, evaluating forces dominates the overall cost, especially in atomistic systems or systems with a large number of bath degrees of freedom. However, the calculation of the correlation functions required scales as $\mathcal{O}(N_{s}^{4})$ and integration of the GQME scales as $\mathcal{O}(N_{s}^{6})$ which could dominate the cost for large $N_{s}$. However, even for the systems studied here which use very simple forms of the bath, the cost of integrating the GQME is negligible.

An additional concern regarding the scalability of the MF-GQME approach comes from amount of memory required to store all elements in $q_{3,mn}^{(j,\pm)}(t)$ and $q_{1,mn}^{(j,\pm)}(t)$. For instance, the number of time-dependent elements in $q_{1,mn}^{(j,\pm)}(t)$ can be calculated as the product of $N_{s}^{2}$ initial electronic conditions $\hat{A}_{m}^{\dagger}$, $N_{sb}$ bath operators $\hat{\Gamma}_{j}$ measured at $t=0$, and $N_{s}^{2}$ electronic operators $\hat{A}_{n}(t)$ and $N_{sb}$ bath operators $\hat{\Gamma}_{j}(t)$ measured at finite times. This corresponds to $\mathcal{O}(N_{s}^{4}\times N_{sb}^{2})$ time dependent elements that need to be stored in memory, and hence one might expect the amount of available memory to be a significant limitation. In practice, however, this problem can be easily circumvented because $q_{3}^{(j,\pm)}(t)$ and $q_{1}^{(jk,\pm)}(t)$ do not themselves need to be stored in memory. Instead, one can evaluate their contributions to $\mathcal{K}_{3}(t)$ and $\mathcal{K}_{1}(t)$ on the fly using Eqs. (7) and (8), so that only $\mathcal{K}_{3}(t)$ and $\mathcal{K}_{1}(t)$ need to be stored. In addition, because not all of the matrix elements are independent, there are only $\frac{1}{2}N_{s}^{2}(N_{s}^{2}+1)$ unique elements in each, with the remainder being complex conjugates of these elements.Shi and Geva (2003) The amount of memory required to store the kernels is therefore proportional to $N_{s}^{2}(N_{s}^{2}+1)\times N_{steps}$ where $N_{steps}$ is the number of timesteps before the kernel decays. In practice, when the memory kernel decays quickly to zero $N_{steps}$ is a small number, meaning that negligible storage is required.

II.2 Accelerating convergence of the GQME using selective sampling of kernel matrix elements

Once the memory kernel has been calculated, solving the GQME, Eq. (5), provides access to the subsystem dynamics for any factorizable initial condition of the subsystem RDM. However, in many cases one is only interested in a particular initial condition, or a set of them. For example, one may be interested in the relaxation dynamics resulting from the initial excitation of a particular chromophore in a molecular system consisting of multiple chromophores. Here it is important to make a distinction between the initial condition of interest at the RDM level, and the electronic initial conditions that need to be sampled in the MFT calculation of the memory kernel. At the memory kernel level, all separable initial conditions (e.g., initial excitation of each chromophore and coherence between chromophores) need to be sampled to construct the full memory kernel. It is this complete sampling of the initial conditions at the memory kernel level that ensures that one can use Eq. (5) to propagate the RDM dynamics from any separable initial condition.

However, as in the case mentioned above, suppose that one is only interested in propagating the RDM starting from one particular initial condition (or a small subset of initial conditions) in a system with a large number of subsystem states. In this case, there may be many pathways that do not play a role in the relaxation dynamics from the RDM initial condition of interest. Hence, it is clearly unnecessary to have accurate knowledge of the memory kernel elements that describe these unimportant relaxation pathways.

One can exploit this realization by focusing resources into converging the elements of the memory kernel that are most important to describe the relaxation process of interest. When the memory kernel elements are generated using MFT, this corresponds to identifying and selectively converging the most important correlation functions from the sets $\{q_{3}^{(j,\pm)}(t)\}$ and $\{q_{1}^{(jk,\pm)}(t)\}$ (defined in Eqs. (10) and (11)), which give rise to the memory kernel. As such, here we present an algorithm that allows one to use fewer trajectories to generate the less important correlation functions while minimizing the resulting loss in accuracy for the RDM.

Any algorithm that allows for the selective convergence of the memory kernel elements should satisfy certain requirements. In particular, after specifying the RDM-level initial condition that one wishes to propagate, the algorithm should provide a prescription for apportioning a fixed total number of quantum-classical trajectories, $N_{tot}$, among the different initial conditions that need to be sampled at the the memory kernel level. Specifically, the algorithm should provide a set of normalized weights, $\left\{w\right\}$, that determine the portion of trajectories $w_{m}N_{tot}$ assigned to a given memory kernel level initial condition, $A_{m}^{\dagger}$. These weights should be chosen such that, when the time-evolved RDM is propagated from the specified initial condition, the error that arises from the statistical error in the memory kernel is minimized. In addition, the calculation of these weights should take into account two important observations. First, correlation functions with certain memory kernel-level initial conditions may contribute more to the relevant relaxation pathways than others. Second, correlation functions arising from different memory kernel-level initial conditions may require different numbers of trajectories to achieve reasonable levels of convergence.

To calculate the optimal weights, some knowledge of the relaxation pathway and current level of convergence of the memory kernel is crucial. Here we suggest an iterative approach, where increasingly improved estimates for the memory kernel and the RDM dynamics are used to re-compute the weights at each iteration. An overview of the algorithm is shown in Fig. 1. We note that while here we present one particular form for the weights, the approach used here is general and does not depend on this particular choice.

To provide a form for the weights, $w_{m}$, we consider how a given realization of the RDM dynamics, which we call the trial RDM dynamics, deviates from its converged form, $\delta\rho_{s}(t)=\rho_{s}^{conv}(t)-\rho_{s}(t)$. Although in practice one does not know the converged dynamics ${\rho}_{s}^{conv}(t)$, considering how an unconverged dynamics may deviate from its converged form will allow us to elucidate how the convergence of the dynamics can be accelerated via selective sampling. Our goal is to reduce the absolute value of this difference from the converged dynamics for a given number of trajectories, $N_{tot}$. Because the dynamics $\rho_{s}(t)$, arise from integrating Eq. (5), one can minimize the absolute difference between the time-derivative of the trial RDM and its converged form, $\delta\rho_{s}(t)=\rho_{s}^{conv}(t)-\rho_{s}(t)$. Because this error metric is time dependent and different for each element of the density, it is convenient to consider the mean absolute deviation, which is the absolute value of $\delta\rho_{s}(t)$, averaged over time and elements of the density matrix:

[TABLE]

where $t_{max}$ is the maximum time that is of interest of the dynamics.

Because there is no error in the initial condition, (i.e., $\delta\rho_{s}(0)=0$), one can minimize the error at subsequent times by minimizing the error in the time-derivative of the RDM, $\delta\dot{\rho}_{s}(t)=\dot{\rho}_{s}^{conv}(t)-\dot{\rho_{s}}(t)=\frac{d}{dt}\delta\rho_{s}(t)$. Again it is convenient to consider the mean absolute deviation in this quantity, which is:

[TABLE]

An approximate upper bound of this quantity is (see Appendix E),

[TABLE]

where the index $m$ runs over all initial conditions. The selective sampling error metric (see also Appendix E), is approximately proportional to the error in the dynamics at time $t$ that arises from finite sampling of correlation functions with initial condition $m$. For the models considered here, this quantity is:

[TABLE]

where $N_{m}$ is the number of trajectories that have been run in all previous iterations of the algorithm with memory kernel level initial condition $m$. The quantity $[\sigma_{K_{1}}(t)]_{m\beta}$ is given by Eq. (59) in Appendix E.

Motivated by the above considerations, the computed weights in our algorithm depend on $\rho_{s}(t)$, as well as the statistical error of the correlation functions that contribute to a specific elements of the memory kernel. By distributing $N_{tot}$ trajectories over $n_{step}$ steps in increments of $\Delta N=N_{tot}/n_{step}$ to different GQME-level initial conditions, we refine our estimate for these weights iteratively at each step, as described below. We repeat this procedure either until the RDM converges to a desired threshold or until $N_{tot}$ trajectories have been run.

Using the definitions provided above, we define the un-normalized weight, $a_{m}$, as

[TABLE]

The normalized weight $w_{m}$ is then

[TABLE]

Using these weights, we then apportion trajectories to each memory kernel level initial condition $\hat{A}_{m}^{\dagger}$ in the subsequent cycle so that,

[TABLE]

where $j$ is the index of the current iteration of the algorithm (for example, in the first iteration of the algorithm $j=1$). $N_{m}^{(j)}$ is therefore the number of trajectories that have been assigned to initial condition $m$ at step $j$. The algorithm can be summarized as follows:

Initialize the algorithm by setting $\mathcal{K}(t)=0$ and $[\sigma_{\mathcal{K}_{1}}]_{m\beta}=1$. 2. 2.

Use $\mathcal{K}(t)$ to intergrate the GQME (Eq. (5)) for some time range $t_{max}$, yielding $\rho_{s}(t)$ 3. 3.

Use the current guess for $\mathcal{K}(t)$ and $\rho_{s}(t)$ to compute the weights, $w_{m}$, corresponding to a given memory kernel-level initial condition $\hat{A}_{m}^{\dagger}$ using Eqs. (15)–(17). 4. 4.

Distribute a number of MFT trajectories, $\Delta N$, among all memory kernel-level initial conditions according to Eq. (18). 5. 5.

Using the correlation functions $q_{3}^{(j,\pm)}(t)$ and $q_{1}^{(jk,\pm)}(t)$ obtained from all MFT trajectories that have been run so far, rebuild $\mathcal{K}(t)$ using Eqs. (6)–(11). 6. 6.

Calculate $\sigma_{\mathcal{K}_{1}}(t)$ using the standard deviations in the correlation functions $q_{1}^{(jk,\pm)}(t)$ (see Appendix E). 7. 7.

Calculate the mean absolute deviation between the previous and the current iteration of $\rho_{s}(t)$. If the maximum number of trajectories, $N_{tot}$ is exceeded, or if the mean absolute deviation falls below the predetermined threshold, stop. Otherwise, return to step 2.

A plot of the number of trajectories assigned to each initial condition for the LHCII model discussed in Section IV is available in the SI.

III Simulation Details

III.1 Hamiltonian

Here, we apply MF-GQME to the Frenkel exciton Hamiltonian. Adolphs and Renger (2006); Ishizaki and Fleming (2009a) This model Hamiltonian forms the basis of many studies of electronic energy transfer in photosynthetic light harvesting systems May and Kühn (2011); Jang and Mennucci (2018). In addition to providing physical insight into these systems, the Frenkel exciton model’s simplicity, i.e., its bilinear coupling to harmonic baths, allows for exact results to be calculated using, for example, the hierarchical equations of motion (HEOM) approach.Ishizaki and Fleming (2009b) Although GQMEs, and in particular MF-GQME, are not limited to any specific form of the Hamiltonian, the availability of exact results for this model system allows us to benchmark the accuracy and efficiency of MF-GQME for multi-state problems. The Frenkel exciton Hamiltonian is of the general form given in Eq. (1), with the subsystem Hamiltonian taking the form,

[TABLE]

where $E^{(k)}$ is the energy of diabatic state $\ket{k}$ and $\Delta_{kk^{\prime}}$ is the electronic coupling between states $\ket{k}$ and $\ket{k^{\prime}}$. Each diabatic state is bilinearly coupled to a bath of $N_{osc}$ independent harmonic oscillators,

[TABLE]

where $j$ indexes the bath degrees of freedom, $c_{j,k}$ is the coupling constant for a given bath oscillator, $\hat{R}_{j,k}$ is the position operator for that oscillator, and in the final equality we define $\hat{\gamma}_{k}$ to be the bath part of the system-bath coupling to state $\ket{k}$. The bath therefore consists of $N_{s}\times N_{osc}$ independent harmonic oscillators:

[TABLE]

where $\hat{P}_{j,k}$ is the momentum operator of the bath mode, and $\omega_{j,k}$ is its frequency, and the Hamiltonian is given in mass-weighted coordinates and momenta. The coupling of each diabatic state to its local bath is characterized by its spectral density,

[TABLE]

For the FMO model, it is common to assume that all spectral densities are equivalent and take the Debye form,

[TABLE]

where $\lambda$ is the reorganization energy that defines the strength of the system-bath coupling and $\omega_{c}$ is the inverse of the characteristic timescale of the bath. For the LHCII model, we adopt a sum of shifted Drude-Lorentz peaks Kreisbeck, Kramer, and Aspuru-Guzik (2014):

[TABLE]

where we use the parameters of Ref. Kreisbeck, Kramer, and Aspuru-Guzik, 2014: $N_{peaks}=7$, $\omega_{c}^{-1}=\left\{30,1400,1000,1400,1000,600,1000\right\}$ fs, $\lambda=\left\{130,6,18,6,16,48,17\right\}\mathrm{cm^{-1}}$ , and $\Omega=\left\{240,297,342,388,518,745,915\right\}\mathrm{cm^{-1}}$.

III.2 Computational Methods

For the FMO model, the spectral density was discretized with frequencies

[TABLE]

For the LHCII model the spectral density was discretized by numerically solving for $\omega_{j}$ that satisfies

[TABLE]

where $\lambda_{tot}=\sum_{k}\lambda_{k}$. For both Eqs. (25) and (26), $\omega_{j}$ is obtained by taking $j=1,2,\ldots,N_{osc}$. For both models, the couplings of the discrete oscillators are given by:

[TABLE]

For Ehrenfest trajectories, the equations of motion were solved using a split evolution method, where at each dynamical timestep, the bath is evolved through a half timestep of $\frac{\Delta t}{2}$ using mean field theory (see Appendix C) with the velocity Verlet algorithm. Swope et al. (1982) The system was then evolved for a full timestep $\Delta t$ by numerically diagonalizing $H_{s}+H_{sb}(R,P)$, holding the bath positions and momenta constant. The bath was then evolved for another half timestep of $\frac{\Delta t}{2}$ using velocity Verlet, keeping the state of the electronic subsystem constant. For the FMO model, the timestep used was $\Delta t=$0.25 fs, while for LHCII the timestep was $\Delta t=$0.05 fs.

The partial memory kernels were constructed using Eqs. (7)-(11), with the necessary correlation functions in Eqs. (10) and (11) approximated using mean field theory as described in Appendix C. The memory kernel was calculated at a lower time resolution (0.5 fs) than the direct MFT dynamics used to generate the underlying correlation functions, Eqs. (10) and (11), which was observed not to change the resulting GQME dynamics. Additional details on the specific form these equations in the context of the Frenkel exciton model are given in Appendix A. Plots of select memory kernel matrix elements that relate to transitions between the states that have significant populations for FMO and LHCII are available in the SI.

Typically the memory kernel decays to zero after a short amount of time and hence can be truncated (set equal to zero) at all times greater than this time. For the FMO model, this kernel cutoff time, $\tau_{c}$, was chosen by increasing the cutoff until the dynamics ceased to change, with $\tau_{c}$ between 125 and 300 fs giving very similar results. In practice, larger values of $\tau_{c}$ require more trajectories to converge. For the dynamics of the FMO model shown in this work we used $\tau_{c}$ = 150 fs. For the LHCII model, the dynamics remain unchanged for $\tau_{c}$ between 50 and 65 fs. To stably integrate the MF-GQME dynamics with cutoff times greater than 65 fs requires a large number of trajectories to achieve stable dynamics. This is because the LHCII complex is in a more strongly coupled regime, which gives rise to highly oscillatory memory kernels that require a very large number of trajectories to converge. Such oscillatory memory kernels have been observed for the dynamics of the spin-boson model with the Debye spectral density in intermediate to strong coupling regimes Montoya-Castillo and Reichman (2016). For all LHCII dynamics shown, we used $\tau_{c}$ = 65 fs. The instabilities observed for LHCII might be able to be alleviated by employing alternative closures,Zhang, Ka, and Geva (2006); Montoya-Castillo and Reichman (2016); Mulvihill et al. (2019) or by using higher level quantum-classical methods to generate the memory kernel.Kelly and Markland (2013); Kelly et al. (2016)

After calculating the partial memory kernels, the GQME, Eq. (5), was integrated using Heun’s method Suli and Mayers (2003), with the integral in Eq. (5) approximated using the trapezoid rule.

Because the sampling of bath initial conditions from which the correlation functions are generated is stochastic, there is variability in how many trajectories are required to converge the memory kernel for a given realization of uniform or selective sampling. To evaluate our selective sampling algorithm on the LHCII model, we performed 50 separate realizations of the selective sampling algorithm, as well as 50 separate realizations of uniform sampling. All GQME dynamics shown are from the realization with the median error in the dynamics, where error in the dynamics is quantified by the mean absolute deviation in the subsystem density matrix from its converged form, Eq. (12). The converged form for LHCII corresponds to 65 million total trajectories, apportioned equally across all initial conditions. A plot that shows the distribution of errors of different realizations of the selective sampling procedure as a function of number of trajectories added can be found in the Supporting Information (SI). The selective sampling procedure was performed by adding a total of $\Delta N=100$ trajectories at each iteration.

IV Results and Discussion

Here we apply the MF-GQME method to electronic energy transfer in Frenkel exciton models of the FMO and LHCII light harvesting complexes. The FMO model consists of $N_{s}=7$ subsystem states whereas $N_{s}=14$ for LHCII. Therefore, for FMO there are $7^{4}=2401$ time dependent matrix elements in the memory kernel, while for LHCII there are $14^{4}=38,416$. However, as described in Sec. II.1, all of the matrix elements can be obtained using very short simulations starting from 28 initial conditions for FMO or 105 distinct initial conditions for LHCII.

Figure 2 compares the direct MFT and MF-GQME dynamics to exact HEOM resultsKreisbeck, Kramer, and Aspuru-Guzik (2014) for the population relaxation following an initial excitation for the FMO model at low temperature in a regime of high nonadiabaticity ($\omega_{c}\gtrsim\Delta_{kk^{\prime}}$ for most pairs of states). As observed in previous work for the FMO model Berkelbach, Markland, and Reichman (2012), direct MFT gives qualitatively correct coherent oscillations at short times, but fails to recover the correct long-time populations. This can be explained by the high nonadiabaticity and the low temperature bath, both of which are problematic for MFT. In contrast, MF-GQME produces highly accurate results, consistent with previous observations that this scheme offers significant improvements over direct MFT in nonadiabatic regimes. Kelly, Brackbill, and Markland (2015); Pfalzgraff, Kelly, and Markland (2015); Montoya-Castillo and Reichman (2016)

For this nonadiabatic problem, the memory kernel decays by a cutoff time of $\tau_{c}=0.15$ ps, which is considerably shorter than the 1-2 ps lifetime of the population dynamics. However, because MF-GQME requires sampling 28 distinct initial conditions to generate the memory kernel, uniform sampling of all initial conditions in this regime, i.e., when trajectories are apportioned equally to each initial condition, requires a total of $10^{5}$ total trajectories ($\sim$3,500 per initial condition) to converge the memory kernel. This is in contrast to the $\sim 10^{4}$ trajectories required to obtain a similar level of convergence with direct MFT. Hence, despite the shorter lifetime of the kernel, the larger number of trajectories required to generate the kernel means that the total cost of generating the MF-GQME dynamics shown in Fig. 2 is $\sim$1.5 times greater than that of MFT. Moreover, because the memory kernel provides the dynamics for all factorizable initial conditions of the subsystem RDM, the dynamics resulting from the two distinct initial subsystem excitations shown in panels (b.) and (d.) of Fig. 2 are generated using the same memory kernel. In contrast, panels (a.) and (c.) of Fig. 2, which correspond to direct MFT dynamics, require separate sets of MFT trajectories initialized from the two different subsystem initial conditions.

An additional advantage of the MF-GQME scheme is that, once the memory kernel is converged, one can obtain dynamics for arbitrarily long times by integrating Eq. (5). In practice, the GQME propagation can be done at negligible computational cost compared to the costs of generating MFT trajectories. Figure 3 shows the long time population dynamics for direct MFT and MF-GQME following initial excitation of state 6 in the same parameter regime as Fig. 2. Because of the prohibitive cost of converging the exact dynamics for such long times, the arrows in Fig. 3 correspond to the Boltzmann equilibrium populations that would be expected based on the subsystem Hamiltonian, $e^{-\beta\hat{H}_{s}}$. As expected, due to the well-known detailed balance problems in methods that rely on a mean field description of the interaction between the electronic subsystem and the nuclei,Parandekar and Tully (2006); Tao and Miller (2010); Kelly and Rhee (2011) MFT incorrectly predicts that all states have nearly equal population at long times, in qualitative disagreement with the expected result. In contrast, MF-GQME accurately captures the expected long-time populations. Furthermore, to converge the MFT dynamics out to the 20 ps shown in Fig. 3 requires significant additional computational effort (at least 20 times more than for 1 ps), while in the GQME approach, Eq. (5) can be integrated for arbitrarily long times using the same memory kernel as Fig. 2. Because integrating the GQME for longer times incurs negligible additional cost, the dynamics shown in Fig. 3 are $\sim 13$ times cheaper when generated with MF-GQME as compared with direct MFT.

A more stringent test of the MF-GQME method is the LHCII model, which has significantly more subsystem states each of which is also more strongly coupled to the bath degrees of freedom. The combination of a large number of states, stronger system-bath coupling, and the complex form of the spectral density of LHCII presents a significant challenge to exact methods. Hence, very limited benchmarks are available using the HEOM method and these required substantial computational effort and some simplification of the spectral density. Kreisbeck, Kramer, and Aspuru-Guzik (2014) Here we use the form and parameterization of the LHCII Hamiltonian of Ref. Kreisbeck, Kramer, and Aspuru-Guzik, 2014 with the J7 spectral density from that reference, and compare against the HEOM results presented therein. These HEOM results relied on the high temperature approximation,Kreisbeck, Kramer, and Aspuru-Guzik (2014) and hence, while they are expected to be highly accurate, it remains unclear whether they are numerically exact.

Figure 4 (a.) shows the dynamics for the sum of the populations of all Chl b chromophores for the LHCII model following initial excitation of the highest energy eigenstate of the subsystem Hamiltonian. For this model the exact HEOM dynamics exhibit an initial fast (sub picosecond) decay in the population of Chl b chromophores followed by a longer time component which decays over tens of picoseconds. Direct MFT fails to capture the magnitude of the short time drop in the population of Chl b chromophores, but qualitatively reproduces the longer timescale of population transfer out of this set of chromophores. As previously observed, modified Redfield dynamics, which is a major workhorse in the modelling and investigation of the excitation dynamics in photosynthetic systems, shows erroneous rapid decay within the first 1-2 psKreisbeck, Kramer, and Aspuru-Guzik (2014). MF-GQME, however, yields accurate short and long time behavior, in good agreement with the HEOM result. This behavior is mirrored in the populations of the Chl 604a and Chl 605b states in Fig. 4 (b.) and (c.) which contribute to the Chl b decay shown in (a.). Direct MFT underestimates the transfer into Chl 604a while overestimating that into Chl 605b, while modified Redfield theory leads to spuriously fast population decay of these states. In contrast, MF-GQME accurately captures the HEOM result.

Figure 5 shows the long-time dynamics obtained from MF-GQME using the same initial condition. In addition to capturing the exact HEOM dynamics at short times MF-GQME also obtains the intermediate timescale (picosecond) trapping in Chl 604a and Chl 605b bottleneck states, and the long timescale (tens of picoseconds) relaxation to the Chl 610a-611a-612a trimer, the Chl 613a-614a dimer, and the Chl 602a-603a dimer. Novoderezhkin et al. (2005)

Because of the stronger system-bath coupling of the LHCII model in comparison to the FMO complex, converging the memory kernel using uniform sampling of all 105 distinct initial conditions would require approximately 65 million total trajectories. Although this number of trajectories is feasible for simple models, it is intractable when more sophisticated atomistic descriptions of the environment are used. The fact that LHCII consists a large number of subsystem states makes it a prime candidate to benefit from using selective sampling to accelerate convergence. As such, we now consider the computational speedups that can be obtained by applying the selective sampling algorithm introduced in Sec. II.2 to this model. To demonstrate the utility of the selective sampling scheme, we thus compare the convergence of the GQME dynamics for a fixed total number of trajectories when using the selective sampling and the uniform sampling approaches to compute the memory kernel. See the SI for a plot showing the average number of trajectories assigned to each memory-kernel level initial condition at each step of the algorithm.

Figure 6 shows the populations for the median realizations of selective and uniform sampling at 100,000 and 200,000 trajectories. Already at 100,000 trajectories the selective sampling yields qualitatively correct dynamics, while the uniformly sampled memory kernel is not numerically stable, with the populations several chromophores becoming negative by $t=1.0$ ps, and the unphysically high frequency and large amplitude oscillations seen in panel (a.) caused by numerical instability. By 200,000 trajectories, selective sampling is already in semi-quantitive agreement to the exact result, while uniform sampling produces negative populations after only $\sim$ 0.5 ps. In fact, the uniform sampling approach fully converges only after $\sim$ 65,000,000 trajectories. Hence, in this regime the selective sampling algorithm is approximately 300 times more efficient than uniform sampling.

V Conclusions

In this paper, we have shown how the MF-GQME approach can be efficiently extended to treat systems where the electronic subspace consists of many states. By analyzing the symmetries of the memory kernel, we have demonstrated that, for a system consisting of $N_{s}$ subsystem states, the memory kernel, despite having $N_{s}^{4}$ matrix elements, can be generated using only $\mathcal{O}(N_{s}^{2})$ quantum-classical trajectories. This scaling analysis is directly applicable to quantum-classical methods using wavefunction-based treatments of the electronic subsystem dynamics, such as MFT and FSSH. The same scaling arguments also apply to methods of generating the kernel using only the original dynamics Kidon, Wilner, and Rabani (2015); Montoya-Castillo and Reichman (2016); Kelly et al. (2016), and to other methods that involve a time convolution master equation such as the transfer tensor method Cerrillo and Cao (2014); Gelzinis, Rybakovas, and Valkunas (2017). This is in contrast to methodsSun, Wang, and Miller (1998); Thoss and Stock (1999); Liao and Voth (2002); Ananth and Miller (2010); Huo and Coker (2011); Hsieh and Kapral (2012); Ananth (2013); Richardson and Thoss (2013) based on the Meyer-Miller-Stock-Thoss (MMST) mapping Meyer and Miller (1979); Stock and Thoss (1997), where the initial conditions necessary to generate all memory kernel elements could be recovered by an appropriate weighting of the variables drawn from a single Gaussian distribution of MMST variables Shi and Geva (2004); Kelly et al. (2016).

To further reduce the cost of generating the memory kernel necessary to propagate an RDM subject to specific initial conditions, we have outlined a general procedure for how one can select an optimal set of quantum-classical trajectories to converge the kernel. Employing a specific realization of a selective sampling procedure on the LHCII model that focuses resources on initial conditions that most contribute to the fast the convergence of the memory kernel, we have demonstrated that one can obtain a further $2$ orders of magnitude speed-up. While the specific procedure used here is one possible choice, alternative ways of apportioning trajectories to specific initial conditions could lead to more efficient selective sampling schemes and represents a promising direction for future work.

By applying MF-GQME to models of excitation energy transfer in the FMO and LHCII complexes, which consist of $7$ and $14$ electronic states, respectively, we have shown that this method can provide accurate results for multi-state systems for both short-time dynamics and long-time populations at the computational cost of low-level quantum classical trajectory-based methods. Further, when memory kernels are short-lived in comparison to the electronic dynamics of interest, the MF-GQME approach usually costs less than a direct application of MFT. Moreover, once the memory kernel is generated up to the required cutoff time, $\tau_{c}$, the GQME can be propagated with linear scaling in time, making long-time dynamics easy to access.

To date, MF-GQME has previously proven successful in treating the dynamics of systems coupled to harmonic Kelly, Brackbill, and Markland (2015); Montoya-Castillo and Reichman (2016, 2017) or fully atomistic Pfalzgraff, Kelly, and Markland (2015) environments, both inMontoya-Castillo and Reichman (2017) and out Kelly, Brackbill, and Markland (2015); Pfalzgraff, Kelly, and Markland (2015); Montoya-Castillo and Reichman (2016) of equilibrium. In the present work, we have demonstrated that this approach is also able to accurately and efficiently address the dynamics of systems where the electronic subspace consists of many states. Indeed, we have shown that MF-GQME can provide quantitatively accurate results in LHCII where commonly invoked methods, such as Redfield theoryBloch (1957); Redfield (1965) and its modified variantChang and Skinner (1993); Zhang et al. (1998), are known to fail Kreisbeck, Kramer, and Aspuru-Guzik (2014). MF-GQME thus provides a powerful tool to investigate charge and energy transfer in the condensed phase.

Supplementary Material

See supplementary material for plots of select memory kernel matrix elements for FMO and LHCII, a plot of the distribution of errors of different realizations of the selective sampling procedure as a function of number of trajectories added, and a plot of the number of trajectories assigned to each initial condition of the LHCII model using the selective sampling approach.

Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0014437. T.E.M also acknowledges support from the Camille Dreyfus Teacher-Scholar Awards Program. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. W.C.P acknowledges support from the Melvin and Joan Lane Stanford Graduate Fellowship. We would also like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that have contributed to these research results.

Appendix A Auxiliary kernels for a multi-level system

In this section, we provide more detail for the derivation of the expressions for the auxiliary kernels $\mathcal{K}_{1}(t)$ and $\mathcal{K}_{3}(t)$ for the GQME that describes the dynamics of the subsystem RDM subject to a spectroscopic initial condition. These expressions are applicable to a general $N_{s}$ state system coupled to an arbitrary bath of the type described in Sec. II. We subsequently specialize our discussion to the Frenkel exciton model presented in Sec. III.1 given by Eqs. (19)-(21).

We start from the general form of the Hamiltonian, Eq. (1), for an $N_{s}$ state electronic system coupled to a bath, which can take a harmonic (e.g., Eq. (20)) or atomistic form, with a general form of the system-bath coupling, given by Eq. (2).

As stated in Sec. II, we are interested in the nonequilibrium RDM dynamics arising from a spectroscopic or Franck-Condon initial condition, $\hat{\rho}(0)=\hat{\rho}_{s}\otimes\hat{\rho}_{b}^{eq}$, where the initial condition of the electronic subsystem can be expanded in terms of the Liouville basis, $\hat{\rho}_{s}(0)=\sum_{n}a_{n}\hat{A}^{\dagger}_{n}$, and the bath is in local equilibrium, $\hat{\rho}_{b}^{eq}=\frac{e^{-\beta\hat{H}_{b}}}{\mathrm{Tr}[e^{-\beta\hat{H}_{b}}]}$. To derive a particular GQME, it is also necessary to choose a projection operator. Here, we choose the $N_{s}$ state system generalization of the Redfield-type projection operator, shown in Ref. Montoya-Castillo and Reichman, 2016 in the context of a two-level system. In the more conventional Nakajima-Zwanzig form of the GQME for nonequilibrium RDM dynamics Nakajima (1958); Zwanzig (1960), this projector corresponds to the commonly used Argyres-Kelley projector Argyres and Kelley (1964), $\mathcal{P}$, which performs a partial trace over the bath and multiples by the canonical density of the bath, i.e., $\mathcal{P}(\hat{O}_{b}\hat{O}_{s})=\hat{\rho}_{b}^{eq}\mathrm{Tr}_{b}[\hat{O}_{b}\hat{O}_{s}]=\hat{\rho}_{b}^{eq}\hat{O}_{s}\mathrm{Tr}_{b}[\hat{O}_{b}]$, where $\hat{O}_{b}$ and $\hat{O}_{s}$ are bath and system operators, respectively.

With this form of the projection operator, and assuming that $\mathrm{Tr}[\hat{\rho}_{b}^{eq}\hat{\Gamma}_{j}]=0$ (which can always be enforced with no loss of generality by redefining the subsystem and bath parts of the Hamiltonian: $\hat{H}_{s}\rightarrow\hat{H}_{s}+\mathrm{Tr}[\hat{\rho}_{b}^{eq}H_{sb}]$ and $\hat{H}_{sb}\rightarrow\hat{H}_{sb}-\mathrm{Tr}[\hat{\rho}_{b}^{eq}H_{sb}]$), it is possible to obtain the general forms for the auxiliary kernels. Focusing first on $\mathcal{K}_{1}(t)$,

[TABLE]

Using similar manipulations, one can show that,

[TABLE]

where the elements of the dynamical matrices, $q_{1}^{(jk,\pm)}$ and $q_{3}^{(j,\pm)}$, are given by Eqs. (11) and (10), the elements of the static transformation matrix $c^{(j,\pm)}$ take the form,

[TABLE]

and we have used the fact that

[TABLE]

Equations (28)-(29) are general, applicable to an $N_{s}$ state system coupled to an arbitrary bath with system-bath coupling given by Eq. (2), and subject to a spectroscopic initial condition. We emphasize that such systems include electronic states coupled to harmonic or atomistic baths and system-bath interactions beyond linear coupling.

We now specialize our discussion to the Frenkel exciton Hamiltonian, Eqs. (19)-(21), where the bath is harmonic and can be decomposed into contributions that are local, in addition to the coupling being bilinear and local. Given the Frenkel exciton form of the Hamiltonian, the matrix elements of the static transformation matrix are given by,

[TABLE]

where $\delta$ is the Kronecker delta.

Appendix B Obtaining $\mathcal{K}(t)$ from $\mathcal{K}_{1}(t)$ and $\mathcal{K}_{3}(t)$ by solving a system of linear equations

Here we describe a procedure for obtaining $\mathcal{K}$ from $\mathcal{K}_{1}$ and $\mathcal{K}_{3}$ by discretizing the self-consistent expansion and solving a system of linear equations, hence avoiding an iterative solution. The self consistency relation for the memory kernel, Eq. (6) can be written in discrete time as,

[TABLE]

where we have discretized the integral using the trapezoid rule. Here $t_{0},t_{1},\ldots,t_{n}$ are equally spaced discrete times, and hence $t_{j+1}-t_{j}\equiv\Delta t$, and $t_{m}-t_{j}=t_{m-j}$ ($m>j$). Equation (33) can be solved to obtain an explicit expression for $\mathcal{K}(t_{n})$ in terms of $\mathcal{K}_{1}$, $\mathcal{K}_{3}$, and $\mathcal{K}(t_{j})$ ($0\leq j\leq n-1$),

[TABLE]

[TABLE]

where $n>0$ and $\mathbf{1}$ is the identity matrix. The above equation allows $\mathcal{K}(t_{n})$ ($n>0$) to be obtained from the $\mathcal{K}(t_{j})$ ($j<n$) using a single, non-iterative step where one solves for the unknown matrix $\mathcal{K}(t_{n})$ in Eq. (34). In practice, the system of linear equations implied by Eq. (34) was solved by computing and storing the Cholesky factorization of $[\mathbf{1}+i\frac{\Delta t}{2}\mathcal{K}_{3}(t_{0})]$, which can be used to solve for $\mathcal{K}(t_{n})$.

Appendix C MFT treatment of the dynamical matrices, $q_{1}^{(jk,\pm)}(t)$ and $q_{3}^{(j,\pm)}(t)$

Here, we outline how one can calculate the dynamical matrices, $q_{1}^{(jk,\pm)}(t)$ and $q_{3}^{(j,\pm)}(t)$ in Eqs. (11) and (10) using the Ehrenfest MFT method. To do this, we first review the partial transformation to Wigner phase space with respect to bath variables of a general correlation function. We then show the general MFT approximation to the dynamical matrices corresponding to a general $N_{s}$-state system, and specialize our discussion to the Frenkel exciton model. Finally, we provide the equations of motion for the subsystem and bath variables.

To treat a general correlation function of a system which can be separated into a subsystem coupled to a bath,

[TABLE]

using the Ehrenfest MFT method, we first perform a partial Wigner transformation Wigner (1932); Hillery et al. (1984) with respect to the bath variables. Here, $\hat{\rho}$ encodes the initial condition of the full system, and $\mathcal{A}$ and $\mathcal{B}$ are arbitrary operators. Upon performing the partial Wigner transformation, we can rewrite the correlation function as follows,

[TABLE]

where $X=(R,P)$ are the classical coordinates, $R$, and momenta, $P$, of the bath, and the Wigner transform of an operator, $\hat{O}$, is defined as Imre (1967); Hillery et al. (1984)

[TABLE]

The Wigner transform of the product of two operators is Imre (1967); Hillery et al. (1984)

[TABLE]

where

[TABLE]

is the Poisson bracket operator and the arrows indicate the directions in which the differential operators act.

The MFT approximation to this correlation function amounts to replacing the time-dependent operator $(\mathcal{B}(t))^{W}=(e^{iHt/\hbar}\mathcal{B}e^{-iHt/\hbar})^{W}$ with $\mathcal{B}^{W}(X_{t},t)$, where the classical bath degrees of freedom evolve according to the Poisson bracket subject to an additional time-dependent mean force arising from the subsystem dynamics and the quantum variables evolve according to the time-dependent Schrodinger equation subject to the time-dependent contribution of the bath. This allows one to rewrite the dynamical matrices, $q_{1}^{(jk,\pm)}(t)$ and $q_{3}^{(j,\pm)}(t)$, corresponding to a general $N_{s}$ state system in the MFT approximation as,

[TABLE]

where

[TABLE]

Here we have explicitly included the dependence of operator $\hat{A}_{n}(t)$ on the time-dependent bath variables, $X_{t}$, to emphasize that the value of $C_{mn}^{(s)}(t,X_{t})$ is dependent on the realization of the bath.

For the Frenkel exciton model, where the bath is harmonic and the bath part of the system-bath coupling is linear in the bath coordinates and given by $\hat{\Gamma}_{j}\rightarrow\hat{\gamma}_{j}=\sum_{k}c_{j,k}\hat{R}_{j,k}$, the Wigner transform of the anticommutator and commutator of the bath canonical distribution and $\hat{\gamma}_{k}$ take simple forms,

[TABLE]

where, using Eq. (38), one obtains

[TABLE]

and

[TABLE]

Application of the Ehrenfest MFT approximation to the dynamics of the Frenkel exciton model leads to the evolution of the subsystem wavefunction subject to the Hamiltonian of the subsystem modulated by the time-dependent system-bath coupling,

[TABLE]

where the subsystem initial condition, $\rho_{s}(0)=\ket{\psi}\bra{\psi}$, is of pure form.

The bath, in turn, follows Hamilton’s equations,

[TABLE]

where the Ehrenfest version of the bath Hamiltonian takes the form,

[TABLE]

where $\overline{S}_{j}(t)=\mathrm{Tr}_{s}[\rho_{s}(t,X_{t})\hat{S}_{j}]$ is the mean potential arising from the subsystem dynamics.

The correlation functions given in Eqs. (40a)-(40d) are generated using MFT and the memory kernel is then generated according to Eqs. (6)-(8).

Appendix D Correlation functions with off-diagonal initial conditions using Ehrenfest mean field theory

In this section, we describe how the correlation functions in Eqs. (10) and (11) are evaluated using MFT for initial conditions that do not correspond to a single wavefunction.

If the subsystem initial condition is diagonal, i.e. $\hat{A}_{j}=\ket{\alpha}\bra{\alpha^{\prime}}$ where $\alpha=\alpha^{\prime}$, then MFT is straightforward to apply because one can use the subsystem initial condition corresponding to pure state $\ket{\alpha}$, i.e. $\ket{\psi(0)}=\ket{\alpha}$. On the other hand, if $A_{j}$ is an off-diagonal (if $\alpha\neq\alpha^{\prime}$), the procedure is less straightforward because $A_{j}$ does not correspond to an initial condition that arises from a single wavefunction. Several possible approaches and subtle issues related to this are discussed in detail in the Appendices of Ref. Montoya-Castillo and Reichman, 2016. A simple and computationally efficient approach for calculating such correlation functions involves an auxiliary wavefunction,

[TABLE]

where $\phi$ is a random number uniformly sampled from the interval $[0,2\pi)$. Then consider the expression

[TABLE]

where $q(t)$ is a general memory kernel correlation function which is partially Wigner transformed over the bath degrees of freedom, such as those defined in Eqs. (40a)-(40d), and we have explicitly used that $\alpha\neq\alpha^{\prime}$. Equation (49) thus provides a way to access two $q$ correlation functions using the same initial wavefunction $\ket{\bar{\psi}}$, in particular those for which $\hat{A}_{m}=\ket{\alpha}\bra{\alpha^{\prime}}$ and $\hat{A}_{m}=\ket{\alpha^{\prime}}\bra{\alpha}$. The explicit procedure for calculating such correlation functions is:

Initialize the wavefunction $\ket{\psi(0)}=\ket{\bar{\psi}}$, sampling a new random number $\phi$ from a uniform distribution on the interval $[0,2\pi)$. Sample the bath initial conditions from the appropriate bath distribution, $\rho_{b,W}$, as in the diagonal case. 2. 2.

At $t=0$, evaluate and store $e^{i\phi}$ and $e^{-i\phi}$. 3. 3.

Run each trajectory using MFT to calculate the desired observables at time $t$ in the same way as for a diagonal initial condition. Multiply the correlation function by $e^{i\phi}$ to obtain $\bar{q}_{+}(t)$ and by $e^{-i\phi}$ to obtain $\bar{q}_{-}(t)$. 4. 4.

Repeat steps 1-4, averaging $\bar{q}_{\pm}(t)$ over an ensemble of trajectories with different bath initial conditions and different initial values of $\phi$. As shown in Eq. (49), the trajectory-averaged $\bar{q}_{+}$ corresponds $q(t)$ with $A_{m}=\ket{\alpha}\bra{\alpha^{\prime}}$ and $\bar{q}_{-}$ corresponds to $q(t)$ with $A_{j}=\ket{\alpha^{\prime}}\bra{\alpha}$.

This approach of phase averaging over a superposition state in order to obtain the correlation function of interest is similar to that employed in the recently proposed CFBT method,Sato, Kelly, and Rubio (2018) which extends Ehrenfest mean field theory.

Appendix E Convergence analysis of the memory kernel

Here, we describe how to motivate our selective sampling algorithm from formal analysis of how the error in the memory kernel arising from lack of convergence translates into error in the GQME dynamics.

Our analysis starts from the GQME corresponding to the fully converged RDM, $\rho_{s}^{conv}(t)$, and the GQME using corresponding to a (possibly unconverged) RDM $\rho_{s}(t)$. The fully converged dynamics arise from a fully converged memory kernel, $\mathcal{K}^{\mathrm{conv}}(t)$, while the unconverged dynamics arise from an unconverged memory kernel, $\mathcal{K}(t)$. A simple metric for the error in $\rho_{s}(t)$ is its deviation from the converged RDM,

[TABLE]

Since this quantity can be positive or negative, it is convenient to consider its absolute value, $\left|\delta\rho_{s}(t)\right|$. This deviation is time dependent and is different for each element of the density. A scalar metric for the total error is the mean absolute deviation of the RDM from its converged value, with the mean taken over time and elements of the density,

[TABLE]

where $t_{max}$ is the total amount of time over which one wants to know the RDM dynamics. Because there is no error in the initial condition, (i.e., $\delta\rho_{s}(0)=0$), one can minimize the error at subsequent times by minimizing the error in the time-derivative of the RDM,

[TABLE]

Equation (5) can be combined with Eq. (52) in order to express the error in the time derivative, $\delta\dot{\rho}_{s}(t)$, in terms of $\delta\rho_{s}(t)$ and $\delta\mathcal{K}(t)=\mathcal{K}^{\mathrm{conv}}(t)-\mathcal{K}(t)$:

[TABLE]

where

[TABLE]

Of the three terms, only $E_{2}$ can be directly controlled by running more trajectories, since more trajectories lead to a smaller $\delta\mathcal{K}(t)$ as the correlation functions that give rise to $\mathcal{K}(t)$ become more converged. In contrast, $E_{1}$ and $E_{3}$ cannot be directly controlled because they depend on $\delta\rho_{s}(t)$, a quantity which arises from $\delta\mathcal{K}(t)$ and previous values of $\delta\rho_{s}(t)$ according to Eqs. (52)-(56). For our choice of selective sampling algorithm we make the ansatz that the mean absolute deviation in the RDM dynamics is proportional to the mean absolute value of $E_{2}$, i.e., $\langle\left|\delta\rho_{s}\right|\rangle\propto\langle|E_{2}|\rangle$. Figure 7a, which is generated using data from LHCII, shows strong correlation between the error in the GQME dynamics and $|E_{2}|$, suggests that this ansatz is valid.

In order to obtain an expression for the selective sampling weights, $w_{m}$, we consider how $\langle|E_{2}|\rangle$ is related to the error in specific elements of $\mathcal{K}(t)$,

[TABLE]

Using this expression to determine the error in the dynamics still requires prior knowledge of the converged kernel. In order to use this expression to obtain an expression for the selective sampling weights, we make the additional approximation that $\left|\delta K(t)\right|\approx\left|\delta K_{1}(t)\right|\approx\sqrt{\frac{\sigma^{2}_{K_{1}}(t)}{N}}$ so that:

[TABLE]

where $\sigma_{\mathcal{K}_{1}}$ is the standard deviation in $\mathcal{K}_{1}$:

[TABLE]

Here it is useful to re-write Eq. (58) as a sum of terms arising from different initial conditions at the memory kernel level:

[TABLE]

where $m$ indexes the initial conditions at the memory kernel level and $[\epsilon(t)]_{m}$ is the sum of all the terms in Eq. (58) arising from initial condition $m$:

[TABLE]

Equations (60) and (61) provide the basis for our choice of weights in our selective sampling algorithm. In particular, the weight $w_{m}$ associated with initial condition $m$ is proportional to the time average of $[\epsilon(t)]_{m}$

[TABLE]

This metric is an upper bound parts of the $|E_{2}|$ error term that arise from initial condition $m$, as shown in Eq. (58). As shown in Fig. 7b, $|E_{2}(t)|$ and $|\epsilon(t)|$ track the error in the GQME dynamics, $|\delta\rho_{s}(t)|$. Hence, a selective sampling algorithm that minimizes this term should also minimize the error in the GQME dynamics. For the models considered here, we note that $[\epsilon(t)]_{m}$ takes a somewhat simpler form because $[c^{(j,\pm)}]_{\alpha m}=0$ if $\alpha\neq m$ and hence

[TABLE]