A Novel Stable Binary BeB2 phase

Potential crystal structures of BeB2 were explored using ab initio evolutionary simulations. A new phase with a Cmcm space group was uncovered. It was determined that the Cmcm phase is mechanically and dynamically stable and has a lower enthalpy, from ambient pressure up to 13 GPa, than any previously proposed phases, as measured using first-principles calculations. The crystal structure, phonon dispersion, phase transitions, and mechanical and electronic properties of this phase were investigated. It was determined that the Cmcm phase may transform into the phase at pressures higher than 13 GPa. The band structures and density of states reveal that the Cmcm phase is metallic. In addition, the Vickers hardness was calculated using three empirical models. To explain the origin of the hardness, charge density difference maps and a Mulliken population analysis were carried out, which demonstrated that there are strong covalent interactions between B atoms. By analyzing the Crystal Orbital Hamilton Population (COHP) diagrams, it was determined that the total interaction of the Be-B bonds is stronger than that of the B-B bonds, indicating a very complex bonding feature in the new phase. It was predicted that the new Cmcm phase is nearly absent of superconductivity.

including the MgB 2 crystal structure (space group P6/mmm, NO.191 20 , one unit per cell, prototype AlB 2 ) and the CaIn 2 and MgGa 2 structures (space group P63/mmc, NO.194, two units per cell 25,26 ); and (3) orthorhombic, including the SrAl 2 structure (space group Imma, NO.74, four units per cell, prototype CeCu 2 27 ), another MgGa 2 structure (space group Pbam, NO.55, eight units per cell 28 ) and the ''CaIn 2 '' structure (space group Pnma, NO.62 20 ). Many phases have already been found. At low pressure (under 160 GPa), the cubic phase (space group F 4 { 3m, NO.216) is the most stable, while the P63/mmc phase is dynamically unstable and can turn into low-symmetry Pnma. At pressures greater than 160 GPa, the P63/ mmc phase is more stable 2 . There is still an open question regarding whether there exist other stable phases with low energy in the low pressure zone. Therefore, further theoretical studies, such as of the crystalline structure, electronic properties and mechanical properties, are important for this potential material.
In the present work, potential crystal structures of BeB 2 were extensively explored using ab initio evolutionary simulations. A new phase with Cmcm symmetry was uncovered. The crystal structures, phonon dispersions, mechanical properties, phase transitions, and electronic structures of this phase as well as other available phases were investigated. The band structures and density of states reveal that the new phase is metallic. Further study on the mechanical properties using three different methods was used to predict the Vickers hardness of the newly discovered Cmcm phase. Charge density difference maps and Mulliken population analysis, as well as COHP diagrams, were also carried out to analyze the character of the bonding between the atoms in Cmcm BeB 2 . The superconductivity of the new Cmcm phase was also investigated.

Results
Crystal Structure. In a previous work, several phases were found both experimentally and theoretically. However, only four structures are stable according to the calculated phonon spectra in the present work. These are cubic (space group F 4 { 3m, NO.216), orthorhombic (space group Pnma, NO.62), orthorhombic (space group Pbam, NO.55) and hexagonal (space group P63/mmc, NO.194). In the present work, a new stable phase with space group NO.63 was uncovered by performing systematic structure searches for BeB 2 . Its optimized equilibrium lattice parameters, axial ratio c/a and unit cell volume are shown in Table 1. The crystal structure information of all other known dynamical stable phases are presented in Table 2 2,20 . Table 1 reveals that the conventional unit cell of the new Cmcm phase contains twelve atoms in total with two inequivalent Wyckoff positions 4c and 8f for Be and B atoms, respectively.
The schematic crystal structure of the Cmcm phase is shown in Figures 1, 2 and 3, where yellow and blue spheres represent Be and B atoms, respectively. Figure 1 shows the unit cell, while Figure 2 shows the super-cell of the same crystal structure. Figure 2(b) is perpendicular to the x-axis. There are two parallel columns formed by boroncentered and beryllium-centered coordination polyhedra that are marked by green and blue colors, respectively. Each boron atom is surrounded by four beryllium atoms that form a distorted tetrahedral configuration. Eight boron atoms embrace a beryllium atom, forming a complex polyhedron composed of 3 rectangles and 6 triangles. Along the z-axis, the structure can be viewed as several layers, as shown in Figure 2 (b) and (c). Six boron atoms connect to each other in a hexagonal ring, where 4 atoms are coplanar (denoted as B and B*) and the other 2 atoms are out of this plane (denoted as A and A*). Be atoms form another layer (denoted as C) between the B and B* layers along the z-axis. The upper image of Figure 2 (a) is taken from the A, B, C, B* and A* layers, whereas the lower image is composed of the B*, A*, C*, A and B layers. The distance between two A layers is 5.120 Å . The configurations of boron-centered and beryllium-centered coordination polyhedra are displayed in Figure 3. The distances of the three types of Be-B bonds are 1.923 Å , 2.021 Å and 2.105 Å , respectively, as shown in Figure 3 (a). The Boron-Boron bond lengths are 1.763 Å (the blue bond) and 1.657 Å (the green bond), as shown in Figure 3 (b).
Phonon Dispersion. The phonon dispersion curves of the new phase as well as of other previously proposed crystalline structures were calculated to determine their dynamical stability, as shown in Figure 4. The new Cmcm phase is dynamically stable, as there are no imaginary phonon frequencies detected in the whole Brillouin zone. For the previously proposed phases 2,20 , the Pbam, Pnma, F43m and P63/mmc phases are dynamically stable according to the same rule. However, the appearance of imaginary frequencies indicates that the P6/mmm, Fd 3 { m and Imma phases mentioned in previously published articles 2,20 are dynamically unstable.
Phase Transition. Figure 5 gives the calculated relative enthalpy as a function of pressure for all stable phases of BeB 2 , relative to the F 4 { 3m phase. The inserted figures show the stable Cmcm phase and the F 4 { 3m phase of BeB 2 , before and after the phase transition at 13 GPa, respectively. As mentioned before, the latter has been considered in previous works to be the most stable phase at atmospheric pressure, having a structure with a diamondoid boron network and beryllium atoms occupying the interstitial tetrahedral sites. The polyhedron representation of the F 4 { 3m phase has also been marked, as shown in the bottom-right corner of Figure 5. As shown in Figure 5, the enthalpy of the newly uncovered phase increases with pressure within the studied pressure range up to 50 GPa, while the enthalpies of the other three stable phases (Pnma phase, Pbam phase, P63/mmc phase) have an inverse trend. It can be concluded that the Cmcm phase is more competitive than the previously known most stable F 4 { 3m phase at pressures up to 13 GPa. When the pressure exceeds 40 GPa, the Pnma phase and P63/mmc phase have almost similar enthalpies, and both of them become more favorable than the new Cmcm phase. The cubic structure (F 4 { 3m phase) continues to have the greatest stability in the pressure range from 13 GPa to 50 GPa. These results are consistent with the results previously found by Andreas Hermann et al. 2 First-principles methods were used to calculate the formation enthalpy of the new phase Cmcm (No.63) at different pressures.
The structural parameters of each phase are listed in Table 3. The enthalpy of formation is calculated as follows: As shown in Table 4, when the pressure is 1 atm, the formation enthalpy of the Cmcm phase is 256 meV/atom, which is within the range of formation enthalpies (240 meV to 2125 meV) of some proposed models for BeB 2.79 phases 20 . This may explain why the BeB 2.79 phase with a very complex structure has been found experimentally, rather than the pure BeB 2 phase with the simple Cmcm structure. It also indicates that the proposed models of BeB 2.79 phases 20 with enthalpies less than 256 meV/atom are more competitive. Upon comparing the formation enthalpy of the new phase (2 0.056 eV/atom) with that of the most stable known phase F 4 { 3m (20.016 eV/atom), it may be observed that the former is more competitive than the latter. In short, the Cmcm phase is more thermodynamically stable than any other known BeB 2 phase.
Using the same approach, the formation enthalpy of the Cmcm phase at 10 GPa and 15 GPa was calculated to be 239.7 meV/atom   and 232.2 meV/atom, respectively. By comparing to the value at 1 atm (256.267 meV/atom), it may be seen that the formation enthalpy value increases continuously with the pressure, indicating that increasing pressure is not beneficial for the synthesis of the Cmcm phase. All of these conclusions agree well with those obtained from the pressure phase diagram 2,20 , and the results will provide guidance for experimental work on BeB 2 synthesis.
In addition, to provide more ''physics'' of the phase transition, we have calculated the solid-solid structural phase transitions between the Cmcm and F 4 { 3m phases at 13 GPa by using the VC-NEB method 29 . Figure 6 illustrate the snapshots from a dynamical trajectory collected from transition path sampling connecting the Cmcm phase and the F 4 { 3m phase. The energy barrier of the transition was calculated to be 0.25 eV/atom.

Discussion
Mechanical Properties. Mechanical properties, including the elastic constants, bulk modulus, shear modulus, Young's modulus, Poisson's ratio, are the main bases for choosing and designing materials. Therefore, all these properties were calculated using the first-principles approach in this work. The elastic constants of all structures under the ambient pressure C ij were calculated with the strain-stress method combined with Hooke's Law, as implemented in the CASTEP code 30 . The results are listed in Table 5. In order to determine the mechanical stability of the new predictions, C ij has to satisfy the elastic stability criteria 31 .
The mechanical stability condition and elastic constants of the orthogonal structure are positive. However, the following inequalities must also be satisfied to indicate stability: As shown in Table 5, these values all meet the mechanical stability criteria, initially confirming that the Cmcm phase is mechanically stable.
Based on the Voigt-Reuss-Hill approximation method 32,33 , we can find the appropriate bulk modulus B and shear modulus G using the elastic constants. In addition, the values of Young's modulus E and Poisson's ratio s can be calculated using the following formula: Usually, the bulk modulus B is used to characterize a material's resistance to volume deformation against external pressure, while the shear modulus G measures a material's ability to resist shear strain. Young's modulus E measures the resistance against longitudinal tensions. Among the 4 mechanical and dynamical stable BeB 2 phases, the Cmcm phase has the lowest value of bulk modulus and the second lowest shear modulus, revealing that it has low resistance to compression and shear strain. In addition, the elastic constant C 22 (282 GPa) is significantly smaller than C 11 (492 GPa) and C 33 (690 GPa), revealing that resistance along the b axis is much smaller than along the a and c axes.
In order to evaluate the ductility of the material, the B/G values were calculated and are listed in Table 5. Higher B/G values greater than 1.75 correspond to a ductile material, while values less than 1.75 correspond to a brittle material 34 . As shown in Table 5, the B/G value of the Cmcm phase is approximately 1.01, suggesting that it is very brittle.
Poisson's ratio reflects the strength of the covalent bond to some extent. A small Poisson's ratio (s50.13) indicates that the Cmcm phase is more intensely covalently linked than the Pbam phase (s50.29) and P63/mmc phase (s50.18). These results show that the Cmcm phase is likely to be a type of potential ultra-incompressible material, despite its bulk modulus and shear modulus not being very high.
Finally, the hardness of the Cmcm phase was calculated in three different ways. Using a recently proposed simple empirical hardness formula H v 52(G 3 /B 2 ) 0.585 -3 35 , the Vickers hardness of these phases was calculated. As the result shows, the Cmcm phase has a hardness value of 36.8 GPa, approaching the critical value of a super-hard material, 40 GPa. It also reveals that the predicted hardnesses of the Pnma phase and the F 4 { 3m phase are approximately 41.4 GPa and 42.3 GPa, respectively, implying that both of them are potential super-hard materials. The hardness value for the Cmcm phase is also calculated to be 25.2 GPa using th he formula of hardness given by Artem R. Oganov 36 (details can be found in the computational methods section).
The hardness of the new phase was also calculated by the microscopic hardness model proposed by F.M. Gao et al. 37,38 (details can be found in the computational methods section). The calculated parameters and hardness of the Cmcm phase calculated by this model are listed in Table 6. It is found that the total hardness of the Cmcm phase is only 13.8 GPa with this method.  Chen's model gives the highest values (36.8 GPa) from these three models followed by Oganov's model (25.2 GPa) while Gao's model provides the lowest values (13.8 GPa), which may be caused by their applicability to the boron compounds and with metals. The scattered values of hardness predicted from these three different models reveal that although it has a simple crystal structure, the electronic structure and bonding characters of the newly uncovered Cmcm phase of BeB 2 are quite complex, as will be analyzed below.
Electronic Properties. The electronic properties of the newly discovered phase, including energy band structures, total and partial density of states (DOS), and charge density maps, were calculated and are shown in Figure 7-9. Figure 7 shows the calculated band structure along high symmetry directions, as well as the total and partial DOS of the optimized Cmcm structure from first-principles calculations within the GGA scheme. The overlapping of the valence bands and conduction bands around the Fermi level suggests that the new phase has a clear metallic character, which is confirmed by the finite value of total DOS at the Fermi level. Figure 7 also plots the partial DOS of the novel Cmcm phase. It reveals that the total DOS of the upper part of the valence bands (from 29.5 eV to the Fermi level) is mainly contributed by the B-p state, while that of the conduction bands come from both the B-p and Be-p states. The Be-s state also contributes to the DOS of the conduction band (above 4 eV), while making a very slight contribution to that of the valence band. The B-s state contributes to the lower part of the valence band (less than 27.5 eV) and to the upper part of the conduction band, especially in the 4 eV to 8 eV region of the total DOS. As discussed above, there is significant hybridization of the s and p states from both Be and B in the region of 4 eV to 8 eV, implying the tendency to form covalent bonds between Be and B atoms. Figure 8 plots the total DOS for all stable phases at ambient pressure, corresponding to the Cmcm, F 4  Figure 8. Figure 9 shows electron density difference maps for the Cmcm phase on the selected slice (100) plane. It can be seen that electrons are gathered at the positions of the B atoms and especially between Boron-Boron bonds, as indicated by the region in red in the top left panel and its enlarged bottom left panel. After removing the atoms, it is clearly seen that electrons are transferred from the Be atoms to the B atoms, as shown in the right two panels. From the electron gathering region between Boron-Boron, one can speculate that there is a strong covalent interaction.
In order to give some insight into these bonding characters, an atomic and bond Mulliken population analysis, which can provide a  good way to quantitatively evaluate the charge transfer chemical bond strength of the studied system, was performed and analyzed by the CASTEP code, with the results listed in Tables 7 and 8. From Table 7, it can be clearly seen that electrons transfer from the Be atoms to the B atoms in the Cmcm phase, as the charges for Be and B are 0.69 and 20.34, respectively. Now let us turn to the bond Mulliken population analysis results. A high nonzero value of overlap population indicates that there is a strong covalent character of the bond, while a small value close to zero shows that there is weak or no interaction between two related atoms and a negative value indicates that the atoms cannot form a bond 39 . From Table 8, it may be seen that the bond population ranges from 0.20 to 1.63 for the Cmcm phase. The maximum number, 1.63, exists between Boron-Boron bonds, indicating their strong covalent character. However, the Be-B bond has very weak covalence, as its Mulliken population is only 0.20. The population of Be-Be is 20.24, suggesting that there are no bonds forming. All of these Mulliken population results are consistent with the conclusions drawn from the electron density difference maps. As discussed in the mechanical properties section of this work, the hardness of the Cmcm phase may be a result of the high population value of the Boron-Boron bonds.
Some techniques such as the crystal orbital overlap population (COOP) 40,41 and its analogous crystal orbital Hamilton population (COHP) 42,43 can provide a straightforward view of orbital-pair interactions; based on these techniques, it is possible to analyze and interpret the bonding situation in solid-state materials. To elucidate the bonding situations in this new BeB 2 phase, we performed crystal orbital Hamilton population (COHP) analysis, which partitions the band structure energy (in term of the orbital pair contributions) into bonding, nonbonding and anti-bonding energy regions within a specified energy range. Figure 10 shows the resulting -pCOHP as a function of energy for the new phase. Positive values of -pCOHP describe bonding energy regions, whereas negative values describe anti-bonding energy regions. As seen in the COHP diagrams in Figure 10, there appear to be obvious Be-B bonding states at the Fermi level, while that is not the case for the B-B combination, indicating that the interactions between Be-B bonds in the unit cell are stronger than those of B-B bonds, despite the fact that a single Be-B bond is weaker than a single B-B bond, as shown by the Mulliken population. Above or below the Fermi level (in the range of 26 to 6 eV), the COHP plots of the Be-B and B-B combinations are clearly dominated by bonding states, which shows that the new phase has a favorable stability performance. Superconductivity Properties. Encouraged by the relatively simple binary MgB 2 having a superconductivity transition temperature of 39 K 44 and the controversy regarding the reported superconductivity properties of BeB 2 and BeB 2.75 19,21 , we also calculated the superconductivity properties of the Cmcm phase of BeB 2 . The T c can be estimated from the Allen-Dynes modified McMillan equation 45 , which has been found to be highly accurate for materials with an EPC constant l,1.5 46 , where v log is the logarithmic average of the phonon frequency and m* is the effective Coulomb repulsion and was assumed to be constant at 0.1. The calculated spectral function a 2 F(v) and integrated l(v) of the Cmcm phase are plotted in Figure 11. Our results reveal that the Cmcm phase exhibits fairly low superconductivity properties, with a T c of only 0.1 K. These results shed light on the controversy regarding the reported superconductivity properties of BeB 2 or BeB 2.75 . The synthesized sample may contain both BeB 2 and BeB 2.75 phases. When the BeB 2 phase dominates, an absence of superconductivity would be observed, as shown in Ref. [21], while when the BeB 2.75 phase dominates, superconductivity appears.  Methods Ab initio evolutionary simulations were run using the USPEX (Universal Structure Predictor: Evolutionary Xtallography) code [47][48][49] . The USPEX code depends on VASP (Vienna ab initio simulation package) 50 to achieve global optimization to calculate the enthalpy of crystal structures and explore the lowest enthalpy phase of a given elemental composition. Here, we used the USPEX code to search for stable compounds and structures with a fixed chemical composition of Be n B 2n (n51 to 5); the Cmcm phase comes from the 25th structure with the stoichiometry of Be 2 B 4 . During the structure search, USPEX selects a whole range of 50 generations to calculate, with each generation containing 50 individuals. The settings used for the variation operators are as follows: 60% of each generation was used to produce the next generation by heredity, 20% comes from soft mutations, 10% is produced randomly from space groups, and the rest is produced through lattice mutations. The minimum length of any lattice vector was defined as 2.0 Å . The cutoff for USPEX relaxation and the k-points for resolution were 318 eV and 2p 3 0.02 Å 21 , respectively.
First-principles 51 calculations were carried out using the density functional theory (DFT) approach by applying a generalized gradient approximation (GGA) for the exchange correlation functional 52,53 . We applied the Ultrasoft pseudo-potential introduced by Vanderbilt 54 , and the k-point samplings in the Brillouin zone were   performed using the Monkhorst-Pack Scheme. The convergence tests used a kinetic energy cutoff of 600 eV and a k-point of 13 3 7 3 8 for the predicted Cmcm phase in the geometry optimization calculations. The self-consistent convergence of the total energy was 5 3 10 27 eV/atom, the maximum force on each atom was below 0.01 eV/ Å and the maximum atomic displacement was below 5 3 10 24 Å . The phonon dispersion curves were plotted using the super-cell calculation method 55 applied in the Phonopy program 56 . The calculation of the elastic constant and Mulliken overlap populations was carried out using the CASTEP code 57 . From the calculated elastic constants C ij , the polycrystalline corresponding bulk modulus B and shear modulus G were calculated using the Voigt-Reuss-Hill approximation 32,33 . In addition, Young's modulus E and Poisson's ratio s were obtained by the equations E5(9G?B)/(3B 1 G), and s5(3B22G)/(6B 1 2G), respectively.
In addition, we used three different methods to calculate hardness. These three methods were proposed by Xing Qiu Chen 35 , Artem R. Oganov 36 and Faming Gao et al. 37,38 , respectively. The description of Chen's model can be found in the text. To use the formula of hardness given by Artem R. Oganov 36 , we need to use the structure file POSCAR (which must contain an element symbol line) and set the parameters for goodBonds, valence and valence electrons. For main group elements, only the outermost electrons are considered as valence electrons under normal circumstances. Hence, the numbers of valence electrons for Be and B atoms are 2 and 3.
Details of Gao's model are described below. The Vickers hardness of complex crystals can be calculated by a geometric average of all bonds as follows, For the Cmcm phase, the total hardness is given by Because no d-orbital valence electrons are involved in the chemical bonds, the hardness of each bond for Cmcm phase can be expressed by: where d u is the length of the bond, N u e is the valence electron density (which can be calculated by where Z Be and Z B are the valence electron numbers of the Be and B atoms constructing Be-B or B-B bonds, N Be and N B are the nearest coordination numbers of the Be and B atoms, N j is the number of j bond in the unit cell, and V is the volume of the unit cell) and f u i is the Phillips ionicity of the bond. According to the generalized ionicity scale, the Phillips ionicity can be obtained from the following formula, where f h is the population ionicity scale of the chemical bond, p is the overlap population of the bonds, and p c is the overlap population of the bonds in a specified pure covalent crystal (here 0.57 is adopted). For the complex crystal compounds, we considered three effects on hardness: the covalent component, the ionic component and the small metallic component. First, we defined a factor of metallicity f m as n m /n e for a simple-structured compound, where n m and n e are the numbers of electrons that can be excited at the ambient temperature and the total number of valence electrons in the unit cell, respectively. According to the electronic Fermi liquid theory, the thermally excited electron number n m can be described by the product of D F and the energy width kT, where k is the Boltzmann constant and T is the temperature. At the ambient temperature, kT is equal to 0.026 eV. Therefore, f m can be written as: When the chemical bonds of a crystal are greater than or equal to two, we refer to it as a complex crystal. For the metallicity of complex crystals, the f m can be calculated by where n u A (or n u B ) is the thermally excited electron number of Be or B atoms in the utype bond, N MA (or N MB ) is the number of chemical bonds with a metallic component around Be or B atoms, and (n u e ) Ã is the number of valence electrons per u-type bond. To elucidate the bonding information in this new phase, we adopted a variant of the familiar COHP approach that stems from a PW calculation and was dubbed ''projected COHP'' (pCOHP) 58,59 . In this approach, all of the projection and analytic methods are implemented in a standalone computer program that processes PAW parameters and self-consistent results from VASP.   The calculation of the electron-phonon coupling (EPC) parameter l was performed using the pseudo-potential plane-wave method within the density functional perturbation theory (DFPT) 60 as implemented in the Quantum Espresso package 61 by using Von Barth-Car type norm-conserving pseudo-potential with cutoff energies of 80 and 320 Ry for the wave functions and the charge density, respectively. A 7 3 4 3 4 q-point mesh in the first Brillouin zone was used in the EPC calculation.

Conclusions
In summary, a new stable phase of BeB 2 with the space group of Cmcm was discovered by using ab initio evolutionary simulations. The Cmcm phase has a lower enthalpy than any previously proposed phase. The new structure is mechanically and dynamically stable, as determined by checking the calculated elastic constants and phonon dispersions, while several previously proposed phases (cubic: Fd 3 { m; hexagonal: P6/mmm; orthorhombic: Imma) were found to be dynamically unstable. The Cmcm phase may transform to the cubic F 4 { 3m phase when the pressure exceeds 13 GPa. The calculated electronic band structure and density of state suggest that the uncovered new phase is metallic. Scattered hardness values calculated from three models suggest the complex electronic and bonding features of the Cmcm phase. The charge density difference maps and the Mulliken population analysis reveal that there are strong covalent interactions between the B atoms. The COHP diagrams show that the total interaction of Be-B bonds is stronger than that of B-B bonds. The Cmcm phase exhibits fairly low superconductivity properties, with a calculated T c of approximately 0.1 K. The current theoretical predictions will most likely promote further experimental and theoretical investigation on the Be-B system.