Introduction

Trains of attosecond pulses generated by high-order harmonic generation in gases typically consist of only odd harmonics of the fundamental radiation of the near-infrared (NIR) driving field1. The absence of even harmonics is imposed by the inversion symmetry of the interaction when using atoms or ensembles of randomly oriented molecules as generating medium. Breaking the inversion symmetry, either by using an asymmetric driving field (e.g., by combining the fundamental field with its second harmonic2,3) or an asymmetric medium (e.g., by orienting asymmetric linear molecules4), produces even and odd harmonics of the fundamental radiation. In both cases, the harmonics are contiguous, i.e., regularly spaced throughout the spectrum. Additional degrees of freedom on the harmonic order generated in the strong-field process can be achieved by combining a driving field consisting of different harmonics with different polarisation states5,6 and orbital angular momentum7.

The experimental characterisation of the relative phase between successive odd harmonics is based on the reconstruction of attosecond beating by interfence of two-photon transition (RABBIT) approach8. In this technique, the relative phase between successive harmonics is encoded in the oscillation of the intensity of the photoelectron peaks generated by the absorption of a photon from the harmonic comb and the exchange of an additional photon with the synchronised NIR field. These peaks are usually referred to as sidebands to distinguish them from the main photoelectron lines due to the absorption of a single XUV photon. While the RABBIT technique is based on the exchange of a single NIR photon with the external field for each pathway contributing to the sideband signal, the effect of the absorption or emission of multiple NIR photons on the reconstructed attosecond pulse trains has also been discussed9. For the application of novel approaches in attosecond quantum metrology, the exchange of multiple NIR photons is also a fundamental point in order to connect the phases of non-adjacent harmonics10.

Recent developments in seeded free-electron lasers (FELs) operating in the extreme ultraviolet (XUV) spectral range11 offer new possibilities for the synthesis and shaping of attosecond pulse trains in the spectral ___domain, and for coherent control experiments in the XUV12,13,14,15,16. Selected harmonics can be generated independently in different undulators, overcoming the constraint of regularly-spaced harmonics. Moreover, the possibility to control the electron bunch dynamics between the different undulators offers the possibility to modify independently the relative phase between the harmonics on an attosecond timescale13.

In this work, we present the generation and control of attosecond pulse trains consisting of non-contiguous harmonics. Using the helium atom as a target, we study interference between multiphoton ionisation paths with differing parity (odd or even number of absorbed photons). While averaging out in angle-integrated measurements, these interferences are visible in the photoelectron angular distributions, which we accessed by collecting the photoelectrons emitted into the upper hemisphere, with positive z-component of the momentum (where the z-axis is defined by the light polarisation direction)17. We study the frequencies of the sideband oscillations in detail, isolating two different components contributed by multiple NIR photon transitions. The temporal reconstruction of the attosecond pulse trains is achieved by reconstructing the phase of the sideband oscillations using a correlation analysis17,18,19 and the attosecond timing tool15 for reordering the single-shot data. The analysis allows us to identify the contribution of a weak spurious harmonic to the oscillations of the sideband intensities.

Results

Temporal characterisation of attosecond pulse trains consisting of non-contiguous harmonics

The main idea of our approach and the scheme of the energy levels of the photoelectrons measured in the experiment are presented in Fig. 1a. The four harmonics H6, H7, H9 and H10 of the fundamental seed laser operating at λUV = 266 nm (angular frequency ωUV) were generated in the six undulators available at FERMI FEL-1 (see Supplementary Note 1: Experimental information). The seed laser was obtained by frequency triplication of a NIR driving field with angular frequency ωNIR. An intense NIR pulse (INIR ≈ 2 × 1012 W/cm2) with the same angular frequency ωNIR was collinearly recombined with the harmonics to generate sidebands of the main photoelectron peaks. The sidebands \({S}_{6,7}^{(\pm )}\) (\({S}_{9,10}^{(\pm )}\)) are expected to present oscillations at frequency 3ωNIR, due to the two interfering pathways contributed by the exchange of NIR photons following the XUV ionisation from the two consecutive harmonics H6 and H7 (H9 and H10). The mechanism leading to the population of these sidebands is exemplified for the sideband \({S}_{6,7}^{(-)}\) in Fig. 1b.

Fig. 1: Energy scheme for four non-consecutive harmonics and intense NIR field.
figure 1

a Scheme of energy levels of the photoelectrons released by the absorption of one photon from the four non-consecutive harmonics H6, H7, H9 and H10 and the absorption or emission of (up to five) NIR photons. Ip indicates the ionisation potential of the atom. b Photoionisation pathways leading to the population of the sideband \({S}_{6,7}^{(-)}\) and contributed by the harmonics H6 and H7. Due to the number of NIR photons involved, the population oscillates at frequency 3ωNIR. c Photoionisation pathways leading to the population of the sideband \({S}_{7,9}^{(--)}\) and contributed by the harmonics H6, H7 and H9. Due to the number of NIR photons involved, the population oscillates at frequencies 3ωNIR  and 6ωNIR as a result of the paths contributed by the harmonics H6-H7 and H7-H9). The subscripts in the sideband notation indicate the harmonic orders closest in energy to the given sideband. The superscripts (− , + ) and (− − , − , 0, + , + + ) are used to distinguish sidebands in order of increasing energy between adjacent and non-adjacent harmonics, respectively. Only transitions involving the exchange of up to five NIR photons are shown for simplicity.

Full size image

Similarly, the sidebands \({S}_{7,9}^{(--)}\) and \({S}_{7,9}^{(-)}\) (\({S}_{7,9}^{(+)}\) and \({S}_{7,9}^{(++)}\)) between the harmonics H7 and H9 are expected to show oscillation frequencies at 3ωNIR due to the pathways contributed by the harmonics H6 and H7 (H9 and H10) and the absorption (emission) of NIR photons, as shown on the left-hand side of Fig. 1c. An additional interference term at 3ωNIR of \({S}_{7,9}^{(--)}\) and \({S}_{7,9}^{(-)}\) (\({S}_{7,9}^{(+)}\) and \({S}_{7,9}^{(++)}\)), contributed by the harmonics H9 and H10 (H6 and H7), is expected to be negligible in the NIR intensity range of the experiment, due to the large number of NIR photons required to populate these sidebands after absorption of one photon of the harmonics H9 or H10 (H6 or H7).

At the same time, all sidebands between the main photoelectron lines H7 and H9 should present an additional term oscillating at 6ωNIR, due to the interference of the pathways involving the harmonics H7 and H9 and the exchange of up to five NIR photons. The two dominant interference mechanisms leading to the oscillation frequencies 3ωNIR and 6ωNIR are shown for the sideband \({S}_{7,9}^{(--)}\) in Fig. 1c. Finally, also the central sideband \({S}_{7,9}^{(0)}\) should present interfering terms oscillating at 3ωNIR and 6ωNIR.

For the temporal reconstruction of the attosecond pulse trains, the relative phase between the harmonics H6 and H7 (H9 and H10) can be determined by investigating the oscillations at frequency 3ωNIR of the sidebands \({S}_{6,7}^{(\pm )}\) (\({S}_{9,10}^{(\pm )}\)). The relative phase between the harmonics H7 and H9 is encoded in the oscillations at frequency 6ωNIR of one of the central sidebands \({S}_{7,9}^{(--,-,0,+,++)}\). By combining this information with an independent measurement of the harmonic intensities we aim to reconstruct the intensity profile of the attosecond pulse train.

XUV harmonic combs and experimental setup

In the experiment, the six undulators available at FERMI (U1–U6) were tuned to the harmonics H10, H9, H7 and H6 of the seed laser. In particular, the first undulator produced harmonic H10, the second H9, the third and fourth H7 and the last two (fifth and sixth) H6 (see Fig. 2a). The relative phase between the harmonics was adjusted by means of phase shifters (PSs) positioned between each pair of successive undulators. In particular, the phase shifter PS1 was used in the experiment to advance the harmonic H10 by the time interval τs1 while PS4 was used to delay the harmonic H6 by the delay τs4 with respect to the other harmonics.

Fig. 2: Harmonic spectra.
figure 2

a Schematic view of the correspondence between the six undulators (U1–U6) of FERMI and the generated harmonics H10-H9-H7-H6. The phase shifters PS1-PS5 placed between the undulators are also indicated, together with the delays introduced by PS1 (τs1) and PS4 (τs4). b Harmonic spectra acquired for two different settings (dataset 1 black dashed line; dataset 2 red dotted line) of the phase shifter PS4, corresponding to a ≈ π phase difference of the phase of the harmonic H6. The intensity of H8 is roughly two orders of magnitude weaker than that of the other harmonics.

Full size image

Two typical XUV spectra acquired with an XUV spectrometer positioned after the end station are shown on a logarithmic scale in Fig. 2b. The difference between the black and red spectra, corresponding to two different experimental settings hereafter referred to as dataset 1 and dataset 2, respectively, was a shift of π rad in the phase of the last harmonic (H6) introduced by the phase shifter PS4. The four harmonics H6-H7-H9 and H10 have comparable intensities (see Table S1 in the Supplementary Information for the corresponding values). It can be observed that a weak harmonic H8 is also present in the experimental spectra. Its origin can be ascribed to the seeding mechanism of the FEL  process (see Methods Section). It should also be noted that in the experiment we could not independently control the phase of this harmonic.

We will show that, even though the presence of the weak harmonic H8 does not significantly affect the reconstructed attosecond pulse train, its presence can significantly modify the oscillations of the sidebands \({S}_{7,9}^{(\pm \pm ,\pm ,0)}\) between the harmonics H7 and H9.

Correlation plots

To reconstruct the intensity profile of the attosecond pulse train, we need to access the relative phase of the harmonics. Due to the timing jitter between the attosecond waveform and the NIR field20, the intensities of all sidebands vary from shot to shot in a correlated way determined by the phases of each harmonic. The typical jitter is in the order of  5 to 6 femtoseconds20, which is greater than the expected periods of the sideband oscillations (T = 885 as and T = 442 as for the 3ωNIR and 6ωNIR oscillation frequencies, respectively).

We start by considering the three lowest harmonics (H6, H7 and H9) and the sidebands produced by the NIR pulse between them. The intensity of the sidebands \({S}_{6,7}^{(\pm )}\) oscillates depending on the relative delay τ between the XUV and NIR field, and on the differences between the phases of the harmonics H7 and H6 (φ7 and φ6) according to the phase term (3ωNIRτ + φ7 − φ6). The oscillatory component of the sideband signal for photelectrons emitted into the upper half space can be isolated by considering the ratio between the difference and the sum of the individual sidebands. The resulting term is given by the expression18:

$${P}_{6,7}=\frac{{S}_{6,7}^{(+)}-{S}_{6,7}^{(-)}}{{S}_{6,7}^{(+)}+{S}_{6,7}^{(-)}}={\alpha }_{6,7}\cos \left[3{\omega }_{{{{\rm{NIR}}}}}\tau +{\varphi }_{7}-{\varphi }_{6}\right],$$
(1)

where α6,7 is a constant.

For the central sideband \({S}_{7,9}^{(0)}\) (corresponding in energy to the position of the weak harmonic H8), the oscillation at 6ωNIR is expected to dominate. According to the strong-field approximation (SFA) model18, the intensity varies according to the relation:

$${S}_{7,9}^{(0)}\approx {a}_{7,9}-{b}_{7,9}\cos (6{\omega }_{{{{\rm{NIR}}}}}\tau +{\varphi }_{9}-{\varphi }_{7}),$$
(2)

where a7,9 and b7,9 are constants. To create a correlation plot of the single-shot photoelectron spectra where both components have the same oscillation frequency with the delay τ, we consider the square of the oscillating component P6,7. The correlation plot \(\left.\right({P}_{6,7}^{2}\) vs \({S}_{7,9}^{(0)}\left.\right)\) results in an ellipse, whose ellipticity depends on the phase difference:

$${\eta }_{L}={\varphi }_{9}-3{\varphi }_{7}+2{\varphi }_{6}.$$
(3)

Similarly, for the group of highest harmonics (H7-H9-H10), we find that the shape of the correlation plot \(\left.\right({P}_{9,10}^{2}\) vs \({S}_{7,9}^{(0)}\left.\right)\) depends on the phase:

$${\eta }_{H}=2{\varphi }_{10}-3{\varphi }_{9}+{\varphi }_{7}.$$
(4)

The reconstructed intensity profile of the attosecond pulse trains (under the assumption of monochromatic harmonics) depends only on the intensity of each harmonic and on the two phases ηL and ηH (see Supplementary Note 2: Temporal reconstruction of the attosecond pulse trains with non-consecutive harmonics).

The phases ηL and ηH are proportional to the group delay dispersion (GDD) of the harmonic groups H6-H7-H9 and H7-H9-H10, respectively. Indeed, neglecting the phase variation within a single harmonic, the GDD for the lowest three harmonics is given by:

$${{{\rm{GDD}}}}=\frac{{d}^{2}\varphi }{d{\omega }^{2}}\propto \frac{1}{{\omega }_{{{{\rm{UV}}}}}}\left[\frac{{\varphi }_{9}-{\varphi }_{7}}{2{\omega }_{{{{\rm{UV}}}}}}-\frac{{\varphi }_{7}-{\varphi }_{6}}{{\omega }_{{{{\rm{UV}}}}}}\right]\propto {\eta }_{L}$$
(5)

The value of these phases can be determined by exploiting the phase tunability of each individual harmonic of the comb and calculating the correlation parameters ρ6,7,9 and ρ7,9,10 (see definition in Methods Section) of the correlation plots obtained by varying the two phases ηL and ηH. According to Eqs. (1) and (2), the correlation parameters oscillate as a function of the second-order phase difference as:

$${\rho }_{6,7,9}=-\cos ({\eta }_{L});\quad {\rho }_{7,9,10}=-\cos ({\eta }_{H})$$
(6)

The value of the correlation parameter allows one to estimate the ellipticity of the correlation plot and thus access the second-order phase difference of the three harmonics.

The experimental results are presented in Fig. 3a, b for the correlation parameters ρ6,7,9 and ρ7,9,10, respectively. For the first correlation plot, the phase of the harmonic H6 was changed using the phase shifter PS4. The delay τs4 introduced by the phase shifter affects the phase ηL through the term ηL = −12ωUVτs4 (see Supplementary Note 2: Temporal reconstruction of the attosecond pulse trains with non-consecutive harmonics). Since the phase of all other harmonics is not changed, the phase ηH remains unaffected. The correlation parameter ρ6,7,9 calculated from the correlation plots for each phase setting shows a clear oscillation with frequency 6ωNIR, which is very well fitted by a sinusoidal curve. By fixing the minima of this fit to the values ηL = 2nπ, we obtain a value of the second-order phase difference for each delay introduced by the PS4.

Fig. 3: Calibration curves for four non-consecutive harmonics.
figure 3

a, b Variation of the correlation coefficients ρ6,7,9 (ρ7,9,10) calculated on the correlation plots \({P}_{6,7}^{2}\) vs \({S}_{6,7}^{(0)}\) (\({P}_{9,10}^{2}\) vs \({S}_{6,7}^{(0)}\)) as a function of the delay τs4 (τs1) introduced by the phase shifter PS4 (PS1). The upper axis indicates the calibrated phase axis for the phase difference ηL (ηH). The error bars are calculated as the standard deviation of the mean values for five sets of data (each with typically 1500 shots) sampled from the total shots of data for each phase-shifter value.

Full size image

A similar procedure was followed for the calibration of the phase ηH using the phase shifter PS1, which only affects the phase of the harmonic H10, leaving the phases of the three lowest harmonics (and therefore the phase ηL) unaffected (see Fig. 3b). With this procedure we can determine the values of the two phases ηL and ηH for each position of the phase shifters. In our previous work we showed that the NIR intensity strongly affects the quality of the correlation plots18. The NIR field must be strong enough to overcome single-shot fluctuations in harmonic amplitudes, and not so strong that it induces strong high-order multiphoton transitions affecting the reconstruction method. These conditions are met for NIR peak intensities between  ≈ 1.4 × 1012 W/cm2 and  ≈ 5 × 1012 W/cm215,18. The NIR intensity used for the acquisition of the experimental data (INIR ≈ 2 × 1012 W/cm2) is well within this range.

Attosecond timing tool and sideband oscillations

An alternative method for reconstructing the second-order phase differences ηL and ηH is to reorder the single-shot experimental data using the attosecond timing tool demonstrated in ref. 15. Using this approach, the sideband oscillations can be reconstructed and the phase differences can be extracted directly from the sinusoidal fits of the sideband oscillations. We use the oscillating components of the sidebands P6,7 and P9,10 with a relative phase difference of about π/2 to implement the attosecond timing tool. This phase difference between the two oscillating components leads to a circular correlation plot \(\left.\right({P}_{9,10}\) vs \({P}_{6,7}\left.\right)\) which allows us to minimise the error in the reconstruction of the sideband oscillations15.

The phase Φ retrieved from the correlation plot is given by:

$$\Phi =3{\omega }_{{{{\rm{NIR}}}}}\tau +\eta$$
(7)

with

$$\eta =\frac{1}{2}\left[({\varphi }_{10}-{\varphi }_{9})+({\varphi }_{7}-{\varphi }_{6})\right]$$
(8)

The oscillating components of the sidebands \({S}_{6,7}^{(\pm )}\) and \({S}_{9,10}^{(\pm )}\) and the oscillation of the central sidebands \({S}_{7,9}^{(0)}\) can be expressed as:

$${P}_{6,7}={\alpha }_{6,7}\cos \left[\Phi -\frac{1}{4}\left({\eta }_{L}+{\eta }_{H}\right)\right]$$
(9)
$${P}_{9,10}={\alpha }_{9,10}\cos \left[\Phi +\frac{1}{4}\left({\eta }_{L}+{\eta }_{H}\right)\right]$$
(10)
$${S}_{7,9}^{(0)}\approx {a}_{7,9}-{b}_{7,9}\cos \left[2\Phi +\frac{1}{2}\left({\eta }_{L}-{\eta }_{H}\right)\right],$$
(11)

where α9,10 is a constant (in general different from α6,7). By fitting the experimental oscillations of the sidebands, and using Eqs. (9) and (11) (or Eqs. (10) and (11)), the values of ηL and ηH can be directly extracted.

As an example, the reconstructed oscillations for the dataset 1 are shown in Fig. 4. As expected, the sidebands between successive harmonics \({S}_{6,7}^{(\pm )}\) (Fig. 4a, b) and \({S}_{9,10}^{(\pm )}\) (Fig. 4c, d) are well fitted by a sinusoidal curve with frequency 3ωNIR (red line). The central sideband between the non-consecutive harmonics \({S}_{7,9}^{(0)}\) (Fig. 4g) oscillates at frequency 6ωNIR due to the two paths contributed by the harmonics H7 and H9 and involving the exchange of (in total) six NIR photons.

Fig. 4: Reconstructed sideband oscillations using the attosecond timing tool for the dataset 1.
figure 4

Reconstructed experimental oscillations of photoelectrons emitted into the upper half space with respect to the light polarization direction (square black points and error bars) of the sideband \({S}_{6,7}^{(-)}\) (a), \({S}_{6,7}^{(+)}\) (b), \({S}_{9,10}^{(-)}\) (c), \({S}_{9,10}^{(+)}\) (d), \({S}_{7,9}^{(--)}\) (e), \({S}_{7,9}^{(-)}\) (f), \({S}_{7,9}^{(0)}\) (g), \({S}_{7,9}^{(+)}\) (h), and \({S}_{7,9}^{(++)}\) (i) as a function of the phase Φ (see Eq. (7)) for the dataset 1. The red curves are the double cosine fit of the experimental oscillations (see Eq. (12)). The error bars in the plots are calculated as the standard deviation of the sideband amplitudes.

Full size image

The other sidebands between harmonics H7 and H9 \({S}_{7,9}^{(--,++)}\) (Fig. 4e, i) and \({S}_{7,9}^{(-,+)}\) (Fig. 4f, h) show a more complex pattern characterised by the combination of two oscillating components at frequencies 3ωNIR and 6ωNIR. The combination of the two oscillations is particularly evident for the sidebands \({S}_{7,9}^{(\pm )}\), while the other two (\({S}_{7,9}^{(--,++)}\)) show a more pronounced oscillation at 3ωNIR. The experimental data for the dataset 2 (see Fig. 2a) are shown in Fig. 5. 

The comparison of the ηL and ηH values for the two datasets is presented in Table S2, showing that the correlation plot and attosecond timing tool methods return consistent values within the error bars.

Strong field approximation

The photoelectron spectra were simulated using the SFA21,22. The model has already been described in detail in refs. 15,18. In the simulations we consider a hydrogenoid transition matrix dipole moment d from the ground state to the continuum state. The pulse duration of each XUV harmonic was set to 50 fs and the duration of the NIR pulse was 65 fs. Furthermore, the total photoelectron yield has been integrated over 2π srad along the polarisation direction of the XUV and NIR fields, to reproduce the experimental conditions.

Time dependent Schrödinger equation

We solve the full-dimensional time-dependent Schrödinger equation (TDSE) for the helium atom, employing the nearly-exact approach thoroughly described in refs. 23,24,25. In brief, we use a numerical description of the two-electron wave function using a finite element discrete variable representation (FEDVR) for the discretization of the radial functions and coupled spherical harmonics for the angular part. The TDSE is solved using a highly efficient short iterative Lanczos algorithm23, and the resulting ionisation probabilities are extracted by a direct projection onto the final scattering states once the external field is turned off. The laser-atom interaction term is written within the dipole approximation, employing linearly polarized laser pulses with a sine-squared envelope to account for the finite duration.

We employ the experimental conditions, i.e., we define the harmonics q−2, q−1, q + 1, and q + 2—for q = 8—of the seed laser with λUV = 266 nm, with 24-fs duration. We then retrieve the RABBIT spectra by performing a time-delay scan relative to a 798 nm NIR field whose duration is set to 28.8 fs, scanning over delays τ ∈ [−887, 887] as (containing two full periods of the lowest-frequency oscillation with frequency 3ωNIR).

Since we are interested in single ionisation, and not too close to the ionisation threshold, a reasonably accurate description can also be obtained by using a single-active electron (SAE) approximation, significantly reducing the computational effort. SAE results are not explicitly shown nor discussed here, but were used to confirm that the sideband phases are essentially independent of the harmonic intensities and quite robust to variations of the NIR intensity close to the experimentally used values, indicating that the lowest-order contribution to each oscillatory component dominates.

Discussion

In order to compare the experimental data with the theoretical predictions, the amplitudes and second-order phase differences of the harmonics H6-H7-H9-H10 (see Tables S1 and S2) were used as input values in the SFA model and the TDSE ab initio simulations. The TDSE simulations here were performed with the same harmonic phases as extracted from the SFA model (which does not account for the phases of the atomic transition matrix elements). The phase shifts between SFA and TDSE results thus give an indication of the magnitude of the atomic phase contribution to the multiphoton ionisation. The experimental and theoretical sideband oscillations as a function of the relative phase Φ were fitted with a double cosine function including the 3ωNIR and 6ωNIR components to take into account the different mechanisms determined by the exchange of multiple NIR photons (see Fig. 1):

$${S}_{{{{\rm{fit}}}}}={A}_{0}+{A}_{3\omega }\cos (\Phi +{\varphi }_{3\omega })+{A}_{6\omega }\cos (2\Phi +{\varphi }_{6\omega }),$$
(12)

where A0, A3ω and A6ω represent the amplitudes of the constant component and those oscillating at frequency 3ωNIR and 6ωNIR, respectively. Similarly, φ3ω and φ6ω represent the corresponding phases in the fit.

The experimental fits with the double cosine are shown in Figs. 4 and 5 by red lines, for the datasets 1 and 2, respectively. The curves agree well with the experimental points, confirming the presence of two oscillatory components.

Fig. 5: Reconstructed sideband oscillations using the attosecond timing tool for the dataset 2.
figure 5

Reconstructed experimental oscillations of photoelectrons emitted into the upper half space with respect to the light polarization direction (square black points and error bars) of the sideband \({S}_{6,7}^{(-)}\) (a), \({S}_{6,7}^{(+)}\) (b), \({S}_{9,10}^{(-)}\) (c), \({S}_{9,10}^{(+)}\) (d), \({S}_{7,9}^{(--)}\) (e), \({S}_{7,9}^{(-)}\) (f), \({S}_{7,9}^{(0)}\) (g), \({S}_{7,9}^{(+)}\) (h), and \({S}_{7,9}^{(++)}\) (i) as a function of the phase Φ (see Eq. (7)) for the dataset 2. The red curves are the double cosine fit of the experimental oscillations (see Eq. (12)). The error bars in the plots are calculated as the standard deviation of the sideband amplitudes.

Full size image

The comparison between experiment and simulations for the two phases φ3ω and φ6ω and the ratio between the amplitude components A3ω/A6ω is shown in Figs. 6 and 7 for dataset 1 and 2, respectively. The corresponding values are reported in Tables S3 and S4 for dataset 1 and 2, respectively. We focus the discussion on the comparison for dataset 1. Similar conclusions apply to the second dataset.

Fig. 6: Comparison of the experimental and simulated amplitudes and phases of the oscillation components of the sidebands for dataset 1.
figure 6

Comparison of the phases φ3ω (a, d), φ6ω (b, e) and ratio of the amplitudes A3ω/A6ω (c, f) derived from Eq. (12) for the experimental dataset 1 (black squares), SFA simulations (blue symbols) and TDSE simulations (red symbols), without (a–c) (circles) and with (d–f) (diamonds) the contribution of the harmonic H8. The error bars were obtained by error propagation considering the errors in the estimation of the phases and of the amplitudes in the two-cosine fit. The data of the XUV field are reported in Tables S1 and S2. The data are summarised in Table S3.

Full size image
Fig. 7: Comparison of the experimental and simulated amplitudes and phases of the oscillation components of the sidebands for dataset 2.
figure 7

Comparison of the phases φ3ω (a, d), φ6ω (b, e) and ratio of the amplitudes A3ω/A6ω (c, f) derived from Eq. (12) for the experimental dataset 2 (black squares), SFA simulations (blue symbols) and TDSE simulations (red symbols), without (a–c) (circles) and with (d–f) (diamonds) the contribution of the harmonic H8. The error bars were obtained by error propagation considering the errors in the estimation of the phases and of the amplitudes in the two-cosine fit. The data of the XUV field are reported in Tables S1 and S2. The data are summarised in Table S4.

Full size image

