Tight Binding Parametrization of Few-layer Black Phosphorus from First-Principles Calculations

Tight Binding Parametrization of Few-layer Black Phosphorus from First-Principles Calculations


W̱e employ a tight-binding parametrization based on the Slater Koster model in order to fit the band structures of single-layer, bilayer and bulk black phosphorus obtained from first-principles calculations. We find that our model, which includes or parameters depending on whether overlap is included or not, reproduces quite well the ab-initio band structures over a wide energy range, especially the occupied bands. We also find that the Inclusion of overlap parameters improves the quality of the fit for the conduction bands. On the other hand, hopping and on-site energies are consistent throughout the different systems, which is an indication that our model is suitable for calculations on multilayer black phosphorus and more complex situations in which first-principles calculations become prohibitive, such as disordered systems and heterostructures with a large lattice mismatch. We also discuss the limitations of the model and how the fit procedure can be improved for a more accurate description of bands in the vicinity of the Fermi energy.


I Introduction

Since the experimental discovery of graphene in 2004, a strong scientific effort has been employed in the study of this material. Many interesting structural and electronic properties were uncovered, with potential applications in the near future Neto et al. (2009); Novoselov et al. (2012). Despite its large carrier mobility in comparison with current silicon based devices, graphene is a zero-gap semiconductor, which undermines its application in electronic devices such as field effect transistors, which require a band gap. As such, a lot of research has been done in finding ways to open a band gap in graphene. One of the most promissing routes so far is the application of external electric fields or doping in multilayer systems Ohta et al. (2006); Lui et al. (2011); Zhang et al. (2010); Guinea et al. (2006).

The success of graphene research and the search for a band gap also led to the study of other materials that share some structural similarities to graphene, the so-called 2D materials. These materials form a single layer of thickness of one up to a few atoms, and can be metallic, semiconducting or insulating, thus being suitable for many different applications Ferrari et al. (2015); Butler et al. (2013). Moreover, the weak van der Waals interactions between such layers allows the possibility of stacking them in many different ways in order to target specific electronic and optical properties, resulting in a set of materials now known as van der Waals heterostructures Geim and Grigorieva (2013); Novoselov et al. (2016). Among such materials, few-layer black phosphorus (BP) is a promissing candidate for future applications. Like graphite, it is a natural layered material that can be exfoliated down to a single layer, which is known as phosphorene in analogy to graphene Li et al. (2014); Lu et al. (2014). Fig. 1 shows the crystal structure of bulk and single layer BP. We can see that each layer is composed of puckered zigzag chains of P atoms, a structure that leads to highly anisotropic transport and optical properties Qiao et al. (2014); Liu et al. (2014); Xia et al. (2014); Wang et al. (2015). Moreover, theoretical predictions and experimental observations have confirmed that few-layer BP is a semiconductor with a band gap that can be tuned by the number of stacked layers, ranging from eV for single-layer down to eV in the bulk, which is a very desirable property for applications in optoelectronic devices Li et al. (2017); Tran et al. (2014); Liu et al. (2014).

Figure 1: Left: Crystal structure of black phosphorus, showing the puckered zigzag chains in each layer. Right: Top view of a single layer, also known as phosphorene. The dashed rectangle indicates an unit cell of the layer.

From a theoretical point of view, the electronic and structural properties of few-layer BP can be accurately determined by first principles calculations, such as those based on Density Functional Theory (DFT) and the GW approximation for quasiparticle corrections Qiao et al. (2014); Tran et al. (2014); Elahi et al. (2015); Dai and Zeng (2014); Yan Li and Li (2014). However, such methods quickly become prohibitive with an increasing number of atoms in the unit cell, as is the case in calculations for transport properties in disordered systems and heterostructures with a large lattice mismatch. In such situations, semi empirical tight-binding (TB) calculations are often preferred, but they require an adequate parametrization, which may be obtained from theoretical calculations such as DFT and GW, or experimental data. There are a few available parametrizations so far, with different degrees of success in describing the electronic properties of this material. In one class of parametrizations, based on projection techniques using Wannier functions and GW calculations, the proposed models include a single -like orbital per atom Rudenko and Katsnelson (2014, 2015). Therefore, they do not fully capture the physics of bonding in this system, due to the staggered nature of each layer. Such models include up to parameters and are successful in describing the low energy properties of few-layer BP, specially near the point of the Brillouin Zone, where the bands are mostly of character. However, they do not fully describe these bands away from the point, since they also contain contributions from , and orbitals. Moreover, they describe only a limited number of bands, not taking into account most of the higher energy states. In a second class of parametrizations, based on atomic orbitals, the proposed models include all and orbitals, but employ a large number of parameters and can become very complex Páez et al. (2016). The reason for that is that the hopping between a pair of sites is parametrized individually (considering the symmetries of the structure), so they do not obey an analytical functional form with distance. Therefore, the number of parameters increases with the number of neighboring interactions considered, reaching up to 58 parameters for 8th neighbors. Such models can describe the valence and conduction bands of few-layer BP with a greater overall accuracy than the previous class of models, but they also fail to fully describe the higher energy bands, since the parameters are tailored to fit the low-energy states.

