Superconductivity of bilayer phosphorene under interlayer compression
Huang Gui-Qin1, Xing Zhong-Wen2, †,
Department of Physics, Nanjing Normal University, Nanjing 210023, China
National Laboratory of Solid State Microstructures and Department of Materials Science and Engineering, Nanjing University, Nanjing 210093, China


† Corresponding author. E-mail:

Project supported by the State Key Program for Basic Researches of China (Grant No. 2014CB921103) and the Natural Science Foundation of Jiangsu Province, China (Grant Nos. BK20141441 and BK2010012).


According to first-principles calculations, it is our prediction that bilayer phosphorene (BLP) will become a quasi-two-dimensional superconductor under a certain degree of interlayer compression. A decreasing interlayer distance may realize the transition in the BLP from a semiconducting phase to a metallic phase. On the other hand, a severe vertical compression may make the BLP lattice become dynamically unstable. It is found that in the stable metallic phase of the BLP, interlayer phonon modes dominate the electron-phonon coupling λ. The obtained λ can be greater than 1 and the superconducting temperature Tc can be higher than 10 K.

1. Introduction

Initiated by the isolation of graphene,[1] research on two-dimensional (2D) materials has aroused wide interests due to the potential for significant and wide-ranging novel applications. By exfoliating a layered crystal into single layers, many 2D materials[24] have been obtained. The structure of black phosphorus (BP) consists of puckered double layers, which are held together by weak van der Waals (vdW) force.[5] Phosphorene (few-layer BP) was a newly found 2D crystal by the exfoliation from a bulk crystal and used to create field-effect transistors.[68] Large current on-off ratios, high mobilities up to 3900 cm2/V·s–4000 cm2/V·s, and quantum oscillations have been demonstrated in few-layer phosphorene devices.[9,10] It was reported that the phosphorenes possess superior mechanical flexibility and can sustain strain up to 30% and even 32%.[11] First-principles calculations suggested that strain could induce a direct-indirect band gap semiconductor transition or a semiconductor–semimetal transition in phosphorene.[1114] Furthermore, Tan et al.[15] showed that the doped phosphorenes with carbon, oxygen, and sulfur are ferromagnetism, which may play an important role in developing spin electric devices.

The investigation about the superconductivity of 2D materials[16,17] has stirred much interest due to possible applications, such as nanoscale superconducting transistors and nanoscale superconducting quantum interference devices. It was reported early that the single crystal of bulk BP showed superconductivity with Tc higher than 10 K under high pressures.[18,19] Recently, the monolayer phosphorene by a theoretically assumed electron doping with a compensating uniform positive background was suggested[20] to be a potential monolayer superconductor by first-principles simulations. Very recently, we predicted that the superconductivity can be induced in bilayer phosphorene (BLP) by Li intercalation, and the transition temperature can be increased up to 16.5 K by first-principles calculations.[21] The early experiment[22] showed that with increasing pressure, the distance between puckered layers decreases much faster than the in-plane lattice constants, for the layered BP is a strong anisotropic structure. As a result, it is interesting to see if few-layers BP can become superconducting under vertical compression (i.e., decreasing the interlayer distance). The BLP is the thinnest few-layer BP system that can be used to study the effect of changing interlayer distance.

In this work, electronic structure, lattice dynamics and electron-phonon (EP) coupling of BLP systems under vertical stress are studied via first-principles density functional theory (DFT). The calculated results show that the vertical stress on the BLP can induce a transition from a semiconductor to a BCS superconductor with the superconducting temperature Tc up to 10 K. The interlayer phonon modes play an important role in enhancing the EP coupling of BLP with decreasing the interlayer distance.

2. Computational details

The calculations are performed in a plane-wave pseudopotential representation through the PWSCF program of the Quantum-ESPRESSO distribution.[23] The ultrasoft pseudopotentials are used to model the electron-ion interactions. The valence electronic wave functions are expanded in a plane-wave basis set with a kinetic energy cutoff of 30 Ry (1 Ry = 13.6056923 eV). The general gradient approximation (GGA-PBE) is used for the exchange and correlation energy functional. The BLP is modeled using a series of slabs, which are separated from each other by a vacuum layer of thickness of about 35 Å. The vdW correction proposed by Grimme (DFT-D2) is included to model the interaction between the two layers of BLP. For the electronic-structure calculations, the Brillouin-zone integrations are performed by using the Gaussian smearing technique with a width of 0.02 Ry. Within the framework of the linear-response theory, the dynamical matrixes and the electron-phonon interaction coefficients are calculated.

3. Results and discussion

