Time-resolved spectroscopy is commonly used to study diverse phenomena in chemistry, biology, and physics. Pump–probe experiments and coherent two-dimensional (2D) spectroscopy have resolved site-to-site energy transfer, visualized electronic couplings, and much more. In both techniques, the lowest-order signal, in a perturbative expansion of the polarization, is of third order in the electric field, which we call a one-quantum (1Q) signal because in 2D spectroscopy it oscillates in the coherence time with the excitation frequency. There is also a two-quantum (2Q) signal that oscillates in the coherence time at twice the fundamental frequency and is fifth order in the electric field. We demonstrate that the appearance of the 2Q signal guarantees that the 1Q signal is contaminated by non-negligible fifth-order interactions. We derive an analytical connection between an *n*Q signal and (2*n* + 1)th-order contaminations of an *r*Q (with *r* < *n*) signal by studying Feynman diagrams of all contributions. We demonstrate that by performing partial integrations along the excitation axis in 2D spectra, we can obtain clean *r*Q signals free of higher-order artifacts. We exemplify the technique using optical 2D spectroscopy on squaraine oligomers, showing clean extraction of the third-order signal. We further demonstrate the analytical connection with higher-order pump–probe spectroscopy and compare both techniques experimentally. Our approach demonstrates the full power of higher-order pump–probe and 2D spectroscopy to investigate multi-particle interactions in coupled systems.

## I. INTRODUCTION

A common technique to investigate ultrafast phenomena is pump–probe (PP) spectroscopy. Here, a first laser pulse excites a sample, and after some time delay *T* (“waiting time” or “population time”), another laser pulse probes the temporal evolution. The transient change in absorption at each population time *T* is obtained by measuring the transmitted spectrum of the probe pulse with and without a prior excitation pulse. These transient maps, available as a function of *T* and of detection wavelength, can be used to investigate ultrafast dynamics.^{1,2} Coherent two-dimensional (2D) spectroscopy can be viewed as an extension of PP spectroscopy that adds frequency resolution for the excitation step. 2D electronic spectroscopy was used to investigate the dynamics of various different systems, such as light-harvesting complexes,^{3,4} reaction centers,^{5} quantum dots,^{6,7} supramolecular aggregates,^{8,9} molecular dimers,^{10} carbon nanotubes,^{11} and 2D materials.^{12,13} Exemplarily, one can extract information on relaxation processes,^{14} homogeneous and inhomogeneous line shapes,^{15} energy transfer,^{16} chemical reaction kinetics,^{17} vibrational and electronic coherences,^{18} and exciton delocalization.^{19} It is worth pointing out that the extraction of state-to-state kinetic rate constants is uniquely determined for 2D spectroscopy under some basic conditions, such as the spectral separation of excitonic states,^{20} unlike for PP spectroscopy where the kinetic modeling is not unique. Note that we concentrate in the present work on 2D electronic spectroscopy, but the underlying concepts apply also to 2D infrared spectroscopy.^{21}

In 2D spectroscopy, a sequence of ultrashort laser pulses interacts with the system with varying time delays.^{21} The total number of pulses for a conventional third-order 2D experiment is often either three or four depending on the beam geometry.^{22,23} Either way, this results in a 2D correlation map for each population time *T* in which one axis is the excitation frequency and the other is the detection frequency. In such a 2D map, excitation and detection of the same transition result in a diagonal peak because the frequency in both steps is the same. Electronic couplings show up as cross peaks for *T* = 0, and energy transfer shows up as the time evolution of cross peaks for *T* > 0, obtained from scanning the population time *T* between the second pump pulse and the probe pulse.^{24} Analysis of cross peaks can be used to determine properties such as exciton delocalization and energy transfer.^{3,19,25}

In the current work, we generalize the analysis to cases in which signals in the 2D spectrum are separated by one or several multiples of the fundamental transition frequency and thus constitute “higher-quantum” peaks, and we also investigate how these signals herald the presence of higher-order contributions to peaks at the fundamental frequency or other, lower, multiple-quantum peaks. In conventional PP and 2D spectroscopies, one seeks to obtain the response in third order of perturbation theory. However, signals are often contaminated with higher-order contributions, and thus, we must extend the analysis to higher orders in the electric field of pump pulses. Higher-order contributions are well known in ultrafast experiments. In a typical experiment, higher-order signals are not isolated, and their influence becomes more relevant with stronger excitation intensity. In intensity-dependent experiments, the influence of higher-order signals was used to investigate 2D materials,^{26} nanorods,^{27} silicon nanocrystals,^{28} and chemical reactions.^{29} For excitonic systems,^{30–36} higher-order signals are sensitive to multi-exciton interactions, such as exciton–exciton annihilation.

Consider, for example, the investigation of photosynthesis with ultrafast spectroscopy.^{14,37–40} Natural light-harvesting systems utilize multiple chromophores to absorb light. These chromophores act as an “antenna” and a “funnel,” directing the absorbed energy to the reaction center of the photosynthetic complex.^{41} If in such a system single-exciton dynamics are investigated but the excitation power is too high, multiple excitons will be generated, leading to distorted kinetics.^{30,42–44} In terms of interaction with the light field, such a regime can be described by additional interactions with the excitation fields, i.e., more than the usual two “pump” interactions in third-order nonlinear experiments. For example, fifth-order contributions add the dynamics of exciton–exciton annihilation to the single-exciton dynamics. The onset of exciton–exciton annihilation depends on the system and on excitation conditions. Thus, in order to reduce artifacts, the accepted approach is to measure time-dependent signals for different excitation intensities. If the intensity change leads to a change in the kinetic evolution, one has to attenuate the excitation power until no such changes are observed.^{45} The problem is that reducing the excitation power also reduces the overall signal such that the signal-to-noise ratio decreases if the acquisition time is held constant. In extended systems with large absorption cross sections, such as molecular aggregates, polymers, or natural-light harvesting complexes, it may become impractical to achieve an adequately high signal-to-noise ratio in the annihilation-free regime, and then, one must accept compromises that may lead to erroneous interpretation of kinetic time constants. In addition, the criterion for a sufficiently small change of kinetics upon power reduction is often subjective.

In recent years, higher-order multidimensional nonlinear spectroscopy gained popularity.^{36,46–58} Higher-order spectroscopy can be used to measure higher-excited states of molecules,^{47,48} energy transfer in light-harvesting complexes,^{46} high-frequency vibronic modes,^{50} exciton–exciton annihilation,^{52,54,59} multiexciton states in quantum dots,^{56,57} coherences between multi-particle collective states,^{58} exciton diffusion in aggregates and polymers via annihilation,^{35,36,53} and non-equilibrium superconducting states in semiconductors.^{60} In coherently detected higher-order spectroscopy, a polarization higher than third order is measured, which requires more than three interactions with excitation fields. In action-detected spectroscopy,^{61–64} higher-order signals are connected to perturbation terms higher than fourth, rather than third, order because of the detection of a population in an excited state.^{54} In this publication, we focus on coherently detected spectroscopy.

In the pump–probe geometry, where three pulses interact with the system (i.e., two pump pulses and one probe pulse), higher-order signals can be isolated by their position along the excitation frequency axis in the 2D spectrum. At low excitation intensities, each pulse interacts once with the system, and the signal is emitted in the phase-matched direction ±**k**_{a} ∓**k**_{b} + **k**_{c} = **k**_{c} since **k**_{a} = **k**_{b} for the first two collinear pump pulses and **k**_{c} is the wavevector of the probe pulse. The “±” and “∓” signs reflect the fact that in the pump–probe geometry, the sum of rephasing and non-rephasing spectrum, i.e., the absorptive part of the spectrum, is directly obtained. The delay between the pump pulses, *τ*, is varied, and a Fourier transform taken with respect to *τ* to obtain the excitation axis of the correlation maps as a function of *ω*_{τ}. The first pump pulse excites a single-quantum coherence (1Q), and therefore, the signal appears at around the central frequency of the pump pulse along the excitation axis. The probe pulse excites a 1Q coherence as well, and therefore, the signal position along the detection axis is also fixed at around the central frequency of the probe pulse. We call this the “1Q1Q signal.” Note that the position of the signal is only approximately located at (in this case, single) multiples of the central laser frequency because the exact position is influenced by the absorption and emission spectra of the system.^{65}

