Introduction

The discovery of graphene laid the foundations for many new areas of research. One of the most important foundation is the intensive investigation of two-dimensional (2D) atomic-layer systems, including graphene itself1,2,3,4,5, transition metal dichalcogenides (TMDCs)5,6,7,8, silicene9,10 and germanane10,11, as candidate materials for future electronics applications12,13,14,15,16,17,18,19. A high-performance device such as a field-effect transistor (FET) requires a moderate electronic band gap, a reasonably high carrier mobility of the channel material and excellent electrode-channel contacts4,5,8,13,14,15,16,17,18,19. Graphene offers extremely high mobilities, due to its very low carrier effective mass, and thus is considered to be a promising candidate for high-speed FET devices, but its intrinsic dispersion is gapless2,3,4,12,13,14. Despite extensive efforts following a wide variety of approaches to the problem of opening a gap in different graphene nanostructures, all devices to date have a relatively large ‘off’ current and thus a low ‘on–off’ ratio13. For this reason the emergence of monolayer TMDs has attracted substantial research interest, and in particular MoS2 has recently been used to fabricate a FET15. Unlike graphene, monolayer MoS2 is a direct-bandgap semiconductor with a carrier mobility of ~200 cm2 V−1 s−1 (ref. 15), improvable up to 500 cm2 V−1 s−1 (ref. 8), a value which is not unreasonable for applications but remains orders of magnitude lower than that of graphene13,14. Thus the entire community working on atomic-layer transport continues to search for a 2D material which is semiconducting, preferably with a direct gap, has high carrier mobility and has the potential to form excellent contacts with known electrode materials.

Here we present the theoretical discovery of a new category of layered direct-bandgap semiconductor, few-layer black phosphorus (BP), with high mobility, high in-plane anisotropy and linear dichroism. BP is an allotrope of phosphorus, whose crystal structure is a strongly folded honeycomb sheet with ‘troughs’ running along the y-axis (b direction). We demonstrate using density functional theory (DFT) calculations that few-layer BP systems, from the monolayer up to five-layer structures, are thermally stable with an interlayer interaction energy of −0.46 eV and a bandgap that falls exponentially with the thickness. Carrier mobilities at room temperature are high, and furthermore exhibit strongly anisotropic behaviour. The BP monolayer has an exceptionally high hole mobility, which we predict to lie in the range from 10,000 to 26,000 cm2 V−1 s−1, and which we explain by visualizing the wavefunctions. We find an explicit linear dichroism in the computed absorption spectra, whereby the positions of the lowest-energy absorption peaks for the two in-plane directions differ strongly. These results demonstrate that few-layer BP is a new category of 2D semiconductor with high potential for novel applications in nanoelectronics and optoelectronics.

Results

Geometric and electronic properties of bulk BP

Because no experimental data are available for the atomic structure and electronic properties of few-layer BP, we begin by considering bulk BP to refine the accuracy of our theoretical predictions. We have performed DFT calculations using a number of different functionals to gauge which provide the best fit to experiment; a full list may be found in Supplementary Table 1 and we summarize our findings here. If one considers only the lattice geometry then the PBE-G06 (refs 20, 21) and optB86b-vdW22 methods produce the best results. However, neither local density approximation-modified Becke–Johnson (LDA-mBJ)23,24 nor HSE06 (refs 25, 26) calculations based on these ‘best’ geometries can provide a satisfactory value for the bulk bandgap, as shown in Supplementary Fig. 1. To optimize both sets of properties simultaneously, we found that the optB88-vdW functional22,27, combined with either the LDA-mBJ or HSE06 method to predict the electronic bandstructure, gave the best fits. Details of choosing the best-fit functional are available in Supplementary Methods. Figure 1a shows the fully relaxed atomic structure of bulk BP, where the equilibrium lattice constants obtained from optB88-vdW are only 1–2% larger than the experiment28,29,30. The accompanying Brillouin zone (BZ) and electronic bandstructures, labelled mBJ (optB88-vdW) and HSE06 (optB88-vdW) are shown in Fig. 1b,c. Both combinations of methods predict that bulk BP is a semiconductor with a direct bandgap at the Z point of 0.31 eV (mBJ) or 0.36 eV (HSE06), both values being fully consistent with the experimental value of 0.31–0.35 eV (refs 31, 32, 33, 34).

