Research ArticleOPTICS

Deconvolution of optical multidimensional coherent spectra

See allHide authors and affiliations

Science Advances  01 Jun 2018:
Vol. 4, no. 6, eaar7697
DOI: 10.1126/sciadv.aar7697


Optical coherent multidimensional spectroscopy is a powerful technique for unraveling complex and congested spectra by spreading them across multiple dimensions, removing the effects of inhomogeneity, and revealing underlying correlations. As the technique matures, the focus is shifting from understanding the technique itself to using it to probe the underlying dynamics in the system being studied. However, these dynamics can be difficult to discern because they are convolved with the nonlinear optical response of the system. Inspired by methods used to deblur images, we present a method for deconvolving the underlying dynamics from the optical response. To demonstrate the method, we extract the many-particle diffusion Green’s functions for excitons in a semiconductor quantum well from two-dimensional coherent spectra.


Optical coherent multidimensional spectroscopy (CMDS) has become an established tool for investigating material properties (112). It has been applied on a wide range of materials including photosynthetic complexes (1, 4), colloidal and epitaxial semiconductor quantum dots (6, 1317), atomic vapors (18), semiconductor quantum wells (7, 11, 12, 19, 20), metal surfaces (21), and two-dimensional (2D) materials (2224). Physical processes that are accessible through the additional spectral dimensions include energy transfer, coherent coupling, relaxation, and dipole-dipole interaction (4, 14, 18, 24, 25). Furthermore, the homogeneous and inhomogeneous linewidths can be determined separately (7, 23, 26, 27), giving access to microscopic dephasing and distributions.

Although CMDS has been successful, it remains a niche technique because of difficulty in understanding the rich spectra and obtaining insight into underlying material properties. Any insight is typically realized by comparison to theoretical results, which must incorporate both the elaborate theoretical tools needed to calculate the spectra (7, 2830) in addition to the particular materials and processes being studied. Furthermore, the information of interest is often obscured by spectral broadening such that interpretation of the spectra is sometimes described as “blobology.” This situation is further exacerbated by the lack of understanding about CMDS outside the spectroscopy community.

To address this challenge, we propose a new paradigm for analyzing coherent multidimensional spectra. Our approach is inspired by methods developed in imaging (3134) to deconvolve the effects of the imaging instrument from the acquired image. In a similar fashion, the effects due to the spectroscopic method can be deconvolved from the acquired spectra to reveal underlying material properties and dynamics. The deconvolution requires a theoretical description of CMDS that accounts for the presence/absence and strength of coherences, incoherent processes, the nature of the eigenstates, optical selection rules, and functions that parametrize the material model. Details inaccessible to the spectroscopy are averaged out. On the basis of this theoretical description, algorithms can be developed that implement deconvolution techniques to extract the functions parametrizing the model for the material. A few general descriptions, and hence algorithms, should cover a wide range of CMDS methods and materials. In the future, when more algorithms are developed and the approach matures, it may be possible to produce standardized programs that can be used by experimentalists as part of the routine data processing. Extracting a material’s properties in this way is easier and more intuitive for a non-expert in CMDS to interpret and allows for direct comparison with theoretical models of the underlying physical phenomena.

The method that we present here as proof of principle does not assume a specific form for the Green’s function, which describes energy flow in the material, and thus is suitable for a continuum, noninvertible cases, and ill-posed, ambiguous cases including noise. A different approach has been proposed for the inversion of 2D spectra for a few coupled pigments (35). In this case, the population transfer matrix is uniquely determined, which assumes a specific form of the Green’s function formulated in terms of relaxation rates and line shape. This assumption restricts the applicability to a few discrete states. Our method is also extendable to the discrete case but is particularly well suited for congested states and even continua or when dark states are important. However, the presented formulas and algorithms will require modifications and extensions for applications/materials beyond the presented example material type, where the assumptions used in the derivation do not hold.


