# Process tensor approaches to modeling two-dimensional spectroscopy

Roosmarijn de Wit, Jonathan Keeling, Brendon Lovett, Alex W. Chin

► To cite this version:

Roosmarijn de Wit, Jonathan Keeling, Brendon Lovett, Alex W. Chin. Process tensor approaches to modeling two-dimensional spectroscopy. Physical Review Research, 2025, 7 (1), pp.013209. 10.1103/PhysRevResearch.7.013209. hal-05406570

HAL Id: hal-05406570

https://hal.sorbonne-universite.fr/hal-05406570v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Process tensor approaches to modeling two-dimensional spectroscopy

Roosmarijn de Wit,¹ Jonathan Keeling,¹ Brendon W. Lovett,¹ and Alex W. Chin²

¹SUPA, School of Physics and Astronomy, University of St Andrews, St Andrews, KY16 9SS, United Kingdom

²Sorbonne Université, CNRS, Institut des NanoSciences de Paris, 4 place Jussieu, 75005 Paris, France

(Dated: December 9, 2025)

Problems in the field of open quantum systems often involve an environment that strongly influences the dynamics of excited states. Here we present a numerical method to model optical spectra of non-Markovian open quantum systems. The method employs a process tensor framework to efficiently compute multi-time correlations in a numerically exact way. To demonstrate the efficacy of our method, we compare 2D electronic spectroscopy simulations produced through our method to Markovian master equation simulations in three different system-bath coupling regimes.

# I. INTRODUCTION

Two-dimensional electronic spectroscopy (2DES) is a valuable tool for probing photophysical processes in light-harvesting systems [1–11]. Generally, these systems do not exist in isolation, and measured quantities are significantly impacted by a vibrational environment. 2DES has yielded key insights into the role of such environment-mediated processes [5–9, 11]. However, extracting the underlying mechanism remains challenging. Theoretically, realistic environments are often non-Markovian, meaning that time-local (Markovian) equations of motion are not sufficient [12, 13]. Experimentally, the complexity of measured signals can make it difficult to accurately interpret observed spectral features. For example, the as-yet-unclear origin of long-lived coherences in biological light absorbing pigments highlights the necessity of accurate theoretical models to explain experimental observations [6–9, 11, 14–16].

Here we present an efficient tensor network method to simulate optical spectra of non-Markovian open quantum systems. To demonstrate its capabilities, we introduce a model describing the essential features of a broad class of molecular chromophores. It consists of a three-level electronic system coupled via a vibrational bath, enabling the intramolecular transfer of energy between system states. We first compare numerical results to master equations in weak and polaronic system-bath coupling regimes to test the validity of our method. Then, we highlight key differences between the 2D spectra obtained via these methods in an intermediate coupling regime. Namely, we find that there are significant differences in peak positions at high temperatures, related to the Lamb shift. At low temperatures, we furthermore find that a master equation underestimates the dephasing time of the optical response, as has been observed in previous work [14].

Mathematically, a 2DES signal is expressed as a sum of four-time correlation functions that encode all possible light-matter interactions. A common approach to calculating multi-time correlation functions is to use the quantum regression theorem, which is based on the Born-Markov approximation [12, 17]. In many realistic scenarios however, a Markovian picture is not sufficient. Even if it accurately represents the reduced dynamics of the sys-

tem, a Markovian description does not necessarily capture multi-time correlations correctly [18]. Some light-harvesting complexes are furthermore known to operate in the intermediate coupling regime, where electronic couplings within the system are comparable in strength to system-bath interactions and a Markovian approximation is not justified. [19–21].

Simulating non-Markovian open quantum systems is a challenging task. Nonetheless, numerous numerical methods have been developed to achieve this [6, 13, 14, 16, 22–41]. Most relevant to this work, one group of such methods uses a process tensor (PT) formalism to capture the influence of the environment on the system. Importantly, the process tensor is cast into a matrix product operator format (PT-MPO), such that only the physically most relevant part of the Hilbert space is efficiently represented [22–31, 41]. Moreover, since PT-MPO methods rely on constructing a reduced density matrix description of the system, they are particularly well-tailored to calculating temporal correlations.

The remainder of this paper is organized as follows. Section II introduces the model, Section III specifies the multi-time correlation functions relevant to our results, and Section IV provides an overview of the PT-MPO computation method. The results are presented in Sections V and VI; Section V compares spectra obtained via PT-MPO methods versus Markovian master equations in different system-bath coupling regimes, while Section VI further analyses the PT-MPO results in an intermediate coupling regime.

# II. THE MODEL

As illustrated in Fig. 1(a), our model consists of a three-level system with excited states linearly coupled to a bosonic bath, $\hat{H} = \hat{H}_S + \hat{H}_B + \hat{H}_I$ with:

$$ \hat{H}_S = (\epsilon + \lambda)(|1\rangle\langle1| + |2\rangle\langle2|) + \Omega(|1\rangle\langle2| + \text{H.c.}), \quad (1) $$

$$ \hat{H}_B = \sum_k \omega_k b_k^\dagger b_k, \quad (2) $$

$$ \hat{H}_I = (|1\rangle\langle1| - |2\rangle\langle2|) \sum_k (g_k b_k + g_k^* b_k^\dagger). \quad (3) $$
Here, the three terms correspond to the system ($\hat{H}_S$), bath ($\hat{H}_B$) and system-bath interaction ($\hat{H}_I$). The energy of the excited states ($|1\rangle, |2\rangle$) is given by the bare energy $\epsilon$ plus the reorganization energy of the bath $\lambda$ (defined below), and electronic coupling $\Omega$. Physically, $|1\rangle$ and $|2\rangle$ could describe two excited states of a chromophore with ground state $|0\rangle$. These excited states could be intramolecular in nature, as in the $S_1, S_2$ states found in many biological pigments, organic dye molecules and nanoparticles [42–44], or the lowest single-particle excitations of coupled chromophores, such as excitonically coupled H-dimers, or charge-transfer complexes [45, 46]. For clarity in later discussions, we imagine our system to be an excitonic dimer, so that the coupling $\Omega$ can be thought of as creating electronic eigenstates (i.e. excitons; $|\pm\rangle = \frac{1}{\sqrt{2}}(|1\rangle \pm |2\rangle)$) that are delocalized across the internal chromophore monomers with energies $E_{\pm} = \epsilon + \lambda \pm \Omega$.

Chromophores in biological systems absorb visible light, which sets the overall energy scale $\epsilon$ of the model. However, they often display photophysics (dynamics) on a picosecond time scale, which we set by choosing $\Omega = 0.2 - 2 \text{ ps}^{-1}$ [7, 9]. Since $\epsilon$ is relatively large, we expect that any effects due to coupling between ground and excited states are primarily associated with absorption of from external light and the emission of the output signal. As such, we do not include any coupling to the ground state in the Hamiltonian in Eq. (1), but instead consider the response to externally applied pulses as discussed below. We will display all results in a rotating frame ($\epsilon = 0$), such that the energy scale set by $\epsilon$ is effectively removed from the response.

The bath is described by vibrational modes with frequency $\omega_k$ and creation (annihilation) operators $b_k^\dagger$ ($b_k$). The system-bath coupling constants $g_k$ are characterized by the spectral density $J(\omega) = \sum_k |g_k|^2 \delta(\omega - \omega_k)$. We set $J(\omega)$ to be Ohmic, as is common in the modeling and interpretation of molecular spectroscopies [1, 2]:

$$J(\omega) = 2\alpha\omega e^{-\frac{\omega}{\omega_c}}, \quad (4)$$

