Pseudogap in a crystalline insulator doped by disordered metals

Key to our understanding of how electrons behave in crystalline solids is the band structure that connects the energy of electron waves to their wavenumber. Even in phases of matter with only short-range order (liquid or amorphous solid), the coherent part of electron waves still has a band structure. Theoretical models for the band structure of liquid metals were formulated more than five decades ago1–15, but, so far, band-structure renormalization and the pseudogap induced by resonance scattering have remained unobserved. Here we report the observation of the unusual band structure at the interface of a crystalline insulator (black phosphorus) and disordered dopants (alkali metals). We find that a conventional parabolic band structure of free electrons bends back towards zero wavenumber with a pseudogap of 30–240 millielectronvolts from the Fermi level. This is wavenumber renormalization caused by resonance scattering, leading to the formation of quasi-bound states in the scattering potential of alkali-metal ions. The depth of this potential tuned by different kinds of disordered alkali metal (sodium, potassium, rubidium and caesium) allows the classification of the pseudogap of p-wave and d-wave resonance. Our results may provide a clue to the puzzling spectrum of various crystalline insulators doped by disordered dopants16–20, such as the waterfall dispersion observed in copper oxides. A back-bending band structure and an emerging pseudogap are observed at the interface between a crystalline solid (black phosphorus) and disordered alkali-metal dopants.

In the 1960s, there was a series of pioneering theoretical studies on the band structure of liquid metals [1][2][3][4][5][6][7][8][9] . Like glassy materials have a well defined complex refractive index, electron waves under the influence of multiple scattering acquire a complex wavenumber shift Δk whose magnitude is sizable at resonance. As shown in Fig. 1a, the band structure of free electrons (energy E ≈ k 2 ) is distorted by the real part of Δk (Re(Δk)) in an unusual sinusoidal form that has no counterpart in crystalline solids. The corresponding imaginary part of Δk (Im(Δk)) represents a spread of k related to the formation of quasi-bound states (QBS) around liquid ions. In the density of states (Fig. 1a, inset), a local minimum is formed at a resonance energy (E r ) 5,[9][10][11][12] , which is the 'pseudogap' coined by Mott 13,14 . This backward-bending band dispersion with pseudogap should be observable by angle-resolved photoemission spectroscopy (ARPES). ARPES has been used to study the electronic structure of a melting lead monolayer on Cu(111) 21 and Si(111) 22 . However, the characteristic k renormalizations and pseudogap induced by resonance scattering depicted in Fig. 1a, despite their fundamental importance, have remained unobserved experimentally.
A key idea in our experiments is to study the interface of disordered dopants (alkali metals) and a crystalline insulator (black phosphorus) [23][24][25][26] . Black phosphorus is a layered material, in which the honeycomb lattice of phosphorus atoms is modulated to form an array of zigzag ridges and valleys (Fig. 1b). The distribution of alkali metals on black phosphorus has been reported to show the radial (and anisotropic) structure factor 25 , which is reproduced by structural simulations in Fig. 1c (see Extended Data Fig. 1 for key assumptions). This radial shape is the key feature of liquid or glassy phases 27,28 , and reflects the presence of the mean interatomic distance (or short-range order) generic to dopants randomly distributed under the repulsive interaction. Each alkali metal donates its valence electron to black phosphorus, and the doped electrons are mostly populated in the topmost layer 23,24 . Then, doped electrons are subject to multiple scattering by the potential of dopant ions with only short-range order, which is the situation required to have resonance effects in theoretical models [1][2][3][4][5][6][7][8][9][10][11][12] (Fig. 1a). Furthermore, the low diffusion barrier of alkali metals on black phosphorus 26 allows us to systematically trace the evolution of ARPES data by varying the density of dopants at low temperature.

