Semiclassical analysis of tunneling through a smooth potential barrier and localized states in graphene monolayer with mass gap Текст научной статьи по специальности «Физика»

SEMICLASSICAL ANALYSIS OF TUNNELING THROUGH A SMOOTH POTENTIAL BARRIER AND LOCALIZED STATES IN GRAPHENE MONOLAYER WITH MASS GAP

1V.V. Zalipaev, 2C.M. Linton

1 Krylov State Research Institute, St. Petersburg, Russia 2 Loughborough University, Loughborough, UK

v.zalipaev@lboro.ac.uk, c.m.linton@lboro.ac.uk PACS 03.65.-w, 03.65.Sq, 04.20Ha, 04.30Db

We present a semiclassical analysis of Dirac electron tunnelling in a graphene monolayer with mass gap through a smooth potential barrier in the ballistic regime. This 1D scattering problem is formulated in terms of a transfer matrix and treated in the WKB approximation. For a skew electron incidence this WKB approximation deals, in general, with four turning points. Between the first and the second, and the third and the fourth, turning points two tunnelling domains are observed. Scattering through a smooth barrier in graphene resembles scattering through a double barrier for the 1D Schrodinger operator, i.e. a Fabry-Perot resonator. The main results of the paper are WKB formulas for the entries of the barrier transfer matrix which explain the mechanism of total transmission through the barrier in a graphene monolayer with mass gap for some resonance values of energy of a skew incident electron. Moreover, we show the existence of modes localized within the barrier and exponentially decaying away from it and its behaviour depending on mass gap. There are two sets of energy eigenlevels, complex with small imaginary part and real, determined by a Bohr-Sommerfeld quantization condition, above and below the cut-off energy. It is shown that total transmission through the barrier takes place when the energy of the incident electron coincides with the real part of one of the complex energy eigenlevels. These facts were confirmed by numerical simulations performed using the finite element method (COMSOL).

Keywords: Graphene monolayer with mass gap, high-energy eigenstates, semiclassical approximation, generalized Bohr-Sommerfeld quantization condition.

1. Introduction

In this paper we study 2D electron transport in a graphene monolayer with mass gap with a smooth potential barrier. The general description of graphene electron transport may be found in the following reviews and publications [1-9]. Theory and experiments on electron transport in graphene and the related phenomenon of Klein tunnelling were described in many reviews and papers, for example, in [3], [10], [11], [12], [13], [14], [15]. It is worth mentioning the papers [16], [17], where analysis of electron-holes conductance oscillation in transport through barriers in graphene nano-ribbons was presented similar to 2D electron gas transport in semiconductors ( [18], [19], [20]).

Graphene, a sheet of single carbon atom honeycomb lattice, is well-known for its unique electron ballistic transport properties. It is regarded as an ideal medium for many applications such as graphene-based electronic devices [1]. However, graphene is a zero-bandgap semiconductor and exhibits semimetallic behaviour. Without bandgap opening it cannot be applied directly to semiconductors devices such as field-effect transistors which cannot operate in the absence of a bandgap in the material [21]. Thus, creating tunable bandgap gives an opportunity of enormous applications of graphene in digital electronics. In

the quest of creating the tunable bandgap, several physical and chemical approaches have been proposed and implemented successfully. This privides the motivation to study electron ballistic transport in graphene in the presence of a tunable bandgap. Generating the tunable bandgap for graphene in studies of electron ballistic transport means an appearance of a mass term in the Dirac system describing electron-hole quasi particle transport dynamics.

In the paper we present a semiclassical analysis of Dirac electron-hole tunnelling in a graphene monolayer with mass gap through a smooth barrier representing electrostatic potential in the ballistic regime. This 1D scattering problem is formulated in terms of transfer matrix and treated in the WKB (adiabatic) approximation for the Dirac system. For a skew electron incidence this WKB approximation deals with four turning points xi, i =1, 2, 3,4. For this scattering problem we assume that incident electron energy belongs to the middle part of the segment [0,Uo], where U0 is the height of the barrier. Thus, between xi and x2, and x3 and x4 we observe two tunnelling strips of total internal reflection similar to the well-known case in electromagnetic optics where the solution exponentially decay and grow. Between the asymptotically small neighbourhoods of xi, i = 1, 2, 3, 4 (boundary layers), we have five domains with WKB type asymptotic solutions, three domains with oscillatory behaviour and two with exponentially decaying and growing asymptotics. A gluing procedure between these five solutions is based on the matched asymptotic techniques (see [22], [23]) applied to so-called effective Schrodinger equation that is equivalent to the initial Dirac system (see [24], [25]). This gluing procedure leads to WKB formulas for the entries of the barrier transfer matrix that give all the transmission and reflection coefficients in this scattering problems.

It is worth mentioning that the scattering through a smooth barrier in graphene resembles scattering through a double barrier for 1D Schrodinger operator, i.e. a 1D Fabry-Perot resonator. From the point of view of physics of graphene, for positive energies close to one-half of the potential height U0 > 0 we observe incident, reflected and transmitted electronic states outside the barrier, whereas under the barrier we have a hole state (n-p-n junction). For negative energies close to one-half of the potential height U0 < 0 we observe incident, reflected and transmitted hole states outside the barrier, whereas under the barrier we have electronic states (p-n-p junction). It was Silvestrov and Efetov [17] who first demonstrated that in graphene a single parabolic barrier is similar to the double barrier potential well in GaAs/AlGaAs. They have shown that for the Dirac chiral relativistic electrons or holes there exist quasibound states. Parts of the left and right slopes of the parabolic potential barrier are acting like tunnelling barriers.

It is important to note that there is a strong difference between rectangular and smooth potential barriers. The presence of stable bound states for a rectangular barrier was first demonstrated by Pereira et al [26]. However in this case the double barrier structure does not arise and therefore the bound states are associated with a trapping into a potential well. Here we show that similar bound states can exist in a more generic situation when there is a single 1D smooth potential barrier formed in a graphene monolayer by some external electrostatic potential.

The main results of the paper are WKB formulas for the entries of the barrier transfer matrix which explain the mechanism of total transmission through an arbitrary smooth barrier in a graphene monolayer with mass gap for some resonance values of the energy of a skew incident electron. Crucially we show the existence of modes localized within the barrier and exponentially decaying away from it. There are two sets of energy eigenlevels, complex (quasibound states with small imaginary part i.e. long lifetime) and real (bound states), determined by a Bohr-Sommerfeld quantization condition, above and below the

