# Extracting Coupling-Mode Spectral Densities with Two-Dimensional Electronic 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. Extracting Coupling-Mode Spectral Densities with Two-Dimensional Electronic Spectroscopy. Journal of Physical Chemistry Letters, 2025, 16 (25), pp.6594-6601. 10.1021/acs.jpclett.5c00928. hal-05406559

HAL Id: hal-05406559

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

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
# Extracting Coupling-Mode Spectral Densities with Two-Dimensional Electronic 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)

**Abstract:** Methods for reconstructing the spectral density of a vibrational environment from experimental data can yield key insights into the impact of the environment on molecular function. Although such experimental methods exist, they generally only access vibrational modes that couple diagonally to the electron system. Here we present a method for extracting the spectral density of modes that couple to the transition between electronic states, using two-dimensional electronic spectroscopy. To demonstrate this, we use a process-tensor method that can simulate two-dimensional electronic spectroscopy measurements in a numerically exact way. To explain how the extraction works, we also derive an approximate analytical solution, which illustrates that the non-Markovianity of the environment plays an essential role in the existence of the simulated signal.

## TOC GRAPHIC

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

Environment-mediated processes play an important role in the photophysics of molecules [1–11]. The electron-vibrational interactions involved in such processes are often characterized by a spectral density that encodes the coupling strength to each mode in the environment [1, 12]. Determination of the spectral density is therefore key to gaining a full description of these interactions. Experimentally, this can be achieved with techniques such as resonance Raman and fluorescence line narrowing spectroscopy [1, 4, 13–16]. These methods are generally sensitive to vibrational modes that tune the electronic energy levels and thus couple diagonally to the electronic eigenstates. However, spectral densities that characterize the off-diagonally coupled modes are challenging to extract [4]. These coupling modes are known to mediate processes such as singlet fission [5], intensity borrowing [17] and efficient excitation transfer in organic chromophores [6–9, 11]. Therefore, the ability to directly extract the coupling-mode spectral density could enhance our understanding of such phenomena, and could provide the potential to leverage the environment to enhance molecular function.

In this Letter, we present a method for extracting the spectral density of coupling modes with two-dimensional electronic spectroscopy (2DES) [1, 18–20]. It relies on directly probing electronic transitions that are coupled vibrationally in the system of interest, resulting in a signal that encodes the spectral density of the coupling modes. We simulate this signal for a three-level chromophore system. To explain how our method works, we derive an approximate analytical solution to the signal, which restricts the vibrational manifold to single phonon interactions. We compare this solution to a numerically exact computation that employs tensor network techniques to capture the influence of the environment.

Using non-linear response theory, the detected signal in a 2DES experiment is usually expressed as a sum of four-time correlation functions [1, 18]. For the 2DES pulse sequence explored in this work, we find that the corresponding multi-time correlation function only has a non-zero solution for a non-Markovian treatment. That is, in the case where the environment can retain a memory of past interactions with the system [12]. In contrast, we find that the signal we study vanishes when applying standard master equation techniques which assume a Markovian environment. Our proposed method for extracting spectral densities may therefore also be relevant to the study of devising experimental protocols for probing non-Markovianity.

Studying non-Markovianity in the context of multi-time correlation functions requires careful consideration. Under a Markovian assumption, one can use the quantum regression theorem to find multi-time correlations straightforwardly [12]. The quantum regression theorem may not always give correct results however, even if the corresponding dynamical map accurately predicts single-time observables [21]. In such cases, computational techniques that can efficiently describe non-Markovian environments are necessary [22–27]. One particular approach is to express the open quantum system as a tensor network that effectively compresses the exponentially large Hilbert space to only the physically most relevant de-

<footnote>* rmadw@st-andrews.ac.uk</footnote>

<footnote>† alex.chin@insp.jussieu.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. a) Sketch of our model, consisting of a three-level system coupled to a bosonic bath. b) Laser pulses in a 2DES experiment, applied at times $\tau_{1,2,3}$. The signal is emitted at time $\tau_4$. Time intervals are denoted by $t_{1,2,3}$. c) Tensor network of seven time steps for calculating a four-time correlation function with a process tensor (red squares) and initial system state $\rho_0^S$ (purple circle). Each dipole operator $\hat{V}$ (orange diamonds) is applied at times $\tau_i$. Propagators related to the system Hamiltonian are denoted by $K$ (blue ovals).</figcaption></figure>

grees of freedom [2, 28–31]. This includes the method employed in this work, which relies on a process tensor (PT) [32] that captures all possible effects of the environment on the system. By casting the process tensor into a matrix product operator (PT-MPO) format, computations can be performed in an efficient and numerically exact way [32–40]. Moreover, since the PT-MPO maps any intermediate control operations onto the final state, this framework lends itself particularly well to the computation of multi-time correlations.

We start by introducing the model we are considering in this Letter, shown in Fig. 1a). It involves a three-level electronic system, linearly coupled to a many-mode vibrational bath. This corresponds to a Hamiltonian of three components, $\hat{H} = \hat{H}_S + \hat{H}_B + \hat{H}_I$, which respectively describe the electronic system, vibrational bath and system-bath interaction:

$$ \hat{H} = (\epsilon - \Omega) |1\rangle\langle 1| + (\epsilon + \Omega) |2\rangle\langle 2| + \sum_k \omega_k b_k^\dagger b_k + (|1\rangle\langle 2| + \text{H.c.}) \sum_k (g_k b_k + g_k^* b_k^\dagger) \quad (1) $$

The average energy of the electronic eigenstates, $|1\rangle$ and $|2\rangle$, is set by $\epsilon$, which defines the overall electronic energy scale, while $\Omega$ defines their splitting. The vibrational modes of the bath with frequencies $\omega_k$ are described by creation (annihilation) operators $b_k^\dagger(b_k)$. The coupling of each mode to the system is set by $g_k$ and characterized by the spectral density $J(\omega) = \sum_k |g_k|^2 \delta(\omega - \omega_k)$.

A model of this form can describe a range of organic chromophore systems, such as carotenoids [41] and methylene blue [42]. The photophysics of such molecules are generally described by three relevant energy levels, $S_0$, $S_1$ and $S_2$, corresponding to the states $|0\rangle$, $|1\rangle$ and $|2\rangle$ in our model. Alternatively, $|1\rangle$ and $|2\rangle$ could represent the lowest excited states of two coupled chromophores in a single-excitation subspace, as found in photosynthetic complexes [6–8, 43]. Organic chromophores generally absorb light in the visible frequency range, while the char-

acteristic timescale of the electron dynamics is typically slower, on the order of femtoseconds to picoseconds [7–9]. We will therefore show our results in a rotating frame, by effectively setting $\epsilon = 0$ in subsequent calculations. Additionally we will set $\Omega = 2 \text{ ps}^{-1}$. We do not consider explicit coupling between the excited and ground state manifolds in the Hamiltonian; instead we will consider coupling between the ground and excited states via the externally applied pulses in a 2DES experiment, as discussed below.

The bath modes represent the intramolecular and solvent vibrations, with which the system can interact. Importantly, the coupling between the system and bath is of an off-diagonal form, meaning that the electronic excited states $(|1\rangle, |2\rangle)$ are connected via the coupling to the bath. To test our proposed extraction method, we will consider two different spectral densities. Firstly, we extract a spectral density of an Ohmic form:

$$ J_1(\omega) = 2\alpha_1\omega e^{-\frac{\omega}{\omega_1}}, \quad (2) $$

where $\alpha_1$ is a dimensionless parameter that sets the system-bath coupling strength and $\omega_1$ denotes the cutoff frequency of the bath. Ohmic spectral densities are commonly used in optical spectroscopy simulations, although realistic spectral densities for chromophore models often contain more structure [4, 11, 14–16]. We therefore secondly investigate a spectral density with some additional structure, containing an Ohmic and super-Ohmic component:

$$ J_2(\omega) = 2\alpha_2\omega e^{-\frac{\omega}{\omega_2}} + 2\alpha_3 \frac{\omega^3}{\omega_{c3}^2} e^{-\frac{\omega}{\omega_3}}, \quad (3) $$

where $\alpha_{2,3}$ and $\omega_{2,3}$ are respectively the system-bath coupling parameters and the cutoff frequencies of the Ohmic and super-Ohmic contributions.

Having introduced the model, we now discuss how the coupling-mode spectral density may be extracted with a 2DES experiment. In general, the response in a 2DES experiment consists of a number of interaction pathways, depending on the experimental phase-matching conditions [1, 18–20]. Using non-linear response theory, each pathway can be expressed as a four-time correlation function [1, 18]. To extract the spectral densities in Eqs. (2) and (3), we employ a particular pulse sequence that outputs a signal corresponding to a correlation function of the form:

$$ R = \text{Tr} \left[ \hat{V}_2(\tau_4) \hat{V}_1(\tau_3) \hat{V}_1(\tau_2) \hat{V}_2(\tau_1) \rho_0 \right]. \quad (4) $$

As shown in Fig. 1a), this pulse sequence involves two different dipole operators that probe two distinct electronic transitions: $\hat{V}_1 = |1\rangle\langle 0| + \text{H.c.}$ and $\hat{V}_2 = |2\rangle\langle 0| + \text{H.c.}$. The operators are evaluated at the times of laser pulse interactions and signal emission, denoted by $\tau_i$. We will use $t_i$ to refer to the time intervals between pulses, as shown in Fig. 1b). We will take the initial system+bath state $\rho_0$ as the tensor product of the system ground state
$|0\rangle\langle 0|$ and the thermal state of the bath at temperature $T = 100$ K.

We next give a brief explanation of why the pulse sequence in Eq. (4) allows for the extraction of the spectral density; this will be explored in more detail below. The sequence starts at time $\tau_1$ by creating an optical coherence between the ground and second excited state $|2\rangle$. Due to the off-diagonal system-bath coupling in Eq. (1), the electronic state can then only evolve into the first excited state manifold $|1\rangle$ by absorbing or emitting vibrational excitations (phonons). Probing the system with the second pulse ($\hat{V}_1$) produces a response that depends on such transitions, and hence encodes the system-bath entanglement accumulated during $t_1$. Measuring the linear response Tr $[\hat{V}_1(\tau_2)\hat{V}_2(\tau_1)\rho_0]$ is however not sufficient to extract the coupling-mode spectral density, as this two-time correlation function exactly vanishes: this pulse sequence has transferred the initial thermal vibrational state to vibrational coherences, which disappear when tracing over the bath at this stage. The two further operations in the 2DES signal can reverse the creation of vibrational coherences, allowing for a non-zero signal. The dependence of this signal on $t_2$ can then be used to directly probe the spectral density of the vibrational modes involved in the $|1\rangle \leftrightarrow |2\rangle$ transitions.

In order to compute the four-time correlation function $R$ in Eq. (4) in a numerically exact way, we use a PT-MPO that can capture any non-Markovian effects on the system's evolution. In particular, we will use the process tensor-time evolving matrix product operator (PT-TEMPO) method, which constructs the PT-MPO from a path integral formalism [32, 34–36, 40, 44]. We provide a brief overview of the computation method here; for a more detailed description, we refer the reader to Refs. [40, 45].

Fig. 1c illustrates how multi-time correlations are calculated with our method. Because the evolution resulting from the bath and system-bath interaction in Eq. (1)) is encoded in the process tensor, the initial input state is the reduced system density matrix $\rho_0^S$ (purple circle). The tensor network is expressed in Liouville space, such that density matrices are represented by vectors and superoperators by matrices. $\rho_0^S$ is therefore drawn as a tensor with a single leg, while the dipole operators (orange diamonds) and system propagators $K$ (two per single time step, blue ovals) have two legs. The red squares represent the PT-MPO, drawn here for seven time steps. To calculate a multi-time correlation function, the system propagators and dipole operators are connected to the PT-MPO at the relevant time steps [45–47]. The bonds in the network are subsequently contracted, which corresponds to tracing over the system subspace (gray semi-circle). Crucially, the process tensor does not depend on system Hamiltonian or any measurement superoperators. It therefore only needs to be constructed once and can be re-used when varying the times at which the dipole operators are evaluated. We provide the convergence parameters and computational resources required in Section 1 of the Supporting Information.

To gain a deeper physical insight beyond the numerical computation, we next derive an approximate analytical expression for the signal in Eq. (4). In contrast to a Markovian master equation, this derivation explicitly tracks the evolution of the (approximate) eigenstates of the full system+bath Hamiltonian in Eq. (1), excited by the $\hat{V}_{1,2}$ operators, and does not assume factorisation between the system and bath at intermediate times. To do this, we restrict the bath Hilbert space to configurations in which only a single phonon is absorbed or emitted after each application of a dipole operator. We start from a particular configuration of the thermal ground state,

$$|0, \mathbf{n}\rangle = |0\rangle \otimes |n_1, n_2, \dots, n_i, \dots\rangle, \quad (5)$$

where $n_i$ corresponds to the occupation of mode $i$. The initial density matrix of the system and bath is therefore $\rho_0 = \sum_{\mathbf{n}} P_{\mathbf{n}} |0, \mathbf{n}\rangle\langle 0, \mathbf{n}|$ where $P_{\mathbf{n}} \propto \exp(-\beta \sum_i \omega_i n_i)$ is the Boltzmann-weighted probability of the occupation vector $\mathbf{n}$. We then apply the $\hat{V}_{1,2}$ operators at times $\tau_{1,2,3,4}$ according to Eq. (4) and track the time evolution of the excited states. This is diagrammatically depicted in Fig. 2a), where we show the evolution of four bath modes in each electronic manifold, represented by the horizontal lines in each sub-diagram. The number of dots correspond to the number of phonon excitations in each mode and thus to an example state $\mathbf{n}$ drawn from the thermal configuration in $\rho_0$. From the form of the system-bath interaction in Eq. (1), it follows that initial excitation of the ground state with $\hat{V}_2$ leads to the formation of a superposition in both excited state manifolds during $t_1$. The ket component of $\rho_0$ will be:

$$|2, \mathbf{n}(t_1)\rangle \approx e^{-iE_{\mathbf{n}}t_1} \left( \alpha(t_1) |2, \mathbf{n}\rangle + \sum_i [\beta_i(t_1) |1, \mathbf{n} + \mathbf{e}_i\rangle + \gamma_i(t_1) |1, \mathbf{n} - \mathbf{e}_i\rangle] \right), \quad (6)$$

neglecting higher order changes of the vibrational number state as noted above. We have denoted the absorption/emission of a phonon in mode $i$ by the vector $\pm \mathbf{e}_i$, such that $[\mathbf{e}_i]_j = \delta_{ij}$. The coefficients $\alpha, \beta_i$ and $\gamma_i$ evolve during the first time interval $t_1$ according to the time-dependent Schrödinger equation. Additionally, the phase arising from the total energy of the initial configuration $E_{\mathbf{n}} = \sum_i \omega_i n_i$ has been factored out. Since the bra-side of the density matrix will evolve with exactly the conjugate of this phase, this term vanishes in the final expression. In Fig. 2a), we illustrate one process contributing to $|2, \mathbf{n}(t_1)\rangle$, corresponding to the creation of a phonon in a single mode (additional blue dot), governed by the evolution of $\beta_i$ in Eq. (6).

At $\tau_2$, we apply the $\hat{V}_1$ operator to obtain

$$\hat{V}_1 |2, \mathbf{n}(t_1)\rangle \approx \sum_i [\beta_i(t_1) |0, \mathbf{n} + \mathbf{e}_i\rangle + \gamma_i(t_1) |0, \mathbf{n} - \mathbf{e}_i\rangle]. \quad (7)$$
<figure><img src="image_4.png" /><figcaption>FIG. 2. a) Diagrammatic representation of one contribution to the system-bath evolution during the spectral density extraction protocol. Each horizontal line represents a vibrational bath mode, depicted within the manifold of the relevant electronic state (|0⟩, |1⟩, |2⟩). For simplicity, only four bath modes are drawn while the model contains a continuum of modes. The dots on each line show the number of phonons in a particular mode. The absorption/emission of a phonon in mode i is denoted by a blue solid dot or a dashed outline respectively. The laser pulses and output signal are represented by orange arrows and the corresponding dipole operators. The double blue arrows represent the build-up of a superposition state in the excited states manifold, resulting from the system-bath coupling $g_i$. b) Double-sided Feynman diagram of the response in Eq. (4). The evolution of the internal state during $t_1$ and $t_3$ is denoted by a dashed line and blue arrow. For clarity, only the bath states resulting from the initial absorption and eventual emission of a phonon in mode i are shown in a) and b).</figcaption></figure>

At this point, we have created a superposition of bath states in the electronic ground state that will evolve during $t_2$ according to the bath Hamiltonian:

$$e^{-iE_n(t_1+t_2)} \sum_i \left[ \begin{aligned} & \beta_i(t_1)e^{-i\omega_i t_2} |0, \mathbf{n} + \mathbf{e}_i\rangle \\ & + \gamma_i(t_1)e^{+i\omega_i t_2} |0, \mathbf{n} - \mathbf{e}_i\rangle \end{aligned} \right]. \quad (8)$$

Note here that if we would now take the trace to calculate the two-time correlation function $\text{Tr}[\hat{V}_1(\tau_2)\hat{V}_2(\tau_1)\rho_0]$, we would find a solution of zero, because the vibrational state in Eq. (8) is never $|\mathbf{n}\rangle$.

We now repeat the same procedure in reverse by applying $\hat{V}_1$ and $\hat{V}_2$ at times $\tau_3$ and $\tau_4$ respectively, described in detail in Section 2 of the Supporting Information. As shown in Fig. 2a), the final two operators enable the backflow of information from environment to system, generating a non-zero signal and allowing the extraction of the system-bath correlations at the final time $\tau_4$. We similarly find a set of coefficients that the describe the evolution of the system-bath states during $t_3$. Through solving the time-dependent Schrödinger equation and employing Weisskopf-Wigner theory, we obtain a set of closed equations of motion [12, 48]. We show the determination of the coefficients as a function of $t_1$ in detail in Section 3 of the Supporting Information. The final solution for $R$ is then found by using the definition of the spectral density and Boltzmann-weighted probabilities of occupations to give:

$$R = \int_{0}^{\infty} d\omega [R_{\beta}(t_1, t_3, \omega)e^{-i\omega t_2} + R_{\gamma}(t_1, t_3, \omega)e^{i\omega t_2}], \quad (9)$$

with

$$R_{\beta}(t_1, t_3, \omega) = J(\omega)(N(\omega) + 1)f_{+\omega}(t_1)f_{+\omega}(t_3); \quad (10)$$

$$R_{\gamma}(t_1, t_3, \omega) = J(\omega)N(\omega)f_{-\omega}(t_1)f_{-\omega}(t_3); \quad (11)$$

$$f_{\pm\omega}(t) = \frac{\left( e^{-\Gamma_2 t - i E'_2 t} - e^{-\Gamma_1 t - i (\pm\omega + E'_1)t} \right)}{(E'_2 - E'_1 - i\Gamma_1 + i\Gamma_2 \mp \omega)}. \quad (12)$$

To obtain the above form of $f_{\pm\omega}(t)$, we have for simplicity assumed that the influence of the environment on the evolution of $\alpha(t_1)$ in Eq. (6) is Markovian, and similar for the evolution during $t_3$. This results in the decay rates $\Gamma_{1,2}$ and Lamb shifts $S(\pm 2\Omega)$ that modify the system eigenvalues in Eq. (1), specified in Section 3 of the Supporting Information. Here the Lamb shifts have been absorbed in the system eigenvalues, such that $E'_1 = \epsilon - \Omega + S(-2\Omega)$ and $E'_2 = \epsilon + \Omega + S(2\Omega)$. Were one to use a fully Markovian approach to calculate $R$, this would neglect the system bath entanglement that exists during the time evolution. Since, as discussed above, that entanglement is crucial to the measured signal, a fully Markovian approach cannot capture this effect.

