Theoretical dissection of superconductivity in two-dimensional honeycomb borophene oxide B2O crystal with a high stability

Atomically thin borophene has recently been synthesized experimentally, significantly enriching the boron chemistry and broadening the family of two-dimensional (2D) materials. Recently, oxides of 2D materials have been widely investigated for next-generation electronic devices. Based on the first-principles calculations, we predict the existence of the superconductivity in honeycomb borophene oxide (B2O), which possesses a high stability and could be potentially prepared by intrinsically incorporating oxygen into the recently synthesized borophene. The mechanical, electronic, phonon properties, as well as electron–phonon coupling of metallic B2O monolayer, have been systematically scrutinized. Within the framework of the Bardeen–Cooper–Schrieffer theory framework, the B2O monolayer exhibits an intrinsic superconducting feature with a superconducting transition temperature (Tc) of ~10.3 K, higher than many 2D borides (0.2–7.8 K). Further, strain can be utilized to tune the superconductivity with the optimal Tc of 14.7 K under a tensile strain of 1%. The superconducting trait mainly originates from the out-of-plane soft-mode vibrations of the system, which are significantly enhanced via the light O atoms’ incorporation compared to other 2D metal-boride superconductors. This strategy would open a door to design 2D superconducting structures via the participation of light elements. We believe our findings greatly bloom the 2D superconducting family and pave the way for future nanoelectronics.

It should be noteworthy that the electron-deficient property of B atoms makes B-B bonds unstable in borophene monolayer, particularly under oxygen-rich conditions 25,26 . This calls for the investigation of the possible 2D "B x O y " materials when exposed to oxygen. By scrutinizing the oxygen adsorption and dissociation on freestanding borophene, O-adsorbed borophene exhibits an enhanced stability via the strong B-O interactions 27 , in line with the characteristics of oxidation as reported in previous experiments 2, 3 . These phenomena are ultimately in favor of the formation of 2D boron oxides under ambient conditions. Besides, oxidation could help to tune the physical and chemical properties of 2D materials as well [28][29][30] . Thus, borophene as electron-deficient monolayer is feasible to be embedded by oxygen via the formation of more stable B-O bonds 31 . It is in this way which the 2D boron oxides could be obtained potentially. Inspired by this, more attentions have been focusing on these 2D boron oxides.  31 . Very recently, a 2D honeycomb boron oxide (B 2 O) with planar monolayer is proposed and exhibits intriguing topological phase transition when subject to external strains 32 . However, superconductivity in 2D boron oxides is yet to be reported experimentally and theoretically [33][34][35][36] . Considering the intrinsic superconductivity in borophene, the successful exploration of superconductivity in 2D boron oxides would not only bloom the 2D boron family but also feasibly facilitate the experimental explorations for realistic applications.
Inspired by prior successes in exploring superconductivity in borophene [16][17][18][37][38][39] , here we systematically study the superconductivity of 2D B 2 O monolayer with the lowest formation energy 31 . Using the first-principles calculations, the thermal and mechanical stability, electron structures, vibration modes, and superconductivity were discussed. The results show that B 2 O is not only a metallic monolayer but also an intrinsic superconductor with superconducting transition temperature (T c ) of 10.3 K that is higher than those of the mostly reported 2D borides. Such superconductivity is attributed to the out-of-plane soft-mode vibrations of O atoms. Moreover, under the tensile strain of 1%, a large softened vibration mode appears along M-X direction, that, leads to an increase of T c by 40%, showing the tunability of superconductivity in B 2 O monolayer.