cut-off energy, respectively. This differs from previous results. It is shown that the total transmission through the barrier takes place when the energy of an incident electron, which is above the cut-off energy, coincides with the real part of a complex energy eigenlevel of the first set of modes localized within the barrier. These facts were confirmed by numerical simulations for the reflection and transmission coefficients performed using the finite element method (COMSOL).

It is worth noting that the Dirac electron scattering by an arbitrary smooth barrier in a graphene monolayer but with zero mass gap in semiclassical approximation was briefly described in [27] stating only the main results. Here we present a more detailed general semiclassical analysis of Dirac electron tunnelling in a graphene monolayer with mass gap through a smooth potential barrier in the ballistic regime.

The paper is organized as follows. First, in section 2, we give a general details of the WKB description of tunnelling through a smooth barrier in a graphene monolayer with mass gap. In the next section we discuss WKB asymptotic solutions for the Dirac system with mass gap. In section, a 4 WKB asymptotic solution uniform with respect to pY for tunnelling through a smooth steps, left and right, is constructed. Building the total barrier transfer matrix is described in section 5. Analysis of quasi-bound and real bound states localized within the barrier is given in section 6. The results of numerical analysis are presented in section 7.

2. Tunnelling through a smooth barrier in graphene monolayer with mass gap

The procedure for deriving a theory describing elementary electronic properties of single layer graphene (see, for example, [9] and [4], section 2) works for electrons whose energy is close to Fermi level when their momenta are close to the Dirac points K and K' of the Brillouin zone. It uses the representation in the tight-binding 2D lattice Hamilton-ian and expanding the operators up to a linear order with respect to momentum within a neighbourhood of the Dirac points, and thus, leading to the effective Hamiltonian with Dirac operator in the first approximation.

Consider a scattering problem for the Dirac operator describing an electron-hole in the presence of a scalar potential representing a smooth localized barrier with the height U0 > 0 (see Fig.1) The problem can be described by the following 2D Dirac system (see, for example, [4])

where x = (x,y). Here u,v are the components of the spinor wave function describing electron density distribution localized on sites of sublattice A or B of the honeycomb graphene structure (see [9]), vF is the Fermi velocity, m = const is the electron effective mass, h is the Planck constant, and a = (ax,ay) with Pauli matrices

If we assume that the potential representing the smooth barrier does not depend on y, i.e. U = U(x), then we can look for a solution in the form

