Pressure-induced zigzag phosphorus chain and superconductivity in boron monophosphide

We report on the prediction of the zinc-blende structure BP into a novel C2/m phase from 113 to 208 GPa which possesses zigzag phosphorus chain structure, followed by another P42/mnm structure above 208 GPa above using the particle-swarm search method. Strong electron-phonon coupling λ in compressed BP is found, in particular for C2/m phase with the zigzag phosphorus chain, which has the highest λ (0.56–0.61) value among them, leading to its high superconducting critical temperature Tc (9.4 K–11.5 K), which is comparable with the 4.5 K to 13 K value of black phosphorus phase I (orthorhombic, Cmca). This is the first system in the boron phosphides which shows superconductivity from the present theoretical calculations. Our results show that pressure-induced zigzag phosphorus chain in BP exhibit higher superconducting temperature TC, opening a new route to search and design new superconductor materials with zigzag phosphorus chains.

T he crucial thermodynamic parameter-high pressure has been emerging as a powerful tool to investigate physical and chemical behaviours of materials, especially to synthesize or design materials with excellent properties such as superhardness and superconductivity [1][2][3][4][5][6] . Recently, B and P have been studied extensively in physical, chemical, and material science fields due to their interesting structural properties when pressure is applied 4,5,[7][8][9] . Boron has always been recognized as a complex element, both structurally and electronically: its crystalline phases are numerous and inevitably complicated, which is related to its electron deficiency and thus the tendency to form multicenter bonds 5,7 . Phosphorus, on the other hand, has six known phases at room temperature under high pressure [10][11][12][13][14][15] . Black phosphorus phase I (orthorhombic, Cmca) is known to be most stable under ambient temperature and pressure 15 . The experiment at low temperature and high pressure showed that the superconducting temperature of phase I increased from 4.5 K to 13 K under pressure 14 . These peculiar physical properties of compressed B and P solids have motivated our attention on high-pressure study of BP. The binary semiconductor boron phosphide (BP) is the zinc-blende crystal structure (space group F-43m) with lattice parameter a 5 4.537 Å 16 , which has extraordinary properties, such as high hardness, high temperature stability, resistance to chemical corrosion, high thermoelectric powers for direct energy conversion, and is regarded as an important candidate for electronic, optical, and other engineering applications [17][18][19][20][21] . In order to understand in detail the pressure-induced structural behaviour, including mainly the phase transitions of BP, further highpressure studies are essential.
Here, on the basis of comprehensive density functional theory (DFT) computations, we report the design of two new high pressure structures of BP. We performed variable-cell structure prediction simulations using CALYPSO 22,23 approach for BP at 0, 50, 100, 150, 200, 250, and 300 GPa, respectively. For comparison, the structures reported in boron nitride (P-3m1, P6 3 mc, P6 3 /mmc) experimentally and theoretically, are also considered in our calculations. Our calculations demonstrate that, F-43m is the most stable phase at ambient conditions, two new high pressure phases C2/m and P4 2 /mnm can be predicted during compression up to 300 GPa. We further elucidate the energetic, mechanical, electronic and superconducting properties of the obtained novel phases, confirming that F-43m is semiconductor, and discovering two high-pressure superconducting phases C2/m and P4 2 /mnm.

Results
The analysis of the predicted structures gives us a list of candidate structures with space groups F-43m, C2/m, P4 2 /mnm, P6 3 mc, Cmmm, Immm, P63/mmc-I, P63/mmc-II, and the previously reported high pressure rock salt structure Fm-3m 20 are depicted in Figure S1 (Supporting Information). The relative enthalpies per chemical formula unit vs. pressure curves of selected structures are plotted in Figure 1(a). Considering that F-43m and P6 3 mc structures hold very close enthalpy, the enthalpy difference vs. pressure curves for the two phases are specially inserted in Figure 1(a), and it indicates that F-43m has lower enthalpy, which is consistent with the experimental results. In order to investigate the phase transition pressures clearly, the enthalpies vs. pressure curves for F-43m, Immm, P4 2 /mnm, C2/m, and P-3m1 are given in Figure 1(b). From Figure 1, one can see that, for compressed BP, the zinc-blende crystal structure F-43m is the most stable structure below 113 GPa, a C2/m phase takes over the pressure range from 113 to 208 GPa, followed by another P4 2 /mnm structure above 208 GPa. The pressure evolution of the unit cell volume of BP in the structures F-43m, C2/m, and P4 2 /mnm is inserted in Figure 1(b). The experimental equation of state data 24 for single-crystal BP with F-43m structure up to 55 GPa were listed for comparison. It can be seen that our calculated data are in good agreement with the experimental results. The zinc-blende crystal structure F-43m is stable on compression to 113 GPa, and then complete change to C2/m is observed. Comparisons with the F-43m cell at 113 GPa shows a 14.5% reduction in volume. Upon further compression into the region of P4 2 /mnm above 208 GPa, the volume difference peaks at only 1.2%.
The optimized structural parameters of zinc-blende crystal structure F-43m are listed in Table S1. For this phase, the equilibrium lattice constants are 4.547 Å and 4.492 Å , for PBE and LDA calculation, respectively. The present results are in good agreement with the previous results 21, [25][26][27] . In this structure, four boron atoms lie in the Wyckoff 4a site and four phosphorus atoms occupy the 4c site, in which three-dimensional (3D) boron and phosphorus network in this structure is formed, as shown in Figure S1. The optimized structural parameters of another two phases C2/m and P4 2 /mnm at high pressure are presented in Table 1. For monoclinic C2/m phase, the equilibrium lattice constants at 120 GPa are a 5 4.275 Å , b 5    3.681 Å , c 5 8.357 Å and b 5 122.627 degree. The structure of C2/m is shown in Figure 2(a) and Figure 2(b), it can be seen that the 1Dzigzag phosphorus chain is formed, which is connected by boron atoms. For hexagonal P4 2 /mnm phase (Figure 2(c)), the equilibrium lattice constants at 210 GPa are a 5 3.631 Å , c 5 3.640 Å . The B-P, B-B, and P-P bonds are formed in this structure, but P-P bonds are not continuous (Figure 2(d)).
The calculated elastic constants for three BP phases with low enthalpy are presented in Table 2. From Table 2, the elastic constants reveal that all BP phases satisfy the mechanical stable criterion at the corresponding pressures. Furthermore, the hexagonal phase P4 2 / mnm can remain stable at ambient pressure, but monoclinic C2/m does not match the mechanical stable criterion at ambient pressure, which only keeps stable at the high pressure range. Also, the calcu-    Figure 3. The maximum optical branch frequencies are 803.5 cm 21 for F-43m phase at atmospheric pressure, 983.7 cm 21 for C2/m phase at 120 GPa, and 1175.2 cm 21 for P4 2 /mnm phase at 210 GPa. One can see that the maximum optical branch frequencies increase with increasing pressure due to an obvious structural difference. From PPHDOS, the lower frequency region is associated with phosphorus atoms, whereas boron atom contributes to the high-frequency region for high-pressure phases (C2/m and P4 2 /mnm) at their corresponding pressure regions. But in the F-43m phase, the boron atom dominates the high-frequency region, whereas phosphorus atom contributes almost equally to the high and low frequency region.
To investigate the possibility of pressure induced metallization and superconducting transition, band structures for F-43m, C2/m and P4 2 /mnm phases have been calculated to examine the electronic properties of these polymorphs of BP at their pressure ranges as shown in Figure 3. Electronic structure calculations indicate that F-43m phase is a semiconductor, whereas C2/m and P4 2 /mnm phases are metallic. From projected density of states (PDOS) plotted in Figure 3, it can be observed that all phases of BP have consistent electronic distributions. B-2p and P-3p electrons dominate in the wide energy range and form the strong covalent bonds as suggested by the matching B-2p and P-3p curve shapes. The contour plots of charge density on the chosen planes are also shown in Figure 4. The strong covalent B-P bonding within in the 3D B-P network is evidenced (Figures 4(a), 4(b) and 4(c)). Moreover, according to Mulliken population analysis, P-P bonding does not exist in F-43m phase, but it exists in C2/m and P4 2 /mnm phases. Interestingly, P-P zigzag chains in phase C2/m transfer to discontinuous bonding in phase P4 2 /mnm. This conclusion is further substantiated by the electron localization function (ELF) shown in Figures 4(d), 4(e) and 4(f), where the calculated ELF patterns have isosurface values of 0.8, 0.65, and 0.7 for F-43m, C2/m, and P4 2 /mnm, respectively.
The high-pressure C2/m phase of BP exhibits interesting features in its electronic band structure. As shown in Figure 3(e), the electronic bands of the C2/m phase crossing the Fermi level along  [28][29][30] . Similar electronic structures have previously been found in niobium that has the highest transition temperature T C of 9.3 K 31 in all element superconductors at ambient pressure, 9.8 K in CaC 2 28 and 220-235 K in CaH 6 at high pressure of 150 GPa 32 . This motivates us to further investigate the superconductivity of BP at high pressures. The calculated spectral function a 2 F(v) of BP phases C2/m and P4 2 /mnm at 120 GPa and 210 GPa were plotted in Figures 3(b) and 3(c). At 120 GPa (Figure 3(b)), the coupling parameter l is 0.61 with the average phonon frequency v ln of 593 K. Using the strong coupling Allen-Dynes equation, an extension of the McMillan theory, with a nominal Coulomb pseudopotential parameter (m*) of 0.12 the estimated superconducting critical temperature T C is 11.5 K. With increasing pressure, the calculated T C become slightly lower as 11.1 K at 150 GPa and 9.4 K at 200 GPa due to slightly smaller l of 0.59 at 150 GPa and 0.56 at 200 GPa, respectively. Here the T C of two high-pressure phases for BP and T C 's dependence on pressure is shown in Figure 5. Calculated coupling parameter l values and logarithmic phonon momentum v log vs. pressure curves are also presented in Figure 5. It can be seen that logarithmic phonon momentum v log increases with pressure, and coupling parameter l decreases with pressure. Among these two high pressure phases, C2/ m phase has the strongest electron-phonon coupling and so has the highest T C . Apparently, the fairly high T C is due to the large electron phonon coupling (l). The origin can be traced from comparing the calculated Eliashberg spectral function (a 2 F(v)/v) with the projected phonon DOS. As shown in Figure 3(b), nearly 80% of the electron phonon coupling is contributed by the low-frequency vibrations in the frequency region from 100 to 600 cm 21 . The low-frequency vibrations in the frequency region from 100 to 600 cm 21 , which are mostly attributed to the P atom. At wide range of pressure, C2/ m phase has comparative superconducting critical temperatures (from 11.5 K at 120 GPa to 9.3 K at 200 GPa) with 11.5 K value of CaC 6 33,34 . However, T C of BP is close to zero above 210 GPa due to the smaller electron-phonon coupling parameter (l , 0.4). We thus further analyze the coupling parameter in each q-point (Figures 3(b) and 3(c)). It is clearly seen that the coupling parameters in C-Z and X-C are very strong. This superconducting feature is typical anisotropy, leading a low superconductivity value.
To summarize, using the PSO technique on crystal structure prediction, we designed two new high-pressure phases of BP, C2/m and P4 2 /mnm, which are stable in the pressure ranges of 113-208 GPa and 208-300 GPa, respectively. Elastic constants and phonon calculations have shown their mechanical and dynamical stability at the dominating pressure ranges. Our calculations reveal that the phosphorus atomic arrangement form from single atoms to 1D zigzag chains to discontinuous bonding. The high superconducting transition temperature of the C2/m phase benefits from the simultaneous presence of the steep bands (highly mobile electrons) and extremely flat bands (highly confined electrons), which is known to favor the electron paring and superconducting behavior.

Computational Methods
The search of high-pressure structures was performed with variable-cell PSO 35 structure prediction simulations using CALYPSO approach 22,23 in combined with Vienna ab initio simulation package (VASP) 36 . CALYPSO was designed to predict stable or metastable crystal structures requiring only chemical compositions of a given compound at specified external conditions. All structures were locally optimized using the DFT method. The 60% structures of each generation with lower enthalpies were selected to generate the structures for the next generation by PSO operation, and the other structures in new generation were randomly generated to increase the structural diversity. The underlying ab initio calculations were performed using density functional theory within the generalized gradient approximation (GGA) 37 , as implemented in the Vienna ab initio simulation package (VASP). The total energy calculations were carried out with the VASP code. The all-electron projector augmented wave (PAW) method 38 was employed with a plane-wave cutoff energy of 600 eV for all phases. The k-point samplings in the Brillouin zone were performed using the Monkhorst-Pack scheme 39 . The total energy convergence tests showed that convergence to within 1 meV/atom was achieved with the above calculation parameters. Single crystal elastic constants were calculated via a strain-stress approach. i.e., by applying a small strain to the equilibrium lattice of the unit cell and fitting the dependence of the resulting stress on the strain. The bulk modulus, shear modulus, Young's modulus, and Poisson's ratio were determined by using the Voigt-Reuss-Hill approximation 40 . The latticedynamical and superconducting properties are calculated using the density functional perturbation theory (DFPT) 41 as implemented in the Quantum Espresso package 42 using with the Troullier-Martins pseudopotentials with cutoff energies of 50 and 500 Ry for the wave functions and the charge density, respectively. Fine k-points mesh of MP was used and the estimated energy error in self-consistency was less than 10 24 a.u. The electron-phonon coupling was convergent with a finer grid and a Gaussian smearing of 0.02 Ry. In order to check the results, the phonon calculations were also carried out by using a supercell approach as implemented in the PHONOPY code 43 .