Mode space approach for tight-binding transport simulations in graphene nanoribbon field-effect transistors including phonon scattering
Roberto Grassi, Antonio Gnudi, Ilaria Imperiale, Elena Gnani, Susanna Reggiani, Giorgio Baccarani
aa r X i v : . [ c ond - m a t . m e s - h a ll ] A p r Mode space approach for tight-binding transport simulations in graphene nanoribbonfield-effect transistors including phonon scattering
R. Grassi, ∗ A. Gnudi, I. Imperiale, E. Gnani, S. Reggiani, and G. Baccarani
ARCES - DEI, University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy (Dated: October 9, 2018)In this paper, we present a mode space method for atomistic non-equilibrium Green’s functionsimulations of armchair graphene nanoribbon FETs that includes electron-phonon scattering. Withreference to both conventional and tunnel FET structures, we show that, in the ideal case of asmooth electrostatic potential, the modes can be decoupled in different groups without any loss ofaccuracy. Thus, inter-subband scattering due to electron-phonon interactions is properly accountedfor, while the overall simulation time considerably improves with respect to real-space, with aspeed-up factor of 40 for a 1 . I. INTRODUCTION
Graphene nanoribbons (GNRs) have been proposed inrecent years as a possible replacement for silicon in thefuture generation of field effect transistors [1, 2]. Despitethe outstanding challenges in producing GNRs with con-trolled width and edges and the resulting degradation ofthe graphene intrinsic mobility, devices with large on-offcurrent ratio have been successfully demonstrated [3].Theoretical performance of GNRFETs have beenwidely investigated by numerous simulation studies (e.gRefs. 4–9). The state-of-the-art approach for modelingthe electronic properties of GNRs is based on an atom-istic tight-binding (TB) Hamiltonian with a p z orbitalbasis set. The transport problem is usually solved withinthe nonequilibrium Green’s function formalism (NEGF)[10], which provides a rigorous framework for includingincoherent scattering processes in the quantum descrip-tion. TB simulations of GNRFETs including phononscattering have been reported too [11–14]. Those simula-tions were performed using a real space (RS) approach.On the other hand, more efficient mode space (MS) meth-ods would be preferable for use in intensive device simu-lations.The MS approach is well established with reference toan effective mass (EM) Hamiltonian [15, 16]. It is basedon the expansion of the Green’s functions in terms ofthe transverse eigenfunctions (modes) and it is most ef-ficient when quantum confinement is relatively strong inthe transverse plane, so that only few of the lowest sub-bands are occupied. The so-called coupled mode space(CMS) is the general method, while a more efficient ver-sion, the uncoupled mode space (UMS), can be adopted ∗ [email protected] when the transverse potential profile is slightly varyingalong the transport direction. The inclusion of electron-phonon scattering in the MS EM model is also well known[17–21].On the other hand, the MS approach is not gener-ally applicable to a TB Hamiltonian. Due to the non-separable Hamiltonian, even when the potential is uni-form along the longitudinal direction (so that the longi-tudinal wavevector k is conserved), the wavefunctions atdifferent k in the same subband are in general different.In other words, by defining the modes as the transversepart of the wavefunctions at a particular k , the wave-functions at the generic k are a linear combination ofthem. This is similar to what occurs in the case of a k · p Hamiltonian [22]. Nevertheless, for the particular case ofthe TB Hamiltonian of GNRs with armchair edges, wehave previously shown by numerical calculation that aMS method is actually possible, since each subband con-tains only few modes [23]. A formal derivation was givenin Ref. 24 using analytically defined modes.In this paper, we propose a novel MS method for theinclusion of graphene acoustic phonon (AP) and opticalphonon (OP) scattering within the NEGF formalism, andtest its accuracy with respect to RS in both the cases ofsmooth and disordered potentials. Two approximationscommonly found in the literature for treating the scatter-ing self-energy are also evaluated. The paper is organizedas follows. Sec. II reviews the standard RS TB methodas well as the MS TB method from Ref. 23, modified soas to include phonon scattering. Simulation results ofGNRFETs are presented in Sec. III and conclusions arefinally drawn in Sec. IV.
FIG. 1. One-dimensional elementary cell of an N a = 13 arm-chair GNR. II. MATHEMATICAL MODELA. Real space formulation
The atomic structure of an armchair GNR unit cell isshown for reference in Fig. 1. The GNR width is denotedby the number of dimer lines N a . We adopt the TBHamiltonian model proposed in Ref. 25, which is basedon a set of orthogonal p z orbitals, one for each atom,with hopping integrals limited to first nearest-neighbororbitals and of value t for internal atom pairs, and t (1+ δ )for atom pairs located along the edges of the GNR ( t = − . δ = 0 . G R , G < and G > , respectively, for eachenergy E (cid:2) ( E + i0 + ) I − H d − Σ R ( E ) (cid:3) G R ( E ) = I (1) G < ( E ) = G R ( E ) Σ < ( E ) G A ( E ) (2) G > ( E ) = G R ( E ) Σ > ( E ) G A ( E ) (3)where H d denotes the restriction of the Hamiltonian ma-trix H to the GNR portion inside the simulation domainand G A = G R † . Besides, Σ R ( E ) = Σ RP ( E ) + Σ RS ( E ) + Σ RD ( E ) and similarly for Σ > ( E ) and Σ < ( E ), where Σ ( ... ) S/D are the self-energies for the semi-infinite source and drainleads, which are assumed in thermodynamic equilibrium,while Σ ( ... ) P account for phonon scattering. Calculationof both (2) and (3) is actually not needed, since eitherone of the two equations can be replaced by the identity G > − G < = G R − G A .Phonon scattering is treated within the self-consistentBorn approximation. The model accounts for both OPand AP scattering, the latter within the elastic and hightemperature limit, resulting in the following expressions of the lesser and greater phonon self energies Σ
P ( E ) = (cid:2) D ap G > ( E ) + D op ( N op + 1) G > ( E − ¯ hω op ) ++ D op N op G > ( E + ¯ hω op ) (cid:3) ◦ I (5)where ◦ indicates element-by-element matrix multiplica-tion, N op is the OP occupation number, and D ap and D op are given by D ap = D ac k B T m c v s , D op = D o ¯ h m c ω op (6)with D ac = 16 eV the AP deformation potential, v s =2 × cm/s the sound velocity in graphene, T the tem-perature, m c the carbon atomic mass, D o = 10 eV/cmand ¯ hω op = 160 meV the zone-boundary OP deformationpotential and energy, respectively. The retarded phononself-energy is computed using the identity Σ RP ( E ) = i2 π + ∞ Z −∞ Σ >P ( E ′ ) − Σ
P ( E ) − Σ
P ( E ′ ) − Σ
In this section the MS approach presented in Ref. 23,and in that paper discussed with reference to coherenttransport simulations of GNRs, is extended to includephonon scattering.The MS formulation starts by defining a set of or-thonormal vectors, or modes, φ ( i ) for each slab i of thedevice. We choose the modes of the generic slab as theeigenvectors computed at k = 0 of a fictitious GNR ob-tained by the periodic repetition (potential energy in-cluded) of that slab along the longitudinal direction. Inthe absence of topological differences between the slabs(such as edge irregularities or internal vacancies), themodes are the solution of the eigenvalue problem (cid:16) H i,i + H i,i +1 + H † i,i +1 (cid:17) φ m ( i ) = ε m φ m ( i ) (8)where H i,j is the Hamiltonian block between slab i and j . Denoting by v ( i ) the matrix whose columns are themodes of layer i , i.e. v ( i ) = [ φ ( i ) · · · φ m ( i ) · · · ], a block-diagonal transformation matrix V satisfying V † V = I is constructed V = v (1) 0 · · · v (2) . . . ...... . . . . . . . . .0 · · · . . . v ( N s ) (9)with N s equal to the number of slabs.In the following, we indicate with a tilde the MS matri-ces to distinguish them from the RS matrices. The CMSmethod consists in approximating G R (and similarly for G < and G > ) as G R ( E ) ≃ V e G R ( E ) V † (10)where e G R is the solution of h ( E + i0 + ) I − f H d − e Σ R ( E ) i e G R ( E ) = I (11)with f H d = V † H d V (12) e Σ R ( E ) = V † Σ R ( E ) V (13)Eq. (11) is the MS version of (1). Analogous consider-ations apply to the other equations. Eq. (10) would beexact if V was a square matrix. In practice, a mode trun-cation is performed so that f H d and e Σ R have a smallersize than their corresponding RS matrices and the solu-tion of (11) instead of (1) is computationally advanta-geous. The computational time of the RGF algorithmscales as O ( N y N s ), where N y is the matrix block size[27], which is equal to the slab size for RS (2 N a accord-ing to our choice of the slab size) and to the number N m of selected modes per slab for CMS. Thus, the speed-up of CMS compared to RS is a factor of the order of(2 N a /N m ) .An algorithm to select the modes to retain in the V matrix was presented in Ref. 23. It is based on identify-ing, among the modes calculated with zero electrostaticpotential, the minimum set of modes that allow to repro-duce with sufficient accuracy the subbands that lie in theenergy range of interest. The mode indexes so identifiedare then used to select the actual modes calculated withthe non-null electrostatic potential. As shown in Fig. 2for a N a = 13 GNR, it turns out that each one of thelowest conduction subbands at zero potential can be well k ( π /3a cc ) -6-4-20246 E ( e V ) RSMS group 1 0 0.2 0.4 0.6 0.8 1RSMS group 2
FIG. 2. Subband structure of a N a = 13 GNR computedwith RS (dashed lines) and MS (solid lines) using two differ-ent groups of modes (left and right). The modes included ineach group are the ones that correspond to the eigenvalues at k = 0 indicated with circles in each figure. The electrostaticpotential is set to zero. 3 a cc is the length of the GNR unitcell, with a cc the carbon interatomic distance. reproduced by just four modes, which correspond to theeigenvalues at k = 0 belonging to ( i ) the considered sub-band, ( ii ) the valence subband symmetrical to it, and( iii-iv ) their respective folded continuations. Symmetri-cal considerations apply to the highest valence subbands.GNRs with different N a behave similarly. Since the setsof four modes described above, hereafter referred to asgroups of modes, are disjoint from each other, a moreefficient MS method, here called uncoupled group modespace (UGMS), can be obtained from (10)-(13) by ne-glecting the coupling between modes belonging to dif-ferent groups. Let v b ( i ) be the matrix whose columnsare the modes in layer i of group b only. A transforma-tion matrix V b similar to (9) is constructed by using thematrices v b as diagonal blocks. The UGMS method isderived by further approximating (10) as G R ( E ) ≃ N g X b =1 V b e G Rb ( E ) V b † (14)where N g is the number of the considered groups and e G Rb is the solution of h ( E + i0 + ) I − f H db − e Σ Rb ( E ) i e G Rb ( E ) = I (15)with f H db = V b † H d V b (16) e Σ Rb ( E ) = V b † Σ R ( E ) V b (17)In the UGMS method, the computational cost of theRGF algorithm scales as O ( N g N y N s ), where N y = 4is the number of modes in each group. Hence, if thesame number of modes N m = 4 N g is used, the speed-upof UGMS compared to CMS scales as N g , which can besizeable for wide ribbons, since N g is proportional to thenumber of subbands that contribute to transport, which,in turn, is proportional to the GNR width.In the following, we refer to the UGMS formulation,since the CMS one can be recovered from Eqs. (14)-(17)by considering all the selected modes as belonging to thesame group and by setting N g = 1. Note that Σ RS and Σ RD can be directly computed in MS without making useof (17). Instead, the calculation of e Σ ( ... ) P is more com-plicated. In particular, in the MS representation thesematrices are no longer diagonal, but just block diagonal.For example, as far as e Σ
All simulated devices have a double-gate structure withSiO N a = 13 (corresponding to a ribbonwidth W ≃ . L G = 17 nm, and dop-ing concentration in the source and drain regions equalto 10 − dopants per carbon atom. Also, unless statedotherwise, the MS method simulations are performed us-ing the uncoupled group approximation (UGMS). For thereference N a = 13 GNR, the 2 groups of 4 modes of Fig. 2are used (8 modes out of a total of 26 that correspond tothe full RS solution). Both conventional FET and tunnelFET (TFET) devices are simulated. The latter have re-ceived great attention in recent years for their potentialin low-power applications [8]. The two structures onlydiffer in the type of doping of the source region. Whilethe drain is always n-type, the source is n-type and p-type for conventional FETs and TFETs, respectively.We consider only ideal GNRs with perfect edges andno internal vacancies. Nevertheless, in the following, wealso study the effect of adding a disorder potential tothe TB Hamiltonian of the ideal GNR. Such disorderpotential can, in principle, mimic the perturbation effectof impurities or the substrate. A. Smooth electrostatic potential
The case of no disorder potential is treated first.The I vs. V GS (“turn-on”) characteristics of a con-ventional n-i-n FET at V DS = 0 . i ) without scattering (i.e. ballistic trans-port), ( ii ) in the presence of only AP scattering, and( iii ) in the presence of both AP and OP scattering. Inthe case of phonon scattering and MS approach, the sim-plified expression of the form-factor in (20) is used (MS1).As regards the effect of phonon scattering and with ref-erence to the RS results, it can be seen that AP scatteringhas only a limited effect on the current at this channellength, resulting in a ballisticity ratio (i.e. ratio betweencurrent in the presence of phonon scattering and ballis-tic current) of 0 . V GS = 0 . V DS values.When also OP scattering is included, the current at high V GS , i.e. on-state current, is only slightly decreased forthe largest V DS value. Similar findings were reported inRef. 12. On the other hand, OP scattering is responsi-ble for an increase of the minimum off-state current bya few orders of magnitudes when V DS is low. This effectis caused by energy relaxation through emission or ab-sorption of optical phonons, which favors band-to-band -9 -6 -3 I ( µ A ) ballisticAPAP+OPN a = 13, L G = 17 nm 036912 I ( µ A ) -0.2 0 0.2 0.4 0.6 0.8 V GS (V) -9 -6 -3 I ( µ A ) I ( µ A ) lines: UGMS1symbols: RS V DS = 0.1 VV DS = 0.8 V FIG. 3. Turn-on characteristics of a n-i-n FET with W =1 . V DS = 0 . V DS = 0 . x (nm) C u rr en t ( µ A ) I (x)I (x)I (x)+I (x) FIG. 4. Current spectrum J b ( x, E ) obtained by UGMS1method for mode group b = 1 (left) and b = 2 (right) for thedevice in Fig. 3 with AP scattering and at the bias indicatedin the figure. Blue (red) color means low (high) density. Thewhite solid (dashed) line is the profile of the first (second) con-duction subband. The source Fermi level is at E = 0. Theinset shows the integrals over energy I b ( x ) = R J b ( x, E )d E and the conservation of the total current I ( x ) + I ( x ). tunneling (BTBT) and shifts toward positive V GS valuesthe onset of the ambipolar conduction by BTBT, similarto what occurs in carbon nanotubes [28].With regard to the MS results, we do not observe anysignificant discrepancy with respect to RS, except for thecurve with OP scattering at V DS = 0 . V GS point of minimum current. This discrepancy will bediscussed later in the text.In order to test the validity of the MS method in de- x (nm) C u rr en t ( p A ) I (x)I (x)I (x)+I (x) FIG. 5. Same as in Fig. 4 but with both AP and OP scatteringand at V GS = − . V DS = 0 . scribing inter-subband scattering, we separately plot inFig. 4 the MS current spectrum for the first and secondgroup of modes (corresponding to the first and secondsubband, respectively) at V GS = 0 . V DS = 0 . I vs. V GS characteristics at V DS = 0 . V GS (V) -9 -6 -3 I ( µ A ) ballisticAPAP+OP 00.40.81.2 I ( µ A ) V GS (V) -6 -4 -2 I ( µ A ) I ( µ A ) solid: UGMS1symbols: RSdashed: UGMS2 N a = 13N a = 40 V DS = 0.1 VV DS = 0.4 VL G = 17 nmL G = 30 nm FIG. 6. Turn-on characteristics of p-i-n FETs. Top: devicewith W = 1 . V DS = 0 . W = 5 nm at V DS = 0 . W = 1 . J(L S ,E) (10 -6 µ A/eV) -0.8-0.6-0.4-0.200.20.4 E ( e V ) UGMS1RSUGMS2AP+OPV GS = 0.2 VV DS = 0.4 V FIG. 7. Left: current spectrum J ( x, E ) obtained by RSmethod for the device in Fig. 6–top with both AP and OPscattering and at the bias indicated in the figure. The verticalline indicates the position L S of the source-channel junction.The two white lines are the profile of the first conduction andvalence subbands. Right: comparison between the currentspectra at x = L S obtained by RS, UGMS1, and UGMS2methods. channel and drain-channel junctions: although the firstconduction and valence subbands do not face each other,electrons can transmit from the valence to the conduc-tion subband by absorption of optical phonons. Theseresults are in agreement with the ones in Ref. 13.Again, the only difference (of about a factor of 1 . -9 -6 -3 I ( µ A ) ballisticAP, ℜ { Σ r }=0AP, ℜ { Σ r } ≠ a = 13, L G = 17 nm, V DS = 0.4 V 00.40.81.2 I ( µ A ) V GS (V) -9 -6 -3 I ( µ A ) ballisticOP, ℜ { Σ r } = ℜ { Σ r } ≠ I ( µ A ) solid: UGMS1symbols: RSUGMS2 FIG. 8. Turn-on characteristics for the same device as inFig. 6–top comparing different approximations for Σ RP . Top:only AP scattering is included, with or without the Hermitianpart of Σ RP . Bottom: only OP scattering is included, with orwithout the Hermitian part of Σ RP . The ballistic curve isshown for reference in both figures and the solution methodfor each curve is indicated in the legend. ence of OP scattering at the point of minimum current.However, the accuracy with respect to RS can be almostcompletely recovered by the MS2 method, which usesthe exact expression of the form-factor in (19). This ismore clearly shown by the comparison in Fig. 7–right be-tween the current spectra at the source-channel junctionobtained with the different methods. The lack of accu-racy of (20) in this bias condition could be related to theneglect of the terms F b,mnb,mn ( i ) with n = m , which are ac-tually of the same size as the terms F b,nnb,mm ( i ) includedin (20). As a drawback, a slow-down of the simulationby about a factor of 1 . V DS is chosen due to the lowerband gap. In addition, we set N a = 40 (correspondingto W ≃ L G = 30 nm, and a doping concentrationin the source and drain regions of 7 · − dopants/atom.4 groups of 4 modes are used in the UGMS simulations.It is worth noticing that the UGMS method is still accu-rate for the wider ribbon, indicating that the decouplingin separate groups is still valid, even though the subbandsare more closely spaced than in the N a = 13 GNR.To evaluate the importance of the Hermitian part of Σ RP , here denoted by ℜ{ Σ RP } , we report in Fig. 8 the I vs. V GS characteristics of the N a = 13 TFET computed withor without ℜ{ Σ RP } , separately for each type of scattering.For AP scattering (Fig. 8–top), the neglect of ℜ{ Σ RP } leads to an underestimation of the on-state current. Thiscan be understood as follows. In general, ℜ{ Σ RP } has the J(E) (10 -9 µ A/eV) -0.2-0.100.10.20.3 E ( e V ) -6 -4 -2 LDOS(x,E) (eV -1 ) E ( e V ) ℜ { Σ r }=0AP, ℜ { Σ r } ≠ GS = 0.2 VV DS = 0.4 V x=0 x=L M x = L M x = 0 FIG. 9. Current spectrum (left) and integral over the slabof the LDOS at x = L M (center) and x = 0 (right) for thedevice in Fig. 8–top at V GS = 0 . V DS = 0 . L M isthe mid-channel position. RS method is used. The ballisticsubband profile is shown in the inset. effect of shifting the Hamiltonian eigenvalues [29]. How-ever, in GNRs, the shift is of opposite sign for energiesabove and below the GNR mid-gap due to the symme-try of the subband structure. The result is a decreaseof the GNR band gap, which favors BTBT, so that alarger current is expected when ℜ{ Σ RP } is included in thesimulation. Interestingly, the two expressions of ℜ{ Σ RP } give almost identical values with respect to the minimumcurrent. To clarify this point, we plot in Fig. 9 the cur-rent spectrum and the local density of states (LDOS)per slab at two positions along the device, for the biaspoint corresponding to the minimum current. First, itcan be seen that the current spectrum with AP scatter-ing and ℜ{ Σ RP } 6 = 0 is larger than the ballistic one: thereason can be ascribed to an enhanced BTBT throughthe channel region due to the band gap narrowing effectmentioned above, which can be appreciated from the log-arithmic plot of the LDOS at the mid-channel positionin Fig. 9–center. Secondly, by looking at Fig. 9–right, itcan be noticed that the peaks of the LDOS in the sourceregion with ℜ{ Σ RP } 6 = 0 are located at the same energypositions as the ones of the ballistic LDOS. On the con-trary, the peaks of the LDOS with ℜ{ Σ RP } = 0 are shiftedup in energy, resulting in a tunneling current larger thanthe ballistic one and similar to the one with ℜ{ Σ RP } 6 = 0(Fig. 9–left). The shift of the LDOS can be attributedto a “loss of charge” when ℜ{ Σ RP } is set to zero [30] andto the combined effect of the electrostatic feedback.For OP scattering instead, no relevant difference is ob-served when including ℜ{ Σ RP } , even at high V GS (Fig. 8–bottom). J(E) ( µ A/eV) -0.3-0.2-0.100.10.20.3 E ( e V ) RSUGMSCMS 0 10 20 30 x (nm) -4-3-2-10123 ( n - p ) Ω ( % ) δ V = 0.5 eV V DS = 0.4 VV GS = 0.6 V FIG. 10. Integral over the slab of the LDOS (left), currentspectrum (center), and average over the slab of the net elec-tron concentration per carbon atom (right) for the device inFig. 6–top at V GS = 0 . V DS = 0 . a cc / B. Disordered electrostatic potential
We focus on the N a = 13 TFET. The simulations areperformed in a non-self-consistent way, by solving theNEGF equations with a fixed electrostatic potential. Wetake the electrostatic potential as the sum of a disorderpotential and the one calculated self-consistently in theabsence of disorder and in the ballistic limit at V GS =0 . V DS = 0 . V i at the atomicsite i , located at position ~r i , is calculated as V i = N X j =1 S j δV exp − (cid:12)(cid:12)(cid:12) ~r i − ~X j (cid:12)(cid:12)(cid:12) l (21)where N , l , δV are parameters, while S j and ~X j are ran-dom variables: S j can take values ± ~X j is uniformly distributed over all the atomicposition ~r i . The second model we study is the Ander-son type of disorder [32], according to which V i = Y i ,where Y i is a random variable uniformly distributed in[ − δV / , δV / N = 0 . N c , l = 5 a cc , and δV = 0 . N c is the total number of carbonatoms inside the device and a cc is the carbon-carbonbond length. The resonant states induced by disorderare clearly visible in the LDOS. In Fig. 10–center and J(E) ( µ A/eV) -0.3-0.2-0.100.10.20.3 E ( e V ) RSUGMSCMS 0 10 20 30 x (nm) -2-1012 ( n - p ) Ω ( % ) δ V = 0.5 eV V DS = 0.4 VV GS = 0.6 V FIG. 11. Same as in Fig. 10 but in the presence of a short-range disorder potential.
Fig. 10–right we compare the current spectrum and netelectron charge along the device, respectively, obtainedwith RS, UGMS, and CMS. The UGMS simulations areperformed with the same 2 groups of 4 modes used inthe case without disorder, while the CMS ones use thesame 8 modes but all coupled in one group. It turns outthat the UGSM method looses accuracy due to disorder-induced mode mixing, while the CMS method leads toresults very close to the RS ones (e.g. the error on thecurrent value is less than 2%).Similar considerations can be made regarding the re-sults obtained with the second model of disorder with δV = 0 . . . µ A). On the other hand, in the case of long-rangedisorder, the current is increased from 0 .
084 to 0 . µ Awhen phonon scattering is included in the simulation,indicating that the localization transport regime [29] thatoccurs in the coherent approximation is broken by thedephasing effect of phonon scattering.
IV. CONCLUSIONS
A mode space method for TB NEGF simulations ofarmchair GNR FETs including phonon scattering hasbeen presented and tested with reference to both conven-tional and tunnel FET structures. When no disorder isincluded in the simulation, an efficient decoupling of themodes in different groups (UGMS) can be employed with
J(L S ,E) ( µ A/eV) -0.200.2 E ( e V ) RSCMS1CMS2 0 10 20 30 x (nm) -4-202 ( n - p ) Ω ( % ) J(L S ,E) ( µ A/eV) -0.200.2 E ( e V ) RSCMS1CMS2 0 10 20 30 x (nm) -2-101 ( n - p ) Ω ( % ) δ V = 0.5 eVV DS = 0.4 VV GS = 0.6 VAP+OPlong-range disordershort-range disorder FIG. 12. Top: current spectrum at the source-channel junc-tion (left) and average over the slab of the net electron con-centration per carbon atom (right) for the device in Fig. 6–topat V GS = 0 . V DS = 0 . excellent accuracy. Despite the decoupling, the methodcorrectly accounts for inter-subband scattering. Simpli-fied expressions of the scattering-self energies have beencompared. The one obtained by neglecting some entriesof the form-factor is found to be accurate except for thebias points where transport occurs by phonon-assistedBTBT. On the other hand, the real part of the scatter-ing self-energy has only a limited effect on the devicecharacteristics, especially for the case of optical phononscattering, where its calculation is most demanding. Inthe presence of a disorder potential, the modes need tobe coupled in a single group (CMS) to account for modemixing, but no additional modes, compared to the onesused to simulate the case without disorder, need to beincluded to achieve accurate results.While the computational advantage of CMS over RS isabout a constant factor with respect to the GNR width(equal to about 30 for the simulation parameters chosenin this paper), the speed-up of UGMS compared to RSis about 40 for a device width of 1.5 nm and increasesproportionally to the second power of the GNR width(for the 5-nm-wide device considered in this paper suchspeed-up is about 360). ACKNOWLEDGMENTS
R. G. would like to thank Dr. E. Baravelli of Universityof Bologna for fruitful discussions on the mode-space ap-proach. This work has been supported by the EU projectGRADE 317839. The authors acknowledge the CINECAAward N. HP10CPFJ69, 2011 for the availability of highperformance computing resources and support.
Appendix A: Calculation of the retarded phononself-energy
We consider here only the RS case, the generalizationof the expressions to MS being straightforward.Replacing (4) and (5) into the last member of (7) andusing the identity analogous to (7) valid for the G ma-trices, one can write Σ RP ( E ) = Σ RAP ( E ) + Σ ROP ( E ) (A1)with Σ RAP ( E ) = D ap I ◦ G R ( E ) (A2)and Σ ROP ( E ) = D op I ◦ (cid:26) ( N op + 1) G R ( E − ¯ hω op ) + N op G R ( E + ¯ hω op ) ++ G < ( E − ¯ hω op ) − G < ( E + ¯ hω op )2 ++iP + ∞ Z −∞ G < ( E ′ − ¯ hω op ) − G < ( E ′ + ¯ hω op )2 π ( E − E ′ ) d E ′ (cid:27) (A3)which does not contain G > . It is worth noticing that theprincipal part integral in (A3) contains only a fraction ofthe Hermitian part of Σ ROP [19]. It is calculated here bymeans of a piecewise constant approximation of G < ( E ′ )over the energy domain [ E min , E max ], which is discretizeduniformly with energy steps ∆ E typically of the order of¯ hω op / + ∞ Z −∞ G < ( E ′ ∓ ¯ hω op ) E − E ′ d E ′ ≃≃ X j G < ( E j ) ln (cid:12)(cid:12)(cid:12)(cid:12) E − E j + ∆ E / ∓ ¯ hω op E − E j − ∆ E / ∓ ¯ hω op (cid:12)(cid:12)(cid:12)(cid:12) (A4) where the summation extends over all discrete energypoints E j ∈ [ E min , E max ]. An expression of Σ ROP analo-gous to (A3) containing G > instead of G < can also bederived. Expression (A3) is preferred for E higher thanthe contact Fermi energies, since the numerical error in-troduced by truncating the upper limit of the integral to E max is minimized due to the decaying nature of G < athigh energies. On the contrary, the alternative expressionof Σ ROP which depends on G > is used for low energies be-low the contact Fermi levels, for analogous reasons. Forintermediate energies an average of the two formulationsis used.In Sec. III A, the results obtained with the full expres-sion of Σ RAP and Σ ROP are compared with approximatesolutions obtained by neglecting the respective Hermitianparts, i.e. Σ RAP ( E ) ≃ D ap I ◦ G > ( E ) − G < ( E )2 (A5)and Σ ROP ( E ) ≃ D op I ◦ (cid:26) ( N op + 1) G > ( E − ¯ hω op ) ++ N op G > ( E + ¯ hω op ) − ( N op + 1) G < ( E + ¯ hω op ) + − N op G < ( E − ¯ hω op ) (cid:27) (A6)Expression (A6) is much more efficient than (A3) due tothe absence of the integral term. [1] M. Y. Han, B. ¨Ozyilmaz, Y. Zhang, and P. Kim, Phys.Rev. Lett. , 206805 (2007).[2] Z. Chen, Y.-M. Lin, M. J. Rooks, and P. Avouris, Phys-ica E , 228 (2007).[3] X. Wang, Y. Ouyang, X. Li, H. Wang, J. Guo, andH. Dai, Phys. Rev. Lett. , 206803 (2008).[4] G. Fiori and G. Iannaccone, IEEE Electron Device Let-ters , 760 (2007).[5] Y. Ouyang, Y. Yoon, and J. Guo, IEEE Trans. ElectronDevices , 2223 (2007). [6] G. Liang, N. Neophytou, M. S. Lundstrom, and D. E.Nikonov, J. Appl. Phys. , 054307 (2007).[7] Y. Yoon, G. Fiori, S. Hong, G. Iannaccone, and J. Guo,IEEE Trans. Electron Devices , 2314 (2008).[8] P. Zhao, J. Chauhan, and J. Guo, Nano Lett. , 684(2009).[9] R. Grassi, A. Gnudi, E. Gnani, S. Reggiani, and G. Bac-carani, J. Comput. Electronics , 441 (2009).[10] S. Datta, Quantum Transport: Atom to Transistor (Cambridge University Press, Cambridge, UK, 2005). [11] Y. Ouyang, X. Wang, H. Dai, and J. Guo, Appl. Phys.Lett. , 243124 (2008).[12] Y. Yoon, D. E. Nikonov, and S. Salahuddin, Appl. Phys.Lett. , 203503 (2011).[13] Y. Yoon and S. Salahuddin, Appl. Phys. Lett. ,263501 (2012).[14] N. D. Akhavan, G. Jolley, G. A. Umana-Membreno,J. Antoszewski, and L. Faraone, J. Appl. Phys. ,094505 (2012).[15] J. Wang, E. Polizzi, and M. Lundstrom, J. Appl. Phys. , 2192 (2004).[16] M. Luisier, A. Schenk, and W. Fichtner, J. Appl. Phys. , 043713 (2006).[17] S. Jin, Y. J. Park, and H. S. Min, J. Appl. Phys. ,123719 (2006).[18] S. Poli, Modelling and simulations of postCMOS devices ,Ph.D. thesis, University of Bologna (2009).[19] A. Esposito, M. Frey, and A. Schenk, J. Comput. Elec-tron. , 336 (2009).[20] D. Nikonov, H. Pal, and G. Bouri-anoff, “Scattering in NEGF: Made simple,” http://nanohub.org/resources/7772 (2009).[21] A. Afzalian, J. Appl. Phys. 110 , 094517 (2011). [22] M. Shin, J. Appl. Phys. , 054505 (2009).[23] R. Grassi, A. Gnudi, E. Gnani, S. Reggiani, and G. Bac-carani, IEEE Transactions on Nanotechnology , 371(2011).[24] P. Zhao and J. Guo, J. Appl. Phys. , 034503 (2009).[25] Y.-W. Son, M. L. Cohen, and S. G. Louie, Phys. Rev.Lett. , 216803 (2006).[26] R. Lake, G. Klimeck, R. C. Bowen, and D. Jovanovic, J.Appl. Phys. , 7845 (1997).[27] A. Svizhenko, M. P. Anantram, T. R. Govindan,B. Biegel, and R. Venugopal, J. Appl. Phys. , 2343(2002).[28] S. O. Koswatta, M. S. Lundstrom, M. P. Anantram, andD. E. Nikonov, Appl. Phys. Lett. , 253107 (2005).[29] S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge University Press, Cambridge, UK, 1997).[30] A. Svizhenko and M. P. Anantram, Phys. Rev. B ,085430 (2005).[31] M. Poljak, E. B. Song, M. Wang, T. Suligoj, and K. L.Wang, IEEE Trans. Electron Devices , 3231 (2012).[32] A. Lherbier, B. Biel, Y.-M. Niquet, and S. Roche, Phys.Rev. Lett.100