Localised States of Fabry-Perot Type in Graphene Nano-Ribbons

This book collects some new progresses on research of graphene from theoretical and experimental aspects in a variety of topics, such as graphene nanoribbons, graphene quantum dots, and graphene-based resistive switching memory. The authors of each chapter give a unique insight about the specific intense research area of graphene. This book is suitable for graduate students and researchers with background in physics, chemistry, and materials as reference.


Introduction
Graphene has been spoken of as a "'wonder material"' and described as paradigm shifting in the field of condensed matter physics [1]. The exceptional behavior of single layer graphene is down to its charge carriers being massless, relativistic particles. The anomalous behavior of graphene and its low energy excitation spectrum, implies the emergence of novel electronic characteristics. For example, in graphene-superconductor-graphene junctions specular Andreev reflections occur [1] and in graphene p-n junctions a Veselago lens for electrons has been outlined [2]. It is clear that by incorporating graphene into new and old designs that new physics and applications almost always emerges. Here we investigate Fabry-Perot like localized states in graphene mono and bi-layer graphene. As one will no doubt appreciate, there are many overlaps in the analysis of graphene with the studies of electron transport and light propagation. When we examine the ballistic regime we see that the scattering of electrons by potential barriers is also described in terms of transmission, reflection and refraction profiles; in analogy to any wave phenomenon. Except that there is no counterpart in normal materials to the exceptional quality at which these occur, with electrons capable of tunneling through a potential barrier of height larger than its energy with a probability of one -Klein tunneling. So, normally incident electrons in graphene are perfectly transmitted in analogy to the Klein paradox of relativistic quantum mechanics. A tunable graphene barrier is described in [3] where a local back-gate and a top-gate controlled the carrier density in the bulk of the graphene sheet. The graphene flake was covered in poly-methyl-methacrylate (PMMA) and the top-gate induced the potential barrier. In this work they describe junction configurations associated with the carrier types (p, for holes and n for electrons) and found sharp steps in resistance as the boundaries between n-n-n and n-p-n or p-n-p configurations were crossed. Ballistic transport was examined in the limits of sharp and smooth potential steps. The PMMA is a transparent thermoplastic that has also been used to great effect in proving that graphene retains its 2D properties when embedded in a polymer heterostructure [4]. The polymers can be made to be sensitive to a specific stimulus that leads to a change in the conductance of the underlying graphene [4] and it is entirely likely that graphene based devices of the future will be hybrids including polymers that can control the carrier charge density. In [5] an experiment was performed to create a n-p-n junction to examine the ballistic regime. Oscillations in the conductance showed up as interferences between the two p-n interfaces and a Fabry-Perot resonator in graphene was created. When there was no magnetic field applied, two consecutive reflections on the p-n interfaces occurred with opposite angles, whereas for a small magnetic field the electronic trajectories bent. Above about 0.3 Tesla the trajectories bent sufficiently to lead to the occurrence of two consecutive reflections with the same incident angle and a π-shift in the phase of the electron. Thus, a half period shift in the interference fringes was witnessed and evidence of perfect tunneling at normal incidence accrued.
Quantum interference effects are one of the most pronounced displays of the power of wave quantum mechanics. As an example, the wave nature of light is usually clearly demonstrated with the Fabry-Perot interferometers. Similar interferometers may be used in quantum mechanics to demonstrate the wave nature of electrons and other quantum mechanical particles. For electrons they were first demonstrated in graphene hetero-junctions formed by the application of a top gate voltage [6]. These were simple devices consisting mainly of the resonant cavity, and with transport channels attached. These devices exhibited quantum interference in the regular resistance oscillations that arose when the gate voltage changed.
Within the conventional Fabry-Perot model [7,8], the resistance peaks correspond to minima in the overall transmission coefficient. The peak separation can be approximated by the condition 2k F L = 2πn. The charge accumulates a phase shift of 2π after completing a single lap (the round-trip) 2L in the resonant cavity, where k F is the Fermi wave vector of the charges, and L is the length of the Fabry-Perot cavity. This is the Fabry-Perot-like resonance condition: the fundamental resonance occurs when half the wavelength of the electron mode fits inside the p-n-p junction representing the Fabry-Perot cavity.
The simplest electron cavity, but still very effective, for the Fabry-Perot resonator may be formed by two parallel metallic wire-like contacts deposited on graphene [9]. There in a simple two terminal graphene structure there are clearly resolved Fabry-Perot oscillations. These have been observed in sub-100 nm devices. With a decrease of the size of the graphene region in these devices, the characteristics of the electron transport changes. Then the channel-dominated diffusive regime is transferred into the contact-dominated ballistic regime. This normally indicates that when the size of the cavity is about 100 nm or less the Fabry-Perot interference may be clearly resolved. The similar Fabry-Perot interferometer for Dirac electrons has been recently developed from carbon nanotubes [10].
Earlier work on the resistance oscillations as a function of the applied gate voltage led to their observation in the p-n-p junctions [6,11]. It was first reported by Young and Kim [6], but the more pronounced observations of the Fabry-Perot oscillations have been made in the Ref. [11]. There high-quality n-p-n junctions with suspended top gates have been fabricated. They indeed display clear Fabry-Perot resistance oscillations within a small cavity formed by the p-n interfaces.
The oscillations arise due to an interference of an electron ballistic transport in the p-n-p junction, i.e. from Fabry-Perot interference of the electron and hole wave functions comprised between the two p-n interfaces. Thus, the holes or electrons in the top-gated region are multiply reflected between the two interfaces, interfering to give rise to standing waves, similar to those observed in carbon nanotubes [12] or standard graphene devices [13]. Modulations in the charge density distribution change the Fermi wavelength of the charge carriers, which in turn is altering the interference patterns and giving rise to the resistance oscillations.
In the present work we consider a simplest model of the Fabry-Perot interferometer, which is in fact the p-n-p or n-p-n junction formed by a one dimensional potential. We develop an exact quasi-classical theory of such a system and study the associated Fabry-Perot interference in the electron or hole transport.
Although graphene is commonly referred to as the "'carbon flatland"' there has been a feeling of discontent amongst some that the Mermin-Wagner theorem appeared to be contradicted. However, recent work shows that the buckling of the lattice can give rise to a stable 3D structure that is consistent with this theorem [14]. In what follows we present the general methodology for analysis of graphene nanoribbons using semiclassical techniques that maintain the assumption of a flat lattice. However, it should be mentioned that the effects found from these techniques are powerful in aiding our understanding of potential barriers and are an essential tool for the developing area of graphene barrier engineering. The natural state of graphene to accommodate defects or charged impurities is important for applications. The p-n interfaces described above may be capable of guiding plasmons and to create the electrical analogues of optical devices to produce controllable indices of refraction [15].
In Part I of this chapter we investigate the use of powerful semiclassical methods to analyze the relativistic electron and hole tunneling in graphene through a smooth potential barrier. We make comparison to the rectangular barrier. In both cases the barrier is generated as a result of an electrostatic potential in the ballistic regime. The transfer matrix method is employed in complement to the adiabatic WKB approximation for the Dirac system. Crucial to this method of approximation for the smooth barrier problem, when there is a skew electron incidence, is careful consideration of four turning points. These are denoted by x i , i = 1, 2, 3, 4 and lie in the domain of the barrier. The incident electron energy in this scattering problem belongs to the middle part of the segment [0, U 0 ], where U 0 is the height of the barrier, and essentially the incident parameter p y should be large enough to allow normal and quasi-normal incidence.
Therefore, between the first two turning points, x 1 and x 2 , and also between the next two, x 3 and x 4 there is no coalescence. Two columns of total internal reflection occur which have solutions that grow and decay exponentially. Looking away from the close vicinity of the asymptotically small boundary layers of x i , there exists five domains with WKB type solutions (See Fig. 2): three with oscillatory behavior and two exhibiting asymptotics that are exponentially growing and decaying. Combining these five solutions is done through applying matched asymptotics techniques (see [16]) to the so-called effective Schrödinger equation that is equivalent to the Dirac system (see [17], [18]). This combinatorial procedure generates the WKB formulas that give the elements of the transfer matrix. This transfer matrix defines all the transmission and reflection coefficients in the scattering problems discussed here.
When the energies are positive around potential height 0.5U 0 , electronic incident, reflected and transmitted states occur outside the barrier. Underneath the barrier a hole state exists (n-p-n junction). The symmetrical nature of the barrier means that we see incident, reflected and transmitted hole states outside the barrier when the energies are negative and close to one-half of the potential height U 0 < 0. Thus, underneath the barrier there are electronic states (a p-n-p junction).
Incorporated into the semiclassical method is the assumption that all four turning points are spatially separated. Consequently, the transverse component of the momentum p y is finite and there is a finite width to the total internal reflection zone. This results in a 1-D Fabry-Perot resonator, which is of great physical importance and may aid understanding in creating plasmonic devices that operate in the range of terahertz to infrared frequencies [19]. Quantum confinement effects are crucial at the nano-scale and plasmon waves can potentially be squeezed into much smaller volumes than noble metals. The basic description of propagating plasma modes is essentially the same in the 2-D electron gas as in graphene, with the notable exception of the linear electronic dispersion and zero band-gap in graphene [20]. Thus, we predicate that the methods applied here are also applicable to systems of 2-D electron gases, such as semiconductor superlattices. Due to the broad absorption range of graphene, nanoribbons as described here, or graphene islands of various geometries may also be incorporated in opto-electronic structures.
In our analysis, if p y → 0 then we have a quasi-normal incidence whereby and first two, x 1 and x 2 , and the second two, x 3 and x 4 , turning points coalesce. In the case of normal incidence, there is always total transmission through the barrier. The vital discovery in this form of analysis is that of the existence of modes that are localized in the bulk of the barrier. These modes decay exponentially as the proximity to the barrier decreases. These modes are two discrete, complex and real sets of energy eigen-levels that are determined by the Bohr-Sommerfeld quantization condition, above and below the cut-off energy, respectively. 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 the complex energy eigen-level of one among the first set of modes localized within the barrier. These facts have been confirmed by numerical simulations for the reflection and transmission coefficients using finite elements methods (Comsol package).
In Part II we examine the high energy localized eigenstates in graphene monolayers and double layers. One of the most fundamental prerequisites for understanding electronic transport in quantum waveguide resonators is to be able to explain the nature of the conductance oscillations (see [25], [26], [27]). The inelastic scattering length of charge carriers is much larger than the size of modern electronic devices and consequently electronic motion is ballistic and resistance occurs due to scattering off geometric obstacles or features (e.g. the shape of a resonator micro or nano-cavity or the potential formed by a defect). It is an interesting area of development whereby defects are engineered deliberately into devices to generate a sought effect. In graphene, defects such as missing carbon atoms or the addition of adatoms can lead to interesting and novel effects, e.g magnetism or proximity effects. In the ballistic regime, conductance is analyzed by the total transmission coefficient and the Landauer formula for the zero temperature conductance of a structure (see monographs [25], [26], [27], and papers [28], [29], [30], for example). The excitation of localized eigenmodes inside a quantum electronic waveguide has a massive effect on the conductivity because these modes could create an internal resonator inside the waveguide. This is a very good reason to research the role of localised eigenmodes for quantum resonator systems and 2-D electronic transport in quantum waveguides. Excitation of some modes could result in the the emergence of stop bands for electronic wave propagation in the dispersion characteristics of the system, whereby propagation through the waveguide is blocked entirely. Other modes will result in total transmission.
In this review, the semiclassical analysis of resonator eigenstates that are localized near periodic orbits is developed for a resonator of Fabry-Perot type. These are examined inside graphene monolayer nanoribbons in static magnetic fields and electrostatic potentials. The first results for bilayer graphene are also presented in parallel to this.
Graphene has generated a fervor throughout the scientific world and especially in the condensed matter physics community, with its unusual electronic properties in tunneling, charge carrier confinement and the appearance of the integer quantum Hall effect (see [31], [33], [34]), [35], [17])). Its low energy excitations are massless chiral Dirac fermion quasi-particles. The Dirac spectrum, that is valid only at low energies when the chemical potential crosses exactly at the Dirac point (see [31]), describes the physics of quantum electrodynamics for massless fermions, except that in graphene the Dirac electrons move with a Fermi velocity of v F = 10 6 m/s. This is 300 times smaller than the speed of light. Graphene is a material that is easy to work with, it has a high degree of flexibility and agreeable characteristics for lithography. The unusual electronic properties of graphene and its gapless spectrum provide us with the ideal system for investigation of many new and peculiar charge carrier dynamical effects. It is also conceivable, if its promise is fulfilled, that a new form of carbon economy could emerge based upon exploitation of graphenes novel characteristics. The enhancements in devices are not just being found at the nano and micron-sized levels, though these hold the most potential (e.g. the graphene transistor, metamaterials etc), but in composites [36], electrical storage [37], solar harvesting [38] and many more applications. Following this train of thought, graphene is also a viable alternative to the materials normally used in plasmonics and nanophotonics. It absorbs light over the whole electromagnetic spectrum, including UV, visible and far-infrared wavelengths and as we have mentioned, it is capable of confining light and charge carriers into incredibly small volumes. Thus, there are a range of applications where band gap engineering is not required and it is satisfactory to directly use nanoribbons of graphene as optical-electronic devices.
In the analysis of graphene one also expects unusual Dirac charge carrier properties in the eigenstates of a Fabry-Perot resonator in a magnetic field. For example, two parts of the semiclassical Maslov spectral series with positive and negative energies, for electrons and holes, correspondingly, with two different Hamiltonian dynamics and families of classical trajectories are apparent. Semiclassical analysis can provide insight into the aforementioned physical systems and good quantitative predictions on quantum observables using classical insights. Application of semiclassical analysis in studying the quantum mechanical behavior of electrons has been demonstrated in descriptions of different nano-structures, electronic transport mechanisms in mesoscopic systems and, as another example, the quantum chaotic dynamics of electronic resonators [25], [26], [27], [39], [40], [41], [42] and many others.
However, it is important to state that the first semiclassical study on two-dimensional graphene systems only recently appeared in [43], [44], [45]. In [43] a semiclassical approximation for the Green's function in graphene monolayer and bilayers was discussed. In [44] and [45] bound states in inhomogeneous magnetic fields in graphene and graphene-based Andreev billiards were studied by semiclassical analysis, accordingly. This was carried out with one-dimensional WKB quantization due to total separation of variables.
In the second half of this review, the semiclassical Maslov spectral series of the proliferation of high-energy eigenstates (see [48], [49] [50]) of the electrons and holes for a resonator formed inside graphene mono and bilayer nanoribbons with zigzag boundary conditions, is specified. These states are localized around a stable periodic orbit (PO) under the influence of a homogeneous magnetic field and electrostatic potential. The boundaries of the nanoribbon act with perfect reflection to confine the periodic orbit to isolation. This system is a quantum electron-hole Fabry-Perot resonator of a type analogous to the "bouncing ball" high-frequency optical resonators found in studies of electromagnetics and acoustics. The asymptotic analysis of the high-energy localized eigenstates presented here is similar to ones used for optical resonators (see [50], [51], [54], and [55]). In this review, the semiclassical methods presented focus upon the stability of POs and electron and hole eigenstates that depend on the applied magnetic field.
We construct a solitary localized asymptotic solution to the Dirac system in the neighborhood of a classical trajectory called an electronic Gaussian beam (Gaussian wave package). In PO theory there are similarities between the asymptotic techniques used here and those used in the semiclassical analysis (see, for example, [27] (chapters 7, 8) or [39] and cited references). Further, the stability of a continuous family of closed trajectories in asymptotic proximity to a PO, confined between two reflecting interfaces, is studied. The classical theory of linear Hamiltonian systems with periodic coefficients gives the basis to study the stability using monodromy matrix analysis. The asymptotic eigenfunctions for electrons and holes are constructed only for the stable PO as a superposition of two Gaussian beams propagating in opposite directions between two reflecting points of the periodic orbit. A generalized Bohr-Sommerfeld quantization condition gives the asymptotic energy spectral series (see [46] and [47], [48], [49], [50], [51] and [55]). This work highlights that the single quantization condition derived herein for the quantum electron-hole graphene resonator fully agrees with the asymptotic quantization formula of a quite general type spectral problem in [51]. It is worth drawing attention to the fact that in a semiclassical approximation for the Green's function in a graphene monolayer and bilayer, the relationship between the semiclassical phase and the adiabatic Berry phase was discussed in the paper [43]. Our asymptotic solutions, for rays and Gaussian beams, possess the adiabatic phase introduced by Berry [64]. The importance of Berry-like and non-Berry-like phases in the WKB asymptotic theory of coupled differential equations and their roles in semiclassical quantization were discussed in [57], [58], [59].
Our results are a special class of POs that occur for graphene zigzag nanoribbons in a homogeneous magnetic field and piece-wise electrostatic potential that is embedded inside the nanoribbon. They are found by giving, to the leading order, a description of the general form of asymptotic solution of Gaussian beams in a graphene monolayer or bilayer. The key point in the asymptotic analysis is the quantization of the continuous one-parameter (energy) family of POs. For one subclass of lens-shaped POs, these localized eigenstates were evaluated against eigenvalues and eigenfunctions that have been computed by the finite element method using COMSOL. For a selectively chosen range of energy eigenvalues and eigenfunctions, agreement between the numerical results and those computed semiclassically is very good. In the graphene Fabry-Perot resonator, the electrostatic potential does not play a role of confinement, it behaves more like an inhomogeneity, but in some cases an electrostatic potential helps to make a family of POs stable.
In this chapter, we describe the tunneling through smooth potential barriers and the asymptotic solutions for a Dirac system in a classically allowed domain. This is done using WKB methods. We then go on to investigate the classically disallowed domain and tunneling through the smooth barrier. The asymptotic WKB solutions are presented for scattering and for quasi-bound states localized within the smooth barrier. The second part of the chapter, goes into detail about high energy localized eigenstates in monolayer and bilayer graphene. The graphene resonator is described when it is subjected to a magnetic field and ray asymptotic solutions outlined. Finally, the construction of periodic orbits, stability analyses and quantization conditions are thoroughly examined. A numerical analysis is given that compares the analytical techniques and results with those found using finite element methods.

