Prediction of superconductivity in bilayer borophenes.

Borophenes and related two-dimensional materials have exhibited many exotic properties, especially for superconductivity, although the superconductivity of single-layer borophene is suppressed by the strains or doping from its substrates. Intriguingly, bilayer (BL) borophenes can be stabilized by appropriate pillar density and hexagonal holes density, rather than being supported by Ag(111) or Cu(111) substrates. Thus, we studied the two most stable structures, namely BL-B8 and BL-B30, stabilized by the above-mentioned two methods. Within density functional theory and Bardeen-Cooper-Schrieffer theory framework, their stability, electron structures, and phonon properties, as well as possible superconductivity are systematically scrutinized. The metallic BL-B8 and BL-B30 exhibit intrinsic superconducting features with superconductivity transition temperatures (T c) of 11.9 and 4.9 K, respectively. The low frequency (below 400 cm-1) consisting of out-of-plane vibrations of boron atoms plays crucial rule in their superconductivity. In particular, a Kohn anomaly appears at the Γ point in BL-B8, leading to substantial electron-phonon coupling. Here, our findings will provide instructive clues for experimentally determining the superconductivity of borophene and will broaden the two-dimensional superconductor family. This journal is © The Royal Society of Chemistry.

Introduction

Inspired by the fascinating properties of graphene,[1] explorations of related two-dimensional (2D) materials have sprung up in recent years.[2-5] Boron, as a neighbor of carbon, also possesses a short covalent radius with sp2 hybridization, leading to the flexibility to form planar boron clusters,[6,7] cage-like boron fullerene,[8,9] and nanotubes.[10] Moreover, similar to graphene, boron atoms can also exist in a 2D graphene-like monolayer, namely borophene. However, different to graphene, the one electron-deficient property of boron atoms results in the biggest challenge to synthesizing borophene. To balance the electrons, hexagonal hole incorporation is an effective method, in which the hexagonal holes can obtain extra electrons from the filled hexagons. In this way, 2D boron sheets (α-, β-,[11,12] g-boron sheet[13,14]) have been predicted to be stable. Moreover, borophene deposited on metal substrates is an alternative method to balance the surplus electron of boron. Following the clues given by first-principle calculation,[15] the single-layer (SL) borophene is stabilized by metal passivation of Cu(111),[16] Ag(111),[17,18] Ir(111),[19] Al(111)[20] and Au(111) surfaces.[21] In addition, bilayer (BL) borophenes, whose stability is related to hexagonal holes and interlayer pillars, have been systematically investigated, greatly enlarging the family of 2D boron sheets (here, a pillar means a direct chemical bond, which connects the upper and lower boron atoms layers).[22] Moreover, BL borophene has been successfully synthesized on atomically flat single-crystal Ag(111).[23] More recently, single-crystalline bilayer borophene with large-size on the Cu(111) surface has been obtained by molecular beam epitaxy.[24] However, their stable phases are supported by substrates, which provide additional electrons for the bonding of additional boron atoms.

As discussed above, both SL and BL borophenes can exhibit various structural polymorphs depending on their supporting substrates, leading to intriguing physical and chemical properties,[25] such as ideal strength, negative Poisson’s ratio[26] and high thermal conductivity.[27] Therefore, borophenes have been explored for desirable applications in many fields, such as as electrode materials,[28-32] sensors,[33-35] and in energy storage,[36,37] flexible devices,[38]etc. Furthermore, superconductivity, as one of the famous macroscopic quantum phenomena, has also been predicted to exist among the family of borophenes, due to its light atomic weight as well as the large density of states (DOS) near the Fermi level. SL β12 and χ3 borophenes with vacancy fractions of 1/6 and 1/5 are predicted to be phonon-mediated superconductors with superconducting transition temperatures Tc of 10–20 K.[39,40] The Tc of SL borophenes as a function of the hexagon hole density shows a V-shape.[41] A honeycomb borophene (h-B2) is also predicted to be a 30 K superconductor without considering the imaginary phonon modes.[41] Besides, super-borophene has a critical temperature Tc of 25.3 K at ambient conditions.[42] Moreover, some boron-based 2D superconductors are also predicted, such as Mo2B2,[43] W2B2,[44] AlB2,[45] Li2B2,[46] Li2B7,[47] B2C,[48] XB6 (X = Al, Ga, In),[49,50] B2O,[51] TiB4,[52] TiB7,[53] MnB3,[54]etc., further indicating the rich superconductivity in boron families.

