# Linear and nonlinear spectroscopy from quantum master equations

###### Abstract

We investigate the accuracy of the second-order time-convolutionless (TCL2) quantum master equation for the calculation of linear and nonlinear spectroscopies of multichromophore systems. We show that, even for systems with non-adiabatic coupling, the TCL2 master equation predicts linear absorption spectra that are accurate over an extremely broad range of parameters and well beyond what would be expected based on the perturbative nature of the approach; non-equilibrium population dynamics calculated with TCL2 for identical parameters are significantly less accurate. For third-order (two-dimensional) spectroscopy, the importance of population dynamics and the violation of the so-called quantum regression theorem degrade the accuracy of TCL2 dynamics. To correct these failures, we combine the TCL2 approach with a classical ensemble sampling of slow microscopic bath degrees of freedom, leading to an efficient hybrid quantum-classical scheme that displays excellent accuracy over a wide range of parameters. In the spectroscopic setting, the success of such a hybrid scheme can be understood through its separate treatment of homogeneous and inhomogeneous broadening. Importantly, the presented approach has the computational scaling of TCL2, with the modest addition of an embarrassingly parallel prefactor associated with ensemble sampling. The presented approach can be understood as a generalized inhomogeneous cumulant expansion technique, capable of treating multilevel systems with non-adiabatic dynamics.

## I Introduction

Non-adiabatic energy transfer in molecular systems represents a problem of broad interest in chemistry, physics, and materials science. In the condensed phase, these processes commonly occur with many comparable energy scales, precluding simple perturbative treatments of the dynamics such as Golden-rule-type rate theories. This class of problems has motivated the important development of accurate numerical techniques capable of evolving the electronic reduced density matrix and offering insight into the population dynamics of multi-level dissipative quantum systems Tanimura and Kubo (1989); Mak and Chandler (1990); Makarov and Makri (1994); Beck et al. (2000); Thoss, Wang, and Miller (2001); Dunkel, Bonella, and Coker (2008); Chen, Cohen, and Reichman (2017). However, with few exceptions, elements of the reduced density matrix are basis-dependent and not directly observable. Instead, the principal experimental probes of energy transfer dynamics are ultrafast time-resolved spectroscopies, such as pump-probe transient absorption and coherent two-dimensional spectroscopy Mukamel (1995); Hybl et al. (1998); Mukamel (2000); Jonas (2003); Brixner et al. (2005); Engel et al. (2007). In this manuscript, we evaluate the accuracy of perturbative, but non-Markovian, quantum master equations for the calculation of linear and nonlinear spectroscopies. In particular, we focus on systems characteristic of protein-protected biological chromophores. This class of problems has been the topic of intense study in the quantum dynamics community Jang et al. (2008); Plenio and Huelga (2008); Mohseni et al. (2008); Huo and Coker (2010); Wu et al. (2010); Kolli, Nazir, and Olaya-Castro (2011); Tiwari, Peters, and Jonas (2013); Tempelaar, Jansen, and Knoester (2014), in part because the environmental fluctuations exhibit long correlation times that challenge conventional theories Ishizaki and Fleming (2009a, b).

The Hamiltonian for multichromophore systems can be generically written as the sum of an electronic (system) Hamiltonian, a nuclear (bath) Hamiltonian, and the interaction between the two, . In the present manuscript, we consider a Frenkel exciton model of coupled chromophores, with the system Hamiltonian

(1) |

where is a mean excitation energy for the excited-state manifold (equal to 10,000–30,000 cm for visible-light absorbing chromophores), is the deviation from this mean excitation energy for site , and is the electronic coupling between sites and . The operators and create and annihilate localized excitations on site and satisfy the commutation relation . The bath Hamiltonian is that of the nuclear degrees of freedom in the electronic ground-state. The system-bath interaction can be generically decomposed into the form where are bath operators and are system operators. For simplicity, we consider the case

(2) |

where is a collective bath operator whose fluctuations act to modulate the energy gap of molecular site . An otherwise generic second-order perturbation theory in the system-bath interaction requires only the equilibrium time correlation function of the nuclear degrees of freedom in the electronic ground state,

(3) |

where . For simplicity, we assume nuclear degrees of freedom belonging to different molecular sites are uncorrelated. This approach is commonly pursued to account for dephasing dynamics via atomistic simulations; in this approach, classical molecular dynamics are used to generate trajectories for the evaluation of the energy gap autocorrelation function in Eq. (3), leading to the second-order cumulant approximation to the lineshape Kubo (1969); Mukamel (1995),

