Electronic Structure of Few-Layer Black Phosphorus from μ-ARPES

Black phosphorus (BP) stands out among two-dimensional (2D) semiconductors because of its high mobility and thickness dependent direct band gap. However, the quasiparticle band structure of ultrathin BP has remained inaccessible to experiment thus far. Here we use a recently developed laser-based microfocus angle resolved photoemission (μ-ARPES) system to establish the electronic structure of 2–9 layer BP from experiment. Our measurements unveil ladders of anisotropic, quantized subbands at energies that deviate from the scaling observed in conventional semiconductor quantum wells. We quantify the anisotropy of the effective masses and determine universal tight-binding parameters, which provide an accurate description of the electronic structure for all thicknesses.

stack was released on Au-coated Si/SiO 2 wafers held at 130 • C during the initial contact before heating to 180 • C (above the glass transition of PC around 150 • C) for the transfer from the stamp to the substrate. Immediately after the transfer, the heterostructures were rinsed in CH 3 Cl (chloroform) for at least 2 hours. In order to reduce the exposure of BP to reactive gases the entire heterostructure fabrication was performed in a N 2 filled glovebox within ∼ 1 h of the exfoliation of BP. Oxygen and water contamination in the glove box were typically below 0.5 ppm.
The as-prepared heterostructures were transported under glovebox atmosphere to the ARPES system using a custom built suitcase, which was pumped down to ultra-high vacuum (UHV) pressure without the samples being exposed to ambient atmosphere at any time. Prior to the measurements, the samples were annealed at 200 • C for a few minutes, which may remove adsorbants and agglomerate the contaminants trapped in between the 2D materials.
The ARPES measurements on bulk BP were performed on single crystals grown with a modified chemical vapor transport method using a mineralizer as reaction promoter [1]. Specifically, a mixture of 300 mg red phosphorus (99.999%), 20 mg Tin powder (99.995%), and 10 mg Tin(IV)iodide powder (99.998%) were weighed and sealed in an evacuated quartz tube with a length of about 10 cm. The tube was placed horizontally in the center of a box furnace.
The temperature of the furnace was increased from room temperature to 650 • C at a rate of 1 • C/min and held at this temperature for 4 hours. Subsequently, the temperature was decreased to 520 • C within 30 hours and was kept at this temperature for 4 hours. Afterwards, the sample was cooled to room temperature over 10 hours. The black phosphorus single crystals were collected from the quartz tube in a glove box. Residual iodine was removed with absolute ethanol, and the crystals were dried in N 2 flow.

II. ARPES MEASUREMENTS
The ARPES measurements of few layer BP were performed in a custom built micro-ARPES system [2]. This setup is based on a continuous wave laser delivering up to 10 15 photons/s in a µeV bandwidth at a wavelength of 206 nm (6.01 eV). The laser beam is first expanded and then focused with an aberration corrected lens mounted in UHV on a 3-axes translator. The samples are mounted on a conventional 6-axes ARPES manipulator and scanned with a stepper-motor driven, bellow sealed xyz-stage. An overall spatial resolution < 2 µm was determined from knife edge scans. ARPES experiments are performed with a hemispherical analyzer (MB Scientific) equipped with a deflector lens permitting the acquisition of 2D k-space maps without rotating the sample. The typical energy and momentum resolution of the measured photoelectrons was 0.003Å −1 and 5 meV. The experiments used linear polarization and a photon flux of approximately 10 14 s −1 (≈ 100 µW), which is around 2 orders of magnitude higher than typical photon fluxes used at synchrotron based micro-ARPES beamlines based on capillary optic focusing. All micro-ARPES data were acquired at a pressure p < 10 −10 mbar and a manipulator temperature of T ∼ 4.5 K.
The ARPES measurements of bulk BP were performed at the SIS beamline of the Swiss Light Source using circularly polarized light with energies of 20 − 40 eV. Bulk crystals were cleaved in UHV at a temperature of 10 K, but measured at 50 K in order to avoid charging.