In this work, we provide a simple, yet reliable, tight-binding parametrization based on the Slater-Koster model for the electronic structure of few-layer BP, obtained from first-principles calculations. Our model includes and orbitals for each site and an exponential decay behavior for the hopping and overlap integrals, resulting in a number of parameters that range from up to , depending on whether overlap parameters are included or not. The number of neighboring interactions can be controlled by a cutoff radius, whose value is chosen after convergence tests. The main advantages of this model, besides its simplicity, are that the associated TB parameters have a clear physical interpretation in terms of chemical bonding and that we can fit a larger number of bands in a wide energy range, as we describe below. Our paper is organized as follows: in the next section we present our methodology, including the first-principles calculations used to obtain the band structures, our tight-binding model used to fit them and the fitting procedure. In Sec. III, we present our results by discussing calculations with and without overlap integrals, different energy ranges and how the resulting TB parameters compare with previous parametrizations. We also discuss the role of the basis choice in the first-principles calculations and the consistency of our results for single-layer, bilayer and bulk BP. Finally, in Sec. IV we present our conclusions and discuss the suitability of our model for multilayer systems and more complex situations. We also discuss how the present model can be improved in order to reach a better fit quality.

Ii Methodology

ii.1 First-Principles Calculations

In the first step of our calculations, we perform first-principles calculations for single-layer, bilayer and bulk BP based on Density Functional Theory (DFT), as implemented on the SIESTA code Hohenberg and Kohn (1964); Kohn and Sham (1965); Soler et al. (2002). We employ a double-zeta-polarized (DZP) pseudoatomic orbital basis to expand the wavefunctions. We have also perfomed calculations with a simpler single-zeta basis (SZ), which bears a closer correspondence to the TB model used for the fits. We have found similar results for both basis sets, so we only discuss the results of the DZP basis here. A PBE exchange-correlation functional is used for the electron-electron interactions and norm-conserving Troullier-Martins pseudopotentials are employed for the ion-electron interactions Perdew et al. (1996); Troullier and Martins (1991). The Brillouin Zone is sampled using a Monkhorst-Pack k-point grid for single-layer and bilayer BP and a for bulk BP Monkhorst and Pack (1976). For single-layer and bilayer BP, a vacuum distance of Å is used in order to isolate the slabs from their periodic images.

Since experimental values for single-layer and bilayer BP are not yet available and theoretical predictions present variations depending on the choice of basis, exchange-correlation functional and inclusion of van-der-Waals interactions Qiao et al. (2014); Elahi et al. (2015); Dai and Zeng (2014); Yan Li and Li (2014), we fix all structural parameters to the experimental values of bulk BP in all cases Brown and Rundqvist (1965). Finally, since DFT underestimates energy gaps, we perform a rigid shift of all conduction bands in order to reproduce the quasiparticle gaps obtained from GW calculations (a procedure also known as scissors shift) Tran et al. (2014). It is found in many cases that such a shift is the main effect of the quasiparticle corrections in semiconductors, so it is enough for the purposes of this work Hybertsen and Louie (1986).

ii.2 Tight-Binding Calculations: Slater-Koster Model

In our TB model, we use an atomic basis with , , and orbitals for each atom, resulting in a hamiltonian in k-space. Following the prescription of the Slater-Koster model within the two-center approximation, any hopping or overlap integral between orbitals centered at different atomic sites can be expressed in terms of their relative orientations and four integrals related to and bonding Slater and Koster (1954).

The simplest case is the hopping between two equivalent orbitals separated by a displacement vector pointing from one to another. It is given simply by and does not depend on the orientation of the bond. The hopping between a orbital and a () orbital separated by a displacement vector from to is given by:


where is the direction cosine of with respect to the corresponding cartesian direction and is the hopping between a orbital and a orbital directed along the bond. Similarly, the hopping between and orbitals () is given by:


where is the hopping between orbitals directed along the bond and is the hopping between orbitals directed perpendicular to the bond and parallel to each other. For two equivalent orbitals centered at different sites, the expression is:


Therefore, all hoppings can be expressed in terms of direction cosines and the amplitudes and calculated at the bond length. These amplitudes decrease as the distance between the orbitals increase. Since the atomic wavefunctions decay exponentially for large distances away from the nuclei, we model the hopping amplitudes in the same way:


In the above equation, is an index that labels the type of hopping. is a reference distance at which the hopping amplitudes are calculated, which we choose as the nearest neighbor distance: Å. is a decay distance and is a cutoff distance. We set to Å, which generates well converged band structures according to our tests.

