Two-dimensional spectroscopy beyond the perturbative limit: The influence of finite pulses and detection modes

The success of 2D spectroscopy derives from its ability to correlate electronic transitions at ultrafast time delays by achieving high resolutions both in frequency and in time.^1–5 By portraying the data in a 2D map, coupling between states may be revealed and dark states can be exposed via their coupling to bright states. Moreover, the two dimensions separate homogeneous broadening and inhomogeneous broadening into perpendicular directions, thereby untangling fast fluctuations from the static or slowly fluctuating energy levels.⁶ By investigating the time evolution of the spectra, information about the dynamics of the system, including vibrations, relaxation, and dephasing, can be obtained.

In multi-pulse spectroscopy, the possibilities explode as the complexity of the experiment increases. However, some designs are more broadly applicable and meaningful and thus establish themselves as routine measurements in labs around the world. As pump–probe became the go-to experiment in two-pulse spectroscopy, 2D spectroscopy is becoming the standard measurement with three pulses.

Pump–probe, which is the simplest spectroscopic technique used to study ultrafast electronic dynamics, measures the frequency spectrum for a range of pump-to-probe times and thus provides information on the transient absorption. The extra pulse in 2D spectroscopy enables two frequency spectra to be resolved, one for excitation and one for detection, which are plotted against each other for a range of pulse delays in order to time-resolve the correlation between excitation and detection frequencies.

2D spectroscopy aims to measure the third-order response of a quantum system to an electric field. By clever design of pulses, this goal can be achieved with a high fidelity. First and foremost, the width of the pulse should be short to enhance the time resolution. Equally important is the phase of the pulse, which is precisely controlled in order to separate the third-order signal from other orders.

Since its inception, 2D spectroscopy has employed heterodyne-detection (HD) to record the amplitude and the phase of the signal.^1–4,7 In HD, the signal is measured by interference with a much stronger electric field, a local oscillator. As HD 2D spectroscopy (HD2D) has matured, so has the understanding of the relationship between the underlying dynamics of the sample and their spectral features. The technique is most widely known for shedding new light on quantum coherent transport in photosynthesis, particularly the hotly debated oscillations in 2D spectra and their origin.^8–12 However, HD2D has been shown to be a broadly applicable method.

In recent years, an alternative approach to measure the third-order signal has gained a lot of interest, namely, fluorescence-detected 2D spectroscopy (FD2D). Until recently, FD2D was primarily a proof-of-concept method, focusing on experimental developments,^13–16 but it is now establishing itself as an incisive tool with real applications.^17–20 In the following, we focus on comparing HD2D and FD2D; for a broader perspective on multidimensional ultrafast spectroscopy, Refs. 7, 21, and 22 are excellent resources.

In this work, we compare the two detection schemes to better understand what they have in common and where they differ. In order to not limit the study to idealized experiments, we chose to go beyond the double-sided Feynman diagrams (DSFD)^5,23,24 and simulate the spectra non-perturbatively.^25–28 This approach enables more aspects of the experiments, such as the pulse duration and the pulse amplitude, to be captured within a self-contained model, without the need for additional calculations within the standard diagrammatic method to study these effects.^29–38 However, other sources of divergence between the detection methods could be interesting to investigate, e.g., dephasing and relaxation mechanisms, transition strengths and quantum yields, the effect of vibrations, and the energy level structure.

A great advantage of HD2D is that the third-order signal is spatially separated into three spots depending on the type of interaction that occurred between the sample and the electric pulse, with similar phase evolutions interfering constructively in these directions. For historical reasons, these contributions are named the rephasing (R), the nonrephasing (NR), and the double-quantum coherence (DQC) signals.

In FD2D, a fourth pulse, similar to the first three, replaces the local oscillator in the conventional experiment, with the effect that the desired third-order signal is encoded into the excited state populations, from which fluorescence is collected during an acquisition time. The modulation in the integrated fluorescence, as a function of the pulse intervals, is then Fourier transformed to give the 2D spectra, which, as for HD2D, can be separated into the R, NR, and DQC contributions.

Other actions, such as photoelectron or photo-ion emission^39,40 or photocurrent,⁴¹ could also be recorded to produce 2D spectra with complementary information, but we devote our attention to FD2D as it is currently more popular and also because it affords a more direct comparison to HD2D.

Detecting the fluorescence instead of the polarization has a number of potential benefits. The most known advantage is that signals from small volumes, in principle, single molecules,^17,42 can be used to generate 2D spectra, whereas the conventional technique requires sample sizes larger than the wavelength of the pulse. This opens up the possibility to pierce through the ensemble and study isolated systems or variations across a sample.¹⁸ 

