PHYSICAL REVIEW B 99, 195406 (2019)

Finite-difference method for Dirac electrons in circular quantum dots

Bartłomiej Szafran, Alina Mreńca-Kolasińska, and Dariusz Żebrowski

AGH University of Science and Technology, Faculty of Physics and Applied Computer Science, Al. Mickiewicza 30, PL-30-059 Kraków, Poland

(Received 29 January 2019; revised manuscript received 15 March 2019; published 3 May 2019)

A simple and reliable finite-difference approach is presented for solution of the Dirac equation eigenproblem
for states confined in rotationally symmetric systems. The method sets the boundary condition for the spinor
wave function components at the external edge of the system and then sweeps the radial mesh in search for
the energies for which the boundary conditions are met inside the flake. The sweep that is performed from the
edge of the system towards the origin allows for application of a two-point finite-difference quotient of the
first derivative, which prevents the fermion doubling problem appearing with spurious solutions and rapidly
oscillating wave functions.

DOI: 10.1103/PhysRevB.99.195406

I. INTRODUCTION be locally introduced to silicene, a graphenelike 2D hexagonal

Si crystal with buckled crystal lattice, using electric fields
Numerical solutions of the Dirac equation for relativistic
vertical to the surface [32,33].
particles have been performed since the 1970s in the context
In this paper we present a very simple and effective finite-
of the problems of nuclear physics and the lattice gauge
difference method for determination of the Dirac Hamiltonian
theories [1]. The interest in computational approaches for the
eigenstates localized in circular quantum dots, applicable to
massless Dirac equation (Weyl equation) has been extended
both the finite flakes and systems with a spatially modulated
to the solid state with the arrival of graphene [2], the two-
energy gap. The method is based on a two-point backward
dimensional (2D) material with gapless and linear dispersion
derivative. Starting from the boundary condition at the exter-
relation near the charge neutrality point.
nal edge of the flake for the two components of the eigen-
Solutions of the discretized version of the Dirac Hamil-
functions one can pin the energies for which the boundary
tonian on a mesh provide both the actual eigenstates and
conditions in the interior of the flake are fulfilled for both
spurious ones that are characterized by rapidly oscillating
the spinor components. The proposed mesh-sweeping method
wave functions, large momenta but low expectation values
resolves only the actual solutions and is free of the spurious
of the energy [1,3–6]. The spurious solutions are degenerate
ones, hence extra terms in the Hamiltonian of a Wilson type
with the regular ones which result in the so-called fermion
or a post-treatment of the eigenstates are not necessary.
doubling problem [1,6–12]. The problem arises in particular
The paper is organized as follows. In Sec. II we present the
with the central, three-point discretization of the first-order
Hamiltonian. The analytical solution given in Sec. III is used
spatial derivative in the Dirac Hamiltonian. The central finite-
for the test calculations. In Sec. IV we illustrate the problem
difference quotient misses the wave function oscillations that
of the spurious solutions, the fermion doubling and the Wil-
occur with the periodicity of the mesh spacing [3,6–8]. The
son procedure using a diagonalization of the finite-difference
spurious solutions are also observed in the finite element
Hamiltonian with the central three-point lattice derivative. The
method [5] despite the exact treatment of shape function
original method reported in this work is presented in Sec. V.
derivatives in this approach. The spurious states can be re-
Examples of applications, including the spatial variation of
moved with the Wilson approach [6,9–11] that introduces
the energy gap and a graphene quantum ring, are presented in
artificial terms in the Hamiltonian depending on the square
Sec. VI. The summary is given in Sec. VII.
of the electron momentum, which removes the spurious states
from the low-energy spectrum. Alternatively, the shift of the
finite-difference quotients can be applied which allows the
lattice derivative to resolve the rapid oscillations of the wave II. HAMILTONIAN
function [3,6–8]. The effective low-energy Dirac Hamiltonian for electrons
In this paper, we focus on confined solutions of the Dirac- in graphene or silicene can be separated into two 2 × 2
type effective Hamiltonians for graphene. The Weyl fermions operators each associated with one of the nonequivalent val-
evade the confinement by the electrostatic potentials with the leys of the Brillouin zone, K or K  . These valleys will be
Klein tunneling phenomenon [13–15]. A way to confine the referred to by the index η = 1 and −1, respectively. The 2 × 2
carriers is to either use a finite flake of graphene [16–28] or Hamiltonian for the valley index η takes the form [2],
introduce a nonzero mass in the region outside of the dot
[29,30]. The nonzero mass along with a finite energy gap is  
introduced to the effective Hamiltonian by the influence of UA (r) h̄v f (kx − iηky )
Hη = , (1)
lattice-matched substrates [31]. The gap or nonzero mass can h̄v f (kx + iηky ) UB (r)