[vf (a • -

hV) + mvFa3 + U(x)]^(x) = E^(x), ^(x) = ( ,

(1)

(a)

Uo+Py U0-p~

(b)

Fig. 1. The three energy zones, shown in green - 1, pink - 2 and yellow - 3 colours are associated with the diverse character of scattering and localization for potential behaviour. The schematic shape of the two barrier potential for 1D Schroodinger, which is equivalent to the smooth potential barrier in graphene presented in the Fig.1(a) when the energy of quasi-particles belongs to 1 and 2 zones, i.e. —pY < E <U0 — pY.

where py means value of the transverse momentum component describing the angle of inci dence. Then, we obtain the Dirac system of two first order ODEs

U(x) — E + mv2F vF [—ihdx — ipy] \ ( u \ ( 0

vF [—ihdx + ipy ] U (x) — E — mvF

0

(2)

We assume that the potential U(x) has just one maximum and vanishes exponentially as x — ±c». Thus, U(x) is being localized within a strip directed along Y axis. It is more convenient to use the dimensionless system

U(x) — E + y —ihdx — ipy \ ( u \ /0

—ihdx + ipy U (x) — E — y

0

(3)

with the dimensionless rescaled variables that we shall be using below, that is, x/D — x with D being a characteristic length scale for the external potential, E/E0 — E with E0 being a characteristic energy scale. Then, we have vFpy/E0 — py, and U(x/D)/E0 ^ U(x). Then the mass term is given by

Y

mv

F

E0

const > 0.

(4)

The WKB solution may be constructed if the dimensionless parameter

hvF

h

E0D

is small (h ^ 1). Typical values of E0 and D are in the ranges 10-100 meV and 100-500 nm. For example, for E0 = 100 meV, D = 66 nm, we have h = 0.1.

Let pY = \Jp2y + Y2 ^ 0. In Fig. 1(a), 3 zones are shown that illustrate different scattering regimes for the smooth barrier scattering problem. These zones are exactly the same as for the rectangular barrier with the height U0. Below zone 1, E < —pY, we have total transmission and exponentially small reflection, asymptotic solutions are of oscillatory type everywhere. In zone 1 (green), —pY < E < pY (E = ±pY are the cut-off energy values),

v

v

2

there is no propagation outside the barrier, however there are oscillatory solutions within the barrier. In the zone 2 (pink), pY < E < U0 — pY, there are oscillatory solutions outside and within the barrier (zone of resonance tunnelling). In zone 3 (yellow), U0 — pY < E <U0 + pY, there is no propagation through the barrier, we have total reflection and exponentially small transmission. Above zone 3 E > U0 + pY, we have total transmission and exponentially small reflection, asymptotic solutions are of oscillatory type everywhere.

In the paper, we study the scattering problem for zone 2 and localized states in zone 1. In the first case, there are 5 domains with different WKB asymptotic solutions:

Q1 = {x : —<x> < x < x1}, n2 = {x : x1 < x < x2}, = {x : x2 < x < x3}, = {x : x3 < x < x4}, = {x : x4 < x <

where the turning points xi with i = 1, 2, 3, 4, are the roots of the equation

(E — U (x))2 — pY = 0. The regions Q1, Q3 and Q5, in which

(E — U(x))2 — pY( > 0, will be referred to as classically allowed domains, whereas and Q4, in which

(E — U(x))2 — pY < 0

are classically forbidden domains. Note that as pY ^ 0 for fixed value of E, the turning points coalesce.

These three energy zones are associated with diverse character of scattering through this potential barrier. The Dirac electron and hole bound states arise in a process of resonance tunnelling through this smooth potential. The quasi-bound states are to be found in the pink zone 2, confined by two tunnelling strips between the turning points x1, x2 and x3, x4, whereas the real bound states are located in zone 1 between x2 and x3. The schematic shape of the two barrier potential for 1D Schrodinger shown in Fig. 1(b) is equivalent to the smooth potential in graphene presented in the Fig.1(a) when the energy of quasi particles belongs to zones 1 and 2, i.e. —pY < E < U0 — pY. Quasi-bound states shown in Fig. 1(a) by dashed lines are confined by two tunnelling strips between x1, x2 and x3, x4.

It is worth remarking that when pY is fixed, and if E moves down from zone 2 to zone 1, the turning points x1 and x4 disappear (x1 ^ — cx>, x4 ^ such that inside zone 1 we have only x2 and x3. When we move down within zone 1, the turning points x2 and x3 get more separated. When E moves up from zone 2 to zone 3, the turning points x2 and x3 coalesce and disappear such that inside zone 3 we have only x1 and x4. When we move up from zone 3, the turning points x1 and x4 coalesce and disappear, and we obtain total transmission. In the paper we assume that U0 > 2pY. Thus, zone 2 should not disappear from the diagram in Fig.1(a).

3. WKB asymptotic solutions for Dirac system

It is convenient to introduce

W = ^ ,V = u——v. (5)

Then, the system (3) reads

(U — E)w + (y + ipy) V — ihw' = 0,

(U — E )V + (y — ipy )W + ihV' = 0. (6)

Eliminating V,

V = ihW' + (E — U )W (7)

Y + p ,

we obtain the so-called effective Schrodinger equation with complex coefficient

h2W" + (f2 — p\ + ihU')W = 0, (8)

where f = U(x) — E. The WKB oscillatory asymptotic solution to (8) in the classically allowed domains is to be sought in the form (see [23]) with real S(x)

W = e * S(x)J2 (—ih)j W3 (x). (9)

j=0

To the leading order we obtain

S = ±S,(x,x,,)= ±j pxdx, px = ^y - pY > n =1,2,3,4, (10)

xn

and up to a constant multiplier w0, for a wave traveling to the right we have

" Sp(x,xn)

W = wo 7~7P~~(1 + O(h)), (11)

yPx V ç + px

while for a wave traveling to the left we have

e-h Sp(X,Xn) Ç I p

W = woe J ç+ppx (1 + O(h)). (12)

VpX y TPj

Here —pY corresponds to solutions referred to x\y4 as U(x1)4) — E = —pY, whereas +pY corresponds to X2,3 as U (x2,3) — E = pY.

We seek the WKB exponentially decaying or growing asymptotic solution to (8) in the classically forbidden domains in the form (see [23]) with real S(x)

W = e h S(x)J^ hj Wj (x). (13)

j=o

To leading order we obtain

x

S = ±Sq(x, xn) = ± J qxdx, Qx = ^p2Y — Ç2 > 0, (14)

xn

and up to a constant multiplier w0, we have

W = wo e-e-2(arcsin 57 ± 2 )(1 + O(h)), (15)

x

or

W = w0-

e- 1 Sq (x,xn)

e 2 (arcsin £ ± 2 )(1 + O(h)),

(16)

and again +n/2 corresponds to solutions referred to x1)4, whereas —n/2 corresponds to x2)3. All these asymptotic solutions break down at the turning points £ = -pY for xi and x4, and £ = pY for x2 and x3.

4. WKB asymptotic solution for tunnelling through a smooth step (uniform asymptotics w.r.t. pY)

Consider the scattering problem for the smooth barrier under the assumption that pY < E < U0 — pY, U0 > 0 and all four turning points are present. From the point of view of scattering through the barrier we observe incident, reflected and transmitted electronic states at x < —a and x > a, whereas under the barrier — a < x < a we have a hole state (n-p-n junction, see Fig. 1). The transfer matrix for the left slope is defined as follows

d

Ta, TL = ( J,

rpL rpL T11 T 12 LL '21 T 22

d

a =

(17)

where we assume that U(-œ) = 0, and U(+œ) = U0 = const. This means that the potential U(x) behaves like a smooth monotone step. The connection coefficients a and d are defined in the asymptotic expansion of the solution at ^œ. Namely, as x — -œ, we have

4

1

\j2p- cos 91

a1e

hpx x+i$o

e-i&-/2 p-e

-, »-,*) + a2e-h - x-i*0( _"-,e-/2

aie-/2 -p~e~

a1ehPx xe1 + a2e

1 h Px xp1

(18)

_ ,- _ J(PX )2 + P

px = Ve2 -pY, p~ = E + Y

and $0 is a constant phase factor. As x — +œ, then

4

1

V2p+ cos 9+

d1ehp+x+i$0

-i -+

e i 2 + i -—

p+el 2

+ d2e

- hpix-i$0

i -+

e 2

+ -i -— , p+e 2

d1 e

h p+xe+ + a2e

-hp+xe+

(19)

and

px+

,- J(pt)2+pi

V<E - U0)2 - *, p+ = Uo^rrr

and 9± = arg (p± + ipy). This asymptotic behaviour at ^^ is very important as the corresponding transfer matrix satisfies canonical properties presented in Appendix A. This normalization leads to the conservation of the probability current (see (74)) written as

M2 — M2 = M2I2 — KI2. (20)

Then, the WKB asymptotic solution valid for x < x1 that matches (18) has to be written as follows

4

1

a1 e

h Sp(x,xi)-i$i

~pY

1

px + £ \ 1 +

Px+g

iPy +Y I + Px+g I +

iPy +Y,

ia2e-hsp/Px±Ç (1 + 2-2] ), (21)

'h sp ,

-pY ^ + 2+j J

where

; S

ar^ e (1 iPv + Y V P-(P- — E ) ^

$2 = arg {ie-i^(1 + P-+EJP-—E}. (22)

iPy + Y y — PyP-

The WKB asymptotic solution valid for x > x2 which corresponds to (19) has to be written as follows

*=M—idie h (1+S) +

dye-*M^)-^/px-+I (1 + g+A ) . (23)

V pY V1 + iPy+1) ' Thus, taking into account (5), for the WKB solution to the effective Schrodinger equation (8) we obtain

W = -L (aie*SrM-^ + iaye-SpM-^ ) (1 + 0(h)), (24)

Vpx Vf + px \I—py

for x < x1 (f < — pY), then,

W = -L ( — idie*+ dye-Sr(x>x2)-i*K ) (1 + O(h)), (25)

vpx V f + px V py

for x > x2 (f> pY). All these asymptotic expansions break down at the turning points x1y2.

According to [22], [23], we construct an asymptotic solution to (8), uniform with respect to pY by means of the comparison equation

1 ¿2

h2yzz + (h(v +2) — j)y = 0. (26)

By gluing this solution with (24) and (25), we derive the slope transfer matrix connecting a1y2 and d1>2. Thus, we seek a uniform asymptotic solution to (8) in the form

W = -2hv/2 ^ ~' (61Dv(h-1/2z) + 62D-v-1(ih-1/2z^ (1 + O(h)), (27) where solutions to (26)

y1,2(z) = Dv(h-1/2z), D-v-1(ih-1/2z)

are parabolic cylinder functions. According to the comparison equation method, after substitution of (27) into (8), we find that the function z(x) is to be determined from the equation

z2

z'2(a2 — zr) = q(x), (28)

where

a2 = h(v +i),q(x)= f2 - p2 + ihU' = q0(x)+ iW, (29)

and a or v are to be determined later. From (28) for x > x2 (£ > pY) we have

if \J q(t)dt = J ^ I —— a2dz,

z2

x2 2a

whereas for x < x1 (£ < -pY ),

xi -2a

il \iq(t)dt = I ^— a2dz,

where x1>2 are complex roots of q(x) = 0. The branches of the complex functions -\Jq(x) and yfz are fixed by the following asymptotic expansions for x >> x2 (|z| — œ)

z2 a2 a2

-hjVWyit - 4h - h '°gz - Th- log^)• <30>

x2

and x << x1

x:i

i I r-r-r, z2 a2, „ s a2

^W)dt - — - h log(-z) - — (1 - log a2). (31)

From (28) we obtain

2a

,/ J-qxdx =1,Ia2 - *dz = na2. (32)

xi -2a

On the other side, using the estimate

x2 x2

h i v-itjdt = h f V-qo(t)dt - + O(h),

xi xi

we have

x2

1 f nh

na2 = nh(- + v) = i / \J-q0(x)dx + —. (33)

xi

Hence, v is given by

v = Q, (34)

where

x2

Q1 = J \/-qo(x)dx. (35)

xi

The following estimates are very important when we glue the solution (27) with asymptotics (24) and (25). For x >> x2 we have

x x

h J vW)dt - h J VÔ0it)dt - ilog2£ - 4 + 2(1 + v) log , (36)