In addition, inherent differences between the two detection schemes affect the selection of interaction pathways, which, in turn, give rise to discrepancies in the spectra. By contrasting HD and FD 2D spectra, it is possible to infer new information, which would otherwise be inconclusive with either detection method. Karki et al. compared the HD and FD DQC spectra of LH2 and were thus able to deduce that the initial excitation is shared between the two bacteriochlorophyll rings, contrary to the generally accepted picture up until that point.²⁰ 

Recently, Malý and Mančal suggested that the acquisition time, which has no analog in HD2D, could be varied in order to isolate specific contributions to the total signal, e.g., the exciton–exciton annihilation.⁴³ This opens up a new window into the underlying dynamics, but long measurement times increase the incoherent mixing of linear signals arising from nonlinear population dynamics.⁴⁴ However, careful analysis can differentiate between true nonlinear signals and incoherently mixed linear signals.⁴⁵ 

Experimentally, FD2D is currently more demanding than its counterpart. This is partly because two pulse delays need to be scanned instead of one. Moreover, labs that have implemented the fluorescence-detection setup are relatively few and do not enjoy the same level of experience and literature to draw upon. However, the increasing rate of articles published suggests that FD2D is leaving infancy and is becoming a powerful addition to the toolbox.

Only a handful of papers have tackled the theory side of FD2D.^43,45–48 Using the double-sided Feynman diagrams derived from perturbation theory, it is readily shown that each diagram from the traditional method has an equivalent diagram in fluorescence-detection, which, in addition, has an extra set of excited state absorption (ESA) diagrams.⁴⁶ Although this knowledge forms an important basis for the interpretation of FD2D spectra, the differences do not end there.

Knowing how challenging it can be to analyze traditional 2D spectra, where the various contributions to the signal are well understood, it is of great interest to explore all possible ways that FD2D spectra can deviate from its counterpart in order to fully unlock the potential of the method.

In order to enhance the comparison between HD and FD and to clarify the interpretation of the 2D spectra, we make some simplifications. Whereas the FD2D experiment requires two pulse intervals to be scanned and the HD2D version requires only one, we chose to scan both in our simulations. In addition, the coherent detection in the HD experiment is performed by taking the instantaneous expectation value of the transition dipole moment; we do not model the local oscillator with a nonequilibrium Green’s function QED approach.⁴⁹ Both of the above-mentioned choices are in line with how the perturbative simulations using double-sided Feynman diagrams are typically performed.

Moreover, as FD2D only requires a single absorber to generate spectra, we do not include a distribution of energy levels as a source of inhomogeneous broadening. For the same reason, we disregard the vector property of the transition dipole moment and treat them as scalars. To promote as identical conditions as possible, we do the same for the HD2D simulations, which do require an ensemble of absorbers—however, only the positions are different.

In our simulations, we use the rotating frame instead of the laboratory frame or the quasi-rotating frame.⁵⁰ We also sample uniformly, opting for simplicity rather than trying to reduce the computational cost with non-uniform sampling.^51,52 We do not investigate the polarization dependency, which is sometimes exploited to select specific pathways,^53–57 but it is known that the early time dynamics suffer from incorrect pulse-ordering artifacts.⁵⁸ 

We assume non-interacting chromophores in our calculations, although the delocalization of excitation energy in many cases is unexpectedly long-range⁵⁹ and its role in energy transfer is a hot topic in the field.^19,60,61

First, we introduce the theoretical and computational methods used throughout this article, starting with the evolution of a quantum system interacting semiclassically with an electric field. We then describe the equations used to detect the third-order signal and construct the 2D spectra for both heterodyne and fluorescence-detection. Finally, we explain the model system that we use in our simulations.

The vast majority of 2D spectroscopy models employ the semiclassical approximation, where the state of the system is described quantum mechanically, but the field, $Er,t$⁠, is described classically,

Here, A_n(t − t_n) describes the envelope of the nth pulse centered at t_n, ω is the frequency of the field, k_n is the wave vector of the nth pulse, and ϕ_n is the phase angle of the nth pulse.

The coupling between the quantum state and the classical field is given by the interaction Hamiltonian

where μ is the transition dipole moment operator, which is assumed to be a scalar for simplicity. For an otherwise isolated quantum system, the dynamics obeys the time-dependent Schrödinger equation

where H₀ is the system Hamiltonian in the absence of the field. In the condensed phase, however, the fluctuations in the quantum system’s environment perturb the system sufficiently that an isolated-quantum system theory is inadequate to describe the dynamics, as it cannot include relaxation and dephasing processes. To accommodate these effects, an open-quantum system approach is needed. For the sake of simplicity, we use the Lindblad equation to propagate the state, which is now represented by the reduced density operator, ρ(t),⁶² 

where an off-diagonal $Lj=am†an$ operator represents spontaneous relaxation or excitation, a diagonal L_j operator represents a dephasing process, and N denotes the number of different decoherence channels. The formal solution to the differential equation in (4), in the case that the Hamiltonian is time-independent, is found by integrating both sides of the equation: $ρ(t)=eLtρ(0)$⁠.

