Partially coherent ultrafast spectrography

Modern ultrafast metrology relies on the postulate that the pulse to be measured is fully coherent, that is, that it can be completely described by its spectrum and spectral phase. However, synthesizing fully coherent pulses is not always possible in practice, especially in the domain of emerging ultrashort X-ray sources where temporal metrology is strongly needed. Here we demonstrate how frequency-resolved optical gating (FROG), the first and one of the most widespread techniques for pulse characterization, can be adapted to measure partially coherent pulses even down to the attosecond timescale. No modification of experimental apparatuses is required; only the processing of the measurement changes. To do so, we take our inspiration from other branches of physics where partial coherence is routinely dealt with, such as quantum optics and coherent diffractive imaging. This will have important and immediate applications, such as enabling the measurement of X-ray free-electron laser pulses despite timing jitter.


MSGPA Fourier space constraint
t' (fs) Supplementary Note 1: Mixture retrieval algorithm General principle. The Principal Components Generalised Projections Algorithm (PCGPA) is a widely used spectrogram inversion algorithm [1]. It iteratively synthesises and then reverses a complex spectrogram s (ω, τ ), each time fitting it to a given model in the least squares sense. The algorithm stops after a given number of iterations or once the error between the original and retrieved spectrograms becomes sufficiently low. It thus extracts the pulse P (t) and the gate G (t) from a measured spectrogram S meas (ω, τ ), with and where F t→ω (resp. F −1 ω→t ) denotes the Fourier Transform (resp. the inverse Fourier Transform). In order to choose the two constraints applied to the synthesised trace in the real and the Fourier domain, one assumes that the measurement is properly described by eq. (1). Consequently, the constraints are: -in the spectral domain, that the intensity |s (ω, τ )| 2 of the synthesised trace matches the measured intensity S meas (ω, τ ), can be summarised by a product of two functions Therefore, in the Fourier domain, the solution simply consists in replacing the amplitude of s (ω, τ ) with the measured one. In the real domain, one rotates the rows of the matrix O (t, t − τ ) to obtain a matrix O (t, t ) having an outer product form. If O (t, t ) is actually an outer product, then its rank equals one and O has only one non-zero singular value. In the general case, O (t, t ) is not an outer product, the solution is therefore to calculate its singular value decomposition and to approximate O as the outer product of the left and right singular vectors associated to the largest singular value, the principal components of O [1]. This is known to be the best outer product fit in the least squares sense. Alternative and less-time consuming techniques, such as the power method [2] or the minimisation of the least-squares error [3,4], also exist to estimate the principal components.
In the presence of a pulse mixture, the spectrogram corresponds to the incoherent superposition of several pure spectrograms: The PCGPA is adapted to partial coherence by parallelising the reconstruction of several pure spectrograms s k , see Supplementary Fig. 1. However, to ensure that these elementary spectrograms are related to each other, the constraints applied to each s k must depend on all the other spectrograms. The constraints in the real and Fourier domain are updated following those proposed in Ref. [5] to adapt the difference map algorithm used in ptychography to partial coherence. The new constraints are now: -in the spectral domain, that the intensity k |s k | 2 of the synthesised partially coherent trace equals S meas , In this new algorithm, called the Generalised Projections Algorithm for Mixed States (MSGPA), each spectrogram s k is weighted by an amplitude function S meas / k |s k | 2 in the spectral domain. As for the time domain constraint, we first rotate the rows of each matrix O k (t, t − τ ) to obtain a series of outer product forms. We then concatenate the matrices O k (t, t ), so that the principal components of the resulting matrix correspond to G (t ) and to a vector where the pulses P k (t) are stitched. This ensures that the new guess for the gate is synthesised from the information contained in every O k . When the MSGPA stops, one obtains a series of pulses P k (t) along with the gate. Each pulse P k (t) being meaningless if treated independently from the rest of the mixture, the retrieved pulses are then converted into statistical quantities, such as the two-time correlation function C (t, t ) or the Wigner distribution W P (t, ω): Initial guesses. For the specific case of FROG-CRAB, the gate is meant to be a pure phase function. The initial guess chosen for the gate is a field having a constant amplitude equal to one and a zero phase. As for the pulse mixture, pulses having various hermite-gauss time profiles are used in order to explore a large variety of pulses.
Convergence. One is able to retrieve partial coherence by giving the MSGPA new degrees of freedom, namely the construction of more than one pure spectrogram. Therefore, the algorithm is less constrained than the PCGPA, and the convergence is observed to require more iterations, typically a few thousands instead of a few hundreds. Moreover, the MSGPA performing the construction of several spectrograms in parallel, each iteration is more time-consuming.
Potential improvements. The various strategies described in Ref. [6] to speed up FROG algorithms, such as the "short-cut generalised projections" or the "over-compensation" technique, should be applicable to mixture retrieval algorithms. Moreover, the retrieval of states mixtures by the construction of several spectrograms in parallel is somewhat independent of the way the spectrogram is generated. Therefore, even though we only treated the case of the PCGP algorithm here, other FROG algorithms can be straightforwardly adapted to state mixtures. As an example, one could imagine using algorithms developed to deal with non-instantaneously nonlinearities [7] to recover mixtures of attosecond pulses with durations approaching the response time of atoms to photoionisation.