(4) |

Strictly, the correlation function in Eq. (3) is a quantum time correlation function and a variety of approximate schemes exist to reconstruct a quantum time correlation function from its classical counterpart Egorov, Everitt, and Skinner (1999); Egorov, Rabani, and Berne (1999). We note that the above formalism neglects energy transfer (or non-adiabatic effects) associated with the intermolecular couplings . In this limit of “pure dephasing,” the second order cumulant approximation is quantum mechanically exact if the operators are linear in the coordinates of a harmonic bath Skinner and Hsu (1986); Reichman, Silbey, and Suárez (1996); this linear coupling model will be adopted below. In the more general case of energy transfer associated with population relaxation, spectroscopy calculations require a more general approach to quantum dynamics, which we discuss in the next section.

## Ii Spectroscopy, correlation functions, and the quantum regression theorem

In order to calculate spectroscopic observables, we augment our Hamiltonian with a light-matter interaction term to be treated via time-dependent perturbation theory Mukamel (1995),

(5) |

here, is a classical electric field and is the dipole operator

(6) |

where, via the Condon approximation, the dipole matrix element is independent of the bath degrees of freedom. For simplicity, here and henceforth we neglect the vectorial nature of the transition dipole matrix element.

Linear absorption spectra are calculated from the equilibrium time-correlation response function,

(7) |

where and the Liouville-space propagator is defined by . Two-dimensional spectra are calculated from the third-order response function,

(8) |

With exact quantum dynamics, the above expressions produce exact spectra that contain the effects of coherence dephasing as well as population relaxation, with explicit treatment of system-bath correlations spanning multiple light-matter interactions.

As a one-time observable, the linear-response function can be calculated exactly from the time evolution of a reduced density-like operator

(9) | ||||

(10) |

Therefore, any theory of reduced dynamics may be used to calculate the linear response function, such as a formally exact time-convolutionless quantum master equation Shibata, Takahashi, and Hashitsume (1977); Chaturvedi and Shibata (1979); Breuer and Petruccione (2007); Yoon, Deutch, and Freed (1975); Mukamel, Oppenheim, and Ross (1978),

(11) | ||||

(12) |

where is the time-ordering operator. Naturally, the lowest-order approximation to the time-dependent relaxation operator may be used; this is the second-order time convolutionless (TCL2) quantum master equation Breuer and Petruccione (2007). Unlike the second cumulant approach described in Sec. I, the TCL2 master equation (detailed below) provides a perturbative description of population relaxation and coherence dephasing on equal footing. Importantly, in the absence of population relaxation, the TCL2 approximation yields uncoupled dephasing dynamics of the (in this case off-diagonal) reduced operator that are identical to those of the second-order cumulant approximation, and thus exact for linear coupling to a harmonic bath. In other words, the dynamical resummation inherent in the TCL2 approximation is exact for this example of the pure-dephasing problem (due to underlying Gaussian statistics). We emphasize that this exactness is independent of the number of bath degrees of freedom, the energy scales of the bath, or the strength of the system-bath coupling.

However, in the presence of population relaxation with , the TCL2 approximation is no longer exact, and its range of validity is commonly understood to be restricted to the nearly-Markovian, weak-coupling limit under which it is typically derived. We will demonstrate that such statements are completely dependent on the observable and not determined exclusively by the energy scales in the Hamiltonian. In particular, we will show that for the same Hamiltonian, TCL2 predicts qualitatively incorrect non-equilibrium population dynamics but quantitatively accurate linear absorption lineshapes.

Unfortunately, the simplicity of the linear response function is not maintained for higher-order correlation functions. If we were to generalize Eqs. (9) and (10) to the third-order response function, we might be led to the tempting approximation

(13) |

which can be shown to be consistent with a quantum version of Onsager’s regression hypothesis known as the quantum regression theorem (QRT) Lax (1963); Gardiner and Zoller (2004). In general, the QRT is not exact Ford and O’Connell (1996); Swain (1999). The exact multi-time correlation function is related to the approximate one by a number of correction terms, which are non-vanishing even to lowest-order in the system-bath interaction Alonso and de Vega (2005, 2007); Goan, Chen, and Jian (2011). To summarize, in the context of the current manuscript, the rigorous calculation of nonlinear response functions requires more information than contained in quantum master equations.

