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 that dominate the signal. Therefore, the use of numerically modeling of 2DES is vitally important, which, in turn, requires approximating the pulsing scheme to some degree. Here, we employ non-perturbative 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 artifacts 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 in which limits more sophisticated approaches are required.

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.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–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 emission39,40 or photocurrent,41 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.18 

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

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.43 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.44 However, careful analysis can differentiate between true nonlinear signals and incoherently mixed linear signals.45 

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.46 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.49 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.50 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.58 

We assume non-interacting chromophores in our calculations, although the delocalization of excitation energy in many cases is unexpectedly long-range59 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,

E(r,t)=nAn(ttn)exp(iω(ttn)+iknriϕn).
(1)

Here, An(ttn) describes the envelope of the nth pulse centered at tn, ω is the frequency of the field, kn 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

Hint(t)=μE(t),
(2)

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

it|Ψ(r,t)=[H0+Hint(t)]|Ψ(r,t)H(t)|Ψ(r,t),
(3)

where H0 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),62 

ddtρ(t)=i[ρ(t),H(t)]+j=1NΓjLjρ(t)Lj12[ρ(t)LjLj+LjLjρ(t)]L[ρ(t)],
(4)

where an off-diagonal Lj=aman operator represents spontaneous relaxation or excitation, a diagonal Lj 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 ksignal = ∓k1 ±k2 + k3, where the upper combination corresponds to the rephasing signal and the lower combination corresponds to the nonrephasing signal. A third (possible) phase–matching direction, kDQC = k1 + k2k3, 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: ≳λ3, where λ is the radiation wavelength.

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

Ej(t)=E0n=13expttn2σ2exp(iωtiknrj),
(5)

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(−ikp · rj) 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, ksignal, 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: −k1 + k2 + k3, k1k2 + k3, and k1 + k2k3. A schematic of the phase-matching method is shown in Fig. 1(a).

FIG. 1.

(a) Heterodyne-detection of two-dimensional electronic spectra employs a local oscillator (LO) to read out the macroscopic polarization 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 polarization, fluorescence is detected up to a cutoff 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.

FIG. 1.

(a) Heterodyne-detection of two-dimensional electronic spectra employs a local oscillator (LO) to read out the macroscopic polarization 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 polarization, fluorescence is detected up to a cutoff 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.

Close modal

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

Pksignal,τ,T,t=2Rejexpiksignalrjα<βμαβ(j)ρβα(j)Ej;t.
(6)

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,

SHD2D(ksignal,ωτ,T,ωt)=i0dτexp(±iωττ)×0dtexp(iωtt)P(ksignal,τ,T,t),
(7)

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,13 which will not discussed here, and phase cycling.63 

1. 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.64Figure 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.

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.45 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.47 It is, however, relevant for other action spectroscopies.39–41 

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,

p(τ,T,t,ϕ1,ϕ2,ϕ3,ϕ4)=α,β,γ,δp̃(τ,T,t,α,β,γ,δ)×exp[i(αϕ1+βϕ2+γϕ3+δϕ4)].
(8)

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: α + β + γ + δ = 0. Consequently, α, β, γ, δ are dependent, which allows us to define ϕ21ϕ2ϕ1, ϕ31ϕ3ϕ1, and ϕ41ϕ4ϕ1. Equation (8) then becomes

p(τ,T,t,ϕ21,ϕ31,ϕ41)=β,γ,δp̃(τ,T,t,β,γ,δ)×exp[i(βϕ21+γϕ31+δϕ41)].
(9)

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,

p̃(τ,T,t,β,γ,δ)=1(2π)302π02π02πdϕ41dϕ31dϕ21×p(τ,T,t,ϕ21,ϕ31,ϕ41)eiβϕ21eiγϕ31eiδϕ41.
(10)

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

p̃(τ,T,t,β,γ,δ)=1LMNn=0N1m=0M1l=0L1×p(τ,T,t,lΔϕ21,mΔϕ31,nΔϕ41)×eilβΔϕ21eimγΔϕ31einδΔϕ41.
(11)

As shown by Tan,64 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

E  lmn(t)=E1(t)+E2  l(t)+E3  m(t)+E4  n(t),
(12)
E1(t)=E0exptt12σ22exp(iωt),
E2  l(t)=E0exptt22σ22exp(iωt+il2π3),E3  m(t)=E0exptt32σ22exp(iωt+im2π3),
(13)
E4  n(t)=E0exptt42σ22exp(iωt+in2π3),

where l, m, n are cycled between 0, 1, and 2 and the resulting populations from the 27 combinations added for each set of {t1, t2, t3, t4}. Note that the kn vector is dropped as the kn · 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

SFD2D(β,γ,δ,ωτ,T,ωt)=i0dτexp(±iωττ)0dt×exp(iωtt)p̃(β,γ,δ,τ,T,t)
(14)

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. Using the fluorescence signal

Instead of using the final populations as proxies45 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

Rel10=0tacqdtΓ10Tr[L10ρS(t)L10],
(15)
L10=a0a1,
(16)

and similarly for relaxations between other states. tacq 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,65 

G±(ω,τf)0τfe±iωtG(t)dt=[±iω1+L]1(e±iωτ1+Lτ1),
(17)

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,

ESA1HD(ω1,t2,ω3)=TrV03G+(ω3)V32G(t2)×V20G(ω1)V01|00|,
(18)
ESA1FD(ω1,t2,ω3;tacq)=0tacqdtγ1Tr|11|G(t)V13G+(ω3)×V32G(t2)V20G(ω1)V01|00|,
(19)
ESA2FD(ω1,t2,ω3;tacq)=0tacqdtγ3Tr|33|G(t)V13G+(ω3)×V32G(t2)V20G(ω1)V01|00|.
(20)
FIG. 2.

Double-sided Feynman diagrams of an identical excited state absorption process but with different read-outs due to the choice of detection method. Traditional HD (left) interferes a local oscillator with the final coherence to produce a signal. FD, on the other hand, interacts with a fourth pulse, which can result in a population of a lower excited state (middle) or a higher excited state (right), from which emission of fluorescence is continually detected. The middle diagram shares the same phase evolution as the HD diagram to the left (both are ESA1) as indeed every HD diagram has an equivalent FD diagram. 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.

FIG. 2.

Double-sided Feynman diagrams of an identical excited state absorption process but with different read-outs due to the choice of detection method. Traditional HD (left) interferes a local oscillator with the final coherence to produce a signal. FD, on the other hand, interacts with a fourth pulse, which can result in a population of a lower excited state (middle) or a higher excited state (right), from which emission of fluorescence is continually detected. The middle diagram shares the same phase evolution as the HD diagram to the left (both are ESA1) as indeed every HD diagram has an equivalent FD diagram. 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.

Close modal

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 tacq 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,

Egr=0,E1=1.46  eV,E2=1.55  eV,E3=E1+E2,
(21)

where Egr is the ground state energy, E1 and E2 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 E3 is the energy of the doubly excited exciton. The optical transitions between the energy levels are all allowed and are equally strong: eE0μ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: L10 = |0⟩⟨1|, L21 = |1⟩⟨2|, and L32 = |2⟩⟨3|, which are scaled by Γ10 = 4.13 μeV, Γ21 = 4.13 meV, and Γ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: L00 = |0⟩⟨0|, L11 = |1⟩⟨1|, L22 = |2⟩⟨2|, and L33 = |3⟩⟨3|, which are all scaled by the same dephasing strength ΓDephasing = 41.3 meV.

FIG. 3.

The model system has four energy levels where all transitions are optically allowed and equally strong, apart from |0⟩ → |3⟩, which is forbidden. The relaxation rates τij indicate an ultrafast relaxation pathway from the highest-excited state and a slightly slower relaxation from the second excited state. The decay from the lowest-excited state to the ground state is orders of magnitude slower. This model system is identical to the one in Ref. 47.

FIG. 3.

The model system has four energy levels where all transitions are optically allowed and equally strong, apart from |0⟩ → |3⟩, which is forbidden. The relaxation rates τij indicate an ultrafast relaxation pathway from the highest-excited state and a slightly slower relaxation from the second excited state. The decay from the lowest-excited state to the ground state is orders of magnitude slower. This model system is identical to the one in Ref. 47.

Close modal

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.

FIG. 4.

The real parts of the zero-waiting time total 2D spectra, calculated with the traditional double-sided Feynman diagrams derived from perturbation theory. (a) The spectrum is simulated with pathways selected by the fluorescence-detection setup, and (b) the heterodyne-detected counterpart is shown on the right. Note that the Fourier transform of the coherence time t1 and signal time t3 is taken for a time interval of 50 fs in order to achieve similar broadening as for the non-perturbative simulation due to the finite Fourier transform.

