Probing Single Vacancies in Black Phosphorus at the Atomic Level.

Utilizing a combination of low-temperature scanning tunneling microscopy/spectroscopy (STM/STS) and electronic structure calculations, we characterize the structural and electronic properties of single atomic vacancies within several monolayers of the surface of black phosphorus. We illustrate, with experimental analysis and tight-binding calculations, that we can depth profile these vacancies and assign them to specific sublattices within the unit cell. Measurements reveal that the single vacancies exhibit strongly anisotropic and highly delocalized charge density, laterally extended up to 20 atomic unit cells. The vacancies are then studied with STS, which reveals in-gap resonance states near the valence band edge and a strong p-doping of the bulk black phosphorus crystal. Finally, quasiparticle interference generated near these vacancies enables the direct visualization of the anisotropic band structure of black phosphorus.

Black phosphorus (BP) is a layered allotrope of phosphorus, which crystallizes into an orthorhombic lattice.[1] Weakly bound phosphorene layers are vertically stacked in an arrangement analogous to graphite; in contrast, however, each phosphorene layer possesses two planes of phosphorus atoms bonded though sp3 orbitals. This leaves two electrons per phosphorus atom paired, yet unbonded.[2,3] The 1D buckling and lone-pair electronic configuration result in an electronic band gap for bulk black phosphorus of 0.3 eV,[1] which is predicted to be dramatically enhanced in the single layer limit (∼1.5–2 eV) due to the suppression of interlayer interactions.[4,5] This strong variation in the electronic structure allows for optical identification of variations in layer thickness,[6,7] which, in addition to the high carrier mobility, make black phosphorus extremely promising for device applications based on 2D layered materials.[8−13] Furthermore, black phosphorus contains an additional distinction: a strongly anisotropic band structure along the principle direction Γ–X compared to Γ–Y, which can be characterized by two free electron-like bands with starkly different effective masses.[14,15] To date, measurable anisotropies have been reported for in-plane carrier mobility, optical extinction, thermal transport, and plasmon dispersion.[14,16−19]

In addition to device applications, it has been proposed that the band structure of black phosphorus can be manipulated by strain or electric fields toward a superconducting phase transition[20,21] or a Lifshitz topological transition.[15,22−24] The first experimental work to this end used surface alkali metals to controllably close the band gap and demonstrate a semimetal transition in which one of the two free electron bands transformed into a Dirac-like band.[25,26] This opens up the possibility of controlling both the topology and band gap of black phosphorus based on surface electric fields.

However, unlike graphene, MoS2, or WS2, black phosphorus is extremely reactive, making its electronic properties sensitive to ambient conditions.[27,28] This motivates a careful understanding of the intrinsic doping of black phosphorus and how native defects modify the band structure.[29−32] Scanning tunneling microscopy/spectroscopy (STM/STS) has been used to characterize the intrinsic band gap of bulk black phosphorus in ultrahigh vacuum, but conflicting values have been reported.[33,34] Here, we characterize the structural and electronic properties of black phosphorus and the influence of inherent single vacancies on those properties, utilizing a combination of low temperature STM/STS with electronic structure calculations. We illustrate that black phosphorus is p-doped with a band gap of ∼0.3 eV; this doping is facilitated by single vacancies, which are characterized by an anisotropic charge density and in-gap resonance states near the valence band edge. Combining tight-binding calculations with atomic-scale STM/STS, we identify the structure of these vacancies as well as depth profile them, correlating their electronic properties near the band gap with their depth from the surface. Moreover, we demonstrate via quasiparticle interference that these vacancies act as strong scattering centers illuminating the anisotropic band structure.

Black phosphorus cleaves along the [001] direction due to weak van der Waals interactions; the STM image in Figure a shows the resultant atomic structure after cleaving in situ, which leaves a macroscopically flat black phosphorus (001) surface. Constant-current STM images reveal zigzag rows (along the [010] or zigzag direction), which correspond to the upper plane of phosphorus atoms (light orange atoms in Figure b) in the topmost layer of the buckled rhombohedral structure. As each phosphorus atom contains a lone pair of nonbonding electrons, the charge density per atom extends quite far into the vacuum, leading to the significant (0.02 nm) atomic corrugation observed in the STM topography. The crystal structure yields two mirror-symmetric sublattices, which we label A and B in Figure b, within each [010]-oriented atomic row. In general, both sublattice sites appear equivalent in STM (Figure c), except in the vicinity of a defect. A close examination of the atomic structure (Figure c) enables the determination of the lattice parameters within the (001) plane: a = 0.32 nm and c = 0.44 nm. These parameters are consistent with previous experimental and theoretical studies of unstrained black phosphorus.[1,35] A FFT of the atomic resolved image in Figure a (Figure d) allows for the identification of the high-symmetry points of the topmost layer, where we adopt the convention denoting the real-space [100] direction as armchair and [010] as zigzag.

