We present a methodology for simulating multidimensional electronic spectra of molecular aggregates with coupling of electronic excitation to a structured environment using the stochastic non-Markovian quantum state diffusion (NMQSD) method in combination with perturbation theory for the response functions. A crucial aspect of our approach is that we propagate the NMQSD equation in a doubled system Hilbert space but with the same noise. We demonstrate that our approach shows fast convergence with respect to the number of stochastic trajectories, providing a promising technique for numerical calculation of two-dimensional electronic spectra of large molecular aggregates.

## I. INTRODUCTION

Modern time-resolved nonlinear optical spectroscopic techniques have expanded our understanding of the photophysics of molecular assemblies.^{1–5} Two-dimensional (2D) electronic spectroscopy, the material response after interacting with three femtosecond laser pulses, is a particularly powerful probe of molecular excitons: 2D spectra provide information about exciton–exciton interactions, dephasing, and relaxation processes.^{6–12} Nevertheless, spectral congestion—even at low temperatures—makes theoretical simulations indispensable for deciphering the dynamics encoded in the 2D spectra.

The key quantity in the simulation of 2D spectra is the third-order optically induced polarization, which is related to third-order nonlinear response functions.^{2} A common theoretical framework for simulating the third-order polarization of molecular aggregates is based on open-quantum system approaches, which propagate the reduced density matrix of the electronic system along different Liouville pathways to obtain nonlinear response functions.^{1,2} Of these approaches, those based on the Redfield or the modified Redfield equations are among the most popular;^{13,14} their applicability, however, is restricted to the case of weak system–bath couplings and the Markovian approximation for the bath. Alternatively, the hierarchy equation of motion (HEOM)^{15,16} and the quasiadiabatic path integral (QUAPI)^{17,18} provide numerically exact descriptions of the non-perturbative and non-Markovian dynamics. However, while both methods have been widely used to simulate exciton dynamics and 2D spectra of molecular assemblies,^{19–27} they become numerically expensive for strong system–bath couplings, low temperatures, and large numbers of pigments.

An alternative to density matrix based methods is the *non-Markovian quantum state diffusion* (NMQSD) formalism, in which stochastic wavefunctions are propagated in the system Hilbert space and the density matrix is obtained from an average of these wavefunctions.^{28,29} Over the years, several approaches have been developed to efficiently solve the NMQSD equation numerically,^{30–34} so that simulations of excitation transport in large molecular aggregates containing thousands of pigments are now tractable. In these propagation schemes, importance sampling via the nonlinear NMQSD equations are essential for efficient convergence with respect to the number of trajectories. Within NMQSD, it is possible to obtain 2D spectra directly by including the femtosecond-pulses explicitly in the time-evolution and extracting the desired signal via phase-cycling, but the convergence with respect to the number of trajectories is slow compared to the calculation of expectation values.^{35}

To overcome this problem, we develop in the present work an NMQSD equation in which the response functions are obtained directly from a perturbative expansion with respect to interactions with the laser field. The crucial point of the scheme is to use the NMQSD formalism to propagate the combined *ket* and *bra* states of the density matrix in a doubled electronic Hilbert space, but having the same noise. Importance sampling via the nonlinear NMQSD equation introduces a coupling between the propagation of the *bra* and *ket* contributions. We refer to NMQSD propagation in the doubled electronic Hilbert space as dyadic NMQSD, consistent with our previous treatment of linear absorption.^{36} Here, we solve the general dyadic NMQSD using a numerically efficient representation known as the hierarchy of pure states (HOPS). Dyadic HOPS exhibits fast convergence with respect to the number of stochastic trajectories and treats singly and doubly excited excitonic states in a unified manner, which is essential to account for ground state bleach (GSB), stimulated emission (SE), and excited state absorption (ESA) contributions to 2D spectra of molecular aggregates.

This paper is organized as follows: In Sec. II, we introduce the details of the molecular system, its interaction with electromagnetic pulses, and the general form of the response functions. In Sec. III, we develop our method to calculate the response function using the NMQSD approach. We particularly emphasize the ability to use the nonlinear NMQSD equation that ensures suitable convergence with respect to trajectories. In Sec. IV, we perform numerical calculation and demonstrate that with only 1000 trajectories, spectra are already well-converged and discuss convergence trends in detail. Finally, we conclude in Sec. V with a summary and a brief outlook. In Appendix A, we connect the first order (linear response) of the present formalism to our previous calculations.^{36} In the supplementary material, we provide further details on the convergence with respect to the number of trajectories and demonstrate the convergence problems associated with the linear NMQSD equation that are resolved using our method.

## II. MOLECULAR SYSTEM, INTERACTION WITH LASER PULSES, AND THE QUANTITIES OF INTEREST

### A. Hamiltonian

We consider a molecular aggregate composed of *N* interacting molecules, where each molecule is described by two electronic levels: the electronic ground state |*g*_{n}⟩ and the electronic excited state |*e*_{n}⟩, *n* = 1, …, *N*. The electronic Hamiltonian $H\u0302ex$ can then be written as

where *ϵ*_{n} is the energy required to excite the *n*th molecule, *V*_{nm} is the electronic coupling between excited molecules *n* and *m*, and $\sigma \u0302n=gnen$. In the case of 2D spectroscopy, we need the common ground state $(|g\u3009=\u220fn=1N|gn\u3009)$, singly excited states $(|n\u3009=|en\u3009\u220fm\u2260n|gm\u3009=\sigma \u0302n\u2020g)$, and doubly excited states ($|nm\u3009=|en\u3009|em\u3009\u220fk\u2260n,m|gk\u3009=\sigma \u0302n\u2020\sigma \u0302m\u2020g$, *n* < *m*) of the molecular aggregate. As is commonly done for molecular aggregates, the Hamiltonian given by Eq. (1) does not couple states with different number of electronic excitations.

For each molecule, there are additional interactions with internal and external nuclear degrees of freedom. In many cases of interest, these interactions can be modeled by (infinite) sets of bosonic modes that couple linearly to the excitonic states. We denote these modes as *environment* or *bath*. In this work, we assume that each molecule has its own set of bath modes so that the Hamiltonian of the bath can be written as

Here, $b\u0302n\lambda $$(b\u0302n\lambda \u2020)$ is the annihilation (creation) operator of *λ*th bath mode of molecule *n* with frequency *ω*_{nλ}. The bath modes couple locally to their respective molecule. The interaction Hamiltonian is then written as

where the coupling operator $L\u0302n$ acts in the system Hilbert space and is given by

and *κ*_{nλ} is the exciton–bath coupling strength of the mode *λ* for molecule *n*. We write the complete matter Hamiltonian as

Initially, there are no electronic excitations in the system and the environment is in thermal equilibrium. Since $gL\u0302ng=0$, there is no initial interaction between electronic and environmental degrees of freedom and we can write the initial state as

It is convenient to introduce the bath spectral density of molecule *n*, *J*_{n}(*ω*) = ∑_{λ}|*κ*_{nλ}|^{2}*δ*(*ω* − *ω*_{nλ}), which is related to the bath-correlation function *α*_{n}(*t*) by

with the inverse temperature *β*.

### B. Interaction with laser field

In a 2D spectroscopy experiment, the system interacts with three laser pulses at controlled inter-pulse delay times. The field–matter interaction Hamiltonian $H\u0302F(t)$ is defined as^{2}

where $\mu \u0302$ is the total transition dipole operator given by

and *μ*_{n} the transition dipole moment of molecule *n*, taken to be real. Note that by writing the dipole operator in terms of local excitation creation and annihilation operators, it connects states that differ by one in the number of electronic excitations and we include an arbitrary numbers of excitations in the Hamiltonian.

The electric field is given by

with **e**_{a}, **k**_{a}, *ω*_{a}, *E*_{a}(*t*), and *t*_{a} denoting the polarization unit vector, the wave vector, the carrier frequency, the envelope, and the central time of the *a*th pulse, respectively. In general, there can also exist a different number of pulses.

We can write the total Hamiltonian $H\u0302tot$ as

with $H\u0302$ given in Eq. (5).

