We present a new data acquisition concept using optimized non-uniform sampling and compressed sensing reconstruction in order to substantially decrease the acquisition times in action-based multidimensional electronic spectroscopy. For this we acquire a regularly sampled reference data set at a fixed population time and use a genetic algorithm to optimize a reduced non-uniform sampling pattern. We then apply the optimal sampling for data acquisition at all other population times. Furthermore, we show how to transform two-dimensional (2D) spectra into a joint 4D time-frequency von Neumann representation. This leads to increased sparsity compared to the Fourier domain and to improved reconstruction. We demonstrate this approach by recovering transient dynamics in the 2D spectrum of a cresyl violet sample using just 25% of the originally sampled data points.

## I. INTRODUCTION

Coherent two-dimensional (2D) electronic spectroscopy^{1–5} has become an established method to study, among others, energy transfer in light harvesting complexes,^{6–8} electronic dynamics in semiconductors,^{9–13} plasmonic coherences,^{14} and photochemical reactions.^{15} Most of the experimental realizations rely on the non-collinear box geometry and heterodyned detection of a coherent optical signal. Another option that has been developed is action-based 2D spectroscopy using a train of four (often collinear) laser pulses and a phase-modulation,^{10,16} or phase-cycling scheme.^{17,18} The action-based ansatz employs the detection of incoherent signals, e.g., fluorescence^{16,17} (see also Fig. 1(a)), (photo)-electrons,^{10,11,13,14} or (photo)-ions.^{19} One advantage of detecting incoherent signals is that it permits the study of low-density samples, even down to the single-molecule limit.^{20} A detailed comparison of various experimental 2D spectroscopy approaches can be found in a recent review.^{21}

Multidimensional spectroscopy in general requires multidimensional scanning of parameters such as time delays or phases. While single-shot techniques have been developed,^{22,23} the question remains whether the amount of data to be acquired can be reduced while retaining the essential information content of the resulting 2D spectrum. In the present work, we focus on the action-based approaches for which the scanning protocol requires the sampling of the coherence time *τ* and the signal time *t* with a rate obeying the Shannon-Nyquist sampling theorem.^{24} A 2D Fourier transform with respect to these two delay times yields a 2D spectrum for a given population time *T*. Especially for gaseous samples, coherences are long-lived^{25} and hence require large time delays for the two indirect dimensions. Furthermore, it may also be useful to perform an additional Fourier transform with respect to the population time *T* and generate a quasi-3D spectrum in order to distinguish the origin of coherences^{26,27} or different channels in molecular reactions.^{15,28} Fully coherent fifth-order 3D spectroscopy methods have also been demonstrated.^{29–32} Thus, in action-based multidimensional spectroscopy, the time-domain parameter space can become very large, reaching easily hundreds of thousands of different configurations. This can be efficiently reduced and hence the acquisition time can be substantially decreased by performing the measurement in a rotating frame,^{33,34} which shifts the coherent oscillations from optical frequencies (close) to the origin of frequency space. Despite this, the acquisition time in these experiments remains large enough to warrant searches for further data acquisition improvements.

Initially developed in the signal processing community, “compressed sensing”^{35,36} (CS) has evolved as a new paradigm for reducing the amount of sampled data. The growing interest in CS can be attributed to its successful broad applicability in several areas of research and engineering such as imaging using a single-pixel detector,^{37,38} nuclear-magnetic resonance imaging,^{39} NMR spectroscopy,^{40,41} and more.^{42} In general, CS can be used to recover certain signals with a significantly reduced number of sampling points. A necessary condition for the successful reconstruction of undersampled data is the existence of a transform domain in which the signal has a sparse representation, i.e., in which only a small percentage of basis coefficients is non-zero. This sparsity, in combination with a random sampling protocol, can be exploited in order to achieve a high-fidelity reconstruction.^{43} First applications of CS to coherent 2D spectroscopy show the suitability of random undersampling to efficiently reduce the amount of time-domain data points without losing spectral information either in experiment^{44} or in theoretical calculations.^{45} Furthermore, it was shown that CS reconstruction can also be used to retrieve the 2D spectrum from a uniformly sampled subset of the full signal.^{46} Other non-uniform sampling protocols, such as exponentially biased, radial, or Poisson-gap sampling, are successfully applied in multidimensional NMR.^{47–50} They benefit from *a priori* assumptions on the approximate signal characteristics and use a dense sampling close to the origin of the time domain where the signal amplitudes are still large.