The rectangular barrier
In a conventional metal or semiconductor there are no propagating states connecting regions either side of the barrier (regions I and I I I). To get through the barrier an electron has to tunnel through the classically forbidden region and the tunneling amplitude depreciates exponentially as a function of the barrier width. Thus, transport between I and I I I is strongly suppressed. However, in each of the three regions of a barrier in a graphene system, the valence and conduction band touches, meaning that there are propagating states connecting I and I I I at all energies. There is no such suppression of the transport at energies incident and below the barrier. At normal incidence transmission is always perfect.
Potential barriers for single quasiparticle tunneling in graphene can be introduced by designing a suitable underlying gate voltage or even as a result of local uniaxial strain [68]. In the following we denote the angle of incidence with respect to the barrier to be θ 1 . We are interested in the dependence of the tunneling transmission on this incidence angle. To illustrate quantum mechanical tunneling one must extract the transmission coefficient from the solution to the graphene barrier problem. The transmission coefficient is the ratio of the flux of the particles that penetrate the potential barrier to the flux of particles incident on the barrier. We demonstrate a rectangular barrier as described in detail in the Reviews by Castro Neto et al [60] and Pereira Jr et al [69]. The problem can be described by the following 2D Dirac system (see, for example, [60]) where x = (x, y) and u, v are the components of the spinor wave function describing electron localization on the sites of sublattice A or B of the honeycomb graphene structure, v F is the Fermi velocity, the symbol (, ) means scalar product,h is the Planck constant andσ = (σ 1 , σ 2 ) with Pauli matrices If we assume that the potential representing the barrier does not depend on y, i.e. U = U(x), then we can look for a solution in the form where p y means value of the transverse momentum component describing the angle of incidence. Then, we obtain the Dirac system of two ODEs The particle incident with energy E < U 0 from the left of the barrier has wavevectors k 1 , q, and k 2 to the left, in the barrier and to the right of the barier, respectively. These regions are denoted I, I I and I I I, respectively. In the symmetric barrier k 1 = k 2 = k. Region I I lies between −L and L, where ±L defines the width of the barrier. The wave functions are defined for each of the three regions below: where we have introduced the wave function, as is done in [31]. The coefficients c 1 , c 2 and a 1 , a 2 are related by means of the transfer matrix, c = Ta. The transfer matrix has unique properties, which are demonstrated in Appendix B. In regions I and I I I the angle of incidence in momentum space is given by, θ 1 = arctan k y /k x and in region I I, θ 2 = arctan k y /q x . In regions I-I I I the valence and conduction bands touch. This allows propagating states to connect the regions at all energies and there is no suppression of transport at the energies below the height of the barrier. There is also perfect transmission at normal incidence. The graphene rectangular barrier can be thought of as a medium with a different refractive index to its surroundings. In an optical analogy, the refractive index of the barrier is 1 − U 0 /E [8]. At the interface of the barrier the incidence angle splits into transmitted and reflected waves with the transmitted wave propagating with angle θ 2 through the barrier. The wave inside the barrier is multiply reflected between −L and L. The parallel wave vector inside the barrier is given by, y and the wave vectors outside the barrier are defined as, The wave functions in regions I and I I are matched at x = −L. Likewise, the wavefunctions between regions I I and I I I are matched at x = L. It is not necessary to match the derivatives, as is done in an analysis using the Schrödinger equation. One requires the wave functions to be continuous at the boundary of each region to generate relationships between the coefficients, a 1,2 , b 1,2 and c 1,2 . We seek solutions such that The elements of the transfer matrix for the rectangular barrier are found to be, where we make the substitutions α = e iθ 1 + e −iθ 2 and β = e −iθ 2 − e −iθ 1 and their complex conjugate forms are denoted byᾱ = e iθ 2 + e −iθ 1 andβ = e iθ 2 − e −iθ 1 . If we assume that the incident wave approaches from the left, then a 1 = 1, a 2 = r 1 and c 1 = t 1 , where r 1 is the reflection coefficient and t 1 is the transmission coefficient. If the incident wave approaches from the right then c 1 = r 2 , c 2 = 1 and a 2 = t 2 . We find that t 1 = t 2 = t and the transmission The reflection coefficients are determined as r 1 = −T 21 /T 22 and r 2 = T 12 /T 22 . The transmission probability is as usual given by |t 1 | 2 with the definition |t 1 | 2 + |r 1 | 2 = 1. At normal incidence the carriers in graphene are transmitted completely through the barrier (Klein tunneling). However, the carriers can be reflected by a potential step when the angle of incidence increases and a non-zero momentum component parallel to the barrier ensues. Thus, the transmission of charge carriers through the potential barrier is anisotropic. When a beam of electrons is fired at an angle into the barrier, it splits up into transmitted and reflected beams, with multiple reflections occurring at the edges of the barrier. As is usual in quantum mechanics, the transmission is found by stipulating that there must be continuity between the wavefunctions. In the above this demand for continuity at the extremities of the barrier allowed us to find the coefficients of the wavefunctions. Thus, using these results and following the work of Castro Neto et al [60], the total transmission as a function of the incident angle is given by T(θ 1 ) = tt * : When the tunneling resonance condition 2Lq x = nπ is met, where n is an integer, T = 1. This statement means that a half-integer amount of wavelengths will fit into the length of the potential barrier. The absolute transmission is the manifestation of Klein tunneling, which is unique for relativistic electrons, and it should occur when an incoming electron starts penetrating through a potential barrier of height, U 0 (which is far in excess of the electrons rest energy). The transport mechanism in a graphene tunneling structure is unique. This perfect transmission at incidence normal to the barrier is due to the pseudo-spin conservation that gives no backscattering. In order to attain an interference effect between the two interfaces an oblique incidence angle is required and it is under this prerequisite that multiple interference effects emerge. Thus, the potential barrier is analogous to two interfaces at −L and L and also a Fabry-Perot interferometer [5]. The analogy of the graphene rectangular barrier to the Fabry-Perot resonator when θ 1 = 0 extends to the potential barrier operating like an optical cavity. In region I I the incoming wave can interfere with itself and with constructive interference, resonances will occur where T(θ 1 = 0) = 1 [5]. The potential barriers for single quasiparticle tunneling in graphene are usually created by suitably changing the underlying gate voltage. In the next section we investigate the smooth barrier and expect that there will be similar scattering behavior as through the rectangular barrier. We seek to explore the similarities and the differences between the two.