Having established the dynamical equations for an open quantum system interacting with a classical field, we proceed to discuss the different detection schemes for the third-order signal. Traditionally, heterodyne-detection has been used, where the emitted light is interfered with a local oscillator and the electric field is measured. More recently, fluorescence-detection has become an increasingly viable and attractive method as technological advances push toward single-molecule 2DES experiments, exploiting the high sensitivity of fluorescence detectors.

Heterodyne-detected 2DES employs phase-matching to isolate the third-order response signal. In short, three pulses impart their unique wave vectors to the system, thus causing the radiated fields from individual emitters to constructively interfere in the direction of the vector sum k_signal = ∓k₁ ±k₂ + k₃, where the upper combination corresponds to the rephasing signal and the lower combination corresponds to the nonrephasing signal. A third (possible) phase–matching direction, k_DQC = k₁ + k₂ − k₃, produces the so-called double-quantum coherence signal, but this will not be considered in our study.

Importantly, a detectable signal in the phase-matched directions relies on an ensemble of oscillating dipoles to coherently add up to a macroscopic polarization in the sample. While this enables high signal-to-noise ratios, it also puts a lower bound on the spatial resolution: ≳λ³, where λ is the radiation wavelength.

Each optically active molecular system, j, picks up a unique phase factor by virtue of its position, r_j, since the electric field is given by

where a Gaussian envelope function is used and the ϕ_n phase has been dropped as it is not necessary for phase-matching. The electric field couples to the transition dipole moment of the system, which makes the exp(−ik_p · r_j) phase factor act as a book-keeper of the transitions. By using linearly independent k-vectors for each of the three pulses, it is possible to differentiate which transitions were caused by which pulses.

However, averaging over the polarization of the individual molecules will cause the third-order signal to vanish due to the random phases, unless the polarization is measured in one of the phase-matching directions. That is, to read out the desired signal, the wave vector of the local oscillator, k_signal, must be chosen such that it simultaneously cancels the phase factors absorbed by all molecular systems for the relevant pathways. For 2DES, we are interested in the part of the response that stems from a single interaction with each of the three pulses. This reduces the number of phase-matching directions to three: −k₁ + k₂ + k₃, k₁ − k₂ + k₃, and k₁ + k₂ − k₃. A schematic of the phase-matching method is shown in Fig. 1(a).

To achieve the phase-matching in a non-perturbative calculation, the dynamics of each individual molecule j must be computed separately using Eq. (4). Here, we assume that the individual molecules evolve independently and that each molecule is initially in its ground state. For simplicity, we employ the identical Hamiltonian and Lindblad operators on each individual quantum system. A generalization to a distribution of Hamiltonians and Lindbladians is straightforward and with no added computational cost; however, this would be a source of both inhomogeneous and homogeneous broadening, which could potentially make it more difficult to interpret (the differences between) the HD and FD 2D spectra.

The polarization in the chosen phase-matching direction is given by

Note that in this equation, the local oscillator is interfered with the full signal, i.e., there is no spatial selectivity to select only the phase-matched one, in contrast to a typical experimental setup where different order signals are spatially separated.

The last step is to perform a two-dimensional Fourier transform from the time domain to the frequency domain,

where “+” is used for the nonrephasing and double-quantum coherence signals and “−” is used for the rephasing signal. In the following, only the real parts of the signals will be investigated, as this information is far more often used in studies—the imaginary part is mostly neglected.

Instead of detecting the third-order polarization with a local oscillator, the integrated fluorescence intensity can be used to report on the nonlinear response of the sample. Fluorescence-detection is extremely sensitive, which allows for studies of single molecules. As fluorescence is an incoherent process, no phase-matching condition can exist, and there is therefore no need to perform the experiment in a noncollinear setup, which is standard for heterodyne-detected 2DES. The advantage of the collinear geometry is that it facilitates rapid data acquisition and that it is inherently phase stable. Two major drawbacks are the need to scan two coherence times and that it is necessary to scan a number of different pulse phases for each unique combination of pulse delays in order to extract the desired rephasing and nonrephasing signals. The latter requirement is achieved by two similar but different methods known as phase-modulation,¹³ which will not discussed here, and phase cycling.⁶³ 

In the collinear setup, the signal is no longer spatially separated according to the order of the light–matter interaction; it now contains all orders. To isolate the rephasing and nonrephasing (and DQC) signals, one must exploit the unique phase evolutions of the third-order signals. This can be done by Fourier transforming the signal with respect to the phase for each pulse interval so that it matches the characteristic phases of the R, NR, and DQC signals. Practical limitations prevent a continuous Fourier transform, but it can be shown that 27 different phase combinations are enough to extract the third-order signals.⁶⁴ Figure 1(b) is a schematic of the FD2D experiment emphasizing the collinear geometry and the phase-cycling approach, where the pulses are tagged with specific phases, which are then rotated or cycled with respect to each other.

