The strain effect on superconductivity in phosphorene: a first-principles prediction

The effects of biaxial and uniaxial strains on electron–phonon coupling and superconductivity in monolayer phosphorene are systematically investigated by first-principles calculations. It is found that the electron–phonon coupling primarily comes from the low frequency optical phonon modes around B 3 g 1 ?> , and the biaxial strain gives rise to more a obvious increase in density of states around the Fermi level and phonon softening in the low frequency regime compared to the other two types of uniaxial strain. Therefore, the electron–phonon coupling is more significantly enhanced by the biaxial strain than the uniaxial strains and the superconducting transition temperature Tc increases sharply from 3 K to 16 K at the typical doping concentration n2D = 3.0 × 1014cm−2 when the biaxial strain reaches 4.0%.


Introduction
Two-dimensional (2D) materials [1], such as graphene based materials and transition metal dichalcogenides (TMDCs), have generated intensive research interest owing to their special physical properties, some of which are not seen in their bulk counterparts. For instance, the massless Dirac fermion of graphene [2] leads to an extremely high room temperature mobility of 10 4 cm 2 V −1 s −1 and other multitudinous electronic properties, but the absence of an energy gap results in a low on/off ratio and limits its application in the mainstream semiconductor industry. In the case of semiconducting MoS 2 [3,4], though it has a moderate gap, a much lower carrier mobility than that of graphene indicates that it may not be a promising candidate for high-performance devices. Very recently, layered black phosphorus, termed phosphorene, has been experimentally isolated (mechanical exfoliation [5][6][7] and plasma-assisted fabrication [8]) and has good prospects in research owing to its tunable energy gap and remarkable electronic properties [9][10][11][12]. The layer-number-dependent direct energy gap [13][14][15][16][17], ranging from 1.51 eV in a monolayer to 0.59 eV in 5-layers [14], leads to its potential application in field-effect transistors (FETs) and photoelectric materials with a wide spectrum. In addition, previous works [6,14,18] have reported that its room temperature mobility is as high as 1000 cm 2 V −1 s −1 [6] with explicit anisotropy, which even reaches 10 4 cm 2 V −1 s −1 [14] in thetheoretical prediction obtained from deformation potential approximation.
In materials science, one important application of 2D materials is in nanoscale superconducting devices, including superconducting quantum interference devices and superconducting transistors [19][20][21]. Previous studies of bulk black phosphorus [22,23] have observed a superconducting phase in the metallic state under high pressures, when T c is above 10 K. Regarding phosphorene, there has been a first-principles calculation predicting superconductivity [24] when T c is about 4 K at a typical doping concentration of × − 3.0 10 cm 14 2 . However, a feasible experimental method to further improve T c still needs to be explored. As a common approach, strain has long been successfully used to tune electronic properties of 2D materials [25][26][27][28][29][30][31][32][33][34][35]. In phosphorene, it is found that its layered structure can sustain a wide range of strains [36], which gives rise to the semiconductor-metal transition and rotates the preferred conducting direction by 90 degrees [18,[36][37][38].
Motivated by the effects of strain on the electronic properties of phosphorene, we employ first-principles calculations in the present work to obtain further insight into the evolution of superconductivity with respect to the magnitude and type of strain. The strain-dependent electron-phonon coupling constant λ and superconducting transition temperature T c are systematically investigated in electron-doped phosphorene. It is found that the low frequency optical phonon modes around B g 3 1 play a key role in the electron-phonon coupling. Through modulating the band structure and phonon spectrum, the strains prominently enhance the electronphonon coupling and the superconductivity, especially in the case of biaxial strain (ε = 4.0%), improving the T c from 3 K to 16  ) illustrates the interaction between the electrons and phonons, expressed as [39]  where| 〉 i k is the Kohn-Sham electronic eigenstate with wavevector k and band index i, and m is the atomic mass.
is the derivative of the Kohn-Sham potential with respect to atomic displacement along the ν phonon mode with wavevector q and the corresponding phonon frequency is ω ν q, . The electron-phonon coupling strength λ ν q, for mode index ν and wavevector q is given by [39], ) where N F is the density of states at the Fermi level, N k is the total number of k points and ϵ i k is the band energy of the Bloch electron.
The dimensionless parameter electron-phonon coupling constant λ is given by the average coupling [39], q q q 0 2 , , α 2 F(ω) is the Eliashberg spectral function, written as [39], where N q is the total number of q points. Based on the isotropic approximation in the Eliashberg theory, which describes the conventional electronphonon coupling superconductors very well, we estimate T c using the Allen-Dynes-modified McMillan equation [39], where μ* is the effective Coulomb repulsion with a reasonable value of 0.1 [24,40] and ω 〈 〉 log is the logarithmically averaged phonon frequency,

