Retrieval of attosecond pulse ensembles from streaking experiments using mixed state time-domain ptychography

The electric field of attosecond laser pulses can be retrieved from laser-dressed photoionisation measurements, where electron wavepackets that result from single-photon ionisation by the attosecond pulse in the presence of a dressing field are produced. In case of fluctuating dressing laser and/or attosecond pulses, e.g. due to pulse-to-pulse fluctuations of the carrier envelope phase of the infrared laser pulse, commonly applied retrieval algorithms result in the erroneous extraction of the pulse fields. We present a mixed state time-domain ptychography algorithm for the retrieval of pulse ensembles from attosecond streaking experiments.


Introduction
Many advances in the field of ultrafast and attosecond science [1] are related to progress in ultrashort pulse laser sources and their metrology [2]. Advances in metrology have led to techniques capable of characterising the electric field of singleand sub-cycle laser pulses [3]. Combining state of the art metrology with advanced laser pulse generation and shaping technology even enabled the synthesis of optical attosecond pulses [4,5]. There has also been a lot of progress towards full spatio-temporal pulse metrology. Reference-based and self-referencing technology has been developed to measure the full spatial field [6,7] and space-time couplings [8], which can be important in the optimisation of next generation laser sources [9].
In all of the above-mentioned methods the laser pulses are assumed to be reproducible from pulse to pulse. Already in Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1969 Fisher and Fleck pointed out that intensity autocorrelation is prone to erroneous results when laser pulses with shot-to-shot spectral phase fluctuations are considered [10]. In this case, the laser pulse duration is not correctly retrieved and in fact underestimated. The autocorrelation measures only the 'coherent artefact', i.e. it only reflects the Fourier-limited rather than the actual pulse duration. More recently, Ratner et al revived the discussion of the 'coherent artefact' and carefully investigated its influence in modern pulse characterisation methods [11,12]. Escoto et al presented updated retrieval algorithms to take into account an incoherent background in FROG measurements and showed that it is possible to overcome the influence of the coherent artefact and correctly retrieve the fluctuating spectral phase [13,14].
Many methods in attosecond metrology are based on laser-dressed photoionisation. One of the common techniques to characterise attosecond pulse trains is reconstruction of attosecond beating by interference of two-photon transitions (RABBITT) [15], whereas isolated attosecond pulses are commonly characterised by attosecond streaking [1]. The two techniques are merged under the common framework of frequency-resolved optical gating for complete reconstruction of attosecond bursts (FROG-CRAB) [16]. Although, the pulses to be measured are not necessarily reproducible from pulse to pulse, this fact is usually ignored in the literature. In this paper we present a retrieval method that can correctly reconstruct distributions of attosecond pulses and dressing laser fields even when pulse-to-pulse fluctuations are present in either pulse.
Our motivation for the development of this method was the characterisation of attosecond pulse trains in experiments using a partially carrier envelope phase (CEP)-stabilised few-cycle laser system. These experiments and the resulting pulse characterisation are described in our companion paper 'generation and characterization of few-pulse attosecond pulse trains at 100 kHz repetition rate' [17], which is also part of the current special issue. In this paper we focus on the retrieval algorithm itself and demonstrate its ability to perform pulse retrieval in the presence of pulse-to-pulse fluctuations by presenting simulations where the dressing field has CEP noise and/or the attosecond pulse has pulse-to-pulse fluctuations of its chirp.