For the increasing intensity of pump pulses, multiple interactions with the excitation field take place, leading to higher-order signals. One possible signal of next-highest order, with the phase-matching condition ±2**k**_{a} ∓ 2**k**_{b} + **k**_{c}, is emitted in the same direction, **k**_{c}, as the 1Q1Q signal (still for **k**_{a} = **k**_{b} in pump–probe geometry). However, due to one more additional interaction with the first pump pulse, the system evolves in a two-quantum (2Q) coherence after the first pump pulse, and the signal appears at around twice the central frequency of the pump spectrum. That is, if data are acquired with sufficiently small spacing of *τ*, the signal from this process will be separated in frequency from the 1Q1Q signal. The probe pulse still only interacts once with the system, and we thus call this the “2Q1Q signal.” We will focus, in this publication, on multiple interactions occurring with the pump pulses and consider the probe pulse to be always weak. Therefore, the second coherence in the signal stays always at the 1Q level, and we will drop this part of the label for the sake of simplicity. We utilized the 2Q signal previously to investigate exciton–exciton annihilation and exciton diffusion in supramolecular aggregates and polymers.^{35,36,53} Other methods to isolate 2Q signals use a noncollinear setup in which the three pulses interact with the sample from three different directions, in which case **k**_{a} ≠ **k**_{b}, and the rephasing (−**k**_{a} + **k**_{b} + **k**_{c}) and non-rephasing (+**k**_{a} − **k**_{b} + **k**_{c}) signals are measured separately.^{55} The 1Q and the 2Q signals are then separated by their distinct phase-matching directions. As we have shown recently, it is also possible to isolate 1Q and 2Q signals in PP spectroscopy.^{66–68} In that method, PP signals are measured at different excitation intensities from which different signals are extracted.

As mentioned above, the excitation intensity has to be chosen carefully in conventional PP or 2D spectroscopy experiments to reduce uncontrolled mixing of higher-order signals. In PP geometry, only odd orders of the nonlinear polarization are, in general, present due to the phase-matching condition. The 1Q signal contains, however, not only the usually desired third order but also all odd orders beginning with the third. The 2Q signal is a higher-order signal; it, likewise, not only contains fifth-order terms but also includes all odd orders beginning with the fifth. Thus, if the dynamics of the 1Q signal do not change with the increasing excitation intensity, the 1Q signal is dominated by a third-order response and the influence of higher-order contributions can be considered minor. The same principle holds for the 2Q signal: If the dynamics of the 2Q signal do not change as the excitation intensity increases, then the 2Q signal is dominated by the fifth-order response in terms of perturbation theory. However, in general, all *n*Q signals inherently contain higher-order contributions, starting at order 2*n* + 1 of perturbation theory, but their strength can be difficult to quantify experimentally. Therefore, even an isolation of different *n*Q signals in 2D and PP experiments, where *n* marks at how many multiples of the fundamental laser frequency the pump-induced coherence oscillates, does not allow for a clean separation of different orders of nonlinearity.

Inspired by isolating clean nonlinear signals in PP spectroscopy,^{66–68} we introduce a method for isolating signals of clean nonlinear orders in 2D spectroscopy in the PP geometry. Our approach works for any *n*Q signals at around *n* times the central frequency of the pump pulse along the excitation axis. As an example, we measure the 1Q, 2Q, and 3Q signals of squaraine oligomers and show how these signals can be used to eliminate fifth- and seventh-order contributions at the 1Q signal and, therefore, obtain a clean third-order signal even at high excitation intensities. This procedure works for *n*Q signals that are integrated over their respective excitation frequency.

This paper is structured as follows. In Sec. II, we discuss the theoretical concepts of isolating higher-order signals focusing on double-sided Feynman diagrams. Starting with general relations between *n*Q signals and nonlinear orders in perturbation theory (Sec. II A), we demonstrate how the *n*Q signals can be isolated at *τ* = 0 in 2D spectroscopy (Sec. II B) and PP spectroscopy (Sec. II C) and then demonstrate how individual nonlinear orders may be extracted from the *n*Q signals (Sec. II D). Section III contains the experiment. After we discuss the experimental setup and the sample in Sec. III A, we show exemplary 1Q, 2Q, and 3Q signals obtained on squaraine oligomers from 2D spectroscopy (Sec. III B) and PP spectroscopy (Sec. III C). In Sec. III D, we demonstrate how the clean third-order signal can be retrieved in both techniques by weighting and adding the different *n*Q signals. In Sec. IV, we summarize the main results and provide an outlook to future experiments.

## II. THEORY OF HIGHER-ORDER SPECTROSCOPY

### A. Relation between multi-quantum signals and perturbative orders

Generally speaking, *n*Q signals, which we will denote as *S*_{nQ}, are those signals that oscillate in coherence time *τ* at roughly *n* times the central frequency, *ω*_{0}, of the pump laser spectrum. In this paper, we take advantage of this separation of *n*Q signals along *ω*_{τ} (see Fig. 1 for an illustration of this effect). *n*Q signals are also often defined by the phase-matching condition that can be used to separate them, when three laser pulses with different wavevectors **k**_{a}, **k**_{b}, **k**_{c} are incident upon a sample.^{66–68} In this case, the rephasing *n*Q signals appear along the −*n***k**_{a} + *n***k**_{b} + **k**_{c} direction, and the non-rephasing *n*Q signals appear along the +*n***k**_{a} − *n***k**_{b} + **k**_{c} direction. In this paper, we focus on the case of PP geometry, where the first two pulses are collinear (**k**_{a} = **k**_{b}). Therefore, all of the rephasing and non-rephasing *n*Q signals propagate with wavevector **k**_{c} and are all measured simultaneously as a change in the absorption of the probe pulse (pulse c). In the PP geometry, *n*Q signals can also be obtained by phase cycling (varying the relative phase between the first two pulses). In this paper, we will use both phase cycling and separation by *ω*_{τ} to study the *n*Q signals.

As shown schematically in Fig. 1, the lowest-order contribution to *n*Q signals derives from a (2*n* + 1)th-order polarization in the sample. However, *n*Q signals always have contributions from higher-order terms as well. As with any nonlinear optical signal, we can expand *S*_{nQ} in orders of electric field amplitudes such that *S*_{nQ} is expressed as a sum over the perturbative terms $SnQ(2r+1)$, where the sum runs over *r* = *n*, …, ∞.

*n*Q signals precise, we write the electric fields as

**e**

_{j}is the polarization vector and $\epsilon jt$ is the complex amplitude of pulse

*j*,

*A*

_{j}is the complex envelope function (which includes spectral phase information of second order or higher in a Taylor expansion),

*ω*

_{j}is the central frequency,

**k**

_{j}is the wavevector,

*ϕ*

_{j}is the absolute phase,

*t*

_{j}is the pulse arrival time, and

**r**is the position. In a 2D experiment, three different time delays can be distinguished. The coherence time

*τ*describes the delay between the first two pulses. The population time

*T*is the delay between the center of the second pump pulse and the probe pulse. The signal time

*t*is the delay between the probe and the emission of the signal. The signal in time

*t*is usually measured implicitly by using a spectrometer, allowing for a direct measurement with respect to the Fourier conjugate variable

*ω*

_{t}. The phase of the emitted signals can be inferred via linear superposition with a known reference pulse, the “local oscillator.” In PP geometry, the probe pulse itself is the local oscillator, and thus, the intensity of the transmitted probe pulse is measured. By constructing the change in optical density, OD [see Eq. (16) in Sec. III A], we extract the transient nonlinear signal, which we refer to throughout this section as

*S*.

_{j}is unitless and $Aj,0t$ is the underlying shape. We focus on the case where the first two pulses are identical except for their arrival times,

*t*

_{j}, and their absolute phases,

*ϕ*

_{j}, so $Aa,0t=Ab,0t$ and λ

_{a}= λ

_{b}= λ. We assume that pulse c is weak with unchanged amplitude and expand

*S*

_{nQ}perturbatively in powers of λ,

*S*

_{nQ}are the

*n*Q signals and $SnQ2m+1$ are the perturbative expansion terms of those

*n*Q signals. We define $IP=\lambda 2$, where the subscript

*P*stands for pump and the letter

*I*indicates that this quantity is related to a unit-less scaling factor for the pump intensity. Therefore,

*n*Q signals can be expanded in powers of the intensity as

_{c}as we do not vary the intensity of the weak probe. Equation (1) includes the arguments

*τ*,

*T*,

*ω*

_{t}to show the explicit dependence of these signals with respect to experimental parameters. The left-hand side of Eq. (1) also includes an argument, showing that

