Abstract
We present an efficient approach to study the carrier transport in graphene nanoribbon (GNR) devices using the nonequilibrium Green's function approach (NEGF) based on the Dirac equation calibrated to the tightbinding πbond model for graphene. The approach has the advantage of the computational efficiency of the Dirac equation and still captures sufficient quantitative details of the bandstructure from the tightbinding πbond model for graphene. We demonstrate how the exact selfenergies due to the leads can be calculated in the NEGFDirac model. We apply our approach to GNR systems of different widths subjecting to different potential profiles to characterize their device physics. Specifically, the validity and accuracy of our approach will be demonstrated by benchmarking the density of states and transmissions characteristics with that of the more expensive transport calculations for the tightbinding πbond model.
Keywords:
graphene nanoribbons; Dirac equation; quantum transport; nonequilibrium Green's function1 Introduction
Recent progress of graphene nanoribbon (GNR) fabrication has demonstrated the possibility of obtaining nanoscale width GNRs, which have been considered as one of the most promising active materials for next generation electronic devices due to their unique properties such as bandgap tunability via controlling of the GNR width or subjecting GNR to external electric/magnetic fields [15]. Device simulations play an important role in providing theoretical predictions of device physics and characteristics, as well as in the investigation of device performance, in order to guide the development of future device designs. Due to the nanoscale structures of GNRs, however, semiclassical treatments of carrier transport [6], which are the mainstay of microelectronics, are no longer valid. As a result, quantum transport formalism based on models incorporating detailed atomic structures, such as the abinitio types [79], is needed for the proper simulation of these materials. Unfortunately, a fullfledge abinitio atomistic model for carrier transport simulation is still very computationally expensive and impractical even with the latest stateoftheart computing resources. In this study, we therefore develop an efficient model in which a tightbinding Dirac equation (TBDE), calibrated with parameters from the tightbinding πbond model (TBπ) [1013], is used together with the nonequilibrium Green's function approach (NEGF) [14] to investigate transport properties of GNRs. We compare the density of states, DOS(E), and the transmission, T(E), of selected GNR devices for our TBDE model with that of the more expensive TBπ model. Good agreement is found within the relevant energy range for flat band, Laplace and single barrier bias condition. We believe that our model and calibrated data for a side selection of GNR widths presented in this article provided researchers in the quantum transport an accurate and practical framework to study the properties, particularly quantum transport in arbitrary bias conditions, of GNRbased devices.
2 Model
The Hamiltonian based on the Dirac equation [15,16] for graphene is given as:
where p_{μ }= iћ∂_{μ }is the momentum for the direction μ = {x,y}, v_{F }is the Fermi velocity of graphene at the Dirac points (fixed at 10^{6 }ms^{1}) and U(x) is the onsite potential. Due to the 1D property of GNRs, the finite difference approach can be used along the transport direction (x) of GNRs and the Hamiltonian (h_{n}) at each site n, and its backward (h_{}) and forward (h_{+}) couplings with its neighbors (separated by a uniform spacing l_{0}) for a particular subband mode k_{y}, can be written as:
where l_{0 }is the effective 1D cell size as a result of the discretized Hamiltonian in (2). Figure 1a shows the schematics for realspace graphene and Figure 1b the 1D GNR model associated with (2). For an infinitely long GNR with uniform U_{0}, the Bloch waves solutions are valid and the dispersion relation, E(k_{x},k_{y}), for (2) is
Figure 1. Schematic representation of mapping of (a) a realspace twodimensional GNR to (b) the onedimensional Dirac Equation model with two degrees of freedom per effective cell of length ℓ_{0}.
where for a fixed k_{y }the positive and negative signs denoting the conduction and valence bands, respectively. In the absence of external potential (U_{0 }= 0) and in the limit of large GNR width at which is small, (3) gives the linear dispersion for graphene . The energy bandgap of a certain width, and hence k_{y}, is given by E_{g }= 2ћv_{F}k_{y }at k_{x }= 0.
For nonequilibrium situations, we have to calculate the device retarded Green's function G(E) for a particular energy E for the Hamiltonian in (2). Assuming the potential energies at the equilibrium source and drain are U_{s }and U_{d}, respectively, and there are N lattice points in the device region, the G(E), of matrix size 2N × 2N, is given by G(E) = [EI_{2N } H  Σ_{s } Σ_{d}]^{1}, where the 'selfenergies' Σ_{s }and Σ_{d }are associated with the effects of the semiinfinitely long source and drain [14]. Consider the selfenergy of the drain (specified by the Hamiltonian H_{d }of size 2M × 2M, where M is an arbitrary number of lattice points with spacing l_{0 }spanning the drain), defined in the NEGF framework [14] by , where the drain Green's function, , is also of the size 2M × 2M, and τ_{ }= (τ_{+})^{† }is the coupling matrix (of size 2M × 2N) between the device and drain, which ends and starts at lattice points n =  1 and 0, respectively. However, the only nonzero component of τ_{± }is that of h_{± }across the n =  1 and 0 interface, and hence only the 2 × 2 drain surface Green's function , makes nontrivial contribution to Σ_{d}, i.e., is the only nonzero 2 × 2 submatrix, associated with lattice point n =  1, of Σ_{d }(of size 2N × 2N). Using the identity for the drain region (n ≥ 0), the system of equations for the dimensionless Green's function G can be written as
where ω^{(0) }= EI_{2 } h_{0 }is independent of sites inside the drain with uniform U_{d}. One can iteratively substitute (second term) in (5) with the same in (4) so that after ℓ ≥ 1 number of iterations, (4) and (5) can be rewritten as [17]
where
The prefactor is such that k_{x }is related to E via (3). The integer m ≥ 0 labels the surviving lattice points with spacing 2 ^{ℓ}l_{0}. The effects of the eliminated nodes after ℓ number of iterations are taken into account in terms of "renormalized" couplings α^{(ℓ) }and β^{(ℓ)}, (which happens to be equal in this model) and site energies (ω^{(ℓ) }at site index 2^{ℓ}m with m ≥ 1 and at m = 0, respectively). The symmetries of h_{0 }and h_{± }in (2) resulted in α^{(ℓ)}, ω^{(ℓ) }and each directly proportional to the "bare" energy ω^{(0) }for all ℓ ≥ 1, with their respective coefficients Λ^{(ℓ)}, Ω^{(ℓ)}, and as scalar functions dependent on λ only. We show by induction that for all ℓ ≥ 1,
uniquely satisfy (11), (12), and (13). Since we are interested in the retarded Green's function (i.e., E → E + iη) for an infinitesimally small energy η > 0, the condition imposed on the propagating waves is such that λ ≈ 1  (l_{0}/ћv_{g})η < 1, where v_{g }≡ ћ^{1 }(∂E/∂k_{x}) > 0 is the relevant group velocity [18,19]. Expanding in terms of λ and taking the limit ℓ → ∞, (14), (15), and (16) give A^{(∞) }= 0, Ω^{(∞) }= (1 + λ^{2})/(1  λ^{2}), and , respectively. The exact value of , in the limit of ℓ → ∞ in (6), is now given by
Similar argument can be applied at the sourcechannel interface where the analog sourceside counterpart of takes the same form as (17) with U_{s }replacing U_{d}. Therefore, the only nonzero 2 × 2 submatrices for Σ_{[s, d] }are
In the past, (6)(13) are evaluated iteratively to calculate , and hence Σ_{[s, d] }[13,17]. In this study, we have shown that (6)(13) can be solved analytically for the Dirac form in (2) and that significant computational saving and accuracy can therefore be achieved by directly using (18) instead of numerically iterating (6)(13). Figure 2 shows that the total computing time to calculate all the relevant modes of (E) for E ∈ [1,1] eV with spacing of 0.001 eV via analytical, i.e., (17), and iterative means, i.e., (6)(13) for a range of GNR width on a typical duo core PC using MATLAB. The time needed to calculate using the iterative method is about 40× larger than that of the analytic method over the entire range of the GNR width considered. In general, it is observed that the computing time increases with the GNR width for both analytical and iterative methods because the number of modes also increases with the width. (See Table 1.) Figure 2 also shows, as a comparison, the corresponding total computing time for calculating the all relevant surface Green's functions (via iterative method) for the same set of GNR width in TBπ model. This time is much larger than that of the TBDE, between about 100× (at 1.1 nm width) and 455× (for 3.8 nm width) that of the analytic method of TBDE. Therefore the computational saving from using our analytic results for the surface Green's function, (17), is compelling. The computing saving will be even more apparent in more realistic quantum transport calculations in which the NEGF and Poisson equation are solved iteratively to achieve selfconsistent solutions.
Figure 2. The total computing time for calculating a series of (E) for all relevant modes in  1 ≤ E ≤ 1 eV with 0.001 eV spacing using analytic (○) and iterative (Δ) methods in the TBDE model for different GNR width. The iterative method takes about 40 × longer than that of analytic method. Included for comparison is the total time to calculate the corresponding surface Green's functions calculated using iterative method in the TBπ model (◊).
Table 1. Results of bestfit l_{0 }(for their respective subbands) to be used for our TBDE model for GNRs of different widths
With G(E) now specified, the DOS(E), T(E), line charge density (ρ_{1D}) and total current (I_{t}) can be obtained, respectively [20], via
where , f_{s,d}(E) is the Fermi function at either the source or drain, Σ_{sb }denotes sum over the subbands, Diag[⋯] and Tr[⋯] denote the diagonal and the trace of a square matrix, respectively.
3 Results and discussions
To incorporate the material details of GNR into the TBπ model, we first fit (3) of different GNR widths with that of the TBπ model, which is widely used to calculate the bandstructures of GNR, for a flat potential (i.e., U = 0). Both real and imaginary parts of (3) are fitted for multiple subbands with different values of l_{0 }for a particular GNR system. Figure 3 shows the comparisons of E(k) for the GNRs with width 1.0 nm and 1.4 nm, labeled as W10 and W14, respectively. At larger k, the E(k) calculated using (3) deviated from the that of the TBπ model. This is expected as the TBDE model for GNR is most accurate near the Dirac points at small k[15]. Since we are interested in semiconductor properties of GNRs, only the wide bandgap armchair GNRs (families with indices of m = 3p and 3p+1) [8,21] are considered here. The GNRs associated with m = 3p + 2 have E_{g }that are too small and are not considered here. Table 1 shows the bestfit l_{0 }at different subbands for the m = 3p and m = 3p + 1 GNRs obtained under this study. With these calibrations, the adequate bandstructure details based on TBπ model can be incorporated in the TBDE model. Figure 4 compare the DOS(E) and T(E) for the same W12 and W14 systems using TBDE model (with the fittedlo values from Table 1) and that of the TBπ model. The very good agreements of results between the two models is a good first step to demonstrate the validity of the TBDE model in tackling quantum transport problems at which accurate T(E) and DOS(E) are the keys.
Figure 3. The E(k) calculated from the TBDE matching that of the TBπ with different best fit l_{0 }for different subbands for the (a) W12 and (b) W14 devices. Only conduction bands for E ≥ 0 are shown. The valance bands are symmetric about E = 0.
Figure 4. The DOS(E) and T(E) of (a) W12 and (b) W14 devices with the bestfit l_{0 }(from E(k) calculations) at U = 0 agreeing to that of the TBπ model. Both the DOS(E) and T(E) are symmetric about E = 0.
To apply the NEGFTBDE to more realistic transport situations, one needs to solve the NEGFTBDE under bias potentials. For a Laplace potential (with a bias of 0.3 V), as shown in Figure 5a, the DOS(E) and T(E) for the W14 GNR are shown in Figure 5b,c, respectively. The corresponding TBπ results and that of TBDE model with U = 0 are also included for reference. As shown in Figure 5a, the 0.3 V bias is achieved by shifting the conduction and valence bands upwards relative to those at the drain. As the highest valence bandedge (E_{v}) (at source) shifted up by 0.3 eV, the onset of DOS(E) for E < 0 also creeped up into the original forbidden zone (with U = 0) by about 0.3 eV as indicated by arrow in Figure 5b. The positions of the DOS(E) associated with the higher subbands have also moved up the energy scale relative to those for U = 0. However, the onset of DOS(E) for E > 0 has not been affected significantly by the Laplace setup because the lowest conduction bandedge, which is at the drain, is still intact at E = E_{g}/2. Although the forbidden zone for DOS(E) has narrowed as indicated in Figure 5b, the forbidden zone for T(E) has actually widen, as shown in Figure 5c, with the onset of nonzero T(E) for E > 0 receding upwards by about 0.3 eV as indicated by the arrow, but unchanged T(E) for E < 0. This is because from carriers are only unhindered sourcetodrain only at E > E_{g}/2 + 0.3 eV and E < E_{g}/2. The newly addition of DOS(E) at the sourceside valence has no state of comparable E to connect to in the channel and drain and hence does not contribute to T(E).
Figure 5. The DOS(E) and T(E) of a simple Laplace Potential. (a) Schematic of a simple Laplace potential profile with a bias of 0.3 V across the GNR channel. (b) The resulting DOS(E) versus E with red arrow indicating the new addition of DOS(E) due to the upward movement of valence bandedge by 0.3 eV. (c) T(E) versus E. Results for U = 0 and that calculated by TBπ are also included for comparison.
Next, we subjected the W14 GNR to a rectangular barrier of 0.1 eV in the channel as shown in Figure 6a. The resulting DOS(E) and T(E) are shown in Figure 6b,c, respectively, with that of TBπ model and U = 0 included for comparison. As expected, the onset of both DOS(E) at the conduction and valence ranges have not changed because the lowest E_{c }and highest E_{v}, at E_{g}/2, and E_{g}/2, respectively, have not been changed by the introduction of the barrier potential compared to that of U = 0. However, it is observed that the magnitude of DOS(E) just above E = E_{g}/2 was reduced significantly due to the lost of states in the channel region dominated by the barrier. The inverted well of depth 0.1 eV at the channel valence bandedge is expected to accommodate some discrete bound states. However, the DOS(E) associated with them may be too sharp to be captured, or partially captured by the E grids being used. This expectation is borne out by the inset of Figure 6b, which shows the logscale of the DOS(E) in the vicinity of E = E_{g}/2. Two discrete bound states, with the heights of their DOS(E) partially captured, are discernible within the inverted well energy range of within 0.1 eV above E_{g}/2. As for T(E), the carriers are unhindered sourcetodrain only for E > E_{g}/2 + 0.1 eV and E < E_{g}/2 eV and hence those boundaries marked the onset of T(E), as shown in Figure 6c. The bound states created by the inverted well in the channel region do not contribute to T(E) as there are no states of comparable energies both at the source and drain to connect to them.
Figure 6. The DOS(E) and T(E) of a rectangular barrier. (a) Schematic of a rectangular barrier of height 0.1 eV across the W14 GNR channel. (b) The DOS(E) with red arrow indicating region of reduced DOS(E) due to the introduction of the barrier at channel. The region near E = E_{g}/2 (as indicated) is magnified as inset with DOS(E) in logscale. Two discrete bound states, created by the inverted well at valence bandedge, as shown in (a). (c) The T(E) with red arrow indicates the receding T(E) away from E_{g}/2 due to the 0.1 eV barrier. Results for U = 0 and that calculated from TBπ are also included for comparison.
In both the Laplace and rectangular barrier potential profiles, the DOS(E) and T(E) for our TBDE model are in satisfactory agreement with that calculated from TBπ model within about 1.5 eV around the midgap. At higher energies, significant deviations in the DOS(E) and T(E) are consistent with the discrepancies we observed in E(k) (as shown in Figure 3b), as discussed earlier. Nonetheless, these deviations are limited to the highenergy range that is of little relevance to the electron transport in GNR devices. Therefore, our TBDE approach is expected to be valid and as a practical and efficient alternative to TBπ for studying carrier transport involving arbitrary selfconsistent electrostatic potentials for device simulations [22,23].
4 Conclusion
We developed a tightbinding Dirac equation for practical and accurate numerical investigation of the electron transport in GNR devices. Based on our knowledge, this is the first time that the surface Green's function arises from applying the Dirac equation in NEGF framework is calculated exactly and hence can be used to achieve significant savings in NEGF calculations. The TBDE model is calibrated, with the appropriate parameters (v_{F }= 10^{6 }ms^{1 }and l_{0}), to match the relevant bandstructure details as that of the TBπ model, especially near the Dirac points. The bestfitted l_{0 }for a selected set of GNR widths are also presented for use. We show that the DOS(E) and T(E) calculated by our calibrated TBDE model can produce very good agreement with those that are calculated by the more expensive TBπ model for the flat, Laplace, and rectangular barrier potentials. These cases validate the accuracy of the TBDE model and provided good confidence that the model can be used as a practical and accurate starting point for quantum transport of GNRbased devices where nonequilibrium and arbitrary electrostatic potentials are involved.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
SKC conceived the possibility of an analytic expression for the surface Green's function in Dirac equation and performed the calculations. He also prepared the manuscript. KTL partially wrote the codes for the Dirac equation and together with DS carried out the numerical studies and analysis. GL conceived of the overall study, designed most of the codes, and participated in the analysis. All authors read and approved the final manuscript.
Acknowledgements
This study was supported by Singapore's Agency for Science, Technology and Research (A*STAR) Public Sector Funding Grant No. 0821010023.
References

Kosynkin DV, Higginbotham AL, Sinitskii A, Lomeda JR, Dimiev A, Price BK, Tour JM: Longitudinal unzipping of carbon nanotubes to form graphene nanoribbons.
Nature 2009, 458(7240):872. PubMed Abstract  Publisher Full Text

Jiao L, Zhang L, Wang X, Diankov G, Dai H: Narrow graphene nanoribbons from carbon nanotubes.
Nature 2009, 458(7240):877. PubMed Abstract  Publisher Full Text

Cai J, Ruffieux P, Jaafar R, Bieri M, Braun T, Blankenburg S, Muoth M, Seitsonen AP, Saleh M, Feng X, Mullen K, Fasel R: Atomically precise bottomup fabrication of graphene nanoribbons.
Nature 2010, 466(7305):470. PubMed Abstract  Publisher Full Text

Choi CY, Lee JH, Koh JH, Ha JG, Koo SM, Kim S: Hightemperature stable operation of nanoribbon fieldeffect transistors.
Nanoscale Res Lett 2010, 5:1795. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kumar SB, Jalil MBA, Tan SG, Liang G: Magnetoresistive effect in graphene nanoribbon due to magnetic field induced band gap modulation.
J Appl Phys 2010, 108:033709. Publisher Full Text

Lundstrom M, Guo J: Nanoscale Transistors: Device Physics, Modeling and Simulations. New York: SpringerVerlag; 2006.

Martins TB, Miwa RH, da Silva AJR, Fazzio A: Electronic and transport properties of Borondoped graphene nanoribbons.
Phys Rev Lett 2007, 98(19):196803. PubMed Abstract  Publisher Full Text

Son YW, Cohen ML, Louie SG: Energy gaps in graphene nanoribbons.
Phys Rev Lett 2006, 97:216803. PubMed Abstract  Publisher Full Text

Yang L, Park CH, Son YW, Cohen ML, Louie SG: Quasiparticle energies and band gaps in graphene nanoribbons.
Phys Rev Lett 2007, 99(18):186801. PubMed Abstract  Publisher Full Text

Saito R, Dresselhaus G, Dresselhaus M: Physical Properties of Carbon Nanotubes. London: Imperial College Press; 1998.

Gunlycke D, White CT: Tightbinding energy dispersions of armchairedge graphene nanostrips.

Saito R, Dresselhaus G, Dresselhaus MS: Trigonal warping effect of carbon nanotubes.
Phys Rev B 2000, 61(4):2981. Publisher Full Text

Wu Y, Childs PA: Conductance of graphene nanoribbon junctions and the tight binding model.

Datta S: Quantum Transport: Atom to Transistor, chap. 9. New York: Cambridge University Press; 2005:217251.

Neto AHC, Guinea F, Peres NMR, Novoselov KS, Geim AK: The electronic properties of graphene.
Rev Mod Phys 2009, 81:109. Publisher Full Text

Ouyang Y, Wang X, Dai H, Guo J: Carrie scattering in graphene nanoribbon fieldeffect transistors.
Appl Phys Lett 2008, 92:243124. Publisher Full Text

Sancho MPL, Sancho JML, Rubio J: Highly convergent schemes for the calculation of bulk and surface Green functions.
J Phys F met Phys 1985, 15:851. Publisher Full Text

Velev J, Butler W: On the equivalence of different techniques for evaluating the Green function for a semiinfinite system using a localized basis.
J Phys Condens Matter 2004, 16:R637. Publisher Full Text

Wang JS, Wang J, Lu JT: Quantum thermal transport in nanostructures.
Eur Phys J B 2008, 62:381. Publisher Full Text

Datta S: Quantum Transport: Atom to Transistor, chap. 11. New York: Cambridge University Press; 2005:305307.

Han MY, Özyilmaz B, Zhang Y, Kim P: Energy bandgap engineering of graphene nanoribbons.
Phys Rev Lett 2007, 98(20):206805. PubMed Abstract  Publisher Full Text

Chin SK, Seah D, Lam KT, Samudra GS, Liang G: Device physics and characteristics of graphene nanoribbon tunneling FETs.

Lam KT, Seah D, Chin SK, Kumar SB, Samudra GS, Yeo YC, Liang G: A simulation study of graphenenanoribbon tunneling FET with heterojunction channel.