Calculating non-linear response functions for multi-dimensional electronic spectroscopy using dyadic non-Markovian quantum state diffusion

We present a methodology for simulating multi-dimensional 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 spectra of large molecular aggregates.

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][20][21][22][23] 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. 24,25 Over the years, sev-eral approaches have been developed to efficiently solve the NMQSD equation numerically, [26][27][28][29][30] 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 non-linear 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 femtosecondpulses 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. 31 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 non-linear 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. 32 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 section II, we introduce the details of the molecular system, its interaction with electromagnetic pulses, and the general form of the response functions. In section III, we develop our method to calculate the response function using the NMQSD approach. We par-ticularly emphasize the ability to use the non-linear NMQSD equation that ensures suitable convergence with respect to trajectories. In section 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 section 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. 32

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 Hamil-tonianĤ ex can then be written aŝ where ε n is the energy required to excite the nth molecule, V nm is the electronic coupling between excited molecules n and m, andσ n = g n e n . In the case of 2D spectroscopy, we need the common ground state (|g = ∏ N n=1 |g n ), singly excited states (|n = |e n ∏ m =n |g m =σ † n g ), and doubly excited states (|nm = |e n |e m ∏ k =n,m |g k =σ † nσ † m g , n < m) of the molecular aggregate.
For each molecule there are additional interactions with internal and external nuclear degrees of freedom. In many cases of interest these interactions can be modelled 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 aŝ Here,b nλ (b † nλ ) 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 aŝ where the coupling operatorL n acts in the system Hilbert space and is given byL and κ nλ is the exciton-bath coupling strength of the mode λ for molecule n, which is specified by the bath spectral density of molecule n, J n (ω) = ∑ λ |κ nλ | 2 δ (ω − ω nλ ). The latter is related to the bath-correlation function α n (t) by with the inverse temperature β = 1/T . We write the complete matter Hamiltonian aŝ

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Ĥ F (t) is defined as 2 where,μ is the total transition dipole operator and µ n the transition dipole moment of molecule n.
The electric field is given by e a E a (t − t a )e ik a ·r−iω a (t−t a ) + H.c. (9) 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 ath pulse, respectively. In general there can also a different number of pulses. We can write the total HamiltonianĤ tot aŝ withĤ given in Eq. (6).

C. Response functions
In this section we present a general notation for non-linear optical response functions that provide a clear connection to the NMQSD formalism. Here, we use a general notation that that is appropriate for arbitrary orders of perturbation theory. In section IV, we specify to 2D spectroscopy.

Perturbation theory for the full density matrix
In the following we denote the time-evolution operator of the system without the electromagnetic field aŝ whereĤ given in Eq. (6). 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) which 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 observableF, 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. (13) 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 row 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 bŷ In Eq. (14) the operators v K j act always on the ket side and the operators v B j always on the bra side ofρ. An important constraint is that for each pair v K j , v B j , with the same index j one of the corresponding operators is the unit operator (which we denote by 1 1). Note that for better readability we will sometimes omit the 'operator hats' on the operators v K/B j , as it has already been done above. Explicitly we have In the following, we introduce a short-hand 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.

A. The NMQSD formalism
For the open quantum system modelĤ =Ĥ ex +Ĥ B +Ĥ int as given in section II A and for a factorized initial stateρ ini = |φ ini φ ini |⊗ρ env , the expectation value of any system operator F can be obtained as 24,25 where M[· · · ] denotes ensemble average over stochastic wave function |φ (t, z * ) obtained by the normalizable (non-linear) NMQSD equation 24,25 where z comprises a set of complex Gaussian stochastic processes z * t,n with mean M[z t,n ] = 0, and correlations M[z t,n z s,m ] = 0 and M[z t,n z * s,m ] = α n (t − s)δ nm . Here α n (t) is the correlation function of the environment, defined in Eq. (5). These processes enter via ζ n (t, z) = z * t,n + t 0 ds α * n (t − s) L † n s , where the expectation values · t are calculated using the normalized state |φ (t, z * ) / φ (t, z * )|φ (t, z * ) .
For completeness we mention that beside the non-linear NMQSD equation (19), there exists also a linear NMQSD formulation where the non-linear terms L † n are dropped in Eq. (19) and in ζ n (t, z) and expectation values are calculated as 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 non-linear NMQSD equation. We first note that in the derivation of the NMQSD equation the finite temperature initial state of the environment has been transformed to the vacuum states by shifting the contributions of temperature to the bath-correlation function. Therefore, the following derivation is done for the initial statê We introduce Φ B (t) and Φ K (t) which represent the evolution of the bra and ket contributions, respectively. We can then write the M th order density matrix aŝ where and the last interval τ M = t − t M contains the time t. In this notation, the response function R(t), an abbreviation for Both Φ B (t) and Φ K (t) can be expanded with respect to coherent states of the bath where φ B/K (t, z * ) are states in the 'system' Hilbert space only and dM(z) = Π nλ d 2 z nλ e −|z nλ | 2 π . Inserting this expansion into Eq. (24) 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 = 0F 0 0 . 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. 33

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 Ψ(t) that lives in the Hilbert space (H S ⊗ H S ) ⊗ H B (note that the 'bath' Hilbert space is not doubled) and obeys where the left hand side is the state introduced in Eq. (28). 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), , and L n = L n 0 0L n .