*S*

_{nQ}is a function of the pump power, whereas the perturbative expansion terms $SnQ2m+1$ do not depend on the pump power (although they do depend on the shape of the pulses).

*n*Q terms are related to one another at

*τ*= 0 as

*r*≤

*n*. Equation (2) is the same as Eq. (S6) of the supplementary material. This relationship represents the key insight of this paper. It allows for the separation of the leading order $SnQ2n+1$ at

*τ*= 0 once the

*n*Q signals are measured. We emphasize that our correction procedure is only valid for

*τ*= 0, i.e., integrated 2D signals. We will explain in detail how to use this result to obtain individual perturbative signals in Sec. II D. We now give some physical intuition for understanding Eq. (2) by comparing Feynman diagrams that contribute to $S1Q5(\tau =0,T,\omega t)$ and $S2Q5(\tau =0,T,\omega t)$. Some example diagrams for both signals are displayed for a visual comparison in Fig. 2.

Double-sided Feynman diagrams are a common tool for visualizing and calculating various excitation pathways of the nonlinear polarization that contribute to a spectroscopic signal.^{69} In these diagrams, time runs from bottom to top. The density matrix is depicted in between two vertical lines. Interactions with electric laser fields are shown as arrows that represent either an excitation or a de-excitation of the system depending on whether an arrow points toward or away from the density matrix, respectively. The nonlinear polarization is emitted in a phase-matching direction that is dictated by the incident wavevectors. It is also possible to isolate different signals based on relative phases between pulses. An arrow pointing to the left contributes with a wavevector −**k**_{j} and a phase of −*ϕ*_{j} and an arrow pointing to the right contributes with +**k**_{j} and +*ϕ*_{j} to the nonlinear polarization. We present this discussion with the possibility of non-collinear pump pulses such that **k**_{a} ≠ **k**_{b}.

The impulsive limit is often assumed when performing calculations or developing intuition for the signals that are measured in a 2D experiment. The impulsive limit corresponds to the assumption that none of the pulses overlap in time. However, in realistic experiments, the first two pulses always overlap when the coherence time is smaller than the pulse duration, which is of critical importance in this article. Thus, we will not assume time ordering between the first two pulses. Furthermore, when the population time is smaller than the pulse duration, no time ordering exists for the interaction with the third pulse either (with respect to the first two pulses), and we must consider all possible temporal orderings of interactions. For population times, further away from zero, pulse c always arrives last and is well separated in time from pulses a and b, which is a situation we will call partial time ordering.

In partial time ordering, sixteen double-sided Feynman diagrams describe the third-order rephasing signal, and a closely related set of sixteen Feynman diagrams describe the non-rephasing signal. The numbers of rephasing and nonrephasing diagrams reduce to six if the rotating-wave approximation (RWA) is applied;^{69} see the supplementary material (Sec. SI). We focus our discussion in the main text on rephasing signals. All the arguments and derivations we provide apply equally to non-rephasing signals, which is discussed in the supplementary material (see especially Sec. SI). When full time ordering holds (that is, none of the pulses overlap), the rephasing signal is described by three diagrams that can be further classified as ground-state bleach (GSB), stimulated emission (SE), and excited-state absorption (ESA). If time ordering between pump pulses and probe pulse does not hold, a total of 16 diagrams contribute to the third-order signal.^{70} We show the full set of double-sided Feynman diagrams for different cases of time ordering in the supplementary material (Sec. SI).

Let us now discuss at which frequency positions and with which amplitudes higher-order signals contribute because we want to derive an exact relationship between the signals that appear at different excitation frequency positions. The fifth-order signal as the next higher-order signal appears at two different positions in a 2D spectrum: (1) at the 1Q excitation coordinate as a contamination of the third-order signal and (2) at the 2Q excitation coordinate (2*ω*_{0}) as a new signal, i.e., the 2Q signal.

The presence of a fifth-order 2Q signal, $S2Q5$, guarantees that the 1Q signal is contaminated by fifth-order signals such that $S1Q(5)$ is no longer negligible. The fifth-order contamination to the 1Q signal must obey the 1Q phase-matching condition, so such contamination occurs with two additional interactions with the same pulse, either a or b, but with oppositely signed wavevector contributions such that the emitted signal direction is not altered. That is, where the 1Q signal occurs at −**k**_{a} + **k**_{b} + **k**_{c}, the fifth-order contamination to the 1Q signal corresponds to interactions at −**k**_{a} + **k**_{a} − **k**_{a} + **k**_{b} + **k**_{c} or, equivalently, for triple interaction with **k**_{b}. In these cases, the contamination leads to a contribution to the 1Q signal from multi-particle dynamics and to a change of the linear scaling of the 1Q signal with higher excitation powers.^{71,72}

At even higher pump-pulse powers, the amplitude of the seventh-order signals will get stronger and exceed the noise threshold of the experiment, at which point they contribute significantly to the 2D spectrum at three different locations (as illustrated in Fig. 1 by the largest light-blue circles): (1) as a contamination of the 1Q signal, (2) as a contamination of the 2Q signal, and (3) as a signal that we call the 3Q signal, and similarly for yet higher orders. In Fig. 1, we show different positions of the 1Q, 2Q, and 3Q signals in the 2D spectrum with representative diagrams below. In Fig. 1, we consider diagrams up to seventh-order contributions. Higher-order contributions, such as the ninth-order signal, would add additional diagrams to Fig. 1 contributing to the 1Q, 2Q, and 3Q signals and diagrams to a 4Q signal. Note that Feynman diagrams do not include state labels as they customarily do because the conclusions we draw from our analysis of time-dependent perturbation theory are valid for any system and any type of states evolving between the pulses. Furthermore, we do not show the arrow representing the emission of the signal for any of the diagrams shown in this publication. Each column below the respective *n*Q signal in Fig. 1 contains the contributing nonlinear orders, each depicted by one exemplary out of many double-sided Feynman diagrams; Sec. SIV of the supplementary material gives the number of diagrams contributing at 1Q, 2Q, and 3Q, up to seventh order. The diagrams on the diagonal correspond to the lowest-order nonlinear signal that contributes to a specific *n*Q signal, and therefore, no diagrams above the diagonal are present. For example, the lowest-order signal at the 2Q position is of fifth order. Let us take a closer look at off-diagonal diagrams that represent the contamination by higher-order contributions. The contaminations must include pairs of interactions with a single pulse to maintain the phase-matching conditions associated with each type of signal. For example, the fifth-order signal contributing to 1Q has two more interactions with pulse a, as compared to third order, and those additional interactions occur with opposite wavevectors, i.e., arrows pointing in opposite directions in the double-sided Feynman diagram.

In the weak-probe limit, fifth-order signals $S1Q(5)$ and $S2Q(5)$ share many features. Both signals report on the same population dynamics and have the same line shapes along the detection frequency axis because the interaction with the probe pulse results in the same state of the system for both contributions. However, the line shapes differ along *ω*_{τ} since 1Q and 2Q evolve differently during *τ*, i.e., the colored red and blue rectangles in Fig. 1, respectively, in the fifth-order row occur after different numbers of prior-interaction arrows. Equation (2) tells us that at *τ* = 0, $S1Q(5)=4S2Q(5)$ if both pump pulses have identical shape. This equivalence can be proven by studying Feynman diagrams that contribute to these signals, and we motivate the proof visually in Fig. 2 using an exemplary set of diagrams.

For the example diagram of $S2Q5$ [Fig. 2(a)], each pump pulse interacts twice with the system, producing the overall phase-matching direction −2**k**_{a} + 2**k**_{b} + **k**_{c}. After the interaction with pulse a, the system evolves in a 2Q coherence, while during the population time, a population in a doubly excited state is reached. Other types of diagrams with a ground-state population or a single-exciton population after the interaction with pump pulses can also contribute to the signal.^{54,66–68} Pulse c only interacts once with the system. In the case where pulses a and b have identical envelopes, $Aat=Abt$, and polarizations, *e*_{a} = *e*_{b}, four fifth-order diagrams at the 1Q position can be constructed [Fig. 2(b)] that are identical to the diagram in Fig. 2(a) if we disregard the pulse labels for pulse a and pulse b. That is, while the 1Q diagrams contribute to the signal at *ω*_{τ} = *ω*_{0} and the 2Q diagram contributes to the signal at *ω*_{τ} = 2*ω*_{0}, when we, instead, consider *τ* = 0, all five of these diagrams produce identical contributions to the signal since there is no distinction between the a and b pulses.