However, although the SL h-B2, α-, β12 and χ3 borophenes have been synthesized in experiments and predicted to be intrinsic superconductors with high Tc by many theoretical works, their superconductivity has not yet been unambiguously observed in experiments.[55] The undetected superconductivity of SL borophene could be attributed to the interference from the substrates, which can introduce strains and/or dopings.[41,45,56] Thus, more effective synthesis methods without the supportive substrates for SL borophenes are urgently required, but current structures of borophenes can only exist on a special substrate. In particular, BL borophene, synthesized by layer-by-layer techniques, could serve as a candidate to experimentally explore the superconductivity of borophene. Besides, some unique properties appear in few-layer borophenes. For example, the in-plane negative Poisson’s ratios in the monolayer become positive in the layered borophene, while out-of-plane negative Poisson’s ratios are preserved;[57] the electronic and magnetic properties of BL borophene can be modulated by transition metal atoms intercalation;[58] P6̄-boron possesses a topologically nontrivial Dirac nodal line, which is protected by the mirror reflection symmetry.[59] BL-α borophene has been synthesized recently and it has a higher crystallinity and local work function, but it can’t maintain a free-standing structure without the supportive substrate.[23] Therefore, further studying the BL borophene will significantly expand the phase space for boron-based nanomaterials.

In the present work, the most stable BL borophenes BL-B8 and BL-30, stabilized by appropriate pillar density and hexagonal holes, respectively, are theoretically dissected via first-principles calculations. Their bonding nature, stability, electronic, phonon vibrational properties and electron–phonon coupling (EPC) are systematically studied. Results indicate that they are not only 2D metals with high stability, but also are 2D intrinsic superconductors with Tc of 11.9 and 4.9 K, respectively. Their superconductivity originates from the coupling between electrons of B-p orbitals and soft or flat phonon modes in the low-frequency region occupied by out-of-plane vibrations. In particular, a Kohn anomaly appears at the Γ point in BL-B8 and provides substantial EPC, leading to a higher Tc than the SL borophenes.

Computational methods

The cell optimizations and electronic structure calculations were carried out using the ab initio plane-wave pseudopotential method as implemented in the Vienna ab initio simulation package (VASP).[60,61] The plane-wave energy cutoff of 500 eV and generalized gradient approximation (GGA)[62,63] formulated by Perdew–Burke–Ernzerhof (PBE)[64] were used for the exchange correlations. A vacuum space of 15 Å along the z direction was adopted to avoid the unphysical interactions within the periodic images. For Brillouin zone (BZ) samplings, the Γ-centered 18 × 16 × 1 and 13 × 13 × 1 k-point meshes were used for rectangular and hexagonal primitive cells, respectively. The ab initio molecular dynamics (AIMD) simulations, with constant number, volume, temperature (NVT) ensemble, were carried out to assess the thermal stabilities, lasting for 5 ps with a time step of 1 fs. Some data post processing steps in the VASP calculations were dealt with using the VASPKIT code.[65] The convergence criteria were set for an energy of 10−5 eV and force of 0.01 eV Å−1.

Within density functional perturbation theory (DFPT)[66] and the Bardeen–Cooper–Schrieffer (BCS) theory,[67] the superconductivity was studied using the open source QUANTUM ESPRESSO (QE) package.[68,69] Here, the phonon spectra were calculated by using the finite displacements method within its primitive cell. The norm-conserving PBE pseudopotentials were chosen to describe the electron–ion interactions. The plane-wave kinetic energy cutoff and the charge–density cutoff were tested to be 70 and 280 Ry, respectively. The self-consistent electron densities for rectangular and hexagonal cells were calculated on 36 × 32 × 1 and 24 × 24 × 1 k-meshes, respectively. Moreover, 9 × 8 × 1 and 6 × 6 × 1 q-meshes were selected for them to calculate the dynamical matrices and the EPC matrix elements.

