Subcycle videography of Landau–Zener–Majorana transitions in graphene

Experimental setup

The subcycle momentum microscopy setup (Extended Data Fig. 1a) is based on two ytterbium-doped potassium gadolinium tungstate (Yb:KGW) amplifiers (CARBIDE, Light Conversion Inc.). Each amplifier generates near-infrared light pulses with a duration of 250 fs, a centre wavelength of 1,030 nm and a pulse energy of 1.5 mJ, at a repetition rate of 50 kHz. In the probe arm, the second harmonic of the output of a four-stage noncollinear optical-parametric amplifier (NOPA; ORPHEUS-OPCPA, Light Conversion Inc.)⁴⁸ yields carrier-envelope phase (CEP)-stabilized pulses centred at 400 nm with a pulse duration of 20 fs. These pulses are focused into a high-pressure argon gas jet to generate high-order harmonics (Extended Data Fig. 1b) in a vacuum beamline⁴⁹. The seventh harmonic at a photon energy of 21.7 eV is spectrally selected by two multilayer mirrors (Extended Data Fig. 1b, grey shaded area). They feature a bandwidth of 0.9 nm around 57 nm and thus support a Fourier-limited pulse duration of less than 10 fs. The XUV pulses enter the experimental ultrahigh-vacuum (UHV) chamber, the momentum microscope³², through a 200-nm-thin aluminium foil, which also filters out the ultraviolet driving pulses.

Because the count rate in a time-of-flight photoemission momentum microscope is limited to only about one electron per laser shot, a high repetition rate is key for keeping measurement times at a reasonable amount while ensuring good statistics. In our experiment, this requirement entails the challenge of pushing the generation of the atomically strong, few-cycle phase-stable MIR pump pulses to new limits for a tabletop setup. Here the generated phase-stable MIR pulses feature, to the best of our knowledge, record-high average powers of up to 1 W and electrical peak field strengths of more than 200 MV cm⁻¹. These field transients are obtained by difference frequency generation in a gallium selenide crystal between the outputs of two parallel (N)OPAs (ORPHEUS-MIR, Light Conversion Inc.)⁵⁰. They feature pulse durations down to sub-50 fs and can be tuned between 4 µm and 15 µm. A small portion of the laser fundamental is picked off to pump a home-built single-stage NOPA. It is optimized to generate 10-fs pulses centred around 825 nm, which are used for electro-optic sampling (EOS) of the MIR waveforms. Extended Data Fig. 1c shows a typical MIR transient detected by EOS. A pair of wire-grid polarizers allows us to precisely control the polarization direction and power of the MIR pulses, which are then coupled into the UHV chamber through a diamond window.

The intense MIR pump (s-polarized) and the XUV probe (p-polarized) are spatially overlapped in the vacuum chamber with a custom-built off-axis parabolic mirror featuring a small aperture for XUV transmission and are focused onto monolayer graphene with an angle of incidence of 70°. A piezo-controlled delay stage in the excitation beam path allows us to set the delay between the MIR and XUV pulses with a temporal precision of 230 as. The emitted photoelectrons are detected by a time-of-flight momentum microscope (lens system: GST mbH; detector: HEX40, RoentDek)⁵¹. Multihit detection ensures linearity of the signal up to a repetition rate of 20 kHz. The combination of the 21.7-eV probe pulses with the large acceptance angle of the momentum microscope allows access to momentum regions up to ±2.14 Å⁻¹, as compared with ±0.16 Å⁻¹ in previous subcycle ARPES experiments^8,20. The energy resolution is limited by the broad spectrum of the probe pulse. The bandwidth of the XUV pulse is tunable, setting a limit for the energy resolution of 180 meV. The momentum resolution is better than 0.01 Å⁻¹. The custom-designed UHV chamber is shielded by a μ-metal liner and maintains a base pressure of 2 × 10⁻¹⁰ mbar. The samples were measured at room temperature.