Figure 1

(a) Constant-current STM image of black phosphorus after in situ cleaving (VS = −0.1 V, It = 200 pA, scale bar = 1 nm). (b) Crystal structure of black phosphorus with side view (upper) illustrating the buckling of each layer. The top view (lower) illustrates the two sublattices (A,B), in comparison to which atoms are imaged with STM (bright orange). (c) Zoomed-in image of black phosphorus, with lattice parameters a = 0.324 nm and c = 0.44 nm (color bar: 0 to 0.02 nm). (d) Fourier transform of the STM image in (a) with green dashed lines showing Brillouin zone boundaries and the principle in-plane reciprocal space locations denoted (scale bar = π/c = 7.14 nm–1). (e) STS of black phosphorus (Vstab = 0.35 V, Istab = 100 pA, Vmod = 5 mV, fmod = 4.2 kHz) illustrating an electronic band gap of 0.317 ± 0.036 eV. The logarithm of the dI/dV tunneling spectra is shown in gray to accentuate the noise floor within the gapped region. The band gap was derived from n = 48 independent spectra.

To understand the electronic properties of black phosphorus, STS was performed to characterize the band gap. Using the analysis shown in Figure S1 for n = 48 individual measurements, we estimate Eg = 0.317 ± 0.036 eV. This band gap is consistent with previous large-area measurements on bulk black phosphorus crystals (0.31–0.35 eV).[1] Scanning tunneling spectroscopy further illustrates that the Fermi level (EF) generally resides within ±30 mV of the valence band edge (Figure S2); considering the band gap, this indicates that black phosphorus is strongly p-doped without deliberate modification. To elucidate the origin of this intrinsic p-doping, we investigate below the structural and electronic properties of the atomic-scale defects in black phosphorus.

Figure a shows a representative large-area constant-current image of black phosphorus near the valence band edge. The surface is decorated with a low density of characteristic, elongated, dumbbell-shaped protrusions, which exhibit enhanced charge density along the [100] (armchair) direction. As these features are present after cleaving in ultrahigh vacuum, we associate them with native defects in the black phosphorus crystal. Line profiles across these defects, shown in the upper left inset of Figure a, reveal that measured heights of the defects range from 0.01 to 0.2 nm above the pristine surface, spanning nearly an order of magnitude greater than the intrinsic corrugation of the atomic lattice. The monotonic changes in apparent height (Figure ) are consistent with STM studies of defects at various heights beneath the surface of III–V semiconductors,[36−39] as well as in various layered compounds.[40,41] We discuss the depth profiling of the defects in more detail later.

Figure 2

(a) Three-dimensional representation of a constant-current STM image with a distribution of single vacancies in black phosphorus (VS = −0.1 V, It = 200 pA, size = 62 nm × 48 nm, color bar = 0–0.2 nm). (Inset) Line profiles taken across the vacancies labeled 1L, 2L, and 3L in (a). Constant-current STM image of a single vacancy at (b) sublattice A and (c) at sublattice B (VS = −0.1 V, It = 200 pA, scale bar = 2 nm). (d) Tight-binding calculations of the charge density of a single vacancy in black phosphorus located at (d) sublattice site A and (e) sublattice site B (scale bar = 1 nm).

Figure 3

(a) Scatter plot of the averaged apparent STM height for vacancy A (see inset). Blue points refer to the mean intensity from both sides of the vacancy (black and gray contours in inset), while the gray squares correspond to the left side with estimated error, and the black diamond shows the height of the right side with corresponding error. (b) Scatter plot of the average apparent height for vacancy B. Red points refer to the mean intensity from both sides of the vacancy, while the gray square shows the left side with estimated error, and the black diamond shows the height of the right side with corresponding error. For all data points: VS = −0.04, −0.05, or −0.1 V, It = 100–200 pA.

Atomic-scale characterization of the observed dumbbell structures is shown in Figure b,c, for two different depths below the surface. Close examination of these images reveals that the two dumbbell structures can be transposed onto one another through a reflection along [010] (see also Figure S3–S4). Utilizing tight-binding calculations, we investigate the resultant charge density associated with midgap states near single vacancies (Figure d,e) in the uppermost layer at both sublattice positions (A, B). The calculations qualitatively reproduce the highly anisotropic charge density seen in STM images (Figure S5), as well as a sublattice-dependent orientation of the charge density. By comparing the curvature near the defect center, each experimentally observed dumbbell structure can be directly assigned to a vacancy at a particular sublattice. Therefore, we refer to the vacancies at sublattice A and B, as vacancy A and B, respectively. Furthermore, tight-binding calculations show that vacancies below the first layer maintain the dumbbell defect shape, but have diminished intensity with respect to vacancies in the surface layer (Figure S6). In reproducing key experimental observations, tight-binding calculations reveal that it is possible to trace each observed defect to a particular depth below the surface as well as the atomic sublattice within that layer.

