Open Access
Add to BibSonomyAbstract
2D electronic spectroscopy is a widely exploited tool to study excited state dynamics. A high density of information is enclosed in 2D spectra. A crucial challenge is to objectively disentangle all the features of the third order optical signal. We propose a global analysis method based on the variable projection algorithm, which is able to reproduce simultaneously coherence and population dynamics of rephasing and non-rephasing contributions. Test measures at room temperature on a standard dye are used to validate the procedure and to discuss the advantages of the proposed methodology with respect to the currently employed analysis procedures.
© 2016 Optical Society of America
1. Introduction
Two dimensional electronic spectroscopy (2DES) is gaining higher recognition worldwide as optical spectroscopic technique. From the experimental point of view, setups and signal acquisition procedures are becoming solid and reliable. On the other hand, fully established and standardized data analysis methods are still missing and this leads to a difficult extraction of complete information. It is thus necessary to develop global and robust data analysis procedures able to capture the complete picture with an increased level of clarity and reliability.
In 2DES, the third order signal is displayed as frequency-frequency 2D maps evolving during the population time (t2), which is the time delay between the second and the third interaction with light [1]. The evolution of the signal as a function of t2 carries information about the dynamics of the excited states, including the possible presence of coherent mechanisms, particularly investigated in the latest years [2–6]. Within the response function formalism, the third order signal can be expressed as a sum of contributions represented graphically by double-sided Feynman diagrams [7]. These contributions can be classified in two groups depending on the evolution of the signal during t2 as shown in Fig. 1. The first group includes non-oscillating pathways, represented by Feynman diagrams where the system reaches a pure state after the first two interactions. The second group consists of oscillating contributions described by Feynman diagrams where, after the first two interactions, the system is in a coherent superposition of states [4].
Fig. 1 Examples of energy level schemes and Feynman diagrams representing (a) non-oscillating and (b) oscillating contributions to the signal. Right panels sketch the dynamics of the corresponding signals as a function of population time.
In the first case, the signal evolves in t2 following the relaxation dynamics of the excited states that can be quantified through the solutions of suitable kinetic differential equations. For example, in the simplest case of parallel relaxation processes, the solutions of the rate equations are real exponential functions [8]. In the second case, the signal oscillates during t2 with a frequency proportional to the energy gap of the states that generates the coherence. These oscillations dampen over time according to their dephasing rates, depending on the nature of states themselves, on the temperature, on the environment etc., and are well described by complex exponential functions.
Several methods have been proposed to analyze 2DES signals as function of t2. Global analysis of population dynamics alone is an established procedure, and real multi-exponential models are the most used to fit the data because they do not require any additional assumption [9, 11]. The coherent dynamics is left in the residuals of the non-oscillatory model, requiring a second step of analysis. Different methods have been proposed to this aim, such as Fourier transform (FT) to achieve FT-maps [12–14], linear prediction Z-transforms [15–17] and time-frequency decomposition techniques. The latter method was proposed for the analysis of ultrafast optical responses by Prior et al. using continuous wavelet transforms [18]. This method and other linear and bilinear approaches were recently analyzed and compared [19]. The methods mentioned above have intrinsic limitations that can hardly be overcome, the most important one is their non-global character.
The procedure proposed here goes beyond the main issues of the standard analyses and it allows extracting all information from a complete 2DES dataset, analyzing simultaneously non-oscillating and oscillating components, without any preliminary subtraction operation. This is particularly relevant when fast decays and strongly damped low-frequency beatings cannot be easily disentangled with FT methods. Using the most general approach, both kinds of contributions can be described with a complex exponential function: in the former case, the imaginary part of the complex function is zero and the real part carries all the information on the decay constants. In the latter case, real and imaginary part are nonzero and the function accounts for amplitude, dephasing rate, frequency and phase of oscillation.
The algorithm consists in a single robust fitting procedure based on complex exponential functions. Moreover, it is also possible to study simultaneously both rephasing and non-rephasing signals. This method is applied to experimental 2DES data collected at room temperature for Zinc(II) phthalocyanine (ZnPc) in tetrahydrofuran (THF).
2. Methods
2.1. Fitting model
Mathematically, 2DES datasets can be described as three dimensional complex valued arrays Xijk in which the third order optical signal is collected as a function of excitation frequency, emission frequency and population time. The signal is visualized as a series of maps with (i, j) frequency indexes at the k-th population time. Hereafter the population time t2 will be simply labeled t with the index k indicating different sampled values.
A simple and general model able to reproduce the oscillating and non-oscillating dynamics of a time-evolving signal is a sum of N independently decaying complex exponential functions. This model is applied to the study of the evolution of the 2DES signal during the population time. Such a superposition of independent decays is often called parallel model [8]. The decomposition of a signal in terms of damped complex exponential functions also recalls the Prony analysis method [20–22].
The n-th exponential component is expressed as . The complex parameter embeds the phase ϕn and amplitude |an| information, while the exponential parameter bn includes decay and frequency properties. In a more explicit form, each exponential component can be re-written as , where the decay constant is the negative inverse of the real part of the exponential parameter, τn = −1/ℜ {bn}, and the angular velocity of the oscillation is the imaginary part, ωn = ℑ {bn}.
In order to ease the formulation of the fitting problem, the dataset dimensionality is reduced. The three dimensional arrays are converted into two dimensional ones using a collective index h which includes all the possible excitation and emission frequency coordinates, Xijk → Ykh. Each column of the newly defined Y matrix represents the decay of the signal at a specific coordinate of the 2DES map. The total number of frequency coordinates is H and the total number of sampled times is K, so that for the running indexes applies that 1 ≤ h ≤ H and 1 ≤ k ≤ K. Moreover, the dataset structure Ykh resembles the one of pump-probe data, thus allowing the application of well established analysis and fitting tools [8, 23].
In this work a global fitting procedure based on the variable projection algorithm [24–26] is developed following an established method for time-resolved spectra analysis [27, 28].
Each complex exponential parameter bn of the n-th component is shared between all the frequency coordinates of the maps. All bn parameters are arranged in a row vector b. The complex amplitudes are conversely resolved in the two frequency dimensions and are arranged in a matrix with elements Anh. Recalling the transformation X → Y, each row of A represents a map of the amplitude of the n-th component. The multi-exponential model is written as
This model can be recast in a more compact form in matrix notation as M = EA, where the matrix E has elements , i.e. each column contains a complex exponential function with unitary amplitude. Note that E is function only of b. The total number of parameters P = N(1+H) of the model can be partitioned in two groups: a small set of N exponential factors and a larger group of N × H complex amplitudes. We can organize all P parameters of the model in an array z of P elements in order to express the model function as M(z) : ℂP → ℂK×H.2.2. Variable projection algorithm
The problem of finding the best z that fits the experimental data Y can be written in term of the following unconstrained minimization of the least squares of the residuals
For a typical 2DES dataset the number of parameters is very large and finding the optimal set which satisfies problem in Eq. (2) is nearly impossible. For example, for a model with N = 10 components to fit maps with a resolution of 256 × 256 pixels, the total number of parameters is P = 10(1 + 2562) = 655370. In order to tackle this major limitation, the mathematical structure of the function M(z) can be exploited. In particular, the partitioning of the parameters in linear Anh and nonlinear bn allows the minimization problem to be separable [28]. One can recognize that the subproblem of minimization
is easy to solve for fixed b. Given E of full rank, i.e. there is no linear relationship between the columns of E, the subproblem can be solved analytically as A = E+Y, where E+ is the Moore-Penrose pseudo inverse [29]. In this way the optimization procedure can operate only on the nonlinear bn parameters and all the linear Anh parameters are analytically computed at each iteration of the minimization algorithm. In other words, for a given set of exponential factors b, the best possible amplitude maps able to fit the experimental data are immediately determined and easily accessible.Given the solution of problem in Eq. (3), the separable optimization problem can be formulated in the reduced space of b alone as
where I is the identity matrix of size K × K. The objective function of this minimization problem is called the variable projection functional and I − E(b)E(b)+ is the projector on the orthogonal complement of the column space of E [28].The elimination of the set of linear amplitudes parameters from the minimization problem has several benefits. The computation of matrix A and its initial guess are not necessary, therefore only the complex array b has to be estimated. Moreover, the number of iterations is drastically reduced with respect to the minimization in the complete space of parameters. The minimization of the variable projection functional is efficient and, more importantly, it has higher probability of finding the global minimum solution instead of a local one [26].
2.3. Tailoring the method for 2DES data analysis
In the specific case of 2DES data analysis, the fitting procedure here reported allows for a convenient global fitting of rephasing and non-rephasing data simultaneously. This can be particularly important in the investigation of the nature (electronic or vibrational) of the coherences excited during the experiments, since electronic and vibrational coherences typically manifest completely different behaviors in rephasing and non-rephasing parts of the signal [30]. To this purpose, the data matrix is built appending blocks of columns relative to the two experiments
where YR and YN are the rephasing and non-rephasing data, respectively. The model described above is then applied to the final matrix Y. As summarized in the scheme reported in Fig. 2, at the end of the fitting procedure, each row of the matrix A contains the two amplitude maps associated with a complex exponential decay component, one for the rephasing and one for the non-rephasing data. Each pair of complex amplitude maps is associated with a damping time and a frequency.Fig. 2 Schematic illustration of the fitting method. Rephasing and non-rephasing data (XR and XN) are subsampled and reshaped into the matrix Y, to which the global fitting procedure is applied. Decay constants, frequencies and matrix A are then recovered. Rephasing and non-rephasing amplitude maps are obtained from matrix A for each complex exponential decay component. Two types of maps can be identified: decay-associated spectra (DAS) and coherence-associated spectra (CAS).
The minimization over two complete rephasing and non-rephasing datasets is computationally intense. In order to lighten the amount of computations involved in each step of the minimization procedure, the data are subsampled in the frequency dimensions. We implemented a simple subsampling algorithm based on the construction of an evenly spaced grid over the 2D maps. The degree of subsampling is controlled by the dimension of the grid step g. The features recorded in typical 2DES maps are usually broad and the information about the evolution of the complete signal is captured using a reasonably small number of points of the map [31]. For the experimental example reported below, we observed that, even retaining less than 5% of the frequency coordinates, the output of the fit is conserved. After the convergence of the minimization problem using the sub-sampled data, the amplitude maps with full resolution are recovered using the solution of the minimization problem Eq. (3) on the complete dataset. The entire procedure on a subsampled dataset takes tens of seconds to converge on a standard computer.
A set of generally reasonable constraints can be implemented in order to reduce the parameter space to be explored. The exponential components are divided into two sets, the decay-set and the coherence-set, each with different boundaries. The former set models the non-oscillating dynamics of the data and it has frequencies fixed at 0 cm−1 and no boundaries are applied to the time constants. The latter set models the oscillating dynamics. Absolute values of frequencies of these components are taken smaller than the Nyquist frequency associated with the sampling over the population time. Damping times are forced to be positive. A coherence associated with a specific process appears with both positive and negative frequencies in complex rephasing and non-rephasing signals [13]; to further reduce the dimensionality of the minimization problem, it is thus possible to consider pairs of oscillating components with the same damping time and with the same frequency but opposite sign. Amplitude maps can be classified according to the set to which they are associated. Maps associated with the decay-set are called decay-associated spectra (DAS), in analogy with the definition previously proposed in similar methods [32, 33], whereas maps associated to the coherence-set are called coherence-associated spectra (CAS), see Fig. 2. Although the definition of DAS and CAS proposed above has been based on the assumption of a set of boundaries, this distinction emerges naturally directly from the unconstrained minimization problem. However, the a priori definition of suitable boundaries prevents solutions with CAS associated to a negative dephasing time and CAS with frequency higher than the Nyquist limit.
To assess if the model is satisfactorily reproducing the experimental data is not an easy task when dealing with multidimensional datasets. Moreover, the use of too many components when working with multi-exponential fits is a well-known issue (see Supporting Info of Ref. [34]). In this context, the singular value decomposition (SVD) is an useful tool to investigate the data in a reduced space [8, 35, 36]. For example, an estimate of the necessary number of components to be included in the fitting model can be obtained as the number of principal values in the singular spectrum of the matrix Y. Moreover the SVD of the residual matrix is helpful to check if the model is capturing all the dynamical features of the data [8]. A second test to assess if the correct number of parameters has been chosen consists in the calculation of the correlation matrix of fitted parameters [37]. The inspection of this matrix allows establishing the possible presence of interdependence between couples of parameters, in which case the number of parameters must be reduced.
3. Experiments
2DES experiment was performed at room temperature on a commercial (Sigma Aldrich®) Zn (II) phthalocyanine dissolved in spectroscopic grade THF (optical density was 0.25 with a pathlength of 1 mm). This dye is used as standard to validate the procedure and to discuss the advantages of the proposed methodology with respect to the currently employed analysis procedures. The experiment was conducted using a 3 KHz Ti:Sapphire Coherent® Libra laser system that pumps a commercial NOPA (Light Conversion® TOPAS White) generating pulses centered at 680 nm with a time duration of about 10 fs. After a pulse shaping stage, the transform-limited pulses enter a fully non-collinear interferometer with BOXCARS geometry inspired by Nemeth et al. [38]. The heterodyne detected third order signal is collected using a double lock-in method proposed by Augulis et al. [39].
As shown in Fig. 3(a) the laser spectrum is tuned to cover the electronic transition to the first excited state of ZnPc. Time-resolved fluorescence measurements ascertained for this state a lifetime in the nanoseconds timescale, confirming that no relaxation of the excited state population should be recorded in the ultrafast regime investigated here [40, 41]. The vibrational properties of the molecule have been characterized by resonant Raman spectroscopy and the relative Raman spectrum is shown in Fig. 3(b). The laser bandwidth used in the 2D experiment can excite simultaneously all the vibrational states within 1000 cm−1 on both red and blue sides of the maximum of the absorption spectrum (14580 cm−1), as highlighted in Fig. 3(c).
Fig. 3 (a) Normalized absorption spectrum of ZnPc in THF (black line) and laser spectrum (yellow area). (b) Raman spectrum of ZnPc powders with 633 nm excitation wavelength. (c) Energy levels diagram for ZnPc.
4. Results
We describe the molecular system in the framework of the displaced harmonic oscillator model (DHO), with the ground and the first excited states having the same vibrational potential energy surface displaced along the mode coordinate. The system consists in a set of independent oscillators associated to the modes coupled to the electronic transition. Within this model the vibrational properties of ground state and excited state are identical, i.e. frequencies and damping times of vibrational coherences are indistinguishable. As already discussed, this assumption, physically meaningful for many systems, simplifies considerably the dimensions of the minimization problem. However, the fitting model can be easily generalized including more components distinguishing between ground state and excited state coherences.
The photophysical data obtained from the preliminary time-resolved and Raman characterization have been used for a first estimate of the number and values of the fitting parameters. The fitting algorithm takes tens of seconds to fit a dataset of 256 × 256 rephasing and non-rephasing maps subsampled with 8 points grid step on a regular laptop computer. The outcoming parameters are listed in table 1.
Table 1. Output parameters of the fitting procedure applied to rephasing and non-rephasing 2D data collected on ZnPc solutions. Confidence intervals obtained from standard errors of the fit are less than 1 cm−1 for frequencies and less than 60 fs for time constants. The estimation of the errors was performed using a procedure based on the analysis of the Jacobian of the residuals as reported in [28].
Rephasing and non-rephasing maps recorded at t = 600 fs are reported in Fig. 4 together with six examples of fitted traces extracted in representative points of both rephasing and non-rephasing datasets. The multi-exponential model function well reproduces the decay and the beating of the experimental data. It is also worthy to notice that the fitting method is able to clearly resolve and distinguish signal features having similar frequencies and damping times, as displayed by the beat between the components 4 and 5 shown in Fig. 4(c).
Fig. 4 Rephasing (upper) and non-rephasing (lower) maps at t = 600 fs. Six traces extracted at representative points in rephasing (red lines) and non-rephasing (green lines) maps are shown (panels a–f).
The first two components (n = 1, 2) are non-oscillating decaying components. The first one has a time constant estimated to be ≫ 2 ps, it is thus almost constant in the investigated time window. It is related to the lifetime of the first excited state. The real DAS for the rephasing and non-rephasing datasets are shown in Figs. 5(a) and 5(c). They effectively account for the contributions to the signal that do not evolve within the investigated time window.
Fig. 5 Real part of DAS of the two non-oscillating components for the rephasing (a,b) and non-rephasing (c,d) signals. (a) and (c) are related to the component n = 1 with a long decay time. (b) and (d) are associated to the component n = 2 with time constant 0.38 ps.
The second non-oscillating component has a time constant of 0.38 ps. The physical origin of this decay is unraveled by its DAS shown in Figs. 5(b) and 5(d). Rephasing real DAS in panel (b), points out that the signal is decaying on the diagonal (red area) and rising on two parallel regions above and below diagonal (blue areas). This DAS represents a broadening of the rephasing peak as a function of t, typically associated with the spectral diffusion phenomenon [9, 10]. Moreover, as expected for spectral diffusion, this contribution is negligible in the non-rephasing real DAS in panel (d) [7].
The five oscillating components identified by the fitting procedure (n = 3–7) have frequencies in agreement with the main vibrational modes detected in the Raman spectrum, except for component 3 that has a frequency of 30 cm−1 lying outside the investigated Raman spectral range. It is convenient to analyze CAS in terms of modulus and phase maps. The modulus of CAS shows where that particular beating frequency contributes more in the 2D spectra, in analogy with the information provided by conventional FT-maps [13]. The phase of CAS displays the phase of the beating component at t = 0 fs. These maps have been demonstrated to be very sensitive to various system and laser pulse parameters. Despite the difficulty in their interpretation, their analysis may be of critical importance in the investigation of the origin of long lived coherences in multichromophoric systems [30]. Fig. 6 summarizes the complete set of information that can be extracted from the amplitude maps Anh for each beating frequency.
Fig. 6 Complete set of information obtained with the fitting procedure for a single oscillating component. The modulus and the phase of the CAS are shown for positive and negative beating frequencies and for rephasing and non-rephasing signals. As an example, the results for the 702 cm−1 component (n = 5) are shown.
The modulus of CAS for the five oscillating components are reported in Fig. 7. CAS of n = 3 component, as expected for a low frequency mode, is mainly localized in the central area of the map close to the maximum of the signal, Figs. 7(a) and 7(f). CAS of the higher frequency components present similar features, as shown in Figs. 7(b)–7(e) and 7(g)–7(j). In particular, in the rephasing CAS the oscillating signals are arranged in a “chair-pattern” as expected for vibrational coherences in the DHO model [42], as shown by the gray dots in Figs. 7(a)–7(e). The features are mainly localized below the diagonal, in particular at excitation frequency ωeg + ωn, where ωeg is the resonant electronic transition frequency. Conversely, in non-rephasing CAS the amplitude is mainly distributed on the upper diagonal.
Fig. 7 Sum of modulus of CAS associated to positive and negative frequency for n = 3–7 components of rephasing (red) and non-rephasing (green) signal. Gray dots identify the coordinates where the oscillating signatures are expected to contribute, according to the DHO model. For an easier comparison, maps (a,f) are scaled by a factor 3.
CAS of n = 7 component in Figs. 7(e) and 7(j) show a slightly different behavior since they present more contributions than what predicted by the DHO model. In fact, the frequency of this component (938 cm−1) is close also to a Raman mode of THF (914 cm−1), which can generate spurious contributions in the amplitude maps. Oscillating signatures resulting from vibrational modes of the solvent are known to contribute in the whole maps because of their non-resonant character [43]. Although the minimization procedure uses a single component to fit simultaneously the system and the solvent contributions, the corresponding CAS is able to clearly identify the presence of both contributions.
5. Discussion
The global fitting method proposed here has several advantages if compared to the currently available methods for the analysis of 2D spectra. The first remarkable feature is the simultaneous access to both the non-oscillating and oscillating dynamics of 2DES datasets in a unique step of analysis, while most of the previously proposed methods need separated and sequential treatments of the decaying part and of the beating part of the signal. Beyond the obvious advantage of reducing time and number of operations, this provides a remarkably higher reliability in the identification of the ultrafast decay and of the early part of the oscillating components. In fact, a standard real multi-exponential fit can easily fail to reproduce a fast decay because, at early times, it could be affected by strongly damped modes and low frequency oscillations. The proposed method naturally overcomes this issue since it fits simultaneously the decaying components and the oscillations. This is a particularly important aspect in the debate about the possible electronic nature of the beating recorded in 2D spectra of multichromophoric systems. Indeed, electronic coherence is usually strongly damped especially at room temperature [44]. The reliable determination of the amplitudes and damping times of oscillations contributing at early times is indeed at the base of a correct interpretation of the physical origins of the recorded beating [45].
A second remarkable advantage is the global character. The methods allows retrieving at the same time the frequencies, damping times and amplitude maps for all the fitted components considering simultaneously real and imaginary parts (i.e. the full complex dataset) of rephasing and non-rephasing signals. Since these features arise from common processes giving rise to real and imaginary, rephasing and non-rephasing signals, the ability of considering the dataset in its completeness makes the final results more reliable and robust.
Moreover, since this method offers the possibility of selecting the number of components of the fitting function, the user can choose to quickly identify only the main components contributing to the overall signal, or analyzing in details also the weaker features. Indeed, as expected for a fitting method based on a least square minimization procedure, the components contributing the most to the overall decay will be identified first, independently on the values of the initial guess. The procedure is very robust and the main components will be always identified with a high degree of reliability. Additional weaker components possibly present are dropped in the residuals and one can choose if further analyze the data increasing the number of components.
A current limitation of the method is that it assumes that both coherence and population dynamics follow an exponential decay. This assumption is often fulfilled in simple systems, and in general it could satisfactorily reproduce any dynamics if enough exponential components are used. In systems affected by more complex dynamics the matrix E should be differently modeled in order to meaningfully capture the desired kinetics, for example associating specific species to the components of the model. Nevertheless, the procedure is quite flexible: any kind of kinetic model can be implemented suitably defining a correct form for the E matrix and then the minimization procedure can be applied as described.
Compared to Fourier transform based methods, the proposed procedure gives a clearer picture if the signal is corrupted by noise. Indeed, Fourier transform uniformly distributes time-domain noise throughout the frequency domain leading to a limitation of the accuracy in extracting peak frequencies, widths and magnitudes. An additional source of spurious signals in FT methods is caused by the finite duration of the experimental time-domain signals. Fourier transform of truncated time-domain signal generates undesirable frequency-domain wiggles (called “Gibbs oscillations”) that hamper the identification of possible weak signals close to more intense peaks [46]. Indeed, if a FT-map at a certain frequency has an intense amplitude, it generates ghosts features in the FT-maps relative to close frequencies, hindering the detection of subtle features. Moreover, the performances of the FT methods are strongly dependent on experimental conditions, such as the time window and the time step used in the experiment, and this has a strong influence on the resolution and on the spectral range of the Fourier transform [46]. The method proposed here overcomes all these issues. The finite duration of the signal is indeed not a limitation since the fitting model reproduces the data within the chosen time window without truncation artifacts. The consequence is that components with close frequencies can be distinguished easier than in FT methods.
A relevant novel aspect of this method, if compared to FT methodologies, is its ability to distinguish components based on their dynamic behavior. Indeed, if the components have different physical origins and thus are characterized by different dephasing times, the fitting method will recognize different components with similar ωn but different τn and will produce CAS, also expected to have different amplitude distribution. It would be thus possible, for example, to distinguish between electronic and vibrational coherences, typically dephasing at room temperature in tens and thousands of femtoseconds, respectively. In the same way the contribution of solvent modes could be distinguished from molecular vibrations contributing in the same spectral region [43].
The disentanglement of beating components with close frequencies is a crucial aspect to assess the electronic or vibrational nature of coherent oscillations and verify the possible interplay between vibrational and electronic degrees of freedom. Time-frequency transform formalism could in principle give access to time resolution and to the same information, but it presents severe limitations due to artifacts. These methods do not present a sufficient time-frequency resolution to perform a robust quantitative analysis for typical 2DES signals [19].
6. Conclusions
We propose a global fitting procedure capable of simultaneously retrieving coherence and population dynamics from a full (real and imaginary, rephasing and non-rephasing) 2DES dataset. The algorithm has proven to be robust and flexible to adapt to the peculiarities of the system under investigation. The advantages with respect to the currently available methodologies have been discussed. In particular, amplitudes and damping times of oscillations contributing at early times are determined with high reliability and coherent dynamics can be disentangled exploiting the different time evolution of its components.
Funding
European Research Council (ERC) (Starting Grant QUENTRHEL 278560); Seventh Framework Programme for Research and Technological Development (FP7) (FET EC project MULTI 317707 and FP7 EU STREP project PAPETS 323901).
References and links
1. D. M. Jonas, “Two-dimensional femtosecond spectroscopy,” Annu. Rev. Phys. Chem. 54, 425–463 (2003). [CrossRef]
2. T. Brixner, T. Mančal, I. V. Stiopkin, and G. R. Fleming, “Phase-stabilized two-dimensional electronic spectroscopy,” J. Chem. Phys. 121(9), 4221–4236 (2004). [CrossRef]
3. M. Cho, Two-Dimensional Optical Spectroscopy (CRC, 2009). [CrossRef]
4. L. Valkunas, D. Abramavicius, and T. Mančal, Molecular Excitation Dynamics and Relaxation: Quantum Theory and Spectroscopy(John Wiley & Sons, Inc., 2013). [CrossRef]
5. E. Collini, “Spectroscopic signatures of quantum-coherent energy transfer,” Chem. Soc. Rev. 42(12), 4932–4947 (2013). [CrossRef]
6. A. Chenu and G. D. Scholes, “Coherence in energy transfer and photosynthesis,” Annu. Rev. Phys. Chem. 66, 69–96 (2015). [CrossRef]
7. S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford University, 1995).
8. I. H. M. Van Stokkum, D. S. Larsen, and R. Van Grondelle, “Global and target analysis of time-resolved spectra,” Biochim. Biophys. Acta Bioenerg. 1657(2-3), 82–104 (2004). [CrossRef]
9. F. V. A. Camargo, H. L. Anderson, S. R. Meech, and I. A. Heisler, “Time-resolved twisting dynamics in a porphyrin dimer characterized by two-dimensional electronic spectroscopy,” J. Phys. Chem. B 119, 14660–14667 (2015). [CrossRef]
10. R. Singh, G. Moody, M. E. Siemens, H. Li, and S. T. Cundiff, “Quantifying spectral diffusion by the direct measurement of the correlation function for excitons in semiconductor quantum wells,” J. Opt. Soc. Am. B 33(7), C137–C143 (2016). [CrossRef]
11. J. Dostál, J. Pšenčík, and D. Zigmantas, “In situ mapping of the energy flow through the entire photosynthetic apparatus,” Nat. Chem. 8, 705–710 (2016). [CrossRef]
12. D. B. Turner, K. W. Stone, K. Gundogdu, and K. A. Nelson, “Three-dimensional electronic spectroscopy of excitons in GaAs quantum wells,” J. Chem. Phys. 131, 144510 (2009). [CrossRef]
13. H. Li, A. D. Bristow, M. E. Siemens, G. Moody, and S. T. Cundiff, “Unraveling quantum pathways using optical 3D Fourier-transform spectroscopy,” Nat. Commun. 4, 1390 (2013). [CrossRef]
14. J. O. Tollerud, S. T. Cundiff, and J. A. Davis, “Revealing and characterizing dark excitons through coherent multidimensional spectroscopy,” Phys. Rev. Lett. 117(9), 097401 (2016). [CrossRef]
15. J. Tang and J. R. Norris, “LPZ spectral analysis using linear prediction and the z-transform,” J. Chem. Phys. 84, 5210–5211 (1986). [CrossRef]
16. G. Panitchayangkoon, D. V. Voronine, D. Abramavicius, J. R. Caram, N. H. C. Lewis, S. Mukamel, and G. S. Engel, “Direct evidence of quantum transport in photosynthetic light-harvesting complexes,” PNAS 108(52), 20908–20912 (2011). [CrossRef]
17. J. R. Caram and G. S. Engel, “Extracting dynamics of excitonic coherences in congested spectra of photosynthetic light harvesting antenna complexes,” Faraday Disc. 153, 93–104 (2011). [CrossRef]
18. J. Prior, E. Castro, A. W. Chin, J. Almeida, S. F. Huelga, and M. B. Plenio, “Wavelet analysis of molecular dynamics: efficient extraction of time-frequency information in ultrafast optical processes,” J. Chem. Phys. 139, 224103 (2013). [CrossRef]
19. A. Volpato and E. Collini, “Time-frequency methods for coherent spectroscopy,” Opt. Express 23(15), 20040–20050 (2015). [CrossRef]
20. M. R. Osborne, “Some special nonlinear least squares problems,” SIAM J. Numer. Anal. 12(4), 571–592 (1975). [CrossRef]
21. D. Kundu, “A modified prony algorithm for sum of damped or undamped exponential signals,” Sankhya 56(B-3), 524–544 (1994).
22. M. R. Osborne and G. K. Smyth, “A modified prony algorithm for exponential function fitting,” SIAM J. Sci. Comput.16(1), 119–138 (1995). [CrossRef]
23. C. Ruckebusch, M. Sliwa, P. Pernot, P. A. de Juan, and R. Tauler, “Comprehensive data analysis of femtosecond transient absorption spectra: a review,” J. Photochem. Photobiol. C 13(1), 1–27 (2012). [CrossRef]
24. G. H. Golub and V. Pereyra, “The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate,” SIAM J. Numer. Anal. 10(2), 413–432 (1973). [CrossRef]
25. G. Golub and V. Pereyra, “Separable nonlinear least squares: the variable projection method and its applications,” Inverse Problems 19(2), R1–R26 (2003). [CrossRef]
26. D. P. O’Leary and B. W. Rust, “Variable projection for nonlinear least squares problems,” Comput. Optim. and Appl. 54(3), 579–593 (2013). [CrossRef]
27. K. M. Mullen, M. Vengris, and I. H. M. Van Stokkum, “Algorithms for separable nonlinear least squares with application to modelling time-resolved spectra,” J. Global Optim. 38(2), 201–213 (2007). [CrossRef]
28. K. M. Mullen and I. H. M. Van Stokkum, “The variable projection algorithm in time-resolved spectroscopy, microscopy and mass spectrometry applications,” Numer. Algorithms 51(3), 319–340 (2009). [CrossRef]
29. G. H. Golub and C. F. Van Loan, Matrix Computations (The Johns Hopkins University, 2012).
30. V. Butkus, D. Zigmantas, L. Valkunas, and D. Abramavicius, “Vibrational vs. electronic coherences in 2D spectrum of molecular systems,” Chem. Phys. Lett. 545, 40–43 (2012). [CrossRef]
31. J. Almeida, J. Prior, and M. B. Plenio, “Computation of two-dimensional spectra assisted by compressed sampling,” J. Phys. Chem. Lett. 3(18), 2692–2696 (2012). [CrossRef]
32. E. E. Ostroumov, R. M. Mulvaney, J. M. Anna, R. J. Cogdell, and G. D. Scholes, “Energy transfer pathways in light-harvesting complexes of purple bacteria as revealed by global kinetic analysis of two-dimensional transient spectra,” J. Phys. Chem. B 117(38), 11349–11362 (2013). [CrossRef]
33. E. E. Ostroumov, R. M. Mulvaney, R. J. Cogdell, and G. D. Scholes, “Broadband 2D electronic spectroscopy reveals a carotenoid dark state in purple bacteria,” Science 340(6128), 52–56 (2013). [CrossRef]
34. T. R. Senty, S. K. Cushing, C. Wang, C. Matranga, and A. D. Bristow, “Inverting transient absorption data to determine transfer rates in quantum dotâǍŞTiO2 heterostructures,” J. Phys. Chem. C 119(11), 6337–6343 (2015). [CrossRef]
35. R. W. Hendler and R. I. Shrager, “Deconvolutions based on singular value decomposition and the pseudoinverse: a guide for beginners,” J. Biochem. Bioph. Methods 28(1), 1–33 (1994). [CrossRef]
36. E. R. Henry and J. Hofrichter, “Singular value decomposition: application to analysis of experimental data,” Methods Enzymol. 210, 129–192 (1992). [CrossRef]
37. H. Motulsky and A. Christopoulos, Fitting Models to Biological Data Using Linear and Nonlinear Regression: A Pratical Guide to Curve Fitting (Oxford University, 2004).
38. A. Nemeth, J. Sperling, J. Hauer, H. F. Kauffmann, and F. Milota, “Compact phase-stable design for single- and double-quantum two-dimensional electronic spectroscopy,” Opt. Lett. 34(21), 3301–3303 (2009). [CrossRef]
39. R. Augulis and D. Zigmantas, “Two-dimensional electronic spectroscopy with double modulation lock-in detection: enhancement of sensitivity and noise resistance,” Opt. Express 19(14), 13126–13133 (2011). [CrossRef]
40. H. Ohtani, T. Kobayashi, T. Ohno, S. Kato, T. Tanno, and A. Yamada, “Nanosecond spectroscopy on the mechanism of the reduction of methylviologen 4431 sensitized by metallophthalocyanine,” J. Phys. Chem. 88, 4431–4435 (1984). [CrossRef]
41. J. Savolainen, D. van der Linden, N. Dijkhuizen, and J. L. Herek, “Characterizing the functional dynamics of zinc phthalocyanine from femtoseconds to nanoseconds,” J. Photochem. Photobiol. A 196(1), 99–105 (2008). [CrossRef]
42. D. B. Turner, R. Dinshaw, K.-K. Lee, M. S. Belsley, K. E. Wilk, P. M. G. Curmic, and G. D. Scholes, “Quantitative investigations of quantum coherence for a light-harvesting protein at conditions simulating photosynthesis,” Phys. Chem. Chem. Phys. 14, 4857–4874 (2012). [CrossRef]
43. K. A. Fransted, J. R. Caram, D. Hayes, and G. S. Engel, “Two-dimensional electronic spectroscopy of bacteriochlorophyll a in solution: elucidating the coherence dynamics of the Fenna-Matthews-Olson complex using its chromophore as a control,” J. Chem. Phys. 137, 125101 (2012). [CrossRef]
44. D. Hayes, G. B. Griffin, and G. S. Engel, “Engineering coherence among excited states in synthetic heterodimer systems,” Science 340(6139), 1431–1434 (2013). [CrossRef]
45. A. Halpin, P. J. M. Johnson, and R. J. D. Miller, “Comment on “Engineering coherence among excited states in synthetic heterodimer systems”,” Science 344(6188), 1099 (2014). [CrossRef]
46. A. G. Marshall and F. R. Verdun, Fourier Transforms in NMR, Optical and Mass Spectrometry(Elsevier, 1990).
47. G. S. Engel, T. R. Calhoun, E. L. Read, T. K. Ahn, T. Mančal, Y. C. Cheng, R. E. Blankenship, and G. R. Fleming, “Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems,” Nature 446(7137), 782–786 (2007). [CrossRef]
48. Y. L. Sheu, H. T. Wu, and L. Y. Hsu, “Exploring laser-driven quantum phenomena from a time-frequency analysis perspective: a comprehensive study,” Opt. Express 23(23), 30459–30482 (2015). [CrossRef]