2469-9950/2019/99(19)/195406(8) 195406-1 ©2019 American Physical Society

BARTŁOMIEJ SZAFRAN et al. PHYSICAL REVIEW B 99, 195406 (2019)

where v f is the Fermi velocity, k = −i∇ + h̄e A, and A is the TABLE I. Zeroes of the J|m+η| Bessel function for η = 1, or
vector potential. This Hamiltonian acts on a wave function of energies of states confined within a circular flake of radius R and
 (r) zigzag termination in the units of vRh̄ .
form (r) = ( 1 ), whose components correspond to the f
2 (r)
A and B graphene sublattices [2], respectively. In Eq. (1), n=1 n=2 n=3
UA and UB stand for potentials at the two sublattices. The
potentials can be made unequal due to the substrate of, e.g., m = −2 3.831706 7.015587 10.173468
m = −1 2.404825 5.520078 8.653727
a hexagonal boron nitride [31]. Vertical electric field applied
m=0 3.831706 7.015587 10.173468
to a graphenelike 2D hexagonal crystal with the buckled
m=1 5.135622 8.417244 11.619841
crystal lattice, e.g., the silicene also introduces unequal UA
and UB values, that introduces the energy gap to the dispersion
relation [32,33].
Hamiltonian for a circular flake of graphene in the ver- [16]. For the flake that is terminated by the A sublattice, the
tical magnetic field (0, 0, B) and symmetric gauge A = second component of the wave function that describes the B
(−By/2, Bx/2, 0) commutes with generalized angular mo- sublattice, needs to vanish at r = R, hence f2 ( Evmnf h̄R ) = 0. This
mentum operator Jz = Lz I + η 2h̄ σz , where Lz is the z com- condition defines the spectrum of the states confined within
ponent of the orbital angular momentum. For the polar co- the flake, where n numbers the eigenvalues for fixed m. For the

ordinates, r = x 2 + y2 and φ = arctan(y/x), the operator is conduction band states we use positive values of n. From the
∂
Lz = −i h̄ ∂φ and σz is the Pauli matrix in the sublattice space. boundary condition we find the nonzero energy eigenvalue,
Therefore, the stationary states can be labeled with magnetic Z|m+η|,|n| v f h̄
quantum number m and have the form, Emn = sgn(n) ,
R
 
f1 (r) where Z|m+η|,|n| is the nth zero of J|m+η| Bessel function [34];
m,η = exp(imφ) . (2)
f2 (r) exp(iηφ) see Table I.
The radial functions f1 (r) and f2 (r) of Eq. (2) are solutions to
the system of eigenequations, IV. DIAGONALIZATION OF THE
  FINITE-DIFFERENCE PROBLEM
i h̄  ieBr
UA (r) f1 + vF −η (m + η) f2 − i h̄ f2 − η f2 = E f1 , For illustration of the problem of spurious solutions, we
r 2
  solve Eqs. (4) and (5) in a finite-difference approach using
i h̄ ieBr diagonalization of the resulting algebraic equations. For that
UB (r) f2 + vF η m f1 − i h̄ f1 + η f1 = E f2 .
r 2 purpose we need a boundary condition at the origin. The
radial functions f1 , f2 , near r = 0 behave as r |m| and r |m+η| ,
(3)
respectively. A natural Dirichlet condition at r = 0 indepen-
dent of the orbital angular momentum can be obtained by
III. ANALYTICAL SOLUTION substitution fi = ψr i for i = 1, 2. Since fi is finite at r = 0,
An analytical solution [16,28] to Eq. (3) without the exter- ψi needs to vanish at the origin independent of m. The system
nal fields will be used for the test calculations. For B = 0 and of eigenequations for functions ψi then reads
UA = UB = 0, Eq. (3) reads  
i h̄ηm 
  vF − ψ2 − i h̄ψ2 = E ψ1 , (8)
i h̄ r
vF −η (m + η) f2 − i h̄ f2 = E f1 , (4)  
r i h̄(ηm + 1)
  vF ψ1 − i h̄ψ1 = E ψ2 . (9)