The smooth barrier
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 U 0 (see Fig.2). It is convenient to use the dimensionless system The generalization of a smooth potential barrier with Gaussian shape (we assume that p y > 0). The Dirac electron and hole states arising in resonance tunneling are shown. The quasibound states are to be found above the green strip, |E| < p y , where bound states are located. Quasibound (metastable) states are confined by two tunneling strips between x 1 , x 2 and x 3 , x 4 , whereas the bound states are located between x 2 and x 3 .
in which we omitted the sign of tilde for brevity. In physical dimensions the energy is U 0 E, the potential is U 0 U(x), the y-component of the momentum is p y U 0 /v F , and the dimensionless Planck constant (small WKB parameter) is given by h =hv F /U 0 D, where U 0 is the height of the potential barrier (|U(x)| < 1 for x ∈ R) and D is a characteristic scale of the potential barrier with respect to the x-coordinate. Typical values of U 0 and D are within the ranges 10-100meV and 100-500nm. For example, for U 0 =100meV, D = 264nm, we have h = 0.025 and also we assume that p y > 0. The six zones are separated by the four red diagonal lines, E = ±p y and E = ±p y + U 0 . We will now discuss the right hand side of this diagram. In zone 1 (orange shading), E < −p y , there is total transmission and exponentially small reflection. The asymptotic solutions are of an oscillatory nature everywhere in this zone. In zone 2 (blue), −p y < E < p y with the cut-off energy at E = ±p y . In zone 2 there is no propagation outside the barrier. However, there are oscillatory solutions within the barrier. Zone 3 (green) is the zone of Klein tunneling. Here, p y < E < U 0 − p y and there are oscillatory solutions everywhere. Zone 4 (hexagons), U 0 − p y < E < U 0 + p y , is devoid of propagation through the barrier. There is total reflection and exponentially small transmission in this zone. The fifth zone (sand color) is limited to E > U 0 + p y and is characterized by total transmission, exponentially small reflection and the asymptotic solutions are oscillatory everywhere. The sixth zone is one of no propagation and exponentially decaying or growing asymptotic solutions, U 0 < E < p y .
In Fig. 3, six zones (horizontal strips in Fig. 3b) are shown that illustrate different scattering regimes for the smooth barrier scattering problem. These six zones are exactly the same as for the rectangular barrier with the height U 0 . In zone 1 E < −p y , we have total transmission and exponentially small reflection, asymptotic solutions are of oscillatory type everywhere. In zone 2, −p y < E < p y (E = ±p y is the cut-off energy), there is no propagation outside the barrier, however there are oscillatory solutions within the barrier. In the zone 3, p y < E < U 0 − p y , there are oscillatory solutions everywhere (zone of the Klein tunneling). In zone 4, U 0 − p y < E < U 0 + p y , there is no propagation through the barrier, we have total reflection and exponentially small transmission. In zone 5 E > U 0 + p y , we have total transmission and exponentially small reflection, asymptotic solutions are of oscillatory type everywhere. Finally, in the zone 6, U 0 − p y < E < p y , there is no propagation, everywhere asymptotic solutions are of exponential type, decaying or growing.
Firstly, we study the scattering problem for zone 3 (see Fig.2). In this case, there are 5 domains with different WKB asymptotic solutions: The regions Ω 1 , Ω 3 and Ω 5 , in which (E − U(x)) 2 − p 2 y > 0, will be referred to as classically allowed domains, whereas Ω 2 and Ω 4 , in which (E − U(x)) 2 − p 2 y < 0, are classically disallowed domains. Note that as p y → 0 for fixed value of E, the turning points coalesce. We exclude this possibility in our analysis.
It is worth to remark that for fixed p y , when E moves down from zone 3 to zone 2, the turning points x 1 and x 4 disappear (x 1 → −∞, x 4 → +∞) such that inside zone 2 we have only x 2 and x 3 . When we move down from zone 2 to zone 1, the turning points x 2 and x 3 disappear. When E moves up from zone 3 to zone 4, the turning points x 2 and x 3 coalesce and disappear such that inside zone 4 we have only x 1 and x 4 . When we move up from zone 4 to zone 5, the turning points x 1 and x 4 coalesce and disappear.

