Theoretical prediction of superconductivity in monolayer h-BN doped with alkaline-earth metals (Ca, Sr, Ba)

Two-dimensional (2D) materials have attracted considerable attention as candidates for materials with novel physical properties [1]. The superconductivity of 2D materials is expected to be different from that in three-dimensional (3D) bulk materials and promising for hosting unique phenomena such as quantum phase transition. Therefore, finding a new 2D superconductor is one of the important topics in this field. Although several types of 2D superconductors have been experimentally and/or theoretically investigated, most of them are metallic/semimetallic materials, such as metal atomic layers on semiconductor surfaces [2–4], layered materials such as FeSe [5, 6] and NbSe₂ [7], and electron-doped graphene [8–15]. The other class of 2D superconductors is doped 2D semiconductors or insulators. In doped 2D semiconductors, superconductivity induced via gate doping using ionic liquid has been recently reported in MoS₂ [16–18] and ZrNCl [19]. The other method for realising superconductivity is via chemical doping using alkali or alkaline-earth metals. Although the tunability of the doping amount in chemical doping is lower than that in gate doping, the former offers additional electron–phonon interaction processes owing to the low-energy phonon modes and electronic states induced by the dopant, thereby increasing the superconducting transition temperature. In addition, the superconductivity of chemically doped MoS₂ [20] and black phosphorus [21] has been reported. However, investigations regarding the possibility of fabricating a superconductor using chemically doped 2D insulators are still limited.

We focussed on hexagonal boron nitride (h-BN) as a candidate for a 2D superconductor. Bulk h-BN has a wide bandgap of approximately 6 eV, and it exists in the monolayer form. Although the application of h-BN as an inert and atomically flat insulating layer to realise novel 2D-material-based electronics [22–25] has been widely investigated, its possibility as a superconductor has been untapped until recently. In our previous study, we reported that a Li-intercalated h-BN bilayer became a phonon-induced superconductor with a Tc of up to 25 K [26], thereby indicating that both the phonon properties and electron–phonon coupling (EPC) in doped h-BN are appropriate for realising superconductivity. However, the fabrication of a superconducting Li-intercalated h-BN bilayer is experimentally challenging, as the stacking order of the h-BN sheets must differ from that of the original bulk structure to attain a high Tc. Therefore, to seek a promising candidate for an h-BN-based superconductor with a more experimentally feasible structure than the bilayer structure, we investigated the properties of monolayer h-BN doped using group-1 (Li, Na, K) and group-2 (Be, Mg, Ca, Sr, Ba) atoms by performing ab initio calculations of the electronic and phononic properties as well as EPC.

Consequently, we reveal that Sr- and Ba-doped monolayer h-BN and Ca-doped monolayer h-BN with 3.5% tensile strain are energetically stable and become superconductors with Tc of approximately 5.83, 1.53, and 10.7 K, respectively. The strain in the Ca-doped case can be realised using the deformation of a polymer substrate or by selecting a substrate to fabricate the h-BN monolayer. Furthermore, the Tc of the Ca-doped h-BN is higher than that of typical BCS superconductors such as Nb (9.29 K) and Pb (7.2 K). Therefore, the superconducting gap in the Ca-doped h-BN can be observed by performing transport measurements and/or low-temperature scanning tunneling spectroscopy [4, 27]. Based on our results, we propose a new approach for the on-surface synthesis of 2D superconductors based on an h-BN monolayer.

The ab initio calculations were performed using the Quantum-ESPRESSO package [28, 29]. The model structure of the doped monolayer h-BN is depicted in figure 1, wherein we assumed a $\sqrt{3}{\times}\sqrt{3}\mathrm{R}3{0}^{{\circ}}$ periodicity, which corresponds to one adsorbate per three h-BN primitive cells. Because the local-density approximation for the exchange–correlation functional can satisfactorily reproduce the structure of h-BN [30, 31], we used the norm-conserving pseudopotential by employing the Troullier–Martins method [32] with LDA. The plane-wave basis set was introduced, and the kinetic-energy and charge cutoffs were set to 60 and 240 Ry, respectively. The Brillouin zone was sampled using a 24 × 24 × 1 k-mesh. The structures were relaxed before the calculations of band dispersions and EPC until the force acting on each atom reduced to less than 1.0 × 10⁻⁴ Ry/a.u. The dynamical matrices and phonon deformation potentials were calculated using a Gamma-centred 12 × 12 × 1 q-mesh. In this study, we used the EPC matrix element that was directly obtained using Quantum-ESPRESSO calculations, as the dispersive band structures above the Fermi level hampered the convergence in disentanglement and Wannierisation procedures to achieve a high-precision Wannier interpolation [33]. To compensate for the effect due to the sparse sampling of the Brillouin zone, we used a finer 120 × 120 × 1 k-mesh to interpolate the Fermi surface in the calculation of the Eliashberg function, defined as follows:

$$\begin{align*}\hfill {\alpha }^{2}F\left(\omega \right)& =\frac{1}{{N}_{k}{N}_{q}N\left(0\right)}\sum _{nk,mq,\nu }{\left\vert {g}_{nk,mk+q}^{\nu }\right\vert }^{2}\delta \left({E}_{nk}\right)\hfill \\ \hfill & \quad {\times}\enspace \delta \left({E}_{mk+q}\right)\delta \left(\omega -{\omega }_{q}^{\nu }\right),\hfill \end{align*}$$

where N(0) denotes the density of states (DOS) at the Fermi level; k and q denote the electron and phonon wavevectors, respectively; N_k and N_q denote the total numbers of k- and q-mesh grid points, respectively; indexes n and m denote the electronic band indexes; ν is the phonon band index; ${\omega }_{q}^{\nu }$ is the phonon energy; E_nk and E_mk+q denote the eigenvalues of the Kohn–Sham wavefunctions with respect to the Fermi level; ${g}_{nk,mk+q}^{\nu }$denotes the electron–phonon matrix element.

Zoom In Zoom Out Reset image size

The superconducting transition temperature was estimated using the McMillan–Allen–Dynes formula [34, 35] as follows:

$$\begin{equation*}{T}_{c}=\frac{{\omega }_{\mathrm{log}}}{1.2}\mathrm{exp}\left[\frac{-1.04\left(1+\lambda \right)}{\lambda \left(1-0.62{\mu }^{{\ast}}\right)-{\mu }^{{\ast}}}\right],\end{equation*}$$

where μ* denotes the Coulomb repulsion pseudopotential, and we set μ* = 0.14 following reference [26]. The EPC constant λ was obtained by setting ω → ∞ in the following expression:

$$\begin{equation*}\lambda \left(\omega \right)=2{\int }_{0}^{\omega }\mathrm{d}{\omega }^{\prime }\frac{{\alpha }^{2}F\left({\omega }^{\prime }\right)}{{\omega }^{\prime }}.\end{equation*}$$

In addition, we decomposed the EPC constant at the Fermi level with respect to the electron wavevector k to determine the contribution of each electronic band. This k-resolved EPC is expressed as follows [10, 36]:

$$\begin{equation*}{\lambda }_{k}=\sum _{q,\nu }\delta \left({\varepsilon }_{k+q}\right)\frac{{\left\vert {g}_{nk,mk+q}^{\nu }\right\vert }^{2}}{{\omega }_{q\nu }}.\end{equation*}$$

First, we determined the stable structure in the doped h-BN and graphene by comparing their total energies, which were obtained by varying the values of the lattice constant, 'a' and the adsorbate–substrate distance, 'h' (see figure 1). We considered the on-top site of N and B, the bridge site between the N–B bond, and the hollow site. The calculation results are summarised in appendix A. For the monolayer h-BN, no minimum was observed in its potential-energy curve as a function of h for Li, Na, K, Be, and Mg for cases where h < 5.5 a.u. This implies that the adsorption of these atoms was unstable. This is further confirmed by the appearance of the phonon mode with an imaginary frequency, which is depicted as a negative frequency in the phonon-dispersion curve (see appendix B). However, Ca, Sr, and Ba were stably adsorbed on the h-BN monolayer at the hollow site. As depicted in figure 2(a), in the Ba case, only a single minimum was observed in the potential-energy curve at approximately 5.0 a.u. By contrast, in both the Ca and Sr cases, the potential-energy curves exhibited two minima, one at approximately 4.7 a.u. and the other at approximately 7.0 a.u. With the equilibrium lattice constant in the Sr case, the closer position of 4.66 a.u. was the global minimum; however, in the Ca case, the further position of 7.0 a.u. was the global minimum and the closer position was a metastable state. To realise a high Tc, a smaller distance between the adsorbate and h-BN sheet is favourable for enhancing the coupling between the low-frequency vibrational mode of the adsorbate and the electronic state in h-BN. Therefore, we investigated the possibility of switching the global-minimum position by introducing a tensile strain on the h-BN monolayer in the Ca-doped case. As shown in figure 2(b), the position of the global minimum depended on the lattice constant, i.e. the strain. If the tensile strain is greater than 3.5%, then the Ca-adatom favors a closer position of 4.66 a.u. The structure parameters, lattice constants, and distances between the adsorbate and layer determined via the calculations above are summarised in appendix C. In addition, we provide the parameters of Ca-, Sr-, and Ba-doped graphene used in the calculations for comparison.

Zoom In Zoom Out Reset image size

We evaluated the superconducting transition temperature of Sr- and Ba-doped h-BN with no strain and that for Ca-doped h-BN with tensile strain. The obtained Tc, DOS at the Fermi level (E_F), and EPC constant (λ) are summarised in table 1. Notably, our calculation results for Tc and λ Ca-doped graphene were slightly lower than those in a previous report [8]. The differences in the calculation results are reasonable because compared with the previous report, we used a slightly larger lattice constant as well as a larger distance between the Ca atom and the graphene layer.

As depicted in figure 3 and table 1, the doped h-BNs yielded a higher EPC constant and Tc than the doped graphene systems (graphene is denoted as G). The higher EPC constant is mainly attributed to the large ${\alpha }^{2}F\left(\omega \right)$ in the region of low-energy phonons, as shown in figure 3. These phonon modes were localised at the adsorbate, and we highlighted them in the phonon dispersion curves in figure 4. Both the in-plane (red) and out-of-plane (blue) vibrations of the adsorbate appeared in the energy range of 0–30 meV. Although similar low-energy phonon modes by the adsorbate were also exhibited in the graphene cases, as compared with the modes in the h-BN cases, these modes contributed less to ${\alpha }^{2}F\left(\omega \right)$, as depicted in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To clarify the difference between h-BN and graphene, we analyzed the electronic band structure depicted in figure 5. The projection of the DOS from the adsorbate s-orbital onto the electronic-band structures indicated the presence of adsorbate bands that crossed the Fermi level in all systems. These bands correspond to the interlayer state reported in previous studies [8, 15]. In doped graphene systems, the interlayer state is crucial for determining Tc. However, the similarity in the interlayer state between doped graphene and h-BN indicates that the Tc in these systems is governed by other factors. Other than the interlayer states, several bands that originated from h-BN and graphene orbitals crossed the Fermi level.

Zoom In Zoom Out Reset image size

These band structures around the Fermi level resulted in Fermi surfaces of different shapes and positions, as depicted in figure 6. In both the h-BN and graphene cases, the Fermi surface around the K point originates from the interlayer state, and those around the Γ point correspond to the band from the h-BN or graphene states. The shape of the Fermi surfaces from the interlayer states resembled each another in both the doped h-BN and graphene cases, except for the size of the pocket. However, the shapes of the inner Fermi surfaces around the Γ point between the doped h-BN and graphene cases differed significantly. In the doped h-BN case, the inner Fermi surface exhibited a snowflake-like shape, but it was hexagram-like in the doped graphene case. The projection of the k-resolved EPC (λ_k) shows that the electron–phonon interaction in the inner Fermi surface is crucial for determining Tc. Compared with doped graphene, the doped h-BNs exhibited greater λ_k in general in the Brillouin zone. The most prominent enhancement in λ_k was observed around the Γ point.

Zoom In Zoom Out Reset image size

The enhancement in λ_k around the Γ point in the Ca-doped h-BN is attributed to the spatial distribution of the wavefunction of the electronic states. In figure 7, we plot the charge-density distribution of each state, which is indicated using crosses (× and +) in figure 5. The distinctive characteristics were shown in the state near the Γ point (marked by +) of the Ca-doped h-BN: the strong intensity of charge distribution connecting the Ca and B atoms. Such a feature was retained but weakened at states that are distant from the Γ point and close to the M point (marked by ×). In the Ca-doped graphene case, the charge-distribution intensity between the Ca and graphene layers was considerably small in both states around the Γ and M points. In addition to the interaction strength between Ca and the substrate electronic state, symmetry is important. Because of the symmetry difference in monolayer h-BN(C_3v) and graphene(D_6h), the number of possible rotation/reflection symmetry operations in the doped h-BN is fewer than that in the doped graphene. Therefore, the restriction on the coupling between intrinsic π* − π* electrons mediated by the in-plane vibration of the Ca adsorbate is reduced in the h-BN case. This explains the coincidence of the peak position in ${\alpha }^{2}F\left(\omega \right)$ and the projected phonon DOS onto the in-plane vibration of the Ca adsorbate, as well as the higher intensity of the peak of ${\alpha }^{2}F\left(\omega \right)$ in the low-energy region in the doped-h-BN compared with that in doped graphene. The abovementioned symmetry restriction and coupling mechanism differ from those in the Li-doped graphene and other graphite intercalation compounds: the coupling between the interlayer state and π* electrons derived via the out-of-plane vibration of the dopant [8, 11, 14]. A similar superconductivity-enhancement mechanism via symmetry lowering was theoretically predicted for black phosphorene [37]. Therefore, instead of tuning the interlayer state, reducing the symmetry and introducing strong adsorbate–substrate interactions are alternative strategies to increase the Tc of doped 2D semiconductors/insulators. In an actual system, the randomness and imperfectness of adsorption would cause additional effects, such as the emergence of localised phonon modes and the inhomogeneity of the electron density, which may affect the superconducting behavior. These issues will be investigated in future studies.