i h̄ r
vF η m f1 − i h̄ f1 = E f2 . (5)
r The boundary conditions are then ψ1 (0) = ψ2 (0) = 0 at
For E = 0 one can eliminate f2 from (5) and plug it into (4) the origin and the zigzag boundary at the edge ψ1 (R) =
to obtain 1, ψ2 (R) = 0. We use N + 1 points on the radial mesh, in-
  cluding the origin r = 0 where the boundary conditions are
2   r2E 2 applied, and the discretization step dr = NR . For discretization
r f1 + r f1 + − m f1 = 0,
2
(6)
v 2f h̄2 we use the difference quotient ψ  = ψ (r+dr)−ψ
2dr
(r−dr)
with the
exception of the last point on the mesh for r = R, where a
which upon introduction of a dimensionless radial coordinate
ρ = h̄v
Er
, produces the Bessel equation, two-point formula is used ψ (R) = ψ (R)−ψdr(R−dr) . The result-
f ing algebraic eigenequation is described by a nonsymmetric
ρ 2 f1 + ρ f1 + (ρ 2 − m2 ) f1 = 0, (7) matrix and is solved with the LAPACK library (procedure
ZGEEV).
where, up to a normalization constant f1 = Jm (ρ), and Jm For the numerical calculations we use v f = 3ta 2h̄
, for the
is the mth Bessel function of the first kind. Using identities silicene tight-binding hopping energy of t = 1.6 eV and the
[34] ρ2 Jm = Jm−1 + Jm+1 , 2Jm = Jm−1 − Jm+1 , one finds f2 = nearest neighbor distance of a = 0.225 nm. For the radius of
iJm+η (ρ). For the test calculations we take a finite flake of the flake we take R = 80 nm. The positive eigenenergies are
radius R and apply a so-called zigzag boundary condition listed in Table II for m = 0 and η = 1, and a varied number

195406-2
FINITE-DIFFERENCE METHOD FOR DIRAC ELECTRONS … PHYSICAL REVIEW B 99, 195406 (2019)

TABLE II. The eigenenergies obtained by diagonalization of the

system of Eqs. (8) and (9) with the finite-difference approach for
m = 0 and η = 1. The eigenenergies of the spurious states that are
put in bold, correspond to the regular solutions but obtained for a
different m; cf. Table I, which results in the fermion doubling.

n = 1 n = 2 n = 3 n = 4 n = 5
N = 200 2.419 3.850 5.548 7.049 8.696
N = 400 2.412 3.841 5.535 7.033 8.675
N = 800 2.409 3.836 5.528 7.024 8.665

of points on the spatial mesh. In Table II we number the

eigenstates by n —the number that counts also the spurious
states—to distinguish from the actual quantum number n used
in Table I and below. For m = 0 and η = 1 one would expect
to obtain the values given by n = 2 and n = 4; see Table I. FIG. 2. Energy levels for η = 1 and m = 0 shifted by the Wilson
The other values are the energies of the spurious states. In term that separates the spurious energy levels from the original ones;
fact, the energies of spurious solutions obtained for m = 0 see Table II. Results were obtained for N = 200 mesh points.
correspond to the regular eigenstates but for a different m.
In particular, the energies for odd n of Table II that were one can numerically remove the spurious states from the low-
calculated with m = 0 correspond to the exact energies for energy part of the spectrum. The procedure [6,9–11] involves
m = −1; see Table I. With the spurious solutions the degener- an extra artificial energy operator of form,
acy of energy levels is artificially increased by a factor of two,
HD = −WD h̄v f ∇ 2 σz dr, (10)
which produces the fermion doubling problem [1,6–12].
The wave functions for n from 1 to 4 are displayed in where WD is a dimensionless Wilson parameter. Figure 2
Fig. 1. With the applied boundary conditions f1 is real and shows the energies of the m = 0, η = 1 states (same as in
f2 is imaginary. The spurious states [Figs. 1(a) and 1(c)] Fig. 1 and Table II), where HD is introduced by the first-order