We can construct all the fifth-order diagrams at the 1Q position by taking each diagram that contributes at the 2Q position and changing one interaction from a to b or vice versa. For each 2Q diagram, four different 1Q diagrams can be constructed, each with identical signal weight to the 2Q diagram at *τ* = 0. We demonstrate this relationship between 1Q and 2Q signals visually for every fifth-order diagram that obeys partial time ordering in Fig. S2 of the supplementary material. The key difference between the 2Q diagram and its four 1Q partners is that the 2Q diagram has two interactions with each pulse, while the 1Q diagrams have a single interaction with one pump pulse and three interactions with the other. Indeed, if we ignore labels “a” and “b” in the diagrams displayed in Fig. 2, we see that each of the diagrams looks identical, i.e., they include the same states and processes during the population time. The contribution of each diagram in Fig. 2 evolves differently as *τ* changes, which is why our argument only applies at *τ* = 0. Note that although we included the state labels in Fig. 2, we do not need to do so because the equivalence of the diagrams does not depend on the specific states that are involved.

The example shown in Fig. 2 contains the key idea of this article: For every diagram contributing to $S2Q5$, we can construct four equivalent diagrams that contribute to $S1Q5$ at *τ* = 0. The diagrams of the 1Q and 2Q signal are equivalent since they differ only by swapping an “a” label to “b” or vice versa, and the “a” and “b” pulses are identical at *τ* = 0. The diagrams contributing to the frequency-resolved 1Q and 2Q signal differ by which specific pulses the interactions occur. We demonstrate how the *τ* = 0 components of *n*Q signals may be extracted from a single 2D spectrum in Sec. II B. When the two pump pulses overlap, but time ordering between the pump and probe pulses holds, $S2Q(5)$ is described by 54 diagrams. If time ordering between a pump and a probe is not fulfilled either, the number of diagrams for $S2Q(5)$ increases to 240. In all cases, for each diagram that contributes to $S2Q(5)$, there are four diagrams that give an equivalent contribution to $S1Q(5)$ at *τ* = 0. The fifth-order contamination at the 1Q position, $S1Q(5)$, is precisely four times the fifth-order signal at the 2Q position, $S2Q(5)$, and can be described by 216 diagrams in the partial time ordering case and 960 diagrams in the case of no time ordering between pump and probe pulses. We show the full set of double-sided Feynman diagrams for the two contributions in the partially time-ordered case in the supplementary material (Sec. SII). All diagrams were generated using an automated diagram generator also delivering the count numbers quoted above.^{70}

We show in this article that we can use the measured *n*Q signals for *n* > *r* to eliminate higher-order contaminations in the *r*Q signal and obtain clean nonlinear signals. For example, the measurement of the 2Q signal allows us to eliminate the fifth-order contamination at the 1Q position, resulting in a clean third-order signal. Since fifth-order contributions at the 1Q and 2Q position have different line shapes along *ω*_{τ}, the relationship between the two signals as a function of *ω*_{τ} is not straightforward, and thus, we leave any discussion of an optional correction of the *ω*_{τ} line shapes for future work. We now turn our attention to describing two different methods of extracting the *τ* = 0 contributions to *n*Q signals. The first method described, in Sec. II B, involves integrating *n*Q signals over their respective frequency regions along the *ω*_{τ} axis. Note that such finite-interval integrations (over each *n*Q region separately) are not equivalent to an integration over all *ω*_{τ}, which would correspond to a measurement at *τ* = 0 according to the projection-slice theorem. The second method, described in Sec. II C, involves extracting different *n*Q signals in PP spectroscopy (i.e., for $\tau =0$ by exploiting their dependence on pump intensity.^{66,68} We again note that *n*Q signals can also be measured individually using phase matching. After describing two methods of extracting *n*Q signals, we demonstrate in Sec. II D how to use the relationship between them, shown in Eq. (2), to extract the individual nonlinear orders.

### B. Isolation of *n*Q signals in 2D spectroscopy

**k**

_{a}=

**k**

_{b}), and therefore, all

*n*Q signals have the same wave vector +

**k**

_{c}, and thus, all are detected as a change in the probe absorption. The total nonlinear signal induced by pump pulses may be written as

*I*

_{P}for this subsection because

*I*

_{P}remains fixed for this analysis. The multi-quantum signals $SnQ\tau ,T,\omega t$ oscillate at roughly

*nω*

_{0}as a function of

*τ*and can thus be separated by taking the Fourier transform with respect to

*τ*. However, we can see by inspecting Eq. (3) that at

*τ*= 0, all

*n*Q signals add together and are not experimentally separable. Therefore, we cannot take a single measurement at

*τ*= 0 and obtain the separated

*n*Q signals. We now demonstrate how we separate the

*τ*= 0 components.

*τ*,

*ω*

_{τ}axis near multiples of the pump pulse center frequency

*ω*

_{0}, as shown schematically in Fig. 1. If the pump pulse has a bandwidth Γ, we expect the signal $S\u03031Q\omega \tau ,T,\omega t$ to be mostly localized within the bounds $\omega 0\u2212\Gamma /2,\omega 0+\Gamma /2$. More generally, we expect the signals $S\u0303nQ\omega \tau ,T,\omega t$ to be mostly localized within the bounds $n\omega 0\u2212\Gamma /2,n\omega 0+\Gamma /2$.

We briefly note that the separation of *n*Q signals along the *ω*_{τ} axis is not guaranteed. Fortunately, if the signals did not separate, it would be visually clear: Instead of localized signals appearing around multiples of *ω*_{0} (as illustrated in Fig. 1), the signals would blend into each other. In such a case, the technique proposed in this paper would not be expected to work; see Sec. SXIII of the supplementary material for two examples where the separation of signals breaks down. We do, in fact, observe a breakdown in signal localization for population times close to 0, which results in errors in extracting *n*Q signals for those population times (see Sec. III D and Secs. SXII and SXIII for further discussion in the supplementary material). For the purposes of this discussion, we will assume that the signals do separate along *ω*_{τ}.

*n*Q signals separate along the

*ω*

_{τ}axis, we can recover, to a good approximation, the

*τ*= 0 component of each individual

*n*Q signal using

_{0}/2. To understand why Eq. (5) gives the

*n*Q signal at

*τ*= 0, it is helpful to recall that in the case where Δ → ∞, we obtain the total signal

*S*at

*τ*= 0 by the projection-slice theorem.

^{21,73}Since

*n*Q signals separate along the

*ω*

_{τ}axis, we can obtain a good approximation of $SnQ\tau =0,T,\omega t$ using Eq. (5). We can understand why the partial integration recovers

*τ*= 0 components of separate

*n*Q signals in a different way by substituting Eq. (4) into Eq. (5). This substitution leads to

*τ*= 0. In the limit of large Δ, we obtain

This limit corresponds to integrating the complete 2D spectrum along *ω*_{τ} and would recover the experimentally measured signal $S\tau =0,T,\omega t$. Setting an appropriate Δ thus leads to a selection of the *n*Q signal for each given *n*, while restricting the interferogram closely to *τ* ≈ 0. In the absence of noise, $\Delta =\omega 02$ is the optimal choice. However, in the presence of noise, smaller values of Δ may be desirable if *n*Q signals are well localized near multiples of *ω*_{0} along the *ω*_{τ} axis. In such a case, using the largest value $\Delta =\omega 02$ may simply add noise to the recovered signal.

We have thus far shown that we can, to good approximation, separate the contributions $SnQ\tau =0,T,\omega t$ via the use of a sinc filter, as shown in Eq. (6), even though those individual *n*Q signals cannot be experimentally directly measured at *τ* = 0 in the collinear PP geometry.

### C. Isolation of *n*Q signals in PP spectroscopy

We now discuss how to, instead, extract the signals $SnQ\tau =0,T,\omega t$ in PP spectroscopy, rather than 2D spectroscopy, using phase cycling. The signals $SnQ\tau ,T,\omega t$ may be extracted for any *τ* by controlling the relative phase between the first two pulses, *ϕ*_{ba} = *ϕ*_{b} − *ϕ*_{a}. One takes repeated measurements of the transient absorption of the probe pulse for different values of *ϕ*_{ba} and then takes the discrete Fourier transform with respect to the varying *ϕ*_{ba}. In this way, the discrete Fourier transform over the phase difference replaces the spatial signal separation that occurs in phase matching when the wavevectors are distinct (**k**_{a} ≠ **k**_{b}).^{74–76} The procedure for carrying out phase cycling is well-known,^{77} and we briefly describe the method here.

