Phase transitions, mechanical properties and electronic structures of novel boron phases under high-pressure: A first-principles study

We have explored the mechanical properties, electronic structures and phase transition behaviors of three designed new phases for element boron from ambient condition to high-pressure of 120 GPa including (1) a C2/c symmetric structure (m-B16); (2) a symmetric structure (c-B56) and (3) a Pmna symmetric structure (o-B24). The calculation of the elastic constants and phonon dispersions shows that the phases are of mechanical and dynamic stability. The m-B16 phase is found to transform into another new phase (the o-B16 phase) when pressure exceeds 68 GPa. This might offer a new synthesis strategy for o-B16 from the metastable m-B16 at low temperature under high pressure, bypassing the thermodynamically stable γ-B28. The enthalpies of the c-B56 and o-B24 phases are observed to increase with pressure. The hardness of m-B16 and o-B16 is calculated to be about 56 GPa and 61 GPa, approaching to the highest value of 61 GPa recorded for α-Ga-B among all available Boron phases. The electronic structures and bonding characters are analyzed according to the difference charge-density and crystal orbital Hamilton population (COHP), revealing the metallic nature of the three phases.

B oron has been recognized as a mystical and important element due to its fascinating chemical and physical properties 1 . However, most of the reported phases are likely to be boron-rich borides [2][3][4] and probably only four phases correspond to the pure element. This might be ascribed to the extreme sensitivity of the element to even small amounts of impurities, given that there are a number of boron-rich compounds with unique icosahedral structures such as YB 66 , B 6 O, AlB 12 , B 13 P, B 50 C 2 5 . Such a complexity may arise from its unique location in the periodic table: situates at the boundary between metals and nonmetals and is the only nonmetal element of group III elements. Extensive studies have been performed to investigate the crystal structures and phase stabilities of this element [6][7][8][9][10][11][12] . Most of the reported boron has complicated crystal structures based on icosahedral B 12 clusters, which can be flexibly linked into rigid frameworks. Such a special structure makes the boron have unique properties among elemental materials, such as a low volatility, a high melting point (2450uC) as well as excellent mechanical properties including high strength and hardness 13 . Now, only three of the reported boron phases have been confirmed to be thermodynamically stable, i.e., arhombohedral boron (a-B) 14 , b-rhombohedral boron (b-B) 15,16 and c-B 17,18 . The a-B was identified in 1958 19 , featured by a single icosahedron in a slightly distorted cubic cell with the C 3 axis of the icosahedron aligned to the c-axis of the unit-cell 20 . In 1957, Sands and Hoard announced the identification of the b-B 21 . Experiments 9,22 and theories 23,24 are thereafter carried out and tested for this phase. Shirai et.al. [25][26][27] studied diagram of a-B and b-B and found that at zero temperature a-B is more stable than b-B. However, the stability of these two phases remains uncertain at low pressure. Oganov et al. 5,8,18 reported an ionic phase of element boron (c-B), which consists of icosahedral B 12 clusters and interstitial B 2 pairs acting as anions and cations, respectively, in a NaCl-type arrangement. The ionicities of boron-boron bonds in B 12 icosahedra of a-B have also been reported in one of our previous work 28 .
In addition, other phases have also been reported. The first crystal structure of a-tetragonal boron called T-50 was identified in 1951 29 , constituted by 50 atoms in the unit cell. It was proposed that b-B is stable up to 30 GPa and 3500 K, and at higher temperatures and pressures, a phase transition to the tetragonal 'T-192' structure occurs 30 . Recently, Pickard et al. 8 obtained a new metastable phase by an ab initio random structure searching method and this new phase can be viewed as a polymorph of a-B, differing in the connectivity of the icosahedral units and just 0.01 eV/atom less stable. He et al. 31 subsequently studied the properties of this phase (termed a*-B) by first-principles calculations. Also, the a-Ga structure and its variants of boron also stimulate great interest in probing the high-pressure phase 5,18,32 .
In this work, we conceived three new phases with or without icosahedral unit of elemental boron. The crystal structure, phase stability, mechanical and electronic properties of these phases have been systematically studied under high-pressure based on the firstprinciples method.