where $\alpha$ is a dimensionless parameter for the system-bath coupling strength and $\omega_c = 3.04 \text{ ps}^{-1}$ represents the bath cut-off frequency. The reorganization energy is then given by $\lambda = \int_0^\infty d\omega \frac{1}{\omega} J(\omega) = 2\alpha\omega_c$. This serves as a characteristic energy scale of energy shifts due to the bath, discussed further below. We note here that realistic spectral densities for biological chromophores generally contain a series of high-frequency underdamped phonon modes, in addition to an Ohmic background [16, 47, 48]. For our PT-MPO method, including such underdamped modes would increase the memory time of the bath and therefore the computational cost of the simulations. Methods to reduce the computational cost associated with highly structured baths do exist however, for example by coarse-graining the environment [49]. As the main purpose of this work is to demonstrate the method, we have chosen to keep the spectral density generic.

<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Sketch of our model, consisting of a three-level system coupled to a bosonic bath. Note that the energy splittings $\epsilon$ and $\Omega$ are not to scale; $\epsilon$ is much larger than $\Omega$, as discussed in the main text. (b) Tensor network of ten time steps for calculating a four-time correlation function with a PT-MPO. The purple circle denotes the initial system state $\rho_0^S$, the orange diamonds the system operators $\hat{V}(\tau_i)$ and the blue ovals the system propagators $K$. The gray half-circle symbolizes a partial trace over the system subspace.</figcaption></figure>

We finally define our initial state $\rho_0$ to be a product state of the system ground state $|0\rangle$ and the thermal state of the environment at temperature $T$.

### III. 2D ELECTRONIC SPECTROSCOPY

2DES involves probing a quantum system with three short laser pulses while systematically varying the time delays between each pulse. For a comprehensive review on 2DES, we refer to Refs. [2, 3]. Here, we will consider the measured third-order optical response, which consists of a series of four-time correlation functions. In order to predict the spectra resulting from this model, we will simulate laser pulses that induce transitions between the ground state and excited state $|2\rangle$, corresponding to the dipole operator $\hat{V}_2 = |0\rangle\langle 2| + \text{H.c.}$. We will assume we are in the semi-impulsive limit, such that we can neglect the temporal width of the laser pulses [1, 2, 50]. We will furthermore make the rotating-wave approximation and assume that the laser pulses applied at times $\tau_i$ are time-ordered [1, 50]. Consequently, only the following four multi-time correlation functions contribute to the 2DES signal:

$$ \begin{aligned} R_1 &= \text{Tr} [\hat{V}_2(\tau_4)\hat{V}_2(\tau_1)\rho_0\hat{V}_2(\tau_2)\hat{V}_2(\tau_3)] \\ R_2 &= \text{Tr} [\hat{V}_2(\tau_4)\hat{V}_2(\tau_2)\rho_0\hat{V}_2(\tau_1)\hat{V}_2(\tau_3)] \\ R_3 &= \text{Tr} [\hat{V}_2(\tau_4)\hat{V}_2(\tau_3)\rho_0\hat{V}_2(\tau_1)\hat{V}_2(\tau_2)] \\ R_4 &= \text{Tr} [\hat{V}_2(\tau_4)\hat{V}_2(\tau_3)\hat{V}_2(\tau_2)\hat{V}_2(\tau_1)\rho_0], \end{aligned} \quad (5) $$

which are commonly grouped together as rephasing ($R_2, R_3$) and non-rephasing ($R_1, R_4$) pathways, accord-
ing to the experimental phase matching conditions [1–3]. The real parts of these pathways in the frequency domain sum up to give the total 2D spectrum; this is a function of an excitation ($\omega_{excite}$) and detection frequency ($\omega_{detect}$) corresponding to Fourier transforms of the first ($\tau_2 - \tau_1$) and third time delay ($\tau_4 - \tau_3$) between pulses. Note that in this work, we will only consider 2D spectra for which the so-called waiting time $\tau_3 - \tau_2$ is zero. For the linear response, there is only a single pathway: $\mathrm{Tr}[\hat{V}_2(\tau_2)\hat{V}_2(\tau_1)\rho_0]$, which is similarly Fourier transformed with respect to $\tau_2 - \tau_1$ to obtain a linear absorption spectrum as a function of frequency $\omega$.

## IV. SYSTEM RESPONSE FUNCTIONS AND PT-MPOS

In order to calculate the multi-time correlation functions in Eq. (5), we will employ the Process Tensor - Time Evolving MPO method (PT-TEMPO) to acquire a PT-MPO that captures any possible non-Markovian effects on the system dynamics [23, 25, 41, 51–53]. In the following section, we will present our method for calculating multi-time correlation functions with PT-MPOs, and give a brief description of the PT-TEMPO method. For the interested reader, more detailed descriptions of process tensors and PT-TEMPO can be found in Ref. [41] (a broad overview) and additionally Refs. [24, 25, 52]. The Python code for constructing process tensors with PT-TEMPO is available in the open-source package OQuPy [54], along with the functionality for calculating multi-time correlations introduced in this work.

### A. Calculating multi-time correlation functions

Let us consider time-ordered multi-time correlations of a general form:

$$R = \left\langle \prod_{p=1}^{P} \hat{V}_p(\tau_p) \right\rangle = \left\langle \prod_{p=1}^{P} \hat{V}_p(M_p \delta t) \right\rangle \quad (6)$$

such that P operators $\hat{V}_p$ are applied at times $\tau_p = M_p \delta t$, where $\delta t$ denotes the time step. A formal approach to calculating R is to consider the time evolution of the total (i.e. system and bath) density matrix $\rho_0$. First, $\rho_0$ is propagated up to $M_1$ time steps according to the formal solution of the von Neumann equation:

$$\rho(M_1 \delta t) = e^{\mathcal{L} M_1 \delta t} \rho_0, \quad (7)$$

with the total Liouvillian $\mathcal{L} = -i[\hat{H}, \cdot]$. Then, after application of the first operator $\hat{V}_1$, the resulting density matrix is propagated for $M_2 - M_1$ time steps, and so on, up to the application of the final operator $\hat{V}_P$. Thus, we

can write Eq. (6) in terms of the time propagators:

$$R = \mathrm{Tr} \left[ \prod_{p=1}^{P} \left( \hat{V}_p^{L,R} \left[ e^{\mathcal{L} \delta t (M_p - M_{p-1})} \right] \right) \rho_0 \right], \quad (8)$$

where $M_0 = 0$. Here the superscript L, R signifies that the super-operator $\hat{V}$ can act to the left ($\hat{V}^L = \hat{V} \cdot$) or to the right ($\hat{V}^R = \cdot \hat{V}$) of the density matrix.

However, applying this approach directly would require direct evolution of the full density matrix $\rho_0$, which is not practical. Instead, we calculate Eq. (8) by constructing a tensor network in which the influence of the bath on the system is encoded in a PT-MPO. The PT is a multi-linear map from operations performed on a system, at a sequence of time steps, to its final state. Crucially, since the PT is completely independent of the system Hamiltonian $\hat{H}_S$, it can be calculated before specifying any system parameters, the initial system state $\rho_0^S$ or the operators $\hat{V}$ in Eq. (8). To separate system and bath propagation, we perform a second-order Suzuki-Trotter splitting [55]:

$$[e^{\mathcal{L} \delta t}]^M = \left[ e^{\mathcal{L}_S \frac{\delta t}{2}} e^{\mathcal{L}_E \delta t} e^{\mathcal{L}_S \frac{\delta t}{2}} \right]^M + \mathcal{O}(\delta t^3), \quad (9)$$

such that $\mathcal{L}_S \cdot = -i[\hat{H}_S, \cdot]$ and $\mathcal{L}_E \cdot = -i[\hat{H}_E + \hat{H}_I, \cdot]$. Once the PT-MPO is constructed (see the next section), multi-time correlation functions are calculated by combining the PT-MPO with the operators $\hat{V}_p$ at times $M_p \delta t$ and system propagators $K = e^{\mathcal{L}_S \frac{\delta t}{2}}$.

