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

Ultra-fast and multi-dimensional spectroscopy gives a powerful looking glass into the dynamics of molecular systems. In particular two-dimensional electronic spectroscopy (2DES) provides a probe of coherence and the flow of energy within quantum systems which is not possible with more conventional techniques. While heterodyne-detected (HD) 2DES is increasingly common, more recently fluorescence-detected (FD) 2DES offers new opportunities, including single-molecule experiments. However in both techniques it can be difficult to unambiguously identify the pathways which dominate the signal. Therefore the use of numerically modelling of 2DES is vitally important, which in turn requires approximating the pulsing scheme to some degree. Here we employ non-pertubative time evolution to investigate the effects of finite pulse width and amplitude on 2DES signals. In doing so we identify key differences in the response of HD and FD detection schemes, as well as the regions of parameter space where the signal is obscured by unwanted artefacts in either technique. Mapping out parameter space in this way provides a guide to choosing experimental conditions and also shows in which limits the usual theoretical approximations work well and which limits more sophisticated approaches are required.


Introduction
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][2][3][4][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. 6 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][2][3][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][9][10][11] 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, [12][13][14][15] but it is now establishing itself as an incisive tool with real applications. [16][17][18][19] In the following we focus on comparing HD2D and FD2D; for a broader perspective on multidimensional ultrafast spectroscopy, references 7, 20 and 21 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 idealised experiments, we chose to go beyond the double-sided Feynman diagrams 5,22 and simulate the spectra non-perturbatively. 23,24 This approach enables more aspects of the experiments to be covered, [25][26][27][28] such as the pulse duration and the pulse amplitude, which are another key focus of our study. However, other sources of divergence bebtween 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.

Background
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 contributios 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, are 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 -ion emission 29,30 or photocurrent, 31 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 polarisation has a number of potential benefits. The most known advantage is that signals from small volumes, in principle single molecules, 16,32 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. 17 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. 19 Recently, Maly and Mancal suggested that the acquisition time, which has no analogue in HD2D, could be varied in order to isolate specific contributions to the total signal, e.g. the exciton-exciton annihilation. 33 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. 34 However, careful analysis can differentiate between true nonlinear signals and incoherently mixed linear signals. 35 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 which have implemented the fluorescence-detection setup are relatively few, and don't enjoy the same level of experience and literature to draw upon. However, the increasing rate of articles published suggest 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. 33,[35][36][37][38] 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. 36 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 analyse 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 only one, we chose to scan both in our simulations. Also, 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. 39 Both of the abovementioned 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. 40 We also sample uniformly, opting for simplicity rather than trying to reduce the computational cost with non-uniform sampling. 41, 42 We do not investigate the polarisation dependency, which is sometimes exploited to select specific pathways, [43][44][45][46][47] but it is known that the early-time dynamics suffer from incorrect pulseordering artifacts. 48 We assume non-interacting chromophores in our calculations, although the delocalisation of excitation energy in many cases is unexpectedly long-range 49 and its role in energy transfer is a hot topic in the field. 18, 50, 51

Theory
First we introduce the theoretical and computational methods used throughout the 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. Lastly, we explain the model system which we use in our simulations.

The quantum dynamics of a system interacting with an electric field
The vast majority of 2D spectroscopy models employ the semiclassical approximation, where the state of the system is described quantum mechanically, but the field, E(r, t), is described classically: Here, A n (t − t n ) describes the envelope of the nth pulse centred 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 timedependent Schrodinger equation where H 0 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 can not 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). 52 where an off-diagonal L j = a † m a n operator represents spontaneous relaxation or excitation, a diagonal L j operator represents a dephasing process, and N denote the number of different decoherence channels. The formal solution to the differential equation in (4) is found by integrating both sides of the equation: ρ(t) = e Lt ρ(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 towards single-molecule 2DES experiments, exploiting the high sensitivity of fluorescence detectors.

Heterodyne-Detected 2DES
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 1 ± k 2 + k 3 , where the upper combination corresponds to the rephasing signal and the lower combination to the nonrephasing signal. A third (possible) phase matching direction, k DQC = k 1 + k 2 − k 3 , 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 polarisation in the sample. While this enables high signal-to-noise ratios, it also puts a lower bound on the spatial resolution: λ 3 , 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 polarisation of the individual molecules will cause the third-order signal to vanish due to the random phases, unless the polarisation 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 which stems from a single interaction with each of the three pulses. This reduces the number of phase-matching Figure 1. (a) Heterodyne-detection of two-dimensional electronic spectra employs a local oscillator (LO) to read out the macroscopic polarisation generated by three interactions with the electric field. In the non-collinear geometry, linearly independent wave vectors imprint unique phases with each photo-excitation and de-excitation. The selection of specific third-order processes, such as rephasing and nonrephasing, is achieved through phase-matching, in which the wave vector of the LO is a combination of the previous wave vectors, but with the necessary signs flipped to reflect the order of (de-)excitation for the respective processes. (b) Fluorescence-detection of two-dimensional electronic spectra is performed in a collinear setup. Instead of measuring the polarisation, fluorescence is detected up to a cut-off acquisition-time. The integrated fluorescence as a function of the pulse delays constitutes the signal. However, the selection of third-order processes requires each pulse train to be executed 27 times with phases cycled. Subsequent Fourier transforms with respect to the distinct phase evolutions of rephasing and nonrephasing processes, produce the 2D spectra which are related to, but not equivalent to, the heterodyne-detected 2D spectra. directions to three: −k 1 + k 2 + k 3 , k 1 − k 2 + k 3 and k 1 + k 2 − k 3 . A schematic of the phase matching method is shown in figure 1(a).
To achieve the phase-matching in a non-perturbative calculation, the dynamics of each individual molecule j must be computed separately using equation 4. Here, we assume that the individual molecules evolve independently and that each molecule is initially in its ground state. For simplicity, we employ identical Hamiltonian and Lindblad operators on each individual quantum system. A generalisation 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 polarisation in the chosen phase-matching direction is given by: 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.

Fluorescence-Detected 2DES
Instead of detecting the third-order polarisation 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, 12 which will not discussed here, and phase cycling. 53

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 is enough to extract the third-order signals. 54 Figure 1(b) is a schematic of the FD2D experiment emphasising 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.

3.3.2.
Example with excited state population 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. 35 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. 37 It is, however, relevant for other action spectroscopies. [29][30][31] The excited state population after interactions with four pulses with phases φ 1−4 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, ρ 0 = |g g|, and the final state is diagonal, the following condition must be fulfilled: 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 byp(β = 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.
p(τ, T, t, β, γ, δ) = 1 (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 1 , t 2 , t 3 , t 4 }. 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 equation 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.

3.3.3.
Using the fluorescence signal Instead of using the final populations as proxies 35 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 undergoes the same phase cycling process as for the population data to construct the rephasing and nonrephasing spectra.

The Diagrammatic Approach to 2DES
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 doublesided 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, 55 as each propagator G(t) ≡ e Lt is independent of other time variables. This obviates the cumbersome process of numerically integrating the state of the system for all realisations of the pulse intervals. More detailed accounts on the subject can be found elsewhere, but we include three diagrams which highlight the difference between HD and FD. Figure 2 shows excited state absorption diagrams which evolve identically until the detection event as is also evident from their corresponding equations: Here, − → V ij denotes the transition dipole moment operator acting on the right, causing a transition from j to i: The rightmost diagram, however, is unique to FD, and is denoted ESA2. Owing to the opposite sign acquired in the last pulse interaction, the middle and the right diagrams largely cancel, which becomes the main source of discrepancy between HD and FD spectra.
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, equation (19), which is equivalent to the HD diagram, equation (18), and the higher excited state, equation (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. Doublequantum 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.

Model
In figure 3 we depict our model system, which we adopted from reference 37. 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 reference 37 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 reference 37: where E 1 and E 2 correspond to the B850 and B800 absorption bands of the inner and outer rings of bacteriochlorophylls in the light-harvesting complex of purple bacteria. [56][57][58][59][60] The optical transitions between the energy levels are all allowed and are equally strong: eE 0 µ ij = 8 meV, with the exception of the transition from ground to the highest excited state, which is forbidden. Electronic relaxation is modelled by the jump operators, the off-diagonal Lindblad operators: L 10 = |0 1|, L 21 = |1 2| and L 32 = |2 3|, which are scaled by Γ 10 = 4.13 µeV, Γ 21 = 4.13 meV, Γ 32 = 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 00 = |0 0|, L 11 = |1 1|, L 22 = |2 2| and L 33 = |3 3|, which are all scaled by the same dephasing strength Γ Dephasing = 41.3 meV.
The model system is initialised with the ground state fully populated and then propagated according to equation 4 with the appropriate (light-matter) interaction Hamiltonian until it is detected as described by equation 6 or 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.

Results
By applying the theory of HD and FD 2DES on the model system, we are able to nonperturbatively 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 nonidealised pulse, we compute the corresponding spectra using the double-sided Feynman diagrams laid out in ref. 37, see figure 4.  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 1 = 50 fs) and last (t 3 = 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.
Before we discuss the response to changes in the pulse amplitude, we note two general differences between FD and HD that we observe in figure 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 favours 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 eE 0 ij = 81 meV eE 0 ij = 56 meV eE 0 ij = 27 meV eE 0 ij = 9 meV Figure 5. The real parts of the zero-waiting time total 2D spectra using heterodynedetection. The insets indicate the peak interaction energies, which were chosen to illustrate failure at low pulse amplitude (linear artifacts), at high pulse amplitude (higher-order effects) and the intermediate regime. The dotted lines guide the eye to the model system energies, the diagonal and antidiagonal.
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-idealised realisations, i.e. real-life experiments and simulations with finite pulses.
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., 37 but with phase-cycling replacing phase-modulation. As can be seen in the appendix, we were able to produce the same spectra with our phasecycling method using the same parameters, specifically a peak interaction energy of 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.

Pulse Amplitude
One of the challenges that arises 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 favourable 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 figures 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 calcaulations. 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 adds 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 61 is counterbalanced by the increased frequency resolution and 300 fs is a good compromise.
4.1.1. The effect on HD spectra By sweeping the pulse amplitude across the thirdorder regime, see figure 5, we gain insight into the behaviour 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 translates 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 figure 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.
In the intermediate regime, see the middle spectra in figure 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 The linear artifacts are strongest for the shortest waiting times such as 0 and 5 fs, with particularly the diagonal peak recovering as the waiting time is increased, but the effect is much less pronounced than for the FD case. Each pair of curves are normalised to the same maximum and are offset for visual clarity.
in steps of 5 fs. Because the linear artifacts are clearly more pronounced in the horizontal axis, we pick the evolution of the ω 3 = E 1 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 on FD spectra
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 figure 7, we do not observe linear artifacts as vertical streaks but instead the spectral eE 0 ij = 27 meV eE 0 ij = 9 meV eE 0 ij = 3 meV eE 0 ij = 1 meV Figure 7. The real parts of the zero-waiting time total 2D spectra using fluorescencedetection. The insets indicate the peak interaction energies, which were chosen to illustrate failure at low pulse amplitude (linear artifacts), at high pulse amplitude (higher-order effects) and the intermediate regime. The dotted lines guide the eye to the model system energies, the diagonal and antidiagonal.
peaks have a higher degree of randomness to them. At high pulse amplitudes, see the rightmost spectrum of figure 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 figure 7) and for an intermediate pulse amplitude (corresponding to the second from the left in figure 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 and 56 meV, see the middle spectra in figure 5. The manner of which the thirdorder 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 upwards 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.

Pulse Durations
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 = 20 = 10 = 5 = 2.5 Figure 9. The real parts of the zero-waiting time total 2D spectra using heterodynedetection. The insets indicate the standard deviation of the Gaussian in fs, and hence the duration of the pulses. The σ values range from 2.5 fs to 20 fs, which is spanning what is currently not possible to what is easily achieved in most 2DES setups. The dotted lines guide the eye to the model system energies, the diagonal and antidiagonal.
ordering becomes, particularly in polarisation-controlled variants of 2DES as shown by Paleček et al. 48 Care should therefore be taken when analysing experiments with considerable pulse overlap, especially of the early-time dynamics where pulse 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. 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.

Conclusions
In our modelling 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 can not 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 behaviour of the two methods as the limits of pulse amplitude and pulse duration are tested. = 20 = 10 = 5 = 2.5 Figure 10.
The real parts of the zero-waiting time total 2D spectra using fluorescence-detection. The insets indicate the standard deviation of the Gaussian in fs, and hence the duration of the pulses. The σ values range from 2.5 fs to 20 fs, which is spanning what is currently not possible to what is easily achieved in most 2DES setups. The dotted lines guide the eye to the model system energies, the diagonal and antidiagonal.
Computationally, the FD simulation is quite cheap because only 1 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.

Acknowledgements
The authors acknowledge support of the Australian Research Council through grant CE170100026. This research was undertaken with the assistance of resources from the National Computational Infrastructure, which is supported by the Australian Government.  Figure A1. A replica of the zero-waiting time fluorescence-detected 2D spectra in ref. 37. The only difference is that here phase-cycling is employed, instead of phasemodulation, as a method to filter the desired third-order signals from the background. Rel ij denotes the relaxation from |i to |j , where Rel 10 is of primary relevance to fluorescence detection. Note that the sign is deliberately changed to promote the similarity to heterodyne-detected spectra, and that the axes are switched with respect to ref. 37