Atomic structural and electronic bandstructure calculations for borophene

Density of states (DOS) and electronic bandstructure diagrams with ε(k) versus k are found for particular allotropes of borophene with much improved accuracy by ab initio quantum calculations using hybrid functionals of several types. The particular types of hybrid functionals are delineated in detail. Varying levels of k-point discretization are utilized to evaluate accuracy. Structural relaxation has been carefully applied prior to electronic bandstructure simulations. Results indicate whether or not one has regions in k-space which display Dirac type non-gapped behavior or parabolic gapped behavior. This work is required in order to determine what types of electronic uses 2D single atomic layer borophene is appropriate for in modern nanoscopic devices.


Introduction
Moving beyond traditional 2D materials is now desirable at this point in time in order to provide alternatives for electronic devices which have various capabilities including switching capabilities. The latter capabilities entail transport for transistor performance. Here we propose using borophene because, as we will show in this letter, obtaining regions of the electronic bandstructure which act as valence and conduction bands, with or without bandgaps, may be obtainable in the foreseeable future, either computationally or experimentally, depending upon how carefully the atomic material is prepared, such that satisfactory devices are obtained for actual electronics, RF or digital.
Just to the left of carbon in the Periodic table is boron. Creation of borophene, consisting entirely of boron atoms in a 2D sheet form [1], makes one recall ideas of using boron nitride (BN) as a material for electronics. However, the intrinsic simplicity of dealing with a monoatomic material is tremendously attractive. One is reminded of the initial fascination with silicon and germanium in the 1950s and 1960s. Silicon grew into the workhorse of electronics for transistors for digital computers and lower frequency RF electronics. In particular, one can think of the use of silicon for affordable cell phones and digital cameras.
However, for high frequency electronic devices, it is well known that antimonide-based compound semiconductors can have very high mobility and offer important advantages for low-power, low-noise RF electronics [2]. Unfortunately, antimonide-based materials suffer several important disadvantages including In and Sb do not have commercially-viable domestic sources [3]; InSb-based electronics can be difficult to monolithically integrate with other electronics technologies (Si-based CMOS digital, GaN-based high-power RF, etc) because of the large lattice-mismatch between the materials families; there is continuing pressure to move to less toxic materials electronics (e.g. ROHS standard); finally, while III-V compound semiconductors have demonstrated f T and f Max cutoff frequencies in the low THz regime, even higher cutoff frequencies are needed for future commercial and DoD systems applications and it is not clear that III-V semiconductors can meet those demands [4]. New 2D materials with band gaps are predicted to have clear advantages for these

First principles calculations
Now look at some results we have obtained. Figure 1 shows atomic structures of the α−1 borophene and β-1 borophene allotropes. Figure 2 gives the density of states (DOS) for these two allotropes of borophene. Figure 3 provides the energy ε(k) versus k bandstructure through symmetry points in the Brillouin Zone for these two borophene allotropes. From the literature, it is known that first-principles calculations using the traditional generalized gradient approximation (GGA) as the exchange correlation functional usually predict no band gap for the borophene allotropes. We performed density functional theory (DFT) calculations using Perdew-Burke-Ernzerhof (PBE) GGA functional for both alpha-1 and beta-1 borophene nanostructures, and find no band gap for both structures, in agreement with previous publications. We show the calculated band structure and density of states (DOS) of alpha-1 borophene in an expanded view in figure 4.
Wu et al [6] demonstrated that when they use hybrid functional PBE0, which is supposed to be more accurate than GGA as the exchange-correlation functional, in their DFT calculations, several alpha-types of borophene allotropes showed band gaps of∼1 eV (figure 3), while other types of allotropes still have no gaps. Because borophene bandgap properties are dialed in through the allotrope crystallographic preparation, one can in principle be able to choose in the material class both the desired band structure and its type. That is, the type of band structure can be selected to be either Dirac type with its relativistic analogy to low mass particles (and no bandgap), or parabolic type with its more massive characteristic of semiconductor action (and a bandgap). In either case the desirable properties should be maintainable over actual 2D monoatomic sheets, and not necessitate the tactic of creating nanoribbons to craft a bandgap as done for graphene, for example. Figure 2. Calculated density of states (DOS) using density functional theory (DFT) for two borophene allotropes α-1 and β-1. The peaks above and below the Fermi level E F at 0 eV suggest conduction and valence band action as in a semiconductor. Non-zero DOS at ε F indicates an imperfect bandgap. Some allotropes have zero here. Used the PBE functional.  One of the most important questions is to identify which allotropes of borophene are mostly likely to present finite reproducable bandgaps, a prerequisite for a useful semiconductor material for nanoscopic and atomic scale electronic devices. In particular, transistor devices acceptable for RF electronics or digital electronics. To better evaluate which allotropes might have potential, various hybrid functionals were utilized to determine the effects on the electronic band structure.

