Sinc-based method for an efficient solution in the direct space of quantum wave equations with periodic boundary conditions
aa r X i v : . [ c ond - m a t . m e s - h a ll ] D ec Sinc-based method for an efficient solution in the direct space of quantum waveequations with periodic boundary conditions
Paolo Marconcini, a) Demetrio Logoteta, and Massimo Macucci Dipartimento di Ingegneria dell’Informazione, Universit`a di Pisa, Via Caruso 16,I-56122 Pisa, Italy.
The solution of differential problems, and in particular of quantum wave equations,can in general be performed both in the direct and in the reciprocal space. However,to achieve the same accuracy, direct-space finite-difference approaches usually in-volve handling larger algebraic problems with respect to the approaches based on theFourier transform in reciprocal space. This is the result of the errors that direct-spacediscretization formulas introduce into the treatment of derivatives. Here, we proposean approach, relying on a set of sinc-based functions, that allows us to achieve anexact representation of the derivatives in the direct space and that is equivalent tothe solution in the reciprocal space. We apply this method to the numerical solutionof the Dirac equation in an armchair graphene nanoribbon with a potential varyingonly in the transverse direction.PACS numbers: 03.65.Ge, 02.60.Lj, 72.80.VpKeywords: wave equations, differential equations, sinc, Dirac equation a) Corresponding author: [email protected] . INTRODUCTION With the progressive downscaling of electron devices, whose size has reached well into thenanoscale, the quantum simulation of charge transport has become increasingly important.This study requires the numerical solution of the quantum wave equation for the chargecarriers in the device.In the case of periodic boundary conditions, the solution in the Fourier domain representsa spontaneous alternative to a study in the direct space. In a Fourier analysis, all thefunctions appearing in the equation are replaced by their Fourier expansions, and the Fouriercoefficients of the wave function become the unknowns of the problem, which in direct spaceapproaches are instead the values of the wave function at the points of a discretization grid.The main advantage of this technique is the correct treatment of the derivatives. Ina reciprocal space analysis, each of the complex exponential functions appearing in theFourier expansions can be derived exactly, without any approximation. On the contrary,in a standard finite difference solution, the derivatives are replaced with finite differenceapproximations involving a certain number of samples of the function. This approximationintroduces a distortion in the dispersion relation, which severely limits the accuracy of high-order solutions. The discretization error decreases if a finer grid is used, at the expenseof an increase of the computational effort. Moreover, in some cases, such as the solutionof the Dirac equation for massless fermions, the discretization of the equation on a directspace grid gives rise to problems such as the appearance of spurious solutions, or fermiondoubling, unless proper, non-standard discretization techniques are applied, while theseissues do not appear using a Fourier solution technique. The solution in the reciprocalspace is especially convenient when a continuum, envelope function description of the deviceis adopted and a slowly varying potential (containing only a limited number of Fouriercomponents) is considered.Here, we illustrate how the use of a particular set of periodic basis functions, obtainedfrom the periodic repetition of sinc (cardinal sine) functions, makes it possible to obtain anexact description of the derivatives in the direct space and a solution approach equivalentto that in the reciprocal space. This approach allows solving the differential problem in avery efficient way in the direct space, without the need to switch to the reciprocal space andvice versa. 2fter a general presentation (in Sec. II) of the method that we propose, in Secs. III andIV we describe how each term of the differential equation is treated, and in Sec. V, wedemonstrate the equivalence of the sinc-based method to a Fourier approach.Finally, in Sec. VI, we apply this method to the solution of the the wave equation in anarmchair graphene nanoribbon (described using an envelope function approach) with a po-tential varying only in the transverse direction. Graphene represents an interesting materialthat, since its isolation in 2004 , has been the focus of a vast theoretical and experimen-tal research activity . Among the many applications that have been proposed for thismaterial, the possibility to use graphene nanoribbons for the implementation of nanoelec-tronic devices is particularly intriguing. The effective wave equation for the envelopefunctions, which in common semiconductors corresponds to the effective mass Schr¨odingerequation (see, for example, Refs. ), in monolayer graphene is represented by theDirac equation . Its numerical solution in graphene ribbons with an arbitrary potentiallandscape has been studied in a few recent papers . The solution of the Dirac equationwithin a slice with a potential varying only in the transverse direction can be exploited toperform a transport analysis for an armchair nanoribbon with a generic potential, if wesubdivide the ribbon into a series of slices in each of which the potential is approximatelyconstant in the longitudinal direction, and then we apply a recursive technique to obtainthe transmission through the overall device . II. MATHEMATICAL METHOD
For the sampling theorem, if a function z ( x ) is band limited with band B ≤ / (2 ∆)(Nyquist criterion), it can be perfectly reconstructed from its samples taken with a sam-pling interval ∆ z ( x ) = + ∞ X n = −∞ z ( n ∆) sinc (cid:18) x − n ∆∆ (cid:19) , (1)where the sinc function is defined as sinc( x ) ≡ sin( πx ) / ( πx ).Moreover, if z ( x ) is periodic with period L and we take N samples within a period (and3 L=N ∆∆ FIG. 1. Four periods of the function g ℓ, ∆ ( x ) for N = 15 and ℓ = 5. thus L = N ∆), we have that z ( x ) = N − X ℓ =0 + ∞ X η = −∞ z (( ℓ + ηN )∆) sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19) = N − X ℓ =0 z ( ℓ ∆) + ∞ X η = −∞ sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19) = N − X ℓ =0 z ( ℓ ∆) g ℓ, ∆ ( x ) , (2)where we have exploited the periodicity of z ( x ) and we have defined the function (periodicwith period L = N ∆) g ℓ, ∆ ( x ) ≡ + ∞ X η = −∞ sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19) (3)(in Fig. 1, we represent 4 periods of the function g ℓ, ∆ ( x ) for N = 15 and ℓ = 5).If we define the scalar product between two functions ψ ( x ) and ψ ( x ) as h ψ ( x ) | ψ ( x ) i = 1∆ Z L ψ ∗ ( x ) ψ ( x ) d x , (4)we can prove (see Appendix A) that the set of functions g ℓ, ∆ ( x ) is orthonormal, i.e. h g j, ∆ ( x ) | g ℓ, ∆ ( x ) i = δ j,ℓ (where δ is the Kronecker delta).Here, we consider a wave equation defined on a domain [0 , L ] and with periodic boundaryconditions. The solutions of this differential equation (i.e., the wave functions), extended byperiodicity with period L , ϕ ( x ) (with ϕ : R → C ) are elements of the infinite-dimensional4ilbert space (with scalar product (4)) L [0 , L ] of the L -periodic functions with finite normand are continuous up to the ( τ − τ beingthe order of the differential equation). In order to solve the problem numerically, here wereformulate it onto the finite-size subspace generated by the N basis functions g that wehave just defined. This allows to preserve the smoothness of the wave functions and, inparticular, to numerically approximate the N slowest-varying solutions, the ones we aremainly interested in. In detail, we approximate each of the terms appearing in the waveequation as a linear combination of the functions g ℓ, ∆ ( x ), with coefficients given (in theirturn) by linear combinations of the samples of the (unknown) wave function ϕ ( x ) at thepoints n ∆ (with n = 0 , . . . , N − N to be an odd number, i.e., N = 2 D + 1(with D an integer). Then, we project the wave equation onto the functions g j, ∆ ( x ), with j = 0 , . . . , N −
1. In this way, exploiting the orthonormality of the g functions, we recastthe wave equation into a system of N linear equations in the N unknowns ϕ ( n ∆) with n = 0 , . . . , N − N × N matrix representation of the wave equationon the chosen set of basis functions). In Sec. V, we will prove that this method is equivalentto a Fourier-Galerkin analysis with cut-off at the spatial frequency D/L (and it is thuscharacterized by the same good convergence properties that have been proved for the Fourier-Galerkin approach, if the involved functions are sufficiently well-behaved: see, for example,Ref. ).In our discussion, we will consider the following terms of the wave equation: a termproportional to the unknown wave function ϕ ( x ), a term proportional to a derivative of thewave function, and a term proportional to the product between the wave function ϕ ( x ) anda function f ( x ) (with f : R → C , in general) which in [0 , L ] represents (or is a functionof) the potential energy and is then extended by periodicity with period L . Other possibleterms appearing in the wave equation will undergo a similar treatment. III. TREATMENT OF THE WAVE FUNCTION AND OF ITSDERIVATIVES
We assume the unknown periodic wave function ϕ ( x ) to be a band limited function witha maximum band B = D/L , such that, sampling it with a step ∆, the Nyquist criterion5 ≤ / (2 ∆) is satisfied.Therefore, the term proportional to the unknown wave function is easily rewritten in thedesired form by applying the sampling theorem with a step ∆ ϕ ( x ) = N − X µ =0 ϕ ( µ ∆) g µ, ∆ ( x ) . (5)Concerning the term proportional to the derivative of the wave function ϕ ( x ), if ϕ ( x ) isband limited with a band equal to B = D/L , also its derivatives have the same property,and thus we can apply the sampling theorem to them, too, obtaining (for the generic s -thderivative): d s ϕ ( x ) d x s = N − X µ =0 (cid:20) d s ϕ ( x ) d x s (cid:21) x = µ ∆ g µ, ∆ ( x ) . (6)We now have to express the values of the derivatives in the points µ ∆ as a function of thesamples of ϕ ( x ). This can be easily achieved by deriving Eq. (5). We have that (cid:20) d s ϕ ( x ) d x s (cid:21) x = µ ∆ = " d s d x s N − X ℓ =0 ϕ ( ℓ ∆) g ℓ, ∆ ( x ) ! x = µ ∆ = N − X ℓ =0 ϕ ( ℓ ∆) + ∞ X η = −∞ (cid:20) d s d x s sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19)(cid:21) x = µ ∆ = N − X ℓ =0 ϕ ( ℓ ∆) α sµℓ , (7)where we have defined as α sµℓ the coefficients with which we have to combine the samples ofthe unknown wave function ϕ ( x ) in order to determine the samples of its s th derivative.In particular, let us develop the calculation for s = 1 (first derivative) and for s = 2(second derivative).For s = 1, we have that dd x sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19) = cos (cid:18) π x − ( ℓ + ηN )∆∆ (cid:19) x − ( ℓ + ηN )∆ − sin (cid:18) π x − ( ℓ + ηN )∆∆ (cid:19) π ∆ ( x − ( ℓ + ηN )∆) (8)and thus (cid:20) dd x sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19)(cid:21) x = µ ∆ = µ − ℓ − ηN = 0( − µ − ℓ − ηN ∆( µ − ℓ − ηN ) if µ − ℓ − ηN = 0. (9)6f we develop the calculation assuming N odd, we have that α µℓ = + ∞ X η = −∞ (cid:20) dd x sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19)(cid:21) x = µ ∆ = µ = ℓ ( − µ − ℓ πN ∆ sin (cid:16) π µ − ℓN (cid:17) if µ = ℓ . (10)Analogously, for s = 2, we have that d d x sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19) = − π ∆ sin (cid:18) π x − ( ℓ + ηN )∆∆ (cid:19) x − ( ℓ + ηN )∆ − (cid:18) π x − ( ℓ + ηN )∆∆ (cid:19) ( x − ( ℓ + ηN )∆) + 2∆ π sin (cid:18) π x − ( ℓ + ηN )∆∆ (cid:19) ( x − ( ℓ + ηN )∆) (11)and thus (cid:20) d d x sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19)(cid:21) x = µ ∆ = − π if µ − ℓ − ηN = 0 − − µ − ℓ − ηN ∆ ( µ − ℓ − ηN ) if µ − ℓ − ηN = 0. (12)Assuming N odd, we obtain that α µℓ = + ∞ X η = −∞ (cid:20) d d x sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19)(cid:21) x = µ ∆ = π N ∆ (1 − N ) if µ = ℓ − π N ∆ ( − µ − ℓ cos (cid:18) ( µ − ℓ ) πN (cid:19) sin (cid:18) ( µ − ℓ ) πN (cid:19) if µ = ℓ . (13)Notably, these expressions coincide with those of the so-called SLAC derivative tech-nique , which were obtained switching to the reciprocal space, operating in that domainand then transforming back (analogously to what we do in Sec. V, where we compare ourmethod with the Fourier one). 7 V. TREATMENT OF THE PRODUCT TERM
The last term that we consider in our analysis is the one proportional to the productbetween the function f ( x ) (the potential energy or a function of it) and the wave function ϕ ( x ).To a first approximation (simplified sinc-based approach), we can proceed as if the prod-uct term were band limited with band B = D/L and thus directly express it as a linearcombination of the functions g ℓ, ∆ ( x ), with coefficients given by the samples of the productterm at the points n ∆ (with n = 0 , . . . , N − n ∆ (analogously to Ref. ). However, since we have assumed a band B for ϕ ( x ), the sameassumption is in general not verified for the product term (unless a constant potential func-tion is considered); and for high-order solutions, this introduces a discrepancy with respectto the solutions from a Fourier analysis.In the following, we discuss a better approximation (advanced sinc-based approach) thatmakes this direct-space method equivalent to the Fourier one.We introduce an odd (for analytical convenience) positive integer number M , whichaccounts for the fact that in general the exploitation of more samples of the potential thanjust those at the positions where we want to evaluate the wave function ϕ can be useful. Inthe calculation we include the samples of the potential function f ( x ) at the N M (= 2 Q + 1,with Q a positive integer) points taken at intervals multiple of ∆ ′ = L/ ( N M ) = ∆ /M . Wereplace the function f ( x ) with the function˜ f ( x ) = MN − X m =0 f ( m ∆ ′ ) + ∞ X η = −∞ sinc (cid:18) x − ( m + ηM N )∆ ′ ∆ ′ (cid:19) = MN − X m =0 f ( m ∆ ′ ) g m, ∆ ′ ( x ) (14)(with a maximum band B ′ = Q/L ), obtained by reconstruction from its
N M samples takenat intervals ∆ ′ .Assuming for the wave function the expression (5), we can therefore write the productterm as: h ( x ) = ˜ f ( x ) ϕ ( x ) = " MN − X m =0 f ( m ∆ ′ ) g m, ∆ ′ ( x ) N − X ℓ =0 ϕ ( ℓ ∆) g ℓ, ∆ ( x ) . (15)8n order to obtain the final matrix representation of our wave equation, we have tocompute only the projections of h ( x ) on the functions g µ, ∆ ( x ) (with µ = 0 , . . . , N − h ( x ) with h ( x ) = N − X µ =0 h g µ, ∆ ( x ) | h ( x ) i g µ, ∆ ( x ) = N − X µ =0 (cid:20) Z L g ∗ µ, ∆ ( χ ) h ( χ ) d χ (cid:21) g µ, ∆ ( x ) . (16)These projections are equal (exploiting Eq. (15)) to1∆ Z L g ∗ µ, ∆ ( χ ) h ( χ ) d χ = N − X ℓ =0 ( ϕ ( ℓ ∆) MN − X m =0 (cid:20) f ( m ∆ ′ ) 1∆ Z L g ∗ µ, ∆ ( χ ) g m, ∆ ′ ( χ ) g ℓ, ∆ ( χ ) d χ (cid:21)) = N − X ℓ =0 ( ϕ ( ℓ ∆) 1 M N MN − X m =0 f ( m ∆ ′ ) ν µℓm ) = N − X ℓ =0 ϕ ( ℓ ∆) β µℓ , (17)where the coefficients β µℓ are given by β µℓ ≡ M N MN − X m =0 f ( m ∆ ′ ) ν µℓm . (18)With some analytical calculations (a possible procedure will be briefly described in thesecond part of Appendix B) it is possible to find a closed form for ν µℓm ≡ M N Z L g ∗ µ, ∆ ( χ ) g m, ∆ ′ ( χ ) g ℓ, ∆ ( χ ) d χ . (19)In detail, if M = 1 we have that ν µℓm = λ ( m, µ ) λ ( m, ℓ ) (20)with λ ( m, µ ) = N if µM = m sin (cid:16) π (cid:16) µ − mM (cid:17)(cid:17) sin (cid:16) πN (cid:16) µ − mM (cid:17)(cid:17) if µM = m . (21)Instead, if M = 1 (and thus ∆ ′ = ∆) we have that ν µℓm = λ ( m, µ ) if ℓ = m ( − ℓ − m ( λ ( ℓ, µ ) − λ ( m, µ ))sin (cid:16) πN ( ℓ − m ) (cid:17) if ℓ = m (22)9ith λ ( m, µ ) = D + 3 D + 1 if m = µ sin (cid:16) πN D ( µ − m ) (cid:17) sin (cid:16) πN ( µ − m ) (cid:17) if m = µ (23)and λ ( m, µ ) = m = µ cos (cid:16) πN ( µ − m ) (cid:17) − ( − µ − m (cid:16) πN ( µ − m ) (cid:17) if m = µ . (24) V. CORRESPONDENCE WITH THE FOURIER ANALYSIS
The differential problem we are interested in can alternatively be solved using a classicalsolution method in the transformed domain. We can replace the wave function ϕ ( x ) andthe potential energy function f ( x ) with their truncated Fourier expansions. In detail, wecan consider only the lowest (in modulus) N spatial frequencies of ϕ ( x ) (i.e., the spatialfrequencies ℓ/L with | ℓ | ≤ D , where N = 2 D + 1) and we can numerically compute thelowest-order M N
Fourier coefficients of f ( x ) by means of a discrete Fourier transform (DFT)of its M N samples within the period L . Then, we can project the resulting equation ontothe N basis functions exp( i πjx/L ) with | j | ≤ D , thereby recasting the problem into an N × N system of linear equations where the Fourier coefficients of the wave function arethe unknowns. In this way, we approximate the infinite-dimensional problem to a finite-sizematrix problem, disregarding the Fourier components corresponding to frequencies greaterin modulus than D/L , both for the unknown wave function and for all the terms appearingin the equation.We will now show that the advanced sinc-based approach described in the previous sec-tions is equivalent to this Fourier-based solution technique.As we have seen, a periodic function ˜ z ( x ) with maximum band D/L can be expressed interms of its N = 2 D + 1 samples at the points n ∆ within the period L as follows:˜ z ( x ) = N − X ℓ =0 ˜ z ( ℓ ∆) g ℓ, ∆ ( x ) . (25)10herefore, we can compute its Fourier series coefficients ˜ Z p noting that the Fourier coeffi-cients [ G ℓ, ∆ ] p of the function g ℓ, ∆ ( x ) are given by (exploiting the definition of g ℓ, ∆ ( x ) as aperiodic repetition of sinc functions )[ G ℓ, ∆ ] p = 1 L F (cid:20) sinc (cid:18) x − ℓ ∆∆ (cid:19)(cid:21) f = pL = 1 N e − i πp ℓN (26)for | p | ≤ D , and 0 otherwise (with F the Fourier transform) and thus˜ Z p = N − X ℓ =0 ˜ z ( ℓ ∆) [ G ℓ, ∆ ] p = 1 N N − X ℓ =0 ˜ z ( ℓ ∆) e − i πp ℓN (27)for | p | ≤ D , and 0 otherwise.It is apparent that Eq. (27) corresponds to the DFT of ˜ z ( x ).Therefore the DFT of a function z ( x ) computed from its N = 2 D + 1 samples coincideswith the exact Fourier series of the function ˜ z ( x ) with maximum band D/L that has thesame samples as z ( x ) in the period L .In particular, this consideration can be applied to the potential function f ( x ): calculatingits DFT on the M N samples (as we do in the Fourier analysis) corresponds exactly tosubstituting f ( x ) with the function ˜ f ( x ) of Eq. (14) (as we do in our advanced sinc-basedapproach).Moreover, as we have seen, when operating in the frequency domain, we disregard thefrequency components outside the interval [ − D/L, D/L ] for all the terms of the differentialproblem. Let us discuss, for each term, the equivalent of this frequency cut-off in the directdomain.Limiting the Fourier components of the unknown wave function ϕ ( x ) to frequencies thathave a modulus less than D/L corresponds exactly to the assumption, we have made whileoperating in the direct space, that ϕ ( x ) is band limited with band D/L (which has allowedus to use the sampling theorem with N = 2 D + 1 samples).The same consideration is valid for the derivatives of ϕ ( x ). The expressions that we havefound for the derivatives can actually be obtained also from a reciprocal space analysis.Indeed, if ϕ ( x ) is periodic and we assume it to be band limited with band B = D/L withFourier coefficients a p , we can write it and its derivatives in the form ϕ ( x ) = D X p = − D a p e i πp xL , (28)11 s ϕ ( x ) dx s = D X p = − D (cid:16) i π pL (cid:17) s a p e i πp xL . (29)The function exp( i πpx/L ) with | p | ≤ D can be seen as a band limited function with band B = D/L and thus expressed in terms of its N samples in the period L : e i πp xL = N − X µ =0 e i πp µN g µ, ∆ ( x ) . (30)Substituting this expression into Eq. (29) we obtain d s ϕ ( x ) dx s = N − X µ =0 " D X p = − D (cid:16) i π pL (cid:17) s a p e i πp µN g µ, ∆ ( x ) . (31)The value of the derivative for x = µ ∆ is the quantity between square brackets. Therefore(if we express a p using Eq. (27)), we have that (cid:20) d s ϕ ( x ) dx s (cid:21) x = µ ∆ = D X p = − D (cid:16) i π pL (cid:17) s N N − X ℓ =0 ϕ ( ℓ ∆) e − i πp ℓN ! e i πp µN = N − X ℓ =0 ϕ ( ℓ ∆) "(cid:18) i π (2 D + 1)∆ (cid:19) s D + 1 D X p = − D p s e i πp µ − ℓ D +1 = N − X ℓ =0 ϕ ( ℓ ∆) α sµℓ . (32)It is easy to verify that the expression for α sµℓ in Eq. (32) yields exactly the same closed-form results that can be found operating in the direct space [for example, for s = 1 and s = 2, those reported in Eqs. (10) and (13)]. Indeed, the expression in Eq. (32) can beeasily obtained also by substituting the Fourier expansion (B8) of g ℓ, ∆ ( x ) into the definition α sµℓ = [ d s g ℓ, ∆ ( x ) /d x s ] x = µ ∆ of Eq. (7).Finally, let us consider the effect of the frequency cut-off on the product term h ( x ).Neglecting the frequency components of this function outside the interval [ − D/L, D/L ]corresponds to limiting its band to the range [ − / (2∆) , / (2∆)], because 1 / (2∆) = (2 D +1) / (2 L ) = ( D/L ) + [1 / (2 L )] but the product term, being periodic with period L , containsonly frequency components multiple of 1 /L . If H ( f ) is the Fourier transform of the function h ( x ) ( f being the spatial frequency), performing this frequency cut-off is equivalent toconsidering, in the Fourier domain, a spectrum H ( f ) = rect(∆ f ) H ( f ), where rect(∆ f ) is12 function equal to 1 between − / (2∆) and 1 / (2∆), and to 0 outside this interval. Theinverse Fourier transform of H ( f ) is equal to h ( x ) = 1∆ Z + ∞−∞ sinc (cid:18) x − χ ∆ (cid:19) h ( χ ) d χ . (33)Since the periodic function h ( x ) is band limited with band B ≤ / (2∆), we can express itin terms of its samples taken with sampling interval ∆ h ( x ) = N − X ℓ =0 (cid:20) Z + ∞−∞ sinc (cid:18) χ − ℓ ∆∆ (cid:19) h ( χ ) d χ (cid:21) g ℓ, ∆ ( x )= N − X ℓ =0 (cid:20) Z L g ∗ ℓ, ∆ ( χ ) h ( χ ) d χ (cid:21) g ℓ, ∆ ( x ) (34)[where we have exploited the fact that the sinc function is even and Eq. (A2)]. This exactlycorresponds to what we have done in the direct space [see Eq. (16)]. Indeed, we can verifythat the expression of the product term obtained operating in the reciprocal domain isequivalent to that achieved in the direct space (see Appendix B).Therefore, with the approach we have described, the solutions in the direct and in thereciprocal space are equivalent and correspond to linear systems of the same size.In order to further clarify this equivalence, let us notice that the periodic and bandlimited function ϕ ( x ) can be equivalently expressed, using both the Fourier expansion andthe sampling theorem, in the following two ways: ϕ ( x ) = D X p = − D ( √ N a p ) e i πp xL √ N = D X ℓ =0 ϕ ( ℓ ∆) g ℓ, ∆ ( x ) . (35)Two different sets of orthonormal basis functions have been used in this equation: thefunctions exp( i πpx/L ) / √ N and the functions g ℓ, ∆ ( x ) [with this choice, the exponentialshave been properly normalized with respect to the scalar product (4)]. Performing thesolution in the reciprocal domain corresponds to using the former basis set, while performingthe solution in the direct space corresponds to using the latter basis set.It is possible to switch from one basis to the other noticing that each of the exponentialfunctions that appear in the Fourier expansion (being band limited with band B ≤ / (2∆))can be expressed in terms of its samples taken with sampling interval ∆ in this way e i πp xL √ N = D X ℓ =0 e i πp ℓN √ N g ℓ, ∆ ( x ) = D X ℓ =0 T ℓp g ℓ, ∆ ( x ) . (36)13he matrix T , made up of the T ℓp elements, can be used to switch from the Fourier basis tothe direct-space one.On the other hand, using Eq. (26), we can write the Fourier expansion of g ℓ, ∆ ( x ) as g ℓ, ∆ ( x ) = D X p = − D [ G ℓ, ∆ ] p e i πp xL = D X p = − D √ N [ G ℓ, ∆ ] p e i πp xL √ N (37)and thus the matrix Θ = T † , made up of the elements Θ pℓ = √ N [ G ℓ, ∆ ] p = exp( − i πpℓ/N ) / √ N = T ∗ ℓp , operates the change from the direct-space basis to the Fourier one.In particular, if the matrix of the direct space linear system is M d and the matrix of thereciprocal-space linear system is M r , we have that M d = T M r T − = T M r T † . VI. APPLICATION TO THE SOLUTION OF THE DIRAC EQUATION INAN ARMCHAIR GRAPHENE NANORIBBON
Transport in graphene is described, within an envelope function approach, by the Diracequation . When considering graphene nanoribbons, at the effective edges of the ribbon(i.e., at the lattice sites just outside the ribbon), we have to enforce Dirichlet boundaryconditions, which couple the envelope functions corresponding to the two inequivalent Diracpoints ~K and ~K ′ . In particular, here we take into consideration the solution of the Diracequation in the case in which the potential in the ribbon depends only on the transverse co-ordinate, i.e., in which the potential does not vary in the longitudinal direction. In this case,the four envelope functions F S ~P ( ~r ) (corresponding to the two Dirac points ~P = ~K, ~K ′ andto the two sublattices S = A, B ) can be written as the product of a confined component inthe transverse y direction Φ S ~P ( y ) and of a propagating wave in the longitudinal x direction(with longitudinal wave vector κ x ): F S ~P ( ~r ) = Φ S ~P ( y ) exp( iκ x x ).As shown in Refs. , the problem can be reformulated into a differential problem withperiodic boundary conditions on a doubled domain (2 ˜ W instead of ˜ W , which representsthe effective width of the ribbon, i.e., the distance between the two effective edges wherethe Dirichlet boundary conditions had to be enforced in the original problem). This newproblem can be expressed in the following form: σ z (cid:18) d ~ϕ ( y ) d y + i ˜ K ~ϕ ( y ) (cid:19) + σ x f ( y ) ~ϕ ( y ) = − κ x ~ϕ ( y ) ~ϕ (2 ˜ W ) = ~ϕ (0) , (38)14here the σ ’s are the Pauli matrices and ~ϕ ( y ) is a two-component function directly relatedto the four envelope functions of graphene in the following way: ~ϕ ( y ) = e − i ˜ Ky Φ A ~K ( y )Φ B ~K ( y ) if y ∈ [0 , ˜ W ] e i ˜ K (2 ˜ W − y ) i Φ A ~K ′ (2 ˜ W − y )Φ B ~K ′ (2 ˜ W − y ) if y ∈ [ ˜ W , W ]. (39)Moreover, K = 4 π/ (3 a ) (where a is the graphene lattice constant) is the modulus of thetransversal co-ordinate of the Dirac points and ˜ K = K − round ( K/ ( π/ ˜ W )) · π/ ˜ W is aquantity such that exp( i K ˜ W ) = exp( i K ˜ W ) and that | ˜ K | < π/ ˜ W . We consider ˜ K in-stead of K , in such a way that slowly varying ~ϕ ( y ) functions correspond to slowly varyingtransverse components of the envelope functions, i.e., those in which we are most interested.Additionally, f ( y ) = ( U ( ˜ W − | ˜ W − y | ) − E ) /γ is the potential term, where U ( y ) is thepotential energy, E is the Fermi energy, γ = ~ v F , ~ is the reduced Planck constant, and v F is the Fermi velocity in graphene. In this formulation, the problem has to be solved overthe domain [0 , W ] and the boundary condition for ~ϕ ( y ) is periodic.The periodic boundary condition makes it possible to solve the problem in the Fourierdomain or equivalently in the basis of g functions. Working in the basis of g functions, wehave rewritten the equation as follows: σ z " D X µ =0 D X ℓ =0 ~ϕ ( ℓ ∆) α µℓ ! g µ, ∆ ( x ) + i ˜ K D X µ =0 ~ϕ ( µ ∆) g µ, ∆ ( x ) + σ x D X µ =0 D X ℓ =0 ~ϕ ( ℓ ∆) β µℓ ! g µ, ∆ ( x ) = − κ x D X µ =0 ~ϕ ( µ ∆) g µ, ∆ ( x ) . (40)Projecting this equation (using the scalar product (4)) onto the generic function g j, ∆ ( x ), weobtain σ z " D X ℓ =0 ~ϕ ( ℓ ∆) α jℓ ! + i ˜ K ~ϕ ( j ∆) + σ x D X ℓ =0 ~ϕ ( ℓ ∆) β jℓ ! = − κ x ~ϕ ( j ∆) (41)for j = 0 , . . . , D , or equivalently D X ℓ =0 h σ z ( α jℓ + i ˜ Kδ j,ℓ ) + σ x β jℓ i ~ϕ ( ℓ ∆) = − κ x ~ϕ ( j ∆) (42)15 (nm) P o t en t i a l ene r g y ( e V ) FIG. 2. Potential energy considered in the armchair ribbon.
Fourier analysisAdvanced sinc−basedanalysis κ x Re( ) ( n m ) − (nm ) −1 κ x I m ( ) −0.3 FIG. 3. Position on the Gauss plane of all the longitudinal wave vectors κ x obtained with theFourier analysis and with our advanced sinc-based approach. for j = 0 , . . . , D . These N = 2 D + 1 equations can be written in matrix form as aneigenproblem with eigenvalues − κ x and eigenvectors containing the values ~ϕ ( j ∆) of theunknown function ~ϕ ( y ) at the N points of the considered grid.In order to adopt the first approximate approach for the product term described at thebeginning of Sec. IV, we just have to replace β jℓ with f ( j ∆) δ j,ℓ in Eq. (42).For example, we have considered a potential U ( y ) (represented in Fig. 2) equal to thesum of two Gaussian functions: U ( y ) = 0 .
15 eV · e − ( y −
300 nm) + 0 .
15 eV · e − ( y −
500 nm) (43)in a 8131 dimer line ( ∼ µ m) wide armchair nanoribbon, considering an electron injectionenergy E = 0 . AK | | ( n m ) − Fourier analysisAdvanced sinc−basedanalysis y FIG. 4. Envelope function Φ
A ~K ( y ) corresponding to the eigenvalue with the largest real part,computed with a Fourier analysis and with our advanced sinc-based approach. ( n m ) − κ I m ( ) x κ x (nm ) −1 Re( )
Simplified sinc−basedanalysis
FIG. 5. Position on the Gauss plane of all the longitudinal wave vectors κ x obtained with thesimplified sinc-based approach. The highest order wave vectors strongly differ from those obtainedwith a Fourier analysis (see Fig. 3). N = 99 and M = 21 with the advanced approach in the direct space, and we compare themwith those achieved by means of a Fourier analysis. In particular, in Fig. 3, we show theposition on the Gauss plane of all the obtained longitudinal wave vectors κ x ; and in Fig. 4,we represent the behavior of the envelope function Φ A ~K ( y ) corresponding to the eigenvaluewith the largest real part. As expected, the two methods yield exactly the same results.Instead, in Fig. 5, we show the values of the longitudinal wave vectors κ x obtained withthe simplified treatment of the product term. As we can see from the highest order wavevectors, this method is not equivalent to the Fourier one, even though it yields quite good17esults for the low-order modes. VII. CONCLUSION
We have proposed a method, based on a set of basis functions resulting from the periodicrepetition of sinc functions, for the numerical solution, in the direct space, of differentialproblems and, in particular, of quantum wave equations, with periodic boundary conditions.This approach is equivalent to a Fourier analysis and (once we have chosen a givendiscretization) requires the solution of a linear system with the same matrix size. It allowsperforming a very efficient solution in the direct space, without the need to Fourier transformthe input functions and finally to inverse transform the solutions.We have applied this method to efficiently solve the Dirac equation in an armchairgraphene nanoribbon in the direct space, verifying its equivalence to a reciprocal spaceapproach.We believe that this treatment also helps in clarifying the exact origin of the differencesin the convergence behavior observed between direct space and reciprocal space approaches.
Appendix A: Orthonormality of the basis functions
The orthonormality of two sinc functions centered on different sampling points can beeasily demonstrated in the following way: Z + ∞−∞ sinc (cid:18) x − j ∆∆ (cid:19) sinc (cid:18) x − ℓ ∆∆ (cid:19) dx = h sinc (cid:16) x ∆ (cid:17) ∗ sinc (cid:16) x ∆ (cid:17)i x =( ℓ − j )∆ = (cid:8) F − (cid:2) ∆ rect(∆ f ) (cid:3)(cid:9) x =( ℓ − j )∆ = ∆ sinc( ℓ − j ) = ∆ δ j,ℓ , (A1)where ∗ is the convolution operator and F − is the inverse Fourier transform.18oreover, if the function z ( x ) is periodic with period L , we have that1∆ Z L g ∗ j, ∆ ( χ ) z ( χ ) dχ = 1∆ Z L g j, ∆ ( χ ) z ( χ ) dχ = 1∆ + ∞ X η = −∞ Z L sinc (cid:18) χ − ( j + ηN )∆∆ (cid:19) z ( χ ) dχ = 1∆ + ∞ X η = −∞ Z ( − η +1) L − ηL sinc (cid:18) ξ − j ∆∆ (cid:19) z ( ξ + ηN ∆) dξ = 1∆ + ∞ X η = −∞ Z ( − η +1) L − ηL sinc (cid:18) ξ − j ∆∆ (cid:19) z ( ξ ) dξ = 1∆ Z + ∞−∞ sinc (cid:18) ξ − j ∆∆ (cid:19) z ( ξ ) dξ , (A2)where we have changed the integration variable from χ to ξ = χ − ηN ∆ = χ − ηL .From the previous observation, we derive the orthonormality of the g functions1∆ Z L g ∗ j, ∆ ( x ) g ℓ, ∆ ( x ) dx = 1∆ Z + ∞−∞ sinc (cid:18) ξ − j ∆∆ (cid:19) g ℓ, ∆ ( x ) dξ = 1∆ + ∞ X η = −∞ Z + ∞−∞ sinc (cid:18) ξ − j ∆∆ (cid:19) sinc (cid:18) x − ( ℓ + ηN )∆∆ (cid:19) = 1∆ + ∞ X η = −∞ ∆ δ j,ℓ + ηN = + ∞ X η = −∞ δ j,ℓ δ η, = δ j,ℓ (A3)(with j, ℓ = 0 , . . . , N − Appendix B: Notes on the product term
Operating in the reciprocal space, if we substitute ϕ ( x ) = D X r = − D a r e i πr xL , ˜ f ( x ) = Q X q = − Q f q e i πq xL (B1)(with f q being the DFT coefficients of f ( x ) computed on its N M = 2 Q +1 samples) into theterm h ( x ) = ˜ f ( x ) ϕ ( x ) and then we consider only the frequency components of the productwithin the interval [ − D/L, D/L ], we obtain that h ( x ) = Q X q = − Q D X r = − D | q + r |≤ D f q a r e i π ( q + r ) xL = D X p = − D D X r = − D | p − r |≤ Q f p − r a r e i πp xL . (B2)19ince the functions exp( i πpx/L ) with | p | ≤ D are band limited with band B ≤ / (2∆),they can be expressed in terms of their samples taken with sampling interval ∆ e i πp xL = N − X µ =0 e i πp µN g µ, ∆ ( x ) . (B3)Moreover, exploiting Eq. (27), the coefficients a r and f p − r can be replaced with a r = 1 N N − X ℓ =0 ϕ ( ℓ ∆) e − i πr ℓN ,f p − r = 1 M N MN − X m =0 f ( m ∆ ′ ) e − i π ( p − r ) mMN . (B4)After these substitutions, we obtain that h ( x ) = N − X µ =0 " N − X ℓ =0 ϕ ( ℓ ∆) β µℓ g µ, ∆ ( x ) (B5)with β µℓ = 1 M N MN − X m =0 f ( m ∆ ′ ) ν µℓm (B6)and ν µℓm = D X p = − D D X r = − D | p − r |≤ Q h e i πN ( µ − mM ) i p h e i πN ( ℓ − mM ) i − r = D X p = − D D X ˜ r = − D | p +˜ r |≤ Q h e i πN ( µ − mM ) i p h e i πN ( ℓ − mM ) i ˜ r (B7)(with ˜ r = − r ).Equations (B5) and (B6) coincide with Eqs. (16)–(18) obtained operating in the directspace. We can prove that Eq. (B7) coincides with Eq. (19), too. In detail, using Eq. (26),we can write the Fourier expansion of the g functions as follows: g µ, ∆ ( x ) = 1 N D X p = − D e i πp x − µ ∆ L , (B8) g m, ∆ ′ ( x ) = 1 M N Q X q = − Q e i πq x − m ∆ ′ L , (B9) g ℓ, ∆ ( x ) = 1 N D X r = − D e i πr x − ℓ ∆ L . (B10)20ubstituting Eqs. (B8)–(B10) into Eq. (19), we find that ν µℓm = 1∆ 1 N D X p = − D Q X q = − Q D X r = − D e i πN ( pµ − qmM − rℓ ) Lδ q,p − r = D X p = − D D X r = − D | p − r |≤ Q h e i πN ( µ − mM ) i p h e i πN ( ℓ − mM ) i − r , (B11)which coincides with Eq. (B7).If M >
1, the condition | p − r | ≤ Q (i.e., | p + ˜ r | ≤ Q ) is always satisfied and thus ν µℓm = D X p = − D h e i πN ( µ − mM ) i p D X ˜ r = − D h e i πN ( ℓ − mM ) i ˜ r = λ ( m, µ ) λ ( m, ℓ ) , (B12)which corresponds to Eq. (20).Instead, if M = 1, the condition | p + ˜ r | ≤ Q has to be taken into consideration. FromEq. (B7), we see that for each value of p , we have to sum over the integers ˜ r for which − D ≤ ˜ r ≤ D and − p − D ≤ ˜ r ≤ − p + D (note that Q = D if M = 1). This set of values of˜ r corresponds for p = 0 to − D ≤ ˜ r ≤ D , for p > − D ≤ ˜ r ≤ − p + D , and for p < − p − D ≤ ˜ r ≤ D . Therefore, Eq. (B7) can be rewritten in the following way: ν µℓm = D X ˜ r = − D h e i πN ( ℓ − m ) i ˜ r + D X p =1 ( h e i πN ( µ − m ) i p − p + D X ˜ r = − D h e i πN ( ℓ − m ) i ˜ r ) + − X p = − D ( h e i πN ( µ − m ) i p D X ˜ r = − p − D h e i πN ( ℓ − m ) i ˜ r ) . (B13)Computing the first sum and noting that the third sum is the complex conjugate of thesecond one, we can rewrite this expression as: ν µℓm = N δ m,ℓ + 2 Re ( D X p =1 (cid:26) h e i πN ( µ − m ) i p − p + D X ˜ r = − D h e i πN ( ℓ − m ) i ˜ r (cid:27)) . (B14)If ℓ = m , this is equal to ν µℓm = N + 2 Re ( D X p =1 nh e i πN ( µ − m ) i p ( N − p ) o) = λ ( m, µ ) . (B15)21f ℓ = m instead, we have that ν µℓm = 2 Re ( D X p =1 (cid:26) h e i πN ( µ − m ) i p ( − ℓ − m i sin (cid:16) πN ( ℓ − m ) (cid:17) h e − i πN p ( ℓ − m ) − i (cid:27)) = ( − ℓ − m sin (cid:16) πN ( ℓ − m ) (cid:17) " Im (cid:26) D X p =1 h e i πN ( µ − ℓ ) i p (cid:27) − Im (cid:26) D X p =1 h e i πN ( µ − m ) i p (cid:27) = ( − ℓ − m ( λ ( ℓ, µ ) − λ ( m, µ ))sin (cid:16) πN ( ℓ − m ) (cid:17) . (B16)Equations (B15) and (B16) correspond to Eq. (22). REFERENCES R. Stacey, Phys. Rev. D , 468 (1982). C. M. Bender, K. A. Milton and D. H. Sharp, Phys. Rev. Lett. , 1815 (1983). L. Susskind, Phys. Rev. D , 3031 (1977). K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grig-orieva and A. A. Firsov, Science , 666 (2004). A. K. Geim and K. S. Novoselov, Nature Mater. , 183 (2007). A. K. Geim, Science , 1530 (2009). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselovand A. K. Geim, Rev. Mod.Phys. , 109 (2009). D. R. Cooper, B. D’Anjou, N. Ghattamaneni, B. Harack, M. Hilke, A. Horth, N. Ma-jlis, M. Massicotte, L. Vandsburger, E. Whiteway and V. Yu, “Experimental Review ofGraphene,” ISRN Condensed Matter Physics , art. number 501686 (2012). S. Mikhailov,
Physics and Applications of Graphene (InTech, Rijeka, 2011). M. R. Connolly, R. K. Puddy, D. Logoteta, P. Marconcini, M. Roy, J. P. Griffiths,G. A. C. Jones, P. A. Maksym, M. Macucci and C. G. Smith, Nano Lett. , 5448 (2012). D. Logoteta, P. Marconcini, M. R. Connolly, C. G. Smith and M. Macucci, “Numericalsimulation of scanning gate spectroscopy in bilayer graphene in the Quantum Hall regime,”in
Proceedings of IWCE 2012 (IEEE Conference Proceedings) (2012), art. number 6242841,DOI: 10.1109/IWCE.2012.6242841. P. Marconcini, A. Cresti, F. Triozon, G. Fiori, B. Biel, Y.-M. Niquet, M. Macucci andS. Roche, ACS Nano , 7942 (2012). 22 P. Marconcini, A. Cresti, F. Triozon, G. Fiori, B. Biel, Y.-M. Niquet, M. Macucci andS. Roche, “Electron-hole transport asymmetry in Boron-doped Graphene Field EffectTransistors,” in
Proceedings of IWCE 2012 (IEEE Conference Proceedings) (2012), art.number 6242844, DOI: 10.1109/IWCE.2012.6242844. H. Raza,
Graphene Nanoelectronics: Metrology, Synthesis, Properties and Applications (Springer, Heidelberg, 2012). F. Schwierz, Nature Nanotech. , 487 (2010). G. Iannaccone, G. Fiori, M. Macucci, P. Michetti, M. Cheli, A. Betti and P. Marconcini,“Perspectives of graphene nanoelectronics: probing technological options with model-ing,” in
Proceedings of IEDM 2009 (IEEE Conference Proceedings) (2009), p. 245, DOI:10.1109/IEDM.2009.5424376. P. Y. Yu and M. Cardona,
Fundamentals of Semiconductors (Springer, Heidelberg, 2003). V. Mitin, V. Kochelap and M. A. Stroscio,
Quantum Heterostructures. Microelectronicsand Optoelectronics (Cambridge University Press, Cambridge, 1999). P. Marconcini, M. Macucci, G. Iannaccone, B. Pellegrini and G. Marola, Europhys. Lett. , 574 (2006). R. S. Whitney, P. Marconcini and M. Macucci, Phys. Rev. Lett. , 186802 (2009). P. Marconcini, M. Macucci, G. Iannaccone and B. Pellegrini, Phys. Rev. B , 241307(R)(2009). P. Marconcini, M. Macucci, D. Logoteta and M. Totaro, Fluct. Noise Lett. , 1240012(2012). P. Marconcini, M. Totaro, G. Basso and M. Macucci, AIP Advances , 062131 (2013). T. Ando, J. Phys. Soc. Jpn. , 777 (2005). P. Marconcini and M. Macucci, “The k · p method and its application to graphene, carbonnanotubes and graphene nanoribbons: the Dirac equation,” La Rivista del Nuovo Cimento , 489 (2011). J. Tworzyd lo, C. W. Groth and C. W. J. Beenakker, Phys. Rev. B , 235438 (2008). A. R. Hern´andez and C. H. Lewenkopf, Phys. Rev. B , 155439 (2012). I. Snyman, J. Tworzyd lo, C. W. J. Beenakker, Phys. Rev. B , 045118 (2008). D. Logoteta, P. Marconcini and M. Macucci, “Numerical simulation of shot noise in dis-ordered graphene,” in
Proceedings of ICNF 2013 (IEEE Conference Proceedings) (2013),art. number 6578889, DOI: 10.1109/ICNF.2013.6578889.23 S. Haykin,
An Introduction to Analog and Digital Communications (Wiley, New York,1989). S. Datta,
Quantum Phenomena (Addison-Wesley, Reading, Massachussets, 1989). C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang,
Spectral Methods: Fundamentalsin Single Domains (Springer, Berlin, 2006). S. D. Drell, M. Weinstein and S. Yankielowicz, Phys. Rev. D , 1627 (1976). J. F¨orster, A. Saenz and U. Wolff, Phys. Rev. E , 016701 (2012). M. Fagotti, C. Bonati, D. Logoteta, P. Marconcini and M. Macucci, Phys. Rev. B ,241406(R) (2011). P. Marconcini, D. Logoteta, M. Fagotti and M. Macucci, “Numerical solution of the Diracequation for an armchair graphene nanoribbon in the presence of a transversally variablepotential,” in