III. TIGHT-BINDING MODEL OF FEW LAYER BLACK PHOSPHORUS
We follow Ref [3] and describe the band structure of few-layer black phosphorus with an effective single orbital tight-binding model with the Hamiltonian: Here, t with the matrix elements: The vectors r ∥ i represented in Fig. S1, are given by: In multilayer BP, sublattices A and B of one layer are coupled to sublattices C and D of the adjacent layer. We describe this coupling with the submatrix: with matrix elements: The vectors r ⊥ i are given by: The tight-binding Hamiltonian for multilayer BP can now be written as: For our analysis of the experimental data, we used the lattice parameters a 1 = 2.22Å, a 2 = 2.24Å, α = 96.5 • , β = 72 • and diagonalized the Hamiltonian in Eq. 8 numerically to obtain the tight-binding band structure of few-layer BP.
The tight-binding hopping elements were determined by fitting the eigenvalues of Eq. 8 to the experimental band dispersion of 2L, 3L, 4L and 5L BP along both high symmetry directions. The fits were additionally constrained to reproduce a bulk band gap of 400 meV, as it was used in Ref [3] -close to the experimentally determined bulk bandgap (about 300 meV) [4]. With this procedure, we find the hopping parameters shown in TABLE 1 of the main text. Figure S2. a) Comparison of the electronic gap from our tight-binding model with ab-initio electronic structure calculations from Refs. [3,[5][6][7][8][9][10][11][12]. b) Comparison with optical gaps reported by Refs. [11,[13][14][15]. The TB bands were obtained by constraining the gap of bulk BP to 400 meV.
As shown in Fig. S2, these hopping elements also provide a fair description of measured optical gaps and of quasiparticle gaps from ab-inito calculations. Note that the quasiparticle-and optical gap differ by the exciton binding energy. We speculate that the two gaps might be similar in our devices because the renormalization of the quasiparticle band gap due to dielectric screening in the encapsulation layer [16] might be comparable to the exciton binding energy.

IV. THICKNESS DETERMINATION OF EXFOLIATED BLACK PHOSPHORUS FLAKES
We determine the thickness of BP flakes from a combination of atomic force microscopy (AFM) measurements, the optical contrast on Si/SiO 2 substrates, tight-binding model, ab-initio calculations, and our ARPES data, as illustrated in Fig. S3. AFM measurements were performed under ambient conditions in a Cypher AFM from Oxford Instruments using the tapping mode. For 2L − 4L flakes, we find a nearly linear optical contrast in the red channel of a calibrated microscope. This allows for a rapid determination of the thickness of very thin flakes prior to the encapsulation in graphene/graphite. For thicker flakes, the contrast begins to saturate and we find it more reliable to determine the thickness by comparing measured subband energies with our tight-binding model fitted to 2L − 5L.
No 1L data is shown throughout the manuscript because we have not succeeded in obtaining high-quality monolayer heterostructures. Density functional theory (DFT) calculations have been carried out using the Quantum ESPRESSO distribution [17].
A pseudopotential approach has been adopted to treat electron-ion interactions, with the P ultrasoft pseudopotential chosen from the pslibrary 1.0.0 [18], which yield results in agreement with all electron calculations. The cutoff energy on wavefunctions was set at 40 Ry, and at 400 Ry for the density. To account for van der Waals interactions between the layers the optB88 exchange-correlation functional has been adopted [19], while spurious interactions between artificially periodic replicas of the multilayers along the vertical direction were suppressed using a cutoff on all relevant interactions [20]. The Brillouin zone was sampled with a Γ-centered 10 × 12 × 1 Monkhorst-Pack grid. The atomic structure and unit cell parameters have been fully relaxed until forces on atoms were smaller than 12.5 meV/Å and stresses below 0.5 kbar. We found that the lattice parameters are strongly dependent on the multilayer thickness and rapidly approach the bulk values.
In Fig. S4, we compare our DFT band structure with the subband dispersion extracted from experiment. Note that our DFT calculations for 9L BP find a metallic state with two Dirac points along the zigzag high-symmetry line. Figure S4. Comparison of the experimental dispersion from the main text with our ab-initio calculations. All energies are relative to the valence band edge E VBM . For clarity, the experimental data points have been symmetrized around k ∥ = 0. The colours symbolize the band index.
In Fig. S5, we compile subband energies and effective masses from published theoretical studies and compare them with the values determined from our experiment and ab-initio calculations. A significant spread of theoretical subband energies is evident from the figure. We note that the agreement between our ab-initio calculations and experiment is good for subbands energies and zigzag effective masses whereas our ab-initio effective masses for the armchair direction are over a factor of two smaller than in experiment. Figure S5. a) Subband energies at the Γ point. Note that all values are given relative to the valence band maximum E VBM . Orange markers represent the experimentally determined subband energies. Theoretical values from Refs. [3,[5][6][7][8][9][10][11] are shown in blue, violet and burgundy. The subband energies determined from the ab initio calculation of our study are shown as hollow red rectangles. Shaded areas group the different subbands. b-c) Effective masses determined from our experiment and from the electronic structure calculations of Refs. [6,9,12].

