Experimental methods

The laser system used in this experiment is the ATTO-PHAROS at KSU, which is based on an industrial-grade CEP-stabilized Yb laser system (PHAROS, Light Conversion). This system delivers 2-mJ pulses with a pulse duration of 170 fs. The details of the CEP-stabilized pulse generation are as follows30: the oscillator operates at a repetition rate of approximately 76 MHz, with the carrier-envelope-offset frequency fCEO locked to frep/4, while the laser repetition rate frep remains unlocked. A Pockels cell reduces frep by an integer multiple of four, lowering the repetition rate to 10 kHz while ensuring that each amplified pulse maintains a constant CEP. We post-compressed the laser pulses to a duration of 3.7 fs (full width at half maximum of the intensity envelope), as characterized by transient-grating FROG and TIPTOE measurements, using cascaded free-space four-pass gas cells36. After post-compression, the pulse energy is 1.4 mJ, and the spectrum spans from 600 to 1250 nm within 25 dB. The single-shot CEP stability is better than 200 mrad (standard deviation) over 24 h. See refs. 42,43 for details of the laser pulse characterization and CEP performance. Note that the IR pulse durations retrieved from the streaking measurements are approximately 6 fs, longer than the values obtained from FROG measurements performed without significant material dispersion. This discrepancy may arise from plasma-induced dispersion of the dressing field in our in-line geometry.

To optimize the pulse dispersion and vary the CEP, we adjusted the insertion depth of a fused-silica wedge pair (3.5∘) using a piezoelectric translation stage (Newport CONEX-SAG-LS32P). As illustrated in Fig. 1 of the main text, three key optical elements were used to construct the stable in-line interferometric beamline. A perforated fused-silica plate (1-mm thick with a 2-mm central hole) serves as the beam splitter. A pair of concentric flat silver mirrors was used to control the time delay between the annular driving IR beam and the inner dressing IR beam, where the diameter of the inner mirror matches the beam size of the inner beam. Before the two beams entered the semi-infinite gas cell, an iris was used to optimize the photon flux and bandwidth of the HHG beam. The HHG gas cell was filled with pure helium at a pressure of 2.0–2.5 bar, and the two beams were focused using a concave mirror with a focal length of 25 cm. The inner dressing pulse arrives earlier than the annular driving pulse but is too weak to generate HHG. After the HHG gas cell, a 1.0-mm pinhole was used to block the residual driving IR beam. Approximately 1 m downstream, a 2.0-mm-diameter metal or carbon filter with variable thickness was mounted on a perforated fused-silica plate identical to the beam splitter. Owing to its larger divergence angle than the soft-X-ray (SXR) HHG beam, the inner dressing IR beam passes through the outer region of the fused-silica plate. Consequently, the temporal delay of the dressing IR pulse relative to the SXR pulse (generated by the annular driving beam), introduced by the beam splitter, is compensated by the fused-silica filter plate. To actively control the SXR-IR time delay during the streaking experiments, we scanned the position of the inner mirror in the concentric mirror assembly before the HHG cell.

The attosecond SXR pulse and the dressing IR field are then focused onto a helium gas target inside a thick-lens high-energy VMI spectrometer45 using a Ni-coated toroidal mirror with a focal length of 500 mm. Our “thick-lens” VMI design is fundamentally different from the Einzel-lens geometry widely used in time-of-flight (TOF) spectrometers. The magnetic bottle and Einzel lens increase the electron count by integrating over the emission angle. However, angle-integrated detection introduces severe artifacts in attosecond streaking measurements and complicates accurate pulse retrieval. In our thick-lens VMI, the term “thick lens” refers to a spectrometer design capable of resolving the full angular distribution even for high-energy photoelectrons. For the VMI voltage settings, the repeller voltage was set to 4–5 kV, while the extractor voltage was maintained at 95% of the repeller voltage. The atomic beam in the VMI was generated by an effusive jet without skimmers. The VMI chamber pressure was better than 1 × 10−8 mbar without gas loading and increased to 5 × 10−6 mbar during operation with the gas load. We operated the VMI in the conventional integration mode without gating the photoelectron time of flight. The exposure time of the VMI CCD camera was set to 5 s per frame, and 12 frames were averaged at each delay position to obtain one image. One delay scan required approximately 2 h, and 10–20 scans were acquired for each filter. No Abel inversion was applied to the VMI images to avoid reconstruction artifacts. For the energy spectra shown in Figs. 5 and 6 of the main text, the transverse momentum Px was fixed at zero.