WKB asymptotic solution for Dirac system in classically allowed domain
The WKB oscillatory asymptotic solution to the Dirac system in the classically allowed domains is to be sought in the form (see [16]) with real S(x) Substituting this series into the Dirac system, and equating to zero corresponding coefficients of successive degrees of the small parameter h, we obtain a recurrent system of equations which determines the unknown S(x) (classical action) and ψ j (x), namely, where I is the identity matrix and S ′ = p x . The Hamiltonian H has two eigenvalues From now on we will omit the dependence on x of U, S, and quantities derived from them. It turns out to be convenient to use different e 1,2 instead with In this way we will be able to solve problems of electron and hole incidence on the barrier simultaneously. Note that, irrespective of whether E > U or E < U, The classical action S(x) is given by the sign indicating the direction of the wave, with + corresponding to a wave traveling to the right.
For electrons and holes one can seek a solution to the Dirac system zero-order problem in the form with unknown amplitude σ (0) . The solvability of the problem (H − EI)ψ 1 = − Rψ 0 requires that the orthogonality condition < e 1 , R(σ (0) e 1 ) >= 0 must hold, written as a scalar product implied with complex conjugation, and from this one obtains the transport equation for σ (0) : It has a solution where a branch of the analytic function √ z is taken that satisfies the condition For higher-order terms, one can seek a solution to where from (15), σ (j) 2 is given by Then, from the orthogonality condition, where c j = const. Below we assume that p x > 0, corresponding to a wave traveling in the positive x-direction. Thus, to the leading order we have This asymptotic approximation is not valid near turning points where S ′ = 0 (see Fig. 1) at x = x j , j = 1, 2, 3, 4 where e iθ = ±i and cos θ = 0, while at x = a, b we have E = U. The WKB asymptotic solution, derived in this section, is valid for the domains Ω i , i = 1, 3, 5.

