Klein tunneling through the trapezoidal potential barrier in graphene: conductance and shot noise
KKlein tunneling through the trapezoidal potential barrier in graphene: conductanceand shot noise
Gheorghe Sorin Paraoanu
QTF Centre of Excellence, Department of Applied Physics,Aalto University School of Science, P.O. Box 15100, FI-00076 AALTO, Finland
When a single-layer graphene sheet is contacted with metallic electrodes, tunnel barriers areformed as a result of the doping of graphene by the metal in the contact region. If the Fermi energylevel is modulated by a gate voltage, the phenomenon of Klein tunneling results in specific featuresin the conductance and noise. Here we obtain analytically exact solutions for the transmission andreflection probability amplitudes using a trapezoidal potential barrier, allowing us to calculate thedifferential conductance and the Fano factor for a graphene sheet in the ballistic regime. We putin evidence an unexpected global symmetry - the transmission probability is the same for energiessymmetric with respect to half of the barrier height. We outline a proposal for the experimentalverification of these ideas using realistic sample parameters.
Graphene is a flat monolayer crystal formed of carbon atoms arranged into a two-dimensional (2D) honeycomb lattice[1, 2]. This material provides a rich experimental platform for testing the predictions of relativistic quantum fieldtheory, due to the fact that its low-energy excitations are massless Dirac fermions [3]. To study these quasiparticles,metallic contacts are fabricated on the graphene sheet, either as backgates or topgates. This results in the formationof tunnel barriers. Since contacts are essential for any electronics application, it is imperative to study the transportof electrons in this type of circuit.In non-relativistic quantum mechanics, tunneling is realized through evanescent waves and consequently the trans-mission probability is lower than unity, being in general dependent on the width and height of the barrier. Incontradistinction, as a consequence of the peculiar energy dispersion, in graphene one expects the occurence of Kleintunneling [4, 5] – where the probability of the particle exiting out on the other side of the barrier can be unity,irrespective to the height and width of the barrier. However, although in non-relativistic quantum mechanics thecalculation of tunneling amplitude probabilities under various barriers is a straightforward exercise, in the relativisticcase solutions have been provided only in some limiting or asymptotic cases. Relatively simple solutions have beenobtained for the case of infinite-slope step potential [6–9] as well as for smooth step potentials [10, 11]. It was alsonoticed that there exists a coincidence between the Fano factor obtained from Klein tunneling near the Dirac pointfor rectangular-barrier models [12] and the value of 0.3 corresponding to diffusion models in disordered metals [13].For the case of linear-slope potentials, Sonin has obtained expressions for transmission and reflection under the as-sumption of an infinitely deep potential under the metallic contacts [14, 15], by building upon the solutions of theDirac equation in an uniform electric field discovered by Sauter [16]. However, in all these works, the contacts arenever explicitly considered, and the wavefunction of the incoming electrons is assumed to extend asymptotically theslope solution. This is obviously not the case in real samples, where the role of contact doping cannot be ignored[17–20].Experimental evidence of Klein tunneling has started first to emerge from resistance measurements of npn graphenesamples at various backgate/topgate voltages [21–24]. Noise data have been acquired shortly after, and they confirmedthe basic phenomenology. Shot noise measurements on single-layer graphene samples, have been first reported indiffusive samples [25] (Fano factor of 0.35-0.38) and ballistic samples [26] (maximum Fano factor of 0.318-0.338).The latter results have confirmed the expected theoretical models based on evanescent waves [12] and large width-to-length ratio (
W/L (cid:29) e / ( π (cid:126) ) and a universal maximumFano factor of 1/3 [7, 12]. The effect of decreasing W/L has been further studied systematically in ballistic-sampleexperiments with graphene flakes on SiO [26–28]. More recently, the asymmetry of shot noise has been thoroughlyinvestigated in samples of suspended graphene [29]. Information about doping can be extracted not only from contactresistance measurements, but also from noise measurements combined with measurements of differential conductance[29].Here we fully solve the paradigmatic case of Klein tunneling under a trapezoidal barrier for propagating waves,providing general analytical results. Based on these, all previous approximate results from the literature can beobtained by taking the appropriate limits. Our results can be readily used to analyze the data in any tunneling ornoise experiment: in order to demonstrate this, we develop a simple and useful electrical model of doping and we showhow to connect it to our parametrization. The solutions provided here not only confirm the previous experimentalresults on tunneling in graphene, but also offer a unified methodology for the design of electronic circuits based ongraphene. As the complexity of these circuits increases, with the use of multi-terminals, different doping substrates,and several sheets of graphene, the earlier approaches to tunneling based on asymptotic approximations will be a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b FIG. 1: Schematic of a typical graphene sample and of the resulting electric potentials. insufficient. Instead, our approach based on analytically-exact solutions provides a straightforward, scalable approachto more sophisticated graphene-based circuits.The paper is organized as follows. In Section I we introduce the problem, the main equations, and the keynotations. In Section II we calculate the transmission and reflection probability amplitudes and we put in evidencetheir symmetries. In Section III we proceed to obtain the conductance and Fano factor for the incoherent and coherentcases. Finally, in Section IV we introduce the doping model and provide experimental predictions that can be testedusing present-day technology. We conclude with an overview of the main results in Section V.
I. TUNNELING UNDER THE TRAPEZOIDAL BARRIER: GENERAL CONSIDERATIONS
A typical graphene sample is presented schematically in Fig. I. The graphene is supported by a dielectric layer (orit can be suspended), and contacted by two metallic leads. A gate voltage can be applied either to the entire graphenelayer (typically as a backgate) or locally (as a top gate). The differential conductance is measured in a standardlock-in setup, with a small bias voltage applied across the structure ib order to obtain a measurable current.Our goal is to provide a complete solution to tunneling in graphene with a trapezoidal-shape potential barrier. Wepresent analytically exact results for the transmission and reflection coefficients, and we calculate the conductanceand Fano factor for realistic experimental setups. The model employed extends the analysis of Sonin [15] to the moregeneral case where the level of doping of graphene under the electrodes is not necessarily very large. It assumes anelectrostatic barrier of trapezoidal form, comprising five distinct regions: the middle region 0 ≤ x ≤ L with constantdoping potential V , two slope regions − d ≤ x ≤ L ≤ x ≤ L + d , where the electric field is finite, and tworegions with zero potential, x ≤ − d and L + d ≤ x . Because we are interested in electrons incident on barriers, thepotential is taken y -independent. This potential is also assumed to vary smoothly on the scale given by the graphenelattice constant of approximately 0 . K and K (cid:48) ,which are separated by | K − K (cid:48) | (of the order of inverse lattice constant [4]). In the left slope region − d ≤ x ≤ V ( x ) = V + sgn [V ] a (cid:126) v F e x, (1)where v F is the Fermi velocity in graphene. The sign sgn [V ] realizes either a positive or negative slope. For clarity, a will denote the modulus of the slope and will therefore be always a positive quantity. The potential Eq. (1) resultsin a Hamiltonian that is translationally invariant in the y -direction, which, as detailed below, allows us to solve thecorresponding Dirac equation by the method of separation of variables. A. Dirac equation in an external potential
Near the Dirac point (valley) K the dispersion relation is linear, and as a result the electrons are described by theequation [ − i (cid:126) v F σ · ∇ + eV ] Ψ AB = E Ψ AB , (2)where σ = ( σ x , σ y ) is a vector of Pauli matrices and v F ≈ m / s ≈ c/ V ( x ), where the wavefunction Ψ AB ( x, y ) = ( ψ A ( x, y ) , ψ B ( x, y )) T is a pseudospin 1/2 corresponding tothe two triangular sublatices A and B, and T stands for transpose. The energy E of the Dirac particle is measuredfrom the same level as the potential V ( x ). Introducing the momentum operator ˆ p = − i (cid:126) ∇ , we can alternativelywrite this equation as ˆ H Ψ AB ( x, y ) = E Ψ AB ( x, y ), where the Hamiltonian is ˆ H = ˆ p · σ + eV . For clarity, in orderto distinguish them from the pseudospin operators, we reserve the hat notation for operators acting on continuousvariables. Since the potential V is invariant in the y -direction, the corresponding y component of the momentum isconserved, i (cid:126) ˙ˆ p y = [ˆ p y , H ] = 0, and we can write its eigenvalue as p y = (cid:126) k y , where k y is the initial-state wavevector.Also the velocity operator is a Pauli matrix, since ˙ˆ r = − ( i/ (cid:126) )[ˆ r , ˆ H ] = σ , where r = ( x, y ). The expression for theelectrical current can be found from the continuity equation as [4] j = ev F ψ † AB σ ψ AB , (3)and owing to the time-independence of the problem we have ∇ j = 0.To proceed, it is convenient to perform a change of basis in Eq. (2), implemented by a rotation R y ( − π/
2) =exp( iσ y π/
4) = ( I + iσ y ) / √ π/ y axis. We then have R y ( − π/ σ x R † y ( − π/
2) = σ z , and weintroduce Ψ( x, y ) = R y ( − π/ AB ( x, y ). In this rotated frame, the Dirac equation reads( − i∂ x σ z − i∂ y σ y ) Ψ( x, y ) = K ( x )Ψ( x, y ) , (4)where we introduced a kernel K ( x ) = [ E − eV ( x )] / (cid:126) v F , and E is the energy of the Dirac particle. Here ψ + ( x, y ) and ψ − ( x, y ) are the two components of the spinor Ψ( x, y ) defined by Ψ( x, y ) = (Ψ + ( x, y ) , Ψ − ( x, y )) T .Furthermore, in this rotating frame the expression for the current in the x direction as obtained from Eq. (3)becomes j x = ev F Ψ † σ z Ψ = ev F ( | ψ + | − | ψ − | ) . (5)Using this equation it immediately follows that the x -component j x of the current is constant along both x - and y - directions. The constancy in the y direction ∂ y j x = 0 is not surprising since it follows immediately from the y -invariance of the potential. However, the condition ∂ x j x = 0 (6)is not trivial. It implies that the x-current is the same in all the regions, a condition which we will verify explicitlythroughout the paper. For convenience, we will work with solutions satisfying | ψ + | − | ψ − | = 1 or j x = ev F ,corresponding to a particle of charge e being transferred across the barrier from the left to right at the Fermi velocity v F .From Eq. (4) we can also put in evidence the three symmetries of the problem:1. Sublattice symmetry:
If Ψ is a solution then σ x Ψ is a solution for the same problem but with a kernel K ofopposite sign. This amounts to the substitutions ψ ± → ψ ∓ and K → −K . In terms of the original lattice, thismeans that ψ A → − ψ A and ψ B → ψ B produces a solution corresponding to a change of sign in energy E → − E , andpotential eV → − eV . The current j changes sign under this transformation.2. Single-valley time-reversal symmetry:
A solution of the complex conjugate of Eq. (4) can be obtained byΨ ∗ → σ y Ψ (or ψ ∗± → ∓ iψ ∓ ) and the same kernel K . This can be seen also from taking the complex conjugate ofthe Hamiltonian generating Eq. (4), ( − i∂ x σ z − i∂ y σ y ) ∗ = σ y ( − i∂ x σ z − i∂ y σ y ) σ y . In terms of the original lattice thismeans ψ ∗ A → − iψ B , ψ ∗ B → iψ A . The current j changes sign by this transformation, as expected from a time-reversedprocess.3. Particle-hole symmetry:
It is obtained by combining the time-reversal with the sublattice symmetry. A solutionof the complex conjugate of Eq. (4) can be obtained by Ψ ∗ → σ z Ψ (or ψ ∗± → ± ψ ± ) and K → −K . This can be seenalso from taking the complex conjugate of the Hamiltonian generating Eq. (4), ( − i∂ x σ z − i∂ y σ y ) ∗ = − σ z ( − i∂ x σ z − i∂ y σ y ) σ z . In terms of the original pseudospin associated with the sublattice A and B this transformation reads ψ ∗ A → ψ B and ψ ∗ B → ψ A . The current j remains invariant under this symmetry.Since in our system the kernel K depends only on x , we can search for a solution of Eq. (4) by writing ψ ± ( x, y ) = ψ ± ( x ) exp( ik y y ). The plane wave in the y -direction is appropriate for the subsequent calculation of the conductanceand Fano factors, for which the standard procedure employs periodic boundary conditions. This results in i ddx Ψ( x ) = −K ( x ) σ z Ψ( x ) − ik y σ x Ψ( x ) . (7)Note at this point that the dimension of a is [length] − , while that of K ( x ) and of k y is [length] − . We have thusobtained a fully adimensional Dirac equation Eq. (7), with all the physical constants embedded in the kernel K .Next, we introduce a few other useful notations for the height of the Fermi level measured from the Dirac pointfor the graphene under the metal, which we denote as ∆ E FM , and the same quantity for the free graphene, which wedenote by ∆ E FG , see Fig. I. We then define adimensional quantities [14, 29], namely the reduced coordinate ξ = x √ a ,the reduced wave-vector component in the y -direction κ = k y / √ a , and the reduced Fermi momentum for the freegraphene, v = ∆ E FG √ a (cid:126) v F . (8)Also we introduce the height of the potential barrier by eV = ∆ E FM − ∆ E FG , (9)and the finite-slope impact parameter p = eV (cid:126) v F √ a . (10)It is also natural to define a reduced kernel K ( ξ ) = K ( x ) / √ a . We can now express Eq. (7) in the variable ξ on the ± components i ddξ ψ + ( ξ ) + iκψ − ( ξ ) = − K ( ξ ) ψ + ( ξ ) , (11) i ddξ ψ − ( ξ ) + iκψ + ( ξ ) = K ( ξ ) ψ − ( ξ ) , (12)obtaining two equations fully written in terms of the reduced parameters. II. TRANSMISSION COEFFICIENTS
With the notations and setup above we are now in a position to solve analytically the problem of tunneling thoughpropagating modes of the Dirac particle across the trapezoidal barrier.
A. General expressions Top of the barrier, x > In this region the reduced kernel reads K ( ξ ) ξ> = v, (13)for both positive and negative ∆ E FG . The solution for a particle moving in the direction of the positive x axis isΨ( ξ ) x> = (cid:20) ψ + ( v, κ ) ψ − ( v, κ ) (cid:21) e i sgn[ v ] √ v − κ ξ , (14)where ψ ± ( v, κ ) = ±| v | + √ v − κ + i sgn[ v ] κ (cid:113) | v |√ v − κ . (15)It is straightforward to calculate the j x current, and we find from the equations above and employing Eq. (5) that j x = ev F . We note that in this region the momentum in the x -direction is p x, FG = sgn[∆ E FG ] (cid:113) (∆ E FG /v F ) − (cid:126) k y and it is a real quantity. This implies | ∆ E FG | ≥ (cid:126) v F k y , or | v | ≥ κ . Region under the metal x < − d For the region of graphene under the metal we have K ( ξ ) x< − d = p + v. (16)The solution can be expressed as a superposition between a transmitted and a reflected waveΨ( x ) x< − d = 1 t Ψ ( t ) ( x ) + rt Ψ ( r ) ( x ) , (17)with r and t denoting respectively the reflection and transmission coefficients.In terms of the variable ξ we write these waves asΨ ( t ) ( ξ ) = (cid:20) ψ + ( p + v, κ ) ψ − ( p + v, κ ) (cid:21) e + i sgn[ p + v ] √ ( p + v ) − κ ξ , (18)and Ψ ( r ) ( ξ ) = − (cid:20) ψ ∗− ( p + v, κ ) ψ ∗ + ( p + v, κ ) (cid:21) e − i sgn[ p + v ] √ ( p + v ) − κ ξ . (19)Note also that in Eq. (19) the reflected component is normalized to a negative current in the x-direction equal to − ev F while the transmitted component is normalized to a positive + ev F . Specifically, j ( t ) x = − j ( r ) x = ev F . (20)Moreover, one can explicitly check that the current corresponding to the total wavefunction Eq. (17) is + ev F , that is j x = ev F , (21)given that that | r | + | t | = 1. Thus the correct normalization is recovered in this region. Similarly to theprevious region, the condition | ∆ E FM | ≥ (cid:126) v F k y must hold, ensuring that the momentum in the x -direction p x, FM = sgn[∆ E FM ] (cid:113) (∆ E FM /v F ) − (cid:126) k y is real. This results in | p + v | > κ . Slope of the barrier, − d < x < In this region we have K ( ξ ) − d We now look at the situation when the slope of the potentialis positive at the left side of the barrier. As a convention, in the following a will be defined as a positive number, a > 0. The height of the electrical potential is positive, V ≡ ∆ E FM − ∆ E FG > 0. A schematic of the potentialexperienced by the Dirac fermions is shown in Fig. 2. The type of junction formed depends on the carrier type inthe “metal” M and “graphene” G region. We get a n − n junction if ∆ E FM > E FG > 0, a n − p junctionif ∆ E FM > E FG < 0, and finally a p − p junction if ∆ E FM < E FG < 0. Note that p − n junctionscannot form for positive slopes. The potential has the form V ( x ) = V + a (cid:126) v F e x, (26)resulting in the useful relation √ ad = p , (27)when V ( − d ) = 0 is used. The kernel on the slope of the barrier is given by K ( x ) = ∆ E FM − eV ( x ) (cid:126) v F = ∆ E FG (cid:126) v F − ax = − a ( x − x (+)0 ) , (28)where the crossing point is x (cid:37) = ∆ E FG /a (cid:126) v F = v/ √ a , or ξ (cid:37) = √ ax (cid:37) = v , where the superscript (cid:37) denotes apositive slope. Thus the crossing point takes positive values for ∆ E FG > E FG < (cid:20) ψ (cid:37) + ( ξ ) ψ (cid:37)− ( ξ ) (cid:21) − d The case of a potential with negative slope has to betreated separated. The reason for this is that the slope appears under square root in the substitutions needed to getthe Sauter solution. With our convention, a is defined as a positive number, a > 0, and the height of the electricalpotential is become negative, V ≡ ∆ E FM − ∆ E FG < 0. A schematic of the potential experienced by the Dirac fermionsis shown in Fig. 3. The type of junction formed depends on the carrier type in the “metal” M and “graphene” Gregion, which is set by the value of the Fermi level. We get a n − n junction if ∆ E FM > E FG > 0, a p − n junction if ∆ E FM < E FG > 0, and finally a p − p junction if ∆ E FM < E FG < 0. For negative slopes, n − p junctions cannot form. The potential has the expression V ( x ) = V − a (cid:126) v F e x, (33)resulting in √ ad = − p . (34)For the kernel on the slope of the barrier we have K ( x ) = ∆ E FM − eV ( x ) (cid:126) v F = ∆ E FG (cid:126) v F + ax = a ( x − x (cid:38) ) , (35)where the crossing point is x (cid:38) = − ∆ E FG /a (cid:126) v F , or ξ (cid:38) = −√ ax (cid:38) = − v , with the superscript (cid:38) denoting a negativeslope. Opposite to the situation of positive slope, here the crossing point has negative values for positive ∆ E FG > E FG < (cid:20) ψ (cid:38) + ( ξ ) ψ (cid:38)− ( ξ ) (cid:21) − d From the results shown in Fig. 4 we can already see a number on interesting features, which are detailed below.First, note that propagation is possible only for κ satisfying − v < κ < v and − p − v < κ < p + v , which implies zerotransmission in regions that do not satisfy these conditions. Geometrically, the wave incident on the barrier in the metalregion forms an angle φ M with the normal, and the transmitted wave in the graphene region forms an angle φ G with thenormal. Here tan φ M = (cid:126) k y /p x, FM = κ/ (sgn[ p + v ]( (cid:112) ( p + v ) − κ ) and tan φ G = (cid:126) k y /p x, GM = κ/ (sgn[ v ] √ v − κ ).We have then an analogue of Snell’s law, v sin φ G = ( p + v ) sin φ M . We notice that the extreme values p + v = ± k , v = ± κ correspond to grazing incidence and transmission, that is φ M = π/ φ G = π/ a. Klein tunneling and pseudo-spin conservation We start by discussing the interesting case in which κ = 0, i.e. the particle is incident on the barrier in the perpendicular x -direction. By analysing the final expression fortransmission or reflection Eqs. (31,32) and Eq. (39,40), we can find immediately that the reflection is zero in thiscase, irrespective to the height or slope of the barrier. This is in contradistinction with the case of tunneling innon-relativistic quantum mechanics. It is instructive to understand what is the origin of this effect.For k = 0, the wavefunction under the metal readsΨ( ξ ) x 0) = exp (cid:0) − i ( ξ − ξ ) / (cid:1) and G ( ξ − ξ , 0) = 0 with solutionsΨ (cid:37) ( ξ ) − d 3. b)Transmission at the left side of the barrier for negative slope and impact parameter p = − . b. Global symmetries The special symmetry of the trapezoidal potential with respect to half of its height leads totwo additional symmetries of this problem. Note that while time-reversal and particle-hole are local symmetries, validfor a given x , the symmetries unearthed below are global symmetries, arising due to the shape of the solution acrossthe slope and the connection points. In terms of energies, this amounts for the substitution p / v ↔ − p / − v or v ↔ − p − v .The first such symmetry is the following one:1 t (cid:37) ( − p − v, κ ) e i sgn[ v ] √ v − κ = 1 t (cid:38) ( v, κ ) e i sgn[ p + v ] √ ( p + v ) − κ , (45) r (cid:37) ( − p − v, κ ) t (cid:37) ( − p − v, κ ) e − i sgn[ v ] √ v − κ = − r (cid:38)∗ ( v, κ ) t (cid:38)∗ ( v, κ ) e i sgn[ p + v ] √ ( p + v ) − κ . (46)Using the expressions for the left side of the barrier, one can immediately obtain the corresponding expressions forthe right side. For example, a particle impinging on the barrier from the left side will experience on the right side astep with the opposite sign. Suppose that the slope at the left side is p : then to find the transmission coefficient atthe right side it is enough to make the substitution v → p + v or v − p → v in the negative-slope formula (if theslope at right side is negative) or in the positive-slope formula (if the slope at the right side is positive).The second symmetry concerns the absolute value of the transmission coefficient, namely in both the case ofpositive slope and in that of negative slope the transmission coefficient is symmetric with respect to the point − p / | t ( v, κ ) | is invariant under the transformation v → − p − v , in other words | t ( v, κ ) | = | t ( − p − v, κ ) | and by consequence of normalization | r ( v, κ ) | = | r ( − p − v, κ ) | . We emphasize that thisproperty is non-trivial: it not a direct consequence of one particular local symmetry of the problem but rather aconsequence of the global inversion symmetry of the potential V ( x ) with respect to the point where the energy − p / t , withthe exception of the case κ = 0 analysed in detail above. In the latter situation, this happens because the phase isanyway independent on v : it is − p / p / III. CONDUCTANCE AND FANO FACTOR FOR THE WHOLE BARRIER With state-of-the-art nanofabrication of graphene junctions, it is possible to realize samples that are in the ballisticregime in the region of interest. Next, we work out in detail how the conductance and Fano factor can be calculated,allowing a straightforward comparison with experiments.0As established above, an incident wave coming from the left of the trapezoidal barrier will be partly reflected, with acomplex reflection coefficient r , and partly transmitted, with transmission coefficient t . These coefficients will dependon p , v and κ , r = r ( v, κ ), t = t ( v, κ ), with each value of k y corresponding to a conduction channel. Therefore, theoverall conductance and noise is a sum of the conductances and respectively noise per channel.Specifically, the conductance is obtained as σ ( v ) = σ (cid:90) min {| p + v | , | v |} dκT ( v, κ ) , (47)and the differential noise reads s ( v ) = σ (cid:90) min {| p + v | , | v |} dκ [1 − T ( v, κ )] T ( v, κ ) . (48)From this, we get the Fano factor as F ( v ) = s ( v ) σ ( v ) , (49)where σ = 4 e w √ a/πh = 2 G w √ a/π , with w the width of the sample and G = 2 e /h = 7 . × − Ω − thequantum of conductance. The factor of 2 that multiplies G in the expression of σ is due to the double-valleydegeneracy of graphene. Note also the upper limits of the integrals in Eqs. (47,49), which come from the conditionof free propagation at the top of the barrier κ < | v | , and at the left side of the barrier κ < | p + v | .An important observation concerns the role played by the barrier impact factor p = √ ad as a fitting parameter.For a given sample geometry and for a certain V g , the Fano factor depends only on p . The gate voltage V g completelydetermines the value of ∆ E FG and of v = k FG / √ a . Then only p remains in the expression of the Fano factor. Notethat p = √ ad , and exactly this quantity appears in the arguments of the hypergeometric functions when we calculate ξ ( − d ). This is not the case for σ – there the integral would depend only on p but σ contains another dependenceon the slope a. A. Ballistic incoherent case In the case of incoherent tunneling, the length of the sample is larger than the coherence length and as a result thephase factors appearing in the transmission is lost as the particle propagates on the top of the graphene barrier. Thisresults in a transmission probability T ( v, κ ) = 12 | t ( v, κ ) | − − , (50)where we have considered a symmetric barrier (same slope at the left and right side) with reduced Fermi level v in theregion G. Using Eqs. (47,49) we can calculate the conductance and the Fano factor for several values of the impactparameter, see Fig. 5 and Fig. 6.A number of interesting features can be noticed from the plots. For large impact parameter | p | (cid:29) | p | (cid:38) | v | smaller that | p | . For | p | = 10 we get the saturation conductance σ/σ ≈ . F ≈ . 33. As the parameter p decreases, the conductance at low voltages tends to decrease, due to the fact that thetransmission decreases. We no longer observe a saturation phenomenon, but instead, in the range −| p | < | v | < −| p | the maximum values for the conductance and Fano factors are reached at v = − p / 2. For values of | p | below 1, themaximum value of the Fano factor is ≈ . 25, and it stays almost independent of p . The slope of the conductanceoutside the range −| p | < | v | < −| p | as well as the value of the Fano factor at v = 0 and v = p also take universalvalues, independent on p . We can find these values by looking at 0 < v (cid:28) p > < − v (cid:28) p < 0, and using F (0 , κ ) = 1 and G (0 , κ ) = 1, to get | t ( v, κ ) | ≈ √ v − κ | v | + √ v − κ , (51)and therefore T ≈ √ v − κ | v | . (52)1 FIG. 5: (color online) Conductance a) and Fano factor b) for positive slope at the left of the barrier. The green lines correspondto p = 10, the blue ones to p = 1 . 3, and the red ones to p = 0 . p = − 10, the blue ones to p = − . 3, and the red ones to p = − . Inserting this expression into Eq. (47, 49) we obtain by direct integration that1 σ dσd | v | = π . , (53)and F ( v = 0) = 1 − π = 0 . . (54) B. Ballistic coherent case If the propagation in the graphene sheet is coherent, then Fabry-P´erot interference fringes will appear in theconductance and in the Fano factor due to interference between multiple reflections at the two slopes of the barrier.The formula for the total transmission in this case is T = | t | | r | − | r | cos[2 ξ L √ v − κ ] , (55)where we used the fact that the wave propagating from the left slope to the right slope in the non-contact regionpicks up a phase √ v − κ ξ L , where ξ L is the reduced length of the top of the barrier, ξ L = √ aL . This result can be2 FIG. 7: (color online) Conductance a) and Fano factor b) for positive slope at the left of the barrier in the coherent case. Thegreen lines correspond to p = 10, the blue ones to p = 1 . 3, and the red ones to p = 0 . employed now in Eq. (47, 49) to calculate the conductance and the Fano factor. As an example, in the case of positiveslope at the left side of the barier the conductance and the Fano factor are plotted in Fig. 7 with ξ L = √ aL = 9.For ballistic and coherent graphene samples, the modulation of conductance and Fano factor around the Dirac point v = 0 has been already observed experimentally in Fabry-P´erot experiments [33]. One notices that in the conductancethe modulation due to multiple reflections is not so pronounced for v > 0. This is another manifestation of the Kleintunneling phenomenon, and it is due to the fact that in that region the transmission is large, therefore there existscomparatively little reflection with respect to the negative- v region. In Eq. (55) the effect of the modulation withrespect to ξ L thus tends to be comparatively small and it is further washed out by the integration over κ . Thus theFabry-P´erot effect tends to be suppressed, and only the Fano factor retains some mild oscillatory dependence. IV. EXPERIMENTAL PROPOSAL To test these prediction, a rather standard experimental setup, together with high-quality clean graphene andmetal-graphene contacts is needed. These are available nowadays, owing to progress in the chemical vapor depositionsynthesis of graphene on metal surfaces [30] and novel techniques for achieving a controlled small contact resistance[31]. The ballistic regime, with electrons having long mean-free paths compared to the sample size, can be achievedby fabrication methods for suspended graphene using either conventional wet etching of the SiO substrate [32] or bydry transfer [33]. Top gating is also possible, for example by using a solid polymer electrolyte [34]. These methodsensure that the electrons have mean free paths larger than the dimensions of the sample [24]. Several calculations haveconfirmed that the details of the contact do not alter significantly the ballistic regime and only slight modificationsof the values of the conductance and Fano factor are to be expected [35, 36]. We have also verified (see AppendixB) that our experimental predictions are largely insensitive to the details of tunneling between metal and graphene.This makes our results robust against uncontrolled fabrication imperfections.We introduce next a doping model that allows us to determine the renormalized Fermi level v and the impactparameter p from the geometric characteristics of the sample and the applied gate and bias voltages. The notationsused are as follows: the workfunction of the pristine graphene is denoted by W G , the workfunction for the metal-covered region is W M , the workfunction of the backgate is W g . The chemical shift [20] of the workfunction of grapheneat the contact with the metal this workfunction is denoted by ∆ c .With these notations, we can now find the shifts of the Fermi level in the region under the metal ∆ E FM and in thenon-contact region ∆ E FG as a function of the applied gate voltage V g . In Fig. 8 b) we present a full equivalent electricschematic of the circuit. From this, all the equations above (and combinations thereof) can be obtained by applyingKichhoff’s laws for electrical circuits. The quantum capacitances corresponding to the graphene under the metal andthe suspended graphene are represented by circles containing the Dirac cone. The differential conductance and thenoise can be measured using a lock-in technique, with a small excitation voltage of a couple of tens of mV [29]. Thisis much smaller compared to the other voltages and energy scales in this circuit, and therefore will be neglected inthe following.3 FIG. 8: a) Schematic of Fermi levels for the doping problem. b) Equivalent electrical circuit corresponding to the electrochemicalpotentials. The circuit equations can be obtained from Kirchhoff’s laws along three independent contours, eV g = W g + ∆ c − W M + eU g M − eU c , (56)0 = ∆ E FM + eU c − W G − ∆ c + W M , (57) eV g = eU g G + ∆ E FG + W g − W G . (58)Here U c /e = n c /C c is the voltage across the surface capacitance C c , with n c the number of electrons per unit areain the contact region, and U g M /e = n g M /C g M is the voltage across C g M , corresponding to an electron distribution n g M on the backgate. In the non-contact region, the surface distribution of electrons n g G produces a voltage drop U g G = en g G /C g G . Denoting the surface electron density in the graphene layer in the metal region by n M and thesurface density on the graphene side by n G , the usual laws of charge conservation yield n M = n g M + n c , (59) n G = n g G . (60)The existence of excess charge carriers in the regions M and G can be related to the values of the correspondingFermi levels by the relations ∆ E FM = (cid:126) v F sgn [n M ] (cid:112) π | n M | , (61)∆ E FG = (cid:126) v F sgn [n G ] (cid:112) π | n G | . (62)From Eqs. (56-58), Eqs. (59-60) and Eqs. (61-62) we can get now expressions for the Fermi level shifts∆ E FM = sgn (cid:20) δV g + χC c eC g M (cid:21) − C c + C g M ζ + (cid:115)(cid:18) C c + C g M ζ (cid:19) + ζ C c (cid:12)(cid:12)(cid:12)(cid:12) χ + C g M C c eδV g (cid:12)(cid:12)(cid:12)(cid:12) , (63)∆ E FG = sgn [ δV g ] − C gG ζ + (cid:115)(cid:18) C gG ζ (cid:19) + eC gG ζ | δV g | . (64)which are directly connected to the experimentally-relevant quantities. Here ζ F = √ π (cid:126) v F /e is a constant called Fermielectric flux [29], and the value v F = 1 . × cm/s for the Fermi velocity of graphene yields (cid:126) v F = 0 . 724 eV · nm and4 FIG. 9: a) Conductance and b) Fano factor for the sample discussed in the text. For this numerical simulation we took d = 10 nm, χ = 0 . w = 1 µ m. For the capacitances per unit area we took C c = 2 . × − F/cm , C g G = 8 . × − F/cm , and C g M = 3 . × − F/cm . ζ F = 1 . × − V · cm. We also defined δV g = V g − ( W g − W G ) /e , and the difference χ = W G + ∆ c − W M betweenthe workfunction of the graphene in the contact region and the workfunction of the metal.Typically the experiments as well the theory [15] have been focused on observing the Dirac point corresponding tothe suspended graphene region, achieved for ∆ E FG = 0, see Eq. (64). We show here that given a right combination ofsample parameters, also the additional Dirac point of graphene under metal corresponding to ∆ E FM = 0 is accessible.This requires C c χ + C g M eδV g = 0 , (65)see Eq. (63). As a rule of thumb, this means that the contribution of the doping and contact coupling in this equationshould not be too large, such that δV g is still within the measurable range.We can then connect the experimental parameters with the barrier parameter p and the effective Fermi level v by eV = ∆ E FM − ∆ E FG . (66)Using p = eV / (cid:126) v F √ a we obtain v = (cid:18) d (cid:126) v F (cid:19) / ∆ E FG (cid:112) | ∆ E FM − ∆ E FG | = (cid:18) √ πdeξ F (cid:19) / ∆ E FG (cid:112) | ∆ E FM − ∆ E FG | , (67) p = sgn[ V ] (cid:18) d (cid:126) v F (cid:19) / (cid:112) | ∆ E FM − ∆ E FG | = sgn[ V ] (cid:18) √ πdeξ F (cid:19) / (cid:112) | ∆ E FM − ∆ E FG | . (68)The slope of the barrier can be expressed as well in terms of the voltage-dependent ∆ E FM and ∆ E FG , a = 1 (cid:126) v F d | ∆ E FM − ∆ E FG | = √ πeξ F d | ∆ E FM − ∆ E FG | , (69)and we also have p = sgn[ V ] d √ a . We have naturally assumed here that the distance d over which the electricpotential drop V occurs is constant (does not change with the gate voltage). However, we have verified numericallythat if one takes the slope a as being constant, there are only minor quantitative changes in the results. This showsthat our results are robust as they are insensitive to the details of the chemistry of the rearrangement of carriers atthe electrostatic barrier.For the experimental observation of the effects predicted before, the sample can be designed following relativelystandard methods. The choice of the contact metal is very important. It is known that some metals (Co, Ni, Ti, Pd)interact very strongly with graphene, resulting in hybridization between the graphene p z orbitals and the electronic d states in the metal. This chemisorption process destroys the Dirac cone. However, for metals such as Au, Ag, Cu,Al and Pt(111) the interaction is weaker (physisorption) and the cone is preserved [20, 37, 38]. For the latter metals,5it was also found by ab-initio calculations that the dangling bonds from either the unsaturated graphene edge orpoint defects further lowers the contact resistance, especially for Au and Ag [39, 40]. From various experiments andtheoretical considerations it is reasonable to expect a value of d of around 10 nm, and this is the value we use in thesimulations below. The contact capacitance C c can be estimated as C c = 2 . × − F/cm by assuming a equilibriumdistance value of about 3 ˚A between the graphene and the metal, as found from DFT calculations [20] for metals suchas Al, Ag, Cu, Au, and Pt. From the same calculations, the workfunction of uncontacted graphene results as W G = 4 . c ≈ . W M of metals are typically also of the order of 1 eV closer to the pristine graphene,it is reasonable to consider for exemplification a metal-graphene workfunction difference χ = W G + ∆ c − W M = 0 . supports (relative permittivity (cid:15) r = 3 . 9) in the contact region. This results in C g G = 8 . × − F/cm and C g M = 3 . × − F/cm .The numerical calculations for the conductance and the Fano factors are plotted in Fig. 9. The conductanceand the Fano factor appear slightly deformed when compared with the representation in terms of the reduced Fermimomentum v presented in the previous sections, due to the fact that Eqs. (63,64) and Eqs. (67,68) are nonlinear, butthe key features are clearly there. The position of the additional Dirac point is obtained numerically at δV g = − . µ eV (0.2 K) [41], much smaller than the other energies ( e.g. the hopping constant in graphene is 2.8 eVand the workfunctions are of the order of a few eV), and therefore can be neglected. Also the presence of evanescentmodes can be neglected for bias voltages larger than ζ F /L = √ π (cid:126) v F / ( eL ), where L is the length of the sample [14].For L = 1 µ m, this corresponds to voltages of the order of 1 . 283 mV, which is easily satisfied experimentally. V. DISCUSSIONS AND CONCLUSIONS In elementary quantum mechanics, tunneling under the trapezoidal barrier is a paradigmatic problem. Yet to date,in the case of graphene this problems has not been completely solved, with only partial solutions being presented inthe literature. Here we present the full solution of the problem of Klein tunneling though the trapezoidal barrier. Wecalculate the transmission and reflection coefficients and describe their symmetries. By using these coefficients, wecan further proceed and calculate the conductance, the noise, and the Fano factor. For realistic experimental setupswith contacts realized from metals with certain workfunctions, we present a model for doping that allows us to makespecific predictions that can be verified experimentally. VI. ACKNOWLEDGEMENTS I acknowledge inspiring discussions with Edouard Sonin on the symmetries of the problem as well as with PerttiHakonen, Bernard Pla¸cais, and Antti Laitinen on possible experimental verifications. This work was supported bythe “Finnish Center of Excellence in Quantum Technology QTF” (project 312296) of the Academy of Finland. Wealso acknowledge the European Microkelvin Platform (EMP, No. 824109). Appendix A: Derivation of the Sauter solutions For completeness, we provide here an explicit derivation of the solutions for the Dirac equation in the positive sloperegion. From Eq. (7) with the kernel K ( x ) = − sgn[ V ] a ( x − x ) we get for each of the pseudospin components i ddx ψ + ( x ) + ik y ψ − ( x ) = sgn[ V ] a ( x − x ) ψ + ( x ) , (A1) i ddx ψ − ( x ) + ik y ψ + ( x ) = − sgn[ V ] a ( x − x ) ψ − ( x ) . (A2)6By employing the notations κ = k y / √ a , ξ = √ ax , and ξ = √ ax we can further write (cid:20) ddξ + i sgn[ V ]( ξ − ξ ) (cid:21) ψ + ( ξ ) + κψ − ( ξ ) = 0 , (A3) (cid:20) ddξ − i sgn[ V ]( ξ − ξ ) (cid:21) ψ − ( ξ ) + κψ + ( ξ ) = 0 , (A4)see also Eq. (11) and Eq. (12). Combining these two equations we obtain (cid:20) d dξ + ( ξ − ξ ) (cid:21) ψ ± ( ξ ) + (cid:0) ± i sgn[ V ] − κ (cid:1) ψ ± ( ξ ) = 0 . (A5)Let us consider the case of positive slope, sgn[ V ] = +1. First, we will find the solution ψ (cid:37) + ( ξ ). We proceed by makingthe substitution z = i ( ξ − ξ ) and ψ (cid:37) + ( ξ ) = exp[ − z/ h ( z ) into Eq. (A5), obtaining z d dz h ( z ) + (cid:18) − z (cid:19) dhdz + iκ h ( z ) = 0 . (A6)We can recognize this as Kummer’s equation, with one solution denoted by M (cid:0) − iκ / , / , z (cid:1) . This produces asolution for ψ (cid:37) + ( ξ ) as the function F defined in the main text, F ( ξ − ξ ) = e − i ( ξ − ξ ) / M (cid:20) − iκ , , i ( ξ − ξ ) (cid:21) . (A7)In order to find the corresponding spinor component ψ (cid:37)− ( ξ ), which we will denote by G , we can use directly Eq. (A3), G ( ξ − ξ ) = − κ (cid:20) ddξ + i ( ξ − ξ ) (cid:21) F ( ξ − ξ ) . (A8)Using the expression for the derivative of the Kummer function Eq. (13.4.8) from Abramowitz and Stegun [42] ddz M ( a, b, z ) = ab M ( a + 1 , b + 1 , z ) (A9)we get ψ (cid:37)− ( ξ ) = G ( ξ − ξ ) = − κ ( ξ − ξ ) e − i ( ξ − ξ ) / M (cid:20) − iκ , , i ( ξ − ξ ) (cid:21) . (A10)As a consistency check, we can now verify the other equation Eq. (A4) with the help of the recurrence relation 13.4.14from Abramowitz and Stegun [42],( b − M ( a − , b − , z ) = ( b − − z ) M ( a, b, z ) + z ddz M ( a, b, z ) . (A11)We have now identified one solution for the case of positive slope, with pseudospin components ψ (cid:37)± ( ξ ) given byEqs. (A7,A10). The other solution for positive slope can be obtained by using the electron-hole symmetry. If we takethe complex conjugate of Eq. (A3,A4) we find the same system of equations but with ( ψ (cid:37)− ) ∗ and ( ψ (cid:37) + ) ∗ playing theroles of ψ (cid:37) + and ψ (cid:37)− respectively. Thus, we have another independent solution, G ∗ ( ξ − ξ ) for the + component and F ∗ ( ξ − ξ ) for the − component.To conclude, for positive slope we have the general form of the solutions (cid:20) ψ (cid:37) + ( ξ ) ψ (cid:37)− ( ξ ) (cid:21) = C (cid:37) (cid:20) F ( ξ − ξ , κ ) G ( ξ − ξ , κ ) (cid:21) + C (cid:37) (cid:20) G ∗ ( ξ − ξ , κ ) F ∗ ( ξ − ξ , κ ) (cid:21) . (A12)Now, for the case of negative slope we can find the solutions by noticing that in Eqs. (A3,A4) the negative-slopesolutions can be obtained from the positive-sign solution by complex conjugation. We then have (cid:20) ψ (cid:38) + ( ξ ) ψ (cid:38)− ( ξ ) (cid:21) = C (cid:38) (cid:20) F ∗ ( ξ − ξ , κ ) G ∗ ( ξ − ξ , κ ) (cid:21) + C (cid:38) (cid:20) G ( ξ − ξ , κ ) F ( ξ − ξ , κ ) (cid:21) . (A13)7 Appendix B: Insensitivity to details of the contact The question to be address in the following concerns the interface between the metal and the graphene sheet. Wewould like to investigate how the tunneling probability depends on the microscopic details of the tunneling betweenthe metal and the graphene sheet.To do so, let us consider that the electrons arriving from the metal side at the metal-graphene contact have only afinite transmission probability T c of crossing this interface into graphene. If this process is incoherent, we have, basedon general considerations (see e.g. Appendix B of Ref. [29] for a pedagogical derivation) that the total transmissionprobability Eq. (50) over the whole structure, becomes T tot ( v, κ ) = 12 | t tot ( v, κ ) | − − , (B1)where t L ( v, κ ) is the one-side transmission probability as seen from the metallic leads,1 | t tot ( v, κ ) | = 1 | t ( v, κ ) | + 1 T c ( v, κ ) − . (B2)We have now encapsulated the details of the interface into an unknown channel transmission probability T c ( v, κ ).In the absence of a detailed microscopic model, by running the simulations with various constant values of T c ( v, κ ) = T c , we can verify that the symmetry properties of the transmission probability are unaltered; the only effect of T c isto reduce the maximum value of σ and increase the maximum value of F in Fig. 9, but otherwise the shape of thesetraces are very similar.In principle, one can now make one more step forward in this analysis and start to include assumptions about theinterface, which may result in an energy-dependent T c . For example, we can consider a simple model for the metal asa graphene sheet with a very low Dirac point, and assume that the electrons are injected into the graphene at somecontact point x c < − d . The wavefunction in the metallic leads then reads1 t c (cid:20) (cid:21) e iv c ξ + r c t c (cid:20) − (cid:21) e − iv c ξ , (B3)where v c is the reduced Fermi vector of the metal, assuming a relatively large value v c (cid:29) κ . By matching Eq. (B3)at x c = ξ c / √ a with the transmitted wave (cid:20) ψ + ( p + v, κ ) ψ − ( p + v, κ ) (cid:21) e i sgn[ p + v ] √ ( p + v ) − κ ξ (B4)at the right side of the contact (corresponding to the graphene sheet on the metal) we obtain1 | t c ( v, κ ) | = 1 T c ( v, κ ) = | ψ + ( p + v, κ ) | = | p + v | (cid:112) ( p + v ) − κ + 12 . (B5)One can immediately see that T c ≤ 1, and also that T c ≈ κ (cid:28) | p + v | . Note also that this contact modeldoes not introduce any additional phenomenological parameters; in particular, owing to the assumption of incoherenttunneling, it does not depend on ξ c . By using Eq. (B5) and Eqs. (B1,B2) we can calculate the conductivity and Fanofactor. We have verified numerically that, even under this energy-dependent model, the traces for conductance andshot noise do not change significantly for samples with parameters as used in this work.Thus, we conclude that the main features of our model of tunneling under the trapezoidal barrier remain remarkablyrobust when the microscopic details of the interface between the metal and graphene are included. [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science ,666 (2004).[2] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevich, S.M. Morozov, A.K. Geim, PNAS , 10451 (2005).[3] A.K. Geim and K. S. Novoselov, Nat. Materials , 183 (2007)[4] P. E. Allain and J. N. Fuchs, Eur. Phys. J. B , 301 (2011).[5] Andrea F. Yound and Philip Kim, Annu. Rev. Condens. Matter Phys. , 101 (2011).[6] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. , 620 (2006). [7] M. I. Katsnelson, Eur. Phys. J. B , 157 (2006).[8] M. I. Katsnelson and K. S. Novoselov, Solid State Commun. , 3 (2007).[9] Markus M¨uller, Matthias Br¨auninger, and Bj¨orn Trauzzettel, Phys. Rev. Lett. , 196801 (2009).[10] V. V. Cheianov and V. I. Fal’ko, Phys. Rev. B , 041403 (R) (2006).[11] J. Cayssol, B. Huard, D. Goldhaber-Gordon, Phys. Rev. B, 075428 (2009).[12] J. Tworzyd(cid:32)lo, B. Trauzzetel, M. Totov, A. Rycerz, and C. W. J. Beenakker, Phys. Rev. Lett. , 246802 (2006).[13] C. W. J. Beenakker, Rev. Mod. Phys. , 1337 (2008).[14] E. B. Sonin, Phys. Rev. B , 233408 (2008).[15] E. B. Sonin, Phys. Rev. B , 195438 (2009).[16] F. Sauter, Z. Phys. , 742 (1931).[17] B. Huard, N. Stander, J. A. Suplizio, and D. Goldhaber-Gordon, Phys. Rev. B , 12402 (2008).[18] E. J. H. Lee, K. Balasubramanian, R. T. Weitz, M. Burghard, and K. Kern, Nat. Nanotechnol. , 486 (2008).[19] F. Xia, V. Perebeinos, Y-m. Lin, Y. Wu, and P. Avouris, Nat. Nanotechnol. , 179 (2011).[20] G. Giovannetti, P. A. Khomyakov, G. Brocks, V. M. Karpan, J. van den Brink, and P.J. Kelly, Phys. Rev. Lett. ,026803 (2008).[21] B. Huard, J. A. Sulpizio, N. Stander, K. Todd, B. Yang, and D. Goldhaber-Gordon, Phys. Rev. Lett. , 236803 (2007).[22] R. V. Gorbachev, A. S. Mayorov, A. K. Savchenko, D. W. Horsell, and F. Guinea, Nano Lett. , 1995 (2008).[23] N. Stander, B. Huard, and D. Goldhaber-Gordon, Rev. Lett. , 025807 (2009).[24] A. F. Young and P. Kim, Nat. Phys. , 222 (2009).[25] L. DiCarlo, J. R. Williams, Yiming Zhang, D. T. McClure, and C. M. Marcus, Phys. Rev. Lett. , 156801 (2008).[26] R. Danneau, F. Wu, M. F. Craciun, S. Russo, M.Y. Tomi, J. Salmilehto, A. F. Morpurgo, and P. J. Hakonen, Phys. Rev.Lett. , 196802 (2008).[27] R. Danneau, F. Wu, M. F. Craciun, S. Russo, M. Y. Tomi, J. Salmilehto, A. F. Morpurgo and P. J. Hakonen, J. LowTemp. Phys. , 374 (2008).[28] R. Danneau, F. Wua, M. F. Craciun, S. Russo, M.Y. Tomi, J. Salmilehto, A.F. Morpurgo, P. J. Hakonen, Solid StateCommun. , 1050 (2009).[29] A. Laitinen, G. S. Paraoanu, Mika Oksanen, Monica F. Craciun, Saverio Russo, Edouard Sonin, and Pertti Hakonen, Phys.Rev. B. , 115413 (2016).[30] J. Wintterlin, and M.-L. Bocquet, Surf. Sci. , 1841–1852 (2009).[31] Teresa Cusati, Gianluca Fiori, Amit Gahoi, Vikram Passi, Max C. Lemme, Alessandro Fortunelli, and Giuseppe Iannaccone,Sci. Rep. , 5109 (2017).[32] Xu Du, Ivan Skachko, Anthony Barker, and Eva Y. Andrei, Nat. Nanotechnol. , 491 (2008).[33] Ying Liu, T. S. Abhilash, Antti Laitinen, Zhenbing Tan, Guan-jun Liu, and Pertti Hakonen, Nanotechnology , 25LT01(2019).[34] A. Das, S. Pisana, B. Chakraborty, S. Piscanec, S. K. Saha, U.V. Waghmare, K.S. Novoselov, H. R. Krishnamurthy, A. K.Geim, A.C. Ferrari, and A. K. Sood, Nat. Nanotechnol. , 210 (2008).[35] H. Schomerus, Phys. Rev. B , 045433 (2007).[36] Ya. M. Blanter and I. Martin, Phys. Rev. B , 155433 (2007).[37] P. A. Khomyakov, G. Giovannetti, P. C. Rusu, G. Brocks, J. van den Brink, and P. J. Kelly, Phys. Rev. B , 195425(2009).[38] Jiaxin Zheng, Yangyang Wang, Lu Wang, Ruge Quhe, Zeyuan Ni, Wai-Ning Mei, Zhengxiang Gao, Dapeng Yu, JunjieShi, and Jing Lu, Sci. Rep. , 2081 (2013).[39] Yuki Matsuda, Wei-Qiao Deng, and William A. Goddard III, J. Phys. Chem. C2010 , 17845–17850 (2010).[40] Bo Ma, Cheng Gong, Yanwei Wen, Rong Chen, Kyeongjae Cho, and Bin Shan, Journal of Applied Physics , 183708(2014).[41] Daniel Huertas-Hernando, F. Guinea, and Arne Brataas, Phys. Rev. B74