Prediction of superconductivity in pressure-induced new silicon boride phases

The crystal structures and properties of boron-silicon (B-Si) compounds under pressure have been systematically explored using particle swarm optimization structure prediction method in combination with ﬁrst-principles calculations. Three new stoichiometries, B 2 Si, BSi, and BSi 2 , are predicted to be stable gradually under pressure, where increasing pressure favors the formation of silicon rich B-Si compounds. In the boron-rich compounds, the network of boron atoms changes from B 12 icosahedron in the ambient phases to the similar buckled graphenelike layers in the high-pressure phases, which crystalize in the same P ¯3 m 1 symmetry but with different numbers of boron layers between adjacent silicon layers. Phonon calculations show that these structures might be retained to ambient conditions as metastable phases. Further electron-phonon coupling calculations indicate that the high-pressure phases of boron-rich compounds might superconduct at 1 atm, with the highest T c value of 21 K from the Allen-Dynes equation in P ¯3 m 1 B 2 Si, which is much higher than the one observed in boron doped diamond-type silicon. Moreover, further fully anisotropic Migdal-Eliashberg calculations indicate that B 2 Si is a two-gap anisotropic superconductor and the estimated T c might reach up to 30 K at 1 atm. On the silicon-rich side, BSi 2 is predicted to be stable in the CuAl 2 -type structure. Our current results signiﬁcantly enrich the phase diagram of the B-Si system and will stimulate further experimental study.