5
The response function (Eq. (29)) is then evaluated as an average over trajectories which obey the nonlinear NMQSD equation (Eq. (19)) in the doubled system Hilbert space. There are many formally equivalent dyadic NMQSD propagation schemes; we will discuss one explicit propagation scheme in Section III D, below. While the bra and ket states are not directly coupled within the dyadic non-linear NMQSD equations, they both contribute to the norm of the dyadic wave function and cannot be propagated independently. For example, note that the expectation values L † n t , that appears in Eq. (19), are calculated using the dyadic wave function 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. 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. (32), which lives in the doubled system Hilbert space and which reads φ ini = φ ini , φ ini T . This state is not normalized. On this state we act with V 1 . We then propagate using the non-linear (but unnormalized) NMQSD during the time interval τ 1 . Then we act with the second interaction V 2 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 for each individual trajectory. According to Eq. (18) 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 un-normalized NMQSD propagation but it includes in addition the norm changes due to the interactions V 1 , . . . , V M . To keep these physically relevant changes of the norm we multiply by these norm changes. In detail: Let us denote the state (in doubled system Hilbert space) before the j-th interaction by ψ(t j , z * ), where t j is the time of the j-th interaction. We define We then have for the 'response function' and obtain the final result by averaging over trajectories We summarize this procedure in Fig. 1(b). We could have also used the normalized non-linear NMQSD equation, with an appropriate change to the normalization factors at the end.

E. HOPS for solving the NMQSD propagation
The NMQSD equation Eq. (19) 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 bathcorrelation function α n (t) is expanded as a finite sum of exponentials, with w n j = γ n j + iΩ n j . In many applications of practical interest, the required number of 'modes' (N modes ) is small. For the interpretation of such modes see for examples. 34,35 A powerful, but approximate, scheme that is based on Eq. (42) is the so called 'zeroth order functional expansion' ZOFE. 26,27 In the present work we employ the numerically exact hierarchy of pure states (HOPS) 28 (43) w = {w 1,1 , · · · , w N,J }, and k = {k 1,1 , · · · , k N,J } with nonnegative integers k n j . The vector e n j = {0, · · · , 1, · · · , 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| ≤ K. More advanced truncation schemes are discussed in Ref. 36. It is also possible to use an adaptive algorithm to reduce the size of the hierarchy 29 or use a matrix product state representation. 30

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. 31. 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 with γ = 0.25 and the vibrational frequency Ω as the unit of energy. Then Eq. (45) can be interpreted as a weakly damped vibrational mode at zero temperature, 34 which requires a non-Markovian treatment. Below, we consider two values for the coupling strength that lead to qualitative 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.

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.
These expressions emerge after applying the rotating wave approximation and phase matching conditions. 1,2 The functions r appearing in Eq. (48) are specific response functions that are evaluated for F = e ·μ − , and contain non-hermitean interaction operatorsV ± j = −μ ± E ± j , wherê 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.  Table I. The operators used in the calculation of the response functions r 1 , . . . r 6 . HereV ± j = −μ ± E ± j withμ ± and E ± j given in Eqs. (49) and (50).

C. Calculations
In Fig. 2 we present the intermediate coupling case, p = 0.5. Fig. 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. 2A, middle row). We see that the HOPS spectra reproduces 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. 2A, 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.
To investigate the convergence with the number of trajectories we introduce a measure E for the integrated difference (details are given in appendix B). This measure of the difference for the GSB, ES, and ESA ( Fig. 2A) are given by E = 0.069, E = 0.076 and E = 0.085, respectively. From this we see that E 0.08 corresponds to quite good agreement. In Fig. 2B, we plot a non-parametric error estimate (calculated via bootstrapping, Appendix B) with respect to the number of trajectories. The SE (orange) and ESA (blue) errors are very similar and both are larger than the GSB error (brown). 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. For the shown waiting times T ≈ 4 the SE and ESA error curves also remain largely unchanged and even decrease slightly. For each waiting time, the error follows the expected 1/ N traj scaling (shown as a solid line), which can be clearly observed in the inset showing a double logarithmic scale. Fig. 3 shows analogous results for the strong coupling case. In particular, there is again the 1/ N traj scaling of the error and the error does not increase for waiting times T ≥ 2.

