Surface buckling of phosphorene materials: determination, origin and influence on electronic structure

TL;DR

This study uses advanced microscopy and calculations to reveal significant surface buckling in phosphorene materials, linking it to surface vacancies and intrinsic hole doping, which impacts their electronic properties.

Contribution

It provides the first detailed measurement of surface buckling in phosphorene and links it to surface vacancies and doping effects using combined experimental and theoretical methods.

Findings

Surface buckling of 0.22 Å in black phosphorus

Surface buckling of 0.30 Å in exfoliated phosphorene

Surface vacancies contribute to buckling and doping

Abstract

The surface structure of phosphorene crystals materials is determined using surface sensitive dynamical micro-spot low energy electron diffraction ({\mu}LEED) analysis using a high spatial resolution low energy electron microscopy (LEEM) system. Samples of (\textit{i}) crystalline cleaved black phosphorus (BP) at 300 K and (\textit{ii}) exfoliated few-layer phosphorene (FLP) of about 10 nm thicknes, which were annealed at 573 K in vacuum were studied. In both samples, a significant surface buckling of 0.22 {\AA} and 0.30 {\AA}, respectively, is measured, which is one order of magnitude larger than previously reported. Using first principle calculations, the presence of surface vacancies is attributed not only to the surface buckling in BP and FLP, but also the previously reported intrinsic hole doping of phosphorene materials.

Tables

Table 1. Table 1: Optimum parameter values for the surface structure of BP crystal and exfoliated BP flake

Model	T	$b_1$ (Å)	$b_2$ (Å)	$z^{\prime }$ (Å) ($\Delta z/z$)	$w^{\prime }$ (Å) ($\Delta w/w$)
Cleaved BP	300 K	0.225	$0.269$	2.287 (+5.3%)	2.825 (-8.0%)
FLP	573 K	0.300	$0.290$	2.381 (+9.9%)	2.877 (-6.3%)
DFT 10	-	$0.02$	-	-	-
DFT 11	-	$0.06$	-	-	-

Equations

$$<u^{2}>_{T}{}=\frac{9\hbar^{2}}{m_{a}k_{B}\theta_{D}}\;\;\Big{\{}\frac{~{}T^{2}}{\theta_{D}^{2}}\int_{0}^{\frac{\theta_{D}}{T}}\frac{xdx}{e^{x}-1}+\frac{1}{4}\Big{\}}$$

Full text

\altaffiliation

Current address: Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA

\altaffiliationCurrent address: Department of Physics and Center for Humanities and Sciences, Kyung Hee University, Seoul 02447, Republic of Korea

\abbreviations(BP), black phosphorus; (FLP) few-layer phosphorene; (LEEM), low energy electron microscopy; ($\mu$LEED) selected area micro-spot low energy electron diffraction; (PEEM) photo-emission electron microscopy

Black phosphorus (BP), which in its monolayer version is denoted as phosphorene, has recently had a rebirth as a new member of the family of the vigorously studied two-dimensional (2D) materials. It has attracted much attention due to its intriguing potential applications for modern electronics 1, 2, 3, 4 and photonics 5, 6. For example, BP exhibits an intrinsic layer-dependent bandgap ranging from 0.3 eV (bulk) to 2 eV (monolayer) 7, which bridges the energy gap between graphene and transition metal dichalcogenides (TMDs) 8. This strong layer-dependent band gap presents a potential for an interesting set of integrated devices on a single supporting platform. Despite this surge of research in the applications of phospherene materials, much remains to be learned of its basic physical properties both from a device and a fundamental physics perspective. For example, its basic surface lattice structure is not well known, but is crucial both for predicting electronic properties and exploring phosphorene based devices.

In fact to date, there is no consensus on the atomic structure of surface region of BP. The crystal structure of BP, as shown in Fig. 1, has a puckered honeycomb structure similar to that of graphene 9. Two previous STM studies of phosphorene 10, 11 have revealed some aspects of the BP surface topography and observed an apparent height difference between the two symmetrically equivalent atoms P1 and P2, as shown in Fig. 1 (d). While these STM measurements were not able to quantify the geometrical height difference between P1 and P2, denoted as surface buckling, these studies did propose very small surface buckling values, 0.02 Å 10 and 0.06 Å 11, based on first-principles calculations. In order to experimentally resolve the surface atomic structure of BP, there are two main challenges for the characterization technique: it has to be (i) non-destructive and sensitive to the 3D atomic structure in the first few layers, and (ii) able to restrict the sampling area to a few $\mu$m because many 2D materials including phosphorene are commonly prepared as very limited area exfoliated flakes. Selected area micro-spot low energy electron diffraction ($\mu$LEED) in a low energy electron microscope (LEEM), combined with dynamical intensity versus incoming electron energy (LEED-IV) calculations, is one of the very few practical techniques to determine the 3D surface structure and composition of 2D materials with atomic resolution 12, 13, 14, 15, 16.