From a purely theoretical point of view, the excited state populations have the same dependence on the interactions with the four pulses as the integrated fluorescence, meaning that the populations can be used as proxies for the measured signal.⁴⁵ In principle, one can distinguish between different excited state populations, but since it is difficult to detect fluorescence from distinct excited states, we neglect this possibility.⁴⁷ It is, however, relevant for other action spectroscopies.^39–41 

The excited state population after interactions with four pulses with phases ϕ₁₋₄ and delay times τ, T, and t can be found by summing all the contributing coherence transfer pathways,

Here, α, β, γ, δ count the number of positive phase factors minus the number of negative phase factors imparted on the system by the respective pulses, which in the language of double-sided Feynman diagrams amounts to the number of arrows pointing to the left minus the number of arrows pointing to the right. $p̃$ denotes that the contributions are specific to the particular set of α, β, γ, δ and that the absorbed phase factors have been taken out.

Because it is assumed that the initial state is diagonal, ρ₀ = |g⟩⟨g|, and the final state is diagonal, the following condition must be fulfilled: α + β + γ + δ = 0. Consequently, α, β, γ, δ are dependent, which allows us to define ϕ₂₁ ≡ ϕ₂ − ϕ₁, ϕ₃₁ ≡ ϕ₃ − ϕ₁, and ϕ₄₁ ≡ ϕ₄ − ϕ₁. Equation (8) then becomes

However, for the 2DES signals, we are not interested in expressing the total population as a sum of its contributions with unique phase factors. Instead, we wish to isolate the individual contributions, in particular, the rephasing, nonrephasing, and double coherence signals, which are given by $p̃(β=1,γ=1,δ=−1)$⁠, $p̃(β=−1,γ=1,δ=−1)$⁠, and $p̃(β=1,γ=−1,δ=−1)$⁠, respectively. This is achieved by Fourier transforming the total population with respect to the characteristic phase evolutions of the third-order signals,

As noted above, a continuous Fourier transform is impractical, so a discrete Fourier transform is used,

As shown by Tan,⁶⁴ a phase-cycling scheme of L × M × N = 3 × 3 × 3 is sufficient when only contributions from fourth order and below are considered, i.e., |α| + |β| + |γ| + |δ| ⩽ 4. This gives $Δϕ21=Δϕ31=Δϕ41=2π3$⁠.

It follows that the appropriate electric field needed to perform the phase-cycling to extract the R, NR, and DQC signals is

where l, m, n are cycled between 0, 1, and 2 and the resulting populations from the 27 combinations added for each set of {t₁, t₂, t₃, t₄}. Note that the k_n vector is dropped as the k_n · r phase factors cancel out in the collinear geometry. The R, NR, and DQC spectra are then found using Eq. (11) and performing the usual

and inserting the appropriate β, γ, δ for the R, NR, and DQC signals. The sign of the Fourier transforms are the same as for the HD2D case.

It should be noted, however, that higher-order processes can have identical phase evolutions as the third-order processes; for example, three interactions with the first pulse can leave the system in the same state as only one interaction with the pulse will. The phase-cycling operation will not be able to filter out these higher-order contributions and the spectra will become distorted as a consequence. Care should therefore be taken to ensure that the third-order contributions are predominant, both in experiments and in simulations.

Instead of using the final populations as proxies⁴⁵ to report on the molecular response, the integrated fluorescence can be recorded. This is also more in line with the actual experiment. Theoretically, the integrated fluorescence can be found as the integral of the spontaneous relaxation, which is given by

and similarly for relaxations between other states. t_acq denotes the acquisition time in which fluorescence is collected. If it is possible to distinguish particular transitions, which is the case for some action spectroscopies, they can give rise to multiple 2D spectra and provide even more incisive information on the studied system. Otherwise, the total integrated fluorescence is just the sum of the individual transitions. Note that the acquisition time can be varied, which opens up another window into the dynamics of the sample.

Once the spontaneous relaxation has been recorded, the data undergo the same phase cycling process as for the population data to construct the rephasing and nonrephasing spectra.

The theoretical framework for both the HD and FD spectra can be simplified substantially by invoking a few approximations, resulting in a set of equations each representing unique spectroscopic pathways. These are commonly depicted as double-sided Feynman diagrams (DSFD), which provide a visual connection to the physical processes and are constructed by following a list of rules dictated by the topology of the model system. A neat property of the DSFD equations is that the Fourier transforms can be calculated analytically in Liouville space,⁶⁵ 

as each propagator $G(t)≡eLt$ is independent of other time variables. This obviates the cumbersome process of numerically integrating the state of the system for all realizations of the pulse intervals.