V. CONCLUSIONS
In the present work, we have developed a framework to simulate multidimensional electronic spectra of molecular aggre- gates using the stochastic non-linear formalism to directly calculate perturbative response functions of arbitrary order. Our approach propagates a dyadic wave function which combines ket and bra states in a doubled system Hilbert space subjected to a common noise.
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. 31 As compared to density matrix based methods, the non-linear NMQSD method propagates vectors instead of matrices, individual simulations with different noise trajectories are trivially parallel and, furthermore, it is consistent with adaptive basis 29 and tensor contraction 30 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 multi-step energy transfer processes. 37 Figure 3. Same as Fig. 2 but for the strong system-bath coupling regime (p = 1.8). The truncation conditions for HOPS are K = 15 for the GSB and SE, and K = 20 for the ESA.

DATA AVAILABILITY
Further 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 R (1) (t) = Tr (e * ·μ − )e −iĤt (e ·μ + )|g g|ρ B e iĤt . (A1)

Expression using the formalism of section III
Within the dyadic NMQSD formalism we write the linear response function as where , and ψ(t 1 , z * ) = g g .
This corresponds to the first-order response function witĥ v K 1 = e ·μ + ,v B 1 = 1 1 and F = e * ·μ − . The numerator of I 1 simplifies to || V 1 ψ(t 1 , z * )|| 2 = 1 + µ 2 eff , with µ 2 eff = ∑ n (e · µ n ) 2 . We calculate the final state ψ(t, z * ) , defined in Eq. (28), following the prescription in Section III D: After the interaction of the initial vector ψ(t 1 , z * ) with the operator V 1 the state becomes g (e ·μ + )|g . During the subsequent time propagation (from t 1 to t), the bra is in the ground state and only acquires a phase, φ B (t) = e −iε g t g . The ket contribution φ K (t) = φ K (t, z * ) can be obtained from propagating the initial state (e ·μ + )|g with the NMQDS equation in the single Hilbert space, where the expectation values of L † n at time s ( L † n s ) are calculated with respect to the norm || ψ(s, z * )|| 2 = (||φ K (s, z * )|| 2 + 1) of the state in the doubled Hilbert space Finally the response function can be written as (A4) where we have introduced ψ ex = 1 µ eff e ·μ + g , to make the connection to our previous result 32 more obvious (see next subsection).

Relation to previous results
In a previous publication 32 we have derived an equation for the perturbative calculation of the linear response function using the non-linear NMQSD equation. In that work the starting point was to treat the non-Hermitean operator (e·μ + )|g g|ρ B as 'initial state', which is then decomposed into a sum of pure states which can be propagated via NMQSD. In Ref. 32 the response was obtained from R decomp (t) = µ 2 eff M ψ ex |χ(t, z * ) 1 2 (||χ(t, z * )|| + 1) e iε g t .
Also here the state χ is propagated in the excited Hilbert space, but expectation values ofL † n are calculated using the normalization with (||χ(t, z * )|| + 1). From this one sees that the only difference to the approach of section 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 (section II C) by scalingV j → V j /m j ,F →F/m and R (M) → m (∏ M j m j ) R (M) . Using m = µ eff and m j = µ eff we have explicitly for the response function (A1) (A6) We then use the same steps as in the previous subsection and arrive at Eq. (A5). with Ω = 1. For two spectra S (1) (ω t , ω τ ) and S (2) (ω t , ω τ ) we then introduce the difference ∆S(ω t , ω τ ) = S (1) (ω t , ω τ ) − S (2) (ω t , ω τ ) Finaly, 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. 38 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.