A home-made HHG spectrometer with a flat-field grating (2400/mm, SHIMADZU 30-003) was used to measure the spectrum of our SXR light pulses, where we used both 10-cm-long MCP and phosphor screen to collect the full spectrum without moving their positions. As shown in main-text Fig. 2, our spectrometer can record the SXR spectrum and its second-order diffraction at the same time. The energy calibration was performed by comparison with the transmission curves of Al, Zr, Sn and C filters. Note that the position of the second-order diffraction provides an independent check of the energy calibration.

For the calibration of the HHG photon flux, we installed a calibrated XUV photodiode (AXUV100G) at the position immediately before the entrance slit of the HHG spectrometer. We measured the photocurrent using four different filters (Sn, Zr, Al, and C) separately. With two 200-nm carbon filters, we measured a photocurrent of 700 nA using a Keithley ammeter, corresponding to an average power of 2.8 μW based on the calibrated responsivity (0.25 A/W for wavelengths below 50 nm). The photon flux at a central photon energy of approximately 150 eV is therefore estimated to be 1.2 × 1011 photons per second. After accounting for the transmission losses of the two carbon filters and the reflectivity of the toroidal mirror, the photon flux at the HHG generation stage exceeds 1 × 1012 photons per second, corresponding to approximately 1 × 108 photons per laser shot. The pulse energy, obtained by dividing the average power by the laser repetition rate (10 kHz), is 2.8 nJ at the generation stage. Using a toroidal mirror with a relatively long focal length (f = 0.5 m), together with a pulse duration of approximately 20 as and an estimated focal spot size of about 100 μm, we estimate the peak SXR intensity at the focus to exceed 1012 W/cm2. Because the ultraviolet transmission of carbon filters is not well characterized, we performed an independent cross-check using Zr and Sn filters. With two 200-nm Zr filters and two 200-nm Sn filters, we measured photocurrents of 450 nA and 370 nA, respectively. These filters effectively suppress low-order harmonics because of their negligible ultraviolet transmission. The photon fluxes derived from the Sn and Zr measurements are in good agreement with the estimate obtained using the carbon filters, as well as with the flux measured by the HHG spectrometer.

Simulation of phase-matched HHG cut-off energy as a function of pulse duration

The phase mismatch Δkq in HHG can be approximated by the simplified expression:

$$\Delta {k}_{q} \sim P\cdot q\left\{(1-\eta )\cdot \delta {n}_{q}\cdot \frac{2\pi }{{\lambda }_{L}}-\eta \cdot {N}_{{{\rm{atm}}}}\cdot {r}_{e}\cdot {\lambda }_{L}\right\},$$

(1)

where P is the gas pressure, q is the harmonic order, η is the ionization fraction, and λL is the laser wavelength57,58. The term δnq = n(λL) − n(λL/q) ≈ n(λL) − 1 represents the difference in refractive index between the fundamental and harmonic wavelengths at a pressure of 1 atm. Since the refractive index at the harmonic wavelength λL/q approaches unity in our regime of extreme ultraviolet or x-ray, δnq is effectively determined by the dispersion at the fundamental wavelength. Natm is the atomic number density at 1 atm, and re is the classical electron radius. The first term describes positive dispersion from neutral atoms, while the second accounts for the negative dispersion due to free electrons. Efficient phase matching occurs when these contributions are balanced at an appropriate ionization fraction.

Following Eq. (1), a fundamental constraint is imposed by the critical ionization level ηc, which sets the maximum ionization fraction below which phase matching is possible:

$${\eta }_{c}={\left\{1+\frac{{\lambda }_{L}^{2}{r}_{e}{N}_{{{\rm{atm}}}}}{2\pi \delta {n}_{q}}\right\}}^{-1}.$$

(2)