Inducing field-driven dynamics in the monolayer requires a field that is parallel to the plane of the monolayer. To maximize the corresponding field component, the MIR pulse is s-polarized. The polarization of the XUV pulse is set to be perpendicular to the MIR pulse to prevent crosstalk between their fields that could introduce coherent artefacts effects in the observed dynamics. Finally, the rotation of the sample is chosen such that the MIR-induced momentum streaking in vacuum is minimized at one of the K points. This requires one K point to be at k_∥ = 0, with k_∥ being the momentum direction parallel to the electric field of the MIR pulse. This alignment also minimizes matrix-element-induced asymmetries along k_∥ = 0 in the observed energy and momentum range. Although an asymmetry along the k_x direction might persist, it is rather small at the observed energies. This is because the dark corridor only begins to set in at relatively high excitation energies, which do not occur along the k_x direction in the experiment. This may also explain the slightly larger broadening observed in the simulations in which the matrix element was neglected.

Graphene samples

For sample preparation, an n-doped 6H-SiC(0001) wafer was used as a substrate for growing an epitaxial monolayer graphene with 0° rotational alignment relative to the substrate lattice (0°-EMLG). The wafer was initially degassed in UHV at 880 °C for 30 min. The temperature was subsequently increased to 1,050 °C and maintained for 30 min under an external flux of Si atoms, resulting in a Si-rich (√3 × √3)R30° surface reconstruction. The sample was then further annealed at 880 °C for an extra 30 min, again under Si flux, resulting in the Si-rich (3 × 3) surface reconstruction. With this procedure, the oxide layer is removed from the SiC surface and atomic cleanness is achieved. In the next step, 0°-EMLG was obtained by annealing the sample at 1,400 °C for 30 min in borazine (B₃N₃H₆) atmosphere with a partial pressure of 1.5 × 10⁻⁶mbar. The borazine molecules act as surfactants on the heated SiC surface. They inhibit the formation of the conventional 30°-rotated graphene layer but instead force the homogeneous graphene layer in an unconventional 0° orientation with respect to the substrate lattice. For more details on the preparation of the 0°-EMLG, see refs. ^38,52,53. Following growth, the surface morphology and crystalline structure were characterized using low-energy electron microscopy and low-energy electron diffraction. The average width of terraces on the SiC surface is approximately 160 nm and the surface is fully covered by 0°-EMLG. Subsequently, the sample was shortly exposed to air, transported in static vacuum (0.1 mbar) to the experimental chambers and degassed at 350 °C for 3 h before the measurements. The photoemission experiments reveal that the 0°-EMLG sample is n-doped, with the Fermi level approximately 200 meV above the Dirac point.

Electric field waveform reconstruction and streaking compensation

The acceleration of photoemitted electrons in vacuum allows for an in situ reconstruction of electric field waveforms²⁰. In addition to acceleration within the sample, the s-polarized field also accelerates electrons after they have left the material. This causes a momentum change of the entire band structure parallel to the surface of the sample, Δk_∥, that follows the vector potential of the field. The electric field waveform can thus be calculated by:

To extract Δk_∥, the extracted photoemission data are shifted stepwise along k_∥ and compared with a reference spectrum without a field. Within a fixed momentum window centred on the Dirac point, spanning +160 meV to −1.5 eV relative to the Fermi level (a range that, at equilibrium, encompasses predominantly the valence band), the sum of absolute pixel-wise differences between the measured data and the reference is computed at each iteration. The iteration with the smallest overall error gives the value for Δk_∥. Repeating this procedure for every delay point, t, allows a reconstruction of the applied field E_∥(t). To refine the retrieved waveform, a second method for extracting Δk_∥(t) is used. Tracking the change in the number of electrons, integrated over all detected energies (+1 eV to −7 eV), with positive or negative momentum along k_∥, yields a unitless quantity directly proportional to Δk_∥(t), which is not subject to discretization caused by the binning of the data along k_∥. The waveform that can be extracted in this way is then scaled by the peak electric field from the first method. The resulting waveform is shown in Fig. 2c.

