1. Introduction and Motivations
In recent years, several studies have been carried out on the structural, electronic and transport properties of beryllium selenide (BeSe). BeSe is a member of the alkaline earth selenides. This wide band-gap semiconductor has attracted great interest for electrical and optoelectronic devices and as a promising base material for blue-green laser diodes and light emitting diodes  . This compound crystallizes in fourfold-coordinated, cubic zinc-blende (B3) structure at ambient temperature. A unique characteristic of this compound is the ratio of the extremely small cation (Be) to a much larger anion (Se). This uniqueness accounts for the high degree of covalent bonding and is similar to the case of boron based group III-V compounds  . Furthermore, the hardness, bonding energy and stability of this compound make it a potentially good material for various technological applications  .
Few experimental studies have been performed on this compound; experimental difficulties include its toxic nature and its instability in air. This compound is difficult to handle experimentally as a single crystal or an epitaxial layer. Yim et al.  prepared crystalline powder of BeSe by passing H2Se vapor over Be metal at 1100˚C, repeatedly, for a total of 12 hours, until an X-ray powder pattern showed sharp K doublets. These authors performed measurements at room temperature; they estimated the band gap to be within the range of 4 - 4.5 eV and clearly stated the need for further investigations. They utilized the optical absorption measurements on cold pressed samples of BeSe powders mixed in KBr. These authors did not specify whether the band gap was direct or indirect because of shallow absorption edges and the lack of a high absorption coefficient. They stated in their work that further studies have to be done on this compound to understand the band structure. In 1999, Wilmers and his group  employed spectroscopic ellipsometry in the UV/VUV region. They investigated the optical properties of various samples of BexZn1-xSe by varying the stoichiometry of beryllium and selenium to obtain BeSe at room temperature. The BexZn1−xSe layers were grown on GaAs in a molecular-beam epitaxy chamber. The thickness of the layers ranged from 200 to 800 nm. This group analyzed the structure of the spectra in the pseudodielectric function to obtain a direct band gap of 5.5 eV. Depending on the thickness of the actual samples of zb-BeSe, the band gap of 5.5 eV could be a slight overestimate for the band gap of bulk zb-BeSe; quantum confinement is known to lead to larger band gaps for thinner films. The above two experimental reports place the band gap of zb-BeSe in the range of 4 eV to 5.5 eV. As shown below, even the lower limit of 4 eV is underestimated by previous, ab-initio theoretical DFT calculations.
Several theoretical studies  -  of electronic and related properties of BeSe have been reported. Recently, Yu et al.  studied the structural and electronic properties of BeSe using the plane-wave pseudopotential method. Their calculations produced an indirect band gap of 2.73 eV. Guo and his colleagues  employed both GGA and LDA potentials to investigate the electronic, optical, and structural properties of zb-BeSe. These authors utilized the plane-wave pseudopotential method in both calculations. They obtain an indirect band gap of 2.787 eV and 2.402 eV, respectively, with GGA and LDA potentials, from this work. The full-potential linearized augmented plane wave plus local orbitals (FP- LAPW + lo) calculations of Allay-e-Abbas et al.  , for the zb-BeSe, led to an indirect band gap of 2.43 eV. The 2011 LDA and GGA work of Al-Douri et al.  respectively, obtained band gap values of 2.397 eV and 2.682 eV; they employed the full-potential linearized augmented plane wave (FP-LAPW) method.
The calculated band gap value of 2.4179 eV was reported in the work of Rached et al.  . They utilized the full-potential linear muffin-tin orbitals (FP- LMTO). The plane-wave pseudopotential method was employed by Srivastava et al.  . They obtained a band gap of 2.43 eV. The LDA study of Khenata et al.  produced a band gap of 2.475 eV. Utilizing the full-potential linearized augmented plane wave (FP-LAPW) method, Hassan and Akbarzadeh  investigated the ground state properties and structural phase transition of zb-BeSe. They employed the LDA and GGA functionals to obtain two different band gaps of 2.33 eV and 2.66 eV, respectively. The full potential linearized augmented- plane wave calculations by Berghout et al.  , with an LDA potential, resulted in a band gap of 2.41 eV. Another LDA study done by the same authors  utilized the plane-wave pseudopotential method to produce a slightly higher band gap value of 2.43 eV. Furthermore, Heciri and his group  performed first- principle calculations to study the electronic structure of BeSe. In this work, the full-potential linearized augmented plane wave plus local orbital (APW + lo) method was used. They reported an indirect band gap of 2.23 eV and 2.51 eV for PW-LDA and PBE-GGA, respectively. The scalar relativistic calculation reported by Okoye  is employed with full-potential linearized augmented plane wave (FP-LAPW) approach. He obtained a band gap of 2.63 eV using the Perdew, Burke and Ernzenhof (PBE) GGA. Table 1 below shows the GGA results are mostly higher than those obtained with LDA potentials. The band gaps obtained from other formalisms such as Hartree Fork, Green function and screened coulomb approximation (GWA) are shown in Table 1.
In addition to the previous LDA and GGA methods used by these authors  , they further employed the screened exchange LDA (sX-LDA) to produce a higher band gap value of 3.455 eV. Also, Alay-e-Abbas et al.  used the modified Becke and Johnson (mBJ) LDA to increase the band gap to 3.53 eV. Another pseudopotential calculation of Yadav and his group  found the band gap to be 3.59 eV. They employed both GGA and GW approximation in their calculations. The Engel-Vosko GGA calculations of Al-Douri et al.  and El Haj Hassan et al.  led to a gap of 3.655 eV and 3.61 eV, respectively. Both authors utilized the full-potential linear augmented plane wave (FP-LAPW) method. In 2000, Fleszar and Hanke  studied the electronic excitation in BeSe by employing the ab-initio GW approximation to produce a gap of 3.66 eV. For the same compound, the early theoretical work performed by Stukel  used both a nonrelativistic self-consistent orthogonalized plane wave (SCOPW) method and the slater's free-electron-exchange approximation to determine the energy band structure. These calculations produced an indirect band gap of 3.61 eV. Seven years after, Sarkar and Chatterjee  obtained a gap of 4.37 eV from Г to K. In their work, they applied the APW method in conjunction with an LCAO interpolation scheme. Gonzalez-Diaz et al.  calculated a band gap of 2.39 eV using the first principle pseudopotential method.
Table 1. Previous, calculated, indirect band gaps of zb-BeSe and room temperature experi mental values.
aReference  bReference  cReference  dReference  eReference  fReference  gReference  hReference  iReference  jReference  kReference  lReference  mReference  nReference  oReference  pReference  qReference  rReference  .
From the theoretical results obtained with ab-initio LDA and GGA potentials, it is evident that the calculated band gap values were underestimated by an average of 1.5 eV or more as compared to the measured ones. While results produced by calculations using DFT-derived potentials, which are not entirely DFT potentials, are closer to the experimental ones, they still underestimate the latter by approximately 0.5 eV or more. The discrepancies between theoretical results and measured band gap values of zb-BeSe are a key motivation for our work. From the theoretical results obtained with ab-initio LDA and GGA potentials, it is evident that the calculated band gap values were underestimated by an average of 1.5 eV or more as compared to the measured ones. While results produced by calculations using DFT-derived potentials, which are not entirely DFT potentials, are closer to the experimental ones, they still underestimate the latter by approximately 0.5 eV or more. The discrepancies between theoretical results and measured band gap values of zb-BeSe are a key motivation for our work.
In this work, the computational method used has been described in detail in previous publications by our group  -  . We performed self-consistent calculations using the Ceperley and Alder local density approximation (LDA) potential  as parameterized by Vosko, Wilk, and Nusair  . The computational package used in this calculation is from the Ames laboratory of the US Department of Energy (DOE), Ames, Iowa   . We implemented the linear combination of atomic orbitals (LCAO), using Gaussian functions in the radial parts of the orbitals. Our calculations are nonrelativistic and were first performed at an experimental lattice constant of 5.152 Å (room temperature)  . In contrast to other calculations, we implemented the Bagayoko, Zhao, and Williams (BZW) method    , as enhanced by Ekuma and Franklin (BZW- EF)  . The difference between the enhanced version of our method (BZW- EF) and the previous method (BZW) is the pattern in which we increase the basis set. That is the addition of orbitals. In BZW, the orbitals representing unoccupied energies are added to the basis set in the order of increasing, excited energies of the atomic or ionic species in the solid. The BZW-EF method add s p, d, and f orbitals, for a given principal quantum number at a given site, if applicable, before adding the corresponding, spherically symmetric s orbital for that principal quantum number. This feature of the BZW-EF method rests on the realization of the primacy of the polarization of p, d, and f orbitals over the spherical symmetry of s orbitals, for the valence electrons of the material under study.
In line with the rules of the BZW-EF method, we start our calculations with a small basis set no smaller than the minimum basis (MB), i.e., the one that is just large enough to account for all the occupied energies in the atomic or ionic species of the material under study. We used the orbitals obtained from the self-consistent calculations for the atomic species of Be2+ and Se2− to construct the basis set for the solid calculations. After the first calculation, we performed a second self-consistent calculation, utilizing a basis set that consists of the one for Calculation I plus an appropriate, unoccupied orbital representing an excited energy level in the ionic species of the system. We compare the occupied energies from Calculations I and II, numerically and graphically. In general, they are found to be different, with some occupied energies from Calculation II being lower than corresponding ones from Calculation I. This lowering means that the initial basis set is not complete in size, angular symmetry or radial function for the description of the ground state of the material. We then perform a third calculation with the basis set of Calculation II plus an appropriate orbital chosen as described above. Again, we compare the occupied energies of Calculations II and III. This process continues until three consecutive calculations, i.e., N, (N + 1) and (N + 2), produce the same occupied energies. This fact is the criterion for the attainment of the absolute minima of the occupied energies, i.e., the ground state. These three calculations, upon the attainment of self-consistency, lead to the same charge density and Hamiltonian content, even though the Hamiltonian matrices have different dimensions. The first of the three calculations, N, is therefore the one providing the DFT description of the material. Its results have the full, physical content of DFT. Basis sets of Calculations (N + 1), (N + 2), and higher produce the occupied energies obtained from Calculation N, whose basis set is known as the optimal basis set. Unoccupied energies from Calculation (N + 1), (N + 2), and higher, that are different from their corresponding values from Calculation N do not belong to the spectrum of the Hamiltonian that is a unique functional of the ground state charge density, a density that did not change from its value obtained with the optimal basis set.
The above explanation for the selection of Calculation N is based on the first DFT theorem, as first provided by Bagayoko  . An equally valid selection of Calculation N is based on the Rayleigh theorem for eigenvalues. According to the Rayleigh theorem  , successive calculations with larger basis sets that contain the optimal one generate increasing numbers of eigenvalues, by virtue of the fundamental theorem of algebra. They do not change the occupied energies. This theorem explains the lowering of some unoccupied energies by Calculations (N + 1), (N + 2) and higher, after attainment of the absolute minima of the occupied energies. This lowering is a mathematical artifact, i.e., the non-trivial basis set and variational effect   totally avoided by the BZW and BZW- EF method as explained below.
The Rayleigh theorem  states that when the same eigenvalue equation is solved with two basis sets containing n and (n + 1) basis functions, respectively, with the smaller basis set totally included in the larger one, then the ordered eigenvalues (from the lowest to the highest) obtained with (n + 1) functions are lower than or equal to their corresponding values obtained with n functions. In the implementation of the BZW-EF method, we avoid the above basis set and variational effect by selecting the outputs from the calculation with the optimal basis set, i.e., the first one to produce the absolute minima of the occupied energies. Larger basis sets, that contain the optimal basis set, produce the charge density and Hamiltonian obtained with the optimal basis set. The changing (i.e. lowering) of an unoccupied energy by these larger basis sets is a mathematical artifact stemming from the Rayleigh    .
Details needed for the replication of our work follow. Beryllium selenide exhibits a cubic lattice in the space of F4-3m. The locations of the ions of Be and Se are at (0, 0, 0) and (1/4, 1/4, 1/4), respectively. We used a room temperature experimental lattice constant of 5.152 Å for the first part of our work. We performed self-consistent calculations for Be2+ and Se2− to obtain the input orbitals required in generating the orbitals used in the LCGO formalism for solid state calculations. A set of even-tempered Gaussian functions were employed in constructing the atomic orbitals of the ionic species. We used 16 even-tempered Gaussian functions with minimum and maximum exponents of 0.24 and 0.9 × 105, respectively, to describe the s and p orbitals of Be2+. The s and p orbitals of Se2− were constructed using 22 even-tempered Gaussian functions, with minimum and maximum exponents of 0.135 and 0.24 × 106, respectively. The convergence for a given self-consistent calculation was attained after 60 iterations, when the potential did not change by more than 10−5 between the last two consecutive iterations. The computational error made in accounting for the valence electrons was 0.00579210 for the 28 electrons or 2.0686 × 10−4 per electron.
Table 2 below contains orbitals employed in the successive solid state calculations. This table also provides the total numbers of orbitals for the description of the valence states, along with band gaps from Г to the minima of the conduction band at high symmetry points and elsewhere.
Our calculated LDA BZW-EF band structure is in Figure 1. This band structure is obtained from the basis set of Calculation II and III. It is observed that the bands from Calculation II cannot be distinguished from those from Calculation III up to 10 eV. In other words, there is a perfect superimposition of not only the occupied ones, but also for the conduction bands up to 10 eV. The calculated indirect, band gap from Γ to a conduction band minimum between Γ and X, is 5.46 eV.
Table 2. Successive, self-consistent calculations of the BZW-EF method for zb-BeSe (Calculations I-IV) and additional, illustrative Calculations (II* & III*). The basis set of Calculation II led to the absolute minima of the occupied energies; it is the optimal basis set. Column 1: Calculation number, Column 2: Valence orbitals for Be2+, Column 3: Valence orbitals for Se2−. Column 4: Total number of valence functions, Columns 5-8: the band gaps (Eg, in eV) from Г to L, Г to Г, Г to X, and from Г to X. The optical band gap is the smallest one, in Column VII, from Г to a conduction band minimum between Г and X, as obtained with the optimal basis set.
Figure 1. Electronic band structure of zb-BeSe as obtained from Calculations II (_) and III (- -), at an experimental lattice constant of 5.152 Å, using our BZW-EF method. The Fermi energy has been set to zero and its position is denoted by the horizontal, dotted lines. The calculated band gap, as obtained with the optimal basis set of Calculation II, is from Γ to a point between Γ and X. This band gap is 5.46 eV.
Our results for the total density of states (DOS) and partial densities of states (pDOS) of beryllium selenide are shown in Figure 2 and Figure 3, respectively. The total valence bandwidth is 13.59 eV. This result is close to those of Okoye (13.8 eV)  and Gonzalez-Diaz (14.16 eV)  . The valence bands are in two groups. The widths of the upper and lower groups are 5.25 eV and 1.81 eV, respectively. Our calculated width for the upper group of valence bands is close to the finding of Rached et al.  of 5.42 eV. The DOS for the valence states has a broad peak between −1 and −2.8 eV, a clear shoulder between −3 and −3.6 eV, and two sharp peaks at −4.4 and 12.0 eV. For the conduction band DOS, we found a mild shoulder at +6.4 eV and a sharp peak at +7 eV, atop a broad one from +6.4 eV to 10 eV.
As shown in Figure 3, the partial densities of states (pDOS) describe the contribution of various valence s, p, d, and states, if applicable, to the band structure. The upper valence band is mainly composed of the p states from Se with slight contributions of p and s states on Be. The lowest group of valence bands consists only of the s state of Se. On the other hand, the minimum of the conduction band, located between Γ and X, is mainly dominated by p and s states of Be, with a minor contribution from p on Se. According to Rached et al.  , the upper valence bands are dominated by p states of Se, while the lower valence bands are mainly from s states on Se. Our results qualitatively confirm this picture. The conduction band is due to the Be 2p. There are no reported experimental data for the DOS and pDOS that are known to us.
Figure 2. Calculated, total density of states (DOS) for zb-BeSe, as obtained from the bands in Figure 1. The vertical, dotted line denotes the Fermi level. The value of the band gap obtained is clearly shown in the insert.
Figure 3. Calculated, partial densities of states (pDOS) for zb-BeSe, as obtained from the bands in Figure 1.
In Table 3, we list the energies of the valence and low laying conduction bands at high symmetry points in the Brillouin zone. These energies are obtained from the self-consistent Calculation II, using the room temperature experimental lattice constant of 5.152 Å. The purpose of listing these energies is to enable detailed comparisons of our results with future experimental measurements, utilizing several techniques ranging from X-ray and UV spectroscopies to optical absorption.
Effective masses of electron and hole are important factors in determining the transport property of a material. In order to determine the accuracy of the shape and curvature of a calculated band, the results of the effective masses from the
Table 3. Calculated, electronic energies (eV) of zb-BeSe at high symmetry points in the Brillouin zone. These energy values are obtained from the optimal basis set of Calculation II at a room temperature experimental lattice constant of 5.152 Å. The Fermi energy is set to zero in the table. The minimum points of the lowest conduction band, at the high symmetry points, are shown with bold values.
theoretical work should be in agreement with those of the measured effective masses. We calculated the electron effective mass at X in the vicinity of the minimum of the conduction band as shown in Figure 1. From X to Г (longitudinal), X to U (transverse) and X to W (transverse) directions, the results for the electron effective mass in units of m0 (free electron mass) are 1.217 m0, 0.303 m0 and 0.302 m0, respectively. Furthermore, we calculated the hole effective masses at the top of the valence band. For heavy hole 1 along these directions, (Г-L)111, (Г-X)100, and (Г-K)110 the effective masses are 1.309 m0, 0.572 m0 and 0.891 m0, respectively. For heavy hole 2, along the same direction, (Г-L)111, (Г-X)100, and (Г-K)110, the effective masses are 1.309 m0, 0.572 m0, and 0.637 m0, respectively. Effective masses of the light hole along (Г-L)111, (Г-X)100, and (Г-K)110, directions are 0.178 m0, 0.285 m0, and 0.235 m0, respectively. According to the calculation done by Stukel  , the heavy hole effective masses are 1.3 m0 for Г-L(111) and 0.6 m0 for Г-X(100) while those of the light hole are 0.2 m0 for Г-L(111) and 0.3 m0 for Г-X(100), also the electron effective mass is 1.2 m0 for Г-X(100). There are no known experimental data available for the effective masses of zb-BeSe to compare our results. We expect future experiments to verify our values. Our results are in agreement with the theoretical works of Stukel  .
The calculated total energy versus the lattice constant is shown in Figure 4. We utilized this curve to obtain the calculated bulk modulus which is a measure of the degree of hardness of the material. The lattice constant corresponding to the minimum value of the total energy is known as the equilibrium lattice constant. The value of our predicted, equilibrium lattice constant is 5.044 Å. The calculated bulk modulus is 92.3 GPa. This result is in excellent agreement with the only experimental value known to us of 92.2 ± 1.8 GPa  . Other calculated          LDA results produced values that are about 1.0 GPa lower than ours. The hardness of zb-BeSe is related to its low ionic nature, unlike some other semiconductors in group II-VI with high ionic bonding.
Figure 4. Plot of the total energy (eV) of BeSe versus the lattice constant (Å). The minimum value of the total energy on the curve is located at the equilibrium lattice constant of 5.044 Å.
A review of the content of Table 1 leads to an obvious question. Indeed, the previous calculations utilizing ab-initio LDA or GGA potentials, unlike ours, have uniformly underestimated the band gap of zb-BeSe. Even the calculation utilizing DFT derived ad hoc potentials mostly underestimated the band gap. The key reasons for these underestimations are thoroughly explained by Bagayoko  . The first of these reasons is that the results of these calculations cannot be expected to possess the full, physical content of DFT, due to the fact that they utilized a single basis set; such calculations lead to a stationary solution among potentially infinite others. The minimization resulting from self-consistent iterations cannot correct for any major deficiency of the selected basis set in terms of size (i.e., number of functions), angular symmetry, and radial characteristics. The successive calculations of the BZW-EF method, with the increase of the size, angular features, and radial components of the basis set, verifiably lead to the absolute minima of the occupied energies, i.e., the ground state, as required by the second Hohenberg-Kohn theorem.
A second reason for the underestimation, also related to the use of a single basis set, stems from the fact that such a basis set, deliberately selected to be adequately large in size, is generally over-complete for the description the ground state. Consequently, unphysical, unoccupied energies, lower than their values obtained with the smallest basis set leading to the ground state, are generally produced by single basis set calculations. The unphysical nature of these lowered, unoccupied energies was explained in the section above on our method, in terms of both the first Hohenberg-Kohn theorem and the Rayleigh theorem for eigenvalues. The explanation for the agreement between our results and experimental ones rests on the fact that the BZW-EF calculations strictly adhere to the conditions inherent to the validity of DFT, including the verified, attainment of the ground state and the avoidance of over-complete basis sets by virtue of the first DFT theorem or the Rayleigh theorem for eigenvalues.
Our calculated, indirect band gap of 5.46 eV calls for further experimental investigations for the following reasons. While this value is clearly in the 4 to 5.5 eV range for available experimental results, it is desirable that new experiments attempt to narrow this rather wide (1.5 eV) range. Additionally, the first experimental report did not specify the direct or indirect nature of the band gap. The second one found a direct band gap while our result is an indirect band gap, as qualitatively found by most of the previous, theoretical calculations. Our calculated, room temperature direct gap of 6.056 eV, at Γ, is only 0.6 eV larger than our calculated, optical band gap of 5.46 eV. All measurement procedures do not find an indirect band gap; this fact is amply illustrated in the case of TiO2 that was considered to be a direct gap material until the work of Ekuma and Bagayoko  found an indirect one of 2.95 eV, very close to the direct one of 3.05 eV. Santara et al.  confirmed the prediction of Ekuma and Bagayoko and explained how measurements with non-polarized light could not find this indirect band gap, but they led to the slightly larger direct gap for TiO2.
In conclusion, we have performed ab-initio, self-consistent calculations of electronic energy bands, total density of states (DOS), partial densities of states (pDOS), effective masses and bulk modulus of zb-BeSe. The distinctive feature of our calculations as compared to previous ab-initio and empirical calculations is the implementation of the Bagayoko, Zhao, and Williams (BZW) method, as enhanced by Ekuma and Franklin (BZW-EF). Our calculated, indirect band gap is 5.46 eV, from Г to a conduction band minimum between Г and X, for a room temperature lattice constant of 5.152 Å. Our bulk modulus is in excellent agreement with experiment within the experimental uncertainties. Our results for the band gap, DOS, pDOS, effective masses and bulk modulus, along with similar ones from this group      , strongly suggest that LDA BZW- EF calculations have the capability to accurately describe and predict electronic and related properties of semi-conductors. Based on this capability, our calculations are expected to inform and to guide the design and fabrication of semiconductor-based devices. Our results and the limited, available experimental ones strongly suggest the need for further measurements of electronic properties of zb-BeSe, including its band gap.
This work was funded in part by the US Department of Energy (DOE), National Nuclear Security Administration (NNSA) (Award No.DE-NA0002630), the National Science Foundation (NSF) (Award No, 1503226), LaSPACE, and LONI- SUBR.