correspond to the real part of f1 and imaginary part of f2 that perturbation Emn,η = Emn,η + m,n,η |HD |m,n,η . The states
rapidly oscillate between the nearest-neighbor mesh points. identified as spurious in the context of Table II and Fig. 1
The solution given by Figs. 1(b) and 1(d) corresponds to the are indeed removed from the low-energy spectrum, while the
analytical eigenstates. actual solutions are only weekly affected by the Wilson term.
In the spurious solutions k2  is large due to the rapid wave
function oscillations [see Figs. 1(a) and 1(c)]. Using this fact V. THE MESH SWEEPING METHOD
We are ready to introduce the method which is the purpose
of this work. The method provides the solution of the Dirac
equation on a finite-difference mesh and is free of the spurious
(a) (b) states, so no discrimination of the rapidly oscillating states is
necessary. The general system of equations given by Eq. (3) is
transformed to a form that allows for sweeping the mesh from
the edge of the flake to its center. For this purpose the radial
derivative is replaced by a two-point finite-difference quotient
f  (r) = f (r)−dr
f (r−dr)
. The sweeping finite-difference scheme is
(c) (d)

idr f2 (r)
f1 (r − dr) = (UB (r) − E )
h̄ vf
  
i h̄ ieBr i h̄ηm
+ − +η + f1 (r) , (11)
dr 2 r
FIG. 1. Radial functions as obtained from diagonalization of 
the finite-difference version of Eqs. (8) and (9) for m = 0, η = idr f1 (r)
f2 (r − dr) = (UA (r) − E )
1, and N = 200 mesh points, radius of the flake R = 80 nm, h̄ vf
and zigzag boundary condition f2 (r = R) = 0. Central, three-point   
i h̄ ieBr i h̄η(m + η)
finite-difference quotient for the derivative is used with the exception − +η + f2 (r) .
of the last point on a mesh N, where a backward quotient is applied. dr 2 r
The results for the eigenvalues n = 1, 2, 3, 4 of Table II are plotted (12)
in (a), (b), (c), and (d), respectively. (a) and (c) Corresponding to
the spurious solutions, and (b) and (d) to the regular ones. The The schematics of the calculations are given in Fig. 3.
eigenfunctions are determined up to a normalization constant so the The procedure starts at the edge of the flake r = R, where
arbitrary units were used. the boundary conditions are applied. For the zigzag boundary

195406-3
BARTŁOMIEJ SZAFRAN et al. PHYSICAL REVIEW B 99, 195406 (2019)

Origin Boundary

FIG. 3. The schematics of the sweeping method. The boundary

(a)
conditions f1 (R), f2 (R) are applied at the external edge of the flake.
Equations (11) and (12) allow for evaluation of the wave function
components to the left for the energy E which is a parameter to be
determined by the boundary condition to be reproduced at the origin.

condition we set f2 (R) = 0 and f1 (R) = 1. The eigenfunc-

tions can be normalized after the entire procedure. With the
starting values of f1 and f2 at r = R one can proceed to
the center of the flake using Eqs. (11) and (12). At the
origin f1 and f2 need to vanish when m = 0 and m + η = 0,
respectively. That, in turn is realized only for discrete values
of the energy. For η = 1, the values of the components of the
wave function for r = 0 are given in Fig. 4 for m = 0 (left (b)
column) and m = 1 (right column). For m = 0 the component
f2 needs to vanish at the origin, since the angular momentum
quantum number corresponding to the second component is
nonzero, m + η = 1. The position of the energy eigenvalues
can be found as the minima of the absolute values of the
f2 (r = 0) [Fig. 4(a)] or equivalently, by zeros of the imaginary
part of f2 (r = 0) [see Fig. 4(c)] as a function of the energy.
The minima [Fig. 4(a)] and the zeros [Fig. 4(c)] correspond to
the actual eigenvalues and the spurious solutions are missing
(cf. Tables I and II). For m = 1 both f1 and f2 components
need to vanish at the origin, and indeed no shift of the minima
of Fig. 4(b) or zeros of Fig. 4(d) is observed on the energy
scale. Besides the solutions of the Bessel form, the results of (c)
Fig. 4 indicate the presence of the zero-energy levels that are
supported by the zigzag edge [16].
The success of the sweeping approach is due to the two-
point derivative that does not allow for sawlike oscillations of
the wave function. To explain this let us look at Fig. 5 where
two functions—a smooth one f and a rapidly oscillating one
g—are plotted on a radial mesh. For even indexes functions
f and g are equal and for odd indexes a shift between the
solutions is present. From the point of view of the central
quotient of the first derivative both the solutions are identical;
in each mesh point the same value of the quotient is obtained,
and there is no relation between the odd and even indexes
on the mesh points, which allows for the spurious states to
appear in the low-energy part of the spectrum. The Wilson (d)
term given by Eq. (10) introduces the second derivative with
FIG. 4. Absolute values (a) and (b), real and imaginary parts of
the quotient f  (i)  f (i+1)+ fdr(i−1)−2
2
f (i)
that links the odd and
the radial functions f1 and f2 (c) and (d) at the origin for the sweeping
even points and removes the rapidly oscillating states from method starting at the edge of the flake for the K valley η = 1.
the low-energy spectrum. The sweeping method that uses the Panels (a) and (c) were calculated for m = 0 and (b) and (d) for
backward quotient f  (i) = f (i)−drf (i−1) and passes through each m = 1.
mesh point consecutively introduces the link between the even
and odd mesh points. The two-point formula for the smooth