More detailed accounts on the subject can be found elsewhere, but we include three diagrams that highlight the difference between HD and FD. Figure 2 shows excited state absorption diagrams that evolve identically until the detection event as is also evident from their corresponding equations,

Here, $Vij⃗$ denotes the transition dipole moment operator acting to the right, causing a transition from j to i, $Vij⃗|j•|=|i•|$⁠. Conversely, $Vij⃖$ acts to the left, causing a transition from i to j: $Vij⃖|•i|=|•j|$⁠. The fluorescence yields are given by γ_i, and t_acq is the acquisition time. The initial state, on the right of the equations, is assumed to be the ground state, hence |0⟩⟨0|, but can be completely general if desired.

In theory, dropping the integral and taking only the projected populations would give perfectly valid third-order signals; however, this would not be in line with the experiment. Additionally, the acquisition time parameter can be exploited to reveal more information about the system.

Using FD, the last pulse interaction can bring the state to two different excited state populations: The lower excited state, Eq. (19), which is equivalent to the HD diagram, Eq. (18), and the higher excited state, Eq. (20), which has no HD counterpart. Moreover, the two FD diagrams have opposite signs resulting in cancellations of these contributions. Depending on the model system, this can lead to peaks appearing with one detection method, but not the other.

For future reference, we note that rephasing (R) designates contributions/diagrams with opposite phase evolutions after the first and third pulse, whereas nonrephasing (NR) contributions/diagrams oscillate with the same sign in these periods. Double-quantum coherence (DQC) contributions/diagrams are in this sense a special case of NR, with the distinction that they are not static or slowly evolving after the second pulse but oscillate with double frequency.

In Fig. 3, we depict our model system, which we adopted from Ref. 47. This four-level system captures a lot of the physics of 2D spectroscopy without complicating the interpretation unnecessarily. It also encompasses the case of an excitonic dimer, which is of fundamental interest and a natural starting point for a comparison between heterodyne-detection and fluorescence-detection. Moreover, the fact that the authors of Ref. 47 simulated the FD2D spectrum using the phase-modulation technique, and not phase-cycling, allows us to validate our calculations and compare the two methods of extracting the desired signal. We use the same values for the energy levels as given in Ref. 47,

where E_gr is the ground state energy, E₁ and E₂ correspond to the B850 and B800 absorption bands of the inner and outer rings of bacteriochlorophylls in the light-harvesting complex of purple bacteria,^66–70 and E₃ is the energy of the doubly excited exciton. The optical transitions between the energy levels are all allowed and are equally strong: eE₀μ_ij = 8 meV, with the exception of the transition from ground to the highest-excited state, which is forbidden. Electronic relaxation is modeled by the jump operators, the off-diagonal Lindblad operators: L₁₀ = |0⟩⟨1|, L₂₁ = |1⟩⟨2|, and L₃₂ = |2⟩⟨3|, which are scaled by Γ₁₀ = 4.13 μeV, Γ₂₁ = 4.13 meV, and Γ₃₂ = 13.78 meV, respectively. This yields relaxation times (=ℏ/Γ_j) of 160 ps, 160 fs, and 48 fs, which is representative of a dimer system. Dephasing is included similarly with diagonal Lindblad operators: L₀₀ = |0⟩⟨0|, L₁₁ = |1⟩⟨1|, L₂₂ = |2⟩⟨2|, and L₃₃ = |3⟩⟨3|, which are all scaled by the same dephasing strength Γ_Dephasing = 41.3 meV.

The model system is initialized with the ground state fully populated and then propagated according to Eq. (4) with the appropriate (light–matter) interaction Hamiltonian until it is detected as described by Eq. (6) or Eq. (15). As FD only requires an individual photoactive system, in contrast to HD which requires an ensemble, we chose to neglect the distributions of energy levels and transition dipole moments as they potentially could obscure the effects we want to study. However, the computational cost of incorporating these details into our model would be negligible.

In order to ensure a random sampling of the position-dependent phase, it is necessary to sample a volume of space with dimensions greater than the wavelength of the incident field. This is achieved by randomly generating the molecular positions and multiplying with a factor, which make intermolecular distances large compared to the wavelength of the pulses.

By applying the theory of HD and FD 2DES on the model system, we are able to non-perturbatively simulate 2D spectra under a range of different conditions and to study effects that will be applicable to a broad range of samples and experiments. Our goal is to investigate under which limits the extraction of the third-order signal breaks down due to pulse intensity or pulse duration and hence to understand the optimal regime for 2D spectra.

For comparison and as a starting point for the discussion of the effect of a non-idealized pulse, we compute the corresponding spectra in the impulsive limit using the double-sided Feynman diagrams laid out in Ref. 47 (see Fig. 4). The structure of our model system means that there are (for both the rephasing and nonrephasing signal) four ground state bleach diagrams, four stimulated emission diagrams, four excited-state absorption diagrams, and for FD four additional excited-state absorption diagrams.

