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.

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.

We consider a molecular aggregate composed of N interacting molecules, where each molecule is described by two electronic levels: the electronic ground state |gn⟩ and the electronic excited state |en⟩, n = 1, …, N. The electronic Hamiltonian Ĥex can then be written as

Ĥex=nϵnσ̂nσ̂n+nmnVnmσ̂nσ̂m,
(1)

where ϵn is the energy required to excite the nth molecule, Vnm is the electronic coupling between excited molecules n and m, and σ̂n=gnen. In the case of 2D spectroscopy, we need the common ground state (|g=n=1N|gn), singly excited states (|n=|enmn|gm=σ̂ng), and doubly excited states (|nm=|en|emkn,m|gk=σ̂nσ̂mg, 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

ĤB=n=1Nλωnλb̂nλb̂nλ.
(2)

Here, b̂nλ(b̂nλ) is the annihilation (creation) operator of λth bath mode of molecule n with frequency ω. The bath modes couple locally to their respective molecule. The interaction Hamiltonian is then written as

Ĥint=n=1NL̂nλκnλ(b̂nλ+b̂nλ),
(3)

where the coupling operator L̂n acts in the system Hilbert space and is given by

L̂n=σ̂nσ̂n
(4)

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

Ĥ=Ĥex+ĤB+Ĥint.
(5)

The detailed derivation of this Hamiltonian can be found in Refs. 2, 3, and 37.

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

ρ̂ini=|gg|ρ̂env,withρ̂env=eĤenv/Tr[eĤenv].
(6)

It is convenient to introduce the bath spectral density of molecule n, Jn(ω) = ∑λ|κ|2δ(ωω), which is related to the bath-correlation function αn(t) by

αn(t)=0dωJn(ω)cothβω2cos(ωt)isin(ωt),
(7)

with the inverse temperature β.

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 as2 

ĤF(t)=μ̂E(r,t),
(8)

where μ̂ is the total transition dipole operator given by

μ̂=n=1Nμn(σ̂n+σ̂n)
(9)

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

E(r,t)=a=13eaEa(tta)eikariωa(tta)+H.c.,
(10)

with ea, ka, ωa, Ea(t), and ta 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 exist a different number of pulses.

We can write the total Hamiltonian Ĥtot as

Ĥtot=Ĥ+ĤF(t),
(11)

with Ĥ given in Eq. (5).

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

Û(τ)=Ûτ=eiĤτ,
(12)

where Ĥ is given in Eq. (5). We also introduce the abbreviation

V̂j=ĤF(tj)
(13)

for the interaction with field at time tj.

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

Rτ1,,τMv1K,,vMKv1B,,vMB=TrF̂ρ̂(M)τ1,,τMv1K,,vMKv1B,,vMB,
(14)

where the expectation value of an observable F̂, in our case the polarization μ, is calculated with respect to a density matrix that is obtained in Mth 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 = ti+1ti between interaction at time ti and the following interaction. The last interval τM contains the final time t of the evolution, i.e., τM = ttM. 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 Mth order density matrix is recursively defined by

ρ̂(j)τ1,,τjv1K,,vjKv1B,,vjB=Û(τj)v̂jKρ̂(j1)τ1,,τj1v1K,,vj1Kv1B,,vj1Bv̂jBÛ(τj)
(15)

and

ρ̂(0)=ϕiniϕiniρenv.
(16)

In Eq. (15), the operators vjK act always on the ket side and the operators vjB always on the bra side of ρ̂. 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

vjK=V̂j,vjB=1,orvjK=1,vjB=V̂j.
(17)

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

ρ̂(M)=ρ̂(M)τ1,,τMv1K,,vMKv1B,,vMB,
(18)

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

For the open-quantum system model as given in Sec. II A, the expectation value of any system operator F̂ can be obtained as28,29

F̂(t)Tr[F̂ρ̂(t)]=Mϕ(t,z*)|F̂|ϕ(t,z*)ϕ(t,z*)|ϕ(t,z*),
(19)

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

t|ϕ(t,z*)=iĤex|ϕ(t,z*)+nL̂nζn(t,z)|ϕ(t,z*)nL̂nL̂nt0tdsαn(ts)δ|ϕ(t,z*)δzs,n*,
(20)

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*]=αn(ts)δnm. Here, αn(t) is the correlation function of the environment, defined in Eq. (7). These processes enter via ζn(t,z)=zt,n*+0tdsαn*(ts)L̂ns, where the expectation values ⟨·⟩t are calculated using the normalized state |ϕ(t,z*)/ϕ(t,z*)|ϕ(t,z*).