195406-4
FINITE-DIFFERENCE METHOD FOR DIRAC ELECTRONS … PHYSICAL REVIEW B 99, 195406 (2019)

Re(f1), Im(f2) [arb.unit.]

Re(f1), Im(f2) [arb.unit.]

TABLE III. First seven rows show the three lowest positive
energy levels for m = 0 and η = 1 in the units of vRh̄ as calculated f2
f
1
with the sweeping method for UA = UB = 0 and B = 0 as a function
of the number of mesh points N. The central row gives the exact 1
results (cf. Table I). Bottom part of the table, row 10 and below, show f1 f2 f1
the difference of the numerical and exact results as a function of N.
0 0
N E (n = 1) E (n = 2) E (n = 3) 0 20 40 60 80 0 20 40 60 80
(a) r[nm] (b) r[nm]
100 3.853094 7.081613 10.304313

Re(f1), Im(f2) [arb.unit.]

Re(f1), Im(f2) [arb.unit.]

200 3.842100 7.046336 10.231384 1 1
400 3.836828 7.030400 10.200578
800 3.8342475 7.022854 10.186564 0
0 f2 f1
1600 3.832970 7.019184 10.179904 f2 f1
-1
3200 3.832339 7.017377 10.176656
-1
6400 3.832023 7.016478 10.175054 -2
Exact 3.831706 7.015587 10.173468
N E (n = 1) E (n = 2) E (n = 3) 0 20 40 60 80 0 20 40 60 80

100 0.021388 0.066026 0.130845 (c) r[nm] (d) r[nm]

R e(f1), Im(f2) [arb.unit.]

200 0.010394 0.030749 0.057916

Re(f1), Im(f2) [arb.unit.]

400 0.005122 0.014813 0.027110 2 f1
3
f2 f1
800 0.002541 0.007267 0.013096
1 2 f2
1600 0.001264 0.003597 0.006436
3200 0.000633 0.00179 0.003188 1
0
6400 0.000317 0.000891 0.001586 0
-1 -1
0 20 40 60 80 0 20 40 60 80
(e) r[nm] (f) r[nm]
and rapidly oscillating solutions produce very different values
of the quotient. FIG. 6. Real and imaginary parts of the radial functions f1 and
Tabel III shows that the convergence of the results of the f2 , respectively, for m = −1 (a), (c), and (e), and for m = 0 (b),
present method to the exact solution is linear as a function of (d), and (f) for the sweeping method starting at the edge of the
dr (or 1/N). Although the convergence is slow, the numerical flake for the K valley η = 1. In (a) and (b) the states are n = 1, in
cost of the calculation increases only linearly with N, and the (c) and (d) n = 2, and in (e) and (f) n = 3. The lines that give the
method needs to keep track of only four complex numbers for analytical solutions are given by f1 = Jm (ρ ), f2 = iJm+η (ρ ), with
the wave function components at adjacent mesh points, so one ρ = h̄v Er
[16,28], with E standing for the energy level of the given
f
can approach the exact solution arbitrarily close at a negligible state. The dots mark the numerical result.
numerical cost.
For completeness, as an additional verification of the
method, Fig. 6 shows the comparison of the analytically cal- The present approach is a simple alternative to the method
culated radial parts of the wave functions, and the numerically proposed in Ref. [35], where an additional finite-difference
obtained ones for 2000 nodes. mesh is introduced to the center of the original 2D mesh.
Next, the values of the wave functions in the additional mesh
are calculated as an arithmetic average of the neighbor lattice
points. A similar approach was proposed in Ref. [8]. The
present approach is limited to the circular symmetry but works
with the original mesh only.

VI. EXAMPLES OF APPLICATIONS