In this letter, we present the first detailed experimental atomic surface structure determination of BP. We use LEEM and dynamical $\mu$LEED-IV analysis to examine in-situ-cleaved bulk BP surface and mechanically-exfoliated few-layer phosphorene (FLP) flakes of about 10 nm thickness. We measure the surface buckling for the two studied systems to be 0.22 Å and 0.30 Å, respectively, which are one order of magnitude larger than two previously reported theoretical values. Finally we use first-principles calculations to identify that the presence of surface vacancies is very likely the origin of not only the surface buckling in BP, but also the intrinsic hole doping of phosphorene that was reported previously 3, 17.

Our experiments were carried out on the Elmitec AC-LEEM system at the Center for Functional Nanomaterials in Brookhaven National Laboratory. The spatial resolution in LEEM mode is better than 3 nm and the electron beam spot size is 2 $\mu$m in diameter in the $\mu$LEED mode. Single-crystal bulk BP was cleaved in ultrahigh vacuum at room temperature. Figure 2 (a) shows the real-space bright field LEEM image of a freshly cleaved BP surface. $\mu$LEED data were acquired at the region denoted by the red 2 $\mu$m circle using a normal incident electron beam. Figure 2 (b) shows the well defined LEED pattern at 35 eV electron energy, indicating a very well-ordered surface. To prepare for the FLP samples, black phosphorous flakes were mechanically exfoliated onto n-doped Si chip with native oxide, using a previously described method 18, 19. The substrate was pre-patterned with gold marks, which allowed us to locate and characterize the flakes of interest using optical microscope; see Fig. 2 (c). This procedure was performed in a Ne atmosphere. Subsequently, the sample was encapsulated and transferred to the LEEM chamber. The total exposure time of the exfoliated sample to air was less than 5 minutes. Even with such a short exposure time, significant surface oxidization and contamination was observed using photo-emission electron microscopy (PEEM). In order to remove the contamination, we annealed the sample at 300*∘C in ultrahigh vacuum for 2 hours. The correct annealing temperature is crucial for successful contamination removal; higher temperatures will lead to accelerated sublimation. As shown in the PEEM and LEEM images, Fig. 2 (d)-(e), the surface was pristine and uniform after successful annealing. $\mu$LEED data were acquired at the region denoted by the red circle using normal incident electrons. Figure 2 (f) shows the sharp LEED pattern at 25 eV electron energy, indicating a very well-ordered layered structure. To fully investigate the surface atomic structure we collected $\mu$LEED-IV spectra for 7 recorded diffraction spots with an energy range of 25 to 135 eV for both sample varieties. The intensities of symmetrically equivalent beams were averaged to minimize intensity anisotropy of the diffraction beam due to small sample titling ($<$ 0.1∘*). Specifically, as shown in Fig. 2 (b) and (f), intensities of spots A were averaged to assign the (01) diffraction beam and beam intensities of spots B were averaged to assign the (11) diffraction beam. The background intensity was then subtracted from the diffraction beam intensity.

Dynamical LEED-IV analysis were carried out to extract the surface atomic structural information for bulk BP and FLP from the corresponding $\mu$LEED-IV curves. In a dynamical LEED-IV analysis, IV curves are calculated for a trial structure and compared with experimental curves. A $\chi^{2}$-based $R_{2}$ reliability factor is used to quantify the difference between calculated and experimental IV curves 20. The surface structural parameters are then adjusted to search for the optimized surface structure that minimizes the $R_{2}$ factor. For electrons with an energy range of 25-135 eV, the mean free path is about 5 to 10 Å. Use of this energy range means that our $\mu$LEED-IV curves are most sensitive to the surface structure of the top few atomic layers of our samples, i. e.  the buckling of the top atomic layer $b_{1}$, the thickness of the first phosphorene layer $z$’, the buckling of the bottom atomic layer $b_{2}$ and the Van der Waals gap between the top and second phosphorene layer w’, as demonstrated in Fig. 1 (d).