Results and discussion

Atomic structures

As reported by previous work,[22] the BL borophene sheets can be obtained by optimal proportions of pillars and hexagonal holes to balance the surplus electrons. Therefore, we focus our study on the most stable crystals, obtained by the above two methods. As shown in Fig. 1(a and b), this BL is constructed with pillar density of 1/4, where pillar density is calculated by (number of atoms forming pillars)/(number of atoms in the unit cell).[22] Based on the number of atoms within its primitive cell, we label it as BL-B8. Compared with other BL sheets, it is the most stable structure with the lowest energy stabilized by pillars.[22] BL-B8 crystallizes in an orthorhombic space group, Pmmm (no. 47), and exhibits three types of B–B bonds with bond lengths of 1.65–1.98 Å, which are comparable to those of the SL borophene (1.62–1.88 Å)[28,32] and bulk borophene (1.73 Å).[70] The seven-coordinated boron atoms, marked as orange atoms in Fig. 1(a and b), serve as pillars to stabilize the top and bottom planes, which balance the surplus electrons and maintain the stability. The high-symmetry paths in its 2D BZ for BL-B8 are shown in Fig. 1(c).

Fig. 1

Top (left panel) and side (middle panel) views for (a and b) BL-B8 and (d and e) BL-B30. Right panel: the 2D Brillouin-zone of (c) BL-B8 and (f) BL-B30 and the points of high symmetry. The rectangle and rhombus with black lines represent the primitive cells. Here, boron atoms in different coordination environments are indicated by different colors in the atom structure.

Another structure of BL borophene, also from ref. 22 and as shown in Fig. 1(d and e), is stabilized by isolated hexagonal holes separated by boron sheets, and has low formation energy among bilayer structures with different hexagonal holes density. For simplicity, we name it as BL-B30. The entitled BL-B30 optimizes in a hexagonal lattice and its high-symmetry paths are along Γ–M–K–Γ as presented in Fig. 1(f). As clearly distinguished in Fig. 1(d and e), there are seven different surroundings of boron atoms and the B–B bonds are in range of 1.65–1.87 Å. Different from BL-B8, here, the pillars are formed by two different B–B bonds: one is 1.72 Å and another 1.78 Å. Moreover, we compare atomic energy of per boron in BL-B8 and BL-B30 with typical SL borophenes (δ6, β12 and χ3). Their corresponding lattice constants and energies are listed in Table 1. Obviously, both our investigated BL borophenes are more stable than the SL borophenes, especially BL-B30. Besides, BL-B30 has more reliable stability when compared with BL-B8, which is in accordance with their bonds’ distributions. Thus, the combination of pillars and hexagonal holes with suitable ratios is beneficial for the stability of BL borophenes.[11,22]

The lattice constants and energies per atom of δ6, β12, χ3 borophenes, BL-B8, and BL-B30

Configuration	Lattice constant (Å)	Energy (eV per atom)
δ 6	a = 2.97; b = 1.62	−6.21
β 12	a = b = 4.45	−6.26
χ 3	a = 2.93; b = 5.00	−6.27
BL-B8	a = 2.85; b = 3.23	−6.37
BL-B30	a = b = 6.56	−7.12

Thermal and mechanical properties

The stabilities of 2D materials are very important. Materials with good stability may be synthesized by experiments. AIMD simulations controlled by the Nosé-Hoover thermostat[71] are adopted to study the thermal properties of BL-B8 and BL-B30. To take periodic boundary conditions into consideration, 3 × 3 × 1 supercells are built. As shown in Fig. 2(a and b), the free energies of BL-B8 and BL-B30 maintain nearly constant fluctuations during the 5 ps simulation period at about 1500 K. At about 1800 K, there may exist structural disruptions with collapsed fluctuations, indicating their theoretical melting temperatures at around 1500–1800 K.

Fig. 2

The variations of free energy of (a) BL-B8 and (b) BL-B30 under the 5 ps AIMD simulations. The orientation-dependent (c) Y and (d) ν as a function of the polar angle for BL-B8 (red lines) and BL-B30 (blue lines).

Moreover, the mechanical properties of BL-B8 and BL-B30 are explored within their rectangle unit cells. Then, the elastic constants C are obtained bywhere Es and ε indicate strain energy and strain, respectively. Here, ε is calculated by , (a and a0 are the lattice constants of the strained and non-strained crystals, respectively.) and ε (ε) is the uniaxial strain along the x(y) direction, and ε is the shear strain. Then, the independent elastic constants (C11, C12, C22 and C66) within the rectangle cell can be calculated from the function of strain energy vs. applied strains.[72] The corresponding elastic constants for BL-B8 and BL-B30, compared with SL borophenes, are listed in Table 2. Clearly, the C11, C22 and C66 of BL-B8 and BL-B30 are much larger than those of β12 and χ3.[73] Additionally, both BL-B8 and BL-B30 satisfy the Born criteria,[74]C11 > 0, C66 > 0 and C11C22 − C122 > 0,[74,75] indicating they are mechanically stable.

Calculated parameters of C (N m−1), Young’s modulus (N m−1), and Poisson’s ratio for BL-B8 and BL-B30. Corresponding results for β12 and χ3 from ref. 73 are also presented for comparison

Compounds	C 11	C 12	C 22	C 66	Y x	Y y	ν x	ν y
β 12	188	36	214	64	182	241	0.17	0.18
χ 3	195	36	188	71	188	181	0.19	0.18
BL-B8	356	17	618	161	355	618	0.03	0.05
BL-B30	509	40	390	146	505	387	0.1	0.08

In addition, BL-B8 and BL-B30 belong to the orthogonal 2D system, their orientation-dependent Young’s moduli (Y) and Poisson’s ratios (ν) can be defined as[51,72,76]andrespectively, where Δ = C11C22–C122, a = cos(θ) and b = sin(θ). Therefore, Y(θ) and ν(θ) can be derived as follows:and

As shown in Fig. 2(c), the variations of Y for BL-B8 and BL-B30 are spindle-like shapes, which are similar to that of the B2O monolayer.[51] However, the distributions of the orientation-dependent Y are opposite. For BL-B8, the in-plane minimal value of Y along the x direction is 355 N m−1, and the maximal value of 618 N m−1 is along y direction. Contrary to BL-B8, for BL-B30 the maximal value of 505 N m−1 and the minimal value of 387 N m−1 are along the x and y directions, respectively. Moreover, their orientation-dependent ν variations show butterfly-like shapes, as shown in Fig. 2(d), further indicating their fully anisotropic traits.

Electron structures

Based on the electron localization function (ELF),[77] plotted in Fig. 3(a–d), the bonding natures in BL-B8 and BL-B30 are studied. In general, a region with ELF = 1 means strong localization of covalent bonding electrons, ELF values close to 0 are a low electron density, and regions with ELF = 0.5 indicates a homogeneous electron gas. Obviously, both in BL-B8 and BL-B30, all the ELF values for B–B bonds are close to 0.8, showing their covalent bonding nature. In particular, large ELF values around 1 are found along the pillars, meaning stronger covalent bonding than the B–B bonds within the planes. Such results indicate again that the pillars in BL-B8 and BL-B30 are important for their stability.[22] Furthermore, the difference charge plots for BL-B8 and BL-B30 are exhibited in Fig. 3(e–h). Evidently, boron atoms in different locations show different electron transfer capabilities, but electron transfer occurs within layers, balancing the surplus electrons and enhancing stability for BL-B8 and BL-B30. Furthermore, the formations of pillars in BL-B8 and BL-B30 are beneficial to the electrons’ transfer and compensation for unsaturated valence electrons in sp2 hybridization.[22] According to the Bader charges analysis,[78] the net charge transfer between B–B covalent bonds in BL-B8 is 0.86|e|, much smaller than that of BL-B30 (2.4|e|). Interestingly, relatively large amounts of electron transfer happen in BL-B8 and BL-B30, hinting at their performances in superconductivity.

Fig. 3

The ELF plots for (a and b) BL-B8 and (c and d) BL-B30. The calculated difference charge density for (e and f) BL-B8 and (g and h) BL-B30. The isosurface value for difference charge density is chosen to be 0.01 a.u.

The orbital-projected band structures and the electronic density of states (DOS) are exhibited in Fig. 4. The same as SL borophenes, BL-B8 and BL-B30 show intrinsic metallic features with large amounts of electron states concentrating near the Fermi level. Clearly, their metallic nature mainly originates from the B-p orbitals with small amounts of s orbitals [Fig. 4(b and f)]. Moreover, the in-plane B-p (B-pσ) and out-of-plane B-p (B-pπ) orbitals have almost equal proportion in the projected DOS of BL-B8 [Fig. 4(c)]. However, as for BL-B30, its metallicity mostly stems from B-pπ orbitals partially hybridized with few B-pσ orbitals [Fig. 4(g)]. Furthermore, the Rashba effect[79,80] is explored, and the band structures calculated with spin–orbit coupling (SOC) interactions are presented in Fig. 4(d and h). Clearly, both in BL-B8 and BL-B30, extremely weak SOC effects happen and no bands split in proximity to the Fermi level. Thus, the Rashba effect will not be considered for the studied systems.

Fig. 4

Orbital projected band structures without SOC (left panel), with SOC (right panel), as well as projected DOS without SOC (middle panel) of (a–d) BL-B8 and (e–h) BL-B30.

Phonon properties and EPC

The phonon properties and possible superconductivity of the metallic BL-B8 and BL-B30 are further investigated. In terms of the displacement directions of B atoms,[81,82] their resolved phonon spectra are given in Fig. 5(a and e). Here, the phonon kp theorem is adopted to sort the phonon branches based on their eigenvectors[83-85]where and Δ are the displacement of the atom j in the eigenvector of (k, σ) vibration modes and a small wave vector, respectively. Clearly, all the phonon modes in BL-B8 and BL-B30 are positive, further suggesting their stability without any dynamical instability. The three acoustic branches are formed by the out-of-plane (ZA), in-plane transverse (TA), and in-plane longitudinal (LA) modes. As shown in Fig. 5(a), in BL-B8, the low-frequency region (below 300 cm−1) is mostly dominated by the in-plane modes of B atoms (B-xy), whilst the LA and the first three optical branches are mainly occupied by the out-of-plane B modes (B-z). In the middle-frequency region (300–1000 cm−1), phonon vibrations also mainly originate from the B-xy. In the high-frequency region (larger than 1000 cm−1), contributions are mostly from B-z modes. As for BL-B30 (Fig. 5(e)), the low-frequency region below 400 cm−1 is related mainly to the B-z modes, along with some B-xy modes in TA, ZA and the first two optical branches. However, phonon frequencies higher than 400 cm−1 mainly consist of the in-plane vibrations of B atoms. The phonon density of states (PhDOS) distributions are in line with above analysis of the vibration modes [Fig. 5(b and f)]. Moreover, the highest vibration frequency in BL-B8 and BL-B30 with light weight reaches to ∼1400 cm−1, comparable to those of δ6 (∼1200 cm−1), β12 (∼1290 cm−1) and χ3 (∼1411 cm−1) borophenes, according well with their strong covalent bond features.[49,51]

Fig. 5

The resolved phonon spectrum in terms of the displacement directions of B atoms and the PhDOS for (a and b) BL-B8 and (e and f) BL-B30. The phonon spectrum weighted by EPC λqν and the distribution of the λ(ω) upon α2F(ω) for (c and d) BL-B8 and (g and h) BL-B30. The magnitude of the λqν is displayed with an color bar scale in (c and g).

