Multiple Unpinned Dirac Points in Group-Va Single-layers with Phosphorene Structure

Emergent Dirac fermion states underlie many intriguing properties of graphene, and the search for them constitute one strong motivation to explore two-dimensional (2D) allotropes of other elements. Phosphorene, the ultrathin layers of black phosphorous, has been a subject of intense investigations recently, and it was found that other group-Va elements could also form 2D layers with similar puckered lattice structure. Here, by a close examination of their electronic band structure evolution, we discover two types of Dirac fermion states emerging in the low-energy spectrum. One pair of (type-I) Dirac points is sitting on high-symmetry lines, while two pairs of (type-II) Dirac points are located at generic $k$-points, with different anisotropic dispersions determined by the reduced symmetries at their locations. Such fully-unpinned (type-II) 2D Dirac points are discovered for the first time. In the absence of spin-orbit coupling, we find that each Dirac node is protected by the sublattice symmetry from gap opening, which is in turn ensured by any one of three point group symmetries. The spin-orbit coupling generally gaps the Dirac nodes, and for the type-I case, this drives the system into a quantum spin Hall insulator phase. We suggest possible ways to realize the unpinned Dirac points in strained phosphorene.


Introduction
Recent years have witnessed a surge of research interest in the study of Dirac fermions in condensed matter systems, ranging from graphene and topological insulator surfaces in two-dimensions (2D) to Dirac and Weyl semimetals in 3D, [1][2][3][4] which possess many intriguing physical properties owing to their relativistic dispersion and chiral nature.Especially, 2D Dirac fermion states have been extensively discussed in honeycomb lattices, commonly shared by group-IVa elements with graphene as the most prominent example, [5][6][7][8][9] for which Dirac points are pinned at the two inequivalent high-symmetry points K and K of the hexagonal Brillouin zone (BZ), around which the dispersion is linear and isotropic.Later on, 2D Dirac points on high-symmetry lines were also predicted in some nanostructured materials, 10 including graphynes 11 and rectangular carbon and boron allotropes. 12,13 owever, the possibility of 2D Dirac points at generic k-points has not been addressed, and such Dirac point has not been found so far.
Meanwhile, the exploration of 2D materials built of group-Va elements (P, As, Sb, and Bi) has just started.5][16][17][18][19][20][21] While 2D allotropes with different lattice structures have been predicted and analyzed for the other group-Va elements, [22][23][24][25] we note that the puckered lattice structure similar to phosphorene has been demonstrated experimentally for Sb [26][27][28][29] and Bi [30][31][32][33] (down to single-layer) grown on suitable substrates, and been predicted for As as well. 22Motivated by these previous experimental and theoretical works, and in view of the ubiquitous presence of the Dirac fermions and the associated interesting physics, one may wonder: Is it possible to have Dirac fermion states hosted in such 2D puckered lattices?
A simple consideration shows that here any possible Dirac point cannot occur at high-symmetry points.The reason is that each Dirac point at k must have a time reversal (TR) partner at −k with opposite chirality, whereas the BZ of the puckered lattice has a rectangular shape, of which all the high-symmetry points are invariant under TR.Therefore, if Dirac states indeed exist in such systems, they must be of a type distinct from those in graphene.
In this work, we address the above question by investigating the electronic structures of group-Va 2D puckered lattices.We find that Dirac fermion states not only exist, but in fact occur with two different types: one type (referred to as type-I) of (two) Dirac points are located on high-symmetry lines; while the other type (referred to as type-II) of (four) Dirac points are located at generic k-points.Depending on their reduced symmetries, dispersions around these points exhibit different anisotropic behaviors.Points of each type can generate or annihilate in pairs of opposite chiralities, accompanying topological phase transitions from a band insulator to a 2D Dirac semimetal phase, and since they are not fixed at high symmetry points, their locations can be moved around in the BZ.Particularly, to our best knowledge, the novel fully-unpinned (type-II) 2D Dirac points are discovered here for the first time.In the absence of spin-orbit coupling (SOC), each Dirac node is protected from gap opening by a sublattice (chiral) symmetry, which can in turn be ensured by any one of three point group symmetries.The inclusion of SOC could gap the Dirac nodes, and in the case of type-I nodes, it transforms the system into a quantum spin Hall (QSH) insulator phase.All these properties make the system distinct from graphene and other 2D materials.We further suggest that the novel unpinned Dirac points can be experimentally realized by the strain engineering of phosphorene.Our discovery therefore greatly advances our fundamental understanding of 2D Dirac points, and it also suggests a promising platform for exploring interesting effects with novel types of Dirac fermions.

Results
A group-Va pnictogen atom typically forms three covalent bonds with its neighbors.As shown in Fig. 1 for a single-layer phosphorene structure, the P atoms have strong sp 3 -hybridization character hence the three P-P bonds are more close to a tetrahedral configuration.This results in two atomic planes (marked with red and blue colors) having a vertical separation of h comparable to the bond length.In each atomic plane, the bonding between atoms forms zigzag chains along y-direction.
The unit cell has a 4-atom basis, which we label as A U , B U , A L , and B L (see Fig. 1(c)), where U and L refer to the upper-and lower-plane respectively.The structure has a nonsymmorphic D 2h (7)   space group which includes the following elements that will be important in our discussion: an inversion center i, a vertical mirror plane σ v perpendicular to ŷ, and two 2-fold rotational axes c 2y and c 2z .Note that due to the puckering of the layer, the mirror planes perpendicular to x and ẑ are broken.With the same valence electron configuration, As, Sb, and Bi possess allotropes with similar puckered lattice structures.
To study the electronic properties, we performed first-principles calculations based on the density functional theory (DFT).The details are described in the Methods.The calculated geometric parameters of group-Va 2D puckered lattices with D 2h (7) symmetry are summarized in the Supplementary Information.The obtained structures agree with the experiments and other theoretical calculations. 17,22,29,33 Th lattice constants a > b, reflecting that the inter-chain coupling is weaker than the coupling along the zigzag chains.The angle θ 2 increases from ∼ 70 • for P to ∼ 85 • for Bi whereas θ 1 remains ∼ 95 • .The inter-plane separation h as well as the bond lengths R 1 and R 3 increase by almost 1 Å; from P to Bi, while R 2 , the distance between sites of neighboring zig-zag chains, increases only slightly, implying that the inter-chain coupling becomes relatively more important with increasing atomic number.
We first examine their corresponding band structures without SOC, whose effect will be discussed later.The results are shown in Fig. 2. The puckered lattice of P is a semiconductor with a bandgap around Γ-point.The band evolution around Γ-point from Sb to Bi and the appearance of type-I Dirac points in Fig. 2 are reminiscent of a band-inversion process.Indeed, by checking the parity eigenvalues at Γ, one confirms that the band order is reversed for Bi around Γ-point (see Supplementary Informa-tion).For a better understanding, we construct a tight-binding model trying to capture the physics around Γ-point.Since the low-energy bands there are dominated with p z -orbital character (Fig. 2), we take one orbital per site, and include couplings along R 1 and R 2 in the same atomic plane (with amplitudes t 1 and t 2 respectively) as well as nearest-neighbor inter-plane hopping along R 3 (with amplitude t ⊥ ) (see Supplementary Information).Written in the basis of (A U , A L , B U , B L ), the Hamiltonian takes the form: where Q(k) is a 2×2 matrix of the Fourier transformed hopping terms (see Supplementary Information).The Hamiltonian (1) can be diagonalized and possible band crossings can be probed by searching for the zero-energy modes, which exist when the condition the system is a band insulator; when λ < 1, it is a 2D Dirac semimetal.The transition occurs at the critical value λ c = 1 when the conduction and valence bands touch at Γ-point and the band order starts to be inverted.This corresponds to a quantum (and topological) phase transition, 34 during which there is no symmetry change of the system.
Model (1) captures the trend observed in DFT results.The overlap between-p z orbitals is larger along the R 3 bond, hence one expects that t ⊥ > t The emergence of low-energy relativistic chiral modes is the most remarkable property of Dirac points. 34To explicitly demonstrate this, we expand Hamiltonian (1) around each Dirac point, which leads to the low-energy Hamiltonian where q is the wave-vector measured from each Dirac point, τ = ±1 for D and D , σ i 's are Pauli matrices for the sub-space spanned by the two eigenstates at the Dirac point (apart from the Bloch phase factor): and v y = b 4(t 1 + t 2 ) 2 − t 2 ⊥ are the two Fermi velocities.The form of (2) may also be argued solely from symmetry.Compared with graphene, these type-I points are unpinned from the highsymmetry points.They can be shifted along Γ-X 2 (and even pair-annihilated) by varying system parameters such as λ, although they cannot go off the line as constrained by the symmetries.
In addition, different from graphene, 6 the dispersion here is anisotropic, characterized by two different Fermi velocities., and TR that map between the two pairs.Expansion to leading order in each wave-vector component q i gives (see Supplementary Information) where q is measured from F (or F ), w, w , m 0 , and m 1 are expansion coefficients.Two Dirac points appear at (0, ±q 0 ) with q 0 = m 0 /m 1 when sgn(m 0 /m 1 ) = 1, corresponding to a local band inversion around q = 0. Further expansion of the Hamiltonian around the Dirac point (0, νq 0 ) This demonstrates that the two points at ν = ±1 are of opposite chirality, as required by σ v .The dispersion is highly anisotropic (at leading order, characterized by three parameters: w, q 0 , and w ) because the Dirac point is at a generic k-point with less symmetry constraint, as compared with type-I Dirac points.
Unlike in 3D systems, Dirac nodes in 2D have a codimension of two hence are generally not protected from gap-opening. 34In the absence of SOC, however, the Dirac nodes here are stable due to the protection by sublattice (chiral) symmetry between {A i } and {B i } (i = U, L) sites, which allows the definition of a winding number 35,36 (i.e., quantized Berry phase in units of π) along a closed loop encircling each Dirac point: , where A k is the Berry connection of the occupied valence bands.And for a 2D Dirac point, the sign of N (or the ±π Berry phase) is also referred to as the chirality. 6Using DFT results, we numerically calculate the Berry phase for each Dirac point and indeed confirm that they are quantized as ±π.The signs are indicated in Fig. 1(d).
More interestingly, in the puckered lattice with a four-atom basis in a non-coplanar geometry, the sublattice symmetry can be ensured by any one of three independent point group symmetries: i, c 2y , and c 2z .The resulting protection of Dirac nodes can be explicitly demonstrated in low-energy models.For example, consider the type-I points described by Eq.( 2).There the representations of i, c 2y , and c 2z (denoted by P, R y , and R z , respectively) are the same, viz., σ x .Then the symmetry requirement R y H τ (q x , q y )R −1 y = H τ (−q x , q y ) by c 2y directly forbids the presence of a mass term mσ z .Meanwhile, since i and c 2z maps one valley to the other, they protect the Dirac nodes when combined with TR (or σ v if it is unbroken), e.g., considering the combined symmetry of c 2z and TR (with representation T = K the complex conjugation operator): H τ (q), which again forbids a mass generation.The underlying reason i, c 2y , and c 2z each protects the Dirac node is that they each maps between the two sublattices hence ensures the sublattice (chiral) symmetry.In comparison, the mirror plane σ v maps inside each sublattice, hence it alone cannot provide such protection.This reasoning is general and applies to type-II points as well.(In model (3), i, c 2y , and c 2z have representations as σ x by construction, and when combined with TR, again each forbids the generation of a mass term ∼ mσ z .See Supplementary Information.)We stress that the three symmetries i, c 2y , and c 2z each protects the Dirac points independent of the other two.For example, we could disturb the system as in Fig. 4 such that only one of the three symmetries survives.The corresponding DFT results confirm that the Dirac nodes still exist.Thus the crystalline symmetries actually offer multiple protections for the Dirac nodes in the current system.
SOC could break the sublattice symmetry.Hence when SOC is included, the Dirac nodes would generally be gapped. 37For type-I points, treating SOC as a perturbation, its leading-order symmetry-allowed form is H SOC = τ ∆σ z s z , where s z is Pauli matrix for real spin.This is similar to the intrinsic SOC term in graphene, 38 which opens a gap of 2|∆| at the Dirac points.For the type-II points, we obtain H SOC = ηq y σ z s z in model ( 3) hence a gap of 2q 0 |η| is also opened at these Dirac points.Gap opening by SOC is closely related to the QSH insulator phase. 1,2,38 H the band topology can be directly deduced from the parity analysis at the four TR invariant momenta. 39This means that only the band inversion at Γ between the two type-I points contributes to a nontrivial Z 2 invariant; whereas that associated with type-II points does not.It follows that Sb is topologically trivial since it has only type-II Dirac points, while Bi is nontrivial since it has additional type-I points.These results are in agreement with previous studies. 33eaking all three symmetries i, c 2y and c 2z can also generate a trivial gap term mσ z at the Dirac points, which competes with the SOC gap.For example, this happens when each atomicplane forms additional buckling structure. 33Nevertheless, as long as the trivial mass term does not close the SOC-induced gap, by adiabatic continuation the band topology will not change.

Discussion
Due to their different locations and the associated symmetries, the two types of Dirac points here exhibit properties distinct from that of graphene.With preserved sublattice symmetry and in the absence of SOC, the Dirac nodes are topologically stable-they can only disappear by pairannihilation between opposite chiralities.This is unlikely for graphene since the Dirac points there are pinned at the high-symmetry points.In contrast, the two types of Dirac points here are less constrained.Pair-annihilation (pair-generation) indeed occurs during the quantum phase transition as observed from the band evolution.
2][13] Meanwhile the type-II points discovered here are completely new.They are fullyunpinned and has highly anisotropic dispersions.With this discovery, now we can have an almost complete picture: 2D Dirac points can occur at high-symmetry points, along high-symmetry lines, and also at generic k-points.
It is possible to have Dirac points, originally sitting at high-symmetry points, to become unpinned when crystalline symmetry is reduced due to structural distortions.However, we stress that the type-II points here are distinct in that they are realized in a native crystalline structure with relatively high symmetry.Only in such a case, we can have a sharp contrast between generic k-points where the group of wave vectors is trivial and the high-symmetry k-points where the group is non-trivial, and accordingly the type-II point can move around (hence fully-unpinned) without any symmetry-breaking.More importantly, it is just because that type-II points occur in a state with high symmetry that the Dirac nodes can be protected (in the absence of SOC): as we discussed, the various crystalline symmetries ensure the protection of the Dirac nodes from gap opening.
It is remarkable that the two different types of Dirac points can coexist in the same 2D material.We emphasize that it is a result of the lattice structure and the valence character of the pnictogen elements.Our DFT result indeed shows that even starting from the P lattice, the two types of Dirac points can be separately tuned to appear or disappear by lattice deformations.
For example, we find that the type-II Dirac points can be generated in phosphorene by applying uniaxial tensile strains along the y-direction.The DFT result in Fig. 5 and Fig. 6 indeed shows the band inversion on Γ-X 1 and the formation of four type-II Dirac points around a strain of 16%.
Since phosphorene has excellent mechanical properties and a critical strains > 25% has been predicted, 40 it is promising that the novel strain-induced topological phase transitions and the appearance of type-II Dirac points can be directly observed in strained phosphorene.The lattice deformation that produces type-I points is discussed in the Supplementary Information.Similar scenarios occur for other group-Va elements as well.
So far, single-layer As, Sb, and Bi in their free-standing form have not been realized yet.Nevertheless, in view of the rapid progress in experimental techniques, we expect that these materials could be fabricated in the near future.7][28][29][30][31][32][33] Besides the topological properties, the presence of Dirac states is expected to endow these 2D materials with many intriguing properties for applications, such as the very high mobility, the half-quantized quantum Hall effect, 41 the universal optical absorption 42 and etc. Due to the highly anisotropic dispersions of these new Dirac points, the electronic transport properties such as the conductivities would show strong direction dependence.In addition, since there is no symmetry connection between the two types of Dirac points, when they are both present, it is possible to independently shift each type of points relative to the Fermi level, e.g. by strain engineering, leading to self-doping and even the interesting scenario with both electron-like and hole-like Dirac fermions in the same system.With the multiple Dirac points with different chiralities, it is possible to further control the carriers near different Dirac points for valleytronic applications.
In conclusion, based on first-principles calculations of 2D allotropes of group-Va elements with puckered lattice structure, we predict the coexistence of two different types of Dirac points: Dirac points on high symmetry lines and at generic k-points.In particular, the 2D Dirac points at generic k-points are fully-unpinned, have highly anisotropic dispersions, and are discovered here for the first time.Combined with low-energy effective modeling, we unveil the low-energy properties of these Dirac points.We show that their appearance is associated with the band inversion process corresponding to a topological phase transition.The topology/symmetry protection of the Dirac nodes is analyzed in detail.Interestingly, because of the unique lattice structure, there is a triple-protection of the nodes by three independent point group symmetries.This also implies versatile methods to control the locations as well as the dispersions around the Dirac points.When SOC is strong, the Dirac nodes are gapped, and in the case of type-I points (such as for Bi), this drives the system into a QSH insulator phase.We further show that the topological phase transition and the novel unpinned Dirac points can be realized in strained phosphorene.Our work represents a significant conceptual advance in our fundamental understanding of 2D Dirac points.The result also suggests a new platform to explore novel types of 2D Dirac fermions both for their fascinating fundamental properties and for their promising electronic and valleytronic applications.

Methods
First-principles calculations.Our first-principles calculations are based on the density functional theory (DFT) implemented in the Vienna ab-initio simulation package. 43The projector augmented wave pseudopotential method is employed to model ionic potentials. 44Kinetic energy cutoff is set to 400 eV and k-point sampling on the rectangular BZ is with a mesh size 20×20.The minimum vacuum layer thickness is greater than 20 Å; which is large enough to avoid artificial interactions with system images.The structure optimization process is performed including SOC with the local density approximation (LDA) for the exchange-correlation energy 45 and with van der Waals corrections in the Grimme implementation. 46The force convergence criteria is set to be 0.01 eV/ Å.
Hybrid functional (HSE06) 47 is used for the band structure calculations.
Correspondence Correspondence should be addressed to Yuhao Lu, Hsin Lin or Shengyuan A. Yang.
Additional Information Supplementary information is available in the online version of the paper.

Figure 6
Band structure of phosphorene with uniaxial tensile strain applied along ydirection.Here the dispersion is plotted along a path crossing point F (as indicated in Fig. 5) and perpendicular to the Γ-X 1 direction.Two type-II Dirac points near F can be clearly observed at strain value > 16%.
Fig.1(d)) shows that they are indeed Dirac points (see Fig.3(a)).Furthermore, along Γ-X 1 line, with two band touching points at (0, ±k D ) where k D = (2/b) arccos(λ).The direct gap at Γ can be obtained as ∆ = 2[t ⊥ − 2(t 1 + t 2 )].Hence this simple model indeed captures the emergence of two Dirac points D and D , along with a transition as parameter λ varies: when λ > 1, Next, we turn to the fully-unpinned type-II Dirac points.The four type-II Dirac points start to appear for Sb in our DFT result, located close to the Γ-X 1 line.They can be more clearly seen for Bi (see Fig.3(b)).Again the band evolution implies a local band inversion near F and F .Here F and F (on Γ-X 1 ) are the mid-points of the lines connecting each pair of the type-II points.The low-energy bands there are mainly of p x -orbital character.To reproduce the fine features using a tight-binding model would require more hopping terms.Instead, we construct a low-energy effective Hamiltonian around point F (F ) based on symmetry analysis.There the Hamiltonian is constrained by σ v which maps inside each pair (labeled by µ = ±1 for F and F ), and by i, c 2y , c 2z

Figure 1
Figure 1 (a,b) Top-and side-view of the 2D puckered lattice structure.The green

Figure 2
Figure 2 Band structures of group-Va elements with 2D phosphorene lattice structures

Figure 3
Figure 3 Dispersion around (a) type-I and (b) type-II Dirac points for 2D puckered Bi.

Figure 4
Figure 4 Perturbations which leave only one of the three symmetries preserved: (a) c 2z

Figure 5
Figure 5Band structure of phosphorene with uniaxial tensile strain applied along y- Supplementary Information).The result shows that (t ⊥ , t 1 , t 2 ) changes from (2.50, 0.77, 0.33) for Sb to (1.86, 0.63, 0.35) for Bi (units in eV).Hence λ crosses the critical value from Sb to Bi, indicating the band inversion at Γ and the appearance of two Dirac points.
1 > t 2 .By fitting the DFT bands around Γ-point, one finds that from P to Bi, t ⊥ decreases a lot, while t 2 increases and becomes rela-tively more important (