The violation of the QRT is intimately linked to the degree of non-Markovianity present in the reduced dynamics. For purely Markovian reduced dynamics, the QRT is exact Swain (1999); Gisin (1993). For weak system-bath coupling leading to nearly Markovian reduced dynamics, the corrections to the QRT are small. In this manuscript, we will study the accuracy of TCL2 dynamics, within the approximation implied by the QRT, Eq. (13), for the simulation of two-dimensional spectroscopy. Unsurprisingly, we find that the results deteriorate with increasing non-Markovianity due to stronger coupling or slower bath degrees of freedom.

In order to maintain the simplicity of quantum master equations while seeking broad applicability to nonlinear spectroscopy, we will evaluate the use of the “frozen modes” approach, recently introduced by one of us and co-workers Montoya-Castillo, Berkelbach, and Reichman (2015). The method will be described in more detail below, but the idea is to simulate generically non-Markovian dynamics by an average over many independent nearly-Markovian trajectories. In each trajectory, the low-frequency bath degrees of freedom are dynamically arrested: they are removed from the master equation’s relaxation operator and treated as a source of static disorder. This frozen-mode approximation is obviously best for low-frequency (slow) degrees of freedom, which are precisely those that contribute to non-Markovian behavior and the inaccuracy of perturbative quantum master equations. Compared to the original problem, each individual trajectory in the frozen modes approach exhibits less non-Markovian behavior and weaker system-bath coupling. Because these are the same effects responsible for the violation of the QRT in nonlinear response calculations, we will demonstrate that the frozen-mode variant of TCL2 dynamics leads to accurate two-dimensional spectra, even in highly non-Markovian regimes.

## Iii Model and Methods

For the remainder of this work, we will adopt the common system-bath model that assumes a harmonic nuclear bath linearly coupled to the system’s excitation number operator, i.e. Eqs. (1) and (2) with

(14) | ||||

(15) |

where and are the mass-weighted momentum and position of mode belonging to the nuclear degrees of freedom of site . For such a simplified form of the system-bath interaction, all properties are determined solely by the equilibrium autocorrelation function in Eq. (3),

(16) |

where we have introduced the spectral density

(17) |

Note that the spectral density can be obtained from the real-part of the cosine-transform of the bath correlation function

(18) |

which provides a way to extract a model spectral density from atomistic simulations of the energy gap autocorrelation function, as done recently for light-harvesting complexes Valleau, Eisfeld, and Aspuru-Guzik (2012); Lee and Coker (2016). In the following, we assume all sites have identical spectral densities and employ an Ohmic spectral density with a high-frequency Lorentzian cutoff

(19) |

which follows from the assumption of an exponentially-decaying autocorrelation function (in the high-temperature limit), . Therefore, the bath relaxation time is given by and the magnitude of fluctuations is related to the reorganization energy . We emphasize that most of our conclusions are not restricted to any particular form of , including those with arbitrary structure (as may be obtained from simulation). Specifically, we expect that the accuracy of our results is only weakly dependent on the form of the spectral density, and can be safely applied to any system-bath Hamiltonian exhibiting only linear coupling to a harmonic bath. However, while the method is entirely applicable to more general system-bath Hamiltonians, the accuracy is harder to assess. For example, we note that even for the pure-dephasing problem, the second cumulant is no longer exact for systems that are quadratically-coupled to harmonic baths Skinner and Hsu (1986); Reichman, Silbey, and Suárez (1996).

### iii.1 Quantum master equations and TCL2

Introduced above, the TCL2 quantum master equation evolves reduced-density-like system operators according to the time-local equation of motion in Eq. (20). This equation of motion assumes a factorized initial condition of the total density-like operator , and we henceforth assume that the initial bath density operator is at equilibrium, . We note that for Hamiltonians with a minimum excited-state gap that is much larger than thermal energy, the total equilibrium density operator is factorized to an excellent approximation, , justifying factorized initial conditions for equilibrium correlation functions. In the basis of eigenstates of the system Hamiltonian, the TCL2 relaxation superoperator is a tensor with elements

(20) |

with

(21) |

and

(22) |