It is evident that the cross peaks in the HD spectrum are weaker than in the FD spectrum as would be expected by investigating the phase evolutions of the diagrams, which reveal cancellations for the HD cross peak diagrams. In fact, the cross peaks in the HD spectrum turn out to be artifacts of finite Fourier transforms of the first (t₁ = 50 fs) and last (t₃ = 50 fs) pulse delay and disappear as these are increased. The reason for using finite Fourier transforms is simply to match the broadening arising from finite Fourier transforms in the non-perturbative simulations below. Potentially, employing a smooth cutoff at the boundary, or nonuniform sampling, could enhance the quality of the spectrum, but for simplicity, we show the bare discrete Fourier transform result without additional smoothing or processing.

Before we discuss the response to changes in the pulse amplitude, we note two general differences between FD and HD that we observe in Fig. 4: (1) The diagonal peaks are equally intense in the HD spectra, whereas the FD counterparts are biased to the lower diagonal peak. The explanation is that our model, as in reality, collects fluorescence emanating from the lowest-excited state; this favors the lower diagonal peak over the upper diagonal, which relies on the appropriate pathways to relax to the lowest-excited state prior to fluorescing. In other words, the upper diagonal grows in as the fluorescence integration time is increased. For the HD version, both the highest and lowest-excited states are detected, and the symmetry is only broken by the relaxation and dephasing processes. (2) The cross peak amplitudes are stronger for FD than for HD. This discrepancy is a direct consequence of the additional pathways selected by the FD scheme, specifically by excited state absorption (ESA) processes where the fourth pulse brings the state to a higher-lying population. The remaining pathways where the frequency evolutions correspond to a cross peak are equivalent for FD and HD and largely cancel out due to pairs of pathways with opposite phase progressions, but the cancellation is incomplete in non-idealized realizations, i.e., real-life experiments and simulations with finite pulses. Moreover, lifting the condition that the transition strengths to the first and second excited state are equal would break further symmetries leading to fewer cancellations.

To calibrate our non-perturbative simulations, we begin by finding appropriate pulse amplitudes for each detection scheme. For FD, we start off by replicating the calculation performed by Damtie et al.⁴⁷ but with phase-cycling replacing phase-modulation. As can be seen in the  Appendix, we were able to produce the same spectra with our phase-cycling method using the same parameters, specifically a peak interaction energy of max(eE(t)μ) = 8 meV, thus validating our approach and finding a suitable starting point for the pulse amplitude.

Intuitively, one might expect the HD and FD versions to operate within the same parameter space and to be similarly affected by changes to the experiment and system variables/parameters. However, when we repeated the simulation with the HD model, keeping all common parameters identical, we found that the two methods do not share the same parameter regime for the extraction of third-order signals. Specifically, the amplitude of the pulses required a sevenfold increase for easy acquisition of the HD signal, i.e., without an exceedingly high number of randomly positioned absorbers. Note that the HD third-order signal is a function of the number of absorbers (and their positions and orientations), so a large number of absorbers can to some degree compensate for low pulse amplitudes. We use 10 000 randomized positions throughout this work, which is sufficient to produce well resolved spectra even when far from the third-order sweet spot. Increasing the number of absorbers further did not result in significant improvements. However, the number of independent molecules in typical experiments will far outnumber this value. Care should also be taken that the volume in which the molecules are dispersed is large enough, as for very small volumes the position-dependent phases are no longer fully randomized.

One of the challenges that arise with explicit simulation of 2D spectra is that of the pulse/field power, which affects how much the state of the sample changes upon each interaction with a pulse. Unlike the pulse durations, which can be taken to be similar to the experimental standard or state-of-the-art, the pulse amplitudes must often be tweaked until a good signal is achieved. If the power is too low or too high, it may be difficult to filter the desired third-order signal from the background of linear or higher orders. The window that is favorable for third-order detection may depend on how the signal is constructed. Given the different nature of the HD and FD schemes, and in particular, how lower and higher orders are suppressed, this gives rise to a different range for the pulse amplitude. To illustrate this, we fix the pulse width by setting σ = 10 fs and vary the pulse amplitude. The resulting total-correlation 2D spectra at zero-waiting time for HD and FD are shown in Figs. 5 and 7, respectively.

It should be noted that the first and third pulse interval are scanned in steps of 10 fs up to 300 fs, which is much longer than the 50 fs interval used for Fourier transforming the corresponding intervals in the double-sided Feynman diagram calculations. This discrepancy stems from the fact that the undersampling, once folded within the Nyquist frequency, appears to sample on a faster time scale. In our case, the folding factor is 3, which means the undersampling at 10 fs corresponds to a normal sampling at $10−3×2πℏωpulse=1.73$ fs. 30 steps of 1.73 fs add up to about 50 fs. It follows that the decay from |2⟩ to |1⟩ will be greater with the longer scanning time. The reduced sensitivity by sampling beyond 1.3 times the decay rate⁷¹ is counterbalanced by the increased frequency resolution and 300 fs is a good compromise.