The extracted Δk_∥ is also used to compensate for the streaking in the photoemission spectra. As well as the trivial shift along k_∥, a shift in energy, Δε, is connected to momentum streaking. The energy of every electron in vacuum is given by:

The change in energy caused by Δk_∥ can therefore be calculated with:

Extended Data Fig. 2a shows band-structure snapshots taken without an applied field (left panel) and with a strong field applied before streaking compensation (middle panel) and after streaking compensation (right panel). The streaking-induced changes of energy and momentum were compensated for in all photoemission spectra before further analysis of the occupation of the band structure.

Keldysh parameter in graphene

A parameter that is frequently used to distinguish between the weak-field and strong-field regimes of light–matter interaction is the Keldysh parameter

It was introduced by Leonid Keldysh in 1965 as the ratio between the carrier frequency of light ω and the tunnelling rate of an electron with mass \({m}_{{\rm{e}}+{\rm{h}}}^{\ast }\) through a barrier Δ, to show that tunnelling and multiphoton ionization in strong electromagnetic fields follow essentially the same physical laws⁵⁴. This ratio scales linearly with frequency ω and the inverse of the electric field strength E₀. Tunnelling dominates for γ ≪ 1, which holds for sufficiently strong fields and/or low frequencies. Equivalently, it can be expressed as \(\gamma =\sqrt{\Delta /{2U}_{{\rm{P}}}}\) by comparing the barrier Δ with the mean kinetic energy of a free electron in an oscillating electric field of amplitude E₀, the ponderomotive potential \({2U}_{{\rm{P}}}={e}^{2}{E}_{0}^{2}/4{m}_{e+h}^{\ast }{\omega }^{2}\). Although similar considerations can be made for quasirelativistic motion of electrons in gapless graphene, it is more realistic for the present experiment to take the full Dirac cone with its distribution of barriers and effective masses into account (Extended Data Fig. 7a). In that case, we obtain a k_x-dependent Keldysh parameter²⁶

because both the gap, Δ = 2ħv_F|k_x|, and the effective mass at the gap, \({m}_{e+h}^{\ast }=\hbar |{k}_{x}|/2{v}_{{\rm{F}}}\), depend on the distance from the Dirac point. By identifying Δ with the MIR photon energy of 124 meV, we get γ = 0.075 for the peak electric fields of 1.57 MV cm⁻¹. If we could switch off intraband acceleration, the relevant transition energy from the valence band to the unperturbed conduction band would increase to two times the doping energy of 200 meV, that is, to Δ = 400 meV. Even then, γ = 0.24 would comfortably place our measurements in the strong-field regime.

Current extraction and Fourier analysis of the phase delay

The microscopic intraband current density for each time instant, j(t), is directly extracted from the individual photoemission measurements, as in Fig. 2c. After compensation of momentum streaking, the photoemission intensity is integrated in two momentum regions with k_y < 0 and k_y > 0, respectively, and normalized by the total number of carriers in the two momentum regions at equilibrium. Subtracting the time-dependent number of counts with positive and negative k_y components, we obtain a dimensionless value proportional to the microscopic current. To reveal the maximum amplitude of the current density, this procedure is performed for all carriers above the Dirac point (E − E_D ≥ 0) within the momentum region |k_x − K| < 0.15 Å⁻¹. The result is then calibrated to physical values by multiplying with the Fermi velocity and carrier density in the conduction band of graphene with a Fermi level located 200 meV above the Dirac point, which can be computed from the density of states²⁷. This yields a maximum current density of j(t) = 38.2 A cm⁻¹. This value is very high, as it is even larger than in previous studies on topological insulators^8,20, in which the scattering times were two orders of magnitude longer. Given the higher Fermi velocity and field amplitude, a higher current could be expected in the presented experiment. Yet, in a Drude-like regime, this would not be sufficient to compensate for the much higher scattering rate. This confirms that, in ballistic regimes, the current does not follow the Drude model. The same procedure is applied to the results of the simulations based on the SBEs.