Our final expression for the signal $R$ consists of two components, $R_\beta$ and $R_\gamma$, that represent the build-up of system-bath entanglement during the first and final time interval. $R_\beta$ encodes the creation and subsequent destruction of a phonon during $t_1$ and $t_3$ respectively, while $R_\gamma$ represents the reverse process. During the waiting time $t_2$, all excitations are contained within the bath, such that the density matrix evolves according to a simple phase (Eq. (9)). Inspired by beating maps which can be extracted from 2DES measurements [49, 50], we can exploit this simple phase dependence to extract the spectral density by Fourier transforming over the waiting
<figure><img src="image_5.png" /><figcaption>FIG. 3. Simulated response for an Ohmic spectral density, Eq. (2), and two system-bath coupling strengths $\alpha_1$. a-b) The modulus of the response in the time domain as a function of $t_2$ and $t_1 = t_3 = 0.01$ ps, calculated with the analytical solution in Eq. (9) (blue) and PT-TEMPO (orange dashed). c-d) Extracted bath information computed with PT-TEMPO (orange dots), compared to the expected result (grey line). The offset has not been subtracted from the results, discussed in the main text. e-f) The extracted spectral density computed with PT-TEMPO and compared to the expected result (grey line). The orange dots represent the raw output without the offset subtracted; green dots correspond to the processed data.</figcaption></figure>

time:

$$ \tilde{R}(t_1, \omega_2, t_3) = \int_0^\infty d\omega \int_0^\infty dt_2 e^{i\omega_2 t_2} \text{Re}[R(t_1, t_2, t_3)]. \quad (13) $$

Here we have chosen to Fourier transform only the real part of $R$, such that when we solve the transform using the relation

$$ \int_0^\infty dt e^{\pm i\omega t} = \pi\delta(\omega) \pm i\mathcal{P}\frac{1}{\omega}, \quad (14) $$

we can discard the Cauchy principle value component, denoted by $\mathcal{P}$, by selecting the real part again after transforming. We then obtain:

$$ \text{Re}[\tilde{R}(t_1, \omega_2, t_3)] = \frac{\pi}{2}(\text{Re}[R_\gamma + R_\beta] + \text{Im}[R_\gamma - R_\beta]). \quad (15) $$

For sufficiently small time intervals $t_1$ and $t_3$, the real part of this Fourier transform is exactly proportional to $J(\omega)(2N(\omega) + 1)$. This follows from Taylor expanding Eq. (12):

$$ f_{\pm\omega}(t) = -it + \mathcal{O}(t^2). \quad (16) $$

If $t_1$ and $t_3$ are sufficiently small compared to the exponents in $f_{\pm\omega}(t)$, Eq. (15) thus reduces to:

$$ \text{Re}[\tilde{R}(t_1, \omega_2, t_3)] \approx -\frac{\pi}{2} t_1 t_3 J(\omega_2)(2N(\omega_2) + 1). \quad (17) $$

Finally, we can divide by $-\frac{\pi}{2}t_1t_3(2N(\omega_2) + 1)$ to obtain the spectral density by itself. Rather than directly setting $t_1$ and $t_3$ to a sufficiently small value, we evaluate Eq. (17) for a range of time intervals and fit the solution to a polynomial with variables $t_1$ and $t_3$. Then, for the term $a_{11}t_1t_3$ in the polynomial, the coefficient $a_{11}$ directly corresponds to $-\frac{\pi}{2}J(\omega_2)(2N(\omega_2) + 1)$. For the results discussed below, we extracted this coefficient from six time intervals, ranging from 0 - 0.05 ps.

In an experimental setting, applying this protocol would thus have to involve the following steps: first, measure the real part of $R(t_0, t_2, t_0)$ for a range of waiting times and for a small $t_1, t_3 = t_0$. Then apply a discrete Fourier transform to the measured signal, and finally select its real part.

In Figs. 3 and 4, we show the results of our extraction method for the spectral densities in Eqs. (2) and (3) respectively. Panel a) in Fig. 3 shows the time domain signal corresponding to Eq. (9), where $t_0 = 0.01$ ps. Here we have chosen to show the modulus of the signal, such that the phase envelope $e^{-i\epsilon(t_1+t_3)}$ disappears. We find that for the chosen parameters ($\alpha_1 = 0.02$), the analytical (blue) and numerical (orange dashed) results agree well.

We display the extracted bath information (orange dots) from the process tensor computation in panel b). Broadly, the extracted data agrees well with the actual $J_1(\omega)(2N(\omega) + 1)$ (grey line). We do find however that the extracted curve does not completely decay to zero at large frequencies $\omega_2$, as expected. This becomes clearly visible in panel c), which displays the extracted spectral density $J_1(\omega)$. This is a consequence of the discrete nature of the numerical Fourier transform we are performing, which is discussed in Section 5 of the Supporting Information. Nonetheless, we can successfully extract the spectral density by removing the offset from zero in panel b), before dividing by $\frac{\pi}{2}(2N(\omega) + 1)$ (green dots).

In panels d)-f) in Fig. (3) we have increased the system bath coupling to $\alpha_1 = 0.4$. For increasing system-bath coupling strengths, we expect that our primary approximation in the analytical solution—limiting each interaction to a single quantum exchange—breaks down. The results in panel d) indeed confirm that the analytical solution does not agree as well with the numerically exact computation when $\alpha_1$ is increased. Panel f) shows that we can still successfully extract the spectral density from the numerical solution for this particular value of $\alpha_1$.
<figure><img src="image_6.png" /></figure>

FIG. 4. Simulated response for the spectral density in Eq. (3), and two sets of system-bath coupling strengths $\alpha_{1,2}$. a-b) The modulus of the response in the time domain as a function of $t_2$ and $t_1 = t_3 = 0.01$ ps, calculated with the analytical solution in Eq. (9) (blue) and PT-TEMPO (orange dashed). c-d) Extracted bath information computed with PT-TEMPO (orange dots), compared to the expected result (grey line). The offset has not been subtracted from the results, discussed in the main text. e-f) The extracted spectral density computed with PT-TEMPO and compared to the expected result (grey line). The orange dots represent the raw output without the offset subtracted; green dots correspond to the processed data.

We repeat the extraction process for the spectral density $J_2(\omega)$ in Eq. (3), as presented in Fig. 4. As shown in panels c) and f), our method is successful in resolving both peaks present in the spectral density.

To conclude, we have presented a method for extracting the spectral density of vibrational modes that couple to electronic transitions, using 2D electronic spectroscopy. The protocol would consist of sequentially probing two distinct electronic transitions as a function of the waiting time. This could for example be achieved by varying the polarisation of the pulses to probe states with