The B-Si system has been the subject of extensive investigation since 1900 [20]. At ambient pressure and high temperature, three approximate compositions, B 3 Si [14,[21][22][23], B 6 Si [24,25], and B n Si (n 14) [26][27][28], have been observed in experiments. α-B 3 Si, with stoichiometries * xieyu@jlu.edu.cn † gaoguoying@ysu.edu.cn ranging from B 2.8 Si to B 4 Si, containing boron-rich icosahedra, was proposed to be rhombohedral, where silicon atoms partially substitute boron atoms in the polar icosahedral sites [22]. This compound forms from direct combination of the elements by isothermal heating at 1225 • C or arc melting in a high-frequency induction furnace. The first fully crystallographically ordered structure, β-B 3 Si [14], was produced in Ga flux at a lower temperature (850-1000 • C), and adopts an orthorhombic structure with B atoms forming regular B 12 icosahedra. Another stoichiometry with an approximate ratio of B:Si = 6:1 (B 6 Si) is also known, which contains boron-rich icosahedra, icosihexahedra, as well as isolated silicon and boron atoms. α-B 3 Si, β-B 3 Si, and B 6 Si are all semiconductors with band gaps of 0.2, 2.0, and 1.2 eV, respectively [14,15]. As an extension of solid solution of silicon in boron, the structure of B n Si (n 14) is similar to pure β boron [26,27], with silicon atoms partially occupying icosahedral boron atom positions or interstitial sites. It is well known that the physical properties of a material are closely related to its crystal structures. The experimentally observed structures of B 3 Si, B 6 Si, and B n Si contain similar B polyhedral frameworks, but there is no report on other structural features formed by B atoms in boron-rich compounds yet. On the silicon-rich side, there is only one theoretical work on BSi 3 by Wang et al. [29], where they predicted three metastable structures, P3 1 21, C2/m, and P2 1 /m, in the pressure range of 0-160 GPa, which are potential hard materials. Small-size B n Si 0/− clusters and two-dimensional boron-silicon compounds were also theoretically investigated [30][31][32][33][34].
Here, we have extensively investigated the crystal structures and properties of B-Si compounds under pressure by using the particle swarm optimization structure prediction method in combination with first-principles calculations. We first explored the already known B 6 Si, B 4 Si, and B 3 Si stoichiometries and then proceed to the analysis of hypothetical B 2 Si, BSi, BSi 2 , and BSi 3 stoichiometries. We studied in detail the crystal structures of the preferred compounds as a function of pressure, together with their dynamical stability and corresponding electronic band structures and superconductivity. We found that three new stoichiometries, B 2 Si, BSi, and BSi 2 , become stable gradually under pressure, where increasing pressure favors the formation of silicon rich B-Si compounds. All the predicted stable high-pressure structures of boron-rich compounds have similar structural characteristics containing buckled graphenelike boron layers. In addition, our calculations indicate that these structures might be preserved to ambient pressure as metastable phases, in which they were predicted to be superconductors with the highest T c value of 30 K in P3m1 B 2 Si.

III. RESULTS AND DISCUSSION
The crystal structure searches of B x Si y (x = 2-4, and 6, y = 1; x = 1, y = 1-3) were carried out with system sizes containing up to 8 formula units (f.u.) per simulation cell over a wide pressure range of 0-400 GPa. Our structure prediction starts from the already known compounds at low pressure. For B 3 Si, an orthorhombic Pmma structure was predicted to be stable against decomposition into elements (B and Si). Its structure is very similar to the experimentally observed Imma β-B 3 Si (Fig. S1 [59]). Both structures are composed of interlinked B 12 cages and Si 4 chains. However, there is a little difference in the relative positions of some Si atoms (atoms 2 and 4) [ Fig. S1(a) [59]]. We made an estimation of the energy barrier between them by shifting Si atoms and the calculated maximum energy barrier reaches 0.12 eV/f.u. at 0 K [ Fig. S1(b) [59]], which prevents the transition from Pmma to Imma. However, this barrier can be easily overcome at high temperatures so that the two structures might be interchangeable.
For B 4 Si and B 6 Si, a series of energetically competitive structures (Figs. S2-3 [59]) were predicted to contain B-rich cages, in agreement with the experimental observations. The predicted R3m B 6 Si is isomorphous with the well-known R3m B 6 O [60]. Cmcm B 6 Si is the most stable phase between 0 and 17.5 GPa and then transforms to a similar Cmca structure, which is the same with the Cmca LiB 6 Si getting rid of Li atoms [61]. Different from B 3 Si, all the predicted structures for B 4 Si and B 6 Si are metastable with positive formation enthalpies, though the enthalpies of P2/c B 4 Si are only 3-11 meV higher than those of 4B+Si at the pressure range of 0-10 GPa. Notably, all the already known B-Si compounds, with the exception of β-B 3 Si, were observed to be disordered and the prediction of disordered structures goes beyond the scope of our research. However, the successful prediction of β-B 3 Si and B 4 Si and B 6 Si cagelike structures confirms the reliability of our method.
To explore the possible new stoichiometries, we focused our structural searches on B x Si y at 1 atm, 10, 50, 100, 200, 300, and 400 GPa. The thermodynamic stabilities were evaluated by calculating their formation enthalpies relative to the products of dissociation into constituent elements. As references, the R3m, Pnnm, Cmca, and P6 3 /mcm structures for B and Fd3m, P6 3 /mcm, and Fm3m for Si are considered, each in its most stable form at the specified pressure [62][63][64]. To account for the enthalpy relationships of each stoichiometry, a convex hull for the B-Si system was calculated at different pressures, as depicted in Fig. 1(a). Formation enthalpies of compounds are located on the hull, showing they are thermodynamically stable with respect to elements or other binary compounds. At 50 GPa, the already known B-rich compounds are predicted to decompose into elemental B and Si, and no other stoichiometries can be formed. With increasing pressure to 65 GPa, B 2 Si, B 3 Si, and B 4 Si become energetically favored over elemental dissociation, although B 3 Si is still unstable against decomposition into B 2 Si and B 4 Si. As pressure increases to 150 GPa, the formation enthalpies of B 2 Si, B 3 Si, and B 4 Si become more and more negative (more stable) and eventually B 3 Si falls on the convex hull. Upon further compression, at 250 and 300 GPa, BSi and BSi 2 , respectively, gradually become stable with respect to disproportionation into other stoichiometries. Pressure clearly favors the formation of B-Si compounds with a high content of silicon. Interestingly, as pressure increases from 65 to 250 GPa, B 2 Si becomes easier to form. Similarly, pressure also plays a positive role in the formation of BSi.
The predicted pressure-composition phase diagram and respective enthalpy curves of corresponding B-Si compounds are shown in Fig. 1(b) and Fig. S4 [59], respectively. An orthorhombic Pnma BSi is predicted to be stable between 233 and 300 GPa and transforms into an energetically favored Cmcm structure above 300 GPa. R3m and P3m1 B 2 Si are enthalpically stable at pressure ranges of 60-170 GPa and 170-400 GPa, respectively. The experimentally observed β-B 3 Si and the similar theoretically predicted Pmma structure are stable up to 18 GPa. Note that B 3 Si decomposes into pure elements and other compounds at the pressure range of 18-150 GPa, while it becomes stable again with a P3m1 structure above 150 GPa. The P3m1 structure for B 4 Si is stable over a wide pressure range of 65-380 GPa and I4/m BSi 2 will be stable at pressures between 282 and 400 GPa.
As a basic requirement for the structural stability, the dynamical stability of various structures for B-Si compounds is also examined. The phonon dispersion curves of all the predicted structures in the corresponding stable pressure ranges are presented in Fig. S13 [59]. The absence of any imaginary phonon frequency in the entire Brillouin zone shows they are dynamically stable. Interestingly, R3m and P3m1 B 2 Si, C2/m and P3m1 B 3 Si, and P3m1 B 4 Si were also calculated to be dynamically stable at 1 atm, which indicates that these predicted high-pressure phases might be also retained to ambient pressure.
The predicted stable structures are shown in Fig. 2. Both Pnma and Cmcm structures of BSi are orthorhombic and have the same configurations as the already known Pnma FeB and Cmcm CrB, respectively [2]. In both structures, B atoms form zigzag chains and Si atoms form a column of trigonal prisms arranging along the B-B chains. With increasing B content, the B-B bonding changes significantly. As shown in Figs. 2(c) and 2(d), the zigzag B-B chains in BSi transform to buckled six-membered rings in the predicted structures for B 2 Si. The P3m1 B 2 Si is similar to the structure of MgB 2 . Si atoms are located at edges of the unit cell and a deformed graphenelike boron layer is located between Si layers. The R3m B 2 Si can be considered a slight modification of the P3m1 structure, containing two types of boron atoms with coordination numbers of 6 and 7, respectively. Interestingly, the predicted stable high-pressure structures of B 3 Si and B 4 Si have similar geometries to that of B 2 Si, which are crystalized in the same P3m1 symmetry. For P3m1 B 3 Si and B 4 Si, B atoms also exist in the form of buckled graphenelike layers [65]. The main difference between P3m1 B 2 Si, B 3 Si, and B 4 Si is the number of boron layers between adjacent silicon layers, which increases with increasing the boron content. However, as boron content increases to B 6 Si, where the predicted structures present more than three boron layers between silicon layers (Fig. S5 [59]), it is unstable over the whole studied pressure range with respect to decomposing into B 4 Si and B [ Fig. S4(e) [59]]. In addition, some metastable structures of BSi, B 2 Si, B 3 Si, and B 4 Si are also shown in Fig. S6 [59], which also presents similar structural characteristics to the predicted stable high-pressure structures.
For Si-rich compounds, the predicted I4/m BSi 2 [ Fig. 2(g)] is similar to the well-known CuAl 2 -type arrangement (some semi-transition-metal borides such as Ti 2 B, Ta 2 B, and Mo 2 B also show the I4/m symmetry [66]), in which B atoms form linear chains, where each B is surrounded by eight Si atoms. BSi 3 was calculated to be unstable and dissociates into pure elements even at 400 GPa. The predicted P-1 and P2 1 /m-2 phases of BSi 3 (Fig. S7 [59]) are energetically more stable than the C2/m phase predicted by Wang et al. [29] in the corresponding pressure range.
The electron localization functions (ELFs) were calculated in order to understand the nature of the chemical bonds. As illustrated in Fig. S8 [59 The calculated electronic band structures and density of states (DOS) of stable structures of B-rich compounds are shown in Fig. 3 and Fig. S9 [59]. Several bands crossing the Fermi level (E f ) reveal the metallic character of these structures at 1 atm. Similar to MgB 2 , the B p orbits give main contributions to the electronic states near the E f . Band projections onto B p x , p y , and p z orbitals are also shown in the band structures of P3m1 B 2 Si, B 3 Si, and B 4 Si. In P3m1 B 2 Si [ Fig. 3(a)], the flat bands along the -A direction and around point near the E f are mainly derived from B p x,y , while the steep bands crossing E f derived from B p z . Moreover, along the -A direction the nearly cylindrical Fermi surfaces also appeared, which arise from p x,y B bands (Fig. S10 [59]). These characters in P3m1 B 2 Si are similar to those in MgB 2 [68], which might imply the superconductivity of P3m1 B 2 Si. On the other hand, the flat bands close to the E f corresponding to localized electronic states lead to a large DOS value at E f , which might also favor an enhancement of the electron-phonon interaction. As we mentioned earlier, B atoms exist in the forms of wrinkled single and double layers in  The insets in the right panels of Fig. 3 show the calculated partial charge densities with isosurface of 0.018 of states around the Fermi level. It can be seen that the calculated charge densities are mainly distributed around B atoms. In P3m1 B 2 Si, B-B σ -bond orbitals form conducting channels for valence electrons. Charge densities of states around the Fermi level come from connected electronic channels in the graphenelike B layer. With increasing the number of B layers in B 3 Si and B 4 Si, the electronic channels are gradually destroyed, and charge densities are mainly distributed around the B atoms bonded with Si atoms, which suggests the charge transfer from Si to B atoms. Moreover, the other B atoms that are not bonded with Si atoms lose electrons and act as donors. According to the calculated partial charge densities, the metallicity of P3m1 B 2 Si is stronger than that of P3m1 B 3 Si and B 4 Si, which is consistent with the calculated DOS. Additionally, the electronic properties of the predicted low-pressure structures of B-rich compounds and high-pressure structures of BSi and BSi 2 are shown in Fig. S11 [59]. β-B 3 Si was estimated to be a semiconductor with an indirect band gap of about 1.3 eV, which is consistent with previous study [14]. The high-pressure phases of BSi and BSi 2 are metallic.
As mentioned above, the predicted high-pressure structures of B-rich compounds are metastable at 1 atm, with structural characteristics similar to MgB 2 . Having in mind that MgB 2 is a good superconductor, with a T c of 39 K, we thus performed phonon dispersion, projected phonon DOS, Eliashberg spectral function α 2 F(ω)/ω, and its integral λ(ω) to further explore the superconducting properties of the predicted metallic phases of B 2 Si, B 3 Si, and B 4 Si at 1 atm (Fig. 4). For B 2 Si, the calculated phonon DOS can be separated into two regions. The low-frequency vibrations (0-12.5 THz) are associated with both Si and B atoms and high-frequency modes are only related to B atoms. The calculated logarithmic average frequency ω log is 393 K and the EPC parameter λ is 0.86. The main contribution to the EPC comes from the low- frequency modes, with 78% of the total λ. To further explore the contribution of different phonon modes, red circles with the area proportional to the EPC strength are also plotted in the phonon dispersion curves. All the phonon modes along -A direction have large contributions to the EPC, as well as the soft modes around point. The cases in P3m1 B 3 Si and B 4 Si are similar to that of B 2 Si. However, the calculated EPC parameter λ decreases from 0.86 in P3m1 B 2 Si to 0.68 and 0.41 in P3m1 B 3 Si and B 4 Si, respectively (Fig. 5), which mainly comes from the decreased contribution of the low-frequency modes from 78%, 69% to 56%, respectively.
The T c values of these predicted phases were calculated by using the Allen-Dynes modified McMillan equation [69]. As shown in Fig. 5, the estimated T c of P3m1 B 2 Si, B 3 Si, and B 4 Si at 1 atm is 21-17.5, 9.8-6.8, and 2.8-1.3 K by applying a Coulomb pseudopotential μ * of 0.1-0.13, respectively, where T c decreases with increasing B content. The trends of ω log and λ as a function of the B content are shown in the inset of Fig. 5. The calculated ω log increases with increasing B content, which can be attributed to the presence of more and stronger B-B bonds. However, λ presents an opposite evolution to ω log . In addition, R3m B 2 Si and Pnma BSi are also estimated to superconduct with T c of 7.7 and 2.4 K at 1 atm and 250 GPa, respectively (Fig. S12 [59]). Our results show that most of the predicted high-pressure phases of B-Si compounds are potential superconductors. P3m1 B 2 Si shows the highest T c value of 21 K at 1 atm, which is much higher than the experimentally observed T c of 0.56 K [19] in boron doped diamond-type silicon.
As we know, MgB 2 is a prototypical anisotropic superconductor. Considering the similarity between B 2 Si and MgB 2 , the superconducting anisotropy of P3m1 B 2 Si was also investigated by combining the fully anisotropic Migdal-Eliashberg theory [52,53] with the electron-phonon interpolation based on maximally localized Wannier functions [57]. Figure 6(a) shows the calculated superconducting energy gap nk as a function of temperature with μ * of 0.1. One can clearly see that the P3m1 B 2 Si displays an anisotropic superconductivity with a two-gap nature. From the momentum-resolved superconducting gap nk on the Fermi surface at 10 K [ Fig. 6(b)], we find that the higher nk values are distributed on the two Fermi surface sheets along the -A-line associated to the p x,y orbitals of boron, while the lower nk values are mapped on the remaining Fermi surface sheets. In comparison, the energy difference of 1.2 meV between two average nk values in B 2 Si is much smaller than that of 6.5 meV in MgB 2 at 10 K [56], indicating that B 2 Si is less anisotropic. In MgB 2 , there are not only strong covalent B-B bonds in the graphenelike planes, but also ionic bonds between Mg and B atoms. Two degenerate σ -bond bands associated to B atoms are highly coupled with the E 2g phonon mode, which is mainly associated with the B atomic motion [70]. However, in the predicted B 2 Si, the interplane B-Si bonds are also strongly covalent. The bands near the Fermi level mainly associated to B orbitals are not only coupled with B-related phonon modes, but also with Si-related ones. As temperature increases, the superconducting gap will decrease. The critical temperature T c can be estimated when the superconducting gap value vanishes. As shown in Fig. 6(a), the critical temperature T c for P3m1 B 2 Si is 30 K, which is higher than that of 21 K derived with the Allen-Dynes formula.

IV. CONCLUSIONS
In summary, we have systematically explored the crystal structures and novel properties of B-Si compounds under pressure. At 1 atm, we reproduced the experimentally observed ordered structure β-B 3 Si and also successfully predicted a series of structures containing B 12 units for B 6 Si and B 4 Si, which are similar to the experimental ones. Upon compression, new compositions, B 2 Si, BSi, and BSi 2 , become gradually stable, indicating that pressure promotes the formation of B-Si compounds with high silicon content. On the B-rich side, from B 2 Si to B 6 Si, all the predicted stable high-pressure structures contain buckled graphenelike boron layers and the number of boron layers between adjacent silicon layers increases with increasing the boron content. Phonon calculations show these structures are dynamically stable, so that they might be preserved as metastable phases at 1 atm. Electronic band structure calculations show that all the stable high-pressure phases are metallic. Moreover, electron-phonon coupling calculations indicate that all the retained high-pressure phases of B-rich compounds are superconductors at 1 atm. P3m1 B 2 Si has the highest T c value of 21 K at 1 atm, which is much higher than that of the observed one in B doped diamond-type Si. Within the fully anisotropic Migdal-Eliashberg formalism, our calculations suggest that P3m1 B 2 Si is a two-gap anisotropic superconductor and the estimated T c might reach up to 30 K at 1 atm. Our studies not only enrich the phase diagram of the B-Si system but also may stimulate experimental exploration of superconductivity in B-Si compounds. Moreover, the B-Si system will be helpful for further understanding structural and physical properties of boron-rich compounds under pressure.