Supplementary Note 2: Decoherence deconvolution
Some sources of decoherence have an impact on the correlation function that can be easily modeled. We here analyse such an impact for two common sources of decoherence: the presence of a temporal jitter between the IR and XUV pulses, and the spectral response of the time-of-flight spectrometer. We then present a compact algorithm to separate the source of decoherence from the XUV pulse. Finally, we give numerical examples to illustrate the extraction of the pulse and of the decoherence from a spectrogram.
Impact on the measurement. Consider the presence of time jitter between a pure XUV pulse and the gate. If the ergodicity condition is fulfilled, i.e. if the acquisition time is sufficiently long so that the temporal mean matches the statistical mean, then the measured spectrogram S (ω, τ ) reads S pure represents the jitter-free spectrogram, and J (δτ ) corresponds to the probability for the IR/XUV delay τ to deviate from the expected value by an amount δτ . Therefore, the impact of time jitter on the spectrogram simply takes the form of a convolution between the pure spectrogram and the jitter envelope J (δτ ). The jitter also has a specific impact on the Wigner function W P (t, ω) and the correlation function C (t, t ) of the XUV pulse mixture. One reformulates eq. (5) to introduce the Wigner distributions, the jitter-free Wigner function. Using eq. (4) linking W P (t, ω) to C (t, t ), one obtains the corresponding two-time correlation function In the presence of jitter, C (t, t ) thus corresponds to the convolution of the pure correlation function with the jitter envelope. Consequently, in the spectral domain, the correlation function conveniently takes the form of a product between C pure (ω, ω ) and J (ω), respectively the Fourier Transforms of C pure (t, t ) and J (t).
The impact of the spectrometer resolution R (ω) on the measurement can be treated equivalently. If the response of the spectrometer is invariant over the spectral range of interest, its effect on the spectrogram can be modeled as a simple convolution of S pure (ω, τ ) with R (ω), which leads to the following expressions for the correlation functions in the spectral and time domains, with R (t) the inverse Fourier Transform of R (ω).
Ideally, one would like to suppress the influence of jitter or of the spectrometer response in order to recover the XUV pulse hidden in the mixture. Before explaining how this can be done numerically, it is fundamental to understand why the information about the pure pulse is still accessible in the correlation function. Consider the examples in Supplementary Fig. 5. In the three depicted scenarios, the pulse ensemble is composed of identical waveforms suffering from different types of temporal jitter. The chosen pulse has a gaussian spectrum (8 eV RMS) and a spectral phase exhibiting second and third order terms, corresponding to a 270 as pulse FWHM. The associated wave packet is centered at 100 eV.
In the jitter-free scenario shown in Supplementary Fig. 5 a), the two-time (resp. two-frequency) correlation function simply reads as P (t) P * (t ) (resp. P (ω) P * (ω )) and exhibits the characteristic symmetries of rank 1 matrices.
In the second case, the arrival time of the pulses in the ensemble can take two possible values, so that the jitter envelope takes the form two Dirac functions, see Supplementary Fig. 5 b). Although unrealistic, this scenario allows one to decompose the impact of jitter on the correlation function. In the time domain, C (t, t ) corresponds to the superposition of identical rank 1 patterns C pure (t, t ) = P (t) P * (t ) translated along the diagonal axis t+t 2 with a probability set by J t+t 2 , hence the convolution C pure (t, t ) ⊗ J t+t 2 . In the spectral domain, the convolution is changed into a product, and C (ω, ω ) is modulated by the oscillatory function J (ω − ω ). Yet, the jitter-free correlation function can easily be guessed, as highlighted by the gaussian function visible in the background.
In the last scenario depicted in Supplementary Fig. 5 c), a jitter envelope with a complex shape (1.4 fs FWHM) is accounted for. Despite the very low purity (ν = 0.09), one can still separate by eye the underlying pure correlation function (grey shaded curves) from the jitter (green shaded curve) in the spectral domain, since C (ω, ω ) must follow the form P (ω) · P * (ω ) · J (ω − ω ). Therefore, the general problem of separating jitter from the XUV pulse in a mixture in fact consists in fitting the two-frequency correlation function to this product form.
Algorithm for the blind deconvolution of the correlation function. One rewrites eq. (8) and affects the indexes i and j respectively to ω and ω to make apparent the number of measurement points: If time jitter is the only source of decoherence, then the correlation function contains a large amount of redundant information. It therefore becomes possible to extract two vectors, the pure XUV pulse P i and the jitter envelope J i , from the matrix C i,j . In fact, we found that it is numerically easier to consider P (ω i ) and P ω j as two distinct functions P (1) i and P (2) j . Therefore, one must split the redundant product between three functions varying along the vertical (i), horizontal (j) and diagonal axes (i − j). This situation is essentially the same as the one faced in the PCGP and MSGP algorithms when applying the product constraint to generate two estimates for P (t) and G (t), as shown in Supplementary Fig. 1. One can thus adopt a similar solution such as the least-squares minimisation. One wants to identify the three functions P One differentiates with respect to these three functions and sets these derivatives to zero. This leads to a set of three coupled equations: One proceeds iteratively to solve this system. Starting from random guesses for P i . The latter estimate is then inserted into eq. (15) along with the initial guess chosen for J j , leading to an estimate for P (2) j . The field associated to J j is finally updated by applying eq. (16). By repeating this process several times, one progressively refines the estimates of the three fields. This algorithm has been tested with experimental data in the context of interferometric measurements, and has proven to be robust against noise and to converge in very few iterations [8].
We now process the two-frequency correlation function C (ω, ω ) shown in Supplementary Fig. 5 c) with this algorithm. As shown in Supplementary Fig. 6 a), the algorithm converges in less than 10 iterations of eqs. (14,15,16), which takes a few seconds in practice. Increasing further the number of iterations does not decrease , the convergence being limited by the noise background originally added to the correlation function.
However, the three fields retrieved with this algorithm are known to suffer from one ambiguity [8]: Each function is thus properly retrieved up to linear phase term e iγ ω and up to an exponential term e γ ω affecting the field modulus. The linear phase term is unimportant since it corresponds to an absolute shift in the time domain, which is a meaningless quantity. However, the term e γ ω needs to be removed. To do so, we take advantage of the fact that P (1) (ω) and P (2) (ω) must correspond to the same function P (ω). In fact, if P (ω) is obtained by estimating the geometric mean of the two pulses, P (ω) = P (1) (ω) · P (2) (ω), then the terms e γ ω naturally cancel out, providing an unambiguous estimate for P (ω). Equivalently, calculating the ratio P (1) (ω)/P (2) (ω) = e γ ω allows one to isolate the exponential term which can then be used to obtain a corrected estimate for J (ω). As a last step, a simple Fourier Transform of P (ω) and J (ω) allows one to obtain the time profile of the XUV pulse P (t) as well as the jitter envelope J (t). The fields retrieved with this procedure along with the corresponding jitter-free correlation function P (t) P * (t ) are shown in Supplementary Fig. 6 (b-d), and agree very well with the original fields.
A similar procedure can be applied to separate the influence of the spectrometer resolution from the XUV pulse. In this scenario, the correlation function adopts a product form P (1) (t i ) · P (2) * t j · R t i − t j in the time domain, see eq. (11). One must therefore replace C ω i , ω j with C t i , t j in eqs. (14,15,16) in order to retrieve both the XUV field P (t i ) and R (t i ), the Fourier Transform of R (ω).
Numerical examples. We now illustrate the full "decoherence removal" procedure starting the processing from the spectrogram, see the three examples depicted in Supplementary Fig. 7. A pure spectrogram is generated and then convolved with a given filter along the delay or spectral axis to simulate the presence of time jitter or the response of the spectrometer (1). The spectrogram is then processed with the MSGPA to retrieve the correlation function C (t, t ) associated to the mixture of XUV pulses. Finally, the correlation function is processed with the help of the blind deconvolution algorithm (2) in order to retrieve the underlying XUV pulse (3).
The XUV pulse is the same 270 as pulse as in Supplementary Figs. 5 and 6. The laser pulse is a 800 nm Fourierlimited gaussian pulse (5 fs intensity FWHM) with an intensity of 8 TW cm −2 . In each case, noise is added to the final spectrogram with a level equal to 1% of the peak intensity of the photoelectron spectrum in the absence of laser field. The spectrograms, having a size of 512 × 512 points, are processed with the MSGPA using 20 modes, i.e. 20 pulses P k (t), to describe the mixture. The algorithm is stopped after 10 4 iterations. The deconvolution is then performed after 5 iterations of eqs. (14,15,16).
In the first scenario, no source of decoherence is added, see Supplementary Fig. 7 a). The pure spectrogram contains fine fringes due to the third-order term in the spectral phase of the attosecond pulse. Consequently, the retrieved correlation function P (t) P * (t ) corresponds to a matrix of rank 1, and the jitter deconvolution algorithm has no effect on it. The amplitude and phase of the pulse are well retrieved, along with the jitter envelope which corresponds to a Dirac distribution.
In the second case depicted in Supplementary Fig. 7 b), a temporal jitter with a gaussian envelope of 355 as FWHM is added to the pure spectrogram in a). As a consequence, the fringe structure is severely blurred and the retrieved correlation function is smoothed by the jitter. The pulse and jitter envelope recovered by the deconvolution algorithm nicely match the original ones.
In order to account for the resolution of the time-of-flight electron spectrometer, the pure spectrogram is convolved along the spectral axis with a gaussian function (3 eV FWHM), resulting in the disappearance of the fringes in the spectrogram, see Supplementary Fig. 7 c). In such a scenario, C (t, t ) reads P (t) P * (t ) R (t − t ), see eq. (11). This product form is equivalent to the one described in the lower panel in Supplementary Fig. 5 c) in the case of jitter. Consequently, the retrieved correlation function exhibits the same structure as C pure near the diagonal axis, as highlighted by the characteristic bounces. This structure rapidly disappears as one moves away from the diagonal due to the product with the gaussian function R (t − t ). The deconvolution algorithm is then applied. Despite a small artifact in the wings of C pure where the relevant signal can no longer be distinguished from the noise, the algorithm efficiently separates the influence of the pulse and of the response of the spectrometer, as shown in Supplementary Fig. 7 c) (3).

