Suppressing Klein tunneling in graphene using a one-dimensional array of localized scatterers

Graphene's unique physical and chemical properties make it an attractive platform for use in micro- and nanoelectronic devices. However, electrostatically controlling the flow of electrons in graphene can be challenging as a result of Klein tunneling, where electrons normally incident to a one-dimensional potential barrier of height V are perfectly transmitted even as V → ∞. In this study, theoretical and numerical calculations predict that the transmission probability for an electron wave normally incident to a one-dimensional array of localized scatterers can be significantly less than unity when the electron wavelength is smaller than the spacing between scatterers. In effect, placing periodic openings throughout a potential barrier can, somewhat counterintuitively, decrease transmission in graphene. Our results suggest that electrostatic potentials with spatial variations on the order of the electron wavelength can suppress Klein tunneling and could find applications in developing graphene electronic devices.

In the following, the basic theory for intravalley multiple scattering in graphene is presented and applied to the problem of a plane wave scattering from an infinite, one-dimensional array of localized scatterers. First, the basic formalism for calculating the scattering solutions for a plane wave scattered from a finite number of localized scatterers in graphene is presented. Next, the scattering wave functions, ψ ± K ( r), for the single and two scatterer cases are explicitly presented. The theory is then extended to the case of scattering from an infinite, one-dimensional array of localized scatterers in graphene where explicit expressions for the reflection and transmission coefficients are provided. Finally, the analogous theory for scattering from an infinite, one-dimensional array of localized scatterers in a two-dimensional electron gas (2DEG) is derived for comparison to the graphene results.

I. INTRAVALLEY MULTIPLE SCATTERING FORMALISM
In the following, we briefly review the basic formalism for intravalley multiple scattering in graphene [1], where the scattering solutions are expanded about either the + K or − K Dirac points, where K = 4π √ 3 9b x and b = 1.42Å is the carbon-carbon bond length in graphene. Let φ ± K inc ( r) be an incident Dirac plane wave spinor of energy E =hν F k 1 = hν F λ normalized to unit flux along the x- with a wave vector given by k 1 = k 1 cos(θ k 1 ) x +k 1 sin(θ k 1 ) y = k X1 x +k Y 1 y for θ k 1 ∈ − π 2 , π 2 , and λ = 2π k 1 is the wavelength of the incident wave. In all calculations,hν F = 1.0558 × 10 −28 J-m. In the following, a scattering potential consisting of a linear arrangement of 2N + 1 identical scatterers along the ydirection is considered with the scatterers' positions indexed by the integer n ∈ [−N, N], r n = nd y where d is the spacing between adjacent scatterers [ Fig. 1].
The total wave function for a Dirac plane wave spinor incident to the linear array of 2N + 1 scatterers, ψ ± K ( r), can be written as: where s l is the scattering amplitude of the l th partial wave, l max + 1 are the number of partial waves that are included in calculation of ψ ± K ( r) in Eq. (1) with l max ≥ 0, T l,± K ψ ± K ( r n ) ≡ T l,± K ψ ± K ( r) r= r n , and where H (1) l (z) is a hankel function of order l, ρ n = | r − r n |, e ±iθ n = ( r− r n )·( x±i y) ρ n , and T l,± K is the l-partial wave t−matrix operator given by: and T 0,± K = 1 is the 2 × 2 identity matrix. In writing Eq. (1), the scattering amplitudes were taken to satisfy the relationship s l = s −(l+1) , which is a consequence of the scattering potential being identical over both inequivalent lattice sites in graphene. For the numerical calculations performed in this work, each scatterer was modeled as a cylindrically symmetric step potential with an effective radius of r s , i.e., the potential for the n th scatterer was given by V 0 Θ( r − r n ) where Θ( r − r n ) = 1 if | r − r n | ≤ r s and Θ( r − r n ) = 0 for | r − r n | > r s . To consider only intravalley scattering and to neglect intervalley scattering, r s b = 1.42Å. For an individual scatterer, the l th partial wave scattering amplitude is given [1,2] by: where k 2 = E−V 0 hν F , and J l (z) is a bessel function of the first kind of order l, respectively. In all simulations, l max was chosen to take into account 99.9% of the total scattering amplitude for an individual scatterer, i.e., ∑ l max l=0 |s l | 2 ≈ 0.999 ∑ ∞ l=0 |s l | 2 . Knowledge of T l,± K ψ ± K ( r n ) for all n and l ∈ [0, l max ] completely determines ψ ± K ( r) in Eq. (1); these can be determined self-consistently as follows: and r n j = | r n − r j | and e ±iθ n j = ( r n − r j )·( x±i y) r n j . In this case, Eq. (5) formally represents a set of 2(2N + 1) × (l max + 1) equations that can be solved self-consistently as follows [1]: Define the following 2(l max + 1) × 1 column vectors related to the total and incident wave functions evaluated at scatterer j ∈ [−N, −N + 1, . . . , N − 1, N]: and the following 2(2N + 1)(l max + 1) × 1 column vectors: Using Eqs. (7) and (8), Eq. (5) for n = −N to n = N and for l = 0 to l = l max can be written compactly as: where 1 is the 2(2N + 1)(l max + 1) × 2(2N + 1)(l max + 1) identity matrix, and TT is a 2(2N + 1)(l max + 1) × 2(2N + 1)(l max + 1) matrix given by: T lmax,± K G 0,± K ( r n , r j , E) T lmax,± K G 1,± K ( r n , r j , E) T lmax,± K G 2,± K ( r n , r j E) . . . T lmax,± K G lmax,± K ( r n , r j , E) and Using Eq. (8) and Eq. (10), ψ ± K ( r) in Eq. (1) can be written compactly as:

II. SCATTERING FROM A SINGLE SCATTERER
We will first consider the problem of scattering from a single scatter located at r 0 = 0. In this case, Eq. (5) becomes T l ,± K ψ ± K (0) = T l ,± K φ ± K inc (0), which gives: Inserting Eq. (15) into Eq. (1) gives: where r = | r − r 0 |. For k 1 r 1, Eq. (16) can be approximated by: As k 1 r → ∞, the scattered wave function asymptotically decays to zero as r − 1 2 leaving only the initial plane wave,

III. SCATTERING FROM TWO IDENTICAL SCATTERERS
Consider two identical scatterers, one located at r 0 = 0 and the other at r 1 = d y. In this case, Eq. (1) becomes: 4ihν F k 1 s l G l,± K ( r, r 0 , E) T l,± K ψ ± K ( r 0 ) + G l,± K ( r, r 1 , E) T l,± K ψ ± K ( r 1 ) where Formally, inverting 1 − TT requires finding the roots to a polynomial of order 4(l max + 1), which must be numerically solved when l max > 0. As in the single scattering case, the scattered wave function in Eq. (18) scales as r − 1 2 for k 1 ρ 1 ≈ k 1 r 1 and k 1 ρ 0 ≈ k 1 r 1, which gives ψ ± K ( r) ≈ φ ± K inc ( r) for k 1 r 1.

IV. SCATTERING FROM AN INFINITE ONE-DIMENSIONAL ARRAY OF SCATTERERS
In the following, the theory for the scattering of plane waves from an infinite, periodic onedimensional array of identical scatterers (N → ∞) with lattice constant d [ Figure 1] is presented.
From translational symmetry, the total wave function satisfies the relation that ψ ± K ( r + nd y) = ψ ± K ( r)e ik Y 1 nd , which means that only ψ ± K ( r) between − d 2 ≤ y ≤ d 2 needs to be calculated. To calculate ψ ± K ( r) in Eq. (1), one can use the fact that translational symmetry implies that T l,± K ψ ± K ( r n ) = e ik Y 1 nd T l,± K ψ ± K ( r 0 ), in which case the total wave function in Eq. (1) can be written as: The various T l,± K ψ ±K ( r 0 ) in Eq. (20) are determined self-consistently by: s l T l ,± K G l,± K ( r 0 , r n , E) e ink Y 1 d T l,± K ψ ± K ( r 0 ) (21) Eq. (21) gives a total of 2(l max + 1) sets of equations, which can be written as: where 1 is a 2(l max + 1) × 2(l max + 1) identity matrix, T ψ ± K ( r 0 ) and T φ ± K inc ( r 0 ) are 2(l max + 1) × 1 column vectors given by: . . .
and TG is a 2(l max + 1) × 2(l max + 1) matrix given by: where e ±iθ 0n = ∓i n |n| and The lattice sum, S l (k 1 d, k Y 1 d) in Eq. (26), can be efficiently calculated using a small modification to a previously published method [3] as: In this work, the expression for S l (k 1 d, k Y 1 d) in Eq. (27) was numerically integrated using MATLAB [4], where the upper limit of the integrals in Eq. (27) was chosen to be a = 4000.
Additionally, the individual sums that comprise G l,±K ( r, nd y, E) in Eq. (20) can be calculated [5] using the following plane wave expansion that is valid for x = 0: , where {z} + corresponds to the smallest integer greater than z, and {z} − corresponds to the largest integer less than z.
In this case, ψ ± K ( r) in Eq. (20) can be written as [for x = 0]: In Eq. (29), ψ ± K ( r) consists of a series of plane waves for n ∈ N that are either transmitted [x > 0] or reflected [x < 0] from the scattering array along with evanescent waves along the x-direction that are freely propagating along the y-direction for n ∈ N . These types of evanescent waves have been predicted to exist in graphene for one-dimensional quantum wells [6] for certain values of k Y 1 , quantum well potential and width. For the scattering array, however, it is the periodicity of the one-dimensional array of scatterers that generates the evanescent waves in ψ ± K ( r).
Finally, it should be noted that calculating ψ ± K ( r) using a plane wave expansion in Eq. (29) is computationally efficient for |x| ≥ d since the number of plane and evanescent waves that signifi- However, the number of evanescent waves that contribute significantly to ψ ± K ( r) increases dramatically as |x| → 0, which renders the plane wave expansion in Eq. (29) as an inefficient method to calculate ψ ± K ( r). In this case, the sum of hankel functions in Eq. (28) should be explicitly evaluated. Note also that although the expression for ψ ± K (x x + y y) in Eq. (29) is valid only for x = 0, the wave function is continuous at Y 1 y for n ∈ N , ψ ± K T ( r), which can be written as: where the sum in Eq. (30) is over all open channels, n ∈ N , with the transmission coefficient for the n th open channel given by: where δ i j is the Kronecker delta (δ i j = 0 for i = j and δ i j = 1 for i = j). The total transmission probability is given by T tot = ∑ n∈N |T n | 2 .

VI. SCATTERED WAVE FUNCTION FROM AN INFINITE ONE-DIMENSIONAL ARRAY OF SCATTERERS IN A TWO-DIMENSIONAL ELECTRON GAS (2DEG)
In the following, we generalize previous work [7] on scattering of plane waves from a onedimensional periodic grating in a 2DEG to include higher partial waves [l > 0] for comparison with the results derived for graphene. Consider an incident wave with effective mass m and energy E = h 2 k 2 1 2m and normalized to unit flux along the x−direction, φ ac inc ( r) = m hk X1 e −i k 1 · r , where k 1 = k X1 x + k Y 1 y = k 1 cos(θ k 1 ) x + k 1 sin(θ k 1 ) y is the wave vector.
In this case, the total wave function, ψ ac ( r), is given by [for x = 0]: where again k (n) and k (n) X1 is real. Furthermore, the various L |l| ± [ψ ac ( r 0 )] are determined from solving the following equation: where 1 is a (2l max + 1) × (2l max + 1) identity matrix, TG ac is a (2l max + 1) × (2l max + 1) matrix with the (m, n) element given by TG ac m,n = s ac |l max +1−n| S m−n (k 1 d, k Y 1 d). In this case, TG ac can be explicitly written as: where S l (k 1 d, k Y 1 d) can be evaluated using Eq. (27), and: For a scattering potential given by V ac ( r, E) = ∑ ∞ n=−∞ V ac (E)Θ(| r − r n |) where Θ(z) = 1 for z ≤ r s and Θ(z) = 0 for z > r s , the scattering amplitudes are given by s ac l = J l (k 1 r s )J l+1 (k 2 r s ) − J l (k 2 r s )J l+1 (k 1 r s ) l+1 (k 1 r s )J l (k 2 r s ) − k 2 J l+1 (k 2 r s )H (1) l (k 1 r s ) where k 1 = 2mĒ h 2 , and k 2 = 2m(E−V ac (E)) Note that s ac l = s ac −l in Eq. (38). For comparing the calculations in a 2DEG to the calculations in graphene, the magnitudes of k 1 and k 2 for the 2DEG were chosen to be the same as in the graphene [details are given in the caption of Figure 3 in the main manuscript].
In this case, the transmitted and reflected wave functions, ψ ac T ( r) for x > 0 and ψ ac R ( r) for x < 0, are determined from Eq.
where the reflection and transmission coefficients satisfy the unitarity condition, ∑ n∈N |T ac n | 2 + |R ac n | 2 = 1, and are given by: