Experimental setup<\/p>\n
The subcycle momentum microscopy setup (Extended Data Fig.\u20091a<\/a>) 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\u2009fs, a centre wavelength of 1,030\u2009nm and a pulse energy of 1.5\u2009mJ, at a repetition rate of 50\u2009kHz. In the probe arm, the second harmonic of the output of a four-stage noncollinear optical-parametric amplifier (NOPA; ORPHEUS-OPCPA, Light Conversion Inc.)48<\/a> yields carrier-envelope phase (CEP)-stabilized pulses centred at 400\u2009nm with a pulse duration of 20\u2009fs. These pulses are focused into a high-pressure argon gas jet to generate high-order harmonics (Extended Data Fig.\u20091b<\/a>) in a vacuum beamline49<\/a>. The seventh harmonic at a photon energy of 21.7\u2009eV is spectrally selected by two multilayer mirrors (Extended Data Fig.\u20091b<\/a>, grey shaded area). They feature a bandwidth of 0.9\u2009nm around 57\u2009nm and thus support a Fourier-limited pulse duration of less than 10\u2009fs. The XUV pulses enter the experimental ultrahigh-vacuum (UHV) chamber, the momentum microscope32<\/a>, through a 200-nm-thin aluminium foil, which also filters out the ultraviolet driving pulses.<\/p>\n 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\u2009W and electrical peak field strengths of more than 200\u2009MV\u2009cm\u22121. 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.)50<\/a>. They feature pulse durations down to sub-50\u2009fs and can be tuned between 4\u2009\u00b5m and 15\u2009\u00b5m. 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\u2009nm, which are used for electro-optic sampling (EOS) of the MIR waveforms. Extended Data Fig.\u20091c<\/a> 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.<\/p>\n 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\u00b0. 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\u2009as. The emitted photoelectrons are detected by a time-of-flight momentum microscope (lens system: GST mbH; detector: HEX40, RoentDek)51<\/a>. Multihit detection ensures linearity of the signal up to a repetition rate of 20\u2009kHz. 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 \u00b12.14\u2009\u00c5\u22121, as compared with \u00b10.16\u2009\u00c5\u22121 in previous subcycle ARPES experiments8<\/a>,20<\/a>. 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\u2009meV. The momentum resolution is better than 0.01\u2009\u00c5\u22121. The custom-designed UHV chamber is shielded by a \u03bc-metal liner and maintains a base pressure of 2\u2009\u00d7\u200910\u221210\u2009mbar. The samples were measured at room temperature.<\/p>\n 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\u2225\u2009=\u20090, with k\u2225 being the momentum direction parallel to the electric field of the MIR pulse. This alignment also minimizes matrix-element-induced asymmetries along k\u2225\u2009=\u20090 in the observed energy and momentum range. Although an asymmetry along the kx 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 kx direction in the experiment. This may also explain the slightly larger broadening observed in the simulations in which the matrix element was neglected.<\/p>\n Graphene samples<\/p>\n For sample preparation, an n-doped 6H-SiC(0001) wafer was used as a substrate for growing an epitaxial monolayer graphene with 0\u00b0 rotational alignment relative to the substrate lattice (0\u00b0-EMLG). The wafer was initially degassed in UHV at 880\u2009\u00b0C for 30\u2009min. The temperature was subsequently increased to 1,050\u2009\u00b0C and maintained for 30\u2009min under an external flux of Si atoms, resulting in a Si-rich (\u221a3\u2009\u00d7\u2009\u221a3)R30\u00b0 surface reconstruction. The sample was then further annealed at 880\u2009\u00b0C for an extra 30\u2009min, again under Si flux, resulting in the Si-rich (3\u2009\u00d7\u20093) surface reconstruction. With this procedure, the oxide layer is removed from the SiC surface and atomic cleanness is achieved. In the next step, 0\u00b0-EMLG was obtained by annealing the sample at 1,400\u2009\u00b0C for 30\u2009min in borazine (B3N3H6) atmosphere with a partial pressure of 1.5\u2009\u00d7\u200910\u22126\u2009mbar. The borazine molecules act as surfactants on the heated SiC surface. They inhibit the formation of the conventional 30\u00b0-rotated graphene layer but instead force the homogeneous graphene layer in an unconventional 0\u00b0 orientation with respect to the substrate lattice. For more details on the preparation of the 0\u00b0-EMLG, see refs.\u200938<\/a>,52<\/a>,53<\/a>. 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\u2009nm and the surface is fully covered by 0\u00b0-EMLG. Subsequently, the sample was shortly exposed to air, transported in static vacuum (0.1\u2009mbar) to the experimental chambers and degassed at 350\u2009\u00b0C for 3\u2009h before the measurements. The photoemission experiments reveal that the 0\u00b0-EMLG sample is n-doped, with the Fermi level approximately 200\u2009meV above the Dirac point.<\/p>\n Electric field waveform reconstruction and streaking compensation<\/p>\n The acceleration of photoemitted electrons in vacuum allows for an in situ reconstruction of electric field waveforms20<\/a>. 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, \u0394k\u2225, that follows the vector potential of the field. The electric field waveform can thus be calculated by:<\/p>\n $${{\\bf{E}}}_{\\parallel }(t)=-\\frac{\\hbar }{e}\\frac{{\\rm{d}}}{{\\rm{d}}t}\\Delta {k}_{\\parallel }(t).$$<\/p>\n To extract \u0394k\u2225, the extracted photoemission data are shifted stepwise along k\u2225 and compared with a reference spectrum without a field. Within a fixed momentum window centred on the Dirac point, spanning +160\u2009meV to \u22121.5\u2009eV 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 \u0394k\u2225. Repeating this procedure for every delay point, t, allows a reconstruction of the applied field E\u2225(t). To refine the retrieved waveform, a second method for extracting \u0394k\u2225(t) is used. Tracking the change in the number of electrons, integrated over all detected energies (+1\u2009eV to \u22127\u2009eV), with positive or negative momentum along k\u2225, yields a unitless quantity directly proportional to \u0394k\u2225(t), which is not subject to discretization caused by the binning of the data along k\u2225. 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.\u20092c<\/a>.<\/p>\n The extracted \u0394k\u2225 is also used to compensate for the streaking in the photoemission spectra. As well as the trivial shift along k\u2225, a shift in energy, \u0394\u03b5, is connected to momentum streaking. The energy of every electron in vacuum is given by:<\/p>\n $$\\varepsilon ({\\bf{k}})=\\frac{{\\hbar }^{2}}{2m}{{\\bf{k}}}^{2}.$$<\/p>\n The change in energy caused by \u0394k\u2225 can therefore be calculated with:<\/p>\n $$\\Delta \\varepsilon ({k}_{\\parallel })=\\varepsilon ({\\bf{k}}+\\Delta {\\bf{k}})-\\varepsilon ({\\bf{k}})=\\frac{{\\hbar }^{2}}{m}{k}_{\\parallel }\\Delta {k}_{\\parallel }+\\frac{{\\hbar }^{2}}{2m}\\Delta {k}_{\\parallel }^{2}$$<\/p>\n Extended Data Fig.\u20092a<\/a> 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.<\/p>\n Keldysh parameter in graphene<\/p>\n A parameter that is frequently used to distinguish between the weak-field and strong-field regimes of light\u2013matter interaction is the Keldysh parameter<\/p>\n $$\\gamma =\\frac{\\omega \\sqrt{{m}_{{\\rm{e}}+{\\rm{h}}}^{\\ast }\\Delta }}{|e|{E}_{0}}.$$<\/p>\n It was introduced by Leonid Keldysh in 1965 as the ratio between the carrier frequency of light \u03c9 and the tunnelling rate of an electron with mass \\({m}_{{\\rm{e}}+{\\rm{h}}}^{\\ast }\\) through a barrier \u0394, to show that tunnelling and multiphoton ionization in strong electromagnetic fields follow essentially the same physical laws54<\/a>. This ratio scales linearly with frequency \u03c9 and the inverse of the electric field strength E0. Tunnelling dominates for \u03b3\u2009\u226a\u20091, 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 \u0394 with the mean kinetic energy of a free electron in an oscillating electric field of amplitude E0, 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.\u20097a<\/a>). In that case, we obtain a kx-dependent Keldysh parameter26<\/a><\/p>\n $$\\gamma =\\frac{\\hbar |{k}_{x}|\\omega }{|e|{E}_{0}}=\\frac{\\omega \\Delta }{2{v}_{{\\rm{F}}}|e|{E}_{0}}$$<\/p>\n because both the gap, \u0394\u2009=\u20092\u0127vF|kx|, 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 \u0394 with the MIR photon energy of 124\u2009meV, we get \u03b3\u2009=\u20090.075 for the peak electric fields of 1.57\u2009MV\u2009cm\u22121. 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\u2009meV, that is, to \u0394\u2009=\u2009400\u2009meV. Even then, \u03b3\u2009=\u20090.24 would comfortably place our measurements in the strong-field regime.<\/p>\n Current extraction and Fourier analysis of the phase delay<\/p>\n The microscopic intraband current density for each time instant, j(t), is directly extracted from the individual photoemission measurements, as in Fig.\u20092c<\/a>. After compensation of momentum streaking, the photoemission intensity is integrated in two momentum regions with ky\u2009<\u20090 and ky\u2009>\u20090, 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 ky 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\u2009\u2212\u2009ED\u2009\u2265\u20090) within the momentum region |kx\u2009\u2212\u2009K|\u2009<\u20090.15\u2009\u00c5\u22121. 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\u2009meV above the Dirac point, which can be computed from the density of states27<\/a>. This yields a maximum current density of j(t)\u2009=\u200938.2\u2009A\u2009cm\u22121. This value is very high, as it is even larger than in previous studies on topological insulators8<\/a>,20<\/a>, 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.<\/p>\n 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)\\}\\):<\/p>\n $$P(\\omega )=\\mathop{E}\\limits^{ \\sim }(\\omega )\\times {\\mathop{j}\\limits^{ \\sim }}^{\\dagger }(\\omega ).$$<\/p>\n The squared modulus of P(\u03c9) 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:<\/p>\n $$\\phi =\\frac{\\sum _{\\omega }|P(\\omega )|\\times \\varphi [P(\\omega )]}{\\sum _{\\omega }|P(\\omega )|},$$<\/p>\n in which \u03c6[P(\u03c9)] is the complex phase of P(\u03c9) and summations only consider frequencies at which |P(\u03c9)| is larger than 50% of its peak value. To compute the phase delay, \u03c4i, we need to divide this quantity by the central frequency of the pump field. We then obtain the mean phase delay, \u03c4\u03d5, and the associated standard deviation of the mean over repeated measurements as:<\/p>\n $${\\tau }_{\\phi }=\\frac{\\sum _{i}{\\tau }_{i}}{N},\\,{\\sigma }_{{\\tau }_{\\phi }}=\\frac{\\sqrt{\\sum _{i}{[{\\tau }_{i}-{\\tau }_{\\phi }]}^{2}}}{N},$$<\/p>\n in which i\u2009=\u20091,\u2026,\u2009N runs over three\u00a0different measurements performed on this sample geometry.<\/p>\n By performing a gated phase delay analysis, it is also possible to calculate a time-dependent phase difference, \u03d5(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\u2009fs 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.<\/p>\n Semiconductor Bloch equations<\/p>\n To simulate the time-dependent MIR-driven electron dynamics in graphene, we numerically solve the SBEs6<\/a>,55<\/a>,56<\/a> for a two-band system:<\/p>\n $$\\hbar \\frac{\\partial {f}_{{\\rm{e}}}}{\\partial t}({\\bf{k}},t)=2\\hbar {\\bf{E}}(t)\\cdot {\\mathfrak{I}}\\{{{\\bf{d}}}_{{\\rm{h}}{\\rm{e}}}({\\bf{k}}){p}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}},t)\\}-e{\\bf{E}}(t)\\cdot {\\nabla }_{{\\bf{k}}}{f}_{{\\rm{e}}}({\\bf{k}},t)+{\\varGamma }_{{\\bf{k}}}^{{\\rm{e}}}$$<\/p>\n $$\\hbar \\frac{\\partial {f}_{{\\rm{h}}}}{\\partial t}({\\bf{k}},t)=-2\\hbar {\\bf{E}}(t)\\cdot {\\mathfrak{I}}\\{{{\\bf{d}}}_{{\\rm{h}}{\\rm{e}}}({\\bf{k}}){p}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}},t)\\}-e{\\bf{E}}(t)\\cdot {\\nabla }_{{\\bf{k}}}{f}_{{\\rm{h}}}({\\bf{k}},t)+{\\varGamma }_{{\\bf{k}}}^{{\\rm{h}}}$$<\/p>\n $$\\begin{array}{l}i\\hbar \\frac{\\partial {p}_{{\\rm{e}}{\\rm{h}}}}{\\partial t}({\\bf{k}},t)=\\left[{\\varepsilon }_{{\\rm{e}}}({\\bf{k}})-{\\varepsilon }_{{\\rm{h}}}({\\bf{k}})-\\frac{i\\hbar }{{T}_{2}}\\right]{p}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}},t)\\\\ \\,+\\,{\\bf{E}}(t)\\cdot {{\\bf{d}}}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}})[{f}_{{\\rm{e}}}({\\bf{k}},t)-{f}_{{\\rm{h}}}({\\bf{k}},t)]+{\\bf{E}}(t)\\cdot [{{\\bf{d}}}_{{\\rm{h}}{\\rm{h}}}({\\bf{k}})\\\\ \\,-\\,{{\\bf{d}}}_{{\\rm{e}}{\\rm{e}}}({\\bf{k}})]{p}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}},t)-ie{\\bf{E}}(t)\\cdot {\\nabla }_{{\\bf{k}}}{p}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}},t)\\end{array}$$<\/p>\n in which k is the crystal momentum, t is the time, \u0127 the reduced Planck constant and e the elementary charge. fe(k,\u2009t) and fh(k,\u2009t) are the electron distributions in the conduction and valence bands, respectively, and peh(k,\u2009t) is the microscopic interband polarization. In these equations, E(t) is the MIR electric field, T2\u2009=\u20095\u2009fs is the phenomenological dephasing time of the microscopic polarization, \u03b5e(k) and \u03b5h(k) are the energy values of the conduction and valence bands and<\/p>\n $$d({\\bf{k}})=\\left[\\begin{array}{cc}{{\\bf{d}}}_{{\\rm{e}}{\\rm{e}}}({\\bf{k}}) & {{\\bf{d}}}_{{\\rm{e}}{\\rm{h}}}({\\bf{k}})\\\\ {{\\bf{d}}}_{{\\rm{h}}{\\rm{e}}}({\\bf{k}}) & {{\\bf{d}}}_{{\\rm{h}}{\\rm{h}}}({\\bf{k}})\\end{array}\\right]=e\\left[\\begin{array}{cc}\\langle {\\varPsi }_{{\\rm{e}}}|i\\nabla |{\\varPsi }_{{\\rm{e}}}\\rangle & \\langle {\\varPsi }_{{\\rm{h}}}|i\\nabla |{\\varPsi }_{{\\rm{e}}}\\rangle \\\\ \\langle {\\varPsi }_{{\\rm{e}}}|i\\nabla |{\\varPsi }_{{\\rm{h}}}\\rangle & \\langle {\\varPsi }_{{\\rm{h}}}|i\\nabla |{\\varPsi }_{{\\rm{h}}}\\rangle \\end{array}\\right]$$<\/p>\n 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.\u20094<\/a> and Extended Data Fig.\u20094<\/a>), we set \\({{\\bf{d}}}_{{\\rm{eh}}}({\\bf{k}})={{\\bf{d}}}_{{\\rm{he}}}^{\\dagger }({\\bf{k}})=0\\).<\/p>\n In the specific case of graphene, the Hamiltonian of a two-level system can be written as27<\/a>:<\/p>\n $$H({\\bf{k}})=\\left[\\begin{array}{cc}0 & \\gamma h({\\bf{k}})\\\\ \\gamma {h}^{\\ast }({\\bf{k}}) & 0\\end{array}\\right],$$<\/p>\n with<\/p>\n $$h({\\bf{k}})=|h({\\bf{k}})|{{\\rm{e}}}^{{\\rm{i}}{\\varphi }_{{\\bf{k}}}}={{\\rm{e}}}^{{\\rm{i}}\\frac{a{k}_{x}}{\\sqrt{3}}}+2{{\\rm{e}}}^{-{\\rm{i}}\\frac{a{k}_{x}}{2\\sqrt{3}}}\\cos \\left(\\frac{a{k}_{y}}{2}\\right)$$<\/p>\n and a\u2009=\u2009\u22122.46\u2009\u00c5 is the nearest-neighbour lattice distance, vF\u2009=\u20091.072\u2009\u00d7\u2009106\u2009m\u2009s\u22121 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:<\/p>\n $$\\begin{array}{cc}{\\varepsilon }_{{\\rm{e}},{\\rm{h}}}({\\bf{k}})=\\pm \\gamma |h({\\bf{k}})|, & {\\varPsi }_{{\\rm{e}},{\\rm{h}}}({\\bf{k}})=\\frac{1}{\\sqrt{2}}\\left[\\begin{array}{c}{{\\rm{e}}}^{{\\rm{i}}{\\varphi }_{{\\bf{k}}}}\\\\ \\pm 1\\end{array}\\right].\\end{array}$$<\/p>\n The out-of-equilibrium dynamics shown in Extended Data Figs. 4<\/a> and 9<\/a>\u00a0were simulated on a 250\u2009\u00d7\u2009250 2D square (kx,\u2009ky) grid in momentum space. The grid is centred on one of the K points with coordinates kx,K\u2009=\u20091.703\u2009\u00c5\u22121 and ky,K\u2009=\u20090 and it has a side length of 0.73\u2009\u00c5\u22121. The simulated time axis ranges from \u2212100\u2009fs to 300\u2009fs, with a step size of 2\u2009fs. The electric field, E(t), is parallel to ky and its temporal profile is analytically fitted to the experimentally measured field. At the initial time instant (t0\u2009=\u2009\u2212100\u2009fs), we assume peh(k,\u2009t0)\u2009=\u20090, whereas the electron distribution corresponds to a Fermi\u2013Dirac distribution. The position of the Fermi level compared with the Dirac crossing is set to \u03b5F\u2009=\u20090.2\u2009eV to match the experiment.\u00a0Simulations comparing the\u00a0dynamics at other driving frequencies, varying field strengths or different\u00a0Fermi levels sometimes\u00a0used analytically modelled waveforms and\u00a0smaller grid sizes to keep up with increased demands for spatial or temporal resolution.\u00a0<\/p>\n Modelling of the scattering rates<\/p>\n In simulations including scattering, the grid was adjusted to 50\u2009\u00d7\u200950 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:<\/p>\n $$\\begin{array}{cc}{\\varGamma }_{{\\bf{k}}}^{{\\rm{e}}}={\\frac{\\partial {f}_{{\\rm{e}}}({\\bf{k}},t)}{\\partial t}|}_{{\\rm{c}}{\\rm{o}}{\\rm{l}}{\\rm{l}}}, & {\\varGamma }_{{\\bf{k}}}^{{\\rm{h}}}={\\frac{\\partial {f}_{{\\rm{h}}}({\\bf{k}},t)}{\\partial t}|}_{{\\rm{c}}{\\rm{o}}{\\rm{l}}{\\rm{l}}}.\\end{array}$$<\/p>\n On these timescales, two physical phenomena mainly contribute to the scattering rates:<\/p>\n $${\\varGamma }_{{\\bf{k}}}^{\\alpha }={\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{p}}}^{\\alpha }+{\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{e}}}^{\\alpha },\\,\\alpha \\in \\{{\\rm{e}},{\\rm{h}}\\},$$<\/p>\n in which \\({\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{p}}}^{\\alpha }\\) and \\({\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{e}}}^{\\alpha }\\) are, respectively, the electron\u2013phonon and electron\u2013electron scattering rates.<\/p>\n Starting from the electron\u2013phonon interaction, the scattering rate can be written as57<\/a>:<\/p>\n $${\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{p}}}^{\\alpha }=\\sum _{v,{\\alpha }^{{\\prime} },{{\\bf{k}}}^{{\\prime} }}\\{{S}_{{\\alpha }^{{\\prime} }\\alpha }^{v}({{\\bf{k}}}^{{\\prime} },{\\bf{k}}){f}_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} },t)[1-{f}_{\\alpha }({\\bf{k}},t)]-{S}_{\\alpha {\\alpha }^{{\\prime} }}^{v}({\\bf{k}},{{\\bf{k}}}^{{\\prime} }){f}_{\\alpha }({\\bf{k}},t)[1-{f}_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} },t)]\\},$$<\/p>\n in which v is the phonon mode considered and \u03b1,\u2009\u03b1\u2032\u2009\u2208\u2009{e,\u2009h} are the band indexes. The first term thus describes transitions from an occupied (k\u2032,\u2009\u03b1\u2032) state to an empty (k,\u2009\u03b1) state owing to the interaction with a v phonon, whereas the second term describes the opposite process. Whereas in the former the state (k,\u2009\u03b1) 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\u2013phonon transition rate:<\/p>\n $$\\begin{array}{c}{S}_{{\\alpha }^{{\\prime} }\\alpha }^{v}({{\\bf{k}}}^{{\\prime} },{\\bf{k}})=\\frac{2{\\rm{\\pi }}}{\\hbar }{|{G}_{{\\alpha }^{{\\prime} }\\alpha }^{v}({{\\bf{k}}}^{{\\prime} },{\\bf{k}})|}^{2}\\\\ \\cdot \\{({n}_{{\\bf{q}}}^{v}+1)\\delta [{\\varepsilon }_{\\alpha }({\\bf{k}})-{\\varepsilon }_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} })+\\hbar {\\omega }_{{\\bf{q}}}^{v}]{\\delta }_{{{\\bf{k}}}^{{\\prime} }={\\bf{k}}-{\\bf{q}}}+{n}_{{\\bf{q}}}^{v}\\delta [{\\varepsilon }_{\\alpha }({\\bf{k}})-{\\varepsilon }_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} })-\\hbar {\\omega }_{{\\bf{q}}}^{v}]{\\delta }_{{{\\bf{k}}}^{{\\prime} }={\\bf{k}}+{\\bf{q}}}\\}\\end{array},$$<\/p>\n in which \\({G}_{{\\alpha }^{{\\prime} }\\alpha }^{v}({{\\bf{k}}}^{{\\prime} },{\\bf{k}})\\) is the electron\u2013phonon coupling constant, \u03b4 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\u2013Einstein distribution.<\/p>\n 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 phonons34<\/a>,47<\/a>. 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\u2013optical phonon coupling constant57<\/a>:<\/p>\n $${|{G}_{{\\alpha }^{{\\prime} }\\alpha }^{{\\rm{LO}}}({{\\bf{k}}}^{{\\prime} },{\\bf{k}})+{G}_{{\\alpha }^{{\\prime} }\\alpha }^{{\\rm{TO}}}({{\\bf{k}}}^{{\\prime} },{\\bf{k}})|}^{2}={|{G}_{{\\alpha }^{{\\prime} }\\alpha }^{{\\rm{OP}}}|}^{2}=1.542\\,{{\\rm{eV}}}^{2}.$$<\/p>\n 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\u2013optical phonon scattering will lead to a rapid spreading of the electron distribution around the Dirac cone.<\/p>\n In the same spirit, the scattering rate for electron\u2013electron interaction can be written as57<\/a>:<\/p>\n $${\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{e}}}^{\\alpha }=\\sum _{{{\\bf{k}}}_{\\ast },{{\\bf{k}}}^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} }}\\sum _{{\\alpha }^{{\\prime} },\\beta ,{\\beta }^{{\\prime} }}\\{{S}_{{\\rm{e}}{\\rm{e}}}({{\\bf{k}}}^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} },{\\bf{k}},{{\\bf{k}}}_{\\ast }){f}_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} }){f}_{{\\beta }^{{\\prime} }}({{\\bf{k}}}_{\\ast }^{{\\prime} })[1-{f}_{\\alpha }({\\bf{k}})][1-{f}_{\\beta }({{\\bf{k}}}_{\\ast })]-{S}_{{\\rm{e}}{\\rm{e}}}({\\bf{k}},{{\\bf{k}}}_{\\ast },{{\\bf{k}}}^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} }){f}_{\\alpha }({\\bf{k}}){f}_{\\beta }({{\\bf{k}}}_{\\ast })[1-{f}_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} })][1-{f}_{{\\beta }^{{\\prime} }}({{\\bf{k}}}_{\\ast }^{{\\prime} })]\\}{\\delta }_{{\\bf{k}}+{{\\bf{k}}}_{\\ast }={{\\bf{k}}}^{{\\prime} }+{{\\bf{k}}}_{\\ast }^{{\\prime} }},$$<\/p>\n 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 (\u03b1\u2032,\u2009k\u2032) and \\(({\\beta }^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} })\\) to the final states (\u03b1,\u2009k) and \\((\\beta ,{{\\bf{k}}}_{\\ast })\\), and the second term refers to the opposite process. The associated transition rate can thus be written as57<\/a>:<\/p>\n $${S}_{{\\rm{e}}{\\rm{e}}}({\\bf{k}},{{\\bf{k}}}_{\\ast },{{\\bf{k}}}^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} })={\\left(\\frac{2{\\rm{\\pi }}}{\\hbar }\\right)}^{2}{|M|}^{2}\\delta [{\\varepsilon }_{{\\beta }^{{\\prime} }}({{\\bf{k}}}_{\\ast }^{{\\prime} })+{\\varepsilon }_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} })-{\\varepsilon }_{\\beta }({{\\bf{k}}}_{\\ast })-{\\varepsilon }_{\\alpha }({\\bf{k}})],$$<\/p>\n in which the matrix element is:<\/p>\n $${|M|}^{2}=\\frac{1}{2}[{|V(q)|}^{2}+{|V({q}^{{\\prime} })|}^{2}-V(q)V({q}^{{\\prime} })]$$<\/p>\n and q\u2009=\u2009|k\u2009\u2212\u2009k\u2032| and \\({q}^{{\\prime} }=|{\\bf{k}}-{{\\bf{k}}}_{\\ast }^{{\\prime} }|\\). The Coulomb potential between electron pairs is:<\/p>\n $$V(q)=\\frac{2{\\rm{\\pi }}{e}^{2}{V}_{{\\rm{e}}{\\rm{f}}{\\rm{f}}}}{{\\epsilon }(q)qA}\\frac{1+\\cos {\\phi }_{{\\bf{k}}{{\\bf{k}}}^{{\\prime} }}}{2}\\frac{1+\\cos {\\phi }_{{{\\bf{k}}}_{\\ast }{{\\bf{k}}}_{\\ast }^{{\\prime} }}}{2}$$<\/p>\n $$V(q{\\prime} )=\\frac{2{\\rm{\\pi }}{e}^{2}{V}_{{\\rm{eff}}}}{{\\epsilon }({q}^{{\\prime} }){q}^{{\\prime} }A}\\frac{1+\\cos {\\phi }_{{\\bf{k}}{{\\bf{k}}}_{\\ast }^{{\\prime} }}}{2}\\frac{1+\\cos {\\phi }_{{{\\bf{k}}}_{\\ast }{{\\bf{k}}}^{{\\prime} }}}{2}$$<\/p>\n with A\u2009=\u20094\u03c02\/\u0394k2. 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<\/a>,58<\/a>, in which:<\/p>\n $${v}_{{\\rm{c}}}(q)=\\frac{2{\\rm{\\pi }}{e}^{2}}{\\kappa q},\\,\\kappa =\\frac{{e}^{2}}{{r}_{{\\rm{S}}}\\hbar {v}_{{\\rm{F}}}},\\,{r}_{{\\rm{S}}}=\\frac{\\hbar {v}_{{\\rm{F}}}}{{\\varepsilon }_{{\\rm{F}}}},\\,D(\\varepsilon )=\\frac{|\\varepsilon |}{2{\\rm{\\pi }}{\\hbar }^{2}{v}_{{\\rm{F}}}^{2}}$$<\/p>\n $$\\mathop{\\Pi }\\limits^{ \\sim }(q)=\\left\\{\\begin{array}{cc}1, & q < 2{k}_{{\\rm{F}}}\\\\ 1+\\frac{{\\rm{\\pi }}q}{8{k}_{{\\rm{F}}}}-\\frac{\\sqrt{{q}^{2}-4{k}_{{\\rm{F}}}^{2}}}{2q}-\\frac{q}{4{k}_{{\\rm{F}}}}\\arcsin \\left(\\frac{2{k}_{{\\rm{F}}}}{q}\\right), & {\\rm{o}}{\\rm{t}}{\\rm{h}}{\\rm{e}}{\\rm{r}}{\\rm{w}}{\\rm{i}}{\\rm{s}}{\\rm{e}}\\end{array}\\right.$$<\/p>\n and rS is the Wigner\u2013Seitz radius of graphene, kF\u2009=\u2009\u03b5F\/(\u0127vF) the Fermi wavevector and D(\u03b5) the density of states close to the Dirac point. This approximation holds for n\u2009\u2265\u2009nth\u2009=\u20091012\u2009cm\u22122, thus if \u03b5F\u2009\u2265\u200937.14\u2009meV, as in our case.<\/p>\n In this form, the scattering rate for electron\u2013electron interaction is computationally extremely expensive. Therefore, we assume that the electron distribution for the second particle can be approximated by the equilibrium Fermi\u2013Dirac distribution of the carriers:<\/p>\n $${f}_{\\beta }({{\\bf{k}}}_{\\ast })\\simeq {f}_{\\beta }^{(0)}({{\\bf{k}}}_{\\ast }),\\,{f}_{{\\beta }^{{\\prime} }}({{\\bf{k}}}_{\\ast }^{{\\prime} })\\simeq {f}_{{\\beta }^{{\\prime} }}^{(0)}({{\\bf{k}}}_{\\ast }^{{\\prime} }).$$<\/p>\n This allows us to write a simplified expression for the electron\u2013electron scattering rate:<\/p>\n $$\\begin{array}{c}{\\varGamma }_{{\\bf{k}},{\\rm{e}}{\\rm{e}}}^{\\alpha }\\simeq \\sum _{{{\\bf{k}}}^{{\\prime} }}\\sum _{{\\alpha }^{{\\prime} },\\beta ,{\\beta }^{{\\prime} }}\\{{S}_{{\\rm{e}}{\\rm{e}},\\beta {\\beta }^{{\\prime} }}^{{\\rm{e}}{\\rm{f}}{\\rm{f}},1}({\\bf{k}},{{\\bf{k}}}^{{\\prime} }){f}_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} })[1-{f}_{\\alpha }({\\bf{k}})]\\\\ \\,-\\,{S}_{{\\rm{e}}{\\rm{e}},\\beta {\\beta }^{{\\prime} }}^{{\\rm{e}}{\\rm{f}}{\\rm{f}},2}({\\bf{k}},{{\\bf{k}}}^{{\\prime} }){f}_{\\alpha }({\\bf{k}})[1-{f}_{{\\alpha }^{{\\prime} }}({{\\bf{k}}}^{{\\prime} })]\\},\\end{array}$$<\/p>\n 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:<\/p>\n $${S}_{{\\rm{ee}},\\beta {\\beta }^{{\\prime} }}^{{\\rm{eff}},1}({\\bf{k}},{{\\bf{k}}}^{{\\prime} })=\\sum _{{{\\bf{k}}}_{\\ast },{{\\bf{k}}}_{\\ast }^{{\\prime} }}{S}_{{\\rm{ee}}}({\\bf{k}},{{\\bf{k}}}_{\\ast },{{\\bf{k}}}^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} }){f}_{{\\beta }^{{\\prime} }}^{(0)}({{\\bf{k}}}_{\\ast }^{{\\prime} })[1-{f}_{\\beta }^{(0)}({{\\bf{k}}}_{\\ast })]{\\delta }_{{\\bf{k}}+{{\\bf{k}}}_{\\ast }={{\\bf{k}}}^{{\\prime} }+{{\\bf{k}}}_{\\ast }^{{\\prime} }}$$<\/p>\n $${S}_{{\\rm{ee}},\\beta {\\beta }^{{\\prime} }}^{{\\rm{eff}},2}({\\bf{k}},{{\\bf{k}}}^{{\\prime} })=\\sum _{{{\\bf{k}}}_{\\ast },{{\\bf{k}}}_{\\ast }^{{\\prime} }}{S}_{{\\rm{ee}}}({\\bf{k}},{{\\bf{k}}}_{\\ast },{{\\bf{k}}}^{{\\prime} },{{\\bf{k}}}_{\\ast }^{{\\prime} }){f}_{\\beta }^{(0)}({{\\bf{k}}}_{\\ast })[1-{f}_{{\\beta }^{{\\prime} }}^{(0)}({{\\bf{k}}}_{\\ast }^{{\\prime} })]{\\delta }_{{\\bf{k}}+{{\\bf{k}}}_{\\ast }={{\\bf{k}}}^{{\\prime} }+{{\\bf{k}}}_{\\ast }^{{\\prime} }}.$$<\/p>\n 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\u2013phonon and electron\u2013electron scattering processes, which are not expected to contribute notably in these conditions34<\/a>. However, for electron\u2013electron 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 Veff\u2009=\u20094.386. It is important to notice that the dependence on the relative angle between k, k\u2032, \\({{\\bf{k}}}_{\\ast }\\) and \\({{\\bf{k}}}_{\\ast }^{{\\prime} }\\) determines the characteristic angular dependence of electron\u2013electron scattering in graphene47<\/a>, also termed collinear scattering41<\/a>. 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.<\/p>\n Relaxation time approximation<\/p>\n 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 approximation6<\/a>,20<\/a>,55<\/a>,56<\/a>. In this context, the scattering rate \\({\\varGamma }_{{\\bf{k}}}^{\\alpha }\\) is assumed to be inversely proportional to a time constant, T1, associated with the mean free time between scattering events. Defining \\({f}_{\\alpha }^{(0)}({\\bf{k}})\\) as the equilibrium electron distribution in the \u03b1 band, we obtain:<\/p>\n $${\\varGamma }_{{\\bf{k}}}^{\\alpha }={\\frac{\\partial {f}_{\\alpha }({\\bf{k}},t)}{\\partial t}|}_{{\\rm{c}}{\\rm{o}}{\\rm{l}}{\\rm{l}}}=-\\frac{{f}_{\\alpha }({\\bf{k}},t)-{f}_{\\alpha }^{(0)}({\\bf{k}})}{{T}_{1}}.$$<\/p>\n 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.<\/p>\n LZM transition region<\/p>\n The original works of Landau22<\/a>, Zener23<\/a>, Majorana24<\/a> and St\u00fcckelberg36<\/a> considered a two-state quantum system with an energy separation \u0394\u03b5\u2009=\u2009\u03b51(t)\u2009\u2212\u2009\u03b52(t)\u2009=\u2009\u03b1t that changes linearly with time and undergoes an avoided crossing with a half gap energy \u03b4. The diabatic transition probability is then given by26<\/a><\/p>\n $${P}_{{\\rm{L}}{\\rm{Z}}{\\rm{M}}}=\\exp (-2{\\rm{\\pi }}\\varGamma )=\\exp \\left(-2{\\rm{\\pi }}\\frac{{\\delta }^{2}}{\\hbar \\alpha }\\right)$$<\/p>\n 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, \u03c4\u2009=\u20092\u03b4\/\u03b5, is short compared with the oscillation period T\u2009=\u20092\u03c0\/\u03c9 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 \u03b4\u2009=\u2009vF\u0127|kx\u2009\u2212\u2009K|, we obtain \\(\\tau =\\frac{\\hbar | {k}_{x}-K|}{|e|{E}_{0}}=\\frac{\\gamma }{\\omega }\\). Because \u03b3\u2009\u2248\u20090.1, the crossing indeed occurs very fast (\u03c4\u2009<\u20091\u2009fs). For the peak electric field E0\u2009=\u20091.57\u2009MV\u2009cm\u22121, 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\/\u03c0 for |kx\u2009\u2212\u2009K|\u2009<\u20090.04\u2009\u00c5\u22122. In this region around the Dirac point, 1\u2009\u2212\u2009e\u22122\u2009=\u200986% 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.\u20098<\/a>).<\/p>\n Injection time analysis<\/p>\n 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.\u20095a<\/a>). 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 interference59<\/a>. 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.\u20099<\/a>). 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.<\/p>\n 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.\u20095b<\/a>). Notably, the weak oscillations in the width of the electron distribution, which are superimposed on the stepwise increase, become more pronounced with scattering (Fig.\u20094c<\/a>). 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 \u03c0\/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.\u20095<\/a>). 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.\u20093c<\/a>). 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.<\/p>\n Data processing<\/p>\n 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.<\/p>\n The raw data from the delay line detector are binned with a spacing of 0.084\u2009eV\u2009\u00d7\u20090.022\u2009\u00c5\u22121\u2009\u00d7\u20090.022\u2009\u00c5\u22121 for the data in Fig.\u20091<\/a> and with a spacing of 0.084\u2009eV\u2009\u00d7\u20090.011\u2009\u00c5\u22121\u2009\u00d7\u20090.011\u2009\u00c5\u22121 for all other photoemission data.<\/p>\n The 3D data in Fig.\u20091c<\/a> are symmetrized to remove the decrease in photoemission intensity caused by the angle of incidence of the XUV pulse.<\/p>\n The 2D momentum maps in Fig.\u20092a<\/a> are obtained by integrating the 3D data within an energy window of 0.33\u2009eV around distinct energy levels (0.25\u2009eV, \u22120.25\u2009eV, \u22120.75\u2009eV, \u22121.25\u2009eV, \u22121.75\u2009eV). 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\u2009\u00c5\u22121) is applied to the data four times. The filtered data are further processed by a 2D curvature filter described by equation\u2009(11) of ref.\u200960<\/a> with C0\u2009=\u20090.2. Finally, a 2D median filter with window size 0.045\u2009eV\u2009\u00d7\u20090.045\u2009\u00c5\u22121 is applied three times to clean up numerical noise.<\/p>\n The band structure cuts in Fig.\u20092b<\/a> are obtained by integrating the 3D data within a momentum window of 0.08\u2009\u00c5\u22121 around kx\u2009\u2212\u2009K\u2009=\u20090. 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\u2009eV and 0.08\u2009\u00c5\u22121, respectively, is applied four times. The filtered data are subject to a 2D curvature filter with C0\u2009=\u20090.2. Finally, a 2D median filter with window size 0.34\u2009eV\u2009\u00d7\u20090.045\u2009\u00c5\u22121 is applied three times. The raw data are shown in Extended Data Fig.\u20092b<\/a>.<\/p>\n The applied field transient and current in Fig.\u20092c<\/a> are extracted from the raw data as laid out in Methods<\/a> sections \u2018Electric field waveform reconstruction and streaking compensation\u2019 and \u2018Current extraction and Fourier analysis of the phase delay\u2019. 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.<\/p>\n The 2D momentum maps in Fig.\u20093a<\/a> are obtained by integrating the 3D data for energies higher than 170\u2009meV. 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\u2009\u00c5\u22121\u2009\u00d7\u20090.02\u2009\u00c5\u22121 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.<\/p>\n To extract the COM motion and change of the width of the carrier distribution shown in Fig.\u20093b,c<\/a>, the plotted 2D maps are integrated along kx and ky. For the data along ky, a Gaussian distribution with a constant offset (accounting for the background) is fitted to the integrated distribution. For the data along kx, the previously constant offset is modified with a linear decrease for growing kx values to account for the decrease in background caused by the photoemission matrix element. The mean of the extracted fit for the ky 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, \u0394\u03c3, of the extracted fits is used to characterize the change of the width of the carrier distribution along kx and ky.<\/p>\n The 2D momentum maps in Fig.\u20094a<\/a> are obtained by integrating the simulated population of the conduction band and applying a Gaussian filter with a small window size of 0.02\u2009\u00c5\u22121\u2009\u00d7\u20090.02\u2009\u00c5\u22121 to the 2D data, as described above for the data shown in Fig.\u20093a<\/a>.<\/p>\n