When the ionization fraction η exceeds ηc, plasma dispersion dominates, increasing the phase velocity of the driving field beyond that of the harmonics and preventing phase matching. This imposes a phase-matching cutoff on the maximum achievable harmonic order, which typically scales as \({I}_{L}{\lambda }_{L}^{1.6-1.7}\). This is in contrast to the single-atom cutoff, which follows the well-known scaling \({I}_{L}{\lambda }_{L}^{2}\).

The ionization fraction η is highly sensitive to the duration of the driving laser pulse, as ionization is a cumulative process. Longer pulses increase the exposure time of atoms to the laser field, resulting in greater electron ionization. Consequently, the ionization fraction can exceed the critical level ηc more rapidly, disrupting phase matching. In contrast, shorter pulses suppress ionization growth, enabling phase matching at higher intensities and thereby increasing the phase-matched cutoff energy.

To evaluate the phase matching cutoff as a function of pulse duration, we calculate the time-dependent ionization fraction η(t) using the ADK model59 and identify the moment when η = ηc. The electric field at that instant determines the maximum return kinetic energy via the classical three-step model, from which the phase-matched cutoff energy is estimated. Fig. 1b in the main text presents the resulting cutoff energies for different gases (helium, neon, argon) across a range of pulse durations. In all cases, shorter pulses yield higher phase-matched cutoffs by delaying the onset of ηc, allowing higher peak fields before phase matching is lost. For example, using a Yb-based post-compressed pulse down to 4 fs in this work, the phase-matched cutoff energy from argon can extend to around 100 eV, while that from neon reaches approximately 220 eV. In helium, the effect is most pronounced, enabling phase-matched harmonics up to 300 eV—well beyond the carbon K-edge (284 eV). These results are in good agreement with experimental observations across all three gases.

According to our simulation shown in main-text Fig. 1b, the phase-matched cut-off energy from helium, neon, and argon can reach 300, 220, and 100 eV, respectively. In the main text, we illustrate the results for helium at 2.0–2.5 bar with ~0.8 mJ single-pulse energy. In  Supplementary Information, we present the CEP-resolved HHG spectra from neon (at 1 bar with ~0.4 mJ single-pulse energy) and argon (at 170 mbar with ~0.2 mJ single-pulse energy) in the Supplementary Figs. 1 and 2, respectively. In each case, the single-pulse energy and gas pressure were iteratively optimized for broad bandwidth. For each gas species, the HHG spectra were measured with the four different filters and presented together with the corresponding filter transmission curve.

Volkov transform quasi-Newton algorithm

In this section, we first formulate pulse reconstruction as a nonlinear optimization problem and then briefly introduce the existing Volkov Transform Generalized Projection Algorithm (VTGPA). By performing an additional QR decomposition, we show that VTGPA is equivalent to a steepest descent algorithm with a fixed step size, whose efficiency is limited by its linear convergence rate. We then replace the steepest descent algorithm with a quasi-Newton optimization scheme that exhibits superlinear convergence, leading to the proposed Volkov Transform quasi-Newton Algorithm (VTQNA). Finally, a set of artificial streaking data is generated to compare the computational efficiency of VTGPA and VTQNA.

For a combination of a weak IR field and an XUV field, the SFA gives the ionization amplitude at momentum p and time delay τ as

$$M(p,\tau )=\int\,d[p+{A}_{{{\rm{IR}}}}(t)]{{{\rm{e}}}}^{{{\rm{i}}}\phi (t;p)}{E}_{{{\rm{XUV}}}}(t-\tau )\,{{\rm{dt}}},$$

(3)

with Volkov phase

$$\phi (t;p)\equiv \int^{t}\left\{\frac{{[p+{A}_{{{\rm{IR}}}}(s)]}^{2}}{2}+{I}_{p}\right\}\,{{\rm{ds}}}.$$

(4)

The ionization potential Ip = 0.9037 a.u. (for helium) and transition amplitude d(p) is taken from ref. 60. Experimentally, only the amplitudes bi ∝ ∣M(pi, τi)∣ over a finite set of points \({\{({p}_{i},{\tau }_{i})\}}_{i=1}^{m}\) are known. To be precise, we have 500 different momentum points and 121 different delays, with m = 500 × 121 = 60500. Our goal is to reconstruct the full information about EXUV(t) and AIR(t) from the measurement. To do so, we first fix the IR field and expand EXUV(t) into a finite basis