Figure 1: Lattice and electronic structures of bulk black phosphorus.
figure 1

(a) Crystal structure of bulk BP marked with coordinate axes (x, y, z), lattice vectors (a, b, c) and structural parameters (R1, R2, θ1 and θ2). (b) Brillouin zone path of BP primitive cell. (c) Electronic bandstructures for bulk BP calculated with the HSE06 functional (red solid line) and the mBJ potential (blue dashed line), together with fitted effective masses along the Z–T′–A′, Z–Q and Z–G directions. At the right of the image, a zoomed-in plot shows the direct bandgap at Z. EVBM is the energy of valence-band maximum.

Full size image

At the bandgap wavevector, which is the Z point, one valence band (VB) and one conduction band (CB) disperse rather strongly along the Z–Q and Z–G directions (Fig. 1c), indicating very small effective masses. A fit of these bands using the nearly-free electron model gives effective carrier masses along the Z–Q direction which are rather small and similar, namely 0.12 m0 for electrons and 0.11 m0 for holes, whereas these along Z–G are slightly larger, taking the respective values 0.15 m0 and 0.30 m0 (all masses taken from HSE06). Considerably larger values are found along Z–T′–A′, where the carrier effective masses are 1.15 m0 and 0.71 m0, respectively. All of these results are fully consistent with the experimental values for bulk BP35. The corresponding mBJ results are within 0.04 m0 of the HSE06 values along Z–Q, Z–G and Z–T′–A′ except for the electron effective mass along Z–T′–A′, which is larger than the HSE06 value by 0.13 m0. Values for the direction Z–Q are similar to other high-mobility semiconductors, such as AlGaAs (0.1 m0)36, and four times smaller than those of MoS2 (ref. 37).

Geometric and electronic properties of few-layer BP

The similarity between the bandstructures of bulk and monolayer BP38, and the low effective masses of bulk BP35, suggest that few-layer BP is likely to be a high-mobility, direct-bandgap 2D semiconductor. We have thus performed direct calculations of the geometric and electronic properties for few-layer BP systems. Table 1 summarizes the changes of geometrical properties as a function of the layer number from 1 to 5. The lattice parameter a increases by 0.11 Å on passing from bulk to monolayer BP, whereas b grows only by 0.02 Å. There is an abrupt reduction of a between the monolayer and the bilayer, which we believe to be the primary consequence of interlayer interactions (0.46 eV) in the bilayer. The significant stretching of a in few-layer BP systems is caused almost exclusively by an increase of the bond angle θ1 but not of the bond lengths.

Table 1 Structural information for few-layer BP.
Full size table

Results from our bandstructure calculations for the five few-layer BP systems are shown in Fig. 2. In the monolayer, the original Z point of the bulk BZ folds back to the Gamma (G) point, so that the original Z–Q and Z–T′–A′ directions of the bulk BZ project onto the G–X and G–Y directions of the monolayer, which correspond respectively to the a and b directions in real space (Fig. 2a,b). Results obtained from the mBJ method are quantitatively the same as those obtained from HSE06. Supplementary Figure 2 and Supplementary Table 2 show that the predicted bandgap is rather less sensitive than the unit-cell dimensions to the functional used for optimizing the atomic structures, with similar values emerging for all few-layer BP systems from the best-performing functionals. However, other physical quantities, notably the location of band maximum (minimum) and the carrier effective mass (and hence the mobility) are in fact very sensitive to the choice of the functional, and hence we have focused primarily on our optB88-vdW results for the discussion of these. Monolayer BP is indeed a direct-bandgap semiconductor, as shown in Fig. 2c, with the gap value of 1.51 eV obtained at the G point. Phosphorus is a heavier element than carbon and therefore should have stronger spin-orbit coupling (SOC) in its 2D forms than graphene. We have considered full SOC effects in calculating the electronic bandstructure of bulk BP. Supplementary Figure 3 shows that inclusion of SOC terms has no appreciable effect on the primary features of the bandstructure, indicating that phosphorus is ‘not sufficiently heavy’ to cause any qualitative changes. Only a very small separation of the formerly fourfold degenerate bands, into two band pairs, can be found from X along N (or M for few-layer BP) to Y, which reaches a maximum of 23 meV around the N (M) point. Along G–X and G–Y, these bands are already separated by bandstructure effects and SOC causes no further splitting. That the spin-orbit splitting energy is so small is a consequence of the small effective nuclear charge (Zeff) of the P atom and the weak variation of the charge gradient in an elemental system such as BP.