We now turn to their potential superconductivity, according to the magnitude of the electron–phonon coupling (EPC) λqν within the whole frequency region, the calculated Eliashberg spectral function α2F(ω) and the cumulative EPC strength λ(ω). Based on the Migdal–Eliashberg theory,[86,87] the λqν can be calculated bywhere γqν, ωqν, and N(EF) are the phonon linewidth, phonon frequency and the electronic DOS near the Fermi level, respectively. Some related superconducting parameters compared with SL borophenes are listed in Table 3. Our calculated results for the three typical SL borophenes are consistent with previous work.[42] According to the λqv of BL-B8 [Fig. 5(c)], the superconductivity mainly originates from the low-frequency vibration modes (below 200 cm−1) with an apparent Kohn anomaly. There are two optical branches that just near 200 cm−1 at the Γ point, intersecting the acoustic branches toward the zone boundary. But for BL-B30, three flat optical branches appear below frequency region of 400 cm−1 in the BZ. As indicated by λqv in Fig. 5(g), the superconductivity of BL-B30 also originates from the low-frequency phonons.

The superconductive parameters of N(EF) (in units of states/spin/Ry/cell), ωlog (in K), λ and Tc (in K) of χ3, δ6, β12 borophenes, BL-B8, and BL-B30

Compounds	N(EF)	ω log	λ	T c
χ 3	6.68	417.57	0.64	11.5
δ 6	2.61	409.29	0.82	20.1
β 12	7.51	371.71	0.71	13.7
BL-B8	7.98	495.51	0.61	11.9
BL-B30	27.01	541.05	0.47	4.9

In terms of α2F(ω), the EPC strength λ(ω) can be obtained byand

Then, the logarithmic average frequency ωlog is derived by

Four clearly distinguished regions in α2F(ω) are indicated by black dots [Fig. 5(d and h)]. Here, region I (below 200 cm−1) and II (200–400 cm−1) are both in the low-frequency region. Tthe intermediate region is from 400 to 1000 cm−1 (region III) and a high-frequency region is above 1000 cm−1 (region IV). Although the total EPC are attributed to all the phonon modes over the entire frequency range, the relative contributions of each region are considerably different.[88] For BL-B8, region I is mainly dominated by the ZA modes and accounts for only about 10% of the total EPC constant (λ = 0.61), and region II with a Kohn anomaly at the Γ point leads to substantial EPC, and provides a contribution of 57% to the total λ. Thus, the low-frequency phonons make the main contribution to the total EPC, approximately 67%. At the same time, region III in the middle-frequency region also has a non-negligible contribution of 27.9% to the total EPC. However, region IV provides a negligible contribution, only about 5%. As mentioned above, the Kohn anomaly in region II dominated by the out-of-plane of boron atoms vibrations is the crucial factor for producing a high value of Tc for BL-B8. Therefore, based upon the above results, the Tc can be calculated by the McMillian–Allen–Dynes formula:[89]where the effective screened Coulomb repulsion constant μ* is selected to be 0.1. Therefore, the Tc of BL-B8 is predicted to be 11.9 K, which is comparable to those of χ3 (11.5 K) and β12 (13.7 K) borophenes. As for BL-B30, the flat optical branches in the low-frequency region consisting of out-of-plane boron atoms vibrations are important for its superconductivity. The low-frequency region (regions I and II) make a major contribution, about 70%, to the total EPC (λ = 0.47). Regions III and IV account for only approximately 19% and 11% to the total EPC, respectively. These results are accordance with its λqν distributions. Thus, the Tc of BL-B30 is finally integrated to be 4.9 K, smaller than the BL-B8. This is understandable since the great coulomb repulsion in BL-B30 can result in the formation of few Cooper pairs.

Conclusions

In summary, we have performed a systematic theoretical study on BL-B8 and BL-B30, including stability, electron, phonon properties and even EPC. They are both intrinsic 2D metals with good stability. Within the BSC microscopic theory, they are further determined to be 2D supercondutors with Tc of 11.9 K and 4.9 K. The EPC in low-frequency occupied by out-of-plane boron vibrations (below 400 cm−1) is crucial for their superconductivity. Here, a Kohn anomaly appears at a Γ point in BL-B8 which provides substantial EPC, while flat optical branches in BL-B30 contribute a little to the EPC with a great quantity of electrons near the Fermi level. In addition, their superconductivity can be tuned by strain and doping engineering with further efforts. Our work would fill in the blanks of superconductivity in BL borophenes, and provide clues for experimental explorations.

Conflicts of interest

There are no conflicts to declare.