Back-bending dispersion and pseudogap
In the simple picture of surface doping that has been widely accepted so far, the conduction bands in the topmost layer of black phosphorus adjacent to alkali-metal ions shift below the Fermi energy E F . Figure 1d shows constant-energy contours at E F (the Fermi surface) 23,24 expected for surface-doped black phosphorus at an electron density n e of 1.0 × 10 14 cm −2 . There is a large, oval-shape contour centred at the Γ point (labelled C1) with a pair of small pockets separated along the y axis (labelled C2). The energy dispersion of the C1 and C2 bands is shown in Fig. 2a-c along three high-symmetry directions (k y , k s and k x as indicated in Fig. 1d), respectively. Owing to the crystal symmetry of black phosphorus, the band dispersion of C1 changes from quadratic in k y (Fig. 2a) to linear in k x (Fig. 2c) 23 . A set of k at E F (or k F ) in every in-plane direction constitutes the oval-shape Fermi contour of C1 in Fig. 1d. c, Fast Fourier transform image taken from the distribution of dopants in b to show their reciprocal lattice. The blue oval is a set of half the reciprocal lattice vectors in every in-plane k direction, which corresponds to k r in the model of liquid metals. d, Constant-energy contours at E F expected for black phosphorus whose surface is doped at n e = 1.0 × 10 14 cm −2 and plotted over the surface Brillouin zone. The grey arrows indicate three high-symmetry k cuts along which a series of ARPES data shown in Fig. 2 are taken.  bottom right of each panel. In data for Na taken along k y (Fig. 2d), there is a clear signature of C2 bands at ±0.6 Å −1 , as expected from Fig. 2a. In contrast, a striking deviation from that expected in Fig. 2a is observed for the C1 band. Unlike the sharp Fermi-Dirac cut-off of C2 bands that can be used as an internal reference of E F , the C1 band shows a gaplike feature with a magnitude of at least 0.2 eV from E F . More noteworthy is that the band dispersion of C1 bends back at around −0.2 eV towards k = 0 or the Γ point with a weak intensity. This backward-bending dispersion of the C1 band can be seen more clearly in k s (Fig. 2e) with a diminishing intensity towards E F , and intriguingly even reaches close to k = 0 in k x (Fig. 2f). However, a difference from the Na case is observed for black phosphorus doped by K, Rb and Cs. In data taken along k y (Fig. 2g, j, m), there is a common pair of C2 bands crossing E F , as observed in Fig. 2d. However, the C1 band shows a well defined energy gap with no sign of the backward-bending band dispersion, as the ARPES intensity diminishes more abruptly towards E F . The magnitude of the energy gaps is about 65-70 meV for K and Rb, and about 30 meV for Cs, which are clearly smaller than 200 meV in the Na case (Fig. 2d). A series of ARPES data taken along k s (Fig. 2h, k, n) and along k x (Fig. 2i, l, o) consistently shows the well defined energy gap with nearly the same magnitudes within the range of ±10 meV, which indicates the isotropy of the energy gap (see Extended Data Fig. 3 for constant-energy maps). This signature of a pseudogap is observed without exception, regardless of the photon energy and samples (Extended Data Fig. 4).

Theoretical model of liquid metals
The backward-bending band dispersion and gap-like features can be naturally explained by the characteristic feature in the band structure predicted in the theory of liquid metals 1-12 . We employ a single-site (structure-independent) model 12 for multiple scattering (Methods), in which each ion is assumed to be a spherical step potential U s whose depth is V 0 and radius is r s (Fig. 3a). This may be a crude approximation, but, as will be shown below, even this simple model is able to capture the essential feature in the band structure of liquid metals. For the partial wave of l ≠ 0 (ref. 9 ), the sum of ionic U s and centrifugal U l forms a potential well around each ion. Electron waves scattered by this effective potential acquire a phase shift δ l as a function of k, such as that in Fig. 3b. If δ l rises through π/2 in the lth partial wave, the magnitude of the scattering amplitude f l peaks at the resonance k, or k r , which is the characteristic feature of resonance 9 (Fig. 3c). The phase difference between incoming and outgoing waves outside of r s is 2δ l , that is, π at k r . Their destructive interference spatially localizes the electron waves in the potential well around each ion, which is the formation of QBS (Fig. 3d).
The variation of complex kf l with δ l is restricted by the unitary relation 29 that can be readily understood by drawing the unitary circle in Fig. 3e. This explains that Re(kf l ) and Im(kf l ) take the simple form of sin(2δ l ) and sin 2 (δ l ), respectively (Fig. 3f). In the theory of liquid metals 1-9 (Methods), a complex Δk of electron waves acquired by the effect of multiple scattering can be approximately written in terms of the predominant f l of the partial wave at resonance as: where n d is the dopant density. It follows that Re(Δk) and Im(Δk) should be directly proportional to sin(2δ l )/k 2 and sin 2 (δ l )/k 2 , respectively. As depicted in Fig. 3g, the typical quadratic band dispersion of free electrons (dotted line) is distorted by Re(Δk) in the sinusoidal form (bold line). The Fourier transform of electrons with complex Δk peaks at k + Re(Δk), but spreads over the range of ±Im(Δk) (grey area in Fig. 3g). This spread of k corresponds to the spatial localization of electrons or QBS in the potential well around each ion (Fig. 3a, d). This is accompanied by changes in the density of states 5,[9][10][11][12] that are roughly proportional to the energy derivative of Re(Δk) (Fig. 3h).
The negative slope of band dispersion near E r (shown in red) forms a . U s is assumed as a step potential whose depth is V 0 and width is r s . b, Phase shift δ l simulated in the form of sigmoid that rises from 0 to π through π/2 at k r . c, Magnitude of the scattering amplitude f l calculated from δ l in b by sin(δ l )/k. The black line is a trace of f l at resonance with k r . d, k series of a probability density calculated for spherical waves in the effective potential in a (Methods). local minimum in the density of states, which is the pseudogap 13,14 . If a residual density of states in the pseudogap is below about one-third of that of free electrons 14 , electronic states in the pseudogap are localized as shown in Fig. 3d, so that conductivity vanishes in the sense of Anderson localization 15 .