RESULTS AND DISCUSSION
The structure properties of B 2 O The artificial B 2 O monolayer in this work was verified to be the global minimum structure by adopting particle-swarm optimization 40,41 . B 2 O has a global minimum of energy of −6.80 eV/atom, 0.48 and 0.50 eV smaller than the configurations with the 2nd and 3rd lowest energies, respectively (Supplementary Fig. 1a, b). However, the latter two structures are not dynamically stable, confirmed by the calculated phonon spectra ( Supplementary Fig.  1c, d). Thus, we next only focus our attentions on the one with the global minimum energy. The title B 2 O monolayer crystallizes in the orthorhombic lattice with space group Cmmm (No. 65), possessing a C 2v symmetry (Fig. 1a). The stability of B 2 O crystal The stability of 2D crystals is very important in predicted structures. Here, the cohesive energy (E coh ) and formation energy (E f ) are performed by and  Supplementary Fig. 2. The average value of the free energy remains nearly constant with small fluctuations during the entire simulation period at about 1500 K ( Supplementary Fig. 2a). After 10-ps simulation, we found that there exists a sign of a structural disruption at about 1700 K ( Supplementary Fig. 2b), producing a calculated melting temperature of 1500-1700 K. This melting temperature is higher than the previous report (1000 K) 32 , suggesting a higher thermal stability and thus the potential applications in extreme high-temperature environment.
To further assess the mechanical stability of our structure, we calculated the elastic constants C ij in the rectangle unit cell by where E s is strain energy, and the tensile strain is defined as ε ¼ aÀa0 a0 , and a and a 0 are the lattice constants of the strained and strain-free structures, respectively. ε xx and ε yy are the strains along the x and y directions, and ε xy is the shear strain. The B 2 O belonging to orthorhombic crystal has four independent elastic constants: C 11 , C 12 , C 22 , and C 66 , corresponding to the second partial derivative of strain energy with respect to the applied strain. In order to calculate the elastic stiffness constants, the E s as a function of ε in the range of − 2% ≤ ε ≤ 2% were calculated. The results of the strain energy curves associated with uniaxial, biaxial, and shear strains are shown in Supplementary Fig. 3a. Then, the C ij can be obtained with the aid of the VASPKIT 43 , a post-processing program for the VASP code. Our calculations estimate C 11 , C 12 , C 22   (123 N/m) 47 , phosphorene (24-103 N/m) and silicene (62 N/m) 48 . The Poisson's ratio reflects the mechanical responses of the system against uniaxial strains and can be calculated as ν x ¼ C12 C22 ¼ 0:29 and ν y ¼ C12 C11 ¼ 1:59, indicating the large anisotropy in mechanical properties. To present a full understanding of the mechanical properties of B 2 O monolayer, we calculated the orientation-dependent Y and ν as a function of the polar angle θ (0-360°). For the orthogogonal 2D system, the strain parallel (ε ∥ ) and perpendicular (ε ⊥ ) to the θ direction induced by the unit stress σ(θ)(|σ| = 1) can be expressed as 49,50 and respectively, where Δ ¼ C 11 C 22 À C 2 12 , a ¼ cos(θ) and b ¼ sin(θ). Then, Y(θ) and ν(θ) are derived as and respectively. Clearly, both the variations of Y and ν show a spindlelike shape, indicating the fully anisotropic traits ( Supplementary  Fig. 3b, c).  The B 2 O monolayer is metallic with two bands crossing the Fermi level (Fig. 2a), which is supported by the Fermi surface distributions along the high paths (Fig. 2b). From the orbitalresolved band structure, the bands around the Fermi level are mainly composed of the B-p and O-p orbitals. In Supplementary  Fig. 4, we can clearly see that the orbital hybridization near the Fermi level mainly stems from B-p x,y and B-p z states, followed by some contributions from O-p x,y , O-p z , and B-s states. Thus, the metallic nature of B 2 O monolayer is essentially dominated by the B-p orbitals. In addition, two dirac cones (DCs) exist around the Fermi level. The orbital-resolved band structures of the B 2 O monolayer are presented in Supplementary Fig. 5 to understand the origin of DPs. The DC1 mainly consists of B-p y and O-p z , and the DC2 involves the dominant contributions from B-p z and O-p y . Nevertheless, the s, p x orbitals of two atoms are not responsible for the formation of the DCs (Fig. 2a; Supplementary Fig. 5). The effects of spin-orbit coupling and Heyd-Scuseria-Ernzerhof (HSE06) 52 are evaluated to play a negligible role in the formation of DCs 32 . Even, the DCs are still well maintained with strain reaching up to 3% ( Supplementary Fig. 8b), indicating that the 2D B 2 O is a robust Dirac material. To assess the electronic transport properties of B 2 O monolayer, we further calculated the Fermi velocity (V F ) within PBE level along Y-M by using the equation V F = ∂E _∂k , where the ∂E ∂k is the slop of the linear band structure and the h is the reduced Planck's constant. The calculated V F closing to DC2 are 9.6 × 10 5 and 6.8 × 10 5 m/s, respectively, larger than and comparable with graphene (8.22 × 10 5 m/s) 53 . This high value of V F suggests B 2 O monolayer to be possess a ultrahigh carrier mobility and would facilitate the future electronics. To probe the effects of the percentage of Hartree-Fock (HF) functionals on electronic properties, screened exchange hybrid density functional of HSE06 were carried out as well (Supplementary Fig. 6). Upon increasing the fraction of HF in the HSE06 calculations, the dispersion of band structure shows a small variation when compared with the PBE results, and the main band traits, such as the two DCs, are maintained. This result suggests that HF plays a negligible role in the electronic structure of system. So, we only calculated and analyzed the superconductivity next on the PBE level.

Phonon and superconductivity of B 2 O crystal
The phonon spectrum of B 2 O monolayer shown in Fig. 3 reveals no imaginary phonon modes, indicating that the rhombic phase is kinetically stable, which is in line with previous result 32 . The phonon k ⋅ p theorem is used to sort the phonon branches based on the continuity of the eigenvectors of vibration modes [54][55][56] , j X e Ã k;σ1 ðjÞ Á e kþ4;σ2 ðjÞj ¼ jδ σ1;σ2 À 0ð4Þj; where e Ã k;σ ðjÞ is the displacement of atom j in the eigenvector of (k, σ) vibration mode, and Δ is a small wave vector. As indicated in Fig. 3a, the out-of-plane (ZA), in-plane transverse (TA), and inplane longitudinal (LA) modes constitute the three acoustic branches for B 2 O.
According to the decomposition of the phonon spectrum with respect to the vibration directions of B and O atoms (Fig. 3a), three acoustic branches dominate the low-frequency region (below 400 cm −1 ), where the main contributions are from in-plane and out-of-plan of modes of O atoms. Meanwhile, the two lowest optical branches consisting of the out-of-plan of modes of B-z also make large contributions to this range. Generally, due to the subtle differences in atomic weight, the out-of-plane modes of B-z and O-z mostly occupy the low-frequency region. The midfrequency region from 400 to 900 cm −1 are related entirely to the in-plane of modes of B-xy and the B-z. Moreover, the vibration frequencies larger than 900 cm −1 originate from the out-of-plan of modes of B-xy and O-xy. The highest vibration frequency, 1471 cm −1 , is much larger than that of Mo 2 B 2 (880 cm −1 ) 21 , W 2 B 2 (920 cm −1 ) 22 , Li 2 B 7 (1120 cm −1 ) 20 , AlB 6 (1150 cm −1 ) 23 , β 12 borophene (1200 cm −1 ) 16,57 , χ 3 borophene (1290 cm −1 ) 16,58 , and B 2 C (1365 cm −1 ) 19 , it is even comparable with that of δ 6 borophene (1411 cm −1 ) 42 . Such a high frequency is consistent with the strong covalent bonding, suggested by the former results of ELF (Fig. 1c) and projected phonon density of states (PhDOS) (Fig. 4a).
The results of the PhDOS, Eliashberg spectral function α 2 F(ω), the EPC constant λ, T c , and the derivatives of T c are presented in Fig. 4. The Eliashberg spectral function exhibits that four major peaks are located at 88.3, 119.3, 210.2, and 288.2 cm −1 , respectively, in the low-frequency region. As shown in Fig. 3b, the low-frequency mode phonons contribute mainly to the EPC, accounting for 78.5% of the total EPC (λ = 0.75). The first peak of α 2 F(ω) is mainly responsible for this part,~51.9% of the total EPC. In the frequency range of 83-120 cm −1 , the large values of λ qν along the M-X − Γ directions are visible and lead to the first two largest peaks of the α 2 F(ω) (Figs.  3b and 4b). As a consequence, λ(ω) increases rapidly in this range. In order to probe the underlying causes, we analyze the vibration modes with the largest value of λ qν in this range (see Fig. 3b). Clearly, , β 12 borophene 16,57 , and χ 3 borophene 16,58 . Using the McMillian-Allen-Dynes formula 59 , the frequency-dependent superconducting transition temperature T c (ω) is obtained from Fig. 4c. Its derivative also exhibits the four main peaks in low frequency, which consistent with distributions of α 2 F (ω). The B 2 O monolayer is a medium-coupling superconductor with λ of 0.75 according to the role proposed by Allen et al 59 . and possesses a T c of 10.35 K, which is higher than those of LiC 6 (5.9 K) 60,61 , C 6 CaC 6 (4.0 K) 62,63 , and Cu-BHT (3.0 K) 64 that their values of T c were determined in experiments. Moreover, the T c of B 2 O is also higher than recently reported 2D boride superconductors using a typical value of μ * = 0.1, such as Li 2 B 7 (6.2 K) 20 , B allotrope (6.7 K) 17 , tetr-Mo 2 B 2 (3.9 K), tri-Mo 2 B 2 (0.2 K) 21 , tetr-W 2 B 2 (7.8 K), hex-W 2 B 2 (1.5 K) 22 , rect-GaB 6 (1.7 K), rect-InB 6 (7.8 K), hex-InB 6 (4.8 K) 24 , and AlB 6 (0.95 K) 23 . This increase of T c can be rationalized by the fact that the vibrations of B atoms are significantly enhanced via incorporating light O atoms into the monolayer, in great contrast to the constraining effect associated with heavier metal atoms within 2D metal-boride superconductors. This provides clues for us to design 2D superconducting systems with light elements and opens the road toward further improvement of 2D superconducting feature. It is true that some intrinsic 2D stable B structures such as δ 6 (27.0 K), χ 3 (24.7 K), and β 12 borophene (18.7 K) 16 show superconductivity with higher T c . This is not in contrast to our designing strategy of 2D superconductors via the participation of light elements: The pure 2D boron structures could be regarded the extreme phase of borides by introducing a more light "B" element into 2D sheet in the form of "B x B" when compared with the O's incorporation in B 2 O monolayer. Here, the T c of borophene is also calculated with μ * = 0.1. To explore the effect of the μ * on the T c of B 2 O monolayer, we calculate the T c by varing μ * from 0.08 to 0.15 ( Supplementary Fig. 7). As expected, T c would decrease monotonically decreasing form 11.9 K to 6.8 K upon the increase of μ * .
The oxidation process of black phosphorene could be wellcontrolled by the assistance of Laser 65 and borophene has been successfully synthesized on the substrates 3,4 . Thus, the B 2 O monolayer may be obtained on borophene substrate by oxidation 32 . To simulate the real samples grown on substrates with different lattice constants, we applied the in-plane biaxial strain to the B 2 O monolayer. The biaxial strain can be calculated by ξ ¼ aÀa0 a0 100% (positive value means tensile strain, while negative one indicates compressive strain). In our calculations, the lattice  constants are changed from −1% to 3%, and atomic coordinates are fully relaxed in each case. The band structures and phonon spectra are plotted in Supplementary Fig. 8. By confirmed the phonon spectra in Supplementary Fig. 8a, B 2 O is stable under the tensile strain from 1% to 3%. As indicated in Supplementary Fig.  8b, the tensile strain has few influence on the band structures around the Fermi level, the same as the N(E F ) ( Fig. 5b; Supplementary Fig. 8b). While, the phonon spectra shift to some extent under applied tensile strain. Especially, compared with freestanding sample, appearing a large soften mode along M-X at strain of 1%, is a good phenomena improving the superconductivity. The variations of superconductive parameters [N (E F ), ω log , λ, and T c ] under series of strain are exhibited in Fig. 5. Along with the increasing tensile strain, the ω log first decrease until strain greater than 1%, and then it goes up conversely. While the λ and T c vary in an inverse way, and the N(E F ) varies a little. When tensile strain equals 1%, the λ and T c can be increased to be 1.04 and 14.7 K, respectively, reaching the maximum values in this strain region. However, the bigger tensile strain suppresses the superconductivity. So, strain engineering offers an effective way to tune the superconductivity and facilitates the potential application in future nanodevices.
In summary, within the framework of the density-functional theory (DFT) and Bardeen-Cooper-Schrieffer (BCS) theory, we have systematically investigated the mechanical and electronic properties, phonon vibrations as well as superconductivity of the proposed 2D honeycomb borophene oxide, B 2 O. The B 2 O monolayer exhibits a high stability and possesses two Dirac cones near the Fermi level, and is an intrinsic BCS-type superconductor with a T c of~10.3 K. This T c is higher than mostly reported boride superconductors. The superconducting trait is attributed to the out-of-plane vibration modes of B and O atoms. Upon applying a tensile in-plane strain of 1%, the T c can achieve the maximum value of 14.7 K, which is associated with the large soften mode appearing along M-X direction in ZA branch. Thus, these interesting results would further trigger efforts on 2D superconducting materials.