### C. Response functions

In this section, we present a general notation for nonlinear optical response functions that provide a clear connection to the NMQSD formalism. Here, we use a general notation that is appropriate for arbitrary orders of perturbation theory. In Sec. IV, we focus on 2D spectroscopy.

#### 1. Perturbation theory for the full density matrix

In the following, we denote the time-evolution operator of the system without the electromagnetic field as

where $H\u0302$ is given in Eq. (5). We also introduce the abbreviation

for the interaction with field at time *t*_{j}.

We are interested in correlation functions (that we loosely call response functions) that depend on the order of interactions with the electric field. We use a condensed notation to track the generic correlation functions of the form

where the expectation value of an observable $F\u0302$, in our case the polarization ** μ**, is calculated with respect to a density matrix that is obtained in

*M*th order of perturbation theory. In Eq. (14), the parameters in the rectangular brackets track the order of interactions influencing the

*ket*(row 2) vs

*bra*(row 3) time-evolution associated with a particular time correlation function. The first row contains the intervals

*τ*

_{i}=

*t*

_{i+1}−

*t*

_{i}between interaction at time

*t*

_{i}and the following interaction. The last interval

*τ*

_{M}contains the final time

*t*of the evolution, i.e.,

*τ*

_{M}=

*t*−

*t*

_{M}. The sequence of operators in the second and third rows identify the operators acting on the

*bra*and

*ket*side (respectively) at each interaction time and thereby define a specific response function. The

*M*th order density matrix is recursively defined by

and

In Eq. (15), the operators $vjK$ act always on the *ket* side and the operators $vjB$ always on the *bra* side of $\rho \u0302$. An important constraint is that for each pair $vjK$, $vjB$, with the same index *j*, one of the corresponding operators is the unit operator (which we denote by $1$). Note that for better readability, we will sometimes omit the “operator hats” on the operators $vjK/B$, as it has already been done above. Explicitly, we have

In the following, we introduce a shorthand notation where we abridge

by omitting all arguments when it is clear which correlation function is being considered or when we consider generic correlation functions.

## III. CALCULATION OF THE NONLINEAR RESPONSE FUNCTION USING NMQSD

### A. The NMQSD formalism

For the open-quantum system model as given in Sec. II A, the expectation value of any system operator $F\u0302$ can be obtained as^{28,29}

where $M[\cdots ]$ denotes ensemble average over stochastic wavefunction |*ϕ*(*t*, **z***)⟩ obtained by the normalizable (nonlinear) NMQSD equation,^{28,29}

where **z** comprises a set of complex Gaussian stochastic processes $zt,n*$ with mean $M[zt,n]=0$ and correlations $M[zt,nzs,m]=0$ and $M[zt,nzs,m*]=\alpha n(t\u2212s)\delta nm$. Here, *α*_{n}(*t*) is the correlation function of the environment, defined in Eq. (7). These processes enter via $\zeta n(t,z)=zt,n*+\u222b0tds\alpha n*(t\u2212s)\u3008L\u0302n\u2020\u3009s$, where the expectation values ⟨·⟩_{t} are calculated using the normalized state $|\varphi (t,z*)\u3009/\u3008\varphi (t,z*)|\varphi (t,z*)\u3009$.

For completeness, we mention that besides the nonlinear NMQSD equation (20), there exists also a *linear* NMQSD formulation where the nonlinear terms $\u3008L\u0302n\u2020\u3009$ are dropped in Eq. (20) and in *ζ*_{n}(*t*, **z**) and expectation values are calculated as $\u3008F\u0302\u3009(t)=M[\u3008\varphi (t,z*)|F\u0302|\varphi (t,z*)\u3009]$. However, the linear NMQSD equation converges slowly with the number of trajectories except for the case of weak system environment coupling or very short propagation times.

### B. Reformulation of the response function equations

Our aim is now to formulate the response function in a way that can be used together with the above nonlinear NMQSD equation. To focus on the main aspects of our treatment, we discuss first the initial state $\rho \u0302(0)=\varphi ini\varphi ini\u229700$, where the environment is in its ground state. In Subsection III D, we briefly comment on different, possibly mixed, states of the environment, with particular attention to the case of a thermal initial state.