Multiple scattering theory and a muffin-tin potential model were implemented to calculate the LEED-IV curves 21, 22. We used computer codes from Adams et al. 20 which were developed from the programs of Pendry 21 and Van Hove and Tong 22. The utilization of the $R_{2}$ factor in Adams’ codes allows the relative intensities of the diffraction beams to be preserved during the optimization, which enhances the reliability of the surface structure determination. The phase shifts (a quantity describing the atomic scattering property 22) were calculated using the Barbieri/Van Hove phase shift calculation package 23. The muffin-tin radii for phosphorus atoms was set to r${}_{P}^{MT}$ = 2.099 a.u. and 12 phase shifts (L = 11) were used for the LEED-IV calculation. The in-plane lattice constant was $a_{1}$ = 3.313 Å and $a_{2}$ = 4.374 Å, the thickness of the phosphorene layer was z = 2.166 Å and the van der Waals distance between phosphorene layers was w = 3.071 Å for the bulk, as indicated in Fig. 1 9.

The mean-square atomic vibrational displacements <$${u^{2}}$$>_{T} for the P atoms were calculated individually according to the relation between Debye temperature $\theta$D and <$${u^{2}}$$>_{T} at the sample temperature of T=300 K for bulk BP and T=573 K for the FLP flakes using the following equation 21:

[TABLE]

where $m_{a}$ is the atomic mass, $\hbar$ is the Planck’s constant and $k_{B}$ is the Boltzmann constant. The Debye temperature $\theta$D was set to 550 K 24. The inner potential, V0+ $i$Vim, was set to be independent of energy. The real part V0 was initially set to 8 eV and adjusted through $\Delta$V0 during the fitting process while the imaginary part Vim was fixed at 6 eV.

The surface structural parameters $b_{1}$, $b_{2}$, w’ and z’ were varied in search for the optimum surface structure. Best-fit structural parameter values are listed in Table 1 and compared with previously reported results 10, 11. The calculated LEED-IV curves using optimized structural parameters match well with the experimental curves for both the BP crystal surface at 300 K and exfoliated FLP flake at 573 K, as shown in Fig. 3 (a)-(d). The minimized $R_{2}$ factors are 0.03 and 0.02, respectively. For comparison, the calculated IV curves (blue dashed lines in Fig. 3 (a)-(b)) using a flat, un-buckled surface, are distinctively different from our experimental results. For the freshly cleaved BP crystal surface, our results show that the top layer surface buckling $b_{1}$ is 0.22 Å and the second phosphorus atomic layer buckling $b_{2}$ is 0.27 Å. The thickness of the top phosphorene layer z’ is expanded by 5.3% from its bulk value of 2.166 Å. The van der Waals gap between the top and second phosphorene layer w’ is contracted by 8% from its bulk value of 3.071 Å. For the mechanically exfoliated flake of FLP at 573 K, the top and second layer buckling are 0.30 Å and 0.29 Å, respectively. The surface bucklings are slightly larger at 573 K than the BP crystal surface at 300 K. We attribute this increase of surface buckling to thermal surface expansion at elevated temperature. For the same reason, the top phosphorene layer z’ and the top van der Waals gap w’ are also slightly increased at 573 K compared to 300 K. z’ shows an expansion of 9.9% and w’ a contraction of 6.3%, with respect to their corresponding bulk values. Figure 3 (e)-(f) show plots of the reliability $R_{2}$ factor as a function of the surface buckling $b_{1}$ and the second atomic layer buckling $b_{2}$ for both of the investigated samples. Well-defined minima were observed for both cases. Along with the good agreement of experimental and calculated IV curves, both results give us strong confidence in our findings.

Our most striking result is that the BP surface buckling $b_{1}$ is one order of magnitude larger than the previously proposed theoretical values 10, 11, for both BP and FLP samples investigated. Note that the buckling extends to second atomic layer. Similar significant surface buckling has also been predicted for other group V thin film such as Bi and other similar elemental 2D materials such as silicene, germanene, by various first-pricinples studies. Specifically, Cahangirov et al. predicted that the buckling height for silicene to be 0.44 Å and 0.64 Å 25; our co-author, Sadowski et al., proposed the buckling of Bi thin film to be 0.5 Å 26, 27. In order to further support our surface structural results, first-principles calculations were utilized to investigate the origin of this surface reconstruction as well as its influence on the electronic properties of BP and FLP.