By sweeping the pulse amplitude across the third-order regime (see Fig. 5), we gain insight into the behavior in the limits of linear and higher orders, as well as the intermediate regime. At the lowest pulse amplitude, linear artifacts appear in HD spectra as vertical streaks, which translate as noise in the excitation frequency for a specific detection frequency. For an experienced experimenter or theoretician observing the spectra, such streaks would be met with suspicion and faced with scrutiny before they would be interpreted as a result of the underlying physics of the sample. In other words, these linear artifacts are not likely to be misinterpreted.

On the opposite end of the third-order regime (see the rightmost spectrum of Fig. 5), we see a perfectly plausible HD2D spectrum. It is only by comparison to the intermediate regime, the middle spectra of the same figure, that it is clear that the high pulse amplitude is introducing higher-order contributions into the spectrum. This demonstrates the importance of a well-calibrated pulse amplitude.

The use of explicit time evolution allows us to track the populations during the pulse sequence. This approach allows us to confirm that for small amplitude (peak interaction energy ≲27 meV), the system only ever absorbs one quantum of energy, i.e., we are in the purely perturbative limit.¹ As the amplitude is increased, increased population of the doubly excited state (the |3⟩ state) indicates the system is absorbing two quanta of energy from a single pulse. This non-perturbative response is seen in the qualitative changes in the spectrum observed in Fig. 5 as higher amplitude but can be uniquely identified as the time-dependent populations are available in the simulation.

In the intermediate regime (see the middle spectra in Fig. 5), we observe slight changes as the amplitude is increased. It can therefore be hard to justify quantitative conclusions based on a single 2D spectrum. Particularly, the cross peaks can be sensitive to the pulse amplitude. These arise from pathways with incorrect time-ordering, which are more likely at short waiting times. The zero-waiting time therefore presents an additional challenge for the explicit simulation, as well as experiments carried out in the lab.

We also investigate the time evolution of the spectra as the waiting time is scanned in steps of 5 fs. Because the linear artifacts are clearly more pronounced in the horizontal axis, we pick the evolution of the ω₃ = E₁ horizontal line to investigate whether the erratic spectrum at zero-waiting time becomes smooth for longer waiting times. Figure 6 shows that there is a slight suppression of the linear artifacts as the waiting time is increased for the lower amplitude case, especially for the diagonal peak, but some noise still remains when compared to the intermediate amplitude regime.

The effect of changes in the pulse amplitude on FD spectra is, similarly to the HD case, slight within the third-order regime and strong at the limits of the regime. Interestingly, the way the spectra break down in these limits is in stark contrast to the HD case. At too low pulse amplitude (see the leftmost spectrum of Fig. 7), we do not observe linear artifacts as vertical streaks, but instead the spectral peaks have a higher degree of randomness to them. At high pulse amplitudes, (see the rightmost spectrum of Fig. 7), we see the lower diagonal peak gaining intensity and the upper diagonal peak losing intensity. The result is a spectrum, which looks reasonable, but can lead to completely wrong analysis.

We investigate the time evolution of the spectra by scanning the waiting time in steps of 5 fs. Figure 8 shows how the diagonal evolves for the low pulse amplitude case (corresponding to the leftmost spectrum in Fig. 7) and for an intermediate pulse amplitude (corresponding to the second from the left in Fig. 7). It is evident that when the waiting time is increased, the low-amplitude case clears up and gradually resembles the intermediate amplitude case as the pulse overlap between the second and third pulse becomes negligible.

Comparing the two detection methods, we note that the HD third-order signal is more robust against higher-order contributions as the power is increased, as evidenced by the spectra acquired with pulse intensities corresponding to peak interaction energies of 27 meV and 56 meV (see the middle spectra in Fig. 5). The manner of which the third-order signal breaks apart in the limits of the third-order regime is very different and stems from the disparate ways of how the third-order signal was constructed.

Note that amplitude affects the number of chromophores needed to phase the signal in HD2D. The comparison to FD2D is therefore not straightforward as only one chromophore is needed to compute its signal. HD tends to prefer a peak interaction energy upward of 9 meV and can go to quite strong laser powers without compromising the third-order signal, whereas FD is more prone to higher order effects, but still produces a clear signal at lower laser powers than HD.

Seeing how influential the pulse overlap can be when the waiting time is varied, we proceed to discuss the effect of the pulse duration, which will increase or decrease not only the overlaps between pulses 2 and 3 but potentially between all pulses.

For a fair comparison of the effect of the pulse duration, the pulse amplitude must be adjusted accordingly to ensure that the pulse power remains the same, thus keeping the population and coherence transfer rates on the same scale.

When the pulse duration is increased, the uncertainty in the excitation and detection frequencies also increases. To limit this effect, one can scan the respective pulse delays for longer. However, it is not possible to correct for the effects caused by pulse overlap. The more the pulses overlap, the stronger the signal from incorrect pulse ordering becomes, particularly in polarization-controlled variants of 2DES as shown by Paleček et al.⁵⁸ Care should therefore be taken when analyzing experiments with considerable pulse overlap, especially of the early time dynamics where pulses 2 and 3 are close together.

Figures 9 and 10 show the total-correlation 2D spectra at zero-waiting time using FD and HD, respectively. Interestingly, HD is more robust against changes in the pulse duration. FD is more easily smeared and also struggles as the duration becomes very short, although this may be of little experimental interest as such pulse durations are not realistic/possible. Additionally, the rotating wave approximation starts to break down and causes more complicated spectra for the shortest pulse as there are only a few optical cycles within one pulse. From a theoretical standpoint, it is interesting to determine why the two detection methods differ. In HD, the phase picked up by each chromophore is solely by virtue of its position, whereas the phases picked up in FD are contained in the waveform and may be distorted as the pulse duration becomes comparable to the period of the pulse frequency. In any case, the linear artifacts seen for FD at the lowest pulse duration disappear when the waiting time is increased, just as we observed for HD at low pulse amplitude. A possible explanation for the smearing of the FD spectrum at long pulse durations is that the error introduced by the construction of the discrete Fourier transform [Eq. (11)] grows with the uncertainty in the interaction times. It should be noted that the phase-cycling equations are derived in the impulsive limit.⁶⁴ It would be interesting to test this hypothesis by experimenting with more elaborate phase-cycling schemes or applying the phase-cycling method to a HD setup, but that is beyond the scope of this work.

In our modeling of FD and HD 2DES, we used pulses with variable amplitudes and durations. Within the semi-classical approximation, this approach can be considered a full model spectroscopy-wise. We observe effects related to the finite pulses, which cannot be accounted for with the commonly employed Feynman diagram method. Of particular interest is the fact that the optimal window for the pulse power is different for the two detection schemes. The underlying reason for this discrepancy is the disparate ways that the third-order signal is constructed from the raw signal data. It is also interesting from a theoretical point of view to observe the dissimilar behavior of the two methods as the limits of pulse amplitude and pulse duration are tested.

Computationally, the FD simulation is quite cheap because only one model system is required to create the 2D spectra, although 27 runs are required for phase-cycling. For HD, the number of simulation runs can be much higher, starting from hundreds but potentially requiring 10 000–100 000, depending on the “phasing” conditions (dephasing, timestep, scanning length, etc).

The investigation of the effect of the pulse amplitude and the pulse duration shows substantial and non-trivial changes within the third-order regime, which calls for great care whenever quantitative analyses are attempted. Conclusions from quantitative analysis should essentially be backed up by non-perturbative simulations to take pulse effects into account. This is even more important when probing short-time dynamics where increased pulse overlapping is detrimental to the selection of spectroscopic pathways.

On top of the complementary information found in FD2D spectra, much of the promise of FD2D is the fact that it can be contrasted against HD2D spectra to provide information that is not contained in either FD2D or HD2D. Therefore, it is imperative that all aspects of the experiment are well understood.

The authors acknowledge support of the Australian Research Council through Grant No. CE170100026. This research was undertaken with the assistance of resources from the National Computational Infrastructure, which is supported by the Australian Government.

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

By adopting the model from Ref. 47, we are able to verify that our calculations, using phase-cycling instead of phase-modulation, are correct. An alternative outcome could be that the two phasing schemes bring about differences in the spectra, but this turns out not to be the case as the resemblance is striking (see Fig. 11 and Fig. 4 in the original work). Apart from replacing phase-modulation with phase-cycling, all parameters and equations are exactly the same.

For completeness and as a starting point for the discussion on differences between explicitly simulated heterodyne-detected and fluorescence-detected 2D spectra, we also compute the HD spectra with the same parameters. However, in calculations of HD spectra, one must also decide on how many chromophores to include (and their positions), which will generally be iterated until convergence is achieved. This number can depend strongly on model parameters, such as dephasing and energy levels, and on the specific realizations of the pulse trains, such as the pulse durations, pulse amplitudes, the scanning times, and the number of scanning steps.

Intuitively, one might expect that the two detection schemes operate in the same field-strength range, meaning that the third-order signal starts to appear at the same laser-power threshold and that higher orders enter at similar intensities. However, the calculated spectra in Fig. 12 show that the pulse amplitude had to be increased sevenfold to achieve optimal third-order signals for the heterodyne-detected simulation.