The approach is suitable for any radial problems involving
the effective Hamiltonian that stems from the Dirac equation
for a given valley. This section provides examples of applica-
FIG. 5. A schematic drawing of a smooth solution f given on a tions.
radial mesh, and a spurious one g, rapidly oscillating from one lattice
point to another. The central quotient for the derivative for point (4)
A. Limit of infinite-mass confinement
can be equal for both solutions f  (4)  f (5)− 2dr
f (3)
= g(5)−g(3)
2dr
 g (5)
(see the black dotted lines). The present sweeping method is based In the precedent section the zigzag-boundary conditions
on an asymmetric formula for which the quotients are very different were used that result from a specific termination of the crystal
for both functions f (5)−
dr
f (4)
= g(5)−g(4)
dr
, so that the spurious solution lattice. The confinement of the particle described by the
is associated with a very different energy. Dirac Hamiltonian can be induced provided by, e.g., a spatial

195406-5
BARTŁOMIEJ SZAFRAN et al. PHYSICAL REVIEW B 99, 195406 (2019)

Re(f1), Im(f2) [arb.unit.]

Re(f1), Im(f2) [arb.unit.]

2 20
f1
1.5 15
1 f1
10
0.5
5 f2
0
f2 0
-0.5
-1 -5
-1.5 -10
40 50 60 70 80 40 50 60 70 80
r[nm] r[nm]
(a) (b)

Re(f1), Im(f2) [arb.unit.]

Re(f1), Im(f2) [arb.unit.]

1.5 2
1 1.5
0.5 1
(a) 0.5 f2 f
0 1
f1 0
-0.5
f2 -0.5
-1 -1
-1.5 -1.5
-2 -2
40 50 60 70 80 40 50 60 70 80
r[nm] r[nm]
(c) (d)

FIG. 8. The radial functions with m = 0 for the K valley (η = 1)

of the (a) and (b) n = +1, (c) and (d) n=+2 state in a ring with
inner (outer) radius Ri = 40 nm (R = 80 nm) and the infinite mass
boundary conditions and B = 0 (a) and (c), B = 5 T (b) and (d).

(b)
rent must vanish (cos φ, sin φ) · j = 0, which implies tan φ =
Re( ∗  )
−η Im(1∗ 22 ) . With Eq. (2) this condition is translated for the
1
Re( f ∗ (R) f (R) exp(iηφ))
radial functions as tan φ = −η Im( f1∗ (R) f22 (R) exp(iηφ)) . To fulfill
1
this condition we set at the end of the flake f1 (R) = 1 and
f2 (R) = i, for which one gets tan φ = η tan ηφ, which is
fulfilled for both the valleys η = ±1.
Let us look at the appearance of the wave function confine-
ment by the mass boundary. We consider a spatial dependence
of the energy gap as introduced by finite UA = −UB potentials.
Figure 7 shows the real part of f1 and the imaginary part
of f2 for η = 1 and m = 0. At the end of the flake the
(c) infinite mass boundary condition is applied. The components
f1 and f2 are—as for the zigzag boundary—purely real and
FIG. 7. Real and imaginary parts of the f1 and f2 radial func- purely imaginary, respectively. Moreover, the infinite-mass
tions, respectively, for an infinite mass boundary condition at the boundary condition implies Re{ f1 (r = R)} = Im{ f2 (r = R)}.
edge of the flake R = 80 nm (a) and R = 100 nm (b) and (c). We The results of Fig. 7(a) were obtained for a flake of radius R =
consider η = 1 state with m = 0. In (a) no external potential is 80 nm in the absence of the external potentials. In Figs. 7(b)
applied. In (b) and (c) the external potential is introduced UA = −UB
with UA = 0.1 eV (b) and UA = 1 eV (c).
(a) (b)
0.1 K K'
dependence of the energy gap [30] of equivalently by an 0.05
external potential due to, e.g., a substrate [31].
E [eV]

For UA = UB an energy gap is opened in the dispersion 0

relation near the Dirac points and the carriers acquire finite
-0.05
masses. For UA and UB that diverge to ∞ and −∞ when
r > R, respectively, the mass in the outside of the confine- -0.1
ment is infinite and the wave functions of finite energies
vanish at r > R [29,35,36]. The boundary condition for this 0 2 4
B [T]
6 8 10 0 2 4
B [T]
6 8 10

type of confinement is derived from vanishing current at the

