7 Density of states of Dirac-Landau levels in a gapped 1 0 graphene monolayer under strain gradient 2 n V O Shubnyi1 and S G Sharapov2 a J 1 DepartmentofPhysics,TarasShevchenko NationalKievUniversity,6 4 AcademicianGlushkovave.,Kiev03680, Ukraine 2 2 BogolyubovInstitute forTheoreticalPhysics,NationalAcademyofScienceof Ukraine,14-bMetrolohichnaStreet,Kiev03680,Ukraine ] E-mail: [email protected] l l a h - Abstract. Westudyagappedgraphenemonolayerinacombinationofuniform s magnetic field and strain-induced uniform pseudomagnetic field. The presence e of two fields completely removes the valley degeneracy. The resulting density of m statesshowsacomplicatedbehaviourthatcanbetunedbyadjustingthestrength . of the fields. We analyze how these features can be observed in the sublattice, t a valleyand fulldensity of states. Theanalytical expression forthe valleyDOSis m derived. - d n o c Submitted to: J. Phys.: Condens. Matter [ 1 1. Introduction v 9 The carbon atoms in monolayer graphene form a honeycomb lattice due to sp2 6 7 hybridisationoftheirorbitals. SincethehoneycomblatticeisnotaBravaislattice,one 6 has to consider the honeycomb lattice as a triangular Bravais lattice with two atoms 0 per until cell. Thus one naturally arrivesat a two-componentspinor wavefunction of 1. the quasiparticle excitations in graphene (see Ref. [1] for some analogies with a real 0 spin). These components reflect the amplitude of the electron wave function on the 7 A and B sublattices. The two-component form of the wave function along with the 1 band structure results in the Dirac form of the effective theory for graphene. : v The Dirac fermions had shown up the celebrated magneto-transport and STS i properties of graphene (see Refs. [2, 3, 4] for the reviews). Recently STM/STS X measurementsallowednotonlytoobserverelativisticLandaulevels,butalsotoresolve r a directly their sublattice specific features. By resolving density of states (DOS) on A and B sublattices of a gapped graphene, it was experimentally confirmed [5] that the amplitude of the wave function of the lowest Landau level (LLL) is unequally distributed between the sublattices depending on its energy sign. In the presence of a gap ∆ driven by inversion symmetry breaking, the LLL splits into two levels with the energy E = η∆sgn(eB), where η = ± distinguishes 0 inequivalent K and K′ points of the Brillouin zone and an external magnetic field B = ∇×A = (0,0,B) is applied perpendicularly to the plane of graphene along Density of states of Dirac-Landau levels 2 the positive z axis [6]. Here e = −|e| is the electron charge and A is the vector electromagnetic potential. The corresponding amplitudes of the wave function of the positive energy electron-like, 0 , and negative energy hole-like, 0 , levels are on A + − and B sublattices. In other words, the individual sublattices are valley polarized for the LLL [7]. An exciting opportunity for manipulating the amplitudes of the wave function on the sublattices opens due to close connection between the impact of deformation and external electromagnetic field on the electron structure of graphene. Change in hopping energybetweenAandB atomsinduced bystraincanbe describedby vector potential A analogousto the vector potential A of the externalmagnetic field (see pm Refs. [8, 9] for a review). The correspondingfield, B =∇×A , is calledpseudomagnetic field(PMF), pm pm because it is directed oppositely in K and K′ valleys. This implies that the LLL breaks the electron-hole symmetry, with the LLL energy, E = ∆sgn(B ) for both 0 pm K and K′ points. Furthermore, the states corresponding to the LLL are sublattice polarized, because they reside exclusively on either A or B sublattice [7]. WhiletheformationoftheLLLisassociatedwithzeromodesanddoesnotrequire ahomogeneousPMF,toformhigherLandaulevelsauniformPMFisneeded[10]. This isinfactthemainchallenge[11]fortheimplementationofstrainedgraphene,although recently there has been some progress both in experiment [12, 13, 14] and in theory [15]. Inthepresenceofeitherexternalmagneticfieldordeformation,thehigherenergy levels from K and K′ points remain degenerate. This degeneracy is lifted when both strainandmagneticfieldarepresent. Oneoftheinterestingconsequencesofthelifting is that for |B | > |eB|, the Hall conductivity is oscillating between 0 and ∓2e2/h pm [16]. The latest experiments [14] show that it is possible to create a homogeneous PMF of order of a few Tesla. Therefore, there is a good chance that the STS/STM measurements of the Dirac Landau levels similar to that done in Ref. [5] are now possible on strained graphene. Thus the purpose of the present work is to study the DOS(including the sublattice resolved)in a combinationofa constantPMFB pm createdbynon-uniformstrainandmagneticfieldB. Inparticular,wewilllookforthe specific effects related to the presence of nonzero gap ∆ and lifting of the degeneracy between K and K′ that can be observed in STS measurements. The paperis organizedasfollows. We beginby presentingin section2 the model describinggappedmonolayergrapheneinthecombinationofPMFandmagneticfield. In section 3 we provide the definitions of the valley, sublattice and full DOS in terms of the Green’s function decomposed over Landau levels. The corresponding DOS are written in section 4 as the sums that in the case of the valley DOS can be calculated analytically. The resultsfor the DOSin the variousregimes arediscussedinsection5 and conclusions are given in section 6. 2. Model We consider gapped monolayer graphene in the continuum approximation described by the effective Hamiltonian HK 0 H = , (1) 0 HK′ (cid:18) (cid:19) Density of states of Dirac-Landau levels 3 The full Hamiltonian (1) acts on wave function with four components ψ• K ψ◦ Ψ= K , (2) ψ◦ K′  ψK•′    where•and◦denote,respectively,AandB sublatticesandwefollowedthe notations of Refs. [6, 3] with exchangingthe sublattices in the K′ valley. Thus the Hamiltonian (1) includes two blocks corresponding to K and K′ valleys e HK =vFτττ −i~∇− A−Apm +τ3∆, (3) c (cid:16) e (cid:17) HK′ = −vFτττ −i~∇− A+Apm −τ3∆. (4) c (cid:16) (cid:17) Here τττ = (τ ,τ ) and τ are Pauli matrices acting in the sublattice space, v is the 1 2 3 F Fermi velocity, the gap ∆ corresponds to the energy difference 2∆ between the A and B sublattices, A and A are the the electromagnetic and strain induced vector pm potentials, respectively. We neglect the spin splitting, because for commonly used strengths of magnetic field the Zeeman splitting is small compared to the distance between the Landau levels. For a fixed direction of external magnetic field, the correspondingtoAtermintheHamiltonianbreakstime-reversalsymmetry,whilethe A termbreakstheinversionsymmetryandleavestime-reversalsymmetryunbroken. pm Withthex-axisalignedinthezigzagdirection,thestrain-inducedvectorpotential reads [17, 18] (see also the reviews [8, 9]) ~βκ u −u A = xx yy , (5) pm 2a0 (cid:18) −2uxy (cid:19) where β = −∂lnt/∂lna| ≈ 3 is the dimensionless electron Gru¨neisen parameter a=a0 for the lattice deformation, t the nearest-neighbour hopping parameter, κ ≈ 1/3 is a parameter related to graphene’s elastic property [17], a is the length of the carbon- carbon bond, (a is the length of the unstrained bond), and u with i,j =x,y is the 0 ij strain tensor as defined in classical continuum mechanics [8, 9]. We also assume that thedeformationisapureshear,sothatu +u =0,andthereisnoscalarpotential xx yy term in the Hamiltonian. The sign of the PMF depends on the valley, and, for example, in K valley, ~βκ 1 B =− ∂ u + ∂ (u −u ) , (6) pm x xy y xx yy a 2 0 (cid:18) (cid:19) whereas it has the opposite sign in K′ valley, because A enters Eqs. (3) and (4) pm with the opposite signs. Eq. (6) illustrates the main problem in this field of research, viz. a uniform PMF can only be created by a non-uniform strain [11]. As was already stated in the Introduction, considering the experimental progress achieved in the field [14], we restrict ourselves to a constant PMF. Thus we arrive at the model with two independent K points characterized by the following combinations of the fields, B = eB/c±B . A more complicated, but analytically intractable ± pm case with a combination of a constant magnetic and inhomogeneous pseudomagnetic fields was considered in Ref. [19], where a circularly symmetric strain is induced by a homogeneous load. Density of states of Dirac-Landau levels 4 3. Green’s function, sublattice and valley resolved DOS Although it is straightforward to obtain the DOS directly from the solution of the corresponding Dirac equation, we rely on the Green’s function (GF) machinery that automatically takes into account the degeneracy of levels and avoids the necessity to work with different directions of fields separately. Since the K points in the model (1) are independent, we will use the GF’s corresponding to the separate K points. In particular, we are interested in the translation invariant part G of the GF that allows to derive both the DOS and the transportcoefficients. Its derivationusing the Schwinger proper-time method and decomposition over Landau-leveel poles has been discussed in many papers (see, e.g., Refs. [20, 21, 22, 23]). Here we begin with the translationinvariantpartforKpointwrittenintheMatsubararepresentation(weset ~=c=k =1 in what follows) B GK(iω,p)=e−|Bp+2| ∞ (−1)n Gn(B+,iω,p), ω =π(2m+1)T, (7) (iω)2−(M+)2 n=0 n X wehere T is the temperature, M± = ∆2+2nv2|B | (8) n F ± are the energies of theqrelativistic Landau levels at K and K′ points (η = ±), respectively, and the function 2p2 G (B ,iω,~p)=(∆τ +iω) (1+τ sgn(B ))L n + 3 3 + n |B | (cid:20) (cid:18) + (cid:19) 2p2 − (1−τ sgn(B ))L (9) 3 + n−1 |B | (cid:18) + (cid:19)(cid:21) 2p2 −4v (p τ +p τ )L1 . F x y y x n−1 |B | (cid:18) + (cid:19) Here Lα(z) are the generalized Laguerre polynomials, and L (z)≡L0(z) (L1 ≡0). n n n −1 WhenderivingGFfromtheknownwave-functions,theLaguerrepolynomialsoriginate fromthe integrationoftwo Hermite polynomials withproper weights. Lookingatthe structure of the GF (7), one can see that the projectors P = (1± τ sgn(B ))/2 ± 3 + take into account that, for example, for B > 0, the states on A and B sublattices + involve L and L , respectively. The most general expression of the propagator in n n−1 the presence of B, B and various types of the gaps is provided in [24]. pm ThecorrespondingcontributionoftheKpointtotheDOSperspinandunitarea on A and B sublattices reads 1 d2p K K K D (ǫ)= [G (ǫ−i0,p)−G (ǫ+i0,p)] (10) A,B 2πi (2π)2 ii ii Z with i=1,2 for A and B sublattices. It feollows from Eq.e(1) that GK′(iω,p)=GK(v →−v ,∆→−∆,B →B ) (11) F F + − and e e 1 d2p DK′ (ǫ)= [GK′(ǫ−i0,p)−GK′(ǫ+i0,p)] (12) B,A 2πi (2π)2 ii ii Z with i = 1,2 for B and A sublattices, ie.e. exchanging tehe sublattices. While the valley resolved DOS presents a theoretical interest and will also be considered below, Density of states of Dirac-Landau levels 5 the STS measurements allow to observe the full DOS involving two valleys on each sublattice K K′ D (ǫ)=D (ǫ)+D (ǫ). (13) A,B A,B A,B We will also consider the valley resolved but summed over sublattices DOS DK,K′(ǫ)=DK,K′(ǫ)+DK,K′(ǫ) (14) A B which presents interest for valleytronics. Finally, the full DOS can also be found by summing the valley resolved DOS K K′ D(ǫ)=D (ǫ)+D (ǫ). (15) 4. Expressions for numerical and analytical calculation of the DOS Using the integral [25] d2p 2p2 p2 (−1)n|B| L exp − = (16) (2π)2 n |B| |B| 2π 2 Z (cid:18) (cid:19) (cid:18) (cid:19) and evaluating the discontinuity of the GF we arrive at the final result DK,K′(ǫ)=DK,K′(0)(ǫ)+DK,K′(n≥1)(ǫ). (17) A,B A,B A,B HereDK,K′(0)(ǫ)istheLLLcontributiontothevalleyandsublatticeresolvedDOSand A,B DK,K′(n≥1)(ǫ) is the corresponding contribution from the Landau levels with n ≥ 1. A,B Explicit expressions for these terms are DK(0)(ǫ)= |B+|θ(B )δ(ǫ−∆), DK′(0)(ǫ)= |B−|θ(−B )δ(ǫ−∆), (18) A 2π + A 2π − DK(0)(ǫ)= |B+|θ(−B )δ(ǫ+∆), DK′(0)(ǫ)= |B−|θ(B )δ(ǫ+∆), (19) B 2π + B 2π − and for the Landau levels with n≥1, DK,K′(n≥1)(ǫ)= ∞ |B±| Mn±+∆δ(ǫ−M±)+ Mn±−∆δ(ǫ+M±) , (20) A 2π 2M± n 2M± n n=1 (cid:20) n n (cid:21) X DK,K′(n≥1)(ǫ)= ∞ |B±| Mn±−∆δ(ǫ−M±)+ Mn±+∆δ(ǫ+M±) , (21) B 2π 2M± n 2M± n n=1 (cid:20) n n (cid:21) X where±signcorrespondstoKandK′ points. Asexpected,presenceofPMFremoves degeneracy of the levels with n≥1 [16] . In section 5 we compute the sublattice and valley resolved DOS numerically on the base of Eqs.(18), (19), (20) and (21) by widening δ-fuction peaks to a Lorentzian shape, viz. 1 1 δ(ǫ−M )→ , (22) n π(ǫ−M )2+Γ2 n n where Γ is the n-thlevelwidth. SuchbroadeningofLandaulevelswith aconstantΓ n wasfoundto be ratheragoodapproximationvalidinnotverystrongmagnetic fields. Density of states of Dirac-Landau levels 6 4.1. The DOS in the zero pseudomagnetic field Setting B =0 werecoverthe well-knownresults thatwereexperimentallyobserved pm in [5]. Then Eqs. (18) and (19) result in the sublattice DOS |eB| |eB| D(0)(ǫ)= δ(ǫ−∆), D(0)(ǫ)= δ(ǫ+∆). (23) A 2π B 2π This confirms that the LLL is valley polarized, because each LLL contribution to the DOS comesfromeither K orK′ valley,as discussedinthe Introduction. This feature has to be contrasted with the valley resolved but summed over the two sublattices DOS DK,K′(0)(ǫ)=DK,K′(0)(ǫ)+DK,K′(0)(ǫ)= |eB|δ(ǫ−η∆sgn(eB)). (24) A B 2π Forn≥1thelevelsatKandK′pointsdescribedbyEqs.(20)and(21)aredegenerate, but the DOS on A and B sublattices differs and this effect is observable [5]. 4.2. The DOS in the zero magnetic field Setting eB = 0 we obtain from Eqs. (18) and (19) that the LLL contribution to the sublattice DOS is |B | |B | D(0)(ǫ)= pm θ(B )δ(ǫ−∆), D(0)(ǫ)= pm θ(−B )δ(ǫ+∆) (25) A 2π pm B 2π pm This confirms that the LLL is sublattice polarized, as discussed in the Introduction. 4.3. Analytical expression for the valley DOS Although the expressions for the sublattice and valley DOS presented in Sec. 3 are sufficient for the numerical study presented in Sec. 5, it is always useful to have a simple analytical expression for the DOS. One can notice that the the valley DOS, Eq. (14) is the sum of delta-functions (or Lorentzians when the the level widening is takeninto account),because the sum ofthe weightfactors (M±±∆)/(2M±) present n n in Eqs. (20) and (21) gives 1. This allows one to use the results of Ref. [26], and calculate the sum over Landau levels analytically 1 Γ Γ DK,K′(ǫ)= −|B |θ(∓B ) −|B |θ(±B ) 2π2 ± ± (ǫ−∆)2+Γ2 ± ± (ǫ+∆)2+Γ2 ( Λ2 ∆2−(ǫ+iΓ)2 +Γln −Im (ǫ+iΓ)ψ . (26) 2|B | 2|B | ± (cid:20) (cid:18) ± (cid:19)(cid:21)) Here ψ is the digamma function, ± sign corresponds to K and K′ points, the width ofalllevelsΓ isassumedtobe the same,andΛisthe cutoffenergythathasthe order of bandwidth. Eq. (26) differs from Eq. (4.15) of Ref. [26] by the first two terms. In the present case they take care of the electron hole asymmetry of the LLL, while in [26] both K and K′ points contribute to the full DOS. The advantage of Eq. (26) is that it allows to consider the low field regime when the direct numerical summation over many Landau levels is consuming. Density of states of Dirac-Landau levels 7 5. Results NowweuseEqs.(18), (19),(20)and(21)tostudythe valley(14),sublattice(13)and thefull(15)DOSnumerically. ForsimplicityweassumethatallLandaulevelshavethe samewidthΓ. InthiscasethevalleyDOSandthenthefullDOScanalsobecalculated using Eq. (26). To fit real experimental data [27] it may be necessary to consider the width, Γ , dependent on the Landau level index. This can be easily done in the n frameworkofnumericalcomputationofthesumoverLandaulevels. However,whenall levelshavethe samewidthandone is interestedinthe valleyDOS,it is moreefficient tocomputeitfromtheanalyticalexpression(26)whichiseasiertouseinthelowfield regime. In all numerical work we take the value of the Fermi velocity, v = 106m/s F that corresponds to the Landau energy scale, ǫ =(~v2B )1/2 = 25.7 B [T] meV. 0 F ± ± The gap∆thatlifts theenergydegeneracyofthe AandB sublattices andbreaksthe p inversion symmetry was observed for a graphene monolayer on top of SiC, graphite [4], and hexagonalboron nitride [28]. Its value ranges from 10 meV to severaltens of meV. Fig. 1 demonstrates how the full density of states, Eq. (15), is formed by the contributionsfromthevalleyresolvedDOS,Eq.(14): leftpanel(a)isfor|eB|<|B | pm and the right panel (b) is for |eB| > |B |. The two curves (thin solid red and pm HaL HbL ÈeBÈ<ÈB È ÈeBÈ>ÈB È pm pm K’ K K K¢ -3 -2 -1 0 1 2 3-3 -2 -1 0 1 2 3 Ε(cid:144)Ε0 Ε(cid:144)Ε0 Figure 1. (Colour online) The full DOS, D(ǫ), (thick solid) and the valley- resolved DOS, DK,K′(ǫ), (thin solid and thin dashed) in arbitrary units as the functions of energy ǫ/ǫ0, where ǫ0=(~vF2B+)1/2 = 25.7pB+[T] meV is the Landau scaleforKvalley. Leftpanel: (a): thefields B=5TandBpm =18T. Rightpanel: (b): thefieldsB=10TandBpm=1T. Thegap∆=50meVand thescatteringrateΓ=10meVinthebothcases. thin dashed blue) in the bottom part of the figure show the valley resolved DOS, DK,K′(ǫ). The curves for K and K′ points have the peaks corresponding to the relativistic Landau levels with the energies ∼ ± n|B |. The position of the peaks ± corresponding the LLL with E =±∆ distinguish the cases (a) the PMF dominated 0 p regime, |eB| < |B |, when both peaks have the same sign of the energy, and (b) pm Density of states of Dirac-Landau levels 8 the magneticfielddominatedregime,|eB|>|B |,whenthe peakshavetheopposite pm energysign. Wecheckedthatthesamecurvesalsofollowfromtheanalyticalexpression (26). Thosearerather trivialconsequencesofhavinga superpositionofmagnetic and PMF. The full DOS D(ǫ) shown by thick black curve in the two panels obviously has twoseriesofpeaks. Onecouldseethatinthespecialcases,thedifferencebetweentwo curves is substantial, and resulting DOS curve has irregular features and/or masked peaks. Depending on the values of the effective fields B , the Landau levels could ± be viewed as a splitting of one level (in case |B | ≈ |B |) or as the two largely + − independent series, as for the case shown in Fig. 1. Let us look closer at the pattern that overlapping Landau levels may create for certain values of B and B . The energies of the Landau levels with indices n , n pm + − for K and K′ points coincide, viz. M = M if there exist some values of B n+ n− + and B satisfying the condition, |B |n = |B |n . This implies that the fraction − + + − − |B |/|B | = a/b has to be rational. In terms of the initial fields B and B this + − pm condition implies that 1−a/b eB = B . (27) pm 1+a/b The corresponding beating patterns for four values of the fraction a/b are shown in Fig. 2. The lowest (green) curve is for the simplest case, B /B =1/2. Eachsecond − + B /B =1/4 - + s e t a St 2/3 f o y t si 1/3 n e D 1/2 0 1 2 3 4 Energy(scaledindependently) Figure2. (Colouronline)ThefullDOS,D(ǫ),inarbitraryunitsasthefunctions of energy ǫ/ǫ0, where ǫ0=(~vF2B+)1/2 = 25.7pB+[T] meV. The green curve is for B−/B+ = 1/2, the blue curve is for B−/B+ = 1/3, the red curve is for B−/B+ =2/3 and the black curve is for B−/B+ =1/4. The gap ∆=10 meV andthescatteringrateΓ=2meV. level with coinciding energies is enhanced. The second from the bottom (blue) curve is for B /B = 1/3. In this case an enhancement occurs for each third level. The − + third from the bottom (red) curve is for B /B =2/3 has even more tricky pattern − + Density of states of Dirac-Landau levels 9 with the highest each third level and each forth level of an intermediate height. The curve on the top (black) is for B /B =1/4. − + It is instructive to represent the dependences of the full DOS, D(ǫ), on the fields B and B employing the density plot. Since a wide range of the fields is involved, pm its consideration in the low field regime may demand summation over many Landau levels. Thus we use Eq. (26), where the summation is done analytically. Figs. 3 (a) and(b) onthe toppanelshowthe full DOS, D(ǫ,B,B )as a functionofenergy ǫin pm meV and magnetic field B in T for B =0 T and B =8 T. Figs.3 (c) and (d) in pm pm the bottom panel show the full DOS, D(ǫ,B,B ) as a function of energy ǫ in meV pm and PMF B in T for B = 0 T and B = 8 T. The density plot is partly overlaid pm with the solid (red) and dashed (blue) curves that show position of the peaks in the DOS originating from the Landau levels at K and K′ points, respectively. Fig. 3 (a) (top left panel) describes unstrained graphene. The Landau levels fan away from the Dirac point at ǫ = 0. One can find a similar DOS map for the STS measurements [4] of graphene on chlorinated SiO . In the real case the spectra are distorted at low 2 fields due to the substrateinduced disorderandarestronglypositiondependent. The density plot [4] allows to observe at higher fields the sequence of broadened Landau levelswithseparatedpeaks. Fig.3(c)(bottomleftpanel)describesstrainedgraphene inzeromagneticfield. ItisalmostidenticaltoFig.3(a)excepttotheLLLthatinthe caseofstrainedgraphenebreakstheelectron-holesymmetry. Fig.3(b)and(d)(right top and bottom panels) describe strained graphene in the external magnetic field. This casewasalsostudiedexperimentallyin[13], whereSMT andSTSmeasurements were made on the deformed by gating graphene drumhead. Comparing all these panels we observe that in the presence of both PMF and magnetic field there exist regions of intersecting Landau levels with the opposite slope that are related to the opposite valleys. In Fig. 3 (b) this is the region with |eB|<|B |, while Fig.3 (d)the correspondingregionisseenfor |B |<|eB|. This pm pm behaviour of Landau levels is almost obvious in the presented case. However, in the caseofpoorlyresolvedLandaulevelsthisfeaturecanberatherhelpfulforprovingthe presence of both PMF and magnetic field. Finally, we illustrate in Fig. 4 how the full DOS is distributed between the sublattices. The DOS on A and B sublattices, D (ǫ), are shown by (dashed) A,B green and solid (orange) curves, respectively. Fig. 4 (a) (top left panel) describes unstrained graphene in the external magnetic field. It corresponds to the situation studied experimentally in [5]. We observe that the positive (negative) energy states resideonA(B)sublattice. Sincethesestatesareassociatedwithdifferentvalleys,the LLLisindeedvalleypolarized. Furthermore,thesublatticeasymmetryisalsoseenfor higher levels, because we took a large value of the gap ∆ = 50 meV. Fig. 4 (b) (top right panel) describes strained graphene in zero magnetic field. As it should be, the LLL is indeed completely sublattice polarized,while higher levels are polarizedin the same fashion as in Fig. 4 (a). The PMF dominated regime, |eB| < |B |, is shown pm in Fig. 4 (c) (bottom left panel). The LLL polarization is similar with Fig. 4 (b). The magnetic field dominated regime, |eB| > |B |, is shown in Fig. 4 (d) (bottom pm right panel) and it is similar to Fig. 4 (a). Fig. 4 (c) and (d) are computed for the same values of the parameters as Fig. 1 (a) and (b), respectively. When both fields are present the asymmetry between the sublattices can be enhanced even for higher levels. We note that in the present work the sublattice asymmetry is directly brought by the inversion symmetry gap ∆. We established that the presence of PMF and Density of states of Dirac-Landau levels 10 20 10 B 0 -10 HaL Bpm = 0T HbL Bpm = 8T -20 20 10 m p 0 B -10 HcL B = 0T HdL B = 8T -20 -200 -100 0 100 200 300 -200 -100 0 100 200 300 Ε, meV Ε, meV Figure 3. Density map of the full DOS D(ǫ,B,Bpm) as a function of energy ǫ inmeV,magneticfieldBinTandPMFinT. Topleftpanel: (a): foraconstant Bpm=0T. Toprightpanel: (b): foraconstant Bpm=8T. Bottom leftpanel: (c): for a constant B =0 T. Bottom right panel: (d): for a constant B = 8 T. The density-map is overlaid with red and blue curves that show the position of thepeaksoriginatingfromKandK′points. Thegap∆=50meV,thescattering rate Γ=5meV and the Landau scale ǫ0=(~vF2B±)1/2 =25.7pB±[T] meV in allcases. magneticfieldfurther enhancesthis effect. Itis shownin[29]thatthe localsublattice