WKB asymptotic solution for Dirac system in classically disallowed domain
The WKB asymptotic solution to the Dirac system in the classically disallowed domain is to be sought in the form with S(x) real. As in section 4, we obtain a recurrent system of equations which determines the unknown S(x) and ψ j (x), namely, where S ′ = q x , and the matrix R is as in (14). The Hamiltonian H is not Hermitian. It has two eigenvalues and not orthogonal eigenvectors Hl 1,2 = h 1,2 l 1,2 , where Thus, the function S(x) in a classically disallowed domain is given by Again, for the sake of simplicity, we shall use different l 1,2 where For electrons and holes one can seek a solution to the Dirac system zero-order problem in the form with unknown amplitude σ (0) . Solvability of the problem requires that the orthogonality condition must hold The vector l 1 is the eigenvector of H, whereas l * 1 is the eigenvector of H * . From the orthogonality condition one obtains the transport equation for σ (0) It has a solution For higher-order terms, we have (H − h 1 I)ψ j = − Rψ j−1 and one should seek solution in the form where σ (j) 2 is given by Then, from the orthogonality condition, < l * 1 , R(σ (j) Below we assume that q x > 0. Thus, to the leading order in classically disallowed domains we have where This asymptotic approximation is not valid near turning points q x = 0. The WKB asymptotic solution, derived in this section, is valid for the domains Ω i , i = 2, 4.

WKB asymptotic solution for scattering through the smooth barrier
Consider a problem of scattering through the smooth barrier (see Fig. 2) under the assumption that |E| > |p y | and all four turning points x i , i = 1, 2, 3, 4 are separated. In this case we have again 5 domains Ω i , i = 1, 2, ..., 5 to describe 5 WKB forms of solution to the leading order. In considering a graphene system, if E > 0 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. 2).
To formulate the scattering problem for transfer matrix T, here we present the WKB solutions in the domains 1 and 5 The barrier is represented by the combination of the left and right slopes. The total transfer matrix T, that is d = Ta, is given by with T R and T L the transfer matrices of the right and left slopes (see formulas (137), (143) in Appendix C), respectively, and The entries of the matrix T read (see formulas (121), (134), (144) in Appendix C) where s i = √ 1 − e −2Q i /h , i = 1, 2. They satisfy the classical properties of the transfer matrix and if a 1 = 1, a 2 = r 1 , d 1 = t 1 , d 2 = 0, then If a 1 = 0, a 2 = t 2 , d 1 = r 2 , d 2 = 1, then The transmission coefficient t = 1/T 22 , looks exactly like the formula (131) in [18]. Total transmission takes place only for a symmetric barrier when However, it is worth noting that r 1 (p y ) = r 2 (−p y ). It is clear that if P(E) + hθ = hπ(n + 1 2 ), n = 0, 1, 2, ... , then we have total transmission |t 1 | = 1.

WKB asymptotic solution for complex resonant (quasibound) states localized within the smooth barrier
Consider a problem of resonant states localized within the smooth barrier (see Fig. 2). In the first case when the energy of the electron-hole is greater than the cut-off energy (E > E c = |p y |), we have five domains Ω i , i = 1, 2, ..., 5 and five WKB forms of solution to the leading order. To determine the correct radiation conditions that are necessary for the localization, we present WKB solutions in the domains 1 and 5 If a 1 = 0, d 2 = 0, then and as a result we obtain Bohr-Sommerfeld quantization condition for complex energy eigen-levels (quasi-discrete) for |p y | < E < U 0 . Solutions to this equation are complex resonances E n = Re(E n ) − iΓ n , where Γ −1 n is the lifetime of the localized resonance. What is important is that the real part of these complex positive resonances is decreasing with n, thus showing off the anti-particle hole-like character of the localized modes. For these resonances we have Γ n > 0. From (45), we obtain the important estimate that is the equivalent of the formula (14) in [35]. Namely, w is the transmission probability through the tunneling strip, ∆t is the time interval between the turning points x 2 and x 3 , and P ′ is the first derivative of P with respect to energy. If p y → 0, then Q → 0, and Γ n → +∞, that is opposite to [35] (to be exact, the estimate for Γ n in [35] works only for a linear potential when p y is not small).  For the second set of real resonances, when the energy of the electron-hole is smaller than the cut-off energy (E < |p y |), we have 2 turning points x 2 and x 3 . Between them there are oscillatory WKB solutions or and outside decaying solutions, By gluing these WKB solutions together through the two boundary layers near x 2 and x 3 , using the results in sections 5.1, 5.2 and the Appendix C, we eliminateā 1,2 andd 1,2 and obtain the homogeneous system ic 1 +c 2 e i h P = 0, Thus, we derive the Bohr-Sommerfeld quantization condition for real energy eigen-levels (bound states) inside the cut-off energy strip for 0 < E < |p y |.