With this model, we have parameters related to hopping: four amplitudes , , and and four decay distances , , , . By setting the on-site energy of the orbitals to zero and including that of the orbitals in the model (), we end up with a minimal model with a total of adjustable parameters. We may also include overlap parameters following the same prescription outlined above, by adding up to four corresponding amplitudes and decay distances, resulting in a model with a maximum of adjustable parameters. In the next section, we present results for both the minimum model and the model with all overlap parameters.

ii.3 Parametrization: Least-Squares Minimization

In order to fit the Tight-Binding band structures to the corresponding DFT results, we use a least-squares minimization. In this method, we consider the TB band structure as a function of the Slater-Koster parameters described above and minimize the following function:


where and are the eigenvalues and and are the corresponding Fermi energies from DFT and TB calculations, respectively. The summation in this equation can be restricted to a given set of bands, k-points or energy range, so the number of data points can be different in each case. For a full-adjustment, in which all bands and k-points are included, . We have considered different sets of k-points and energy ranges, as we discuss in the next section.

The minimization of was carried on using Powell’s minimization algorithm Press et al. (1996). The iterative procedure stops when the relative value of differs by less than between two consecutive iterations. We find in our calculations that the minimized value of this function depends on the initial choice of parameters, especially when overlap is included. This is an indication that this functions may contain different local minima, so we have tested several different starting points and looked for the best fit.

Iii Results

iii.1 Full Adjustment

In Fig. 2, we show a comparison between the DFT band structures of single-layer, bilayer and bulk BP and the corresponding TB full adjustments without (left) and with (right) overlap parameters in the model. For bulk BP, we have used a conventional rectangular unit cell with four atoms in order to facilitate the comparison with the single-layer and bilayer Brillouin Zones. As we can see, the minimum TB model without overlap parameters gives a good overall description of the DFT band structure, especially the occupied bands. The inclusion of overlap parameters in the model improves the description of the unoccupied band and fixes a few avoided crossings, leading to a better quality fit.

Figure 2: Comparison between DFT band structures (black) and corresponding Tight-Binding fits (red) for single-layer (top), bilayer (middle) and bulk Black Phosphorus (bottom). The left (right) panels correspond to a model without (with) overlap parameters. The Fermi energy is set to zero in all cases (horizontal dashed line).

Table 1 shows the fit parameters for each case and the corresponding standard deviation of the adjusment, which is the minimized value of . The first three (next three) rows of the table correspond to a model without (with) overlap parameters, which we label TB1 (TB2). The last row contains reference parameters from Harrison’s model, which describes the hopping between and orbitals as a decay Harrison (1999). We can see that the hopping parameters are consistent throughout single-layer, bilayer and bulk BP, which is an indication that this model is suitable for calculations on multilayer BP and more complex situations. When overlap parameters are included in the model, the adjusted values for the hopping amplitudes and decay constants change, but they remain consistent throughout the three cases. Overlap amplitudes are found to be very small, with the exception of , which may be an indication that the fit procedure is less sensitive to these parameters. Therefore, a model with all hopping parameters and only and might give a good description of the electronic structure.

\arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash
Without Overlap (TB1)
1 \arraybackslash-1.59 \arraybackslash2.39 \arraybackslash4.03 \arraybackslash-1.14 \arraybackslash-8.80 \arraybackslash0.33 \arraybackslash0.53 \arraybackslash0.58 \arraybackslash0.53 \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash0.54
2 \arraybackslash-1.56 \arraybackslash2.43 \arraybackslash3.91 \arraybackslash-1.18 \arraybackslash-8.59 \arraybackslash0.36 \arraybackslash0.54 \arraybackslash0.59 \arraybackslash0.44 \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash0.55
Bulk \arraybackslash-1.54 \arraybackslash2.37 \arraybackslash3.63 \arraybackslash-1.25 \arraybackslash-8.44 \arraybackslash0.35 \arraybackslash0.59 \arraybackslash0.58 \arraybackslash0.41 \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash0.58
With Overlap (TB2)
1 \arraybackslash-2.02 \arraybackslash2.76 \arraybackslash3.99 \arraybackslash-1.07 \arraybackslash-7.34 \arraybackslash0.43 \arraybackslash0.55 \arraybackslash0.76 \arraybackslash0.51 \arraybackslash0.06 \arraybackslash0.01 \arraybackslash0.04 \arraybackslash0.14 \arraybackslash0.04 \arraybackslash6.37 \arraybackslash0.00 \arraybackslash0.36 \arraybackslash0.45
2 \arraybackslash-2.01 \arraybackslash2.63 \arraybackslash3.91 \arraybackslash-1.09 \arraybackslash-8.13 \arraybackslash0.48 \arraybackslash0.53 \arraybackslash0.71 \arraybackslash0.36 \arraybackslash \arraybackslash0.01 \arraybackslash0.00 \arraybackslash0.12 \arraybackslash0.02 \arraybackslash0.20 \arraybackslash0.00 \arraybackslash0.30 \arraybackslash0.47
Bulk \arraybackslash-1.75 \arraybackslash2.75 \arraybackslash3.61 \arraybackslash-1.07 \arraybackslash-7.38 \arraybackslash0.48 \arraybackslash0.56 \arraybackslash0.73 \arraybackslash0.35 \arraybackslash0.04 \arraybackslash0.01 \arraybackslash0.00 \arraybackslash0.15 \arraybackslash \arraybackslash0.00 \arraybackslash0.00 \arraybackslash0.40 \arraybackslash0.50
Har \arraybackslash-2.03 \arraybackslash2.19 \arraybackslash3.42 \arraybackslash-0.97 \arraybackslash-9.68 \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash-
Table 1: Tight-Binding paramters for a full adjustment of the DFT band structure, as shown in Fig. 2. Hopping and overlap amplitudes are calculated at nearest neighbor distance ( Å) and are in eV. Decay distances are in Å and is the standard deviation of the adjusment in eV. The last line contains reference parameters from Harrisons’ model Harrison (1999). Zero values mean that they are smaller than in the corresponding units.
single-layer bilayer bulk
\arraybackslashDFT \arraybackslashTB1 \arraybackslashTB2 \arraybackslashRef. \arraybackslashDFT \arraybackslashTB1 \arraybackslashTB2 \arraybackslashRef. \arraybackslashDFT \arraybackslashTB1 \arraybackslashTB2
\arraybackslash0.13 \arraybackslash0.40 \arraybackslash0.37 \arraybackslash0.15 \arraybackslash0.10 \arraybackslash0.37 \arraybackslash0.32 \arraybackslash0.15 \arraybackslash0.15 \arraybackslash0.34 \arraybackslash0.26
\arraybackslash0.13 \arraybackslash0.44 \arraybackslash0.40 \arraybackslash0.17 \arraybackslash0.10 \arraybackslash0.44 \arraybackslash0.32 \arraybackslash0.18 \arraybackslash0.17 \arraybackslash0.37 \arraybackslash0.25
\arraybackslash3.42 \arraybackslash2.37 \arraybackslash1.58 \arraybackslash6.35 \arraybackslash1.71 \arraybackslash1.90 \arraybackslash1.35 \arraybackslash1.81 \arraybackslash- \arraybackslash- \arraybackslash-
\arraybackslash1.14 \arraybackslash2.37 \arraybackslash- \arraybackslash1.12 \arraybackslash1.14 \arraybackslash4.73 \arraybackslash- \arraybackslash1.13 \arraybackslash- \arraybackslash- \arraybackslash-
\arraybackslash2.00 \arraybackslash2.83 \arraybackslash2.62 \arraybackslash1.53 \arraybackslash1.30 \arraybackslash2.48 \arraybackslash2.12 \arraybackslash1.00 \arraybackslash0.59 \arraybackslash1.88 \arraybackslash1.29
Table 2: Electron () and hole () effective masses along (x) and (y) and energy gap at the point () from DFT and TB fits. TB1 is the model which does not include overlap parameters, while the TB2 model contains all parameters. Theoretical reference values are taken from Ref. Qiao et al. (2014). Masses are in units of the bare electron mass and gap energies are in eV.

For transport properties, the most relevant bands are those in the vicinity of the Fermi energy. Therefore, if the model is intended to be used for such applications it is desirable that the fit reproduces well these bands, especially for k-points in the vicinity of the band-edge states. In order to further check the description of these bands, we have calculated the band gaps and effective masses along and for electrons and holes for both DFT and TB calculations. The results are reported in Table 2 together with reference values from DFT calculations from Ref. Qiao et al. (2014). We recall that the unoccupied bands from DFT were shifted in order to match the quasiparticle gap obtained in the reference GW calculation Tran et al. (2014). In particular, our DFT calculations predicts bulk-BP to be a semimetal, so the shift is specially important in this case. With the energy shift, bulk-BP has an indirect gap with the top of the valence band at the line and the bottom of the conduction band at the line, so we do not report the corresponding effective masses. Similarly, we can see that the TB model with overlap gives a very flat conduction band at the line in all cases, so we do not report the corresponding effective mass. As we can see, both TB models overestimate the gaps and effective masses in all cases. Nevertheless, we see that the full adjustment provides a good description of the whole band structure, with consistent effective masses for single-layer, bilayer and bulk BP and standard deviations eV in an energy range of about eV. The fit also captures the high effective mass anisotropy along and lines, but at a qualitative level.