Technical details
All of the first-principles calculations in the present work were carried out with the QUANTUM ESPRESSO package [41], using the norm-conserving pseudopotentials [42]. The Perdew, Burke and Ernzerhof parameterized generalized-gradient approximation (PBE-GGA) [43] was used to describe the exchangecorrelation effect. The kinetic energy cutoff of 40 Ry was used in all the calculations. The energy and force convergence criteria for fully relaxing atomic positions were 10 −5 Ry and 10 −4 Ry a.u. −1 , respectively. A Monkhorst-Pack k-mesh of 42 × 30 × 1 was used for electronic calculations. In our slab model, a vacuum layer of 15 Å was set to avoid interactions between the adjacent atomic layers. In order to be a metallic state, the electron-doped phosphorene was achieved by increasing the valence charge while simultaneously introducing the same amount of uniform background charge. The strains were introduced by adjusting the lattice constants with uniaxial strain capacity ε = − × a a a ( ) 100% and biaxial strain capacity ε ε ε = = x y , where a 0 (b 0 ) is the lattice constant along the it x (y) direction for the strained-free structure. The phonon spectrum and λ were calculated on a 14 × 10 × 1 q-grid by using the density functional perturbation theory (DFPT) [44].

Results and discussion
3.1. Band structure and phonon spectrum Unlike graphene, the lattice structure of monolayer phosphorene has a puckered structure (figure 1), because of the five valence electrons of phosphorus the atom and the sp 3 hybridization [38]. In our calculations, the inplane lattice constants a 0 and b 0 are optimized to be 3.30 Å and 4.62 Å, respectively, similar to previous works [5,36]. In our calculations for different doping concentrations, we have used fixed lattice constants of a unit cell and allowed the atomic positions to be determined by full relaxation for every calculation. As summarized in table 1, R 2 increases by 0.033 Å from = × n 0.26 10 2D 14 cm −2 to = × n 2.6 10 2D 14 cm −2 , whereas R 1 grows only by 0.002 Å. The bond angles, α and β, decrease by 0.15°and 0.26°, respectively.
In the band structure, there is an energy gap of 0.91 eV at the Γ point, which is in agreement with previous theoretical calculations [36]. We define the first and second conduction band minima at the Γ point to be Γ 1 and Γ 2 , respectively (figure 2(a)). There are other two valleys near the X point, X 1 and X 2 , with close energies. When increasing the doping concentration, Γ 1 slightly drops relative to X 1 and in the meantime one band falls from above to somewhere between Γ 1 and Γ 2 in the high doping level (figures 2(b) and (c)). Similar behavior-the shift of the conduction band with doping-has also been found in other layered materials [45][46][47]. Monolayer phosphorene has twelve branches of phonons and belongs to the non-symmorphic space group D . h The eight irreducible representations at the Γ point help to denote the vibration modes and figure 2(e) shows the three main optical vibration modes in the following discussions. Comparing the phonon spectra at different doping concentrations (figure 2(d)), one finds that the optical phonons around the Γ point are greatly changed with doping. In the low frequency branches of the phonon spectrum, the it y direction stretching optical mode [24,48] B g 3 1 softens evidently and other low frequency modes are nearly invariable. Another remarkable softening appears in the high frequency branches, which are the outplane optical modes [24,48] stretching A 1 g and interlacing B g 3 2 . The significant phonon softening in the doping condition may be attributed to the changes of bond lengths (table 1), the Kohn effect, and the nesting effect around the Fermi surface [24], which are also found in graphene and MoS 2 [35,49].
When strains are turned on, they show significantly different effects on the different electronic eigenstates owing to their different bonding natures [36]. The application of biaxial strain ε brings Γ 1 and Γ 2 close to each other and also leads to the apparent decline of the conduction band whose bottom is along the Γ-Y line. As a  result, the five conduction band bottoms are simultaneously located in a small energy window around the Fermi level (figures 3(a)-(c)). Similar behaviors are also observed in the case of strain along the zigzag (x) direction (ε x = 4.0%) ( figure 3(d)), but with a greater energy difference between Γ 1 and Γ 2 than under of biaxial strain (ε = 4.0%) ( figure 3(c)). The effect of strain along the armchair (y) direction on the conduction band (figure 3(e)) is far less obvious than the other two types of strain.
In the effect of strains on the phonon spectrum, the biaxial strain weakens the bond energy which, overall, lowers the atomic vibration frequencies, as shown in figures 4(a) and (b). In particular, the B g 3 1 mode exhibits conspicuous softening, which plays a crucial role in enhancing the superconducting T c via equation 2. In the case of uniaxial strains (figures 4(c) and (d)), both uniaxial strains affect the phonon spectrum in the low frequency branches less than the biaxial strain. In addition, due to the direction of vibration of the B g 3 1 mode, i.e. the y direction (figure 2(e)), the strain along the armchair (y) direction has more influence on the B g 3 1 mode than that of the zigzag direction. By comparison, in the high frequency branches of the phonon spectrum there is more obvious phonon softening under strain along the zigzag direction.