First-principles studies of the surface structure were carried out based on the framework of density functional theory (DFT) with projector augmented plane-wave (PAW) potential 28 as implemented in the Vienna ab initio simulation package (VASP) 29, 30, 31. The plane-wave functions expanded with an energy cutoff of 400 eV were employed throughout calculations. The exchange-correlation energy was described by generalized gradient approximation (GGA) in Perdew, Burke, and Ernzerhof (PBE) form 32. The k points in two-dimensional Brillouin zone (BZ) of the 1$\times$1 unit cell of monolayer BP containing 4 phosphorous atoms were sampled on a 16$\times$12 mesh. The van der Waals (vdW) interactions were also incorporated within the Tkatchenko-Scheffler method 33. In addition, we employed the Heyd, Scuseria, and Ernzerhof (HSE06) hybrid functional 34, 35 for the band structure calculations.

The structure of monolayer phosphorene and the top phosphorene layer of bulk BP (a six-layer supercell) were calculated and compared. Only very small structural differences were observed, $<$ 0.001 Å, between the atomic positions and and bond lengths of the monolayer phosphene and that of the top layer of bulk. This is expected for layered materials with weak van der Waals bonding in between adjacent layers. In order to simplify our calculations, we focus on single layer phosphorene. The thickness of the vacuum layer in each slab structure is more than 15 Å.

First, the defect-free monolayer phosphorene with different sizes were investigated. The lattice structure was optimized until the atomic force, both Hellmann-Feynman and vdW terms included, on each relaxed atom was less than 1 meV/Å. In an up to 8$\times$4 supercell, no buckling was found within the accuracy of the calculation. This result is reasonable since both BP bulk and monolayer structures have the insulating electronic structures with band gaps, and exposed surfaces do not bring about the electronic mismatch or additional dangling bonds. Surface reconstruction is thus not necessary in such a stable structure.

However, if an impurity, such as vacancy defect 11 or doping 3, 17 is induced on the surface, the situation changes completely. In fact, Liang et al. have recently observed vacancy defects on their freshly cleaved surfaces of BP crystals 11 using STM. Here, we introduced a single point defect into the monolayer phosphorene by removing one atom. Several supercells were calculated with their sizes ranging from 2$\times$4 to 8$\times$4. After the structure optimization, deviations of the atoms along out-of-plane direction were observed in all of these structures. As shown in the top view (middle panel) of Fig. 4, each supercell has 8 zig-zag rows, and the defect is located on the upper layer of row 3. The magnitude of buckling in each row is summarized in the bottom panel of Fig. 4. by calculating the standard deviations of the phosphorus atoms’ z components for each entire row. It is seen that the buckling is maximized in rows around the defect, and the maximum buckling ranges from 0.15 Å to 0.33 Å through all the supercell sizes under investigation. These calculations agree well with our experimental values of 0.22 Å to 0.30 Å. Although the buckling magnitude decays rapidly along the armchair direction, away from the row, on which the defect is located, no significant decay in the buckling magnitude was found in the zig-zag direction. Based on these results, it is concluded that the buckling is significantly enhanced near the point defect. It is anisotropic and long-range along the zig-zag direction while it is short-range along the armchair direction. The defect-induced buckling cannot be maintained in the armchair direction. This interesting insight agrees well with the recent experimental observation of an anisotropy in the surface density-of-state (DOS) on the BP surface by STM 11.

Intuitively one would expect that such long-range buckling would be reflected in the band structures as well. Thus the electronic structure of the 4$\times$4 supercell with a single vacancy was investigated and compared with that of a clean monolayer. According to the density-of-state results shown in Fig. 5 (a)-(b), the clean monolayer BP is insulating with a band gap of 1.5 eV, while an impurity state is present in the defect containing supercell across the Fermi level close to the top of valence states. A similar state was also observed by Zhang et al. 10 in their STM dI/dV measurement. This indicates the existence of the defect-induced hole-doping electronic structure in these defect structure, i.e. each phosphorus vacancy generates three dangling bonds that need to be saturated by more electrons. This suggests that the distortion of the lattice, such as buckling, appears in order to eliminate this instability of the electronic structure.