For dataset 1, the agreement between the phases of the 3ωNIR components (Fig. 6a) extracted from the experimental results and the SFA simulations for the sidebands \({S}_{6,7}^{(\pm )}\) and \({S}_{9,10}^{(\pm )}\) is excellent. On the contrary, large deviations between experiment and SFA simulations are observed for the phases φ3ω for the five sidebands between harmonics H7 and H9.

The agreement between the SFA model and the experiment for the phases φ6ω of the central sidebands is very good (see Fig. 6b). We do not present the comparison for the phases φ6ω for the other sidebands, since the oscillations at these frequencies are not clearly visible in the experimental data and cannot be extracted reliably.

The agreement for the ratios is also reasonable, although the SFA model tends to underestimate the ratios for the central sidebands (see Fig. 6c).

The evolution of the phases obtained with the SFA model and the TDSE simulations show some differences. As mentioned above, the TDSE simulations here were performed with the harmonic phases as extracted from the SFA model, such that the difference in the obtained sideband phases gives an indication of the contribution of the atomic transition phases. This conclusion is supported by the observation that, in the case of the φ3ω phase, the differences decrease with increasing photoelectron energy. The differences between the φ6ω phases show an almost constant offset due to the multiple continuum-continuum phase terms contributing to the overall phase of these sideband oscillations26,27,28. Instead of using the SFA-based approach for extracting the harmonic phases discussed above, it is also possible to extract these parameters by nonlinear fitting of all pulse parameters to achieve good agreement between the experimental results and the full TDSE simulations. This approach does not rely on perturbative approximations and includes the atomic phases, but requires many computationally expensive TDSE simulations (performed within the SAE here). Below we show that this gives a very similar reconstructed attosecond pulse train, indicating that the SFA-based extraction method is accurate as long as the target atomic phases are sufficiently smooth as a function of energy.

The term oscillating at frequency 3ωNIR for the central sidebands between the harmonics H7 and H9 can be populated through the interference of paths involving the absorption of an XUV photon from each harmonic H6 and H7 (for the sidebands \({S}_{7,9}^{(--)}\) and \({S}_{7,9}^{(-)}\)) or from each harmonic H9 and H10 (for the sidebands \({S}_{7,9}^{(+)}\) and \({S}_{7,9}^{(++)}\)) and the exchange of several NIR photons (see Fig. 1c). The discrepancies observed for the phases φ3ω for these sidebands indicate that these pathways are not sufficient to reproduce the experimental phases.

These discrepancies can be resolved by taking into account that in our experimental conditions, the weak harmonic H8 also contributes to the oscillations at frequency 3ωNIR with an additional path. This path corresponds to the interference mechanism involving harmonics H7-H8 for the lowest central sidebands and harmonics H8-H9 for the two highest central sidebands. The results including the effect of the weak harmonic H8 are shown in Fig. 6d–f. The amplitude of the field associated with harmonic H8 has been derived from the XUV spectrometer measurements (see Fig. 2a and Table S1). The value of the phase of harmonic H8 was derived by optimising the agreement between the experimental values φ3ω and the results of the SFA simulations, and keeping the phases ηL and ηH fixed at the experimental values.

We can observe that by including the contribution of H8, the comparison of the experimental and simulated phases φ3ω is very good for all the sidebands except \({S}_{7,9}^{(0)}\) (Fig. 6d). The presence of the harmonic H8 has a very weak effect (of the order of a few attoseconds) on the simulated phases of the oscillations of the sidebands \({S}_{6,7}^{(\pm )}\) and \({S}_{9,10}^{(\pm )}\), while its effect is significant for the sidebands \({S}_{7,9}^{(++,--)}\), \({S}_{7,9}^{(\pm )},\) and \({S}_{7,9}^{(0)}\).

As shown in the Supplementary Discussion, in general four different terms oscillating at the frequency 3ωNIR contribute to the signal of the central sidebands, when considering the five harmonics H6-H10. The magnitude of these terms depends on the product of the amplitudes of the harmonics contributing to them and the product of Bessel functions, whose order depend on the number of NIR photons exchanged with the NIR field. In our experimental conditions, the most relevant terms are from the pathways contributed by the harmonics H6-H7 and H7-H8 for the sidebands \({S}_{7,9}^{(--)}\) and \({S}_{7,9}^{(-)}\), and by the harmonics H8-H9 and H9-H10 for the sidebands \({S}_{7,9}^{(+)}\) and \({S}_{7,9}^{(++)}\), respectively. In particular, the sidebands closer to the harmonic H8 (\({S}_{7,9}^{(\pm )}\)) are dominated by the contribution of the interference H7-H8 for the sideband \({S}_{7,9}^{(-)}\) and that of H8-H9 for \({S}_{7,9}^{(+)}\). These contributions are about one order of magnitude higher than the other terms (see Table S5). For the sidebands \({S}_{7,9}^{(--)}\) and \({S}_{7,9}^{(++)}\), the contributions of the term H7-H8 and H8-H9 are still the strongest, but other terms have comparable magnitude (H6-H7 for \({S}_{7,9}^{(--)}\) and H9-H10 for \({S}_{7,9}^{(++)}\)) (see Table S5).

For the central peak \({S}_{7,9}^{(0)}\), all four terms present comparable magnitude (see Table S5) and the phase of the oscillation of this photoelectron peak is extremely sensitive to their ratios. An experimental error in the estimation of the amplitude of the harmonics could be responsible for the limited agreement between the experimental data and the theoretical simulations in Fig. 6d.

The inclusion of the harmonic H8 in the simulations only marginally affects the phase of the 6ωNIR oscillations for the central sidebands, as shown in Fig. 6e. This trend is expected considering that the 6ωNIR contributions due to the harmonic H8 (H6-H8 for the lowest sidebands and H8-H10 for the highest sidebands) are negligible due to the low amplitude of H8. This is confirmed also by the quantitative analysis considering the three terms oscillating at frequency 6ωNIR (see Supplementary Discussion and Table S6).

Finally, the inclusion of harmonic H8 also leads to a better agreement of the ratio of the amplitudes of the A3ω and A6ω components (see Fig. 6f), at least for the TDSE simulations. For the sake of completeness, we also report the contributions of the components oscillating at frequencies 9ωNIR and 12ωNIR in Tables S7 and S8, respectively. These components, however, do not contribute significantly to the sideband yield under our experimental conditions. The values shown in Fig. 6 are presented in Table S3.

In Fig. 7 we report the same analysis obtained for dataset 2, where only the phase of the harmonic H6 was changed by π rad with respect to the phase of the harmonics used for the sideband measurements presented in Fig. 6. The comparison between experiment and theory qualitatively reproduces the conclusions derived for the dataset 1 with a good agreement for the phases φ3ω for the sidebands between successive harmonics and a poor agreement for the central sidebands when H8 is not included (Fig. 7a). Similar conclusions for the phases φ6ω and the ratios A3ω/A6ω remain valid also for dataset 2 (Fig. 7b, c). Also in this case, the inclusion of the contribution of the harmonic H8 significantly improves the agreement for the phases φ3ω of the central sidebands (Fig. 7d), without significantly affecting the phases φ6ω (Fig. 7e). The harmonic H8 also leads to an overall improvement in the experimental-theoretical comparison for the ratios (Fig. 7f). The values shown in Fig. 7 are presented in Table S4.

We observe that SFA simulations performed with a lower intensity of the harmonic H8 (between 0.1% and 0.5% of the intensity of the harmonic H7) show a slightly better agreement with the experimental data for both datasets (not shown), suggesting that the intensity of H8 might be slightly overestimated. This could be due to a different scaling of the harmonic signals measured with the XUV spectrometer between the intense harmonics (H6-H7-H9-H10) and the less intense one (H8), due to their large intensity difference (between two and three orders of magnitude).

For the temporal reconstruction of the attosecond pulse trains consisting of the harmonics H6-H7-H9-H10, we used the expression described in the Supplementary Note 2: Temporal reconstruction of the attosecond pulse trains with non-consecutive harmonics and the amplitudes and phases contained in Tables S1 and S2. In order to include the contribution of the harmonic H8, we followed the approach demonstrated in refs. 17,29, which requires the second-order phase differences Δφq−1,q,q+1 = φq+1−2φq + φq−1 between three consecutive harmonics.

To account for the photoionisation phases included in the TDSE simulations in the temporal reconstruction, we performed an optimisation procedure to fit the phases of all the sideband oscillations and extracted the second-order phase differences from the harmonic combination that best fits the experimental data (see Supplementary Discussion). We remark that for these TDSE simulations, we did not use SFA-extracted phases, but the phases were obtained by fitting the TDSE simulations directly to the experimental data. The second-order phase differences obtained with this approach are presented in Table S9 and differ by about  ≈ 0.06–0.2 rad from those obtained from the SFA model. This result is in agreement with the previously reported contribution of the photoionisation phase17.

The reconstructed intensity profiles of the attosecond pulse train using the SFA model (with and without a contribution of the harmonic H8) and the TDSE fitting procedure are shown in Fig. 8 for the two datasets. We observe that the harmonic H8 does not contribute significantly to the intensity profile due to its low amplitude. Furthermore, the TDSE reconstruction differs only slightly from the SFA-based reconstruction, due to the small differences in the second-order phase differences. Finally, we observe that the π phase change of the harmonic H6 between the two datasets (Fig. 8a,.b for dataset 1 and Fig. 8c, d for dataset 2) strongly affects the intensity profile, as expected.

Fig. 8: Attosecond pulse reconstruction.
figure 8

Comparison of the reconstructed attosecond pulse using the reconstruction formulas (see Eq. S3 and refs. 29) (black dots) and the TDSE best fitting procedure (red dotted lines) (see Table S9) for the dataset 1 (a, b) and dataset 2 (c, d), including (F8 ≠ 0 a and c) or excluding (F8 = 0 b and d) the contribution of the harmonic H8. The blue shaded areas are the errors in the reconstruction estimated as the standard deviation of the reconstructions obtained for 500 input amplitudes and phases of the harmonics randomly chosen within their respective error bars. F8 stands for the amplitude of H8.

Full size image

The double-peak structure observed in the reconstructed attosecond intensity profile is due to the second-order phase difference between the low (H6-H7) and high (H9-H10) pair of harmonics, which is close to  −π/2 (panels a and b) and π/2 (panels c and d).

The typical uncertainty of our technique for the temporal duration of an attosecond pulse in the train is about 2 to 3%, corresponding to 4 to 5 as. This value is smaller than the typical uncertainty obtained with the angular streaking method used for the temporal characterisation of isolated attosecond pulses30,31,32, which is of the order of 20% or larger (see Supplementary Information in ref. 33). We attribute this difference mainly to the fact that in angular streaking the pulse duration is obtained from each single-shot measurement, whereas in our approach, thanks to the shot-to-shot reproducibility of the harmonic properties, the characterisation is based on several (typically thousands) of single-shot experimental spectra. For this reason, the correlation approach leads to a significant reduction in the uncertainty of the pulse duration.

Conclusion

We have shown that seeded FELs offer the unique possibility of generating selected harmonics of the driving field, including non-equally spaced harmonics. For the temporal reconstruction of the attosecond pulse train, we have taken advantage of the technique already demonstrated for consecutive harmonics29. The two different approaches (correlation analysis and sideband reordering) provide consistent results for the second-order phase difference of the harmonic combs and thus for the attosecond pulse reconstruction.

A non-contiguous harmonic comb, like the one demonstrated in this work, could be useful, for example, in studying attosecond time delays in photoionisation in molecules. In these systems, the close energy spacing between the different electronic states of the cation typically leads to congested photoelectron spectra in which the contributions of the different photoionisation channels cannot be easily disentangled. The use of non-consecutive and tunable harmonics offers two important experimental tools not only to reduce the complexity of the photoelectron spectra, but also to address specific resonances in the molecular systems. The generation of attosecond pulses at FELs will pave the way for the study of ultrafast electronic dynamics in the soft X-ray34 and for the implementation of attosecond pump-attosecond probe experimental approaches35. The combination of these elements with shaping capabilities will provide the international community with unique tools for the study and control of ultrafast electronic dynamics.

Methods

Experimental methods

The interaction between the electron beam and the UV external laser induces a periodic density modulation in the beam at the seed laser wavelength and its higher harmonics, up to H10 and beyond. The radiator undulators of the FEL are then tuned to a subset of these harmonics, providing strong coherent emission at these photon energies. In this experiment, the selected harmonics were H6, H7, H9 and H10, producing approximately a few microjoules (μJ) per color.

However, the bunched electron beam can coherently emit at all harmonics in any magnetic element defining the electron beamline, such as the terminations along the FEL amplifier or the dipole that deflects the beam to the beam dump. This can lead to a spurious signal, several orders of magnitude weaker than the emission produced in the resonant undulators, and is the source of the signal at H8 detected during the experiment.

The XUV harmonics were transported along the PADReS system36, where they were also characterised in terms of intensity and spectrum. Moreover, they were focused into the experimental end station in a helium gas target by means of the KAOS active optics system37 serving the LDM beamline. The focusing was optimised monitoring the focal spot size and shape using a scintillator crystal positioned in the interaction region. As the photon beam arriving on the sample during the experiment was composed by different harmonics, with the harmonics’ effective sources displaced along the undulator chain, the ‘source—refocusing optics’ distances were also different. As a consequence, it was impossible to optimally focus all the harmonics at the same time. However, as the request for the experiment was to have a non-tightly focused spot, a good compromise was found in order to produce similar spots for the separate harmonics in the interaction region, with their sizes being about 60 × 60 μm2 FWHM. In any case, simulations showed that the focal spot size variation upon moving the photon source from the first to the sixth undulator was not greater than 20%. The XUV harmonics were collinearly recombined with a NIR field using a drilled mirror. The photoelectrons emitted in the two-color photoionisation process were collected in a solid angle of 2π srad along the common polarisation direction of the fields by a magnetic bottle electron spectrometer. Photoelectron spectra were measured for each individual laser shot in the experiment and used to reconstruct the correlation plots. After the interaction region, the XUV radiation was analysed using an XUV spectrometer38.

Theoretical methods

To ensure numerical convergence, the TDSE simulations have been carried out using an expansion in coupled spherical harmonics including total angular momenta up to \({L}_{\max }={l}_{\max }=12\), both for the single and total angular momenta, limiting the individual angular momenta to \(\min ({l}_{1},{l}_{2})\le 5\). The radial box used for the calculations is relatively large, \({r}_{\max }=2248\) a.u. to avoid unphysical reflections during the propagation, and with a grid density defined by finite elements of 4 a.u., and DVR functions of order 11. Since only single ionisation processes in the 1s channel are relevant, we use an L-shaped grid restricted to \(\min ({r}_{1},{r}_{2})\le 32\) a.u.

For the SAE simulations, we have employed the model potential for He from ref. 39. The SAE implementation also uses a FEDVR numerical representation of the wave function, but solves the TDSE using a split-operator propagator.

Data analysis methods

The correlation coefficient ρ6,7,9 is defined as:

$${\rho }_{6,7,9}=\frac{{{{\rm{cov}}}}(x,y)}{{\sigma }_{x}{\sigma }_{y}}=\frac{\langle xy\rangle -\langle x\rangle \langle y\rangle }{\sqrt{\langle {x}^{2}\rangle -{\langle x\rangle }^{2}}\sqrt{\langle {y}^{2}\rangle -{\langle y\rangle }^{2}}}\quad x\equiv {P}_{6,7}^{2}(i)\quad y\equiv {S}_{7,9}^{(0)}(i)$$
(13)

where:

$$\langle x\rangle =\frac{\mathop{\sum }_{i = 1}^{N}x(i)}{N}\quad i=1,2,3,...N$$
(14)

and i indicates the ith single-shot photoelectron spectra and N the total number of single-shot photoelectron spectra used in the analysis. A similar definition was used for 〈y〉.

The correlation coefficient ρ7,9,10 is defined in a similar way by replacing \(x={P}_{6,7}^{2}\) with \(x={P}_{9,10}^{2}\).