# Third nearest neighbor parameterized tight binding model for graphene nano-ribbons

Van-Truong Tran, Jérôme Saint-Martin, Philippe Dollfus, Sebastian Volz

► To cite this version:

Van-Truong Tran, Jérôme Saint-Martin, Philippe Dollfus, Sebastian Volz. Third nearest neighbor parameterized tight binding model for graphene nano-ribbons. AIP Advances, 2017, 7 (7), 10.1063/1.4994771. [hal-01631290](https://hal.science/hal-01631290)

HAL Id: hal-01631290

https://hal.science/hal-01631290

Submitted on 24 Jul 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Third nearest neighbor parameterized tight binding model for graphene nano-ribbons

Cite as: AIP Advances **7**, 075212 (2017); https://doi.org/10.1063/1.4994771
Submitted: 09 March 2017. Accepted: 14 July 2017. Published Online: 24 July 2017

Van-Truong Tran ID, Jérôme Saint-Martin, Philippe Dollfus, and Sebastian Volz

## COLLECTIONS

Paper published as part of the special topic on **Chemical Physics, Energy, Fluids and Plasmas, Materials Science and Mathematical Physics**

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

View Online

<figure><img src="image_4.png" /></figure>

Export Citation

<figure><img src="image_5.png" /></figure>

CrossMark

## ARTICLES YOU MAY BE INTERESTED IN

* A generic tight-binding model for monolayer, bilayer and bulk MoS₂

AIP Advances **3**, 052111 (2013); https://doi.org/10.1063/1.4804936

* The quantum spin Hall effect and topological insulators

Physics Today **63**, 33 (2010); https://doi.org/10.1063/1.3293411

* Band gap engineering in finite elongated graphene nanoribbon heterojunctions: Tight-binding model

AIP Advances **5**, 087121 (2015); https://doi.org/10.1063/1.4928450

NEW!

Sign up for topic alerts

New articles delivered to your inbox

AIP
Publishing
# Third nearest neighbor parameterized tight binding model for graphene nano-ribbons

Van-Truong Tran,<sup>1,ª</sup> Jérôme Saint-Martin,<sup>2</sup> Philippe Dollfus,<sup>2</sup> and Sebastian Volz<sup>1,3</sup>

(Received 9 March 2017; accepted 14 July 2017; published online 24 July 2017)

The existing tight binding models can very well reproduce the *ab initio* band structure of a 2D graphene sheet. For graphene nano-ribbons (GNRs), the current sets of tight binding parameters can successfully describe the semi-conducting behavior of all armchair GNRs. However, they are still failing in reproducing accurately the slope of the bands that is directly associated with the group velocity and the effective mass of electrons. In this work, both density functional theory and tight binding calculations were performed and a new set of tight binding parameters up to the third nearest neighbors including overlap terms is introduced. The results obtained with this model offer excellent agreement with the predictions of the density functional theory in most cases of ribbon structures, even in the high-energy region. Moreover, this set can induce electron-hole asymmetry as manifested in results from density functional theory. Relevant outcomes are also achieved for armchair ribbons of various widths as well as for zigzag structures, thus opening a route for multi-scale atomistic simulation of large systems that cannot be considered using density functional theory. © 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4994771]

## I. INTRODUCTION

Density functional theory (DFT) and tight binding (TB) method are widely used to investigate and predict various properties of materials, from electronic to phononic, thermoelectric and optical properties.<sup>1-6</sup> While the former technique does not require any empirical input parameters as it is derived directly from first principles, the latter needs several parameters such as onsite energy, hoping energy and eventually overlap terms to construct the Hamiltonians.<sup>5,7</sup> Although DFT usually provides relevant results compared to experimental data, it remains computationally very expensive and therefore its use is limited to small-size structures from few to few hundreds of atoms.<sup>8</sup> In contrast, TB models do not require self-consistent procedures to get the band structures, it hence requires much less computational resources. Consequently, TB models can be implemented to investigate large structures with up to millions of atoms. Additionally, in some specific cases TB calculation can lead to analytical expressions which are very convenient to deepen the analysis of materials properties.<sup>9-11</sup> Thus DFT and TB methods have their own advantages according to the desired level of accuracy and the size of the system.

The empirical TB parameters are usually generated by fitting TB calculations with DFT<sup>12,13</sup> or experimental data.<sup>14</sup> Although the first study of graphene band structure was done by Wallace in 1947,<sup>15</sup> the concept of parameterized TB model was just introduced in 1954 by Slater and Koster.<sup>7</sup>

<footnote n="1"><sup>1</sup>EM2C, CentraleSupélec, Université Paris Saclay, CNRS, 92295 Châtenay Malabry, France</footnote>

<footnote n="2"><sup>2</sup>C2N, CNRS, Université Paris-sud, Université Paris Saclay, 91405 Orsay, France</footnote>

<footnote n="3"><sup>3</sup>LIMMS, Institute of Industrial Science, University of Tokyo, CNRS-IIS UMI2820, 4-6-1 Komaba Meguro-Ku, Tokyo 153-8505, Japan</footnote>

<footnote><sup>ª</sup>vantruongtran.nanophys@gmail.com</footnote>
To calculate the energy bands of 2D graphene in low energy regions, it has been shown that the $2p_z$ orbital TB model is relevant and allows us to describe important electronic properties of this material.¹⁶ In the case when only the nearest neighbor interaction (1NN) is taken into account, the nearest neighbor hoping energy to be considered is $t_1 = 2.7$ eV.¹⁷ This simple 1NN TB model is basically accurate around the Dirac points (K points), but it poorly mimics DFT or experimental data in high-energy regions.¹²

In 2002, Reich *et al*¹² have shown that the match between TB and DFT results for 2D graphene can be substantially improved by introducing a third nearest neighbor (3NN) TB model, including overlap terms. Indeed, the set of TB parameters proposed by Reich fits DFT predictions over a large range of energies. However, the solution for such a set is certainly not unique and in this model the values of the parameters seem to be fitted mathematically rather than physically as hoping and overlap terms of the third nearest neighbors are even larger than that of the second ones. Realizing this problem, Kundu¹⁶ proposed to re-fit DFT data of Reich's paper and introduced a new set in which hoping and overlap terms decay as the neighbor distance increases.

For the ribbon form of graphene, the simple 1NN TB model predicts that armchair graphene nanoribbons (AGNRs) of the $3p + 2$ group are semi-metallic, while AGNRs of other groups $3p + 1$ and $3p$ are semi-conducting.¹⁸⁻²⁰ However, DFT calculations²⁰ and experiments²¹ have shown that all AGNRs are semiconducting. To theoretically explain the semiconducting behavior of AGNRs, Son *et al*²⁰ have introduced an edge deformation (ED) effect in the edges of armchair ribbons into the 1NN TB model (we note this model as 1NN + ED) and an additional term was introduced to increase the hoping energy between edge atoms up to about 12%. The added effect indeed corresponds to the underlying physics of the bandgap opening in group $3p + 2$ although it still can not replicate accurately the width of the bandgap and the slope of energy bands. In 2008, Gunlycke and White²² have improved the 1NN + ED model by introducing an additional term as third nearest neighbor hoping parameter $t_3$ and thus constructed a 3NN + ED model. This model can accurately reproduce the bandgap of AGNRs in most cases but still has discrepancies with the DFT data in the high-energy region of the conduction band. More important, both Son's and Gunlycke's models always present a symmetry between conduction and valence bands (electron-hole symmetry), while DFT results show that electron-hole pair is asymmetrical. When we used Reich's or Kundu's sets for ribbon band structure calculation, the bandgap was underestimated because these sets are not optimized for ribbon structures. It is worth noting that Hancock *et al*²³ have also introduced a Hubbard model up to third nearest neighbor interactions that replicate DFT results satisfyingly. However, Hubbard's models always require a self-consistent procedure and are only necessary for spin-based studies. If the Hubbard term in this model is ignored, the results no longer fit with DFT data.

Hence, it is still required to build a robust set of TB parameters that can adequately reproduce not only bandgaps but also the band shapes of nanoribbons predicted by DFT, using a simple approach without the need of a self-consistent procedure. Since the bandgap is directly linked with on/off states of electronic devices²⁴ and the slope of the bands defines the group velocities and the effective masses of carriers,²⁵,²⁶ the accuracy of the band structure calculation may thus affect the results and conclusions in a large range of transport problems. In addition, the impact of each TB parameter on the band structure is still indeterminate and needs to be revealed.

In the present work, by implementing both DFT and TB calculations for AGNRs, we deduce and propose a new set of parameters for 3NN TB models that presents excellent agreement with DFT results in most cases, even in high-energy regions. Furthermore, although the new set was adapted for narrow AGNRs, we show that this set precisely describes large armchair ribbons and zigzag structures as well. This modeling scheme thus opens a route for accurate atomistic simulation of large systems.

## II. MODELING AND METHODOLOGY

The sketch of both armchair and zigzag graphene nanoribbons is shown in Fig. 1, where *M* characterizes the width of the ribbon and refers to the number of dimer lines along the width of the armchair ribbons and the number of chain lines in zigzag ribbons. The red rectangles indicate
<figure><img src="image_6.png" /><figcaption>FIG. 1. Two typical ribbon structures with (a) armchair edges and (b) zigzag edges. The width of the ribbons is characterized by the parameter M referring to the number of dimer lines in armchair ribbons and to the number of chain lines in zigzag ribbons.</figcaption></figure>

primary cell in each structure and the circles indicate the ranges over which an atom interacts with its neighbors. $a_0$ denotes the nearest distance between two carbon atoms ($a_0 \sim 0.142$ nm), the distance to the first, second and third nearest neighbor atoms are referred to as $r_1 = a_0$, $r_2 = \sqrt{3}a_0$, and $r_3 = 2a_0$, respectively.

In this work, DFT computations were implemented within the QUANTUM ESPRESSO suite of codes²⁷ in the framework of the plane wave basis set while TB predictions were obtained with our house-made code.

In DFT calculations, ribbon edges were passivated by hydrogen atoms to avoid any unexpected states generated inside the bandgap due to charge transfer induced by edge dangling bonds. In all calculations, we have used the Perdew-Zunger (PZ) exchange-correlation functional²⁸ and a norm-conserving Hartwigsen-Goedeker-Hutter pseudopotential²⁹ within the local density approximation (LDA).³⁰ A kinetic energy cutoff of 90 Ry was chosen to safely converge total energies. To mesh the Brillouin zone for integrals, a Monkhorst-Pack 40×1×1 was used.³¹ All structures were relaxed until the force on each atom was less than 0.001 Ry·au⁻¹.

Regarding TB calculations, we started with general Bloch wave functions in periodic structures:⁷

$$|\Psi_{\vec{k}}(\vec{r})\rangle = \frac{1}{\sqrt{N}} \sum_{\beta=-\infty}^{+\infty} \sum_{j=1}^{P} e^{i\vec{k}\cdot\vec{R}_{\beta}} \langle c_j | \varphi_j(\vec{r}-\vec{R}_{\beta}) \rangle, \quad (1)$$

where N is the number of unit cells in the crystal, P is number of atoms in a unit cell and $|\varphi_j(\vec{r}-\vec{R}_{\beta})\rangle$ is the $2p_z$ orbital of the j-th atom in the β-th unit cell.

The time-independent Schrödinger equation is commonly written as

$$H |\Psi_{\vec{k}}(\vec{r})\rangle = E |\Psi_{\vec{k}}(\vec{r})\rangle \quad (2)$$

Substituting (1) into (2) and multiplying both sides by $\langle\varphi_i(\vec{r}-\vec{R}_\alpha)|e^{-i\vec{k}\cdot\vec{R}_\alpha}$ yields

$$\sum_{\beta=-\infty}^{+\infty} e^{i\vec{k}\cdot(\vec{R}_{\beta}-\vec{R}_{\alpha})} \sum_{j=1}^{P} H_{i\alpha,j\beta} c_j = E \sum_{\beta=-\infty}^{+\infty} e^{i\vec{k}\cdot(\vec{R}_{\beta}-\vec{R}_{\alpha})} \sum_{j=1}^{P} S_{i\alpha,j\beta} c_j, \quad (3)$$

where $H_{i\alpha,j\beta} = \langle\varphi_i(\vec{r}-\vec{R}_\alpha)|H|\varphi_j(\vec{r}-\vec{R}_\beta)\rangle = -t_{i\alpha,j\beta}$ is the Hamiltonian element directly associated to the hoping coupling between atom i of the α-th unit cell and atom j of the β-th unit cell and $S_{i\alpha,j\beta} = \langle\varphi_i(\vec{r}-\vec{R}_\alpha)|\varphi_j(\vec{r}-\vec{R}_\beta)\rangle = s_{i\alpha,j\beta}$ is the overlap of the two wave functions. In a general 3NN TB parameterized model, $\{t_{i\alpha,j\beta}, s_{i\alpha,j\beta}\}$ will be fitted to $\{t_1, s_1\}, \{t_2, s_2\}$ or $\{t_3, s_3\}$ depending on the distance between atoms i and j. In the case where $\alpha \equiv \beta$ and $i \equiv j$, $H_{i\alpha,i\alpha} = \langle\varphi_i(\vec{r}-\vec{R}_\alpha)|H|\varphi_i(\vec{r}-\vec{R}_\alpha)\rangle = E_{2p}$ is the $2p_z$ onsite energy of a carbon atom, and obviously $S_{i\alpha,i\alpha} = \langle\varphi_i(\vec{r}-\vec{R}_\alpha)|\varphi_i(\vec{r}-\vec{R}_\alpha)\rangle = 1$.
Combining the P equations constructed from equation (3) (as i = 1:P), a matrix equation can be formed as

$$ \left( H_{\alpha\alpha} + \sum_{\beta \neq \alpha} H_{\alpha\beta} e^{i\vec{k} \cdot (\vec{R}_{\beta} - \vec{R}_{\alpha})} \right) \phi_0 = E \cdot \left( S_{\alpha\alpha} + \sum_{\beta \neq \alpha} S_{\alpha\beta} e^{i\vec{k} \cdot (\vec{R}_{\beta} - \vec{R}_{\alpha})} \right) \phi_0, \quad (4) $$

where $H_{\alpha\beta} = \{H_{i\alpha,j\beta}\}$, $S_{\alpha\beta} = \{S_{i\alpha,j\beta}\}$ are the matrices containing all interactions of atoms between the two $\alpha$-th and $\beta$-th cells and $\phi_0 = (c_1 c_2 \dots c_p)^T$.

Setting $H = H_{\alpha\alpha} + \sum_{\beta \neq \alpha} H_{\alpha\beta} e^{i\vec{k} \cdot (\vec{R}_{\beta} - \vec{R}_{\alpha})}$ and $S = S_{\alpha\alpha} + \sum_{\beta \neq \alpha} S_{\alpha\beta} e^{i\vec{k} \cdot (\vec{R}_{\beta} - \vec{R}_{\alpha})}$ leads to the Eigenvalue problem which provides the band structure:

$$ (S^{-1}H) \phi_0 = E \phi_0 \quad (5) $$

To obtain the density of states (DOSs) in the frame of the TB method, we used the Gaussian smearing of the delta function, i.e.:³²

$$ D(E) = \sum_{n} \sum_{\vec{k} \in BZ} \delta(E - E_n(\vec{k})) = \sum_{n} \sum_{\vec{k} \in BZ} \frac{1}{\eta\sqrt{\pi}} e^{-\frac{(E-E_n(\vec{k}))^2}{\eta^2}}, \quad (6) $$

where n refers to the band index and η to a small positive number.

To reach a high resolution, we used in both DFT and TB calculations a 1000×1×1 k-mesh grid for DOS calculations.

### III. RESULTS AND DISCUSSION

#### A. Assessing existing sets of TB parameters for ribbon calculations

We first examine the merit of existing sets of TB parameters for ribbon structures. In Fig. 2 we display the band structure and the DOS calculated by DFT and different TB models. Two different groups of sets of TB parameters were distinguished: in Fig. 2(a), we employed three different sets of TB parameters which have been fitted for a 2D graphene sheet, while in Fig. 2(b), two sets of TB parameters proposed for AGNRs by Son and Gunlycke, respectively, were used. The parameter values for each set are reported in Table I.

In Fig. 2(a), the simple 1NN TB model with only one nearest-neighbor hoping $t_1$ leads to energy bands (solid pink lines) with a zero bandgap in an armchair ribbon of width $M = 5$. This gapless characteristic has been also predicted for all ribbons of group $M = 3p + 2$ as reported in Refs. 18,19 In contrast, results from the 3NN TB models of both Reich's¹² and Kundu's¹⁶ sets indicate that the structure is a semiconductor, which is thus in agreement with conclusions of DFT data.²⁰ However, both sets lead to smaller gaps compared to that deduced from DFT, which is reflected clearly in the inset of the panel of the DOSs.

Kundu's set was shown to give results in agreement with the DFT band structure of 2D graphene, and particularly to be more physically grounded than the set of Reich, with hoping and overlap parameters decaying for longer neighbor distances. However, observing the energy bands and the DOS shows that Reich's set leads to results that conform to DFT results better than Kundu's set in the armchair ribbon structure with $M = 5$. The weakness of Kundu's set when applied to ribbons may be related to the fact that rule of decay versus distance of hoping and overlap parameters are different for atoms inside the ribbons and for the ones near the edges. We hence should only consider the average effect of this rule.

We also observe a discrepancy in the vicinity of the bandgap when applying these two sets to ribbons of different widths (not shown). In fact, these sets cannot accurately reproduce the bandgap together with the shape of the bands in ribbon structures because they have been fitted for 2D graphene sheet only. As a consequence, they do not include properly the finite size effects of ribbons which strongly influence the bandgap.

It is also worth noting that the highest peaks appearing in the DOS are due to flat bands resulting from the simple 1NN TB model.
<figure><img src="image_7.png" /><figcaption>FIG. 2. Comparison of energy bands and density of states (DOSs) for an armchair ribbon of width *M* = 5 calculated by DFT and TB method. (a) Using sets of TB parameters fitted for 2D graphene. (b) Using sets of TB parameters fitted for graphene ribbons.</figcaption></figure>

To overcome the failure of the 1NN TB model in reproducing the bandgap of AGNRs, Son²⁰ has proposed to introduce a deformation effect at the edges of the ribbon and added a term in the model to describe this effect. In this corrected model, the nearest hoping coupling at the edge increases by 12% compared to the one in bulk. As a result, this effect induces a finite bandgap in ribbons of group

<table><caption>TABLE I. Different sets of tight binding parameters for 2D and ribbon graphene structures.</caption><thead><tr><th>Fit for structure</th><th>Set</th><th><i>E</i><sub>2p</sub> (eV)</th><th><i>t</i><sub>1</sub> (eV)</th><th><i>t</i><sub>2</sub> (eV)</th><th><i>t</i><sub>3</sub> (eV)</th><th><i>s</i><sub>1</sub></th><th><i>s</i><sub>2</sub></th><th><i>s</i><sub>3</sub></th><th>Δ<i>t</i><sub>1</sub> (eV)</th></tr></thead><tbody><tr><td rowspan="4">2D Graphene</td><td>1NN</td><td>0</td><td>2.7</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>3NN (Reich2002)</td><td>-0.28</td><td>2.97</td><td>0.073</td><td>0.33</td><td>0.073</td><td>0.018</td><td>0.026</td><td></td></tr><tr><td>3NN (Kundu2011)</td><td>-0.45</td><td>2.78</td><td>0.15</td><td>0.095</td><td>0.117</td><td>0.004</td><td>0.002</td><td></td></tr><tr><td>1NN+ED (Son2006)</td><td>0</td><td>2.7</td><td></td><td></td><td></td><td></td><td></td><td>0.12 * <i>t</i><sub>1</sub></td></tr><tr><td rowspan="3">Graphene ribbons</td><td>3NN+ED</td><td>0</td><td>3.2</td><td></td><td>0.3</td><td></td><td></td><td></td><td>0.0625*<i>t</i><sub>1</sub></td></tr><tr><td>(Gunlycke2008)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>3NN (this work)</td><td>-0.187</td><td>2.756</td><td>0.071</td><td>0.38</td><td>0.093</td><td>0.079</td><td>0.070</td><td></td></tr></tbody></table>
$M = 3p + 2$ as confirmed in the left and the right panels of Fig. 2(b) (yellow lines). Essentially, this correction turns a semi-metallic ribbon into a conducting one in agreement with DFT and experimental results. However, the accuracy of the bandgap and the band shape still needs to be improved. As can be clearly observed in the inset of Fig. 2(b), the DOS arising from the 1NN TB model shows a mismatch with DFT results around the bottom of the first conduction band ($E_{1c}$) and the top of the first valence band ($E_{1v}$).

In the model of Gunlycke,²² the calculation of bandgap was greatly improved as the authors introduced an additional term corresponding to the third nearest neighbor hoping energy. As it will be shown later, this added term is a pertinent choice in comparison with solely adding ED effects, and makes TB predictions closer to the DFT ones. From the DOS lines, we note that the 3NN TB model using Gunlycke's set gives a value of bandgap in very good agreement with that of DFT in this armchair structure. However, some discrepancies in conduction bands remain noticeable even in the low energy range from 0.4 eV to 1eV.

More important, both TB models for ribbons yield a symmetry between conduction and valence bands with respect to the middle of the bandgap. However, DFT clearly renders an electron-hole asymmetry, which cannot be reproduced by these sets of TB parameters.

A better set is therefore needed not only to reproduce the bandgap but also the electron-hole asymmetry and the shape of bands that plays a significant role in high-energy transport problems.

## B. Impact of TB parameters on the energy bands

To understand which term must be introduced to make predictions of TB method closer to the DFT ones, we first examine the difference between, on one hand, conduction bands and valence bands resulting from DFT and, on the other hand, the outcomes of the simple 1NN TB model (with only $t_1$ involved). Due to the multiplicity of the bands, we simplify the calculations by only considering the difference in the lowest conduction band ($\Delta E_{1v}$) and the highest valence band ($\Delta E_{1c}$) as they are the bands that most contribute to transport.

As it can be seen in Fig. 3(a), the first conduction band obtained from DFT calculation (red line with open circles) is higher than the one given by the 1NN TB model in a short range of energy in the vicinity of the Gamma point. But it becomes lower than its counterpart of the 1NN TB model for

<figure><img src="image_8.png" /><figcaption>FIG. 3. (a) Difference of energy in the first conduction ($\Delta E_{1c}$) and the first valence ($\Delta E_{1v}$) bands between DFT and simple 1NN TB calculations. (b) The quantities $\Delta E_{1c}$ and $\Delta E_{1v}$ were calculated to show the impact of each single TB parameter, as the difference between results of [1NN + 1 parameter] and 1NN TB calculations. All calculations were performed for the armchair ribbon of width $M = 5$.</figcaption></figure>
higher values of k-points. An inverse result is observed for the first valence band. These considerations can also be deduced from Fig. 2(a).

By adding different terms into the 1NN TB model, we may understand which parameters are the most relevant to mimic properly the outcomes of DFT calculation.

We consequently introduced six new sets of two parameters with set 1 being constructed by introducing an edge relaxation $\Delta t_1 = 0.12t_1$ which is nothing but Son's model. Other sets correspond to the addition of $t_2, t_3, s_1, s_2$ or $s_3$, respectively. It is worth to note that by comparing the first and third nearest neighbour distances in the TB description with the radial cut-off of the p orbitals in the localized orbital DFT descriptions,¹² it can be understood that the addition of each TB parameter when introducing further neighbour interactions reflects the longer (and more complete) radial cut-off of the p orbitals considered. The impact of overlap terms is obvious as the orbitals spread over longer regions than the distance between the two nearest lattice sites.

We have calculated again $\Delta E_{1v}$ and $\Delta E_{1c}$ to compare the results induced by these new sets of parameters to the simple 1NN TB model. The results are presented in Fig. 3(b).

At the first glance, a finite gap is obtained when introducing either edge relaxation $\Delta t_1$ (Fig. 3(b1)) or coupling to the third nearest neighbors $t_3$ (Fig. 3(b3)). Adding another term such as $s_1$ or $s_3$ brings almost no change in the bandgap as $E_{1v}$ and $E_{1c}$ are preserved at the Gamma point but strongly alters the bands at high k-point values. Although the adding of $t_2$ leads to a shift of all bands, the equality of $\Delta E_{1v}$ and $\Delta E_{1c}$ at $k=0$ indicates that the bandgap remains the same (equal to 0) as in the 1NN model. The introduction of $t_2$ is actually equivalent to adding a potential energy on each lattice site, which shifts the whole band structure.

Adding $s_2$ also retains the bandgap as $\Delta E_{1c} = \Delta E_{1v} = 0$ at $k=0$. However, for other structures of different widths such as $M=6$ or $M=7$, we observe that adding $s_2$ leads to an opening of the energy gap (not shown).

It has been shown in Ref. 26 that the effects of adding either the edge relaxation term $\Delta t_1$ or the hopping energy $t_3$ are equivalent regarding the opening of a bandgap. However, they are actually equivalent only in the vicinity of $k=0$. Indeed, they both lead to a gap opening, but away from the Gamma point, adding $\Delta t_1$ yields a shift up (down) of the conduction (valence) bands. In contrast, adding $t_3$ shifts down (up) the conduction (valence) bands at high k-point values, with $\Delta E_{1v} < 0$ ($\Delta E_{1c} > 0$). Adding $t_3$ finally appears more relevant to match the results of DFT than introducing $\Delta t_1$. Actually, in the presence of the third neighbor coupling $t_3$ and, as discussed below, of the overlap terms, the edge distortion effect becomes a minor effect that can be safely ignored.

It is also very important to note that Figs. 3(b1), 3(b3), 3(b5) depicting the effect of introducing the ED term $\Delta t_1, t_3$, and $s_2$, respectively, reveals that these parameters retain the symmetrical property between conduction and valence bands, as in the simple 1NN TB model. In contrast, the introduction of $t_2, s_1$ and $s_3$ as presented in Figs. 3(b2), 3(b4), and 3(b6), respectively, can induce an asymmetrical behavior of conduction and valence bands.

As it can be seen from the DFT energy bands in Fig. 2, electron-hole symmetry should not be expected. As a consequence, adding $t_3$ or/and $\Delta t_1$ is not enough to precisely mimic the band shapes predicted by DFT. This outcome justifies why both Son's and Gunlycke's models are failing in reproducing the bands of DFT although they successfully explain the semiconducting behavior of the ribbons of group $M=3p+2$. Accordingly, it is mandatory to also introduce overlap terms in order to optimize the accuracy of the TB technique against DFT calculation.

### C. The new 3NN TB model for graphene ribbons

Following the idea that adding $t_3$ is more relevant than introducing $\Delta t_1$, we use here a general 3NN TB model without $\Delta t_1$ to aim at the best fit between TB calculation and DFT data.

We started with a reasonable guess set of 3NN TB parameters, $E_{2p} = -0.187$ eV being qualitatively accurate to reproduce the DFT bandgap of the structure $M=5$. Then we varied the other six parameters including $t_1, t_2, t_3, s_1, s_2$ and $s_3$, around the guess values, i.e, we selected eleven values around the guess value of each parameter and constructed $11^6$ sets to be scanned. We then calculated the energy error of each band provided by the TB model with each set compared to DFT data. Due to a large number of bands, we considered the error for only the first conduction and valence bands, which are the most relevant to transport properties. The error of each set was taken as the largest error over all k-points
and the set providing the lowest error compared to the other ones was then selected. An advanced step was then added to optimize the value of each term by varying it independently and successively around the obtained value but with smaller increments and we searched for the lowest error compared to the DFT data. The final optimized values are shown in Table I.

Though it may be considered as surprising, the fact that $t_3 > t_2$ is consistent with previous results of Reich ($t_2 = 0.073$ eV, $t_3 = 0.33$ eV) and Gunlycke ($t_2 = 0$). These TB parameters are just fitting parameters that implicitly include complex effects as the influence of other bands and many-body effects, which can make their values counterintuitive.

We employed the parameters of the selected set and plotted the corresponding energy bands together with the DFT ones for comparison. The results are displayed in Fig. 4 where the DOSs were also reported for further comparison. As clearly manifested from the band structure and the DOS panels, TB calculations (dashed red lines) obtained from the set proposed in this work show an excellent agreement with DFT results for energies ranging from -3 eV to 3 eV. Around the region of the bandgap, the inset of the DOSs indicates the high accuracy of TB data as both the bandgap and the DOSs are identical to the DFT results. This new set of parameters seems to be well optimized compared to the existing TB ones proposed by Reich¹² and Kundu¹⁶ fitted for 2D graphene and Son,²⁰ Gunlycke²² fitted for ribbons.

To further validate this new set of parameters, we compare other results of band structure obtained from DFT, 3NN TB with Gunlycke's set and 3NN TB model with the set of this work. First, results for other groups of ribbons, i.e. $M = 6$ (group 3p) and $M = 7$ (group 3p + 1), are shown to support the robustness of the new set.

In Figs. 5(a) and 5(b), the results obtained from our set of parameters are obviously better than the ones of Gunlycke's set and, overall, they fit very accurately with DFT outcomes although both Gunlycke's set and the new set give bandgaps slightly smaller than the one of DFT in the case of $M = 6$. Though the new set was constructed from a fit in a narrow structure with $M = 5$, Figs. 5(c) and 5(d) displaying the results obtained for $M = 11$ and $M = 19$, respectively, show that this set still has strong relevance for larger ribbons. Band structures resulting from our model (red lines) still fit very precisely DFT data (black circles) even in the high energy regions. In the large ribbons, the results of Gunlycke's set (blue) are good in the low energy regions around the gap and in a large part of the valence band. However, it still exhibits a substantial inconsistency with DFT results in the high-energy regions of the conduction bands.

Additionally, we have also used the set for zigzag structures. The band structure for a zigzag ribbon of width $M = 11$ chain lines along the width is displayed in Fig. 6, including both DFT and TB results. The bands resulting from the new TB set show a very good agreement with that of DFT, particularly in the first conduction and valence bands and at the energy points at the boundary of the

<figure><img src="image_9.png" /><figcaption>FIG. 4. (a) Energy bands and (b) Density of states (DOSs) calculated for the armchair ribbon $M = 5$ by DFT and 3NN TB model with the new set of parameters proposed in this work.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 5. Comparison of energy bands calculated by DFT (black filled circle), the 3NN TB model using the set of Gunlycke (blue) and the 3NN TB model using the set proposed in this work (red) for different ribbons of width: (a) $M = 6$, (b) $M = 7$, (c) $M = 11$, and (d) $M = 19$.</figcaption></figure>

Brillouin zone. Gunlycke's set also matches well with DFT predictions for the first valence bands but obviously poorly fits in the high-energy region of the conduction bands. These results reinforce the relevance of the new set of TB parameters proposed in the present work.

Although the set proposed in this work is devoted to GNRs, we have also checked the merit of the set for 2D graphene sheets. The DFT calculation for 2D graphene was performed using the same method as for ribbons (LDA). Then the DFT data were compared with the results obtained from different TB models. The outcome is presented in Fig. 7(a).

At first glance, we see that the bands obtained with the Reich's set are more accurate overall than the others are, which is consistent with the fact that this set was optimized for 2D layers. However, by zooming in on the energy bands around the Dirac point (K point), as shown in Fig. 7(b), we see that the TB results obtained from the set proposed in this work fit accurately the DFT data and are comparable with Reich's results. Thus, in the energy range from -3 eV to 3 eV that is the most relevant for many applications, our set of TB parameters can be used to get accurate results for both 2D and ribbon graphene structures.
<figure><img src="image_11.png" /><figcaption>FIG. 6. Comparison of energy bands for a zigzag graphene ribbon of width $M = 11$ with DFT (black with filled circle), the 3NN TB model using the set of Gunlycke (blue) and the 3NN TB model using the set proposed in this work (red).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 7. (a) The band structure of 2D graphene sheet: Comparison between the DFT data and the ones obtained by different TB models. (b) Zoom in on the low energy region around the Dirac point (K point).</figcaption></figure>

## IV. CONCLUSION

In conclusion, we have analyzed the relevance of commonly used sets of TB parameters for graphene structures. We have shown that although 3NN TB parameters such as Reich's or Kundu's ones have been demonstrated to accurately fit DFT bands of 2D graphene, they are not efficient in reproducing DFT predictions for ribbon structures. Other parameters of Son and Gunlycke were optimized for armchair ribbons but these sets of parameters cannot mimic the electron-hole asymmetry of the energy bands. By implementing both DFT and TB calculations and then introducing a fit, we have shown the high accuracy of a new set of TB parameters including up to 3NN plus overlap terms for graphene ribbons. The new set has been demonstrated to be in excellent agreement with DFT results in most cases, even in high energy regions and for large ribbons. Although the set has been fitted for armchair ribbons, it has been shown to also very well describe zigzag structures and even the energy bands around the Dirac points of 2D graphene.

## ACKNOWLEDGMENTS

This work was supported by the TRANSFLEXTEG EU project.
¹ R. G. Parr and Y. Weitao 352 (1994).

² C. D. Simão, J. S. Reparaz, M. R. Wagner, B. Graczykowski, M. Kreuzer, Y. B. Ruiz-Blanco, Y. García, J.-M. Malho, A. R. Goñi, J. Ahopelto, and C. M. Sotomayor Torres, *Carbohydr. Polym.* **126**, 40 (2015).

³ N. M. Harrison, *Technology 2*, 1 (1995).

⁴ C. Bena and G. M. Baux, *New J. Phys.* **11**, 95003 (2009).

⁵ D. A. Papaconstantopoulos and M. J. Mehl, *J. Phys. Condens. Matter* **15**, R413 (2003).

⁶ V.-T. Tran, J. Saint-Martin, and P. Dollfus, *Nanotechnology* **26**, 495202 (2015).

⁷ J. C. Slater and G. F. Koster, *Phys. Rev.* **94**, 1498 (1954).

⁸ L. E. F. F. Torres, S. Roche, and J. C. Charlier, *Introduction to Graphene-Based Nanomaterials: From Electronic Structure to Quantum Transport* (Cambridge University Press, 2014).

⁹ Y. W. Son, M. L. Cohen, and S. G. Louie, *Phys. Rev. Lett.* **97**, 216803 (2006).

¹⁰ D. Gunlycke and C. T. White, *Phys. Rev. B - Condens. Matter Mater. Phys.* **77**, 1 (2008).

¹¹ N. A. Pike and D. Stroud, *Phys. Rev. B - Condens. Matter Mater. Phys.* **89**, 1 (2014).

¹² S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejón, *Phys. Rev. B* **66**, 1 (2002).

¹³ Y. Tan, M. Povolotskyi, T. Kubis, Y. He, Z. Jiang, G. Klimeck, and T. B. Boykin, *J. Comput. Electron.* **12**, 56 (2013).

¹⁴ M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damjanović, A. Bosak, M. Krisch, and C. Thomsen, *Phys. Rev. B* **76**, 35439 (2007).

¹⁵ P. R. Wallace, *Phys. Rev.* **71**, 622 (1947).

¹⁶ R. Kundu, *Mod. Phys. Lett. B* **25**, 163 (2011).

¹⁷ A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, *Rev. Mod. Phys.* **81**, 109 (2009).

¹⁸ K. Nakada, M. Fujita, G. Dresselhaus, and M. Dresselhaus, *Phys. Rev. B. Condens. Matter* **54**, 17954 (1996).

¹⁹ H. Zheng, Z. Wang, T. Luo, Q. Shi, and J. Chen, *Phys. Rev. B* **75**, 165414 (2007).

²⁰ Y. W. Son, M. L. Cohen, and S. G. Louie, *Phys. Rev. Lett.* **97**, 216803 (2006).

²¹ W. X. Wang, M. Zhou, X. Li, S. Y. Li, X. Wu, W. Duan, and L. He, *Phys. Rev. B - Condens. Matter Mater. Phys.* **93**, 1 (2016).

²² D. Gunlycke and C. White, *Phys. Rev. B* **77**, 115116 (2008).

²³ Y. Hancock, A. Uppstu, K. Saloriutta, A. Harju, and M. J. Puska, *Phys. Rev. B - Condens. Matter Mater. Phys.* **81**, 3 (2010).

²⁴ V. T. Tran, J. Saint-Martin, and P. Dollfus, *Appl. Phys. Lett.* **105**, 73114 (2014).

²⁵ R. Grassi, S. Poli, E. Gnani, A. Gnudi, S. Reggiani, and G. Baccarani, *Solid. State. Electron.* **53**, 462 (2009).

²⁶ P. Zhao, M. Choudhury, K. Mohanram, and J. Guo, *Nano Res.* **1**, 395 (2008).

²⁷ P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. Fabris, G. Fratesi, S. de Gironcoli, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, *J. Phys. Condens. Matter* **21**, 395502 (2009).

²⁸ J. P. Perdew and A. Zunger, *Phys. Rev. B* **23**, 5048 (1981).

²⁹ C. Hartwigsen, S. Goedecker, and J. Hutter, *Phys. Rev. B* **58**, 3641 (1998).

³⁰ W. Kohn and L. J. Sham, *Phys. Rev.* **140**, A1134 (1965).

³¹ H. J. Monkhorst and J. D. Pack, *Phys. Rev. B* **13**, 5188 (1976).

³² C. Pickard and M. Payne, *Phys. Rev. B* **59**, 4685 (1999).