Semi-Empirical Pseudopotential Method for Graphene and Graphene Nanoribbons

We implemented a semi-empirical pseudopotential (SEP) method for calculating the band structures of graphene and graphene nanoribbons. The basis functions adopted are two-dimensional plane waves multiplied by several B-spline functions along the perpendicular direction. The SEP includes both local and non-local terms, which were parametrized to fit relevant quantities obtained from the first-principles calculations based on the density-functional theory (DFT). With only a handful of parameters, we were able to reproduce the full band structure of graphene obtained by DFT with a negligible difference. Our method is simple to use and much more efficient than the DFT calculation. We then applied this SEP method to calculate the band structures of graphene nanoribbons. By adding a simple correction term to the local pseudopotentials on the edges of the nanoribbon (which mimics the effect caused by edge creation), we again obtained band structures of the armchair nanoribbon fairly close to the results obtained by DFT. Our approach allows the simulation of optical and transport properties of realistic nanodevices made of graphene nanoribbons with very little computation effort.


Introduction
The electronic properties of low-dimensional materials are unique due to unprecedented properties that are unparalleled by those in bulk counterparts [1,2]. A scientific breakthrough occurred in 2004 in the isolation of monolayer graphene by mechanical exfoliation. After the discovery of graphene, the research on two-dimensional (2D) materials stimulated a great deal of interest due to their promising optical, physical, and chemical properties. These materials are layered, weakly coupled materials that can exist in a few or single-layer forms. Many 2D materials can be easily fabricated due to recently developed cutting-edge technology.
Ab initio formalism such as density-functional theory (DFT) [3] has been used extensively in calculating the electronic structures [4] and thermoelectric [5] and optical properties of solids [6]. In most calculations, three-dimensional (3D) plane waves are used as the basis functions. Plane wave (PW) basis is easy to implement, and the convergence of the calculation can be checked systematically. DFT methods are still affected by the band-gap problem, which requires numerically expensive GW calculation [7] or hybrid functionals to correct the band gap and effective mass near the Fermi level. The other drawback of this basis is the requirement of a large number of PWs needed when the unit cell is large. their semiconducting behavior [20], meaning they can exhibit a band gap in their electronic structure. The width of AGNRs determines the size of the band gap, with narrower ribbons having larger band gaps in general.
In electronics, AGNRs offer tunable bandgaps and unique edge states, making them potential candidates for high-performance transistors and logic devices [26]. AGNRs also exhibit attractive optical properties, with efficient light absorption and emission capabilities, making them suitable for optoelectronic applications such as photodetectors [27] and light-emitting diodes [28]. Furthermore, AGNRs have shown potential in topological physics [29][30][31][32], with the ability to host topological edge states and phenomena such as the quantum spin Hall effect. This opens up possibilities for the development of topological quantum devices and spintronics. Additionally, AGNRs can be tailored to exhibit magnetic properties [33,34] by introducing magnetic dopants or adatoms, enabling applications in spintronics, magnetic sensors, and information storage. AGNRs show promise in thermoelectric applications [35][36][37][38] due to their low thermal conductivity and tunable bandgap, enabling efficient energy conversion. The multifunctionality of armchair graphene nanoribbons makes them a highly versatile platform for advancing electronic, optical, topological, thermoelectric, and magnetic technologies, with ongoing research aimed at uncovering their full potential.
In this paper, we implement a more sophisticated semi-empirical pseudopotential method (SEPM) within the planar basis to study the atomic and electronic structure of 2D systems, using graphene and armchair graphene nanoribbons as examples. The applications of the currently developed SEP include optical and transport properties of graphene and related structures such as graphene nanoribbons and various junction devices.

B-Splines
To calculate the band structures of 2D materials in the semi-empirical pseudopotential (SEP) model, it is more efficient to adopt a planar basis. The basis consists of localized finite-element functions, i.e., B-splines [39] in the z coordinate multiplied by plane waves in the in-plane coordinates (x and y). B-splines of order κ consist of positive piecewise polynomials of z with degree κ − 1. These polynomials vanish everywhere outside the sub-inervals t i ≤ z < t i+k . The B-spline basis set with order κ defined by a knot sequence {t i } is generated by the following relationship: with B i1 (z) = 1 if t i ≤ z < t i+1 0 otherwise , where i = 1, 2, 3 . . . up to the number of knot sequence.
The first derivative of the B-spline of order κ is given by the following: The derivative of B-splines of order κ is a linear combination of B-splines of order κ − 1, which is also a simple polynomial and is continuous across the knot sequence. Using the polynomial expansion we obtain the Fourier transform of B iκ (z) as the following: where I i,j n (g) can obtain by the following recursion relationship: with I i,j Here, L c is the period length used along z, and I n (g) can be obtained by the recursion relationship. The B-spline functions used for the current calculation are shown in in Figure 1.

Kinetic and Overlap Matrix Elements within Planar Basis
Throughout this paper, we work in atomic units, where energy is in units of Ry, and distance is in units of Bohr. The Hamiltonian of the system in the ultrasoft pseudopotentia approach [16] is given by the following: Here, and denote the local pseudopotential and non-local pseudopotential Planar basis, used to expand , is defined as the following: where denotes the in-plane reciprocal lattice vector in 2D, and = ( , ) is the projection of r in the x-y plane.
is the in-plane wave vector. A is the surface area of the sample. The overlap matrix elements between two planar basis functions are given by the following: where Figure 1. B-spline consists of an exponential-type of a 17-point knot sequence of order κ = 4. Here, (•) denotes the knot points in the sequence.

Kinetic and Overlap Matrix Elements within Planar Basis
Throughout this paper, we work in atomic units, where energy is in units of Ry, and distance is in units of Bohr. The Hamiltonian of the system in the ultrasoft pseudopotential approach [16] is given by the following: Here, V loc andV nl denote the local pseudopotential and non-local pseudopotential. Planar basis, used to expand ∑ i , is defined as the following: where G denotes the in-plane reciprocal lattice vector in 2D, and ρ = (x, y) is the projection of r in the x-y plane. k is the in-plane wave vector. A is the surface area of the sample. The overlap matrix elements between two planar basis functions are given by the following: where is the overlap integral between two B-spline functions. The kinetic energy matrix elements are given by the following: where

Implementation of the Semi-Empirical Local Pseudopotential for Graphene
The local pseudopotential of the crystal is given by the following: where the first term describes the ionic local potential with in which R denotes a bulk lattice vector, and τ σ denotes the position of different atoms within the bulk unit cell. We note that V σ L (r) consists of a long-range term that decays like Z σ /r for large r. For charge-neutral systems, there is a counter long-range term in V H (r) due to the valence charges, and the sum of V ion (r) + V H (r) will be short-ranged.
For fitting empirical pseudopotentials of bulk materials, one adjusts the Fourier transform of the potential of V loc (r) (called "form factors" ∼ V loc (G)) at a small number of reciprocal lattice vectors until the band structure agrees with the experimental data or first-principal calculations [14]. To understand the nature of potential, a set of form factors for the first few shells at a large number of additional reciprocal lattice vectors is needed. A variety of algebraic forms have been used in the past for bulk materials [40][41][42]. We have not found an existing form that has sufficient flexibility to obtain the correct band structure for 2D materials. In our case, we try to mimic the full local potential by an analytic expression of the form.
where G denotes an in-plane 2D reciprocal lattice vector. The first term in Equation (14) denotes the main part of the local pseudopotential (denoted V 0 ), which is parametrized in terms of three spherical Gaussian functions, with C s and α s being fitted parameters. Namely, V 0 (r) = ∑ 3 s=1 C s e −α s r 2 . C s and α s are related to the spatial average of V loc (r) in the 2D plane as a function of z. We can rewrite V 0 (r) in terms of 2D plane waves as the following: Let ∼ V loc (z, G) denote the Fourier transform of V loc (r) in the 2D reciprocal space. We have for the τ σ = 0 term at G = 0. Here, N a denotes the number of atoms per unit cell, and A c is the area of the 2D unit cell. We determine the fitting parameters D s and α s by fitting the corresponding ∼ V loc (z, 0), determined by a DFT calculation within the same mixed basis as defined in (7). The implementation of the DFT package within this basis for graphene and AGNRs was reported in [12]. The DFT results obtained by using the package developed in [12] were checked against results obtained with the VASP package [43], and the calculated results for graphene and related nanoribbons obtained by using both packages are essentially the same. The exchange-correlation functional used is deduced from the Monte Carlo results calculated by Ceperley and Alder [44] and parametrized by Perdew and Zunger [45]. The fitted coefficients for D s can be directly converted to C s through the relationship C s = D s A c α s /(N a π). The best-fit parameters for graphene are given in the first row of Table 1, and the quality of the fit for ∼ V loc (z, 0) is shown in Figure 2a. Table 1. Fitting parameters for the main part of local potential V 0 (r).

Exponents
Coefficients fitted coefficients for can be directly converted to through the relationship = /( ). The best-fit parameters for graphene are given in the first row of Table 1, and the quality of the fit for ( , 0) is shown in Figure 2a.  The second term of Equation (14) denotes the difference , which is expressed in reciprocal space. Since = is already well fitted by ( , 0), we only have to consider the contribution. Due to point-group symmetry, the reciprocal lattice vector can be sorted into many shells with the magnitude of vectors being the same in each shell, and ∆ ( , ) is the same for all vectors in the same shell. We found that for vectors with magnitude > 4 a.u., ∆ ( , ) can be well fitted by a short-range correction function of the form: ( ) = ∑ • denotes the structure factor of graphene. The second term of Equation (14) which is expressed in reciprocal space. Since G = 0 is already well fitted by ∼ V loc (z, 0), we only have to consider the G = 0 contribution. Due to point-group symmetry, the reciprocal lattice vector G can be sorted into many shells with the magnitude of G vectors being the same in each shell, and ∆ ∼ V loc (z, G) is the same for all G vectors in the same shell. We found that for G vectors with magnitude G > 4 a.u., ∆ ∼ V loc (z, G) can be well fitted by a short-range correction function of the form: S 1 (G) = 1 A c ∑ σ e −iG·τ σ denotes the structure factor of graphene. Note that the chosen form of ∆ ∼ V S (z, G) goes to zero at G = 0, so the fitting does not affect ∼ V loc (z, 0). We can transform ∆ ∼ V S (z, G) back to real space analytically and obtain the following: where . The best-fit parameters for short-range correction functions defined above are listed in the first row of Table 2. Finally, for the three shells with 0 < G < 4 a.u., we need to fit the by another long-range correction function of the form The best-fit parameters for fitting f L (z) and ∼ D L (G) for the three shells are shown in the second row of Table 2. The Fourier transforms of the effective local pseudopotential, i.e., ∼ V loc (z, G) for G shells with magnitude G 1 , G 2 , G 3 = 1.56, 2.702, and 3.12 a.u., respectively, are shown in Figure 2b-d. Similarly, we can transform ∆ ∼ V L (z, G) back to real space analytically and obtain the following: with D L (ρ) = ∑ σ,G ∼ D L (G)e iG·(ρ−τ σ ) . Then, we obtain the following relationship: where γ = S or L labels the short-range or long-range term. Table 2. Fitting parameters for the short-range and long-range shape function f S (z) and f L (z) for the correction terms to ∼ V loc (z, G) used in this work.

Type
Exponent Combining all above, the total local potential is given by the following: For convenience, D γ (ρ) can be expressed as a polynomial multiplied by a Gaussian function. where b γ m can be determined by the following relationship: which leads to the following relationship between b γ m and the set of parameters P To check if the real space form of the local pseudopotential V loc (r) obtained by the current SEPs as given by Equation (21) can truthfully represent the DFT results, we compared the results for V loc (r) obtained by DFT and by the current method in Figure 3. In Figure 3a, we show the x-dependence of V loc (r) at z = −0.5∆z and y = a/ 2 √ 3 (a line passing through the center of a C atom). Here, ∆z = W/104 is the grid size for the fast Fourier transform used in the DFT calculation, and W is the width of the domain along the z-axis used to define the B-spline basis.
the center ( 0). We also noticed that the x-and y-dependence almost coincide for | | < 0.2 , indicating the potential has nearly cylindrical symmetry near the center of the hexagon cell. We then took the difference of the local potential obtained by DFT and by SEP and plotted the difference function ∆ ( , 0) in Figure 3c. It turns out that such a minor correction will still be important for the states derived from the -bonds, but it has a negligible effect on states derived from the -bonds since the -orbitals have negligible amplitude at the center of the hexagon cell. This minor correction is caused by the fact that all terms included in ( ) in Equation (21) are centered at atomic sites, while the effect of charge redistribution due to valence charges in the solid cannot be fully absorbed by the atom-centered terms. Thus, we should also consider contributions described by functions localized at the center of hexagon cells in the graphene lattice. We shall call such a term the "bond-charge contribution".
To determine this contribution, we fit the difference of the local potential obtained by DFT and by SEP, i.e., ∆ ( , 0), in Figure 3c by the following analytic function: The best-fit parameters { ; = 0, … ,4} and are listed in Table 3. We approximate the net local potential near the center of the hexagon by a separable form: where ( ) = ( , )/ ( , 0) describes the z-dependence of the bond-charge contribution. Here, ( ) describes the variation of ( , )/ ( , 0) along the z-axis at = , i.e., the center of the hexagon cell. By taking the difference between ( ) obtained by DFT and by SEP, we obtain the following: where we assume that the z-dependent functions obtained by DFT and SEP are the same. , which passes through a row of C atoms in the graphene sheet. (b) Along two perpendicular lines (along the x-axis and y-axis) both passing through the origin, which is the center of the hexagon cell. (c) The difference of V loc (ρ, z) obtained by DFT and SEP plotted along the x-axis (green) or y-axis (red). The blue curve indicates the average of the red and green curves. (d) The best-fit results of the average ∆V b (ρ, z).
It is seen that the result obtained by SEP matches the corresponding DFT results very well. In Figure 3b, we show both the x-and y-dependence of V loc (r) at z ≈ 0 along lines in x-and y-directions, with both passing through the center of the hexagon cell. The agreement with the DFT results is still very good, except there is a small deviation near the center ( ρ ≈ 0). We also noticed that the x-and y-dependence almost coincide for |ρ| < 0.2a, indicating the potential has nearly cylindrical symmetry near the center of the hexagon cell. We then took the difference of the local potential obtained by DFT and by SEP and plotted the difference function ∆V b (ρ, 0) in Figure 3c. It turns out that such a minor correction will still be important for the states derived from the σ-bonds, but it has a negligible effect on states derived from the π-bonds since the π-orbitals have negligible amplitude at the center of the hexagon cell. This minor correction is caused by the fact that all terms included in V loc (r) in Equation (21) are centered at atomic sites, while the effect of charge redistribution due to valence charges in the solid cannot be fully absorbed by the atom-centered terms. Thus, we should also consider contributions described by functions localized at the center of hexagon cells in the graphene lattice. We shall call such a term the "bond-charge contribution".
To determine this contribution, we fit the difference of the local potential obtained by DFT and by SEP, i.e., ∆V b (ρ, 0), in Figure 3c by the following analytic function: The best-fit parameters { C b n ; n = 0, . . . , 4 and α b are listed in Table 3. We approximate the net local potential near the center of the hexagon by a separable form: i.e., the center of the hexagon cell. By taking the difference between V loc (r) obtained by DFT and by SEP, we obtain the following: where we assume that the z-dependent functions obtained by DFT and SEP are the same. Table 3. Fitting parameters for the bond-charge contribution in pseudopotential localized at the center of the hexagon cell as described by Equations (25)- (27). The net local potential of graphene, V loc (ρ, z), evaluated at ρ = 0 and obtained by DFT is shown in Figure 4. For convenience, we fit its normalized z-dependent function f b (z) = V loc (ρ, z)/V loc (0, z) by the following expression:

Fitting of the Non-Local Pseudopotential for Graphene
The semi-empirical pseudopotential contains both local and non-local term non-local pseudopotential contains angular momentum and energy-dependen [15]. Incorporating non-local in addition to local terms can provide a more accurate range for the valence band edge and best matches to the experimental data. The n potential is given by a separable form [16]: where ( ) denotes the projector functions for an atom at position . For eac the best-fit β functions take the following form:

Fitting of the Non-Local Pseudopotential for Graphene
The semi-empirical pseudopotential contains both local and non-local terms. The non-local pseudopotential contains angular momentum and energy-dependent terms [15]. Incorporating non-local in addition to local terms can provide a more accurate energy range for the valence band edge and best matches to the experimental data. The nonlocal potential is given by a separable form [16]: where β nσ lm (r) denotes the projector functions for an atom at position τ σ . For each atom, the best-fit β functions take the following form: for l, m = 00, 10, 11, 20, 21, 22. We fit the beta function for 2S and 2P orbitals for C atom from the Vanderbilt ultrasoft pseudopotential (USPP) [16]. The fitting potential used has the form for r < R s , where R s is the cut-off radius.
For the second 2P orbital (2P 2 ), we break the β function into two segments with a dividing radius at R s = 0.9228 (indicated by a vertical dashed line in Figure 5d). We fit the first segment (seg 1) for r < 0.9228 by expression (32). For the second segment (seg 2), we fit A n l (r) with the following expression: where r s = r − 0.9228. The bet-fit parameter obtained is given in Table 4 The fitting result for each orbital is shown in Figure 5. The quality of fitting to the input of USPP [16] is excellent.

Matrix Elements of Local and Nonlocal Pseudopotential
The local pseudopotential ( ) is given by Equation (12). It consists of the atomiclike term ( ), which is spherical in 3D space, and the correction term, which has a shortrange part and a long-range part. We define the following:  Table 4. Fitting parameters (a.u.) for β functions used in the non-local pseudopotential of C atom. The fitting result for each orbital is shown in Figure 5. The quality of fitting to the input of USPP [16] is excellent.

Matrix Elements of Local and Nonlocal Pseudopotential
The local pseudopotential V loc (r) is given by Equation (12). It consists of the atomiclike term V 0 (r), which is spherical in 3D space, and the correction term, which has a short-range part and a long-range part. We define the following: B i (z) is the B-spline basis functions with the subscript κ dropped for brevity. Within the mixed basis, the matrix elements for the atomic-like term can be written as follows: within the B-spline basis. Here, ∆G = G − G and C s denote the fitting parameters for a given C atom in the unit cell, as given in Table 1.
denotes the structure factor for graphene. Due to the inversion symmetry of graphene, we choose the origin to be the center between two C atoms in the unit cell, and τ 1 is the position of one C atom in the unit cell. Thus, S 1 ( ∆G) is real. Similarly, we obtain the matrix elements for the remaining terms for the local potential as follows: where denotes the in-plane Fourier transform of the bond-charge contribution, ∆V b (ρ, 0), given in Equation (25).

Nonlocal Corrections in Overlap and Potential
The matrix elements for nonlocal potential read as follows [12]: Here, lm denotes the projection of the β functions (with the center shifted to the origin) into our basis. A detailed description of the evaluation of P in lm (K) is given in Appendix A. For calculating the band structure, we also need to evaluate the correction to the overlap matrix elements. The correction to the overlap matrix can be written as follows [12]: This expression is identical to (39) except that the energy parameters E nn lm are replaced by the overlap parameters q nn lm . Here, we adopt the same parameters for E nn l and q nn l as in the ab initio input data [12,16] to calculate the matrix elements of non-local pseudopotentials. The input parameter used for E nn l and q nn l are shown in Table 5. The formulas derived above can be extended to other two-dimensional materials beyond graphene, such as transition-metal dichalcogenides. Due to the inversion symmetry in graphene, all matrix elements become real. Furthermore, there is a mirror symmetry for the z-axis. By adopting the B-spine basis along the z-axis and taking symmetric and antisymmetric combinations of these basis functions, we can decouple the eigenstates of the Hamiltonian within the symmetric and antisymmetric basis sets for any wave vector in the Nanomaterials 2023, 13, 2066 13 of 30 2D plane. This significantly improves the speed of computation for the band structures. We found that using a direct solver for diagonalization becomes even faster than the iterative solver based on the conjugate-gradient (CG) approach [46] for this system with both inversion and mirror symmetry. The cut-off used for plane waves used in wavefunctions is 30 a.u. and that in pseudopotential is 180 a.u.

Band Structure of Graphene
The all-electron (AE) Bloch states of graphene are written as follows: where ϕ n,k (ρ, z) denotes the pseudo-wavefunction, ∼ B i (z) denotes the orthogonalized Bspline function, andÔ is the overlap operator defined in (A5). The all-electron Kohn-Sham equation [47] reads as follows: whereĤ AE is the all-electron Hamiltonian operator, and E ν (k) denotes the energy of the ν-th band at wavevector k. Substituting Equation (41) into Equation (42) gives rise to the following generalized eigenvalue problem: whereĤ is the Hamiltonian operator given in Equation (6). We apply the current SEP to calculate the band structure of graphene. The realspace structure of graphene can be described by a 2D unit cell with primitive lat-  Figure 6 shows the geometry of graphene in real space and reciprocal space.
In Figure 7, the band structure of graphene calculated by SEP with best-fitted parameters deduced in Section 2 is shown and compared with that obtained from a self-consistent calculation based on DFT [12]. The overall band structures obtained by the present method are in close agreement with the ab initio calculation. Since we solve the Bloch states with even and odd parities (with respect to the mirror symmetry about the z-axis) separately, the symmetry characteristics of the bands can be easily distinguished. Here, the even (odd) parity states are presented in green (red) color. It is noted that the red curves (with odd parity) are related to the p z orbitals of the carbon atoms, and they form π-bonded bands, whereas the green curves (with even parity) are related to the s, p x , andp y orbitals.
As shown in Figure 7, the agreement between the SEP and DFT results obtained by using the method described in [12] for the lowest 20 bands is excellent. The advantage of the SEP is that it is easy to use and requires no self-consistent calculation to establish the charge density. The time needed to calculate the band structures via a direct solver on a laptop PC is less than 6 s, which is only a fraction (~1/5) of the time needed for the CG calculation used in the DFT package (after the self-consistent density is established). The saving is mainly due to the use of the inversion and mirror symmetry properties of graphene.
We apply the current SEP to calculate the band structure of graphene. The real-spa structure of graphene can be described by a 2D unit cell with primitive lattice vectors (1, − √3) a/2 and = (1, √3) a/2. The positions of carbon atoms are (1, −1/√3)a/2 a (−1,1/√3)a/2. The basis vectors in reciprocal space are = √ √3, −1 , = √ (√3, Figure 6 shows the geometry of graphene in real space and reciprocal space. In Figure 7, the band structure of graphene calculated by SEP with best-fitted parameters deduced in Section 2 is shown and compared with that obtained from a selfconsistent calculation based on DFT [12]. The overall band structures obtained by the present method are in close agreement with the ab initio calculation. Since we solve the Bloch states with even and odd parities (with respect to the mirror symmetry about the zaxis) separately, the symmetry characteristics of the bands can be easily distinguished. Here, the even (odd) parity states are presented in green (red) color. It is noted that the red curves (with odd parity) are related to the orbitals of the carbon atoms, and they form -bonded bands, whereas the green curves (with even parity) are related to the , , and orbitals.
As shown in Figure 7, the agreement between the SEP and DFT results obtained by using the method described in [12] for the lowest 20 bands is excellent. The advantage of the SEP is that it is easy to use and requires no self-consistent calculation to establish the charge density. The time needed to calculate the band structures via a direct solver on a laptop PC is less than 6 s, which is only a fraction (~1/5) of the time needed for the CG calculation used in the DFT package (after the self-consistent density is established). The saving is mainly due to the use of the inversion and mirror symmetry properties of graphene.

Band Structure of Armchair Graphene Nanoribbon
We consider a graphene nanoribbon with armchair edges with a width of , where is the lattice constant of graphene. The supercell of the armchair graphene nanoribbon (AGNRs) contains 4 + 2 atoms, as illustrated in Figure 8a. To adopt the 2D plane-wave basis, we introduce a vacuum region with the width of ( − ) inside the supercell such that the dimension of the rectangular supercell is √3 × . The reciprocal lattice vectors of a superlattice consisting of a 1D periodic of AGNRs along the -axis can be written as = + , where = and = are the basis vectors in Figure 7. Band structure of graphene obtained by the present SEP with best-fitted parameters (solid curves). For comparison, the band structure obtained by self-consistent calculation based on DFT with Vanderbilt USPP is also included (dotted curves).

Band Structure of Armchair Graphene Nanoribbon
We consider a graphene nanoribbon with armchair edges with a width of Na, where a is the lattice constant of graphene. The supercell of the armchair graphene nanoribbon (AGNRs) contains 4N + 2 atoms, as illustrated in Figure 8a. To adopt the 2D plane-wave basis, we introduce a vacuum region with the width of (M − N)a inside the supercell such that the dimension of the rectangular supercell is √ 3a × Ma. The reciprocal lattice vectors of a superlattice consisting of a 1D periodic of AGNRs along the x-axis can be written as g n 1 n 2 = n 1 Max and ∼ b 2 = π aŷ are the basis vectors in reciprocal space, and n 1 and n 2 are arbitrary integers. Let g i (i = 1, . . . , 2M) denote those non-equivalent g vectors falling within the first Brillouin zone of graphene. (See Figure 8b.) Then, all g vectors of the AGNRs superlattice can be expressed as g i + G, where G denotes the reciprocal lattice vectors of graphene. edges. To describe such a change, it is computationally more efficient to start with the zone-folded band structures of graphene that can be calculated by using the graphene SEP code developed above at wave vectors of + ( = 1, … ,2 ), where is within the first mini zone of the AGNR superlattice. A similar approach was applied to the Si 7 × 7 surface, and it was demonstrated that such a method works very well for large superstructures [48]. If the vacuum region introduced in the superlattice is thick enough, the coupling between adjacent AGNRs can be negligible, and of interest will be along the long axis of the AGNR (taken to be parallel to here). The low-lying pseudo-Bloch states of graphene (with band labeled by ν) obtained at + (denoted as , ) can then be used as a set of contracted basis functions for calculating the band structures of the AGNR. This set of contracted bases is a nearly complete basis for AGNR if we use a large number of graphene bands until the convergence is established for the calculated results.
denotes the sum of local potentials of C atoms defined in Equation (12) for atomic sites (indexed by ) in the nanoribbon region [with | | ≤ /2] inside the supercell used for AGNRs.
denotes the nonlocal pseudopotential of the nanoribbon. The AGNRs Hamiltonian matrix within the contracted basis is then given by the following: The Bloch states of the AGNRs superlattice can be written as linear combinations of basis functions that are products of 2D plane waves (labeled by wave vectors k + g j + G ) and B-splines in z (i.e., the mixed-basis introduced in [11]). For wide AGNRs, this basis set can be quite large. The band structure of the AGNRs is expected to be quite close to the zone-folded band structure of graphene, with deviation mainly coming from the quantum confinement effect from the vacuum region and the effect due to the creation of edges. To describe such a change, it is computationally more efficient to start with the zone-folded band structures of graphene that can be calculated by using the graphene SEP code developed above at wave vectors of k + g j (j = 1, . . . , 2M), where k is within the first mini zone of the AGNR superlattice. A similar approach was applied to the Si 7 × 7 surface, and it was demonstrated that such a method works very well for large superstructures [48]. If the vacuum region introduced in the superlattice is thick enough, the coupling between adjacent AGNRs can be negligible, and k of interest will be along the long axis of the AGNR (taken to be parallel toŷ here). The low-lying pseudo-Bloch states of graphene (with band labeled by ν) obtained at k + g j (denoted as ϕ ν,k+g j ) can then be used as a set of contracted basis functions for calculating the band structures of the AGNR. This set of contracted bases is a nearly complete basis for AGNR if we use a large number of graphene bands until the convergence is established for the calculated results.
Thus, we can write the Hamiltonian of the AGNR as H = −∇ 2 + U 1 +Û 2 , where denotes the sum of local potentials of C atoms defined in Equation (12) for atomic sites (indexed by σ) in the nanoribbon region A NR [with |x| ≤ W/2] inside the supercell used for AGNRs.Û 2 denotes the nonlocal pseudopotential of the nanoribbon. The AGNRs Hamiltonian matrix within the contracted basis is then given by the following: where the first term describes the kinetic energy, and the second term consists of the semi-empirical local ( U 1 ) and nonlocal (Û 2 pseudopotentials for AGNR. The detailed expressions for these matrix elements are given in Appendix B.
In addition to the z-mirror symmetry, the AGNR also has x-mirror symmetry for any wave vector, i.e., k, along the y-axis. Thus, we can define symmetrized basis states with respect to the x-mirror as follows: where g j = −g jx , g jy , ϕ s ν,k+g j (x, y, z) = ϕ s ν,k+g j (−x, y, z) (corresponding to flipping the sign of G x in the plane-wave basis). f j = 1/ √ 2 for g jx = 0, and f j = 1 otherwise. The symbol "±" denotes even/odd parity states with respect to the x-mirror, while the superscript s = e, o is the label for the even/odd parity states with respect to the z-mirror. The Hamiltonian matrix elements between symmetrized states of the same parity read as follows: The overlap matrix is given by the same expression above, withĤ GNR replaced byÔ GNR .
Here, we divide the basis set into four different subsets according to different symmetry types labeled by { π + , π − , σ + , σ − }. π + and π − states represent π-bonded (p z -like) states with even and odd parity, respectively, under the x-mirror operation. σ + and σ − states represent σ-bonded (p x , p y , or s-like) states with even and odd parity, respectively, under the x-mirror operation. Eigenstates ofĤ GNR with different symmetry types are decoupled based on group theory. Thus, we can block-diagonalize the Hamiltonian matrix into four diagonal sub-blocks, and each diagonal sub-block can be diagonalized separately. By taking advantage of these symmetry properties, the electronic structures of graphene nanoribbons can be solved very efficiently with the present SEP via a direct diagonalization method. Furthermore, the symmetry characteristics of different bands can be easily identified in the band structure by using different colors. This can help sort out the complex characteristics in AGNR electronic states and improve the understanding of the electronic properties of AGNR.
The graphene bulk basis adopted above can describe the band folding due to the quantum confinement effect in the nanoribbon well, but it requires a large number of basis states to fully capture the effect of localized edge states. The AGNRs band structures obtained by solving Equation (43) are shown in Figure 9a,b to compare with the DFT results. The number of graphene bands included in the calculation is 10 for the π + (π − ) states and 60 for the σ + (σ − ) states. Adding more bands produces no significant effect. For the odd states, the change of energy on the AGNR band structure is less than 0.005 eV within the energy window of interest. For the σ-bonded (even in z) states, most minibands remain nearly unchanged when the number of graphene bands included varies from 14 to 60, while the two pairs of minibands closest to the band-gap region change significantly and still do not reach the desired convergence level even with 60 graphene bands included in the basis. These edge-induced states correspond to the bonding and antibonding states of the dimers on two edges. Since there will be an edge-induced correction in the local potential, it is unnecessary to spend a significant effort to obtain fully convergent results for the preliminary investigation in this section. sufficient, we also performed the DFT calculation of the AGNR, in which the outermo four rows of atoms are allowed to relax, and we found that the major displacement occu on the edge atoms, while the other inner rows only relax slightly, and the resulting ban structure is quite similar to the one with relaxation only on the edge atoms. F comparison, we shifted the DFT band structure rigidly in energy, so the valence-ban maximum is aligned with our SEP calculation. (a) Our SEP results without t relaxation of the edge atoms. (b) Our SEP results with edge relaxation but without modifying t pseudopotentials on the edge atoms. (c) DFT results were obtained by using the method describ in [12] for the AGNRs with relaxed atomic positions for the edge atoms. (d) Our SEP results inclu the modification of pseudopotentials on the edge atoms. In SEP results (a,b,d), the bands in bl ( ) and red (with symmetry) are derived from -bonded states (odd wi respect to the z-mirror), while the bands in green are derived from the -bonded states (even wi respect to the z-mirror). Here, we do not distinguish the from states since the importa states are edge states, and they are essentially degenerate. Our SEP results with edge relaxation but without modifying the pseudopotentials on the edge atoms. (c) DFT results were obtained by using the method described in [12] for the AGNRs with relaxed atomic positions for the edge atoms. (d) Our SEP results include the modification of pseudopotentials on the edge atoms. In SEP results (a,b,d), the bands in blue (with π − symmetry ) and red (with π + symmetry) are derived from π-bonded states (odd with respect to the z-mirror), while the bands in green are derived from the σ-bonded states (even with respect to the z-mirror). Here, we do not distinguish the σ + from σ − states since the important states are edge states, and they are essentially degenerate.
We also performed DFT calculations of the same AGNR within the basis constructed by 2D plane waves multiplied by B-spline functions of z, as described in [12]. The edge atoms are allowed to relax to minimize the total energy of the AGNR. In the relaxed geometry, the displacement of the edge atom in the left and upper corner of the AGNR supercell (the green box) as shown in Figure 8a is described by ∆τ = (∆x, ∆y), with ∆x = 0.2107 a.u. and ∆y = 0.1788 a.u. The displacement of three other edge atoms can be deduced by applying the x-mirror and y-mirror symmetry. The resulting band structures are shown in Figure 9c. To make sure that the relaxation of only the edge atoms is sufficient, we also performed the DFT calculation of the AGNR, in which the outermost four rows of atoms are allowed to relax, and we found that the major displacement occurs on the edge atoms, while the other inner rows only relax slightly, and the resulting band structure is quite similar to the one with relaxation only on the edge atoms. For comparison, we shifted the DFT band structure rigidly in energy, so the valence-band maximum is aligned with our SEP calculation.
We found that, as shown in Figure 9a, the band structure of the unrelaxed AGNR looks qualitatively similar to the corresponding DFT results (with edge relaxation) in Figure 9c. However, the band gap obtained is about 0.45 eV, which is substantially smaller than the DFT result of~0.6 eV. Furthermore, we found a pair of mid-gap bands (in green) that are related to the σ-bonding states of the dimer atoms on the edge. There is also a pair of σ-antibonding states of the edge dimers with energy near 0 eV at the Γ point (not shown). These edge-dimer-related minibands are nearly doubly degenerate, corresponding to dimer states located at the left and right edges of the AGNR. Since the width of the ribbon is much larger than the spread of the wavefunctions localized on edges on opposite sides, these states form a closely spaced pair of minibands. In the DFT band structure, the corresponding σ-bonding states are lowered by~2.5 eV at the Γ point, and there are no mid-gap states left. Since the DFT calculation did not separate the even and odd states, some of them even states look artificially entangled with the odd states.
The SEP band structures of the AGNR with the same relaxed atomic positions as determined by DFT but with the pseudopotential (PP) given by Equation (29) are shown in Figure 9b. With the relaxation, the band gap becomes enlarged to~0.5 eV, and the overall band structure agrees better with the DFT results. The pair of σ-bonding states of edge-dimers (in green) is also lowered but only by~0.5 eV at Γ point, and there are no states left in the mid gap. We note that for device applications, the π-bonded states (near the band gap region) play a dominant role. The band structure for these π-bonded states is not significantly changed due to the relaxation of atomic positions at the edges. However, there are still noticeable differences between the band structures obtained by SEP and by DFT, which are mainly caused by the deviation of local potential near the edge atoms as the graphene is cleaved to form a nanoribbon. We discuss the consequence of this effect and its remedy below. Figure 10 shows the contour plot of the local potential of the 9 × 2 AGNR with edge relaxation. The~5% shift of atomic positions near the edges is noticeable in this plot. To examine the deviation of local potential (V loc ) calculated by SEP from the DFT results, we plot the line cuts of V loc (x, y, z) as functions of x at the four different values of y ( y = −L/2, −3L/17, 0, and 6L/17) and z ≈ 0 in Figure 11. These lines are indicated by black lines in Figure 10. Here, L = √ 3a is the length of the AGNR supercell along the y-axis.

Modification of Pseudopotential for Edge Atoms of Armchair Graphene Nanoribbon
As shown in Figure 11, we found that the local potentials obtained by SEP agree very well with the DFT results in the interior region of the AGNR, with noticeable deviation from the DFT results only near the two edges (with |x|> 4.5a ). It implies that the SEP potential centered at an atom follows the position of the relaxed atom very well, while the cleavage of graphene to form an AGNR causes some charge redistribution that can lead to the change of the local potential near the edges. Figure 10 shows the contour plot of the local potential of the 9 × 2 AGNR with edge relaxation. The ~5% shift of atomic positions near the edges is noticeable in this plot. To examine the deviation of local potential ( ) calculated by SEP from the DFT results, we plot the line cuts of ( , , ) as functions of x at the four different values of ( = − 2 ⁄ , − 3 17 ⁄ , 0, and 6 /17) and 0 in Figure 11. These lines are indicated by black lines in Figure 10. Here, = √3 is the length of the AGNR supercell along the y-axis. Figure 10. Contour plot of net local potential in one supercell of the AGNR obtained by DFT [11].

Modification of Pseudopotential for Edge Atoms of Armchair Graphene Nanoribbon
Here, is the lattice constant of graphene, and = √3 is the the length of supercell along theaxis.
As shown in Figure 11, we found that the local potentials obtained by SEP agree very well with the DFT results in the interior region of the AGNR, with noticeable deviation from the DFT results only near the two edges (with | | > 4.5 ). It implies that the SEP potential centered at an atom follows the position of the relaxed atom very well, while the cleavage of graphene to form an AGNR causes some charge redistribution that can lead to the change of the local potential near the edges.  Figure 10), (c) = 0 (along the middle black line going through the centers of nine bonds in Figure 10), and (d) = 6 /17 (along the black line going through the centers of ten atoms in Figure  10).
To investigate the effect of local potential induced by the edge creation, we took the difference between the DFT (red curve) and SEP results (green curve) in Figure 11 for between 3 and 6 at = − /2, −3 /17, 0, and 6 /17. The results for the differences at the selected values of are shown as dotted lines in Figure 12. These difference curves (for one of the two edges) can be reasonably fitted by the following functional form: denotes the -coordinate of the relaxation edge atom, including the relaxation ∆ . The AGNR potential also has a y-mirror symmetry, so we duplicate the terms at and and add the same potentials at − and − . We note that the potentials at = − 2 ⁄ (equivalent to = 2 ⁄ due to periodicity) and = 0 will be mapped to themselves by the y-mirror. Here, = 5 denotes the bond-charge location at = − /2 , and = 4.65 denotes the common bond-charge location near the edges for other values. The Figure 11. The net local pseudopotential in one supercell of the AGNRs obtained by the current SEP (green line) and by DFT (red line) along the lines with z ≈ 0 (near the AGNRs plane) and some selected values of y. (a) y = −L/2 (along the bottom black line going through the centers of ten bonds in Figure 10), (b) y = −3L/17 (along the black line going through the centers of nine atoms in Figure 10), (c) y = 0 (along the middle black line going through the centers of nine bonds in Figure 10), and (d) y = 6L/17 (along the black line going through the centers of ten atoms in Figure 10).
To investigate the effect of local potential induced by the edge creation, we took the difference between the DFT (red curve) and SEP results (green curve) in Figure 11 for x between 3a and 6a at y = −L/2, −3L/17, 0, and 6L/17. The results for the differences at the selected values of y are shown as dotted lines in Figure 12. These difference curves (for one of the two edges) can be reasonably fitted by the following functional form: where ρ a = (x a , −L/2), ρ b = (x b , −3L/17), ρ c = (x b , 0), and ρ d = (x b , y d ) denote the four locations of effective bond charges near the AGNR right edge. y d = y e + ∆y ≈ 6L/17 denotes the y-coordinate of the relaxation edge atom, including the relaxation ∆y. The AGNR potential also has a y-mirror symmetry, so we duplicate the terms at y b and y d and add the same potentials at −y b and −y d . We note that the potentials at y a = −L/2 (equivalent to y a = L/2 due to periodicity) and y c = 0 will be mapped to themselves by the y-mirror. Here, x a = 5a denotes the bond-charge location at y = −L/2, and x b = 4.65a denotes the common bond-charge location near the edges for other y values. The positions of these bond charges are related to the minimum of the difference in the local potentials shown in Figure 12. S a , S b , S c , and S d , describing the variation of potential strength along the y-axis. Here, we choose α a = 0.4 and α b = 1.1, which approximately describe the width of the potential wells shown in Figure 12. The best-fit results for ∆V e (ρ, 0) are shown as green curves in Figure 12, with best-fit values of S a , S b , S c , and S d listed in Table 6. We approximate the net three-dimensional local potential, including the bond-charge redistribution near the edges by a separable form: where ( ) = ( , )/ ( , 0) describes the z-dependence of the bond-charge contribution near the edges. ( ) can be extracted from the net local potential of AGNR, ( , 0, ), obtained by DFT with = and = 0. The result is shown in Figure 13, and it can be well fitted by the sum of two Gaussian functions. The normalized zdependence of the local potential is given by the following: where = 1 − due to the normalization requirement similar to Equation (28). We obtain = 0.979 , = 1.0926, and = 0.1026 . Thus, the correction to local pseudopotential due to the bond-charge redistribution near the edges can be approximately given by the following: where ∆ ( , 0) is given in Equation (48) and ( ) in Equation (50). Finally, we add the edge-induced correction in the local pseudopotential, ∆ ( , ), to the Hamiltonian of the AGNR in our semi-empirical pseudopotential model. The matrix elements of ∆ ( , ) within the contracted basis functions derived from the selected graphene eigenstates at special reciprocal lattice vectors for AGNR ( ) are given by the following:  We approximate the net three-dimensional local potential, including the bond-charge redistribution near the edges by a separable form: where f e (z) = V loc (0, z)/V loc (0, 0) describes the z-dependence of the bond-charge contribution near the edges. f e (z) can be extracted from the net local potential of AGNR, V loc (x b , 0, z), obtained by DFT with x = x b and y = 0. The result is shown in Figure 13, and it can be well fitted by the sum of two Gaussian functions. The normalized z-dependence of the local potential is given by the following: where C e 2 = 1 − C e 1 due to the normalization requirement similar to Equation (28). We obtain C e 1 = 0.979, α e1 = 1.0926, and α e2 = 0.1026. Thus, the correction to local pseudopotential due to the bond-charge redistribution near the edges can be approximately given by the following: where ∆V e (ρ, 0) is given in Equation (48) and f e (z) in Equation (50). Here, ∆ ( , 0) is the edge-induced correction in pseudopotential for the right edge, as given in Equation (48), and c.c. denotes the term contributed from the left edge, which is the complex conjugate of the previous term due to the inversion symmetry in the 2D plane of the AGNR. With the use of the analytic fitting functions introduced above, ( ) can be written as follows: After adding the correction term ∆ ( , ) , the calculated band structure for the AGNR is shown in Figure 9d. To improve the convergence for the edge-dimer-related states, we included 90 graphene bands in our basis for the calculation of even states, while only 10 graphene bands are needed to achieve the desired convergence for odd states. We found that by adding more graphene bands to the calculation, the results remain nearly the same. Comparing with Figure 9c, where we did not include the edge-induced correction in pseudopotential, ∆ ( , ), we found that the main effect of ∆ ( , ) is to give a significant improvement of the miniband structures for the -bonded states (in red or blue) of the AGNR, which agree very well with the DFT results as shown in Figure 9c, and the band gap obtained by SEP (~0.7 eV) also matches the DFT result well. The other important feature is the lowering of the -bonding states of edge dimers at the zone center from −4.3 eV in Figure 9b to −6.3 eV in Figure 9d, while the -antibonding states of edge dimers are also lowered by about 2 eV, and they appear in Figure 9d at −2.3 eV at the Γ point. The energy spacing between these edge-dimer bonding and anti-bonding states is about 4 eV at k = 0 (the zone center), and it reduces to about 3.3 eV at = / (the zone boundary). This feature also agrees reasonably well with the DFT result shown in Figure  9b. The energy position of these edge-dimer bonding states (near −6.75 eV at the zone boundary) is close to the DFT result (−6.8 eV), while the edge-dimer antibonding states sit around −3.46 eV, which is higher than the corresponding DFT results by about 0.3 eV. This discrepancy is likely due to other corrections in the pseudopotential that are not yet included by the current SEP. On the other hand, the band structure of other even states Finally, we add the edge-induced correction in the local pseudopotential, ∆V e (ρ, z), to the Hamiltonian of the AGNR in our semi-empirical pseudopotential model. The matrix elements of ∆V e (ρ, z) within the contracted basis functions derived from the selected graphene eigenstates at special reciprocal lattice vectors for AGNR (g j ) are given by the following: where I(α; i, , i ) is defined in Equation (34), ∆g jj + ∆G = g j − g j + G − G, and Here, ∆V e (ρ, 0) is the edge-induced correction in pseudopotential for the right edge, as given in Equation (48), and c.c. denotes the term contributed from the left edge, which is the complex conjugate of the previous term due to the inversion symmetry in the 2D plane of the AGNR. With the use of the analytic fitting functions introduced above, ∼ v e (q) can be written as follows: with I 0 (α, q) = π α e −q 2 /4α . Due to the x-mirror symmetry, there are four corresponding bond charges at the left edge.
After adding the correction term ∆V e (ρ, z), the calculated band structure for the AGNR is shown in Figure 9d. To improve the convergence for the edge-dimer-related states, we included 90 graphene bands in our basis for the calculation of even states, while only 10 graphene bands are needed to achieve the desired convergence for odd states. We found that by adding more graphene bands to the calculation, the results remain nearly the same. Comparing with Figure 9c, where we did not include the edge-induced correction in pseudopotential, ∆V e (ρ, z), we found that the main effect of ∆V e (ρ, z) is to give a significant improvement of the miniband structures for the π-bonded states (in red or blue) of the AGNR, which agree very well with the DFT results as shown in Figure 9c, and the band gap obtained by SEP (~0.7 eV) also matches the DFT result well. The other important feature is the lowering of the σ-bonding states of edge dimers at the zone center from −4.3 eV in Figure 9b to −6.3 eV in Figure 9d, while the σ-antibonding states of edge dimers are also lowered by about 2 eV, and they appear in Figure 9d at −2.3 eV at the Γ point. The energy spacing between these edge-dimer bonding and anti-bonding states is about 4 eV at k = 0 (the zone center), and it reduces to about 3.3 eV at k = π/L (the zone boundary). This feature also agrees reasonably well with the DFT result shown in Figure 9b. The energy position of these edge-dimer bonding states (near −6.75 eV at the zone boundary) is close to the DFT result (−6.8 eV), while the edge-dimer antibonding states sit around −3.46 eV, which is higher than the corresponding DFT results by about 0.3 eV. This discrepancy is likely due to other corrections in the pseudopotential that are not yet included by the current SEP. On the other hand, the band structure of other even states (in green) not related to the edge-dimer states are in very good agreement with the DFT results. Since the edge-dimer antibonding states have a minimum at the zone boundary, and their energies are higher than the conduction band minimum at the zone center, they may not play a significant role in certain device applications. In case they are needed, a "scissor operation" can be applied to shift these minibands rigidly to match the DFT results, and phenomena such as the Gunn effect and negative differential resistance [47,48] can be simulated with this simple approach, and it remains a reasonable approximation to use the wavefunctions obtained by SEP with the current set of basis to calculate the carrier-scattering process in the nano Gunn-diode application by using AGNRs. Since these edge-dimer states are very sensitive to the chemical modification of the AGNR edges, our model can be used to simulate the effect of edge modification on the I-V characteristics by imposing a model potential on the edge atoms.

Comparison with Experiments
Many experimental studies on graphene and AGNR-related devices have been reported in the last decade [49][50][51][52][53][54][55][56][57][58][59][60][61], which validate and complement the theoretical predictions, providing critical empirical evidence. Techniques such as scanning tunneling microscopy (STM) [51][52][53][54][55][56][57], electrical transport measurements [58][59][60], and angle-resolved photo emission spectroscopy (ARPES) [61] are commonly used in AGNR research. Experimental measurements verify the existence of energy bandgaps, electronic states, and other electronic and transport phenomena predicted by theory. The band gap of AGNR is determined by the energy difference between the conduction band and the highest valence band at the Γ point. We assume that the edge relaxation and the edge-induced correction in pseudopotential remain unchanged for AGNRs with N = 3 ∼ 12. The corresponding numbers of dimer lines in these AGNRs are N d = 2N + 1 = 7 ∼ 25. The calculated band energies for the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) for AGNRs as functions of N d are shown in Figure 14. Here, we only consider the cases with odd N d so that the x-mirror symmetry holds. The computation time is less than 10 s for each case shown.
According to simple models that impose rigid boundary conditions on the edges of AGNRs, the band gap becomes zero when N d = 3m + 2, where m is a positive integer [20,62]. Our calculations show that the band gaps indeed shrink to relatively small values at N d = 11, 17, 23. However, at N d = 11 and 17, the band gap is still appreciable (~0.14 eV and 0.13 eV). This is because the boundary conditions become not so rigid in the realistic situation. Our calculated value of~0.14 eV is consistent with the small band gaps (~0.18 eV) observed in AGNRs of similar dimensions [51]. The band gap obtained by our calculation for the N d = 7 case is 1.43 eV, which is close to the result of 1.47 eV obtained by a previous DFT calculation [63] and also consistent with the experimental result of~1.4 eV based on Fourier-transformed scanning tunneling spectroscopy [52]. However, our result is 0.9 eV, lower than that obtained by the dI/dV measurements for the 7-AGNR, which shows a band gap~2.3 eV [52]. Since our calculation does not include the many-body effect for quasiparticle excitation (nor does the DFT calculation without GW correction), our result is expected to be significantly lower than the DFT-GW calculation [54], which predicts a band gap of~3.94 eV when the many-body effect is included via the GW approximation [7]. However, most dI/dV measurements were taken for AGNR samples placed on Au substrate [51][52][53][54][55][56][57], and there is a strong screening effect on the many-body effect, which explains the large difference in band gap predicted by DFT-GW calculation and experimental results [63]. For the case of 9-AGNR, the band gap obtained by the current SEPM is 0.68 eV, which is 0.7 eV lower than the 9-AGNR band gap measured by dI/dV measurements (~1.4 eV) [54] due to the many-body effect. For the 13-AGNR, our calculation predicts a band gap of 0.9 eV, and the LUMO level is ~ 0.8 eV lower than the 7-AGNR (See Figure 14). This energy difference in LUMO levels between 7-and 13-AGNRs is close to the result (~0.7 eV) obtained by the dI/dV measurements on a 7-13 AGNR heterojunction [22]. Since we adopted a B-spline basis along the z-axis, our calculated band energies are absolute values with respect to the vacuum level. Therefore, it is meaningful to make such a comparison. For the case of 15-AGNR, the band gap obtained by the current SEPM is 0.49 eV, which is ~0.5 eV lower than the dI/dV measurements of 1.03 eV [56]. For the case of 21-AGNR, the SEPM band gap is 0.32 eV, which is ~0.4 eV lower than the value of 0.7 eV obtained by dI/dV measurements [57]. Therefore, in general, the band gaps for AGNRs predicted by the current SEPM are fairly close to those obtained by other DFT calculations (without including the GW correction), and the values are consistently lower than the dI/dV measurements for AGNRs on Au substrate by an amount that varies from 0.4 eV at = 21 to 0.9 eV at = 7. Including the GW correction tends to obtain much larger band gaps in comparison to the experimental results [52], which can be attributed to the screening effect from the metal substrate on the GW correction. Thus, it is necessary to carry out DFT-GW calculations for AGNRs on Au substrate, as reported in [63], in order to determine the amount of band-gap correction due to the many-body effect.

Conclusions
We developed a semi-empirical pseudopotential (SEP) method that is easy to implement and capable of obtaining accurate band structures for graphene and armchairedged graphene nanoribbons. For the -bonded states (with odd symmetry with respect to the z-axis mirror), our SEP method can nearly reproduce all salient features of the DFT results with good accuracy. The time needed to compute the whole band structure associated with -bonded states of graphene is only about a few seconds on a personal For the case of 9-AGNR, the band gap obtained by the current SEPM is 0.68 eV, which is 0.7 eV lower than the 9-AGNR band gap measured by dI/dV measurements (~1.4 eV) [54] due to the many-body effect. For the 13-AGNR, our calculation predicts a band gap of 0.9 eV, and the LUMO level is~0.8 eV lower than the 7-AGNR (See Figure 14). This energy difference in LUMO levels between 7-and 13-AGNRs is close to the result (~0.7 eV) obtained by the dI/dV measurements on a 7-13 AGNR heterojunction [22]. Since we adopted a B-spline basis along the z-axis, our calculated band energies are absolute values with respect to the vacuum level. Therefore, it is meaningful to make such a comparison. For the case of 15-AGNR, the band gap obtained by the current SEPM is 0.49 eV, which is~0.5 eV lower than the dI/dV measurements of 1.03 eV [56]. For the case of 21-AGNR, the SEPM band gap is 0.32 eV, which is~0.4 eV lower than the value of 0.7 eV obtained by dI/dV measurements [57]. Therefore, in general, the band gaps for AGNRs predicted by the current SEPM are fairly close to those obtained by other DFT calculations (without including the GW correction), and the values are consistently lower than the dI/dV measurements for AGNRs on Au substrate by an amount that varies from 0.4 eV at N d = 21 to 0.9 eV at N d = 7. Including the GW correction tends to obtain much larger band gaps in comparison to the experimental results [52], which can be attributed to the screening effect from the metal substrate on the GW correction. Thus, it is necessary to carry out DFT-GW calculations for AGNRs on Au substrate, as reported in [63], in order to determine the amount of band-gap correction due to the many-body effect.

Conclusions
We developed a semi-empirical pseudopotential (SEP) method that is easy to implement and capable of obtaining accurate band structures for graphene and armchair-edged graphene nanoribbons. For the π-bonded states (with odd symmetry with respect to the z-axis mirror), our SEP method can nearly reproduce all salient features of the DFT results with good accuracy. The time needed to compute the whole band structure associated with π-bonded states of graphene is only about a few seconds on a personal computer (PC), and it takes around 100 s to compute the whole band structure for π-bonded states of an AGNR with supercell size of 16a × √ 3a. Thus, it will be a highly efficient method for modeling AGNR-related devices when the π-bonded states play the primary role. For the modeling of AGNR devices that involve the indirect valley minimum at the zone boundary, the edge-dimer antibonding states (which possess even symmetry with respect to the z-axis mirror) will be needed. For such a case, the current method becomes less efficient. It will take about one hour on a PC to obtain the nearly correct energy and dispersion of the localized σ-bonding and σ-antibonding states of the edge dimmers of the AGNR, while the corresponding DFT calculation will take more than one day.
The SEP method with the same edge-induced correction potential given in Equation (51) is also applicable for AGNRs of sizes different from the example used here since the dimer formation on AGNR edges should not be significantly affected by the width of the interior region. Thus, the current model is suitable for modeling any size of AGNR (as long as the width is at least 3a) for device applications.
The current approach can still be improved by adding localized basis functions at the edges to make the convergence of σ-bonding and antibonding states at AGNR edges much faster. Furthermore, an SEP for the other popular graphene nanoribbons with zigzag edges will also be useful. The current approach can also be extended to develop similar SEPs for transition-metal dichalcogenides (TMDs) and related moiré superlattices. All these improvements are worthy topics for future research. For application to TMDs, some complications will occur due to the more complicated structure factors, which cannot be made real. Therefore, finding semi-empirical potentials to fit both the real and imaginary parts of the complex quasi-2D form factors, ∼ V loc (z, G), for two kinds of atoms will require more tedious procedures, although it can still be done. Once it is done, the more interesting application is to develop an SEPM for TMD moiré superlattices considering their high scientific impact. Here, we have demonstrated the success of using graphene eigenstates at different g j points enclosed in the graphene Brillouin zone (BZ) (see Figure 8b) as a contracted basis set to calculate the AGNR band structures efficiently. The same idea can be applied to deal with twisted bilayer TMDs with a large supercell. Since there are no edge states to deal with in the moiré superlattices, we expect that only a small number of TMD bands need to be included in the contracted basis set. By using this contracted basis, the miniband structures of TMD moiré superlattices can be calculated very efficiently with good accuracy (similar to the π-bonded states in AGNR). The wavefunctions of AGNRs calculated by SEPM are very close to the DFT results. Therefore, they can be used to calculate the electron-phonon scattering with good accuracy. Furthermore, since all our AGNR wavefunctions (including edge states) are written as linear combinations of graphene Bloch states, we can relate the electron-phonon scattering matrix elements to those for graphene, and it will be convenient to model the transport properties in AGNR-related optoelectronic devices. The software developed here will be valuable for IC designs of 2D material-based nanoelectronics devices that will be of interest to the semiconductor industry.
The close interplay between theoretical and experimental studies fosters a deeper understanding of AGNRs. Moreover, experimental data provide feedback to refine and improve theoretical methods. The collaboration between theorists and experimentalists allows for the identification of new phenomena and the validation of theoretical models. Together, theoretical and experimental studies offer a comprehensive perspective on AGNRs' properties and pave the way for potential applications in nanoelectronics and beyond. Here, γ = S, L, b represent the short-range (S), long-range (L) range, and bond-charge (b) contributions to ∆V, and D γ (ρ) is related to ∼ D γ (q) by the Fourier transform. S γ (q) = S N R (q) for γ = S or L, and the structure factor for the bond-charge contribution ( γ = b) is as follows: where τ h denotes the center positions of hexagon cells inside the AGNR. Similarly, the matrix elements of the nonlocal pseudopotential are as follows: where Z ν,k+g j iG denotes the eigenvectors obtained by solving Equation (43) for the ν-th band at wavevector k + g j . The kernelV nl in the above has a separable form, as given by Equation (A4), with K = k + g j + G and K = k + g j + G and the index σ running through atoms in the nanoribbon region, A NR . Thus, the sum can be evaluated efficiently. We obtain the following: where ∑ σ∈Gr indicates a sum over two graphene atoms closest to the origin. ∼ Z νn σlm k + g j = ∑ i,G Z ν,k+g j iG P in lm k + g j + G e −i(k+g j +G)·τ σ (A20) and P in lm k + g j + G is given in Equation (A2). ∼ S NR (q) = 1 A SC ∑ s∈A NR e iq·R s , with the lattice vectors R s running through all graphene unit cells within the area A NR . Here, R s = 0 corresponds to the origin marked by a red cross in Figure 8a.