Appearance of Negative Differential Conductivity in Graphene Nanoribbons at High-Harmonics

We theoretically study current dynamics of graphene nanoribbons subject to bias dc and ac driven fields. We showed that graphene nanoribbons exhibit negative high-harmonic differential conductivity. Negative differential conductivity appears when bias electric filed is in the neighborhood of applied ac filed amplitude. We also observe both even and odd high-harmonic negative differential conductivity at wave mixing of two commensurate frequencies. The even harmonics are more pronounced than the odd harmonics. A possible use of the present method for generating terahertz frequencies at even harmonics in graphene is suggested.


I. INTRODUCTION
Graphene has continued to surprise scientists since its discovery in 2004 by Geim and his team 1 . Theoretically, the carrier transport properties are fantastic. Especially, its high carrier mobility of 44000cm 2 V −1 s −12 . But attempts to utilize these astonishing properties in graphene devices are posing some difficulties. The limitation is probably due to several factors including; lack of band gap in graphene sheets, edge defects, disorder, among others. To overcome some of these obstacles, particularly the lack of band gap, the dimension of graphene sheet has to be reduced. After all, new physics (quantization) emerge when dimensions of materials are reduced. An infinite 2D graphene could become a one dimensional plus quantization along one other direction, the consequence is the opening of a gap in the its dispersion spectrum. The resulting material from the 2D infinite sheet form a graphene nanoribbon (GNR). Depending on the nature of the ribbon edges, one gets two symmetry groups; Armchair Graphene Nanoribbon (aGNR) and Zigzag Graphene Nanoribbon (zGNR). Electron dynamics of both aGNR and zGNR have different properties, mostly due to the berry phase and pseudo spin 9 . Edge states have significant contribution to graphene properties, because in a nanometer size ribbon, massless Dirac fermions can reach the edges within a femto-second before encountering any other lattice effects, like electronelectron interaction, electron-phonon interaction, etc.
In this paper, we study the phenomenon of negative differential conductivity (NDC) in GNRs. In conventional semiconductor devices, a negative differential conductive behavior is known to offer great potential for high frequency applications as Bloch oscillators, frequency multipliers, and fast switching devices. For this reason, a) Corresponding author. Email: rabpeace10gh@gmail.com. b) Electronic mail: profsymensah@yahoo.co.uk the NDC effect has been greatly explored and discussed in several graphene nanostructures, particularly in 6 . NDC can also be observed in other graphene allotropes; carbon nanotubes (CNT) 12 . The unique energy spectrum of holes and electrons in GNRs, especially its narrow gapless nature leads to nontrivial features such as Negative Differential Conductivity (NDC) in the THz regime 6 .
In fact, we must emphasize that the phenomenon of NDC and Bloch oscillations in a material is a possible signature for THz generations in the material 7,11 , since NDC occurs in the THz spectral range. Most of the methods of producing THz frequencies are experimental and only very few analytical (without computer numerics) approaches are known. Though Green function techniques have also been employed in some cases. Motivated by the fact that a rigorous analytical approach is necessary for studying NDC in GNRs, we adopt a semiclassical method used in references 5,12 for armchair and zigzag CNTs. We predict that GNRs should also reproduce similar results as in reference 5 because both CNTs and GNRs have almost the same full tight binding (TB) nonlinear complex band structure.
The rest of this paper is organized as follows; In Ssection II, we derived the current density of aGNR and zGNR and imposed certain conditions to reduce the equations to forms appropriate for our models. The results obtained in Section III are plotted and discussed in Section IV. The paper finally concludes in Section V where some recommendations for future applications are made.

II. THE THEORY
As it is usually done in semi-classical treatment of quantum systems, we assume that the dynamics of the free π-electrons in graphene satisfies the time-dependent Boltzmann transport equation (BTE) in zero magnetic arXiv:1202.3569v2 [cond-mat.mtrl-sci] 11 Mar 2012 field. That is, We are also assuming relaxation time approximation (RTA) and a spatial uniform graphene nanoribbon. Also, the inverse of the relaxation time Γ is momentum independent. For the case of energy varying Γ, see 3 . In Eq.1, f 0 (k) and f (t, k) are the equilibrium and non-equilibrium Fermi electron distribution functions, respectively. e is the electronic charge, k is the electron wave vector and is the reduced plank constant. We consider an external applied field as a superposition of n harmonic waves polarized along one direction with the angular frequency ω. The phase difference between the (j + 1)th and jth wave being α j+1 − α j = α is arbitrary, j is an integer. E j are the field amplitudes. We require that ω 0 = α 0 = 0. In the following section we will look at current density for GNRs.
(4) λ = ±1. (+) for conduction band and (-) for valence band. l = √ 3a/2 and l = a/2, a is the lattice spacing with value 0.246nm, γ 0 ∼ 3.0eV is the overlap integral and θ is the phase perpendicular to the quasi-momentum k. The 1BZ of aGNR is bounded by kl = [−π/2, π/2] and the zGNR is kl = [0, π]. k is parallel to the edge and has translational symmetry along this direction. For aGNR, the transverse wave vector (phase) is quantized according to the rule 9 , θ s = s∆θ with ∆θ = π N +1 and s = 1, 2, . . . , N . Unlike aGNR, the nature of transverse wave vector quantization is complicated in zGNR, depending on both k and θ as ∆θ s = (πj +Λ(k, θ))/(N +1). However, for simplicity we assume Λ is constant, say π/2, so that ∆θ s = (2j+1) π 2 N +1 . Except this little subtlety for zGNRs, all that will be discussed in the following for aGNR are equally applied to the zGNR. Now, employed translational invariance of the graphene ribbon in the reciprocal space and expand in Fourier series functions f , f 0 and E along the edge having the periodicity in k. i.e, The Fourier coefficient f r is expressed as s in Eq.(8) counts the number of dimers N in GNRs.
The factor Φ r in Eq. (9) is a central point in this paper and so has to be determined. r is an integer and not equal to zero. We consider a classical limit in which energy levels could be excited due to thermal fluctuations, i.e ∆E << K B T << E C . This condition is also necessary for large enough field, so that charge carriers can escape low energy scattering 16 . The energy level spacing ∆E = πW γ 0 l/A, E C is the charging energy, K B is the Boltzmann constant, T is the lattice temperature and W is the graphene width. In what follows next, we will find the form of Φ(t). To do this Eqs.(5), (6), (7) are substituted in Eq.(1) to yield where Ω(t) = Ω 0 + el n j=1 E j e iωj t+αj is the modulation degree of anharmonicity in electron motion. The solution of eq.(10) is a straight forward one, using the boundary conditions t = 0, Φ r (0) = 1, one has Φ r (t) = Γ dte Γt+i n j=1 βj e iω j t+α j +iβ0t e Γt+i n j=1 βj e iω j t+α j +iβ0t , β j = erlE j / ω j and β 0 = Ω 0 = erlE 0 . One can introduce product notation in eq.(11) as Eqs. (12), (13) are connected with the well known Bessel functions via Jacobi-Anger expansion J m (β) is the m th order Bessel function. Using the expansions above, Eq.(13) becomes Applying the integration in the preceding equation and letting j = j = 1, 2, 3, . . . and m j , n j = ±1, ±2, ±3, . . .
where ν j = n j − m j .

B. Sheet current density
The sheet current density of graphene can be determined from the relation The sheet area A = W L, with W the width. g s , g v are the spin and valley degeneracies respectively. For the aGNR, The velocity of Dirac fermions in graphene is defined as v(k) = ∂E/ ∂k. In terms of the Fourier coefficients, giving with j 0,r = 2g s g v eγ 0 πl ∆θ n s=1 rE rs f rs , j * 0,r = −j 0,−r .
Substitute Φ r in Eq. (19) to get Note the r dependence of β j , β 0 and the summation over the index. Using the formalism by Litvinov and Manasson 14 , Eq. (20) can be put in a taylor-like expan- where is the differential dc conductivity (for ν j = 0), and is the large-signal dynamic nonlinear conductivity at ν j harmonic with drive frequency ω j .

A. Pure dc limit
To see immediately that Eq.(22) demonstrates NDC, we consider a pure dc limit where ω j → 0. The Bessel functions except the n = 0 term will vanish. The real part of the differential conductivity σ(0) = lim ωj →0 ∂j/∂E 0 becomes so that if β 0 > τ −1 , the differential conductivity is negative and NDC is manifest in GNRs. Electron dynamics may be come more complicated in the presence of high-frequency components in addition to the static electric fields. High negative differential conductivity thus may result in GNRs if an external drive force is applied.

B. Monoharmonics
If one component of an ac filed in Eq.(1) is applied, then n = 1 and Eq.(22) simplifies to after dropping the subscripts on n. Here, it not clear immediately how NDC can be seen. To observe it, we plot j versus E 0 for ωτ << 1 is shown in Fig.1 (left) for armchair ribbon and in Fig.1 (right) for zigzag ribbon,

C. Biharmonics
One can also allow the graphene nanoribbon subject to two ac fields. In that case, we let j = 1, 2 in our general formalism in Section II . This case has been studied in literature, especially in 15 for superlattices. Eq.(20) takes the form from which Eq.(23) follows. We have eliminated the time dependence by averaging over the period of the fields to find the time-independent current j. In the left hand side, we replaced j(t) = j, and in the right hand side a delta function emerges which ensures that ν 1 = − ω2 ω1 ν 2 . If ω 1 = ω 2 , then one must put α 1 = 0 and α 2 = α so that ν 1 = ν 2 . However, we shall generalized this to a case of commensurate frequencies. We exemplified the case by biharmonic having frequencies which can be periodic ω 2 = µω 1 or non-periodic, ω 2 = µω 1 with µ = 1, 2, . . .. These two cases were studied in 17,18 for semiconductor superlattices. Defining j 0,r = 2gvgseγ0 √ 3 (n+1)a n s=1 rE rs f rs , Eq.(26) assumes the form Simplifying further, we linearize with respect to one of the field amplitudes (say, E 2 ). For a week field, β 2 << 1 and J n (β) ≈ (β/2) 2 /n!, which allows us to take n 2 (n 2 − ν 2 ) = ±1 (0, ±1). We obtained with β 1,2 = erlE 1,2 / ω 1,2 . Where l = √ 3a/2 for armchair and l = a/2 for zigzag graphene nanoribbons respectively. Finally, the current density becomes which reduces to the monoharmonic case when µ = 1.
The nature of the NDC is observed for a simultaneously varying harmonic field and phase difference in a three dimensional plot shown in Fig.3

IV. DISCUSSION
In Fig.1, normalized current density j/j 0 is plotted against reduced static electric field E 0 /E cr for aGNR (left) and zGNR (right) for an applied ac field. At low fields up to E 0 (j max ), the quantum derivative of the j − E 0 characteristic yields a positive slope. A negative slope results for E 0 > E 0 (j max ). The whole of the region E 0 > E 0 (j max ) gives what is called Negative Differential Conductivity (NDC). A consequence of NDC in GNRs is a formation of electric field domains that impedes a continuous motion of electric field waves and thus blocks high frequency generation in these nanoribbons. NDC disappears quite faster in aGNR as E → ∞ as compared to zGNR which is more rubust at this limit.
The curves in Fig.2 demonstrate NDC, they are obtained at wave mixing of two commensurate frequencies, ω 2 = µω 1 . Fig.2 (left) µ = odd and Fig.2 (right) µ = even. The onset of NDC in odd-harmonics occurs around E 0 ∼ E 1 , and in even-harmonics it starts at E 0 ≤ E 1 . In both cases, as in the previous NDC graphs, ωτ << 1.
The combined effect of phase shift and ac amplitude on NDC is depicted in Fig.3. There are three peaks at low bias fields at points (E ∼ E cr , α = 0), (E ∼ E cr , α = π) and (E ∼ E cr , α = 2π). For now, it is not clear what these crests and throughs represents, they might be associated with field domains along one direction (for α = 0, 2π) and others along the opposite direction (for α = π) and vice-versa.

V. CONCLUSION
We have demonstrated that graphene nanoribbons exhibit NDC regions in its j − E 0 characteristics at low bias field when ωτ << 1. NDC is observed either in the presence of bias field alone or by superimposing ac field amplitudes on the bias field. For one ac field, NDC occurs around ωτ ∼ 0.1. When two ac fields at commensurate frequencies are applied, high-harmonic NDC emerge for both even and odd harmonics at rather very low frequencies ωτ ∼ 0.01. The even-series gives pronounced high-harmonic NDC than the odd-series. The presence of high-harmonic NDC means that it is possible for high-frequency generation in graphene nanoribbons when electric field domains are suppressed at high enough applied frequencies ωτ >> 1 and E 0 > E cr . We therefore suggest this approach for the study of terahertz generation in graphene.