To extract the relative timing between the microscopic current in the Dirac cone, j(t), and the electric field waveform, E(t), we apply a phase delay analysis. We compute the product between the Fourier transform of the field waveform, \(\mathop{E}\limits^{ \sim }(\omega )=F\{E(t)\}\), and the complex conjugate of the Fourier transform of the current, \(\mathop{j}\limits^{ \sim }(\omega )=F\{j(t)\}\):

The squared modulus of P(ω) peaks at the frequency of the pump field, whereas its phase corresponds to the timing difference (or phase delay) between the two oscillating quantities. We thus define the average phase difference as:

in which φ[P(ω)] is the complex phase of P(ω) and summations only consider frequencies at which |P(ω)| is larger than 50% of its peak value. To compute the phase delay, τ_i, we need to divide this quantity by the central frequency of the pump field. We then obtain the mean phase delay, τ_ϕ, and the associated standard deviation of the mean over repeated measurements as:

in which i = 1,…, N runs over three different measurements performed on this sample geometry.

By performing a gated phase delay analysis, it is also possible to calculate a time-dependent phase difference, ϕ(t), which indicates the evolution of the scattering processes throughout the duration of the pulse. Here, however, the delay remains relatively stable over the first 100 fs of the driving pulse. We therefore use a single average phase shift, which we associate with a characteristic scattering time. In principle, however, this method could open new avenues for investigating relaxation dynamics across different samples with unprecedented temporal resolution.

Semiconductor Bloch equations

To simulate the time-dependent MIR-driven electron dynamics in graphene, we numerically solve the SBEs^6,55,56 for a two-band system:

in which k is the crystal momentum, t is the time, ħ the reduced Planck constant and e the elementary charge. f_e(k, t) and f_h(k, t) are the electron distributions in the conduction and valence bands, respectively, and p_eh(k, t) is the microscopic interband polarization. In these equations, E(t) is the MIR electric field, T₂ = 5 fs is the phenomenological dephasing time of the microscopic polarization, ε_e(k) and ε_h(k) are the energy values of the conduction and valence bands and

are the dipole moments. \({\varGamma }_{{\bf{k}}}^{{\rm{e}}}\) and \({\varGamma }_{{\bf{k}}}^{{\rm{h}}}\) are the time-dependent and momentum-dependent scattering rates for electrons and holes, respectively, and will be discussed in the next section. In simulations without interband coupling (Fig. 4 and Extended Data Fig. 4), we set \({{\bf{d}}}_{{\rm{eh}}}({\bf{k}})={{\bf{d}}}_{{\rm{he}}}^{\dagger }({\bf{k}})=0\).

In the specific case of graphene, the Hamiltonian of a two-level system can be written as²⁷:

with

and a = −2.46 Å is the nearest-neighbour lattice distance, v_F = 1.072 × 10⁶ m s⁻¹ is the experimental value of the Fermi velocity and \(\gamma =-2\hbar {v}_{{\rm{F}}}/(\sqrt{3}{\rm{a}})\) is the hopping integral. From the Hamiltonian, it is then possible to compute the equilibrium eigenenergies and eigenfunctions:

The out-of-equilibrium dynamics shown in Extended Data Figs. 4 and 9 were simulated on a 250 × 250 2D square (k_x, k_y) grid in momentum space. The grid is centred on one of the K points with coordinates k_x,K = 1.703 Å⁻¹ and k_y,K = 0 and it has a side length of 0.73 Å⁻¹. The simulated time axis ranges from −100 fs to 300 fs, with a step size of 2 fs. The electric field, E(t), is parallel to k_y and its temporal profile is analytically fitted to the experimentally measured field. At the initial time instant (t₀ = −100 fs), we assume p_eh(k, t₀) = 0, whereas the electron distribution corresponds to a Fermi–Dirac distribution. The position of the Fermi level compared with the Dirac crossing is set to ε_F = 0.2 eV to match the experiment. Simulations comparing the dynamics at other driving frequencies, varying field strengths or different Fermi levels sometimes used analytically modelled waveforms and smaller grid sizes to keep up with increased demands for spatial or temporal resolution. 