First-principles calculations
The first-principles calculations based on the DFT were performed through the Vienna ab initio simulation package (VASP) 66,67 and the Quantum-ESPRESSO code. After the full convergence tests, the exchange correlation interaction was simulated within the generalized gradient approximation (GGA) 68,69 with the Perdew-Burke-Ernzerhof (PBE) 70 -type pseudopotential. The electronic wave functions were expanded via the plane wave basis set with a energy cutoff of 600 eV. The Γ-centered 15 × 15 × 1 k-point mesh were adopted using the Monkhorst-Pack method. To avoid the interaction between adjacent monolayers, the vacuum thickness was set to be 15 Å along the z direction. The structure were fully relaxed until the total energy and force on per atom were <10 −5 eV and 0.01 eV/Å, respectively 71,72 . The EPC and superconductivity were calculated by the QE within the densityfunctional perturbation theory (DFPT) 73 and the BCS theory 74 . The optimized norm-conserving Vanderbilt pseudopotentials 75 were used to model the electron-ion interactions. The kinetic energy cutoff and the charge density cutoff of the plane wave basis were chosen to be 80 and 320 Ry, respectively. Self-consistent electron density was evaluated by employing 24 × 24 × 1 k mesh with a Methfessel-Paxton smearing width of 0.02 Ry. The phonon calculations were carried out on the 6 × 6 × 1 q mesh. Meanwhile, the convergences of λ and T c are tested and shown in Supplementary Fig. 9, verifying that the q mesh of 6 × 6 × 1 is large enough to used in our calculations.