BLP consists of two puckered layers, which are interacted by weak vdW forces and are stacked following an AB stacking order. The top and side views of the atomic structure of BLP are presented in Figs. 1(a) and 1(b), respectively. Each layer may be viewed as consisting of armchair and zigzag chains along the x and y directions, respectively. For bulk BP, the early experiment[22] showed that the vdW bonds are shortened while the covalent bonds change slightly with increasing pressure. As a result, the BLP under vertical pressure is modeled by decreasing the interlayer distance. We first take the interlayer distance D as fixed, and then relax the unit cell and other atomic positions to obtain their optimized values. With the decrease of D, the calculated lattice constants, puckered distance dz, bond length and bond angles are shown in Figs. 1(c) and 1(d). It is found that lattice constant b, bond length L, and bond angle θ1 within the layer hardly change with D, reflecting the rigidity of the strong covalent bonding along the zigzag direction. For lattice constant a and bond angle θ2, however, they show different responses to the vertical compression due to the anisotropic geometric structure. For 0.75D0DD0, the lattice constant a and the bond angle θ2 hardly change with D, but for D ≤ 0.75D0, they gradually increase with decreasing the interlayer distance. The puckered distance dz shows little change, indicating that the puckered character of BLP remains basically unchanged in the range of interlayer distance studied here.

Fig. 1. Top and side views of the atomic structure for BLP [(a) and (b)], the variation of structural parameters of BLP under interlayer compression [(c) and (d)].

The calculated band structures of BLP for several interlayer distances are shown in Fig. 2. For the BLP, it has a puckered structure and its bonding is formed via sp3 hybridization. The top of the valence band has π-like character, while the bottom of the conduction band has π*-like character. For the free BLP (Fig. 2(a)), the maximum of the valence band and the minimum of the conduction band are both at the point, so that it is a direct band gap semiconductor with a band gap of about 0.4 eV. Under vertical pressure, the BLP soon becomes an indirect band gap semiconductor, as shown in Fig. 2(b) with D = 0.95D0. The minimum of the conduction band moves from point to a point between points Ȳ and , exhibiting the character of an indirect band gap. With decreasing D, the indirect band gap decreases gradually and it goes to zero at D = 0.85D0 (see Fig. 2(c)). A further decrease in D will result in the transition from a semiconductor to a semimetal (see Fig. 2(d)). The density of states (DOS) at the Fermi level (N(EF) listed in Table 1) increases with the decrease of interlayer distance D.

Fig. 2. Calculated electronic band structures of BLP under interlayer compression.

It is well known that DFT underestimates the band gap of semiconductors. Other methods incorporated many-body corrections, such as hybrid functionals,[24] GW approximation[25] seem to correct the systematic DFT underestimation. We use the hybrid functional Gau-PBE[26] to calculate the electronic structure of free BLP, the band gap is estimated to be about 1.0 eV, which is larger than the band gap in the conventional DFT. Although the hybrid functional approach results in an appreciably larger band gap, the observed trend of electronic structure under vertical compression with DFT investigations above is in line with the hybrid functional approach.

Table 1.

The calculated DOS at Fermi level N(EF), EP coupling constant λ, the logarithmically averaged frequency ωln, and the estimated Tc for BLP at three different interlayer distances.


In order to study lattice dynamical stability of the metallic BLP, we calculate its phonon spectra. The calculated phonon dispersions and the phonon DOS for several D are shown in Fig. 3. Roughly speaking, the acoustic modes and the interlayer modes for opposite vibrations of two puckered layers are in the low-frequency range, the sublayer modes for opposite vibrations of two sublayers for each puckered layer are in the intermediate-frequency range, and the P–P bond-stretching modes are in the high-frequency range due to the strong covalent bond. For D = 0.8D0, 0.75D0, and 0.7D0, their phonon spectra and phonon DOS are very similar, reflecting that the change of vdW force has only a little effect on the lattice dynamics. For D = 0.65D0, however, there appear imaginary frequencies in the two acoustic branches, as shown by the dotted lines of Fig. 3(d). For such a small interlayer distance, there may appear covalent-bond coupling between the puckered layers, which would destroy the stability of the puckered bilayer structure. If D is smaller than 0.65D0, the bilayer structure will become more unstable. As a result, the range of D for the stable metallic BLP should be within the range of 0.65D0 < D < 0.85D0.

Fig. 3. Calculated phonon dispersions, phonon densities of states, and the Eliashberg spectral function for BLP under interlayer compression.

In what follows we discuss the EP interaction. According to the Migdal–Eliashberg theory,[26] the Eliashberg spectral function α2F(ω) is given by

where N(EF) is the electronic DOS at Fermi level, is the EP matrix element, which can be determined self-consistently by the linear response theory.[23] By plotting α2F(ω), we can estimate the relative strength of the EP coupling. The frequency-dependent EP coupling is given by

Within the range of stable metallic BLP, the calculated Eliashberg spectral function α2F(ω) and frequency-dependent EP coupling λ(ω) are shown in Fig. 3 for three different interlayer distances. One can see that their phonon structures are very similar, and the larger EP coupling for smaller interlayer distance may mainly arise from the enhancement of metallicity. Furthermore, for smaller interlayer distance, we notice that α2F(ω) has a significantly enhanced peak at the frequency below 55 cm−1, which will dominate the EP coupling. For example, at D = 0.7D0, λ (55 cm−1) ≃ 1.08 is beyond 70% of the total EP coupling λ(∞)≃ 1.46. The characteristic phonon modes which dominate the EP coupling are mainly from the interlayer modes. According to the previous first-principles calculation by Du et al.,[27] the interlayer interaction is the vdW Keesom force, and the permanent dipole is formed within each puckered layer. It is speculated that this interlayer attractive interaction plays an important role in enhancing the EP coupling of BLP for smaller interlayer distance.

For three different D, the obtained EP coupling constant λ and logarithmically averaged frequency ωln, as well as the estimated Tc with μ* = 0.1 taken are summarized in Table 1. Most remarkably, one sees that λ increases dramatically with the decrease of D. The increase of both N(EF) and λ is favorable to enhancing Tc. When D = 0.7D0, λ reaches as high as 1.46 and the estimated Tc is about 10 K. The superconductivity of bulk BP was reported early under high pressure.[18,19] Our results suggest that BLP can become a BCS superconductor under a certain degree of vertical compression. Due to the decrease of the van der Waals gap, the conduction band and the valence band are overlapped, resulting in the semiconductor-metal transition. In our previous work,[21] we predicted that the superconductivity could be induced in BLP by Li intercalation, and the transition temperature could be increased up to 16.5 K. The intercalated Li atoms act as a charge reservoir of the BLP and the charge transfer from Li to P atoms leads to a semiconducting to metallic transition. The prediction of superconductivity in BLP induced by vertical compression or by Li intercalation is an interesting subject for future experimental study.

4. Summary

In conclusion, we predict that the BLP can become superconducting by adjusting the interlayer distance. A decreasing interlayer distance may realize the transition in the BLP from a semiconducting phase to a metallic phase. Our calculations of phonon structures show that the stable metallic BLP under vertical compression is within the range of 0.65D0 < D < 0.85D0. The interlayer phonon modes dominate the EP coupling, and so the interlayer vdW attractive interaction may play an important role in enhancing the EP coupling of BLP at smaller interlayer distance. The obtained EP coupling λ can be greater than 1 and the superconducting temperature Tc can be higher than 10K, suggesting that the phosphorene is a good candidate for low-dimensional superconductor.

1Novoselov KGeim A KMorozov SJiang DZhang YDubonos SGrigorieva IFirsov A 2004 Science 306 666
2Zeng FZhang W BTang B Y 2015 Chin. Phys. B 24 097103
3Wang Y YQuhe R GYu D P J 2015 Chin. Phys. B 24 087201
4Feng B JLi W BQiu J LCheng PChen LWu K H 2015 Chin. Phys. Lett. 32 037302
5Hultgren RGingrich N SWarren B E 1935 J. Chem. Phys. 3 351
6Li L KYu Y JYe G JGe Q QOu X DWu HFeng D LChen X HZhang Y B 2014 Nat. Nanotech. 9 372
7Liu HNeal A TZhu ZLuo ZXu X FTománek DYe P D 2014 ACS Nano 8 4033
8Koenig S PDoganov R ASchmidt HCastro Neto A HÖzyilmaz B 2014 Appl. Phys. Lett. 104 103106
9Li L KYe G JTran VyFei R XChen G RWang H CWang JWatanabe KTaniguchi TYang LChen X HZhang Y B 2015 Nat. Nanotechnol. 10 608
10Gillgren NWickramaratne DShi Y MEspiritu TYang J WHu JWei JLiu XMao Z QWatanabe KTaniguchi TBockrath MBarlas YLake R KLau C N 2015 2D Materials 2 011001
11Wei QPeng X 2014 Appl. Phys. Lett. 104 251915
12Hu THan YDong J M 2014 Nanotechnology 25 455703
13Rodin A SCarvalho ACastro Neto A H 2014 Phys. Rev. Lett. 112 176801
14Huang G QXing Z W 2015 J. Phys.: Condens. Matter 27 175006
15Tan X YWang J HZhu Y YZuo A YJin K X2014Acta Phys. Sin.63207301(in Chinese)
16Xing YWang J 2015 Chin. Phys. B 24 117404
17Zhang W HSun YZhang J Set al. 2014 Chin. Phys. Lett. 31 017401
18Kawamura HShirotani ITachikawa K 1984 Solid State Commun. 49 879
19Wittig JMatthias B T 1968 Science 160 994
20Shao D FLu W JLv H YSun Y P 2014 Europhys. Lett. 108 67004
21Huang G QXing Z WXing D Y 2015 Appl. Phys. Lett. 106 113107
22Chang K JCohen Marvin L 1986 Phys. Rev. B 33 6177
23Giannozzi P 2009 J. Phys.: Condens. Matter 21 395502
24Chai J DHead-Gordon M 2008 Phys. Chem. Chem. Phys. 10 6615
25Faleev S Vvan Schilfgaarde MKotani T 2004 Phys. Rev. Lett. 93 126406
26Allen P BMitrovic B1982Solid State PhysicsNew YorkAcademic371
27Du Y LOuyang C YShi S QLei M S 2010 J. Appl. Phys. 107 093718