Modelling of the scattering rates

In simulations including scattering, the grid was adjusted to 50 × 50 while maintaining the size in momentum space, to save resources. To reproduce the time-dependent electron distribution observed in momentum space, it is imperative to realistically model the scattering rates:

On these timescales, two physical phenomena mainly contribute to the scattering rates:

in which \({\varGamma }_{{\bf{k}},{\rm{e}}{\rm{p}}}^{\alpha }\) and \({\varGamma }_{{\bf{k}},{\rm{e}}{\rm{e}}}^{\alpha }\) are, respectively, the electron–phonon and electron–electron scattering rates.

Starting from the electron–phonon interaction, the scattering rate can be written as⁵⁷:

in which v is the phonon mode considered and α, α′ ∈ {e, h} are the band indexes. The first term thus describes transitions from an occupied (k′, α′) state to an empty (k, α) state owing to the interaction with a v phonon, whereas the second term describes the opposite process. Whereas in the former the state (k, α) is empty and becomes occupied (thus, the rate is positive), in the latter, the rate must be negative. The probability of these transitions is associated with the electron–phonon transition rate:

in which \({G}_{{\alpha }^{{\prime} }\alpha }^{v}({{\bf{k}}}^{{\prime} },{\bf{k}})\) is the electron–phonon coupling constant, δ is the Kronecker delta imposing energy and momentum conservation, q and \({\omega }_{{\bf{q}}}^{v}\) are the momentum and angular frequency of a v phonon, respectively, and \({n}_{{\bf{q}}}^{v}\) is the number of v phonons with a given frequency and momentum, which is described by the Bose–Einstein distribution.

To simplify our description, we assume that: (1) Umklapp processes are negligible and (2) the number of phonons \({n}_{{\bf{q}}}^{v}\) is constant and equal to its value at equilibrium (phonon bath). Moreover, previous works explained the dominating contribution of LO and TO phonons^34,47. We will thus limit our model to LO and TO phonons and, because they have similar phonon energies (\({\omega }_{{\bf{q}}}^{{\rm{LO}}}\simeq {\omega }_{{\bf{q}}}^{{\rm{TO}}}\simeq 164.6\,{\hbar }^{-1}\,{\rm{meV}}\)) and matrix elements but their angular dependencies compensate each other, we will treat them together with an effective electron–optical phonon coupling constant⁵⁷:

The effective matrix element is constant in momentum space and its value is adjusted to reproduce the experimental data. Owing to the isotropic angular dependence of this term, electron–optical phonon scattering will lead to a rapid spreading of the electron distribution around the Dirac cone.

In the same spirit, the scattering rate for electron–electron interaction can be written as⁵⁷:

in which \(\alpha ,{\alpha }^{{\prime} },\beta ,{\beta }^{{\prime} }\in \{{\rm{e}},{\rm{h}}\}\) are the band indices. The first term, therefore, describes the scattering between a pair of particles in the initial states (α′, k′) and \(({\beta }^{{\prime} },{{\bf{k}}}_{\ast }^{{\prime} })\) to the final states (α, k) and \((\beta ,{{\bf{k}}}_{\ast })\), and the second term refers to the opposite process. The associated transition rate can thus be written as⁵⁷:

in which the matrix element is:

and q = |k − k′| and \({q}^{{\prime} }=|{\bf{k}}-{{\bf{k}}}_{\ast }^{{\prime} }|\). The Coulomb potential between electron pairs is:

with A = 4π²/Δk². The screening of the Coulomb interaction is given by the dielectric function \({\epsilon }(q)=1+{v}_{{\rm{c}}}(q)D({\varepsilon }_{{\rm{F}}})\mathop{\Pi }\limits^{ \sim }(q)\)^57,58, in which:

and r_S is the Wigner–Seitz radius of graphene, k_F = ε_F/(ħv_F) the Fermi wavevector and D(ε) the density of states close to the Dirac point. This approximation holds for n ≥ n_th = 10¹² cm⁻², thus if ε_F ≥ 37.14 meV, as in our case.

In this form, the scattering rate for electron–electron interaction is computationally extremely expensive. Therefore, we assume that the electron distribution for the second particle can be approximated by the equilibrium Fermi–Dirac distribution of the carriers:

This allows us to write a simplified expression for the electron–electron scattering rate:

in which we use that \({S}_{{\rm{ee}}}({\bf{k}},{{\bf{k}}}_{\ast },{{\bf{k}}}^{{\prime} },{{\bf{k}}}_{\ast }^{{\prime} })={S}_{{\rm{ee}}}({{\bf{k}}}^{{\prime} },{{\bf{k}}}_{\ast }^{{\prime} },{\bf{k}},{{\bf{k}}}_{\ast })\) and define:

In our simulations, we only focus on a small momentum space region around one of the six equivalent K points. In doing so, we neglect intervalley electron–phonon and electron–electron scattering processes, which are not expected to contribute notably in these conditions³⁴. However, for electron–electron scattering, it is possible that an electron scatters from one state to another within the same valley by interacting with a charge carrier in a different one. For this reason, we effectively account for these processes by multiplying \({S}_{{\rm{ee}},\beta \beta {\rm{{\prime} }}}^{{\rm{eff}},1}\) and \({S}_{{\rm{ee}},\beta \beta {\rm{{\prime} }}}^{{\rm{eff}},2}\) by the number of equivalent K valleys, that is, six. To reproduce the experimental data, we assumed V_eff = 4.386. It is important to notice that the dependence on the relative angle between k, k′, \({{\bf{k}}}_{\ast }\) and \({{\bf{k}}}_{\ast }^{{\prime} }\) determines the characteristic angular dependence of electron–electron scattering in graphene⁴⁷, also termed collinear scattering⁴¹. Notably, at variance with the relaxation time approximation, our approach does not enforce any specific electron distribution at any time instant, while matching our experimental results.

Relaxation time approximation

To reduce the computational costs when relating the phase delay to an effective scattering time, a customary approach when solving the SBEs is to invoke the relaxation time approximation^6,20,55,56. In this context, the scattering rate \({\varGamma }_{{\bf{k}}}^{\alpha }\) is assumed to be inversely proportional to a time constant, T₁, associated with the mean free time between scattering events. Defining \({f}_{\alpha }^{(0)}({\bf{k}})\) as the equilibrium electron distribution in the α band, we obtain:

It is important to notice that, even though this approach allows us to extract an effective time constant, it is not enough to reproduce the experimental data. In fact, in this approximation, the scattering rate at a given point k in momentum space will only depend on the instantaneous electron distribution at that point and not on any other point in momentum space. This makes it impossible to reproduce any carrier redistribution dynamics in momentum space owing to, for example, the interaction with optical phonons.

LZM transition region

The original works of Landau²², Zener²³, Majorana²⁴ and Stückelberg³⁶ considered a two-state quantum system with an energy separation Δε = ε₁(t) − ε₂(t) = αt that changes linearly with time and undergoes an avoided crossing with a half gap energy δ. The diabatic transition probability is then given by²⁶