Figure 2: Electronic structures of few-layer BP.
figure 2

(a,b) Top view of the atomic structure of the monolayer and the associated Brillouin zone. (d,e) Side views of the atomic structure of the bilayer. (c,f) Bandstructures of monolayer and bilayer BP calculated with the HSE06 functional (red solid lines) and the mBJ potential (blue dashed lines), respectively. Two valence (VB1 and VB2) and two conduction states (CB1 and CB2) are marked in panel f. (g) Spatial structure of wavefunctions for the four marked states illustrated in the xz and yz planes using an isosurface of 0.0025 e Å−3. (h) Evolution of the direct bandgaps as a function of the sample thickness. Functionals used for structural optimization are shown in parentheses. Bandgap values are marked for the monolayer system, for the extrapolation of our results and for real bulk BP.

Full size image

When two monolayers are combined to form a bilayer (Fig. 2d,e), the gap is reduced to 1.02 eV and two additional bands emerge around the gap at the G point. Together with the original VB and CB, we denote these bands as VB1, VB2, CB1 and CB2, as shown in Fig. 2f. In real space, the states from these four bands with k close to G are extended throughout the bilayer, as shown in Fig. 2g. States VB1 and VB2 differ in the interlayer region although they share the same origin in terms of the atomic orbitals. A clear bonding-like feature is visible in the interlayer region (marked by red rectangles) for VB1, which lies lower in energy, whereas VB2 shows an anti-bonding feature. Similar bonding and anti-bonding features are also found in CB1 and CB2, where they are observed not between the layers but across the troughs. These features indicate that wavefunction overlap, rather than van der Waals effects, plays the primary role in mediating the interlayer interaction, which explains the abrupt reduction of the lattice constant a from monolayer to bilayer. This interlayer interaction introduces band dispersions of VBs and CBs in the direction normal to the layers, which leads to the reduction of the bandgap (by 0.5 eV) from monolayer to bilayer.

The thicker the few-layer BP system, the stronger is the interlayer interaction and thus the larger is the dispersion of VBs (CBs), resulting in a smaller bandgap. Thus the gap falls continuously on adding more layers to the bilayer, reaching 0.59 eV in five-layer BP, as illustrated in Fig. 2h. We fitted these values by an exponential decay relation and found that the corresponding bulk gap obtained by extrapolation is 0.53 eV, a value 0.17 eV larger than the one we calculated for bulk BP. We ascribe this difference to the elongated lattice parameter a in few-layer BP and results calculated with a series of constrained a values confirm our expectation; these results imply that the bandgap is very sensitive to the lateral strain along the a (x) direction and could be modulated by varying interatomic separations and angles. During the preparation of this manuscript, we became aware that a strain-induced bandgap modulation has recently been demonstrated theoretically39.

The small effective masses remain in all of the few-layer BP systems (Table 2). For the G–X direction in the monolayer, equivalent to Z–Q in the bulk, the carrier effective masses are 0.15 m0 (hole) and 0.17 m0 (electron), only 0.04–0.05 m0 larger than those in bulk BP. It is remarkable that in the G–Y direction, equivalent to Z–T′–A′ in the bulk, the VB appears to be nearly flat close to the G point, with an effective mass of 6.35 m0 (see Fig. 2c), nine times its bulk value of 0.71 m0. By contrast, the effective mass for the CB is 1.12 m0, very close to its bulk value of 1.15 m0. Unlike the bandgap, only the hole effective mass along G–Y shows a strongly layer-dependent evolution, in that it decreases from 6.35 m0 for the monolayer to 1.81 m0 for the bilayer and eventually to 0.89 m0 in the five-layer system, close to the bulk value of 0.71 m0.

Table 2 Predicted carrier mobility.
Full size table

Carrier mobility