Supplementary Methods
Two-color photoionisation. The two-color XUV-IR photoionisation can be efficiently described under the Strong-Field Approximation (SFA). In the frame of the first order perturbation theory with the single active electron approximation, the SFA gives the following expression in atomic units for the amplitude a p of the transition from the ground state to the state |p of momentum p in the continuum [9]: In this equation, τ stands for the delay between the IR and XUV fields, I p is the ionisation potential of the atomic species, p (t) = p + A (t) represents the instantaneous momentum, A (t) is the vector-potential of the laser field E IR (t) = −∂A/∂t, E XU V (t) corresponds to the XUV electric field, and d p the dipole matrix element of the transition from the ground state to the state |p . The accessible quantity in practice is the square modulus of a p , that is the energy spectrum of the photoelectrons. Introducing the kinetic energy of the photoelectrons W = p 2 /2, one rewrites eq. (17) to make the spectrogram equation appear: As in the main text, denotes the averaging of the spectrogram over an ensemble. The pulse P (t) corresponds to the electron wave packet d p E XUV (t), and the gate G (t) equals exp (iφ (p, t)), where φ (p, t) = − +∞ t dt p · A (t ) + A (t ) 2 /2. Note that we do not apply the "central momentum approximation" here: φ (p, t) is evaluated for each momentum p, and the final amplitude a p is obtained by integrating numerically the function P (t − τ )·e iφ(p,t) ·e i(W +Ip)t over t. This procedure is used to calculate all the photoelectron spectra presented in the article. Moreover, the dipole matrix element d p has been assumed to be constant (flat cross-section and constant or negligible atomic phase) across the XUV pulse bandwidth. This is equivalent to saying that the electron wave packet is an exact replica of the XUV pulse, up to an energy shift I p . If the dipole matrix element is known beforehand, this procedure could be improved by dividing the retrieved wave packet in the frequency domain d p E XU V (W ) with d p in order to recover the optical wave packet E XU V (ω + I p ) [10]. Finally, we consider that the electrons are collected along the direction of polarisation of the laser field, which enables one to replace p · A with pA in the gate expression.
Simulation of the RABBIT measurement in Fig. 2. The mixture used to simulate the RABBIT measurement was formed by the incoherent averaging of the attosecond pulse trains P k (t) in Supplementary Fig. 2 a), i.e. C (t, t ) = P k (t) P * k (t ) k . All the pulses in the mixture have a common spectrum composed of gaussian harmonics (0.6 eV FWHM, 1.55 eV fundamental frequency) contained in a 9 eV FWHM gaussian envelope centered on harmonic 31. In order to illustrate the ability of Mixed-FROG to retrieve highly mixed states, we used exaggerated variations of the group delay dispersion (GDD) and of the absolute delay of the pulse trains. The GDD (resp. the absolute group delay) of the pulse trains varies linearly with k from 10 5 to 2 · 10 5 as 2 /rad (resp. from −1 to +1 fs). This results in a mixture with a very low purity ν of 0.26.
The photoelectron spectra were then simulated using argon as detection gas and a 8 fs (intensity FWHM) 800 nm gaussian dressing pulse. The traces consisted of 256 × 256 traces. In the case where the laser intensity equals 0.02 TW cm −2 , the spectrum and group delay retrieved with RABBIT and with FROG-CRAB after 10 4 iterations of the PCGPA appear in Supplementary Fig. 2 b). In order to identify the optimal mode number for the mixed state reconstruction, we processed the same RABBIT trace using various numbers of modes, see Supplementary Fig. 3 a). Following the definition of the R-factor used in ptychography to monitor the convergence of the difference map algorithm [11], we define the error R between the original S O and retrieved S R spectrograms as When using a single mode, which corresponds to the FROG-CRAB reconstruction in Supplementary Fig. 2 b), the error rapidly stagnates, highlighting that partial coherence is not properly accounted for. Using two modes dramatically reduces the error, which confirms the presence of a pulse mixture. Convergence can be slightly improved with 5 modes. However, using 10, 20 or 30 modes makes hardly any difference. The same conclusion can be drawn from the evolution of the retrieved temporal marginals shown in Supplementary Fig. 3 b). |P (t)| 2 no longer changes when more than 10 modes are used. Using 10 modes thus constitutes the best tradeoff between computation time and quality of the reconstruction of the mixture. This also shows that increasing the number of modes does not enable one to properly recover the original mixture shown in Fig. 2 in the main text. This confirms that the complete information about the mixture is fundamentally not present in the spectrogram when the RABBIT trace is obtained with a low intensity dressing field. The missing information can only be obtained by imprinting a stronger modulation onto the electron wave packet, that is by increasing the laser intensity. Therefore, the traces obtained for various laser intensities in Fig. 2 d) in the main text were processed with the MSGPA, using 10 modes and 10 4 iterations, and starting from the same initial state to ensure a relevant comparison.
Simulation of the attosecond pulse in Fig. 3. The photoelectron spectra were simulated using a replica of the driving pulse as the dressing beam (intensity 5 TW cm −2 ) and argon as detection gas. The spectrograms were first calculated at each radius r and then averaged to account for the spatial integration S (W, τ ) = 2π S (W, τ, r) rdr.
The resulting 512 × 512 spectrogram was then processed using 10 modes and 10 4 iterations of the MSGPA. The computation time was about 2 hours on a standard computer. The convergence is depicted in the supplementary movie.
Simulation of the Free-Electron Laser pulse in Fig. 4. The two-color photoionisation process was simulated using the ionisation potential of argon. The IR pulse has a gaussian spectrum centered at 800 nm (19 fs Fourier-limit) and a third order term (5 · 10 4 fs 3 /rad 2 ) was added on the spectral phase, resulting in a duration of 20 fs FWHM. The pulse intensity in the gas jet equals 2 TW cm −2 . The final spectrogram consisted in a 1024 × 1024 trace. In order to simulate the absence of CEP stabilisation of the dressing laser, the trace was first averaged over 20 values of the carrier-envelope-phase of the IR pulse equally distributed between 0 and 2π, resulting in the trace shown in the left panel of Fig. 4 a) in the main text. Using eq. (5), the trace was then convolved with a 20 fs FWHM jitter gaussian envelope to account for the laser/XUV jitter, leading to the trace in the right panel of Fig. 4 a) in the main text. The reconstruction was performed over 30 modes and 5 · 10 4 iterations of the MSGPA. 25 iterations of eqs. (14,15,16) were then necessary to deconvolve the obtained correlation function and retrieve the FEL pulse and jitter envelope. We depict in Supplementary Fig. 4 the evolution of the error R defined in eq. (19) with respect to the iteration number in the case of the Mixed-FROG reconstruction (MSGPA with 30 modes) and of the FROG-CRAB reconstruction (monomode PCGPA) depicted in Fig. 4 (d-e) in the main text. The same observation as in the case of the RABBIT reconstruction in Supplementary Fig. 3 a) can be made. While R rapidly stagnates at a value of ∼ 0.15 in the single mode case, it continuously decreases down to ∼ 5·10 −3 after 5·10 4 iterations when 30 modes are used. This highlights that in the presence of jitter the spectrogram is much more accurately reproduced if partial coherence is accounted for.