Finite difference method for Dirac electrons in circular quantum dots

TL;DR

This paper introduces a finite difference method for solving the Dirac equation in circular quantum dots, effectively handling boundary conditions and avoiding fermion doubling issues.

Contribution

It presents a simple, reliable finite difference approach that sets boundary conditions at the system edge and sweeps inward, preventing spurious solutions.

Findings

Method effectively solves Dirac eigenproblems in circular geometries.

Prevents fermion doubling and spurious oscillations in solutions.

Applicable to rotationally symmetric quantum systems.

Abstract

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 and the appearance of the spurious solutions with rapid oscillations of the wave functions in space.

Figures

Tables

Table 1. Table 1: Zeroes of the J | m + η | subscript 𝐽 𝑚 𝜂 J_{|m+\eta|} Bessel function for η = 1 𝜂 1 \eta=1 , or energies of states confined within a circular flake of radius R 𝑅 R and zigzag termination in the units of R v f ​ ℏ 𝑅 subscript 𝑣 𝑓 Planck-constant-over-2-pi \frac{R}{v_{f}\hbar} .

	$n=1$	$n=2$	$n=3$
$m=-2$	3.831706	7.015587	10.173468
$m=-1$	2.404825	5.520078	8.653727
$m=0$	3.831706	7.015587	10.173468
$m=1$	5.135622	8.417244	11.619841

Table 2. Table 2: The eigenergies obtained by diagonalization of the system of Eqs. ( 8 ) and ( 9 ) with the finite difference approach for m = 0 𝑚 0 m=0 and η = 1 𝜂 1 \eta=1 . The eigenergies of the spurious states that are put in bold, correspond to the regular solutions but obtained for a different m 𝑚 m – cf. Table 1 , which results in the fermion doubling.

	$n^{\prime }=1$	$n^{\prime }=2$	$n^{\prime }=3$	$n^{\prime }=4$	$n^{\prime }=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

Table 3. Table 3: First seven rows show the three lowest positive energy levels for m = 0 𝑚 0 m=0 and η = 1 𝜂 1 \eta=1 in the units of R v f ​ ℏ 𝑅 subscript 𝑣 𝑓 Planck-constant-over-2-pi \frac{R}{v_{f}\hbar} as calculated with the sweeping method for U A = U B = 0 subscript 𝑈 𝐴 subscript 𝑈 𝐵 0 U_{A}=U_{B}=0 and B = 0 𝐵 0 B=0 as a function of the number of mesh points N 𝑁 N . The central row gives the exact results (cf. Table 1 ). Bottom part of the table: rows 10th and below show the difference of the numerical and exact results as function of N 𝑁 N .

$N$	$E\left(n=1\right)$	$E\left(n=2\right)$	$E\left(n=3\right)$
100	3.853094	7.081613	10.304313
200	3.842100	7.046336	10.231384
400	3.836828	7.030400	10.200578
800	3.8342475	7.022854	10.186564
1600	3.832970	7.019184	10.179904
3200	3.832339	7.017377	10.176656
6400	3.832023	7.016478	10.175054
exact	3.831706	7.015587	10.173468
$N$	$\Delta E\left(n=1\right)$	$\Delta E\left(n=2\right)$	$\Delta E\left(n=3\right)$
100	0.021388	0.066026	0.130845
200	0.010394	0.030749	0.057916
400	0.005122	0.014813	0.027110
800	0.002541	0.007267	0.013096
1600	0.001264	0.003597	0.006436
3200	0.000633	0.00179	0.003188
6400	0.000317	0.000891	0.001586

Equations

$$H_{\eta}=\left(\begin{array}[]{cc}U_{A}({\bf r})&\hbar v_{f}(k_{x}-i\eta k_{y})\\ \hbar v_{f}(k_{x}+i\eta k_{y})&U_{B}({\bf r})\end{array}\right),$$

$$\Psi_{m,\eta}=\exp(im\phi)\left(\begin{array}[]{c}f_{1}(r)\\ f_{2}(r)\exp(i\eta\phi)\end{array}\right).$$

