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^[2–4] 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.^[6–8] Large current on-off ratios, high mobilities up to 3900 cm²/V·s–4000 cm²/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.^[11–14] 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 T_c 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 T_c up to 10 K. The interlayer phonon modes play an important role in enhancing the EP coupling of BLP with decreasing the interlayer distance.

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.

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 d_z, 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 θ₁ 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 θ₂, however, they show different responses to the vertical compression due to the anisotropic geometric structure. For 0.75D₀ ≤ D ≤ D₀, the lattice constant a and the bond angle θ₂ hardly change with D, but for D ≤ 0.75D₀, they gradually increase with decreasing the interlayer distance. The puckered distance d_z shows little change, indicating that the puckered character of BLP remains basically unchanged in the range of interlayer distance studied here.

Fig. 1.

	Figure Option ViewDownloadNew Window
	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 sp³ 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.95D₀. 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.85D₀ (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(E_F) listed in Table 1) increases with the decrease of interlayer distance D.

Fig. 2.

	Figure Option ViewDownloadNew Window
	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.

Table 1.

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. . D/D0 N(EF)/(states/eV) λ ωln/K Tc/K 0.80 0.68 0.32 183.2 0.2 0.75 0.86 0.85 107.9 5.7 0.70 1.08 1.46 96.1 10.6

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.8D₀, 0.75D₀, and 0.7D₀, 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.65D₀, 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.65D₀, 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.65D₀ < D < 0.85D₀.

Fig. 3.

	Figure Option ViewDownloadNew Window
	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 α²F(ω) is given by

For three different D, the obtained EP coupling constant λ and logarithmically averaged frequency ω_ln, as well as the estimated T_c 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(E_F) and λ is favorable to enhancing T_c. When D = 0.7D₀, λ reaches as high as 1.46 and the estimated T_c 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.

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.65D₀ < D < 0.85D₀. 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 T_c can be higher than 10K, suggesting that the phosphorene is a good candidate for low-dimensional superconductor.