To demonstrate this paradigm for analyzing CMDS, we select the example of spectral diffusion of the exciton distribution in a disordered semiconductor quantum well. For this example, a theoretical model and experimental data already exist (25). Since incoherent exciton relaxation dominates the evolution of the spectra, only the spectral line shapes and relaxation Green’s functions enter the simplified model (see the Supplementary Materials), making it a good candidate for demonstrating this concept. The extracted line-shape function contains the inhomogeneous distribution and energy-dependent homogeneous line shape. Prior efforts to analyze CMDS peak shapes only extracted constants characterizing inhomogeneous and homogeneous broadening, which assumes certain functional forms (26, 27).

An exciton is an electron-hole pair bound by the Coulomb attraction but free to move as a unit. Confining them in a quantum well increases their oscillator strength, and hence their optical nonlinearity, which provides a strong signal in a CMDS experiment (5, 12, 19, 20). Real quantum wells always have some degree of disorder, primarily due to fluctuations in the well thickness, which results in localization of some states and a mobility edge marking the gradual transition from localized to delocalized states (25, 36). The corresponding variation in energy produces inhomogeneous broadening in the optical spectrum of the excitonic resonance. Spectral diffusion results from spatial migration of excitons among the states, often mediated by acoustic phonons, including across the mobility edge.

The specific 2D coherent spectra presented here are produced by exciting a sample with a sequence of three cocircularly polarized pulses with wave vectors k1, k2, and k3, as illustrated in Fig. 1A. Their interaction gives rise to a signal in the direction kI = −k1 + k2 + k3, which corresponds to a photon echo if the pulse with wave vector k1 arrives first. 2D spectra can be generated by measuring the amplitude and phase of the signal as a function of both τ, the time between pulses k1 and k2, and τ’, the time over which the signal is emitted, and by performing a 2D Fourier transform. During the delay T between pulses k2 and k3, exciton relaxation processes can occur. Spectral diffusion from the initial absorption energy to the final emission energy of the excitons can be tracked through the evolution of the 2D spectra as function of T. Cocircularly polarized pulses are used to avoid the participation of bound biexcitons. Furthermore, coherences during T and exciton-exciton interactions do not significantly influence the spectra. Since the exciton resonance is spectrally narrow, we neglect effects of finite pulse bandwidth. Furthermore, the phonon system is assumed to be a bath with constant temperature.

Fig. 1 Pulse sequence of a 2D photon echo and visualization of optical transformation.

(A) The pulse sequence applied in the 2D photon echo spectroscopy. (B) Visualization of the transformation from an object O(x′, y′) to an image I(x, y) using the convolution with the PSF cf. Eq. 2.

Considering these assumptions and approximations, the 2D spectrum isS(Ω1,T,Ω2)=dω1dω2L*(Ω1ω1,ω1)G(ω1,ω2;T)L(Ω2ω2,ω2)(1)based on the sum-over-states approach (28), where frequencies Ω1 and Ω2 result from Fourier transformation with respect to τ and τ’ (see the Supplementary Materials for the derivation of Eq. 1 and eq. S33 for a detailed discussion of the validity range of Eq. 1).

The line-shape function, L(Δω, ω), depends on ω the exciton frequency and Δω, the detuning from ω. The line-shape function describes the two-time correlation function of the absorption and emission processes (37). The relaxation Green’s function, G1, ω2; T), is the probability that an excitation absorbed at ω1 is emitted at ω2 after time T. Extracting G1, ω2; T) is the main goal because it captures the exciton redistribution dynamics, which give rise to the spectral diffusion. Equation 1 shows that G1, ω2; T) is convolved in two dimensions with L(Δω, ω), so we must find a way to deconvolve them.

The problem of 2D deconvolution has been addressed in image processing. Specifically, the image of an object O(x′, y′) can be represented asI(x,y)=dxdyPSF(xx,yy;x,y)O(x,y)(2)where the point spread function (PSF) describes the effect of the optical apparatus (see Fig. 1B). The PSF is often extracted by using the image of a point source to enable reconstruction of the original O(x′, y′) from the image.