We introduce $\Phi B(t)$ and $\Phi K(t)$, which represent the evolution of the *bra* and *ket* contributions, respectively. We can then write the *M*th-order density matrix as

where

and the last interval *τ*_{M} = *t* − *t*_{M} contains the time *t*. In this notation, the response function *R*(*t*), an abbreviation for *R*(*t*) = *R*^{(M)}(*τ*_{1}, …, *τ*_{M}), becomes

Both $\Phi B(t)$ and $\Phi K(t)$ can be expanded with respect to Bargmann coherent states^{38} of the bath,

where $\varphi B/K(t,z*)$ are states in the “system” Hilbert space only and $dM(z)=\Pi n\lambda d2zn\lambda e\u2212|zn\lambda |2\pi $. Inserting this expansion into Eq. (25),

and using the “reproducing property” of coherent states, we obtain the important result

where the *bra* and the *ket* now evolve with the same coherent states **z**.

Introducing a state

in a doubled “*system*” Hilbert space, we can write

with $F\u0303=0F\u030200$. These formulas are the starting point for the dyadic NMQSD approach similar to the doubling used for the case of a quantum state diffusion unravelling of Lindblad equations.^{39}

### C. Dyadic NMQSD

The dyadic NMQSD equations use the construction of the response functions in the doubled “system” Hilbert space to time-evolve the combined *bra* and *ket* states. We can introduce a state $\Psi \u0303(t)$ that lives in the Hilbert space $(HS\u2297HS)\u2297HB$ (note that the “bath” Hilbert space is not doubled) and obeys

where the left-hand side is the state introduced in Eq. (29). For the corresponding time-evolution, we can write

with an initial state

a time-evolution operator

and a Hamiltonian

where the *κ*_{nλ} are the same as in Eq. (3),

The response function [Eq. (30)] is then evaluated as an average over trajectories that obey the nonlinear NMQSD equation [Eq. (20)] in the doubled system Hilbert space. There are many formally equivalent dyadic NMQSD propagation schemes; we will discuss one explicit propagation scheme in Sec. III E, below.

While the *bra* and *ket* states are not directly coupled within the dyadic nonlinear NMQSD equations, they both contribute to the norm of the dyadic wavefunction and cannot be propagated independently. For example, note that the expectation values $\u3008L\u0303n\u2020\u3009t$, appearing in Eq. (20), are calculated using the dyadic wavefunction in the doubled system Hilbert space, where

and

We note that for the linear dyadic NMQSD equation, the bra and ket states are not coupled; however, these equations show poor convergence in parameter regimes beyond weak system–bath coupling.

### D. Comment on different initial states and finite temperature

The derivation above has been done assuming the environment initially in its ground state $(00)$. The “thermofield” approach^{40} is a well-known method for mapping the dynamics of a finite temperature environment onto a formally enlarged environment in its ground state. Already in early works, the thermofield technique has been applied to NMQSD (see Appendix C of Ref. 29), and in Ref. 41, the thermofield approach is discussed in detail for the case of molecular aggregates. For Hermitian coupling operators $Ln=Ln\u2020$, such as we have here, this approach leads to the NMQSD equation (20), where the initial thermal state is encoded in the temperature dependent bath-correlation function given by Eq. (7). A non-Hermitian coupling operator *L*_{n} introduces independent noise processes for *L*_{n} and $Ln\u2020$ with different correlation functions.^{29,41} Recent developments include alternatives to the thermofield approach, such as the P-representation of the thermal state^{42} that introduces a second noise term associated with the system Hamiltonian.

### E. Summary of the numerical propagation scheme

In the present work, the response function is calculated in the following way: We start with the system part of the initial state Eq. (33), which lives in the doubled system Hilbert space and which reads $\varphi \u0303ini=\varphi ini,\varphi iniT$. This state is not normalized. On this state, we act with $V\u03031$. We then propagate using the nonlinear (but unnormalized) NMQSD during the time interval *τ*_{1}. Then, we act with the second interaction $V\u03032$ and continue the propagation during time interval *τ*_{2}. We repeat this until the end of the last time interval and then calculate the expectation value of $F\u0303$ for each individual trajectory.

According to Eq. (19), we normalize each trajectory before taking the average. Here, some care is necessary. The normalization should take care of the change of norm caused by the unnormalized NMQSD propagation, but it includes in addition the norm changes due to the interactions $V\u03031,\u2026,V\u0303M$. To keep these physically relevant changes of the norm, we multiply by these norm changes. The details are as follows: Let us denote the state (in doubled system Hilbert space) before the *j*th interaction by $\psi \u0303(tj,z*)$, where *t*_{j} is the time of the *j*th interaction. We define

We then have for the “response function”

and we obtain the final result by averaging over trajectories,

We summarize this procedure in Fig. 1(b). We could have also used the *normalized* nonlinear NMQSD equation, with an appropriate change to the normalization factors at the end.

### F. HOPS for solving the NMQSD propagation

The NMQSD equation Eq. (20) is not particularly suitable for a direct numerical implementation because of the functional derivative with respect to the stochastic processes. Numerically convenient schemes can be derived when the bath-correlation function *α*_{n}(*t*) is expanded as a finite sum of exponentials,

with *w*_{nj} = *γ*_{nj} + *i*Ω_{nj}. In many applications of practical interest, the required number of “modes” (*N*_{modes}) is small. For the interpretation of such modes, see, for example, Refs. 37 and 43. A powerful, but approximate, scheme that is based on Eq. (43) is the so-called “zeroth order functional expansion,” ZOFE.^{30,31} In the present work, we employ the numerically exact hierarchy of pure states (HOPS),^{32}

Here, $w=w1,1,\u2026,wN,J$ and $k=k1,1,\u2026,kN,J$ with non-negative integers *k*_{nj}. The vector $enj=0,\u2026,1,\u2026,0$ is one at the (*n*, *j*)th position and is zero otherwise. The relevant contribution to perform calculations of expectation values is the zeroth order element, i.e.,

The HOPS consists of an infinite set of coupled equations, which must be truncated at a finite hierarchy for numerical calculations. In this work, we use a simple triangular truncation condition for the hierarchy: $|k|\u2264K$. More advanced truncation schemes are discussed in Ref. 44. It is also possible to use an adaptive algorithm to reduce the size of the hierarchy^{33} or use a matrix product state representation.^{34}

## IV. EXAMPLE CALCULATIONS

One reason for developing the present perturbative approach is the large number of trajectories required to converge the non-perturbative approach of Ref. 35. Therefore, in the following exemplary calculations, we focus in particular on the convergence with respect to the number of trajectories.

### A. Model system

Here, we perform calculations for a dimer (*N* = 2) consisting of identical monomers (i.e., *ϵ*_{n} = *ϵ*) with parallel transition dipoles (*μ*_{n} = ** μ**). For the bath-correlation functions, we choose a single exponential given by

with *γ* = 0.25 and the vibrational frequency Ω as the unit of energy. Then, Eq. (46) can be interpreted as a weakly damped vibrational mode at zero temperature,^{37} which requires a non-Markovian treatment. Below, we consider two values for the coupling strength that lead to qualitatively different 2D spectra: the intermediate coupling case and the strong coupling case (*p* = 0.5 and 1.8, respectively, in units of Ω^{2}). For the interaction between the monomers, we use *V* = 0.3. We note that the linear HOPS equation has severe convergence problems (see Sec. II of the supplementary material), as expected for our parameter regime.

### B. The various 2D spectra

In the following, we consider the third-order response functions that contribute to 2D electronic spectroscopy. For the three time intervals, we adopt the commonly used notation

We present plots for the ground state bleaching (GSB), stimulated emission (SE), and excited state absorption (ESA) signals as follows:

with

These expressions emerge after applying the rotating wave approximation and phase matching conditions.^{1,2} The functions *r*_{ℓ} appearing in Eq. (49) are specific response functions that are evaluated for $F=e\u22c5\mu \u0302\u2212$ and contain non-Hermitian interaction operators $V\u0302j\xb1=\u2212\mu \u0302\xb1Ej\xb1$, where

