Ultrafast dynamics of the surface photovoltage in
potassium-doped black phosphorus

Geoffroy Kremer, Maxime Rumo, Changming Yue, Aki Pulkkinen,
Christopher Nicholson, Thomas Jaouen, Fabian von Rohr, Philipp Werner,
Claude Monney

► To cite this version:

Geoffroy Kremer, Maxime Rumo, Changming Yue, Aki Pulkkinen, Christopher Nicholson, et al. Ultrafast dynamics of the surface photovoltage in potassium-doped black phosphorus. Physical Review B, 2021, 104 (3), pp.035125. 10.1103/PhysRevB.104.035125. hal-03285155v2

HAL Id: hal-03285155

https://hal.science/hal-03285155v2

Submitted on 16 Feb 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Ultrafast dynamics of the surface photovoltage in potassium-doped black phosphorus

G. Kremer,¹,*, M. Rumo,¹ C. Yue,¹ A. Pulkkinen,¹,² C. W. Nicholson,¹ T. Jaouen,¹,³ F. O. von Rohr,⁴
P. Werner,¹ and C. Monney¹

¹Département de Physique and Fribourg Center for Nanomaterials, Université de Fribourg, CH-1700 Fribourg, Switzerland

²School of Engineering Science, LUT University, FI-53850, Lappeenranta, Finland

³Univ Rennes, CNRS, (IPR Institut de Physique de Rennes), UMR 6251, F-35000 Rennes, France

⁴Department of Chemistry, University of Zurich, CH-8000 Zurich, Switzerland

<figure><img src="image_2.png" /></figure>

(Received 23 February 2021; revised 29 April 2021; accepted 23 June 2021; published 12 July 2021)

Black phosphorus is a quasi-two-dimensional layered semiconductor with a narrow direct band gap of 0.3 eV. A giant surface Stark effect can be produced by the potassium doping of black phosphorus, leading to a semiconductor to semimetal phase transition originating from the creation of a strong surface dipole and associated band bending. By using time- and angle-resolved photoemission spectroscopy, we report the partial photoinduced screening of this band bending by the creation of a compensating surface photovoltage. We further resolve the detailed dynamics of this effect at the pertinent timescales and the related evolution of the band structure near the Fermi level. We demonstrate that after a fast rise time, the surface photovoltage exhibits a plateau over a few tens of picoseconds before decaying on the nanosecond timescale. We support our experimental results with simulations based on drift-diffusion equations.

DOI: 10.1103/PhysRevB.104.035125

I. INTRODUCTION

Semiconducting materials are the building block of modern
technologies such as electronics, optoelectronics or com-
puting devices. Their band gap is a fundamental property
governing their optical and electrical abilities, leading in re-
cent decades to an increasing amount of research aiming to
enhance the control of both its size and its nature (direct or
indirect). In this context, quasi two-dimensional (2D) ma-
terials have appeared as a new promising emergent family
[1]. Indeed, 2D semiconducting materials form a rich class
of compounds whose electronic band gap can vary from
0 to more than 7 eV, i.e., ranging from the tera-hertz to
the infrared in monolayer Xenes [2–10], to the visible with
transition metal dichalcogenides (TMDCs) [11] and finally
to the ultraviolet with hexagonal boron nitride [12] and 2D
oxides [13–16]. In addition to covering a large part of the
electromagnetic spectrum, the main advantage of quasi 2D
materials is the possibility to tune their band gap with nu-
merous approaches such as strain [17,18], electric fields [19],
surface gating [20], chemical substitution [21,22], dimen-
sionality reduction [23,24], or by the creation of complex
heterostructures [25]. Besides the band-gap tunability, this
opens a broad range of possibilities to create and engi-
neer new quantum materials with unprecedented physical
properties.

Among all these materials, black phosphorus (BP) is one
of the most prominent because of its high carrier mobilities
[26,27] and its tunable band gap over the whole infrared range,
giving it a singular place between graphene and the TMDCs
[28]. Black phosphorus is a layered semiconductor composed

of only phosphorus atomic planes which are weakly inter-
acting via out-of-plane van der Waals forces. Its band gap
is highly dependent on the thickness, ranging from 0.3 eV
in the bulk (BP) to 1.5 eV in the monolayer case (phospho-
rene). In addition to this intrinsic variety, external stimuli
such as mechanical strain or surface chemical gating have
been demonstrated as powerful tools to modify its electronic
properties [29].

Chemical gating has attracted a lot of attention since the BP surface is very sensitive to external electric fields [30,31] and induced surface dipoles [32,33]. Indeed, by doping the surface of BP with alkali atoms, it is possible to create a strong dipole generating an electric field responsible for the continuous reduction of the band gap at the surface. At a sufficient alkali coverage the band gap closes, as recently demonstrated by angle resolved photoemission spectroscopy (ARPES), creating a remarkable semiconductor to semimetal transition [32]. The electric field of alkali-doped BP is related to the large electron transfer from the alkali atoms to the BP. It creates an excess of free carriers, leading to the formation of a space charge region (SCR) and band bending (BB), a generic effect in semiconductors physics [34]. There is also a charge redistribution between the valence band (VB) and the conduction band (CB) whose wave functions are respectively pushed to the bulk and confined at the surface of the material [30–33]. In particular, the CB forms a two-dimensional electrons gas (2DEG) [32,35] and experiences a larger BB, which can be decomposed into the same BB feeled by the VB (φ<sub>BB</sub>) and a second one we name Stark potential (φ<sub>S</sub>), localized at the extreme surface. At a sufficient coverage, the total CB bending is enough to close the surface band gap but the bulk one remains the same. In analogy with the atomic Stark effect created by an electric field, this is called the surface giant Stark effect.

<footnote>*Corresponding author: geoffroy.kremer@unifr.ch</footnote>
In recent works [36,37], it has been shown that $\phi_{BB}$ in alkali-doped BP at this critical coverage equals a few hundreds of meV and can be compensated via the illumination of the surface. Under the preexisting electric field, photogenerated electrons and holes are redistributed in the SCR and create a compensating field called surface photovoltage (SPV). Hedayat *et al.* reported that this is felt as an energy shift of the whole surface photoemission spectrum but also as a spectral shape modification of the VB [37]. However, their work mainly focused on the transient dynamics up to 5 ps after the photoexcitation.

Here we report on the temporal evolution of the SPV and the band structure near the Fermi level ($E_F$) of potassium- (K) doped BP, from a few ps before to a few hundreds ps after photoexcitation. To do so, we use time-resolved ARPES (Tr-ARPES) which is the ideal tool to study the ultrafast dynamics of the electronic structure of materials [38,39]. Tr-ARPES is particularly well suited since it allows the generation and the study of the SPV [40–43] with energy and time resolutions. Whereas it appears quasi instantaneously, the SPV is remarkably stable over 60 ps and decreases on the nanosecond (ns) timescale. We demonstrate that the characteristic timescales of the SPV are correlated to the spectral shape modification of the VB, which is highly affected by the potential profile in the SCR and as a consequence by the transient screening of $\phi_{BB}$. Furthermore, we solve the drift-diffusion equations (DDEs) for K-doped BP and compare the simulated temporal evolution of the BB to the measured SPV dynamics. Our simulations show a good qualitative agreement and unambiguously confirm the ultrafast transient evolution of $\phi_{BB}$ after the photoexcitation of the pump. Finally, the calculations demonstrate the importance of the electron-hole recombination processes for the long-time evolution.

## II. METHODS

Large single crystals of BP were synthesized by a modified chemical vapor transport technique similar to the method first reported in Ref. [44]. Dry red phosphorus (150 mg, purity 99.999%), dry $SnI_4$ (5 mg, purity 99.999%), tin shots (40 mg, purity 99.98%-purity), and gold powder (20 mg, purity 99.999%) were sealed in an evacuated quartz glass tube. The mixture was placed in a box furnace at 680°C and held at this temperature for 24 h. It was then cooled with 0.1 K/h to 500°C and quenched to room temperature (RT) by removal from the furnace. The crystals were manually removed from the metals. Residual $SnI_4$ was removed from the crystals by sealing them in an evacuated quartz tube and heating them to 400°C for 5 h in a tube furnace. Black phosphorus samples were cleaved at RT in a pressure of $1 \times 10^{-8}$ mbar. ARPES measurements were performed using a ScientaOmicron DA30 photoelectron analyzer with energy and angular resolutions better than 10 meV and 0.1°, respectively. The base pressure in the photoemission chamber was better than $3 \times 10^{-11}$ mbar. For the static measurements, we used monochromatized $HeI_\alpha$ radiation ($h\nu = 21.2$ eV, SPECS UVLS with TMM 304 monochromator). For the Tr-ARPES measurements, half of the power of a femtosecond laser (Pharos, Light Conversion, operating at 1030 nm) is converted into 780-nm light with an optical parametric amplifier, which

is then frequency-quadrupled to 6.3 eV in $\beta$-BaB$_2$O$_4$ crystals to generate probe pulses [45]. The pump-probe measurements were performed at 200 kHz and the pump fluence was 160 $\mu$J/cm$^2$. The intrinsic resolution of the probe laser is 47 meV, as determined by the fit of the Fermi edge of a polycrystalline metal. The remaining half of the fundamental laser power is directed into a collinear optical parametric amplifier (Orpheus, Light Conversion) to generate pump pulses at 780 nm (1.55 eV) with a duration of about 50 fs. The temporal resolution was estimated to be better than 100 fs by measuring the width of the cross-correlation between the pump and the probe pulses. We used $s$ and $p$ polarization for the probe and an incident angle of 55° with respect to the normal of the surface to measure along the armchair (AC) and zigzag (ZZ) directions, respectively. The Tr-ARPES data have been acquired with a negative bias of -5 V. All the photoemission measurements were realized at 30 K if not further specified. The K evaporation was performed at 30 K from a SAES getter source in a pressure better than $2 \times 10^{-10}$ mbar.

Density-functional theory (DFT) calculations were performed using the projector augmented wave method implemented in the Vienna ab-initio software package [46–51], with a kinetic energy cutoff set to 400 eV and a k point spacing <0.2 Å$^{-1}$. Exchange and correlation effects were included within the strongly constrained and appropriately normed (SCAN [52]) meta-generalized gradient approximation. The structure of BP was relaxed until the residual forces on atoms were smaller than 1 meV/Å. The lattice parameters obtained after the relaxation of the structure with DFT calculations in the SCAN approximation are $a = 4.53$ Å, $b = 3.29$ Å, and $c = 10.91$ Å, in good agreement with x-ray diffraction measurements [53].