In order to thoroughly understand the relationship between the experimentally determined apparent height and the vacancy depth below the surface, over 100 vacancies were characterized according to their sublattice location and constant-current height profile. This analysis, with the index sorted according to apparent height, is presented in Figure for both vacancy types. From the total number of characterized vacancies, we observe vacancy B nearly twice as frequently as vacancy A. The plots in Figure display the mean intensity extracted from line profiles of the apparent height (illustrated in the insets) for (1) both sides of the vacancy (red/blue), (2) the left side (gray square), (3) and the right side (black diamond). Clear trends emerge in the apparent height for vacancy A, seen in the step-like features at approximately 0.01, 0.04, 0.07, and 0.10 nm. The plateaus are attributed to vacancies residing in discrete layers below the surface, which we label by increasing number with increasing depth underneath the surface, with the label 0 referring to vacancies in the surface layer. Furthermore, a majority of A vacancies exhibited an enhanced intensity on the left region of the charge density compared to the right, with asymmetry between regions decreasing with increasing depth below the surface. Height plateaus for vacancy B are not readily identifiable in Figure b; however, there is a height dependent trend in the dI/dV spectra, as we show below (Figure a and Figure S6), similar to vacancy A, which we believe is related to the specific layer assignment. The preference of vacancy B compared to vacancy A, combined with the differences in their height profiles indicates that the symmetry of the A–B sublattices is broken near the surface at the vacancy site. We attribute this trend to surface relaxation effects, which may relax the two sublattices asymmetrically at the vacancy site near the surface. While the symmetry breaking of the A–B sublattices may be due to vertical relaxation, we note that we cannot rule out a small vertical relaxation for deeper vacancies, which show A–B sublattice symmetry.

Figure 4

Tunneling spectroscopy of the depth dependence and spatial variation for both vacancy types taken at red circles. (a,b) Representative tunneling spectra sorted by decreasing apparent height (h) according to progressive blue scale in the middle for (a) vacancy A and (b) vacancy B taken at the position of the red dot as denoted in the inset. All apparent heights are within ±0.01 nm, except for the h = 0.165 nm, which is ±0.02 nm. (c,d) Spatially resolved tunneling spectra taken for vacancy B (h = 0.150 nm) along: (c) [100] (armchair) direction over a total length of 20 nm and (d) along [010] (zigzag) direction over a total length of 10 nm. Stabilization conditions for all spectra: Vs = 0.35 V, It = 200–400 pA. The red spectra were taken at the red point in the inset.

To correlate the electronic properties of each vacancy type with their depth from the surface, spatially resolved tunneling spectroscopy was performed for vacancies at various depths from the surface (Figure ). All spectra for both vacancy types taken near the center of the vacancy share two key departures from the pristine band gap spectrum (gray curve at the bottom of Figure a,b): (a) spectral resonance(s) emerge in all curves near EF and (b) nonzero tunneling conductance is measured throughout most of the bulk band gap in the proximity of the vacancy (see Figure S7). Moreover, the spectra (Figure a,b) clearly show a phenomenological trend with respect to the resonance states near EF. The upper curves (h > 0.06 nm) show two distinct spectral features split (not necessarily equivalently) about the Fermi level. For both vacancy types, these split peaks are seen for vacancies near the surface, whereas deeper vacancies converge to a single spectral resonance slightly above EF (see also Figure S8). Furthermore, while the absolute energy of the double-resonance structure for vacancy B varies, the separation between the resonance features is nearly constant with a value of approximately 27 ± 3 mV. The rigid separation likely originates from the monovalent nature of single vacancies, while the absolute position of spectral features can be influenced by the depth from the surface and the surrounding environment. Additionally, the relative intensity of the spectral peaks varies with the precise measurement position (see Figure S9). It is important to note that we do not observe any significant effects of tip-induced band bending (Figure S10), most likely due to the pinning of the valence band. At h < 0.05 nm, the spectra for vacancies A and B merge to a single resonance at 9 ± 2 mV. This further indicates that the symmetry between both sublattices may be broken by surface relaxation for vacancies at or near the surface (e.g., 0L–2L), and this symmetry is restored in the bulk. These spectral features reveal that the atomic vacancies host shallow donor states, which contribute to the inherent p-doping and lead to electrical conduction within the band gap.

