Two-Dimensional Pnictogen Honeycomb Lattice: Structure, On-Site Spin-Orbit Coupling and Spin Polarization

Because of its novel physical properties, two-dimensional materials have attracted great attention. From first-principle calculations and vibration frequencies analysis, we predict a new family of two-dimensional materials based on the idea of octet stability: honeycomb lattices of pnictogens (N, P, As, Sb, Bi). The buckled structures of materials come from the sp3 hybridization. These materials have indirect band gap ranging from 0.43 eV to 3.7 eV. From the analysis of projected density of states, we argue that the s and p orbitals together are sufficient to describe the electronic structure under tight-binding model, and the tight-binding parameters are obtained by fitting the band structures to first-principle results. Surprisingly large on-site spin-orbit coupling is found for all the pnictogen lattices except nitrogen. Investigation on the electronic structures of both zigzag and armchair nanoribbons reveals the possible existence of spin-polarized ferromagnetic edge states in some cases, which are rare in one-dimensional systems. These edge states and magnetism may exist under the condition of high vacuum and low temperature. This new family of materials would have promising applications in electronics, optics, sensors, and solar cells.

Exploration of unknown phases of materials has been the scientific endeavor for the past decades, and the discovery of new materials can spurt new field of study for both experimentalists and theorists. The low dimensional systems are gaining increasing attention in recent years. Low dimensional systems generally refer to systems in zero-, one-or two-dimension. In particular, various two-dimensional (2D) materials are found in recent years 1,2 . The novel properties of 2D materials make them important in both fundamental research and applications. For example, graphene has remarkable mechanical, thermal, and electronic properties. 2D transition metal dichalcogenides have peculiar optical properties 3,4 . They have promising applications in optoelectronics, spintronics, catalysts, chemical and biological sensors, supercapacitors, and solar cells 5 . In other words, the study of 2D materials is of great importance.
Successful exfoliation of graphene marked the beginning of the study of 2D materials [6][7][8] . However, the application of graphene is limited by the lack of intrinsic gap. The planer structure of graphene also renders it impossible to control the electronic structure by perpendicular electric field. Recently, black phosphorus was synthesized and studied [9][10][11] , which has a direct band gap and high carrier mobility 10,11 . Its puckered structure also enables the modification of band structure by perpendicular field. Recently, another allotrope of phosphorus called blue phosphorus was predicted. It has buckled structure similar as silicene and germanene, and its band gap is close to black phosphorus 12 . Afterwards, two additional phosphorus allotropes were predicted by first-principle study 13 . Accordingly, the phosphorus allotropes show great importance in both theory and experiment. Phosphorus is pnictogen, and belongs to the nitrogen group (VA). Besides the study of phosphorus, 2D systems involving other VA elements are also carried out. Quantum spin hall effect is predicted to exist in honeycomb X-hydride/halide (X = N-Bi) monolayers 14 and BiX/SbX (X = H, F, Cl and Br) monolayers 15 . Arsenic, antimony and bismuth honeycomb monolayers are also studied in topological aspects, and they are found to be topological insulators [16][17][18][19] .
Despite all these studies, it is still unknown whether all pnictogens can have honeycomb lattice structure, especially nitrogen. Secondly, the orbitals involved in the electronic structures near Fermi level are unclear, thus no tight-binding models applicable to these systems have been proposed. The spin-orbit couplings (SOC) are believed to play important role in the electronic structure, which have not been studied seriously. Furthermore, the properties of nanoribbons based on these materials are not known well, such as their magnetism. These questions are highly nontrivial and worth systematic study.
In this paper, we use first-principle calculation to systematically explore the unknown phase of VA elements in the periodic table. Both the structures of bulk and 2D monolayers are optimized and found to be stable. The electronic structures of monolayers are intensively studied, and the evolution of band structures is clearly shown. A tight-binding Hamiltonian is constructed, and the hopping parameters are obtained by fitting to band structures from first-principle calculations. We find large SOC for all elements except nitrogen, in which the on-site contribution is the leading term. Unlike other types of SOC, on-site SOC is solely determined by atomic orbitals of lattice sites; and on-site SOC does not break time-reversal symmetry, therefore these materials can be good candidate for topological superconductors. Electronic structures of both armchair and zigzag nanoribbons of these materials are also studied. The conducting edge states appear in nitrogen and bismuth zigzag nanoribbon. More interestingly, ferromagnetism shows up in the zigzag nanoribbons of P, As and Sb, and in the armchair nanoribbons of As and Sb. The dangling bonds on the edges are mainly responsible for the edge states and the magnetism. The ferromagnetic edge states may contribute to useful applications, such as spin-polarized angular-resolved photoemission spectroscopy (ARPES), spin-polarized scanning tunneling microscopy (SP-STM), and spin-polarized field emission. On the other hand, the appearance of conducting edge states on nitrogen zigzag nanoribbon is a novel and counter-intuitive phenomenon that worth deep investigation.