In the above, and is given in Eq. (16). The relevant dimensionless parameter that controls the accuracy of weak-coupling master equations is roughly given by ; in particular, very slow bath relaxation dynamics are a challenge as they violate the quasi-Markovian approximations inherent in such master equations. As such, we will also consider a quantum-classical hybrid scheme that microscopically treats slow nuclear degrees of freedom as static variables sampled from their respective canonical distributions, discussed next.

### iii.2 TCL2 master equation with frozen modes

The spectral density defined in Eq. (17) can formally be separated into fast and slow components under the condition that . In the ideal partitioning, the slow part is comprised of quasi-adiabatic modes that evolve much more slowly than the system and can safely be treated classically, while the fast part contains modes which evolve much more quickly than the system and therefore can be treated with nearly-Markovian weak-coupling theories Thoss, Wang, and Miller (2001); Berkelbach, Reichman, and Markland (2012); Berkelbach, Markland, and Reichman (2012). The partitioning between slow and fast parts of can be done with a switching function,

(23a) | ||||

(23b) |

where is a function which goes from a value of 1 at to a value of 0 at a specified splitting frequency . Following previous work Berkelbach, Reichman, and Markland (2012); Berkelbach, Markland, and Reichman (2012); Montoya-Castillo, Berkelbach, and Reichman (2015) we use the smooth switching function

(24) |

which avoids long-time oscillatory features in the bath correlation function. The free parameter determines the set of modes to be treated classically. Two example partitionings of an Ohmic-Lorentz spectral density are shown in Fig. 1.

While the slow modes could be treated with several different semiclassical techniques, including mean-field Ehrenfest dynamics Thoss, Wang, and Miller (2001); Berkelbach, Reichman, and Markland (2012); Berkelbach, Markland, and Reichman (2012), here we make the simplest approximation and treat the modes described by as completely frozen. In this case, positions of the slow degrees of freedom are sampled from the Boltzmann distribution of a harmonic oscillator and used to alter the energy levels of the system Hamiltonian,

(25) |

where is a renormalized coupling constant due to the partitioning enforced by the switching function . This new system Hamiltonian is used with the fast part of the spectral density to perform a single realization of TCL2 dynamics. This process is repeated and averaged. Further details and discussion of the method can be found in Ref. Montoya-Castillo, Berkelbach, and Reichman, 2015.

In the spectroscopic context, the frozen modes are most physically interpreted as a source of inhomogeneous broadening. If the nuclear degrees of freedom are extremely slow, so as to appear effectively frozen during the decay of the dipole autocorrelation function, then this is the correct result. We emphasize that this inhomogeneous broadening is entirely microscopic because it originates from degrees of freedom in the total Hamiltonian; it is not an artificial, unidentified source of inhomogeneous broadening. When applied to linear absorption spectra of systems without non-adiabatic effects, this approach is equivalent to the inhomogeneous cumulant expansion described in Refs. Fried and Mukamel, 1993; Mukamel, 1995, but differs for higher-order spectroscopies.

The combination of TCL2 dynamics with frozen-mode sampling (TCL2-FM) leads to many trajectories that are each more Markovian than the original problem. This property mitigates corrections to the QRT, which are not considered explicitly in this work. Furthermore, the TCL2-FM approach prevents the application of perturbation theory beyond its regime of validity (through the use of a faster, more weakly-coupled bath in the quantum master equation). As we show below, these two effects collectively produce semiquantitative accuracy in the prediction of nonlinear spectroscopy.

## Iv Results

### iv.1 Model parameters and methodological details

For simplicity, we present results for a system of two chromophores that has four accessible electronic states: the ground state, two states where each chromophore is separately excited, and one state where both chromophores are simultaneously excited. The absence of a quartic exciton-exciton interaction in Eq. (1) implies that the energy of the doubly-excited state is simply the sum of the excitation energies of the individual chromophores. In all results, we use the electronic parameters cm, cm and cm, along with the bath parameters fs and K. The reorganization energy will be varied. The bath frequency and temperature lead to energy scales cm and cm, i.e. all energy scales are comparable, with the bath frequency being the smallest.

Our TCL2 dynamics are generated with a fourth-order Runge-Kutta integrator. For TCL-FM calculations, the total spectral density was discretized into 300 modes, and those with were discarded. We performed realizations of frozen-mode sampling to converge the dynamical observables. For all frozen-mode calculations, we use a splitting frequency that characterizes the timescale of isolated electronic dynamics, . Other choices that differ by factors of order one give similar qualitative results, and occasionally better quantitative results, but we find that the above form is reasonable over a broad range of parameters and observables.