x2 x2

x

z

x

and for x << x1

X1

i h

VW)dt

h

i———, 1 2f 1 1 . , , V +1

vqitjdt — 2 log ^ — 4 + 2(1 + V)log -—.

(37)

Then, after taking into account (30) and (31), it follows that for x >> x2 (|z| >> |a|)

x

2f

and for x << x1

z a i i 1—TT 7 1 1 2f

4h — h l0g z ~ hj^q0(t)dt — 2 log pY + C'

z2 a2

4h — h l0g(—z)

x2

xi

i h

2f

— ' Vq0(t)dt — ^log —'p- + C

(38)

(39)

where

1

1

V +1

Z = 2h(1 — log a2) + 2(v + 2)log V

1 = 1(V + 2)(1 — log (hV)) — 1.

In a case of linear potential with constant U' we may use the substitute

z = V%U'ein/4(x — a), (40)

(see [25]). In this case equation (8) transforms exactly into (26). For a more general case of U(x), considered as a perturbation of the case of linear potential, our z belongs to a finite neighbourhood of the complex plane based on the line (40), where U' = U'(a) (E = U(a)). Thus, we assume that for large | z| the asymptotic expansions of the parabolic cylinder functions in (27) is applied in a way similar to the case of the linear potential.

Using the asymptotic expansions of the parabolic cylinder functions for large argument (see the appendix C), we obtain for x >> x2

W

M 2

(f2 — pI )1/4

