Superhard and Superconducting Bilayer Borophene

Two-dimensional superconductors, especially the covalent metals such as borophene, have received significant attention due to their new fundamental physics, as well as potential applications. Furthermore, the bilayer borophene has recently ignited interest due to its high stability and versatile properties. Here, the mechanical and superconducting properties of bilayer-δ6 borophene are explored by means of first-principles computations and anisotropic Migdal–Eliashberg analytics. We find that the coexistence of strong covalent bonds and delocalized metallic bonds endows this structure with remarkable mechanical properties (maximum 2D-Young’s modulus of ~570 N/m) and superconductivity with a critical temperature of ~20 K. Moreover, the superconducting critical temperature of this structure can be further boosted to ~46 K by applied strain, which is the highest value known among all borophenes or two-dimensional elemental materials.


Introduction
Superconductivity in two-dimensional (2D) materials has attracted perennial and ever-increasing attention over past decades, owing both to fundamental scientific interest and tantalizing applications [1][2][3][4][5].One of the most important goals in superconductivity research is raising the superconducting transition temperature, T c .In general, superconducting materials can be classified into two categories: low-temperature superconductors, characterized by critical temperatures (T c ) below 30 K, such as Nb-Ge films [6], and hightemperature superconductors with T c exceeding 30 K, exemplified by Cu-based oxides [7][8][9].Within the framework of the conventional Bardeen-Cooper-Schrieffer (BCS) theory [10], it is reasonably anticipated that metals composed of light elements have a better chance to induce a high T c , because the Debye temperatures within such metals are usually high enough to trigger a strong phonon-mediated superconducting pairing.More specifically, according to the celebrated McMillan-Allen-Dynes (MAD) formula [11,12]: T c should be elevated by increasing the log-averaged characteristic phonon frequency (ω log ) and the electron-phonon coupling (EPC) parameter, λ.The light elemental materials typically have high frequency phonon modes, enlarging ω log and thus increasing T c .Furthermore, large phonon frequency and EPC potential V = λ N(E f ) induced by strong covalent bonding in light elemental material, together with N E f , and the electronic density of states (DOS) at the Fermi level [13], should result in a higher T c .In fact, the metals with strong covalent bonding can be grouped as "covalent metals" [13] and their potential to harbor a high T c in 2D form have been confirmed, e.g., hydrogenated monolayer MgB 2 (67 K) [14], doped graphene (8.1 K~30 K) [15][16][17], hydrogenated monolayer borophene (32.4 K) [5] and especially 2D metal borides (1.4 K~72 K) [18][19][20][21][22].Among covalent metals, 2D boron sheets (borophenes) recently came to the fore of superconducting research, motivated by their inherent metallicity, light weight and significant experimental progress on their synthesis [23][24][25].Based on the first-principles calculations and MAD formula, Penev et al. [26] reported early on the T c of the (experimentally synthesized δ 6 , β 12 and χ 3 borophenes) to be in the range of 11.5 K~20.5 K and Gao et al. [27] found them to span 18.7 K-24.7 K. Using a more accurate and sophisticated anisotropic Migdal-Eliashberg (ME) equation, Zhao et al. [28] found the T c in those borophenes in 26.2 K-33 K. Apart from intrinsic metallicity and light weight, another impetus for searching high-temperature superconductivity in borophenes is their vast polymorphs and superior mechanical properties, which provide a great benefit to potential superconductors, suggesting effective ways to modulate the superconductivity [29].
The monolayer borophenes have been intensively studied; however, the investigations are scarce [30][31][32] for bilayer or multilayer borophenes, where more intriguing properties and tunability could be manifested, compared to their monolayer counterparts.Very recently, the realization of bilayer borophenes were reported [33,34].Unfortunately, the first reported bilayer α-borophene is expected to be unstable if peeled off from the metal substrate [35], and the β 12 -like bilayer synthesized on Cu(111) surface appears complicated due to ambiguity (with more than three hundred atoms in its unit cell and the atomic structure still in debate), preventing an exploration of its properties or further application.It can be anticipated that more stable and versatile bilayer borophenes are still to be unveiled for research.
In this work, we investigate the mechanical and superconducting properties of a bilayer composed of two covalently bonded δ 6 borophene monolayers, "(BL)-δ 6 " in Figure 1a-c.The compact atomic structure endows BL-δ 6 with remarkable stability and tantalizing properties compared to its monolayer counterparts.Its computed high Young's modulus of 570 N/m is even higher than 346 N/m for graphene (or near 280 N/m "per layer", for a fair comparison).Based on the anisotropic ME equation, we find its T c = 20 K and can be boosted to 46 K by strain effect, reaching the highest T c found among all borophenes or elemental 2D materials known to date.The coexistence of superhard (in basal plane direction) and superconducting properties in one material is rare, making BL-δ 6 a fascinating superconductor for emergent nanoscale devices such as quantum interferometers, superconducting transistors, superconducting qubits, and wear-resistant parts of superconducting devices [3,36,37].metals with strong covalent bonding can be grouped as "covalent metals" [13] and their potential to harbor a high Tc in 2D form have been confirmed, e.g., hydrogenated monolayer MgB2 (67 K) [14], doped graphene (8.1 K~30 K) [15][16][17], hydrogenated monolayer borophene (32.4 K) [5] and especially 2D metal borides (1.4 K~72 K) [18][19][20][21][22].Among covalent metals, 2D boron sheets (borophenes) recently came to the fore of superconducting research, motivated by their inherent metallicity, light weight and significant experimental progress on their synthesis [23][24][25].Based on the first-principles calculations and MAD formula, Penev et al. [26] reported early on the Tc of the (experimentally synthesized δ , β and χ borophenes) to be in the range of 11.5 K~20.5 K and Gao et al. [27] found them to span 18.7 K-24.7 K. Using a more accurate and sophisticated anisotropic Migdal-Eliashberg (ME) equation, Zhao et al. [28] found the Tc in those borophenes in 26.2 K-33 K. Apart from intrinsic metallicity and light weight, another impetus for searching high-temperature superconductivity in borophenes is their vast polymorphs and superior mechanical properties, which provide a great benefit to potential superconductors, suggesting effective ways to modulate the superconductivity [29].
The monolayer borophenes have been intensively studied; however, the investigations are scarce [30][31][32] for bilayer or multilayer borophenes, where more intriguing properties and tunability could be manifested, compared to their monolayer counterparts.Very recently, the realization of bilayer borophenes were reported [33,34].Unfortunately, the first reported bilayer α-borophene is expected to be unstable if peeled off from the metal substrate [35], and the β -like bilayer synthesized on Cu(111) surface appears complicated due to ambiguity (with more than three hundred atoms in its unit cell and the atomic structure still in debate), preventing an exploration of its properties or further application.It can be anticipated that more stable and versatile bilayer borophenes are still to be unveiled for research.
In this work, we investigate the mechanical and superconducting properties of a bilayer composed of two covalently bonded δ borophene monolayers, "(BL)-δ " in Figure 1a-c.The compact atomic structure endows BL-δ with remarkable stability and tantalizing properties compared to its monolayer counterparts.Its computed high Young's modulus of 570 N/m is even higher than 346 N/m for graphene (or near 280 N/m "per layer", for a fair comparison).Based on the anisotropic ME equation, we find its Tc = 20 K and can be boosted to 46 K by strain effect, reaching the highest Tc found among all borophenes or elemental 2D materials known to date.The coexistence of superhard (in basal plane direction) and superconducting properties in one material is rare, making BL-δ a fascinating superconductor for emergent nanoscale devices such as quantum interferometers, superconducting transistors, superconducting qubits, and wear-resistant parts of superconducting devices [3,36,37].,c) side views of the atomic structure of BL-δ 6 .The yellow shaded region is its unit cell, and the green shaded region is the unit cell of δ 6 .There are three irreducible boron atoms colored as B1 (green), B2 (blue), and B3 (pink).The direction dependence of (d) Young's modulus and (e) Poisson's ratio of BL-δ 6 .(f) Stress-strain curves of BL-δ 6 , the vertical blue and red dashed lines denote the fractured stress along the a and b directions, respectively.1a), and BL-δ 6 can be viewed as AB stacking of two δ 6 monolayers bonded via interlayer covalent bonds (side view in Figure 1b,c).The unit cell of BL-δ 6 is a rectangle with lattice constants a = 3.243 Å, b = 2.883 Å (Table 1), which is obtained with a Vienna ab initio Simulation Package [38,39] (computational details can be found in Supplementary Note S1).The total energy of BL-δ 6 is computed to be lower than all experimental synthesized monolayers, due to the interlayer bonding: −6.38 eV/atom, which is 0.171 eV/atom, 0.125 eV/atom and 0.113 eV/atom lower than that of the synthesized δ 6 , β 12 and χ 3 , respectively (see Table 1, for atomic structures see Figure S2).Moreover, the slight dynamical instability of δ 6 is also eliminated by the interlayer bonding of BL-δ 6 .The energetic stability of BL-δ 6 was also confirmed by an ab initio evolutionary global structure search in our prior work [40] and reported by Zhou et al. [41].Given that the monolayer δ 6 has been fabricated on a Ag(111) substrate, the thermal stability of BL-δ 6 is also tested on Ag(111) substrate by an AIMD simulation at 500 K, and one can observe that its whole structure is well maintained after 5 ps with a timestep of 1 fs under such a high temperature (Figure S3).Therefore, it is well-expected that BL-δ 6 will be accessed experimentally, supported by its outstanding energetic and thermal stability, and the progress in bilayer borophene synthesis [42,43].We obtain the four independent elastic constants for BL-δ 6 : C 11 = 570 N/m, C 22 = 322 N/m, C 12 = 22 N/m, and C 44 = 146 N/m, which apparently satisfy the Born-Huang criteria [46]: 12 > 0 and C 44 > 0, demonstrating the mechanical stability of BL-δ 6 .The θ dependence of in-plane Young's modulus Y and Poisson's ratio ν of BL-δ 6 are plotted with polar coordinates in Figure 1d,e (computational details can be found in Supplementary Note S2 and Figures S5 and S6).We find that the in-plane Young's moduli and the Poisson's ratios of BL-δ 6 are highly anisotropic.For comparison, the values of elastic constants, Young's moduli, Poisson's ratios of δ 6 , β 12 , χ 3 and graphene are summarized in Table 1.The introduced B3-B3 bonds in BL-δ 6 boost the Y a and Y b to 570 N/m and 321 N/m, respectively.Although graphene is well-known as the strongest 2D material, along the armchair direction Y a = 570 N/m exceeds the Y a = 364 N/m of graphene, partially due to the strong directional B-B σ-bonds, making it stand out among 2D materials.In contrast, the Y b = 321 N/m (still comparable to graphene) is lower, due to the much weaker multicenter bonds involved in the zigzag direction.
We find the fracture strain is 13% along the a and 8% along the b direction, according to the strained phonon spectra (Figure S4), before reaching the elastic maximum of 16% and 14% (Figure 1f).Thus, the failure mechanism of BL-δ 6 is phonon instability, for both directions, which is different from that of δ 6 (elastic instability in the zigzag direction and phonon instability in the armchair direction) [47].The fracture strengths of BL-δ 6 are 43 and 23 N/m along the a and b directions, much above that of δ 6 (20.26 and 12.98 N/m) [47], black phosphorene (9.99 and 4.44 N/m) [48], and even graphene (40.41 and 36.74N/m) [49] along the zigzag direction.Considering the ultrahigh Young's moduli and large fracture strength, one can judge the BL-δ 6 being in-plane superhard.