Hybrid functional theory summary
To simulate the 2D planar nanostructure, we set up a large-size supercell, with one single atomic layer of boron atoms in the middle. We use period boundary conditions along all the directions in the simulation, thus the borophene layers in adjacent cells are separated by vacuum of ∼30 angstroms, and the interactions between them are essentially negligible, as shown in figure 5.
Hybrid functional [36] is a class of approximations to the exchange-correlation energy functional in DFT that incorporate a portion of exact exchange from Hartree-Fock theory with the rest of the exchangecorrelation energy from other sources (ab initio or empirical), which provides a simple scheme for improving the calculation of many materials properties.
Hybrid functionals, such as B3LYP or PBE0, containing some amount of exact exchange, as in Hartree-Fock theory (spin-restricted case for simplicity): There are three major types of hybrid functional energy E , xc hyb known as PBE0 [37], HSE [38], and B3LYP [39], Here the first three terms of equation (1) are identified as, T s the kinetic energy part of the Hamiltonian, V(r) the potential energy part of the Hamiltonian at spatial location r with an electron density per unit volume of n(r), and E H is the Hartree part of the Hamiltonian. Fourth term of (1) is E , x HF the exchange Hartree-Fock part of the Hamiltonian. Generally, the E xc HF is the Hartree-Fock exchange and correlation energy part of the Hamiltonian (the more complete form implicitly includes the Pauli exclusion principle with spin indices), but in our above expression, only the exchange part is included and the correlation part is taken from DFT calculations. In the spirit of including DFT calculations in the first equation, exchange DFT Hamiltonian E x DFT is mixed in as a proportion (1−α x ), with the Hartree-Fock part Hamiltonian E x HF weighted proportionally as α x . The hybrid DFT Hamiltonian part including the last three terms of (1) is labeled as E , xc hyb and is equal to one of the three possibilities shown in (2). That is, , .  Finally, the Hartree-Fock Hamiltonian is given by the integral Returning to the suite of equations (2), in (2b), SR and LR in the superscripts of energy stand for, respectively, short range and long range Coulomb interactions in the exchange part, where ω is an adjustable parameter. In (2c), the superscripts LDA and GGA in the energies denote, respectively, local density approximation and generalized gradient approximation (discussed above). For LDA contributions, Equations (5) gives either the total exchange energy or the total correlation energy, depending upon the subscript, based upon weighting the contributions per electron particle for a homogeneous electron gas (HEG) model. Here ( ) n r is the electron density. Particularization to spin density orientation upgrades (5a) to (5b). For the HEG model, for exchange, the energy density for a particle is known analytically, and results in a formula weighted by charge density when non-uniformity is the case for total energy.
For correlation energy, again analytical expressions exist for HEG. In the low and high density limits (extremely low and high correlation), the expressions are of the forms [40], or for a single formula [41], with a and b arising from constraint that high density limit be matched, In equations (7), r s is the Wigner-Seitz parameter, a dimensionless quantity. It represents a sphere which encompasses exactly one electron charge, and is normalized to the Bohr radius. In a 3D system, it is given by To improve upon the LDA which works only with electron density n(r), the GGA method includes gradient information on the electron density too, with the form shown handling spin too.
n n n n n n n r r , , , Here this GGA form incorporates the spin information in a formalistic way, in parallel with the LDA method upgraded to the LSDA in equations (5). Simplest form of the GGA is the Chachiyo exchange-correlation functional (CECF), given by [42,43] The CECF is comparable in accuracy to other well known functionals such as B3LYP, PBE and TPSS. Here in equations (10a) and (10b), respectively, x, t and h represent two scaled gradients of the electron density, and h = 0.06672632 Hartree.
Several different methods to construct the hybrid functional have been developed recently. The PBE0 functional mixes the PBE exchange energy and Hartree-Fock exchange energy in a set 3:1 ratio, along with the full PBE correlation energy. The HSE (Heyd-Scuseria-Ernzerhof) exchange-correlation functional uses an error-function-screened Coulomb potential to calculate the exchange portion of the energy in order to improve computational efficiency. The popular B3LYP (Becke, 3-parameter, Lee-Yang-Parr) exchange-correlation functional combines the Becke 88 exchange functional, the correlation functional of Lee, Yang and Parr for B3LYP, and GGA and local density approximations (LDA) to the correlation functiona.

