Effects of high pulse intensity and chirp in two-dimensional electronic spectroscopy of an atomic vapor

The effects of high pulse intensity and chirp on two-dimensional electronic spectroscopy signals are experimentally investigated in the highly non-perturbative regime using atomic rubidium vapor as clean model system. Data analysis is performed based on higher-order Feynman diagrams and non-perturbative numerical simulations of the system response. It is shown that higher-order contributions may lead to a fundamental change of the static appearance and beating-maps of the 2D spectra and that chirped pulses enhance or suppress distinct higher-order pathways. We further give an estimate of the threshold intensity beyond which the high-intensity effects become visible for the system under consideration.


I. INTRODUCTION
Coherent two-dimensional (2D) femtosecond (fs) spectroscopy [1][2][3][4][5] is a powerful nonlinear spectroscopic technique which allows for the study of complicated dynamics and couplings in various quantum systems. Usually, 2D spectroscopy experiments are performed with laser intensities in the regime where perturbation theory holds (within the χ (3) limit) and hence the acquired signals can be described by the third-order response function formalism [6]. However, several new schemes have been developed, for example action-detected variants based on photoionization to study dilute gas-phase samples [7][8][9], 2D spectroscopy combined with microscopy to reach high spatial resolution [10][11][12] and higher-order techniques capable of providing information not accessible by third-order schemes [13][14][15][16][17][18][19][20][21]. In many of these implementations, higher laser intensities than usual may be reached (e.g. due to tight focusing in spatially resolved experiments), or may be of advantage (e.g. to enhance weak signals from dilute samples) or are even required (e.g. in higher-order spectroscopy schemes).
To the best of our knowledge, there is only limited literature available on the influences of high intensities on the appearance of 2D spectra and possible distortions/artifacts that may be induced by higher-order effects and saturation. In the theoretical study of Brüggemann et al. [22], non-perturbative calculations were performed which reproduce experimental 2D spectra of the Fenna-Matthews-Olson (FMO) photosynthetic antenna complex in the weak-field limit, i.e. where only a fraction of <10% of the FMO complexes is excited by every pulse. In the same study, they also performed simulations with higher field strengths, corresponding to the case of populating mostly the two-exciton state. They noticed a slightly broader and less structured spectrum, but the overall shape was conserved, from which they concluded that changing the laser intensity has only little influence on the spectra for such systems. In Ref. [23], 2D degenerate four-wave mixing measurements on atomic rubidium vapor are presented. In this work, three unexpected features that cannot be explained by third-order perturbation theory have been observed, which might be produced by fifth-order or cascaded third-order signals, but no clear conclusion about the origin of these peaks could be given. Furthermore, nonlinear two-phonon and two-photon interband coherences in InSb were investigated with 2D terahertz experiments [24]. As the optical interband dipole of InSb is exceptionally large, these experiments have been performed in the highly non-perturbative regime which has manifested itself in the impulsive off-resonant excitation of the interband coherences.
Very recently, a systematic theoretical study was published exploring the effects of intense laser fields on 2D electronic spectroscopy (2DES) signals using the concept of nonperturbative response functions [25]. These simulations based on a vibronic model system indicate the occurrence of peak shape modulations and phase shifts, as well as the enhancement of weak features. On the other hand, these results imply that experiments in the high-intensity regime could easily lead to a misinterpretation of 2D spectroscopic data.
Likewise, it is well known from coherent control experiments that chirped laser pulses can have a significant impact on coherent excitation schemes and population transfer [26][27][28][29][30][31][32][33][34][35][36][37]. Yet, the effect of chirped laser pulses on 2DES has not been much explored. Tekavec et al. studied peak shape distortion of 2DES for a two-level system interacting with a solvent caused by chirped pulses, both experimentally and via calculations of the third-order response of the system [38]. They introduced a chirp-correction scheme for 2D experiments with a continuum probe [39]. The impact of chirp on peak shapes in 2D spectra was studied analytically in Ref. [40].
In the current work, we systematically investigate, in a combined experimental and theoretical study, the effects of high laser intensities and chirp on 2DES. Our study is based on action-detected 2DES using phase-modulation/phase-cycling, but our results can also be transferred to non-collinear phase-matching based 2DES. As spectroscopic system, we chose atomic rubidium (Rb) vapor which provides us with a clean model system with well separated, narrow absorption features to facilitate the direct observation of peak distortions, amplitude modulation and unexpected features. This is in contrast to previous studies on systems with strongly broadened and overlapping features where subtle effects are much more difficult to identify [22]. The study of Rb vapors has recently also gained more interest due to the application of coherent 2D spectroscopy to ultracold samples [7]. Furthermore, the influence of the intensity of the pump and the probe pulses on transient absorption of atomic Rb vapor beyond the χ (3) limit has been investigated very recently in Ref. [41].