*τ*and

*T*, the change in the probe absorption is measured for $\varphi ba=\varphi p=\pi pN$, where

*p*= 0, 1, …, 2

*N*− 1. We denote the nonlinear transient signal obtained as $S\tau ,T,\omega t,\varphi p$. The discrete Fourier transform over

*ϕ*

_{p}is

*k*=

*N*.

^{78}For

*k*=

*n*<

*N*, we have $S\u0303n=SnQ$ and $S\u03032N\u2212n=SnQ*$. For the degenerate point

*k*=

*N*, Eq. (7) gives $S\u0303N=Re[SNQ]$. Phase cycling requires precise control of the relative phase between the two pump pulses, for example, using a pulse shaper. We solely focus on

*τ*= 0. We demonstrate elsewhere

^{66–68}that at

*τ*= 0, phase cycling is, in fact, intensity cycling of a single pump pulse since at

*τ*= 0, there is, in fact, only one pulse (with the two pump pulses precisely overlapping in space and time). To understand this relationship, consider the sum that the first two pulses have the same envelope and let their sum be $EPt=Ea(t)+Eb(t)$. We let $EP=\epsilon P+c.c.$, where

**k**

_{P}is the common wavevector of both pulses,

*ω*

_{0}is the carrier frequency,

**r**is the position (which is integrated out when taking into account the phase-matching condition from propagation through the sample), and c.c. denotes the complex conjugate of the previous term. The temporal intensity $Ipt\u2261\epsilon p2$ at

*τ*= 0 is then given by

^{77}but making use of Eq. (8), we can rewrite Eq. (7) to find the signals

*S*

_{nQ}at

*τ*= 0 as

*n*≤

*N*and

*p*is the intensity-cycling step index. Note that we have replaced $S\u0303k$ in Eq. (7) with

*S*

_{nQ}in Eq. (9) because we have restricted

*n*≤

*N*and because

*S*

_{nQ}are real-valued at

*τ*= 0 so that $S\u0303N=Re[SNQ]$ is no longer a special case. The argument

*I*

_{0}on the left-hand side of Eq. (9) indicates that all the

*n*Q signals are obtained at a single power, specifically the base power. This contrasts with the terms inside the sum on the right-hand side of Eq. (9), which are the PP signals measured at varying powers

*I*

_{p}. The signal

*S*$\tau =0,T,\omega t,Ip$ is the same as in Eq. (1).

*τ*≠ 0, phase cycling involves 2

*N*separate phase-cycling steps according to Eq. (7). However, when

*τ*= 0, phase cycling reduces to the intensity cycling of Eq. (8). In this case, we find that

*N*− 1 values of the 2

*N*values of intensity,

*I*

_{p}, are redundant because

*I*

_{p}for

*p*= 1, 2, …,

*N*− 2,

*N*− 1 and for

*p*= 2

*N*− 1, 2

*N*− 2, …,

*N*+ 2,

*N*+ 1, respectively, are the same, and therefore, the signals collected at those intensity pairs must be the same. We omit the value

*p*=

*N*, which corresponds to zero intensity, and note that

*p*= 0 has no redundant partner. Given that $e\u2212in\varphi p+ein\varphi p=2cos(n\varphi p)$, we may combine the redundant terms in Eq. (9) to obtain

*n*Q signal via linear combination. The Kronecker delta

*δ*

_{p,0}arises because the

*p*= 0 value had no redundant partner in Eq. (10).

*N*= 3. The 1Q signal can be isolated by setting

*n*= 1 in Eq. (11) and using Eq. (8) to find the values of

*I*

_{p}, which results in

*N*= 3 dataset, we can isolate the 2Q signal with

*n*= 2 by

*n*= 3 by

*n*Q signals for

*n*> 3, the number of intensity-cycling steps must be increased to

*N*=

*n*, resulting in smaller steps between the intensities.

### D. Extraction of nonlinear orders from *n*Q signals

We have described two methods for obtaining the signals $SnQ\tau =0,T,\omega t,I0$. Regardless of the method used for obtaining them, we may use Eq. (2) to isolate the individual nonlinear orders. We now demonstrate how to perform this extraction procedure.

*τ*= 0. As discussed in Sec. II A, for each fifth-order diagram at the 2Q position, four diagrams corresponding to a fifth-order contribution at the 1Q position can be constructed, which is the

*n*= 2,

*r*= 1 case of Eq. (2). The fifth-order contribution at the 1Q position can therefore be eliminated by subtracting four times the 2Q from the 1Q signal. If seventh-order corrections are negligible, the result is the clean third-order signal,

*τ*, which is always zero in this method. Note that in previous work, we already discussed contamination correction using higher-order signals in 2D spectroscopy.

^{35}However, in that work, we had derived an incorrect factor of six instead of the correct factor of four for the relation between the fifth-order signal at 2Q and the fifth-order signal at the 1Q position because we had considered only time-ordered diagrams, whereas in the present work, we take into account all diagrams including those without time ordering between the two pump pulses.

*n*= 3 and

*r*= 2, resulting in a ratio of six between the two contributions. If the ninth-order responses are negligible, the clean fifth-order contribution to the 2Q signal can be obtained by subtracting six times the 3Q signal from the 2Q signal. For obtaining a clean third-order signal, we have to eliminate both the fifth order (by subtracting four times the corrected 2Q data) and the seventh order [by subtracting 15 times the 3Q signal, obtained from Eq. (2) with

*n*= 3,

*r*= 1]. Since the 3Q signal is also used to correct the 2Q signal, it turns up twice in the overall correction procedure: to eliminate the seventh-order contribution from 1Q and to eliminate the seventh-order contribution from 2Q. The overall correction procedure thus results in

*S*

^{(2n+1)}$T,\omega t$ with the measured

*n*Q signals $SnQT,\omega t,I0$,

*n*Q signal has order 2

*n*+ 1, the matrix is upper triangular. Note that the correction factors are the same regardless of whether time ordering between the pump and probe pulses is fulfilled or the rotating-wave approximation (RWA) holds, as we show in the supplementary material (Sec. SIV). An important point is that the procedure corrects up to a certain nonlinear order. For example, in Eq. (14), the 1Q signal is corrected for fifth- and seventh-order contributions using the 2Q and 3Q signals. However, if the 3Q signal is contaminated by non-negligible ninth-order contributions, the correction using only the 2Q and 3Q signal is not sufficient, and the 4Q signal has to be taken into account. For any given pump intensity, there is a high-order signal that is small enough to be neglected, and orders lower than that must be corrected. In our recent publication, we used an independent approach using this intensity dependence to derive Eq. (15).

^{68}In summary, Eq. (15) allows us to obtain

*N*clean nonlinear signals, given that we have first obtained

*N*multi-quantum (

*n*Q) signals.

## III. RESULTS AND DISCUSSION

### A. Experiment

We now illustrate the general concept derived in Sec. II exemplarily by describing measurements and the associated analysis on a squaraine oligomer (oSQB8) dissolved in toluene. The absorption spectrum is shown in Fig. 3(a) (green) together with the laser spectra of the pump (red shaded area) and probe pulses (black dashed line). The oligomer is made from eight repeating units of the cisoid indolenine squaraine monomer (SQB) [Fig. 3(b)]. A thorough analysis of the fluorescence and absorption spectra on a series of squaraine oligomers with varying lengths made out of SQB revealed that spectroscopic properties in toluene can be explained by a structurally disordered linear chain.^{79} The linear structure in toluene results in a J-type coupling and a red-shifted absorption maximum compared to the monomer. In our measurements, we excited the blue part of the main absorption peak.

^{19,53}consists of a Ti:sapphire laser amplifier (SpitfirePro, Spectra Physics, 1 kHz, 800 nm, 4 mJ, 35 fs) whose output attenuated to about 1 mJ pulse energy was focused into a hollow-core fiber (Ultrafast Innovations) filled with neon (∼1 bar). The resulting broadband white light was compressed of a set of chirped mirrors and then split into a pump beam and a probe beam by a pair of wedges. The population time

*T*was controlled by a mechanical delay stage (M-IMS600LM, Newport) within the probe beam. The pump beam was guided for further compression through a GRISM compressor (Fastlite) and an acousto-optic programmable dispersive filter (AOPDF) pulse shaper (Dazzler, Fastlite). For PP measurements, the AOPDF was used to control the intensity of the pump pulse. For 2D experiments, the AOPDF was used to create a double pulse with variable delay

*τ*. The pump and probe beams were focused into the sample via two focusing mirrors, resulting in beam sizes of 59

*μ*m full width at half maximum (FWHM) for the main axis of the probe (at an eccentricity of 0.75) and 335

*μ*m FWHM for the pump (eccentricity of 0.83). The pulse duration of a single pump pulse was measured by pulse-shaper-assisted collinear frequency-resolved optical gating (cFROG)

^{80}and found to be 17 fs (FWHM of the intensity). The pump pulse characterized by cFROG was used to determine the pulse duration of the probe pulse in a second step by cross-correlation FROG.

^{80}The reconstructed duration of the probe pulse was ∼50 fs (FWHM), slightly longer than the pump pulse due to some remaining dispersion because the pulse shaper acted only on the pump beam. To obtain the signal (in TA and in 2D), we used a double-chopping scheme in which both beams were chopped with different frequencies. In the pump beam, sequences of two consecutive pulses were blocked and then transmitted; in the probe beam, alternating pulses were blocked and transmitted. This resulted in four combinations, all measured as powers of the probe beam: (1) both beams blocked corresponding to the background only (

*I*

_{B}), (2) both beams open (

*I*

_{Pu−Pr}), (3) pump beam open and probe beam blocked (

*I*

_{Pu}), and (4) vice versa (

*I*

_{Pr}). The transient absorption signal as the change of the optical density was constructed by

*μ*m beam path.

### B. 2D spectroscopy

We begin with the isolation of *n*Q signals, starting with 2D spectroscopy in this section. Carrying out the procedure derived in Sec. II B requires scanning of the coherence time in small enough steps to separate different *n*Q signals along the excitation axis in the 2D spectrum. We want to demonstrate our correction procedure at different excitation intensity regimes, i.e., where different higher-order signals are contributing. Therefore, we measured 2D spectroscopy at three excitation intensities and isolated the 1Q, 2Q, and 3Q signal for each excitation intensity. In practice, we kept all pulse parameters constant and varied the pulse energy, which is proportional to intensity, using the pulse shaper. We refer to these measurements by the pulse energy of the excitation pulses at temporal overlap (*τ* = 0). The probe pulse was always kept weak such that it effectively interacted only once with the system.

The 2D spectrum at a population time of *T* = 220 fs, measured with an excitation pulse energy of 220 nJ, is shown in Fig. 4(a). The signal near the excitation frequency of *ω*_{0} (corresponding to $\u223c13000$ cm^{−1}) represents the 1Q signal (blue rectangle). It is dominated by a negative signal, which is elongated along the detection axis. At around twice the excitation frequency (2*ω*_{0}), a positive 2Q peak is visible (red rectangle), called “EEI2D signal” in previous publications.^{35,36,54} Here, we adopt the “*n*Q” nomenclature to enable generalization to arbitrary orders *n* and non-excitonic systems.

In Fig. 4(b), we show the individual regions of 1Q, 2Q, and 3Q regions with rescaled color bars such that we can clearly distinguish the three different *n*Q signals. Interestingly, for higher-order signals (that we define as all *n*Q signals with *n* > 1), a pronounced elongation along the excitation axis can be observed. All *n*Q signal maxima are located slightly below their respective diagonals ($\nu \u0304t=\nu \u0304\tau $ for 1Q, $2\nu \u0304t=\nu \u0304\tau $ for 2Q, and $3\nu \u0304t=\nu \u0304\tau $ for 3Q). The reason for the peak position could be energy relaxation and stimulated emission contributions, leading to a shift to smaller detection wavenumbers. The line shapes of the 2Q signal were investigated theoretically for molecular dimers and trimers^{59} and directly compared to experimental results on the example of a squaraine trimer.^{55} We leave an analysis of *n*Q line shapes in the present oligomer to future work and instead concentrate on the signals integrated over $\nu \u0304\tau $ for each *n*Q region shown in Fig. 4(b) (1Q: 11 258–15 010 cm^{−1}, 2Q: 24 183–27 936 cm^{−1}, and 3Q: 37 109–40 862 cm^{−1}).

We show in Fig. 5(a) the integrated 2D signals of the respective regions for three different pulse energies. The $\nu \u0304\tau $ integration, while keeping the $\nu \u0304t$ coordinate, leads to effective transient absorption maps. We adjusted small differences of the temporal overlap between pump and probe pulses (*T* = 0, *τ* = 0) between the measurements by shifting the population time axis using the coherent artifact in the 1Q signal, as described in the supplementary material (Sec. SV). This shift only affected the population time axis and not the coherence time. We show the 1Q, 2Q, and 3Q signals (from left to right) for pulse energies of 15, 120, and 220 nJ (from top to bottom). The 1Q signal at all pulse energies is dominantly negative, which corresponds to ground-state bleach and stimulated emission in our sign convention. Now, we focus on higher-order signals. For the lowest pulse energy (15 nJ), only a 1Q signal is visible. At the next higher pulse energy (120 nJ), a 2Q signal can be clearly seen, and the 3Q signal just begins to emerge from the noise floor. At the highest pulse energy (220 nJ), all three signals (1Q, 2Q, and 3Q) are clearly visible.

In Fig. 5(b), we show the average (along $\nu \u0304t$ over the marked regions in the maps) transient signals for 1Q, 2Q, and 3Q. These transients are normalized to the average signal from 200 to 300 fs to facilitate a direct comparison of *n*Q kinetics. Thus, the transients are all positive for long delay times and no longer reflect the sign alternation from Fig. 5(a). The 1Q signals at 220 nJ [Fig. 5(b), blue solid line] and 15 nJ (blue dotted line) both rise at around *T* = 0. The differences between the two measurements are visible in the first 100 fs after *T* = 0. During the first 100 fs, the 1Q signal at 220 nJ rises to a maximum above 1.0 and then slightly decreases. The 1Q signal at 15 nJ first initially agrees with the signal at 220 nJ but rounds off without building up a maximum. The change in the dynamics between the measurements at different pulse energies is a clear indicator that at high pulse energies, higher-order contributions contaminate the 1Q signal.

With the *n*Q analysis described in the present work, contamination of 1Q is trivial to confirm because we can directly measure the higher-order signals, i.e., at 2Q [Fig. 5(b), red solid line] and 3Q positions [Fig. 5(b), yellow solid line] at 220 nJ. These higher-order signals correspond to negations of the lower-order signals and SE and ESA from higher-excited states, such as biexciton states. While the third-order signal reports on single-exciton dynamics, the fifth-order signal adds the dynamics of exciton–exciton annihilation, the seventh-order signal adds the annihilation of three excitons, and analogously for higher-order signals. Note that each higher-order signal also contains reductions of lower-order signals, just as the ground-state bleach at third order is a reduction of the linear absorption. For example, the fifth-order signal includes the dynamics of states that begin with two excitons and single- and zero-exciton dynamics.

As previously discussed, in extended systems, higher-order signals are dominated by the signal belonging to the maximum number of excitons that can be probed in the given order of nonlinearity.^{53,68} The 2Q signal at 220 nJ rises more slowly than the 1Q signal at 220 nJ during the first 100 fs although the difference is small. The 2Q signal rises slower since at short time, the SE from the two-exciton states opposes the contributions from ground state and single-exciton states.^{54} As the excitons annihilate, the contribution of the two-exciton states to the overall signal becomes weaker, allowing the total signal to grow in magnitude. In order to annihilate, the excitons have to arrive in close proximity to each other, and therefore, a period of exciton diffusion takes place, which is reflected by a rise of the 2Q signal. Compared to larger excitonic systems that were investigated previously, such as polymers and aggregates,^{35,53} the rise of the 2Q signal is faster here since the studied oligomers are smaller, and therefore, the excitons do not need to propagate for large distances before annihilation. For a better comparison, we show the 1Q and 2Q signals at 220 nJ on top of each other in Fig. S6 of the supplementary material (Sec. SVI). At an excitation pulse energy of 120 nJ, the 2Q signal [Fig. 5(b), red dashed line] rises more slowly than for 220 nJ (red solid line). The slower 2Q rise for lower pulse energies is an analogous effect as that observed for 1Q, i.e., higher-order signals contaminate the signal. In this case, at 220 nJ, the seventh-order signal appears at the 2Q position. We can observe the seventh-order signal directly at the 3Q position [Fig. 5(b), yellow solid line]. A comparison of the measured dynamics of the 3Q signal with the 1Q and 2Q signals is difficult because of the low signal-to-noise ratio. While the 2Q signal is dominated by the dynamics of biexcitons, the 3Q signal additionally contains the dynamics of triexcitons, i.e., the annihilation of three excitons.

### C. PP spectroscopy

As discussed in Sec. II C, *n*Q signals can also be extracted from PP experiments using a systematic variation of the excitation intensity. We show in Fig. 6(a) the 1Q (left), 2Q (middle), and 3Q signals (right) measured via PP. In order to obtain the 1Q, 2Q, and 3Q signals, a set of three different pulse energies is needed as discussed in Sec. II C. We constructed three different sets of 1Q, 2Q, and 3Q signals using Eq. (11), which are labeled by the highest power in each sequence: 200 nJ (constructed from measurements at 200, 150, and 50 nJ), 50 nJ (constructed from measurements at 50, 37.5, and 12.5 nJ), and 12.5 nJ (constructed from measurements at 12.5, 9.375, and 3.125 nJ).

We show the extracted 1Q, 2Q, and 3Q signals for the three different values of 4*I*_{0} in the supplementary material (Sec. SVII). Let us now focus on the transients arising after averaging over the marked regions in Fig. 6(a) and normalization to the average signal between 200 and 300 fs [Fig. 6(b)]. The difference in the dynamics for the measurements at different excitation powers is again mostly visible in the population times from *T* = 0 to *T* = 100 fs. For the 1Q signal at 200 nJ (4*I*_{0}), the signal rises rapidly at around time zero, reaches a maximum above 1.0, and then decays [Fig. 6(b), blue solid line]. The 1Q signal of 4*I*_{0} = 12.5 nJ (blue dotted line) does not reach a value above 1.0 and has a more rounded shape. The differences of the dynamics between these two curves are again a clear sign that at 200 nJ, the 1Q signal is contaminated by higher-order contributions. The change in the dynamics is also visible for the 2Q signal at 200 nJ (red solid line) compared to 50 nJ (red dashed line), which can be explained by the seventh-order contribution, clearly visible as a 3Q signal (yellow solid line). The 2Q signal at 200 nJ [Fig. 6(b), red solid line] rises more slowly than 1Q, which is better visible by a direct comparison of the traces, as shown in the supplementary material (Sec. SVIII).

The difference for the measurements between low and high powers in PP spectroscopy is similar to the differences between low and high excitation powers in 2D spectroscopy. The 1Q signal in PP spectroscopy at low power (4*I*_{0} = 12.5 nJ) has a similar roundish shape as the 1Q signal at lower powers in 2D measurements [Fig. 5(b), blue dashed line]. In addition, the 2Q signal for PP and 2D spectroscopies (dark gray solid line) has similar dynamics for *T* > 0. For *T* near and less than zero, *n*Q signals from PP and 2D methods do not agree, with the most obvious differences in the 3Q signal. The problem arises in 2D signals, where the *n*Q signals cease to be well-separated from each other along the *ω*_{τ} axis. The windowed integration along *ω*_{τ} then does not capture the full *n*Q signal and therefore does not well reproduce the *τ* = 0 contribution to the signals. We illustrate this problem in the supplementary material Sec. SXIII, showing that it occurs when *T* becomes smaller than the pulse duration. The PP signals do not rely on the windowed *ω*_{τ} integration, so those signals are more accurate in the small- and negative-*T* regions.

In summary, the higher-order signals in 2D and PP spectroscopies reveal similar dynamics except for the population times close to *T* = 0. The contamination of higher-order signals in both methods is visible as a change in the dynamics of the 1Q signal by 2Q and 3Q contributions and of the 2Q signal by 3Q contributions, i.e., faster rises in each signal are seen for higher excitation pulse energy.

### D. Extraction of contamination-free third-order signal

We now focus on the correction of the 1Q signal with respect to higher-order contributions using both 2D and PP data. In Fig. 7, we show 1Q transients from 2D (a) and PP (b) measurements, averaged over $\nu \u0304t$ (from 12 200 to 12 800 cm^{−1}). In the case of 2D spectroscopy, the measurements of 220 and 15 nJ are already shown in Fig. 5(b). The transients extracted from a set of three PP measurements with 4*I*_{0} of 200 nJ are shown in Fig. 6(b). For the high-power measurement (220 nJ for 2D, 4*I*_{0} of 200 nJ for PP spectroscopy), the 1Q signal is contaminated by higher-order contributions, i.e., third-, fifth-, and seventh-order signals are mixed together. The spectrally resolved maps are shown in the supplementary material (Sec. SIX).