Fig. 1(b) shows a tensor network diagram of a PT-MPO for ten time steps (red squares) and illustrates how multi-time correlations are calculated. The tensor network is expressed in Liouville space, such that density matrices are represented by vectors and superoperators by matrices. Diagrammatically, the initial system state $[\rho_0^S]^j$ therefore corresponds to a tensor with a single leg (purple circle), where the index j runs from 1 to the squared dimension of the system Hilbert space. Similarly, the dipole operators $[V]_j^{j'}$ (orange diamonds) and system propagators $[K]_j^{j'}$ (two per single time step, blue ovals) have two legs. Thus, multi-time correlation functions are computed by applying the system propagators and dipole operators as a set of interventions at the relevant time steps and subsequently tracing over the bonds in the network to obtain R [41, 53, 56]. Since the final dipole operator $\hat{V}_P$ is applied at the end of the chain, we only have to run the simulation once to calculate R over a range of the final times $M_P \delta t$. When varying some earlier time arguments $M_p \delta t$, e.g. over a range of i time steps, the simulation is repeated i times, moving the position of $\hat{V}_p$ in the tensor network (Fig. 1(b)) with each repetition. For further details on the required computational resources, see Appendix A and Ref. [41].
## B. PT-MPO construction

The method described above is general, and there exist multiple algorithms to construct the PT-MPO [22–31, 41]. In this work, we employ the PT-TEMPO method to construct a PT-MPO, which we will summarize here [24, 25, 41, 54]. For readers already familiar with PT-MPOs, this section can be skipped.

The general TEMPO framework expresses the (non-Markovian) impact of the bath on the system in terms of a discretized Feynman-Vernon influence functional [52, 57–59]. In PT-TEMPO, the parts of the TEMPO network corresponding to the influence functional are contracted into a PT-MPO that does not depend on the system propagators. To illustrate how the tensor network is constructed, let us consider the propagation of the system density matrix $\rho^S$ for M time steps, such that the total propagation consists of M propagators over short time steps $\delta t: [e^{\mathcal{L}\delta t}]^M$. We start with the Suzuki-Trotter splitting in Eq. (9). Then, resolutions of identity are added between each system and bath propagator. By tracing over the environment, we obtain a path sum over system states composed of the discretized influence functional and system propagators:

$$ \rho_{jM}^S (M\delta t) = \sum_{\substack{j_0, \dots, j_{M-1} \\ j'_0, \dots, j'_{M-1}}} \left[ \prod_{m=0}^{M-1} K_m(j_{m+1}, j'_m) \times \left( \prod_{k=0}^{m} I_k(j'_m, j'_{m-k}) \right) K_m(j'_m, j_m) \right] \rho_{j_0}^S. \quad (10) $$

In this expression, the system propagators are given by $K(j, j') = [\exp(\mathcal{L}_S \frac{\delta t}{2})]_{j, j'}$. Because each propagator evolves $\rho_j^S$ by half a time step, we require two indices: $j_m$ runs over each full step, while $j'_m$ connects the pairs of propagators within a step. The influence functions $I_k(j, j')$ capture the effect of the bath on the system, and connect system states separated by k time steps. For the exact form of $I_k$, see Ref. [52].

Fig. 2(a) illustrates how Eq. (10) can be constructed as a tensor network. Per time step, we apply two system propagators (blue ovals, each evolving by half a time step), represented by the rank-2 tensors $K_j^{j'}$. The influence functions are incorporated as the following bath tensors:

$$ [b_k]_{\alpha', j'}^{\alpha, j} = \delta_{j'}^{j} \delta_{\alpha'}^{\alpha} I_k(\alpha, j), \quad (11) $$

which are drawn as red squares (labeled with index k) in Fig. 2(a). At the left and top edges of the network, we instead require tensors that lack the $\alpha'$ or $j'$ leg. These are obtained from Eq. (11) by respectively omitting the $\delta_{\alpha'}^{\alpha}$ or $\delta_{j'}^{j}$ Kronecker deltas. To illustrate how Eq. (10) maps to a tensor network, consider the propagation of $\rho^S$

<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) The TEMPO tensor network. The bath tensors $b_k$ are represented by red squares, labeled with index k. (b) The network for a single time step. (c) The bath tensors of the TEMPO network are horizontally contracted into a process tensor.</figcaption></figure>

for a single time step ($M=1$):