Partial-wave analysis and simulations
The resonance scattering, which is characterized by δ l that rises through π/2 ( Fig. 3b), occurs in one of the partial waves (l ≠ 0) for a given k r (Extended Data Fig. 5) 9 . Which partial wave is at resonance depends on the depth of the scattering potential V 0 (Fig. 3a). As shown in Fig. 4a, there are three distinct V 0 ranges, each of which has a partial wave of l = 1 (p-wave), l = 2 (d-wave) or l = 3 (f-wave) at resonance. The variation of δ l with k is calculated for three different V 0 values to have p-wave, d-wave and f-wave resonance at the same k r (dotted line in Fig. 4a), as compared in Fig. 4b. The rise of δ l at d-wave resonance is steeper (narrower in k) than that at p-wave resonance. The corresponding density of states is calculated based on thin-slab approximations (Methods) 12 , which is shown in Fig. 4c. Indeed, the p-wave pseudogap is greater in magnitude than the d-wave pseudogap. Another remarkable difference in Fig. 4b is that δ l at p-wave resonance rises through π/2 but turns back without reaching close to π (incomplete resonance, Extended Data Fig. 6). This makes the p-wave pseudogap less clear with a residual density of states, compared with the d-wave pseudogap (Fig. 4c). The band structure renormalized by resonance scattering and the calculated density of states are taken to simulate a series of ARPES spectra (Methods), as shown in Fig. 4e-j. The simulated APRES spectra for p-wave resonance in Fig. 4e-g reproduce the backward-bending dispersion of Na-doped black phosphorus with a diminishing intensity towards E F (Fig. 2d-f). This could have been observed owing to the finite density of states in the p-wave pseudogap (red curve in Fig. 4c). On the other hand, the sharp dip of the d-wave pseudogap (yellow curve in Fig. 4c) renders the simulated ARPES spectra in Fig. 4h-j to exhibit a clear gaplike feature with little spectral weight on the backward-bending part of the band dispersion, as observed for K, Rb and Cs in Fig. 2g-o (see Extended Data Fig. 7 for a fuller set of simulations). Therefore, our spectral simulations ( Fig. 4e-j) not only collectively reproduce key aspects of experimental observations ( Fig. 2d-o) but also naturally explain the difference between the Na case and the K, Rb and Cs cases by the orbital character of the pseudogap, p-wave or d-wave. This is further supported by the δ l of partial waves calculated with the realistic screened scattering potential of alkali ions 30 (Extended Data Fig. 8), where K, Rb and Cs show d-wave resonance, whereas Na shows p-wave resonance, as exactly observed in our experimental results.
For the quantitative analysis on the magnitude of the pseudogap, the ARPES spectral weight, which is proportional to the density of states, is plotted as a function of energy in Fig. 4d for Na and K (Methods). They commonly show the diminishing spectral weight towards E F over different energy scales. The magnitude of the pseudogap is defined by the energy at which the spectral weight drops by half relative to E F . A fit to this pseudogap with that in the calculated density of states (black curves in Fig. 4d, Methods) yields 235 ± 53 meV for Na, 65 ± 23 meV for K, 65 ± 17 meV for Rb and 33 ± 12 meV for Cs. They are in the range of p-wave (200−250 meV) and d-wave (30−80 meV) pseudogaps at their corresponding k r , respectively (Fig. 5a).