In the present work, we use optimal control concepts to optimize the sparse sampling matrix for best performance at a given undersampling percentage. In combination with rotating-frame data acquisition, this allows us to decrease the measurement time substantially without significant loss in information content of the resulting 2D spectra. The experimental realization is made using a customized pulse shaper that is able to scan an arbitrary sequence of sampling points in the *τ*-*t* parameter space on a shot-to-shot basis. In contrast to setups based on mechanical delay stages, our pulse-shaper approach allows us to realize solely the desired sampling points while retaining the regular duty cycle. Additionally, we use a compact representation of 2D spectra generated by the von Neumann (vN) basis^{51} which increases the sparsity in the transform domain of the CS reconstruction algorithm.

## II. MATERIALS AND METHODS

### A. Concept

The general idea of our new approach is illustrated schematically in Fig. 1. As the first step (Fig. 1(a)), we measure a reference data set at a fixed population time *T* using coherent 2D fluorescence spectroscopy. For this, we apply the full time-domain sampling protocol in a rotating frame and retrieve the signal contribution of interest after phase cycling.^{18} Next, we generate a random sampling matrix based on the temporal grid defined by the previous measurement with the desired undersampling percentage. This is shown in Fig. 1(b), with the sampled data points being indicated by white squares. This matrix is used to undersample the reference data set *a posteriori*. We then employ a genetic algorithm to optimize this reduced sampling matrix such that the reconstruction of the undersampled reference 2D spectrum agrees as well as possible with the fully sampled reference 2D spectrum. In the present work, we particularly focus on rephasing (photon-echo) and nonrephasing (reverse transient grating) signals; however, treatment of others such as two-quantum coherence is also conceivable. Data reconstruction is performed using the established “two-step iterative shrinkage/thresholding” (TwIST) algorithm^{52} and a new, joint time-frequency von Neumann representation of 2D spectra. Finally (Fig. 1(c)), we use the optimal sampling matrix to measure and reconstruct 2D spectra at all other population times. The savings in data acquisition time arise from the reduced amount of required data for all population times other than the reference set. In this paper, we benchmark our approach experimentally using the laser dye cresyl violet.

### B. Rapid-scan 2D fluorescence spectroscopy

As a basis for the experiment, we employ our recently developed setup for action-based fluorescence 2D spectroscopy using shot-to-shot femtosecond pulse shaping.^{53} Briefly, we take 1 mJ of the output of a commercial Ti:Sa amplifier system (Spitfire, Spectra Physics) with 120 fs pulse duration and 1 kHz repetition rate to pump a noncollinear optical parametric amplifier (NOPA) which delivers pulses centered at 590 nm (25 nm FWHM). After the NOPA, the pulses pass a single-prism compressor^{54} to account in part for the positive group-delay dispersion introduced by the following commercial pulse shaper (Dazzler, Fastlite). The shaped output pulses are focused into a 250 $\mu m$ glass capillary with an f = 150 mm focusing mirror. Using a standard 0.25 NA microscope objective (CVI Melles Griot) we collect the sample fluorescence and detect it with an avalanche photodiode (A-Cube S5003, Laser Components). The photodiode signal is recorded by a high-speed digitizer card (ADQ14, Signal Processing Devices).

Utilizing the pulse shaper we generate a sequence of four collinear laser pulses and scan the first coherence time *τ* and the second coherence time *t* (“signal time”) from 0 to 84 fs with 6 fs steps in a rotating frame centered at 590 nm. In addition, the population time *T* is varied from 0 fs to 300 fs with 10 fs steps. Furthermore, we cycle the phases of the second, third, and fourth pulses in the sequence through $23\pi $, $43\pi $, and $2\pi $ with respect to the first pulse on a shot-to-shot basis. This means that the relative phase between the pulses is changed with every laser shot and that all possible phase combinations are completed after 27 laser shots. This 1 × 3 × 3 × 3 phase-cycling scheme allows us to disentangle the different fourth-order contributions to the fluorescence signal.^{18} During data analysis the time-domain maps of different phase-cycling steps are summed with weights according to the contribution of interest and Fourier-transformed to generate the 2D spectrum. We emphasize that phase cycling can also be performed for randomly undersampled data since it only consists of a sum and no change of basis (e.g., Fourier transformation) is involved in this analysis step. In total, the configuration space for the complete sampling including all coherence and population time steps and phase cycling consists of 15 × 31 × 15 × 27 = 188 325 different pulse shapes in our specific case. The acquisition of one complete map takes only about 3 min because of the shot-to-shot pulse-shaping capability. We average over 473 maps. The studied sample consists of 0.1 mM cresyl violet perchlorate (Radiant Dyes GmbH) dissolved in ethanol.