## III. RESULTS

### A. Static characterization of black phosphorus

The crystal structure of BP is shown in Fig. 1(a). It crystallizes in an orthorhombic crystal structure (space group *Cmca*) with superimposed buckled layers of phosphorus atoms along the z axis. The layered structure of BP and the weak interlayer van der Waals interactions favor mechanical cleaving of the material with the scotch-tape technique. The strongly anisotropic nature of BP in the basal plane is visible in its two high-symmetry directions along x and y, namely the AC and ZZ directions. The corresponding three dimensional Brillouin zone (BZ) and its 2D projection on the (001) plane (green), i.e., the cleave plane, are shown in Fig. 1(b).

The static electronic band structure of pristine BP along the AC direction measured by ARPES is shown in Fig. 1(c) (note that the signal has been symmetrized with respect to the L point). At this particular photon energy, we almost probe the in-plane Z-L-Z high-symmetry line of the BZ [55] [thick red line in Fig. 1(b)]. The ARPES intensity exhibits all the spectroscopic signatures of a clean BP surface [56]. In these energy and momentum ranges, the band structure is composed of six dispersive bands with the VB top presenting a maximum at the Z point and located about 0.1 eV below $E_F$. In Fig. 1(d), we compare our experimental data with bulk projected DFT calculations, i.e., the superposition of the calculated band
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Crystal structure of BP plotted with the VESTA visualization program [54] showing the two inequivalent in-plane directions, namely the ZZ and AC directions. (b) Corresponding three dimensional BZ and its 2D projection on the (001) plane (green). (c) ARPES spectrum along the Z-L-Z high-symmetry line [thick red line in (b)] measured at T = 300 K and $h\nu = 21.2$ eV. (d) Bulk projected electronic band structure [red plane in (b)] calculated from DFT in the SCAN approximation.</figcaption></figure>

structure of BP along the $k_z$ direction, from the $\Gamma$ to the Z point of the 3D BZ [red plane in Fig. 1(b)]. The overall agreement between ARPES and DFT is excellent. Furthermore, the intensity variations and the energy broadening of the bands are also well reproduced.

## B. Dynamics of the surface photovoltage in K-doped black phosphorus

Since the static ARPES only gives access to the occupied part of the band structure, we now turn our attention to Tr-ARPES measurements in order to obtain the band gap of BP. Figures 2(a) and 2(b) show the corresponding data taken at a pump-probe delay of -1 ps and +2 ps along the ZZ direction, respectively. We are using this orientation in the following because of detrimental matrix element effects along the AC direction. At negative delays, we only see the occupied part of the band structure as in the static case (shifted by the SPV as discussed below). At positive delay, the CB is now transiently populated and can be probed. Thus, we can directly estimate the band-gap amplitude of the material to be 0.3 eV, in good agreement with previous measurements [57,58] and our DFT calculations from Fig. 1(d) which give a band-gap value of 307 meV.

As visible in Fig. 2(b), $E_F$ is closer to the VB than to the CB. Since the effective masses of these bands are almost the same in pure BP [26], it means that our sample is p doped. Recent scanning tunneling microscopy (STM) experiments

have demonstrated that p-doped BP samples have an important density of surface acceptor states [58]. This is confirmed by our measurements showing that photoexcitation of the pump induces a global shift of -13 meV of the photoemission spectrum (see Fig. 7), i.e., a shift to low kinetic energy (KE), originating from a negative SPV as expected by the existence of an upward BB at the surface of our sample [34]. This negative SPV compensates the initial upward BB at the surface of our p-doped BP sample. We illustrate the situation in Fig. 3(a), where we plot the energy band diagrams of the p-doped BP surface. We consider the nonequilibrium, the equilibrium and finally the photoexcited configurations.

At a negative pump-probe delay, the photoelectrons traveling from the surface to the photoemission analyzer are decelerated/accelerated (depending on the initial BB curvature) by the electric field generated by the SPV due to the pump pulse arriving afterwards [42,43]. The experimental consequence is a rigid KE shift of the entire photoemission spectrum. In other words, at negative delays we probe the spatial distribution of the electric field generated outside the sample by the SPV but not the photoexcited band structure of the material. At positive delays, we are sensitive to the modification of the band structure induced by the SPV.

Subsequently, we have deposited K atoms on the BP surface until reaching the semiconductor to semimetal transition, as depicted in Fig. 2(c). This panel shows the Tr-ARPES intensity taken at -1 ps pump-probe delay. Compared to the pristine BP surface, the VB has shifted down by 200 meV relative to $E_F$, whereas the CB is now visible in the occupied states and touches the top of the VB. By doping the surface with electrons from donor atoms, we switch from an upward to a strongly downward BB configuration. This is the reason why we now observe a positive SPV at the origin of a global shift to high KE of the whole photoemission spectrum (see Fig. 8 for measurements with and without pump). Both VB and CB feel this common BB, namely $\phi_{BB}$, but the CB also experiences a Stark potential at the extreme surface and forms a 2DEG. As a consequence, its total bending is larger than the VB one and, at a sufficient K coverage, lead to a collapse of the surface band gap. At the low fluence used in our experiment, we do not screen this strong Stark potential which is confined in the first layers of BP and, consequently, the 2DEG of the CB persists. The coexistence of a long range BB and a strong surface potential appears to be a common effect in layered semiconducting materials [59]. Similarly to the VB, the 2DEG is only energy shifted by the SPV and the associated screening of $\phi_{BB}$ [see Fig. 3(b)].

Let us now discuss the measurement taken at +2 ps and presented in Fig. 2(d). The most striking effect relates to the top of the VB which is shifted to low KE: the gap between the CB and the VB is apparently reopened. This is even more evident by looking at the energy-distribution curves (EDCs) taken at -1 ps and +2 ps [see Fig. 2(e)] or the difference in ARPES intensity between these two pump-probe delays [see Fig. 2(f)]. It is worth to mention that this effect is independent of the crystal orientation, as shown in Fig. 9. Nevertheless, it is more convenient to work along the ZZ line because the CB has no spectral weight at the normal emission along the AC direction due to matrix element effects [60].
<figure><img src="image_4.png" /><figcaption>FIG. 2. ARPES spectra recorded at pump-probe delays of [(a) and (c)] −1 ps and [(b) and (d)] +2 ps for bare BP and K-doped BP. In panel (b), the intensity at high KE has been multiplied by a factor 100 with respect to the low KE part in order to visualize the CB and the VB on the same linear color scale. The data have been acquired along the ZZ direction, i.e., the *M*-Z-*M* high symmetry direction of the BZ. (e) EDCs at −1 ps (red) and +2 ps (blue) integrated ±2° around the normal emission for K-doped BP. We highlight the position of the leading edge inflexion point of the VB, which is the marker we use in the following to obtain the evolution of the energy shift induced by the SPV. (f) Difference between the ARPES spectra of panels (c) and (d). Red and blue colors respectively correspond to an increase and a depletion of photoemission intensity. (g) and (h) are sketches of K-doped BP photoemission spectra from panels (c) and (d) at negative and positive pump-probe delays.</figcaption></figure>

The EDCs in Fig. 2(e) clearly show that the spectral shape of the VB is unambiguously affected by the photoexcitation at positive pump-probe delay. Indeed, the leading edge of the VB is strongly sharpened, leading to a depletion of intensity at 6.95 eV and an increase of intensity at 6.8 eV KE, visible as a blue/white/red contrast in the difference ARPES intensity plot in Fig. 2(f). In contrast, the CB spectral shape is marginally affected. This confirms that the effect of the SPV on the VB is not limited to a global energy shift of the band structure; it also modifies its spectral shape [37].

At negative delay, we probe the equilibrium state of the system for which the common BB, $\phi_{BB}$, is not yet screened. Since the VB is delocalized in the bulk of the material and the 6.3 eV laser is bulk sensitive (several nm) due to the extraction of very low kinetic energy electrons [61], we integrate both the contributions of the VB from the surface and the bulk which are energy separated by the amplitude of $\phi_{BB}$, i.e., almost 200 meV [see Fig. 3(b)]. The corresponding sketch of the band structure is represented in Fig. 2(g). It explains why the VB width is larger at negative delay, as observed in Fig. 2(e). On the contrary, the CB forms a 2DEG confined at the surface and is visible as a single band we call $CB_{surface}$, which aligns with the bulk contribution of the VB.

At positive delay, we now probe the photoexcited system at a pump fluence of 160 µJ/cm². At this fluence, the SPV is saturated [36] and we recover a flat band regime for the VB: surface and bulk contributions of the VB are now aligned [see Fig. 3(b)]. In contrast, the 2DEG persists because the Stark potential is not screened at this fluence. The consequence is that the bulk contribution of the VB shifts relatively to the surface one and to $E_F$, leading to a decrease of the VB width and to a sharpening of its leading edge, as visible in Fig. 2(e). This is also the reason why $CB_{surface}$ and the bulk contribution of the VB are now energy separated and why we observe an apparent transient band-gap reopening. This second situation is schematized in Fig. 2(h). It is apparent because, in reality, the surface band gap was not closed since $CB_{surface}$ and $VB_{surface}$ are energy separated by 200 meV, as revealed once $\phi_{BB}$ is compensated by the SPV [Fig. 2(d)]. With our K coverage, the surface band gap has not yet collapsed and has only been reduced by 100 meV with respect to the bulk band gap of BP.

Nevertheless, it is an additional evidence that the surface band gap can be reduced with respect to the bulk by alkali doping. By suppressing the photoemission artifact related to the integration of the VB along the out-of-place direction, we access the true surface band-gap amplitude. For this reason,
<figure><img src="image_5.png" /><figcaption>FIG. 3. Schematic diagram of electron energy levels near the surface of (a) p-doped BP sample and (b) K-doped BP with and without illumination. For the p-doped (K-doped) BP surface, an upward (downward) long-range BB and consequently a negative (positive) SPV are expected. In the K-doped case, the CB is confined at the surface and experiences an additional Stark potential which is not screened by the pump excitation in our fluence regime. This Stark potential is responsible of the progressive collapse of the surface band gap as a function of K coverage. Note that the scales in panels (a) and (b) are not respected.</figcaption></figure>

the claim of a semiconductor to semimetal transition has to be used with care since it is only true (i) in the first layers of BP where the Stark potential exists and (ii) it is dependent on the depth integration of the technique, which varies with the photon energy in ARPES.

Thus, in this scenario the spectral shape evolution of the VB is directly correlated to the photoinduced SPV at positive delay. To confirm this, it is important to show that the SPV dynamics is unambiguously related to the dynamics of the spectral shape of the VB. First, we characterized the energy shift induced by the SPV both at negative and positive delays over a few hundreds of ps. To do so, we estimated the energy shift of the photoemission spectra as extracted from the EDCs analysis as a function of time [see a few EDCs at different pump-probe delays in Fig. 10(a)]. The energy shift has been evaluated as the evolution of the position of the leading edge inflexion point of both the VB and the CB. The corresponding results are presented in Figs. 4(a) and 4(b). It is important to mention that the energy position of the inflexion point is independent of the evolution of the spectral shape and consequently only reflects the evolution of the experimentally measured SPV (in other words the position of CB<sub>surface</sub> and VB<sub>surface</sub>).

In both cases, the SPV is increasing at negative delays up to a maximum around time zero. Then, it exhibits a remarkable plateau over 60 ps before slowly decaying over a few hundreds of ps. We have fitted the recovery of the SPV with an exponential fit. The extracted time constants for the CB and the VB are >1 ns [see Figs. 11(a) and 11(b)]. The temporal evolution of the photoemission cutoff position shows the same features

<figure><img src="image_6.png" /><figcaption>FIG. 4. Temporal evolution of the surface photovoltage (SPV) at the surface of K-doped BP for (a) the CB and (b) the VB. Panels (c) and (d) represent EDCs taken at normal emission, respectively for the CB and the VB, corrected by the SPV determined from panels (a) and (b) for a few pump-probe delays. (e) Temporal evolution of the integrated intensity $I_{CB}$ under the EDCs in panel (c) (gray area). The inset corresponds to a zoom over the first ps. (f) Same for $I_{VB}$ [red area in panel (d)].</figcaption></figure>
[see Fig. 10(b)]. Thus, the energy shift induced by the SPV is the same for all the photoemission spectrum and not restricted to a particular KE range.

Now, we address the temporal evolution of the spectral shape of the bands. In order to remove the influence of the energy shift induced by the SPV, we have plotted the EDCs corrected by the above determined SPV for both the CB and the VB, as shown in Figs. 4(c) and 4(d) for selected pump-probe delays. To understand the temporal evolution of the spectral shape, we have plotted the variation of the photoemission intensity above the energy position of the inflexion points for both the CB and the VB (it would be the complementary of taking an integration region below these points). The integration regions are highlighted by the gray and red rectangles. These plots are represented in Figs. 4(e) and 4(f) and respectively called $I_{CB}$ and $I_{VB}$.

The temporal evolution of $I_{CB}$ shows a strong increase around time zero before decaying over a few ps to its equilibrium value [see inset in Fig. 4(e)]. It can be reproduced by an exponential fit with a time constant of 2 ps [see Fig. 11(c)] which is expected for the transient dynamics of the electronic temperature [62]. Thus, the effect of the SPV on the CB is limited to an energy shift and its spectral shape is not modified as expected by the strong surface confinement of electrons from the CB in the above mentioned scenario.

On the contrary, the evolution of $I_{VB}$ is more complex. Over the first ps after photoexcitation, $I_{VB}$ is drastically reduced: It corresponds to the intensity depletion we discussed in Fig. 2. This decrease can be reproduced by an exponential fit with a time constant of 520 fs [see Fig. 11(e)]. During the next 60 ps, $I_{VB}$ shows a plateau before slowly recovering with a time constant of 1.1 ns [see Fig. 11(d)]. The plateau and the recovery timescales of $I_{VB}$ are in very good agreement with the ones observed for the SPV. Thus, it unambiguously demonstrates that the spectral shape evolution of the VB in K-doped BP is directly related to the magnitude of the SPV in the material. We can further note that, as expected and contrary to the experimentally measured SPV, $I_{VB}$ is not affected at negative delays because we probe the band structure of the unexcited system, i.e., before the pump pulse arrival.

### C. Discussion and simulations

We now summarize the above discussed evolution of the band structure, both in the static and dynamic cases. By doping the surface of a p-doped BP sample with K, we create a strong surface dipole which is the origin of the surface confinement and the shift of the CB down to the occupied states. At a sufficient coverage, the $CB_{surface}$ and the bulk contribution of the buried VB superimpose, visible as a band-gap closing in static ARPES [Figs. 5(a) and 5(b)]. By photoexciting the system with an intense infrared laser pulse, we induce and reveal the dynamics of the SPV, which is responsible for the shift of the whole surface photoemission spectrum and for the transient modification of the VB spectral shape.

The pertinent timescales are given in Fig. 5(c) and can be decomposed into four main steps. (i) At negative delays, since the probe induces photoemission before the pump arrival, we are sensitive to the nonphotoexcited band structure of the system. Nevertheless, the photoemission spectrum is shifted

to higher KE by the electric field generated afterwards by the pump pulse. The smaller the pump-probe delay, the larger the field experienced by the outgoing photoelectrons. At a sufficient time before the pump excitation, the photoelectrons are not sensitive to this electric field anymore and the band structure looks the same as in the static case (see the extreme left panel). (ii) During the first picoseconds after photoexcitation, i.e., at positive delays, two processes are occuring. First, a SPV quickly builds up and modifies the energy level diagram of the system. The VB is largely affected because it is delocalized in the bulk of BP while the surface confined $CB_{surface}$ is not altered. The consequence is the transient modification of the VB spectral shape and the apparent reopening of the band gap. Second, the pump pulse drives the system to a highly excited state. The electronic temperature increases, progressively equilibrates with the lattice and decreases: The consequence is a transient broadening of the Fermi-Dirac

<figure><img src="image_7.png" /><figcaption>FIG. 5. [(a) and (b)] Cartoon representation of the band structure before and after K deposition on BP as measured by ARPES. The strong surface dipole originating from alkali deposition leads to a progressive surface semiconductor to semimetal transition named giant surface Stark effect. (c) Cartoon of the temporal evolution of K-doped BP band structure at the pertinent timescales as measured by Tr-ARPES: ① Photoelectrons acceleration before time zero, ② SPV build up, band structure modification and excited electrons thermalization, ③ Plateau, and ④ SPV relaxation and band structure recovery. (d) Schematic representation of the evolution of the SPV only, as experimentally measured with Tr-ARPES (red) and as expected from theory (blue).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Nonequilibrium dynamics of photoexcited carriers obtained by solving the DDEs. (a) The SPV as a function of time after the pump for different recombination rates $B_r$. When $B_r = 0$ (blue dashed line), there is no recombination. The inset shows the fast initial dynamics in the range $0 \le t < 10$ fs. (b) The band energy as a function of $x$ at different times. The upper panel is for $B_r \ne 0$ and the lower panel is for $B_r = 0$ [similarly in panels (c) and (d)]. The red-dashed line in the lower panel shows the solution of the steady state after the pump when $B_r = 0$. (c) The internal electric field as a function of $x$ for different $t$. (d) The net excess carrier density ($\delta n - \delta p$) as a function of $x$ for different $t$. The parameters used in the simulation are shown in Table I.</figcaption></figure>

distribution in the first ps after the photoexcitation. (iii) Once the electronic temperature equilibrates and the band gap apparently reopens, the system exhibits a remarkable quasi steady state over 60 ps. (iv) Then, the SPV progressively recovers with a time constant of the order of ns. The consequence is a progressive shift to lower KE this time and a recovery of the VB spectral shape. At a sufficient time after the pump excitation, the SPV has fully recovered and the band structure is similar to the one observed a long time before the photoexcitation or in the static case (see the extreme right panel).

In Fig. 5(d), we give a schematic representation of the SPV time evolution as experimentally measured with Tr-ARPES (red curve), and as we would expect from theory (blue curve). Before time zero, the progressive increase of the experimentally extracted SPV is an artifact of the pump-probe approach, as discussed above. As schematized in the step 1, we are theoretically expecting a zero SPV in this time regime (see blue curve). At time zero, the system is photoexcited and the SPV builds up (step 2). As discussed below, it occurs on the timescale of the pump pulse ($\approx 50$ fs). In theory, the SPV starts from zero before building up to its maximum after this characteristic time (without step 1). This is in contrast with the experimentally measured SPV which is already at an almost maximum value around time zero due to the step 1, further complicating its experimental resolution. Steps 3

and 4 are similar in both cases with a plateau and relaxation processes.

In order to simulate the nonequilibrium dynamics of photoexcited carriers in K-doped BP at $T = 30$ K we solve the DDEs with realistic parameters (see Appendix). Since we are looking at the effect of the photoexcitation of the pump on the common long-range BB, i.e., $\phi_{BB}$, we assume that the depth profile potential of the CB and the VB are the same. The results are shown in Fig. 6. While electron-hole recombination can be neglected on the sub-ps timescales of the simulation, we also show results with an artificially short life time of the photocarriers. In these simulations, the band-to-band recombination [63] between excess electron-hole pairs results from a term $B_r n \delta p$ ($n = n_0 + \delta n$ with $n_0$ the static free electron distribution) in the DDEs, where $B_r$ is the recombination rate and $\delta n$ ($\delta p$) is the excess density of electrons (holes) generated by the laser pump. The SPV amplitude $\phi_{SPV} = -e\phi(x=0)$ is plotted as a function of time in Fig. 6(a). An ultrafast increase of $\phi_{SPV}$ is observed for $t < 0.1$ fs, as a result of two main factors. On the one hand, there is a large discrepancy between the high mobility of the holes, $\mu_p \approx 6000$ cm$^2$/Vs, and the smaller mobility of electrons $\mu_n \approx 540$ cm$^2$/Vs at $T = 30$ K [26]. On the other hand, a strong internal electric field ($E \sim 10^7$ V/m) exists in the surface region due to the spatial separation between the K ions and the free electrons (see Fig. 12). The drift velocity of the holes and electrons can
<table><caption>TABLE I. Parameters used in the simulation. The diffusion coefficients are obtained by the Einstein relation $D_{n,p} = \frac{\mu_{n,p}k_B T}{e}$. The values of $\mu_n$, $\mu_p$, and $\epsilon_r$ are from Ref. [66]. The coefficient $\alpha$ in $B_r$ is a unitless number.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Temperature</td><td>T</td><td>30 K</td></tr><tr><td>Relative permittivity</td><td>ε<sub>r</sub></td><td>8.3</td></tr><tr><td>Electron mobility</td><td>μ<sub>n</sub></td><td>540.0 × 10<sup>-4</sup> m<sup>2</sup>/(Vs)</td></tr><tr><td>Electron diffusion coefficient</td><td>D<sub>n</sub></td><td>1.396 × 10<sup>-4</sup> m<sup>2</sup>/s</td></tr><tr><td>Hole mobility</td><td>μ<sub>p</sub></td><td>6000 × 10<sup>-4</sup> m<sup>2</sup>/(Vs)</td></tr><tr><td>Hole diffusion coefficient</td><td>D<sub>p</sub></td><td>15.51 × 10<sup>-4</sup> m<sup>2</sup>/s</td></tr><tr><td>Potassium coverage</td><td>σ<sub>0</sub></td><td>6.3 × 10<sup>17</sup> m<sup>-2</sup></td></tr><tr><td>Photoexcited charge density</td><td>σ<sub>p</sub></td><td>~4σ<sub>0</sub></td></tr><tr><td>Length unit</td><td>l<sub>0</sub></td><td>10<sup>-10</sup> m (1 Å)</td></tr><tr><td>Sample length</td><td>L</td><td>1000l<sub>0</sub></td></tr><tr><td>Number of x-grid points</td><td>N</td><td>20 001</td></tr><tr><td>Electric field unit</td><td>E<sub>u</sub> = σ<sub>0</sub>e / (2ε<sub>r</sub>ε<sub>0</sub>)</td><td>6.867 × 10<sup>8</sup> V/m</td></tr><tr><td>Density unit</td><td>n<sub>u</sub> = σ<sub>0</sub> / (l<sub>0</sub>α)</td><td>6.3 × 10<sup>27</sup> m<sup>-3</sup></td></tr><tr><td>Radiative recombination rate</td><td>B<sub>r</sub> = 0.01fs·n<sub>u</sub></td><td>1.587α × 10<sup>-9</sup> m<sup>3</sup>/s</td></tr></tbody></table>

reach $v_p = \mu_p E \sim 56$ Å/fs and $v_n = \mu_n E \sim 5$ Å/fs, respectively [64]. As a result, the holes flow rapidly into the bulk (the direction of the net field), while the electrons flow toward the surface. The time evolution of the excess carrier density, $\delta n - \delta p$, plotted in Fig. 6(d), illustrates how excess electrons rapidly accumulate near the surface while an excess of holes rapidly builds up slightly away from the surface. The effect of this reshuffling of charge is a rapid weakening of the internal electric field, as shown in Fig. 6(c), and a rapid flattening of the energy band $-e\phi(x)$, i.e., $\phi_{BB}$, as shown in Fig. 6(b).

The diffusion term becomes important when the density profile is sharp. It causes particles to flow from the high density region to the lower density region. As we can see from Fig. 6(d), the surface peak of $\delta n - \delta p$ decreases with time and a broader peak appears which moves toward the bulk. At $t = 0.1$ fs, the internal electric field $E$ stops weakening and $\phi_{SPV}$ reaches its maximum value as shown in Fig. 6(a). For $t > 0.1$ fs, the surface peak in the net charge density begins to disappear slowly and shifts further into the bulk region. As a result, the internal electric field $E$ starts to increase again and some BB reappears, as shown in Fig. 6(b). We should note that these simulations are for an instantaneous photo-doping at $t = 0_+$. In our experiment, the pump pulse has a duration of about 50 fs. Hence, in the experiment, the flattening of the band should occur on a timescale controlled by the pump pulse.

Without recombination between electrons and holes ($B_r = 0$), the system will relax to a new steady state controlled by the different electron and hole mobilities. As shown by the lower panel of Fig. 6(d), the time evolution of $\delta n - \delta p$ approaches this steady-state solution (red dashed line) rather slowly after $t = 0.1$ fs. This solution corresponds to an almost flat band, which is distinctly different from the initial equilibrium solution [green line in Fig. 6(b)]. Hence, as shown by the blue dashed line in Fig. 6(a), $\phi_{SPV}$ will eventually reach a plateau much higher than its initial value at $t = 0$.

If we switch on recombination processes ($B_r > 0$) (with a rate much larger than in the realistic system), the $\delta n - \delta p$ peak starts to decay even on the timescale of our simulation, see the upper panel in Fig. 6(d). Now, the system will return

to the initial equilibrium state and we observe a continuous recovery of the BB, as shown by the green and red lines in Fig. 6(a) and also in the upper panel of Fig. 6(b). The higher the recombination rate, the faster the return to the initial state.

These simulations present a good qualitative agreement with our Tr-ARPES measurements. Very quickly after the system is excited by the pump pulse, the photogenerated charge carriers are redistributed in the SCR under the effect of the pre-existing electric field and diffusion. This charge redistribution leads to a transient compensation of the electric field at the surface, and as a consequence to a compensation of the initial BB. As experimentally extracted from the transient evolution of the VB [see Fig. 4(f)], the SPV builds up in the first ps. According to our quench calculations, it happens on the 0.1 fs timescale [see Fig. 6(b)] [64]. There are two reasons for this discrepancy. First, in the experiment the photoexcitation occurs with a finite duration of about 50 fs and the total time resolution of our experiment is 100 fs, well above this timescale. Second, it has been demonstrated in the pristine BP case that a fast transient broadening of the VB happens in the first ps after the photoexcitation [60]. This effect is still present in the K-doped case [37] and explain the very fast increase of intensity in Fig. 4(f), which is quickly vanishing and counterbalanced by the depletion caused by the SPV. Consequently, the finite pulse duration and transient broadening explain why the time trace in Fig. 4(f) is broadened and why the measured time constant of the SPV build up has a finite duration >500 fs.

Then, the delayed recombination [the plateau on Fig. 5(d)] regime is experimentally seen, but not captured by our theoretical model. Such features have already been observed in other systems such as silicon oxide [40] or GaAs [41] but, to the best of our knowledge, no model has been able to explain it. We speculate that the origin of this delayed recombination could be due to an additional effect, not captured by our theoretical model, that prevents the recombination of electrons and holes immediately after the pump, like spatial separation or the trapping by defects. Another possibility is that, in this temporal regime, recombination takes place but is exactly compensated by another effect which is unknown to us. We
leave the explanation of the origin of this feature for future
works.

Finally, the calculations unambiguously demonstrate that
the relaxation regime can only be explained by taking into
account a recombination term of the photoexcited charge car-
riers, otherwise the system would remain in a photo-doped
nonequilibrium state with almost zero BB indefinitely, which
is not observed experimentally.

Overall, our dynamic simulations based on the solution
of the DDEs equations satisfactorily capture the physics of
the problem [as expected in Fig. 5(d)], although simulations
with realistically slow recombination times were not possible.
Since simulations up to ns timescales are beyond our capabil-
ities we had to choose extremely fast recombination times (a
few fs against ns). Still these calculations demonstrate that the
relaxation of the band structure on a timescale is controlled by
the electron-hole recombination.

IV. CONCLUSIONS

To summarize, the ultrafast dynamics of the SPV in pho-
toexcited K-doped BP has been measured with Tr-ARPES and
compared to numerical simulations of the DDEs. We analysed
the detailed temporal evolution of the SPV and the band
structure near the Fermi level, from a few ps before to a few
hundreds ps after the photoexcitation. Our results demonstrate
that the timescales of the SPV and the spectral shape modica-
tion of the VB are correlated, and support the scenario where
the transient screening of φBB is responsible for the apparent
reopening of the band gap. Our simulations confirm the ultra-
fast transient evolution of φBB (on the timescale of the pump
pulse) and the importance of the electron-hole recombination
processes for the long-time evolution. Without recombination,
the system would remain trapped in a nonequilibrium steady
state with almost no BB.

Our experimental observations and our theoretical ap-
proach are very general and should have implications also for
other Tr-ARPES measurements of 2D semiconducting mate-
rials with significant BB. In particular, care has to be taken
in the interpretation of transient spectral shape modications,
since these may reflect the evolution of the SPV rather than
qualitative changes in the band structure.

ACKNOWLEDGMENTS

This project was supported by the Swiss National Sci-
ence Foundation (SNSF) Grants No. P00P2_170597, No.
200021-165539, and No. 200021-19696 and European Re-
search Council (ERC) Consolidator Grant No. 724103. The
drift-diffusion calculations have been performed on the Beo05
cluster at the University of Fribourg. We are very grateful to
P. Aebi for sharing with us his photoemission setup. Skillful
technical assistance was provided by J. L. Andrey, M. Andrey,
F. Bourqui, B. Hediger, and O. Raetzo. A.P. acknowledges the
Osk. Huttunen Foundation for financial support, and CSC-IT
Center for Science, Finland, for computational resources.

APPENDIX

1. SPV characterization in pristine black phosphorus

Figure 7 displays the Tr-ARPES measurements of BP at
-1 ps pump-probe delay without (red) and with (blue) the

<figure><img src="image_9.png" /></figure>

FIG. 7. Tr-ARPES measurements of BP along the ZZ direction at -1 ps pump-probe delay for a pump fluence of (a) 0 µJ/cm² and (b) of 160 µJ/cm². (c) Corresponding EDCs taken at normal emission. A negative SPV of -13 meV is observed due to the compensation of the upward φBB as depicted in Fig. 3(a).

pump. From this, it is possible to evaluate the magnitude and
the sign of the SPV. We observe that the blue spectrum is
shifted by 13 meV to low KE with respect to the red one.
Since the pristine BP surface is p doped due to the presence of
surface acceptor states, an upward BB of the VB is expected as
illustrated in the middle panel of Fig. 3(a). This is confirmed
by the negative value of the SPV which is compensating the
upward BB.

2. SPV characterization in K-doped black phosphorus

As discussed above in the p-doped BP case, Fig. 8 dis-
plays the Tr-ARPES measurements of K-doped BP at –1 ps
pump-probe delay without (red) and with (blue) the pump.
In contrast to the pristine BP case, we observe that the blue
spectrum is shifted by 200 meV to high KE with respect to
the red one. Due to the strong surface dipole, a long-range
downward BB of both the CB and the VB is expected as
illustrated in the middle panel of Fig. 3(b). This is confirmed
by the large positive value of the SPV which is compensating
this BB.

3. Evidence of spectral shape modification along
the armchair direction

In Fig. 9, we show that the spectral shape modification of
the VB is also visible along the AC direction. As discussed in
<figure><img src="image_10.png" /><figcaption>FIG. 8. Tr-ARPES measurements of K-doped BP along the AC direction at −1 ps pump-probe delay for a pump fluence of (a) 0 µJ/cm² and (b) of 160 µJ/cm². (c) Corresponding EDCs taken four degrees off normal emission. A positive SPV of +200 meV is observed due to the compensation of the downward φBB as depicted in Fig. 3(b).</figcaption></figure>

the case of the ZZ direction in Fig. 2 of the main text, we
observe the same blue/white/red contrast in the difference
plot [see Fig. 9(c)] before and after the photoexcitation of

<figure><img src="image_11.png" /><figcaption>FIG. 10. (a) EDCs of the K-doped BP ARPES spectra taken at normal emission for a few pump-probe delays. These EDCs show both a time dependent energy shift and VB spectral shape evolution. (b) Temporal evolution of the low energy electron cutoff position.</figcaption></figure>

the pump. The main difference along the ZZ direction is the
absence of spectral weight at the bottom of CB, as already
discussed in Ref. [60].

<figure><img src="image_12.png" /><figcaption>FIG. 9. Tr-ARPES measurements of K-doped BP along the AC direction using a pump fluence of 160 µJ/cm² for pump-probe delays of (a) −1 ps and (b) +2 ps. (c) Intensity difference between the panels (b) and (a). Red and blue colors respectively correspond to an increase and a depletion of photoemission intensity. The transient dynamics is the same at the one observed in Fig. 2 of the main text.</figcaption></figure>

4. Raws EDCs and photoemission cutoff dynamics

Figure 10(a) displays EDCs of the K-doped BP sample taken at normal emission showing both a time-dependent energy shift and VB spectral shape evolution. In contrast, the CB is only energy shifted. This is a selection of EDCs we used to produce Figs. 4(a) and 4(b). Figure 10(b) shows the temporal evolution of the low energy electron cutoff position. Its shape is similar to the time-dependent SPV of the VB and CB [see Figs. 4(a) and 4(b)].

5. Time constants as obtained by exponential fits

In Fig. 11, we show a systematic fit of the pertinent time
constants discussed in the main text, for both the VB and the
CB of photoexcited K-doped BP.
<figure><img src="image_13.png" /><figcaption>FIG. 11. Temporal evolution of the SPV at the surface of K-doped BP for (a) the CB and (b) the VB. (c) Temporal evolution of $I_{CB}$ and [(d) and (e)] $I_{VB}$ as defined in the Fig. 4 of the main text. For all these spectra, solid lines are the exponential fit used to extract the characteristic time constants discussed in the text.</figcaption></figure>

## 6. Simulation of the carrier dynamics in photoexcited black phosphorus

### a. Model and parameters

Pristine bulk BP is an intrinsically p-doped semiconductor with a small band gap of $E_g = 0.3$ eV [26,60]. The intrinsic carrier density $p_{intr}$ can be as large as $O(10^{20})$ m⁻³ at T = 200 K [65], but at very low temperature, $p_{intr}$ (and also $n_{intr}$) decays exponentially with decreasing T. At the temperature of the experiment, T = 30 K, the density is as small as 0.016 m⁻³ and thus negligible. In order to electron dope the system, K atoms are deposited on the surface of BP. These K atoms release their s electrons into the sample and create a 2D electron gas near the (positively charged) surface. The spatial separation between the ionized K atoms and the electron gas creates a strong electric field [32]. When electron-hole pairs are created by the pump, their dynamics is affected by this field and by the different mobilities in the VB and CB.

The BB along the z axis is mainly controlled by the carrier dynamics along the z direction. Assuming a homogeneous system in the x-y plane, we consider a one-dimensional drift-diffusion problem along the z axis. The mobilities of the electrons and holes in BP are taken from Ref. [66]. At T = 30 K, $\mu_n = 540$ cm²/Vs, and $\mu_p = 6000$ cm²/Vs, which makes the drift velocity $v = \mu E$ as large as $O(10)$ Å/fs. Thus the simulation needs a timestep of the order of attoseconds to resolve the ultrafast carrier dynamics.

In the simulation, we assume a sample thickness of L = 1000 Å, which is large compared to the atomic scales, but

much smaller than the actual sample thickness (~1 mm). There are two reasons for this choice. On the one hand, the strong electric field requires a real-space grid spacing as small as 0.1 Å. Thus a realistic thickness requires huge matrices of dimension $O(10^8)$ which is computationally very demanding. On the other hand, the penetration depth of the pump pulse is only $O(10)$ nm, and as we will show, the net difference in the electron and hole population is significant only in a narrow region near the surface. Hence, the relevant short-time carrier dynamics takes place near the surface.

Furthermore, we assume that the photoexcited electron-hole pairs, with densities $\delta n(x, t = 0_+)$ and $\delta p(x, t = 0_+)$, are generated instantaneously at time $t = 0_+$, with an exponential profile (see Fig. 12)

$$ \delta n(x, t = 0_+) = \delta p(x, t = 0_+) = \gamma_0 e^{-x/\lambda_{ph}}, \quad (A1) $$

where the surface photoexcited carrier density $\gamma_0 = 5.7 \times 10^{25}$ m⁻³ and the penetration depth $\lambda_{ph} = 500$ Å are chosen in accordance with the experiment. The total density (per area) of photoexcited electron-hole pairs is denoted by $\sigma_p = \int_0^L \delta n(x, t = 0_+)dx = \int_0^L \delta p(x, t = 0_+)dx$. The K atoms are assumed to be deposited near the surface with a distribution that decays exponentially into the bulk, $N_d(x) = n_u \lambda_d e^{-\lambda_d x/\lambda_0}$, where $n_u = \sigma_0/l_0$ is the unit of density and $l_0/\lambda_d \approx 5$ Å represents the doping depth, while $\sigma_0$ is the total K coverage per unit area. In the case of a $\delta$-function distribution of K ions, the thickness of the 2D electron gas is less than 1 Å, which is unrealistic. A more accurate quantum mechanical
<figure><img src="image_14.png" /><figcaption>FIG. 12. The static equilibrium-state solution before the laser pump and the photoexcited carriers distribution at $t = 0_+$.</figcaption></figure>

treatment based on a Schrödinger-Poisson solver would yield a broader charge distribution, and the above exponential profile for $N_d(x)$ is a simple way of mimicking this effect.

As for the carrier densities $n$ and $p$, we neglect the intrinsic densities [67]

$$n_{\text{intr}} = p_{\text{intr}} = \sqrt{\left(\frac{m_e^* m_h^* k_B^2 T^2}{\pi^2 \hbar^4}\right)^{3/2}} e^{-\frac{E_g}{2k_B T}} \\ \approx 1.547 \times 10^{21} T^{3/2} e^{-\frac{1740.7}{T}} \text{m}^{-3}, \quad (\text{A2})$$

where the averaged effective masses are $m_e^* \approx 0.413 m_e$ and $m_h^* \approx 0.335 m_e$ [26] and the band gap is $E_g \approx 0.3$ eV. At temperature $T = 30$ K, $n_{\text{intr}} \approx 0.016 \text{ m}^{-3}$ which is much smaller than the excess carrier density [65] [see also Fig. 3 of Ref. [65]]; even at $T = 200$ K, $n_{\text{intr}} \sim O(10^{20}) \text{ m}^{-3} \ll 10^{25} \text{ m}^{-3}$. We only consider the photoexcited electron-hole densities ($\delta n$, $\delta p$) and the free electron density [$n_0(x)$] resulting from the fully ionized K atoms with surface density $\int_0^L n_0(x) dx = \sigma_0$. Hence, the total electron and hole densities are

$$n(x, t) = n_0(x) + \delta n(x, t) \quad (\text{A3})$$

and

$$p(x, t) = \delta p(x, t). \quad (\text{A4})$$

All the simulation parameters are listed in Table I.

### b. Drift-diffusion equations

The continuity equations for the electrons and holes read

$$-\nabla \cdot \mathbf{J}_n + G_n - R_n = \frac{\partial n}{\partial t}, \quad (\text{A5})$$

and

$$-\nabla \cdot \mathbf{J}_p + G_p - R_p = \frac{\partial p}{\partial t}, \quad (\text{A6})$$

respectively. Here the particle flux densities are

$$\mathbf{J}_n = -D_n \nabla n - n \mu_n \mathbf{E} \quad (\text{A7})$$

for the electrons and

$$\mathbf{J}_p = -D_p \nabla p + p \mu_p \mathbf{E} \quad (\text{A8})$$

for the holes. $\mu_{n,p}$ are the drift mobilities and $D_{n,p}$ the diffusion coefficients which are related to $\mu_{n,p}$ by the Einstein relations $D_{n,p} = \frac{\mu_{n,p} k_B T}{e}$. $G$ and $R$ denote the generation rate and recombination rate, respectively. To illustrate the effect of recombination, we make the simple ad hoc approximation

$$G - R \approx B_r n \delta p. \quad (\text{A9})$$

Here we neglect the thermal generations and recombinations and consider only inter-band recombinations [63].

Combining Eqs. (A3) to (A9), we arrive at the drift-diffusion equations

$$D_n \nabla^2 n + \mu_n (\mathbf{E} \cdot \nabla n + n \nabla \cdot \mathbf{E}) - B_r n \delta p = \frac{\partial \delta n}{\partial t}, \quad (\text{A10})$$

$$D_p \nabla^2 \delta p - \mu_p (\mathbf{E} \cdot \nabla \delta p + \delta p \nabla \cdot \mathbf{E}) - B_r n \delta p = \frac{\partial \delta p}{\partial t}. \quad (\text{A11})$$

Furthermore, the total electric field $\mathbf{E}$ and total charge density satisfy the Poisson equation

$$\nabla \cdot \mathbf{E} = \frac{e}{\kappa} (\delta p - n_0 - \delta n + N_d). \quad (\text{A12})$$

where $\kappa \equiv \epsilon_r \epsilon_0$ with $\epsilon_r$ the relative permittivity and $\epsilon_0$ the vacuum permittivity.

For the stability of the simulation, it is advantageous to reformulate these equations in terms of $\delta n$ and $\delta p$. The free electron distribution $n_0(x)$ is obtained by solving the static (equilibrium) problem before the pump ($\delta n = \delta p = 0$ for $t < 0$)

$$D_n \nabla^2 n_0 + \mu_n (\mathbf{E}_0 \cdot \nabla n_0 + n \nabla \cdot \mathbf{E}_0) = 0, \quad (\text{A13})$$

$$\nabla \cdot \mathbf{E}_0 = \frac{e}{\kappa} (-n_0 + N_d). \quad (\text{A14})$$

The numerical method used for computing the static $n_0$ and $\mathbf{E}_0$ is described in Appendix 6d, where we have to set $\sigma_p = 0$. Equation (A10) can then be rewritten as

$$\frac{\partial \delta n}{\partial t} = D_n \nabla^2 \delta n + \mu_n (\delta \mathbf{E} \cdot \nabla n_0 + \mathbf{E} \cdot \nabla \delta n) \\ + \mu_n (n_0 \nabla \cdot \delta \mathbf{E} + \delta n \nabla \cdot \mathbf{E}) - B_r n \delta p, \quad (\text{A15})$$

where $\mathbf{E} = \mathbf{E}_0 + \delta \mathbf{E}$ and

$$\nabla \cdot \delta \mathbf{E} = \frac{e}{\kappa} (\delta p - \delta n). \quad (\text{A16})$$

For given $\delta n$ and $\delta p$, we calculate $\delta \mathbf{E}$ and $\mathbf{E}$ as

$$\delta E(x) = \frac{e}{\kappa} \int_{0}^{x} [\delta p - \delta n](x') dx', \\ E(x) = \frac{e}{\kappa} \int_{0}^{x} [\delta p - \delta n - n_0 + N_d](x') dx'. \quad (\text{A17})$$

### c. Discretization and semi-implicit Euler method

The density and electric field are continuous functions of the coordinate $x$ and time $t$. Numerically, we discretize $x$ and $t$ as $\{0 = x_0 < x_1 < \dots < x_N = L\}$ and $\{0 = t_0 < t_1 < \dots < t_M = T\}$, and define $\tau_m \equiv t_{m+1} - t_m$. Let $f_i^m$ denote the value of the function $f(x, t)$ at $x_i$ and $t_m$, and $\tau_m \equiv t_{m+1} - t_m$, $f_0^m, f_1^m, \dots, f_N^m$. The first- and second-order derivative with respect to $x$ of a given function $f(x)$ can be represented in a
matrix-vector multiplication form, $f'_i \approx \sum_j D_{1;i,j} f_j$ and $f''_i \approx \sum_j D_{2;i,j} f_j$, where $D_1$ and $D_2$ denote the finite difference matrices for the first and second-order derivative, respectively. Here, we will use the central difference approximation to $D_1$ and $D_2$,

$$f'_i \approx \frac{f_{i+1} - f_{i-1}}{2dx}, \quad f''_i \approx \frac{f_{i+1} - 2f_i + f_{i-1}}{dx^2}. \qquad (A18)$$

One needs to properly treat the zero-flux boundary conditions $J_{n,p}|_{x=0,L} = 0$ [Eqs. (A7) and (A8)] which are equivalent to

$$\nabla n|_{x=0,L} = \nabla p|_{x=0,L} = 0, \qquad (A19)$$

since $E|_{x=0,L} = 0$ due to charge neutrality. By introducing two auxiliary points $x_{-1}$ and $x_{N+1}$, we can write (using $n$ as an example)

$$n'_{0} \approx \frac{n_{1} - n_{-1}}{2dx} = 0, \quad n'_{N} \approx \frac{n_{N+1} - n_{N-1}}{2dx} = 0, \qquad (A20)$$

which gives $n_{-1} = n_1$ and $n_{N+1} = n_{N-1}$. Thus

$$\begin{align}
n''_{0} &\approx \frac{n_{1} - 2n_{0} + n_{-1}}{dx^{2}} \approx \frac{2n_{1} - 2n_{0}}{dx^{2}}, \nonumber \\
n''_{N} &\approx \frac{n_{N+1} - 2n_{N} + n_{N-1}}{dx^{2}} \approx \frac{2n_{N-1} - 2n_{N}}{dx^{2}}. \tag{A21}
\end{align}$$

Hence, we can write the finite-difference matrices for the density as

$$D_1 = \frac{1}{dx} \begin{bmatrix} 0 & 0 & \cdots & 0 \\ -\frac{1}{2} & 0 & \cdots & \frac{1}{2} \\ \vdots & \vdots & \ddots & \vdots \\ -\frac{1}{2} & 0 & \cdots & \frac{1}{2} \\ 0 & 0 & \cdots & 0 \end{bmatrix}, \qquad (A22)$$

and

$$D_2 = \frac{1}{dx^2} \begin{bmatrix} -2 & 2 & \cdots & 1 \\ 1 & -2 & \cdots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ 1 & -2 & \cdots & 1 \\ 2 & 2 & \cdots & -2 \end{bmatrix}. \qquad (A23)$$

The coupled Eqs. (A11), (A12) (A15) and (A16) are solved by using a semi-implicit Euler method:

$$\begin{align}
\frac{\delta \boldsymbol{n}^{m+1} - \delta \boldsymbol{n}^m}{\tau_m} &\approx D_n D_2 \delta \boldsymbol{n}^{m+1} - B_r \delta \boldsymbol{p}_{\text{diag}}^m \boldsymbol{n}^m \nonumber \\
&\quad + \mu_n \left[ \delta \boldsymbol{E}_{\text{diag}}^m D_1 \boldsymbol{n}_0 + \boldsymbol{E}_{\text{diag}}^m D_1 \delta \boldsymbol{n}^{m+1} \right] \nonumber \\
&\quad + \mu_n \left[ (\nabla \cdot \delta \boldsymbol{E})_{\text{diag}}^m \boldsymbol{n}_0 + (\nabla \cdot \boldsymbol{E})_{\text{diag}}^m \delta \boldsymbol{n}^{m+1} \right], \tag{A24}
\end{align}$$

and

$$\begin{align}
\frac{\delta \boldsymbol{p}^{m+1} - \delta \boldsymbol{p}^m}{\tau_m} &\approx D_p D_2 \delta \boldsymbol{p}^{m+1} - B_r \boldsymbol{n}_{\text{diag}}^m \delta \boldsymbol{p}^m \nonumber \\
&\quad - \mu_p \left[ \boldsymbol{E}_{\text{diag}}^m D_1 \delta \boldsymbol{p}^{m+1} + (\nabla \cdot \boldsymbol{E})_{\text{diag}}^m \delta \boldsymbol{p}^{m+1} \right]. \tag{A25}
\end{align}$$

Here $[f]_{\text{diag}}$ represents a diagonal matrix with the elements of the vector $f$. $\delta \boldsymbol{E}_{\text{diag}}^m$ and $\boldsymbol{E}_{\text{diag}}^m$ are calculated from Eq. (A17) using the density at $t_m$ and similarly for $\nabla \delta \boldsymbol{E}^m$ and $\boldsymbol{E}^m$.

Equations (A24) and (A25) can be reduced to the linear equations

$$\begin{align}
&[\mathbb{I} - \tau_m D_n D_2 - \tau_m \mu_n (\boldsymbol{E}_{\text{diag}}^m D_1 + (\nabla \cdot \boldsymbol{E})_{\text{diag}}^m)] \delta \boldsymbol{n}^{m+1} \nonumber \\
&\approx \delta \boldsymbol{n}^m - \tau_m B_r \delta \boldsymbol{p}_{\text{diag}}^m \boldsymbol{n}^m \nonumber \\
&\quad + \tau_m \mu_n [\delta \boldsymbol{E}_{\text{diag}}^m D_1 \boldsymbol{n}_0 + (\nabla \cdot \delta \boldsymbol{E})_{\text{diag}}^m \boldsymbol{n}_0] \tag{A26}
\end{align}$$

and

$$\begin{align}
&[\mathbb{I} - \tau_m D_p D_2 + \tau_m \mu_p (\boldsymbol{E}_{\text{diag}}^m D_1 + (\nabla \cdot \boldsymbol{E})_{\text{diag}}^m)] \delta \boldsymbol{p}^{m+1} \nonumber \\
&\approx [\mathbb{I} - \tau_m B_r \boldsymbol{n}_{\text{diag}}^m] \delta \boldsymbol{p}^m. \tag{A27}
\end{align}$$

The numerical simulation consists of the following steps:

(i) Find the static solution ($\boldsymbol{n}_0$ and $E_0$) before the pump (see Appendix 6d); turn on the pump $\delta \boldsymbol{n}^{m=0} = \delta \boldsymbol{p}^{m=0}$ as in Eq. (A1) at $t = 0_+$ ($m = 0$).

(ii) Update the electric field $\delta \boldsymbol{E}_{\text{diag}}^m$ and $\boldsymbol{E}_{\text{diag}}^m$ by Eq. (A17) using the density at $t_m$; update $(\nabla \cdot \boldsymbol{E})^m$ and $(\nabla \cdot \delta \boldsymbol{E})^m$ by Eqs. (A12) and (A16) using the density at $t_m$.

(iii) Update the densities $\delta \boldsymbol{n}^{m+1}$ and $\delta \boldsymbol{p}^{m+1}$ at $t_{m+1}$ by solving the linear Eqs. (A26) and (A27).

(iv) Repeat steps (ii) and (iii) until $t = t_N$.

For any intermediate time $t_m$, we can calculate the potential profile as

$$\phi(x, t_m) = \phi(L, t_m) + \int_x^L E(x', t_m) dx', \qquad (A28)$$

where the potential reference is chosen as $\phi(L, t_m) = 0$. The band bending profile is given by

$$-e\phi(x, t_m) = -e \int_x^L E(x', t_m) dx', \qquad (A29)$$

and the surface photovoltage is given by $-e\phi(0, t_m)$.

### d. Steady-state solution

Without recombination term in Eqs. (A10) and (A11), the photoexcited carriers drift and diffuse inside the material and eventually reach a steady state. In this steady state, the current is zero everywhere

$$0 = \mathbf{J}_n = -D_n \nabla n - n \mu_n \boldsymbol{E}, \qquad (A30)$$

$$0 = \mathbf{J}_p = -D_p \nabla p + p \mu_p \boldsymbol{E}. \qquad (A31)$$

To find the steady-state solution, we need to solve the coupled Eqs. (A12), (A30) and (A31).

**Dimensionless formulation.** For simplicity, we split $\boldsymbol{E}$ into three components, $E = E_n + E_p + E_d$, with $E'_n = -\frac{e}{\kappa} n$, $E'_p = \frac{e}{\kappa} \delta p$ and $E'_d = \frac{e}{\kappa} N_d$. Furthermore, we introduce the unit $E_u = \frac{\kappa \delta p}{2\kappa}$ for the electric field and $n_u = \frac{\sigma_0}{l_0}$ for the density, and introduce the unitless functions $y$, $h$, and $f$ for $E_n$, $E_p$, and $E_d$, respectively. $E_d(x)$ can be easily calculated from the known $N_d(x)$ as $E_d(x) = E_u f(x)$, with $f(x) = 1 - 2e^{-\lambda_d \frac{x}{l_0}}$. We can thus simplify Eqs. (A30) and (A31) as

$$D_n E_u y''(x) + \mu_n E_u y'(x) E_u [y(x) + h(x) + f(x)] = 0, \qquad (A32)$$

and

$$-D_p E_u h''(x) + \mu_p E_u h'(x) E_u [y(x) + h(x) + f(x)] = 0. \qquad (A33)$$
Now we introduce the unitless coefficient

$$cl_0 = \frac{\mu_n E_u l_0}{D_n} = \frac{\mu_p E_u l_0}{D_p} = \frac{e \sigma_0 e}{k_B T 2\kappa} l_0, \quad (A34)$$

which allows to simplify Eqs. (A32) and (A33) as

$$y''(x) + \frac{cl_0}{l_0} y'(x)[f(x) + y(x) + h(x)] = 0, \quad (A35)$$

and

$$-h''(x) + \frac{cl_0}{l_0} h'(x)[f(x) + y(x) + h(x)] = 0. \quad (A36)$$

We furthermore introduce the two auxiliary functions $P = y + h$, $Q = y - h$ which simplify the implementation of the boundary conditions, as shown below. The sum of Eq. (A35) and Eq. (A36) gives

$$P'' + \frac{cl_0}{l_0} Q'[f + P] = 0, \quad (A37)$$

while the difference of Eq. (A35) and Eq. (A36) gives

$$Q'' + \frac{cl_0}{l_0} P'[f + P] = 0. \quad (A38)$$

In the solution of the coupled second-order nonlinear differential Eqs. (A37) and (A38) on the discretized grid points $0 = x_0 < x_1 < x_2 < \dots < x_{N-1} < x_N = L$, we need to impose the proper boundary conditions for $P$ and $Q$. In the following, we represent a function on the grid points by a vector $f$ with the $i$th element labeled as $f_i = f(x_i)$ for simplicity. The conditions to be imposed are the following:

(a) Charge neutrality requires that the electric field vanishes at the boundaries, i.e., $E|_{x=0} = E|_{x=L} = 0$. Since $E = E_u(P + f)$, we thus have $P_0 = -f_0 = 1.0$ and $P_N = -f_N \approx -1.0$ if $L \gg l_0$.

(b) The electric field component arising from $\delta p(x)$, i.e., $E_p(x) = E_u h(x)$, satisfies $E_{p|_{x=L}} = E_{p|_{x=0}} + \frac{e}{\kappa} \int_0^L \delta p(x') dx' = \frac{e}{\kappa} \sigma_p = E_u \frac{2\sigma_p}{\sigma_0}$. Furthermore, the Gauss law requires that $E_{p|_{x=0}} = -E_{p|_{x=L}}$. Thus we have $h_0 = -\frac{\sigma_p}{\sigma_0} = -h_N$, which gives $Q_0 = y_0 - h_0 = -f_0 - 2h_0$ and $Q_N = y_N - h_N = -f_N - 2h_N$.

(c) $P''_0 = Q''_0 = P''_N = Q''_N = 0$, since $E \sim f + P$ is zero at the boundaries.

**Newton-Raphson iterative method.** Solving Eqs. (A37) and (A38) for $P$ and $Q$ is equivalent to solving a multivariable root problem $F(X) = 0$. The Newton-Raphson iterative method [68] is based on the first order Taylor expansion of $F(X)$ as $F(X + \delta X) \approx F(X) + \mathcal{J} \cdot \delta X$, where $\mathcal{J}$ denotes the Jacobian matrix with elements $\mathcal{J}_{i,j} = \frac{\partial F_i}{\partial X_j}$. By setting $F(X + \delta X) = 0$, we obtain a set of linear equations for $\delta X$:

$$\mathcal{J} \cdot \delta X = -F. \quad (A39)$$

$\delta X$ is the correction to $X$ and we update the solution as

$$X_{\text{new}} = X_{\text{old}} + \delta X \quad (A40)$$

until convergence.

To utilize the Newton-Raphson method, we need to combine Eqs. (A37) and (A38) into a matrix form. Since $P_{0,N}$, $Q_{0,N}$, and $Q''_{0,N}$ are constants given by the boundary conditions, we only need to define the finite-difference matrices

$D_{1,2}$ on the internal grid points $\{x_1, x_2, \dots, x_{N-1}\}$ with

$$D_1 = \frac{1}{dx} \begin{bmatrix} 0 & \frac{1}{2} & & & \\ -\frac{1}{2} & 0 & \frac{1}{2} & & \\ & \ddots & \ddots & \ddots & \\ & & -\frac{1}{2} & 0 & \frac{1}{2} \\ & & & -\frac{1}{2} & 0 \end{bmatrix}, \quad (A41)$$

and

$$D_2 = \frac{1}{dx^2} \begin{bmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & \ddots & \ddots & \ddots & \\ & & 1 & -2 & 1 \\ & & & 1 & -2 \end{bmatrix}. \quad (A42)$$

Let us introduce the vectors

$$\vec{P}^T = [P_1, P_2, \dots, P_{N-1}], \quad (A43)$$

$$\vec{Q}^T = [Q_1, Q_2, \dots, Q_{N-1}],$$

and

$$\begin{align*}
\vec{v}_1^T &= \frac{1}{dx^2} [P_0, 0, \dots, 0, P_N], \\
\vec{v}_2^T &= \frac{1}{dx} \left[ -\frac{1}{2} P_0, 0, \dots, 0, \frac{1}{2} P_N \right], \\
\vec{v}_3^T &= \frac{1}{dx^2} [Q_0, 0, \dots, 0, Q_N], \\
\vec{v}_4^T &= \frac{1}{dx} \left[ -\frac{1}{2} Q_0, 0, \dots, 0, \frac{1}{2} Q_N \right]. \tag{A44}
\end{align*}$$

$(\forall i, \dim \vec{v}_i^T = N - 1)$. The first and second derivatives are given by

$$\begin{align*}
\vec{P}'' &= D_2 \vec{P} + \vec{v}_1, \\
\vec{P}' &= D_1 \vec{P} + \vec{v}_2, \\
\vec{Q}'' &= D_2 \vec{Q} + \vec{v}_3, \\
\vec{Q}' &= D_1 \vec{Q} + \vec{v}_4.
\end{align*} \quad (A45)$$

Now we define the vector $F$ from Eqs. (A37) and (A38) as

$$\begin{align}
F &= \begin{pmatrix} F_1 \\ F_2 \end{pmatrix} = \begin{pmatrix} \vec{P}'' + c[\vec{f} + \vec{P}]_{\text{diag}} \vec{Q}' \\ \vec{Q}'' + c[\vec{f} + \vec{P}]_{\text{diag}} \vec{P}' \end{pmatrix} \nonumber \\
&= \begin{bmatrix} D_2 \vec{P} + \vec{v}_1 + c[\vec{f} + \vec{P}]_{\text{diag}} (D_1 \vec{Q} + \vec{v}_4) \\ D_2 \vec{Q} + \vec{v}_3 + c[\vec{f} + \vec{P}]_{\text{diag}} (D_1 \vec{P} + \vec{v}_2) \end{bmatrix}. \tag{A46}
\end{align}$$

The Jacobian matrix reads

$$\begin{align}
\mathcal{J} &= \begin{bmatrix} \frac{\partial F_1}{\partial \vec{P}} & \frac{\partial F_1}{\partial \vec{Q}} \\ \frac{\partial F_2}{\partial \vec{P}} & \frac{\partial F_2}{\partial \vec{Q}} \end{bmatrix} \nonumber \\
&= \begin{bmatrix} D_2 + c(D_1 \vec{Q} + \vec{v}_4)_{\text{diag}} & c[\vec{f} + \vec{P}]_{\text{diag}} D_1 \\ c(D_1 \vec{P} + \vec{v}_2)_{\text{diag}} + c[\vec{f} + \vec{P}]_{\text{diag}} D_1 & D_2 \end{bmatrix}. \tag{A47}
\end{align}$$

According to Eq. (A40), $\vec{P}$ and $\vec{Q}$ are updated as

$$\begin{pmatrix} \vec{P} \\ \vec{Q} \end{pmatrix}_{\text{new}} = \begin{pmatrix} \vec{P} \\ \vec{Q} \end{pmatrix}_{\text{old}} + \begin{pmatrix} \delta \vec{P} \\ \delta \vec{Q} \end{pmatrix} \quad (A48)$$
after solving the linear equation

$$ \mathcal{J} \begin{pmatrix} \delta \vec{P} \\ \delta \vec{Q} \end{pmatrix} = -F. \qquad (A49) $$

### e. Static solution

The static equilibrium state before the pump can be obtained in a similar manner as described in Appendix 6d, by

setting the parameter $\sigma_p$ to zero. In Fig. 12, the distribution of $N_d(x)$ is shown as the red line, and the distribution of free electrons from ionized K is shown as the black dashed line. The inset shows the difference between $N_d(x)$ and $n_0(x)$ near the surface region. The sharp surface peak and the small broad peak in $N_d - n_0$ at $x \approx 40$ Å imply a spatial separation between the K ions and the free electrons, and this gives rise to a strong internal electric field in this region, as shown by the green line ($t=0$) in Fig. 6(c).

[1] A. Chaves, J. G. Azadani, H. Alsalman, D. R. da Costa, R. Frisenda, A. J. Chaves, S. H. Song, Y. D. Kim, D. He, J. Zhou et al., Npj 2D Mater. Appl. 4, 29 (2020).

[2] G. Ak and K. Novoselov, Nat. Mater. 6, 183 (2007).

[3] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Phys. Rev. Lett. 108, 155501 (2012).

[4] L. Li, S.-z. Lu, J. Pan, Z. Qin, Y.-q. Wang, Y. Wang, G.-y. Cao, S. Du, and H.-J. Gao, Adv. Mater. 26, 4820 (2014).

[5] J. P. Ji, X. F. Song, J. Z. Liu, Z. Yan, C. X. Huo, S. L. Zhang, M. Su, L. Liao, W. H. Wang, Z. H. Ni, Y. F. Hao, and H. B. Zeng, Nat. Commun. 7, 13352 (2016).

[6] Z. Zhu, X. Cai, S. Yi, J. Chen, Y. Dai, C. Niu, Z. Guo, M. Xie, F. Liu, J.-H. Cho et al., Phys. Rev. Lett. 119, 106101 (2017).

[7] V. Kochat, A. Samanta, Y. Zhang, S. Bhowmick, P. Manimunda, S. A. S. Asif, A. S. Stender, R. Vajtai, A. K. Singh, C. S. Tiwary et al., Sci. Adv. 4, e1701373 (2018).

[8] J. Deng, B. Xia, X. Ma, H. Chen, H. Shan, X. Zhai, B. Li, A. Zhao, Y. Xu, W. Duan et al., Nat. Mater. 17, 1081 (2018).

[9] R. Wu, I. K. Drozdov, S. Eltinge, P. Zahl, S. Ismail-Beigi, I. Božović, and A. Gozar, Nat. Nanotechnol. 14, 44 (2019).

[10] G. Bihlmayer, J. Sassmannshausen, A. Kubetzka, S. Blügel, K. von Bergmann, and R. Wiesendanger, Phys. Rev. Lett. 124, 126401 (2020).

[11] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nat. Rev. Mater. 2, 17033 (2017).

[12] C. Elias, P. Valvin, T. Pelini, A. Summerfield, C. J. Mellor, T. S. Cheng, L. Eaves, C. T. Foxon, P. H. Beton, S. V. Novikov et al., Nat. Commun. 10, 2639 (2019).

[13] H. Zhang, M. Holbrook, F. Cheng, H. Nam, M. Liu, C.-R. Pan, D. West, S. Zhang, M.-Y. Chou, and C.-K. Shih, ACS Nano 15, 2497 (2021).

[14] L. Lichtenstein, M. Heyde, S. Ulrich, N. Nilius, and H.-J. Freund, J. Phys.: Condens. Matter 24, 354010 (2012).

[15] G. Kremer, J. C. Alvarez Quiceno, S. Lisi, T. Pierron, C. González, M. Sicot, B. Kierren, D. Malterre, J. E. Rault, P. Le Fèvre, F. Bertran, Y. J. Dappe, J. Coraux, P. Pochet, and Y. Fagot-Revurat, ACS Nano 13, 4720 (2019).

[16] G. Kremer, J. C. Alvarez Quiceno, T. Pierron, C. González, M. Sicot, B. Kierren, L. Moreau, J. E. Rault, P. Le Fèvre, F. Bertran, Y. J. Dappe, J. Coraux, P. Pochet, and Y. Fagot-Revurat, 2D Mater. 8, 035021 (2021).

[17] H. Guo, N. Lu, L. Wang, X. Wu, and X. C. Zeng, J. Phys. Chem. C 118, 7242 (2014).

[18] A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Phys. Rev. Lett. 112, 176801 (2014).

[19] A. Ramasubramaniam, D. Naveh, and E. Towe, Phys. Rev. B 84, 205325 (2011).

[20] M. Kang, B. Kim, S. H. Ryu, S. W. Jung, J. Kim, L. Moreschini, C. Jozwiak, E. Rotenberg, A. Bostwick, and K. S. Kim, Nano Lett. 17, 1610 (2017).

[21] Y. Gong, Z. Liu, A. R. Lupini, G. Shi, J. Lin, S. Najmaei, Z. Lin, A. L. Elías, A. Berkdemir, G. You et al., Nano Lett. 14, 442 (2014).

[22] Q. Ma, M. Isarraraz, C. S. Wang, E. Preciado, V. Klee, S. Bobek, K. Yamaguchi, E. Li, P. M. Odenthal, A. Nguyen et al., ACS Nano 8, 4672 (2014).

[23] R. A. B. Villaos, C. P. Crisostomo, Z.-Q. Huang, S.-M. Huang, A. A. B. Padama, M. A. Albao, H. Lin, and F.-C. Chuang, Npj 2D Mater. Appl. 3, 2 (2019).

[24] M.-K. Lin, Rovi Angelo B. Villaos, J. A. Hlevyack, P. Chen, R.-Y. Liu, C.-H. Hsu, J. Avila, S.-K. Mo, F.-C. Chuang, and T.-C. Chiang, Phys. Rev. Lett. 124, 036402 (2020).

[25] K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, 943 (2016).

[26] A. Morita, Appl. Phys. A 39, 227 (1986).

[27] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nat. Nanotechnol. 9, 372 (2014).

[28] X. Ling, H. Wang, S. Huang, F. Xia, and M. S. Dresselhaus, Proc. Natl. Acad. Sci. U.S.A. 112, 4523 (2015).

[29] X. Mu, J. Wang, and M. Sun, Mater. Today Phys. 8, 92 (2015).

[30] K. Dolui and S. Y. Quek, Sci. Rep. 5, 11699 (2015).

[31] Y. Liu, Z. Qiu, A. Carvalho, Y. Bao, H. Xu, S. J. Tan, W. Liu, A. Castro Neto, K. P. Loh, and J. Lu, Nano Lett. 17, 1970 (2017).

[32] J. Kim, S. S. Baik, S. H. Ryu, Y. Sohn, S. Park, B.-G. Park, J. Denlinger, Y. Yi, H. J. Choi, and K. S. Kim, Science 349, 723 (2015).

[33] S.-W. Kim, H. Jung, H.-J. Kim, J.-H. Choi, S.-H. Wei, and J.-H. Cho, Phys. Rev. B 96, 075416 (2017).

[34] Z. Zhang and J. T. Yates Jr, Chem. Rev. 112, 5520 (2012).

[35] B. Kiraly, E. J. Knol, K. Volckaert, D. Biswas, A. N. Rudenko, D. A. Prishchenko, V. G. Mazurenko, M. I. Katsnelson, P. Hofmann, D. Wegner et al., Phys. Rev. Lett. 123, 216403 (2019).

[36] Z. Chen, J. Dong, C. Giorgetti, E. Papalazarou, M. Marsi, Z. Zhang, B. Tian, Q. Ma, Y. Cheng, J.-P. Rueff et al., 2D Mater. 7, 035027 (2020).

[37] H. Hedayat, A. Ceraso, G. Soavi, S. Akhavan, A. Cadore, C. Dallera, G. Cerullo, A. C. Ferrari, and E. Carpene, 2D Mater. 8, 025020 (2021).

[38] U. Bovensiepen and P. S. Kirchmann, Laser Photon. Rev. 6, 589 (2012).

[39] C. L. Smallwood, R. A. Kaindl, and A. Lanzara, Europhys. Lett. 115, 27001 (2016).
[40] W. Widdra, D. Bröcker, T. Gießel, I. Hertel, W. Krüger, A. Liero, F. Noack, V. Petrov, D. Pop, P. Schmidt *et al.*, Surf. Sci. **543**, 87 (2003).

[41] S. Tokudomi, J. Azuma, K. Takahashi, and M. Kamada, J. Phys. Soc. Jpn. **77**, 014711 (2008).

[42] S.-i. Tanaka, J. Electr. Spectrosc. **185**, 152 (2012).

[43] S.-L. Yang, J. A. Sobota, P. S. Kirchmann, and Z.-X. Shen, Appl. Phys. A **116**, 85 (2014).

[44] S. Lange, P. Schmidt, and T. Nilges, Inorg. Chem. **46**, 4028 (2007).

[45] J. Faure, J. Mauchain, E. Papalazarou, W. Yan, J. Pinon, M. Marsi, and L. Perfetti, Rev. Sci. Instrum. **83**, 043109 (2012).

[46] G. Kresse and J. Hafner, Phys. Rev. B **47**, 558 (1993).

[47] G. Kresse and J. Hafner, Phys. Rev. B **49**, 14251 (1994).

[48] G. Kresse and J. Furthmüller, Comput. Mater. Sci. **6**, 15 (1996).

[49] G. Kresse and J. Furthmüller, Phys. Rev. B **54**, 11169 (1996).

[50] P. E. Blöchl, Phys. Rev. B **50**, 17953 (1994).

[51] G. Kresse and D. Joubert, Phys. Rev. B **59**, 1758 (1999).

[52] J. Sun, A. Ruzsinszky, and J. P. Perdew, Phys. Rev. Lett. **115**, 036402 (2015).

[53] A. Brown and S. Rundqvist, Acta Crystallogr. **19**, 684 (1965).

[54] K. Momma and F. Izumi, J. Appl. Crystallogr. **44**, 1272 (2011).

[55] T. Takahashi, N. Gunasekara, H. Ohsawa, H. Ishii, T. Kinoshita, S. Suzuki, T. Sagawa, H. Kato, T. Miyahara, and I. Shirotani, Phys. Rev. B **33**, 4324 (1986).

[56] E. Golias, M. Krivenkov, and J. Sánchez-Barriga, Phys. Rev. B **93**, 075207 (2016).

[57] B. Kiraly, N. Hauptmann, A. N. Rudenko, M. I. Katsnelson, and A. A. Khajetoorians, Nano Lett. **17**, 3607 (2017).

[58] Z. Qiu, H. Fang, A. Carvalho, A. Rodin, Y. Liu, S. J. Tan, M. Telychko, P. Lv, J. Su, Y. Wang *et al.*, Nano Lett. **17**, 6935 (2017).

[59] E. Papalazarou, L. Khalil, M. Caputo, L. Perfetti, N. Nilforoushan, H. Deng, Z. Chen, S. Zhao, A. Taleb-Ibrahimi, M. Konczykowski *et al.*, Phys. Rev. Materials **2**, 104202 (2018).

[60] Z. Chen, J. Dong, E. Papalazarou, M. Marsi, C. Giorgetti, Z. Zhang, B. Tian, J.-P. Rueff, A. Taleb-Ibrahimi, and L. Perfetti, Nano Lett. **19**, 488 (2018).

[61] M. P. Seah and W. Dench, Surf. Interface Anal. **1**, 2 (1979).

[62] J. A. Sobota, S. Yang, J. G. Analytis, Y. L. Chen, I. R. Fisher, P. S. Kirchmann, and Z.-X. Shen, Phys. Rev. Lett. **108**, 117403 (2012).

[63] S. Li, Semiconductor Physical Electronics (Springer-Verlag New York, 2016).

[64] We caution here that the sub-fs hopping timescale is obtained by assuming that the experimental mobilities are still valid for such a very high electric field. In practice, this might not be the case, resulting potentially a longer rise-time. As long as it remains shorter than the pump duration, our conclusions are not qualitatively affected.

[65] P. Bhaskar, A. W. Achtstein, M. J. W. Vermeulen, and L. D. A. Siebbeles, J. Phys. Chem. C **120**, 13836 (2016).

[66] Y. Akahama, S. Endo, and S.-i. Narita, JPSJ **52**, 2148 (1983).

[67] J. H. Davies, *The Physics of Low-Dimensional Semiconductors* (Cambridge University Press, Cambridge, UK, 1997).

[68] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes: The Art of Scientific Computing*, 3rd ed. (Cambridge University Press, Cambridge, UK, 2007).