edge [29] of the confinement area. The probability density FIG. 9. The energy levels of a ring with inner (outer) radius Ri =
current that is derived from the Hamiltonian [29] is j = 40 nm (R = 80 nm) and with infinite mass boundary conditions, with
2vF [Re(1∗ 2 ), ηIm(1∗ 2 )]. For the infinite mass at the m = 0, ±1, ..., ±4. The black curve shows the m = −4 level, and the
outside of the dot, the orthogonal component of the cur- brighter curves have higher m.

195406-6
FINITE-DIFFERENCE METHOD FOR DIRAC ELECTRONS … PHYSICAL REVIEW B 99, 195406 (2019)

(a) IM (b) ZZ of the radius of R = 80 nm we remove the central disk of a

20 radius of Ri = 40 nm. The infinite mass boundary condition
in the inner edge of the flake reads f2 (Ri )/ f1 (Ri ) = −i. We
15 start from the external edge R as in the precedent subsection.
The sweep at the finite-difference mesh stops at Ri . From
the inner boundary condition we construct a function F (E ) =
10 f2 (Ri ; E ) + i f1 (Ri ; E ) and we look for its zeros as a function
of the energy.
5 The m = 0 wave functions for the K valley (η = 1) found
in this way are displayed for n = 1 in Figs. 8(a) and 8(b)
0 and for n = 2 in Figs. 8(c) and 8(d). The energy spectra
for η = ±1 are given in Fig. 9 for |m|  4. In the lowest
energy states of the conduction and the valence bands one ob-
-5
serves the characteristic angular momentum transitions with
the Aharonov-Bohm periodicity of B  0.385 T [37], that
-10 corresponds to a flux quantum threading a ring of an effective
radius 59 nm—close to the central radius of the ring Rc =
-15 (R + Ri )/2.

-20 C. Comparison with analytical results in the

0 5 10 15 20 25 0 5 10 15 20 25 external magnetic field
In the absence of the external potential UA = UB = 0, or
for a massless particle, analytical results were obtained for
FIG. 10. The energy spectrum calculated for graphene quantum the Weyl equation in the external magnetic field in Ref. [16]
dot of radius R = 70 nm with the present approach, to be compared for both zigzag and infinite mass boundary conditions. We
with Fig. 1 of Ref. [16]. In (a) the infinite mass boundary condition is adopted the graphene parameters of this work [16] t = 2.7 eV,
applied and in (b) the zigzag boundary. The dashed lines correspond the lattice constant 0.142 nm, and the radius of the dot R = 70
to K  valley and the solid ones with the K valley. The energies and nm. The results for the 500 mesh points are displayed in
the magnetic field are expressed in dimensionless units,
= hv ER
and
F Fig. 10 with a perfect agreement with Fig. 1 of Ref. [16].
eBR2
β= 2h̄
. Energy levels for m ∈ [−4, 4] are displayed.
VII. SUMMARY

and 7(c) we extended the flake to R = 100 nm but introduced We have presented the finite-difference method, applicable
the energy gap beyond r > R = 80 nm with UA = −UB . In to Dirac carriers in circular confinement, that sweeps the mesh
Fig. 7(a), UA = 0.1 eV and UA = 1 eV in Fig. 7(c). We can see from the external edge of the system, where the boundary con-
that for larger UA the functions penetrate only weakly into the ditions are defined, to the center. The energies of the confined
gapped region for r > R = 80 nm, and the values Re{ f1 } and states are pinned by the internal boundary conditions. The
Im{ f2 } become equal at r = R , where the gap is introduced. method is simple, free of spurious solutions and the fermion
In this way the infinite boundary condition is found in the limit doubling problem. Applications to circular flakes, quantum
of large UA at r = R . rings, and the energy gap modulated in space were presented.

B. The quantum ring spectra ACKNOWLEDGMENT

The method with a slight modification can be applied to This work was supported by the National Science Centre
a problem of a graphene quantum ring [28]. From the flake (NCN) according to Grant No. DEC-2016/23/B/ST3/00821.