To better understand the relation between hole-doping and the surface structure of BP, the hole doped 2$\times$1 clean supercell structures with a tunable total electron number was optimized. As shown in Fig. 5 (c), buckling appears when the hole number exceeds 0.5 per 8 phosphor atoms. The magnitude of the buckling as well as the energy difference between the buckled and ideal structures increases rapidly with the rise of the hole number. In particular, the buckling reached 0.2 Å when the hole number is 0.6 per 8 atoms. Our first principle calculations thus show that the presence of defect induces hole doping on the clean BP surface, which in turn leads to lattice distortion and the surface buckling. It was confirmed experimentally that both undoped bulk BP 17 and FLP 3 are p-type semiconductors, but the origin of intrinsic p-type doping is unclear so far. Recently, Osada proposed that the edge state of finite bi-layer phosphorene might be the origin of the intrinsic hole-doping around the edge 36. Our DFT calculations, together with the experimental observation of BP surface reconstruction, strongly indicate that surface defect induced buckling is a likely explanation for the intrinsic hole-doping for both bulk BP and FLP.

To summarize, we determined the surface atomic structure of black phosphorus using dynamical selected area $\mu$LEED-IV in a high spatial resolution LEEM system. Two samples were studied at two different temperatures; freshly cleaved bulk BP at room temperature and mechanically exfoliated few-layers phosphorene at 573 K. The results show that for both samples, the surface buckling is one order of magnitude larger than previously proposed. The surface buckling on black phosphorus was found to be profound, and is close to the buckling of other similar elemental 2D materials such as silicene and germanene 25 and other group V thin film materials such as Bi 27. Specifically, the top surface buckling $b_{1}$ for the cleaved BP crystal surface at 300 K is 0.22 Å and for exfoliated FLP of about 10 nm thickness at 573 K is 0.30 Å. The slight increase in the buckling magnitude from 300 K to 573 K is most likely due to the thermal expansion. Similar phenomenon has been observed on other 2D materials, such as MoS2 12. Furthermore, due to the high sub-surface sensitivity of $\mu$LEED-IV, we found a similar buckling for the second phosphorus atomic layer. We denoted this buckling as $b_{2}$=0.27 Å and 0.29 Å for T=300 K and 573 K, respectively. Besides the surface buckling, our results show that the BP surfaces have an expansion-contraction surface relaxation. The thickness of the top phosphorene layer z’ is expanded by 5.3% in the cleaved BP crystal at 300 K and expanded by 9.9% in the exfoliated FLP at 573 K. The van der Waals gap between the first and second phosphorene layer w’ is contracted by 8% in cleaved BP crystal at 300 K and ccontracted by 6.3% in exfoliated FLP at 573 K.

We further confirmed our surface structural results using first-principles calculations and identified a vacancy defect driven mechanism as the cause of the surface buckling. The surface vacancy defect also introduces an impurity state in the band gap, which suggests that the vacancy defect induced surface buckling is most likely the reason why phosphorene materials are intrinisically p-type.

We have shown that $\mu$LEED-IV via LEEM, with its high surface sensitivity and selected area investigation ability, is a powerful technique to study atomic lattice of 2D materials. We believe that our results provide a clear picture of the fundamental surface properties of BP, including the detailed surface structure, the origin of its surface deformation, the electronic properties of phosphorene materials, and the origin of its intrinsic hole-doping property. We expect the significant surface buckling to potentially be even more pronounced in the monolayer phosphorene form, and having a significant impact on the electronic properties of monolayer phosphorene based devices.

{acknowledgement}

The experimental work presented here was carried out at the Center for Functional Nanomaterials, which is a U.S. DOE Office of Science Facility, at Brookhaven National Laboratory under Contract No. DE-SC0012704. The work of Z. D. and K. P. were supported by NSF DMR 1006863. First principles analysis by J.X. Y. and J. Z. was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award No. DE-SC0016424. In addition, the work of R.M. O. and W. J. were financially supported by the U.S. Department of Energy under Contract No. DE-FG 02-04-ER-46157. The work of J. H. and Y. D. K. were financially supported by NSF MRSEC program through Columbia in the Center for Precision Assembly of Superstratic and Superatomic Solids (NSF DMR-1420634).