Electron-phonon coupling
The above obvious changes of electronic band structure and phonon spectra indicate a strong dependence of electron-phonon coupling on the strains. To further substantiate this statement, we have calculated the Eliashberg spectral function α ω F ( ) 2 and the electron-phonon coupling constant λ under different strains (figure 5).
As shown in figure 5(a), the α ω F ( ) 2 spectrum spreads through all phonon frequencies with three major parts lying on different frequency ranges. The stronger peaks of the spectrum lie in the low frequency regime around 150 cm −1 and primarily come from the coupling between the electrons and phonons, with frequencies around that of the B g  The doping-dependent total electron-phonon coupling constant λ under the strains is shown in figure 5(c) and (d). It is usually expected that the increase of doping concentration strengthens the electron-phonon coupling as verified by the overall increasing trend in figure 5(c), but a local decline of λ appears in a small range around = × n 1.5 10 2D 14 cm −2 under the low strains ε ⩽ 2.0% (we shall explain this in the next paragraph). It is important that the biaxial strain (ε = 4.0%) markedly enhances the electron-phonon coupling, improving λ from 0.4 to 1.6 at = × n 3.0 10 2D 14 cm −2 . In the case of uniaxial strains ( figure 5(d)), due to the obvious changes of electronic band structure and the high frequency branches of the phonon spectrum under the strain along the zigzag direction, it has more influence on the electron-phonon coupling than the uniaxial strain along the armchair direction. But both effects of the uniaxial strains are far less obvious than those of the biaxial strain. It is  also found that the increase of strain does not always enhance the electron-phonon coupling. As shown in figure 5(d), the electron-phonon coupling of ε x = 5.0% is slightly below that of ε x = 4.0%.
In order to more clearly decipher the relative contributions of phonons with different wave vectors to the electron-phonon coupling, we introduce the distribution map of the wave-vector-resolved electron-phonon coupling parameter λ λ = ∑ ν ν q q , , as plotted in figure 6. We define the region around Γ point as I and its adjacent highlight region as II ( figure 6(b)). In the strain-free condition, the strong electron-phonon coupling in area I is mainly attributed to the appearance of X 1 and X 2 near the Fermi level (figures 2(b)), while the electron-phonon coupling in area II (figure 6(a)) mainly originates from the electrons on the lowest conduction band around the Γ point. With the increase of doping concentration, the lower of the conduction band minima, Γ 1 , and the band which falls from above to somewhere between Γ 1 and Γ 2 (figure 2) lessen the electron population on the X 1 and X 2 valleys. This significantly weakens the electron-phonon coupling in area I but hardly affects it in area II (figure 6(c)), so that the decline of λ appears around = × n 1.5 10 2D 14 cm −2 . The application of biaxial strains obviously strengthens the electron-phonon coupling in area I, while in area II it is enhanced slightly (figures 6(d) and (e)). The reason is that the biaxial strain, on the one hand, impacts the conduction band to increase the density of states near the Fermi level (figures 3(b) and (c)) and, on the other hand, clearly softens the phonons in the low frequency regime (figures 4(a) and (b)). The two effects are advantageous for strengthening the electron-phonon coupling. Based on the study for the changes of electronphonon coupling caused by different phonon modes, one can mainly attribute the enhancement of the electronphonon coupling in area I to the increase of the contribution from the low frequency optical phonon modes around B g 3 1 , which sharply soften with increasing biaxial strain. In the case of uniaxial strains (figures 6(f)-(g)), the increase of electron-phonon coupling in area I is far less significant than in the case of biaxial strains because of the weaker phonon softening in the low frequency regime under uniaxial strain. In addition, the effect of strain along the armchair direction on the electron-phonon coupling in area I is slightly stronger than that along the zigzag direction, consistent with the softening of the B g 3 1 phonon mode above (figures 4(c) and (d)). When increasing the strain along the zigzag direction to 5% (figures 6(h)), the electron-phonon coupling is weaker than in the case of ε x = 4.0%, because X 1 and X 2 are far away from the Fermi level and lower the density of states (figure 3(f)).