Electronic Properties
The electronic structure of BL-δ 6 is summarized in Figure 2. Three bands cross the Fermi level and form the Fermi surface, of which one forms a small electron pocket centered around the corner of the first BZ and the other two contribute to the rest of Fermi surface (Figure 2b).The atomic-and orbital-resolved band structures suggest that the main contribution near the Fermi level is from the p y orbital of B1/B3 atoms and the p z orbital of all B atoms (Figure 2f,g, more details in Figure S7), which is also confirmed by the charge distribution in the states within ±0.3 eV from the Fermi level (Figure 2c-e).Using Wan-nier90 (v 3.1.0)code [50], the Wannier-interpolated band structures are also provided and show excellent agreement with those by first-principles calculations (Figure 2a), which lays a clear ground for subsequent anisotropic superconductivity calculations, as implemented in electron-phonon Wannier (EPW v 5.4) code [51].phonon instability in the armchair direction) [47].The fracture strengths of BL-δ are 43 and 23 N/m along the a and b directions, much above that of δ (20.26 and 12.98 N/m) [47], black phosphorene (9.99 and 4.44 N/m) [48], and even graphene (40.41 and 36.74N/m) [49] along the zigzag direction.Considering the ultrahigh Young's moduli and large fracture strength, one can judge the BL-δ being in-plane superhard.