All approximate results will be compared to numerically exact results obtained with the hierarchical equations of motion (HEOM) Tanimura and Kubo (1989); Ishizaki and Tanimura (2005); Xu and Yan (2007); for the bath parameters used in our results, the HEOM calculations required zero Matsubara frequencies (i.e. the high-temperature approximation with ), but as many as levels in the hierarchy. For spectroscopy calculations, we use transition dipole matrix elements satisfying and spectra will be presented with arbitrary units. All calculations were performed with our open-source quantum dynamics package pyrho pyr ().

### iv.2 Linear spectroscopy and population dynamics

In the frequency domain, we present the imaginary (absorptive) part of the linear-response susceptibility,

(26) |

Again we emphasize that for the model Hamiltonian adopted here, the TCL2 master equation is identical to the second-cumulant approximation and exactly solves the pure-dephasing lineshape problem. In the case of excited-state electronic coupling , the TCL2 generalizes the second cumulant approximation. While it is not exact, the results are remarkably accurate, as shown in the right-hand column of Fig. 2, even for significant coupling to a slow bath. The TCL2-FM approach yields very minor quantitative improvements, which confirms that those local-frequency modes treated as frozen are properly interpreted as giving rise to inhomogeneous broadening.

Unlike the non-Markovian TCL2-based approaches, the Markovian Redfield theory predicts incorrect spectra, with an accuracy that deteriorates for increasing reorganization energy or increasing bath relaxation times. This can be understood simply from the Markovian limit of the relevant pure-dephasing term in the TCL2 tensor: where is the Bose-Einstein distribution function. For the Ohmic-Lorentz spectral density employed here, one finds . As shown in Fig. 2, this Markovian theory predicts linewidths that are much too large.

This large discrepancy between two weak-coupling theories (Markovian Redfield theory and non-Markovian TCL2-based theories) is quite surprising given the perturbative nature of both approaches. In the left-hand column panels of Fig. 2, we show the population relaxation dynamics generated by the same methods, using the initial condition . Clearly, Redfield theory and conventional TCL2 fail in a similar manner as the perturbation becomes large. The TCL2-FM approach now yields results that are quite different than vanilla TCL2, and the former predicts population dynamics that are in good agreement with the numerically exact HEOM results. These observations demonstrate one of our main conclusions: the accuracy of a given dynamical technique depends on the observable. More specifically, TCL2-based approaches are well-suited to the evolution of electronic coherences; we believe that this is because TCL2 reduces to the exact second-cumulant solution for pure-dephasing problems. Population dynamics are less accurate.

### iv.3 Third-order nonlinear spectroscopy

Compared to linear response considered above, third-order nonlinear spectroscopy is a more challenging test, due to the simultaneous importance of coherence and population dynamics. As discussed in Sec. II, nonlinear spectroscopy presents an additional complication for non-Markovian quantum master equations in the form of violations to the QRT.

We present two-dimensional electronic spectra, obtained by Fourier transforming over the pump () and probe () time delays, resolved by the population waiting time . In particular, we simulate the photon echo spectrum, generated by six terms in the rotating wave approximation associated with rephasing (rp) and non-rephasing (nr) pathways,

(27) |

In the above,

(28a) | ||||

(28b) |

and

(29a) | ||||

(29b) | ||||

(29c) | ||||

(29d) |

For a moderate reorganization energy of cm, our results are shown in Fig. 3 for three values of the waiting time . Although TCL2 dynamics are qualitatively correct at , this accuracy degrades with increasing waiting times. In particular, the spectra become artificially broadened along the axis – as discussed in Ref. Ishizaki and Tanimura, 2008 – and completely fail to describe the four-peak structure and regions of negativity. In contrast, the TCL2-FM two-dimensional spectra are in qualitative agreement with the exact results, reproducing the appropriate peak shapes, spectral features, and timescales. The TCL2-FM results show an unphysical diagonal elongation at long waiting times due to the completely frozen modes that are unable to relax. We have found that using a smaller splitting frequency (i.e. freezing fewer modes) does improve the long-time results at the expense of short-time results. This observations suggests an interesting time-dependent frozen-modes scheme – where modes are successively unfrozen at times longer than their characteristic relaxation time – though we do not pursue this approach here. We conclude that the computationally efficient TCL2-FM approach yields accurate coherence and population dynamics, while mitigating correction terms due to violation of the QRT, leading to qualitatively correct two-dimensional spectra. Although we only show one parameter set for simplicity, we have tested other parameter sets and find that the TCL2-FM accuracy is robust over a wide range of Hamiltonians.

For a simpler representation of the third-order response, we also present the time-resolved pump-probe spectra, obtained by setting (or equivalently, integrating over ). We note that the correlation functions contributing to pump-probe spectroscopy are two-point nonequilibrium correlation functions, which should depend sensitively on the accuracy of excited-state population dynamics and exhibit nontrivial corrections to the QRT Alonso and de Vega (2005, 2007). For the same parameters as in Fig. 3, we present the pump-probe spectra in Fig. 4. Consistent with our findings for the full two-dimensional spectra, we see that TCL2 dynamics are accurate at zero waiting time, but completely wrong for nonzero waiting times; the spectra are much too broad and lacking multi-peaked structure. TCL2-FM dynamics cures these deficiencies and predicts accurate pump-probe spectra.

## V Conclusions

We have investigated the accuracy of the second-order time-convolutionless (TCL2) quantum master equation for the prediction of linear and third-order time-resolved spectroscopy. We have argued that TCL2 dynamics generalizes the second-cumulant approximation to problems exhibiting excited-state electronic coupling leading to non-adiabatic dynamics and population relaxation. For many contemporary problems concerning the dynamics of multichromophore systems, the bath timescales are too slow to permit classic Markovian theories of lineshapes. For such systems, the TCL2 approach generates nearly-exact linear-response spectra, even though simulations of population dynamics – using the same Hamiltonian parameters – are grossly in error.

In contrast, the TCL2 approximation breaks down away from its perturbative regime when used to simulate nonlinear spectroscopies. We have attributed this failure to two distinct effects: the greater importance of population dynamics and the violation of the quantum regression theorem (QRT). By partitioning microscopic bath degrees of freedom into fast and slow sets, and treating the latter as frozen variables sampled from a thermodynamic ensemble, we have argued that the TCL2-FM approach alleviates both of these problems. Physically, the TCL2-FM formalism treats modes responsible for highly non-Markovian dynamics as a source of inhomogeneous broadening; the resulting Hamiltonian exhibits dynamics which are consequently more Markovian and thus well-described by weak-coupling, TCL2-type quantum master equations. It will be interesting to implement the second-order corrections to the QRT, in order to assess their relative importance in the prediction of nonlinear spectra. Work along these lines is currently in progress.

Importantly, unlike many reduced quantum dynamics techniques, the formalism of perturbative quantum master equations is not tied to specific forms of the bath Hamiltonian or the system-bath interaction, as long as multi-point equilibrium correlation functions of isolated bath operators can be computed. Nonetheless, the accuracy of such approaches remains to be assessed for more generic Hamiltonians, but the absence of numerically exact results makes this challenging.

Our findings have implications for other quantum dynamics techniques. In particular, we recall that the hierarchical equations of motion (HEOM) reduce to time-convolutionless and time-convolution quantum master equations when the hierarchy is truncated at low order Xu et al. (2005, 2017). These low-order HEOM calculations will exhibit the same features described here. For linear spectroscopy, this observation argues strongly for the use of the time-convolutionless closure Xu et al. (2005) of the hierarchy, in support of the numerical observations made in Refs. Chen et al., 2009, 2010. For nonlinear spectroscopy, low-order approximations will also violate the QRT, and higher levels in the hierarchy will be needed in order to eliminate these (neglected) corrections.

With regards to computational cost, the TCL2-FM approach is extremely attractive. In its conventional form, presented here, it scales only as , where is a parallel prefactor associated with frozen-mode ensemble sampling and is size of the system Hilbert space; this scaling arises from the action of the relaxation tensor on the reduced density matrix. Because the frozen-modes approach leads to more Markovian dynamics, it is tempting to explore a stochastic unraveling of the reduced density matrix Breuer and Petruccione (2007); Zoller, Marte, and Walls (1987); Gisin and Percival (1992); Carmichael (1993); Gardiner and Zoller (2004). This latter approach replaces a single reduced density matrix simulation by an average over wavefunction dynamics with stochastic relaxation events. This procedure would have a scaling of , where is a parallel prefactor associated with the stochastic dynamics. We thus anticipate an accurate method which only requires many parallel simulations of system-wavefunction dynamics, each scaling like for generically non-sparse system Hamiltonians; it is hard to imagine a more attractive computational cost with comparable qualitative accuracy over an extremely broad range of Hamiltonian parameters.

## Acknowledgments

We thank Dr. Roel Tempelaar and Prof. Jim Skinner for helpful discussions. This research was supported by start-up funds from the University of Chicago. Calculations were performed with resources provided by the University of Chicago Research Computing Center (RCC) and we thank Dr. Jonathan Skone of the RCC for assistance.

## References

- Tanimura and Kubo (1989) Y. Tanimura and R. K. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989).
- Mak and Chandler (1990) C. H. Mak and D. Chandler, Phys. Rev. A 41, 5709 (1990).
- Makarov and Makri (1994) D. E. Makarov and N. Makri, Chem. Phys. Lett. 221, 482 (1994).
- Beck et al. (2000) M. Beck, A. Jackle, G. A. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000).
- Thoss, Wang, and Miller (2001) M. Thoss, H. Wang, and W. H. Miller, J. Chem. Phys. 115, 2991 (2001).
- Dunkel, Bonella, and Coker (2008) E. R. Dunkel, S. Bonella, and D. F. Coker, J. Chem. Phys. 129, 114106 (2008).
- Chen, Cohen, and Reichman (2017) H.-T. Chen, G. Cohen, and D. R. Reichman, J. Chem. Phys. 146, 054105 (2017).
- Mukamel (1995) S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford University Press, 1995).
- Hybl et al. (1998) J. D. Hybl, A. W. Albrecht, S. M. Gallagher Faeder, and D. M. Jonas, Chem. Phys. Lett. 297, 307 (1998).
- Mukamel (2000) S. Mukamel, Annu. Rev. Phys. Chem. 51, 691 (2000).
- Jonas (2003) D. M. Jonas, Annu. Rev. Phys. Chem. 54, 425 (2003).
- Brixner et al. (2005) T. Brixner, J. Stenger, H. M. Vaswani, M. Cho, R. E. Blankenship, and G. R. Fleming, Nature 434, 625 (2005).
- Engel et al. (2007) G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. T. Mancal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming, Nature 446, 782 (2007).
- Jang et al. (2008) S. Jang, Y. C. Cheng, D. R. Reichman, and J. D. Eaves, J. Chem. Phys. 129, 101104 (2008).
- Plenio and Huelga (2008) M. B. Plenio and S. F. Huelga, New J. Phys. 10, 113019 (2008).
- Mohseni et al. (2008) M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, J. Chem. Phys. 129, 174106 (2008).
- Huo and Coker (2010) P. Huo and D. F. Coker, J. Chem. Phys. 133, 184108 (2010).
- Wu et al. (2010) J. Wu, F. Liu, Y. Shen, J. Cao, and R. J. Silbey, New J. Phys. 12, 105012 (2010).
- Kolli, Nazir, and Olaya-Castro (2011) A. Kolli, A. Nazir, and A. Olaya-Castro, J. Chem. Phys. 135, 154112 (2011).
- Tiwari, Peters, and Jonas (2013) V. Tiwari, W. K. Peters, and D. M. Jonas, Proc. Natl. Acad. Sci. 110, 1203 (2013).
- Tempelaar, Jansen, and Knoester (2014) R. Tempelaar, T. L. C. Jansen, and J. Knoester, J. Phys. Chem. B 118, 12865 (2014).
- Ishizaki and Fleming (2009a) A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234110 (2009a).
- Ishizaki and Fleming (2009b) A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009b).
- Kubo (1969) R. Kubo, Adv. Chem. Phys. 15, 101 (1969).
- Egorov, Everitt, and Skinner (1999) S. A. Egorov, K. F. Everitt, and J. L. Skinner, J. Phys. Chem. A 103, 9494 (1999).
- Egorov, Rabani, and Berne (1999) S. A. Egorov, E. Rabani, and B. J. Berne, J. Phys. Chem. B 103, 10978 (1999).
- Skinner and Hsu (1986) J. L. Skinner and D. Hsu, J. Phys. Chem. 90, 4931 (1986).
- Reichman, Silbey, and Suárez (1996) D. Reichman, R. J. Silbey, and A. Suárez, J. Chem. Phys. 105, 10500 (1996).
- Shibata, Takahashi, and Hashitsume (1977) F. Shibata, Y. Takahashi, and N. Hashitsume, J. Stat. Phys. 17, 171 (1977).
- Chaturvedi and Shibata (1979) S. Chaturvedi and F. Shibata, Z. Phys. B Condens. Matter 35, 297 (1979).
- Breuer and Petruccione (2007) H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, 2007).
- Yoon, Deutch, and Freed (1975) B. Yoon, J. M. Deutch, and J. H. Freed, J. Chem. Phys. 62, 4687 (1975).
- Mukamel, Oppenheim, and Ross (1978) S. Mukamel, I. Oppenheim, and J. Ross, Phys. Rev. A 17, 1988 (1978).
- Lax (1963) M. Lax, Phys. Rev. 129, 2342 (1963).
- Gardiner and Zoller (2004) P. Gardiner and P. Zoller, Quantum Noise (Springer-Verlag, Berlin, 2004).
- Ford and O’Connell (1996) G. W. Ford and R. O’Connell, Phys. Rev. Lett. 77, 798 (1996).
- Swain (1999) S. Swain, J. Phys. A. Math. Gen. 14, 2577 (1999).
- Alonso and de Vega (2005) D. Alonso and I. de Vega, Phys. Rev. Lett. 94, 1 (2005).
- Alonso and de Vega (2007) D. Alonso and I. de Vega, Phys. Rev. A - At. Mol. Opt. Phys. 75, 052108 (2007).
- Goan, Chen, and Jian (2011) H. S. Goan, P. W. Chen, and C. C. Jian, J. Chem. Phys. 134 (2011), 10.1063/1.3570581.
- Gisin (1993) N. Gisin, J. Mod. Opt. 40, 2313 (1993).
- Montoya-Castillo, Berkelbach, and Reichman (2015) A. Montoya-Castillo, T. C. Berkelbach, and D. R. Reichman, J. Chem. Phys. 143, 194198 (2015).
- Valleau, Eisfeld, and Aspuru-Guzik (2012) S. Valleau, A. Eisfeld, and A. Aspuru-Guzik, J. Chem. Phys. 137, 224103 (2012).
- Lee and Coker (2016) M. K. Lee and D. F. Coker, J. Phys. Chem. Lett. 7, 3171 (2016).
- Berkelbach, Reichman, and Markland (2012) T. C. Berkelbach, D. R. Reichman, and T. E. Markland, J. Chem. Phys. 136, 034113 (2012).
- Berkelbach, Markland, and Reichman (2012) T. C. Berkelbach, T. E. Markland, and D. R. Reichman, J. Chem. Phys. 136, 084104 (2012).
- Fried and Mukamel (1993) L. E. Fried and S. Mukamel, Adv. Chem. Phys. 84, 435 (1993).
- Ishizaki and Tanimura (2005) A. Ishizaki and Y. Tanimura, J. Phys. Soc. Japan 74, 3131 (2005).
- Xu and Yan (2007) R.-X. Xu and Y. Yan, Phys. Rev. E 75, 031107 (2007).
- (50) “pyrho: a python package for reduced density matrix techniques,” https://github.com/berkelbach-group/pyrho.
- Ishizaki and Tanimura (2008) A. Ishizaki and Y. Tanimura, Chem. Phys. 347, 185 (2008).
- Xu et al. (2005) R. X. Xu, P. Cui, X. Q. Li, Y. Mo, and Y. Yan, J. Chem. Phys. 122, 041103 (2005).
- Xu et al. (2017) M. Xu, L. Song, K. Song, and Q. Shi, J. Chem. Phys. 146, 064102 (2017).
- Chen et al. (2009) L. Chen, R. Zheng, Q. Shi, and Y. Yan, J. Chem. Phys. 131, 094502 (2009).
- Chen et al. (2010) L. Chen, R. Zheng, Q. Shi, and Y. Yan, J. Chem. Phys. 132, 024505 (2010).
- Zoller, Marte, and Walls (1987) P. Zoller, M. Marte, and D. Walls, Phys. Rev. A 35, 198 (1987).
- Gisin and Percival (1992) N. Gisin and I. C. Percival, J. Phys. A. Math. Gen. 25, 5677 (1992).
- Carmichael (1993) H. J. Carmichael, Phys. Rev. Lett. 70, 2273 (1993).