Doping dependence and phase diagram
The magnitude of the pseudogap is plotted in Fig. 5b as a function of n d for different kinds of alkali metal (see Extended Data Fig. 9 for a doping series of raw ARPES data). In the wide range of n d shown in grey, the pseudogap was persistently observed at E F . k r is inversely proportional to the mean interatomic distance of dopants, and its inverse square corresponds to n d , that is, n d ∝ k r 2 . For an insulator doped by monovalent dopants, n d is equal to n e , which is proportional to the square of k F , that is, n e ∝ k F 2 . Therefore, the rate of increase in k r as a function of n d should be the same as that in k F , which means that k r follows k F (E r is always located at E F ) as n d increases. Once this subtle balance between the rate of increase in k F and k r is broken by the C2 band crossing E F , n d is no longer equivalent to the n e of C1, and the magnitude of the pseudogap reduces as observed in Fig. 5b. This can also naturally Energy (eV) Article explain the absence of a pseudogap in C2 by the deviation of k F from k r (0.2-0.3 Å -1 ), which is clearly greater than half the Δk of the pseudogap (at most 0.08 Å −1 for Na).
It is also remarkable that the pseudogap observed in the wide range of n d is nearly isotropic in magnitude. The pseudogap can be isotropic if the spatial distribution of dopants is made to locate a set of k r (shown by the blue oval in Fig. 1c) along the Fermi surface of C1 (Fig. 1d), which is supported by the radial and anisotropic structure factor reported in a previous scanning tunnelling microscopy (STM) study 25 at n d = 1.8 × 10 13 cm −2 . This is a natural extension of electronic interactions between dopants mediated by substrates 25,27 , with extra energy gains arising from the opening of a pseudogap in the band structure, which is connected to the concept of glassy charge-density waves 31 .

Discussion
It is difficult to distinguish the ARPES spectra of the d-wave pseudogap ( Fig. 4h-j) from those of actual energy gaps caused by any kind of symmetry breaking in a crystalline solid unless there are other signatures, such as band folding 32,33 and Bogoliubov bands 34 . Our systematic observation of both a p-wave pseudogap and a d-wave pseudogap allows us to clearly identify their origin as k renormalizations induced by resonance scattering, which were predicted in theoretical models for the band structure of liquid metals. This renormalized band structure is a common characteristic regardless of dimensionality, but it may affect the transport of QBS as in Anderson localization 15 , depending on how well confined the QBS are by the potential of ions. It should be noted again that the backward-bending band dispersion of Na-doped black phosphorus in Fig. 2f even reaches close to k = 0 (|Δk| is comparable to the magnitude of k), as predicted in the self-consistent model proposed by Anderson and McMillan 5 . A potential pairing instability in this band structure might be an interesting topic for future studies.
Our results may provide a clue to the puzzling spectrum of other crystalline insulators doped by disordered dopants if their correlation function peaks at a mean interatomic distance. A possible example is the waterfall dispersion widely observed and debated in the study of copper oxides [16][17][18] . Extended Data Fig. 10a shows ARPES spectra taken from Bi 1.5 Pb 0.5 Sr 2 CaCu 2 O 8+δ (critical temperature T c = 88 K), where there is a clear signature of waterfall dispersion. n d is estimated from T c (ref. 35 ), the pseudogap (32 meV) 36 and the Luttinger count to be 1.2 × 10 14 −1.4 × 10 14 cm −2 , corresponding to k r = 0.29 ± 0.01 Å −1 . Our spectral simulation for resonance scattering at k r = 0.29 Å −1 in Extended Data Fig. 10b (Methods) reproduces the waterfall dispersion remarkably well. At this stage, given that the dopants in copper oxides are likely to be randomly distributed, we cannot rule out that the effect of electron correlation plays a key role in the waterfall dispersion. In any case, however, the physics revealed here is worth considering to account for the puzzling spectrum of copper oxides [16][17][18] and other crystalline insulators doped by disordered dopants 19,20 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-021-03683-0. b, Magnitude of pseudogap estimated from ARPES data of black phosphorus doped by Na, K, Rb and Cs (Extended Data Fig. 9), which is plotted as a function of n d over the range of 0−1 ML. The grey area represents the pseudogap phase. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