Spatially resolved STS can also be utilized to map the lateral variation of the observed spectral features. The line spectra in Figure c show that the spectral resonances near the Fermi level are localized to the topographically bright dumbbell regions (red point in inset and red curve, see also Figure S7). Such resonance states are observable as far as 5 nm from the vacancy site along [100]. In the direction of [010] (Figure d), however, the resonance states decay much faster, generally within 2 nm of the vacancy site. The strong wave function anisotropy likely stems from the drastic differences in the effective masses of the valence band along Γ–X and Γ–Y. The full spatial profile of the differential tunneling conductance is shown in Figure S7, at energies both off and on resonance.

The anisotropic nature of the black phosphorus band structure is strongly pronounced in quasiparticle interference (QPI) near the vacancy sites.[42,43]Figure a shows a constant-current dI/dV image acquired at VS = −0.6 V near several vacancies with varying depth from the surface. Clear oscillations are observed around each vacancy, with a much stronger intensity along [100] compared to [010]. This asymmetric scattering pattern is a signature of the anisotropy of the valence bands near the Fermi edge, giving rise to nesting of the quasiparticle scattering vectors. This asymmetry can be clearly traced in the Fourier transform of the dI/dV map (Figure b), where higher intensity is seen for vectors oriented along the Γ–X direction. The edge of the Brillouin zone (BZ) is also highlighted with a green rectangle for reference. As seen in the calculated constant energy contours in Figure d (progressive ellipses around Γ, see Figure S9 for cuts on the band structure), the shape of the valence band is strongly ellipsoidal near EF. This ellipticity causes an imbalance of available scattering vectors in the Γ–X versus Γ–Y directions leading to more prominent carrier scattering along [100]. As the calculations reveal, the constant energy contours become more circular at increasing negative energies, which weakens the nesting condition and results in stronger intensity fringes along [010], as shown in Figure c.

Figure 5

(a) STS map (taken concurrently with constant-current topography) at Vs = −0.6 V showing interference patterns surrounding point defects in BP (Vs = −0.6 V, It = 600 pA, scale bar = 8 nm). (b) FFT of STS image in (a) (scale bar = π/a = 7.14 nm–1). The open-ended ellipse around Γ corresponds to valence band carrier scattering at defect sites. (c) Dispersion of scattering behavior at indicated energies (Vs = varying, It = 600 pA, scale bars = 4 nm). (d) Theoretical calculation of valence band constant energy contours in the vicinity of Γ, shown with the BZ for comparison. From dark to light red (outer to inner): E = −0.9, −0.7, −0.5, and −0.3 eV. For (a) and (c), Vmod = 5 mV and fmod = 4.2 kHz.

In conclusion, we have characterized the structural and electronic properties of inherent single vacancies in black phosphorus. We illustrate that vacancies exhibit a highly anisotropic and delocalized charge density. Moreover, intrinsic vacancies possess in-gap resonance states near the valence band edge, which strongly perturb the bulk band gap. Therefore, we conclude that intrinsic vacancies contribute to the observed strong p-doping of black phosphorus. This work motivates further studies of the bulk doping to conclude if intrinsic vacancies are solely responsible for the observed strong p-doping. We further reveal that these single vacancies serve as strong scattering centers, as seen in quasiparticle interference near the observed vacancy sites. This study provides fundamental insight into the influence of intrinsic point defects on the electronic properties of black phosphorus, which is crucial for future developments of black phosphorus technologies. Moreover, the recent findings of magnetic vacancies in graphene and indications of magnetism in black phosphorus vacancies,[32] make this an interesting platform to study magnetic ordering in sp-driven systems.[44]

Experimental Methods

STM/STS measurements were carried out on a commercial Omicron low-temperature STM with a base temperature of 4.6 K, operating in ultrahigh vacuum (<1 × 10–10 mbar), with the bias applied to the sample. Etched W tips were utilized for these measurements and were treated in situ by electron bombardment, field emission, as well as dipped and characterized on clean Au surface. Scanning tunneling spectroscopy was performed using a lock-in technique to directly measure dI/dV, with a modulation frequency of fmod = 4.2 kHz and amplitude of Vmod = 5 mV. Black phosphorus crystals were provided by HQ graphene and stored in vacuum at a temperature less than 25 °C, cleaved under ultrahigh vacuum conditions at pressures below 1 × 10–9 mbar, and immediately transferred to the microscope for in situ characterization.

Theoretical Calculations

Theoretical calculations were based on the tight-binding model proposed for multilayer BP in ref (5). To model structural defects, we considered a missing atom in a three-layer supercell with dimensions (16a × 12c) ≈ (52.5 × 53.0) Å and performed exact diagonalization of the tight-binding Hamiltonian. Wave functions of a defect state were represented as Ψ(r) = Σici · φ(r), where ci is the contribution of the ith atom given by the corresponding eigenvector of the Hamiltonian, and φ(r) is the cubic harmonic representing the 3p-like orbital of P atoms, being the basis functions of the Hamiltonian. Constant energy contours appearing in Figure d were extracted directly from first-principles GW calculations performed in ref (5).