VI. SPECTRAL FUNCTION OF A FILLED BAND COUPLED TO EINSTEIN PHONONS
At second-order in ionic displacements and neglecting the Debye-Waller correction, electron-phonon coupling is described by the Hamiltonian [21]: H = k,n,σ ξ k,n c † n,k,σ c n,k,σ + q,ν ℏω q (b † q,ν b q,ν + 1/2) + k,q,m,n,ν,σ g mnν (k, q)c † m,k+q,σ c n,k,σ (b † q,ν + b q,ν ) Here g mnν (k, q) are the electron-phonon coupling matrix elements, ξ k,n is the dispersion of band n relative to the chemical potential µ with corresponding creation and annihilation operators c † m,k,σ , c † m,k,σ for each spin orientation σ. ω q,ν is the phonon frequency in branch ν for crystal momentum q and b † q,ν , b q,ν are the corresponding phonon creation and annihilation operators. To capture the essence of the spectral function of an electron-phonon coupled system, it is sufficient to consider a single band and a single, dispersionless phonon branch, Ω 0 . We further assume that this phonon mode couples to the electrons with a momentum independent electron-phonon matrix element, which we denote by g 2 . Using these simplifications, the self-energy Σ is given at second order by: where ρ(ξ) is the density of states per spin of the single-particle band considered here, f (ξ) and b(ℏΩ 0 ) are the Fermi-Dirac and the Bose-Einstein distribution, respectively and γ is a broadening factor that can be interpreted as impurity scattering. Note that within this approximation the self-energy is local, i.e. it does not dependent on momentum.
From Eq. 10, one obtains for the mass enhancement λ at zero temperature: Here, ρ 0 is the bare density of states per spin at the subband edge. We point out that this expression differs by a factor of 2 from the equivalent expression for metals [21]. This is because the second term in Eq.10 vanishes at zero temperature in our case, contrarily to the case of metals. [ω−ξ k −ReΣ(ϵ)] 2 +[ImΣ(ϵ)] 2 obtained from this model for a realistic phonon frequency ℏΩ 0 = 60 meV and a weak to moderate coupling strength λ = 0.32. The calculation reproduces a non-dispersive spectral feature at the same energy observed in the experiment. However, the weight of this feature is smaller than in our data and decays more rapidly away from the main band. The same behavior is found within the momentum average approximation of electron-phonon coupling [22]. Stronger effects of electron-phonon coupling on the spectral function are obtained using the momentum average approximation in 1D [22,23] but this can be traced back to the divergence of the bare density of states at the band edge, which is absent in 2D.