Numerical results
Based upon the analytical descriptions in the preceding sections for the smooth barrier, we present the results for the energy eigenvalues and eigenfunctions. These are shown in Fig's. 4-6 and compare favorably with those obtained through finite difference methods, as detailed in [71]. The energy dispersion curves, E n (p y ), are shown for the complex resonant and real bound states for h = 0.1 and potentials of different widths. In Fig. 4(a), the energy levels are illustrated for the potential, U = 1/coshx, with n = 0, 1, ...., 15. For complex resonant states the real parts are shown. It must be emphasized that in zone "3", which is restricted by E < U 0 − p y and E = p y with p y > 0, the complex quasibound states reside. The bound states are located in zone "2", which lies between E = ±p y and E = U 0 − p y . In zone "3" there are nine complex resonances. In Fig. 4(b), the results for a narrower potential of U = 1/cosh2x can be seen (all other parameters being the same as in Fig. 4 (a)). In this case, there are four complex resonances in zone "3" and n = 0, 1, ...., 9. The lifetimes of the local resonances, Γ n , are shown in Fig. 5 (a) and (b) for the same two potentials as described in Fig. 4. The complex quasi-bound states that are witnessed in zone "3" in Fig. 4 are shown in Fig. 5 for Γ n . The thinner potential allows less complex bound states. The bound states have infinite lifetimes. Both types of states are confined within the barrier in the x-direction, while the motion in the y-direction is controlled by the dispersion relations. In Fig. 6 we present the transmission probabilities |t| 2 for the two potentials. There are nine tunneling resonances, i.e. complex quasi-bound states within potential barrier defined as U = 1/coshx that correlate with those shown in Fig. 4 in zone "3". Likewise, for the thinner barrier there are four complex quasi-bound states. These resonance states are a clear indication of the Fabry-Perot multiple interference effects inside the barrier.

Graphene resonator in a magnetic field
We consider a spectral problem for the Dirac operator describing the electron-hole quantum dynamics in a graphene monolayer without a gap, in the presence of a homogeneous magnetic field A and arbitrary scalar potential U(x) (see [31]) where x = (x 1 , x 2 ), and u, v are the components of the spinor wave function that describes electron localization on the sites of sublattice A or B of a honeycomb graphene structure.
Here, e is the electron charge, c is the speed of light, A is magnetic potential in axial A = B/2(−x 2 , x 1 , 0) or Landau gauge A = B(−x 2 , 0, 0) (magnetic field is directed along the x 3 axis), and v F is the Fermi velocity. The symbol <, > means a scalar product, andh is the Planck constant (which is a small parameter (h → 0) in semiclassical analysis). The vector σ = (σ 1 , σ 2 ) with Pauli matrices corresponds to the K Dirac point of the first Brillouin zone (see [31]). The case of the second K ′ Dirac point withσ * = (σ 1 , −σ 2 ) is treated similarly and is not considered here. Figure 7. A periodic orbit inside the graphene nanoribbon resonator with magnetic field and electrostatic potential (electronic trajectory). Magnetic field is directed along the x 3 axis, the electrostatic field is piece-wise linear U(x 2 ) = β|x 2 |.
We study the high energy spectral problem, using the semiclassical approximation, for a vertical graphene nanoribbon confined between two flat reflecting interfaces L 1,2 (see Fig.7). It is assumed that the spinor wave function satisfies zigzag boundary conditions on the interfaces L 1,2 : u| L 1 = 0, v| L 2 = 0. It will be discussed later that the electrostatic field U(x 2 ) = β|x 2 | makes the orbit shown in Fig. 7 periodic and stable. In the gener al case, as it was noted earlier for the Schrödinger operator (see [55]), if high-energy localised eigenstates are sought, which decay exponentially away from the resonator axis AB, the separation of variables will not help construct an exact solution due to the difficulty of satisfying the boundary conditions.

Ray asymptotic solution
The WKB ray asymptotic solution to the Dirac equation is sought through consideration of the eigenvalue problem associated with Hφ = Eφ. The magnetic vector potential A = B/2(−x 2 , x 1 , 0) is given in terms of the axial gauge for a magnetic field. The Hamiltonian of the Dirac system (see equations (2) and (10)) takes the form for monolayer graphene, In contrast, for bilayer graphene the Hamiltonian takes the form, where g ≈ 0.4eV/υ F is the interlayer coupling. We consider the case when bilayer graphene has Bernal stacking as shown in Fig. 6. Bernal-stacked bilayer graphene occurs with half of the carbon atoms in the second layer sitting on top of the empty centers of hexagons in the first layer. An external electric field can tune its bandgap by up to 250meV [32]. This form of structure of bilayer graphene can be experimentally created using chemical vapor deposition [32]. Considering the low energy states of electrons, we can reduce the 4 × 4 matrix describing the bilayer graphene to a form similar to that of monolayer graphene [56]. The only difference from the monolayer form is the squaring of the off-diagonal entries and the inclusion of a band mass for bilayer electrons.
It is now convenient to introduce some dimensionless variables. The coordinate system x ⇒ xD, where D is a characteristic scale associated with a change in the potential (and correspondingly U(x) ⇒ U(xD)). Then, we write U = U(xD)/E 0 ; where we define the characteristic energy scale as E ⇒ E 0 E. For single layers of graphene the small parameter, h << 1, is h = υ Fh /U 0 D. In double layer graphene it is slightly different; h = αD/ √ 2mU 0 , with the magnetic length, as a function of the applied magnetic field, given to be α = αD/ √ 2mE 0 with α = eB/c. We now write the dimensionless forms of the one and two layer graphene systems as (with the tildes omitted for brevity), and The solution for monolayer graphene is sought in the same form as equation (11). Substituting this series into the Dirac system, and equating to zero the corresponding coefficients of successive degrees of the small parameter h, we obtain a recurrent system of equations which determines the unknown S(x) and ψ j (x).
The Hamiltonian H has two eigenvalues. In the domain Ω e = {x : E > U(x)}, the Hamiltonian eigenvalue h 1 = U(x) + p on the level set h 1 = E describes the dynamics of electrons. The corresponding classical trajectories can be obtained from the Hamiltonian systemẋ = H e p ,ṗ = −H e x , x = (x 1 , x 2 ), p = (p 1 , p 2 ), with an equivalent Hamiltonian (see [48]) on the level set H h = 0. The Hamiltonian dynamics with h 1,2 or with H e,h are equivalent (see [48]).Classical action S(x) satisfies the Hamilton-Jacobi equation in the monolayer case to be Likewise, in the case of bilayers, and The Hamiliton-Jacobi equation is satisfied in the two-layers of graphene case by, The requires that the orthogonality condition with complex conjugation < e 1 , R(σ (0) (x)e 1 ) >= 0 must hold, where where, Using the basic elements of the techniques described in [48], from the orthogonality condition one obtains the transport equation for σ (0) (x). The geometrical spreading for an electron or hole with respect to the Hamiltonian system with h 1,2 = U ± v F p has a solution Monolayer : Bilayer : where where we have introduced θ, which is the angle made by the velocity of the particle trajectory with the x 1 axis: Here −θ/2 is the adiabatic phase for monolayer graphene, as introduced by Berry [64]. Chirality results in a Berry phase of θ in bilayer graphene and the confinement of electronic states. Conservation of chirality in monolayer graphene means that the particles cannot backscatter and this leads to normal incidence transmission equal to unity. This is not the case in bilayer graphene and backscattering can occur.