with the diabaticity parameter \(\varGamma =\frac{{\delta }^{2}}{\hbar \alpha }\). This simple formula can be used to estimate the region in momentum space in which transitions occur in our experiment. The model is applicable under the condition that the characteristic sweep time through the transition, τ = 2δ/ε, is short compared with the oscillation period T = 2π/ω of the driving electric field. In this case, the field can be considered quasistatic during the transition and we simply have \(\alpha =2{v}_{{\rm{F}}}\hbar |{k}_{x}-K|=2{v}_{{\rm{F}}}|e|{E}_{0}\). With δ = v_Fħ|k_x − K|, we obtain \(\tau =\frac{\hbar | {k}_{x}-K|}{|e|{E}_{0}}=\frac{\gamma }{\omega }\). Because γ ≈ 0.1, the crossing indeed occurs very fast (τ < 1 fs). For the peak electric field E₀ = 1.57 MV cm⁻¹, the diabaticity parameter \(\varGamma =\frac{\hslash {v}_{{\rm{F}}}}{2|e|{E}_{0}}{|{k}_{x}-K|}^{2}\,\cong \) \(210\,{{\rm{\mathring{\rm A} }}}^{2}{|{k}_{x}-K|}^{2}\) becomes smaller than 1/π for |k_x − K| < 0.04 Å⁻². In this region around the Dirac point, 1 − e⁻² = 86% of the transitions from the valence band to the conductions band are expected to take place as long as Pauli blocking can be neglected. This region roughly coincides with the momentum area in which the conduction band of our doped graphene samples is occupied. This shows that transitions can only take place after the electrons have been accelerated away from the Dirac point. After several driving cycles, the observed dynamics of the doped and a hypothetical undoped sample are expected to be similar (Extended Data Fig. 8).

Injection time analysis

For an empty upper Dirac cone, we expect the carrier injection rate owing to LZM tunnelling to closely follow the magnitude of the electric field and the electron population to increase stepwise, with the highest injection rates occurring at the largest absolute field values. In an n-doped system and for high injection rates, as in the present experiment, Pauli blocking modifies this simple picture (Extended Data Fig. 5a). Electrons can only be injected when there is sufficient empty space in the conduction band near the K point. Consequently, intraband transport plays a decisive role in determining the injection timing. Carrier injection happens most efficiently during the intervals between the extrema of the electric field and the carrier displacement. Our simulations clearly show that this results in a phase delay between the population jumps and the field, which grows as more population has been transferred. Beyond the effects of Pauli blocking and scattering, the injection dynamics can be strongly influenced by quantum interference⁵⁹. Here we find that they can cause a deviation of a step-like increase in the conduction band population (or electron distribution width) by transforming the plateau regions, following a strong injection phase, into local minima through destructive interference of electron wavepackets (Extended Data Fig. 9). Although the signal-to-noise ratio of the measured conduction band population does not allow us to resolve the jumps in the population during carrier injection experimentally, the dynamics of the increase of the carrier density can be well reproduced with the theory.

The simulations also show that the elongation of the conduction band electron distribution in momentum space is closely tied to the carrier population (Extended Data Fig. 5b). Notably, the weak oscillations in the width of the electron distribution, which are superimposed on the stepwise increase, become more pronounced with scattering (Fig. 4c). This effect is probably caused by a shortening of the effective injection time window. In the absence of scattering, the electron displacement follows the vector potential and is therefore phase-shifted by π/2 relative to the electric field. With scattering, the phase shift between the electric field and the displacement of the electrons shrinks, such that periods of rapid population increase occur only a few femtoseconds before an extremum in the displacement (Extended Data Fig. 5). Experimentally, the oscillation of the width of the electron distribution at twice the frequency of the driving field is clearly visible during the first three cycles of the electric field, that is, at times when there is still substantial carrier injection (Fig. 3c). The observed relative phase with respect to the electric field transient and the COM displacement agrees well with the simulations that include scattering. The data thus support the picture of a well that is periodically driven by the electric field and opened and closed by the oscillating carriers in the conduction band.

Data processing

The new experimental platform presented here was developed to give direct access to strongly driven electron dynamics without relying heavily on theory or models. This philosophy also translates to data processing: our conclusions only rely on features present in the raw data. The previously described streaking compensation is the first and only step in the data analysis pipeline that is required to ensure that the observed dynamics can be attributed to processes that are happening in the solid and not the vacuum. Following that, the general strategy is that filtering is only used to remove technical noise and signals related to the matrix elements before the data are analysed or to enhance visibility in the figures. To ensure full transparency, the following contains a summary of all notable data processing steps for the presented figures.

