Six-fold symmetry origin of Dirac cone formation in two-dimensional materials

Dirac materials possess many excellent electrical properties, resulting that the search and design of Dirac materials have become a hot research area. Revealing the formation conditions of Dirac cone (DC) can provide theoretical guidance for the search and design of Dirac materials. To obtain the necessary conditions for the formation of DC of two-dimensional (2D) materials with six-fold symmetry (SFS), the DC formation mechanism was analyzed by the ‘divide-and-couple’ approach in the framework of tight-binding theory, confirmed by the subsequent density functional theory calculations. The simple ‘6n + 2’ rule was proposed to determine whether the 2D materials with SFS have DCs, i.e. when the number of atoms in a unit cell is 6n + 2, the systems would possess DCs at the vertex of Brillouin zone for the 2D materials composed of the elements of the IV main group. Moreover, the ‘3n + 1’ rule was derived as the condition for the DC formation in graphene-like silagraphene with SFS and used to design a silagraphene Si6C8 with DCs. Understanding the DC formation mechanism of 2D materials with SFS not only provides theoretical guidance for designing novel Dirac materials but also sheds light on the symmetry origin of the formation mechanism of DC.

However, the 2D materials with DCs are still rare and become the searching targets for potential applications. The mechanism discussions on the DC formation claim the following requirements: (a) specific symmetry; (b) appropriate tight-binding (TB) parameters; (c) there are few bands other than DC at Fermi level [25].
Understanding the formation mechanisms of DC band structures is of great significance to design Dirac materials. Some researchers [9,10,26,27] use the real space renormalization group (RSRG) [28] scheme to explain the formation mechanisms of the DC band structures of 2D materials. In the RSRG scheme, the atomic cluster between two vertex atoms is reduced into an effective hopping term between the two atoms, thus the complex system can be approximated by a simpler system whose origin of the DC is readily explained. For a detailed introduction to the RSRG scheme, see a recent review article [29].
We have previously explored the formation origin of the DC band structures of some 2D systems using 'divide-and-couple' scheme in the framework of TB, in which, according to the characteristics of the system, atomic wave functions are divided into predefined groups, and then the couplings between atomic wave functions are carried out step by step conceptually. The formation origin of DCs can be revealed by the formation process of the band structure [16,17,30,31]. And through the analysis of the process, we can derive the conditions that the TB parameters (hopping parameters and onsite energy parameters) should be satisfied. For example, in the study about t1-SiC, by using 'divide-and-couple' scheme, the conditions that the TB parameters should be satisfied to make the materials have DC are easily derived [16]. The clear physical pictures of the 'divide-and-couple' scheme can provide insights for the design of DC materials.
Because the 2D materials with DCs always possess certain symmetry [25], one can use the 'divide-and-couple' scheme to explain the formation origin of the DCs of 2D materials by classification according to the symmetry of materials. Graphene holds six-fold symmetry (SFS) as the first experimental Dirac material, however, not all 2D materials with SFS have DCs. To find the conditions of DC formation of 2D materials with SFS, in this work, we used the 'divide-and-couple' scheme combining with TB to explain how SFS led to the formation of DC. A simple general '6n + 2' rule was proposed to judge whether the 2D material with SFS has DCs. Moreover, according to the '6n + 2' rule, a '3n + 1' rule was derived as the condition for the DC formation in graphene-like silagraphene with SFS and used to design a silagraphene Si 6 C 8 with DCs that has only SFS but no mirror symmetry.

Results and discussion
2.1. '6n + 2' rule 2.1.1. Restricted conditions To reduce the complexity, this work only considered the 2D materials containing the elements of the IV main group. All atoms of the system are located in the same plane, and the bands near the Fermi level are formed by the p z orbits of the atoms without the participation of other orbitals. So, if there is no special explanation, the wave functions mentioned in this work refers to the p z orbital wave functions.
Because there are four electrons in the outermost shell of the IV main group elements, they can occupy half of the s orbital and the p orbitals (four orbitals in total). The onsite energies corresponding to these four orbitals are very close, so the Fermi level is usually near the onsite energies after couplings. The p z orbitals are not coupled with the other three orbitals, so half of the bands formed by p z orbitals are usually above Fermi level and half below Fermi level. So, if the p z orbitals are to form DC, there should be even numbers of p z bands. Therefore, this work only discussed the case with even numbers of p z orbitals (in the case of odd numbers of p z orbitals, DC would not form at Fermi level usually), and it is assumed that the number of p z bands above is equal to that below Fermi level.

Wigner-Seitz cell and the possible positions of atoms
The angle between the two lattice vectors of a 2D system with SFS is 120 • and its Wigner-Seitz cell is a regular hexagon ( figure 1(a)). The Brillouin zone of the reciprocal space is also a regular hexagon ( figure 1(b)). The Dirac point is the degeneracy point of bands, so it is usually located at the high symmetry point of the Brillouin zone. The high symmetry points of 2D materials with only SFS are K point and Γ point ( figure 1(b)). We first discuss the possibility that forming a DC at the K point.
The atoms of the system can be located at different positions in the Wigner-Seitz cell, such as the center, the vertexes, the interior (except for the center), and the edge of the regular hexagon ( figure 1(a)). We will discuss these various situations.

Origin of DC by 'divide-and-couple' approach
Next, we adopt the 'divide-and-couple' approach to analyze the origin of DC. This approach includes the following steps generally: (1) for each type of position, the atom wave functions are combined into the eigenfunctions of the symmetry operator which change the K point into either the K point itself or an integer period of difference when the DC at K point in Brillouin zone is discussed. The symmetry operator is the operation of counter-clockwise rotation of 2π/3, i.e. For the case that there is an atom in the center of the Wigner-Seitz cell, the coordination number of the atom is six due to the SFS of the system. However, this work only discussed the systems composed of the elements of the IV main group, whose coordination numbers generally cannot be six, so this situation is not in the scope of this work. Now, it is considered that the atoms are located inside the Wigner-Seitz cell (except for the center). If there is an atom at B 1 , the SFS will produce another five atoms B 2 -B 6 , as shown in figure 1(a). These six atoms are divided into two groups: B 1 , B 3 , B 5 and B 2 , B 4 , B 6 , each of which possesses three-fold symmetry. Within each group, the wave functions can be combined to form the eigenfunctions of C 1 3 : where |B m n is the atomic wave function at B m site of the nth Wigner-Seitz cell. The possible values of l are 1, 2, or 3. |B odd − l n and |B even − l n are all the eigenfunctions of C 1 3 : Then, the Bloch functions corresponding to|B odd − l n and |B even − l n are written as where R n is the vector of the center point of the corresponding Wigner-Seitz cell, and N is the total number of cells in the system. When the system is implemented with the operation of C 1 3 , the Bloch function|B odd − l k becomes: So |B odd − l k is not the eigenfunction of C 1 3 at a general point of Brillouin zone, but at K point (2π/3a, 2π/ √ 3a), it is still the eigenfunctions of C 1 3 : Similarly, at the K point, Because the possible values of l are 1, 2, or 3, the corresponding eigenvalues are e i 2π 3 , e −i 2π 3 , and 1. There would be no coupling between two Bloch functions at the K point if they possess different eigenvalues.
Next, the vertex atoms are considered. When an atom is located at a vertex of the regular hexagon, there will be an atom of the same type at each vertex due to the SFS of the system (figure 1(a), A 1 -A 6 ). Since a vertex is shared by three hexagons, each Wigner-Seitz cell contains two atoms. Here, let A 1 and A 2 belong to the hexagon with the center point vector being 0, then A 3 belongs to the hexagon with the center point vector being −a, A 4 and A 5 belong to the hexagon with the center point vector being −a−b, and A 6 belong to the hexagon with the center point vector being −b. The coordinate of A 1 is (a/2, a/2 √ 3) and the coordinates of A 2 is (0, a/ √ 3). When A 1 is acted by C 1 3 , A 3 and A 5 (twice actions) can be acquired but they are belonging to different unit cells, so it is impossible to obtain the eigenfunctions of C 1 3 in one unit cell for A 1 atom. The situation is similar for A 2 . Then we write their corresponding Bloch functions directly: where |A 1 n is the atomic wave function at A 1 position of the nth Wigner-Seitz cell, R n is the vector of the center point of the corresponding Wigner-Seitz cell, and N is the total number of cells in the system. The other symbols have similar meanings.
When the system is implemented with the operation of C 1 3 , the Bloch function becomes: At the K point (2π/3a, 2π/ √ 3a)of the Brillouin zone, Similarly, At the K point, there is no coupling between |A 1 k and |A 2 k due to the different eigenvalues corresponding to the eigenfunctions of C 1 3 . Now, we only consider the case that there are two atoms at the vertex of the hexagon of the Wigner-Seitz cell and six atoms in the interior of the hexagon (eight atoms in total). The band structure of the system is formed by the couplings among the Bloch functions |A 1 k , |A 2 k , |B odd − l k , and |B even − l k . At K point, these Bloch functions are all the eigenfunctions of C 1 3 , and they can be divided into three groups according to the eigenvalues. The eigenfunctions with the eigenvalue e i 2π 3 are of group 1 (G 1 ), and they are The eigenfunctions with the eigenvalue e −i 2π 3 are of group 2 (G 2 ), and they are The eigenfunctions with the eigenvalue 1 are of group 3 (G 3 ), and they are Through the above grouping, we can get the following conclusions: (1) at K point, there are no intergroup couplings due to the different eigenvalues of C 1 3 . (2) The intragroup couplings of G 1 and G 2 are exactly same. The reason is that if the functions in G 1 are implemented with the united operations of counter-clockwise rotation of π, i.e. C 3 6 , and time reversal T, they can be transformed into the corresponding functions in G 2 , i.e.
So, SFS and time-reversal symmetry guarantee that the intragroup couplings of G 1 and G 2 are exactly same. Please refer to appendix A for a detailed demonstration.
To explain the formation process of the band structure clearly, we set the TB parameters, then implemented TB calculations only considering the nearest-neighbor hopping to obtain the various intermediate band diagrams through step-by-step couplings as follows. During TB parameter setting, the atoms at A sites are regarded as carbon atoms, and the atoms at B sites are regarded as silicon atoms. Then the hopping energy between the two adjacent B sites is set to be t B = −1.12 eV as the V ppπ value in silicon crystal [32], the hopping energy between a vertex and the nearest B sites and the onsite energy difference are set to be t = −1.64 eV and ΔE AB = −2.63 eV as the values in the 2D materials h-SiC [16], respectively.
First, when none of the couplings are considered, the bands were located at the onsite energies of A atoms and B atoms. As shown in figure 2(a), the orange lines are six bands from the onsite energies of B atoms, and green lines are two bands from the onsite energies of A atoms, where, G 1 contains two orange lines and one green line, G 2 also contains two orange lines and one green line, and G 3 contains only two orange lines.
Second, the intragroup couplings of each group are considered, and the bands formed are shown in figure 2(b). Two orange bands and one green band of G 1 in figure 2(a) form three blue bands in figure 2(b), one of which is located near the Fermi level, and the other two are separated on both sides of the Fermi level. Since the intragroup couplings of G 1 and G 2 are the same exactly, the bands from G 2 are also three blue lines completely coincident with the lines of G 1 . Two orange bands of G 3 in figure 2(a) form two red bands in figure 2(b), which are separated on both sides of Fermi level. The coupling parameters of each group are shown in table 1. The coupling parameters in G 1 and G 2 are same by complex conjugate transformation except for a constant phase factor, so the bands by the intragroup couplings of G 1 and G 2 are the same exactly.
Third, the intergroup couplings are considered among G 1 , G 2 , and G 3 . There is a band gap between the middle band in G 1 and the middle band in G 2 (these two bands are conduction band (CB) and valence band (VB), respectively). However, at K point, because there are no intergroup couplings among G 1 , G 2 , and G 3 , CB and VB are still in contact, then a DC is formed at Fermi level, as shown in figure 2(c).
The above step-by-step coupling analysis is similar to the analysis in our previous work [17]. From the formation process of band structure above, the addition of six atoms to the interior of the Wigner-Seitz cell results in two more bands in each group of G 1 , G 2 , and G 3 . Thus, there are still odd bands in G 1 and G 2 groups, respectively, and the middle band (CB or VB) in each group of G 1 and G 2 is located near Fermi level after considering intragroup couplings. Further couplings lead to a band gap between CB and VB, but they keep in contact at the K point, thus forming a DC. If there are no atoms at the vertexes of the Wigner-Seitz cell, G 1 and G 2 would contain even bands that are separated on both sides of the Fermi level after considering intragroup couplings. Further couplings cannot lead to the formation of DC at the Fermi level. It can be seen that the key to the formation of DC is to make each of G 1 and G 2 contains odd bands.
Next, we discussed the case that there are atoms on the edge of the Wigner-Seitz cell that can be divided into two cases: (1) at the midpoint of the edge; (2) at the other positions. For the second case [the atoms are located at C 1-12 ( figure 1(a))], the couplings are analogous to that of the atoms inside the Wigner-Seitz cell: the number of atoms is multiple of six, and every six atoms in the edge add two Bloch functions to each group of G 1 , G 2 , and G 3 , so as not to affect the formation of DC. See appendix B for details.
Finally, for the case of atoms at the midpoint of the edge, this will lead to an increase of three atoms in the cell (the cell has six edges, each of which is shared by two cells), which results that the whole cell contains odd atoms due to the atoms in various possible positions discussed before being in the form of even numbers. Thus, a cell includes odd p z orbitals, and DC cannot be formed usually, so this case is not within the scope of this work.

Presentation and verification of '6n + 2' rule
From the above analysis, it can be seen that the possible positions of atoms in the Wigner-Seitz cell are vertexes (A), interior except for the center (B), and the edges except for the midpoint (C) of the Wigner-Seitz cell. In a unit cell, A sites contain two atoms, and the numbers of atoms at B and C sites are all a multiple of 6. The condition that the system has a DC at K point is that each of G 1 and G 2 contains an odd number of Bloch functions, which requires that there must be atoms at A sites, but there is no requirement for the existence of atoms at B or C sites. Therefore, when the system has a DC at the K point, the total number N of atoms in a Wigner-Seitz cell is as follows: Equation (21) can be regarded as the criterion of whether the system with SFS has DC at the K point in the Brillouin zone, dubbed a '6n + 2' rule. When the rule is satisfied, the system has a DC at the K point, and when the rule is not satisfied, the system has no DC at the K point. For example, graphene [2] (N = 2), α-graphyne [8,30] (N = 8), δ-graphyne [11] (N = 20), g-SiC 3 silagraphene [15,17] (N = 8), and g-Si 3 C silagraphene [15,17] (N = 8) with DC at the K point all satisfy the '6n + 2' rule. However, γ-graphyne [33] (N = 12) and β-graphyne [8] (N = 18) do not have DC at the K point because they do not satisfy the '6n + 2' rule. Although β-graphyne has a DC in the GM path, it does not violate our conclusion on how SFS leads to DC at the K point, because the DC of β-graphyne originates from mirror symmetry [8,31].
It should be noted that the premise of the '6n + 2' rule as a criterion for the existence or absence of DC is that the bands near the Fermi level are completely formed by the p z orbitals without the participation of other orbitals. For example, we insert a carbon atom into each carbon chain of δ-graphyne 11 , and the obtained system satisfies the '6n + 2' rule. However, because there are some bands participated by non-p z orbitals near the Fermi level, the system has no clean DCs at the Fermi level. See appendix C for details.

Influence of TB parameters on DC band structure
When the 'divide-and-couple' approach is used to analyze the band structure of the system, it is based on the symmetry of the system, which limits the choices of TB parameters. The 'divide-and-couple' approach reflects the phenomenon that the different sets of TB parameters can give rise to different band structures. Of course, under the limitation of symmetry, the choice of different TB parameters can also have a slight impact on the band structure. For example, in the 'divide-and-couple' approach analysis above, increasing the difference of the onsite energies between B and A a large extent (larger than 5.37 eV) can make the value of the CB at Γ less than that at K so that the DC point deviates from the Fermi-level (see appendix D for details), but this does not affect the existence of DC. Since there are various materials satisfying the '6n + 2' rule, we cannot give the specific range of TB parameters when DC does not deviate from the Fermi level.
However, for the fourth main group materials, when the '6n + 2' rule is satisfied, there will be DC at the Fermi level under the normal case.

Influence of buckled geometry on DC band structure
When the system shows buckled geometry, the bands near the Fermi-level would be coupled with non-p z orbitals. However, if the system maintains six-fold rotation symmetry, the wave functions in G 1 can still be translated to the wave functions in G 2 by the united operations of counter-clockwise rotation of π, i.e. C 3 6 , and time reversal T, so the intragroup couplings of G 1 and G 2 are exactly same, and the VB and CD still degenerate after the intragroup couplings; then the intergroup couplings would open a gap between VB and CB, but VB and CB are still in contact at the K point because there are no intergroup couplings among G 1 , G 2 , and G 3 due to the different eigenvalues of C 1 3 , and the DC is formed. Furthermore, the buckled geometry may break six-fold rotation symmetry, however, if the system possesses three-fold rotation symmetry and inversion symmetry (e.g. silicene and germanene), the wave functions in G 1 can still be translated to the wave functions in G 2 by the united operations of inversion and time reversal T, so the DC still can be formed.

Possibility of the bands at Γ point being degenerated at the Fermi level
So far we discussed the formation of DC at the K point in the Brillouin zone. As for the Γ point, we can conclude that there are usually no degenerated bands at the Fermi level by using a similar analysis. See appendix E for detailed analysis.

'3n + 1' rule
As an example of the application of the '6n + 2' rule, this work designed a silagraphene DC material with SFS and honeycomb structure. The materials can be obtained by replacing some C atoms with Si atoms in graphene. First, the lattice vectors of graphene are redefined, one of which is l 1 a + l 2 b (a, b are the two original lattice vectors of graphene, and l 1 and l 2 are integers). Then, some C atoms are replaced by Si atoms while maintaining the SFS of the system, and then the atomic structure of the system is optimized by density functional theory (DFT).
Next, we discussed the conditions that l 1 and l 2 should satisfy so that the designed silagraphene material has DCs at K point. Assuming that the number of atoms in a unit cell of the new system is N. The area occupied by a single atom in the no optimized system is equal to that of graphene, so |l 1 a + l 2 b| 2 : N = |a| 2 : 2 (22) i.e. N = 2(l 2 1 + l 2 2 − l 1 l 2 ).
According to equation (24), dubbed '3n + 1' rule, the possible values of integer l 1 and l 2 are not arbitrary, so not all graphene-like silagraphene with SFS have DC at the K point. On the other hand, the possible values of integer n are not also arbitrary, which shows that the possible numbers of atoms of graphene-like silagraphene with DC at the K point cannot cover all the numbers of atoms satisfying the '6n + 2' rule.
Next, we discussed possible DC materials by using the '3n + 1' rule. When n = 0, then l 1 = 1, l 2 = 0 (or other equivalent solutions), and the corresponding material is graphene, which can be regarded as a special kind of silagraphene. When n = 1, then l 1 = 2, l 2 = 0 (or other equivalent solutions), and the corresponding material is g-SiC 3 or g-Si 3 C, which was theoretically predicted to have DC [15,17]. When n = 2, then l 1 = 1, l 2 = 3 (or other equivalent solutions), the Dirac silagraphene material Si 6 C 8 was designed in this work. When n = 3, equation (24) has no integer solutions; and when n = 4, equation (24) has integer solutions (such as l 1 = 4, l 2 = 1). Therefore, the silagraphene with SFS and honeycomb structure does not have DC at the K point in the Brillouin zone if its atomic number N in a unit cell is in the open interval 14 < N < 26.

DFT calculations
For confirmation, we optimized the atomic structure of Si 6 C 8 by DFT ( figure 3(a)) and calculated the band structures and density of states ( figure 3(b)). It was confirmed that Si 6 C 8 does have DC at the K point. Si 6 C 8 has only SFS, but no mirror symmetry, which indicates that the formation of DC does originate from SFS. The stability of Si 6 C 8 is proved by its phonon spectrum without imaginary frequency. See appendix F for details.

Conclusions
In this work, by using the 'divide-and-couple' scheme [17], the condition of DC at the K point in Brillouin zone caused by the SFS of 2D materials formed by the IV main group elements is analyzed, that is, the total number of atoms in a unit cell should satisfy '6n + 2' rule, with the conditions of the restriction that there are no atoms at the center and the midpoint of the edge of the Wigner-Seitz cell. This rule holds not only when the system presents a flat geometry, but also even if the system presents a buckle geometry as long as the system possesses six-fold rotation symmetry. Furthermore, when the buckle geometry breaks six-fold rotation symmetry, but as long as the system still has three-fold rotation symmetry and inversion symmetry, the '6n + 2' rule still holds.
The '6n + 2' can be used to explain the origins of DC formation in various graphyne and silagraphene materials, and can also guide the designs for new DC materials. According to this rule, we further derived the condition that the graphene-like silagraphene with SFS has a DC at the K point, '3n + 1' rule, by which, the Dirac material Si 6 C 8 is designed. This material has only SFS, but no mirror symmetry, which indicates that the formation of DC does originate from SFS.
Also, the focus of this study was limited to the fourth main group elements. If the electron filling of bands was carefully considered, the analytical method discussed in this work can be extended to analyze whether the 2D materials formed by non-fourth main group elements exhibit DCs.

Computation methods
In the DFT calculations, the Vienna ab initio simulation package (VASP) [34,35] was used at the GGA-PBE [36] level, and the projector augmented-wave method [37] was adopted for the pseudopotentials. The energy cutoff for plane-wave basis expansion was 750 eV. The k-point sampling adopted an (11 × 11 × 1) Monkhorst-Pack mesh for most self-consistent field (SCF) iteration calculations. A 5.0 × 10 −7 eV per atom energy tolerance was adopted in SCF iteration. In the geometry optimization with a conjugated gradient method, the distance between two adjacent atom layers was fixed to be 20 Å to avoid mirror image interaction at periodic boundary conditions, the other cell parameters and atom parameters were all relaxed, and the convergence criteria was that the remaining force on each atom is lower than 0.001 eV Å −1 .

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Demonstration on why the intragroup couplings of G 1 and G 2 are are exactly same
The system possesses SFS and time-reversal symmetry, resulting that the intragroup couplings of G 1 and G 2 are exactly same. The reason is that if the functions in G 1 are implemented with the united operations of counter-clockwise rotation of π, i.e. C 3 6 , and time reversal T, they can be transformed into the corresponding functions in G 2 , i.e. Similarly:

Appendix B. Influence of the edge atoms excluding the midpoint atoms of Wigner-Seitz cell on the formation of Dirac cone
Due to the SFS and translational symmetry, if there is one atom at C 1 on the edge (except for the midpoint), an additional 11 atoms will be generated at C 2 -C 12 , as shown in figure 1(a) in the manuscript. Among the 12 atoms, C 7 -C 12 can be obtained by translating C 1 -C 6 into integer lattice vectors, so each Wigner-Seitz cell contains six atoms. Here, let C 1 -C 6 belong to the hexagon with the center point vector being 0. Similar to the treatment of atomic wave functions at the position B 1 -B 6 in the manuscript, these six atomic wave functions are combined into the eigenfunctions of the operator C 1 3 (which counter-clockwise rotates the wave function by 2π/3): The six Bloch functions expressed in equations (B.5) and (B.6) can be divided into the G 1 , G 2 , and G 3 groups described in the manuscript according to their eigenvalues for C 1 Figure D1. (a) Band structure with the difference between the onsite energies of B and A is increased to 5.37 eV. (b) Band structure with the difference between the onsite energies of B and A is increased to 7.00 eV. ΔE BA = 7.00 eV, the value of CB at Γ is 0.262 eV lower than the value of DC point, which moves the DC point 0.127 eV upward from the Fermi level (as shown in figure D1(b)). Our previous articles have similar discussions [17]. no matter whether there are atoms at A points, each of G 1 and G 2 contains even numbers of Bloch functions, then there are usually no degenerated bands at the Fermi level at Γ point.

Appendix F. Phonon spectrum of Si 6 C 8
To investigate the stability of Si 6 C 8 designed in the manuscript, its phonon spectrum is calculated by PHONOPY using the finite displacement method [38]. A supercell of 2 × 2 × 1 is built first, and the scf-cycle is run with VASP [34,35] for the fore calculations. The k-point sampling for the supercell calculations adopts a 7 × 7 × 1 Monkhorst-Pack mesh. Other calculation parameters can be found in the chapter 'Computation methods'. The phonon spectrum of Si 6 C 8 is shown in figure F1, and there is no virtual frequency found, proving its stability.