Attosecond streaking with fluctuating pulses
In order to characterise attosecond laser pulses, laser-dressed photoemission spectra are usually recorded as a function of the delay between an NIR dressing pulse and the attosecond XUV pulse or pulse train. Let us consider the transition amplitude of photoelectrons from the ground state to their final continuum state in the framework of the strong-field approximation (SFA) [18,19]: In this expression atomic units are used. p(t) is the instantaneous momentum of the free electron and v is the canonical momentum (i.e. the momentum that is recorded on the detector), related to each other by p = v + A(t). W is the kinetic energy of the photoelectron W = v 2 /2. τ is the time delay between the NIR and XUV pulses, I p the ionisation potential of the target atom, and d p the dipole transition matrix element. The phase modulation imparted by the vector potential A(t), which is related to the dressing laser field by E NIR (t) = −∂ A(t)/∂t, onto the photoelectron wavepacket produced by the XUV pulse E XUV (t) is given by We define the photoelectron wavepacket produced by the XUV as the unknown pulse P(t) to be determined, and the phase modulation term caused by the NIR dressing field as gate G(t): with φ(t) given by equation (2). This allows writing the acquired signal as the Fourier transform of the product of the unknown pulse with a time-shifted gate: with F {} denoting forward Fourier transformation. The formulas (1)-(4) as defined above are applicable in case both the XUV and NIR fields are fully reproducible from pulse to pulse, however this is often not the case. In many experiments there exists a CEP jitter, that leads to shot-toshot fluctuations in the shape of both the NIR and XUV fields. This leads to a situation where the experimentally measured streaking trace (or RABBITT spectrogram) differs from the trace that results when using the averages of both the pulses and the gates. In other words, the experimentally measured streaking trace is not equal to the trace that would be obtained when the coherent sum of all pulses would be streaked by the coherent sum of all gates. Mathematically this is expressed by the inequality where k sums over the K different combinations of pulses and gates encountered in the experiment. The consequence of this inequality is as follows: if either the attosecond pulses or the dressing fields do not reproduce from pulse to pulse, then the retrieval will lead to the determination of an attosecond pulse and a dressing field that will be incorrect and that might underestimate their respective pulse durations. An example of this is shown in figure 1. Here, an RABBITT spectrogram is simulated for the case of an attosecond pulse train consisting of one dominant attosecond pulse and two satellites that each have an intensity that is about 40% of the intensity of the dominant attosecond pulse (grey curve in figure 1(b). The RABBITT spectrograms are calculated assuming a 7 fs, 0.5 TW cm −2 dressing field, with a CEP jitter of 0.25π rad (curves in figure 1(a)). Figure 1(c) shows an RAB-BITT spectrogram for one selected NIR pulse and figure 1(d) shows the incoherent sum of the spectrograms calculated for all NIR pulses with their different CEPs. Retrieval with a standard (single pulse) time-domain ptychography algorithm based on the sum-spectrogram (shown in figure 1(d)) leads to the attosecond pulse shown as the blue curve in figure 1(b), revealing clear deficiencies in the retrieval. Compared to the attosecond pulse train that was used to calculate the RAB-BITT spectrograms, the intensities of the satellites are severely understimated, suggesting that the attosecond pulse train is close to an isolated attosecond pulse. This is in fact precisely the situation that we encountered in the experiment that we report in our companion paper [17], where the initial pulse retrieval produced an attosecond pulse train that was substantially shorter than expected on the basis of the duration of the NIR driver laser used in the high harmonic generation.
The situation encountered in figure 1 is somewhat reminiscent of the coherent artefact in femtosecond metrology [10]. We conclude that conventional attosecond pulse retieval algorithms will fail to correctly retrieve the attosecond pulses and electric fields of the dressing laser in case of fluctuations in either pulse. (e) Residuals map, calculated as the difference between the measured and the retrieved traces. The residuals map reveals that the algorithm failed to retrieve some periodic structures that were part of the measured trace (d). This is also reflected in the high normalised rms error F = 2.7% (commonly called FROG-error). The pulse retrieval converges to a flatter spectral phase (not shown here), and thus underestimates the intensity of the satellites in the XUV pulse temporal structure (b) (blue line).

Retrieval algorithm
Recently, a principal components generalised projections algorithm (PCGPA) has been extended to mixed states and has been used to retrieve ensembles of XUV pulses from streaking traces [20]. Shortcomings of the PCGPA algorithm for attosecond pulse retrieval have been pointed out in the past [21]. For the present paper we opted to use a time-domain ptychography algorithm. The concept of ptychography was originally proposed by Hoppe for phase retrieval from diffraction patterns in x-ray crystallography [22]. Spatial domain ptychography is now widely used in lensless imaging [23]. With progress in the processing speed of modern computers, fast retrieval algorithms exist for spatial imaging, e.g. the extended ptychographic iterative engine (ePIE) presented by Maiden and Rodenburg [24].
The concept of ptychography is readily transferred to ultrafast optics when considering the unknown time-domain field P(t) as 'spatial domain object', and measured spectra S(ω, τ ) resulting from a non-linear interaction with a probe pulse G(t) as 'k-space diffraction patterns'. Time-domain ptychography was first demonstrated by Feurer [25][26][27] and subsequently also used to characterise near single-cycle laser pulses [28] and short UV laser pulses [29]. Lucchini et al have introduced time-domain ptychography to the XUV domain and have reconstructed attosecond XUV and short NIR pulses from streaking measurements [30].
Here we demonstrate an extension of the time-domain extended ptychographic iterative engine (td-ePIE) algorithm presented in [28] to ensembles of pulses. We refer to our algorithm as time-domain ensemble extended ptychographic iterative engine (td-e 2 PIE). The argumentation follows the algorithm of Thibault and Menzel, which was presented for the retrieval of microscopic objects under partially coherent illumination [31].
The starting point of the retrieval procedure is the definition of two ensembles that contain K XUV pulses P k (t) and L gates G l (t) that are vectors of size N t . Within the algorithm, the number of points N t used for representing the pulse and the gate in the time domain is equal to the number of points N ω that is used for their spectral representation. These two ensembles represent a finite sampling of the probability distributions of the pulses and gates present in the experiment. The purpose of the retrieval is to determine the ensembles that best reproduce the experimentally measured streaking or RABBITT trace, i.e.
In the above expression the equality l = k is a consequence of the fact that the unknown pulse and the gate are part of a timeordered sequence of corresponding pulses and gates. However, the retrieval method in principle also allows the analysis of experiments where this is not the case. The retrieval consists of an iterative procedure of updating the guesses of P and G that is continued until (7) is satisfied. The ensembles are initialised as follows: the pulse ensemble is built by constructing K vectors of complex numbers P k (t) with random amplitudes and phases: Similarly all L gates G l (t) are initialised with random phases according to In these expressions R [a, b] indicates a vector of size N t with random numbers in the interval a to b. We found initialisation with random numbers led to better results compared to using e.g. the result of a single pulse algorithm as a seed. Random initialisation also avoids any bias. The ensembles P initial and G initial are then fed into the algorithm. Within the algorithm the ensembles are iteratively adapted until they converge to a result that reproduces the measured trace. The operating principle of the algorithm is similar to that of an FROG algorithm [32] or the td-ePIE algorithm in [28]. However, whereas there typically exists only a single pulse and a single gate in single pulse algorithms, the current retrieval algorithm needs to update K pulses and L gates. In the standard FROG/td-ePIE algorithm the single pulse and the single gate are updated once per iteration by applying both a frequency domain and a time domain constraint. In the current algorithm this needs to be done for each pair of pulses and gates. Accordingly, in each iteration n we loop over all M time delays τ m , and for each τ m we loop over all combinations of the K XUV pulses P (n−1) k (t) and L gates G (n−1) l (t) in the respective ensembles. All L gate pulses G (n−1) l (t) in the ensemble from the previous iteration (or the initial guess) are shifted to the current delay position τ m by means of a linear phase ramp in the spectral domain For each pair of pulses and gates the signal spectrum S (n−1) k,l (ω, τ m ) at time-delay τ m is calculated using pulse P (n−1) k (t) and the shifted gate G (n−1) l (t − τ m ) resulting from the previous iteration n − 1. The first step is then to apply the Fourier domain or intensity constraint [32]. During this step the calculated signal spectrum S (n−1) k,l (ω, τ m ) is replaced by a signal spectrum S (n) k,l (ω, τ m ), in which the phase is kept, but the amplitude is replaced by the experimentally measured amplitude: with and Note that, just like in equation (4) the signal spectrum is evaluated as Fourier transform of the time-domain signal field, which in turn is the product of pulse and gate. But now this is evaluated for each combination of K pulses and L gates: One other key difference to the standard single pulse algorithm is the denominator of equation (11), which now contains an incoherent sum over all possible combinations of K pulses and L gates. Taking the inverse Fourier transform of S (n) k,l (ω, τ m ) yields the modified complex-valued signal field in the time domain: With this updated time-domain signal field, the pulses and gates can be updated. This is done using the update rules presented in [28], which are now modified to take into account the ensemble: The scaling parameters β P and β G are chosen as random numbers in the interval [0.01, 0.1], and [0.1, 0.5], respectively. These value ranges have been determined by assessing the convergence performance for simulated traces. As one can recognise from equations (15) and (16), the new estimate of the K pulses P (n) k (t) involves the evaluation of a sum over the L gates, whereas the new estimate for the L gates G (n) l (t) involves the averaging over the K pulses. Note that the sums in the denominators are incoherent sums, whereas the sums on the right-hand sides of (15) and (16) are coherent sums. The * denotes complex conjugation. Finally, all updated gates are shifted back to the time origin by a phase ramp in the spectral domain that reverts the time shift provided by (10). After all τ m have been processed the iteration n is completed. In our implementation the individual streaked spectra S(ω, τ m ) are processed in random order within each iteration. I.e. the indices m = {1, . . . , M} are randomly permuted, when picking the delays τ m in each iteration n. For ensembles of size N ω = N t = 512, K = 10, L = 10, and M = 100, convergence is usually achieved after N ≈ 100 iterations. On a contemporary i7 CPU one iteration takes about 1.5 s.

Results
We demonstrate the td-e 2 PIE algorithm for various scenarios shown in figures 2-5. RABBITT and streaking traces are simulated using the SFA model of equations (1) and (2). First, we consider the characterisation of an isolated attosecond pulse with a Gaussian spectrum, and with a Fourier limited FWHM pulse duration of 316 as. Because in case of fluctuations a single pulse algorithm fails, and retrieves a pulse close to the Fourier limit (compare figure 1), a chirp of 0.1 fs 2 is applied, which stretches the attosecond pulse to 943 as. The attosecond pulse is characterised by means of a streaking measurement with a dressing pulse (λ 0 = 800 nm, Δt = 5 fs, and an intensity of 2 TW cm −2 ) with a CEP jitter. A measured trace is constructed by sampling 120 CEP values from a Gaussian distribution with σ orig. = 0.235π rad (see black line in figure 2(h)). With these gate pulses the streaking traces S k,l (ω, τ m ) are calculated. An example trace is shown in figure 2(c). These individual traces are then incoherently summed to construct the 'measured' trace (figure 2(e)). Next, the td-e 2 PIE algorithm is used to retrieve the ensembles of gates and pulses from this 'measured' trace. Given the assumption of an isolated attosecond pulse that is reproducible from pulse to pulse, the ensemble of pulses P contains only a single member, whereas the choice L = 100 was made for the gate pulse ensemble. Note that the choice of L determines the ability of the algorithm to reproduce the assumed distribution of CEPs. On the one hand, if L is chosen too small the Gaussian distribution of CEP values is sampled too coarsely. On the other hand if L too large, convergence is compromised. After 150 iterations the streaking trace figure 2(f) is retrieved, which very closely resembles the 'measured' streaking trace shown in figure 2(e). The good retrieval is demonstrated by the low normalised rms error (commonly called FROG-error) F = 0.2%, which compares the differences of measured and retrieved traces. The original chirped XUV pulse is perfectly retrieved (see figures 2(b) and (d)), as well as the ensemble of NIR pulses with varying CEPs (see figure 2(a)). From the retrieved ensemble of gate pulses the CEP values of the 100 individual pulses in the ensemble are obtained and plotted in the histogram in figure 2(h). The rms of the retrieved distribution is σ retr. = 0.227π rad, which is within 3% of the original value. All this shows that the algorithm has successfully retrieved both the stable isolated attosecond pulse and the fluctuating dressing field. To further underscore this point, in (g) the average of the original (grey line) and the average of the retrieved ensemble of NIR pulses (red line) are shown together with the original and retrieved XUV pulse intensity envelopes (grey and blue lines, respectively). In addition to the commonly used FROG error F the agreement between two laser pulses can be quantified by an rms error for normalised fields NF (equation (6) in [33]). Within this formalism good agreement corresponds to NF < 0.02. For our original and retrieved attosecond pulses NF = 0.0076, showing that the attosecond pulse has been retrieved to very high accuracy.
Next, a stable 5 fs NIR dressing pulse in combination with an ensemble of XUV pulses with a group velocity dispersion (GVD) jitter is considered. The ensemble of 120 XUV pulses used to calculate the 'measured' streaking trace is constructed by sampling GVDs from a Gaussian distribution with σ = 0.037 fs 2 . This stretches the XUV pulses from the Fourierlimited duration of 314 as to an average FWHM duration of 755 as. This ensemble of 120 XUV pulses is displayed as grey lines in figure 3(b). A streaking trace for one of the XUV pulses is displayed in figure 3(c), whereas the incoherent sum of all 120 traces, representing the 'measured' spectrogram, is shown in figure 3(e). The td-e 2 PIE algorithm is now run for 150 iterations with L = 1, and K = 100. The retrieved trace is displayed in figure 3(f). Again, the rms error is small F = 0.2% indicating a very good convergence. Both the NIR dressing field and the ensemble of XUV pulses are correctly retrieved. The NIR electric field is plotted in figure 3(a). All retrieved XUV intensity envelopes |E XUV (t)| 2 are shown in figure 3(b). The original pulses are plotted as grey lines, the retrieved pulses as blue lines. All 100 pulses from the retrieved ensemble are shown. For better comparison the means are displayed in figure 3(g). The mean of the NIR streaking pulses as well as the mean of the XUV pulse ensemble agree perfectly. The normalised field error [33] for the comparison of the mean of the retrieved XUV pulse to the original XUV pulse mean is NF = 0.0055 indicating a high accuracy retrieval.  In figure 4 the case of an NIR gate with a fluctuating CEP in combination with an isolated attosecond XUV pulse with a fluctuating GVD is considered. In this case we run the tde 2 PIE algorithm with ensemble sizes of K = 10 pulses and L = 10 gates. The algorithm manages to retrieve the ensembles of both the NIR and XUV pulses. The retrieved trace closely resembles the measured trace (compare figures 4(e) and (f)). The rms trace error is F = 0.3%. In figures 4(a) and (b) the ensembles of the original NIR and XUV pulses are shown as grey lines. The ensemble of the retrieved gate pulses is shown as red lines in (a), whereas the ensemble of XUV pulses is plotted as blue lines in (b). In (d) all XUV pulse spectra are displayed (original in grey, retrieved in blue). The original phases are shown as grey dashed lines, the retrieved phases as blue dashed lines. In the subplot (g) the real part of the average of the original gate fields (grey) is compared to the real part of the average of the retrieved ensemble of gate pulses R{1/L L l=1 E l,NIR (t)} (red), as well as the modulus squared of the coherent average |1/K K k=1 P k (t)| 2 of the original and retrieved XUV pulse fields (blue line). The error for normalised fields [33] for the XUV pulses is NF = 0.0057. The mean GVD for the original ensemble is equal to the retrieved mean GVD μ orig. = μ retr. = 0.100 fs 2 , and the original and retrieved standard deviations of the GVD distributions are σ orig. = 0.023 fs 2 , and σ retr. = 0.035 fs 2 . The standard deviation of the retrieved CEPs is σ retr. = 0.16π rad compared to the standard deviation σ orig. = 0.14π rad of the original ensemble.
In stark contrast to the situation that we started with in figure 1, the average of the original XUV pulse fields that were used to calculate the streaking trace and the average retrieved XUV pulses agree very well with each other in the simulations presented in figures 2-4). Indeed, the td-e 2 PIE method introduced in this paper can also be used to correctly retrieve the attosecond pulse trains underlying the simulation shown in figure 1. Whereas the pulse retrieval in figure 1 led to the determination of an attosecond pulse train with satellite pulses that were dramatically weaker than the satellite pulses with which the RABBITT trace had been calculated, the retrieval in figure 5 shows an attosecond pulse train with the correct ratio between the intensities of the main peak and the satellites. The field error for this td-e 2 PIE retrieval is NF = 0.0737 compared to NF = 0.5252 for the retrieval using the single pulse algorithm in figure 1. It should be noted, that the simulations presented here have not included noise. We confirmed the td-e 2 PIE algorithm's robustness to noise by further simulations. Even for traces with a signal-to-noise level as low as 1, retrievals of an attosecond pulse with a satellite pulse resulted in pulse field errors below NF < 0.2. This is in agreement with Lucchini et al who showed that a pytchographic retrieval algorithm outperforms PCGPA and LSGPA in the correct retrieval of satellite pulses for noisy streaking traces [30].
It is interesting to note that the dressing field intensity in the examples shown in figures 1 and 5 was only 0.5 TW cm −2 . In [20] it was mentioned that an RABBITT measurement at low intensity is not a complete set of measurements that suffices to fully describe the quantum state [34], and therefore pulse ensembles cannot be retrieved correctly. We believe that the unsuccessful retrieval of APTs at low gate pulse intensity in [20] was rather caused by the way the APT was constructed using identical spectral components for each pulse in the train. This leads to a trace, which has a pure 2ω modulation, leading to a π ambiguity in the CEP of the dressing field. If we construct the APT using the more realistic assumption of varying XUV spectra for each constituent pulse, with the more intense pulses in the attosecond pulse train having higher energy cut-offs, then the trace not only has a 2ω, but crucially also a 1ω delay-dependent modulation, thus lifting the CEP ambiguity, as shown in our companion paper [17] on the retrieval of short APTs generated using 7 fs driver laser pulses.

Conclusion
In conclusion, we have introduced a retrieval algorithm based on mixed state time-domain ptychography, that allows to retrieve ensembles of dressing laser and attosecond pulse fields. This algorithm will be particularly useful for attosecond streaking measurements with an unstabilised or partially stabilised carrier envelope phase, or other experimentally present fluctuations.