[b1e-z2/4h zv h-v/2+

(41)

b2^ez2/4h-i 2 (v+1)z-v-1hv/2+1/2 + e-z2/4h-inv/2zv ^/2.

for x << x1 W

r(V +1)J

(bi [e-z2/4hZvh-v/2 - ez2/4h-nv -v-1hv/2+1/2 V^ ]+ (42) (f 2 — pi )1/4^D1[e zh e z h r( —V)]+ (42)

(—z)1/2

where r(z) is the Gamma function. Matching these two asymptotic expansions with the WKB asymptotics (24), (25), correspondingly, and using (38) and (39), leads to the following system

a1e

—i$i

= b1(—1)ve

-Z

ia2e d2e-i

(—b1e-inv rft + b2e-i n (v+1) )hv+1/2 (—1)-v-1eZ

r(-v)

b2eZ-i n (v+1)hv+1/2, (b1 + b2e-i2v)e-Z.

2

r^J

Let us introduce new notation:

Ho = -(-1)-v = -einv = -e-Q, (44)

= ieJny+2C ^ hv+1/2 U2 = -e~2C v/2n fi-v-1/2 (45)

ie r(-v)h , H2 e r(1 + v)h . (45)

Then, the system (43) reads

a^-^ = -b1x, 1 1 № '

ia2e-i*2 = b1 ^e-z - ib2ez-invhv+1/2Ho, -id1e-i®2 = -ib2ez-i n vhv+1/2,

d2e-i^i = bxe-z + b2ez+i n v hv+1/2 g.

(46)

Now eliminating b1 and b2 from the system (46), we obtain the relations determining the transfer matrix TL

d1 = -a1eia^ - a2 —, 1 1 № 2 № '

d2 = -a1(Ho - ^) + e-iaa2 fo,

that is

Ho Ho Ho Ho Ho

where a — $2 — $1. The expressions for H1 can be simplified as follows

H1 = —ivexp (inv + (v + 1)(1 — log (hv)) — 1) hv+1/2

- 11/ I-1 V » -

/ Q1 + .Q1 (1 ! (Q1 )) 1, \ V^ _iv ex^ + ^(1 - log( )) log ^.

2h 2 & yr(1 — v)'

Using the properties of the Gamma function (see [29])

|r(1 T v)| ' ™ ' 2Q1

sin (nv) V h(eQi/h - e-Qi/h) '

we derive where

(47)

rpL ia rpL 1 rpL ,, | ^1^2 1 rpL -ia № (

J11 = -e -, T 12 =--, T 21 = -Ho +--=--, T 22 = e -, (48)

H1 = ei0V 1 - e-2Qi/h, (49)

91 = 9(Q1) = Qh(1 - log (Qh)) - 4 - arg r(1 - iQh). (50)

Similarly, taking into account that

argr(1 + v) = — argr(1 — v),

we obtain _

H2 = —e-i9W 1 — e-2Qi/h. (51)

Hence, the left slope transfer matrix reads

TL(a,Q1)= — 1 \ t t

ei(e1+a)+Q1/h^1 _ e-2Qi/h eQ

e% e-i(8l+a)+Q1/h^1 _ e—2Q\/h

It is clear that the transfer matrix for the left slope satisfies all the properties in Appendix A, namely

Tl = (Tl)*, Tl = (T^y, det TL = -1.

One can easily understand that the quantities ri>2 mean the corresponding reflection coefficients, and t is the transmission coefficient.

It is worth remarking that due to the asymptotics as x ^

n 1

Imlog (r(—ix)) = — + x(1 — logx) + O(—),

if h (the turning points £ = ±pY do not coalesce), we observe that

tvi Qi\ i Qi\ i TV Qi\ n i i

argr(1 — ihn) = arg (—ihn) + argr(—ihn) = — 4 + hn(1 — log hn

and, consequently from (50), we obtain that 9i ^ 0. Thus, up to small exponential errors the transfer matrix (52) coincides with the corresponding non-uniform (wrt pY) asymptotic representation (121) obtained in the Appendix D.

For the right slope, uniform asymptotics for the entries of the transfer matrix may be obtained similarly to the derivation of the non-uniform asymptotic result (124) in the Appendix D, namely,

tr( Q ) i e^-a)+Q2/hV1 — e-2q2/n eq

1 (a,Q)=! eQ- e-i(e2-a)+q2/h^1 — e-2q2/h ). (53)

5. WKB asymptotic solution for scattering through a smooth barrier

Consider a problem of scattering through the smooth barrier (see Fig. 1(a)) under assumption that pY < E < U0 — pY and all four turning points xi} i =1, 2, 3, 4 are separated. In this case we have again 5 domains Qj, i = 1, 2,..., 5 to describe 5 WKB forms of solution to leading order. From the point of view of physics of graphene, we observe incident, reflected and transmitted electronic states at x < a and x > b, whereas under the barrier a < x < b we have a hole state (n-p-n junction, see Fig. 1(a)).

We formulate the problem for scattering through the barrier for the transfer matrix T, connecting the coefficients a\,2 of the asymptotics

^ = e h Sp(x,xi) ai6i + e-h sP(x,Xl)a2 e2, (54)

for x < x1, and d1y2 of the asymptotics

^ = e h Sp(xx)diei + e-h Sp(xx)d2e2, (55)

for x > x4, such that d = Ta, where e1>2 are corresponding eigenvectors (see (21)). Taking into account the transfer matrices obtained for scattering through the left and right slopes of the barrier (TL and TR in (52), (53)), the total transfer matrix T may be obtained as the product

T = TR(ehP 0i TL, (56)

\ 0 ehP J

where

x3 _

P = j \J{U(x) — E)2 — p21 dx.

x2

Then, the entries of the matrix T read

Tn = eh+h [s1s2ei(0ip) + e-iP ], (57)

T22 = eh+h [sls2e-i(ei+e2+p) + eiP], (58)

T12 = e-ia+h+h [s2ei(02+p) + s1e-i№+p) ], (59)

T21 = eia+h+h [s2e-i(02+p) + s^P)], (60)

where

si = \A — e-2Qi/h, i = 1, 2. They satisfy the classical properties of a transfer matrix (see Appendix B)

T22 = T*1, T21 = T*2, det T = 1. If a1 = 1, a2 = r1, d1 = t1, d2 = 0, then

, =-L = _ Ti

'1 rp , r1 rri

T 22 T 22

If a1 = 0, a2 = t2, d1 = r2, d2 = 1, then

t1 = 7^, H = - J21, |t112 + |n|2 = 1.

t2 = t1 = t, T2 = ^, |t2 |2 + |r2|2 = 1. T 22

Correspondingly, the unitary scattering matrix, defined as

a2 \ = S f a1 d1 d2

may be written as follows

M it)-

The transmission coefficient is given by

1

t = 7~ -T 22

Total transmission takes place only for symmetric barrier when Q2 = Q1 = Q (92 = 91 = 9). Then, we have

/ P 2Q P \-1