FIG. 4.

The real parts of the zero-waiting time total 2D spectra, calculated with the traditional double-sided Feynman diagrams derived from perturbation theory. (a) The spectrum is simulated with pathways selected by the fluorescence-detection setup, and (b) the heterodyne-detected counterpart is shown on the right. Note that the Fourier transform of the coherence time t1 and signal time t3 is taken for a time interval of 50 fs in order to achieve similar broadening as for the non-perturbative simulation due to the finite Fourier transform.

Close modal

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 (t1 = 50 fs) and last (t3 = 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.47 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.

FIG. 5.

The real parts of the zero-waiting time total 2D spectra using heterodyne-detection. 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 at the intermediate regime. The dotted lines guide the eye to the model system energies, the diagonal and antidiagonal.

FIG. 5.

The real parts of the zero-waiting time total 2D spectra using heterodyne-detection. 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 at the intermediate regime. The dotted lines guide the eye to the model system energies, the diagonal and antidiagonal.

Close modal

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 103×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 rate71 is counterbalanced by the increased frequency resolution and 300 fs is a good compromise.

1. The effect on HD spectra

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.1 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 ω3 = E1 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.

FIG. 6.

The intensity along the horizontal line ω3 = E1 of a HD2D spectrum is plotted as a function of the waiting time in steps of 5 fs, with 0 fs waiting time at the top and 50 fs at the bottom. The full line is calculated using an intermediate pulse amplitude, corresponding to a peak interaction energy of 27 meV, and the dotted line is calculated using a low pulse amplitude with a peak interaction energy of 9 meV. 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 is normalized to the same maximum and is offset for visual clarity.

FIG. 6.

The intensity along the horizontal line ω3 = E1 of a HD2D spectrum is plotted as a function of the waiting time in steps of 5 fs, with 0 fs waiting time at the top and 50 fs at the bottom. The full line is calculated using an intermediate pulse amplitude, corresponding to a peak interaction energy of 27 meV, and the dotted line is calculated using a low pulse amplitude with a peak interaction energy of 9 meV. 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 is normalized to the same maximum and is offset for visual clarity.

Close modal

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

FIG. 7.

The real parts of the zero-waiting time total 2D spectra using fluorescence-detection. 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.

FIG. 7.

The real parts of the zero-waiting time total 2D spectra using fluorescence-detection. 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.

Close modal

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.

FIG. 8.

The intensity on the diagonal line ω1 = ω3 of a FD2D spectrum is plotted as a function of the waiting time in steps of 5 fs, with 0 fs waiting time at the top and 50 fs at the bottom. The full line is calculated using an intermediate pulse amplitude with a peak interaction energy of 3 meV, and the dotted line is calculated using a low pulse amplitude corresponding to a peak interaction energy of 1 meV. When the waiting time becomes large compared to the pulse width, the noise vanishes as a result of a diminishing contribution from incorrect time-ordering pathways. Each pair of curves is normalized to the same maximum and is offset for visual clarity.

FIG. 8.

The intensity on the diagonal line ω1 = ω3 of a FD2D spectrum is plotted as a function of the waiting time in steps of 5 fs, with 0 fs waiting time at the top and 50 fs at the bottom. The full line is calculated using an intermediate pulse amplitude with a peak interaction energy of 3 meV, and the dotted line is calculated using a low pulse amplitude corresponding to a peak interaction energy of 1 meV. When the waiting time becomes large compared to the pulse width, the noise vanishes as a result of a diminishing contribution from incorrect time-ordering pathways. Each pair of curves is normalized to the same maximum and is offset for visual clarity.

Close modal

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.58 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.64 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.

FIG. 9.

The real parts of the zero-waiting time total 2D spectra using heterodyne-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 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.

FIG. 9.

The real parts of the zero-waiting time total 2D spectra using heterodyne-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 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.

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

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

Close modal

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.

FIG. 11.

A replica of the zero-waiting time fluorescence-detected 2D spectra in Ref. 47. The only difference is that here phase-cycling is employed, instead of phase-modulation, as a method to filter the desired third-order signals from the background. Relij denotes the relaxation from |i⟩ to | j⟩, where Rel10 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. 47.

FIG. 11.

A replica of the zero-waiting time fluorescence-detected 2D spectra in Ref. 47. The only difference is that here phase-cycling is employed, instead of phase-modulation, as a method to filter the desired third-order signals from the background. Relij denotes the relaxation from |i⟩ to | j⟩, where Rel10 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. 47.

Close modal

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.

FIG. 12.

Zero-waiting time heterodyne-detected 2D spectra. Left: The rephasing and nonrephasing signals were calculated using the same parameters as in Ref. 47. Right: The same parameters were used, apart from the pulse amplitude which was increased from 8 to 56 meV in peak interaction energy. 10 000 randomized positions were employed to converge the phasing of the third-order signal.

FIG. 12.

Zero-waiting time heterodyne-detected 2D spectra. Left: The rephasing and nonrephasing signals were calculated using the same parameters as in Ref. 47. Right: The same parameters were used, apart from the pulse amplitude which was increased from 8 to 56 meV in peak interaction energy. 10 000 randomized positions were employed to converge the phasing of the third-order signal.

Close modal
1.
J.
Yuen-Zhou
,
J. J.
Krich
,
I.
Kassal
,
A. S.
Johnson
, and
A.
Aspuru-Guzik
,
Ultrafast Spectroscopy
(
IOP Publishing
,
2014
), pp.
2053
2563
.
2.
D. M.
Jonas
, “
Two-dimensional femtosecond spectroscopy
,”
Annu. Rev. Phys. Chem.
54
(
1
),
425
463
(
2003
).
3.
A. M.
Branczyk
,
D. B.
Turner
, and
G. D.
Scholes
, “
Crossing disciplines—A view on two-dimensional optical spectroscopy
,”
Ann. Phys.
526
(
1-2
),
31
49
(
2014
).
4.
A.
Gelzinis
,
R.
Augulis
,
V.
Butkus
,
B.
Robert
, and
L.
Valkunas
, “
Two-dimensional spectroscopy for non-specialists
,”
Biochim. Biophys. Acta, Bioenerg.
1860
(
4
),
271
285
(
2019
).
5.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
, Optical and Imaging Sciences Series (
Oxford University Press
,
1995
).
6.
J. A.
Davis
,
C. R.
Hall
,
L. V.
Dao
,
K. A.
Nugent
,
H. M.
Quiney
,
H. H.
Tan
, and
C.
Jagadish
, “
Three-dimensional electronic spectroscopy of excitons in asymmetric double quantum wells
,”
J. Chem. Phys.
135
(
4
),
044510
(
2011
).
7.
P.
Hamm
and
M.
Zanni
,
Concepts and Methods of 2D Infrared Spectroscopy
(
Cambridge University Press
,
2011
).
8.
G. S.
Engel
,
T. R.
Calhoun
,
E. L.
Read
,
T.-K.
Ahn
,
T.
Mančal
,
Y.-C.
Cheng
,
R. E.
Blankenship
, and
G. R.
Fleming
, “
Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems
,”
Nature
446
,
782
(
2007
).
9.
M. B.
Plenio
,
J.
Almeida
, and
S. F.
Huelga
, “
Origin of long-lived oscillations in 2D-spectra of a quantum vibronic model: Electronic versus vibrational coherence
,”
J. Chem. Phys.
139
(
23
),
235102
(
2013
).
10.
H.-G.
Duan
,
P.
Nalbach
,
V. I.
Prokhorenko
,
S.
Mukamel
, and
M.
Thorwart
, “
On the origin of oscillations in two-dimensional spectra of excitonically-coupled molecular systems
,”
New J. Phys.
17
(
7
),
072002
(
2015
).
11.
H.-G.
Duan
,
V. I.
Prokhorenko
,
R. J.
Cogdell
,
K.
Ashraf
,
A. L.
Stevens
,
M.
Thorwart
, and
R. J. D.
Miller
, “
Nature does not rely on long-lived electronic quantum coherence for photosynthetic energy transfer
,”
Proc. Natl. Acad. Sci. U. S. A.
114
(
32
),
8493
8498
(
2017
).
12.
J.
Lim
,
C. M.
Bösen
,
A. D.
Somoza
,
C. P.
Koch
,
M. B.
Plenio
, and
S. F.
Huelga
, “
Multicolor quantum control for suppressing ground state coherences in two-dimensional electronic spectroscopy
,”
Phys. Rev. Lett.
123
,
233201
(
2019
).
13.
P. F.
Tekavec
,
G. A.
Lott
, and
A. H.
Marcus
, “
Fluorescence-detected two-dimensional electronic coherence spectroscopy by acousto-optic phase modulation
,”
J. Chem. Phys.
127
(
21
),
214307
(
2007
).
14.
A. K.
De
,
D.
Monahan
,
J. M.
Dawlaty
, and
G. R.
Fleming
, “
Two-dimensional fluorescence-detected coherent spectroscopy with absolute phasing by confocal imaging of a dynamic grating and 27-step phase-cycling
,”
J. Chem. Phys.
140
(
19
),
194201
(
2014
).
15.
S.
Draeger
,
S.
Roeding
, and
T.
Brixner
, “
Rapid-scan coherent 2D fluorescence spectroscopy
,”
Opt. Express
25
(
4
),
3259
3267
(
2017
).
16.
S.
Goetz
,
D.
Li
,
V.
Kolb
,
J.
Pflaum
, and
T.
Brixner
, “
Coherent two-dimensional fluorescence micro-spectroscopy
,”
Opt. Express
26
(
4
),
3915
3925
(
2018
).
17.
V.
Tiwari
,
Y. A.
Matutes
,
A.
Konar
,
Z.
Yu
,
M.
Ptaszek
,
D. F.
Bocian
,
D.
Holten
,
C.
Kirmaier
, and
J. P.
Ogilvie
, “
Strongly coupled bacteriochlorin dyad studied using phase-modulated fluorescence-detected two-dimensional electronic spectroscopy
,”
Opt. Express
26
(
17
),
22327
22341
(
2018
).
18.
V.
Tiwari
,
Y.
Acosta Matutes
,
A. T.
Gardiner
,
T. L. C.
Jansen
,
R. J.
Cogdell
, and
J. P.
Ogilvie
, “
Spatially-resolved fluorescence-detected two-dimensional electronic spectroscopy probes varying excitonic structure in photosynthetic bacteria
,”
Nat. Commun.
9
(
1
),
4219
(
2018
).
19.
M.
Schröter
,
T.
Pullerits
, and
O.
Kühn
, “
Using fluorescence detected two-dimensional spectroscopy to investigate initial exciton delocalization between coupled chromophores
,”
J. Chem. Phys.
149
(
11
),
114107
(
2018
).
20.
K. J.
Karki
,
J.
Chen
,
A.
Sakurai
,
Q.
Shi
,
A. T.
Gardiner
,
O.
Kühn
,
R. J.
Cogdell
, and
T.
Pullerits
, “
Before Förster. Initial excitation in photosynthetic light harvesting
,”
Chem. Sci.
10
,
7923
7928
(
2019
).
21.
T. A. A.
Oliver
, “
Recent advances in multidimensional ultrafast spectroscopy
,”
R. Soc. Open Sci.
5
(
1
),
171425
(
2018
).
22.
Y.
Song
,
X.
Li
, and
J. P.
Ogilvie
, “
Two-dimensional electronic spectroscopy
,” in
Encyclopedia of Modern Optics
, 2nd ed., edited by
B. D.
Guenther
and
D. G.
Steel
(
Elsevier
,
Oxford
,
2018
), pp.
184
196
.
23.
J.
Lim
,
D. J.
Ing
,
J.
Rosskopf
,
J.
Jeske
,
J. H.
Cole
,
S. F.
Huelga
, and
M. B.
Plenio
, “
Signatures of spatially correlated noise and non-secular effects in two-dimensional electronic spectroscopy
,”
J. Chem. Phys.
146
(
2
),
024109
(
2017
).
24.
P. A.
Rose
and
J. J.
Krich
, “
Automatic Feynman diagram generation for nonlinear optical spectroscopies
,”
J. Chem. Phys.
154
,
034109
(
2021
).
25.
B.
Brüggemann
,
P.
Kjellberg
, and
T.
Pullerits
, “
Non-perturbative calculation of 2D spectra in heterogeneous systems: Exciton relaxation in the FMO complex
,”
Chem. Phys. Lett.
444
(
1
),
192
196
(
2007
).
26.
B.
Brüggemann
and
T.
Pullerits
, “
Nonperturbative modeling of fifth-order coherent multidimensional spectroscopy in light harvesting antennas
,”
New J. Phys.
13
(
2
),
025024
(
2011
).
27.
P. A.
Rose
and
J. J.
Krich
, “
Numerical method for nonlinear optical spectroscopies: Ultrafast ultrafast spectroscopy
,”
J. Chem. Phys.
150
,
214105
(
2019
).
28.
P. A.
Rose
and
J. J.
Krich
, “
Efficient numerical method for predicting nonlinear optical spectroscopies of open systems
,”
J. Chem. Phys.
154
,
034108
(
2021
).
29.
M.
Richter
,
R.
Singh
,
M.
Siemens
, and
T. S.
Cundiff
, “
Deconvolution of optical multidimensional coherent spectra
,”
Sci. Adv.
4
(
6
),
eaar7697
(
2018
).
30.
C. L.
Smallwood
,
T. M.
Autry
, and
S. T.
Cundiff
, “
Analytical solutions to the finite-pulse Bloch model for multidimensional coherent spectroscopy
,”
J. Opt. Soc. Am. B
34
(
2
),
419
429
(
2017
).
31.
V.
Perlík
,
J.
Hauer
, and
F.
Šanda
, “
Finite pulse effects in single and double quantum spectroscopies
,”
J. Opt. Soc. Am. B
34
(
2
),
430
439
(
2017
).
32.
T. N.
Do
,
M. F.
Gelin
, and
H.-S.
Tan
, “
Simplified expressions that incorporate finite pulse effects into coherent two-dimensional optical spectra
,”
J. Chem. Phys.
147
(
14
),
144103
(
2017
).
33.
X.
Leng
,
S.
Yue
,
Y.-X.
Weng
,
K.
Song
, and
Q.
Shi
, “
Effects of finite laser pulse width on two-dimensional electronic spectroscopy
,”
Chem. Phys. Lett.
667
,
79
86
(
2017
).
34.
H.
Li
,
A. P.
Spencer
,
A.
Kortyna
,
G.
Moody
,
D. M.
Jonas
, and
S. T.
Cundiff
, “
Pulse propagation effects in optical 2D Fourier-transform spectroscopy: Experiment
,”
J. Phys. Chem. A
117
(
29
),
6279
6287
(
2013
).
35.
P. F.
Tekavec
,
J. A.
Myers
,
K. L. M.
Lewis
,
F. D.
Fuller
, and
J. P.
Ogilvie
, “
Effects of chirp on two-dimensional Fourier transform electronic spectra
,”
Opt. Express
18
(
11
),
11015
11024
(
2010
).
36.
M. K.
Yetzbacher
,
N.
Belabas
,
K. A.
Kitney
, and
D. M.
Jonas
, “
Propagation, beam geometry, and detection distortions of peak shapes in two-dimensional Fourier transform spectra
,”
J. Chem. Phys.
126
(
4
),
044511
(
2007
).
37.
N.
Belabas
and
D. M.
Jonas
, “
Fourier algorithm for four-wave-mixing signals from optically dense systems with memory
,”
Opt. Lett.
29
(
15
),
1811
1813
(
2004
).
38.
S. M.
Gallagher Faeder
and
D. M.
Jonas
, “
Two-dimensional electronic correlation and relaxation spectra: Theory and model calculations
,”
J. Phys. Chem. A
103
(
49
),
10489
10505
(
1999
).
39.
L.
Bruder
,
M.
Binz
, and
F.
Stienkemeier
, “
Efficient isolation of multiphoton processes and detection of collective resonances in dilute samples
,”
Phys. Rev. A
92
,
053412
(
2015
).
40.
L.
Bruder
,
M.
Mudrich
, and
F.
Stienkemeier
, “
Phase-modulated electronic wave packet interferometry reveals high resolution spectra of free Rb atoms and Rb*He molecules
,”
Phys. Chem. Chem. Phys.
17
,
23877
23885
(
2015
).
41.
K. J.
Karki
,
J. R.
Widom
,
J.
Seibt
,
I.
Moody
,
M. C.
Lonergan
,
T.
Pullerits
, and
A. H.
Marcus
, “
Coherent two-dimensional photocurrent spectroscopy in a Pbs quantum dot photocell
,”
Nat. Commun.
5
,
5869
(
2014
).
42.
M.
Liebel
,
C.
Toninelli
, and
N. F.
van Hulst
, “
Room-temperature ultrafast nonlinear spectroscopy of a single molecule
,”
Nat. Photonics
12
(
1
),
45
49
(
2018
).
43.
P.
Malý
and
T.
Mančal
, “
Signatures of exciton delocalization and exciton–exciton annihilation in fluorescence-detected two-dimensional coherent spectroscopy
,”
J. Phys. Chem. Lett.
9
(
19
),
5654
5659
(
2018
).
44.
P.
Grégoire
,
A. R.
Srimath Kandada
,
E.
Vella
,
C.
Tao
,
R.
Leonelli
, and
C.
Silva
, “
Incoherent population mixing contributions to phase-modulation two-dimensional coherent excitation spectra
,”
J. Chem. Phys.
147
(
11
),
114201
(
2017
).
45.
A. A. S.
Kalaee
,
F.
Damtie
, and
K. J.
Karki
, “
Differentiation of true nonlinear and incoherent mixing of linear signals in action-detected 2D spectroscopy
,”
J. Phys. Chem. A
123
(
19
),
4119
4124
(
2019
).
46.
A.
Perdomo-Ortiz
,
J. R.
Widom
,
G. A.
Lott
,
A.
Aspuru-Guzik
, and
A. H.
Marcus
, “
Conformation and electronic population transfer in membrane-supported self-assembled porphyrin dimers by 2D fluorescence spectroscopy
,”
J. Phys. Chem. B
116
(
35
),
10757
10770
(
2012
).
47.
F. A.
Damtie
,
A.
Wacker
,
T.
Pullerits
, and
K. J.
Karki
, “
Two-dimensional action spectroscopy of excitonic systems: Explicit simulation using a phase-modulation technique
,”
Phys. Rev. A
96
,
053830
(
2017
).
48.
O.
Kühn
,
T.
Mancal
, and
T.
Pullerits
, “
Interpreting fluorescence detected two-dimensional electronic spectroscopy
,”
J. Phys. Chem. Lett.
11
,
838
(
2020
).
49.
C. A.
Marx
,
U.
Harbola
, and
S.
Mukamel
, “
Nonlinear optical spectroscopy of single, few, and many molecules: Nonequilibrium Green’s function QED approach
,”
Phys. Rev. A
77
,
022110
(
2008
).
50.
P. L.
Kramer
,
C. H.
Giammanco
,
A.
Tamimi
,
D. J.
Hoffman
,
K. P.
Sokolowsky
, and
M. D.
Fayer
, “
Quasi-rotating frame: Accurate line shape determination with increased efficiency in noncollinear 2D optical spectroscopy
,”
J. Opt. Soc. Am. B
33
(
6
),
1143
1156
(
2016
).
51.
S.
Roeding
,
N.
Klimovich
, and
T.
Brixner
, “
Optimizing sparse sampling for 2D electronic spectroscopy
,”
J. Chem. Phys.
146
(
8
),
084201
(
2017
).
52.
Z.
Wang
,
S.
Lei
,
K.
Jung Karki
,
A.
Jakobsson
, and
T.
Pullerits
, “
Compressed sensing for reconstructing coherent multidimensional spectra
,”
J. Phys. Chem. A
124
,
1861
(
2020
).
53.
T.
Zhang
,
I.
Kuznetsova
,
T.
Meier
,
X.
Li
,
R. P.
Mirin
,
P.
Thomas
, and
S. T.
Cundiff
, “
Polarization-dependent optical 2D Fourier transform spectroscopy of semiconductors
,”
Proc. Natl. Acad. Sci. U. S. A.
104
(
36
),
14227
14232
(
2007
).
54.
K. W.
Stone
,
K.
Gundogdu
,
D. B.
Turner
,
X.
Li
,
S. T.
Cundiff
, and
K. A.
Nelson
, “
Two-quantum 2D FT electronic spectroscopy of biexcitons in GaAs quantum wells
,”
Science
324
(
5931
),
1169
1173
(
2009
).
55.
S.
Mueller
,
S.
Draeger
,
X.
Ma
,
M.
Hensen
,
T.
Kenneweg
,
W.
Pfeiffer
, and
T.
Brixner
, “
Fluorescence-detected two-quantum and one-quantum-two-quantum 2D electronic spectroscopy
,”
J. Phys. Chem. Lett.
9
(
8
),
1964
1969
(
2018
).
56.
E.
Thyrhaug
,
R.
Tempelaar
,
M. J. P.
Alcocer
,
K.
Žídek
,
D.
Bína
,
J.
Knoester
,
T. L. C.
Jansen
, and
D.
Zigmantas
, “
Identification and characterization of diverse coherences in the Fenna–Matthews–Olson complex
,”
Nat. Chem.
10
(
7
),
780
786
(
2018
).
57.
T.
Kramer
and
M.
Rodríguez
, “
Effect of disorder and polarization sequences on two-dimensional spectra of light-harvesting complexes
,”
Photosynth. Res.
144
(
2
),
147
154
(
2020
).
58.
D.
Paleček
,
P.
Edlund
,
E.
Gustavsson
,
S.
Westenhoff
, and
D.
Zigmantas
, “
Potential pitfalls of the early-time dynamics in two-dimensional electronic spectroscopy
,”
J. Chem. Phys.
151
(
2
),
024201
(
2019
).
59.
L.
Bruder
,
A.
Eisfeld
,
U.
Bangert
,
M.
Binz
,
M.
Jakob
,
D.
Uhl
,
M.
Schulz-Weiling
,
E. R.
Grant
, and
F.
Stienkemeier
, “
Delocalized excitons and interaction effects in extremely dilute thermal ensembles
,”
Phys. Chem. Chem. Phys.
21
,
2276
2282
(
2019
).
60.
J.
Strümpfer
,
M.
Şener
, and
K.
Schulten
, “
How quantum coherence assists photosynthetic light-harvesting
,”
J. Phys. Chem. Lett.
3
(
4
),
536
542
(
2012
).
61.
A. S.
Márquez
,
L.
Chen
,
K.
Sun
, and
Y.
Zhao
, “
Probing ultrafast excitation energy transfer of the chlorosome with exciton–phonon variational dynamics
,”
Phys. Chem. Chem. Phys.
18
,
20298
20311
(
2016
).
62.
G.
Lindblad
, “
On the generators of quantum dynamical semigroups
,”
Commun. Math. Phys.
48
(
2
),
119
130
(
1976
).
63.
P.
Tian
,
D.
Keusters
,
Y.
Suzaki
, and
W. S.
Warren
, “
Femtosecond phase-coherent two-dimensional spectroscopy
,”
Science
300
(
5625
),
1553
1555
(
2003
).
64.
H.-S.
Tan
, “
Theory and phase-cycling scheme selection principles of collinear phase coherent multi-dimensional optical spectroscopy
,”
J. Chem. Phys.
129
(
12
),
124501
(
2008
).
65.
A.
Anda
,
D.
Abramavičius
, and
T.
Hansen
, “
Two-dimensional electronic spectroscopy of anharmonic molecular potentials
,”
Phys. Chem. Chem. Phys.
20
,
1642
1652
(
2018
).
66.
G.
McDermott
,
S. M.
Prince
,
A. A.
Freer
,
A. M.
Hawthornthwaite-Lawless
,
M. Z.
Papiz
,
R. J.
Cogdell
, and
N. W.
Isaacs
, “
Crystal structure of an integral membrane light-harvesting complex from photosynthetic bacteria
,”
Nature
374
(
6522
),
517
521
(
1995
).
67.
A.
Anda
,
T.
Hansen
, and
L.
De Vico
, “
Multireference excitation energies for bacteriochlorophylls a within light harvesting system 2
,”
J. Chem. Theory Comput.
12
(
3
),
1305
1313
(
2016
).
68.
A.
Anda
,
L.
De Vico
, and
T.
Hansen
, “
Intermolecular modes between LH2 bacteriochlorophylls and protein residues: The effect on the excitation energies
,”
J. Phys. Chem. B
121
(
22
),
5499
5508
(
2017
).
69.
L.
De Vico
,
A.
Anda
,
V. A.
Osipov
,
A. Ø.
Madsen
, and
T.
Hansen
, “
Macrocycle ring deformation as the secondary design principle for light-harvesting complexes
,”
Proc. Natl. Acad. Sci. U. S. A.
115
(
39
),
E9051
E9057
(
2018
).
70.
A.
Anda
,
T.
Hansen
, and
L.
De Vico
, “
Qy and Qx absorption bands for bacteriochlorophyll a molecules from LH2 and LH3
,”
J. Phys. Chem. A
123
(
25
),
5283
5292
(
2019
).
71.
M. W.
Maciejewski
,
M.
Mobli
,
A. D.
Schuyler
,
A. S.
Stern
, and
J. C.
Hoch
, “
Data sampling in multidimensional NMR: Fundamentals and strategies
,”
Top. Curr. Chem.
316
,
49
77
(
2012
).