$$\left\{\begin{array}[]{c}U_{A}(r)f_{1}+v_{F}\left[-\eta\frac{i\hbar}{r}(m+\eta)f_{2}-i\hbar f^{\prime}_{2}-\eta\frac{ieBr}{2}f_{2}\right]=Ef_{1},\\ U_{B}(r)f_{2}+v_{F}\left[\eta\frac{i\hbar}{r}mf_{1}-i\hbar f^{\prime}_{1}+\eta\frac{ieBr}{2}f_{1}\right]=Ef_{2}.\end{array}\right.$$

$$v_F\left[-\eta \frac{i\hslash }{r}(m+\eta )f_2-i\hslash f_2^{\prime }\right]=E{f}_1,$$

$$v_F\left[\eta \frac{i\hslash }{r}m{f}_1-i\hslash f_1^{\prime }\right]=E{f}_2.$$

$$r^2{f}_1^{\prime \prime }+r{f}_1^{\prime }+\left(\frac{r^2{E}^2}{v_f^2{\hslash }^2}-m^2\right)f_1=0,$$

$${\rho }^2{f}_1^{\prime \prime }+\rho f_1^{\prime }+({\rho }^2-m^2)f_1=0,$$

$$E_{mn}=\mathrm{sign}(n)\frac{Z_{|m+\eta |,|n|}{v}_f\hslash }{R}$$

$$v_F\left[-\frac{i\hslash \eta m}{r}{\psi }_2-i\hslash {\psi }_2^{\prime }\right]=E{\psi }_1,$$

$$v_F\left[\frac{i\hslash (\eta m+1)}{r}{\psi }_1-i\hslash {\psi }_1^{\prime }\right]=E{\psi }_2.$$

$$H_D=-W_D\hslash v_f{\nabla }^2{\sigma }_{z}dr,$$

$$f_1(r-dr)$$

$$f_2(r-dr)$$

Full text

I Introduction

Numerical solutions of the Dirac equation for relativistic particles have been performed since the 1970’s in the context of the problems of nuclear physics and the lattice gauge theories kogut . The interest in computational approaches for the massless Dirac equation (Weyl equation) has been extended to the solid state with the arrival of graphene neto , the two-dimensional material with gapless and linear dispersion relation near the charge neutrality point.

Solution of the discretized version of the Dirac Hamiltonian on a mesh is commonly pested with spurious solutions that are characterized by rapidly oscillating wave functions, large momenta but low expectation values of the energy kogut ; stacey ; bender ; fermdi ; poszl . The spurious solutions are degenerate with the regular ones which results in the so-called fermion doubling problem kogut ; stacey ; tanimura ; bottcher ; wilson ; susskind ; tworzydlo ; spl4 . The problem arises in particular with the central, three-point discretization of the first order spatial derivative in the Dirac Hamiltonian. The central finite difference quotient misses the wave function oscillations that occur with the periodicity of the mesh spacing stacey ; bender ; bottcher ; tworzydlo . The spurious solutions are also observed in the finite element method poszl despite the exact treatment of shape functions derivatives in this approach. The spurious states can be removed with Wilson approach stacey ; tanimura ; wilson ; susskind that introduces artificial terms in the Hamiltonian depending on the square 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 function stacey ; bottcher ; bender ; tworzydlo .

In this paper, we focus on confined solutions of the Dirac type effective Hamiltonians for graphene. The Weyl fermions evade the confinement by the electrostatic potentials with the Klein tunneling phenomenon k1 ; k2 ; k3 . A way to confine the carriers is to either use a finite flake of graphene gru ; zarenia ; flake1 ; flake2 ; flake4 ; flake5 ; mori ; Yamamoto ; Zhang ; Guclu ; Ezawaf ; Potasz1 ; tho or introduce a non-zero mass in the region outside of the dot berry ; nori . The non-zero mass along with a finite energy gap is introduced to the effective Hamiltonian by the influence of lattice-matched substrates sachs . The gap or non-zero mass can be locally introduced to silicene, a graphene-like 2D hexagonal Si crystal with buckled crystal lattice, using electric fields vertical to the surface ni ; Drummond12 .

In this paper we present a very simple and effective finite difference method for determination of the Dirac Hamiltonian eigenstates localized in circular quantum dots, applicable to both the finite flakes and systems with spatially modulated energy gap. The method is based on a two-point backward derivative. Starting from the boundary condition at the external edge of the flake for the two components of the eigenfunctions one can pin the energies for which the boundary conditions in the interior of the flake are fulfilled for both the spinor components. The proposed mesh-sweeping method resolves only the actual solutions and is free of the spurious ones, hence extra terms in the Hamiltonian of a Wilson type or a post-treatment of the eigenstates are not necessary.

The paper is organized as follows. In Section II we present the Hamiltonian. The analytical solution given in Section III is used for the test calculations. In Section IV we illustrate the problem of the spurious solutions, the fermion doubling and the Wilson procedure using a diagonalization of the finite difference Hamiltonian with the central three point lattice derivative. The original method reported in this work is presented in Section V. Examples of applications, including the spatial variation of the energy gap and a graphene quantum ring are presented in Section VI. The summary is given in Section VII.

II Hamiltonian

The effective low-energy Dirac Hamiltonian for electrons in graphene or silicene can be separated into two $2\times 2$ operators each associated with one of the non-equivalent valleys of the Brillouin zone, $K$ or $K^{\prime}$. These valleys will be referred to by the index $\eta=1$ and $-1$, respectively. The $2\times 2$ Hamiltonian for the valley index $\eta$ takes the form neto

[TABLE]

where $v_{f}$ is the Fermi velocity, ${\bf k}=-i\nabla+\frac{e}{\hbar}{\bf A}$, and ${\bf A}$ is the vector potential. This Hamiltonian acts on a wave function of form \Psi({\bf r})=\left(\begin{array}[]{c}\Psi_{1}({\bf r})\\ \Psi_{2}({\bf r})\end{array}\right), whose components correspond to the A and B graphene sublattices neto , respectively. In Eq. (1), $U_{A}$ and $U_{B}$ stand for potentials at the two sublattices. The potentials can be made unequal due to the substrate of e.g. a hexagonal boron nitride sachs . Vertical electric field applied to a graphene-like 2D hexagonal crystal with the buckled crystal lattice, e.g. the silicene also introduces unequal $U_{A}$ and $U_{B}$ values, that introduces the energy gap to the dispersion relation ni ; Drummond12 .

Hamiltonian for a circular flake of graphene in the vertical magnetic field $(0,0,B)$ and symmetric gauge ${\bf A}=(-By/2,Bx/2,0)$ commutes with generalized angular momentum operator $J_{z}=L_{z}{\bf I}+\eta\frac{\hbar}{2}\sigma_{z}$, where $L_{z}$ is the z-component of the orbital angular momentum. For the polar coordinates: $r=\sqrt{x^{2}+y^{2}}$ and $\phi=\arctan(y/x)$, the operator is $L_{z}=-i\hbar\frac{\partial}{\partial\phi}$ and $\sigma_{z}$ is the Pauli matrix in the sublattice space. Therefore, the stationary states can be labeled with magnetic quantum number $m$ and have form

[TABLE]

The radial functions $f_{1}(r)$ and $f_{2}(r)$ of Eq. (2) are solutions to the system of eigenequations

[TABLE]

III Analytical solution

An analytical solution gru ; tho to Eq. (3) without the external fields will be used for the test calculations. For $B=0$ and $U_{A}=U_{B}=0$, Eq. (3) reads

[TABLE]

For $E\neq 0$ one can eliminate $f_{2}$ from (5) and plug it into (4) to obtain

[TABLE]

which upon introduction of a dimensionless radial coordinate $\rho=\frac{Er}{\hbar v_{f}}$, produces the Bessel equation

[TABLE]

where, up to a normalization constant $f_{1}=J_{m}(\rho)$, and $J_{m}$ is the $m$th Bessel function of the first kind. Using identities bee $\frac{2}{\rho}J_{m}=J_{m-1}+J_{m+1}$, $2J_{m}^{\prime}=J_{m-1}-J_{m+1}$, one finds $f_{2}=iJ_{m+\eta}(\rho)$. For the test calculations we take a finite flake of radius $R$ and apply a so-called zigzag boundary condition gru . For the flake that is terminated by the $A$ sublattice, the second component of the wave function – that describes the $B$ sublattice, needs to vanish at $r=R$, hence $f_{2}(\frac{E_{mn}R}{v_{f}\hbar})=0$. This condition defines the spectrum of the states confined within the flake, where $n$ numbers the eigenvalues for fixed $m$. For the conduction band states we use positive values of $n$. From the boundary condition we find the non-zero energy eigenvalue

[TABLE]

where $Z_{|{m+\eta}|,|n|}$ is the $n$th zero of $J_{|m+\eta|}$ Bessel function bee – see Table 1.

IV Diagonalization of the finite difference problem

For illustration of the problem of spurious solutions, we solve Eqs. (4,5) in a finite difference approach using diagonalization of the resulting algebraic equations. For that purpose we need a boundary condition at the origin. The radial functions $f_{1}$, $f_{2}$, near $r=0$ behave as $r^{|m|}$ and $r^{|m+\eta|}$, respectively. A natural Dirichlet condition at $r=0$ independent of the orbital angular momentum can be obtained by substitution $f_{i}=\frac{\psi_{i}}{r}$ for $i=1,2$. Since $f_{i}$ is finite at $r=0$, $\psi_{i}$ needs to vanish at the origin independent of $m$. The system of eigenequations for functions $\psi_{i}$ then reads

[TABLE]

The boundary conditions are then $\psi_{1}(0)=\psi_{2}(0)=0$ at the origin and the zigzag boundary at the edge $\psi_{1}(R)=1$, $\psi_{2}(R)=0$. We use $N+1$ points on the radial mesh, including the origin $r=0$ where the boundary conditions are applied, and the discretization step $dr=\frac{R}{N}$. For discretization we use the difference quotient $\psi^{\prime}=\frac{\psi(r+dr)-\psi(r-dr)}{2dr}$ with the exception of the last point on the mesh for $r=R$, where a two point formula is used $\psi(R)=\frac{\psi(R)-\psi(R-dr)}{dr}$. The resulting algebraic eigenequation is described by a non-symmetric matrix and is solved with the ARPACK library (procedure ZGEEV).

For the numerical calculations we use $v_{f}=\frac{3ta}{2\hbar}$, for the silicene tight-binding hopping energy of $t=1.6$ eV and the nearest neighbour distance of $a=0.225$ nm. For the radius of the flake we take $R=80$ nm. The positive eigenenergies are listed in Table 2 for $m=0$, and $\eta=1$ and a varied number of points on the spatial mesh. In Table 2 we number the eigenstates by $n^{\prime}$ – the number that counts also the spurious states – to distinguish from the actual quantum number $n$ used in Table 1 and below. For $m=0$ and $\eta=1$ one would expect to obtain the values given by $n^{\prime}=2$ and $n^{\prime}=4$ – see Table 1. The other values are the energies of the spurious states. In fact, the energies of spurious solutions obtained for $m=0$ correspond to the regular eigenstates but for a different $m$. In particular, the energies for odd $n^{\prime}$ of Table 2 that were calculated with $m=0$ correspond to the exact energies for $m=-1$ – see Table 1. With the spurious solutions the degeneracy of energy levels is artificially increased by a factor of two, which produces the fermion doubling problem kogut ; stacey ; tanimura ; bottcher ; wilson ; susskind ; tworzydlo ; spl4 .

The wave functions for $n^{\prime}$ from 1 to $4$ are displayed in Fig. 1. With the applied boundary conditions $f_{1}$ is real and $f_{2}$ is imaginary. The spurious states [Fig. 1(a,c)] correspond to the real part of $f_{1}$ and imaginary part of $f_{2}$ that rapidly oscillate between the nearest-neighbor mesh points. The solution given by Fig. 1(b,d) corresponds to the analytical eigenstates.

In the spurious solutions $\langle{\bf k}^{2}\rangle$ is large due to the rapid wave function oscillations (see Fig. 1(a,c)). Using this fact one can numerically remove the spurious states from the low-energy part of the spectrum. The procedure stacey ; tanimura ; wilson ; susskind involves an extra artificial energy operator of form

[TABLE]

where $W_{D}$ is a dimensionless Wilson parameter. Figure 2 shows the energies of the $m=0$, $\eta=1$ states (same as in Fig. 1 and Table 2), where $H_{D}$ is introduced by the first order perturbation $E_{mn,\eta}^{\prime}=E_{mn,\eta}+\langle\Psi_{m,n,\eta}|H_{D}|\Psi_{m,n,\eta}\rangle$. The states identified as spurious in the context of Table 2 and Fig. 1 are indeed removed from the low-energy spectrum, while the actual solutions are only weekly affected by the Wilson term.

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 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^{\prime}(r)=\frac{f(r)-f(r-dr)}{dr}$. The sweeping finite difference scheme is

[TABLE]

The schematics of the calculations given in Fig. 3. The procedure starts at the edge of the flake $r=R$, where the boundary conditions are applied. For the zigzag boundary condition we set $f_{2}(R)=0$ and $f_{1}(R)=1$. The eigenfunctions can be normalized after the entire procedure. With the starting values of $f_{1}$ and $f_{2}$ at $r=R$ one can proceed to the center of the flake using Eqs. (11) and (12). At the origin $f_{1}$ and $f_{2}$ need to vanish when $m\neq 0$ and $m+\eta\neq 0$, respectively. That, in turn is realized only for discrete values of the energy. For $\eta=1$, the values of the components of the wave function for $r=0$ are given in Fig. 4 for $m=0$ (left column) and $m=1$ (right column). For $m=0$ the component $f_{2}$ needs to vanish at the origin, since the angular momentum quantum number corresponding to the second component is nonzero, $m+\eta=1$. The position of the energy eigenvalues can be found as the minima of the absolute values of the $f_{2}(r=0)$ [Fig. 4(a)] or equivalently, by zeroes of the imaginary part of $f_{2}(r=0)$ [see Fig. 4(c)] as a function of the energy. The minima [Fig. 4(a)] and the zeroes [Fig. 4(c)] correspond to the actual eigenvalues and the spurious solutions are missing [cf. Tables 1 and 2]. For $m=1$ both $f_{1}$ and $f_{2}$ components need to vanish at the origin, and indeed no shift of the minima of Fig. 4(b) or zeroes of Fig. 4(d) is observed on the energy scale. Besides the solutions of the Bessel form, the results of Fig. 4 indicate the presence of the zero-energy levels that are supported by the zigzag edge gru .

The success of the sweeping approach is due to the two-point derivative that does not allow for saw-like 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 are 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 term given by Eq. (10) introduces the second derivative with the quotient $f^{\prime\prime}(i)\simeq\frac{f(i+1)+f(i-1)-2f(i)}{dr^{2}}$ that links the odd and even points and removes the rapidly oscillating states from the low-energy spectrum. The sweeping method that uses the backward quotient $f^{\prime}(i)=\frac{f(i)-f(i-1)}{dr}$ and passes through each mesh point consecutively introduces the link between the even and odd mesh points. The two-point formula for the smooth and rapidly oscillating solutions produce very different values of the quotient.

Tabel 3 shows that the convergence of the results of the present method to the exact solution is linear as a function of $dr$ (or $1/N$). Although the convergence is slow, the numerical cost of the calculation increases only linearily with $N$, and the method needs to keep track of only four complex numbers for the wave function components at adjacent mesh points, so one can approach the exact solution arbitrarily close at a negligible numerical cost.

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 applications.

VI.1 Limit of infinite-mass confinement

In the precedent section the zigzag-boundary conditions were used that result from a specific termination of the crystal lattice. The confinement of the particle described by the Dirac Hamiltonian can be induced provided by e.g. a spatial dependence of the energy gap nori of equivalently by an external potential due to e.g. a substrate sachs .

For $U_{A}\neq U_{B}$ an energy gap is opened in the dispersion relation near the Dirac points and the carriers acquire finite masses. For $U_{A}$ and $U_{B}$ that diverge to $\infty$ and $-\infty$ when $r>R$, respectively, the mass in the outside of the confinement is infinite and the wave functions of finite energies vanish at $r>R$ berry ; huang ; nipre . The boundary condition for this type of confinement is derived from vanishing current at the edge berry of the confinement area. The probability density current that is derived from the Hamiltonian berry is ${\bf j}=2v_{F}\left[\Re(\Psi_{1}^{*}\Psi_{2}),\eta\Im(\Psi_{1}^{*}\Psi_{2})\right]$. For the infinite mass at the outside of the dot, the orthogonal component of the current must vanish $(\cos\phi,\sin\phi)\cdot{\bf j}=0$, which implies $\tan\phi=-\eta\frac{\Re(\Psi_{1}^{*}\Psi_{2})}{\Im(\Psi_{1}^{*}\Psi_{2})}$.

With Eq. (2) this condition is translated for the radial functions as $\tan\phi=-\eta\frac{\Re(f_{1}^{*}(R)f_{2}(R)\exp(i\eta\phi))}{\Im(f_{1}^{*}(R)f_{2}(R)\exp(i\eta\phi))}.$ To fulfill this condition we set at the end of the flake $f_{1}(R)=1$, and $f_{2}(R)=i$, for which one gets $\tan\phi=\eta\tan\eta\phi$, which is fulfilled for both the valleys $\eta=\pm 1$.

Let us look at the appearance of the wave functions confinement by the mass boundary. We consider a spatial dependence of the energy gap as introduced by finite $U_{A}=-U_{B}$ potentials. Figure 6 shows the real part of $f_{1}$ and the imaginary part of $f_{2}$ for $\eta=1$ and $m=0$. At the end of the flake the infinite mass boundary condition is applied. The components $f_{1}$ and $f_{2}$ are – as for the zigzag boundary – purely real and purely imaginary, respectively. Moreover, the infinite-mass boundary condition implies $\Re\{f_{1}(r=R)\}=\Im\{f_{2}(r=R)\}$. The results of Fig. 6(a) were obtained for a flake of radius $R=80$ nm in the absence of the external potentials. In Fig. 6(b,c) we extended the flake to $R=100$ nm but introduced the energy gap beyond $r>R^{\prime}=80$ nm with $U_{A}=-U_{B}$. In Fig. 6(a), $U_{A}=0.1$ eV and $U_{A}=1$ eV in Fig. 6(c). We can see that for larger $U_{A}$ the functions penetrate only weakly into the gapped region for $r>R^{\prime}=80$ nm, and the values $\Re\{f_{1}\}$ and $\Im\{f_{2}\}$ become equal at $r=R^{\prime}$, where the gap is introduced. In this way the infinite boundary condition is found in the limit of large $U_{A}$ at $r=R^{\prime}$.

VI.2 The quantum ring spectra

The method with a slight modification can be applied to a problem of a graphene quantum ring tho . From the flake of the radius of $R=80$ nm we remove the central disk of a radius of $R_{i}=40$ nm. The infinite mass boundary condition in the inner edge of the flake reads $f_{2}(R_{i})/f_{1}(R_{i})=-i$. We start from the external edge $R$ as in the precedent subsection. The sweep at the finite difference mesh stops ar $R_{i}$. From the inner boundary condition we construct a function $F(E)=f_{2}(R_{i};E)+if_{1}(R_{i};E)$ and we look for its zeroes as a function of the energy.

The $m=0$ wave functions for the $K$ valley ($\eta=1$) found in this way are displayed in Fig. 7 for $n=1$ (a,b) and $n=2$ (c,d). The energy spectra for $\eta=\pm 1$ are given in Fig. 8 for $|m|\leq 4$. In the lowest energy states of the conduction and the valence bands one observes the characteristic angular momentum transitions with the Aharonov-Bohm periodicity of $\Delta B\simeq 0.385$ T pers , that corresponds to a flux quantum threading a ring of an effective radius $59$ nm – close to the central radius of the ring $R_{c}=(R+R_{i})/2$.

VI.3 Comparison with analytical results in the external magnetic field

In the absence of the external potential $U_{A}=U_{B}=0$, or for a massless particle, analytical results were obtained for the Weyl equation in the external magnetic field in Ref. gru for both zigzag and infinite mass boundary conditions. We adopted the graphene parameters of this work gru $t=2.7$ eV, the lattice constant $0.142$ nm and the radius of the dot $R=70$ nm. The results for the 500 mesh points are displayed in Fig. 9 with a perfect agreement with Fig. 1 of Ref. gru .

VII Summary

We have presented the finite difference method, applicable to Dirac carriers in circular confinement, that sweeps the mesh from the external edge of the system, where the boundary conditions are defined, to the center. The energies of the confined states are pinned by the internal boundary conditions. The method is simple, free of spurious solutions and the fermion doubling problem. Applications to circular flakes, quantum rings, and the energy gap modulated in space were presented.

Acknowledgments

This work was supported by the National Science Centre (NCN) according to decision DEC-2016/23/B/ST3/00821.