Construction of eigenfunctions, periodic orbit stability analysis and quantization conditions
Let x 0 = (x 1 (s), x 2 (s)) be a particle (electron or hole) classical trajectory, where s is the arc length measured along the trajectory. Consider the neighborhood of the trajectory in terms of local coordinates s, n, where n is the distance along the vector, normal to the trajectory, such that where e n (s) is the unit vector normal to the trajectory. Introducing ν = n/ √ h = O(1), we seek an asymptotic solution to the Dirac system related to (2) where S 0 (s) and S 1 (s) are chosen similar to [55], [66] as they give a linear approximation for solution to the Hamilton-Jacobi equation (55) (see [55], [66])). The parameter for monolayers a(s) = E − U 0 (s) and for Bernal bilayers, In the following γ i (s), i = 1, 2 are the Cartesian components of e n (s). Following [54], [55], [65], and [66], we apply the asymptotic boundary-layer method to the Dirac system (2). We allow that the width of the boundary layer is determined by |n| = O( √ h) ash → 0. We assume that we deal with a continuous family of POs symmetric with respect to both axes (see Fig. 4). Thus, the trajectory of the PO consists of two symmetric parts between two reflection points A and B. We seek the asymptotic solution of the eigenfunction for electrons and holes localized in the neighborhood of a PO as a combination of two Gaussian beams where we have defined the Berry phase to be, e iθ = γ 2 − iγ 1 , and e 1,2 = (1, e iθ )/ √ 2. Here, each beam is related to the corresponding part of the periodic orbit. Namely, ψ 1 is determined by z = z 1 (s), p = p 1 (s) for 0 < s < s 0 , describing the electrons propagation along the lower part of the orbit from A to B, whereas ψ 2 is determined by z = z 2 (s), p = p 2 (s) for s 0 < s < 2s 0 , for the electrons propagation along the upper part of the orbit from B to A. The complete derivation of the electronic Gaussian beams for monolayer graphene can be found in the work of Zalipaev [66]. Following the methodology developed in the previous work [66], we state the problem in terms of the function z(s) and write the Hamiltonian in its terms, with the Hamiltonian, The above are the same for both mono and bilayer graphene, but with different d(s) (and a(s), as mentioned above), where ρ(s) is the radius of curvature of a trajectory. Thus, (z 1 (s), p 1 (s)) and (z 2 (s), p 2 (s)) define (z(s), p(s)) for 0 < s < 2s 0 . The asymptotic localized solution of Gaussian beam ψ(s, n) is constructed in an asymptotically small neighbourhood of the PO. This solution is to be periodic with respect to s ∈ R with the period 2s 0 , and satisfies the zigzag boundary conditions. The reflection coefficient R is derived in the short-wave approximation, and given by where γ is the angle of incidence, and δ 1 = θ(s 0 + 0) − θ(s 0 − 0). In the bilayer graphene system the reflection coefficient is, The localized solution can be constructed if z(s), p(s) is a complex (in the complex phase space C 2 z,p ) quasi-periodic Floquet solution of Hamiltonian system (65) with periodic coefficients (see [50], [54]). Namely, for the monodromy 2 × 2 matrix M, describing the mapping for the period 2s 0 , a Floquet solution for arbitrary s is defined as The structure of the monodromy matrix M is given by the following product where M 1 = M 1 (s 0 ) and M 2 = M 2 (2s 0 ) are fundamental matrices of the system (65) describing the evolution (z(s), p(s)) for 0 < s < s 0 and s 0 < s < 2s 0 , correspondingly.
The reflection matrices at points A and B (see Fig. 6), R A and R B are given by , where γ is the angle of incidence of the trajectory at the points A, B. To attain R A and R B , the classical action S of the phase function at the reflecting boundary requires continuity to be set between the incident and the reflected beams (see [50], [54], [51]).
In a general case, the entries of M 1,2 are to be determined numerically as the Hamiltonian system has variable coefficients. It is worth to remark that all the multipliers in (71) are symplectic matrices. Thus, the monodromy matrix M is symplectic.
The classical theory of linear Hamiltonian systems with periodic coefficients states that, if |TrM| < 2, we have a stable PO (elliptic fixed point, for example, see [27]), and ||M n || < const for arbitrary n ∈ N. Then, there exist two bounded complex Floquet's solutions for −∞ < s < +∞, namely, (z(s), p(s)) and (z(s),p(s)) with Floquet's multipliers λ 1,2 = e ±iϕ (0 < ϕ < π), which are solutions of These solutions (z(s), p(s)) and (z(s),p(s)) may be obtained as follows. If |TrM| < 2, the monodromy matrix has complex eigenvectors w = (w z , w p ) andw = (w z ,w p ) Then, the first solution is determined by Satisfaction of the solution to the zigzag boundary condition at the interface L 1 (u| L 1 = 0), as well as the requirement that the solution should be periodic with respect to s ∈ R with the period 2s 0 lead to a generalized Bohr-Sommerfeld quantization condition determining semiclassical asymptotics of the high energy spectrum. Namely, after the integration around the closed loop of PO, the total variation of the classical action S and the phase of the amplitude of ψ 0 must be equal to 2πm 1 . Thus, we obtain the quantization condition for electrons and holes in the form Monolayer : ∆ = π − γ 2 , where m 1,2 ∈ N are the longitudinal and the transversal quantization indexes, and for electrons we have +, for holes −. The index m 2 and factor ∆ appear due to the variation of the phase of σ m 2 (s, ν) (see the formulas in [66]). Here in the left-hand side in (72)

Numerical results
In this section we concentrate on the example for monolayer graphene with piece-wise linear potential U(x 2 ) = β|x 2 |. The numerical techniques used in this section are described in [55]. We deal with the Dirac system (2) by incorporating the Landau gauge A = B(−x 2 , 0, 0). Thus, using dimensionless U, E, α and dimensionless coordinates, the Dirac system is written in the following form The energy in eV is given by U 0 E, where U 0 = v Fh /hD = 6.59meV/h is the characteristic scale of the potential U. Here we assume that D = 10 −7 m. A new small dimensionless parameter h (0 < h << 1) is supposed to be predetermined. The magnetic induction amplitude is given by B = αcU 0 /v F eD = α/h 6.5910 −2 T. Consider, as an example, a family of continuous POs which are symmetric with respect to both axes, with two reflection points A, B (see Fig. 1). The formulas describing electronic POs as solutions of the corresponding integrable system with the Hamiltonian in the Landau gauge on the level set H = 0, are easily obtained and given by for the lower part 0 < t < t 0 , and for the upper part 0 < t < t 0 . Here π 1 and π 2 are the initial values of the components of momentum p 1 and p 2 at the point A, and Ω = α 2 − β 2 . It is important that α > β.
In this case a drift motion of electrons and holes takes place in the positive direction of the x 1 axis, from the point A to the point B (see Fig. 6). This fact helps to construct POs. We assume that everywhere in a domain, in which we construct asymptotic solutions for electronic eigenfunctions, that the inequality E > U(x 2 ) holds. In equations (75-76) π 1 is a fixed parameter, whereas π 2 and t 0 as functions of π 1 are determined uniquely by the equations f 1 (t 0 , π 1 , π 2 , β) = D, f 2 (t 0 , π 1 , π 2 , β) = 0. The formulas, describing PO holes with Hamiltonian on the level set H = 0, are given by for the upper part 0 < t < t 0 , and for the lower part 0 < t < t 0 . It is worth to remark that holes move along a clockwise direction of PO whereas electrons run counter-clockwise around the PO contour. Thus, we have for electrons and holes the continuous family of POs with respect to parameter π 1 which look like lens-shaped contours. As soon as the parameter π 1 has been determined from the generalized Bohr-Sommerfeld quantization condition (72), the semiclassical energy levels for electrons and holes are computed by For the lens-shaped class of POs the high-energy semiclassical localized eigenstates were tested successfully against the energy eigenvalues and the eigenfunctions computed by the finite element method using COMSOL (see [70]). The boundary conditions were used in the following numerical experiments.
In Fig's. 9 and 10, the electronic eigenfunction density component |u| 2 is shown that was computed semiclassically for the the states m 1 = 27, m 2 = 0, 1, 2 with E = 2.2538, 2.2812, 2.3078 for h = 0.025, α = 1.0, β = 0.5. It is worth to remark that in this case the localization of eigenfunction density components takes place in close neighborhood of PO. In all shown figures computed by semiclassical analysis one can easily see a white contour of PO. In Fig. 9 (a) dependence of TrM/2 on π 1 -(a) and dependence of Im(Γ) on s -(b), for the state m 1 = 27, m 2 = 0 with E = 2.2538 α = 1, β = 0.5, h = 0.025 are shown.