### C. The von Neumann basis as a sparse representation of 2D spectra

We use the von Neumann (vN) basis—also known as the Gabor basis in signal processing^{55,56}—as a sparsifying transform for the CS reconstruction algorithm. In the vN basis, a time or frequency-domain signal with a support of *N* data points is transformed into a joint $N\xd7N$ time-frequency space. As we have shown previously, this results in a compact representation of ultrashort laser pulses^{57–59} and a reduced parameter space for quantum control experiments.^{60,61}

Dealing with two time and two frequency axes, we are going to denote in the following the first coherence time *τ* and its Fourier-associated frequency *ω* with Greek characters and the second coherence time *t* and its associated frequency *w* with Latin characters. The vN basis functions $|gtm\mathit{w}n\u27e9$ are Gaussian pulses located at time *t*_{m} with a center frequency *w*_{n} which can be represented in the time domain as

wherein *α* is chosen as $\alpha =\sigma t2\sigma \mathit{w}$ with $\sigma t$ and $\sigma \mathit{w}$ being the full widths at half maximum in time and frequency, respectively.^{57} Any signal *r*(*t*) can be expanded in the vN basis by

with the expansion coefficients $atm\mathit{w}n$ and the inverse overlap matrix element $S(m,n),(k,l)\u22121$. The latter arises due to the non-orthogonality of the vN basis, leading to a basis-function overlap

for (*m*, *n*) ≠ (*k*, *l*).^{57} Thus, the coefficient $atm\mathit{w}n$ also has a contribution from the overlap of the basis function $|gtk\mathit{w}l\u27e9$ with the signal. This non-locality can be removed by using a basis which is bi-orthogonal to $|gtm\mathit{w}n\u27e9$, as shown by Shimshovitz and Tannor.^{56,62} In the bi-orthogonal basis, the inverse overlap element $S(m,n),(k,l)\u22121$ becomes part of a new basis function

Hence, the signal can be expressed in the new basis as

where the coefficients $ctm\mathit{w}n$ are now localized near the times *t*_{m} and frequencies *w*_{n}.

Since we deal with 2D spectroscopy based on two coherence time axes, we now proceed to generalize the vN representation for a 2D signal $r(\tau ,t)$. In that case, we require a vN basis set expressed in terms of four parameters $(\tau i,\omega j,tm,\mathit{w}n)$ that we will call the 2DvN representation. It can be obtained by numerically calculating the tensor product of two 1D vN basis sets,

A few examples of these 2D basis functions are presented in the supplementary material. The vN coefficients for a 2D time-domain signal are then calculated in analogy to the one-dimensional case (see Eq. (4)) by projecting $r(\tau ,t)$ onto the 2DvN basis,

Thus, the 2D time-domain signal is represented in a four-dimensional parameter space, which we flatten to two dimensions for purposes of graphic representation. This allows us to visualize the spectral as well as the temporal properties of the signal simultaneously.

### D. Genetic algorithm and signal reconstruction

As we aim to undersample the 2D time-domain data non-uniformly, we require methods of sparse recovery to reconstruct the 2D spectrum. Equation (5) can be recast in vector notation,

where $r\u21920$ denotes the time-domain signal and $c\u2192$ the vN coefficients. B is the matrix containing the time-domain representation of the bi-orthogonal basis functions. Choosing a non-uniform sampling scheme and measuring just a subset $r\u2192$ of the full time-domain data $r\u21920$ is equivalent to multiplication with a sampling matrix Ψ. Hence, the undersampled time-domain data can be expressed as a matrix multiplication