Supercomputer results using hybrid functionals
Wu et al [6] used only PBE0 as the hybrid functional in their research on borophene nanostructures. It is important to know how different types of hybrid functional might affect the results, so this article addresses all three types of hybrid functionals, PBE0, HSE, and B3LYP to examine borophene nanostructures, with the technical objectives to evaluate and identify the accuracy and efficiency of recently-developed hybrid functionals and determine effects on the borophene electronic band structure.
We use quantum espresso, an integrated suite of open-source computer codes for electronic-structure calculations and materials modeling at the nanoscale. Quantum espresso is based on density-functional theory, plane waves basis set, and pseudopotentials, with many features to examine the ground-state calculations, including structural optimization, molecular dynamics, potential energy surfaces, electrochemistry and special boundary conditions, response properties, spectroscopic properties, quantum transport, and platforms. Quantum espresso implements all different types of popular hybrid functional: PBE0, B3LYP, HSE hybrid GGA, and works with both the normal-conserving pseudopotential (NCPP) and ultrasoft-pseudopotential (USPP) with some limitations. Hybrid functional DFT calculation is the newly developed and evolving feature in quantum espresso, which we exploit here.
In typical DFT calculations, in order to obtain the band structure of a given material, one usually first performs a Self-Consistent Field (SCF) calculation to approach convergence, and then a non-SCF calculation using the SCF charge density but on a much finer k-point mesh to obtain the band energies on enough k-points in the Brillouin zone. Unfortunately, such non-SCF DFT calculations using hybrid functional are still not supported in the latest version of quantum espresso. In order to obtain the band structure from hybrid functional DFT calculations, it is necessary to perform SCF calculations on much finer k point mesh, which are usually very computationally demanding. Since our supercell has relatively large size (∼30 angstrom) along the direction perpendicular to the borophene plane, we only need one k point along that direction. The k point sets we tested include a list of five items: Hybrid functional calculations using 16x16x1 k-point mesh turns out to be too computationally demanding, which the tests are not able to approach SCF convergence after running large-scale parallel calculations using hundreds of cores for 120 h on the supercomputer Thunder at AFRL DSRC. For smaller k-point mesh such as 12 × 12 × 1, we are able to approach SCF convergence using hundreds of cores with a reasonable amount of time, namely, about 70 h. The tests and computational work have been performed on supercomputers Thunder and Mustang at AFRL DSRC, Excalibur at ARL DSRC, and Polar at NRL-DC. See figures 6 (a) and (b) for SCF convergence time versus grid case number on list above. Figure 6(b) highlights the nonlinear nature of the points by plotting a rough guide through the actual points. When the horizontal scale is converted to logarithm of the k-point area grid in normalized momentum grid space (just the calculated value from the products in each listed case above), another interesting plot is found, figure 6(c), which demonstrates since the curve looks very similar to that in figure 6(b), that the grid point case number and the logarithm of the area in k-point space are proportional to each other. Figure 7 shows typical ε(k) versus n, the grid areal integer number, from runs on the supercomputer Thunder at AFRL DSRC, for three energy bands, shown in (a), (b) and (c). Here three symmetry points are shown for each band: the Gamma (or Γ), Y and S points. Each symmetry point shows a decline in energy value until the last n with a slight or noticeable uptick in magnitude. The last increase in value may indicate an oscillatory behavior to a converged value. Since we report our results in eV, that is what is used as the vertical scale.
For alpha-1 borophene nanostructure shown in the left side of figure 2, we use 12 × 12 × 1 k point mesh in the hybrid functional DFT calculations, which produces totally 49 k-points in the irreducible part of the Brillouin zone. In figure 8 we show the calculated band structures of alpha-1 borophene using PBE0, HSE and B3LYP as the hybrid functional, respectively. B3LYP predicts a slightly different band structure from PBE0 and HSE. The band structures from PBE0 and HSE are essentially identical, although the exchange-correlation functionals used are quite different, as shown in the quantum espresso output files,    Table 1. Calculated band negative energies (eV) of α-1 borophene nanostructure using PBE0.  Table 3. Calculated band negative energies (eV) of α-1 borophene nanostructure using B3LYP. might account for the discrepancies here. First, different pseudopotentials have been used. Second, and the most likely reason, is due to the structural relaxation. We performed relaxations using the structure from Wu et al [6] as the initial structure, and we notice some changes in the supercell dimensions after the structural relaxation.
Borophene is known to have many different metastable allotropes where the energies of different allotropes are pretty close. We also notice that for all our hybrid functional calculations for borophene alpha-1 structure, there is only one band crossing near the Fermi level, and along that band there is only one point (Y point in the Brillouin zone) that has energy below the Fermi level. Slight change in the a direction on the borophene plane (corresponding to the Y direction in the Brillouin zone) might change the band energy around the Y point. To make further detailed comparisons of the calculated band structure from different hybrid functionals, we shift the Fermi level of all the calculated band structures to zero and plot the calculated band structures together, as shown in figure 9. It can be seen that the band structures of PBE0 and HSE are essentially identical for both the valence band and the conduction band. The calculated valence bands from B3LYP are in close agreements with PBE0 and HSE, but display some noticeable discrepancies in the conduction band regime, especially close to the Y and S points in the Brillouin zone. Detailed tabulated values for respectively PBE0, HSE, and B3LYP functionally determined energies at three of the symmetry points S, Γ and Y are shown in tables 1, 2, and 3, for alpha-1 borophene. We show the calculated band structures of beta-1 borophene nanostructure using PBE0, HSE and B3LYP as the hybrid functional in figure 10. There are totally 74 k points in the irreducible part of the Brillouin zone due to the symmetry f of the crystal structure. Similar to the alpha-1 borophene results, the band structures calculated from PBE0 and HSE are essentially identical. To make further detailed comparisons of the calculated band structure from different hybrid functionals, we shift the Fermi level of all the calculated band structures to zero and plot the calculated band structures together, as shown in figure 11. The calculated band structure of beta-1 borophene using B3LYP agree well with HSE and PBE0 in the valence band regime, but shows a noticeable difference only in the highest conductance band energy, near the Y point in the Brillouin zone. All the hybrid functional calculations show that there are more than one band crossing the Fermi level, indicating the metallic nature of beta-1 borophene, in agreement with previous research. Again, detailed tabulated values for respectively PBE0, HSE, and B3LYP functionally determined energies at three of the symmetry points S, Γ and Y are shown in tables 4, 5, and 6, for beta-1 borophene.

Conclusions and summary
In summary, we tested and evaluated various forms of hybrid functionals in quantum espresso for borophene nanostructures. For both alpha-1 and beta-1 borophene nanostructures, the calculated band structures from PBE0 and HSE are essentially identical, and the band structure calculated using B3LYP are in close agreements with HSE and PBE0 results in the valence band regime, but show some modest differences in the conductance band regime near one symmetry point. All the calculated band structures of alpha-1 borophene nanostructure show no band gap, in contrast to Wu et al [6] who predicted a small band gap using PBE0, which is probably caused by the structural relaxation and the availability of many allotropes of borophene nanostructures. Hybrid functional calculations in quantum espresso take much longer computational time (∼50-100 times) than traditional local density approximation and general gradient approximation as the exchange-correlation functional, and thus it requires significantly more supercomputer resources if one wants to perform hybrid functional DFT calculations using quantum espresso.
Our main conclusion is that it is not a foregone conclusion that electronic structure bandgap can be assumed to occur in the 2D borophene allotropes. Use of several hybrid functionals all agree on the fact that obtaining a finite reproducible bandgap in grown material sheets of single atomic borophene will be a challenge. To make borophene single atomic sheets competitive with what silicon did for more macroscopic bulk-like transistor devices, will require an integrated research program of materials growth, first principles calculations, new materials growth, and careful experimental evaluation of the grown borophene sheets.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).