Conclusion
In this review we have outlined our work on the semiclassical analysis of graphene structures and introduced some new results for monolayer and bilayer graphene. We have outlined a range of new asymptotic methods and a semiclassical analysis of Dirac electron-hole tunneling through a Gaussian shaped barrier that represents an electrostatic potential. We started by analyzing the rectangular barrier and have found some important differences between it and the smooth barrier. Namely, the smooth barrier exhibits a quasi-discrete spectrum and complex bound states that do not exist for the rectangular barrier (in zone "3" in Fig. 3). In both types of barrier Klein tunneling occurs. The WKB approximation deals with the asymptotic analysis of matched asymptotic techniques and boundary layers for the turning points in the barrier. The main results of this work are eloquent WKB formulas for the entries in a smooth barrier transfer matrix. This matrix explains the mechanism of total transmission through the barrier for some resonance values of energy of skew incident electrons or holes. Moreover, it has been shown that the existence of modes localized within the barrier, and exponentially decaying away from it, for two discrete complex and real sets of energy eigenlevels can be determined by the Bohr-Sommerfeld quantization condition. It was shown that the total transmission through the barrier takes place when the energy of an incident electron or hole coincides with a real part of the complex energy eigenlevel of one among the set of localized modes. These facts were confirmed by numerical simulations done by finite elements methods and have been found to also be in excellent agreement with the results found using finite difference methods as in [71].
We have also applied the Gaussian beam methods, originated by Popov [51] and expanded by Zalipaev [66] to quantum problems, to describe monolayer and bilayer graphene. We have constructed eigenfunctions and defined the stable periodic orbit conditions and the quantization conditions. The reflection and transmission coefficients of monolayer and bilayer graphene have been derived within the context of semiclassical physics in full. It is clear that these methods can offer good insights into the behavior of the graphene Fabry-Perot resonator.
Such systems will find applications in plasmonics, and nanoribbon heterostructures made from graphene are promising to emerge. The kind of bilayer structure analyzed here can be created by chemical vapor deposition [32] and this opens up the road to a flurry of geometrically optimized graphene resonator systems, whether acting in isolation or as part of a composite, or array.

Appendix A. Transfer and scattering matrix properties for a smooth step
Let us formulate this scattering problem in terms of transfer matrix T for the left slope of the entire barrier (see [67]) and Taking into account the conservation of the x-component of the probability density current (see equation (8) in [17] or (18) in [18]) we obtain that Thus, for the slope transfer matrix T it holds that or As a result we have |T 11 | = |T 22 |, |T 12 | = |T 21 |, |det(T)| = 1. For the scattering matrix we have and From (90) we obtain that thus, the scattering matrix is unitary. If the entries of S are known, then, T = −S 11 /S 12 1/S 12 S 21 − S 11 S 22 /S 12 S 22 /S 12 , Time-reversal symmetry in scattering through the graphene barrier would mean that are both asymptotic solutions to the Dirac system, and In what follows that Thus, S 12 = S 21 , and If a 1 = 1, a 2 = r 1 , d 1 = 0, d 2 = t 1 , then If a 1 = 0, a 2 = t 2 , d 1 = 1, d 2 = r 2 , then

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.
The definition of T is given by (81),(82), and looks the same Ta = d. However, for the barrier we have and For the scattering matrix S we have and From (105) we obtain that If the entries of S are known, then, Taking into account the time-reversal symmetry in scattering through the graphene barrier, we obtain S = S T , and 16. Appendix C. WKB asymptotic solution for tunneling through a smooth step.

Left slope tunneling
Let us assume that E > E c , where E c = |p y | is the cut-off energy. In the case |E| < E c there is no wave transmission through the barrier. On the other side, we assume that E < U 0 − δE, and δE is chosen such as to avoid coalescence of all four turning points. Consider a scattering problem through a smooth step that is the left slope of the barrier. Assume that the right slope in Fig. 1 does not exist, that is U(x) = U 0 if x > x max . In this case we have three domains Ω i , i = 1, 2, 3 with the only difference for Ω 3 = {x : x 2 < x < +∞}. Thus, to the leading order, in the domain Ω 1 we have a superposition of waves traveling to the left and to the right.
In the domain Ω 2 we have exponentially decaying and growing contributions. In the domain Ω 3 we have where d = T L a. It is worth remarking that the electron state in x < a transfers into a hole state for x > a.
To determine the unknown entries of the transfer matrix T L (see Appendix A), we have to match the principal terms of all asymptotic expansions by gluing them through asymptotically small boundary layers at x = x 1 and x = x 2 . To perform matching asymptotics techniques in this case we introduce a new variable U(x) − E = ξ and derive an effective Schrödinger equation. Then, we have where α = α(ξ) = dξ/dx > 0. Changing u, v as follows we obtain the system of −ihαW ξ + ip y V + ξW = 0, and ihαV ξ − ip y W + ξV = 0.
Next, differentiating the first equation with respect to ξ and eliminating V, we obtain a second order ODE for W Then, after we have found W, we have Both boundary layers for two turning points ξ = −|p y | and ξ = |p y | are determined by following scale, well-known in WKB asymptotics for turning points in 1D Schrödinger equation as h → 0 (see for example [16]), On the other side, this scattering problem for the equation (112) written as effective Schrödinger equation may be represented to leading order as follows w = 1 2(ξ 2 − p 2 y ) 1/4 |p y |α 2 a 1 sgn(p y ) e i h Φ − (ξ) for ξ < −|p y |, where for ξ > |p y |. According to the method of comparison equations described in [61] and [62], we seek asymptotic solutions, uniform with respect to |p y |, as follows and the function z(ξ) is determined by z ′2 (a 2 − z 2 /4) = q(ξ). where The asymptotics include the parabolic cylinder function D ν (z) that is a solution to h 2 y zz + (h(ν + 1/2) − z 2 /4)y = 0.

Right slope tunneling
Now let us formulate the scattering problem with transfer matrix T R for the right slope. Taking into account that α = | dξ dx |, the problem for transfer matrix written in terms of solution to the effective Schrödinger equation may be represented as follows for ξ > |p y |, for ξ < −|p y |. If w is a solution to (138), then w * is the solution to (114). Thus, the coefficients from (139), (140) are connected by Hence, we have Since d = T R a, we obtain where It is worth to remark that Q 1 and Q 2 differ as the function α(ξ) behave differently for the same segment ξ → (−|p y |, |p y |) for left and right slopes of non-symmetric barrier.