The raw data from the delay line detector are binned with a spacing of 0.084 eV × 0.022 Å⁻¹ × 0.022 Å⁻¹ for the data in Fig. 1 and with a spacing of 0.084 eV × 0.011 Å⁻¹ × 0.011 Å⁻¹ for all other photoemission data.

The 3D data in Fig. 1c are symmetrized to remove the decrease in photoemission intensity caused by the angle of incidence of the XUV pulse.

The 2D momentum maps in Fig. 2a are obtained by integrating the 3D data within an energy window of 0.33 eV around distinct energy levels (0.25 eV, −0.25 eV, −0.75 eV, −1.25 eV, −1.75 eV). Before plotting, the data are further processed to enhance the visibility in the figure. First, a 1D simple moving average along both momentum axes (window size: 0.08 Å⁻¹) is applied to the data four times. The filtered data are further processed by a 2D curvature filter described by equation (11) of ref. ⁶⁰ with C₀ = 0.2. Finally, a 2D median filter with window size 0.045 eV × 0.045 Å⁻¹ is applied three times to clean up numerical noise.

The band structure cuts in Fig. 2b are obtained by integrating the 3D data within a momentum window of 0.08 Å⁻¹ around k_x − K = 0. Before plotting, the data are further processed to enhance the visibility in the figure. First, a 1D simple moving average along the energy and momentum axes with window sizes of 0.59 eV and 0.08 Å⁻¹, respectively, is applied four times. The filtered data are subject to a 2D curvature filter with C₀ = 0.2. Finally, a 2D median filter with window size 0.34 eV × 0.045 Å⁻¹ is applied three times. The raw data are shown in Extended Data Fig. 2b.

The applied field transient and current in Fig. 2c are extracted from the raw data as laid out in Methods sections ‘Electric field waveform reconstruction and streaking compensation’ and ‘Current extraction and Fourier analysis of the phase delay’. Before plotting, a 4-THz high-pass filter is applied to the data to remove slowly varying offsets. The data showing the increase of the carrier population in the conduction band were acquired under almost identical experimental conditions in which special care has been taken to align the MIR polarization such that the electric field exhibits no perpendicular component leading to coherent oscillations of the integrated photoemission intensity.

The 2D momentum maps in Fig. 3a are obtained by integrating the 3D data for energies higher than 170 meV. This energy level is chosen to minimize contributions of the valence band caused by the limited energy resolution of the setup. Before extracting the COM displacement and width, a Gaussian filter with a small window size 0.02 Å⁻¹ × 0.02 Å⁻¹ has been applied to the 2D data. We applied the same filter to the experimental 2D data for a better comparison between experiment and theory. For the theory data, this is necessary to reduce numerical or binning artefacts. Owing to its small window size, the filter does not change or remove any of the features retrieved for data analysis.

To extract the COM motion and change of the width of the carrier distribution shown in Fig. 3b,c, the plotted 2D maps are integrated along k_x and k_y. For the data along k_y, a Gaussian distribution with a constant offset (accounting for the background) is fitted to the integrated distribution. For the data along k_x, the previously constant offset is modified with a linear decrease for growing k_x values to account for the decrease in background caused by the photoemission matrix element. The mean of the extracted fit for the k_y distribution is used to characterize the COM motion. Before plotting, a 15-THz high-pass filter is applied to these data to remove a slowly varying offset. The change of the standard deviation, Δσ, of the extracted fits is used to characterize the change of the width of the carrier distribution along k_x and k_y.

The 2D momentum maps in Fig. 4a are obtained by integrating the simulated population of the conduction band and applying a Gaussian filter with a small window size of 0.02 Å⁻¹ × 0.02 Å⁻¹ to the 2D data, as described above for the data shown in Fig. 3a.

The extraction of the simulated current density and change of the width of the carrier distribution shown in Fig. 4b,c follows the procedures outlined for the experimental data.