The electronic properties of few-layer BP are governed to a large extent by the carrier mobilities, which are in turn strongly influenced (but not solely determined) by their effective masses. Thus we also provide theoretical predictions for the carrier mobilities, in both x and y directions, in few-layer BP systems by applying a standard 2D model as discussed under Methods below. We apply a phonon-limited scattering model, in which the primary mechanism limiting carrier mobility is scattering due to phonons36,40,41,42. Two properties of the few-layer BP lattice, namely the deformation potential E1 and the elastic modulus C2D in the propagation direction of the longitudinal acoustic wave, are then the most relevant factors determining the mobility43,44. These quantities were calculated using the optB88-vdW functional and the data are presented in Table 2. If few-layer BP is cut into ribbons, a one-dimensional model should be adopted for predicting the carrier mobilities. We present the results and discussion for a 10-nm-width BP ribbon in the Supplementary Discussion and Supplementary Table 3.

Our predicted mobilities for the five few-layer BP systems (Table 2) are in general moderately large (hundreds to thousands of cm2 V−1 s−1), moderately anisotropic and asymmetric between electrons and holes, with the holes being more mobile in both directions. Concerning directional anisotropy, electron and hole mobilities for the x direction are nearly twice the values along y for holes and four times as large for electrons. However, the monolayer is exceptional in every way. The electron mobility along x is nearly 14, rather than four, times the value along y, that is, 1,100–1,140 cm2 V−1 s−1 compared with ~80 cm2 V−1 s−1, which is consistent with the trend of the other few-layer systems. To our great surprise, however, the hole mobility along x is 16–38 times smaller than that along y, that is, 640–700 cm2 V−1 s−1 versus 10,000–26,000 cm2 V−1 s−1, making y the direction of higher hole conductivity. This extraordinarily large value for monolayer BP is a consequence of the extremely small deformation potential, E1y=0.15±0.03 eV, and occurs despite the fact that the carriers are very heavy (6.35 m0). This ‘monolayer exception’ should be valuable in applications such as the separation of electrons and holes.

The value of E1y for holes in a monolayer is very striking in that it is an order of magnitude smaller than typical values of E1, which are 5.0 eV for graphene44, 3.9 eV for MoS2 (ref. 37), and 3.7 eV for h-BN40 and 7 eV for AlGaAs quantum wells36. In few-layer BP systems, E1y for holes has a conventional value 2.97±0.18 eV in a five-layer structure, decreasing smoothly to 1.63±0.16 eV in a bilayer. These results can be explained rather well by the form of the VB wavefunctions shown in Fig. 2g, where interlayer and stacking-induced intra-layer overlap increases with layer number. These effects are completely absent in monolayer BP, making the VB wavefunction quite isolated along y, as it is also for the VB2 (anti-bonding) state of a bilayer highlighted by the green rectangle in Fig. 2g. Thus small structural deformations due to longitudinal phonon oscillations along y have almost no effect on this wavefunction and cause little change in its energy, resulting in the very small deformation potential. The situation is quite different for the CB (electron) wavefunction, which has E1x=2.72±0.02 eV in the monolayer and 5.02±0.02 eV in the bilayer. The in-plane bonding and anti-bonding features of the CB wavefunctions in few-layer BP systems cause a substantial enhancement of the energy-level perturbations due to structural deformations along y. The monolayer exception aside, all E1 values for CB states are rather large (5.0–7.7 eV), for both x and y directions, and are almost independent of the number of layers, a result which can be explained by the strong CB wavefunction overlap, especially along y (Fig. 2g). For the VB states, the potentials are somewhat smaller (2–3 eV), precisely because their wavefunctions have less overlap between different P atoms.

Optical absorption spectra and linear dichroism

We have also predicted the optical absorption spectra of few-layer BP systems by computing the dielectric function. Two absorption spectra are shown in Fig. 3a,b for light incident along the z direction and linearly polarized in the x and y directions, respectively; while the results for incident light polarized in the z direction may be found in Supplementary Fig. 4. These results demonstrate a strong linear dichroism: for a dielectric polarization in the x direction, the band edge of the first absorption peak is found at the bandgap and thus falls rapidly with the thickness of the sample (Fig. 3a). By contrast, with y-polarized light this peak is found at 3.14 eV in the monolayer and its position falls only slightly with thickness, remaining at 2.76 eV in bulk BP (Fig. 3b). From the symmetries of the wavefunctions in Fig. 2g, clearly the (odd) dipole operator connects the VB and CB states for x-polarization, allowing the direct-bandgap process, but this is symmetry-forbidden for y-polarization and the transition occurs between VB and CB states elsewhere in the BZ.