Electronic Properties
The electronic structure of BL-δ is summarized in Figure 2. Three bands cross the Fermi level and form the Fermi surface, of which one forms a small electron pocket centered around the corner of the first BZ and the other two contribute to the rest of Fermi surface (Figure 2b).The atomic-and orbital-resolved band structures suggest that the main contribution near the Fermi level is from the  orbital of B1/B3 atoms and the pz orbital of all B atoms (Figure 2f,g, more details in Figure S7), which is also confirmed by the charge distribution in the states within ±0.3 eV from the Fermi level (Figure 2c-e).Using Wannier90 (v 3.1.0)code [50], the Wannier-interpolated band structures are also provided and show excellent agreement with those by first-principles calculations (Figure 2a), which lays a clear ground for subsequent anisotropic superconductivity calculations, as implemented in electron-phonon Wannier (EPW v 5.4) code [51].

Isotropic and Anisotropic Superconducting Properties
In the phonon spectrum of BL-δ (Figure 3a), we find no negative frequencies, demonstrating its dynamic stability.This is in stark contrast to δ , a layer whose phonon dispersion displays small imaginary frequency near the Γ point [24].Mechanically intuitive, the BL-δ is stabilized by the covalent B3-B3 bonds.The () mode-and (q) momentum-resolved EPC  is also given, calculated by the following equation:

Isotropic and Anisotropic Superconducting Properties
In the phonon spectrum of BL-δ 6 (Figure 3a), we find no negative frequencies, demonstrating its dynamic stability.This is in stark contrast to δ 6 , a layer whose phonon dispersion displays small imaginary frequency near the Γ point [24].Mechanically intuitive, the BL-δ 6 is stabilized by the covalent B3-B3 bonds.The (ν) mode-and (q) momentum-resolved EPC λ qν is also given, calculated by the following equation: where N k represents the total number of k points in the k space, g nm k,qν is an EPC matrix element, n/m and ν represent the indices of electronic bands and phonon mode, are the eigenvalues with respect to the Fermi level, and ω qν is the phonon frequency.We found three modes (B 2u and two A g ) have comparatively large EPC near the Γ point (left panel of Figure 3a).In the B 2u shear mode at 39.3 meV, the B1 and B2 atoms move in one direction and the B3 atoms move in the opposite (Figure 3b).The first A g mode at 102.6 meV corresponds to a bond stretching contributed by B3 atoms (Figure 3c).The second A g mode at 168.7 meV is also a stretch mode, but predominantly contributed by B1 atoms (Figure 3d).The vibrational contribution of these three modes can also be inferred from the phonon density of states (right panel of Figure 3a).The Eliashberg spectral function α 2 F(ω) is a key quantity, which determines the total EPC λ by the following equation: where Nk represents the total number of k points in the k space,  , is an EPC matrix element, n/m and  represent the indices of electronic bands and phonon mode,  and  are the eigenvalues with respect to the Fermi level, and  is the phonon frequency.We found three modes (B2u and two Ag) have comparatively large EPC near the Γ point (left panel of Figure 3a).In the B2u shear mode at 39.3 meV, the B1 and B2 atoms move in one direction and the B3 atoms move in the opposite (Figure 3b).The first Ag mode at 102.6 meV corresponds to a bond stretching contributed by B3 atoms (Figure 3c).The second Ag mode at 168.7 meV is also a stretch mode, but predominantly contributed by B1 atoms (Figure 3d).The vibrational contribution of these three modes can also be inferred from the phonon density of states (right panel of Figure 3a).The Eliashberg spectral function  () is a key quantity, which determines the total EPC  by the following equation: Accordingly, the  () and cumulative EPC () are calculated (middle panel of Figure 3a).We observe three sharp peaks mainly contributed by the three phonon modes analyzed above.Correspondingly, the EPC strength λ calculated by Equation ( 3) is 0.59 and the Tc of BL-δ within the MAD approximation calculated by Equation ( 1) is 10.1 K. Akin to structure and mechanics anisotropy, the Fermi surface of BL-δ , formed by multiple bands with different orbital contributions (py and pz orbitals, Figure 2f,g) is also significantly anisotropic.For such a system with a complicated Fermi surface, using the anisotropic ME equation is essential, in order to obtain an anisotropic EPC, and an accurate Tc, compared with the isotropic superconductivity calculated from the MAD formula using an isotropic EPC [28].As shown in Figure 4a,b, the variation in momentum-dependent EPC parameter  and the superconducting gap Δ at 10 K displays similar Accordingly, the α 2 F(ω) and cumulative EPC λ(ω) are calculated (middle panel of Figure 3a).We observe three sharp peaks mainly contributed by the three phonon modes analyzed above.Correspondingly, the EPC strength λ calculated by Equation ( 3) is 0.59 and the T c of BL-δ 6 within the MAD approximation calculated by Equation (1) is 10.1 K.
Akin to structure and mechanics anisotropy, the Fermi surface of BL-δ 6 , formed by multiple bands with different orbital contributions (p y and p z orbitals, Figure 2f,g) is also significantly anisotropic.For such a system with a complicated Fermi surface, using the anisotropic ME equation is essential, in order to obtain an anisotropic EPC, and an accurate Tc, compared with the isotropic superconductivity calculated from the MAD formula using an isotropic EPC [28].As shown in Figure 4a,b, the variation in momentum-dependent EPC parameter λ k and the superconducting gap ∆ k at 10 K displays similar anisotropy, i.e., the maximum along the Γ − X direction and the minimum along the Γ − Y direction.We find the superconducting gap ratio ∆ aniso = (∆ max − ∆ min )/∆ ave = (2.499− 0.444)/1.470≈ 140%, a measure of its strong anisotropy at the Fermi surface, also signifying that the anisotropic ME formula is indispensable for predicting T c . Figure 4c shows the evolution of the superconducting gap as a function of temperature, based on the ME equations solved in either isotropic or fully anisotropic approximations.The T c is identified as the lowest temperature at which the vanishing gap is observed.Under the isotropic approximation T c = 12 K, comparable to the value (10 K) obtained by the MAD formula but is much lower than the value calculated with the fully anisotropic approximation (20 K), further attesting to the anisotropic superconducting nature of BL-δ 6 , omitted previously [32].
Materials 2024, 17, x FOR PEER REVIEW 6 of 11 anisotropy, i.e., the maximum along the Γ − X direction and the minimum along the Γ − Y direction.We find the superconducting gap ratio Δ = (Δ − Δ )/Δ = (2.499− 0.444)/1.470≈ 140%, a measure of its strong anisotropy at the Fermi surface, also signifying that the anisotropic ME formula is indispensable for predicting Tc. Figure 4c shows the evolution of the superconducting gap as a function of temperature, based on the ME equations solved in either isotropic or fully anisotropic approximations.The Tc is identified as the lowest temperature at which the vanishing gap is observed.Under the isotropic approximation Tc = 12 K, comparable to the value (10 K) obtained by the MAD formula but is much lower than the value calculated with the fully anisotropic approximation (20 K), further attesting to the anisotropic superconducting nature of BL-δ , omitted previously [32].The mechanical robustness of BL-δ permits a consideration of whether the tensile strain could enhance the Tc, motivated mainly by two reasons.First, according to Equation (2), it can be seen that EPC  is inversely proportional to the phonon frequency  , which would be helpful to lower, perhaps by tension, weakening the atomic force constants and thus softening the phonon modes.Second, the atomic orbitals overlap is reduced so that the electronic bands become less dispersive, enlarging the ( ), which provides more electrons susceptible to pairing interactions, mediated by dynamic phonons, to certainly contribute to the Tc rise.
To examine the effect of strain on the BL-δ superconductivity, the evolution of EPC  , the log and maximum frequencies  and  , as well as Tc, are all calculated within three approximations (i.e., MAD, isotropic ME, anisotropic ME), as a function of tensile strain along the a direction, and presented in Figure 5a,b (for details see Figures S8-S11).As expected, the increment in Tc is accompanied by a decrease in frequency-dependent  / and an increase in EPC  (Figure 5b).The phonon spectrum under a 13% tension along a (near fracture strain) weighted by  , the Eliashberg spectral function  () and EPC () are shown in Figure 5c,d.The red shift of phonon frequency The mechanical robustness of BL-δ 6 permits a consideration of whether the tensile strain could enhance the T c , motivated mainly by two reasons.First, according to Equation (2), it can be seen that EPC λ qν is inversely proportional to the phonon frequency ω qν , which would be helpful to lower, perhaps by tension, weakening the atomic force constants and thus softening the phonon modes.Second, the atomic orbitals overlap is reduced so that the electronic bands become less dispersive, enlarging the N(E f ), which provides more electrons susceptible to pairing interactions, mediated by dynamic phonons, to certainly contribute to the T c rise.
To examine the effect of strain on the BL-δ 6 superconductivity, the evolution of EPC λ, the log and maximum frequencies ω log and ω max , as well as T c , are all calculated within three approximations (i.e., MAD, isotropic ME, anisotropic ME), as a function of tensile strain along the a direction, and presented in Figure 5a,b (for details see Figures S8-S11).As expected, the increment in T c is accompanied by a decrease in frequency-dependent ω log /ω max and an increase in EPC λ (Figure 5b).The phonon spectrum under a 13% tension along a (near fracture strain) weighted by λ qν , the Eliashberg spectral function α 2 F(ω) and EPC λ(ω) are shown in Figure 5c,d.The red shift of phonon frequency can be clearly observed and the enhanced EPC λ qν due to the phonon softening can be seen in all frequency ranges (Figure 5d).In addition, the contribution of low-frequency part becomes more dominant.In particular, four largely softened phonon Kohn anomalies with considerable EPCs around 15 meV appear on the path of X-M-Y.Eventually, these combined effects lead to a significant increase, by strain, in λ from 0.59 to 1.33 and T c from 20 K to 46 K.It should be emphasized that 46 K is the highest T c among all borophenes and elemental 2D materials (Table 2).It should be noted that for bilayer or multilayer materials with positive Poisson ratios, such as BL-δ 6 , tensile strain would lead to vertical shrinking (Figure 5a), suggesting that the vertical pressure may also be an effective way to enhance the T c of BL-δ 6 , as the tensile strain does.can be clearly observed and the enhanced EPC  due to the phonon softening can be seen in all frequency ranges (Figure 5d).In addition, the contribution of low-frequency part becomes more dominant.In particular, four largely softened phonon Kohn anomalies with considerable EPCs around 15 meV appear on the path of X-M-Y.Eventually, these combined effects lead to a significant increase, by strain, in  from 0.59 to 1.33 and Tc from 20 K to 46 K.It should be emphasized that 46 K is the highest Tc among all borophenes and elemental 2D materials (Table 2).It should be noted that for bilayer or multilayer materials with positive Poisson ratios, such as BL-δ , the tensile strain would lead to vertical shrinking (Figure 5a), suggesting that the vertical pressure may also be an effective way to enhance the Tc of BL-δ , as the tensile strain does.[26] 1.10 20.5   [28] 0.82 27.0   [26] 0.80 16.1   [27] 0.89 18.7   [28] 1.01 33.0   [26] 0.60 11.5   [27] 0.95 24.7   [28] 0.79 26.2 Graphene (Li deposition) [15] 0.61 8.1 BL-  0.59 10.1 20 BL-  [32] 0.61 11.9 BL-  (13%) 1.33 30.78 46   [15] and BL-δ 6 .δ 6 [26] 1.10 20.5 δ 6 [28] 0.82 27.0 β 12 [26] 0.80 16.1 β 12 [27] 0.89 18.7 β 12 [28] 1.01 33.0 χ 3 [26] 0.60 11.5 χ 3 [27] 0.95 24.7 χ 3 [28] 0.79 26.2 Graphene (Li deposition) [15] 0.61 Before concluding, a few remarks need to be addressed.Superhard materials are usually semiconductors with large bandgaps or insulators, because the majority of electrons are bound into covalent bonds, depleting any excess electrons from electron transport.For instance, diamond is an insulator and the hardest material known.It may seem that metallicity and superhardness cannot coexist in one covalent material.However, this dilemma can be resolved in 2D borophenes, because the electron deficiency nature of boron atom enables strong localized covalent bonds and delocalized multicenter metal bonds in one system, which endows the system with superb mechanical performance and excellent metallicity, as evidenced by the high Young's moduli and metallicity in BL-δ 6 .

System
It is worth mentioning research on BL borophenes is still in its infancy and has just been ignited by two recent independent experimental investigations.We would also like to point out that in addition to the striking mechanical properties and superconductivity studied in this work, BL-δ 6 is expected to possess other interesting properties, such as anisotropic plasmonics [52] and ultrahigh thermal conductivity [53], awaiting more research.On the other hand, the excellent mechanical stability is also beneficial for BL-δ 6 transfer to or even growth on some inert substrates [54], for convenience of characterization.Our results also show that the interlayer bonds, while strengthening the bilayer, do not destroy the σ-bond resonance nor the conducting π-bonds [55], adding a fine tuning-knob to the properties.This also highlights the potential of covalent metals in the quest for high T c superconductors.
In our investigation, the tensile strain is applied to enhance the superconducting temperature of bilayer borophene.In addition to strain, the phenomenon of disorder, particularly correlated disorder, can also significantly enhance the superconductivity in a material [56,57].In bilayer borophene, manipulating disorder could pave new paths to enhance superconducting performance.On the other hand, a caveat common for all ≤2D materials, due to the Mermin-Wagner theorem, must be mentioned, since it restricts the stability, both merely structural and of the superconducting phase.Nevertheless, in any BL borophene realization, the sample finite size (and support by a 3D-substrate) should mitigate these concerns, although quantifying the size limitation is far beyond the scope of the present study.