Results
Layered Structure and Monolayers. The three-dimensional (3D) AB stacking structures and corresponding 2D structures of VA elements are optimized by first-principle calculations. We calculated vibration frequencies for monolayers of all VA elements, and no imaginary frequency was found, suggesting the stability of all these structures. Pnictogens have five outer electrons, and they need three additional electrons to achieve octet stability. Their hybridized sp 3 orbitals account for their similar buckled structures, as shown in Fig. 1, and the optimized structure parameters of nitrogen (N), phosphorus (P), arsenic (As), antimony (Sb) and bismuth (Bi) with and without van de Waals (vdW) correction, are listed in Table 1. We use l bond to denote bond length, and l′ to denote nearest distance between atoms in adjacent layers.d and Δ z stand for inter-layer distance and buckling. a is the lattice constant. As shown in Fig. 1, each layer in these 3D layered structures, has silicene-like structure. The calculated a and l bond for P are 3.32 Å/3.28 Å (with/without vdW) and 2.27 Å for bulk, respectively, and 3.28 Å and 2.27 Å for monolayer, respectively, in agreement with previous study by Z. Zhu et. al. 12 , although their inter-layer distance is larger than our value.
As shown in Table 1, the l bond 's are much smaller than the l′ 's. The structure parameters of monolayer differ only slightly from that of the bulk. It is clear that the structure parameters for nitrogen are much smaller than those of other VA elements. The inclusion of the vdW correction mainly affects d while having only minor effect on the a's, Δ z's and l bond 's. When vdW is included, all the l′ 's are about twice the vdW radii. In other words, the inclusion of vdW interaction does not dramatically change the structure parameters within a layer, even though there are noticeable changes in the d's, such as for P and As. All these phenomena suggest that layers couple with each other by weak vdW interaction and thus can easily exfoliated to form monolayer 20 . This theoretical prediction may provide guidance to experimental synthesis and applications.
Even though there is no monotonous increase of d's, the monotonous increases of a, Δ z's, l bond 's and l′ 's with respect to increasing atomic number from N to Bi are clear. This is in agreement with the periodic law, in that the atoms become larger with increasing atomic number. Interestingly, if we compare the bond length with twice the covalent radius, the bond length is a little bit longer.
It is necessary to point out that, during the optimization process of three-dimensional (3D) structure of N, we found the tiny jittering of total energy as the inter-layer distance changes. This may suggest that the interaction between N layers are extraordinarily weak (if they have any), and the calculations find it difficult to optimize the structure to minimize the total energy. This means that these layers can detach from each other very easily, therefore the 3D structure of N may not exists, while the 2D structure of N can exist. This is supported by the calculation of phonon dispersion for the nitrogen monolayer, which shows non-negative frequencies for all vibration modes. Moreover, due to the interactions between the sheet and the substrate, the presence of substrate can enhance the stability of the monolayer. Experimentally, the 2D nitrogen atomic sheet has been grown epitaxially on GaAs(001) surface 21 .
Electronic Structure of Monolayer. The band structures obtained from our first-principle calculations shown in Fig. 2 (without SOC) and Fig. 3 (with SOC), and the band gaps for all the monolayers are listed in Table 2, with and without SOC. The projected density of states (PDOS) are presented in Fig. 4.
From the analysis of PDOS, it is clear that the electronic structures of all the monolayers are mainly made up of s and p orbitals. N has no 2d orbital; P has 3d orbitals, but they are vacant. Therefore, in the case of N and P monolayer, the contribution of the d orbitals to the total density of states (DOS) is exactly zero. For As, Sb and Bi, SOC in d orbitals has non-zero but only minor contribution to the total DOS. Based on this analysis, we argue that the band structures of all these monolayers can be sufficiently well-depicted by s and p orbitals, and d orbitals are not necessary. Meanwhile, the gap between the lower branch and the middle branch enlarges as atomic number increases due to the relativistic effect.
From the band structures shown in Fig. 2 and Fig. 3, three branches of bands are evident, and they are represented by three different colors in the figures. For all the five monolayers, the lower branch mainly consists of s orbitals while the other two branches mainly consist of p orbitals. For N monolayers, the band structure resembles that of silicene 22 ; PDOS shows that the lower branch and the middle branch are entangled. As atomic number increases, the lower branch gradually disentangles from the middle one, becoming lower in energy, and formed a graphene-like band structure. From analysis of PDOS presented in Fig. 4, as atomic number increases, the s component of the lower branch increases with decreasing p character, but the middle and upper branches show an opposite trend, i.e., the p character   increases as s character decreases. It can be seen from both the band structures and PDOS's, the ranges of all these three bands shrink, indicating the weakening of interactions between orbitals as the covalent bonds become longer. As shown in Table 2, the band gap decreases from 3.70 eV of nitrogen to just 0.55 eV/0.42 eV(without/ with SOC) of Bi as atomic number increases, representing a transition from insulator to semiconductor. This is consistent with the periodic law: the metal characteristic increases as atomic number increases in VA elements. We obtained a slightly smaller band gap (1.97 eV) than previous study of blue phosphorus monolayer (~2 eV) 12 . In the absence of SOC, the Bi monolayer has direct band gap while all others have indirect band gaps. The valence band maximum (VBM) of As, Sb and Bi, and the conducting band minimum (CBM) of Bi are located in Γ point while others are not. When SOC is included, all band structures have indirect band gaps, and the VBM of Bi shifts away from Γ point. We also apply tensile strain to the bismuth monolayer, and calculate the band gap with various strains in the presence of SOC. We find that the band gap decreases with increasing tensile strain, then it closes at approximately 5% tensile strain, and reopens if tensile strain keeps increasing. This is consistent with previous studies that focus on the topological aspect of bismuth monolayer 17 ; they claimed that bismuth monolayer is non-trivial 2D topological insulator.
By comparing the band structure with and without SOC, as shown in Fig. 2 and Fig. 3, respectively, SOC breaks energy degeneracy of certain points in the band structures. As mentioned before, the lower branch is mainly consisted of s orbitals, while the other two branches are mainly from p orbitals, thus SOC has nonzero contribution to p orbitals, but vanishes for s orbitals because the angular moment of s orbital is zero. Therefore, SOC influences the middle and upper branches of bands while having no effect on the lower branch.
Tight-Binding Approach to Monolayer. Based on the analysis of PDOS, it is clear that a tight-binding (TB) model with just s and p orbitals is sufficient to describe the system. Zolymoi et al. have constructed a TB model for silicane and germanane with only s and p orbitals 23 . Similarly, we construct a TB model with SOC.
Here, a † (R i ) and a(R i ) are the creation and annihilation operators of s electron in lattice site R i , and i are the creation and annihilation operators of p electrons in lattice site R i . ε s and ε p in  0 are the on-site energies of the s and p orbitals.  1 contains all the nearest neighbor hoppings and  2 contains all the next nearest neighbor hoppings. Both  1 and  2 can be divided into three parts: hoppings between two s orbitals (h s and ′ h s ), hoppings between two p orbitals (h p and ′ h p ), and hoppings between an s orbital and a p orbital (h sp and ′ h sp ). The corresponding hopping integrals are denoted by γ ss , γ ′ ss , γ ppσ and so on. The first term in  SO represents intrinsic SOC and the second term is Rashba SOC 24 . α and β in  SO are spin indices. (d i × d j ) z corresponds to the z component of d i × d j . μ ij = ± 1 for A(B) sites. d i and d j are the two nearest bonds connecting the next nearest neighbors. −  On Site represents the on-site SOC that comes from the relativistic effect, and it serves as a correction term for the non-relativistic approximation of Dirac equation in central field. It can be proven that S·L does not commute with , therefore the contribution of the SOC to the system is nontrivial. The non-zero matrix elements of the on-site SOC between the orbitals and spin components are listed in Table 3.
All together there are 13 parameters in this TB Hamiltonian, including three parameters of SOC. We have obtained these parameters by fitting TB band structure to the result from first-principle calculations, as shown in Table 4. The comparison between first-principle band structures and TB band structures are presented in Fig. 5 (without SOC) and Fig. 6 (with SOC). It can be seen from Table 4 that these parameters obey the periodic law, with some abnormality. γ sp and γ ppσ are dominant over all other parameters because of their better overlap between orbitals. The hoppings between nearest neighbors are greater than those between next nearest neighbors. Since the SOC in N monolayer is negligibly small,  Zigzag and Armchair Nanoribbons. To be more comprehensive, we investigate nano-ribbons and their properties. Nano-ribbons demonstrate novel physical phenomenon different from the 2D counterpart, and the properties of nano-ribbons of VA elements worth studying. Using the structure parameters obtained from the studies on the monolayers, we have calculated zigzag nanoribbons (ZNR) and armchair nanoribbons (ANR) with various widths. Figure 7(a) shows a particular case of ZNR and ANR of 20 atoms wide. The outermost three atoms on each side are relaxed before the calculation of electronic structures. As shown in Fig. 7(b-f), the edges of nitrogen and bismuth ZNR, and the edges of nitrogen, phosphorus and bismuth ANR are found to be rearranged, while for other cases only minute atom displacements are observed. The nitrogen atoms on the edge of nitrogen ZNR and ANR moves towards the central plane ( Fig. 7(b,d)). The edges of bismuth ZNR are obviously bent (Fig. 7(c)).The outermost bond of phosphorus and bismuth ANR are rotated, making the buckling on the edges larger than that in the center of ANR (Fig. 7(e,f)). There is even a possible reconstruction of bismuth ANR on the edges in such a way that each atom is bonded to three neighboring ones. Compared with the band structure of monolayer, the correspondence of the branches of bands are obvious. A particular example of phosphorus is shown in Fig. 8. The lower branch of ZNR and ANR mimics the counterpart of graphene. Any bands of ZNR and ANR with this correspondence can be obtained by this projection. And the 2D band gap is 2.05 eV for phosphorus ZNR, slightly larger than that of the monolayer.
It is interesting to notice that there are two bands of phosphorus monolayer that do not correspond to any of the bands in the monolayer. The two bands in ZNR are doubly degenerate in energy and spin, but their spins are opposite to each other. There is a gap between these two bands, which is about 0.28 eV, and the Fermi level lies inside the gap. In other words, this system is ferromagnetic, and demonstrates a magnetic moment about 2.02μ B per supercell. It is reasonable since the band below Fermi level is doubly degenerate, and is filled with two electrons from each unit cell. This is different from the antiferromagnetic edge state of graphene ZNR 26,27 . On the other hand, the two bands in ANR are four-fold degenerate, but not spin-polarized. In other words, the edges of ANR of phosphorus are anti-ferromagnetic, possessing a zero magnetic moment. We also calculate the band structures of phosphorus ZNR and ANR that are passivated by hydrogen, and the edge states disappear. Therefore it is possible that dangling bonds are responsible for edge states. One possible explanation for the magnetism is that, the electrons in the dangling bonds are slightly mixed with other states inside the ribbon, and there are slightly delocalized and thus can overlap with the wave functions of the dangling bonds on the other edge.
Similar band structures can be obtained for other cases, as shown in Fig. 9. Surprisingly, the edge states of N and Bi ZNR, and the edge states of N and P ANR are not spin-polarized. And there seems to be no edge states for Bi ANR. Interestingly, the edge states of N ZNR are conducting, which is in contrast   Analysis of charge density reveals the existence of edge states. We found that the charges of this band are mainly distributed on both edges of ZNR. Figure 10 shows the charge distribution of edge states of phosphorus ZNR with 20 atoms wide. The charge distributions for other elements with various widths are very similar, and are not presented here. Interestingly, for bismuth ANR, it shows no edge states because of possible reconstruction in which no dangling bonds exists that contribute to edge states.
The dangling bonds on the edges are mainly responsible for the magnetism. The dangling bonds exist on the outermost atoms on the edges, and they are chemically active in the presence of air. However they can be stable under ideal conditions, such as high vacuum and low temperature. Recently, single dangling bond on Si surface has been observed 28 . The conclusion about the edge states and magnetism are obtained with the most ideal condition. We are not sure if the edge magnetism is robust or not in the presence of air or thermal excitation. However, they are at least stable under the ideal condition, without the presence of disturbance.
The spin-polarized edge states of ZNR are very useful in applications. It can be used in spin-polarized angular-resolved photoemission spectroscopy (ARPES), a necessary technique in the study of topological insulators. It can also be used in spin-polarized scanning tunneling microscopy (SP-STM) as a tool to detect the spin texture of material surfaces. Likewise, spin-polarized field emission can take advantage of this polarized spin as well. The graphene ZNR may not be useful in all these cases, since the edge state of graphene ZNR is not spin-polarized, as mentioned above.

Discussion
Using first-principle calculations, we demonstrate the existence of two-dimensional honeycomb monolayer of VA elements, and determine the corresponding structure parameter. From the three-dimensional structure, it is clear that the systems have AB stacking layered structure (with some uncertainty about N), and each layer is similar to silicene. The geometry of monolayer is very similar to that of the bulk, indicating that layers are bound together by weak van de Waals interaction. Bond length, lattice constant, and buckling increase monotonously with atomic number from nitrogen to bismuth, obeying the periodic law.  We find that the electronic structure obeys the periodic law, as atomic number goes from nitrogen to bismuth. Three branches of bands are obvious in the band structure, and they are similar for all VA elements with minor variations. Spin-orbit coupling (SOC) plays important role in the band structures, and its effect becomes stronger as atomic number increases, however its influence is still limited. It is shown that SOC influences p orbitals, which can be seen from the changes of bands with major characters of p orbitals. There is no noticeable effect of SOC on s orbitals. The band gap decreases monotonously, from 3.7 eV for nitrogen to only 0.55 eV (without SOC)/0.43 eV(with SOC).
From projected density of states, it is clear that the s and p characters of bands changes monotonously. When the atomic number is small, s and p orbitals are mixed, such as for cases of nitrogen and phosphorus. As the atomic number increases, the s character shifts downwards in energy, and becomes the dominant component of lower branch of bands; on the other hand, the p character shifts upwards in energy, becoming the dominant component of the middle branch and upper branches of bands. In all these cases, the effect of d orbital to band structure is negligible. Based on this analysis, we have constructed a tight-binding (TB) Hamiltonian using only one s orbital and three p orbitals from each atom in the unit cell, and the hopping parameters are obtained by fitting the band structures. Large on-site SOCs are found for the monolayers of all these elements except nitrogen.
Lastly, we study the atomic and electronic structure of zigzag nanoribbons (ZNR) and armchair nanoribbons (ANR). Rearrangements in nitrogen and bismuth ZNR and nitrogen, phosphorus and bismuth ANR are found, and there is possible reconstruction in bismuth ANR. In the study of electronic structure, we found an obvious correspondence of band structure to that of monolayer. Band structures of ZNR and ANR show spin-polarized ferromagnetic edge states for phosphorus, arsenic and antimony ZNR and arsenic and antimony ANR. The conducting edge states of nitrogen ZNR is also unusual since it contradicts our convention that materials made from nitrogen are usually insulators.
All these findings contribute to the study of future use of these materials. Firstly, these materials can be used in electronic devices. Secondly, because of their large on-site SOC, they are candidates for topological superconductors. Thirdly, the spin-polarized edge states may be useful in spin-polarized angular-resolved photoemission spectroscopy (ARPES), Spin-polarized scanning tunneling microscopy (SP-STM) and spin-polarized field emission.

Methods
In this paper, we use first-principle calculations to explore the unknown phase of VA elements in the periodic table. Our calculations are based on Plane Augmented Wave (PAW) with Perdew-Burke-Ernzerh (PBE) of exchange-correlation 29 as implemented in the Vienna Ab initio Simulation Package (VASP) code 30 . To ensure reliable result, we use large cut-off energy (400eV) for bulk and monolayers, however we use default cut-off energy for nanoribbons. The systems are restricted to periodic boundary conditions. A vacuum at least 15 Å thick is inserted to eliminate the interaction between different units in the study of system with low dimensions (monolayer and nanoribbon). For optimization, ions are relaxed using conjugate-gradient algorithm until the total force is less than 0.01 eV/ Å. Grimme's scheme for van de Waals (vdW) correction 31 is used to optimize bulk structure, and results are compared with those obtained without vdW correction. Unfortunately, vdW correction is not included in the case of bismuth because of lack of parameters in Grimme's scheme. We believe that vdW has no substantial effect on the layer structure of bismuth monolayer. Some vdW scheme with higher precision 32 could have been used, but getting accurate vdW energy is not the main interest of this paper, and a simple time-saving scheme of Grimme is sufficient to prove our statements. Spin-orbit coupling (SOC) is included only in the calculation of band structure of monolayers.