Figure 3: Optical absorption spectra.
figure 3

(a,b) Optical absorption spectra of few-layer BP for light incident in the c (z) direction and polarized along the a (x) and b (y) directions, respectively. Black dashed lines show an approximate linear fit used to estimate the band edges for the first absorption peak, which are highly anisotropic between x and y. The band edge drops rapidly with sample thickness for x (from 1.55 eV for the monolayer (red) to 0.60 eV for five layers (orange)) but only slightly for y (from 3.14 eV to 2.90 eV). (c) Schematic illustration of a proposed experimental geometry to determine the orientation of few-layer BP structures using optical absorption spectroscopy, and thus to utilize the anisotropic electronic properties of BP. The light is linearly polarized in a chosen orientation and near-normally incident on the sample. The sample should be rotated to identify the a or b direction by monitoring the absorption signal, after which the source and drain electrodes, denoted as ‘S’ and ‘D’ respectively, may be deposited to fabricate an FET device.

Full size image

Discussion

High carrier mobilities would be expected in BP systems due to the small effective masses of 0.1–0.2 m0. However, our calculations reveal that the key issue is not the effective mass but the very small elastic modulus C2D: all the BP systems we have studied are ‘very soft’ in the x direction, due to the fact that bond angles are changed easily (with little effect on bond lengths). When compared to other monolayer systems, the modulus of BP for the x direction is 29.0 J m−2, a value respectively 2, 10 and 11 times smaller than MoS2 (ref. 37), h-BN40 and graphene44. These results indicate that much higher mobilities could be most likely achieved if few-layer BP systems could be modified to become stiffer, that is, to have a larger elastic modulus for the x direction, for example by external stress, doping or substrate confinement. However, possible changes of the deformation potential as a result of modifying the elastic moduli are not yet clear and require further investigation. The modulus for the y direction is significantly larger than that for x, but the higher carrier effective masses, result in carrier mobilities being smaller than those along x. Both the bandgap and the mobility depend strongly on the structural properties and this offers a number of possibilities for further control over the electrical and mechanical properties of few-layer BP.

Further, the linear dichroism we find allows an electronic determination of the sample orientation. We propose an experimental setup to deduce the orientation of few-layer BP systems using optical spectroscopy, as illustrated in Fig. 3c. If the incident light is linearly polarized in a chosen orientation and is near-normally incident, the absorption spectrum should vary when the sample is rotated and such that the a and b directions can be identified by monitoring the absorption signal. Once the sample orientation is determined, it is far easier to fabricate electrodes or gates utilizing the highest-mobility direction of a few-layer BP sample in an FET-type device.

In summary, we have shown theoretically that few-layer BP is a novel category of 2D semiconductor offering a direct bandgap, high carrier mobility and high transport anisotropy, all of which are tunable by controlling the layer thickness. The bandgap decreases from 1.51 eV for a monolayer to 0.59 eV for five-layer BP, a set of values fitting well into the gap between graphene nanoribbons and TMDCs. The anisotropy of electric and optical properties is another unique feature distinguishing BP from other 2D materials, such as silicene9,10, h-BN5,45 and germanane10,11. Higher conductivity is generally found in the direction oriented perpendicular to the troughs (the x direction) and holes are usually more mobile than electrons. For holes, the mobility in the x direction increases from 600 cm2 V−1 s−1 for a monolayer to over 4,000 cm2 V−1 s−1 for five-layer BP. We found that monolayer BP is a very special material, in which we predict an extremely high hole mobility of 10,000–26,000 cm2 V−1 s−1. All of our mobility results can be understood rather well from the form of the real-space wavefunctions, which show a much stronger interlayer coupling in few-layer BP than in graphene or TMDCs. All these results make few-layer BP a very promising candidate material for future applications in electronics and optoelectronics.