Conclusions
In conclusion, we comprehensively investigated mechanical and superconducting properties of a bilayer borophene (BL-δ 6 ), which can be viewed as AB stacking, with interlayer covalent bonding, of earlier realized δ 6 borophene.The original good stability of δ 6 and the introduced covalent bond endows BL-δ 6 with very prominent energy-thermodynamic, thermal, and mechanical stability, suggesting that it has a high chance of occurrence in experiments.The strong directional σ B-B bonds and very compact atomic configuration give BL-δ 6 a 2D Young's modulus along the a direction (569.4N/m) higher than graphene (346.1 N/m), while ensuring that BL-δ 6 can sustain an ultimate 13% and 8% tensile strain along the a and b directions, respectively.Furthermore, according to a fully anisotropic solution of the ME equation, we predict BL-δ 6 is a conventional phonon-mediated 2D superconductor with T c as high as 20 K. We also highlighted the effects of tensile strain on EPCs and found that the T c can be boosted to 46 K at a 13% tensile strain along the a direction.To the best of our knowledge, it is the highest value currently known for borophene or elemental 2D materials.Our findings also consolidated the justification that covalent metals, such as borophene, should benefit the search for high T c superconductors.The concurrent superior mechanical performance and excellent superconductivity is scarce; thus, promising BL-δ 6 many potential applications, such as quantum interferometers, superconducting transistors, superconducting qubits, and even wear-resistant parts for superconducting devices.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/ma17091967/s1.The results of BL-δ 6 under 6% tensile strain along a direction.The meaning of (a-c) are the same as Figure S8. Figure S10: The results of BL-δ 6 under 10% tensile strain along a direction.The meaning of (a-c) are the same as Figure S8. Figure S11: The results of BL-δ 6 under 13% tensile strain along a direction.The meaning of (a-c) are the same as Figure S8.