t = ei0( cos (- + 9)(2eh - 1) + i sin (- + 9) j , (61)

2 cos (p + 9)eh+i(0+a)V1 - e-2Q/h

11 = + ^Q ...... „ . (62)

cos (P + 9)(2eh - 1) + i sin (P +

It is clear that if

P(E) + h9 = hn(n + 1), n = 0,1,2,... , (63)

than we have total transmission ^ = 1.

It is worth noting that the formula for the transmission coefficient (61) was first

obtained in [25] for a single layer of graphene without gap.

6. WKB asymptotic solution for complex resonant and real bound states localized within the smooth barrier

Consider a problem of resonant states localized within the smooth barrier (see Fig. 1). For the sake of simplicity consider the case of symmetric barrier, that is, Q1 = Q2 = Q, 91 = 62 = 9. Now, let the turning points be —x2, —x1,x1,x2, and 0 < x1 < x2. Again, when the energy of electron-hole is greater than the cut-off energy (E > Ec = pY), we have 5 domains Qi, i =1, 2,..., 5 and 5 WKB forms of solution to leading order. To determine the correct radiation conditions as x ±c» that are necessary for the localization, we present WKB solutions in the domains 1 and 5 in the form

4 = e-h Sp(x'-x2)a2 e2, 4 = e h Sp(xx2)d1e1,

(64)

respectively, where e1)2 are the corresponding eigenvectors (see (18) in section 4). If a1 = 0, d2 = 0, a2 = 0, then

T22 (E ) = 0, (65)

and as a result we obtain a Bohr-Sommerfeld quantization condition for complex energy eigenlevels

i —2Q 1

h-1P(E) + 9 + 2 log (1 - e — ) = *(n + 2),

n

0,1, 2,

N1

(66)

for pY < E < Uo. Solutions to this equation are complex resonances En = Re(En) — iYn, where r-1 is the lifetime of the localized resonance. What is important is that the real part of these complex positive resonances decrease with n, thus showing off the anti-particle hole-like character of the localized modes. For these resonances we have rn > 0. From (66), we obtain the important estimate

r

hw 2Ât

w

- log (1 - e-2Q/h), Ât = -

dP

dE

(67)

E—En

This is the equivalent of the formula (14) in [17]. Namely, w is the transmission probability through the tunnelling strip, and At is the time interval between the turning points —x1 and x1. If pY ^ 0, then Q ^ 0, and rn ^ that is opposite to [17] (to be exact, the estimate for rn in [17] works only for linear potential when pY is not small).

For the second set for the bound states, when the energy of electron-hole is smaller than the cut-off energy (E < pY), we have 2 turning points —x1 and x1. Between them we have got oscillatory WKB solutions

4 = e h Sp(x'-xi)d1e1 + e-h Sp(x>-xi)d2e2,

or

^ = e h SP(X,X1) a1e1 + e-hh ^^ e2,

and outside decaying

^ = ehSq(x,-xi)c2l2, x < —x1, 4 = e-hhSq(x'xi)c1h, x > x1, where /1>2 are corresponding eigenvectors

(68)

(69)

h:

±2 (arcsin - 2) /1 T

y/ox

Y+iPy

1 iqx±l

Y+iPy ,

2

(a)

(b)

N 0.5 0.4 0.3

Fig. 2. The dependences of |i(E)| are shown for py = 0.^ - a and py = 0.2 -b for the potential U =1/ cosh2x, and the other parameters are h = 0 .1, Y = 0 . 1.

By gluing these WKB solutions through the two boundary layers near —xi and xi and using the techniques described in the Appendix D, we eliminate ai>2 and di>2 and come to the homogeneous system

ici + C2e h P = 0,

• - - --p ^ ici — c2e h = 0.

Thus, we derive the Bohr-Sommerfeld quantization condition for real energy eigenlevels (bound states) inside the cut-off energy strip for 0 < E < pY.

P(E) = hn(n +2), n = Ni + 1,...N2.

(71)

7. Numerical results

In this section we illustrate the effectiveness of the semiclassical formulae for the transmission coefficients (61) and the energy spectrum obtained from Bohr-Sommerfeld quantization conditions (66) by comparing their data with the results obtained by means of the finite element method (COMSOL). All the computations were done for the symmetric barrier

U

1

cosh 2x

(72)

The other dimensionless parameters are h = 0.1, y = 0.1. In Fig. 2a and 2b the dependences of |i(E)| are shown for py = 0.1 and py = 0.2. One can see a very good agreement between semiclassical results and the data obtained by means of finite element method. In Fig. 2a there are 3 resonances of total transmission (|i(E)| = 1) that take place for E = Re(E12,3) where E123 are obtained from (66) and E1 = 0.607 — i0.035, E2 = 0.371 — i0.016, E3 = 0.205 — i0.0044. In Fig. 2b there are only 2 resonances that take place for E = Re(E12) and E1 = 0.572 — ¿0.01, E2 = 0.343 — i0.002. Note that rn = Im(En) decrease with n what is in complete agreement with the thickness of the resonances shown in Fig. 2a and 2b.

In Fig. 3 the real parts of En for the first seven eigenvalues of quasibound (shown in the triangular domain only) and the next six En for the bound states are computed

(a)

(b)

0 -0.2 -0.4 -0.6

Fig. 3. The real parts of En for the first seven eigenvalues of quasibound (shown in the triangular domain only) and the next six En for the bound states are computed semiclassically with respect to py using Bohr-Sommerfeld quantization conditions (66) and (71) with mass gap y = 0.1 - a and y = 0.0 - b. The lines depicted with stars represent E = ±pY and E = U0 — pY.

semiclassically for the symmetric barrier

U

1

h (73)

coshx

with respect to py using Bohr-Sommerfeld quantization conditions (66) and (71) with mass gap y = 0.1 and y = 0.0. The structure of the energy eigenlevels has been deformed due to the presence of the mass gap. It is clearly seen that in the case of normal incidence py = 0 there are six bound states within the gap opposite to the gapless monolayer.

8. Appendix A. Transfer and scattering matrix properties for a smooth step

Let us come back to the scattering problem in terms of transfer matrix T (see [28]) for the left slope of the entire barrier formulated in (17). Taking into account the conservation of the ж-component of the probability density current (see equation (8) in [24] or (18) in [25])

Jx = Vu + Uv, (74)

we obtain that

\ai\2 — \a2\2 = \d2\2 -\d±\2. (75)

Thus, for the slope transfer matrix T it holds that

T21T22 — T11T12 = 0, \T21 \2 — \Tn\2 = 1, \T22\2 — \Tu\2 = —1

or

T+f 01 1 IT

As a result we have \T11\

we have

1 0 0 —1

\T22\, \T12\ = \T21\, \det(T)\ = 1. For the scattering matrix 5*11 S12

S =

S

S21 S22

(76)

(77)

(78)

(79)

E

E

barrier in graphene and

M2 + KI2 = H2 + |d2|2.

From (80) we obtain that

= = I,

thus, the scattering matrix is unitary. If entries of 5 are known, than,

T

— S11/S12 1/S12 S21 - S11S22/S12 S22/S

det(T) = -S21/S12.

741 (80) (81)

(82)

Time-reversal symmetry in scattering through the graphene barrier would mean that

(a34)* = e-hPx xa1e- + e hPx a*2e-, x G Q1,

(a34)* = e-hp+ d1e+ + ehpXd*2e+, x G Q3,

where

&3

1 0 01

are both asymptotic solutions to the Dirac system. Thus, we have

S

a2 a1

d1

(83)

(84)

(85)

In what follows that

T

(86)

Thus, S12 = S21,

and

s = st , i° i! t h = t*.

T

1 0/ v 1 0

detT = 1,

T11 T12

rp* rp* 1 12 T11

If a1 = 1, a2 = r1, d1 = 0, d2 = t1, then

1 TL t1 = Tl, i1 = - JL, N2 + |t1|2 = 1.

112 112

If a1 = 0, a2 = t2, d1 = 1, d2 = r2, then

1 TL

t2 = W, 12 = H, N2 + |t2|2 =1. 1 12 112

(87)

(88)

(89)

(90)

d

d

a

9. Appendix B. Transfer and scattering matrix properties for a smooth barrier

Let us formulate this scattering problem in terms of transfer matrix T for the entire barrier (see (54), (55), Ta = d). However, for the barrier we have

M2 — | a212 = |d1|2 — Id2l2, (92)

and

T +(1 M T = CiM. (93)

0 —1 0 —1

For the scattering matrix 5 we have

5 da21 = da12 , (94)

and

M2 + |d2|2 = |a2|2 + |d1|2. (95)

From (95) we obtain that

5+5 = 55+ = I. (96)

If entries of 5 are known, than,

rp (521 — 511522/512 522/512\ j c in (n.r,\

T = 1 — 511/512 1/512) , ) = 521/512. (97)

Taking into account the time-reversal symmetry in scattering through the graphene barrier, we obtain 5 = ST, and

T ( rp * rp *) , detT 1.

T12 T11

(98)

10. Appendix C

For the equation (8) (effective Schrodinger equation)

h2W" + (£2 - pY + ihU')W = 0, (99)

we shall have

W = — (a1ehSP(x,xi)-i*i JOËL + ia2e-SAx,xi)-i*2 ] (1 + O(h)), (100)

V Ç + Px y Py

for x < x1 (£ < -pY), then,

W = — (de-hSq(x'xi)e2(arcsin PV+2) + ^ehSq(x'xi)e-2(arcsin PV+2) )(1 + O(h)), (101)

(102)

for x1 < x < x2 (—pY < C < pY), and

_L ( — id1ehSP(xx2)-i*2 ljPl_ + d2e-Sp(x>x2)-i*i /Lt^x )(1 + 0(h)),

y 'px V C + px V pY

for x > x2 (C> pY).

To leading order, the uniform asymptotic representation within the neighbourhood of x1 is given by (see ( [23], chapter 4, section 3.3))

W = h-1^^ cosh^^6(11)Ai(h-2/3z) + b^Bi^-2^ + (103)

+ih1/6]J—Z sinh ^b^Ai'(h-2/3z) + b2)Bi'(h-2/3z^ , for x ^ x1 (£ ^ —pY), and

W = h-1/6jn cos b((1) Ai(h-2/3z) + b21]Bi(h-2/3z )) +

—ih1/e\[~n' sin ^ (b(11)Ai'(h-2/3z) + b^ Bi'(h-2/3z)^J , for x ^ x1 (£ ^ —pY). Here we have introduced

(104)

z = — ^2 JPxdx'^ < 0, x ^ xi (£ ^ —Py), (105)

z = Q JPxdx^ > 0, x ^ xi (£ ^ —p7),

2 ' ix

x1

(106)

^i=log w£ _+px, vi = 2(arcsin —+2). (107)

y pY 2 pY 2

Note that for this neighbourhood z'(x) > 0. According to ( [23], chapter 4, section 3.3), the functions A0 = cosh B0 = — sinh fii are the solutions to the following systems for x ^ xi (£ ^ — py ),

2z'^—zB'0 + Ao = 0,

2z'^—zA0 + Bo = 0, (108)

where

z'yf—z = px.

For x ^ xi (£ ^ — pY), A0 = cos vi, B0 = — i sin ^ are the solutions to

2z'^B'0 + iAo = 0,

2z'^zA'0 + iB0 = 0, (109)

where

z y/z = qx.

To leading order, the uniform asymptotic representation within the neighbourhood of x2 is similarly given by

W = cosh b(2)Ai(h-2/3z) + b(22)Bi(h-2/3z+ (110)

—ihi/6\—z\ Sinhb{2)Ai'(h-2/3z) + b{2)Bi'(h-2/3z)) for x ^ x2 (£ ^ pY ), and

W = h-i/6\[v\ C°s bf^Ai(h-2/3z) + b(2)Bi(h-2/3z(111)

ihi/^ y^zyz^T Sin ^b(2)Ai'(h-2/3z) + b(2)Bi'(h-2/3z)

for x ^ x2 (C ^ pY ). Here we have introduced

z = - (2 J Pxdx') / < 0, x ^ x2 (C ^ Py), (112)

x2

z = J Pxdx^ > 0, x ^ x2 (C ^ Py ), (113)

x

^2 = log WC +Px , v2 = 2(arcsin — - 72). (114)

PY 2 PY 2

Note that for this neighbourhood z'(x) < 0. According to ( [23], chapter 4, section 3.3), the functions A0 = cosh¡i2, B0 = sinh¡i1 are the solutions to the systems (108) for x ^ x2 (C ^ PY ), where

-z'y-z = Px.

For x ^ x2 (C ^ PY), A0 = cos v2, B0 = i sin v2 are the solutions to (109) with

z 'yfz = qx.

ia2e-i*2 = i(b()e-iп/4 + b^e^4). (116)

Now, using the asymptotic expansions of the Airy functions Ai(x) and Bi(x) (see Appendix C) and matching the uniform asymptotic expansions (103) and (104) within the neighbourhood of x1 with WKB asymptotics (100) for x ^ x1, we obtain

a1 e-^ = 1(b(11) ein/4 + b(21)e-in/4 ), (115)

1 2

Matching the uniform asymptotic expansions (103) and (104) within the neighbourhood of x1 with WKB asymptotics (101) for x ^ x1, we obtain

ein/4C1 = 1 b11)ein/4, e-in/4C2 = 1 b(2)e-in/4. (117)

Matching the uniform asymptotic expansions (110) and (111) within the neighbourhood of x2 with WKB asymptotics (101) for x ^ x2, we obtain

e-Qi/h+in/2Ci = b22), 2eQi/h-in/2c2 = bf0, (118)

where

x2

Q1 = J qx(x)dx.

xi

Matching the uniform asymptotic expansions (110) and (111) within the neighbourhood of x2 with WKB asymptotics (102) for x ^ x2, we obtain

-id1e-i*3 = 1(b12)e-in/4 + b22)ein/4), (119)

d2e-i*4 = 1(b(2)ein/4 + b22) e-in/4). (120)

Now we derive the entries of the transfer matrix TL connecting d and a

TL(a,Q,) = I Vf> ^ I, (12!)

e h + e-ia(eQ — ^)

where a = $2 — It is clear that the transfer matrix for the left slope satisfies all the properties (75)-(91).

For the right slope we have the following scattering problem

W = ( — iaiehJIK1 + a2e- SP(x,x3)-i*i ) (1 + O(h)), (122)

vpx V£+Px \jPi

for x < x3 (£ > pY), and

W = (diehSr(xx4)-^ HEEL + zd2e-SP(xx)-i*2 )(i + o(h)), (123)

Vpx V £ + Px V — py

for x > x4 (£ < —pY). By performing complex conjugation for (122), (123), and using (24), (25) and (121), derived for the left slope, we obtain

'—aie^ _ rrLf ( die1®1

^ a*2ei*1

and thus, we have d = TRa with

let*2} = TL(a Q ere*1}

ei*i ) = 1 ^^—dlel*2) '