As we are interested in the complete set of coefficients $c\u2192$, this equation needs to be inverted. In general, a unique inversion is not possible since the inverse problem is under-determined. However, it was shown by Donoho^{36} and Candès *et al.*^{35} that the signal can be recovered with a high probability if it is sufficiently sparse in the coefficients. Several solvers have been developed^{63} which aim to find the coefficient vector $c\u2192$ with the smallest number of non-zero entries yielding the best representation of the observed time-domain signal $r\u2192$. In this paper, we use a Matlab implementation of TwIST^{52} to recover the 2DvN coefficients $c\u2192$. Technically, this is done by minimizing the objective function

where $\u2225\u22c5\u2225p$ is the *l*_{p}-norm of the vector. The regularization parameter *λ* controls the weight of the sparsity term with respect to an exact reconstruction of the measured data points.

We incorporate this reconstruction method into the framework of a genetic algorithm in order to find the sampling matrix Ψ with the best reconstruction results. Prior to optimization we need to determine the regularization parameter *λ*. For this we sample the data with a high density (35%) and apply a ternary search algorithm in order to find *λ* such that it yields the best reconstruction. Next, we start the optimization with a random Poisson-gap sampling matrix^{49} with the desired undersampling percentage, which is in our case 25% of the data points of the full sampling matrix. As depicted in the flowchart of the algorithm (Fig. 2) the matrix is mutated, i.e., a small fraction of the sampling points (“mutation rate”) is moved to a different, random position. The data are sampled using the new matrix and subsequently used as the input for the TwIST reconstruction algorithm. A fitness value is calculated comparing the reconstructed and the original coefficients. If the fitness did not increase with the mutation, the matrix is discarded. Otherwise, it is selected as the best matrix and the procedure is repeated. If there was no beneficial mutation in a certain number of generations, the mutation rate is lowered in order to explore the vicinity of the best matrix in more detail. The algorithm is stopped when the mutation rate has decreased to zero.

The total fitness *F*_{tot} is calculated as a weighted sum over the individual pixels,

where the individual pixel fitness is given by

comparing the amplitude of corresponding pixels *i* of the original 2DvN coefficients $\mu O(i)$ and the reconstructed 2DvN coefficients $\mu R(i)$. The weights

are calculated from the amplitudes of the original2DvN coefficients. This procedure is done for the real as well as the imaginary part of the rephasing and nonrephasing signals. All four fitness values are averaged to yield a single value for the reconstruction quality.

## III. RESULTS AND DISCUSSION

### A. 2D spectra in 2D von Neumann representation

The rephasing (photon-echo) signal is of particular interest in coherent 2D spectroscopy since it allows to observe molecular couplings and to extract the homogenous and inhomogeneous linewidths of the sample under study. For this reason, we use this signal contribution in order to demonstrate the properties of the 2DvN representation. In Fig. 3, we compare the real part of the measured rephasing 2D signal at a population time of *T* = 300 fs in three different representations, the time domain (Fig. 3(a)), the 2DvN domain (Fig. 3(b)), graphically flattened to two dimensions by plotting a 2D grid of 2D spectra, and the frequency domain (Fig. 3(c)).

Whereas the 1D vN transformation maps a 1D time- or frequency-domain signal into a joint 2D time-frequency representation $(\tau i,\omega j)$, the 2DvN transformation maps 2D time- or 2D frequency-domain signals into a joint 4D space with basis parameters $(\tau i,\omega j,tm,\mathit{w}n)$. We use a grid size of *N* = 25 for the time axes with the same temporal step size as in the experiment. This yields a $N\xd7N$ = 5 × 5 time-frequency grid for each of the two temporal dimensions, such that the same total number of 5 × 5 × 5 × 5 data points arises as in the Fourier representations. By choosing an odd number of frequency points we explicitly include zero frequency and thus non-oscillating contributions. Since the experimental data consist of 15 time steps, we zero-pad it accordingly prior to transforming it. We visualize the 2DvN space $(\tau i,\omega j,tm,\mathit{w}n)$ by plotting a 2D spectrum $(\omega ,\mathit{w})$ for each time pixel in Fig. 3(b) with the axis labels denoting the center position $(\tau ,t)$ of the 2DvN basis functions in the time domain. The 2D spectrum at the first time pixel $(\tau ,t)$ = (14 fs, 14 fs) is shown enlarged in the upper right of the figure (blocking graphically the underlying eight 2D spectra that are close to zero anyway, see below). It has the same general structure as the 2D spectrum in the regular Fourier domain, as can be deduced by comparison with Fig. 3(c). In addition, the decay of the amplitudes along the 2DvN time grid reflects the evolution of the signal in the time domain. This illustrates that the temporal as well as the spectral structure of the signal is preserved in the 2DvN representation (the preservation of information in the vN representation was proven previously^{57}). Thus, the 2DvN representation essentially behaves as a 2D short-time Fourier transformation with a Gaussian window. It is apparent that the 2DvN representation contains the essential information about the signal but is concentrated on a much smaller number of coefficients as compared to the Fourier representation, i.e., there are many more coefficients with vanishing amplitude in the whole plane of Fig. 3(b) than in Fig. 3(c). Therefore, mapping a 2D spectrum into 2DvN space increases the sparsity of the signal, and this will enable an optimized performance of compressed sensing reconstruction as shown below.

In order to further elucidate how the 2DvN basis compares to Fourier space with respect to sparsity, we subsequently remove 2DvN/Fourier coefficients starting with the smallest absolute values and transform the signal back into the time domain for comparison with the original signal. In Fig. 4 we plot the error which is introduced by this procedure in the real and imaginary parts of the time-domain data as a function of the percentage of retained non-zero basis coefficients. Defining an error of one percent as the threshold for a reliable reconstruction, the 2DvN representation requires 4.5% (16.8%) of the basis functions for reconstructing the real (imaginary) part. In contrast, using the Fourier representation, the threshold is reached with 45.8% (57.4%) percent of the basis coefficients.

### B. Optimized sampling and signal recovery

We applied the acquisition concept introduced above to measured data of the laser dye cresyl violet. The 2D spectrum exhibits pronounced oscillations in the crosspeaks as a function of the *T* delay, consistent with the literature^{64–67} and originating from vibrational coherences.^{68} Anticipating a decay of the signal for larger waiting times *T* and thus a greater sensitivity to the influence of experimental noise, we have chosen the 2D spectrum with the smallest signal amplitude at *T* = 300 fs as a reference data set and use this to optimize the sampling matrix. However, other times *T* are also feasible. We tested our approach with *T* = 90 fs as a reference data set and obtained similar results as the ones discussed below. Prior to reconstruction, we zero-pad the sampled data to fill the chosen *N* = 25 grid for each temporal dimension. For the calculation of the reconstructed 2D spectra arising from optimal sampling and their comparison with the regularly sampled 2D spectra, we truncate the reconstructed time-domain signal to the size of the original data.

The evolution of the fitness value during the optimization of the sampling matrix as a function of the number of generations is shown in Fig. 5(a). Within 7000 generations, the fitness has increased from an initial value of 0.78 to a final value of 0.99. Accordingly, the sampling matrix has evolved, as shown in Fig. 5(b). The sampling is now mainly concentrated at the origin of the *τ*-*t* parameter space, since most of the information about the signal is present in this region. However, sampling in intermediate and more distant regions is also necessary to get the best reconstruction result. In order to test if sampling in distant regions is really needed, we started with an initial sampling matrix having the sampling points distributed uniformly along *t* and *τ* in 6 fs steps from 0 fs to 42 fs, which corresponds to the bottom left quarter of the original time grid and thus also to a 25% coverage of the full sampling matrix. The genetic algorithm evolved this starting matrix into a similar structure as the optimized matrix shown here, demonstrating the need to sample in the outer regions for a better reconstruction. Applying a specifically designed semi-uniform starting distribution similar to the optimized sampling matrix obtained here might lead to faster convergence of the genetic algorithm and could be used if optimization time is of critical importance.

In Fig. 5(c) we correlate the 2DvN coefficients of the original and the recovered data for both the rephasing and the nonrephasing contributions on a double-logarithmic scale. For this we take the absolute values of the real and imaginary parts and show them in a joint scatter plot. The largest coefficients (down to the ∼5% level) line up very well along the diagonal (black line, indicating perfect correlation) and therefore prove the excellent reconstruction quality using just 25% of the original sampling points. The deviations in the smaller coefficients only play a minor role, since they do not contain the main information about the signal structure, as shown in Section III A and displayed in Fig. 4. This is also demonstrated in Fig. 5(d), comparing the rephasing and nonrephasing 2D spectra which were obtained from transforming the (reconstructed) 2DvN representation via the time domain into frequency space.