orthogonal dipole moments, as is the case for the $Q_x$ and $Q_y$ states in chlorophylls [51–53]. An alternative method for selecting different states could be to vary the color of the pulses [54]. By keeping the first and final time delay sufficiently small, the measured signal in the Fourier domain is directly proportional to the spectral density. We simulated this pulse sequence for a model chromophore coupled to a vibrational environment with a process tensor method that can compute multi-time correlations in a numerically exact way. To gain physical insight into the evolution of the system and environment during the pulse sequence, we compared the numerical results to an approximate analytical derivation. These results show that both an Ohmic spectral density and one with additional structure can be accurately extracted with our method.

We moreover found that the simulated response is zero under a standard Markovian master equation approach, indicating that non-Markovianity plays a significant role. The analytical derivation sheds more light on this: it shows that the response arises from the transfer of excitations, which leave and re-enter the system at a later time during the pulse sequence. Previous studies have suggested that an experimental witness of non-Markovianity should involve the measurement of multi-time correlations [55, 56]. Additional work has shown how a 2D NMR spectroscopy setup could be used to witness initial correlations by sequentially disentangling en entangling the system and environment [57]. Beyond extracting spectral densities, the method introduced in this work could thus be an interesting step towards designing a protocol for a measure of non-Markovianity in open quantum systems.

**Supporting Information:** Information on the computational resources and convergence parameters used for the PT-TEMPO calculations; a detailed description of the time evolution during $t_3$, the determination of the time-dependent coefficients, and of taking the weighted average over thermal states in the analytical derivation; further explanation on the error in the discrete Fourier transform of the data (PDF).

## ACKNOWLEDGMENTS

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 support from ANR Project ACCEPT (Grant No. ANR-19-CE24-0028).

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

[2] 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 coher-
ence and recoherence in pigment-protein complexes, Nat.
Phys. **9**, 113 (2013).

[3] 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).

[4] I. Gustin, C. W. Kim, D. W. McCamant, and I. Franco, Mapping electronic decoherence pathways in molecules, Proc. Natl. Acad. Sci. U.S.A. **120**, e2309987120 (2023).

[5] F. A. Schröder, D. H. Turban, A. J. Musser, N. D. Hine, and A. W. Chin, Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation, Nat. Commun. **10**, 1062 (2019).

[6] J. M. Womick and A. M. Moran, Vibronic enhancement of exciton sizes and energy transport in photosynthetic complexes, J. Phys. Chem. B **115**, 1347 (2011).

[7] E. A. Arsenault, Y. Yoneda, M. Iwai, K. K. Niyogi, and G. R. Fleming, The role of mixed vibronic Qy-Qx states in green light absorption of light-harvesting complex II, Nat. Commun. **11**, 6011 (2020).

[8] S. Sil, R. W. Tilluck, N. Mohan TM, C. H. Leslie, J. B. Rose, M. A. Domínguez-Martín, W. Lou, C. A. Kerfeld, and W. F. Beck, Excitation energy transfer and vibronic coherence in intact phycobilisomes, Nat. Chem. **14**, 1286 (2022).

[9] R. Moya, A. C. Norris, T. Kondo, and G. S. Schlau-
Cohen, Observation of robust energy transfer in the
photosynthetic protein allophycocyanin using single-
molecule pump-probe spectroscopy, Nat. Chem. **14**, 153
(2022).

[10] K. E. Hunter, Y. Mao, A. W. Chin, and T. J. Zuehlsdorff, Environmentally driven symmetry breaking quenches dual fluorescence in proflavine, J. Phys. Chem. Lett. **15**, 4623 (2024).

[11] S. Sohoni, P.-J. E. Wu, Q. Shen, L. T. Lloyd, C. MacGregor-Chatwin, A. Hitchcock, and G. S. Engel, Resonant vibrational enhancement of downhill energy transfer in the c-phycocyanin chromophore dimer, J. Phys. Chem. Lett. **15**, 11569 (2024).

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

[13] L. A. Pachón and P. Brumer, Direct experimental determination of spectral densities of molecular complexes, J. Chem. Phys. **141**, 174102 (2014).

[14] M. Rätsep and A. Freiberg, Electron–phonon and vibronic couplings in the fmo bacteriochlorophyll antenna complex studied by difference fluorescence line narrowing, J. Lumin. **127**, 251 (2007).

[15] M. Rätsep, J. Pieper, K.-D. Irrgang, and A. Freiberg, Excitation wavelength-dependent electron-phonon and electron-vibrational coupling in the cp29 antenna complex of green plants, J. Phys. Chem. B **112**, 110 (2008).

[16] J. Pieper, M. Rätsep, K.-D. Irrgang, and A. Freiberg, Chromophore-chromophore and chromophore-protein interactions in monomeric light-harvesting complex ii of green plants studied by spectral hole burning and fluorescence line narrowing, J. Phys. Chem. B **113**, 10870 (2009).

[17] P.-A. Yin, Q. Wan, Y. Niu, Q. Peng, Z. Wang, Y. Li, A. Qin, Z. Shuai, and B. Z. Tang, Theoretical and Ex-

perimental Investigations on the Aggregation-Enhanced
Emission from Dark State: Vibronic Coupling Effect,
Avd. Electron. Mater. **6**, 2000255 (2020).

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

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

[20] E. Fresch, F. V. Camargo, Q. Shen, C. C. Bellora, T. Pul lerits, G. S. Engel, G. Cerullo, and E. Collini, Two- dimensional electronic spectroscopy, Nat. Rev. Methods Primers **3**, 84 (2023).

[21] G. Guarnieri, A. Smirne, and B. Vacchini, Quantum re-
gression theorem and non-Markovianity of quantum dy-
namics, Phys. Rev. A - At. Mol. Opt. Phys. **90**, 022110
(2014).

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

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

[24] 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).

[25] L. Chen, D. I. G. Bennett, and A. Eisfeld, Calculating nonlinear response functions for multidimensional electronic spectroscopy using dyadic non-markovian quantum state diffusion, J. Chem. Phys. **157**, 114104 (2022).

[26] T. Gera, L. Chen, A. Eisfeld, J. R. Reimers, E. J. Taffet, and D. I. G. B. Raccah, Simulating optical linear absorption for mesoscale molecular aggregates: An adaptive hierarchy of pure states approach, J. Chem. Phys. **158**, 174103 (2023).

[27] 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).

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

[29] 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).

[30] 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).

[31] 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).

[32] 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).

[33] 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).

[34] A. Strathearn, B. Lovett, and P. Kirton, Efficient real-time path integrals for non-Markovian spin-boson models, New J. Phys. **19**, 93009 (2017).
[35] 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).

[36] 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).

[37] 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).

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

[39] 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).

[40] Fux, Gerald E. and Fowler-Wright, Piper and Beckles, Joel and Butler, Eoin P. and Eastham, Paul R. and Gribben, Dominic and Keeling, Jonathan and Kilda, Dainius and Kirton, Peter and Lawrence, Ewen D. C. and Lovett, Brendon W. and O'Neill, Eoin and Strathearn, Aidan and de Wit, Roosmarijn, OQuPy: A Python package to efficiently simulate non-Markovian open quantum systems with process tensors, J. Chem. Phys. **161**, 124108 (2024).

[41] D. Manawadu, T. N. Georges, and W. Barford, Photoexcited state dynamics and singlet fission in carotenoids, J. Phys. Chem. A **127**, 1342 (2023).

[42] 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**, 144112 (2021).

[43] 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).

[44] The TEMPO collaboration, OQuPy: A Python package to efficiently simulate non-Markovian open quantum systems with process tensors. (2024).

[45] R. de Wit, J. Keeling, B. W. Lovett, and A. W. Chin, Process tensor approaches to modeling two-dimensional spectroscopy, Phys. Rev. Res. **7**, 013209 (2025).

[46] 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).

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

[48] V. Weisskopf and E. Wigner, Berechnung der natürlichen Linienbreite auf Grund der Diracschen Lichttheorie, Z. Phys. **63**, 54 (1930).

[49] A. A. Bakulin, S. E. Morgan, T. B. Kehoe, M. W. Wilson, A. W. Chin, D. Zigmantas, D. Egorova, and A. Rao, Real-time observation of multiexcitonic states in ultrafast singlet fission using coherent 2d electronic spectroscopy, Nat. Chem. **8**, 16 (2016).

[50] V. Petropoulos, P. S. Rukin, F. Quintela, M. Russo, L. Moretti, A. Moore, T. Moore, D. Gust, D. Prezzi, G. D. Scholes, E. Molinari, G. Cerullo, F. Troiani, C. A. Rozzi, and M. Maiuri, Vibronic coupling drives the ultrafast internal conversion in a functionalized free-base porphyrin, J. Phys. Chem. Lett. **15**, 4461 (2024).

[51] S. Westenhoff, D. Paleček, P. Edlund, P. Smith, and D. Zigmantas, Coherent picosecond exciton dynamics in a photosynthetic reaction center, J. Am. Chem. Soc. **134**, 16484 (2012).

[52] G. S. Schlau-Cohen, A. Ishizaki, T. R. Calhoun, N. S. Ginsberg, M. Ballottari, R. Bassi, and G. R. Fleming, Elucidation of the timescales and origins of quantum electronic coherence in LHCII, Nat. Chem. **4**, 389 (2012).

[53] Y. Song, A. Schubert, E. Maret, R. K. Burdick, B. D. Dunietz, E. Geva, and J. P. Ogilvie, Vibronic structure of photosynthetic pigments probed by polarized two-dimensional electronic spectroscopy and ab initio calculations, Chem. Sci. **10**, 8143 (2019).

[54] B. S. Prall, D. Y. Parkinson, G. R. Fleming, M. Yang, and N. Ishikawa, Two-dimensional optical spectroscopy: Two-color photon echoes of electronically coupled phthalocyanine dimers, J. Chem. Phys. **120**, 2537 (2004).

[55] S. Milz, M. S. Kim, F. A. Pollock, and K. Modi, Completely positive divisibility does not mean markovianity, Phys. Rev. Lett. **123**, 040401 (2019).

[56] G. E. Fux, *Process tensor networks for non-Markovian open quantum systems*, PhD thesis, University of St Andrews, St Andrews (2022), available at https://doi.org/10.17630/sta/259.

[57] P. Jatakia, S. Vinjanampathy, and K. Saha, Detecting initial correlations via correlated spectroscopy in hybrid quantum systems, Sci. Rep. **11**, 20718 (2021).

[58] Note that while we make this Markovian approximation for the dissipation acting on $β_i(t_1)$, our overall treatment is not Markovian, in that we rely on system-environment entanglement.
# SUPPORTING INFORMATION: EXTRACTING COUPLING-MODE SPECTRAL DENSITIES WITH TWO-DIMENSIONAL ELECTRONIC SPECTROSCOPY

## I. COMPUTATIONAL RESOURCES

The numerical computations discussed in the main text were performed with the PT-TEMPO method, which is part of the open source package OQuPy [40, 44]. The precision and efficiency of this method relies on the following three parameters: the time step $\delta t$, the maximum bath memory length $\Delta K_{max}$ and the maximal relative error in the singular value cutoff $\epsilon_{rel}$. The meaning and function of these parameters are discussed in detail in Ref. [40]. For this work, the computations were performed with $\delta t = 0.01$ ps, $\Delta K_{max} = 1000$ and $\epsilon_{rel} = 10^{-7}$. We show in Ref. [45] that these parameters produce numerically converged results.

For the results shown in the main text, the computational cost depends on the system-bath coupling strength $\alpha$ and the particular form of the spectral density $J(\omega)$. Generally, stronger coupling and more structured spectral densities are computationally more expensive [40]. For an Ohmic spectral density and $\alpha = 0.02$ in Fig. 3a), it took 78 mins to compute the four-time correlation function $R(t_2; t_1 = t_3 = 0.01$ ps) for 300 time steps on a single core of a Intel i5 (8th Generation) processor, with a pre-computed process tensor. The construction of the process tensor took 43 mins. In contrast, the computation for the structured spectral density with $\alpha_2 = 0.6$ in Fig. 4b) took 17 core hours for 300 time steps, plus 16 core hours to construct the process tensor.

## II. EVOLUTION DURING THE FINAL INTERVAL

Picking up from Eq. (8) in the main text, we apply $\hat{V}_1$ to excite the state into the $|1\rangle$ manifold at $\tau_3$:

$$e^{-iE_n(t_1+t_2)} \sum_i \left[ \begin{aligned} & \beta_i(t_1)e^{-i\omega_i t_2} |1, \mathbf{n} + \mathbf{e}_i\rangle \\ & + \gamma_i(t_1)e^{+i\omega_i t_2} |1, \mathbf{n} - \mathbf{e}_i\rangle \end{aligned} \right]. \quad (S1)$$

During $t_3$, both terms will evolve back into a superposition in both excited state manifolds:

$$|1, \mathbf{n} \pm \mathbf{e}_i(t_3)\rangle \approx e^{-iE_n t_3} \alpha_i^{\pm}(t_3) |1, \mathbf{n} \pm \mathbf{e}_i\rangle \\ + e^{-iE_n t_3} \sum_j \left[ \begin{aligned} & \beta_{ij}^{\pm}(t_3) |2, \mathbf{n} \pm \mathbf{e}_i + \mathbf{e}_j\rangle \\ & + \gamma_{ij}^{\pm}(t_3) |2, \mathbf{n} \pm \mathbf{e}_i - \mathbf{e}_j\rangle \end{aligned} \right], \quad (S2)$$

where $\alpha_i^{\pm}$, $\beta_{ij}^{\pm}$, $\gamma_{ij}^{\pm}$ are now new coefficients that will have similar forms to $\alpha$, $\beta_i$ and $\gamma_i$ respectively.

After $t_3$, we apply the final operator $\hat{V}_2$ and take the trace. Since we are applying all four dipole operators to the ket-side of the density matrix, the bra-side will still

be in the ground state after $t_3$, i.e. $\langle 0, \mathbf{n}|$. To obtain a non-zero trace, the only terms in the sum over $i$ and $j$ that contribute to the final solution are those for which $i = j$ and such that the final number state is $\mathbf{n}$. For example:

$$\mathrm{Tr} \sum_{ij} \beta_{ij}^{-}(t_3) |0, \mathbf{n} - \mathbf{e}_i + \mathbf{e}_j\rangle \langle 0, \mathbf{n}| = \sum_{ij} \delta_{ij} \beta_{ij}^{-} \quad (S3)$$

Thus, in Eq. (S2), only $\gamma_{ii}^{+}$ and $\beta_{ii}^{-}$ contribute to the final solution. After taking the trace, the correlation function will then be of the form:

$$\sum_i [\beta_i(t_1)e^{-i\omega_{ii}t_2}\gamma_{ii}^{+}(t_3) + \gamma_i(t_1)e^{i\omega_{ii}t_2}\beta_{ii}^{-}(t_3)], \quad (S4)$$

We discuss the determination of the coefficients in the above expression in the next Section.

## III. WEISSKOPF-WIGNER APPROACH

Here we describe in detail the derivation of the coefficients in Eq. (S4), focussing on the first time interval $t_1$. First, we insert the coefficients $\alpha(t_1)$, $\beta_i(t_1)$, $\gamma_i(t_1)$ into the time-dependent Schrödinger equation:

$$i \frac{d\alpha(t_1)}{dt} = E_2\alpha(t_1) + \sum_i \left[ g_i^* \sqrt{n_i + 1} \beta_i(t_1) + g_i \sqrt{n_i} \gamma_i(t_1) \right], \quad (S5)$$

$$i \frac{d\beta_i(t_1)}{dt} = (E_1 + \omega_i)\beta_i(t_1) + g_i \sqrt{n_i + 1}\alpha(t_1) \\ + \sum_m \left[ g_m^* \sqrt{n_m + 1} b_{im}^+(t_1) + g_m \sqrt{n_m} c_{im}^+(t_1) \right], \quad (S6)$$

$$i \frac{d\gamma_i(t_1)}{dt} = (E_1 - \omega_i)\gamma_i(t_1) + g_i^* \sqrt{n_i} \alpha(t_1) \\ + \sum_m \left[ g_m^* \sqrt{n_m + 1} b_{im}^{-}(t_1) + g_m \sqrt{n_m} c_{im}^{-}(t_1) \right], \quad (S7)$$

where $E_{1,2}$ are the eigenvalues of $\hat{H}_S$ in Eq. (1) and $g_i$ the system-bath coupling of mode $i$ in $\hat{H}_I$. $n_i$ refers to the occupation of mode $i$. In order to ensure dissipation of the same order in $g_i$ as in the evolution of $\alpha(t_1)$, the equations of motion for $\beta_i(t_1)$ and $\gamma_i(t_1)$ additionally include coupling terms to states involving bath modes $m$, which
evolve according to the coefficients $b_{im}^{\pm}(t_1)$ and $c_{im}^{\pm}(t_1)$:

$$i \frac{d b_{im}^{+}(t_1)}{dt} = (E_2 + \omega_i + \omega_m) b_{im}^{+}(t_1) + g_m \sqrt{n_m + 1} \beta_i(t_1) \quad (S8)$$

$$i \frac{d c_{im}^{+}(t_1)}{dt} = (E_2 + \omega_i - \omega_m) c_{im}^{+}(t_1) + g_m^* \sqrt{n_m} \beta_i(t_1) \quad (S9)$$

$$i \frac{d b_{im}^{-}(t_1)}{dt} = (E_2 - \omega_i + \omega_m) b_{im}^{-}(t_1) + g_m \sqrt{n_m + 1} \gamma_i(t_1) \quad (S10)$$

$$i \frac{d c_{im}^{-}(t_1)}{dt} = (E_2 - \omega_i - \omega_m) c_{im}^{-}(t_1) + g_m^* \sqrt{n_m} \gamma_i(t_1). \quad (S11)$$

Now, following a Weisskopf-Wigner approach [12, 48], we first solve the equations of motion for $b_{im}^{\pm}(t_1)$ and $c_{im}^{\pm}(t_1)$ in terms of $\beta_i(t_1)$ and $\gamma_i(t_1)$. For example:

$$b_{im}^{+}(t_1) = -i g_m \sqrt{n_m + 1} \int_0^t dt' e^{-i(E_2 + \omega_i + \omega_m)(t-t')} \beta_i(t') \quad (S12)$$

$$c_{im}^{+}(t_1) = -i g_m^* \sqrt{n_m} \int_0^t dt' e^{-i(E_2 + \omega_i - \omega_m)(t-t')} \beta_i(t'), \quad (S13)$$

and similar for $d_{im}^{-}(t_1)$ and $c_{im}^{-}(t_1)$. We can now insert these back into the equations of motion for $\beta_i(t_1)$ and $\gamma_i(t_1)$. We will explicitly show the next steps for $\beta_i(t_1)$. Inserting $b_{im}^{+}(t_1)$ and $c_{im}^{+}(t_1)$:

$$\frac{d\beta_i(t_1)}{dt} = -i(E_1 + \omega_i)\beta_i - ig_i\sqrt{n_i + 1}\alpha(t_1) - \int_0^t dt' f_m^+(t-t')e^{-\omega_i(t-t')}\beta_i(t'), \quad (S14)$$

where

$$f_m^+(\tau) = \sum_m |g_m|^2 e^{-E_2\tau} [(n_m + 1)e^{-i\omega_m\tau} + n_m e^{i\omega_m\tau}]. \quad (S15)$$

In Eq. (S14), the complex exponential dependent on $\omega_i$ has been explicitly separated from $f_m^+(t-t')$ to emphasize that this term will cancel in the next step. In order to solve for $\beta_i(t_1)$, we first move into the interaction picture $\beta_i(t_1) = \tilde{\beta}_i(t_1)e^{-i(E_1+\omega_i)t}$. For simplicity, we assume the influence of the bath modes $m$ on the time evolution of $\beta_i(t_1)$ is Markovian [58], though we note that an exact solution could be found via Laplace transforms [12]. We then obtain:

$$\frac{d\tilde{\beta}_i(t_1)}{dt} = -i g_i \sqrt{n_i + 1} e^{i(E_1 + \omega_i)t} \alpha(t_1) - \Lambda_1 \tilde{\beta}_i(t_1). \quad (S16)$$

with

$$\Lambda_1 = \int_0^\infty d\tau' f_m^+(\tau) e^{iE_1\tau}. \quad (S17)$$

$\Lambda_1$ is composed of a real part $\Gamma_1$ leading to exponential decay and an imaginary part corresponding to the negative Lamb shift $S(-2\Omega)$:

$$\Gamma_1 = \pi N(2\Omega)J(2\Omega); \quad (S18)$$

$$S(\mu) = \mathcal{P} \int_{-\infty}^{\infty} \frac{J(\omega)(N(\omega) + 1) + J(-\omega)N(-\omega)}{\mu - \omega} \quad (S19)$$

We can follow exactly the same procedure for $\gamma_i(t_1)$, which leads to the same dissipation term. Thus, for the equations of motion for $\beta_i(t_1)$ and $\gamma_i(t_1)$ we now have:

$$i \frac{d\beta_i(t_1)}{dt} = (E_1 + \omega_i - i\Lambda_1)\beta_i(t_1) + g_i \sqrt{n_i + 1}\alpha(t_1) \quad (S20)$$

$$i \frac{d\gamma_i(t_1)}{dt} = (E_1 - \omega_i - i\Lambda_1)\gamma_i(t_1) + g_i^* \sqrt{n_i}\alpha(t_1). \quad (S21)$$

We now repeat the same Weisskopf-Wigner procedure on the differential equation for $\alpha(t_1)$: we solve Eqs. (S20) and (S21) in terms of $\alpha(t_1)$, and insert $\beta_i(t_1)$ and $\gamma_i(t_1)$ into Eq. (S5). By moving into the interaction picture, making a Markov approximation and solving the differential equation we obtain:

$$\alpha(t_1) = e^{-i(E_2 + S(2\Omega))t - \Gamma_2 t}, \quad (S22)$$

where $S(2\Omega)$ now corresponds to a positive Lamb shift with a dissipation rate $\Gamma_2 = \pi(N(2\Omega) + 1)J(2\Omega)$. Finally, substituting $\alpha(t_1)$ into the Eqs. (S20) and (S21) gives us the final expressions for $\beta_i(t_1)$ and $\gamma_i(t_1)$:

$$\beta_i(t_1) = g_i \sqrt{n_i + 1} \left( \frac{e^{-i(E'_1 + \omega_i)t_1 - \Gamma_1 t_1} - e^{-iE'_2 t_1 - \Gamma_2 t_1}}{E'_2 - E'_1 - \omega_i - i\Gamma_2 + i\Gamma_1} \right), \quad (S23)$$

$$\gamma_i(t_1) = g_i^* \sqrt{n_i} \left( \frac{e^{-i(E'_1 - \omega_i)t_1 - \Gamma_1 t_1} - e^{-iE'_2 t_1 - \Gamma_2 t_1}}{E'_2 - E'_1 + \omega_i - i\Gamma_2 + i\Gamma_1} \right), \quad (S24)$$

where the Lamb shifts have been absorbed into the system eigenvalues, i.e. $E'_1 = E_1 + S(-2\Omega)$ and $E'_2 = E_2 + S(2\Omega)$.

We now repeat the entire procedure for the $t_3$-dependent coefficients, keeping in mind that only the coefficients with indices $i = j$ contribute to the final solution, as shown in Eq. (S3). This results in the following solutions:

$$\gamma_{ii}^{+}(t_3) = g_i^* \sqrt{n_i + 1} \left( \frac{e^{-i(E'_1 + \omega_i)t_3 - \Gamma_1 t_3} - e^{-iE'_2 t_3 - \Gamma_2 t_3}}{E'_2 - E'_1 - \omega_i - i\Gamma_2 + i\Gamma_1} \right), \quad (S25)$$

$$\beta_{ii}^{-}(t_3) = g_i \sqrt{n_i} \left( \frac{e^{-i(E'_1 - \omega_i)t_3 - \Gamma_1 t_3} - e^{-iE'_2 t_3 - \Gamma_2 t_3}}{E'_2 - E'_1 + \omega_i - i\Gamma_2 + i\Gamma_1} \right). \quad (S26)$$
<figure><img src="image_7.png" /><figcaption>FIG. S1. Magnitude of the offset from zero of Re[ $\tilde{R}(t_1, t_2, t_3)$ ], when $\omega_2 \to \infty$, for different values of the time step $\delta t$. Data points were computed using the Ohmic spectral density $J(\omega)$ with $\alpha = 0.4$.</figcaption></figure>

IV. OBTAINING THE FINAL SOLUTION

We can now insert the solutions for the coefficients in Eqs. (S23)-(S26) into Eq. (S4). To obtain the final solution for the signal $R$ stated in the main text, we lastly take the thermal average over all vibrational configurations. Since the solution is linear in configurations, this corresponds to a weighted sum over occupation vectors $n$:

$$R = \sum_{n,i} P_n \left[ \beta_i(t_1) e^{-i\omega_{ii}t_2} \gamma_{ii}^{+}(t_3) + \gamma_i(t_1) e^{i\omega_{ii}t_2} \beta_{ii}^{-}(t_3) \right], \quad (S27)$$

where $P_n$ is the Boltzmann-weighted probability of $n$, stated in the main text. Additionally, since the environment is composed of a continuum of modes, we can replace the sum over modes $i$ by an integral over vibrational frequencies, resulting in the spectral density dependence. For example:

$$\sum_{n,i} P_n |g_i|^2 (n_i + 1) = \int_0^\infty d\omega J(\omega) (N(\omega) + 1), \quad (S28)$$

using the definition of the spectral density, $J(\omega) = \sum_i |g_i|^2 \delta(\omega - \omega_i)$, and the Bose-Einstein occupation number $N(\omega)$. This results in the final solution of $R$, described in the main text in Eqs. (9)-(12).

V. ERROR IN FOURIER TRANSFORM AS FUNCTION OF TIME STEP

In the main text, we mentioned that the extracted data has a finite offset from zero at large frequencies compared to the exact form $J(\omega_2)(2N(\omega_2)+1)$, which decays to zero. In Fig. S1 we show the magnitude of this offset, calculated by performing a discrete Fourier transform

on the analytical solution to Re[$R(t_1, t_2, t_3)$] in Eqs. (9) and (13), for different values of the time step $\delta t$. The largest time step displayed, $\delta t = 0.01$ ps, corresponds to the value used in the main results in Figs. 3 and 4. By extrapolation, we can see that the offset tends to zero for an infinitesimally small time step. Since both the analytical and numerical PT-TEMPO results decay to a finite value rather than to zero, we can conclude that this error is likely a consequence of the discrete nature of the numerical Fourier transform. Because the time step sets the finite length of the frequency window, the signal may be distorted due to aliasing effects, which cause replicas of the signal in the Fourier domain to overlap.

Thus we find that, if it were computationally feasible to perform the PT-TEMPO calculation with a smaller time step, the offset after performing the discrete Fourier transform would disappear.