TR(a,Q2)=| e Q h — JOT) e \ + ~ I, (124)

where

e—h" ia, e—h ■ e h + ^i- e (e h —

Q2 = J Qx(x)dx.

x3

11. Conclusion

Semiclassical analysis of Dirac electron with mass gap tunnelling through a smooth Gaussian shape barrier representing an electrostatic potential in the ballistic regime has been presented. The corresponding 1D scattering problem is formulated in terms of a transfer matrix and treated in the WKB (adiabatic) approximation. For skew electron incidence the WKB approximation deals with the asymptotic analysis of matched asymptotic techniques and boundary layers for four turning points. Scattering through a smooth barrier in graphene resembles scattering through a double barrier for 1D Schrodinger operator, i.e. a 1D Fabry-Perot resonator. The main results of the paper are WKB formulas for the entries of the transfer matrix. They explain the mechanism of total transmission of a Dirac electron in graphene with mass gap through a smooth barrier for some resonance values of the energy of skew incident electrons. Moreover, we have showed the existence of modes localized within the barrier and exponentially decaying away from it for two discrete complex and real sets of energy eigenlevels determined by a Bohr-Sommerfeld quantization condition. It has been shown that total transmission through the barrier takes place when the energy of the incident electron coincides with the real part of one of the complex energy eigenlevels. These facts have been confirmed by numerical simulations performed using the finite element method (COMSOL). It is necessary to note that some of the details of the comparison equation method of WKB analysis used in the paper require more rigorous treatment. These part of the work along with further detailed discussion on other scattering regimes through the barrier will be carried on in future publications.