As a reference for 2D measurements, we use the signal measured at *I*_{ref,2D} = 15 nJ. For PP spectroscopy, we use the signal measured at *I*_{ref,PP} = 3.125 nJ, which is the lowest intensity-cycling step for the 12.5 nJ set. We call these two signals “reference signals.” The transients are scaled using two factors. First, each signal is divided by the absolute value of the transient (averaged from *T* = 200 to *T* = 300 fs) of the reference signal. In addition, as shown in Eq. (11), the extraction procedure of the *n*Q signals from the PP data yields the quantities $SnQT,\omega t,I0$, which is to say that the signals are extracted at a particular synthetic power *I*_{0}, even if the intensity-cycling procedure does not happen to include *I*_{0}. This means that to leading order, the *S*_{nQ} signals scale as $I0n$. Similarly, according to Eq. (15), the individual orders that are obtained are scaled as $S2n+1T,\omega tI0n$. Therefore, all 1Q and third-order signals are also normalized by the ratio of *I*_{0}/*I*_{ref}. Therefore, the 1Q signal extracted from the set with 4*I*_{0} corresponding to 200 nJ, 3*I*_{0} corresponding to 150 nJ, and *I*_{0} corresponding to 50 nJ is divided by the ratio (50/3.125). Analogously, the 1Q signal extracted from measurement of 50 nJ, 37.5, and 12.5 nJ is divided by (12.5/3.125). As seen in Fig. S10 of the supplementary material, the raw PP data at 3.125 nJ are practically indistinguishable from the corrected 1Q signal constructed from the full set of three intensity cycling steps, giving us confidence that the original PP data taken at 3.125 nJ can serve as a reasonable reference for a pure third-order signal. If the 1Q signal dominated by a negative signal scaled linearly with excitation intensity, i.e., if no higher-order contaminations contributed at higher pulse energies, the minimum signal of each trace would be at −1.0 with this normalization procedure because we divide the negative signal by the absolute of the averaged reference signal and the ratio of powers, which are both positive. Saturation of the signal for higher excitation intensity leads to lower signal amplitude and, therefore, a change of the signal from −1.0 to, for example, −0.8.

In Fig. 7, uncorrected transients are shown as solid lines. For the lowest-power measurements in the 2D experiment, no higher-order signals could be detected. Therefore, this signal should be free from contamination from higher orders [Fig. 7(a), light blue line]. For higher pulse energies, such as 120 nJ [Fig. 7(a), dark blue solid line] and 220 nJ [Fig. 7(a), purple solid line], the signal saturates, which leads to a signal of about −0.6 for the measurement at 120 nJ and −0.4 for the measurement at 220 nJ. As stated above, the change of the signal dynamics by higher-order contributions is small in this particular sample, but changes can be dramatic in the general case of other excitonic systems.^{30,32,71} We therefore use the signal strength as an additional indicator for whether a clean third-order signal can be recovered.

We now apply the correction procedure, Eq. (13), from Sec. II D to the measurements of 120 nJ by subtracting four times the 2Q signal from the 1Q signal, resulting in a clean third-order signal, which is shown in Fig. 7(a) (dark blue dashed line). For the measurement at 220 nJ, we use both the 2Q and 3Q signals for correction according to Eq. (14) and show the result in Fig. 7(a) (purple dotted line). The corrected curves show similar dynamics as the low-power measurement, and the scaling is also correct, i.e., their minimum signals are around −1.0. At close inspection, it can be seen that the corrected curves are both slightly below the reference measurements. This suggests that there is a small contribution of fifth-order contamination present even for the lowest-power measurement. There are larger discrepancies at negative *T*, as discussed in Sec. III C.

For the PP experiments shown in Fig. 7(b), the situation is similar to 2D experiments. Since the pulse energies are slightly lower in PP measurements than in 2D measurements, the saturation for the power of 50 nJ [Fig. 7(b), dark blue solid line] is not as strong as in the 2D measurement. The change in the dynamics is also not as clear for the power of 50 nJ [Fig. 7(b), dark blue solid line] as for 2D spectroscopy. However, for the highest power of 200 nJ [Fig. 7(b), purple solid line], the signal starts to lose its roundish shape. We compare the signal to the annihilation-free PP measurement at 3.125 nJ [Fig. 7(b), orange line]. The corrections of the 50 nJ data up to 2Q (dark blue dashed) and the 200 nJ data up to 3Q (purple dotted) eliminate higher-order contaminations. We show that no seventh-order contamination is present in the 50 nJ measurement in the supplementary material (Sec. SXI). Both of the corrected curves and the low-power measurement in Fig. 7(b) agree well, and only at close inspection, a slight systematic deviation can be noticed. For the highest power, the corrected curve [Fig. 6(b), purple dotted line] is slightly above the lowest-power measurement. This is a hint that at 200 nJ, the 3Q signal is contaminated by a (small) ninth-order contribution that we did not extract with our measurement scheme but could be obtained using one more measurement step in the intensity-cycling procedure. In the 2D experiment, the signal at highest pulse energy seems to be corrected reasonably well using the 2Q and 3Q signals and does not require any 4Q corrections. Slightly different beam sizes or a variation in the spatial overlap between pump and probe beams between the 2D and PP measurements might also be reasons for the small deviations between the two experiments.

In addition to the transients of Fig. 7, we show in Fig. 8 cuts along $\nu \u0304t$ for the 1Q signal with and without correction for higher-order contamination at two population times, 40 and 300 fs. The transient spectra for *T* = 0 fs are shown in Sec. SXII of the supplementary material. For a population time of 300 fs, the 2D *n*Q signals are well localized in *ω*_{τ}, so the windowed integration proceeds without issue. However, at a population time of 40 fs, the 2D *n*Q signals are not well localized in *ω*_{τ}, and thus, the windowed integration procedure does not fully capture the *τ* = 0 component of the *n*Q signals, which causes the discrepancies between the 15 nJ measurement and the corrected measurements, seen in the top panel of Fig. 8(a), as well as the visual differences between 2D and PP techniques. We directly compare the data from 2D spectroscopy [Fig. 8(a)] and from PP spectroscopy [Fig. 8(b)]. The signals are scaled in a similar way as in Fig. 7: First, the signal for each specific *T* step is divided by the absolute value of the minimum of the reference signal, and then, each signal is divided by the ratio between the pulse energy of the measurement and the reference pulse energy.

Both methods shown in Fig. 8 show signs of saturation of the 1Q signal as the pump power increases, as previously discussed. Figure 8(a) also shows further evidence that the reference 2D measurement taken at 15 nJ includes some fifth-order contamination since the corrected curves are slightly below −1.0. For the larger population time of *T* = 300 fs, the spectra from both methods are similar. For both population times for PP (and for *T* > 40 fs for 2D), the corrections rebuild the shape of the reference signal. We emphasize that our procedure allows one to judge very easily the excitation regime, i.e., which higher-order contaminations are present. Besides the correction, one also obtains higher-order contributions that contain additional information about the interaction of multiple (quasi)particles and higher-excited states.

## IV. CONCLUSION

In this paper, we demonstrated a general procedure that can be applied in coherent two-dimensional (2D) spectroscopy and in pump–probe (PP) spectroscopy to obtain clean nonlinear signals of a particular perturbative order. The first step in both methods is to isolate multi-quantum (*n*Q) signals. In 2D spectroscopy, different signals can be isolated by their position along the excitation-frequency axis. In PP spectroscopy, *n*Q signals are isolated by an intensity-cycling procedure where linear combinations of PP signals measured at specified excitation intensities produce various *n*Q signals. General procedures for these linear combinations were derived that are also applicable to higher orders. The second step involves taking appropriate linear combinations of *n*Q signals in order to isolate the individual nonlinear orders. We derived the analytical connection between the *n*Q signals by studying the structure of the Feynman diagrams that contribute to each nonlinear order of each *n*Q signal and provided visual motivation for the relationship by inspecting the fifth-order diagrams that contribute to the 1Q and 2Q signals. Based on the knowledge of how much each nonlinear order contributes to different *n*Q signals, we developed a correction protocol using the isolated *n*Q signals in PP and 2D spectroscopies. The clean nonlinear orders were extracted by adding *n*Q signals together with specific correction factors. This result is remarkable because the correction procedure is independent of the explicit system properties.

We experimentally confirmed our theoretical analysis by isolating the 1Q, 2Q, and 3Q signals of squaraine oligomers in 2D spectroscopy and demonstrated that these signals can be isolated in a PP experiment as well. We showed that at high excitation pulse energies, fifth- and seventh-order contributions contaminate the 1Q signal, leading to a change in observed dynamics. In our example, these contaminations arose from multi-exciton processes, such as exciton–exciton annihilation. Additionally, a saturation of the signal with higher pulse energies was observed. Following a full analysis of all contributing Feynman diagrams from all relevant higher orders, the measured 2Q and 3Q signals were utilized to correct the 1Q measurement for contaminations from higher orders, resulting in a clean third-order signal.

The experimental *n*Q results and their corrections were found to be equivalent for 2D and PP spectroscopies, apart from some differences at negative and small positive population times that are related to streaking of the 2D data along the excitation frequency axis. Using the correction procedure, an annihilation-free third-order signal was obtained even for high pulse energy. We confirmed our approach by comparing the corrected 1Q (i.e., third-order) signals with measurements at low pulse energies, i.e., in the absence of higher-order contamination, and found the dynamics to be the same.

Our technique poses a simple solution to the longstanding problem to avoid annihilation in ultrafast spectroscopic experiments. The extraction of annihilation-free signals is especially interesting for the investigation of natural light-harvesting complexes or other “chromophore-crowded” systems where annihilation is challenging to avoid. While in this publication we focused on one specific case of higher-order signals, namely, multi-exciton annihilation, the problem of uncontrolled mixing of higher-order signals is more general, and our approach can be applied to other systems where the higher-order effects report on phenomena other than annihilation. In any system, this procedure extracts clean third-order signals.

Furthermore, one gains simple access to higher-order signals, allowing us to investigate multi-particle interactions and probe their dynamics.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the double-sided Feynman diagrams of the third-order 1Q signal, double-sided Feynman diagrams of the fifth-order 1Q and 2Q signals, analytic proof of the correction procedure, ratio of double-sided Feynman diagrams, shifting of the population time, 1Q and 2Q signals obtained from two-dimensional spectroscopy, pump–probe signals for different pulse energies, comparison of pump–probe signals, corrected 1Q signals for different pulse energies, corrections of the pump–probe signal at a pulse energy of 12.5 nJ, corrections of the pump–probe signal at a pulse energy of 50 nJ, correction at T = 0 fs, and separation of overlap diagrams along the ωτ axis.

## ACKNOWLEDGMENTS

We acknowledge the funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant No. 423942615 (T.B.) and the Canadian Natural Sciences and Engineering Research Council (J.J.K. and P.A.R.). P.M. was supported by the Alexander von Humboldt Foundation. J.L. was supported by a scholarship of the Cusanuswerk.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author contributions

J.L. and P.A.R. contributed equally to this work.

**Julian Lüttig**: Conceptualization (equal); Formal analysis (lead); Investigation (equal); Methodology (supporting); Software (supporting); Validation (lead); Visualization (lead); Writing – original draft (equal). **Peter A. Rose**: Conceptualization (equal); Formal analysis (supporting); Investigation (equal); Methodology (lead); Software (lead); Validation (supporting); Visualization (supporting); Writing – original draft (equal). **Pavel Malý**: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting). **Arthur Turkin**: Resources (lead). **Michael Bühler**: Formal analysis (supporting); Investigation (supporting). **Christoph Lambert**: Funding acquisition (supporting); Supervision (supporting). **Jacob J. Krich**: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal). **Tobias Brixner**: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

## REFERENCES

*Concepts and Methods of 2D Infrared Spectroscopy*

*Molecular Mechanisms of Photosynthesis*

*Principles of Nonlinear Optical Spectroscopy*

*Frequency-Resolved Optical Gating: The Measurement of Ultrashort Laser Pulses*