iii.2 Adjustments for Low Energy Bands

We now discuss how the fit procedure can be improved in order to get a better description of the band structure near the Fermi energy. For that end, we focus our attention only on single-layer BP, since the model has already shown to be consistent for multilayer systems. As we are mainly interested in an accurate description of the energy gaps and effective masses near the point, we narrow down the energy range and k-point grid of the adjustment, which reduces the number of data points included in the summation in Eq. 5.

For this fit, we have considered k-points only along the path and included energy levels only within a eV window from the Fermi energy. Moreover, since the most relevant features of the valence (conduction) band are near the top (bottom) at the point and local maxima (minima) at nearby points, we have sampled these regions with a higher density of k-points. The TB fits for the models with (TB3) and without (TB4) overlap parameters are shown in Fig. 3. As we can see, the full model improves the fit of the valence and conduction bands, at the expense of the deeper occupied bands. The corresponding adjusted parameters are shown in Table 3. Note that some of these parameters are very different from those reported in Table 1, especially those involving hoppings and overlaps to orbitals. In fact, such parameters are more relevant to the description of deeper energy bands, not included in this fit. As such, the fit procedure becomes less sensitive to these parameters and greatly depends on their initial choice. This problem is particularly bad in the full model, where the calculation usually does not converge for energy ranges smaller than eV because the overlap parameters assume unphysical values during the minimization.

Figure 3: Comparison between DFT (black) and Tight-Binding band structures (red) for single-layer Black Phosphorus for bands in the vicinity of the Fermi Energy, which is set to zero. The top (bottom) panels correspond to a model without (with) overlap parameters.
Model \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash
TB3 \arraybackslash-3.58 \arraybackslash4.39 \arraybackslash3.78 \arraybackslash-1.43 \arraybackslash-7.34 \arraybackslash0.22 \arraybackslash0.54 \arraybackslash0.82 \arraybackslash0.68 \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash0.33
TB4 \arraybackslash-2.78 \arraybackslash3.93 \arraybackslash3.97 \arraybackslash-1.07 \arraybackslash-7.88 \arraybackslash0.42 \arraybackslash0.48 \arraybackslash0.88 \arraybackslash0.53 \arraybackslash0.32 \arraybackslash-0.05 \arraybackslash0.17 \arraybackslash0.06 \arraybackslash1.03 \arraybackslash5.18 \arraybackslash0.21 \arraybackslash2.27 \arraybackslash0.28
Table 3: Tight-Binding parameters for the adjustments shown in Fig. 3. Hopping and overlap amplitudes are calculated at nearest neighbor distance ( Å) and are in eV. Decay distances are in Å and is the standard deviation of the adjusment in eV.

Another feature observed in our calculations is that both models fit better the occupied states, as is evident in Figs. 2 and 3. In fact, the fit of these states can be further improved by completely leaving the unnocupied states out of the adjustment, as verified by our test calculations. On the other hand, even if the occupied states are completely left out of the fit, the adjustment of the conduction bands does not improve. In order to further understand this issue, we have performed projected density of states (PDOS) calculations within DFT. The results for the single-layer are shown in Fig. 4. Bilayer and bulk BP yield similar results, so we do not discuss them here. As we can see, the orbitals give an important contribution to the conduction band states. Since these orbitals are not included in our TB model, it is expected that the fit quality is not the same for these bands. However, it should also be pointed out that DFT calculations with the SZ basis, which does not include orbitals, give similar results when we use the same TB model to fit the resulting band structures, even tough the standard deviations are slightly smaller. This may be an indication that our simple model has a limitation in the description of the conduction bands, as also seen in similar models for other semiconductors Vogl et al. (1983). Another important feature of the PDOS is that the orbitals contribute mainly to the deeper bands, in agreement with our discussion above in terms of the TB parameters. There is some hybridization with the orbitals in the low energy range, but the latter still dominate. Moreover, our calculations indicate that the low-energy bands do not have a full character. Therefore, TB parametrizations which include only a single -like orbital per site are expected to not fully capture all features of these bands, as seen in the proposed models so far Rudenko and Katsnelson (2014, 2015). However, they are still very successful in fitting these bands in the vicinity of the point (and some higher energy states as well), thus correctly reproducing band gaps and effective masses.

Figure 4: Projected density of states (PDOS) of single-layer BP as given by DFT calculations with the DZP basis. The energy range is the same as the band structures in Fig. 2 and the Fermi energy is set to zero. We have also performed a scissors shift in order to reproduce the correct quasiparticle gap

Even with the limitations of our model, we may also get an improved description of the gap and effective masses near the point. To this end, we perform a new fit with k-points in a narrow vicinity of the point and a eV energy range. Only k-points with are included, where Å is the experimental lattice constant of bulk BP. The results are shown in Fig. 5 and the corresponding parameters are reported in Table 4 (TB5 and TB6 without and with overlap, respectively). Note that the standard deviation of the adjustment is now greatly reduced, ranging from to eV, due to the reduced number of data points. For the same reason, the adjusted parameters are much more sensitive to the fit conditions than those of the full adjustment reported in Table 1.

Figure 5: Comparison between DFT (black) and TB band structures (red) for single-layer black phosphorus for bands in a narrow vicinity of the point. The top (bottom) panels correspond to a model without (with) overlap parameters.
Model \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash \arraybackslash
TB5 \arraybackslash-0.78 \arraybackslash3.61 \arraybackslash3.46 \arraybackslash-1.32 \arraybackslash-16.4 \arraybackslash2.15 \arraybackslash0.36 \arraybackslash0.49 \arraybackslash0.01 \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash- \arraybackslash0.08
TB6 \arraybackslash-2.13 \arraybackslash2.64 \arraybackslash4.76 \arraybackslash-1.51 \arraybackslash-10.8 \arraybackslash0.84 \arraybackslash0.45 \arraybackslash0.49 \arraybackslash0.60 \arraybackslash0.19 \arraybackslash-0.07 \arraybackslash0.05 \arraybackslash0.14 \arraybackslash0.03 \arraybackslash0.27 \arraybackslash0.23 \arraybackslash0.18 \arraybackslash0.02
Table 4: Tight-Binding parameters for the adjustments shown in Fig. 5. Hopping and overlap amplitudes are calculated at nearest neighbor distance ( Å) and are in eV. Decay distances are in Å and is the standard deviation of the adjusment in eV.

The effective masses and the gap energies at the point are shown in Table 5. We can see that the narrow fit greatly improves the description of these quantities in comparison with the full adjustment shown in Table 2, especially the gap and the effective masses along . For the direction, the TB model without overlap parameters gives a conduction band minimum outside , so is not reported. On the other hand, the model with overlap gives a good description of both masses. We point out, however, that since the valence band is almost flat along this line, numerical values for the corresponding effective mass can be imprecise.

\arraybackslashDFT \arraybackslashTB5 \arraybackslashTB6
\arraybackslash0.13 \arraybackslash0.21 \arraybackslash0.16
\arraybackslash0.13 \arraybackslash0.22 \arraybackslash0.15
\arraybackslash3.42 \arraybackslash1.18 \arraybackslash2.68
\arraybackslash1.14 \arraybackslash- \arraybackslash1.07
\arraybackslash2.00 \arraybackslash2.08 \arraybackslash1.99
Table 5: Electron () and hole () effective masses along (x) and (y) and energy gap at the point () from DFT and TB fits in a narrow vicinity of the point. TB5 is the model which does not include overlap parameters, while TB6 contains all parameters. Masses are in units of the bare electron mass and gap energies are in eV.

Iv Conclusions

In summary, we have developed a Tight-Binding parametrization based on the Slater Koster model in order to describe the band structures of multilayer black phosphorus systems as obtained from first principles calculations. The model includes and orbitals and we model the behavior of the hopping and overlap amplitudes with distance as an exponential decay, which reduces the total number of parameters. We find that our model can fit the DFT band structures over an energy range of about eV with standard deviations of eV. It fits especially well the occupied bands, but we find it has limitations in the description of the conduction bands, even tough the inclusion of overlap parameters improves their fit. The adjusted hopping parameters are consistent throughout single-layer, bilayer and bulk BP, which is an indication the model is suitable for multilayer systems. On the other hand, overlap parameters are in general less consistent and more dependent on the conditions of the fit.

For the valence and conduction bands, we find that the full adjustment, in which all bands are included, does not give a sufficiently accurate description of the energy gaps and effective masses, but it captures the mass asymmetry at a qualitative level. For an improved description of these quantities, the fit procedure needs to be carried on bands in a narrow energy range and k-points in the vicinity of the gap, leading to standard deviations of eV.

The description of the conduction bands may be improved with the inclusion of fictitious orbitals which lead to band repulsion, as is usually done in TB models for semiconductors Vogl et al. (1983). The dependence of the hopping and overlap amplitudes with distance can also be modified to more complex expressions, such as the product of polynomials and exponentials. One can also go beyond the two center approximation and include more types of hoppings and overlaps. Finally, we point out that for the study of impurities or heterostructures, local modifications of hopping and on-site energies may also be required. For the least-squares fit, the use of methods that perform a search over all the parameter space, such as genetic algorithms, may lead to improved adjustments. All these modifications will lead to more complex models with an increased number of parameters and longer computation times, so we believe our model provides an important, yet simple starting point for tackling these situations.

This work was supported by the following Brazilian funding agencies: CNPq, CAPES, FAPERJ and INCT-Nanomateriais de Carbono. We also thank the computational support from DIMAT-INMETRO.


  1. preprint:


  1. A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,  and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. 81, 109 (2009).
  2. K. S. Novoselov, V. I. Falko, L. Colombo, P. R. Gellert, M. G. Schwab,  and K. Kim, “A roadmap for graphene,” Nature 490, 192 (2012).
  3. Taisuke Ohta, Aaron Bostwick, Thomas Seyller, Karsten Horn,  and Eli Rotenberg, “Controlling the electronic structure of bilayer graphene,” Science 313, 951 (2006).
  4. Chun Hung Lui, Zhiqiang Li, Kin Fai Mak, Emmanuele Cappelluti,  and Tony F. Heinz, “Observation of an electrically tunable band gap in trilayer graphene,” Nature Physics 7, 944 (2011).
  5. Fan Zhang, Bhagawan Sahu, Hongki Min,  and A. H. MacDonald, “Band structure of abc-stacked graphene trilayers,” Phys. Rev. B 82, 035409 (2010).
  6. F. Guinea, A. H. Castro Neto,  and N. M. R. Peres, “Electronic states and landau levels in graphene stacks,” Phys. Rev. B 73, 245426 (2006).
  7. Andrea C. Ferrari, Francesco Bonaccorso, Vladimir Fal’ko, Konstantin S. Novoselov, Stephan Roche, Peter Boggild, Stefano Borini, Frank H. L. Koppens, Vincenzo Palermo, Nicola Pugno, Jose A. Garrido, Roman Sordan, Alberto Bianco, Laura Ballerini, Maurizio Prato, Elefterios Lidorikis, Jani Kivioja, Claudio Marinelli, Tapani Ryhanen, Alberto Morpurgo, Jonathan N. Coleman, Valeria Nicolosi, Luigi Colombo, Albert Fert, Mar Garcia-Hernandez, Adrian Bachtold, Gregory F. Schneider, Francisco Guinea, Cees Dekker, Matteo Barbone, Zhipei Sun, Costas Galiotis, Alexander N. Grigorenko, Gerasimos Konstantatos, Andras Kis, Mikhail Katsnelson, Lieven Vandersypen, Annick Loiseau, Vittorio Morandi, Daniel Neumaier, Emanuele Treossi, Vittorio Pellegrini, Marco Polini, Alessandro Tredicucci, Gareth M. Williams, Byung Hee Hong, Jong-Hyun Ahn, Jong Min Kim, Herbert Zirath, Bart J. van Wees, Herre van der Zant, Luigi Occhipinti, Andrea Di Matteo, Ian A. Kinloch, Thomas Seyller, Etienne Quesnel, Xinliang Feng, Ken Teo, Nalin Rupesinghe, Pertti Hakonen, Simon R. T. Neil, Quentin Tannock, Tomas Lofwander,  and Jari Kinaret, “Science and technology roadmap for graphene, related two-dimensional crystals, and hybrid systems,” Nanoscale 7, 4598–4810 (2015).
  8. Sheneve Z. Butler, Shawna M. Hollen, Linyou Cao, Yi Cui, Jay A. Gupta, Humberto R. Gutiérrez, Tony F. Heinz, Seung Sae Hong, Jiaxing Huang, Ariel F. Ismach, Ezekiel Johnston-Halperin, Masaru Kuno, Vladimir V. Plashnitsa, Richard D. Robinson, Rodney S. Ruoff, Sayeef Salahuddin, Jie Shan, Li Shi, Michael G. Spencer, Mauricio Terrones, Wolfgang Windl,  and Joshua E. Goldberger, “Progress, challenges, and opportunities in two-dimensional materials beyond graphene,” ACS Nano 7, 2898–2926 (2013).
  9. A. K. Geim and I. V. Grigorieva, “Van der waals heterostructures,” Nature 499, 419–425 (2013).
  10. K. S. Novoselov, A. Mishchenko, A. Carvalho,  and A. H. Castro Neto, “2d materials and van der waals heterostructures,” Science 353 (2016).
  11. Likai Li, Yijun Yu, Guo Jun Ye, Qingqin Ge, Xuedong Ou, Hua Wu, Donglai Feng, Xian Hui Chen,  and Yuanbo Zhang, “Black phosphorus field-effect transistors,” Nature Nanotechnology 9, 372–377 (2014).
  12. Wanglin Lu, Haiyan Nan, Jinhua Hong, Yuming Chen, Chen Zhu, Zheng Liang, Xiangyang Ma, Zhenhua Ni, Chuanhong Jin,  and Ze Zhang, “Plasma-assisted fabrication of monolayer phosphorene and its raman characterization,” Nano Research 7, 853–859 (2014).
  13. Jingsi Qiao, Xianghua Kong, Zhi-Xin Hu, Feng Yang,  and Wei Ji, “High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus,” Nature Communications 5 (2014).
  14. Han Liu, Adam T. Neal, Zhen Zhu, Zhe Luo, Xianfan Xu, David Tománek,  and Peide D. Ye, “Phosphorene: An unexplored 2d semiconductor with a high hole mobility,” ACS Nano 8, 4033–4041 (2014).
  15. Fengnian Xia, Han Wang,  and Yichen Jia, “Rediscovering black phosphorus as an anisotropic layered material for optoelectronics and electronics,” Nature Communications 5 (2014).
  16. Xiaomu Wang, Aaron M. Jones, Kyle L. Seyler, Vy Tran, Yichen Jia, Huan Zhao, HanWang, Li Yang, Xiaodong Xu, ,  and Fengnian Xia, “Highly anisotropic and robust excitons in monolayer black phosphorus,” Nature Nanotechnology 10, 517–521 (2015).
  17. Likai Li, Jonghwan Kim, Chenhao Jin, Guo Jun Ye, Diana Y. Qiu, Felipe H. da Jornada, Zhiwen Shi, Long Chen, Zuocheng Zhang, Fangyuan Yang, Kenji Watanabe, Takashi Taniguchi, Wencai Ren, Steven G. Louie, Xian Hui Chen, Yuanbo Zhang,  and Feng Wang, “Direct observation of the layer-dependent electronic structure in phosphorene,” Nature Nanotechnology 12, 21 (2017).
  18. Vy Tran, Ryan Soklaski, Yufeng Liang,  and Li Yang, “Tunable band gap and anisotropic optical response in few-layer black phosphorus,” Phys. Rev. B 89, 235319 (2014).
  19. Mohammad Elahi, Kaveh Khaliji, Seyed Mohammad Tabatabaei, Mahdi Pourfath,  and Reza Asgari, “Modulation of electronic and mechanical properties of phosphorene through strain,” Phys. Rev. B 91, 115412 (2015).
  20. Jun Dai and Xiao Cheng Zeng, “Bilayer phosphorene: Effect of stacking order on bandgap and its potential applications in thin-film solar cells,” Journal of Phys. Chem. Lett. 5, 1289 (2014).
  21. Shengxue Yang Yan Li and Jingbo Li, “Modulation of the electronic properties of ultrathin black phosphorus by strain and electrical field,” Journal of Phys. Chem. C 118, 23970 (2014).
  22. A. N. Rudenko and M. I. Katsnelson, “Quasiparticle band structure and tight-binding model for single- and bilayer black phosphorus,” Physical Review B 89, 201408 (R) (2014).
  23. A. N. Rudenko and M. I. Katsnelson, “Toward a realistic description of multilayer black phosphorus: From gw approximation to large-scale tight-binding simulations,” Physical Review B 92, 085419 (2015).
  24. Carlos J. Páez, Kursti DeLello, Duy Le, Ana L. C. Pereira,  and Eduardo R. Mucciolo, “Disorder effect on the anisotropic resistivity of phosphorene determined by a tight-binding model,” Phys. Rev. B 94, 165419 (2016).
  25. P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136, B864 (1964).
  26. W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev. 140, A1133 (1965).
  27. J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón,  and D. Sánchez-Portal, “The siesta method for ab initio order-n materials simulation,” J. Phys.: Cond. Matt. 14, 2745 (2002).
  28. J. P. Perdew, K. Burke,  and M. Ernzerhof, “Generalized gradient approximation made simple,” Phys. Rev. Lett. 77, 3865 (1996).
  29. N. Troullier and Jose Luis Martins, “Efficient pseudopotentials for plane-wave calculations,” Phys. Rev. B 43, 1993 (1991).
  30. H. J. Monkhorst and J. D. Pack, “Special points for brillouin-zone integrations,” Phys. Rev. B 13, 5188 (1976).
  31. A. Brown and S. Rundqvist, “Refinement of the crystal structure of black phosphorus,” Acta Cryst. 19, 684 (1965).
  32. Mark S. Hybertsen and Steven G. Louie, “Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies,” Phys. Rev. B 34, 5390 (1986).
  33. J. C. Slater and G. F. Koster, “Simplified lcao method for the periodic potential problem,” Phys. Rev. 94, 1498 (1954).
  34. William H. Press, Saul A. Teukolsky, William T. Vetterling,  and Brian P. Flannery, Numerical Recipes in Fortran 90 (2Nd Ed.): The Art of Parallel Scientific Computing (Cambridge University Press, New York, NY, USA, 1996).
  35. W.A. Harrison, Elementary Electronic Structure, Elementary Electronic Structure (World Scientific, 1999).
  36. P. Vogl, Harold P. Hjalmarson,  and John D. Dow, “A semi-empirical tight-binding theory of the electronic structure of semiconductors,” J. Phys. Chem. 44, 365 (1983).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description