Acknowledgments

The authors would like to thank Prof S.Yu. Dobrokhotov, Prof J. Ferapontov, Prof F. Kusmartsev, Dr D. Maksimov and Dr A.Vagov for constructive discussions and valuable remarks.

References

[1] A. K. Geim and K. S. Novoselov. The rise of graphene. Nature Mater. 6, 183, 2007.

[2] M. I. Katsnelson. Graphene: carbon in two dimensions. Mater. Today, 10:20, 2007.

[3] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim. Chiral tunnelling and the Klein paradox in graphene. Nature Physics, 2, 620, 2006.

[4] Castro Neto, A.H., Guinea, F., Peres, N.M.R., Geim, A.K., and Novoselov, K.S. The electronic properties of graphene. Rev. Mod. Phys. 81, Jan-March, 109-162, 2009.

[5] C. W. J. Beenakker. Colloquium: Andreev reflection and Klein tunnelling in graphene. Rev. Mod. Phys, 80:1337, 2008.

[6] A. K. Geim. Graphene: Status and prospects. Science, 324:1530, 2009.

[7] N. M. R. Peres. Colloquium: The transport properties of graphene: An introduction. Rev. Mod. Phys, 82:2673, 2010.

[8] M. A. H. Vozmediano, M. I. Katsnelson, and F. Guinea. Gauge fields in graphene. Phys. Rep., 496:109, 2010.

[9] Semenoff G.W. Condensed-matter simulation of the three-dimensional anomaly. PRL 53, 26 (1984).

10] O. Klein. Die Reflexion von Elektronen an einem Potentialsprung nach der relativistischen Dynamik von Dirac. Z. Phys, 53:157-165, 1929.

11] A. Calogeracos and N. Dombey. History and physics of the Klein paradox. Contemp. Phys., 40:313, 1999.

12] N. Dombey and A. Calogeracos. Seventy years of the Klein paradox. Phys. Rep., 315:41, 1999.

13] S. V. Vonsovsky and M. S. Svirsky. The Klein paradox and the Zitterbewegung of an electron in a field with a constant scalar potential. Uspekhi Fiz. Nauk, 163(5):115, 1993.

14] R.-K. Su, G. G. Siu, and X. Chou. Barrier penetration and Klein paradox. J. Phys. A, 26:1001, 1993.

15] Cheianov, V.V., Fal'ko, V.I. Selective transmission of Dirac electrons and ballistic magnetoresistance of n-p junctions in graphene. Phys.Rev. B. 74, 041403(R), 2006.

16] Fistul, M.V., Efetov, K.B. Electromagnetic-Field-Induced Suppression of Transport through n-p Junctions in Graphene. PRL 98, 256803, 2007.

17] Silvestrov, P.G., Efetov, K.B. Quantum dots in graphene. PRL 98, 016802, 2007.

18] Datta, S., 1995 Electronic transport in mesoscopic systems, Cambridge University Press, Cambridge.

19] Mello, P. A., Kumar, N., 2004 Quantum Transport in Mesoscopic systems, Oxford University Press, New York.

20] Stockmann, H. J., 2000, Quantum Chaos. An Introduction, Cambridge University Press, Cambridge.

21] K. S. Novoselov et al.. Electric Field Effect in Atomically Thin Carbon Films. Science. 306, 666, 2004.

22] F. W. J. Olver. Asymptotics and Special Functions, Academic Press, 1974.

23] M. V. Fedoryuk. Asymptotic Analysis: Linear Ordinary Differentia,l Equations. Springer, 1993.

24] Sonin, E.B. Effect of Klein tunneling on conductance and shot noise in ballistic graphene. Phys.Rev. B. 79, 195438, 2009.

25] Tudorovskiy, T., Reijnders, K.J.A., Katsnelson, M.I. Chiral tunneling in single-layer and bilayer graphene. Phys. Scr. T146, 014010, 17pp, 2012.

26] Pereira, J.M., Mlinar, Jr.V., Peeters, F.M., Vasilopoulos, P. Confined states and direction-dependent transmission in graphene quantum wells. Phys.Rev. B. 74, 045424, 2006.

27] V. Zalipaev, D.Maximov, C.M.Linton and F.V.Kusmartsev. Quantum states localized in a smooth potential barrier on single layer graphene. Phys. Let. A, 377, 216-221, 2012.

28] P. Markos and C Soukoulis. Wave Propagation. From electrons to photonic crystals and left-handed materials. Priceton University Press, 2008.

29] M. Abramowitz and I. Stegun, editors. Handbook on Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York, 1965. Chapter 19 on parabolic cylinder functions.

30] The Wolfram special functions site: http://functions.wolfram.com/Hypergeometrical Functions/Parabolic CylinderD/.