Results
Structural features. The schematic ambient crystal structure of the first new phase is shown in Figure 1, suggesting an orthorhombic structure with the space group Pmna. The unit cell is composed of 24 atoms, which can be classified into 4 non-equivalent atomic sites stated as B1, B2, B3, and B4. For the convenience of further discussion, we refer to it as o -B 24 . Structurally, o -B 24 is another phase of boron consisting of pure slightly distorted icosahedral B 12 clusters (except a-boron and a*-boron). The B 12 cluster units are equivalent, built by 4 B1 atoms, 4 B2 atoms, 2 B3 atoms, and 2 B4 atoms. The centres of the B 12 icosahedra projected along the [010] direction (overlapped by B4 sites) form a rhombic packing2see Figure 1. In o -B 24 , the four inequivalent boron atoms can be  divided into two types: the B2 and B4 atoms have five intraicosahedral B-B bonds and another single inter-icosahedral B-B  bond while the B1 and B3 atoms possess five intra-icosahedral B-B  bonds and double inter-icosahedral B-B bonds. See Table S1 in the Supporting Information for more details. The relation between intraicosahedral and inter-icosahedral bond has been intensively studied. For example, it has been generally thought that the inter-icosahedral bonding is characterized as the covalent bond and is stronger than the intra-icosahedral bonding since Longuet-Higgins disclosed the nature of the icosahedral bonding of bond molecules 33 . However, apparent contradictions have been found between such basic bonding concepts and experimental results, which may be solved by considering the special geometrical effects, as introduced to explain the elastic responses of boron carbides 34 . In table S1, we found that the shortest bonds are inter-icosahedral B 2 -B 2 bonds (1.65 Å ), indicating they are stronger than each intra-icosahedral B-B bonds. There are also some weak inter-icosahedral bonds, like B 1 -B 1 (1.92 Å ) and B 1 -B 3 (2.04 Å ), as well as some inter-icosahedral bonds comparable to those intra-icosahedral bonds, like B 4 -B 4 (1.70 Å ). To summarize, the relation between intra-icosahedral and intericosahedral bonds in the o -B 24 phase is quite complicated, implying it may have some special elastic responses. In addition, the coordination number N c of a selected boron atom (defined by the number of surrounding bonds with a bond length less than 2 Å ) has also been investigated in this work. For o -B 24 , a coordination description of B1 [6] B2 [6] B3 [5] B4 [6] giving an average coordination numbers 5.8. In addition, each icosahedral B 12 cluster is surrounded by eight neighbouring B 12 clusters connected by different B-B bonds. Like a*-boron, o -B 24 may also be considered as another twinned polymorph of a-boron which contains pure icosahedral B 12 clusters without linking atoms. The slightly difference in the atomic positions within an icosahedral B 12 cluster and the difference in the packing way among B 12 clusters in a-boron, a*-boron, and o -B 24 might explain their distinct symmetries. It is noted that two other known phases (c -B and b -B) of boron also contain icosahedral B 12 clusters. In c-B, each icosahedral B 12 cluster connects to ten equivalent neighbouring B 12 clusters directly with inter-B 12 B-B bonds and two secondly adjacent B 12 clusters indirectly   by interstitial B 2 pairs 31 . The connecting between icosahedral B 12 clusters and interstitial atoms in the b-B phase is more complicated.
The second new structure2denoted as c-B 56 2belongs to cubic crystal system and contains 56 atoms in the conventional unit cell (space group Ia3, No. 206). In this structure, the ordered B atoms are clearly split into two categories, B1 and B2, which occupy Wyckoff 48e and 8a positions, respectively. The B1 atoms form a cubic area and the B2 atoms occupy its centre when the structure is projected along the [100] direction, which is denoted by the dashed box shown in the left panel of Fig. 2. The middle panel of Fig. 2 is the dashed box view along [010] direction. One can find that each B1 atom is coordinated by one B2 and five B1 atoms and each B2 atom by six B1 atoms (forming a corrugated hexagon) 2 see the right panel of Fig. 2. In short, the repeated cubic unit denoted by the dashed box in Fig. 2 can be considered as stacking congruent distorted hexagons linked by 2 types of B1-B1 bonds with bond lengths of 1.70 Å (B 1 -B 1 ) and 1.76 Å (B 1 -B 1 9), see Table S2. There is another type of B1-B1 bond length with 1.82 Å (B 1 -B 1 0), forming the edge of the distorted hexagon. In the unit cell, one find the central three B2 atoms are bonded with surrounding twelve B1 atoms with a single bond length of 1.74 Å . Regarding the inter-cubic areas, only two types of B1-B1 bonds of 1.70 Å (B 1 -B 1 ) and 1.76 Å (B 1 -B 1 9), are found to connect them of each other, see Table S2 in the Supporting Information for details.
The third new phase is a monoclinic structure that contains 16 atoms in the unit cell (space group C2/c, No. 15, hereafter denoted by m-B 16 ). B atoms have two inequivalent sites in this structure, both of which occupy Wyckoff 8f positions. The m-B 16 structure is constructed by two layers, named A and B in Figure 3 Table S3 in the Supporting Information. From Table S3 and Fig. 3, the A layer is bonded to B layer by two different B1-B1 bonds with bond length of 1.73 Å and 2.13 Å , respectively, while only one type of B2-B2 bond with bond length of 1.78 Å between B layer and the next A layer. It has also been found there are six different types of intra-layer bonds from Table S3. In short, each B1 atom is bonded with four B1 atoms and four B2 atoms (seven types of B-B bonds) while each B2 atom is bonded with four B1 atoms and three B2 atoms (six types of B-B bonds). Moreover, the average coordination number is calculated to be 6.0 in the m-B 16 phase according to the aforementioned definition. Interestingly, when the pressure exceeds 68 GPa, the m-B 16 phase transforms into another new phase. The new phase2denoted as o-B 16 2 belongs to orthorhombic crystal system (space group Imma, No. 74, see Figure 3b) and also contains 16 atoms in the unit cell. In this structure, the B atoms are also packed in a layer structure with the stacking order of KLKLKL… along the crystallographic c axis. Here both K and L denote the hexagonal boron networks, i.e., the    Phase transformations. According to our first-principles calculation results, the relative enthalpies for the chosen structures compared to that of a-B as a function of pressure up to 120 GPa are presented in Figure 4. It is confirmed that the a-B is more favorable than any other boron phases at ambient pressure. The relative enthalpy of a*-B has a constant trend within the pressure range studied, that is, the enthalpy of a-B is only about 0.01 eV/atom lower than that of a*-B during the whole pressure range. When pressure increases above 19 GPa, the c -B has a lower enthalpy than that of a-B and becomes stable. The c -B keeps stable up to 93 GPa and then the a-Ga-B prevails. Our prediction results are in accordance with the earlier study results 18,31 . The present calculation also shows that the enthalpy curve of the m-B 16  It has been widely known that pressure and temperature are both pivotal factors that determine the states of materials. During past   Space group  Space group decades, the pressure limit for laboratory experiments has been progressively enhanced. In 1978, the pressure in the experiments with diamond2window pressure cell exceeded 178 GPa 35 . In 1986, a diamond2anvil, high-pressure apparatus was used to extend the upper limit from 210 to 550 GPa 36 , which is beyond the maximum pressure (360 GPa) of the earth9s core. In present work, pressure as high as 500 GPa have been applied in predicting the stability region of the c-B 56 and o-B 24 phases, unfortunately, no new stable phases have been found. We also try to investigate the high temperature stabilities of the new phases by performing quasiharmonic free energy calculations. As showed in Figure S4 of the supplementary file, it is found that the m-B 16 phase transforms to the o-B 24 phase when the temperature up to 575 K at ambient pressure, suggesting temperature plays a very important role during phase transition. However, in the case of o-B 24 , there are 200 phonon dispersions are required for a pressure point, such computational resources are not available in our group now, and will be explored in the future work.
To study the dynamical stability of the new high-pressure phases of boron, the phonon properties are investigated by phonon package 37 . The calculated phonon dispersion curves are shown in Fig. 5 (a),(d), respectively. We can see that there are no imaginary frequencies for them, indicating that they are all dynamically stable.
Mechanical properties. Several fundamental solid-state properties, such as equation of state (EOS), specific heat, thermal expansion, Debye temperature, Grüneisen parameter, melting point and many others are closely related to elastic properties of solids and, thus, the knowledge of elastic constants C ij is essential for investigating the mechanical and thermodynamic properties of a system. In calculating the elastic constants, different types of deformations are adopted for different phases according to the space group symmetry as implanted in the CASTEP code 38 . Table 5 summarizes the C ij for the new boron phases. There are three, nine or thirteen independent elastic constants for the new cubic, orthorhombic and monoclinic crystal systems, respectively. The criteria for mechanical stability 39 of cubic phases are given by For the orthorhombic crystal, the corresponding mechanical stability criterion is as follow: C ii w0 i~1, 2, 3, 4, 5, 6 ð Þ C 11 zC 22 zC 33 z2 C12zC 13 zC 23 ð Þ w0 For the monoclinic structure, the mechanical stability under ambient pressure can be judged by where V 0 is the volume per formula unit at ambient pressure, V is the volume per formula unit at pressure P given in GPa, B 0 is the isothermal bulk modulus, and B 0 9 is the first pressure derivative of the bulk modulus. The values of the bulk modulus B 0 and its pressure derivative B 0 9 are listed in Table 6. The bulk modulus value is in excellent agreement with that calculated from elastic constants, confirming the reliability of our calculations. Figure 6 (a) and (b) plot the pressure dependence of the lattice constants a, b, and c for o-B 24 and m-B 16 up to 50 GPa. It can be seen that the b axis is the most compressible crystallographic direction for both structures. One can also notice that the pressure has an even stronger effect on the b axis direction of o-B 24 than on that of m-B 16 . Concerning the o-B 24 , the c axis is the most incompressible crystallographic direction and followed by the a axis. Whereas in the m-B 16 structure, a and c axes exhibit similar compressibility at relatively lower pressure, the c axis becomes rigid with a lower compressibility than a axis as the pressure further increases. Interestingly, the c axis direction has an even lower compressibility than the b axis at over 47 GPa. The pressure dependence of cell volume is also shown in the insets. Comparatively, the rate of the volume shrinkage for o-B 24 (about 19.7%) is larger than that of m-B 16 and c-B 56 (by 15.3% and 16.5%, respectively), resulting in a remarkably smaller bulk modulus than those of other two new structures (see Table 6).
Based on the Voigt-Reuss-Hill approximation 41 , the corresponding bulk and shear moduli (B and G) are obtained from the calculated elastic constants. The Young's modulus (E) and the Poisson's ratio (s) can also be calculated from B and G. The values of B, G, E, s and B/G for the new and known boron phases are all illustrated in Table 6. The bulk modulus is a measure of the resistance against volume change imposed by the applied pressure, while the shear modulus denotes the resistance against the reversible deformations upon shear stress 42 . The calculated shear modulus for a-Ga-B (297 GPa) is higher than the others, indicating it can withstand higher shear stress than other structures, followed by the new phase (o-B 16 ) with a shear modulus of 287 GPa. On the contrary, the o-B 24 phase has the lowest shear modulus. Young's modulus is a measure of the stiffness of the solid. From Table 6, the a-Ga-B is also of the stiffest structure, followed by o-B 16 Table 6 | B 0 (GPa) (is the isothermal bulk modulus, obtained by fitting a third-order Birth-Murnaghan equation of state), B 0 9, V 0 (Å 3 /f.u.), bulk modulus B (GPa) (is the adiabatic bulk modulus, obtained by the DFT calculation), shear modulus G (GPa), Young9s modulus E (GPa), Poisson9s ratio s and B/G ratio for the whole structure types of boron at 0 GPa and 0 K   204 GPa) 48,49 , m-B 16 has a very similar bulk modulus, but its shear modulus is 26.5% higher. Therefore, it is conceivable that m-B 16 is also a superhard material considering its higher shear modulus than that of B 6 O. Thus we analyzed the hardness of different phases by adopting the empirical scheme 50 which correlates the Vicker9s hardness and the Pugh9s modulus ratio via the formula According to Eq. (5), the obtained values of Vicker9s hardness are illustrated in Table 7. From the data, m-B 16 has a hardness of 56.2 GPa, which is higher than the hardness of other boron phases except for that of a-Ga-B 51 . Note that the hardness of the orthorhombic o-B 16 is 60.7 GPa, suggesting the high-pressure phase transition may result in high hardness values. To give a range of possible values, we also use the Lyakhov-Oganov method 52 to evaluate the hardness of the whole boron phases. From Table 7, one can find that there is about 20 , 30% difference between the calculated hardness values from the two approaches. Interestingly, we find that both m-B 16 and o-B 16 as well as a-Ga-B are potential superhard materials. Our results are in good agreement with existing experimental and theoretical data in relevant references as shown in Table 7.
Electronic properties. To probe the electronic properties of the four boron phases, we calculated the electronic band structure and the partial electronic density of states (DOS) for o-B 24 , c-B 56 , and m-B 16 at 0 GPa and o-B 16 at 68 GPa, as shown in Figure 7. It can be clearly seen that the four phases are metallic as the energy bands crossing over the Fermi level (E F ), which are quite different from most known boron phases. It is generally known that boron has many hypothetical structures, and some of them are also predicted to be metallic [53][54][55] . However, none of these phases have been confirmed experimentally. Many scholars have made unremitting efforts to understand difficulties of boron crystals by the experimental and theoretical methods. Such difficulties maybe attributed to the   defect states, which are found commonly in complicated boron crystals 24,56,57 . According to Ref. 57, boron crystals show another freedom 2 compositional freedom for complicated crystals, which causes deviation from stoichiometry, such effect has also been found in boron compounds like BeB 2 [58][59][60][61][62] . A possible consequence is that doping is countered by change of defect density in a way that the system is kept semiconducting, such as the known b-phase. In case of proposed structures, since the host structures are band crossing metal, above scenario does not apply. Instead, introduction of defects may lead to the shift Fermi level in a way that density states at the Fermi level are minimized.
The total electronic DOS at the Fermi level for the m-B 16 and o-B 16 phases is 0.55 and 0.21 eV 21 , respectively, both of which are smaller than that (0.64 and 1.49 eV 21 , respectively) of the o-B 24 and c-B 56 phases. Moreover, we can see that both conduction and valence band for the whole phases are contributed from the p states, especially around the Fermi level. Figure 8 plots the difference charge density maps of the four new structures along selected planes at ambient pressure. In the present work, the charge density difference is defined as: where r sc is the surface charge density obtained after self-consistent calculations. r atom is the surface charge density obtained by corresponding non-self-consistent calculations, namely, it is a superposition of atomic charge densities of the original geometry configuration. Therefore, the charge density difference maps can be used to analyze the charge transfer before and after the electronic structure relaxation as well as the charge transfer during the corresponding bonds forming procedure. Before carrying out such calculations, the validity of our calculation procedure has been confirmed by reproducing the differential charge density maps as shown in Fig. 4 of Ref. 31.
For the convenience of discussion, the locations of boron atoms are also marked. For the phase o-B 24 containing icosahedra, one observes that there are electrons distribute on the centre of boron pentagon2see the centre of Figure 8 (a), indicating that the charge accumulation occurs inside the icosahedral B 12 units. It is obvious that a smaller charge accumulation around the B atoms (the pentagon forming ones), and some charges move from adjacent atoms forming the pentagon (the green area) to the center of the pentagon (the red area). It is worth mentioning that the nearest B-B distance is 1.65 Å , which is also the nearest distance between two adjacent icosahedra. Significant charge depletion (about 0.045 e/Å 3 ) between the shortest B-B bonds in the differential charge density maps has been observed, revealing some charges are moved out of the bonds during relaxation. Such results contradict to those of a-Boron and a*-Boron, where the bonding electrons prefer to distribute on the inter-B 12 B-B bonds. Such abnormal electronic response of the o -B 24 phase may also be attributed to the complicated relation between intra-icosahedral and inter-icosahedral bonds as we discussed before. In addition, some high charge accumulation in the interstitial space (marked by 1 , 5 in Fig. 8(a)) among icosahedra has been observed, implying gaining electrons during bonding therein.
For the c-B 56 structure, the charge density difference on the (111) atomic plane is obtained. The charge depletion regions form regular hexagons, indicating the neighbouring B atoms perpendicular to the plane is losing electrons during relaxation. Compared to the other B atoms, the B atoms located on the edge of the slice show substantial charge accumulation. A small amount of electron distribution in the interstitial sites around the centre of this map can be also observed.
For the m-B 16 and o-B 16 structures, the charge density difference maps for the same slices are obtained as shown in Figure 8  Some techniques such as the crystal orbital overlap population (COOP) 63,64 and its analogous crystal orbital Hamilton population (COHP) 65,66 can provide a straightforward view for oribital-pair interactions. Also, based on the techniques, one can analyze and interpret the bonding situation in solid-state materials. To elucidate the bonding situation in these four new boron phases, we also perform crystal orbital Hamilton population (COHP) analysis, which is a bond-detecting tool for solids and molecules. COHP partitions the band structure energy (in term of the orbital pair contributions) into bonding, nonbonding and antibonding energy regions within a specified energy range. In Figure 9    the crystal orbital Hamilton population (COHP) analysis one find that the o-B 24 phase has stronger bonds than that of the c-B 56 phase and the interaction changes during the m-B 16 phase transforming to the o-B 16 phase. The present work strongly suggests further study is needed to explore the behaviours of boron under high pressure conditions.

Methods
Although evolutionary simulation method like USPEX [67][68][69] have been adopted in successfully predicting potential crystal structures of Boron, all new structures studied in this work were conceived and constructed by hand, which is activated by Chaoyu He9s work 31 . The details of procedures in arriving these structure are described in the supplementary file.
For each crystal structure, the structural relaxations and electronic properties calculations were performed in the framework of density functional theory 70 , as carried out within the Vienna ab initio simulation package (VASP) 71,72 , with the projector augmented wave (PAW) method 73 . The 2s 2 2p 1 electrons were treated as valence electrons. The generalized gradient approximation with the Perdew-Burke-Ernzerhof (PBE) functional 74 for the exchange correlation was employed. A planewave basis with a cutoff energy of 500 eV was used to expand the wave functions. To ensure that the obtained structures are dynamically stable, phonon frequencies were calculated throughout the Brillouin zone using the phonon package 37 with the forces calculated from VASP. The reliability of the pseudopotential approach has also been confirmed by the full-potential linearized augmented plane waves approach. The calculation of the elastic constants by the strain-stress relations was carried out using the CASTEP code 38 .
Some techniques such as the crystal orbital overlap population (COOP) 63,64 and its analogous crystal orbital Hamilton population (COHP) 65,66 can provide a straightforward view onto oribital-pair interactions, based on these techniques, one can analyze and interpret the bonding situation in solid-state materials. In the present work, to elucidate the bonding information in these four new boron phases, we adopted a variant of the familiar COHP approach that stems from a PW calculation and was dubbed ''projected COHP'' (pCOHP) 75,76 . In this approach, all the projection and analytic methods have been implemented in a standalone computer program which processes PAW parameters and self-consistent results from VASP.