[1] J. P. Kogut, Rev. Mod. Phys. 55, 775 (1983). [8] J. Tworzydło, C. W. Groth, and C. W. J. Beenakker, Phys. Rev.
[2] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, B 78, 235438 (2008).
and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009). [9] Y. Tanimura, K. Hagino, and H. Z. Liang, Prog. Theor. Exp.
[3] C. M. Bender, K. A. Milton, and D. H. Sharp, Phys. Rev. Lett. Phys. 2015, 073D01 (2015).
51, 1815 (1983). [10] K. G. Wilson, Phys. Rev. D 10, 2445 (1974).
[4] J. C. Wells, V. E. Oberacker, A. S. Umar, C. Bottcher, M. R. [11] L. Susskind, Phys. Rev. D 16, 3031 (1977).
Strayer, J.-S. Wu, and G. Plunien, Phys. Rev. A 45, 6296 (1992). [12] F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk, Comput.
[5] W. Poeschl, D. Vretenar, and P. Ring, Comput. Phys. Commun. Phys. Commun. 183, 1403 (2012).
99, 128 (1996). [13] A. V. Rozhkov, G. Giavaras, Y. P. Bliokh, V. Freilikher, and F.
[6] R. Stacey, Phys. Rev. D 26, 468 (1982). Nori, Phys. Rep. 503, 77 (2011).
[7] C. Bottcher, M. R. Strayer, A. S. Umar, and P.-G. Reinhard, [14] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys.
Phys. Rev. A 40, 4182 (1989). 2, 620 (2006).

195406-7
BARTŁOMIEJ SZAFRAN et al. PHYSICAL REVIEW B 99, 195406 (2019)

[15] A. F. Young and P. Kim, Nat. Phys. 5, 222 (2009). [27] P. Potasz, A. D. Güçlü, A. Wöjs, and P. Hawrylak, Phys. Rev. B
[16] M. Grujic, M. Zarenia, A. Chaves, M. Tadic, G. A. Farias, and 85, 075431 (2012).
F. M. Peeters, Phys. Rev. B 84, 205441 (2011). [28] M. R. Thomsen and T. G. Pedersen, Phys. Rev. B 95, 235427
[17] M. Zarenia, A. Chaves, G. A. Farias, and F. M. Peeters, Phys. (2017).
Rev. B 84, 245403 (2011). [29] M. V. Berry, and R. J. Mondragon, Proc. R. Soc. London A 412,
[18] J. Fernandez-Rossier and J. J. Palacios, Phys. Rev. Lett. 99, 53 (1987).
177204 (2007). [30] G. Giavaras and F. Nori, Phys. Rev. B 83, 165427 (2011).
[19] B. Wunsch, T. Stauber, and F. Guinea, Phys. Rev. B 77, 035316 [31] B. Sachs, T. O. Wehling, M. I. Katsnelson, and A. I.
(2008). Lichtenstein, Phys. Rev. B 84, 195414 (2011).
[20] P. Potasz, A. D. Güçlü, and P. Hawrylak, Phys. Rev. B 81, [32] Z. Ni, Q. Liu, K. Tang, J. Zheng, J. Zhou, R. Qin, Z. Gao, D.
033403 (2010). Yu, and J. Lu, Nano Lett. 12, 113 (2012).
[21] W. L. Wang, O. V. Yazyev, S. Meng, and E. Kaxiras, Phys. Rev. [33] N. D. Drummond, V. Zolyomi, and V. I. Fal’ko, Phys. Rev. B
Lett. 102, 157201 (2009). 85, 075423 (2012).
[22] S. Moriyama, Y. Morita, E. Watanabe, and D. Tsuya, Appl. [34] M. Abramowitz and I. A. Stegun (eds.), in Handbook of Math-
Phys. Lett. 104, 053108 (2014). ematical Functions with Formulas, Graphs, and Mathematical
[23] T. Yamamoto, T. Noguchi, and K. Watanabe, Phys. Rev. B 74, Tables, 10th printing (Dover, New York, 1972), pp. 355–389.
121409(R) (2006). [35] X. Ni, L. Huang, Y.-C. Lai, and C. Grebogi, Phys. Rev. E 86,
[24] Z. Z. Zhang, K. Chang, and F. M. Peeters, Phys. Rev. B 77, 016702 (2012).
235411 (2008). [36] L. Huand, H.-Y. Xu, C. Grebogi, and Y.-C. Lai, Phys. Rep. 753,
[25] A. D. Güçlü, P. Potasz, and P. Hawrylak, Phys. Rev. B 82, 1 (2018).
155445 (2010). [37] S. Viefers, P. Koskinen, P. Singha Deo, and M. Manninen,
[26] M. Ezawa, Phys. Rev. B 77, 155411 (2008). Physica E 21, 1 (2004).

195406-8