Prediction of nodal-line semimetals in two-dimensional black phosphorous films

Semimetals are a new kind of quantum materials, in which the conduction and valence bands cross each other near the Fermi level. Based on density-functional theory calculations and symmetry analysis, we propose nodal-line semimetals in layered stacked black phosphorus (BP) films which are designed to have a mirror symmetry lying in the BP layer plane and thus rendering them different from the BP film systems previously studied. A closed nodal-line degenerate band can appear around the Fermi level in the BP films after a biaxial compressive strain is applied. The calculated Z2 number of Z2 =  − 1 indicates the robustness of the nodal-line semimetals obtained in the BP films, protected by the in-plane mirror symmetry. Intriguingly, with the increase of the film thickness, a smaller biaxial compressive strain is required to produce the nodal-line semimetals, more accessible in experiments. Our results provide a promising route to carrying out the nodal-line semimetals based on various two-dimensional materials.

Black phosphorus (BP) is a layered structure and has potential applications in nanoelectronics owing to its unique behaviors like high carrier mobility, large anisotropy, and negative Poisson's ratio etc [21][22][23][24][25] . 3D BPs, the most stable allotrope of phosphorus, were successfully synthesized in 1914 26 with a layered structure. Each layer of the BP is a 2D hexagonal and puckered along the armchair direction. The specially puckered atomic structure of the BP can result in the negative Poisson's ratio in the BP film 22 . Due to the weak interlayer interactions, monolayer or few-layer BPs can be cleaved from the 3D bulk black phosphorus. Monolayer BP flakes with the thickness of ~ 0.7 nm were exfoliated by Wang et al. in 2015 23 . Based on BP films, field-effect transistors with a high mobility of ~ 1000 cm 2 /V/s were fabricated in experiments 27 . Dirac semimetals were reported in experiments in BPs by applying an externally vertical electric field 28 . In theories, NLSs 29 were predicted in 3D BPs by applying hydrostatic pressure etc. For 2D BP films, even though there are numerous previous researches, most of them have focused on digging out the Dirac semimetals 30,31 and NLSs have not been explored in 2D BP films yet. For future applications, it is meaningful to explore various semimetal behaviors in 2D black phosphorus due to the need for the integration and miniaturization of the devices in nanoelectronics.
In this work, based on density function theory (DFT) calculations and symmetry analysis, we obtain NLSs in 2D layered BP ultrathin films with different thickness for the first time. For a bilayer BP built with a mirror symmetry, an elliptic closed nodal-line degenerate band is found being located around the Fermi level (E F ) if a

Results and discussion
We focus on N-layer BP films (N = 2, 4, 6) to produce the 2D NLS. Since a mirror symmetry possibly leads to nodal-line degenerate bands 9,10 , we built a bilayer BP structure with a mirror symmetry as shown in the right panels of Fig. 1a. The weak interlayer interactions can lead to various stacking orders for the bilayer BPs, which were studied by Zhang et al. 32 and Dai et al. 33 . The bilayer BP with the same stacking order as the 3D bulk BP has been found owning a quantum spin Hall state after certain controlling schemes are applied 32 . As shown in Fig. 1a, the bilayer BP we studied can be obtained from the bilayer BP, in which the top layer stacks vertically on the bottom layer (in the left panels of Fig. 1a), by shifting the 1 and 2 atoms of the top layer to the positions right above the 3′ and 4′ atoms of the bottom layer, respectively. The designed bilayer BP structure (right panels in Fig. 1a) has an orthorhombic lattice with a layer group 41 (space group 51, Pmma) 34 . The M z mirror symmetry owned in the structure is illustrated in Fig. 1b. The structures of four-and six-layer BPs can be similarly built, with the same layer group and M z mirror symmetry. The first Brillouin zone (BZ) and the projected one dimensional (1D) BZ of the bilayer BP film built are displayed in Fig. 1c.
The calculated structural parameters of the N-layer BP films (N = 2, 4, 6) are listed in Table 1. The parameters of the bilayer BP with a 5.5% biaxial compressive strain are also given. For the pristine bilayer BP, the optimized lattice constants are a 1 = 3.298 Å, a 2 = 4.628 Å. As shown in Table 1, with the increase of the layer thickness N, a 1 increases slightly and a 2 decreases gradually. R 1 and R 2 do not change much with the variation of the layer thickness. The interlayer distance D, however, increases obviously with the increase of N, showing the weakening of the interlayer interactions between this two BP layers. The applied biaxial compressive strain does not affect the symmetries the system owns. Thus, the crystal with the biaxial strain applied still possesses the mirror symmetry (M z ). For the bilayer BP with a 5.5% compressive strain, the distance D between the two BP monolayers decreases much compared to that of the pristine bilayer BP, ascribed to the obvious decrease of the angle θ between the R 1 and R 2 bonds under the compressive strain. To study the structural stability of the built bilayer BP films, we construct a 4 × 4 × 1 supercell and calculate the frequency dispersion of the bilayer BP along high  Table 1. The optimized lattice constants (a 1 and a 2 ), the bond lengths (R 1 and R 2 ), the angle between the R 1 and R 2 bonds (θ), and the distance (D) of the layer interval having the mirror symmetry for the multiple layer BP films (N = 2, 4, 6). The second line gives the results of the bilayer BP with a 5.5% biaxial compressive strain. www.nature.com/scientificreports/ symmetry lines in the BZ. As shown in Fig. 2a, no imaginary frequencies appear, indicating that the bilayer BP is dynamically stable at ambient pressure. The band structure of the pristine bilayer BP is presented in Fig. 2b. The band dispersions along the Г-X and Г-Y directions are so different, indicating the large anisotropy in the BP 21 . There is a direct band energy gap of 0.57 eV between the valence band maximum (VBM) and conduction band minimum (CBM) near the Г point. The band gap as a function of the biaxial compressive strain is shown in Fig. 2c. The strain is defined as (l' − l) × 100% /l, where l' stands for the lattice constant with a strain applied and l stands for the relaxed lattice constant without the strain. With the increase of the biaxial compressive strain, the band gap decreases gradually. When a 5.2% compressive strain is applied, the band gap is closed. Note that the band gap is defined as E B1 -E B2 , where E B1 and E B2 are the minimum and the maximum energies of the B 1 and B 2 bands, respectively (see Fig. 2b). Thus, the band gap becomes negative when the compressive strain is larger than 5.2%. Hence, the bilayer BP undergoes a phase transition from a semiconductor (SC) to a metal or semimetal after a compressive strain with a certain magnitude is applied. Besides the strain, applying a pressure is also an effective tuning tactics in experiments. With the extraction method mentioned in Ref. 30 , the band gap as a function of the external pressure is calculated. As displayed in Fig. 2c, when the external pressure reaches 11.2 GPa, the band gap becomes zero and the phase transition occurs.
The mechanism of above phase transition is now analyzed. The band gap of the monolayer BP is about 1.51 eV 21 , while the bilayer BP gives a much smaller band gap (0.57 eV) (Fig. 2b), showing that the interlayer interactions between the two BP monolayers tend to reduce the band gap. For few-layer BP films, the interlayer interactions contain two aspects: the traditional van der Waals (vdW) interaction and the electronic hybridization 35 . The latter comes from the lone electron-pairs of P atoms due to the three sp 3 hybridization orbitals formed on each P atom. This extra interaction can lead to some exceptional behaviors in the BP film, such as the thicknessdependent vibrational properties etc., not owned by graphene films 35 . Thus, the interlayer electronic hybridization and the vdW interaction together influence the electronic states when the interlayer distance D is varied. When a biaxial compressive strain is applied, the D as a function of the strain is displayed in Fig. 2d. The greater the strain strength is, the smaller the distance D is. Thus, the interlayer interactions are strengthened, in favor of closing the band gap. On the other hand, the bands around the E F in the bilayer BP are primarily made up of the p z orbitals (see Fig. 3a). Therefore, the band gap of the bilayer BP must be formed due to the bonding and anti-bonding states of the p z orbitals, similar to the case in the monolayer BP 24 . As listed in Table 1, when a 5.5% compressive strain is applied, the bond length R 2 decreases slightly (0.001 Å), indicating the p z bond between the two P atoms (such as the 2 and 3 (or 1′ and 4′) atoms in Fig. 1a) in one BP monolayer becomes strong. This trend results in a large band gap between the bonding and anti-bonding states. The R 2 change is, however, very www.nature.com/scientificreports/ small with the increase of the compressive strain (Table 1). Hence, the enlarging effect of the band gap from the decrease of R 2 is negligible, compared to the reduction effect of the band gap from the D decrease in the process. Hence, the main cause of the phase transition from a SC to a metal or semimetal (Fig. 2c) is the increase of the interlayer interactions in the bilayer BP, as illustrated in the inset of Fig. 2d. As displayed in Fig. 2c, when the compressive strain is greater than 5.2%, the band gap of the system is closed. To show the bilayer BP entering the NLS state, we take the bilayer BP with a 5.5% biaxial compressive strain as an example and analyze concretely its electronic states in following. Its structural stability is first explored. Figure 2e gives the phonon spectrum of the bilayer BP under the 5.5% compressive strain. Similar to the case without the strain (Fig. 2a), no imaginary frequencies appear in Fig. 2e, inferring the dynamical stability of the structure under the compressive strain. Besides, the thermal stability of the bilayer BP under the 5.5% compressive strain is also investigated by an ab initio molecular dynamics (MD) simulation. A 4 × 4 × 1 supercell is adopted in the calculations. As illustrated in Fig. 2f, no structural reconstruction happens and the structure remains intact at 300 K for 12 ps, indicating the thermal stability of the bilayer BP under a large compressive strain.
The orbitally resolved band structure for the bilayer BP under the 5.5% compressive strain is displayed in Fig. 3a. The result clearly reveals that a crossing of the conduction and valence bands, composed of the p z orbitals, appears around the E F . This band inversion may lead to interesting band structures. To check whether the bilayer BP forms a NLS near the E F under a 5.5% compressive strain, we plot its 3D band dispersion diagram and the corresponding Fermi surface. As illustrated in Fig. 3b, the highest occupied band and the lowest unoccupied band form a circle of band degenerate intersection near the Г point in the BZ. The Fermi surface plotted in the Fig. 3c also displays a nodal loop around the Г point. Thus, a NLS is formed in the bilayer BP with a biaxial compressive strain larger than 5.2%. We check the electronic structures of the bilayer BP film with hybrid Heyd-Scuseria-Emzerhof (HSE06) functional 36 since the GGA functional may underestimate the band gap. The band gap obtained from the HSE06 is about 1.12 eV, larger than that (0.57 eV) from the GGA. And a larger biaxial compressive strain (> 5.2%) is expected to produce the NLS in the system. Fortunately, as illustrated in Fig. 2c, the phase transition can be induced by applying an external pressure with a not very large magnitude (11.2 GPa). Even if the HSE06 functional is employed, the critical pressure increases to 16.4 GPa, still easily carried out in experiments. Besides, as discussed below, the band gap decreases with the increase of the BP thickness. Thus, the NLS may be easily achieved in experiments in the thicker BP films. Due to the weak spin-orbit coupling (SOC) strength of the P atoms, the SOC effect is ignored in above discussion. If the SOC is considered, all the 2D crossing points do not exist anymore. Small band gaps appear along the nodal loop, consistent with the trend from the symmetry analysis that the two inverted bands have the same M z eigenvalues instead of different ones. The band gaps opened along the Г-X and Г-Y directions are of 1.5 and 6.6 meV, respectively, also reflecting the large anisotropic behaviors of the BP films. Since the band gaps are very small, they are neglected in the following. www.nature.com/scientificreports/ To characterize the 2D NLS obtained, we employ mirror eigenvalues of the M z to define the Z 2 number, as performed in Ref. 10 . The Z 2 number is calculated through 10 where ξ a = N occ n=1 ξ an ( ξ a is the mirror eigenvalue and N occ is the number of the occupied states). For the above bilayer BP system, Z 2 = −1 is obtained. Since there is only one nodal line in the BZ (Fig. 3c), we prove that Z 2 = (−1) N , where N is the total number of the nodal lines in a 2D system, i.e., N = 1 and Z 2 = −1 are obtained for the bilayer BP film. Besides, as illustrated in Fig. 3d, the different mirror eigenvalues of ±1 for the two crossing bands also indicate that the degeneracy of the two bands along the nodal loop is protected by the mirror symmetry. Therefore, the formed 2D NLS in the bilayer BP is very robust, protected by the mirror symmetry M z . The isosurfaces of the charge densities of the B 1 and B 2 bands (Fig. 3a) on the yz plane are displayed in Fig. 3e. A chemical bonding-like character from the B 1 band is obviously observed between the two P monolayers, manifesting the interlayer quasi-covalent (hybridization) interaction, which together with the traditional vdW interaction can tune effectively the band gap under the strain.
To fully understand the origin of the nodal line, we shift the A and B atoms (Fig. 4a) in the bilayer BP with a 5.5% compressive strain along the z axis by 0.01 and − 0.01 Å, respectively, to artificially break the mirror symmetry M z . Note that this bilayer BP without the M z symmetry is less stable than that of the bilayer BP with the M z symmetry (Fig. 1b) by 435 meV per unit cell. The corresponding band structure is given in Fig. 4b, from which we find that all the nodes around the Г point are annihilated and the degenerate points are all gapped after the mirror symmetry (M z ) of the lattice is broken. An indirect band gap of 64 meV is formed in the band (Fig. 4b). This evidence shows that the nodal-line semimetal obtained in the bilayer BP with a large (> 5.2%) compressive strain is protected by the mirror symmetry M z. Thus, the NLSs may appear in the 2D systems if the M z exists in the structures. If there is a random shift across the film and the mirror symmetry M z of the structure is broken, the nodal line does not exist anymore. Due to the relatively weak interlayer interactions of the BP film, a low temperature condition is expected to observe the nodal line behavior in the BP thin film (Fig. 3a).
The electronic states of the N-layer BP films with N = 4 and 6 are also explored. Since the thicker the BP film is, the stronger the interlayer interactions will be, which makes the dispersion of the valence and conduction bands larger, resulting in a smaller band gap. We, thus, expect that the thicker BP films can also possess robust NLSs as the bilayer BP, under maybe smaller biaxial compressive strains. This speculation is confirmed by our DFT results. The N-layer BP films (N = 4 and 6) also have geometrical structures with a layer group 41 Pmma (space group 51) and M z mirror symmetry (Fig. 5a). The dispersion of the low-energy bands of the 4-layer BP film is found to be similar to that of the bilayer. The band gap becomes 0.38 eV, smaller than that (0.57 eV) of the bilayer BP. This smaller band gap of the pristine 4-layer BP leads to smaller compressive strains required to close the band gap and trigger the phase transition. Our calculations show that the band gap of the 4-layer BP can be closed with a 3% biaxial compressive strain. With the further increase of the strain, the band inversion occurs near the E F and a NLS is formed in the system. The band structure of the 4-layer BP with a 3.5% biaxial compressive strain is displayed in Fig. 5b. Similar to the bilayer case, the two crossing bands around the E F also have opposite mirror eigenvalues ±1 along the high symmetry lines of the Г-X and Г-Y directions and N = 1 and Z 2 = −1 are obtained, indicating the robustness of the nodal-line semimetal achieved. For the 6-layer BP, the nodal line is found appearing under a biaxial compressive strain larger than 1.7% (see Fig. 5c). Thus, a NLS can be more easily acquired in the thicker BP films, favorable for experimental observation.

Conclusion
Based on first principles calculations and symmetry analysis, we find two-dimensional nodal-line semimetals in the built multiple layer BP films with certain biaxial compressive strains applied. The nodal-line appears around the Fermi level and comes from the inversion of the two p z bands with opposite mirror eigenvalues ±1 . The Z 2 = −1 shows the robustness of the nodal-line semimetals obtained. The thicker the BP film is, the smaller the compressive strain required to produce the phase transition from an ordinary semiconductor to a nodal-line semimetal will be. Our results provide a route to designing the two-dimensional nodal-line semimetals and promote promising applications of the black phosphorus films in future nanoelectronics.

Methods
Our first-principles calculations are performed based on density functional theory (DFT) using the projector augmented wave method as implement in the Vienna ab initio simulation package (VASP) 37 . The Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE) is adopted for the exchange-correlation functional 38 . The cutoff energy is set as 400 eV for the plane-wave basis and the Brillouin zone is sampled with k meshes of 14 × 10 × 1 by using Monkhorst-Pack method 39 . Forces on the ions are calculated according to the Hellmann-Feynman theorem. The convergence threshold for the total energy is set to 1 × 10 −5 eV. The forces on all atoms are optimized to be less than 0.01eVÅ −1 . The phonon state calculations are carried out by using the supercell approach implemented in the PHONOPY code 40 . Molecular dynamics simulations are carried out to determine the thermal stability of the bilayer BP. The maximally localized Wannier functions (MLWFs) are constructed by employing the WANNIER90 code 41 , from which the Fermi surface is performed by WannierTools package 42,43 . www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.