$$ \rho_{j_1}^S = \sum_{j_0, j'_0} K_0(j_1, j'_0) I_0(j'_0, j'_0) K_0(j'_0, j_0) \rho_0^S, \quad (12) $$

which can be written as the tensor contraction

$$ [\rho^S]^{j_1} = \sum_{j_0, j'_0, \alpha} [K_0]_{\alpha}^{j_1} [b_0]_{j'_0}^{\alpha} [K_0]_{j_0}^{j'_0} [\rho^S]^{j_0}, \quad (13) $$

shown in Fig. 2(b).

The number of bath tensors—and thus the bath memory—builds up with each time step, such that at time $m\delta t$, m bath tensors are added to the tensor network (Fig. 2(a)). We can introduce a maximal memory time of the bath, set by the parameter $\Delta K_{max}$, which ensures that the number of bath tensors added to the network per time step stops increasing after $\Delta K_{max}$ time steps.

With PT-TEMPO, the bath tensors of the TEMPO network shown in Fig. 2(c) are contracted into a PT-MPO using a sequence of standard tensor contraction and compression techniques [22]. This results in an array of rank-4 tensors, where the unconnected legs have the physical dimension of the system Hilbert space squared. The dimensions of the bonds that connect the tensors quantify the degree of non-Markovianity of the bath [60]. Through the truncation of singular values, the compression steps ensure that the tensor bond dimensions (and thus the overall complexity of the network) are kept to a minimum. PT-TEMPO is a numerically exact algorithm, meaning that all errors are controlled by a set of numerical convergence parameters (Appendix B). The truncation threshold is determined by the parameter $\epsilon_{rel}$,
which sets the upper bound to the singular values to be discarded, relative to the largest singular value in the tensor.

Overall, the applicability of the method largely depends on the size of the bond dimension required to accurately simulate the environment. The efficiency furthermore strongly depends on the size of the system, making the algorithm primarily useful for simulating few-level systems. For a wide variety of physically relevant non-Markovian problems, however—including the model discussed in this work—it is possible to efficiently perform the computations with reasonable resources [41].

## V. MASTER EQUATION COMPARISON

To test the validity of our method, we compare spectra obtained with PT-TEMPO against a Markovian weak-coupling master equation (WCME) and polaron master equation (PME). To do this, we present the results in two different system-bath coupling regimes. We characterize these regimes by comparing the frequency scales of the system (set by the electronic coupling $\Omega$) against the frequency scales of the bath (set by the reorganization energy $\lambda$). We first consider a regime where bath frequency scales are much smaller than those of the system ($\lambda \ll \Omega$). Here the WCME is expected to give accurate results [12]. Conversely, the PME is valid when the system-bath coupling strength is much stronger than the electronic coupling ($\lambda \gg \Omega$) [61–64]. As with PT-TEMPO, the master equation spectra were simulated with OQuPy [54], which has a Lindblad master equation solver implemented that works in tandem with the multi-time correlations module.

### A. Weak-coupling master equation

Starting with the WCME, we obtain the following master equation in Lindblad form using standard techniques [12]:

$$ \frac{d}{dt} \rho_S = -i[\hat{H}'_S, \rho_S] + \sum_n^2 \gamma_n \left( \hat{L}_n \rho_S \hat{L}_n^\dagger - \frac{1}{2} \{ \hat{L}_n^\dagger \hat{L}_n, \rho_S \} \right), \quad (14) $$

where $\hat{L}_1 = |-\rangle\langle +|$ and $\hat{L}_2 = \hat{L}_1^\dagger$. The rates $\gamma_n$ are given by $\gamma_1 = 2\pi J(2\Omega)(N(2\Omega) + 1)$ and $\gamma_2 = 2\pi J(2\Omega)N(2\Omega)$, where $N(\omega)$ is the Bose-Einstein occupation number and $J(\omega)$ is the spectral density in Eq. (4). The electronic energies are furthermore modified by an energy shift due to the bath—a form of Lamb shift—so that

$$ \hat{H}'_S = (E_- + S(-2\Omega)) |-\rangle\langle -| + (E_+ + S(2\Omega)) |+\rangle\langle +|, \quad (15) $$

in the eigenbasis $|\pm\rangle$ with $E_{\pm} = \epsilon + \lambda \pm \Omega$, and

$$ S(\nu) = \mathcal{P} \int_{-\infty}^{\infty} d\omega \frac{J(\omega)(N(\omega) + 1) + J(-\omega)N(-\omega)}{\nu - \omega}. \quad (16) $$

Here $\mathcal{P}$ denotes the Cauchy principal value and we define $J(\omega) = 0$ for $\omega < 0$. Note that at zero temperature one has $|S(0)| = \lambda$, the reorganization energy. $S(+2\Omega)$ corresponds to a positive shift to $E_+$, while $S(-2\Omega)$ corresponds to a negative shift to $E_-$. The total energy splitting between the excited eigenstates therefore becomes:

$$ \delta E = 2\Omega + S(+2\Omega) - S(-2\Omega) \quad (17) $$

$$ = 2\Omega \left( 1 + 2 \int_{0}^{\infty} d\omega \frac{J(\omega) \coth\left(\frac{\omega}{2k_{B}T}\right)}{(2\Omega)^2 - \omega^2} \right). \quad (18) $$

As illustrated in Figs. 3(a-c), the PT-TEMPO simulations of linear absorption and 2D spectra agree well with the WCME in a weak-coupling regime. Here we have set $\Omega = 2 \text{ ps}^{-1}$ and $\alpha = 0.005$, such that $\lambda = 0.03 \text{ ps}^{-1}$. In the 2D spectra, the two peaks along the diagonal of the spectra reflect the transition frequencies of the model, which coincide with the peak frequencies in the linear absorption spectrum. For $\alpha = 0.005$, the Lamb shift, Eq. (16), is small compared to $\epsilon$ and $\Omega$, such that the splitting between peaks is close to the bare splitting $2\Omega$. The 2D spectra furthermore contain two cross-peaks that correlate the transition frequencies and are a signature of the coupling between the electronic excited states.

For some problems with weakly-structured environments, i.e. those with non-constant $J(\omega)$, one can improve on the Lindblad master equation by restoring non-secular terms from Redfield theory [65, 66]. For our particular choice of dipole operators and initial conditions, the non-secular terms in the master equation do not affect the outcome. To show this, consider the individual components of the reduced system density matrix

$$ \rho_S = \begin{pmatrix} \rho_{00} & \rho_{0-} & \rho_{0+} \\ \rho_{-0} & \rho_{--} & \rho_{-+} \\ \rho_{+0} & \rho_{+-} & \rho_{++} \end{pmatrix}, \quad (19) $$

expressed in terms of system eigenstates. The non-secular terms in the master equation would have the form:

$$ \sum_n^2 \gamma_n (\hat{L}_n \rho_S \hat{L}_n + \hat{L}_n^\dagger \rho_S \hat{L}_n^\dagger), \quad (20) $$

where $\hat{L}_1 = |-\rangle\langle +|$, $\hat{L}_2 = \hat{L}_1^\dagger$ and $\gamma_n$ are the transition rates stated in the main text. Note here that $\hat{L}_n \hat{L}_n^\dagger = \hat{L}_n^\dagger \hat{L}_n^\dagger = 0$, and terms involving these are therefore not included in Eq. (20). For the individual components of $\rho_S$ in Eq. (19), we find that only $\rho_{-+}$ and $\rho_{++}$ would be modified by these non-secular terms. The dipole operator $\hat{V}_2 = |0\rangle\langle 2| + \text{H.c.}$, on the other hand, only creates $\rho_{0\pm}$ and $\rho_{\pm 0}$ components when applied to the initial state $|0\rangle\langle 0|$. Thus, for the absorption spectra presented in this work,
<figure><img src="image_4.png" /><figcaption>FIG. 3. Comparison of linear absorption and 2D spectra in (a-d) a weak system-bath coupling regime, (e-h) a stronger coupling (polaronic) regime and (i-l) an intermediate coupling regime, calculated using PT-TEMPO (blue solid in linear spectra), a weak-coupling master equation (WCME; orange dashed in in linear spectra) or a polaron master equation (PME; green dot-dashed in linear spectra). All spectra were calculated at $T = 13 \text{ ps}^{-1}$ (100 K), using dipole operator $\hat{V}_2$. For (j) and (k), the blue dashed lines mark the diagonal peak maxima.</figcaption></figure>

there would be no difference between secular and non-secular master equations. For the 2D spectra, $\rho_{+-}$ and $\rho_{-+}$ components would be created after application of the second dipole operator at time $\tau_2$, Eq. (5). Since in this work we are only considering multi-time correlations for which $\tau_2 = \tau_3$ (the third operator is applied immediately after the second), there would be no difference for the 2D spectra presented here either.

## B. Polaron master equation

To derive a master equation for our three-level system, Eq. (1), in a regime where bath frequency scales are much

larger than characteristic system frequencies, we first apply a polaron transformation to the Hamiltonian [61–64]:

$$ \tilde{H} = e^G \hat{H} e^{-G}; \quad G = \sum_k \frac{g_k}{\omega_k} (b_k^\dagger - b_k) (|1\rangle\langle 1| - |2\rangle\langle 2|), \qquad (21) $$

to obtain:

$$ \tilde{H}_S = \epsilon(|1\rangle\langle 1| + |2\rangle\langle 2|) + \Omega'(|1\rangle\langle 2| + \text{H.c.}) $$

$$ \tilde{H}_B = \hat{H}_B $$

$$ \tilde{H}_I = \Omega D_c |1\rangle\langle 2| + \text{H.c.} \qquad (22) $$

In the polaron frame, the interaction Hamiltonian $\hat{H}_I$ contains the coupling operator $D_c$, which depends on the
displacement operator

$$ \mathcal{D} = \exp\left\{\sum_k \frac{2g_k}{\omega_k}(b_k^\dagger - b_k)\right\}, \quad (23) $$

such that $\mathcal{D}_c$ is given by:

$$ \mathcal{D}_c = \mathcal{D} - \langle \mathcal{D} \rangle. \quad (24) $$

Here, the expectation value of $\mathcal{D}$

$$ \langle \mathcal{D} \rangle = \exp \left\{ -2 \sum_k \frac{g_k^2}{\omega_k^2} \coth \left( \frac{\beta \omega_k}{2} \right) \right\} \quad (25) $$

has been subtracted from $\mathcal{D}$ and is instead added to $\hat{H}_S$, which now contains the renormalized electronic coupling:

$$ \Omega' = \Omega \langle \mathcal{D} \rangle. \quad (26) $$

Note here that for our choice of an Ohmic spectral density in Eq. (4), $\langle \mathcal{D} \rangle$ (and therefore $\Omega'$) are infinitesimally small. Following the same standard techniques as for the WCME, we can now perturbatively expand $\hat{H}_I$ to find the polaron master equation:

$$ \frac{d}{dt} \tilde{\rho}_S = -i \left[ \tilde{H}'_S, \tilde{\rho}_S \right] + \sum_n^3 \Lambda_n \left( \hat{L}_n \tilde{\rho}_S \hat{L}_n^\dagger - \frac{1}{2} \left\{ \hat{L}_n^\dagger \hat{L}_n, \tilde{\rho}_S \right\} \right), \quad (27) $$

with Lindblad operators $\hat{L}_1 = |-\rangle\langle +| = \hat{L}_2^\dagger$ and $\hat{L}_3 = |-\rangle\langle -| - |+\rangle\langle +|$. The transition rates $\Lambda_n$ depend on the four bath correlation functions $\langle \mathcal{D}_c(t) \mathcal{D}_c(0) \rangle$, $\langle \mathcal{D}_c^\dagger(t) \mathcal{D}_c^\dagger(0) \rangle$, $\langle \mathcal{D}_c^\dagger(t) \mathcal{D}_c(0) \rangle$ and $\langle \mathcal{D}_c(t) \mathcal{D}_c^\dagger(0) \rangle$. For example:

$$ \langle \mathcal{D}_c(t) \mathcal{D}_c^\dagger(0) \rangle = \langle e^{i\hat{P}(t)} e^{-i\hat{P}(0)} \rangle - \langle \mathcal{D} \rangle^2, \quad (28) $$

where $\hat{P}(t) = -i \sum_k \frac{2g_k}{\omega_k} (b_k^\dagger(t) - b_k(t))$. To solve Eq. (28), we can add the exponents together using the Baker-Campbell-Hausdorff formula:

$$ \langle \mathcal{D}_c(t) \mathcal{D}_c^\dagger(0) \rangle = \left\langle e^{i\hat{P}(t)-i\hat{P}(0)+\frac{1}{2}[i\hat{P}(t), -i\hat{P}(0)]} \right\rangle - \langle \mathcal{D} \rangle^2. \quad (29) $$

Then, since the problem is Gaussian in the bosonic operators, we can shift the expectation to the exponent. In our case $\langle \hat{P} \rangle = 0$, so only the second order term in $\hat{P}$ survives:

$$ \langle e^{\hat{P}} \rangle = \exp \left[ \frac{1}{2} \langle \hat{P}^2 \rangle \right]. \quad (30) $$

This leaves us with the following expression:

$$ \langle \mathcal{D}_c(t) \mathcal{D}_c^\dagger(0) \rangle = \exp \left[ \frac{1}{2} \left\langle \hat{P}(0)^2 \right\rangle + \left\langle \hat{P}(t) \hat{P}(0) \right\rangle \right] - \langle \mathcal{D} \rangle^2, \quad (31) $$

which we can further simplify by defining the phonon propagator

$$ \phi(t) = \int_0^\infty d\omega \frac{J(\omega)}{\omega^2} \left( \cos(\omega t) \coth\left(\frac{\omega}{2k_B T}\right) - i \sin(\omega t) \right). \quad (32) $$

Note here that for an Ohmic spectral density, the integrand of $\phi(0)$ diverges as $\omega \to 0$, meaning that $e^{-4\phi(0)} = 0$. Applying this limit to all four bath correlation functions leaves us with:

$$ \begin{aligned} \langle \mathcal{D}_c(t) \mathcal{D}_c^\dagger(0) \rangle &= \langle \mathcal{D}_c^\dagger(t) \mathcal{D}_c(0) \rangle = e^{4(\phi(t)-\phi(0))}; \\ \langle \mathcal{D}_c^\dagger(t) \mathcal{D}_c^\dagger(0) \rangle &= \langle \mathcal{D}_c(t) \mathcal{D}_c(0) \rangle = 0. \end{aligned} \quad (33) $$

Finally, taking the time integral over these correlations and inserting the coupling $\Omega$ gives us the transition rates $\Lambda_n$ in the master equation in Eq. (27):

$$ \begin{aligned} \Lambda_1 &= \Lambda_2 = \int_0^\infty dt \, \Omega^2 e^{4(\phi(t)-\phi(0))}, \\ \Lambda_3 &= 2\Lambda_1, \end{aligned} \quad (34) $$

which can be solved numerically. As for the WCME, the imaginary parts of these integrals correspond to the Lamb shift, and are incorporated into $\hat{H}'_S$. For the PME however, we find that the Lamb shift is about two orders of magnitude smaller than system frequency scales, and thus has a limited impact on the peak locations in Fig. 3.

As shown in Fig. 3(e-g), the PME results agree well with PT-TEMPO when $\lambda \gg \Omega$. In particular, we have set $\alpha = 0.4$, $\lambda = 2.4$ ps$^{-1}$ and $\Omega = 0.2$ ps$^{-1}$. Unlike the WCME spectrum (Fig. 3(g)) in which the peak is homogeneously broadened in all directions, the PT-TEMPO and PME results capture correlations between excitations at $t_1$ and $t_3$, leading to an inhomogeneously broadened peak stretched along the diagonal axis. For the PME, the transition rates $\Lambda_n$ in Eq. (34) are negligibly small. Rather, the peak broadening in Fig. 3(h) comes from the polaron transformation in Eq. (21), which encodes optical dephasing between the ground and excited states.

Because $\hat{H}_I$ is proportional to $\Omega$ (see Eq. (22)), larger values of $\Omega$ would not give correct results. This is illustrated in Fig. 3(a)(d), which shows the 2D and linear spectrum obtained with the PME in the weak-coupling regime ($\lambda = 0.03$ ps$^{-1}$, $\Omega = 2.0$ ps$^{-1}$). Since $\Omega'$ is infinitesimally small, only a single peak (and therefore no cross-peaks) can be resolved. Additionally, because the rates in Eq. (34) are proportional to $\Omega^2$, this peak is strongly broadened compared to the PT-TEMPO and WCME solutions.

## VI. DIFFERENCES IN THE INTERMEDIATE COUPLING REGIME

We will now shift our focus to the intermediate coupling regime found in photosynthetic complexes [19–21], where neither the weak-coupling nor polaron master
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Splitting between diagonal peaks in 2D spectra as a function of temperature, as predicted by PT-TEMPO and the WCME (Lamb shift). Data is shown for system-bath couplings $\alpha = 0.1$ (green solid, blue dots) and $\alpha = 0.02$ (red dashes and orange cross markers). Note that due to increasing computational resources with decreasing $T$, data points for PT-TEMPO below $3.9 \text{ ps}^{-1}$ (30 K) are not shown. (b-d) Diagonal slices through the 2D spectra for $\alpha = 0.1$, corresponding to the data points $T = 3.9 \text{ ps}^{-1}$ (30 K), $T = 6.5 \text{ ps}^{-1}$ (50 K) and $T = 25 \text{ ps}^{-1}$ (190 K) in panel (a). For all: $\Omega = 2 \text{ ps}^{-1}$.</figcaption></figure>

equation are expected to give accurate results. For this purpose, we will set $\Omega = 2 \text{ ps}^{-1}$ and $\alpha = 0.1$, such that $\lambda = 0.6 \text{ ps}^{-1}$. As in the weak-coupling regime, the 2D spectrum obtained with PT-TEMPO for these parameters, Fig. 3(i), contains two diagonal peaks and cross-peaks, although the peaks are more broadened due to the stronger system-bath coupling. Firstly, we observe that the PME spectra, Fig. 3(i)(l), cannot replicate the double peak structure found in the numerically exact result, because the electronic coupling $\Omega$ is renormalized to zero (see Section V B). The WCME spectra, Fig. 3(i)(k), on the other hand do accurately reflect the presence of all peaks. For the chosen parameters however, the Lamb shift in the WCME spectra, Eq. (18), overestimates the splitting between the peaks compared to the PT-TEMPO results. This difference is illustrated in more detail in Fig. 4(a): for both PT-TEMPO and the WCME, the relative splitting between peaks ($\delta E - 2\Omega$, Eq. (18)) increases with temperature. For the WCME, this increase is linear for large temperatures ($\frac{\omega_c}{2k_B T} \ll 1$). Although the splitting predicted by PT-TEMPO agrees well with the Lamb shift below $T \approx 4 \text{ ps}^{-1}$, an increasing discrepancy arises when the temperature is increased. Since the Lamb shift is directly proportional to $\alpha$ for our choice of spectral density, Eq. (4), this discrepancy disappears when $\alpha$ is decreased, depicted by the red dashes and orange cross markers in Fig. 4(a) for $\alpha = 0.02$.

While the WCME gives more accurate predictions for peak splitting in the intermediate coupling regime as

temperature is decreased, its prediction for peak height breaks down. As shown in Figs. 4(b-d), the ratio between diagonal peak amplitudes increases with decreasing temperature. While this is also the case for the numerically exact results, the WCME predicts a sharper and higher amplitude peak for the $|-\rangle$ transition at sufficiently low $T$ (Fig. 4(b)), and similarly a broader peak for $|+\rangle$ transition. From the WCME in Eq. (14), we can observe that the transition rate $\gamma_1$ is always finite, while $\gamma_2$ tends to zero as temperature decreases. As a result of this and the absence of a pure dephasing rate, the optical response predicted by the WCME decays more slowly at low temperatures compared to PT-TEMPO, leading to a more pronounced peak corresponding to the $|-\rangle$ state [14]. As temperature increases on the other hand, Fig. 4(c), $\gamma_1$ and $\gamma_2$ become more comparable in magnitude, decreasing the relative difference between diagonal peak amplitudes.

## VII. CONCLUSIONS

In this work, we have presented a tensor network method for the efficient calculation of multi-time correlation functions, based on a process tensor framework. We have employed our technique to simulate linear absorption and 2D spectra in three different system-bath coupling regimes, comparing the results to a weak-coupling and polaron master equation. Here we observed that in an intermediate coupling regime, both master equations break down in the following ways: for an Ohmic spectral density, a polaron master equation fails to resolve the two energy transitions probed in our model. A weak-coupling master equation on the one hand overestimates peak splitting (given by the Lamb shift) at high temperatures. At low temperatures on the other hand, it overestimates the decay time of the optical response, leading to discrepancies in peak amplitude ratios compared to our numerically exact results.

We furthermore note that several methods were recently developed that exploit the time translational invariance of the process tensor. These improve the scaling of the PT-MPO construction algorithm with the number of time steps [30, 31]. Employing such algorithms could further reduce the computational effort of calculating multi-time correlation functions with process tensors. This could be particularly useful for more computationally expensive problems than explored in this work.

Finally, we hope that beyond light-harvesting, the methods employed in this work may aid in the study of other non-Markovian open quantum systems, such as coherent dynamics in semiconductor quantum dots and the polaronic wave functions of excitons in perovskite materials [67–70].
ACKNOWLEDGMENTS

For insightful comments on an earlier version of this
paper, we would like to thank Martin Plenio, Susana
Huelga, Nicola Lorenzoni and Jaemin Lim. R. d. W.

acknowledges support from EPSRC (EP/W524505/1).
B. W. L. and J. K. acknowledge support from EPSRC
(EP/T014032/1). A. W. C. wishes to acknowledge sup-
port from ANR Project ACCEPT (Grant No. ANR-19-
CE24-0028).

[1] S. Mukamel, *Principles of nonlinear optical spectroscopy* (Oxford University Press, Oxford, 1999).

[2] M. Cho, Coherent two-dimensional optical spectroscopy, Chem. Rev. **108**, 1331 (2008).

[3] E. Collini, 2D Electronic Spectroscopic Techniques for Quantum Technology Applications, J. Phys. Chem. C **125**, 13096 (2021).

[4] G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mančal, Y. C. Cheng, R. E. Blankenship, and G. R. Fleming, Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems, Nature **446**, 782 (2007).

[5] K. L. Lewis and J. P. Ogilvie, Probing photosynthetic energy and charge transfer with two-dimensional electronic spectroscopy, J. Phys. Chem. Lett. **3**, 503 (2012).

[6] S. H. Yeh, R. D. Hoehn, M. A. Allodi, G. S. Engel, and S. Kais, Elucidation of near-resonance vibronic coherence lifetimes by nonadiabatic electronic-vibrational state character mixing, Proc. Natl. Acad. Sci. U. S. A. **116**, 18263 (2019).

[7] L. Wang, M. A. Allodi, and G. S. Engel, Quantum co-
herences reveal excited-state dynamics in biophysical sys-
tems, Nat. Rev. Chem. **3**, 477 (2019).

[8] E. Bukartė, A. Haufe, D. Paleček, C. Büchel, and D. Zigmantas, Revealing vibronic coupling in chlorophyll c1 by polarization-controlled 2D electronic spectroscopy, Chem. Phys. **530**, 110643 (2020).

[9] J. Cao, R. J. Cogdell, D. F. Coker, H. G. Duan, J. Hauer, U. Kleinekathöfer, T. L. Jansen, T. Mančal, R. J. Dwayne Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H. S. Tan, R. Tempelaar, M. Thorwart, E. Thyraug, S. Westenhoff, and D. Zigmantas, Quantum biology revisited, Sci. Adv. **6**, eaaz4888 (2020).

[10] Y. Kim, F. Bertagna, E. M. D'souza, D. J. Heyes, L. O. Johannissen, E. T. Nery, A. Pantelias, A. S. P. Jimenez, L. Slocombe, M. G. Spencer, J. Al-Khalili, G. S. Engel, S. Hay, S. M. Hingley-Wilson, K. Jeevaratnam, A. R. Jones, D. R. Kattnig, R. Lewis, M. Sacchi, N. S. Scrutton, S. R. P. Silva, and J. McFadden, Quantum biology: An update and perspective, Quantum Reports **3**, 80 (2021).

[11] J. S. Higgins, L. T. Lloyd, S. H. Sohail, M. A. Allodi, J. P. Otto, R. G. Saer, R. E. Wood, S. C. Massey, P. C. Ting, R. E. Blankenship, and G. S. Engel, Photosynthesis tunes quantum-mechanical mixing of electronic and vibrational states to steer exciton energy transfer, Proc. Natl. Acad. Sci. U. S. A. **118**, e2018240118 (2021).

[12] H. P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Oxford University Press, 2007).

[13] I. de Vega and D. Alonso, Dynamics of non-Markovian open quantum systems, Rev. Mod. Phys. **89**, 15001 (2017), arXiv:1511.06994.

[14] A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, The role of non-

equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes, Nat. Phys. **9**, 113 (2013).

[15] J. Lim, D. Palecek, F. Caycedo-Soler, C. N. Lincoln, J. Prior, H. Von Berlepsch, S. F. Huelga, M. B. Plenio, D. Zigmantas, and J. Hauer, Vibronic origin of long-lived coherence in an artificial molecular light harvester, Nat. Commun. **6**, 7755 (2015), arXiv:1502.01717.

[16] F. Caycedo-Soler, A. Mattioni, J. Lim, T. Renger, S. F. Huelga, and M. B. Plenio, Exact simulation of pigment-protein complexes unveils vibronic renormalization of electronic parameters in ultrafast spectroscopy, Nat. Commun. **13**, 2912 (2022), arXiv:2106.14286.

[17] S. Talkington and M. Claassen, Linear and Non-Linear Response of Quadratic Lindbladians (2024).

[18] G. Guarnieri, A. Smirne, and B. Vacchini, Quantum regression theorem and non-Markovianity of quantum dynamics, Phys. Rev. A - At. Mol. Opt. Phys. **90**, 022110 (2014), arXiv:1406.5340.

[19] A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen, and G. R. Fleming, Quantum coherence and its interplay with protein environments in photosynthetic electronic energy transfer, Phys. Chem. Chem. Phys. **12**, 7319 (2010).

[20] F. Fassioli, R. Dinshaw, P. C. Arpin, and G. D. Scholes, Photosynthetic light harvesting: Excitons and coherence, J. R. Soc. Interface **11**, 20130901 (2014).

[21] C. C. Jumper, S. Rafiq, S. Wang, and G. D. Scholes, From coherent to vibronic light harvesting in photosynthesis, Curr. Opin. Chem. Biol. **47**, 39 (2018).

[22] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. (N. Y). **349**, 117 (2014), arXiv:1306.2164.

[23] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Non-Markovian quantum processes: Complete framework and efficient characterization, Phys. Rev. A **97**, 012127 (2018), arXiv:1512.00589.

[24] M. R. Jørgensen and F. A. Pollock, Exploiting the Causal Tensor Network Structure of Quantum Processes to Efficiently Simulate Non-Markovian Path Integrals, Phys. Rev. Lett. **123**, 240602 (2019), arXiv:1902.00315.

[25] G. E. Fux, E. P. Butler, P. R. Eastham, B. W. Lovett, and J. Keeling, Efficient Exploration of Hamiltonian Parameter Space for Optimal Control of Non-Markovian Open Quantum Systems, Phys. Rev. Lett. **126**, 200401 (2021), arXiv:2101.03071.

[26] M. Cygorek, M. Cosacchi, A. Vagov, V. M. Axt, B. W. Lovett, J. Keeling, and E. M. Gauger, Simulation of open quantum systems by automated compression of arbitrary environments, Nat. Phys. **18**, 662 (2022).

[27] N. Ng, G. Park, A. J. Millis, G. K.-L. Chan, and D. R. Reichman, Real-time evolution of Anderson impurity models via tensor network influence functionals, Phys. Rev.
B **107**, 125103 (2023), arXiv:2211.10430.

[28] J. Thoenniss, M. Sonner, A. Lerose, and D. A. Abanin, Efficient method for quantum impurity problems out of equilibrium, Phys. Rev. B **107**, L201115 (2023), arXiv:2211.10272.

[29] J. Thoenniss, A. Lerose, and D. A. Abanin, Nonequilibrium quantum impurity problems via matrix-product states in the temporal domain, Phys. Rev. B **107**, 195101 (2023), arXiv:2205.04995.

[30] M. Cygorek, J. Keeling, B. W. Lovett, and E. M. Gauger, Sublinear Scaling in Non-Markovian Open Quantum Systems Simulations, Phys. Rev. X **14**, 011010 (2024), arXiv:2304.05291.

[31] V. Link, H.-H. Tu, and W. T. Strunz, Open quantum system dynamics from infinite tensor network contraction, Phys. Rev. Lett. **132**, 200403 (2024).

[32] Y. Tanimura and R. Kubo, Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath, J. Phys. Soc. Japan **58**, 101 (1989).

[33] Y. Tanimura, Stochastic Liouville, Langevin, Fokker-Planck, and master equation approaches to quantum dissipative systems, J. Phys. Soc. Japan **75**, 82001 (2006).

[34] Y. Tanimura, Numerically "exact" approach to open quantum dynamics: The hierarchical equations of motion (HEOM), J. Chem. Phys. **153**, 20901 (2020).

[35] J. Prior, A. W. Chin, S. F. Huelga, and M. B. Plenio, Efficient Simulation of Strong System-Environment Interactions, Phys. Rev. Lett. **105**, 050404 (2010).

[36] A. W. Chin, Á. Rivas, S. F. Huelga, and M. B. Plenio, Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials, J. Math. Phys. **51**, 10.1063/1.3490188 (2010), arXiv:1006.4507.

[37] D. Tamascelli, A. Smirne, J. Lim, S. F. Huelga, and M. B. Plenio, Efficient Simulation of Finite-Temperature Open Quantum Systems, Phys. Rev. Lett. **123**, 090402 (2019), arXiv:1811.12418.

[38] A. D. Somoza, O. Marty, J. Lim, S. F. Huelga, and M. B. Plenio, Dissipation-Assisted Matrix Product Factorization, Phys. Rev. Lett. **123**, 100502 (2019).

[39] F. Mascherpa, A. Smirne, A. D. Somoza, P. Fernández-Acebal, S. Donadi, D. Tamascelli, S. F. Huelga, and M. B. Plenio, Optimized auxiliary oscillators for the simulation of general open quantum systems, Phys. Rev. A **101**, 052108 (2020), arXiv:1904.04822.

[40] L. S. Sirkina and E. A. Muljarov, Impact of the phonon environment on the nonlinear quantum-dot-cavity qed: Path-integral approach, Phys. Rev. B **108**, 115312 (2023).

[41] G. E. Fux, P. Fowler-Wright, J. Beckles, E. P. Butler, P. R. Eastham, D. Gribben, J. Keeling, D. Kilda, P. Kirton, E. D. C. Lawrence, B. W. Lovett, E. O'Neill, A. Strathearn, and R. de Wit, OQuPy: A Python package to efficiently simulate non-Markovian open quantum systems with process tensors (2024), arXiv:2406.16650.

[42] M. B. Oviedo, C. F. Negre, and C. G. Sánchez, Dynamical simulation of the optical response of photosynthetic pigments, Phys. Chem. Chem. Phys. **12**, 6706 (2010).

[43] A. J. Dunnett, D. Gowland, C. M. Isborn, A. W. Chin, and T. J. Zuehlsdorff, Influence of non-adiabatic effects on linear absorption spectra in the condensed phase: Methylene blue, J. Chem. Phys. **155**, 10.1063/5.0062950 (2021), arXiv:2107.06587.

[44] M. Dufour, V. Steinmetz, E. Izquierdo, T. Pons, N. Lequeux, E. Lhuillier, L. Legrand, M. Chamarro, T. Barisien, and S. Ithurria, Engineering bicolor emission in 2D core/crown CdSe/CdSe1-xTex nanoplatelet heterostructures using band-offset tuning, J. Phys. Chem. C **121**, 24816 (2017).

[45] C. Wang, W. Chi, Q. Qiao, D. Tan, Z. Xu, and X. Liu, Twisted intramolecular charge transfer (TICT) and twists beyond TICT: From mechanisms to rational designs of bright and sensitive fluorophores, Chem. Soc. Rev. **50**, 12656 (2021).

[46] S. Rafiq, B. Fu, B. Kudisch, and G. D. Scholes, Interplay of vibrational wavepackets during an ultrafast electron transfer reaction, Nat. Chem. **13**, 70 (2021).

[47] J. Adolphs and T. Renger, How proteins trigger excitation energy transfer in the fmo complex of green sulfur bacteria, Biophysical Journal **91**, 2778 (2006).

[48] M. Rätsep and A. Freiberg, Electron-phonon and vibronic couplings in the fmo bacteriochlorophyll a antenna complex studied by difference fluorescence line narrowing, Journal of Luminescence **127**, 251 (2007), proceedings of the Ninth International Meeting on Hole Burning, Single Molecule, and Related Spectroscopies: Science and Applications.

[49] N. Lorenzoni, N. Cho, J. Lim, D. Tamascelli, S. F. Huelga, and M. B. Plenio, Systematic coarse graining of environments for the nonperturbative simulation of open quantum systems, Phys. Rev. Lett. **132**, 100403 (2024).

[50] P. Hamm and M. Zanni, Concepts and Methods of 2D Infrared Spectroscopy (Cambridge University Press, 2011).

[51] A. Strathearn, B. Lovett, and P. Kirton, Efficient real-time path integrals for non-Markovian spin-boson models, New J. Phys. **19**, 93009 (2017).

[52] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett, Efficient non-Markovian quantum dynamics using time-evolving matrix product operators, Nat. Commun. **9**, 1 (2018), arXiv:1711.09641.

[53] G. E. Fux, D. Kilda, B. W. Lovett, and J. Keeling, Tensor network simulation of chains of non-Markovian open quantum systems, Phys. Rev. Res. **5**, 033078 (2023), arXiv:2201.05529.

[54] The tempo collaboration, OQuPy: A python 3 package to efficiently compute non-Markovian open quantum systems (2022), 10.5281/zenodo.4428316.

[55] M. Suzuki, General theory of higher-order decomposition of exponential operators and symplectic integrators, Phys. Lett. A **165**, 387 (1992).

[56] D. Gribben, A. Strathearn, G. E. Fux, P. Kirton, and B. W. Lovett, Using the Environment to Understand non-Markovian Open Quantum Systems, Quantum **6**, 847 (2022), arXiv:2106.04212.

[57] R. P. Feynman and F. L. Vernon, The theory of a general quantum system interacting with a linear dissipative system, Ann. Phys. (N. Y.) **24**, 118 (1963).

[58] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory, J. Chem. Phys. **102**, 4600 (1995).

[59] N. Makri and D. E. Makarov, Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory, J. Chem. Phys. **102**, 4600 (1995).

[60] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Operational markov condition for quantum processes, Phys. Rev. Lett. **120**, 040405 (2018).
[61] T. Holstein, Studies of polaron motion: Part I. The molecular-crystal model, Ann. Phys. **8**, 325 (1959).

[62] S. Jang, Theory of coherent resonance energy transfer for coherent initial condition, J. Chem. Phys. **131**, 101104 (2009).

[63] F. A. Pollock, D. P. McCutcheon, B. W. Lovett, E. M. Gauger, and A. Nazir, A multi-site variational master equation approach to dissipative energy transfer, New J. Phys. **15**, 75018 (2013), arXiv:1212.5713.

[64] A. Nazir and D. P. McCutcheon, Modelling exciton-phonon interactions in optically driven quantum dots, J. Phys. Condens. Matter **28**, 103002 (2016), arXiv:1511.01405.

[65] P. R. Eastham, P. Kirton, H. M. Cammack, B. W. Lovett, and J. Keeling, Bath-induced coherence and the secular approximation, Phys. Rev. A **94**, 012110 (2016), arXiv:1508.04744.

[66] R. Hartmann and W. T. Strunz, Accuracy assessment of perturbative master equations: Embracing nonpositivity, Phys. Rev. A **101**, 012103 (2020), arXiv:1906.02583.

[67] E. Cassette, J. C. Dean, and G. D. Scholes, Two-Dimensional Visible Spectroscopy for Studying Colloidal Semiconductor Nanocrystals, Small **12**, 2234 (2016).

[68] S. T. Cundiff, Coherent spectroscopy of semiconductors, Opt. InfoBase Conf. Pap. **12**, 1800171 (2016).

[69] E. Collini, H. Gattuso, L. Bolzonello, A. Casotto, A. Volpato, C. N. Dibenedetto, E. Fanizza, M. Striccoli, and F. Remacle, Quantum Phenomena in Nanomaterials: Coherent Superpositions of Fine Structure States in CdSe Nanocrystals at Room Temperature, J. Phys. Chem. C **123**, 31286 (2019).

[70] W. Tao, Y. Zhang, and H. Zhu, Dynamic Exciton Polaron in Two-Dimensional Lead Halide Perovskites and Implications for Optoelectronic Applications, Acc. Chem. Res. **55**, 345 (2022).
# APPENDIX A: COMPUTATIONAL RESOURCES

For this work, the required PT-MPOs were constructed using the open source package OQuPy [41, 54]. For the Ohmic spectral density used in this work, $J(\omega) = 2\alpha\omega e^{-\omega/\omega_c}$ ($\alpha = 0.1$, $\omega_c = 3.04$ ps$^{-1}$), it took 8.7 mins to construct a PT (150 time steps) at a temperature $T = 13.09$ ps$^{-1}$ on a single core of an Intel i5 (8th Gen) processor. Computing a single four-time correlation function as a function of $\tau_4$ (e.g. $R(\tau_4; \tau_1 = \tau_2 = \tau_3 = 0)$ took 101 s on a single CPU core. This computation time scales linearly with the number of additional time steps when one of the earlier time arguments ($\tau_{1,2,3}$) is varied. For example, the 2D spectra in Fig. 2 were computed over 50 time steps each in $\tau_1$ and $\tau_4$. Taking into account that it is composed of four correlation functions (Eq. (5)), the total 2D spectrum in Fig. 2(a) required 5.6 core hours with a pre-computed PT.

For a set number of time steps, constructing the PT at lower temperatures reduces the computation time. However, since multi-time correlations generally take longer to decay to zero with decreasing temperature, a longer PT has to be constructed to capture the full signal. For example, the PT constructed for Fig. 4(a) and $T = 3.9$ ps$^{-1}$ was 300 time steps in length and took 12 mins to compute on a single CPU core.

# APPENDIX B: CONVERGENCE OF MULTI-TIME CORRELATIONS

The PT-TEMPO algorithm relies on three computational parameters: the time step $\delta t$, the maximum memory length $\Delta K_{max}$ and the maximal relative error in the singular value cutoff $\epsilon_{rel}$ [41]. The product $\delta t \Delta K_{max}$ corresponds to the maximal memory time of the bath captured by the computations. Therefore, the value of $\delta t \Delta K_{max}$ should be larger than the time it takes for the bath autocorrelation function to decay to zero:

$$C(\tau) = \int_{0}^{\infty} d\omega J(\omega) \left[ \cos(\omega\tau) \coth\left(\frac{\omega}{2k_{B}T}\right) - i \sin(\omega\tau) \right]. \quad (A1)$$

For an Ohmic spectral density (Eq. (4)) with $\alpha = 0.1$ and $\omega_c = 3.04$ ps$^{-1}$, we find that $C(\tau)$ decays to a value $10^{-3}$ times smaller than its maximum when $\tau = 10.3$ ps. In this work, we have set $\epsilon_{rel} = 10^{-6}$, $\Delta K_{max} = 1000$ and $\delta t = 0.05-0.1$ ps, which places the cutoff well beyond the memory time of the bath. To justify our choice of $\delta t$ and $\epsilon_{rel}$, we test for numerical convergence by computing the four-time correlation function $R_4(\tau_4; \tau_1 = \tau_2 = \tau_3 = 0)$ in Eq. (5):

$$R_4(\tau_4) = \mathrm{Tr} [\hat{V}_2(\tau_4)\hat{V}_2(0)\hat{V}_2(0)\hat{V}_2(0)\rho_0]. \quad (A2)$$

Fig. A1(a) shows $R_4(\tau_4)$ for a range of $\epsilon_{rel}$ values and a constant time step $\delta t = 0.1$ ps. To better illustrate the convergence, the absolute difference $\Delta R_4$ between the

<figure><img src="image_6.png" /><figcaption>FIG. A1. Numerical convergence of the four-time correlation function $R_4(\tau_4)$ (Eq. (A2)) with respect to (a-b) the truncation error $\epsilon_{rel}$ for a constant time step $\delta t = 0.1$ ps and (c-d) with respect to $\delta t$ and $\epsilon_{rel} = 10^{-6}$. Figures (b) and (c) show the absolute difference between $R_4$ calculated with the most precise convergence parameter ((b) $\epsilon_{rel} = 10^{-8}$, (c) $\delta t = 0.025$ ps) and with each lower precision.</figcaption></figure>

curve with the highest precision ($\epsilon_{rel} = 10^{-8}$) and each subsequent lower precision curve is plotted in Fig. A1(b). Similarly, Fig. A1(d) shows $\Delta R_4$ for different $\delta t$, taking $\delta t = 0.025$ ps as the baseline.