$${E}_{{{\rm{XUV}}}}(t)={\sum }_{j=1}^{n}{x}_{j}{f}_{j}(t),$$

(5)

where {xi} are complex coefficients and {fi(t)} are Hermite-Gaussian functions up to n = 96. Inserting Eq. (5) into Eq. (3), the estimated ionization amplitude takes the form

$$M({p}_{i},{\tau }_{i})={\sum }_{j=1}^{n}{W}_{ij}{x}_{j},$$

(6)

in which

$${W}_{ij}\equiv \int\,d[{p}_{i}+{A}_{{{\rm{IR}}}}(t)]{{{\rm{e}}}}^{{{\rm{i}}}\phi (t;{p}_{i})}{f}_{j}(t-{\tau }_{i})\,{{\rm{dt}}}.$$

(7)

Eq. (6) is referred to as the Volkov transformation40. Our aim thus is to minimize the target function

$${{\min }_{{{\boldsymbol{x}}}}} \; G({{\boldsymbol{x}}})={{\min }_{{{\boldsymbol{x}}}}}\frac{1}{2}{\sum }_{i=1}^{m}{\left(\left| {\sum }_{j=1}^{n}{W}_{ij}{x}_{j}\right| -{b}_{i}\right)}^{2},$$

(8)

to obtain the best estimate of EXUV(t). For convenience, the measurement data {bi} are normalized such that \(\frac{1}{2}{\sum }_{i=1}^{m}{b}_{i}^{2}=1\). The value of minimized G(x) is used as the optimization merit value in the main text.

If the phase information of

$${z}_{i}\equiv \mathop{\sum }_{j=1}^{n}{W}_{ij}{x}_{j}$$

(9)

is known, the problem becomes a standard linear least-squares problem61. Any full-rank matrix W can be decomposed into W = QR with an m × n unitary matrix Q and an n × n upper triangular matrix R. The formal solution is given by

$${{\boldsymbol{x}}}={{{\bf{R}}}}^{-1}{{{\bf{Q}}}}^{{\dagger} }({{\boldsymbol{b}}}{{{\rm{e}}}}^{-{{\rm{i}}}\arg {{\boldsymbol{z}}}}).$$

(10)

Although we do not know \(\arg {{\boldsymbol{z}}}\) in the beginning, Eq. (9) and Eq. (10) provide an iterative procedure for solving Eq. (8):

  1. (1)

    Assume \(\arg {{\boldsymbol{z}}}=0\) initially;

  2. (2)

    Compute \({{\boldsymbol{y}}}={{{\bf{Q}}}}^{{\dagger} }({{\boldsymbol{b}}}{{{\rm{e}}}}^{{{\rm{i}}}\arg {{\boldsymbol{z}}}})\), where we introduce y ≡ Rx to simply the iteration;

  3. (3)

    Obtain a new estimation of z = Qy and its phase;

  4. (4)

    Check whether the change in z is below the iteration tolerance. If not, return to step (2);

  5. (5)

    Compute x = R−1y after convergence.

This idea originates from phase retrieval algorithms for Fourier62,63 and short-time Fourier64 transformations. Now it is referred to as the generalized projection algorithm (GPA) in the literature65, and has been applied to the reconstruction of XUV pulses in attosecond streaking experiments40,66, which is referred to as VTGPA.

Actually, Equation (8) can also be viewed as a general nonlinear optimization problem. In this case, we also replace x with y ≡ Rx

$${{\min }_{{{\boldsymbol{y}}}}} \; H({{\boldsymbol{y}}})={{\min }_{{{\boldsymbol{y}}}}}\frac{1}{2}{\sum }_{i=1}^{m}{\left(\left| {\sum }_{j=1}^{n}{Q}_{ij}{y}_{j}\right| -{b}_{i}\right)}^{2},$$

(11)

and compute its derivative

$$\begin{array}{lll}{\partial }_{{y}_{k}^{*}}H({{\boldsymbol{y}}})&=&{\sum }_{i=1}^{m}\left(\left\vert {z}_{i}\right\vert -{b}_{i}\right)\frac{{z}_{i}{Q}_{ik}^{*}}{| {z}_{i}| }\\ &=&{y}_{k}-{\sum }_{i=1}^{m}{b}_{i}\frac{{z}_{i}}{| {z}_{i}| }{Q}_{ik}^{*}\end{array}$$

(12)

The simplest nonlinear optimization algorithm is the steepest descent method. It gives y at the next iterative step with step size α

$$\begin{array}{lll}{{{\boldsymbol{y}}}}^{{\prime} }&=&{{\boldsymbol{y}}}-\alpha {\partial }_{{y}_{k}^{*}}H({{\boldsymbol{y}}})\\ &=&(1-\alpha ){{\boldsymbol{y}}}+\alpha {{{\bf{Q}}}}^{{\dagger} }({{\boldsymbol{b}}}{{{\rm{e}}}}^{{{\rm{i}}}\arg {{\boldsymbol{z}}}}).\end{array}$$

(13)

One finds that, by choosing α = 1, the steepest descent method yields the same iterative procedure as in the VTGPA. The residue ϵ of the steepest descent method converges linearly, i.e. \({{\lim }_{k\to \infty }}| {\epsilon }_{k}/{\epsilon }_{k-1}|=(\kappa -1)/(\kappa+1)\), with κ being the condition number of the problem, which is generally very high for such a high-dimensional problem with m = 60500. This linear convergence leads to a slow exponential decrease of the residue. However, other optimization algorithms with superlinear convergence \({{\lim }_{k\to \infty }}| {\epsilon }_{k}/{\epsilon }_{k-1}|=0\) are generally more numerically efficient than the steepest descent method. In this work, we adopt the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method together with Wolfe-Powell inexact line search67. The computational cost of each iteration step is similar for the two methods, but the number of iteration steps to reach the convergence for a given IR field is reduced dramatically.

With a good estimate of XUV field, one can then update the IR field to minimize the target function given in Eq. (8). Again we expand AIR(t) into basis functions

$${A}_{{{\rm{IR}}}}(t)={\sum }_{i=1}^{k}{u}_{i}{g}_{i}(t),$$

(14)

here {gi(t)} are taken to be trigonometric functions up to k = 41 and the coefficients {ui} are assumed to be real. We adopt the same BFGS quasi-Newton algorithm for the optimization, which requires the first-order derivatives

$${\partial }_{{u}_{i}}G={\sum }_{j=1}^{m}\left(\left| {z}_{j}\right| -{b}_{j}\right){{\rm{Re}}}\,{{{\rm{e}}}}^{-{{\rm{i}}}\arg {z}_{j}}{\partial }_{{u}_{i}}{z}_{j},$$

(15)

with

$$\begin{array}{lll}{\partial }_{{u}_{i}}{z}_{j}&=&\int\,\left\{\frac{{d}^{{\prime} }[{p}_{j}+{A}_{{{\rm{IR}}}}(t)]}{d[{p}_{j}+{A}_{{{\rm{IR}}}}(t)]}{f}_{i}(t)+{{\rm{i}}}\int^{t}[{p}_{j}+{A}_{{{\rm{IR}}}}(s)]{f}_{i}(s)\,{{\rm{ds}}}\right\}\\ &&\times d[{p}_{j}+{A}_{{{\rm{IR}}}}(t)]{{{\rm{e}}}}^{{{\rm{i}}}\phi (t;{p}_{{{\rm{j}}}})}{E}_{{{\rm{XUV}}}}(t-{\tau }_{j})\,{{\rm{dt}}}.\end{array}$$

(16)

Although the computational cost of evaluating the above equations is relatively high, only tens of iterations are needed to obtain a converged IR field because its behavior is quite simple.

By alternating between the two steps—optimizing the IR and XUV fields—it is relatively easy to reconstruct both the XUV and IR fields within 20 minutes on a laptop. With the GPA-type algorithm, the same task takes hours to days to find a barely acceptable result.

To test the accuracy and robustness of our reconstruction method, we generated artificial streaking traces using the strong-field approximation (SFA). Poisson-distributed noise was added to mimic the statistical fluctuations in the photoelectron counts, as shown in Supplementary Fig. 4a. The maximum number of events in each pixel was set to 1000. Both the input XUV and IR pulses were assumed to have Gaussian temporal envelopes with linear chirp, as shown by the blue solid curves in Supplementary Figs. 4c-g. The retrieved temporal and spectral profiles agree very well with the input pulse parameters.

We first considered the case in which the IR field was assumed to be fully known and iteratively retrieved only the XUV field using both VTGPA and VTQNA, starting from the same random initial guess. This configuration allows a direct comparison of the convergence behavior and reconstruction accuracy of the two algorithms under identical conditions. In Supplementary Fig. 5, we further evaluate the performance of VTQNA for the more challenging case of double attosecond pulses. The satellite pulse is separated from the main pulse by half an optical cycle of the dressing IR field while remaining spectrally overlapped with it. The retrieved temporal and spectral profiles agree well with the input pulse parameters, demonstrating that VTQNA can reliably reconstruct both the main and satellite pulses. Because of the spectral overlap between the two pulses, interference fringes appear in the spectrum. The oscillatory feature associated with the satellite pulse, marked by the red arrows in Supplementary Fig. 5a, b, appears in the opposite direction from that of the main pulse. Note that both the satellite and main attosecond pulses are positively chirped in this simulation.

The convergence histories of the two algorithms are shown in Fig. 7. For both methods, the residual decreases rapidly within the first few hundred iterations. However, as the number of iterations increases, the residual of VTGPA approaches a straight line in the logarithmic plot, corresponding to an exponential decrease of \({\epsilon }_{k} \sim \exp (-k/1{0}^{5})\) as expected, while the VTQNA continues to exhibit superlinear convergence. We then remove the prior knowledge of the IR field, retain an initial guess derived from the center-of-mass momentum of the streaking trace at each time delay, and apply the full VTQNA procedure. Our retrieval algorithm converges within 5 minutes for this test case, resulting in a perfect match between input and retrieved pulses.

In the main text, we also present a histogram of the FWHM values of the retrieved pulses. These values are obtained following the procedure described in ref. 21: (i) first, we alternately optimize the XUV and IR pulses using the full trace; (ii) we then fix the IR pulse and, at each XUV-IR delay, determine the optimal XUV pulse and its corresponding FWHM; and (iii) after scanning all 121 delays, we obtain a distribution of estimated FWHM values, which is used to estimate the uncertainty.

Auto-correlation representation of attochirp

In Supplementary Fig. 3, we present both the measured and retrieved streaking traces obtained without a filter and with 200-, 400-, 500-, and 600-nm carbon filters, respectively. It is worth noting that, in addition to the sub-cycle asymmetry observed in the streaking traces, some of the present authors previously introduced an autocorrelation representation to provide a more intuitive assessment of the attochirp54. Here, we present the autocorrelation and streaking representations side by side for comparison.

The autocorrelation is defined as Q(τ1, τ2) = ∫ S(E, τ1)S(E, τ2), dE, where S(E, τ) is the retrieved streaking trace as a function of the photoelectron energy E and the SXR-IR delay τ. In the main text, we demonstrate that the attochirp can be visualized through the left-right asymmetry of the streaking trace within one optical cycle. Based on this feature, the attochirp is clearly overcompensated in panel E1 (600-nm carbon filter) of Supplementary Fig. 3. Here, we show that both the attochirp and its overcompensation can also be intuitively visualized in the autocorrelation representation. A square diamond-shaped pattern corresponds to a chirp-free pulse. As the linear chirp increases, the square pattern gradually distorts and skews along one diagonal. A linear chirp with the opposite sign causes the pattern to skew along the opposite diagonal. Higher-order chirps generally smooth the square pattern, making it more rounded. For the result obtained without a filter (panel A3), the autocorrelation pattern appears more rounded, indicating a more complex chirp. With the 200- and 400-nm carbon filters, the autocorrelation pattern expands toward the upper-right corner, as indicated by the arrows in panels C2 and C3. In contrast, with the 600-nm carbon filter, the pattern expands toward the bottom-left corner. For our systematic experimental data, both the original streaking representation and the autocorrelation representation clearly reveal the attochirp and its progressive compensation with increasing filter thickness.