In the next step, we apply the optimized sampling pattern to the data at other population times, which we would hypothetically measure after the optimization using the result of the genetic algorithm (in the present work, the full data set was measured to facilitate a comparison). As the genetic algorithm takes about 20 min on a standard workstation, this can be incorporated into the usual lab routine without much loss of measurement time or changes in the experimental conditions. In Fig. 6(a) we show the absorptive spectrum of cresyl violet at a population time of 300 fs, which is calculated from the experimentally measured rephasing and nonrephasing signal contributions.^{34} The spectrum is presented with lab-frame frequency axes by adding the rotating-frame center frequency to the measured frequencies. In order to extract and analyze the vibrational coherence over the population time *T*, we select four regions of interest (ROIs) in the 2D spectrum. They are indicated by colored, filled circles. Their evolution as a function of population time *T* is plotted in Fig. 6(b) (same colors) in comparison with the full data set (gray circles), normalized to the ROI value of the reference data (*T* = 300 fs). Especially for the crosspeak (ROI 4) we can resolve the expected oscillation with a frequency of 580 cm^{−1}. The compressed sensing reconstruction closely follows the values from the full set, and the coherent oscillations in ROI 4 can be accurately resolved. A decreased modulation depth is attributed to deviations in the reconstructed 2DvN coefficients caused by the influence of experimental noise which poses a challenge to CS reconstruction algorithms. In the 2DvN representation, only a few coefficients carry the full information on the signal. Hence, deviations in the reconstructed values can alter the 2D spectrum in frequency space. Thus, by using only 25% of the original sampling points, we can recover the almost complete information about both the peak shape (Fig. 5(d)) and its temporal evolution (Fig. 6(b)).

## IV. CONCLUSION

We have suggested and implemented a new scheme for employing compressed sensing reconstruction in coherent two-dimensional (2D) spectroscopy. The basic concept is to optimize a reduced sampling matrix for best data recovery at one particular population time using a genetic algorithm and then apply that matrix for data acquisition at all other population times. The lower limit to the sampling density is given by the sparsity of the signal. Due to inhomogeneous broadening, liquid-phase 2D spectra measured in the rotating frame are not considered as sparse. In this case, the von Neumann representation of 2D spectra is useful to reduce the number of non-zero coefficients and thus to increase sparsity, which is beneficial for sparse recovery algorithms as shown here. Further applications of the 2DvN representation can be envisioned. For example, designing a pulse train capable of directly measuring the 2DvN coefficients would significantly speed up acquisition times, since only about 5% of the largest coefficients are sufficient to fully describe the signal. The 2DvN representation introduced here can also be useful for the analysis of conventionally recorded 2D spectra by condensing the data to few relevant coefficients while retaining the full information content.

We demonstrated that the approach of the present work can decrease measurement times in action-based 2D spectroscopy by a factor of four, especially when shot-to-shot pulse shaping is available and thus the number of collected data points can be reduced while keeping the full duty cycle of data acquisition. We have demonstrated the fidelity with experimental 2D fluorescence data obtained on cresyl violet. Optimal sparse sampling in connection with 2D spectroscopy will be especially interesting, for example, for the detailed study of photophysical processes and photochemical reactions, for which a high time resolution and a large range of population times *T* are required. We anticipate even larger savings in acquisition time for molecular systems that exhibit long-lived coherences such as in 2D spectra of atomic vapor^{16,17,44} or molecular beams,^{19} and thus feature sparse signals in the frequency and likewise in the von Neumann domain.

## SUPPLEMENTARY MATERIAL

See supplementary material for time-domain representations of the 2DvN basis functions centered at time pixel (*τ*,*t*) = (14 fs, 14 fs).

## ACKNOWLEDGMENTS

This work was funded by the European Research Council (ERC) within the Consolidator Grant “MULTISCOPE.” We further thank D. Tannor for fruitful discussions on the implementation of the 2DvN basis and S. Draeger for performing the measurements on the laser dye.