In Table I, we summarize the operators that enter into the calculation of the response functions *r*_{1} to *r*_{6}. We take all fields to be identical and polarized parallel to the transition dipole moments of the molecules.

$v1K$ | $v1B$ | $v2K$ | $v2B$ | $v3K$ | $v3B$ | |

r_{1} | $V\u03021+$ | $1$ | $1$ | $V\u03022+$ | $1$ | $V\u03023\u2212$ |

r_{2} | $1$ | $V\u03021+$ | $V\u03022+$ | $1$ | $1$ | $V\u03023\u2212$ |

r_{3} | $1$ | $V\u03021+$ | $1$ | $V\u03022\u2212$ | $V\u03023+$ | $1$ |

r_{4} | $V\u03021+$ | $1$ | $V\u03022\u2212$ | $1$ | $V\u03023+$ | $1$ |

r_{5} | $1$ | $V\u03021+$ | $V\u03022+$ | $1$ | $V\u03023+$ | $1$ |

r_{6} | $V\u03021+$ | $1$ | $1$ | $V\u03022+$ | $V\u03023+$ | $1$ |

$v1K$ | $v1B$ | $v2K$ | $v2B$ | $v3K$ | $v3B$ | |

r_{1} | $V\u03021+$ | $1$ | $1$ | $V\u03022+$ | $1$ | $V\u03023\u2212$ |

r_{2} | $1$ | $V\u03021+$ | $V\u03022+$ | $1$ | $1$ | $V\u03023\u2212$ |

r_{3} | $1$ | $V\u03021+$ | $1$ | $V\u03022\u2212$ | $V\u03023+$ | $1$ |

r_{4} | $V\u03021+$ | $1$ | $V\u03022\u2212$ | $1$ | $V\u03023+$ | $1$ |

r_{5} | $1$ | $V\u03021+$ | $V\u03022+$ | $1$ | $V\u03023+$ | $1$ |

r_{6} | $V\u03021+$ | $1$ | $1$ | $V\u03022+$ | $V\u03023+$ | $1$ |

### C. Calculations

In Fig. 2, we present the intermediate coupling case, *p* = 0.5. Figure 2(a) shows in the upper row, the GSB, SE, and ESA spectra obtained for T = 0 using *N*_{traj} = 1000 trajectories. We compare our HOPS spectra to reference calculations, performed with the HEOM method^{15,16} [Fig. 2(a), middle row]. We see that the HOPS spectra reproduce the relevant features from HEOM. To see in detail the deviations of the HOPS spectra from the HEOM ones, we show the point-wise difference between the HOPS spectrum and the reference HEOM spectrum [Fig. 2(a), bottom row]. From this, one sees that the maximal differences are around 10% of the peak signal. For GSB and ESA, the fluctuations are spread around a large region in the vicinity of the signal. For the SE spectrum, the fluctuations are largest along the diagonal. We note that the shown HOPS spectra and, therefore, the difference plots also depend on the specific realization of the *N*_{traj} trajectories.

To investigate in more detail the difference between the spectra, we introduce a measure $E$ for the integrated difference (details are given in Appendix B). This measure of the difference for the GSB, SE, and ESA [Fig. 2(a)] is given by $E=0.069$, $E=0.076$, and $E=0.085$, respectively. From this, we see that $E\u22720.08$ corresponds to fairly good agreement. Using this measure, we now investigate convergence with respect to the number of trajectories. In a first step, we construct the distribution of $E$ for different realizations of the *N*_{traj} using bootstrapping and then estimate the mean value and variance of the error (see Appendix B for details). In Fig. 2(b), we plot the mean and standard deviation of the error with respect to the number of trajectories. For each waiting time, the mean error follows the expected $1/Ntraj$ scaling (shown as a solid line), which can be clearly observed in the inset showing a double logarithmic scale. Furthermore, the standard deviation also decreases as $1/Ntraj$ and has values smaller than the mean error. Thus, the HOPS spectra associated with a large truncation depth of the hierarchy we have used here converge to the HEOM result and the associated error is dominated by statistical noise.