Superconductivity
The superconducting transition temperature T c under different strains was estimated by using equation 5 with a reasonable μ* of 0.1 [24,40], as shown in figure 7. Table 2 also summarizes The results of T c in a reasonable range of μ* (e.g., 0.1 ∼ 0.15). In the strain-free condition, we got T c = 3 K with λ = 0.4 at the typical doping concentration 3.0 × 10 14 cm −2 . More importantly, T c can be enhanced to 16 K at the same doping concentration with a biaxial strain of 4.0%. However, the uniaxial strain can only enhance the T c to 10 K (8 K) for ε x (ε y ) = 4.0%   3.0 10 2D 14 cm −2 , which is below the value in the case of biaxial strain. For this reason, the uniaxial strains cannot bring about the enhancement of density of states around the Fermi level and obvious phonon softening simultaneously, which are the major factors in enhancing the T c of phonon-mediated superconductivity. Compared to graphene [50], there are both similar and different aspects of the enhancing effect of biaxial tensile strain on superconductivity. In both systems, strain gives rise to the softening of the phonon mode, while in graphene the mode is the shear horizontal optical in-plane C-C stretching mode, in phosphorene it is the y direction stretching optical mode. However, the strain in graphene has little effect on the conduction band, which is in contrast to phosphorene. According to the enhancement of T c under strain, our results effectively serve as guidance for experimentalists to improve the superconductivity of phosphorene.

Summary
In the present work, we investigate the effects of strain on electron-phonon coupling and superconductivity in monolayer phosphorene using first-principles calculations. We found that the low frequency optical phonon modes around B g 3 1 mainly contribute to electron-phonon coupling. Our important discovery is that the biaxial strain enhances the electron-phonon coupling significantly and hence raises the superconducting transition temperature from 3 K to 16 K at the typical doping concentration × 3.0 10 14 cm −2 used here. The boosting of T c is attributed to the simultaneous increases of the density of states around the Fermi level and the softening of the phonons. Furthermore, the anisotropy in phosphorene leads to different effects on the electron-phonon coupling under uniaxial strain along different directions. Unlike the high T c in the case of biaxial strain, the uniaxial strains of ε x = 4.0% and ε y = 4.0% enhance T c to 10 K and 8 K, respectively, at a doping concentration of 3.0 × 10 14 cm −2 .