Structure screening and AIMD calculations
The particle-swarm optimization (PSO) scheme, as implemented in the CALYPSO code 40,41 , was adopted to search for the global minimum structure for the B 2 O compound. In the PSO calculations, both planar and buckled structures of B 2 O were considered, and the population size and the number of generations were set to be 30. Through the highthroughput calculations, thousands of different structures of B 2 O were generated and ranked by CALYPSO code in order of enthalpy from low to high. The electronic structure calculations were performed through the VASP. We also performed the AIMD simulations at a series of temperatures (300, 500, 700, 900, 1100, 1300, 1500, and 1700 K) with constant number, volume, temperature (NVT) ensemble, lasting for 10 ps with a time step of 1 fs to assess the thermal stabilities of B 2 O monolayer. The 3 × 3 × 1 supercell was adopted and the temperature was controlled using the Nosé-Hoover thermostat 76 .

Superconductivity calculations
Within the BCS and Migdal-Eliashberg theories 77,78 framework, to examine the contribution to λ from individual phonon modes, the magnitude of the EPC λ qν can be calculated by where γ qν , ω qν , and N(E F ) are the phonon linewidths, the frequency of a lattice vibration with crystal momentum q in the branch ν and the density of states (DOS) at the Fermi level, respectively 79 . In addition, the phonon linewidths γ qν can be estimated by γ qν ¼ 2πω qν Ω BZ X k;n;m jg ν kn;kþqm j 2 δðϵ kn À ϵ F Þδðϵ kþqm À ϵ F Þ; where Ω BZ is the volume of Brillouin zone (BZ); ϵ kn and ϵ k+qm are the Kohn-Sham energy, ϵ F is the Fermi energy, and g ν kn;kþqm denotes the EPC matrix element. Moreover, according to the linear response theory 59 the g ν kn;kþqm can be determined self-consistently. Subsequently, the Eliashberg electron-phonon spectral function α 2 F(ω) and the cumulative frequencydependent EPC function λ(ω) can be calculated by and λðωÞ ¼ 2 respectively. The logarithmic average frequency ω log and the superconducting transition temperature T c can be calculated as follows: and T c ¼ ω log 1:2 exp À 1:04 ð1 þ λÞ λ À μ Ã ð1 þ 0:62λÞ ; where μ* = 0.1, a typical value of the effective screened Coulomb repulsion constant 16,57,[80][81][82][83][84] .