Figure 1 .
Figure 1.(a) Planar and (b,c) side views of the atomic structure of BL-δ6.The yellow shaded region is its unit cell, and the green shaded region is the unit cell of δ6.There are three irreducible boron atoms colored as B1 (green), B2 (blue), and B3 (pink).The direction dependence of (d) Young's modulus and (e) Poisson's ratio of BL-δ6.(f) Stress-strain curves of BL-δ6, the vertical blue and red dashed lines denote the fractured stress along the a and b directions, respectively.

Figure 1 .
Figure 1.(a) Planar and (b,c) side views of the atomic structure of BL-δ 6 .The yellow shaded region is its unit cell, and the green shaded region is the unit cell of δ 6 .There are three irreducible boron atoms colored as B1 (green), B2 (blue), and B3 (pink).The direction dependence of (d) Young's modulus and (e) Poisson's ratio of BL-δ 6 .(f) Stress-strain curves of BL-δ 6 , the vertical blue and red dashed lines denote the fractured stress along the a and b directions, respectively.

Figure 2 .
Figure 2. (a) Electronic band structures of BL-δ6 by DFT (red solid lines) and Wannier90 (blue dashed lines).(b) Fermi surface, (c-e) top and side views of the charge density corresponding to the energy range shaded green in (a), and electronic band structures weighted by the py orbital of B1 and B3 atoms (f) and the pz orbital contributions of all boron atoms (g) of BL-δ6.

Figure 2 .
Figure 2. (a) Electronic band structures of BL-δ 6 by DFT (red solid lines) and Wannier90 (blue dashed lines).(b) Fermi surface, (c-e) top and side views of the charge density corresponding to the energy range shaded green in (a), and electronic band structures weighted by the p y orbital of B1 and B3 atoms (f) and the p z orbital contributions of all boron atoms (g) of BL-δ 6 .