For completeness, we mention that besides the nonlinear NMQSD equation (20), there exists also a linear NMQSD formulation where the nonlinear terms L̂n are dropped in Eq. (20) and in ζn(t, z) and expectation values are calculated as F̂(t)=M[ϕ(t,z*)|F̂|ϕ(t,z*)]. 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.

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 ρ̂(0)=ϕiniϕini00, 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 ΦB(t) and ΦK(t), which represent the evolution of the bra and ket contributions, respectively. We can then write the Mth-order density matrix as

ρ̂(M)(t)=ΦK(t)ΦB(t),
(21)

where

ΦB(t)=(eiĤτMv̂MB)(eiĤτ1v̂1B)ϕini0,
(22)
ΦK(t)=(eiĤτMv̂MK)(eiĤτ1v̂1K)ϕini0,
(23)

and the last interval τM = ttM contains the time t. In this notation, the response function R(t), an abbreviation for R(t) = R(M)(τ1, …, τM), becomes

R(t)=Tr{F̂ΦK(t)ΦB(t)}
(24)
=ΦB(t)F̂ΦK(t).
(25)

Both ΦB(t) and ΦK(t) can be expanded with respect to Bargmann coherent states38 of the bath,

ΦB/K(t)=dM(z)zz|ΦB/K(t)=dM(z)zϕB/K(t,z*),
(26)

where ϕB/K(t,z*) are states in the “system” Hilbert space only and dM(z)=Πnλd2znλe|znλ|2π. Inserting this expansion into Eq. (25),

R(t)=dM(z)dM(z)z|zϕB(t,z*)F̂ϕK(t,z*),
(27)

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

R(t)=dM(z)ϕB(t,z*)F̂ϕK(t,z*),
(28)

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

Introducing a state

ψ̃(t,z*)=ϕB(t,z*)ϕK(t,z*)
(29)

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

R(t)=dM(z)ψ̃(t,z*)F̃ψ̃(t,z*),
(30)

with F̃=0F̂00. 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 

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 (HSHS)HB (note that the “bath” Hilbert space is not doubled) and obeys

ψ̃(t,z*)=dM(z)zz|Ψ̃(t),
(31)

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

Ψ̃(t)=Ũ(τM)ṼMŨ(τ1)Ṽ1Ψ̃(0),
(32)

with an initial state

Ψ̃(0)=ϕiniϕini0,
(33)

a time-evolution operator

Ũ(τ)=eiH̃τ,
(34)

and a Hamiltonian

H̃=H̃S+ĤB+nL̃nλκnλ(b̂nλ+b̂nλ),
(35)

where the κ are the same as in Eq. (3),

H̃S=ĤS00ĤS,Ṽj=v̂jB00v̂jK,andL̃n=L̂n00L̂n.
(36)

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 L̃nt, appearing in Eq. (20), are calculated using the dyadic wavefunction in the doubled system Hilbert space, where

L̃nt=ψ̃(t,z*)L̃nψ̃(t,z*)ψ̃(t,z*)|ψ̃(t,z*)
(37)
=ϕK(t,z*)L̂nϕK(t,z*)+ϕB(t,z*)L̂nϕB(t,z*)ϕK(t,z*)2+ϕB(t,z*)2
(38)

and

ϕ(t,z*)2=ϕ(t,z*)|ϕ(t,z*).
(39)

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.

The derivation above has been done assuming the environment initially in its ground state (00). The “thermofield” approach40 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, 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 Ln introduces independent noise processes for Ln and Ln with different correlation functions.29,41 Recent developments include alternatives to the thermofield approach, such as the P-representation of the thermal state42 that introduces a second noise term associated with the system Hamiltonian.

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 ϕ̃ini=ϕini,ϕiniT. This state is not normalized. On this state, we act with Ṽ1. We then propagate using the nonlinear (but unnormalized) NMQSD during the time interval τ1. Then, we act with the second interaction Ṽ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. (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 Ṽ1,,ṼM. 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 jth interaction by ψ̃(tj,z*), where tj is the time of the jth interaction. We define

Ij=Ṽjψ̃(tj,z*)2ψ̃(tj+1,z*)2.
(40)

We then have for the “response function”

R(M)(z)=j=1MIjψ̃(t,z*)F̃ψ̃(t,z*),
(41)

and we obtain the final result by averaging over trajectories,

R(M)=MzR(M)(z).
(42)

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.

FIG. 1.

(a) Sketch of the time sequence of the interaction with three short laser pulses, centered at times t1, t2, and t3. In the present work, we focus on the impulsive limit, where the pulses are described by δ-functions. After the third pulse, a specific signal is detected at time tS. The resulting expressions can be readily applied also to pulses with finite width.2 Time intervals between two successive pulses are labeled by τj and indexed according to the index of the first of the two pulses. (b) Symbolic representation of the calculation of the response function according to the formalism of Sec. III B.

FIG. 1.

(a) Sketch of the time sequence of the interaction with three short laser pulses, centered at times t1, t2, and t3. In the present work, we focus on the impulsive limit, where the pulses are described by δ-functions. After the third pulse, a specific signal is detected at time tS. The resulting expressions can be readily applied also to pulses with finite width.2 Time intervals between two successive pulses are labeled by τj and indexed according to the index of the first of the two pulses. (b) Symbolic representation of the calculation of the response function according to the formalism of Sec. III B.

Close modal

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,

αn(t)j=1Nmodespnjewnjt,
(43)

with wnj = γnj + iΩnj. In many applications of practical interest, the required number of “modes” (Nmodes) 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 

t|ψ(k)(t,z*)=iĤexkw+nL̂nζn(t,z)|ψ(k)(t,z*)+nL̂njknjpnj|ψ(kenj)(t,z*)nL̂nL̂ntj|ψ(k+enj)(t,z*).
(44)

Here, w=w1,1,,wN,J and k=k1,1,,kN,J with non-negative integers knj. The vector enj=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.,

ϕ(t,z*)=|ψ(0)(t,z*).
(45)

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. 44. It is also possible to use an adaptive algorithm to reduce the size of the hierarchy33 or use a matrix product state representation.34 

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.

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

αn(t)=α(t)=peiΩtγt,t0,
(46)

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.

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

τ1τ,τ2T,τ3t.
(47)

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

SGSB(ωτ,T,ωt)=S3()(ωτ,T,ωt)+S4(+)(ωτ,T,ωt),SSE(ωτ,T,ωt)=S2()(ωτ,T,ωt)+S1(+)(ωτ,T,ωt),SESA(ωτ,T,ωt)=S5()(ωτ,T,ωt)+S6(+)(ωτ,T,ωt),
(48)

with

S(±)(ωτ,T,ωt)=Re00dtdτr(τ,T,t)e±iωττeiωtt.
(49)

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μ̂ and contain non-Hermitian interaction operators V̂j±=μ̂±Ej±, where

μ̂+=nμnσ̂n,μ̂=nμnσ̂n,
(50)
Ej+=ej,Ej=ej*.
(51)

In Table I, we summarize the operators that enter into the calculation of the response functions r1 to r6. 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 r1, …, r6. Here, V̂j±=μ̂±Ej±, with μ̂± and Ej± given in Eqs. (50) and (51).

 v1K v1B v2K v2B v3K v3B 
r1 V̂1+ 1 1 V̂2+ 1 V̂3 
r2 1 V̂1+ V̂2+ 1 1 V̂3 
r3 1 V̂1+ 1 V̂2 V̂3+ 1 
r4 V̂1+ 1 V̂2 1 V̂3+ 1 
r5 1 V̂1+ V̂2+ 1 V̂3+ 1 
r6 V̂1+ 1 1 V̂2+ V̂3+ 1 
 v1K v1B v2K v2B v3K v3B 
r1 V̂1+ 1 1 V̂2+ 1 V̂3 
r2 1 V̂1+ V̂2+ 1 1 V̂3 
r3 1 V̂1+ 1 V̂2 V̂3+ 1 
r4 V̂1+ 1 V̂2 1 V̂3+ 1 
r5 1 V̂1+ V̂2+ 1 V̂3+ 1 
r6 V̂1+ 1 1 V̂2+ V̂3+ 1 

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 Ntraj = 1000 trajectories. We compare our HOPS spectra to reference calculations, performed with the HEOM method15,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 Ntraj trajectories.

FIG. 2.

Two-dimensional spectra and their convergence properties for a dimer with two identical monomers and intermediate system–bath coupling (p = 0.5). The other parameters are γ = 0.25 and V = 0.3, where the energies are in units of Ω. For p = 0.5, the truncation conditions for HOPS are K=10 for the GSB and SE and they are K=11 for the ESA. (a) GSB, SE, and ESA spectra at waiting time T = 0. Upper row: HOPS calculations with Ntraj = 1000 trajectories. Second row: HEOM calculations as reference. Third row: point-wise difference SHOPS(ωτ, ωt) − SHEOM(ωτ, ωt) between the HOPS and the HEOM spectra. Note that for convenience, we have normalized all spectra to the maximal value a = max[S(ωτ, ωt)], which is indicated in each panel. This allows us to use a common colorbar for all plots and to quickly assess the maximal error. (b) Convergence with the number of trajectories for different waiting times T. In each panel, the error measure, defined in  Appendix B, is shown as function of number of trajectories. The crosses are the numerical data obtained from a bootstrapping procedure (described in  Appendix B) and the solid lines are curves that follow a 1/Ntraj scaling. GSB: brown, SE: orange, ESA: blue.

FIG. 2.

Two-dimensional spectra and their convergence properties for a dimer with two identical monomers and intermediate system–bath coupling (p = 0.5). The other parameters are γ = 0.25 and V = 0.3, where the energies are in units of Ω. For p = 0.5, the truncation conditions for HOPS are K=10 for the GSB and SE and they are K=11 for the ESA. (a) GSB, SE, and ESA spectra at waiting time T = 0. Upper row: HOPS calculations with Ntraj = 1000 trajectories. Second row: HEOM calculations as reference. Third row: point-wise difference SHOPS(ωτ, ωt) − SHEOM(ωτ, ωt) between the HOPS and the HEOM spectra. Note that for convenience, we have normalized all spectra to the maximal value a = max[S(ωτ, ωt)], which is indicated in each panel. This allows us to use a common colorbar for all plots and to quickly assess the maximal error. (b) Convergence with the number of trajectories for different waiting times T. In each panel, the error measure, defined in  Appendix B, is shown as function of number of trajectories. The crosses are the numerical data obtained from a bootstrapping procedure (described in  Appendix B) and the solid lines are curves that follow a 1/Ntraj scaling. GSB: brown, SE: orange, ESA: blue.

Close modal

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 E0.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 Ntraj 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.

FIG. 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.

FIG. 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.

Close modal

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 basis33 and tensor contraction34 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 

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.

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.

The authors have no conflicts to disclose.

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).

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

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*μ̂)eiĤt(eμ̂+)|gg|ρ̂BeiĤt.
(A1)

1. Expression using the formalism of Sec. III

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

R(1)(t)=MI1ψ̃(t,z*)|F̃|ψ̃(t,z*),
(A2)

where

F̃=0(e*μ̂)00,I1=Ṽ1ψ̃(t1,z*)2ψ̃(t,z*)2,
Ṽ1=100(eμ̂+),andψ̃(t1,z*)=gg.

This corresponds to the first-order response function with v̂1K=eμ̂+, v̂1B=1, and F̂=e*μ̂. The numerator of I1 simplifies to Ṽ1ψ̃(t1,z*)2=1+μeff2, with μeff2=n(eμn)2. We calculate the final state ψ̃(t,z*), defined in Eq. (29), following the prescription in Sec. III E: After the interaction of the initial vector ψ̃(t1,z*) with the operator Ṽ1, the state becomes g(eμ̂+)|g. During the subsequent time propagation (from t1 to t), the bra is in the ground state and only acquires a phase, ϕB(t)=eiϵgtg. 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̂ns) are calculated with respect to the norm ψ̃(s,z*)2=(ϕK(s,z*)2+1) of the state in the doubled Hilbert space,

L̂ns=ϕK(s,z*)L̂nϕK(s,z*)/(ϕK(s,z*)2+1).
(A3)

Finally, the response function can be written as

R(1)(t)=M1+μeff2ϕK(t,z*)2+1μeffψex|ϕK(t,z*)eiϵgt,
(A4)

where we have introduced ψex=1μeffeμ̂+g, to make the connection to our previous result36 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μ̂+)|gg|ρ̂B 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

Rdecomp(t)=μeff2Mψex|χ(t,z*)12(χ(t,z*)+1)eiϵgt.
(A5)

Moreover, here, the state χ is propagated in the excited Hilbert space, but expectation values of L̂n 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̂jV̂j/mj, F̂F̂/m, and R(M)m(jMmj)R(M). Using m = μeff and mj = μeff, we have explicitly for the response function (A1), the expression

R(1)(t)=μeff2Tr(e*μ̂)μeffeiĤt(eμ̂+)μeff|gg|ρ̂BeiĤt.
(A6)

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

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

S(ωt,ωτ)S(ωt,ωτ)/S,
(B1)

where

