Flat bands without twists: periodic holey graphene

\textit{Holey Graphene} (HG) is a widely used graphene material for the synthesis of high-purity and highly crystalline materials. In this work, we explore the electronic properties of a periodic distribution of lattice holes, demonstrating the emergence of flat bands with compact localized states. It is shown that the holes break the bipartite sublattice and inversion symmetries, inducing gaps and a nonzero Berry curvature. Moreover, the folding of the Dirac cones from the hexagonal Brillouin zone (BZ) to the holey superlattice rectangular BZ of HG with sizes proportional to an integer $n$ times the graphene's lattice parameter leads to a periodicity in the gap formation such that $n \equiv 0$ (mod $3$). Meanwhile, it is shown that if $n \equiv \pm 1$ (mod $3$), a gap emerges where Dirac points are folded along the $\Gamma-X$ path. The low-energy hamiltonian for the three central bands is also obtained, revealing that the system behaves as an effective $\alpha-\mathcal{T}_{3}$ graphene material. Therefore, a simple protocol is presented here that allows obtaining flat bands at will. Such bands are known to increase electron-electron correlated effects. This work provides an alternative system, much easier to build than twisted systems, to obtain highly correlated quantum phases.

Furthermore, the size and distribution of the holes can also modify the electronic properties.Using DFT calculations, Barkov et al. [52] investigated holey graphene's transmission function, T (E).They found that introducing periodic holes can modulate the values of T (E) by changing the distances between holes and maintaining the hole size; thus, the conductance tends to increase farther apart.
Another technique used to fabricate HG is template-assisted chemical vapor deposition [48,51,53,62].In this method, HG is grown on specific inorganic templates.This technique allows for the bottom-up growth of HG and offers control over nanopores' size, shape, and distribution [62].However, these fabrication techniques have certain drawbacks because they require careful control of various parameters such as temperature, pressure, and duration, making the process lengthy, cumbersome, and costly.
A third approach to synthesizing HG is using block copolymer lithography techniques on CVD-grown graphene [60,65].This technique involves patterning the graphene with a copolymer material and using reactive ion etching to create nanopores in the desired arrangement.One advantage of this technique is that it allows for the control of the neck width, which can be important in determining the properties and applications of HG.However, this approach still requires further optimization to achieve large-scale production of HG.
Nevertheless, laser-based methods have recently emerged as a promising technique for synthesizing HG [54][55][56][57][58][59].These methods offer a precise means of controlling the synthesis process of the holey structure.Researchers can adjust parameters like laser intensity and target material to create holes with specific sizes, shapes, and distributions.This level of control enables the production of HG with tailored properties for various applications.Moreover, this method allows for the fabrication of high-purity and highly crystalline HG without the need for catalysts or chemicals.Laser ablation methods are also scalable, making them ideal for large-scale HG production.Furthermore, these techniques permit the customization of HG by introducing dopants and functional groups.This versatility contributes to the development of graphene-based devices with enhanced performance.However, it's worth noting some notable disadvantages, including the cost and complexity of the equipment, limited scalability for large-area synthesis of pristine graphene, and challenges in achieving structural uniformity.
In this work, we explore different configurations of holey graphene which we show are able to produce flat-bands.From there, we study its electronic properties of the resulting lattices.The structure of this work is as follows: in Section 2, we present the tight-binding model and the unit cells of the holey graphene primarily under discussion.In Section 3 we present the electronic properties of the systems introduced in Section 2. Therein, the emergence of flat bands with compact localized states and energy gaps is discussed in terms of the Brillouin zone folding and bipartite sublattices site imbalance.Finally, our conclusions and future perspectives are outlined in Section 4.

Tight-binding holey graphene model
Our model is based on periodically distributed holes within a graphene lattice, maintaining translational symmetry.This enables us to simplify the treatment by leveraging Bloch's theorem.As seen in Fig. 1 a), we use a 4-atom superlattice unit cell containing the usual graphene non-equivalent sites A and B, where A and B refer to carbon atoms on each of the graphene's bipartite lattices, distinguished in the figure by colors blue and red, respectively.
Additionally, we use two orthogonal lattice vectors: l 1 = ( √ 3a 0 , 0) and l 2 = (0, 3a 0 ), where a 0 = 1.42 Å is the carbon-carbon distance.As visually depicted in Fig. 1 a) and b), this specific configuration ensures the alignment of the horizontal axis with the zigzag orientation and the perpendicular axis with the armchair orientation.
After obtaining the graphene lattice, we selectively remove atoms located within circles of radius R (Å) (refer to Fig. 1c).We define the unit cell by employing lattice vectors a 1 = nl 1 and a 2 = ml 2 (as shown in Fig. 1b), where n ∈ N to ensure periodicity in the x direction, and similarly, m ∈ N is employed to establish periodicity in the y direction.It is worth noting that, without loss of generality, the circle delineating the atoms for removal is centered at (a 1 + a 2 )/2.As depicted in Fig. 1, we obtain a lattice featuring holes; and cutting exclusively in the horizontal or vertical direction results in nanoribbons with zigzag or armchair boundaries, respectively (refer to Fig. 1d and e).We will denote the unit cell of graphene with a hole of radius R and vectors Figure 1.a) Unit cell of graphene with a 4-atom basis and lattice vectors l 1 = ( √ 3a0), l 2 = (0, 3a 0 ) arranged such that the zigzag and armchair directions align with the horizontal and vertical axes, respectively.b) Unit cell for hollow graphene, featuring lattice vectors a 1 = nl 1 and a 2 = ml 2 , where n, m ∈ N. Additionally, the atoms to be removed are highlighted in red, situated within a red circle of radius R. For notation purposes, we will designate this unit cell as UCHG(n, m, R). c) Two-dimensional holey graphene.d) and e) Holey graphene nanoribbons with zigzag and armchair edges, respectively, denoted as ZHGN(n, m, R) and AHGN(n, m, R).Except in the case of b), sites A are denoted in blue, while sites B are represented in red UCHG(n, m, R), and correspondingly, the nanoribbons with zigzag and armchair directions as ZHGN(n, m, R) and AHGN(n, m, R), respectively.
Under such considerations, our investigation is based on a tight-binding model with only first-neighbors hopping transfer integral, which yields a Hamiltonian given by where ⟨ij⟩ represents the sum over the neighbors with positions r i and r j that satisfy the hopping integral between the i-th and j-th sites.Additionally, we numerically construct the Hamiltonian operator in reciprocal space k, which depends on the number of atoms in the unit cell and is denoted as ĤT (k).The eigenvalues and eigenfunctions were thus numerically found by using python dedicated libraries [66][67][68].
In the following section we will discuss the resulting electronic and optical properties.

Electronic properties of bidimensional Holey Graphene
In this section, we will discuss the electronic and optical properties of systems with unit cells of type UCHG(n, n, R), focusing primarily on those depicted in Figure 2, namely, UCH(5, 5, 3), UCHG(7, 7, 3), UCHG(9, 9, 3), UCHG(7, 7, 5.2), UCHG(9, 9, 5.2), and UCHG (9,9,5).These examples were chosen to include the effects of different radius , it is noted that as n increases while maintaining R, the localization decreases; however, this does not hold true for case f).This behavior can be understood in connection with the inbalance density ρ(n, R), as discussed in the text.On the other hand, a consequence of maintaining dangling bonds distributed with C 3 symmetry is the creation of a gap.This gap remains strongly protected even when increasing the size of the cell.
and edge terminations.Observe that among our examples, in Fig. 2 f) we include the case of a hole with dangling bonds.In Fig. 3 we illustrate the bands obtained for each of the systems shown in Fig. 2. As we can see, gaps are observed in some cases while Dirac cones are seen in others.The opening of these gaps will be discussed later on in this section.Meanwhile, we observe that at energy E = 0 flat bands are obtained in all cases, corresponding to the Fermi energy (E F = 0 eV) at half-filling.Recent research suggests that the formation of flat bands is intricately linked with Compact Localized States (CLS) [10], also known as confined states [70,71], signifying the existence of a non-trivial localization behavior.
In order to ascertain whether there exists differences in the localization properties of states within the flat bands and those outside of them, we employ the Inverse Participation Ratio (IPR) [44,[72][73][74][75] where ψ k,s (r) represents the wave function in the s-th bands and characterized by momentum k and energy E(k, s), and due to the periodic condition of the lattice, we perform the integral over the unit cell (UC).The IPR of an extended state goes as 1/N where N is the number of extended states, while for localized states does not depend on N .For exponentially localized states, it can be proven that IP R(E) ∼ λ −2 where λ is the localization length,.
In Figure 3, the IP R(E(k, s)) values for each E k,s are projected onto the band structure using a color code.As anticipated, the flat band states exhibit a significantly high degree of localization in contrast to the majority of other states.Flat band states preferentially spatially localize in the zig-zag boundary regions that form around the hole, as can be observed in Figure 4 a) and b).However, it is worth noting that these localized states are not the sole examples; in addition, we find states for which the bands exhibit quasi-flat behavior within specific regions of the Brillouin zone, notably those with energies at E = |t| and momenta between the X and M points.Such behavior is due to the presence of Van Hove singularities in which a dimerization effect is seen [8,29,77].
To further confirm the nature of localization, in Fig. 3 a)-b) we present the local density of states (LDOS) for a flat band state as a function of the sites for two different examples.Therein it can be seen that localization mainly occurs at the edges of the holes as expected.Also, from Figs. 3 a)-d) and 4 a) and b), it is observed that if n increases while R remains constant, or vice versa, the localization in the flat band decreases.However, for n = 9 and r = 5 such rule does not to hold since in addition to a increased localization, a gap opens.To gain a better understanding on the formation of flat bands, the gap size, ∆ and the density of states at E = 0, we introduce two parameters.The first one is the imbalance number, denoted as δN (n, R), given by which indicates the difference between number of sites A and B in the unit cell UCHG(n, n, R).The second is the imbalance density, ϱ(n, R), given as where N is the total number of atoms in each unit cell.
As established in Ref. [78], sublattice site imbalance in bipartite lattices can induce flat bands of confined states, and in certain cases, the formation of energy gaps and pseudogaps.Bipartite sublattice site imbalance will occur here only at hole edges.To see if such effect is in play , in Table 1 we present the values for the site imbalance number, the site imbalance density, the energy gap, the density of states (DOS) at E = 0 and the observed number of flat bands.[10,76].f)g)The density of states, ρ(E), is presented for the unit cells f)UCHG(5, 5, 3) (black line), UCHG(7, 7, 3) (blue line), UCHG(9, 9, 3) (orange line), g) UCHG(7, 7, 5.2) (black line), and UCHG(9, 9, 5.2) (red line).It is observed that ρ(E = 0) decreases proportionally to the size of the cell n, which is a consequence of the imbalance density ϱ(n, R) as established in Table 1 and discussed previously.
From Table 1, we conclude that the number of flat bands formed at energy E = 0 eV coincides with the imbalance number for all the studied cases.In recent studies, it has been demonstrated that the formation of degenerate flat bands is also a consequence of a broken path-exchange symmetry [10,76].Upon breaking, this symmetry introduces a complex phase [10,76], leading to the lifting of the degeneracy while simultaneously preserving the flatness of the band.In the periodic holey system presented here, it can be demonstrated that the number of imbalances is related to the path exchange symmetry, and consequently the number of flat bands is equal to the imbalance number.Specifically, when δN (n, R) ̸ = 0 such symmetry is broken and an isolated flat band emerges.Therefore, such symmetry distinguishes the sublattice with the highest number of sites.This is illustrated in Fig. 4, panels d) and e), for the case of HG(9, 9, 5.2), where the nonzero energy bands result from equal contributions from sublattices A and B. Additionally, it is worth noting that the density of states (DOS or ρ(E)), at E = 0 is directly proportional to the imbalance density (see Table 1).In cases where there are no dangling bonds, we observe that ρ(E = 0) ∼ 10.5ϱ(n, R), while in cases with dangling bonds, we find that ρ(E = 0) ∼ ϱ(n, R) (see Fig. 4 f) and g)).This can be understood by considering the DOS for each unit cell, in this case, ρ(E = 0) ∼ n st /N , where n st is the number of sites supporting the state with E = 0. Therefore, ρ(E)/ϱ(n, R) ∼ n st /δN (n, R).In the studied cases without dangling bonds, n st fluctuates between 9 − 12 sites (see for example Fig. 4 a) and b)), and δN (n, R) = 1, maintaining an average ratio of 10.5 (see Table 1 and Fig. 4 f) and g) ).For cases with dangling bonds, the E = 0 state is strongly localized on these dangling sites, and n st corresponds to this number of sites, hence ρ(E = 0) ∼ ϱ(n, R).
Let us know discuss the origin of the gaps.In Fig. 5 we present the gap sizes for three distinct graphene series with hole radii of R = 3, 5, 5.2 Å, respectively.For each of these series, the size of the unit cell,n, is considered within the range of values n = 6, . . ., 15, with the exception of the first series where n = 4, . . ., 15.The initial observation of note in Fig. 5 is that, for the series featuring hole radii of R = 3 and R = 5.2 Å, a periodicity with the property ∼ n ≡ 0, (mod3) emerges for n sufficiently large values relative to the specific R value at which the gap reaches zero.This phenomenon is particularly evident in Figure 3 c) and e), where the emergence of Dirac cones at the Γ points is clearly observed.To understand the origin of the gaps, we rely on Fig. 5, and excluding the semimetallic cases, it is obtained that in the three series the gap size, ∆, follows a power law as a function of the imbalance density ∼ ϱ(n, R) p , where p is approximately 0.5166, 0.5797, and 0.6384 for the respective series.These values, close to p = 1/2, are comparable to the case where there is a density of impurities breaking the sublattice symmetry, and the imbalance density is given by x imp .In a previous work [44], it was found analytically that |∆| ≈ |t 0 | √ 6x ≈ 6.85857x 0.5 .The nearly similar exponents indicate that in the end, big holes enter as effective impurities and that the mechanism of underlying frustration due to sublattice imbalance is in play here [44].
On other hand, we will now demonstrate that gap periodicity arises due to the folding of the Dirac cones from the pristine graphene hexagonal Brillouin zone to the holey superlattice rectangular Brillouin zone (see Fig. 4 c) and from broken symmetries.To do this, consider the original Dirac points on each of graphene's valley, in the hexagonal Brillouin zone of graphene and ξ is the valley pseudospin index.On other hand, the reciprocal vectors of holey graphene with UCHG(n, m, R).It can be demonstrated that if n = 3l + η with l ∈ Z and η = −1, 0, 1, the following equation is obtained, Therefore, K ξ folds onto the point κ ξ,η , i.e., along Γ to X = b 1 /2 or Γ to X ′ = −b 1 /2 path and its corresponding valley pseudospin ξ and showing the periodicity of ∆n = 3 through the η index, and is independent of the size index m or hole radii R. Fig. 6) presents such folding sequence from n = 1 to n = 6.
It is noteworthy that when η = 0, K ± folds to the Γ point, leading to valley degeneracy.Upon introducing holes in graphene (without dangling bonds), sublattice and inversion symmetries are broken, resulting in gaps at the Dirac points for η = ±1 [79][80][81][82][83][84].However, Γ is a point protected by other symmetries, ensuring the preservation of the two degenerate cones when η = 0. To induce a gap at Γ, dangling bonds are necessary, which breaks bond symmetry [85,86].Therefore, the imbalance density, the folding of the Brillouin zone, and symmetry breaking explain the behavior of the gap size shown in Fig. 5.
Additionally, Fig. 7 shows the Berry curvature for the unit cells shown in Fig. 2 as a function of k and the band s.The Berry curvature is denoted as Ω s (k) and defined as, where A s (k) is the Berry connection and u s,k are the periodic functions in the Bloch states that solve the Schrödinger equation.Figure Fig. 7 illustrates how the symmetry breaking induces a nonzero Berry curvature at the folding points K ± and in the vicinity of E = 0.
Finally, near the high-symmetry points κ ηξ of the reciprocal folded lattice and without dangling bond three symmetric bands around E = 0 are formed, one for valence, the flat band, and the conduction band, denoted by s = −1, 0, 1, respectively.Folding the Dirac points at κ ξ,η would, in principle, yield two bands with linear dispersion and a flat band at E = 0. Therefore, an appropriate Hamiltonian for this description is the α − T 3 model [87], defined by where Ŝξ (α) = ξ Ŝx (α), Ŝy (α) with 0 ≤ α ≤ 1 a parameter, and Ŝx (α), Ŝy (α), and v F are the pseudospin operators and Fermi velocity, respectively; given by with ) a constant that depends on unit cell size n and hole size R.Such hamiltonian is the most simple one that has flat bands and under electromagnetic radiation, behaves as a two-or three-level Rabi system with clear optical signatures of flat bands [87].
As discussed previously, breaking inversion and sublattice symmetry results in the opening of a gap.Therefore, the low-energy Hamiltonian around κ ξ,η will exhibit an additional mass-like term to the α − T 3 Hamiltonian, thus the low-energy hamiltonian is where Σ z is a pseudospin operator defined as and M τ (k) is the mass-like term that depends on k and pseudospin index τ = ξη = −1, 0, 1.It can be constructed from the expansion of energies for the bands s, given by where |E 1 (κ ξ,η )|, a 2 (n, R) can be numerically determined.Then M τ (k) is expressed as, In Fig. 8, the electronic bands obtained from the Hamiltonian ĤT (dashed black, blue and red lines) and the low-energy approximation with the Hamiltonian H ξ,η (Eq.11) (solid black, blue and red lines) are depicted.In general a excellent agreement is obtained.

Concluding remarks
In this paper, we investigated the formation of flat bands in two-dimensional periodic holey graphene (2D HG).Our findings reveal that these flat band states exhibit significantly higher localization when compared to other states, particularly in the zigzag edge regions surrounding the hole.This highlights a connection between flat band formation and Compact Localized States (CLS), as investigated in prior works [10,70,71].Moreover, we establish that sublattice site imbalance, achieved through the breaking of path-exchange symmetry [10,76], induces the formation of flat bands, while the breaking of sublattice and inversion symmetry leads to the creation of energy gaps.
Furthermore, we have discussed the influence of unit cell size and hole dimensions on the gap size and density of states (DOS) at E = 0, through the imbalance density ϱ(n, R).Additionally, we have explored the existence of semimetallicity in 2D HG with a periodicity of ∆n = 3, arising from the folding of K ± points to the Γ point, protected by bond symmetry, except when introducing dangling bonds that break this symmetry, resulting in a gap.We demonstrated that the breaking of the discussed symmetries (inversion, sublattice, and bond) generates a non-zero Berry curvature.
Finally, we establish a continuous Hamiltonian model for three bands in cases where the hole radius R ̸ = 0 and there are no dangling bonds.This model consists of an α − T 3 type Hamiltonian and a mass-like term that opens the gap at κ ξη points.
Our work gives a protocol that allows to obtain flat bands at will and thus gives a possible new 2D material to obtain highly correlated quantum phases without twists.

Figure 3 .
Figure 3.The band structure for systems with unit cells of a) UCHG(5, 5, 3), b) UCHG(7, 7, 3), c) UCHG(9, 9, 3), d) UCHG(7, 7, 5.2), e) UCHG(9, 9, 5.2), and f) UCHG(9, 9, 5.0) are shown.It is observed that the states in the flat bands exhibit characteristics of Compact Localized States (CLS).Furthermore, for cases a) to d), it is noted that as n increases while maintaining R, the localization decreases; however, this does not hold true for case f).This behavior can be understood in connection with the inbalance density ρ(n, R), as discussed in the text.On the other hand, a consequence of maintaining dangling bonds distributed with C 3 symmetry is the creation of a gap.This gap remains strongly protected even when increasing the size of the cell.

Figure 7 .
Figure 7.As discussed in Eq. 8, the symmetries broken induce a nonzero Berry curvature at the folding points K ± and in the vicinity of E = 0.

Table 1 .
Table displaying the number of imbalances, δN (n, R), the imbalance density, ϱ(n, R), gap size, ∆, state density, ρ(E), and the number of flat bands for graphene with a gap of size R [Å] and unit cells UCHG(n, n, R).