Non-equilibrium Green's function treatment of phonon scattering in carbon nanotube transistors

We present the detailed treatment of dissipative quantum transport in carbon nanotube field-effect transistors (CNTFETs) using the non-equilibrium Green's function formalism. The effect of phonon scattering on the device characteristics of CNTFETs is explored using extensive numerical simulation. Both intra-valley and inter-valley scattering mediated by acoustic (AP), optical (OP), and radial breathing mode (RBM) phonons are treated. Realistic phonon dispersion calculations are performed using force-constant methods, and electron-phonon coupling is determined through microscopic theory. Specific simulation results are presented for (16,0), (19,0), and (22,0) zigzag CNTFETs that are in the experimentally useful diameter range. We find that the effect of phonon scattering on device performance has a distinct bias dependence. Up to moderate gate biases the influence of high-energy OP scattering is suppressed, and the device current is reduced due to elastic back-scattering by AP and low-energy RBM phonons. At large gate biases the current degradation is mainly due to high-energy OP scattering. The influence of both AP and high-energy OP scattering is reduced for larger diameter tubes. The effect of RBM mode, however, is nearly independent of the diameter for the tubes studied here.

large gate biases the current degradation is mainly due to high-energy OP scattering. The influence of both AP and high-energy OP scattering is reduced for larger diameter tubes.
The effect of RBM mode, however, is nearly independent of the diameter for the tubes studied here.

I. INTRODUCTION
Since the first demonstration of carbon nanotube (CNT) field-effect transistors in 1998 [1,2], there has been tremendous progress in their performance and the physical understanding [3]. Both electronic as well as optoelectronic devices based on CNTs have been realized, and the fabrication processes have been optimized. Ballistic transport in CNTs has been experimentally demonstrated for low-bias conditions at low temperatures [4,5]. High-performance CNT transistors operating close to the ballistic limit have also been reported [6,7,8]. The experimentally obtained carrier mobilities are of the orders 10 4~1 0 5 cm 2 /Vs [9,10] so exceptional device characteristics can indeed be expected.
Current transport in long metallic CNTs, however, is found to saturate at ~ 25 µA at high biases, and the saturation mechanism is attributed to phonon scattering [11]. On the other hand, for short length metallic tubes, the current is found not to saturate but to increase well beyond the above limit [12,13].
Nevertheless, carrier transport in these shorter tubes is still influenced by phonon scattering, and warrants a detailed physical understating of the scattering mechanisms due to its implications on device characteristics for both metallic as well as semiconducting CNTs.
There have been many theoretical studies on the calculation of carrier scattering rates and mobilities in CNTs using semiclassical transport simulation based on the Boltzmann equation [14,15,16,17,18,19,20]. Similarly, phonon mode calculations for CNTs are also performed with varying degrees of complexity: continuum and forceconstant models [21,22,23] to first-principles based methods [24,25,26]. The determination of electron-phonon (e-ph) coupling strength is performed using tightbinding calculations [27,28,29] as well as first-principles techniques [30]. It has been shown, however, that the influence of phonon scattering on device performance depends not only on the phonon modes and e-ph coupling, but also on the device geometry [31,32]. Therefore, in order to ascertain the impact of phonon scattering on the device performance, aforementioned calculations should be done in the context of specific device geometry. To that end, phonon scattering in CNT transistors has been treated using the semiclassical Boltzmann transport to determine its effects on device characteristics [31,33]. Semiclassical transport, however, can fail to rigorously treat important quantum mechanical effects, such as band-to-band tunneling, that have been deemed important in these devices [34,35,36]. Therefore, a device simulator based on dissipative quantum transport that rigorously treats the effects of phonon scattering will be essential for the proper assessment of CNT transistor characteristics, and to gain a deeper understanding of carrier transport at nanoscale.
The non-equilibrium Green's function (NEGF) formalism has been employed to describe dissipative quantum transport in nanoscale devices [37,38,39]. It has been used to treat the effects of phonon scattering in CNT Schottky barrier transistors (SBFETs) [40,41]. It has also been successfully used to investigate the impact of phonon scattering, and to explore interesting transport mechanisms such as phonon-assisted inelastic tunneling, in CNT metal-oxide semiconductor field-effect transistors (MOSFETs) with doped source and drain contacts (hereafter, simply referred to as CNTFETs) [32,34,35,36]. The NEGF simulation of ballistic transport in CNTFETs is reported in [42]. Here, we extend the previous work [42], and present the detailed simulation technique employed for the treatment of phonon scattering in them. Section II describes the tight-binding scheme, the self-consistent electrostatics, and the treatment of e-ph coupling for NEGF modeling of the CNTFETs. Section III summarizes the numerical procedures used for the simulation of phonon scattering in the self-consistent Born approximation. Section IV, followed by the conclusion in section V, has the detailed simulation results, and discusses the impact of phonon scattering on CNTFET characteristics. It compares the diameter dependence of the effect of phonon scattering in (16,0), (19,0), and (22,0) zigzag CNTs (i.e.: mod(n-m,3) = 1 type) that are in the experimentally useful diameter range (1.2nm ~ 1.8nm), below which the contact properties degrade, and above which the bandgap is too small for useful operation [43].

A. Treatment of Transport by NEGF
A detailed description of the NEGF modeling of ballistic transport in CNTFETs is described in [42]. Here we present a brief overview of that device model for the sake of completeness. The device Hamiltonian used in this study is based on the atomistic nearest-neighbor p z -orbital tight-binding approximation [21]. The device geometry, shown in Fig. 1(a), is a CNT MOSFET with doped source and drain regions (L SD ) and a cylindrical wrap-around metallic gate electrode over the intrinsic channel region (L ch ).
The gate oxide with thickness t OX covers the full length of the tube. We employ artificial heavily doped extension regions, L ext . They do not influence the transport in the working part of the transistor, but useful for better numerical convergence purposes when phonon scattering is present (however are not necessary for ballistic simulations). The cylindrical geometry of this device ensures symmetry in the angular direction thus drastically simplifying the mode-space treatment of electron transport [42,44]. It also permits the treatment of self-consistent electrostatics using 2D finite difference method [42]. The source and drain electrodes are treated as quasi-continuum reservoirs in thermal equilibrium and are modeled by the contact self-energy functions as in [42].
The NEGF model of the CNTFET used for transport simulations is shown in Fig.   1 where η + is an infinitesimal positive value, and I the identity matrix [37].
The self energy contains the contributions from all mechanisms of relaxation; the source and drain electrodes, and from scattering [37] (2) Note that in Eq. (2) the self-energy functions are, in general, energy dependent.
In the mode-space treatment of an (n,0) zigzag CNT, the dependence of the electronic state on the angle along the tube's circumference, ϕ , is expanded in a set of circular harmonics exp( ) imϕ with the angular quantum number, m. It spans the integer values of 1 to 2n, or, equivalently, -n+1 to n. Integer values on m outside this range would produce equivalent harmonics at the crystal lattice sites. The total Hamiltonian splits into independent matrices for subbands associated with each value of m [42], giving rise to a 1D Hamiltonian with two-site unit cell, as schematically shown in Fig. 1(c The diameter of the zigzag nanotube is [21] 3 cc t n a d π = (4) The mode-space transformation procedure of the real-space atomistic tight-binding Hamiltonian is well described in [42], and is not repeated here. The two-site unit cell, as expected, gives rise to two subbands corresponding to the conduction and the valence band. The Hamiltonian matrix for the subbands with angular quantum number m in an (n,0) zigzag CNT is then given by [42], 1 is the nearest neighbor hopping parameter, and N is the total number of carbon rings along the device. Here, the diagonal elements U j correspond to the on-site electrostatic potential along the tube surface. All electronic subbands in a CNT are four-fold degenerate: due to two spin states and the valley degeneracy of two [21]. The valley degeneracy comes from the two subbands with the same energy dispersion, but different m-values. Each subband can be represented as a cut of the graphene 2D Brillouin zone by a line with a constant momentum y k . In this paper we equate momentum with wavevector, having the dimension of inverse length. The cuts closest to the K-points of graphene correspond to lowest-energy conduction subbands as well as highest-energy valence subbands, and correspond in zigzag tubes to angular momenta m L1 = round(2n/3) and m L2 = round(4n/3).
Level broadening is defined as follows and can be shown [37] to be † are the in/out-scattering functions (see below). The same relations apply separately to each mechanism of relaxation.
For a layered structure like the carbon nanotube, the source self-energy function has all its entries zero except for the (1,1) element. That is [ and, Similarly, has only its (N,N) element non-zero and it is given by equations similar to (7) and (8) where f(E) is the Fermi distribution, and where the energy dependence of the Green's function and in/out-scattering functions is suppressed for clarity. The spectral function is [37] ( Note that the electron and hole correlation functions, , are matrices defined in the basis set of ring numbers i,j and subbands m (we will imply the last index in the rest of the paper where summation is performed over the spin and subband variables, and produces the degeneracy factor of 4 (for each non-equivalent subband). In the view of Eq. (13) one recognizes that the spectral function is proportional to the density of states which is traditionally defined [45] to include the spin summation, but is taken separately for each subband , 1 Finally, the current flow from site z j to z j+1 in the nearest-neighbor tight-binding scheme can be determined from [38,39], with the transmission coefficient, T(E), given by † Eq. (19) is the famous Landauer equation widely used in mesoscopic transport [37].
One can better understand the bandstructure of carbon nanotubes in by solving for the eigenvalues of the Hamiltonian (5) for zero external potential, and thereby obtaining [42] the energy dispersion relations, , versus the momentum along the length of the ( ) z E k tube, for each subband. For the lowest conduction and the highest valence subbands, close to the K-points the graphene band edge is approximately conic, thus and the distance to the K-point of The velocity of carriers in the band is Far enough from the band edge, the velocity tends to the constant value The 1D density of states including spin summation but only one subband (valley) can thus be expressed as or, in other terms ( )

B. Poisson's Equation
This section summarizes the implementation of self-consistent electrostatics in our simulation. The diagonal entries of the Hamiltonian in Eq. (5) contain the electrostatic potential on the tube surface, which thereby enters the NEGF calculation of charge distribution in Eqs. (14) and (15) Here, ρ(r,z) is the net charge density distribution which includes dopant density as well.
At this point, it should be noted that even though Eqs. (14) and (15) give the total carrier densities distributed throughout the whole energy range, what we really need for determining the self-consistent potential on the tube surface, , is the induced charge density ( = CNT radius). This can be determined by performing the integrals in Eqs.
CNT r (14) and (15) in a limited energy range defined with respect to the local charge neutrality energy, E N [42,46]. In a semiconducting CNT, due to the symmetry of the conduction and valence bands, E N is expected to be at the mid-gap energy. Finally, the induced charge density at site z j can be calculated from [42], where the first and second terms correspond to the induced electron and hole densities, respectively, with charge of the electron e.
Knowing the induced charge Q ind , the net charge distribution ρ(r,z) is given by, where, D N + and A N − are ionized donor and acceptor concentrations, respectively. Here, it is assumed that the induced charge and the dopants are uniformly distributed over the CNT surface. Finally, Eq. (27) is solved to determine the self-consistent electrostatic potential U j along the tube surface. The finite difference solution scheme for the 2D Poisson equation is described in [42]. The calculated potential, new j U , gives rise to a modified Hamiltonian (Eq. (5)), eventually leading to the self-consistent loop between electrostatics and quantum transport (Fig. 3).
Even though the self-consistent procedure we have just outlined appears conceptually straightforward, it has poor convergence properties. Therefore, a non-linear treatment of the Poisson solution is used in practice, as explained in [38,47], in order to expedite the electrostatic convergence. The convergence criterion used in this process is to monitor the maximum change in the potential profile between consecutive iterations, i.e.: where the tolerance value is normally taken to be 1meV.

C. Phonon Modes
The parameters of the phonons are obviously determined by the structure of the nanotube lattice. The one-dimensional mass density of an (n,0) nanotube is, where is the mass of a carbon atom. The energy of a phonon of momentum (in the The index of the phonon subband l is implicitly combined with the momentum index here. The half-amplitude of vibration for one phonon in a tube of length L is [45], For the reservoir in a thermal equilibrium at temperature T , the occupation of modes is given by the Bose-Einstein distribution As discussed earlier, the electron states in semiconducting CNTs have two-fold valley degeneracy with the lowest-energy subbands having angular quantum numbers m L1 and m L2 . Electron-phonon scattering is governed by energy and momentum conservation rules. Thus, as shown in Fig. 2(a) electrons can be scattered within the same subband (intra-valley) where they do not change their angular momentum, and, such scattering is facilitated by zone-center phonons having zero angular momentum (l = 0). As shown in where d CNT is the CNT diameter in nanometers [24,25,30].
The Hamiltonian of electron-phonon interaction in a general form is [45] ( where are the creation and annihilation operators for phonons in the mode q . The summation over momenta is generally defined via an integral over the first Brillouin zone, † , where is the number of unconfined dimensions. For carbon nanotubes and the limits of the integral are as follows from (3).
Electron-phonon (e-ph) coupling calculations have also been carried out, as described in [27], in conjunction with the dispersion calculations in order to account for the mode polarization effect on e-ph coupling value [47]. We find that only a few phonon modes that effectively couple to the electrons. As highlighted in Fig. 4(a) .
In this paper, we take the matrix elements as inputs and describe the general method of treatment of electron-phonon interaction in nanotubes for both optical and acoustic phonon modes.

D. Electron-Phonon Scattering
As derived in Appendix B, the in/out-scattering functions for electron-phonon scattering in a ring from subband to subband The imaginary part of self-energy is The real part of self-energy is manifested as a shift of energy levels and is computed by using the Hilbert transform [37] ' ( ') In this paper we neglect the real part of electron-phonon self energy in order to simplify For acoustic phonon scattering, the coupling constant is In Appendix B, we provide the justification for using only diagonal terms of the selfenergy and in/out-scattering functions. We also make the connection (in Appendix C) between the in/out-scattering functions in the coordinate space and the traditionally considered scattering rates in the momentum space.

III. NUMERICAL TREATMENT OF DISSIPATIVE TRANSPORT
Here, we summarize the overall simulation procedure used in this study.
Throughout this work we encounter many energy integrals such as Eqs. (17) and (28).
The use of a uniform energy grid becomes prohibitive when sharp features such as quantized energy states need to be accurately resolved. Therefore, an adaptive technique for energy integrations is used based on the quad.m subroutine of Matlab® programming language. The treatment of phonon scattering is performed using the self-consistent Born approximation [38,39]. In that, we need to treat the interdependence of the device Green's function, Eq. (1), and the scattering self-energy, Eq. (2), self-consistently. The treatment of OP scattering is presented first, followed by that for AP scattering.

A. Treatment of Optical Phonon Scattering
The determination of in/out-scattering self-energy functions, Eqs. (3) and (4) 6) Repeat steps 1 through 5 until convergence criterion is satisfied. We use the convergence of the induced carrier density, Eq. (28), as the criterion.
In the above calculations, there is a repetitive need for the inversion of a large matrix, Eq. (1), which can be a computationally expensive task. However, we only need a few diagonals of the eventual solution such as the main diagonal of G n/p for the calculation of scattering and carrier densities, and the upper/lower diagonals of G n for the calculation of current in Eq. (17). The determination of these specific diagonals, in the nearest-neighbor tight-binding scheme, can be performed using the efficient algorithms given in [49]. A Matlab® implementation of these algorithms can be found at [50].
Finally, it should be noted that the overall accuracy of the Born convergence procedure described above is confirmed at the end by observing the current continuity throughout the device, Eq. (17).

B. Treatment of Acoustic Phonon Scattering
Similar to the above method, AP scattering is treated using the following procedure, For the case of AP scattering we have introduced an additional convergence loop (step 5 above) since, unlike in inelastic scattering, here the self-consistent Born calculation at a given energy is decoupled from that at all other energy values. Similar to OP scattering, we use the efficient algorithms of [49] for numerical calculations, and confirm the overall accuracy of the convergence procedure by monitoring current continuity throughout the device.

IV. RESULTS AND DISCUSSION
Dissipative transport simulations are carried out as explained in the previous sections, and the results are compared to that with ballistic transport. Here, we first study the effects of phonon scattering on CNTFET characteristics using a (16,0) tube as a representative case. Then, we compare the diameter dependence using (16,0), (19,0) and (22,0) tubes, that belong to the mod(n-m,3) = 1 family. The device parameters ( Fig. 1(a) (Table I). The intra-LA mode is treated under AP scattering separately.  Fig. 6. Here, it is seen that up to moderate gate biases AP scattering causes a larger reduction in the device current compared to OP scattering. Even in this case, the current reduction seen for OP scattering is mainly due to the low-energy RBM mode [32]. At large gate biases, however, the effect of OP scattering becomes stronger, reducing the current by ~ 16% from the ballistic level at V GS = 0.7V. Previous studies have shown that the strong current degradation at larger gate biases is due to high-energy OP scattering processes becoming effective (mainly inter-LO/TA and intra-LO modes) [31,32]. Nevertheless, the importance of AP and low-energy RBM scattering should be appreciated since these might be the relevant scattering mechanisms under typical biasing conditions of a nanoscale transistor.
The relative behavior of OP and AP scattering can be understood by studying Fig.   7. It shows the energy-position resolved current spectrum, which is essentially the integrand of Eq. (17), under ballistic transport and OP scattering. In Fig. 7(a), it is seen that under ballistic conditions, carriers injected from the source reaches the drain without losing energy inside the device region. There exists a finite density of current below the conduction band edge (E C ) which is due to quantum mechanical tunneling. In the presence of OP scattering, however, it is seen that the carriers near the drain end relaxes to low energy states by emitting phonons (Fig. 7(b) [31,32]. Figure 8 shows the energy-position resolved electron density spectrum, which is essentially the integrand of Eq. (14). Examining Fig. 8(a), one can see that electrons are filled up to the respective Fermi levels in the two contact regions. In these regions, a characteristic interference pattern in the distribution function is observed due to quantum mechanical inference of positive and negative going states [42]. Quantized valence band states in the channel region are due to the longitudinal confinement in this effective potential well [42]. In the presence of OP scattering, few interesting features are observed in Fig. 8(b). The interference pattern seen in the contact regions are smeared due to the   [28,29].
Here we compare the ballisticity of tubes, defined as the ratio between current under scattering and the ballistic current (I scat /I ballist ), vs. FS η , defined in Fig. 7(b). Positive FS η corresponds to the on-state of the device at large positive gate biases, and negative FS η is for the off-state. The characteristic roll-off of ballisticity under OP scattering is seen in Fig. 9(a) [32]. In that, the roll-off is due to high-energy OP scattering becoming effective at large gate biases. The ballisticity reduction at small gate biases is due to the lowenergy RBM scattering [32]. In Fig. 9(a) it is seen that the impact of high-energy OP scattering decreases for larger diameter tubes. This can be easily understood by noting that the e-ph coupling parameter for these modes (intra-LO and inter-LO/TA) monotonically decreases with increasing diameter (Table I). On the other hand, the impact of the RBM mode at low gate biases seems to be nearly diameter independent for the tubes considered here, even though there is a similar decrease in e-ph coupling for larger diameter tubes (Table I). This behavior is due to the concomitant reduction of energy of the RBM mode at larger diameters that leads to an increased amount of scattering events, which ultimately cancels out the overall impact on device current.
Diameter dependence of AP scattering is shown in Fig. 9(b). The ballisticity for larger tubes is higher due to the corresponding reduction of the e-ph coupling parameter shown in Table I. They all show a slight increase in the ballisticity at larger gate biases due to majority of the positive going carriers occupying states well above the channel conduction band edge [32]. The backscattering rate is a maximum near the band edge due to increased 1D density of states and decays at larger energies [14,16,18]. It is seen that for all the tubes on Fig. 9, the impact of AP scattering is stronger compared to OP scattering until the high-energy modes become effective. Under typical biasing conditions for nanoscale transistor operation, FS η will be limited ( FS η ≤ 0.15eV) and the transport will be dominated by AP and low-energy RBM scattering [ ]. 51

V. CONCLUSION
In conclusion, we present here the detailed self-contained description of the NEGF method to simulate transport of carriers in carbon nanotube transistors with the account of both quantum effects and electron-phonon scattering. This capability is especially necessary, since it provides the rigorous treatment in the practically important where the first term in the expressions corresponds to emission of a phonon, and the second one -to absorption of a phonon. The electron-phonon coupling operator contains the sum over the phonon momentum that operates on the factors to the right of it (62) One needs to insert the factor of 4, for the number of rings in the period, to obtain that the electron-phonon coupling is a constant factor (44) and the expression for the in/out-scattering functions (37) and (38).Also, a very important conclusion is that the self-energy and the in/out-scattering functions can be treated as diagonal in this case. This significantly simplifies the problem and permits the use of various algorithms of solution of the matrix equations only applicable to 3-diagonal matrices, such as the recursive inversion method [38].
Second case, for elastic scattering, when one can neglect the energy of a phonon compared to characteristic energy differences. This is approximately fulfilled for acoustic phonons. For this case, the dependence on the momentum is typically ( ) By going beyond the assumption of a constant product of the coupling factor and the phonon occupation, we can determine how good the approximation of a diagonal selfenergy is. By representing it as a Taylor series (and we know that it is an even function) 1 ... q q q q K a n K a n q and examining the second term, we obtain ( ) This can be restated as: the off-diagonal terms of the self energy and the in/out-scattering functions have the order of magnitude of the variation of the product (67) over the first Brillouin zone. By doing an inverse Fourier transform of (67), we recognize the parameter as the inverse characteristic radius of electron-phonon interaction. Thus the alternative formulation of the above criterion is: the self-energy is diagonal if the corresponding interaction radius is much less than the crystal lattice size. (2) q Appendix C. Connection between self-energy, scattering rates, and mean free path.
In this section we draw the correspondence between the in/out-scattering functions and the scattering rates, which researchers typically deal with in the classical description of transport. The probability of scattering between two specific momentum states and of carriers is calculated according to Fermi's "golden rule" where the upper sign corresponds to absorption of a phonon and the lower sign -to emission of a phonon. The total scattering rate for carriers with momentum p is where summation is performed only over momentum variables but not the spin variables.
In other words, the spin state is assumed unchanged in scattering. For isotropic scattering, such as deformation potential of acoustic of optical phonons, the scattering rate (70) is equal to the momentum relaxation rate [45]. Also in this case, the momentum summation can be replaced with the help of Eq. (36) by the integral over energies this yields A general expression for the scattering rate (for one-dimensional structures) is For electron-phonon scattering, the constant in this expression is related to the constant in the in/out-scattering functions as Similarly one obtains for elastic scattering (both with emission and absorption of phonons) with a similar relation between the constant in the scattering rate and in the in/out- which does, in fact, coincide with the first term in (37).
The mean free path for carriers of certain energy is given by the product their velocity and scattering time By substituting the scattering rate (73) and the density of states (25), we obtain for the mean free path relative to scattering with emission of a phonon as This expression simplifies in the limit of high enough energies, i.e., far from the band edge, according to (26). We also take the limit of phonon occupation number . (82) The above nominal mean free path (81) is the upper limit over all energies. In semiconducting nanotubes, velocity is smaller for energies closer to the band edge, and the density of states is larger. Therefore the specific mean free path is shorter for energies closer the band edge, and likewise, the mean free path averaged over the carriers' distribution can be orders of magnitude shorter than (81). Therefore scattering can be significant in a 20nm-channel transistor even if the nominal mean free path is close to 1 micrometer. However the value for the nominal mean free path is sometimes used as a parameter in experiments. Note that it would provide a good estimate for the mean free path in metallic nanotubes, which have zero band gap and linear energy dispersion. Table Captions   c) e-ph coupling for acoustic phonons is determined according to Eq. (45).        η is defined as the energy difference between the source Fermi level and the channel barrier (see Fig. 8(b)).      η is defined as the energy difference between the source Fermi level and the channel barrier (see Fig. 8(b)).