Zoom In Zoom Out Reset image size

Finally, we discuss the 3.5% tensile strain in the Ca-doped h-BN case. Experimentally, high-quality h-BN thin films have been fabricated via the exfoliation of bulk crystals or chemical vapor deposition on metal surfaces [38–41]. In the former method, a strain can be applied by deforming the polymer substrate, and experiments with several strain percentages up to 5.5% have been reported [42, 43]. However, in the latter method, the mismatch in lattice constant between h-BN and the substrate induces a certain strain in the h-BN sheet. Notably, on Au and Ag substrates, the strain induced by the lattice mismatch might be greater than 3.5% [44]. Therefore, a tensile strain of 3.5% on the h-BN monolayer might be experimentally feasible.

We investigated the possibility of superconductivity in a monolayer h-BN that was doped using group-1 (Li, Na, K) and group-2 (Be, Mg, Ca, Sr, Ba) atoms by performing ab initio calculations. We confirmed that each Ca, Sr, and Ba atom was stably adsorbed on h-BN, unlike the other atoms (Li, Na, K, Be, and Mg). In addition, in both the Ca and Sr cases, the potential-energy curve as a function of the adsorbate–surface distance exhibited two minima. In particular, in the Ca-doped case, the Ca-adsorption position at the global minimum with an equilibrium lattice constant was far from the h-BN surface; however, a tensile strain greater than 3.5% enabled the adsorption of Ca atoms at a considerably closer distance from the h-BN surface. In all the Ca, Sr, and Ba cases, electron doping made h-BN metallic. Furthermore, we calculated the Tc of both Sr- and Ba-doped h-BNs with an equilibrium lattice constant and that of Ca-doped h-BN with a tensile strain of 3.5%. The obtained Tc values were 5.83, 1.53, and 10.7 K, respectively. All the Tc values were considerably higher than those in the doped-graphene case. A significant difference between doped h-BN and graphene appeared in the k-resolved EPC projected onto the Fermi surface. In the h-BN case, the intrinsic π* orbitals formed a Fermi surface around the Γ point. The EPC in this π* orbital, mediated by the adsorbate in-plane vibration, improved owing to the strong interaction between Ca and B atoms and the lower symmetry in h-BN compared with that in graphene, thereby explaining the origin of the high Tc in Ca-doped h-BNs.

The following financial support is acknowledged: JST PRESTO JPMJPR17I7 (EM), and MEXT KAKENHI 17H05215, 19H02544 (SW and EM). The calculations were performed using the following computer facilities: Institute for Solid State Physics (Kashiwa, Japan), Information Technology Centre (University of Tokyo, Tokyo, Japan), Centre for Computational Materials Science (Tohoku University, Sendai, Japan), and Research Centre for Computational Science (Okazaki, Japan).

In this appendix, figure A1 shows the calculation results of the potential-energy curves obtained for four different adsorption configurations: N-on top, B-on top, bridge, and hollow sites. For Li, Na, K, Be, and Mg, all the adsorption sites indicate no minimum for h < 5.5 a.u.. For Sr and Ca, the hollow site is energetically favourable. Only for Ba, in addition to the hollow site, the B-on top site is a possible adsorption site. For simplicity, we focus on the hollow site herein.

Zoom In Zoom Out Reset image size

In this appendix, figure B1 shows the calculation results for the phonon dispersion of Li-, Na-, K-, Be-, and Mg-doped h-BNs, which demonstrate the dynamical instability of the adsorption structure. In all cases, the dopant was adsorbed at the hollow site.

Zoom In Zoom Out Reset image size

In this appendix, we summarise the structure parameters used in this study for the hollow-site adsorption case in tables C1 and C2. The adsorption energies for the Ca-doped h-BN with 3.5% strain, Sr-doped h-BN, and Ba-doped h-BN with the equilibrium constant were 1.14, 1.41, and 1.70 eV, respectively.