Paolo Marconcini, Demetrio Logoteta, Maurizio Fagotti, Massimo Macucci, "Numerical solution of the Dirac equation for an armchair graphene nanoribbon in the presence of a transversally variable potential", Proceedings of the 14th International Workshop on Computational Electronics, Pisa, Italy, 26-29 October 2010, edited by Giovanni Basso, Massimo Macucci, IEEE Conference Proceedings, IEEE Catalog Number: CFP10462-PRT, ISBN: 978-1-4244-9381-4, Ed. Plus, 53 (2010).

Numerical solution of the Dirac equation for an armchair graphene nanoribbon in the presence of a transversally variable potential ∗ Dipartimento

† Dipartimento

P. Marconcini∗‡ , D. Logoteta∗ , M. Fagotti† , M. Macucci∗ di Ingegneria dell’Informazione, Universit`a di Pisa, Via Caruso 16, I-56122 Pisa, Italy di Fisica and INFN, Universit`a di Pisa, Largo Bruno Pontecorvo 3, I-56127 Pisa, Italy ‡ e-mail: [email protected]

Abstract—We compute, using a continuum model based on the Dirac equation, the transverse modes and the longitudinal wave vectors in an armchair ribbon, in the presence of a transversally variable external potential. We show that the application of a standard finite difference method to this problem is not effective, since it can lead to the appearance of spurious solutions and to poor precision. We show how an approach based on Fourier analysis allows to overcome this difficulty and to drastically speed up the calculations, using well-established FFT algorithms.

I. I NTRODUCTION

W

y=0

Transport in armchair graphene nanoribbons (Fig. 1(a)) can be studied either with an atomistic approach (for example using a tight-binding model) or with a continuum, envelope function method. This latter approach, albeit missing the atomic details of the structure, allows us to treat relatively large ribbons with a still reasonable numerical effort: the transport properties of the nanoribbon can then be investigated subdividing it into transverse regions with a longitudinally constant potential and using, for example, a scattering matrix approach. This method has been successfully applied to the study of transport in large nanoribbons with an external potential U that varies only in the transport direction [1]. In this case, the potential is assumed constant in each section and thus an analytical form of the envelope functions is known. However, in general the potential inside each section varies in the transverse direction (Fig. 1(b)) and an analytical form for the envelope functions is not available. As a first step for a future transport analysis, here we describe a procedure to determine the envelope functions of the armchair nanoribbon in the general case of an external potential depending also on the transverse coordinate, pointing out the problems we have met with a finite difference approach and describing an alternative method, based on the Fourier analysis, which simplifies the solution of the problem. II. T HE NUMERICAL PROBLEM The overall electron wavefuction in a graphene nanoribbon can be written [2]    j )ϕ(r − R  j) , ψj (R (1) ψ(r) = j=A,B R j

y=W

a/2

W

(a)

a/2

(b)

Fig. 1. a) Sketch of an armchair nanoribbon. Since we consider a relatively large ribbon, we show the atomic details only for the edge regions. The hydrogen atoms on the edge rows are for passivation. b) A possible (in this case Lorentzian-shaped) external potential U .

 A and R B where ϕ(r) is the 2pz atomic orbital of carbon, R are the positions of the atoms of the two sublattices A and B of graphene, while ψA (r) and ψB (r) are given, in an armchair nanoribbon, by      ψA (r) = eiK·r FAK (r) − ieiK ·r FAK (r) . (2)     ψB (r) = ieiK·r FBK (r) + eiK ·r FBK (r)   K   and S = A, B) are The functions FSP (r) (with P = K, the envelope functions relative to the two unequivalent Dirac  and K   and to the two sublattices A and B of points K graphene, which vary slowly around the corresponding Dirac point. They have to satisfy the Dirac equation:    U (r)I γσ · (−i∇) F (r) = E F (r) (3)  γσ · (−i∇) U (r)I

978-1-4244-9382-1/10/$26.00 ©2010 IEEE

with the boundary conditions: ⎧   K ⎪ ΦK ⎪ A (0) = iΦA (0) ⎪   ⎨ K ΦB (0) = iΦK B (0) ˜ K   ˜ i2K W ˜ ⎪ ΦK ΦA (W ) ⎪ A (W ) = ie ⎪ ⎩ K  ˜  ˜ ˜ ΦB (W ) = iei2K W ΦK ( W ) B

(5)

(where K=4π/(3a) specifies the position of the Dirac points:  = −K yˆ, while K   = K yˆ). K III. F INITE DIFFERENCE APPROACH To solve the previous problem, we have first adopted a finite difference approach. We consider a uniform mesh of N points yi along the transverse direction, with y1 = 0 and ˜ . We write the first and the second equation of the yN = W system (4) in each of the points yi with i = 1, . . . , N − 1 and the third and the fourth in each of the points yi with i = 2, . . . , N , using a 5-points discretization formula for the first derivatives (in its symmetric form inside the ribbon, in its asymmetric forms at the borders) [4]. In order to exploit the 4 boundary conditions, we do not consider explicitly the     K K K 4 quantities ΦK A (yN ), ΦB (yN ), ΦA (y1 ) and ΦB (y1 ), but, wherever they appear in the discretized equations, they are replaced by the quantities given by Eq. (5). In this way, for each given injection energy E, the problem reduces to an algebraic eigensystem, where the κx s are the eigenvalues, while the vectors containing the discretized form     K K K of ΦK A , ΦB , ΦA and ΦB are the eigenvectors. Since we are interested in the low-order Φ (those having slow variations in space), it is not necessary to consider a very large number N of discretization points in the calculation.

0.8 0.04

0.6

0

0.4 Im {κx } (nm −1)

(where F (r) is the column vector     K K K K (FA (r), FA (r), FB (r), FB (r)), I is the 2 × 2 identity matrix and σ is the column vector of Pauli spin matrices (σx , σy )). If y is the transverse coordinate, while x is the longitudinal one, at the edges of the nanoribbon we have to enforce that for every value of x the functions ψA (r) ˜ (where and ψB (r) vanish for y = 0 and y = W + a ≡ W W is the width of the nanoribbon, a is the graphene lattice constant, y = 0 is taken a/2 below the lower edge of the ˜ is taken a/2 above the upper edge of the ribbon, and y = W ribbon: see Fig. 1(a)) [3]. If we restrict ourselves to a transverse section of the nanoribbon where U depends only on y, the envelope   functions can be factorized as FSP (r) = ΦP S (y) exp(iκx x). In this condition, the previous relations reduce to the following system ⎧    E−U (y) K ⎪ − ddy ΦK ΦB (y) = κx ΦK ⎪ A (y) + A (y) γ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪    E−U (y) K ⎪ K ⎪ ΦA (y) + ddy ΦK ⎨ B (y) = κx ΦB (y) γ (4) ⎪    E−U (y) K ⎪ d K K ⎪ Φ (y) + Φ (y) = κ Φ (y) x A ⎪ B dy A γ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪    ⎩ E−U (y) ΦK d K K A (y) − d y ΦB (y) = κx ΦB (y) γ

−0.1716 −0.1708

−0.04 −0.17

0.2 0 −0.2 −0.4 −0.6 −0.8 −0.2 −0.15 −0.1 −0.05

0

0.05

0.1

0.15

0.2

Re {κx } (nm −1)

Fig. 2. Eigenvalues computed, using the finite difference approach, for a 1 μm wide ribbon, with an injection energy of 0.1 eV and no external potential. The larger, empty dots represent the physically meaningful κx s that we consider, while the smaller, filled dots are those we neglect. The inset shows, for the most negative real eigenvalues, the agreement between the numerical results (empty dots) and the values obtained with an analytical calculation (crosses).

Moreover, since the matrix of the eigenvalue problem is very sparse, we use a routine optimized for sparse matrices, which exploits the Implicitly Restarted Arnoldi Method and is part of the ARPACK package. Unfortunately, as a consequence of the discretization, from the numerical solution of the eigensystem we obtain, besides the physically meaningful eigenvalues, many complex spurious eigenvalues. Since, as one could expect, the physically meaningful eigenvalues correspond to the most slowly varying Φs, while the spurious ones are associated with fast-varying envelope functions, for which the discretization procedure shows its limits, the problem can be solved ordering the Φs on the basis of their frequency spectrum and taking into consideration only the first Φs (and the corresponding κx s). These solutions correspond to the propagating modes and to the lowest-order evanescent modes, which are exactly those we need to compute the transmission inside the structure. In Fig. 2 we report the numerical results obtained for a 1 μm wide ribbon and an injection energy of 0.1 eV, in the absence of an external potential. The larger empty dots represent the physically meaningful κx s that we consider, while the smaller filled dots are higher-order physical eigenvalues and spurious eigenvalues. The inset shows, for the more negative real eigenvalues, the agreement between the numerical results (empty dots) and the values obtained with an analytical calculation, that is available for null external potential (crosses). As we can see, the physically meaningful eigenvalues that come out from our calculation agree very well with the analytically expected ones. IV. F OURIER APPROACH In the attempt to avoid the presence of spurious solutions, we have then applied a different technique, based on Fourier analysis. Defining f (y) ≡ (U (y) − E)/γ and introducing the two 2-

dimensional vectors: 

 ΦK (y) A , ϕ  ↑ (y) =  ΦK B (y)

ϕ  ↓ (y) =



iΦK A (y)  K iΦB (y)

the system (4) can, more compactly, be written as:  (f (y)σy − iκx σz − iI ddy ) ϕ↑ (y) = 0 d ϕ↓ (y) = 0 (f (y)σy − iκx σz + iI d y )

 ,

(6)

potentials. As a consequence, if we choose a positive integer D in such a way that   π   + K   max |f | , (13) (n0 ± D) ˜  W

we can approximately write the matrix M as the direct sum of three matrices: M− (containing the blocks Mn, (7) with n, < n − D), M (with the blocks M 0 0 n, with n0 − D ≤ n, ≤ n0 + D) and M+ (with the blocks Mn, ˜ ]) and the boundary conditions (5) become: with n, > n0 + D). While in the submatrices M− and M+ (with y ∈ [0, W the blocks An dominate on the blocks Bn, and thus M− ϕ  ↑ (0) = ϕ  ↓ (0) (8) and M+ are approximately diagonal, inside the submatrix M0 ˜ ˜ ˜) . ϕ  ↑ (W ) = ei2K W ϕ  ↓ (W both the An and the Bn, blocks have to be considered. The eigenvalues of the matrices M− and M+ are approximately Now let us introduce the function ϕ  (y), defined for y ∈ ˜ equal to −κ x = ±i((πn/W ) + K) and (being n far from ˜ ] in this way: [0, 2W ) represent imaginary eigenvalues with large modulus; thus n 0 ˜] for y ∈ [0, W ϕ  ↑ (y) they correspond to high-order evanescent modes, which we ϕ  (y) = ˜ ˜ − y) for y ∈ [W ˜ , 2W ˜ ] . (9) can neglect in a transport analysis. Therefore, for a sufficiently ei2K W ϕ  ↓ (2W large cut-off threshold D the submatrix M0 provides all of After some algebra, it can be proved that the problem (7)-(8) the information we are interested in. The longitudinal wave ˜ ]: is equivalent to solving for y ∈ [0, 2W  can be vectors κx are the eigenvalues of M0 . The functions ϕ ⎧  by computing the reconstructed from the eigenvectors of M d 0 ⎨ − σ ˜ − |W ˜ − y|)σx ϕ  (y) = κx ϕ + f (W  (y) z . Finally, the envelope (cut) series with Fourier coefficients  a n dy (10) ˜ ⎩ ˜ ) = ei2K W functions are completely determined by the ϕ  s. ϕ  (0) ϕ  (2W The Fourier coefficients f are numerically obtained using the (we have halved the number of equations of the system Discrete Fourier Transform of the samples of the function doubling the solution domain). ˜ − |W ˜ − y|). In order to gain in terms of computaf (W If we consider the function e−iKy ϕ  (y), we see that, due to tional speed, the calculation has been performed using the the boundary condition, it has the same value in y = 0 and in Fast Fourier Transform, with some well established versions ˜ . Therefore we can extend it with continuity outside y = 2W optimized for the case of real even functions [5]. An FFT ˜ ˜. [0, 2W ] making it periodic with period 2W algorithm has been used also for the final step, in which we ˜ the Analogously, we can extend by periodicity with period 2W obtain the function ϕ(y) from the Fourier coefficients a . ˜ − |W ˜ − y|), which has the same value in real function f (W This method does not suffer from the spurious eigenvalue ˜ , obtaining in this way an even periodic y = 0 and in y = 2W problem that we have observed (for this particular problem) function. using the finite difference technique. Moreover, the previously We can rewrite the first equation of (10) in terms of the described numerical tricks have allowed to drastically reduce coefficients a and fm of the Fourier series expansion of the computational times. ˜ − |W ˜ − y|),  (y) and f (W the periodic extensions of e−iKy ϕ In Fig. 3 we show the position on the Gauss plane of the respectively, obtaining for the generical index n: numerically obtained eigenvalues, using the Fourier approach, +∞ for a 1 μm wide ribbon, with an injection energy of 0.1 eV,

    πn i + K σz δ,n + f|n−| σx a = −κxan , (11) in the presence of a Lorentzian-shaped potential. We note ˜ W the presence around the origin of values of κx which are =−∞ neither real nor purely imaginary. These eigenvalues, which where δ,n is the Kronecker delta. Since for each considered appear also using the finite difference method and correspond value of n we obtain an equation of the form (11), globally we to quite slowly-varying modes, are not spurious, but physically have a (in principle infinite-dimensional) eigensystem Ma = meaningful values. −κxa, where a is a vector given by the juxtaposition of the 2-dimensional vectors an , while M is a matrix, the 2 × 2 In Fig. 4 we show, for the same case of Fig. 3, and in particular for the value of κx with the greatest real part, the square subblock in position (n, ) of which is equal to modulus of the corresponding slowly-varying envelope func  πn   K K + K σz δ,n + f|n−| σx = An δ,n + Bn, . (12) tions ΦA (y) and ΦA (y), and of the fast-varying transverse Mn, = i ˜ W component of ψA (r), as a function of y. As we can see, in The blocks An give a small contribution near n0 = this case the probability density has the greatest values in ˜ /π), while have a great value for indices n far the regions where the potential has low values, as one could −round(K W from n0 . Instead the contribution of the blocks Bn, vanishes expect. for large values of |n − | (i.e. sufficiently far away from In Fig. 5 instead we show what happens if we apply an external the diagonal of the matrix) for not too fast-varying external potential with the opposite sign with respect to Figs. 3-4.

0.018

0.6

0.0135

K 2 |ΦA | (nm −1)

0.8

0.2

0

−0.2 −0.4 −0.6 0

0.05

0.1

0.15

0.2

Re {κx } (nm −1)

|ΦAK | 2 (nm −1) |ΦAK’ |2 (nm −1)

0.0045 −0.25

0.056

−0.2

0.042

−0.15

0.028

−0.1

0.014 0

−0.05 0

200

400

600

800

0 1000

y (nm)

0.0035 0.0028 0.0021 0.0014 0.0007 0

Fig. 5. For a potential with opposite sign with respect to that used for Fig. 3 and 4 (represented on a reversed scale by a dashed line in the bottom panel) and for the eigenvalue with the largest real part: square modulus (computed   K using the Fourier approach) of the ΦK A (y), ΦA (y), and the transverse component of ψA ( r), as a function of position.

0.0035 0.0028 0.0021 0.0014 0.0007 0 0.24

0.016

0.192

0.012

0.144

0.008

0.096

0.004

U (eV)

0.02

0

0.009

U (eV)

Fig. 3. Eigenvalues obtained using the Fourier approach for a 1 μm wide ribbon, with an injection energy of 0.1 eV, in the presence of a Lorentzianshaped potential with a half-width of 50 nm and a peak value of 0.2 eV, off-center with respect to the longitudinal axis of the ribbon.

0.0135

0 0.07 |ψA |2 (nm −1)

−0.8 −0.2 −0.15 −0.1 −0.05

|ψA |2 (nm −1)

0.009 0.0045

0.018

0 |ΦAK’|2 (nm −1)

Im {κx } (nm −1)

0.4

0.048 0

200

400 600 y (nm)

800

0 1000

Fig. 4. For the same potential as in Fig. 3 (represented by a dashed line in the bottom panel) and for the eigenvalue with the largest real part: square   K modulus (computed using the Fourier approach) of the ΦK A (y), ΦA (y), and the transverse component of ψA ( r ), as a function of position.

This problem is, somewhat unexpectedly, much more challenging than the solution of the Schr¨odinger equation in a similar domain. A finite difference approach produces a large number of spurious eigenvalues and leads to very poor computational efficiency. On the other hand, an approach based on a modified formulation of the problem and on operating in a Fourier-transformed domain is shown to be effective and suitable for integration into a code for the simulation of transport in graphene nanoribbons. ACKNOWLEDGMENT We thank Claudio Bonati for useful discussions and suggestions. Support from the European Commission through the project GRAND (Contract No. 215752) is gratefully acknowledged. R EFERENCES

In this case the probability density has its largest values at the bottom of the well. These are results that agree with intuition; for different values of the particle injection energy a counterintuitive behavior of the probability density may appear due to the particular properties of the Dirac equation. V. C ONCLUSION We have presented a study of the numerical problem consisting in the solution of the Dirac equation in the transverse section of a graphene armchair nanoribbon.

[1] J. Wurm, M. Wimmer, I. Adagideli, K. Richter, and H. U. Baranger, New J. Phys. 11, 095022 (2009). [2] T. Ando, J. Phys. Soc. Jpn. 74, 777 (2005). [3] L. Brey, H. A. Fertig, Phys. Rev. B 73, 235411 (2006). [4] L. Brugnano, D. Trigiante, Solving Differential Problems by Multistep Initial and Boundary Value Methods, Gordan and Breach, Amsterdam, 1998. [5] B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in Fortran 77: The Art of Scientific Computing, Cambridge University Press, Cambridge, 1992.