The structural similarity between Eqs. 1 and 2 suggests that methods used to reconstruct images might be applicable to 2D spectra. However, there is no equivalent to a point source for a 2D spectrum. Thus, we need a different strategy for determining L(Δω, ω) from a spectrum. For zero waiting time T ≈ 0, G1, ω2; T ≈ 0) = δ(ω1 − ω2)/D1), where δ(x) is the Dirac delta distribution and D(ω)is the density of states, and thus, the spectrum isSL(Ω1,0,Ω2)=dωL*(Ω1ω,ω)L(Ω2ω,ω)(3)where L(ω)=L(ω)/D(ω). The spectrum corresponding to Eq. 3 has a diagonal form, as shown in Fig. 2A, with the inhomogeneous width along the diagonal and the homogeneous width in the cross-diagonal direction. Since Eq. 3 only depends on L, we can extract L using an optimization algorithm, as described in the Supplementary Materials. If we can also extract D(ω), then L(ω) can be determined.

Fig. 2 Photon echo spectra and data extracted from T = 0 ps.

(A and B) Normalized experimental photon echo spectra (absolute value) for T = 0 ps and T = 20 ps at 5 K. (C) Reconstructed relaxation Green’s function for T = 0 ps at 5 K. (D) Absolute value of line-shape function L(Δω, ω) for 20 K. (E and F) Rescaled line-shape function L(Δω, ω)/|L(0, ω)| for 5 and 20 K, respectively, with the corresponding oscillator strength L(0, ω) and D(ω) given as inset. In (D), (E), and (F), the gray lines mark the area with low reconstruction error.

To extract D(ω), we need to use a different spectroscopic measurement that also depends on D(ω) in conjunction with L. One possibility is the linear absorptionαD(ω)=dωL(ωω,ω)D(ω)(4)from which we extract D(ω) using an optimization algorithm.


Examples of the input spectra and extracted functions used in the deconvolution are given in Fig. 2. Figure 2D shows a reconstructed L(Δω,ω). The absolute error between the calculated spectrum and the experimental data is minimized. As a result, the quality of the extracted line-shape function is only good in areas with large signals. In areas with lower signal strength, the reconstructed line-shape function may have random phase jumps and oscillations, resulting later in artifacts in the reconstructed Green’s function. L(Δω, ω) includes the line shape along Δω and the oscillator strength distribution multiplied by the density of states along ω. In Fig. 2 (E and F), the line shape L(Δω, ω)/|L(0, ω)|, the oscillator strength L(0, ω)/D(ω), and density of states D(ω) are plotted separately for temperatures 5 and 20 K. We focus on high oscillator strength areas with low reconstruction error ranging from 1543.8 to 1545.2 meV, marked by the gray lines. For 5 K, the linewidth increases with increasing energy since an increased number of scattering states are reachable for higher energies due to the increasing D(ω). For 20 K, the linewidth is broader than for 5 K, as expected, and the width stays almost constant inside the trusted area. This broadening results from the higher bath temperature that opens most scattering channels for lower exciton states, which do not contribute at 5 K. All states inside the distribution have a similar lifetime. The oscillator strength distributions are very similar for both temperatures, as expected.

After extracting the line-shape function L(Δω,ω) and the density of states D(ω), we are now ready to extract the relaxation Green’s function, G1, ω2; T), from Eq. 1 using an optimization algorithm. The parts of the Green’s function connected to bright states are successfully extracted, whereas those connected to dark states do not contribute to the signal. Thus, only part of the full Green’s function is successfully extracted, and the overall probability is not conserved since relaxation involving the dark states and exciton recombination occurs. In the following discussions, we focus on the area with sufficient oscillator strength for valid reconstruction (the area between 1543.8 and 1545.2 meV). For energies lower than 1543.8 meV, no excitons and therefore no oscillator strength exist in the quantum well, and thus, many spurious features appear in the Green’s functions in this spectral region. For energies higher than 1545.2 meV, a continuum of excitons with smoothly decreasing oscillator strength exists; therefore, distortion above 1543.8 meV is expected to be smaller but can still lead to false results. For T = 0 ps (see Fig. 2C), a perfect reconstruction would lead to a strict diagonal shape. The deviations from the expected diagonal shape will be used for T > 0 ps as an indicator for problems such as spurious features or noise that are not reliable.

In Fig. 3 (A and B), the Green’s function for the 1s exciton relaxation is shown for T = 10 and 20 ps at a temperature of 20 K. The reconstructed G1, ω2; T) is compared to the simulated result using the theory from the study of Singh et al. (25). The details of the reconstruction of the Green’s function are ambiguous, so it is possible that multiple Green’s functions reproduce the experimental spectrum equally well. The ambiguity represents the resolution limit and is influenced by the width of the line shapes, as well as the discretization and resolution of the experimental data. This ambiguity causes visible (oscillatory) noise in the reconstructed Green’s functions (examples of the ambiguous reconstruction can be found in the Supplementary Materials). Starting at T = 10 ps, off-diagonal contributions (around the horizontal line) in the Green’s function show exciton redistribution (spectral diffusion), almost covering a weak diagonal contribution. After 10 ps, the excitons are broadly distributed over more localized states with larger oscillator strength and lower energy, closer to the initially excited energy and nonequilibrium temperature. After longer delay times, the maxima of the distributions move toward higher energy until the maxima converge at the same final energy for different initial energies (visible as a horizontal feature parallel to the abscissa), which reflects the quasi-equilibrium distribution at the lattice temperature.

Fig. 3 Reconstructed and simulated relaxation Green’s function at 20 K.

(A and B) Reconstructed relaxation Green’s functions G1, ω2; T) for the initial energy ℏω1 and the final energy ℏω2 for T = 10 ps and T = 20 ps at 20 K. (C and D) Corresponding simulated relaxation Green’s functions. (The exciton–acoustic-phonon scattering in the second order Born-Markov approximation overestimates the relaxation times by a factor of 2 at 20 K (25); therefore, we use T = 5 and 10 ps from theory for the comparison.) In (A) and (B), the green contour line shows off-diagonal contribution at T = 0 ps, its presence indicates areas that may contain large artifacts and spurious features. Diagonal and horizontal lines provide visual guidance.

The simulated Green’s function shown in Fig. 3 (C and D) shows qualitatively the same behavior and prominent features, such as the disappearing diagonal and the horizontal off-diagonal contribution moving toward higher final energies. It includes scattering of exciton with acoustic phonons and radiative recombination in a disorder potential (25, 36). Since we know that the model with exciton–acoustic-phonon scattering in the second order Born-Markov approximation overestimates the relaxation times by a factor of 2 at 20 K (25), we use T = 5 and 10 ps from theory for the comparison. Overall, the Green’s function from simulation is much smoother than the reconstructed Green’s function.

At 5 K in Figs. 2C and 4 (A and B), the diagonal is more dominant than for 20 K. It is clearly visible in the Green’s function for T = 0 and T = 10 ps, since only a few excitons have recombined and scattered. It is much sharper in the Green’s function than in the spectra in Fig. 2 (A and B), highlighting the success of the deconvolution. The decay of the diagonal contribution is slower, as in the high temperature case with no redistribution of the off-diagonal contribution to higher energies for longer delay times. Instead, the distribution moves toward lower temperatures for longer delay times compared to the higher initial excitation, reflecting the lower bath temperature. Limitations in reconstructing the off-diagonal distribution with lower amplitude near the high-amplitude diagonal appear in the plot. The high-amplitude diagonal masks part of the low-amplitude contribution and generates echoes of the diagonal visible along the off-diagonal, visible above and below the diagonal. Again, the simulated relaxation Green’s function in Fig. 4 (C and D) shows qualitatively the same behavior as the reconstructed, for 5 K; also, the quantitative agreement is better compared to 20 K, since the second order Born-Markov approximation is more suitable for lower temperatures. However, we observe a stronger contribution above the diagonal (relaxation toward higher energies) in the extracted data than in the simulated data. We believe that this is caused partially by the larger reconstruction error from low oscillator strength in the area, which is seen as false signal at T = 0 ps as well. Higher (hot)–phonon temperature in the experiment caused by the excitation may be another reason, but we believe that it is mainly caused by reconstruction errors.

Fig. 4 Reconstructed and simulated relaxation Green’s function at 5 K.