S=1Ωωminωmaxdωtωminωmaxdωτ|S(ωt,ωτ)|,
(B2)

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,ωτ).
(B3)

Finally, we define the integrated difference

E=ωminωmaxdωtωminωmaxdωτ|ΔS(ωt,ωτ)|.
(B4)

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 × 104 trajectories. For each value of Ntraj, we then construct 500 ensembles by randomly choosing Ntraj trajectories from the original 4 × 104 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.

1.
D.
Abramavicius
,
B.
Palmieri
,
D. V.
Voronine
,
F.
Šanda
, and
S.
Mukamel
,
Chem. Rev.
109
,
2350
2408
(
2009
).
2.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
1995
).
3.
V.
May
and
O.
Kühn
,
Charge and Energy Transfer Dynamics in Molecular Systems
, 3rd ed. (
Wiley-VCH
,
2011
).
4.
H.
van Amerongen
,
R.
van Grondelle
, and
L.
Valkunas
,
Photosynthetic Excitons
(
World Scientific
,
2000
).
5.
L.
Chen
,
P.
Shenai
,
F.
Zheng
,
A.
Somoza
, and
Y.
Zhao
,
Molecules
20
,
15224
15272
(
2015
).
6.
M.
Cho
,
Chem. Rev.
108
,
1331
1418
(
2008
).
7.
M.
Cho
,
Two-Dimensional Optical Spectroscopy
(
CRC Press
,
New York
,
2009
).
8.
J.
Dostál
,
F.
Fennel
,
F.
Koch
,
S.
Herbst
,
F.
Würthner
, and
T.
Brixner
,
Nat. Commun.
9
,
2466
(
2018
).
9.
G. D.
Scholes
,
G. R.
Fleming
,
L. X.
Chen
,
A.
Aspuru-Guzik
,
A.
Buchleitner
,
D. F.
Coker
,
G. S.
Engel
,
R.
van Grondelle
,
A.
Ishizaki
,
D. M.
Jonas
,
J. S.
Lundeen
,
J. K.
McCusker
,
S.
Mukamel
,
J. P.
Ogilvie
,
A.
Olaya-Castro
,
M. A.
Ratner
,
F. C.
Spano
,
K. B.
Whaley
, and
X.
Zhu
,
Nature
543
,
647
656
(
2017
).
10.
H.-G.
Duan
,
V. I.
Prokhorenko
,
R. J.
Cogdell
,
K.
Ashraf
,
A. L.
Stevens
,
M.
Thorwart
, and
R. J. D.
Miller
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
8493
8498
(
2017
).
11.
J.
Cao
,
R. J.
Cogdell
,
D. F.
Coker
,
H. G.
Duan
,
J.
Hauer
,
U.
Kleinekathöfer
,
T. L. C.
Jansen
,
T.
Mancǎl
,
R. J. D.
Miller
,
J. P.
Ogilvie
,
V. I.
Prokhorenko
,
T.
Renger
,
H. S.
Tan
,
R.
Tempelaar
,
M.
Thorwart
,
E.
Thyrhaug
,
S.
Westenhoff
, and
D.
Zigmantas
,
Sci. Adv.
6
,
eaaz4888
(
2020
).
12.
G. S.
Schlau-Cohen
,
A.
Ishizaki
, and
G. R.
Fleming
,
Chem. Phys.
386
,
1
22
(
2011
).
13.
A. G.
Redfield
,
IBM J. Res. Dev.
1
,
19
(
1957
).
14.
W. M.
Zhang
,
T.
Meier
,
V.
Chernyak
, and
S.
Mukamel
,
J. Chem. Phys.
108
,
7763
7774
(
1998
).
15.
Y.
Tanimura
,
J. Phys. Soc. Jpn.
75
,
082001
(
2006
).
16.
Y.
Tanimura
,
J. Chem. Phys.
153
,
020901
(
2020
).
17.
N.
Makri
and
D. E.
Makarov
,
J. Chem. Phys.
102
,
4600
4610
(
1995
).
18.
N.
Makri
and
D. E.
Makarov
,
J. Chem. Phys.
102
,
4611
4618
(
1995
).
19.
A.
Ishizaki
and
G. R.
Fleming
,
Proc. Natl. Acad. Sci. U. S. A.
106
,
17255
17260
(
2009
).
20.
C.
Kreisbeck
,
T.
Kramer
,
M.
Rodríguez
, and
B.
Hein
,
J. Chem. Theory Comput.
7
,
2166
2174
(
2011
).
21.
C.
Kreisbeck
and
T.
Kramer
,
J. Phys. Chem. Lett.
3
,
2828
2833
(
2012
).
22.
L.
Chen
,
R.
Zheng
,
Q.
Shi
, and
Y.
Yan
,
J. Chem. Phys.
132
,
024505
(
2010
).
23.
Y.
Tanimura
,
J. Chem. Phys.
137
,
22A550
(
2012
).
24.
J.-J.
Ding
,
J.
Xu
,
J.
Hu
,
R.-X.
Xu
, and
Y.
Yan
,
J. Chem. Phys.
135
,
164107
(
2011
).
25.
A. G.
Dijkstra
and
Y.
Tanimura
,
J. Chem. Phys.
142
,
212423
(
2015
).
26.
S.
Kundu
and
N.
Makri
,
J. Phys. Chem. Lett.
11
,
8783
8789
(
2020
).
27.
X.-T.
Liang
,
J. Chem. Phys.
141
,
044116
(
2014
).
28.
L.
Diósi
and
W. T.
Strunz
,
Phys. Lett. A
235
,
569
573
(
1997
).
29.
L.
Diósi
,
N.
Gisin
, and
W. T.
Strunz
,
Phys. Rev. A
58
,
1699
(
1998
).
30.
T.
Yu
,
L.
Diósi
,
N.
Gisin
, and
W. T.
Strunz
,
Phys. Rev. A
60
,
91
(
1999
).
31.
J.
Roden
,
A.
Eisfeld
,
W.
Wolff
, and
W. T.
Strunz
,
Phys. Rev. Lett.
103
,
058301
(
2009
).
32.
D.
Suess
,
A.
Eisfeld
, and
W. T.
Strunz
,
Phys. Rev. Lett.
113
,
150403
(
2014
).
33.
L.
Varvelo
,
J. K.
Lynd
, and
D. I. G.
Bennett
,
Chem. Sci.
12
,
9704
(
2021
).
34.
X.
Gao
,
J.
Ren
,
A.
Eisfeld
, and
Z.
Shuai
,
Phys. Rev. A
105
,
L030202
(
2022
).
35.
P.-P.
Zhang
and
A.
Eisfeld
,
J. Phys. Chem. Lett.
7
,
4488
4494
(
2016
).
36.
L.
Chen
,
D. I. G.
Bennett
, and
A.
Eisfeld
,
J. Chem. Phys.
156
,
124109
(
2022
).
37.
J.
Roden
,
W. T.
Strunz
,
K. B.
Whaley
, and
A.
Eisfeld
,
J. Chem. Phys.
137
,
204110
(
2012
).
38.
M.
Combescure
and
D.
Robert
,
Coherent States and Applications in Mathematical Physics
, 2nd ed. (
Springer
,
2012
).
39.
H.-P.
Breuer
,
B.
Kappler
, and
F.
Petruccione
,
J. Phys. A: Math. Gen.
31
,
L147
L151
(
1998
).
40.
G. W.
Semenoff
and
H.
Umezawa
,
Nucl. Phys. B
220
,
196
(
1983
).
41.
G.
Ritschel
,
D.
Suess
,
S.
Möbius
,
W. T.
Strunz
, and
A.
Eisfeld
,
J. Chem. Phys.
142
,
034115
(
2015
).
42.
R.
Hartmann
and
W. T.
Strunz
,
J. Chem. Theory Comput.
13
,
5834
5845
(
2017
).
43.
D. W.
Schönleber
,
A.
Croy
, and
A.
Eisfeld
,
Phys. Rev. A
91
,
052108
(
2015
).
44.
P.-P.
Zhang
,
C. D. B.
Bentley
, and
A.
Eisfeld
,
J. Chem. Phys.
148
,
134103
(
2018
).
45.
Y.
Ke
and
Y.
Zhao
,
J. Chem. Phys.
149
,
014104
(
2018
).
46.
J. T.
Stockburger
and
H.
Grabert
,
Chem. Phys.
268
,
249
(
2001
).
47.
J. T.
Stockburger
and
H.
Grabert
,
Phys. Rev. Lett.
88
,
170407
(
2002
).
48.
Z.
Zhang
,
P. H.
Lambrev
,
K. L.
Wells
,
G.
Garab
, and
H.-S.
Tan
,
Nat. Commun.
6
,
7914
(
2015
).
49.
A. C.
Davison
and
D. V.
Hinkley
,
Bootstrap Methods and Their Application
, Cambridge Series in Statistical and Probabilistic Mathematics (
Cambridge University Press
,
1997
).

Supplementary Material