Note added in proof: during the review process for this manuscript, we became aware that few-layer BP samples have very recently been exfoliated from bulk BP and fabricated into electronic devices46,47. Measurements on these samples verify the high hole mobilities46,47 and transport anisotropy46 we predict. Although preliminary optical characterizations have been performed46,47, linear dichroism has not yet been explored.

Methods

DFT calculation

DFT calculations were performed using the generalized gradient approximation for the exchange-correlation potential, the projector augmented wave method48,49 and a plane-wave basis set as implemented in the Vienna ab-initio simulation package50. The energy cutoff for the plane-wave basis was set to 500 eV for all calculations. Two k-meshes of 10 × 8 × 4 and 10 × 8 × 1 were adopted to sample the first BZ of the conventional unit cell of bulk and few-layer BP, and the mesh density of k points was kept fixed when calculating bandstructures using primitive cells. In optimizing the system geometry, van der Waals interactions were considered by the vdW-DF level with the optB88 exchange functional (optB88-vdW)22,27. The shape and volume of each supercell were optimized fully and all atoms in the supercell were allowed to relax until the residual force per atom was <0.001 eV Å−1. Electronic bandstructures were calculated by the mBJ23,24 and hybrid functional (HSE06)25,26 methods based on the atomic structures obtained from the full optimization by optB88-vdW. The charge gradient used in the mBJ calculations, c=1.1574, was extracted from the bulk value.

Carrier mobility calculation

In 2D the carrier mobility is given by the expression40,41,42

where is the effective mass in the transport direction and is the average effective mass determined by . The term E1 represents the deformation potential constant of the valence-band minimum for hole or conduction-band maximum for electron along the transport direction, defined by . Here ΔVi is the energy change of the ith band under proper cell compression and dilatation (calculated using a step of 0.5%), l0 is the lattice constant in the transport direction and Δl is the deformation of l0. The elastic modulus C2D of the longitudinal strain in the propagation directions (both x and y) of the longitudinal acoustic wave is derived from (E−E0)/S0=C(Δl/l0)2/2, where E is the total energy and S0 is the lattice volume at equilibrium for a 2D system. All structural properties in the calculation of carrier mobilities were obtained from optB88-vdW and properties related to the electronic structure were computed with the HSE06 functional. The temperature used for the mobility calculations was 300 K.

Uncertainty of deformation potential

In calculations of the electron and hole mobilities, the deformation potential is derived from linear fits to the respective energies of the conduction-band maximum and the valence-band minimum as functions of the lattice dilation or compression. Typical fitting results for E1y (electron state) and E1x (hole state) are shown in Supplementary Fig. 5a–c. The slope represents the deformation potential, which usually has an unavoidable uncertainty accounted for in the standard fitting error. Supplementary Figure 5d summarizes the relative errors for the deformation potentials of electrons and holes in few-layer BP for the x and y directions. For a given direction and type of carrier, the error is generally smallest for monolayer BP, while that for five-layer BP is the largest, except for the hole deformation potential along y. The error in E1y for hole states appears rather larger than the others, especially for monolayer BP where it exceeds 20%, presumably due to the small absolute value of E1y for this case. Supplementary Fig. 5c shows details of the fit. The absolute value of E1y for hole states is only 0.15 eV, an order of magnitude smaller than the other energies; an error of 1 meV, believed to be the accuracy limit of DFT, may then lead to a fitting error of more than 10%. In this respect, a fitting error of 20% is definitely reasonable, and suggests that our calculation is highly accurate.

Absorption spectra calculation

Absorption spectra were calculated from the dielectric function using the expression A(ω)=1−e−α(ω)·Δz, where is the absorption coefficient, is the index of refraction, ε1 and ε2 are the real and imaginary parts of the dielectric function, ω is the light frequency, c is the speed of light in vacuo and Δz represents the unit-cell size in the c direction. The electronic structures were obtained from HSE06 results and the k-mesh was doubled in calculating dielectric functions. Excitonic contributions were not considered in our calculations. The total number of bands considered was set to be twice that used in the total-energy and bandstructure calculations. Because the dielectric function is a tensor, the absorption spectra along the three directions a (x), b (y) and c (z) were obtained separately.

Additional information

How to cite this article: Qiao, J. et al. High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus. Nat. Commun. 5:4475 doi: 10.1038/ncomms5475 (2014).