The topology and robustness of two Dirac cones in S-graphene: A tight binding approach

Present work reports an elegant method to address the emergence of two Dirac cones in a non-hexagonal graphene allotrope S-graphene (SG). We have availed nearest neighbour tight binding (NNTB) model to validate the existence of two Dirac cones reported from density functional theory (DFT) computations. Besides, the real space renormalization group (RSRG) scheme clearly reveals the key reason behind the emergence of two Dirac cones associated with the given topology. Furthermore, the robustness of these Dirac cones has been explored in terms of hopping parameters. As an important note, the Fermi velocity of the SG system (vF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\simeq $$\end{document}≃ c/80) is almost 3.75 times that of the graphene. It has been observed that the Dirac cones can be easily shifted along the symmetry lines without breaking the degeneracy. We have attained two different conditions based on the sole relations of hopping parameters and on-site energies to break the degeneracy. Further, in order to perceive the topological aspect of the system we have obtained the phase diagram and Chern number of Haldane model. This exact analytical method along with the supported DFT computation will be very effective in studying the intrinsic behaviour of the Dirac materials other than graphene.

Chern insulators. It is to note that the Chern insulators exhibit non-zero Hall conductance in the absence of any global magnetic flux 21,22 . The topological phase transition connecting two regimes with a distinct Chern numbers (0 to ±1). The Chern number is a topological index that reveal the information about the phase evolution of wavefunctions around the Brillouin zone (BZ) torus. Moreover, remarkable experimental advances in the field of ultracold atoms [23][24][25][26] and photonic systems [27][28][29][30] have successfully evinced the realization of non-trivial topological phases. One of the key advantages of Haldane model is certainly the experimental accessibility of its phase diagram in case of the non-interacting ultracold systems 25 . Besides, the lattice dynamics can be efficiently attainable via manipulating any single site 31 . Furthermore, an extended version of Haldane model 32 exhibits that the long range interaction generates additional satellite Dirac points. These Dirac points play a crucial role in explaining the sudden quenching dynamics of the model. The nonequilibrium dynamics of the edge current of the semiopen Haldane model has also been reported 33 . Moreover, these symmetry protected edge states are extremely important in the study of the superconductors with Majorana Fermions 34 , spin Hall insulators [35][36][37] , three dimensional topological insulators 38 etc. The paper is organized as follows. In the next section (section-II) we have discussed the exact RSRG scheme on S-Graphene like tetragonal lattice. In section-III we have described an analytical description on the emergence of two Dirac cones in the IBZ of the S-graphene. In the subsequent section-IV the robustness of the Dirac cones has been well explored. In addition, we have proposed a method in section-V to obtain the Haldane lattice of S-graphene for describing its topological phase with non-zero Chern number. Finally, the conclusion has been drawn at the end of this paper.

the RSRG scheme
Our primary aim is to address the emergence of two Dirac cones and their robustness associated with the lattice shown in Fig. 1 keeping the topology intact. This type of lattice has been included in the carbon family by Xu et al. 16 as a stable allotrope of graphene called S-graphene. The tetragonal unit cell of S-graphene is defined by the lattice vectors → =â ai and → =b bj . In this paper we have extensively used the real space renormalization group (RSRG) scheme [39][40][41][42] to compute the dispersion (E-k) relation associated with the system through tight binding Hamiltonian. In particular, we have decimated out the preferred subset of atomic sites of the original lattice to achieve a scaled version of it. The scaled lattice, however, carries all the information of the original lattice. We can perform the process several times without losing any generality. The RSRG process has been initiated by evaluating the set of difference equations for the original lattice using Eq. 1 given below.
Here, ε, ψ i and t ij denote the uniform on-site potential, amplitude of the wave function at i-th site and the hopping matrix element between the i-th and the j-th sites. The electron's hopping is restricted between nearest neighbour (NN) only. It is clear from Fig. 2(a) that we can chose five different hopping parameters of the original lattice keeping the topology intact. Therefore, the difference equations for the section of the original lattice depicted in Fig. 2(a) are given below.  Now, we want to decimate the red shaded atoms shown in Fig. 2(a). For that purpose we have adopted two elegant and legitimate approximations described as follows. First we argue that, since we are solely interested in the formation of Dirac cones near the Fermi energy (E F = ε) we can substitute ε ≈ E . Thereafter, we have chosen a hopping relation = t t t 2 2 1 4 that decouples the whole structure into two identical and non-interacting systems. The above relation is also valid for the uniform hopping parameter of the complete system. With these set of approximations we can substitute the following values to perform the decimation process. 1 . The effect of the first RSRG step on the original lattice has been schematically shown Fig. 3 for better understanding. It is clear from the Fig. 3 that the new scaled lattice consists of only 6 atoms in the unit cell despite 8 atoms in the original system. Hence, the RSRG process reduces the degree of difficulty in obtaining the band dispersion near E F . In a similar vein, we have repeated the RSRG scheme for the output lattice of the first RSRG process in a more rigorous manner. In accord with the first RSRG scheme, we have chosen another section as shown in Fig. 4(a) to initiate the RSRG scheme once again. The red shaded atomic sites have further been decimated using the following difference equations As an example, we have decimated the atomic sites b and e using the following relations We can use similar relations for c and f site also for the corresponding decimation and the renormalized lattice is depicted in Fig. 4(b). Moreover, the decimation process essentially turns the system into a structure topologically equivalent to two interpenetrating identical and non-interacting graphene like hexagonal sheets as depicted in

emergence of two Dirac points
It is evident that the unit cell of the renormalized lattice consists of two atoms and there are only two types of hopping parameters t 3 and t″ present in the system. Such a two-band structure can generically be treated as a two-level system that corresponds to a two-dimensional Hilbert space (   H k 2 ) at each point → q q q ( , ) x y of the Brillouin torus 2 . The components of the momentum vector → q can be defined as  www.nature.com/scientificreports www.nature.com/scientificreports/ two-level system is a 2 × 2 Hermitian matrix. Therefore, the Hamiltonian can be rewritten on the basis of Pauli matrices as follows Here, σ 0 is simply the 2 × 2 identity matrix. The energy levels of the two-level system will thus be described as follow 0 where h k ( ) has the following expression It is worth mentioning that the only role of the term h 0 is to equally shift the energy values of both the bands. Therefore, it is insignificant in determining the band dispersion relation. As a consequence we can straight away substitute h 0 = 0 in above Eq. 7.
In case of the completely renormalized two-level lattice given in Fig. 5 we have obtained the following values of the h parameter , 0 x x y y y z 3 3 Therefore, we can write  www.nature.com/scientificreports www.nature.com/scientificreports/ We now start with the most trivial case, i.e. all the hopping parameters of the original lattice are uniform say, τ. As a consequence, two hopping parameters of the renormalized lattice (t 3 = τ and t″ = t 2 /t 1 = τ 2 /τ = τ) have the same value similar to the original lattice. Therefore, the dispersion relation has the form Moreover, we have plotted the dispersion relation to realize the E-k spectra over the entire tetragonal BZ. As already anticipated we have observed the occurrence of two distinct Dirac cones ' A' and 'B' at specific k-points lying between Γ (0.0, 0.0) → X (0.5, 0.0) and M (0.5, 0.5) → Y (0.0, 0.5) respectively. Therefore, it is obvious that near the first Dirac point ' A' , k y is invariably zero. Whereas, k y has a constant value 0.5 near the second Dirac point 'B' . The dispersion relations between different symmetry points has been described in Table 1.
It is worthy to note that the dispersion relations are symmetric about E F . This is a signature of particle hole symmetry in the system. Since, we have only considered the NN hoppings in a two-level system, this symmetry is well expected 17 . A pictorial representation of all the dispersion relations has been provided in Fig. 6. The positions of these Dirac points have been determined from the relation h(k) = 0. It turns out that the first Dirac point is at A (1/3, 0) and the second Dirac point is at B (1/6, 1/2). We have then calculated the Fermi velocity (v F ) by expanding the dispersion equations in Taylor series near Dirac points. Neglecting the quadratic and higher order terms of the series we can compare the result with In the above expression q x = 2πk x /a where, = . a 5 30 is the lattice constant in angstrom along x direction. We have considered all the bond lengths hence, the hopping parameters are equivalent to that of graphene 2 , τ = .
2 7 eV. Therefore, the carriers near both the Dirac points possess same Fermi velocity v F = ± τ a 3 / = 3.7 × 10 6 m/s  c/80. It is to note that the Fermi velocity of graphene 2,43 has previously been reported to be ≈1 × 10 6  c/300. The results suggest that the Fermi velocity of the S-graphene network is significantly (≈3.75 times) large compared to that of graphene.
In order to support the band structure obtained analytically, we have performed the density functional theory (DFT) computation with the help of SIESTA package 44 . During the computation, the Perdew-Burke-Ernzerhof (PBE) parameters are employed to treat the exchange-correlation part of density functional. Besides, the generalized gradient approximation (GGA) with double ζ plus polarized basis sets has been chosen with 32 × 32 × 1 Monkhorst-Pack. The band structure of the self consistent ground state of S-graphene computed from DFT also exhibits two Dirac cones in the band structure as described in Fig. 6(b). It is to note that, while achieving the self-consistent ground state of S-graphene we have not imposed any symmetry restriction. As a result, the bond lengths and the corresponding hopping parameters of the system become non-uniform. However, the Figure 5. The effect of the second RSRG process on the lattice obtained after first RSRG process. The decimated atoms of the initial lattice are marked with shaded red. The final lattice is basically two interpenetrating graphene like hexagonal systems. New allowed hoppings are marked by different colours (red and green) for two distinct hexagonal rings. The unit cells have been marked with a box.

BZ path constraint
Dirac cone E-k relation

Robustness of the Dirac cones
Now, we want to deal with the more generic case i.e. the hopping parameters of the renormalized lattice t 3 and t″ are distinct. In this case, the dispersion relation is given by Eq. 10. We aim to attain the relation between t 3 and t″ that allows the band gap opening in the S-graphene system. In this context, we would like to define a dimensionless parameter r that satisfies the relation = ″ = r t t t t t / / 3 2 1 3 . It is evident from Eq. 10 that the robustness of the Dirac point ' A' is protected by the relation π . This essentially gives us the condition π = − r sec k (2 )/2 x . On the other hand, the robustness of the Dirac point 'B' is protected by the condition π = r sec k (2 )/2 x . There is an inherent constraint > r 0 applied in the above two conditions. Therefore, = r 1/2 is the lower limit for the survival of the Dirac cones. In other words, > t t t 2 1 3 2 essentially removes the degeneracy and opens a gap in the system. It is evident that the relation = t t t 2 2 1 4 that decouples the system into two non-interacting distinct parts, forbids us to tune the Dirac cone individually. In brief, the simultaneous implementation of the conditions = t t t 2 2 1 4 and > t t t 2 1 3 2 induces band gap even in the presence of inversion symmetry.
Next, we are interested to investigate the robustness of the individual Dirac cone. The underlying objective is obviously to discard the relation = t t t 2 2 1 4 . We have availed the NN tight-binding (NNTB) model for the p z orbital owing to the fact that these orbitals solely contribute in the formation of Dirac cones. The NNTB Hamiltonian in the Wannier basis can thus be expressed as follows Here, ε i , † c i , c i and + t i i , 1 represent the on-site potential, creation, annihilation operator and NN hopping parameter between i-th to (i + 1)-th site. For the consistency check, we have obtained the band spectra corresponding to the Hamiltonian given in Eq. 12 for uniform onsite energy and hopping parameter i.e. ε ε = = 0 i and τ = 1 eV. The system exhibits the existence of two Dirac points under above conditions as depicted in Fig. 7(a). Furthermore, we have obtained the band structure for different set of hopping parameters without any constraint. Initially, we want to inspect the role of individual variation of hopping parameters in determining the fate of Dirac cones. We start with the sole variation t 2 . Let the other hopping parameters are uniform say τ (>t 2 ). It has been observed that the position of the Dirac points move towards the higher symmetry points with decreasing t 2 . There is a critical value t 2 = τ/ 2 below which the degeneracies break and the band spectra exhibit a band gap in the system. This is well expected as this condition satisfies the hopping relation = t t t 2 2 1 4 . On the other hand, τ < t 2 can not break the degeneracies. In a similar way, the sole variation of t 3 , however, can not split the Dirac cones for t 3 < τ. This observation is also valid for t 3 = 0. It can be argued that both the Dirac points are situated at the BZ path parallel to k x axis. Hence, breaking of the periodicity along the y-axis (i.e. t 3 = 0) transforms the sheet into 1D nanoribbon of particular width. This ribbon will also exhibit Dirac cones. In contrary, the condition t 3 > τ, induces a band gap in the system above t 3 = 2τ. The above criteria is also valid for t 4 . In addition, we have considered the asymmetrical hopping parameters and obtained that the Dirac point A is robust to the parameter t 4 upto the relation t 4 = 2t′. Here, we have considered, t 1 = t 3 = t′ and t, t 2 ≥ t′.
Additionally, we have explored two distinct conditions that lift the degeneracy of one Dirac point but other one remain unaltered. The analytical calculation reveals that the condition 0 < v ≤ 4 3 protects the Dirac point A ≤ v < 0 protects the Dirac point B. Here, v is the ratio between the non-zero vertical renormalized hopping parameter (between 'a' and 'd' site of Fig. 4) and the rest of the uniform hopping parameter. In particular, t 4 < 2τ/3 is a special criterion that splits the Dirac cone B but not A. Whereas, another relation, t 1 > t 2 (=t 4 ) ≤ t and t 3 ≤ 2t 2 removes the degeneracy of A but not B. These results are listed in Table 2 and depicted in Fig. 7(b,c) for better understanding. It is to note that, the tuning of the Dirac cones alters the slope of the bands near E F giving rise to the significant change in their carrier mobility. In order to realize the feasibility of such hopping relations in band gap tuning we have applied stress in the system along x-axis using DFT. We have observed (not shown here) that 1% strain can split the Dirac cones. The result is consistent as the hopping parameter strongly depends on the separation between the neighbouring atoms.
Subsequently, the on-site symmetry within the unit cell has been broken, while the hopping parameters remain uniform. Therefore, the atom with on-site energy +ε is connected with the three neighbouring atoms with on-site potential −ε and vice versa. We have observed that even a small anisotropy in the on-site potential can induce band opening in the system as depicted in Fig. 7(d). This result is highly expected as in the renormalized lattice this offset in the on-site energy will break the inversion symmetry of the system although, the time-reversal symmetry is present in the system. It is worth mentioning that the simultaneous presence of these two symmetry breaking terms will provide us the rich physics associated with its topological phase.

the topological phase Diagram
The hexagonal system is a non-Bravais lattice. It comprises two distinct sublattices A and B per unit cell as shown in Fig. 8(a). If we consider single tight binding orbital at each lattice site and allow the NN hopping (t 1 ) only, then the system exhibits semimetallic band structure and the particle-hole symmetry is restored. This unreal particle-hole symmetry can be destroyed in presence of next nearest neighbour (NNN) hopping ( = t t /3 2 1 ) but the valence and conduction bands still touch each other at two inequivalent symmetry points K and K′, depicted in Fig. 8(b). Moreover, the degeneracy can be lifted on breaking the inversion and time-reversal symmetry of the system. In the Haldane model 21 , the spatial inversion symmetry has been broken by considering two different Semenoff mass terms M and −M for the two sublattices. On the other hand, the local time-reversal symmetry has been broken by means of staggered magnetic flux ( ( ) φ Φ ). It is to note that the global time reversal symmetry must be conserved in this model. Therefore, the staggered magnetic flux must be so chosen that the total magnetic flux through each hexagonal plaquette remains zero as depicted in Fig. 8(a). As a consequence the NN hopping parameters are real while the NNN hopping parameters are complex. Thus, the appropriate two-level Hamiltonian of the Haldane model can be written as

and NNN vectors
Complex second NN hopping has also been shown. Besides, two distinct sublattices are marked with red and blue colours. (b) The BZ for hexagonal lattice with reciprocal lattice vectors → d 1 and → d 2 . Here, K and K′ are two inequivalent lattice points. (c) Chern phase diagram for the Haldane model of the hexagonal system, plotted as a function of the ratio of the on-site energy M also known as the sublattice symmetry breaking Semenoff mass term to the NNN hopping term against the time-reversal symmetry breaking staggered flux φ. The background white region is topologically trivial phase with Chern number ν = 0 and the coloured region is topologically nontrivial Chern phase with Chern number ν = ±1. is characterized by the non-zero Chern number ν = ±1 21,33 . Under such condition the charge conducting edge modes occur for any finite sized system. The Chern phase diagram of the hexagonal lattice has been described in Fig. 8(c). Henceforth, we have proposed an elegant method to reveal the Chern phase diagram of the non-hexagonal system considered in this work. Therefore, we start with the completely renormalized two-level system with uniform lattice parameters. To achieve the Haldane lattice we have introduced some imaginary sites with zero probability amplitude of the tight binding orbital. The sites are marked black in the Fig. 9(a). In particular, the hopping to and from these imaginary sites are restricted. Such consideration provides the hexagonal plaquette to break the time reversal symmetry locally as depicted in Fig. 9(a). As a consequence, the NN hopping remain unperturbed while, the NNN hoping parameter gains an additional Aharonov-Bohm (AB) phase 21 . Besides, we have considered two distinct onsite potential of the sites in the unit cell to break the sublattice symmetry of the system. These conditions invariably satisfy all the basic requirements of Haldane model. It is evident that the Haldane lattice is also a two-level system. Therefore, we have calculated the Haldane Hamiltonian using Eq. 13 for S-graphene system. Furthermore, the 2 × 2 matrix has been expanded in terms of Pauli matrices as described in Eq. 6. It is to note that near both the Dirac points the relation holds good. However, unlike the previous case h z has a non zero value. As a result, the non-zero energy gaps near Dirac points are given as below It is well known that the closure of either one gap indicates the topological phase transition. Therefore, the relation 2 φ | ± | = defines the boundary between normal topological insulation phase. The Chern number (ν) has the following form www.nature.com/scientificreports www.nature.com/scientificreports/ The region of the topological phase (ν = ±1) in the phase diagram of the Haldane model has been compared between the tetragonal and hexagonal lattice as depicted in Fig. 9(b).

conclusion
In summary, we have introduced an elegant analytical description to address the surprising emergence of two Dirac cones in the irreducible BZ (IBZ) of the tetragonal lattice identical to S-graphene. The exact RSRG scheme efficiently reduces the degree of difficulty of the model from eight atoms per unit cell to the two atoms per unit cell. Moreover, the NNTB Hamiltonian of such two-level tetragonal system essentially exhibits two distinct Dirac points in the (IBZ). As an important note, the Fermi velocity of the SG system (  v c/80 F ) is remarkably larger (3.75 times) than that of graphene. Besides, we have attained a particular relation between the hopping parameters that decouples the system into two identical interpenetrating hexagonal parts. The above condition also restrict us to tune the individual Dirac cones. The position of the Dirac cones can be shifted along the momentum axis for different set of hopping parameters. Besides, the relation between the momentum and hopping parameter that protects the degeneracies has also been explored. Furthermore, we have obtained the conditions for tuning the individual Dirac cones. The non-uniform hopping parameters essentially breaks the symmetry of slope of the linear bands near Dirac points A and B. Therefore, we can control the Fermi velocity by tuning the individual Dirac cones. In addition, we have proposed a scheme to obtain the topological Haldane lattice of non-hexagonal lattice S-graphene. Moreover, the Chern number has also been evaluated in the topologically trivial and non-trivial states to achieve the complete description of topological phase transitions in S-graphene.