At waiting time *T* = 0, all three errors are comparable; upon increasing the waiting time, the error of the GSB signal remains essentially unchanged, while the error of the SE and ESA signal increases. In general, the SE (orange) and ESA (blue) errors are very similar and both are larger than the GSB error (brown). For the shown waiting times *T* ≈ 4, the SE and ESA error curves also remain largely unchanged and even decrease slightly. Details of the dependence of the error as a function of waiting time *T* are given in Sec. I of the supplementary material.

Figure 3 shows analogous results for the strong coupling case. In particular, there is again the $1/Ntraj$ scaling of the error and the error does not increase for waiting times *T* ≥ 2.

## V. CONCLUSIONS

In the present study, we have developed a framework to simulate multidimensional electronic spectra of molecular aggregates using the stochastic nonlinear formalism to directly calculate perturbative response functions of arbitrary order. Our approach calculates nonlinear response functions by propagating a pure state that combines *ket* and *bra* states in a doubled electronic Hilbert space subjected to a common noise, which corresponds to a dyadic equation. Importantly, our scheme enables us to use the nonlinear HOPS equation, which we have shown has superior convergence properties for strong system–bath coupling, compared to the linear version. The use of the nonlinear NMQSD and the appearance of a common noise distinguishes our scheme, e.g., from that of Ref. 45, where the 2D signal is also calculated using a kind of dyadic HOPS equation. In contrast to our pure state approach, there also exist stochastic schemes where the time-evolution of the reduced density matrix is calculated using dyadic propagation with correlated noise.^{46,47} It would be interesting to see how response functions obtained from the propagation scheme of Refs. 46 and 47 relate to the ones derived in the present work.

Numerical simulations for a dimer system coupled to a structured environment demonstrate that our theory has favorable convergence properties with respect to the number of stochastic trajectories. It should be noted that the new formalism developed here needs many fewer stochastic trajectories to obtain converged 2D spectra when compared to those calculated by a non-perturbative phase-cycling scheme.^{35} As compared to density matrix based methods, the nonlinear NMQSD method propagates vectors instead of matrices, individual simulations with different noise trajectories are trivially parallel, and it is consistent with adaptive basis^{33} and tensor contraction^{34} approaches recently developed for HOPS. Our theory, thus, offers a promising technique to simulate 2D spectra of large molecular aggregates. This is especially the case for the description of the excited state absorption contribution to 2D spectra where a large number of doubly excited excitonic states are involved. Furthermore, it is straightforward to account for the effect of the static disorder induced by the inhomogeneity of the solvent environment by simply sampling excitonic parameters from a certain distribution for each stochastic trajectory. Our theory can be readily applied to simulate higher-order response functions, for example, fifth-order 3D signals, which are a powerful tool to reveal multistep energy transfer processes.^{48}

## SUPPLEMENTARY MATERIAL

In the supplementary material, we provide further details on the convergence of the nonlinear HOPS with respect to the number of trajectories and demonstrate the convergence problems associated with the linear HOPS that are resolved using our method.

## ACKNOWLEDGMENTS

We thank Jacob K. Lynd for proofreading. L.C. acknowledges support from the Max-Planck Gesellschaft via the MPI-PKS visitors program. A.E. acknowledges support from the DFG via a Heisenberg fellowship (Grant No. EI 872/10-1). D.I.G.B. acknowledges support from the Robert A. Welch Foundation (Grant No. N-2026-20200401) and the US National Science Foundation CAREER Award under Grant No. CHE-2145358.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Lipeng Chen**: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Project administration (equal); Validation (lead); Writing – review & editing (equal). **Doran I. G. Bennett**: Conceptualization (equal); Investigation (supporting); Methodology (supporting); Project administration (supporting); Writing – review & editing (equal). **Alexander Eisfeld**: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (lead); Supervision (lead); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: LINEAR RESPONSE

It is instructive to also consider linear response within the present formalism to elucidate the normalization with respect to the state in doubled Hilbert space. The linear response function as it appears in the calculation of absorption is defined as

#### 1. Expression using the formalism of Sec. III

Within the dyadic NMQSD formalism, we write the linear response function as

where

This corresponds to the first-order response function with $v\u03021K=e\u22c5\mu \u0302+$, $v\u03021B=1$, and $F\u0302=e*\u22c5\mu \u0302\u2212$. The numerator of $I1$ simplifies to $\Vert V\u03031\psi \u0303(t1,z*)\Vert 2=1+\mu eff2$, with $\mu eff2=\u2211n(e\u22c5\mu n)2$. We calculate the final state $\psi \u0303(t,z*)$, defined in Eq. (29), following the prescription in Sec. III E: After the interaction of the initial vector $\psi \u0303(t1,z*)$ with the operator $V\u03031$, the state becomes $g(e\u22c5\mu \u0302+)|g\u3009$. During the subsequent time propagation (from *t*_{1} to *t*), the *bra* is in the ground state and only acquires a phase, $\varphi B(t)=e\u2212i\u03f5gtg$. The *ket* contribution $\varphi K(t)=\varphi K(t,z*)$ can be obtained from propagating the initial state $(e\u22c5\mu \u0302+)|g\u3009$ with the NMQDS equation in the *single* Hilbert space, where the expectation values of $L\u0302n\u2020$ at time *s*$(\u3008L\u0302n\u2020\u3009s)$ are calculated with respect to the norm $\Vert \psi \u0303(s,z*)\Vert 2=(\Vert \varphi K(s,z*)\Vert 2+1)$ of the state in the doubled Hilbert space,

Finally, the response function can be written as

where we have introduced $\psi ex=1\mu effe\u22c5\mu \u0302+g$, to make the connection to our previous result^{36} more obvious (see next subsection).

#### 2. Relation to previous results

In a previous publication,^{36} we have derived an equation for the perturbative calculation of the linear response function using the nonlinear NMQSD equation. In that work, the starting point was to treat the non-Hermitian operator $(e\u22c5\mu \u0302+)|g\u3009\u3008g|\rho \u0302B$ as “initial state,” which is then decomposed into a sum of pure states that can be propagated via NMQSD. In Ref. 36, the response was obtained from

Moreover, here, the state $\chi $ is propagated in the excited Hilbert space, but expectation values of $L\u0302n\u2020$ are calculated using the normalization with (‖*χ*(*t*, **z***)‖ + 1). From this, one sees that the only difference to the approach of Sec. III is that one starts the excited state propagation of the *ket* with a different normalized state. The change in initial condition leads to different trajectories even for the same noise realization. Nevertheless, both methods result in equivalent average response functions. Numerically, we have found that for our examples, there is little difference in the convergence of the two approaches.

To derive Eq. (A5) within our present formalism, we redefine the response function (Sec. II C) by scaling $V\u0302j\u2192V\u0302j/mj$, $F\u0302\u2192F\u0302/m$, and $R(M)\u2192m(\u220fjMmj)R(M)$. Using *m* = *μ*_{eff} and *m*_{j} = *μ*_{eff}, we have explicitly for the response function (A1), the expression

We then use the same steps as in the previous subsection and arrive at Eq. (A5).

### APPENDIX B: THE ERROR MEASURE

To quantify the difference between different 2D spectra, we introduce an error measure in the following way: First, we normalize each spectrum according to

where

with Ω = 1. For two spectra *S*^{(1)}(*ω*_{t}, *ω*_{τ}) and *S*^{(2)}(*ω*_{t}, *ω*_{τ}), we then introduce the difference

Finally, we define the *integrated difference*

To obtain a detailed analysis of the statistical error due to a finite number of trajectories shown in Figs. 2 and 3, we employ bootstrapping.^{49} We first calculate 4 × 10^{4} trajectories. For each value of *N*_{traj}, we then construct 500 ensembles by randomly choosing *N*_{traj} trajectories from the original 4 × 10^{4} trajectories. For each ensemble, we calculate the averaged *integrated difference* and finally obtain the *error* as the mean of the *integrated difference* over the 500 ensembles.