Figure 3 .
Figure 3. (a) From left to right: phonon band structure of BL-δ6 weighted by EPC λqν with purple circles, isotropic Eliashberg function α 2 F and EPC λ(ω), and phonon density of states contributed by different kinds of boron atoms.(b-d) Vibrational modes of the B2u mode and the two Ag modes at Γ, as labeled in (a).

Figure 3 .
Figure 3. (a) From left to right: phonon band structure of BL-δ 6 weighted by EPC λ qν with purple circles, isotropic Eliashberg function α 2 F and EPC λ(ω), and phonon density of states contributed by different kinds of boron atoms.(b-d) Vibrational modes of the B 2u mode and the two A g modes at Γ, as labeled in (a).

Figure 4 .
Figure 4. (a) Momentum-dependent EPC λk across the full BZ and (b) momentum-dependent superconducting band gap ∆k at 10 K and projected onto the Fermi surface.(c) Variation in the superconducting gap ∆k with temperature, calculated by solving the ME equations in the isotropic approximation (yellow dots and dashed line interpolated) and with the fully anisotropic solution, where the purple shadowed regions indicate the magnitude distribution of the ∆k and the light green dots connected with the dashed line represents the average value of the entire anisotropic ∆k.

Figure 4 .
Figure 4. (a) Momentum-dependent EPC λ k across the full BZ and (b) momentum-dependent superconducting band gap ∆ k at 10 K and projected onto the Fermi surface.(c) Variation in the superconducting gap ∆ k with temperature, calculated by solving the ME equations in the isotropic approximation (yellow dots and dashed line interpolated) and with the fully anisotropic solution, where the purple shadowed regions indicate the magnitude distribution of the ∆ k and the light green dots connected with the dashed line represents the average value of the entire anisotropic ∆ k .

Figure 5 .
Figure 5. (a) Evolution with tensile strain along the a direction, of the Tc at three approximations (MAD, circles; isotropic ME, triangles; anisotropic ME, squares) and the distance between δ planes measured as B3-B3 bond length d.(b) Evolution of EPC λ (solid crosses), the log-weighted frequency ωlog, (circles) and maximum phonon frequency ωmax (triangles) versus the strain.(c) Phonon band structure weighted by EPC ωqν under 13% strain.(d) Isotropic Eliashberg α2F and EPC λ with 13% tensile strain along the a direction.Data with no strain are shown in gray, for comparison.

Figure 5 .
Figure 5. (a) Evolution with tensile strain along the a direction, of the T c at three approximations (MAD, circles; isotropic ME, triangles; anisotropic ME, squares) and the distance between δ 6 planes measured as B3-B3 bond length d.(b) Evolution of EPC λ (solid crosses), the log-weighted frequency ω log , (circles) and maximum phonon frequency ω max (triangles) versus the strain.(c) Phonon band structure weighted by EPC ω qν under 13% strain.(d) Isotropic Eliashberg α 2 F and EPC λ with 13% tensile strain along the a direction.Data with no strain are shown in gray, for comparison.
Figure S1: The evolution of T c with different k and q meshes.
Figure S2: The atomic structures of (a) δ 6 , (b) β 12 and (c) χ 3 .Figure S3: (a) The total potential energy fluctuation of BL-δ 6 during AIMD simulation at 500 K.(b) The top view and (c) side view of the atomic configuration of BL-δ 6 on Ag (111) surface after 5 ps of AIMD simulation at 500 K.

Figure S5 :
The strain energy per area under different kinds of strains for BL-δ 6 .FigureS6:The strain energy per area under different kinds of strains for (a) δ 6 , (b) β 12 , (c) χ 3 and (d) graphene.
Figure S7:The atomic and orbital resolved band structures of BL-δ 6.
Figure S8: (a) The electronic band properties of BL-δ 6 under 3% tensile strain along a direction, red line and blue dashed line denote the results calculated from DFT and Wannier90, respectively.(b) Phonon band structure of BL-δ 6 weighted by EPC λ qν with purple circles, isotropic Eliashberg function α 2 F and EPC λ(ω), PHDOS from different kinds of boron atoms' contribution.(c) Evolution of the superconducting gap ∆ k a function of temperature, calculated by solving the ME equations in the isotropic approximation (yellow dots and dashed line interpolation) and with a fully anisotropic solution where the purple shadowed regions indicate the magnitude distribution of the ∆ k and the light green dots connected with dashed line represents the average value of the entire anisotropic ∆ k .Figure S9: