Flat bands without twists: periodic holey graphene

In recent years, the concept of electronic flat bands has gained prominence in the fields of materials science and condensed matter physics [1–10]. This heightened interest can be attributed to the electronic, geometrical and topological properties inherent to flat bands [10–13] and their potential to give rise to novel phenomena [7, 14–20]. Specifically, the discovery of superconductivity in twisted bilayer graphene in 2018 fueled considerable attention in the field [7]. However, the primary challenge lies in the experimental control of the rotation angle. For this reason, alternative systems are currently being explored to preserve flat bands without relying solely on the rotation angle. Examples of these systems include graphene subjected to mechanical deformations [21–29] , multilayer graphene twisted (MTLG) under pressure [30–34], cyclicgraphyne, cyclicgraphdiyne [35] and holey graphene (HG) [36–39], where the theoretical possibility of chiral superconductivity emergence has been demonstrated [40].

HG, also called graphene nanomesh, refers to the presence of holes distributed in the material [37, 41–43]. Introducing holes in graphene provides an additional degree of freedom to control optoelectronic and mechanical properties [44–46]. For example, holes create a mobility gap in graphene making it a narrow gap semiconductor [44]. This opens opportunities for exploring electronic, optical and correlated dependent phenomena and designing graphene-based devices that can exploit the unique properties of HG [42, 47–50]. Furthermore, the size and distribution of the holes can also modify the electronic properties. Using DFT calculations, Barkov and Glukhova [51] investigated HGs charge transmission in periodic arrays of holes.

Several techniques are employed to fabricate HG, including electron/ion beam lithography [47, 48, 52], plasma etching [47, 48], chemical etching [37], laser-based methods [53–58], and others [42, 43, 47–49, 59–63] . These techniques offer an experimental advantage over MTLG. Among them, the most common techniques used to fabricate HG are:

• Top-down lithography involves the selective removal of graphene sheets to create nanopores with uniform sizes and spatial arrangement. Although successful in producing graphene sheets with uniform nano-size pores, this technique faces challenges in large-scale fabrication and requires specialized equipment and expertise, making it expensive [47, 48, 52].
• Template-assisted chemical vapor deposition [42, 48, 52, 61]. 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 [61]. 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.
• Block copolymer lithography techniques on CVD-grown graphene [59, 64]. 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 [53–58]. 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 is 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 HG which we show are able to produce flat-bands. From there, we study its electronic properties. The structure of this work is as follows: in section 2, we present the tight-binding model and the unit cells of the HG primarily under discussion. In section 3 we present the electronic properties of the systems introduced in section 2. Therein, the emergence of flat bands and energy gaps is discussed in terms of the Brillouin zone (BZ) folding and bipartite sublattices site imbalance. In section 4 we study the topological properties and present an effective low-energy model. Finally, our conclusions and future perspectives are outlined in section 5.

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 figure 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.

Zoom In Zoom Out Reset image size

Additionally, we use two orthogonal lattice vectors: $\boldsymbol{l}_{1} = ( \sqrt{3} a_0,0)$ and $\boldsymbol{l}_{2} = (0,3a_0)$, where $a_0 = 1.42$ Å is the carbon-carbon distance. As visually depicted in figures 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. We define the unit cell by employing lattice vectors $\boldsymbol{a}_{1} = n \boldsymbol{l}_{1}$ and $\boldsymbol{a}_{2} = m \boldsymbol{l}_{2}$ (as shown in figure 1(b), where $n \in \mathbb{N}$ to ensure periodicity in the x direction, and similarly, $m \in \mathbb{N}$ is employed to establish periodicity in the y direction.

After obtaining the graphene lattice, we selectively remove atoms located within circles of radius R (Å) (refer to figure 1(c). It is worth noting that, without loss of generality, the circle delineating the atoms for removal is centered at $(\boldsymbol{a}_{1}+\boldsymbol{a}_{2})/2$ in the unit cell. As depicted in figure 1, we obtain a lattice featuring holes. We will denote the unit cell of graphene with a hole of radius R and vectors $\boldsymbol{a}_{1} = n \boldsymbol{l}_{1}, \boldsymbol{a}_{2} = m \boldsymbol{l}_{2}$ as UCHG$(n,m,R)$.

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

$$\begin{equation} \mathcal{H} = \sum_{\langle ij \rangle} t_{0} \hat{c}^{\dagger}_{\boldsymbol{r}_{i}}\hat{c}_{\boldsymbol{r}_{j}} + h.c \end{equation} \tag{ 1 }$$

where $\langle ij \rangle$ represents the sum over the neighbors with positions r _i and r _j that satisfy $|\boldsymbol{r}_{j}-\boldsymbol{r}_{i}| = a_0$; $\hat{c}^{\dagger}_{\boldsymbol{r}_{i}} (\hat{c}_{\boldsymbol{r}_{i}})$ is the creation (annhilation) operator and $t_0 = -2.8$ eV is the hopping integral between the ith and jth 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 $\hat{\mathcal{H}}_{T}(\boldsymbol{k})$. The eigenvalues and eigenfunctions were thus numerically found by using Python dedicated libraries [65–67]. In the following section we will discuss the resulting electronic properties.

This section is divided in three subsections, the first is devoted to present the band structure, the second analyzes the mechanism for flat band formation and the third explains the observed periodicity of gap formation.

3.1. Band structure and localization properties

In this section, we will discuss the electronic 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 and edge terminations. Observe that among our examples, in figure 2(f) we include the case of a hole with many sites of coordination one, akin to have dangling bonds, that will allow to discuss the origin of flat bands.

Zoom In Zoom Out Reset image size

In figure 3 we illustrate the bands obtained for each of the systems shown in figure 2. The BZ and high symmetry points are detailed in figure 4(c). 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. More important for our purpose is that flat bands appear at energy E = 0, corresponding to the Fermi energy ($E_\mathrm{F} = 0$ eV) at half-filling. Additional plots are provided in the appendix, showing that such flat bands are also seen for bigger sizes of the system. Recent research suggests that the formation of flat bands is intricately linked with compact localized states (CLS) [10], also known as confined states [68, 69], 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, 70–73]

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

$$\begin{equation} \mathrm{IPR}\left(E\left(\boldsymbol{k},s\right)\right) = \int_{\mathrm{UC}} |\psi_{\boldsymbol{k},s}\left(\boldsymbol{r}\right)|^{4} d^{2} r \end{equation} \tag{ 2 }$$

where $\psi_{\boldsymbol{k},s}(\boldsymbol{r})$ represents the wave function in the sth band, characterized by wave vector k and energy $E(\boldsymbol{k},s)$. Due to the periodic conditions of the lattice, the integral that appears in equation (2) is performed over the unit cell (UC). The IPR of an extended state goes as $1/N$ where N is the number of sites, while for localized states does not depend on N. For exponentially localized states, it can be proved that $\mathrm{IPR}(E) \sim \lambda^{-2}$ where λ is the localization length.

In figure 3, the $\mathrm{IPR}(E(\boldsymbol{k},s))$ values for each $E_{\boldsymbol{k,s}}$ are projected onto the band structure using a color code. In the appendix, the IPR for bigger system sizes is also presented. As anticipated, the flat band states seem to exhibit a lesser participation in contrast to the majority of other states. For example, in the case UCHG(5,5,3) the systems has N = 90 sites. An extended state will have $IPR = 1/90 = 0.01$. Such value is in perfect agreement with the numerical results seen in figure 3(a) for non-flat bands. Flat-bands in figure 3(a) have $\mathrm{IPR} \approx 0.03$. How we interpret this value? First we remark here that the system is periodic even if holes are present, thus all states are in principle extended and will scale as $1/N$, but with a slightly different participation. Secondly, we observe that flat band states preferentially spatially localize in the zig-zag boundary regions that form around the hole, as can be observed in figures 4(a) and (b). This allows to estimate the $\mathrm{IPR}$ as follows. The number of sites in the hole edge seen in figure 3(a) is $N_\mathrm{hole} = 22$, from where the fraction $N_\mathrm{hole}/N = 22/90 = 0.24$. Taken this number as the fraction of active sites, we estimate that $\mathrm{IPR} \approx N_\mathrm{hole}/N = 0.24$ in reasonable agreement with the numerically observed value. The difference from the numerical value arises due to the fact that not only the sites at the hole's edge contribute; there are also sites close to the edge that contribute. In spite of this, the message here is that the flat-bands IPR shows a tendency for a decrease in their spatial participation when compared with others.

It is worth noting that these flat-band states are not the sole examples of increased IPR; in addition, we find quasi-flat bands within specific regions of the BZ that show a similar behavior, 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, 75].

To further confirm the nature of the states, in figures 4(a) and (b), we present the site contribution of the local density of states (DOS) for a flat band state in two different examples. It can be seen therein that participation in the flat band mainly occurs at the edges of the holes. From figures 3(a)–(d) and 4(a) and (b), it is observed that if n increases while R remains constant, or vice versa, the participation in the flat band decreases. However, for n = 9 and r = 5 such rule does not to hold since in addition, a gap opens.

3.2. Flat band formation mechanism: sublattice imbalance number and broken path-exchange symmetry

To gain a better understanding on the formation of flat bands, we introduce two parameters. The first one is the imbalance number, denoted as $\delta N(n,R)$ and defined by,

$$\begin{equation} \delta N\left(n,R\right) = |N_{A}-N_{B}| \end{equation} \tag{ 3 }$$

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, $\varrho(n,R)$,

$$\begin{equation} \varrho\left(n,R\right) = \delta N\left(n,R\right)/N, \end{equation} \tag{ 4 }$$

where N is the total number of atoms in each unit cell. The reason to define such quantities is that sublattice site imbalance in bipartite lattices is known to induce flat bands of confined states in other systems [76]. Such confinement is due to destructive interference [10]. Bipartite sublattice site imbalance can only occur in our periodic HG 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 DOS at E = 0 and the observed number of flat bands. Notice that in the last row of table 1 we included a case in which the imbalance number was made big by performing a tailor-made cut that maximizes the number of atoms in one sublattice along the hole rim. As seen in figure 2(f), this is equivalent to have dangling bonds or sites with coordination one in the rim.

Table 1. table displaying the number of imbalances, $\delta N(n,R)$, the imbalance density, $\varrho(n,R)$, gap size, Δ, state density, $\rho(E)$, and the number of flat bands for graphene with a gap of size R [Å] and unit cells UCHG$(n,n,R)$. Notice that in the last row, the imbalance number is big as the cut was chosen to increase the number of atoms in one sublattice along the hole rim. This is equivalent to have dangling bonds or sites with coordination one as seen in figure 2(f).

n	R[Å]	$\delta N(n,R)$	$\varrho(n,R)$	Δ [eV]	$\rho(E = 0)$	No. Flat Bands
5	3	1	${\sim}1.15\times 10^{-2}$	${\sim}0.55$	${\sim}0.119$	1
7	3	1	${\sim}5.46 \times 10^{-3}$	${\sim}0.37$	${\sim}0.060$	1
9	3	1	${\sim}3.21 \times 10^{-3}$	${\sim}0.00$	${\sim}0.0364$	1
7	5.2	1	${\sim}6.30 \times 10^{-3}$	${\sim}0.37$	${\sim}0.0684$	1
9	5.2	1	${\sim}3.50 \times 10^{-3}$	${\sim}0.00$	${\sim}0.0402$	1
9	5	7	${\sim}2.39 \times 10^{-2}$	${\sim}0.53$	${\sim}0.024$	7

From table 1 we observe that the number of flat bands perfectly matches the imbalance number for all the studied cases. Noteworthily, the last row of table 1, which contains the tailor-made site imbalance, produces the same number of flat bands. Therefore, we conclude that flat bands are produced by the sublattice imbalance and thus we can generate flat bands at will. Moreover, it has been demonstrated that flat bands appear as a consequence of brakening the path-exchange symmetry [10, 74]. Upon breaking, this symmetry introduces a complex phase [10, 74], 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 imbalance number is related to the path-exchange symmetry, and consequently the number of flat bands is equal to the imbalance number. The broken path-exchange symmetry distinguishes the sublattice with the highest number of sites. This is illustrated in figure 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 DOS or $\rho(E)$, at E = 0 is directly proportional to the imbalance density (see table 1). In cases with no dangling bonds, we observe that $\rho(0) \sim 10.5 \varrho(n,R)$, while $\rho(0) \sim \varrho(n,R)$ (see figures 4(f) and (g)) for the lattice with dangling bonds. This can be understood by considering $\rho(0) \sim n_\mathrm{st}/N$, where $n_\mathrm{st}$ is the number of sites supporting the state with E = 0. Therefore, $\rho(0)/\varrho (n,R) \sim n_\mathrm{st}/ \delta N(n,R)$. In the studied cases without dangling bonds, $n_\mathrm{st}$ fluctuates between $9-12$ sites (see for example figures 4(a) and (b)), and $\delta N(n, R) = 1$, maintaining an average ratio of 10.5 (see table 1 and figures 4(f) and (g)). For cases with dangling bonds, the states at E = 0 are strongly localized on sites with coordination one. Then $n_\mathrm{st}$ corresponds to such number of sites and hence $\rho(0) \sim \varrho(n,R)$.

3.3. Gap formation mechanism: periodicity of the BZ folding

Let us know discuss the origin of the gaps. In figure 5 we present the gap size (Δ) evolution for three distinct graphene series with hole radii of $R = 3, 5, 5.2$ Å, in which the standard error of the gap, $\sigma_{M}(\Delta)$, and the standard error of the imbalance density, $\sigma_{M}(\varrho(n,R))$, for the three data series shown are $0.02, 0.01, 0.02$ eV and $5 \times 10^{-3}, 2\times 10^{-3}, 1 \times 10^{-3}$, respectively. Also, for each of these series, the size of the unit cell n is varied within the range $n = 6, \ldots, 15$, with the exception of the first series where $n = 4, \ldots,15$. The initial observation in figure 5 is that for the series featuring hole radii of R = 3 and R = 5.2 Å, a periodicity with the property $\sim n \equiv 0 , (\mbox{mod} 3)$ emerges for n sufficiently large, relative to the specific R value at which the gap reaches zero. This phenomenon is particularly evident in figures 3(c) and (e) where the emergence of Dirac cones at the Γ points is clearly observed.

Zoom In Zoom Out Reset image size

The second observation that comes out from figure 5 is that excluding the semimetallic cases, Δ follows a power law for all cases,

$$\begin{equation} \Delta \sim \varrho\left(n, R\right)^p \end{equation} \tag{ 5 }$$

where p is approximately 0.5166, 0.5797, and 0.6360 for the respective series and their statistical coefficient of determination $\overline{R}^{2}$ are 0.9997, 0.9947, and 0.9997.

In a previous work [44], it was found analytically and numerically that for graphene with random sustituitional disorder $|\Delta| \approx |t_{0}| \sqrt{6 x} \approx 6.85857 x^{1/2}$. The nearly similar exponents indicate that in the end, big holes enter as nearly effective impurities. There are two reasons for this coincidence. The first is that impurities in the split band limit act as holes in the host lattice bands as wave functions are not able to penetrate such energy barriers [77]. The hole array acts as a third sublattice which produces the flat band as explained in section 4. It is akin to graphene with strongly biased circular patches [78]. The second reason for the square root scaling of the gap is that the chiral, electron–hole symmetric Hamiltonian when squared, removes one bipartite sublattice. Then E = 0 is the ground state of a triangular lattice with some underlying frustration [44, 76]. Sublattice imbalance leads to a perturbation on the Dirac cone mechanism of frustration in the triangular lattice [44]. Here the system is periodically perturbed and with some effective self-energy. This mechanism works in many other chiral systems [69].

However, Δ also oscillates. This is due to the periodic hole array perturbation perturbation and its relation with the unit cell size and geometry. In particular, now we will demonstrate that the gap periodicity arises due to the folding of the pristine graphene's BZ into the holey superlattice rectangular BZ (see figure 4(c)). To do this, consider the original Dirac points on each of the graphene's valley,

$$\begin{equation} \boldsymbol{K}_{\xi} = \xi \frac{4 \pi}{3 \sqrt{3}a_{0}} \left(1,0\right), \,\, \xi = \pm 1 \end{equation} \tag{ 6 }$$

in the hexagonal BZ of graphene and ξ is the valley pseudospin index. On the other hand, consider the reciprocal vectors,

$$\begin{equation} \boldsymbol{b}_{1} = \frac{2 \pi}{\sqrt{3} n a_{0}}\left(1,0\right), \hspace{1mm} \boldsymbol{b}_{2} = \frac{2 \pi}{3m a_{0}} \left(0,1\right) \end{equation} \tag{ 7 }$$

of HG with UCHG$(n,m,R)$. It can be demonstrated that if $n = 3 l + \eta$ with $l \in \mathbb{Z}$ and $\eta = -1,0,1$, the following equation is obtained,

$$\begin{align} \boldsymbol{K}_{\xi} & = \xi\left(2l+\eta\right) \boldsymbol{b}_{1} - \xi \frac{\eta}{3} \boldsymbol{b}_{1} = \boldsymbol{G}_{\xi \left(2l+ \eta\right),0}\nonumber\\ & \quad + \boldsymbol{\kappa}_{\xi,\eta}; \,\, \, \boldsymbol{\kappa}_{\xi, \eta} \equiv - \frac{\xi \eta}{3} \boldsymbol{b}_{1} \end{align} \tag{ 8 }$$

where

$$\begin{equation} \boldsymbol{G}_{p,q}\equiv p \boldsymbol{b}_{1}+q \boldsymbol{b}_{2} \,\,\mbox{and } \boldsymbol{\kappa}_{\xi, \eta} \equiv - \frac{\xi \eta}{3} \boldsymbol{b}_{1}, \end{equation} \tag{ 9 }$$

and therefore, $\boldsymbol{K}_{\xi}$ folds onto the point $\boldsymbol{\kappa}_{\xi,\eta}$. In other words, $\boldsymbol{K}_{\xi}$ folds along the path $\boldsymbol{K}_{-}-\boldsymbol{\Gamma}-\boldsymbol{K}_{+}$ with a destination point determined by the η index and a period $\Delta n = 3$. Figure 6) presents, from n = 1 to n = 6, the folding sequence that results from equation (8). Notice that this result is independent of the size index m or hole radii R.

Zoom In Zoom Out Reset image size

It is noteworthy that when η = 0, $K_{\pm}$ 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 $\eta = \pm 1$ [79–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 to break the path-exchange bond symmetry [85, 86]. Therefore, the imbalance density and the folding of the BZ explain the behavior of Δ seen in figure 5.

As is well known, flat bands in Moire systems are often topological. Thus it is interesting to explore the flat band topology of our HG system, and if possible, to build an effective low-energy model. Figure 7 shows the Berry curvature for the unit cells shown in figure 2 as a function of k and the band s. The Berry curvature is denoted as $\Omega_{s}(\boldsymbol{k})$ and defined as,

$$\begin{align} \boldsymbol{\Omega}_{s}\left(\boldsymbol{k}\right) & = \nabla_{\boldsymbol{k}} \times \mathcal{A}_{s} \left(\boldsymbol{k}\right) \nonumber \\ \mathcal{A}_{s}\left(\boldsymbol{k}\right) & = i \langle u_{s \boldsymbol{k}}| \nabla_{\boldsymbol{k}}| u_{s, \boldsymbol{k}}\rangle \end{align} \tag{ 10 }$$

where $\mathcal{A}_{s}(\boldsymbol{k})$ is the Berry connection and $u_{s, \boldsymbol{k}}$ are the periodic functions in the Bloch states that solve the Schrödinger equation. Figure 7 illustrates how in some cases the symmetry breaking induces a nonzero Berry curvature at the folding points $K_{\pm}$ and in the vicinity of E = 0. However, in figures 7(c) and (e) we see that although the Dirac cones are preserved, the bands are non topological. This is due to the protection at the Γ point. Therefore, to have topological bands, it is necessary that n must not be divisible by three. Another option is to add dangling bonds. In summary, the flat bands discovered in this system are generally non topological, and are thus different from the moiré flat bands, which are often topological.

Zoom In Zoom Out Reset image size

As the electronic, optical and topological properties are determined mainly by electrons at the Fermi level, it is very convenient to develop a low-energy model to the system, akin to the Dirac effective equation that occurs in graphene. We start by observing that without dangling bonds and near the high-symmetry points $\boldsymbol{\kappa}_{\eta \xi}$ of the folded reciprocal lattice, 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. In other words, holes can be considered as a third sublattice producing an extra band. Folding the Dirac points at $\boldsymbol{\kappa}_{\xi,\eta}$ 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 $\alpha-\mathcal{T}_{3}$ model [87], defined by

$$\begin{equation} \hat{H}_{0}\left(\alpha\right) = \hbar v_{F} \boldsymbol{\hat{S}}_{\xi}\left(\alpha\right) \end{equation} \tag{ 11 }$$

where ${\hat{\boldsymbol{S}}_{\boldsymbol{\xi}}}(\alpha) = ( \xi \hat{S}_{x}(\alpha), \hat{S}_{y}(\alpha))$ with $0 \unicode{x2A7D} \alpha \unicode{x2A7D} 1$ a parameter, and $\hat{S}_{x}(\alpha), \hat{S}_{y}(\alpha)$, and $v_\mathrm{F}$ are the pseudospin operators and Fermi velocity, respectively; given by

$$\begin{align} v_\mathrm{F} & = \mathfrak{c}\left(n,R\right) \frac{t_{0}a_{0}}{\hbar}\nonumber \\ \hat{S}_{x}\left(\alpha\right) & = \left( \begin{array}{ccc} 0& C_{\alpha} & 0 \\ C_{\alpha} & 0 & S_{\alpha}\\ 0 & S_{\alpha } & 0 \end{array} \right) \nonumber \\ \hat{S}_{y}\left(\alpha\right) & = \left( \begin{array}{ccc} 0& -iC_{\alpha} & 0 \\ iC_{\alpha} & 0 & -iS_{\alpha}\\ 0 & i S_{\alpha } & 0 \end{array} \right) \end{align} \tag{ 12 }$$

with $C_{\alpha} = 1/\sqrt{1+ \alpha^{2}}$, $S_{\alpha} = \alpha/ \sqrt{1+ \alpha^{2}}$ and $\mathfrak{c}(n,R)$ 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 $\boldsymbol{\kappa}_{\xi,\eta}$ will exhibit an additional mass-like term to the $\alpha-\mathcal{T}_{3}$ Hamiltonian, thus the low-energy hamiltonian is

$$\begin{equation} \mathcal{H}_{\eta \xi}\left(\boldsymbol{k}\right) = \hbar v_{F} \hat{\boldsymbol{S}}_{\xi}\left(|\xi|\right) \cdot \boldsymbol{k}+ M_{\tau}\left(\boldsymbol{k}\right) \Sigma_{z} \end{equation} \tag{ 13 }$$

where $\Sigma_{z}$ is a pseudospin operator defined as

$$\begin{equation} \Sigma_{z} = \left( \begin{array}{ccc} 1&0&0 \\ 0&0&0\\ 0&0&-1 \end{array}\right) \end{equation} \tag{ 14 }$$

and $M_{\tau}(\boldsymbol{k})$ is the mass-like term that depends on k and pseudospin index $\tau = \xi \eta = -1,0,1$. It can be constructed from the expansion of energies for the bands s, given by

$$\begin{align} E_{s}\left(\boldsymbol{\kappa}_{\xi,\eta}+\boldsymbol{q}\right) &\approx \langle u_{s, \boldsymbol{\kappa}_{\xi,\eta}}|\hat{\mathcal{H}_{T}}\left(\boldsymbol{\kappa}_{\xi,\eta}\right)| u_{s, \boldsymbol{\kappa}_{\xi,\eta}}\rangle \nonumber\\ & + \langle u_{s, \boldsymbol{\kappa}_{\xi,\eta}}\left|\left.\boldsymbol{q} \cdot \left( \nabla_{\boldsymbol{k}} \hat{\mathcal{H}}_{T}\left(\boldsymbol{k}\right)\right)\right|_{\boldsymbol{k} = \boldsymbol{\kappa}_{\xi,\eta}} \right|u_{s, \boldsymbol{\kappa}_{\xi,\eta}}\rangle\nonumber\\ & + \mathcal{O}\left(|\boldsymbol{q}|^{2}\right) = s\left( |E_{1}\left(\boldsymbol{\kappa}_{\xi,\eta}\right)|+ \mathfrak{a}_{2}\left(n,R\right) q_{x}\right.\nonumber\\ & \left.+ \mathcal{O}\left(|\boldsymbol{q}|^{2}\right)\right). \end{align} \tag{ 15 }$$

where $|E_{1}(\boldsymbol{\kappa}_{\xi,\eta})|, \,\mathfrak{a}_{2}(n,R)$ can be numerically determined. Then $M_{\tau}(\boldsymbol{k})$ is expressed as,

$$\begin{equation} M_{\tau}\left(\boldsymbol{k}\right) = |E_{1}\left(\boldsymbol{\kappa}_{\xi,\eta}\right)|- \tau |\mathfrak{a}_{2}\left(n,R\right)| k_{x}. \end{equation} \tag{ 16 }$$

In figure 8, the electronic bands obtained from the Hamiltonian $\hat{\mathcal{H}}_{T}$ (dashed black, blue and red lines) and the low-energy approximation with the Hamiltonian $\mathcal{H}_{\xi, \eta}$ (equation (13)) (solid black, blue and red lines) are depicted. In general, an excellent agreement is obtained between the numerical and the low-energy effective model.

Zoom In Zoom Out Reset image size

In this paper, we investigated the formation of flat bands in two-dimensional periodic HG (2D HG). Our findings reveal that these flat band states exhibit an increased IPR when compared to other states, particularly their have maximum density along the zigzag edges surrounding the hole. This highlights a connection between flat band formation and CLS, as investigated in prior works [10, 68, 69].

One of the main results established in this work is that flat bands arise as a consequence of an induced sublattice site imbalance, i.e. by having more sites in one of the graphene's bipartite sublattice than in the other. This is equivalent to break the path-exchange symmetry [10, 74]. Therefore, the results derived for the flat bands is completely general for bipartite sublattice site imbalance, but not in general for different folding schemes and hole shapes. Also, it was established that the breaking of both sublattice and inversion symmetries 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 DOS at E = 0, through the imbalance density $\varrho(n,R)$. Additionally, we have explored the existence of semimetallicity in 2D HG with a periodicity of $\Delta n = 3$, arising from the folding of $K_{\pm}$ 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) can generate a non-zero Berry curvature. However, in many cases the flat bands discovered here can also be non topological, and are thus different from the moire flat bands, which are often topological.

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 $\alpha-\mathcal{T}_{3}$ type Hamiltonian and a mass-like term that opens the gap at $\boldsymbol{\kappa}_{\xi \eta}$ 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.

This work was supported by CONAHCyT Project 1564464 and UNAM DGAPA Project IN101924. Abdiel de Jesús Espinosa-Champo is supported by a CONAHCyT PhD fellowship (No. CVU 1007044). The authors acknowledge and express gratitude to Carlos Ernesto López Natarén from Secretaria Técnica de Cómputo y Telecomunicaciones for helping with the high performance computing infrastructure at Instituto de Física, UNAM, where we run our calculations and his valuable support.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix we present the bands for bigger systems than those studied in the main text. Figure 9 presents three representative cases and its corresponding IPR. Flat bands are clearly seen as expected as site imbalance and not size determine its existence. To further confirm this point, in figure 10 we plot the contribution of each bipartite sublattice to the bands. The flat band is clearly sublattice polarized as expected form site imbalance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size