II. EXPERIMENTAL SETUP AND METHOD
A scheme of our experiment is depicted in Fig. 1. We employ a collinear action-detected type of 2DES for our experiments, that is, we are measuring the nonlinear population of our sample after the interaction with four fs laser pulses. More precisely, we exploit a phasemodulation approach combined with fluorescence detection and lock-in demodulation, as established by Marcus and coworkers [42]. Our experimental setup and measurement scheme are described in more detail elsewhere [7,8]. Very briefly, the output of a noncollinear optical parametric amplifier (NOPA) with 200 kHz repetition rate is converted into a collinear fourpulse sequence with adjustable delays τ , T and t. The laser pulses resonantly excite room temperature 87 Rb vapor contained in a glass cell and its fluorescence is detected with a photo-multiplier tube (PMT) perpendicular to the laser propagation direction.
The carrier-envelope phase (CEP) of each of the four pulses is modulated on a shot-toshot basis at distinct frequencies via acousto-optic modulators (AOM). Combined with the high repetition rate of the laser, this leads to a quasi continuous modulation of the relative CEPs between the laser pulses in the low kHz regime. Similar to phase-cycling approaches [43], the rephasing (RP) and non-rephasing (NRP) signal contributions (S RP and S NRP ) can be extracted from the total signal due to their specific phase signatures. Explicitly, the phase cycling conditions for the RP and the NRP signals are given by: Here, it is φ 21 = φ 2 − φ 1 and φ 43 = φ 4 − φ 3 . We choose the phase-modulation condition in a way that φ RP = 5 kHz and φ NRP = 13 kHz. Lock-in demodulation of the total fluorescence yield with respect to these two frequencies extracts then S RP (τ, T, t) and S NRP (τ, T, t). A Fourier transformation (FFT) of S RP (τ, T, t) and S NRP (τ, T, t) with respect to the time delays τ and t yields the complex-valued RP and NRP 2D frequency spectrã S RP (ω τ , T, ω t ) andS NRP (ω τ , T, ω t ), respectively. The sum of these two spectra gives the 2D frequency-correlation spectrumS C (ω τ , T, ω t ). Note that our definition of the NRP signal phase-signature is different from the standard definition: We account for the different convention in our data analysis by using accordingly adapted signs in the Fourier transformation of our data to obtain the proper rephasing, non-rephasing and correlation spectra. Throughout this work, we use the convention to plot the excitation frequency (ω τ ) on the x-axis and the detection frequency (ω t ) on the y-axis in the 2D spectra.
Data is recorded by scanning τ and t from 0 fs to 2960 fs in 80 fs steps. Furthermore, we apply zero-padding (factor 2) and multiply the time-data by a Gaussian window prior the Fourier transformation to reduce Fourier transform artifacts.
Measurements in the low-intensity regime are performed with the unfocused laser having a beam diameter (1/e 2 ) of ≈ 2.3 mm and a pulse energy of ≈ 33 nJ per excitation pulse.
Acquiring spectra beyond the perturbative regime is achieved by focusing the same laser pulses with a lens (f = 100 mm) down to a focal beam diameter of ≈ 44 µm. Furthermore, our pulses are stretched to ≈ 185 fs which corresponds to an estimated quadratic chirp of ≈ +1510 fs 2 . With this parameter, the peak intensity I 0 for the low-intensity measurement is about 8 MW/cm 2 and that for the high-intensity measurement ≈ 22 GW/cm 2 . For comparison, the intensities of the pump and probe pulses used in Ref. [41] to reveal high-intensity effects in transient absorption signals of Rb vapor were within 10-60 GW/cm 2 .
To collect the fluorescence and guide it to the PMT, we use a 4-f lens mapping with a spatial filter (SF) (iris diaphragm, opening diameter a ≥ 0.8 mm) implemented (Fig. 1).  Figure 2. Level-scheme, transitions and laser spectrum. (a) Energy level diagram for the relevant transitions in Rb. For the spectral bandwidth of the optical pulses, the system behaves as a fivelevel system composed of the ground state 5 2 S 1/2 , two lower-lying electronically excited states 5 2 P n (n = 1/2, 3/2) plus two higher-lying electronically excited states The reasons for using a spatial filter are (i) the suppression of stray light and (ii) to restrict the fluorescence detection volume to a small region around the focal spot. Theoretically, for a perfectly aligned detector, the detection volume should have a diameter of a along the beam propagation axis (see Fig. 1). This masking of the fluorescence is important to suppress signals originating from low intensity laser excitation occuring outside the laser focus volume.
We deliberately choose atomic Rb vapor as a target system. This gas-phase system provides well-known sharp spectral lines, while the level structure providing first and second excited state manifolds serves as an adequate model system for common molecular systems.
The relevant levels and allowed transitions of Rb atoms along with a typical laser spectrum used in the experiments are shown in Fig. 2. Note that the resolution of the present measurements does not allow to resolve the energy splitting between the two D m states (m = 3/2, 5/2). The transition wavenumbers and transition dipole moments used as input parameters for the numerical simulations of the 2D spectra are listed in Table I.  Table I. Transition wavenumbers and transition dipole moments of the Rb atom.ν ij is the transition wavenumber (taken from Ref. [44]) and µ ij the transition dipole moment in units of ea 0 , calculated using the transition probabilities listed in Ref. [45].