ARPES experiments and analysis
We conducted ARPES measurements at Beamline 7.0.2 (MAESTRO) at the Advanced Light Source. The microARPES end-station is equipped with a hemispherical spectrometer. We collected the data at a sample temperature of 15−35 K with a photon energy of 104−125 eV and linear-horizontal polarization. At these settings, the energy and k resolutions are better than 20 meV and 0.01 Å −1 , respectively. We carefully aligned the in-plane orientation of black phosphorus along the k y , k s or k x directions with respect to the analyser slit oriented in the direction perpendicular to the scattering plane for symmetry. Then, we repeated the cycle of in situ deposition and data acquisition to take a doping series of ARPES data shown in Fig. 2d-o, Extended Data Fig. 9. After converting the emission angle to k, a series of k distribution curves at a constant energy was fit with the standard Lorentzian function. The spectral weight of C1 integrated in the first Brillouin zone is plotted as a function of energy in Fig. 4d. The energy at which the spectral weight drops by half relative to E F is taken as the magnitude of pseudogap in Fig. 5b. The ARPES data for copper oxides (Extended Data Fig. 10a) are taken at Beamline I05 at the Diamond Light Source with a sample temperature of 5 K and a photon energy of 20-140 eV.

Sample preparation and surface doping
The single-crystal samples of black phosphorus (99.995%, HQ graphene) were glued on sample holders by conductive epoxy. The sample holders were transferred in the ARPES chamber through a semi-automated sample highway system. The black phosphorus samples were cleaved in the ultrahigh vacuum chamber with a base pressure better than 5 × 10 −11 torr. We scanned over the surface with a photon beam of 80 μm in diameter to carefully optimize the sharpness of the ARPES spectra ( Fig. 2d-o). The in situ deposition of alkali metals (Na, K, Rb and Cs) was carried out using commercial dispensers (SAES) at an adsorption rate of over 3.6 × 10 11 cm −2 s −1 . n d is assumed equivalent to n e estimated from the density of surface unit cells (6.9 × 10 14 cm −2 ) multiplied by the ratio of area enclosed by constant-energy contours at E F (Fig. 1d) to that of the surface Brillouin zone. We extrapolate k F in x (k F,x ) and k F in y (k F,y ) from the band dispersion of C1, and estimate the area of Fermi contours as πk F,x k F,y . The unit of ML in Fig. 5b is defined as 4.6 × 10 14 cm −2 , which is the density of closely packed K atoms in two dimensions. Single-crystal samples of Bi 1.5 Pb 0.5 Sr 2 CaCu 2 O 8+δ overdoped at T c = 88 K were used for Extended Data Fig. 10a.