(A and B) Reconstructed relaxation Green’s functions G1, ω2; T) for the initial energy ℏω1 and the final energy ℏω2 for T = 10 ps and T = 20 ps at 5 K. (C and D) Corresponding simulated relaxation Green’s functions. In (A) and (B), the green contour line shows off-diagonal contribution at T = 0 ps, and its presence indicates areas that may contain large artifacts and spurious features. Diagonal and horizontal lines provide visual guidance.

In conclusion, we have extracted Green’s functions from CMDS using deconvolution methods inspired by those developed for image processing. This approach allows a direct comparison between theoretical and simulated Green’s functions, which can be calculated by specialists in materials physics theory without requiring that they also become experts in the spectroscopic method. We illustrated this procedure for spectral diffusion inside the exciton manifold of a semiconductor quantum well, producing good agreement and enhanced insight into the processes occurring in the system. We were able to extract the energy-dependent homogeneous line shape and the oscillator strength.



The photon echo signal was generated by a sequence of three actively phase-stabilized cocircularly polarized excitation pulses with wave vectors k1, k2, and k3 (cf. Fig. 1A). The photon echo signal was collected in the direction kI = −k1 + k2 + k3. The signal was heterodyned with a reference pulse and detected through spectral interferometry to measure both amplitude and phase. The signal was recorded as delay τ was scanned, and delay T was kept constant. The signal was then Fourier-transformed with respect to τ. The sample was a four-period 10-nm-wide GaAs quantum well with 10-nm-wide Al0.3Ga0.7As barriers. The excitation was restricted to the heavy-hole exciton resonance with ~150-fs-long pulses. The experiment was carried out for different sample temperatures 5 K to 20 K using a sample-in-vapor helium flow cryostat [cf. study of Singh et al. (25) for more details].


The simulation used a sum over state treatment analogous to the study of Abramavicius et al. (28) for calculating the spectra. The exciton wave functions were obtained from a numerical solution of the 2D Schrödinger equations in relative and in center of mass coordinates. The calculation of the wave functions included Coulomb interaction and a random disorder potential caused by quantum well width fluctuations (36). The exciton wave function was used for calculating radiative and exciton-phonon scattering rates in the second-order Born-Markov approximation (36). The density matrix equations of motion were then solved numerically using the Portable, Extensible Toolkit for Scientific Computation (PETSC) library (38, 39) for obtaining the relaxation Green’s functions. In the end, the resulting quantities were averaged over several random realizations [cf. study of Singh et al. (25) for more details].

Extraction algorithms

For extracting a quantity f such as the line-shape or the Green’s function, first, a cost function C(f) was defined. The major contribution to the cost function C(f) is the error between the measured quantity (for example, a spectrum) calculated from f compared to the experimental data. Other contributions to the cost function C(f) ensured specific features of f, such as a specific functional form, smoothness, etc.

The cost function was then minimized using the TAO (Toolkit for Advanced Optimization) package from PETSC (3840) applying suitable constraints. More details about the used cost functions and the extraction procedure are provided in the Supplementary Materials.


Supplementary material for this article is available at

Supplementary Text: Reconstruction algorithms

Supplementary Text: Derivation of material model

fig. S1. Different reconstructed Green’s function G1, ω2; T) at 5 K and T = 10 ps.

fig. S2. Contributing pathways to the photon echo signal.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: This work was inspired, in part, by the suggestions of R. Merlin (U. Michigan). Funding: The work at University of Michigan and JILA was primarily supported by the Chemical Sciences, Geosciences, and Energy Biosciences Division, Office of Basic Energy Science, Office of Science, U.S. Department of Energy under award no. DE-FG02-02ER15346 and no. DE-SC0015782. The work at Technische Universität Berlin was supported by the Deutsche Forschungsgemeinschaft through SFB 951 B12 and GRK 1558 A4. Author contributions: S.T.C. conceived the experimental concept. R.S. and M.S. ran the experiments. M.R. designed and calculated the simulation and the extraction algorithms. S.T.C. and M.R. wrote the manuscript. All authors discussed the results and commented on the manuscript. The discussions of all authors lead to the idea of the extraction algorithms. Competing interests: S.T.C. is an inventor on a patent application related to this work filed by the University of Michigan (no. 20180073856, 15 September 2017). All the other authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article