Edge states and sublattice imbalance of rectangular graphene nanoflakes

The energy states of π-electrons in a rectangular graphene flake with a pair of zigzag edges and a pair of armchair edges were studied using a tight-binding method. It is demonstrated that the known exact solution can be interpreted as the results derived from the analytical solutions for graphene nanoribbons on a brick-type lattice by considering the wavefunctions and sublattice imbalance. Edge states in rectangular graphene flakes with a symmetric structure, whose energies are not exactly zero but approach zero as the size increases, are characterized by discrete allowed modes, whereas the zero-energy state appears in rectangular graphene flakes with an asymmetric structure. By applying the edges states of rectangular graphene flakes, zero-energy states were examined in Y-shaped graphene flakes, which involve the same three zigzag edges as in the triangular graphene flakes. Sublattice imbalance in the model was found to force the edge states, including the undamped mode, into zero energy. Configurations of the allowed modes in the zero-energy state are revealed.

In this study, we deal with rectangular GFs, which is one of systems discussed by Ryu and Hatsugai [8] in their investigation of the edge states in particle-hole symmetric Hamiltonians. Following an analytical study of the spectrum of π-electrons in rectangular GFs by Malysheva and Onipko [13], we considered the wavefunctions and energy levels on the basis of the analytical solutions of two types of graphene nanoribbons by Wakabayashi et al [7] to investigate the relationship between edge states and structural asymmetry. Then, we treated a model derived from zigzag-edged triangular GFs attached to a distinguished sublattice imbalance, and clarified the characteristics of edge states.
We employed the simplest tight-binding Hamiltonian defined by  The allowed wave vectors (k x , k y ) corresponding to the ordinary states are indicated by circles. The independent wave vectors are depicted by red closed circles. Numerals inside the BZ are the values of j y , while the ones outside are those of j x in equation (7). Imaginary wavenumbers ¢ p y are created at ( ) p x j x for j x = 5 and j x = 7 in the yellow region. K and K′ are the Dirac points where the energy is zero. Thin lines indicate constant energy lines with E = ±1.
where † c i and c i are the creation and annihilation operators for a π-electron at site i, respectively; t is the hopping integral between the nearest-neighbor sites 〈i, j〉; and ò is the site energy. Hereafter, the energies will be measured from the value of ò in units of |t|. In this study, zero energy corresponds to ò = 0 for the Fermi energy. This is a common approach to study the behavior of π-electrons near the Fermi energy in graphene. Bery and Fertig [10] demonstrated that the results obtained from tight-binding calculations are comparable to the solutions of Dirac equations for graphene nanoribbons. However, a report by Zarenia et al [19] on triangular and hexagonal GFs showed that the results of both approaches are qualitatively similar but quantitatively different in some cases.

Results and discussion
2.1. Rectagular graphene flakes A rectangular GF consisting of N x × N y hexagons [13,14] is shown in figure 1. Each site of carbon atoms is assigned by a pair of indices (m x , m y ) in a brick-type lattice, where m x = 1, 2, L ,M x and m y = 1, 2, L ,M y , in the same manner as [6,7]. The rectangular GF can also be specified by The total number of sites is We define two independent sublattices, A and B, for m x + m y = odd and m x + m y = even, respectively. Malysheva and Onipko [13] proved that the energy states of rectangular GFs can be described by two wavenumbers perpendicular to each other, k x and k y . For simplicity, these are made dimensionless as follows: where a 0 is the graphene lattice constant. The Brillouin zone of a rectangular GF is in the shape of a rectangle ( < < where integer j x = 1, 2, L ,M x within 0 < p x < π. For each ( ) p x j x , the wavenumbers p y in the y-direction along an armchair-shaped boundary are obtained by solving the following equation [7,13,14]: , equation (8) has M y real solutions: y ) within the range from 0 to π. However, when The other one in the latter case is an imaginary solution = ¢ p ip y y satisfying equation (9), which corresponds to an edge state. Although N wave vectors, including the ordinary states and edge states, are allowed in the Brillouin zone, they are reduced to half their number because one wave vector can describe both energy states in the valence and conduction bands. Consequently, we have only to consider that x j x The existence of the edge states (or imaginary wavenumbers ¢ p y ) depends on the existence of integer j x to satisfy which is derived from = 0 dF dp y according to [7]. Equation (12) can be rewritten as Clearly, the number of edge states increase with increasing M x . For example, those in the rectangular GF with M x = 11 and M y = 8 appear only when j x = 5 as 4.2408 < j x < 6 from equation (14). Those with M x = 21 and M y = 10 appear when j x = 8, 9, 10 as 7.6956 < j x < 11. When M y → ∞ , equation (14) becomes 1 5 x j x which corresponds to the left half of the yellow region in figure 2. The number of j x satisfying equation (15) is It should be noted that an edge state corresponding to ( ) p x j x closer to π/2 (or j x closer to (M x + 1)/2) appears preferentially when M x or M y increase.
The energies and wavefunctions of the ordinary or edge states can be obtained directly by applying the wave vectors ( ) to the π-electron band energy of graphene. The energy and wavefunctions of the ordinary states can be expressed as Those of the edge states are where C and C e are the normalization factors, and σ = − 1 and +1 distinguish between the valence and conduction bands, respectively. C and C e are expressed as follows: , corresponding to the M-point in figure 2. In fact, it can be demonstrated that by taking into account equation (7). On the other hand, the summations of y-components are not always reduced to simple constants because ( ) p y j y and ¢ p y are solutions that satisfy equations (8) and (9), respectively. Furthermore, we can demonstrate that S x is divided into two halves for sublattices A and B regardless of whether M x is even or odd. To distinguish between sublattices A and B, the sites at m y are divided into two groups with m x = even and m x = odd, whose numbers are N and ¢ N , respectively. Therefore, The contribution from A-sites to S x is the same as that from B-sites despite the sublattice imbalance in a row of sites with any m y . From the above, we can derive the various characteristics of the energy spectrum of rectangular GFs. It is significant that there is no wave vector satisfying E = 0 or E e = 0 within the region of equations (10) and (11), as pointed out by Ryu and Hatsugai [8]. Moreover, the existence of peaks at E = ±1 is considered as one of the characteristics of the energy spectrum of clusters made from a graphene sheet [15,30]. The degeneracy at E = ±1 is also evident in an N x × N y rectangular GF, which can be interpreted as the number of wave vectors ( ) figure 2. This is because when p x = π/2 or p y = π − p x , then E = ±1 from equation (17). Therefore, we can predict that the degree of degeneracy is N y + GCF − 1, which is N y at a minimum, or 2N y at a maximum for N x = N y , where GCF is the greatest common factor of N x + 1 and N y + 1.

Decay ratios of edge states
For edge states of graphene ribbons with zigzag edges, the intensity decays at a rate of 2 from the edge [4][5][6][7], that is, ( ) µ I m D y k m y . In rectangular GFs, the rate is discrete using j x in equation (14). However, the decay in rectangular GFs is subjected to which is normalized by the value at the edge, the decay rate at each m y can be interpreted using equation (9) as follows:

Size-dependence of edge state energies
The edge states of rectangular GFs are associated with j x satisfying equation (14). The energy levels expressed in equation (19) do not appear at zero energy. In the rectangular GF with N x , the first nearest energies to E = 0 are associated with j x = N x ; the second nearest ones with j x = N x − 1, and so on. That is, the n-th nearest energies to E = 0, each of which are in the valence and conduction bands, are associated with Clearly, the energies approach zero as the size N x × N y increases. To describe this behavior, we estimated the approximate relation between the energy and size using the energies obtained in numerical calculations for approximately 1500 rectangular GFs, with sizes ranging from N x = 3 to N x = 50 and N y = 2 to N y = 50. The relation is as follows: where D j x involves n and N x through equations (7) and (38). For example, a rectangular GF with N x = 35 and N y = 20 has N R = 11 edge states corresponding to j x from 25 to 35 satisfying equation (14) in the valence or conduction band. Provided that the edge states degenerate substantially when the energy difference is within 10 −8 , the rectangular GF is deemed to have seven degenerate edge states at E = 0, as seen in figure 4. However, the degeneracy is substantial but not essential, as noted previously. It is more significant that the order of the edge state energies n and j x are directly connected to equation (38) in rectangular GFs.

Zero-energy state
We have hitherto considered N x × N y rectangular GFs, where M x is odd and M y is even from equations (2) and (3). In this section, we deal with an asymmetric structure: rectangular GFs with M x = 2N x + 1 = odd and M y = odd, where the number of states N = M x M y is odd differently from symmetric cases. In other words, we consider the case of ΔN = 1, where the sublattice imbalance ΔN = |N A − N B | using the numbers of A-and B-sites N A and N B , respectively. The former is also a pure edge state corresponding to p x = π/2, which cannot satisfy equation (19). The amplitudes are normalized by the maximum or minimum values of the respective wavefunctions. Sites depicted by transparent circles indicate zero probability. Whereas N states are separated into two halves in the valence and conduction bands when N = even, a zeroenergy state must appear inevitably along with two (N − 1)/2 states in the valence and conduction bands, when N = odd. The zero-energy state, which cannot be described by equation (17) or (19), is a pure edge state corresponding to  Figure 6. A Y-shaped GF, where an equilateral triangular GF with zigzag edges connecting three N x × N y rectangular GFs is a model involving three zigzag edges with a size of N x , which are kept away from the sides of the triangular part at a distance of N y .

Y-shaped graphene flakes
To investigate zero-energy states due to sublattice imbalance or an asymmetric structure, we consider a Y-shaped GF consisting of three N x × N y rectangular GFs, and an equilateral triangular GF with a zigzag edge length of N x , as shown in figure 6. This model, which is a triangular GF with three zigzag edges when N y = 0, holds the sublattice imbalance ΔN = N x − 1, and can be proved to have ΔN-fold degenerate states at zero energy, whose wavefunctions produce zero amplitude on sublattice B sites in the same manner as [18]. Therefore, the charge density distribution at zero energy is completely localized on sublattice A sites, and has a tendency to be localized near the zigzag edges. The features can be confirmed in  (7) and (31)). However, this was solved by applying them to the cases with small N x shown in figure 8, where the intensities I(m y ) for the rectangular parts at zero energy in the Y-shaped GFs with N y = 16 are depicted. It is clear that the cases with N x = 2, 3, 4 have a single mode j x , while the cases with N x = 5, 6, 7 consist of two modes. We have estimated that the allowed modes in the cases with larger N x satisfy ) for j x = N S in the case of Y-shaped GFs wit N x = 3ℓ − 1 (ℓ = 1, 2, 3, L). The decomposition results in the cases of N T = 1, 2, and 3 are shown in figure 8 and table 1. The examples of larger N T can be confirmed in figure 9 and table 2. In the tables, it can be seen that the value of I tri is mainly occupied by the component of mode j , because the wavefunctions with smaller D j x are localized Table 1. Decomposition of the number of zero-energy states in Y-shaped GFs with small N x and N y = 16. The size of a rectangular part is denoted by N x × N y . j x and c j x are the mode and coefficient, respectively, in equation (49), I rec in equations (46) is the number of zeroenergy states for a rectangular part, which consists of N T components corresponding to j x and ( ) I j rec x in equation (47). The contribution of a triangular part I tri = ΔN − 3I ( rec) . . nearer the zigzag edges. By compensating the contribution from the triangular part I tri , we can determine how to distribute the number of zero-energy states ΔN to N T allowed modes represented by j x , ranging from N S to N x : The configurations of the allowed modes in the zero-energy states can be represented by the notation: N j x ( j x = N x , L ,N S in descending order) for following the electron configuration of an atom, where numbers of the zero-energy states along with the allowed modes are symbolized by the main part of ( ) p x j x in equation (7). For example, the configurations are (2/3) 1 , (3/4) 2 , (4/5) 3 for N x = 2, 3, 4 (N T = 1); (5/6) 3 (4/6) 1 , (6/7) 3 (5/7) 2 , (7/8) 3 (6/8) 3 for N x = 5, 6, 7 (N T = 2); and (8/9) 3 (7/9) 3 (6/9) 1 , (9/10) 3 (8/10) 3 (7/10) 2 , (10/11) 3 (9/11) 3 (8/11) 3 for N x = 8, 9, 10 (N T = 3). It is significant that there is a priority in the distribution of zero-energy states to the modes. Finally, we turn to the dependence of the decomposition on N y . The cases of N x = 8 and smaller N y can be compared in table 3. We considered that the configuration of the zero-energy states given by equation (55) holds true independent of N y , while the least squares method requires a larger M y to obtain more accurate resultseven more so for a Y-shaped GF with larger N x . The coefficient c j x was affected by N y : c j x for  D 1 j x was more dependent than c j x for < D 1 j x because the wavefunction extended from the zigzag edges in the rectangular parts to the triangular part. The coefficient can be roughly estimated as ( )