III. PERTURBATIVE DESCRIPTION OF ELECTRONIC 2D SPECTROSCOPY
A convenient way to keep track of the various signal contributions in 2D spectroscopy are double-sided Feynman diagrams. This diagrammatic description is based on perturbation theory and hence is only valid for adequately weak laser intensities. Nonetheless, it is helpful to discuss the basic features of our signals using Feynman diagrams. In Fig. 3 we show an exemplary selection of diagrams contributing to the fourth-order signal in collinear population-detected 2D spectroscopy using four excitation pulses and the phase-cycling con-    The factors behind the signs indicate the relative magnitude of the specific contribution (see text).
Pathways that beat with respect to T are marked by red dashed boxes.
radiative de-excitation of excited states can be neglected in the studied dilute sample. Hence, ESA features should appear as negative signals in the 2D spectra. Note that for better con- scheme has been multiplied by an additional factor of -1.
In Fig. 3, we also provide a rough estimate of relative amplitudes with which each pathway contributes to the fourth-order population-detected signal. To this end, we assumed a flat laser spectrum and used the dipole moments from Table I. All amplitudes are relative to the first RP pathway (labeled D 2 D 2 ), which is normalized to 1. Adding up all possible pathways, the strongest ESA feature is roughly a factor of 40-50 weaker than the D 2 D 2 feature. Hence, the ESA peaks should have negligible amplitudes in the measurements performed in the low-intensity regime.
It is also possible to make a statement about the T dependence of the signal with the help of the Feynman diagrams depicted here. Obviously, pathways that are in a coherent state (red dashed boxes in Fig. 3) during the time interval T , exhibit a coherent beating with the energy/frequency difference of the excited states P 1/2 and P 3/2 . This is the case for pathways contributing to off-diagonal spectral features for RP signals, and diagonal spectral features for NRP signals, respectively.