Theoretical model for liquid metals
We employ the theory of multiple scattering 1-12 , in which each ion is assumed to have a spherical step potential whose depth is V 0 and radius is r s (Fig. 3a). This boundary-condition problem is solved based on spherical wavefunctions and partial-wave expansion to obtain the phase shift of partial waves δ l (Fig. 3b). Then, the scat tering amplitude of partial waves f l can be simply written in terms of δ l as: l l δ i l Re(kf l ) and Im(kf l ) take the simple mathematical form of sin(2δ l ) and sin 2 (δ l ) (Fig. 3f), and their relation can be straightforwardly shown by the unitary circle in the complex plane 29 (Fig. 3e). The magnitude of f l in Fig. 3c is given by sin(δ l )/k with its maximum possible value of 1/k r at resonance. By taking the square of spherical wavefunctions at a constant k, the probability density is calculated and plotted as a function of r and k in Fig. 3d after the normalization of the first peaks outside of r s . In the thin-slab approximation (k ≫ Δk, forward scattering) 9,12 , Δk can be written in terms of the predominant f l of the partial wave at resonance as: l d It follows that Re(Δk) and Im(Δk) are directly proportional to sin(2δ l )/k 2 and sin 2 (δ l )/k 2 with respect to k r , respectively. The quadratic band structure of E ≈ k 2 is renormalized by Re(Δk) as k′ = k + Re(Δk) with the spread of k′ over the range of ± Im(Δk), which is the band structure of liquid metals depicted in Fig. 3g.
The resonance scattering, which is characterized by δ l passing through π/2, occurs in one of the partial waves for a given k r . Which partial wave is at resonance scattering depends on the depth of the scattering potential V 0 in Fig. 3a. For a constant V 0 (set to 16.4 eV for d-wave resonance), a partial-wave series of δ l and Δk is calculated as described above and shown in Extended Data Fig. 5. We systematically calculated p-wave, d-wave and f-wave resonances as a function of V 0 with constant r s (set to 2.1 Å for liquid Na 38 ), as shown in Fig. 4a. The representative set of δ l at p-wave, d-wave and f-wave resonances in Fig. 4b is calculated for V 0 = 7.4 eV, V 0 = 16.4 eV and V 0 = 27.8 eV to have the same k r of 0.36 Å −1 (see Extended Data Fig. 6 for the fuller set of U l , δ l and Δk). The density of states in Fig. 4c is calculated using equation (3.13) in ref. 12 with the input parameters of V 0 = 7.15 eV for p-wave resonance and V 0 = 16.26 eV for d-wave resonance at constant r s = 2.1 Å. The half-width at half-maximum of Im(Δk), which corresponds to half the dip width in the energy derivative of Re(Δk) in Fig. 3h, is taken as the magnitude of pseudogap in Fig. 5a.

Spectral simulations
The band structure of black phosphorus doped by disordered metals is simulated based on the model of k renormalizations. The k distribution of ARPES intensity is modelled in the form of a Lorentzian, the integration of which over k is a constant function with energy. The energy distribution of spectral weight, which is proportional to the density of states, is given by n E , which accounts for the diminishing spectral weight towards E F by the formation of pseudogap as: where k E is the dispersion of non-interacting bands in a form of E ≈ k p with k F and the bottom energy taken from ARPES data. p is set to 1.8-2.0 (quadratic) in y and 1.2-1.4 (nearly linear) in x to reflect the well known armchair-zigzag anisotropy of black phosphorus 23 . Re(Δk) and Im(Δk) are taken in a form of sin(2δ l )/k 2 and sin 2 (δ l )/k 2 , respectively, with respect to k r located near k F . The magnitude of Re(Δk) and Im(Δk) is scaled to fit the spectral peak position and width of ARPES spectra. η is the offset k broadening in the range of 0.03−0.14 Å −1 and f FD is the Fermi-Dirac distribution function. n E is the density of states calculated using equation (3.13) in ref. 12 and convoluted by the error function to account for finite experimental broadening. The input parameters are V 0 and r s , the latter of which is set to 2.1 Å. V 0 is optimized to reproduce the magnitude of a pseudogap in the spectral weight obtained by integrating the ARPES intensity of C1 over all of k space in the first Brillouin zone in Fig. 4d, which yields 7.15 eV for Na and 16.26 eV for K at experimental broadening of 0.06−0.14 eV. The simulated spectra are shown in Fig. 4e-j, and directly compared with a fuller set of data in Extended Data Fig. 7.
As for Extended Data Fig. 10b, k E is assumed to be the quadratic band dispersion shown by the grey curve. The n d of our samples is estimated from T c (88 K) 35 , pseudogap (32 meV) 36 and the Luttinger count to be 1.2 × 10 14 −1.4 × 10 14 cm −2 , corresponding to a k r of 0.29 ± 0.01 Å −1 . Re(Δk) and Im(Δk) are calculated based on the model of liquid metals with input parameters of V 0 = 9 eV and r s = 2 Å for the p-wave resonance at k r = 0.29 Å −1 . η is set to be 0.01 Å −1 at E F , from which it monotonically increases with the binding energy up to 0.08 Å −1 . , no structure is expected in its FFT (f). Considering that each dopant is ionized by donating its valence electron, there should be Coulomb repulsion between the ionized dopants, which prevents them from being located close to each other in a certain range corresponding to the radius of hard spheres. If the position of dopants is randomized under this hard-sphere assumption with no further interactions considered (c), one can always find the radial structure factor in its FFT (g). The radial shape itself in the structure factor reflects the presence of a mean interatomic distance (or short-range order) generic to ionic dopants randomly distributed under the repulsive interaction. This is the structure factor required for the resonance effect in the theory of liquid metals [1][2][3][4][5][6][7][8][9][10][11][12] , which is clearly distinguished from both crystalline (a, e) and fully random (b, f) cases.
In the case of monovalent dopants on a crystalline insulator, the area of a circle whose radius is half the radial peak is the same as that enclosed by the Fermi surface in any shape due to charge conservation. To account for the isotropic magnitude of a pseudogap over the anisotropic Fermi surface of C1, the anisotropy of the structure factor is additionally considered based on a recent STM study 25 , where the radial and anisotropic structure factor was found at the density of 1.8 × 10 13 cm −2 . When the anisotropy factor of 0.42 taken from the anisotropic k F values in the Fermi surface is applied to the hard-sphere model (d), the radial structure factor becomes anisotropic, as seen in h.  A series of U l , δ l , Re(Δk) and Im(Δk) is calculated as described in Methods for the partial waves from l = 0 (top row) to l = 3 (bottom row). The depth of potential V 0 is set to 16.4 eV, which corresponds to d-wave resonance at k r = 0.36 Å −1 , as shown in Fig. 4a. The d-wave resonance can be identified by δ 2 passing through π/2, which is accompanied by Re(Δk) and Im(Δk) taking the form of sin(2δ l )/k 2 and sin 2 (δ l )/k 2 , respectively. For l = 0, there is no potential well, in which electron waves can be trapped, and it is impossible to have s-wave resonance.  Fig. 4a). The variation of δ l , Re(Δk) and Im(Δk) is narrower in k or E width for the higher number of l, which accounts for the smaller magnitude of pseudogap in Fig. 5a.