IV. NON-PERTURBATIVE SIMULATIONS OF THE 2D SPECTRA
To compare our high-intensity measurements with theory, we perform non-perturbative calculations of 2DES, which has the advantage that all orders of interaction are included automatically and hence the simulated spectra are not restricted to the 4th-order pathways presented above.
Our simulation procedure is similar to that used in Refs. [46,47]. It is briefly summarized as follows. The system Hamiltonian H involves five relevant Rb levels ( Fig. 2) and is thus The values of the transition frequenciesν ij = (E i − E j )/ and transition dipole moments Table I. The spectra of the laser pulses (in RWA) used in the simulations are approximated by the expression where t FWHM is the full width at half maximum (FWHM) pulse duration, while χ 2 is the quadratic chirp and ω L = 2πc/(787 nm) denotes the laser central angular frequency (see Fig.   2 for the comparison of the experimental and simulated E(ω)). The complex pulse envelopes For Rb vapor at room temperature, no relaxation/dephasing occurs on the time scale (3 ps) of the experiment [48]. Hence the time-dependent Schrödinger equation is appropriate for the description of the photo-induced system dynamics. In matrix-vector notation, it reads In our simulations, Eq. (4) is solved numerically for all given inter-pulse delay times and phase angles by a fourth-order Runge-Kutta integrator with a time step in the range of 0.5 fs to 10 fs depending on the used pulse intensity (see [46,47] for the simulation details). The atoms are assumed to be in the ground state before interaction with the laser pulses. Hence the initial condition is Here τ 0 is a time moment before the arrival of the first laser pulse (in the simulations, we The detected signal is then proportional to the weighted sum of the excited-state asymptotic populations, Here τ ∞ is a time moment after the arrival of the fourth laser pulse (in the simulations, we take τ ∞ = τ 4 + 10t FWHM ), while 0 ≤ Γ ≤ 2 is the weighting factor which quantifies contributions of the higher-excited states to the signal (see Refs. [42,[49][50][51] for an in-depth discussion). In the present case, no non-radiative processes occur and we set Γ = 2 in all the simulations.
In the present experiment, the RP and NRP contributions S RP and S NRP to the signal are discriminated by the phase-modulation procedure [42]. This procedure can be explicitly implemented in perturbative computer simulations of 2DES signals [52,53]. The explicit realization of the phase-modulation in non-perturbative simulations of 2DES signals is, however, numerically inefficient since it requires a long-time propagation of the driven Schrödinger equation (4). Instead, we extract S RP and S NRP from the total signal S by using a discrete Fourier transformation (this [54,55] and similar [56,57] numerical methods have widely been used in the literature).
The procedure goes as follows. Due to isotropy of space, only relative phases of the pulses matter. Hence all φ a can be shifted by e.g. φ 1 to yield the modified wave vectors φ a ≡ φ a − φ 1 and the appropriately modified phase-matching conditions (1) and (2). The phase-matched signals S RP and S NRP can therefore be computed as follows: Here the upper (lower) signs in the exponential function correspond to the RP (NRP) contributions andφ n ≡ 2πn/N φ . Eq. (7) becomes exact when N φ → ∞ and summations are replaced by integrals. If the system-field interaction is not too strong, much smaller values of N φ are sufficient. In particular, N φ = 3 yields the so-called 3 × 3 × 3 phase cycling scheme, which is known to be exact for weak-field fourth-order signals [43]. In the context of action-detected 2DES spectroscopy, this procedure has been employed by several experimental groups [58,59], while a 5 × 5 × 5 phase cycling scheme was used to extract sixth-order contributions [21]. The simulations of the present work are carried out via Eq. The radial intensity profile in the laser focus follows a Gaussian distribution, which means that the signal arises from a superposition of different laser intensities. To account for this effect, we discretize the expected radial intensity profile at the focus position by calculating the intensity at 15 radial positions r 1 to r 15 (Fig. 4). We choose these positions to be equally spaced by ∆r and range from r 1 = 0 µm to r 15 = 2.1ω 0 , with ω 0 being the 1/e 2 -radius of the laser beam in the focus. For each of these intensities, a numerical simulation of the 2D signal is performed yielding the RP spectraS RP,i (ω τ , T, ω t ) and NRP spectraS NRP,i (ω τ , T, ω t ) with i = 1, 2, . . . , 15. These spectra are then added up considering the weighting factors A i (see Fig. 4) to obtain the intensity weighted spectrã  andS with and for i = 2, 3, . . . , 15.
We neglect variations of the intensity along the laser propagation direction (pointing into the drawing plane of Fig. 4) due to the use of a spatial filter in the fluorescence detection unit (see Sec. II).

V. RESULTS AND DISCUSSION
In the following we present and discuss the effects of high pulse intensity we observed in 2DES experiments. Section V A deals with the change of the static appearance of 2D spectra when going to the highly non-perturbative intensity regime. It is shown that due to higher-order contributions new peaks can appear and peak shape distortions of the GSB/SE features are possible. In section V B we show by non-perturbative simulations that chirped pulses have a strong influence on 2D spectra in the high-intensity regime as distinct higherorder pathways can be enhanced or suppressed. We further discuss the effects of high laser intensities on the temporal dynamics reflected in the 2D signals in section V C. We show that higher-order processes can significantly influence the coherent beat behavior of the 2D maps.
Finally, in section V D, we present how we define an estimate of the threshold intensity at which high-intensity effects start to play a role for the system under investigation.
A. Perturbative regime vs. highly non-perturbative regime In this section, we study the changes observed in our 2D spectra when going from intensities assigned to the perturbative regime to intensities far beyond this regime. For this purpose, Fig. 5 shows experimental spectra, for a fixed population time of T = 360 fs, acquired using the unfocused laser beam (exp. low I) and data measured by focusing the laser onto the sample (see Sec. II) (exp. high I). The data of the high-intensity experiment are compared with a numerical simulation including the pulse chirp estimated to be present in the experiment (sim. high I). Note that, in the experimental data, saturation of the detector and detection electronics may also lead to the occurrence of intensity-dependent saturation effects stemming from the nonlinear mixing of the linear signal response. The absence of such artifacts was carefully checked. The dashed lines in the spectra indicate the transition frequencies corresponding to all possible transitions in our target system (see Table I). Since the frequency resolution of the present measurements is not sufficient to resolve the frequency difference between the two transitions P 3/2 ↔ D 3/2 and P 3/2 ↔ D 5/2 , a single reference line is plotted at the arithmetic mean of the two transition frequencies.
Furthermore, to speed-up data acquisition and accompanied simulations (see below), large increments (80 fs) for τ and t were chosen. As a consequence, the highest frequency features in the spectra (P 1/2 ↔ D 3/2 ) are aliased and projected onto lower frequencies (green dashed lines). Table II lists all crossing points and hence the positions where signals may appear.
We categorize these signals into 4th-order and higher-order signals according to the scheme shown in Fig. 6.
As expected, the low intensity spectra (upper row Fig. 5  probe trans. label pump trans. probe trans. 4th-order Feynman diagrams (GSB/SE) with absent ESA contributions due to their much smaller amplitude. To confirm that this measurement can be adequately described by 4thorder perturbation theory (and does not necessarily require non-perturbative treatment), we performed a perturbative simulation (not shown here) using the response functions formalism analogous to the procedure described in [42]. The sign (positive) and the shapes of the  Figure 6. Sketch to categorize expected peaks for 4th-order (yellow) and 6th-order (red) processes.
The diagram shows the expected peak positions and peak categorization. Peaks are labeled A-P, dashed lines correspond to the expected transition frequencies (see Fig. 5). Circles: 4th-order GSB/SE peaks. Squares: 4th-order ESA peaks. Red triangles: peaks allowed only in ≥6th-order.
Yellow and red shaded symbols: peaks with overlapping 4th-order and ≥6th-order contributions.
features in these simulated 2D spectra match very well with the experimental spectra. In the reminder we use this intensity condition as the reference for 2D spectra in the perturbative regime.
The spectra obtained for high laser intensities (middle row in Fig. 5) show a different behavior. The most prominent differences are highlighted by black circles. One consequence of the high intensity is a change in the amplitude ratio of the diagonal peaks G and M. This is visible in the RP, NRP and the correlation spectrum. Furthermore, in the RP spectrum, the off-diagonal feature E changes sign whereas O does not. Hence, an asymmetry is introduced, which is not present in the perturbative regime. For the NRP spectrum this effect is absent.
Additionally, new features appear at positions J and B in the RP spectrum and at A, B and F in the NRP spectrum. For other population times T (not shown here) also a feature at position A can be clearly identified in the RP spectrum. The correlation spectrum C combines all the above mentioned changes, as it is the sum of the RP and NRP spectra.
However, the additional features, only appearing at high intensities, are less pronounced in the correlation spectrum.
The appearance of the additional peaks when using intense laser pulses can qualitatively be described by Feynman diagrams beyond 4th-order. Such pathways include more than one interaction with a single laser pulse. For example, 6th-order pathways fulfilling the phase-cycling conditions of Eq. (1) and Eq. (2) can be constructed by three interactions with one pulse and one interaction with each of the other pulses in the following way: Hence, these higher order signals are not filtered by phase-cycling and at high intensities may "leak" with significant contributions into the 4th-order signal detection. In analogy, >6th-order signals may be constructed that also leak into the 4th-order detection. Due to their nonlinear intensity dependence their amplitudes are expected to be small and are not considered here for simplicity. It is important to note that, although demonstrated here for phase-cycling, the analog case accounts for phase-matching based 2D spectroscopy experiments. Fig. 7 depicts several representative 6th-order Feynman diagrams to explain the highintensity spectra. Many more 6th-order pathways exist, but are omitted here for simplicity. (b), (e) and (f). This is attributed to a possibly smaller net amplitude of these contributions due to the involved combinations of dipole moments as well as chirp effects (see below). We want to further point out that for peaks F and P, the system has to be in a different P n state during the time interval τ than during t. This is possible for 6th-order NRP diagrams (b) but not for 6th-order RP diagrams (a) (here only ≥8th-order diagrams contribute to F and P). This nicely matches to the fact that the NRP spectrum shows a peak at F and the RP spectrum does not. Diagrams (c) and (g) correspond to peak positions (A, C, I, K) where 6th-order contributions can overlap with 4th-order ESA contributions (yellow and red shaded squares). Hence, with the present data, we can not distinguish if the appearance of peak A in the NRP spectrum is caused by 4th-order, 6th-order or by both contributions. Due to the different sign, the overlap of 6th-order contributions (diagrams (d) and (h)) with 4th-order contributions can in prinicple lead to sign and shape changes of the GSB/SE features (E, G, M, O). We observe this effect for peak E in the RP high-intensity spectrum. Using this argumentation, in principle also the off-diagonal peak O should change its sign and shape and the same behavior should be visible in the NRP spectrum (see, e.g. the simulated 2D signals in Ref. [25]). However, this is not the case for our data. We think that this asymmetry is due to the fact that our laser pulses exhibit a significant chirp. Numerical simulations (see below) show that chirp has a strong influence on the 2D spectra in the high-intensity regime, especially on higher-order pathways.
While Feynman diagrams offer an intuitive understanding of 2D spectra, they are limited to a perturbative description of nonlinear signals. Therefore, we performed non-perturbative numerical simulations of the high-intensity 2D signals (lower row in Fig. 5). We obtain a good agreement between experiment and simulation. The most prominent features are reproduced, although the magnitude and appearance of the high-intensity effects differ to some extent from the experimental data. All additional features corresponding to >4th-order pathways are reproduced in the simulations. Also the sign change of peak E is apparent in the calculated RP spectrum. However, the magnitudes of the peaks B, F, J are somewhat overestimated by the simulations. The most striking discrepancy between simulated spectra and experimental spectra is the absence of the strong change in the relative ratio of the diagonal peaks M and G in the simulation.
These differences are attributed to the non-trivial spatial intensity distribution in the laser focus and temporal chirp of the laser pulses (see below). Although we account for these effects in the simulation, it is difficult to match these parameters precisely with the experiment.
The chirp was estimated from the measured time-bandwidth product of the pulses without directly measuring the spectral phase of the pulses. The dispersion of the focusing lens and glass walls of the vapor cell were included analytically. The laser intensity can be misgauged if the 4-f lens system is not precisely aligned to the laser focus. Furthermore, the spectral intensity of the laser at the atomic resonances was estimated by assuming a Gaussian spectral profile which deviates to some extent from the experimental laser spectrum (Fig. 2 (b)).
To set the low and high intensity measurements into perspective of saturation and chirp effects, we have calculated the excitation probability for a single laser pulse exciting the five-level system (5LS) of the Rb atom (Fig. 8). The high intensity measurements (chirped pulses) are performed at a peak intensity of ≈ 22 GW/cm 2 (black dashed line in Fig. 8 (b)).
In contrast, the low intensity measurements (chirped pulses) were performed at a much smaller peak intensity of ≈ 8 MW/cm 2 (not shown). For a two-level system (2LS) excited by a transform-limited laser pulse, it is well known that the population probability of the excited state, after the interaction with the laser pulse, exhibits an oscillatory behavior as a function of the laser intensity which is the result from Rabi oscillations [60]. Similarly, in a more complex system, such as the 5LS studied here, the transitions exhibit a clear oscillatory behavior which is, however, changing in frequency and amplitude as a function of the laser intensity ( Fig. 8 (a)  the behavior is completely different for chirped pulses, similar to what was found in Refs. [61,62]. Here, we observe some strong damping of the excited-state population oscillation by the chirp. Note that Fig. 8 shows the population after the interaction with a single laser pulse. For the case of a four-pulse interaction, the dependency on the intensity and the chirp is expected to be more complicated. We therefore expect that errors in the estimates of the experimental parameters chirp and laser intensity are responsible for the mismatch between experiment and simulation.
In this analysis, the Feynman diagrams provided us with an intuitive explanation of the origin of the additional peaks and peak shape distortions appearing at high intensities. In addition, the numerical non-perturbative simulations provide us a quantitative confirmation of these assumptions and show that we can reasonably well simulate the high-intensity effects in the 2D experiments.

B. Influence of chirp in the high-intensity regime
To further investigate the influence of chirp on the 2DES signals, we performed simulations at low and high laser intensities for chirped and transform-limited pulses. At low intensities, we do not observe a significant impact of chirp on the 2D spectra, despite the fairly large amount of the applied chirp (not shown). In contrast, at high laser intensities we observe several differences between the cases of chirped and unchirped pulses (Fig. 9).
Note that instead of sampling the τ and t delay with a 80 fs step size, as done for the data presented in Fig. 5, here a step size of 30 fs is used to avoid aliasing of the P 1/2 ↔ D 3/2 features to lower frequencies (fully-sampled transition freq. in Fig. 9: ≈ 13122 cm −1 , aliased freq. in Fig. 5: ≈ 12700 cm −1 , indicated by green dashed lines). Due to the smaller step size, Fig. 9 shows also a considerably larger frequency range.
Comparing the two simulations it becomes apparent that several spectral features are amplified/damped when introducing chirp (e.g. peaks A, B, C, D, I, J, K, L, N and F).
Consequently, it is possible to enhance or suppress distinct higher-order pathways to a certain degree by choosing properly chirped pulses. In combination with the beating analysis (below) our study indicates that the influence of chirp increases for >4th-order pathways. This is expected, since all >4th-order signals observed with our phase-cycling condition involve multiphoton transitions (multiple interactions with a single laser pulse). For such processes, the strong effect of chirp is in principle well-known from coherent control and population transfer schemes [26][27][28][29][30][31][32][33][34][35][36][37]. As a second observation, the asymmetry in shape, sign and amplitude, induced between the two off-diagonal peaks E and O using intense  to its fully-sampled frequency (≈ 13122 cm −1 ). In all spectra, the most intense peak is normalized to one. chirped pulses, is not present when unchirped pulses are employed. Hence, this strong chirp dependence could explain why peaks B, F and J are stronger and the asymmetry between the off-diagonal peaks E and O is more pronounced in the simulated spectra than observed in the experimental data (Fig. 5).

C. Investigation of coherent beatings
Besides the alternation of peak shapes and amplitudes, it is also important to consider the effects of high laser intensities on the temporal dynamics reflected in 2D spectra. To this end, we studied experimentally and theoretically the dependence of the 2D signals on the population time T in the range from 0-360 fs in 30 fs steps. Exemplarily, the time evolution of the amplitude of the diagonal (G, M) and off-diagonal (E, O) peaks of the real part of the RP and NRP spectra are shown in Fig. 10. The plotted signal corresponds to the amplitude of the pixel nearest to the theoretical transition frequency in the experimental/simulated 2D spectra.
The low-intensity case follows the predictions from 4th-order Feynman diagrams (Fig.   3). An amplitude beating of the off-diagonal peaks is present in the RP contribution, whereas the diagonal peaks exhibit no time dependence. The opposite behavior is true for the NRP spectrum. The frequency of the present beating matches to the predicted value corresponding to the energy difference between the P 1/2 state and the P 3/2 state.
In the high-intensity regime, the dependence of the signals on the population time changes. Here, also the diagonal peaks show an oscillation with respect to T in the RP spectrum. The time period of this beating is the same as for the off-diagonal peaks (note that the T -evolution in both low-and high-intensity regime is governed by the field-free system Hamiltonian), but appears with a ≈ π-phase shift. We attribute this new beating to 6th-order contributions as qualitatively represented by 6th-order Feynman diagrams (Fig.   11). Here, RP pathways exist, leading to signal at the diagonal positions M and G, in which the system is in a coherence between the P 1/2 state and the P 3/2 state during T . A typical diagram for feature G is depicted in Fig. 11 (a). We further note that the change of the prefactor of the perturbative expansion from i 4 = 1 for 4th-order pathways to i 6 = −1 for 6th-order pathways implicates a phase shift of π [21]. Hence, also the observed phase shift between the diagonal and off-diagonal peak beatings nicely agrees with our assumption of 6th-order contributions overlaying with 4th-order pathways. Note that intensity-dependent phase shifts in the time dynamics have been also detected for transient absorption signals in Ref. [41].
For the NRP spectrum, an analogous argumentation applies. Here, 6th-order contributions (e.g. Fig. 11 (b)) induce oscillatory pathways for the off-diagonal peaks E and O which are not allowed in 4th-order. However, just one of the two off-diagonal features (E) shows the predicted oscillation. The phase of this beat is again roughly the opposite of the features which oscillate due to the 4th-order contributions, at least compared to the diagonal feature G. Our interpretation is confirmed by numerical simulations (Fig. 10 (e,f)) where we find exactly the same oscillatory behavior and phase shift.   Figure 11. 6th-order contributions influencing the T -dynamics. Exemplary 6th-order diagrams being in a coherent state during T (highlighted red) causing a coherent beating of the diagonal peaks in the RP (a) and the off-diagonal peaks in the NRP (b) spectrum at high intensity.
We want to point out that we observe a non-negligible phase shift between the beatings of the two diagonal peaks (blue and cyan curves) in the low-and high-intensity measurements (Fig 10 (b), (c), (d)). We attribute this to be an effect of the data analysis. In the analyzed RP and NRP spectra, peaks comprise of positive and negative parts oscillating out of phase.
The phase of the deduced beat thus shows a dependence on the position of the finite-size peak integration area. To reduce this effect, we minimized the integration area to a single pixel of the the 2D maps. On the other hand, increasing the integration area to cover the whole peak would avoid this issue. However, in this case we observe a significant influence on the beat behavior induced by the tails of neighbouring peaks. Our beat analysis thus implies, that even for sparse 2D frequency spectra as analyzed here, precise retrieval of beat phases is difficult. The situation improves for truly absorptive peaks in 2D correlation spectra, yet, flanking ESA features may lead there to a similar phase issue.
Above, we discussed that chirped pulses cause an asymmetry between the off-diagonal peaks E and O in the high-intensity regime. This suggests that the difference in the time dynamics of these two features observed here may be also due to chirp effects. To check this, we repeated the simulation of the T -scan with high intensity, but this time assuming transform-limited laser pulses ( Fig. 10 (g,h)). The pulse energy is kept the same as for the chirped case. As expected, now also a beat of the feature O shows up in the NRP spectrum, having the same frequency, amplitude and phase as the other off-diagonal feature E.
We hence find, that higher-order processes superimposing with 4th-order signals can significantly alter the beat behavior of 2D spectra in amplitude and phase. Chirp can in addition change the beat properties, especially if high-order signal contributions are involved.

D. Intensity threshold for higher-order effects
In the present work, we performed 2DES measurements and simulations of a model system at two distinct laser intensities, on the one hand, at a very low intensity in the linear excitation regime and on the other hand, at a much higher intensity where the excitation probability clearly behaves in a nonlinear way and shows Rabi oscillations (Fig. 8). We find good agreement between experiment and simulation. This allows us to use our numerical simulation to explore in more detail the threshold intensity at which distortions due to higher-order effects occur in 2D spectra.
To this end, we performed numerical simulations of 2D spectra for several laser intensities for a fixed population time of T = 300 fs. To focus on the transition with the highest transition dipole moment where nonlinearities should be visible first, we shift the carrier frequency from the experimental value (≈ 787 nm) to the S 1/2 → P 3/2 resonance (780 nm).
To provide insight how the chirp and the complexity of the target system influence the onset of high-intensity effects, the simulations have been performed for the following cases: (i) the complete five-level system (5LS) excited by transform-limited pulses, (ii) the 5LS excited by pulses with a quadratic chirp of +1510 fs 2 , (iii) just a two-level system (2LS) consisting of ground state S 1/2 and excited state P 3/2 excited by unchirped pulses, (iv) the 2LS excited by chirped pulses. We have chosen the input intensities for this simulation series such that the population probability of the P 3/2 state, in a pure 2LS, interacting with a single transform limited pulse, exhibits the following values: 1% and 10-100% in 10% increments ( Fig. 12 (a)). These values are calculated with the following formula for the excited-state population which is derived from first order perturbation theory adopting the rotating-wave approximation (RWA) and assuming a resonant laser pulse described by Here, µ eg is the transition dipole moment from the ground (S 1/2 ) to the excited state (P 3/2 ), have been used for the Rabi and the perturbative calculation (Eq. (16) and (14)), but this time plugging in the pulse duration and peak intensity of the chirped pulse. Obviously, these simple formulas fail to estimate the population for strongly chirped pulses, even for very low intensities. Fig. 12 reveals the intensity regime where perturbation theory starts to deviate from the exact numerical simulation, which implies at which intensity the experiment cannot be described by low-order perturbation theory anymore. Note that measuring the fluorescence saturation curve would be an experimental way to estimate this onset regime.
We further note, that the values calculated here are valid for the specific system studied here.
Yet, our simulations may provide a rough general estimate for the threshold intensity below which perturbation theory is valid. This conclusion is also corroborated by the observation that the breakdown of the perturbative description of transient absorption of pump-probe signals in the displaced harmonic oscillator model starts at µ eg 2I 0 /( 0 c) > 0.01 eV [55].
For the transition dipole moment of the S 1/2 → P 3/2 transition, this yields the threshold field intensity of 0.3 GW/cm 2 , in good agreement with Fig. 12. Note also the different behavior for a 2LS and a 5LS at higher intensities, implying that for more complex systems (e.g. molecular systems) the situation may change. However, at low intensities, both systems behave very similar, suggesting that our conclusions may be generalized to more complex systems. We also point out that we consider here population probabilities (two field-matter interactions exciting a population), which corresponds to the linear absorption of the system, a quantity that is readily experimentally accessible. In 2D spectroscopy it is also common to think in terms of single field-matter interactions. For a 2LS this probability corresponds to the square-root of the population probability.
While Fig. 12 considers the interaction with a single laser pulse (i.e. two field-matter interactions), for 2DES experiments it is more relevant at which intensities deviations are observed in the nonlinear system response generated by four field-matter interactions. To define the threshold intensity at which higher-order effects appear, we investigated the distortions in the corresponding 2D spectra. To do so, we normalized, for all the simulated cases (i)-(iv), the real parts of the different contributions (RP, NRP, C) for each intensity to have the maximum amplitude of one. The spectra simulated for the intensity corresponding to a population probability of 1% after one pulse is chosen as reference spectra for the perturbative regime. Our simulations show that distortions are best visible in NRP spectra, for which the onset of distortions are shown in Fig. 13   for the 5LS.
For the simple model system of Rb atoms yielding clean, well-separated and sharp spectral features, distortions due to high pulse energies can be directly identified by a qualitative comparison of the spectra with the low-intensity reference spectrum. Here, we find that first visible changes appear for the off-diagonal features. For the given pulse energy in Fig. 13 (≈ 0.15 nJ), these changes are more pronounced for the case of chirped pulses, although the population probability (numerically exact calculation) of the strongest transition induced by one pulse is nearly the same for both cases (≈ 18% assuming unchirped pulses and ≈ 17% assuming chirped pulses). For population probabilities of 10%, no distortions can be identified. Hence, we conclude that the limit of the perturbative intensity regime lies in the range of 10-20% excitation probability of the P 3/2 state, corresponding to an intensity of ≈ 0.4 − 0.8 GW/cm 2 for transform-limited pulses. Note, that the overall signal scaling (integrated absolute spectra) deviates already at lower intensities from the I 2 dependency expected for 4th-order processes. However, we did not analyze this in detail as our focus was more on distortions of the spectral appearance.  For comparison, Fig. 14 shows corresponding simulations for a 2LS (cases (iii) and (iv)).
Here, we find that even for the highest tested pulse energy (≈ 0.73 nJ) in this study, which delivers about 71% excited-state population for transform-limited pulses and 65% for chirped pulses (numerical exact calculations), the feature in the 2D spectrum does not change much.
Hence, for a pure 2LS, distortions due to high laser intensities and chirp seem to be absent and we conclude that distortions arise from the mixing of multiple states during the nonlinear multi-pulse excitation in 2DES experiments.
For a more quantitative analysis we calculated the deviation of the 2D spectra by subtracting the reference spectrum (normalized to the maximum absolute amplitude of one) from the higher intensity spectra, normalized in the same way. The maximum absolute value of the difference spectrum as a function of pulse energy for the cases (i)-(iv) is shown in Fig. 15. As already indicated in our qualitative analysis, the deviation in the NRP spectrum increases more strongly than in the RP spectrum and chirping the pulses results in an increased maximal deviation for the 5LS. Furthermore, the maximal difference to the reference spectrum assuming a 2LS does not show such a drastic increase with intensity. the results for the 5LS assuming TFL pulses, the red squares for the 2LS assuming TFL pulses, the dark blue points for the 5LS assuming chirped pulses and the dark red points for the 2LS assuming chirped pulses.
In conclusion, our analysis implies that observable intensity-induced distortions of the 2D spectra start to appear at excitation probabilities of 10-20% for a single pulse. This is consistent with other work and may be generalized to more complex systems as a rule of thumb. We note, that at a population probability of 20%, the linear absorption of the sample hardly deviates from a linear curve (≈ 10% deviation for the case of no chirp) (see Fig. 12 (a)). Such small deviations are usually difficult to determine in the experiment.
Therefore, care has to be taken in determining the critical intensity in each experiment.

VI. CONCLUSIONS
In this work, we investigated effects caused by high intensities and chirped laser pulses on 2D electronic spectroscopy. To this end, we performed collinear phase-modulated 2DES experiments of Rb atoms, which serve as a well-defined model system with well-known and sharp features, in different intensity regimes and compared the results with non-perturbative numerical simulations. The advantage of having sharp and well-separated spectral features in the atomic gas-phase system facilitated the observation of subtle changes in the spectra, such as amplitude modulations and additional peaks appearing at high intensity. On the other hand, the absence of dephasing on the observation time scale (3 ps) required the application of a window function (here Gaussian) to the data to avoid Fourier transform artifacts [1], resulting in artificial lineshapes with no physical meaning. Hence, the investigation of lineshape effects [38] was not possible with our data.
We found that even though we use a high-repetition rate laser featuring low pulse energies (≈ 30 nJ), high intensities in the regime of multiple Rabi-cycles are already reached by focusing the laser onto the target system with a moderately short focal length of f = 100 mm (Fig. 8). As a consequence, new features appear in the 2D spectra which can be explained by 6th-order contributions that cannot be discriminated by phase-modulation/cycling or phase-matching. Moreover, higher-order contributions also influence the peak shape, sign and amplitude of the features associated with GSB/SE pathways (diagonal and off-diagonal peaks). We further found that also the beating behavior of these peaks changes due to superimposing higher-order signals. The different beating behavior for diagonal and off-diagonal features in RP and NRP 2D spectra, respectively is commonly used to separate coherent from incoherent dynamics [65]. We find, that at high laser intensities the beat behavior of RP and NRP 2D spectra intermix, which in principle compromises such analysis. Besides the change in the relative ratio of the diagonal GSB/SE features, all the experimentally observed effects are qualitatively reproduced by non-perturbative simulations. Although shown here for phase-cycling, analogous effects are expected for phase-matching experiments.
To investigate the influence of chirp in the high-intensity regime, we compared simulations for chirped (+1510 fs 2 ) and transform-limited pulses. Here, we found that it is possible to enhance or suppress distinct higher-order pathways using chirped pulses. Furthermore, an asymmetry between the off-diagonal peaks, as observed in the experimental data, is present when using chirped pulses, but vanishes for transform-limited pulses.
As we observe all these effects already at rather low pulse energies and moderate focusing conditions, this raises the question at which intensity regime the onset of high-intensity effects occurs in 2D spectra. To this end, we performed simulations for several laser intensities. We found that the onset of high-intensity effects is best visible in the NRP contribution and manifests itself primarily in changes of the off-diagonal features. For a given pulse energy, the changes are more pronounced in the case of chirped pulses. For transform-limited pulses, we conclude that the high-intensity limit lies in the range of ≈ 0.4 − 0.8 GW/cm 2 , corresponding to an excitation probability per pulse of 10-20% of the P 3/2 state (transition with the highest transition dipole moment). While these estimates are valid for a specific system (Rb atom as five-level system), they may be used as a general qualitative estimate of the threshold intensity below which the spectral appearance of 2D spectra can be described by perturbation theory.