Data
Sim. Cs Sim. -0.4 Sim. The experimental band structure is taken from black phosphorus whose surface is doped by different kinds of alkali metal (Na, K, Rb and Cs), as marked at the bottom right of each panel. This series of ARPES spectra was taken along k y (left column), k s (middle column) and k x (right column). Each dataset is directly compared with corresponding spectral simulations at p-wave resonance for Na and d-wave resonance for K, Rb and Cs (Methods). The bottom energy of spectral simulations is limited by that of the non-interacting (bare) band, which is the reason that those for Rb and Cs are slightly cut at a binding energy of 0.42-0.45 eV.  Fig. 8 | Partial-wave analysis with the screened scattering potential. a, Phase shift δ l of the partial waves calculated by the screened scattering potential for the liquid phase of Cs ions, which clearly shows d-wave resonance 9 . The dotted line is a fit with that calculated by the spherical step potential at V 0 = 12.96 eV and r s = 2.1 Å. b, Phase shift δ l of partial waves calculated by the screened scattering potential of Na, K, Rb and Cs ions 30 . K, Rb and Cs ions favour d-wave resonance, whereas the Na ions favour p-wave resonance, which is exactly as observed in our experiments. Our results reproduce other details as well. (1) δ l crosses π/2 clearly for K, Rb and Cs ions, whereas that for Na ions is relatively incomplete. A doping series of APRES data are taken from black phosphorus whose surface is doped by Na, K, Rb and Cs, marked at the bottom right of each panel. Shown together in number is the dopant density n d in units of 10 14 cm −2 . It can be clearly seen that the magnitude of the pseudogap progressively decreases to zero with increasing n d , as summarized in the phase diagram (Fig. 5b). The hole-like states developing inside of the C1 band (pronounced for Na at n d = 3.3 × 10 14 cm −2 ) are the valence band of black phosphorus 23 . The additional feature near E F and at the Γ point (pronounced for K at n d = 2.9 × 10 14 cm −2 ) is another conduction band of black phosphorus reproduced in first-principles band calculations 24 .