Finite-temperature Coulomb Excitations in Extrinsic Dirac Structures
Andrii Iurov, Godfrey Gumbs, Danhong Huang, Ganesh Balakrishnan
FFinite-temperature Coulomb Excitations in Extrinsic Dirac Structures
Andrii Iurov ∗ Center for High Technology Materials, University of New Mexico,1313 Goddard SE, Albuquerque, NM, 87106, USA
Godfrey Gumbs
Department of Physics and Astronomy, Hunter College of the City University of New York,695 Park Avenue, New York, NY 10065, USA andDonostia International Physics Center (DIPC), P de Manuel Lardizabal,4, 20018 San Sebastian, Basque Country, Spain
Danhong Huang
Air Force Research Laboratory, Space Vehicles Directorate,Kirtland Air Force Base, NM 87117, USA andCenter for High Technology Materials, University of New Mexico,1313 Goddard SE, Albuquerque, NM, 87106, USA
Ganesh Balakrishnan
Center for High Technology Materials, University of New Mexico,1313 Goddard se, Albuquerque, NM, 87106, USA (Dated: July 28, 2017)We have derived algebraic, analytic expressions for the chemical potential without any restrictionon temperature for all types of doped, or extrinsic, gapped Dirac cone materials including gappedgraphene, silicene, germanene and single-layer transition metal dichalcogenides. As an importantintermediate step of our derivations, we have established a reliable piecewise-linear model for cal-culating the density-of-states in molybdenum disulfide, showing good agreement with previouslyobtained numerical results. For the spin- and valley-resolved band structure, we obtain an addi-tional decrease of the chemical potential due to thermally induced doping of the upper subband atfinite temperature. It has been demonstrated that since the symmetry between the electron and holestates in
MoS is broken, the chemical potential could cross the zero-energy level at sufficiently hightemperature. These results allow us to investigate the collective properties, polarizability, plasmonsand their damping. Emphasis is placed on low temperatures, when initial electron doping plays acrucial role. We clearly demonstrated the contribution of the initial doping to the finite-temperaturecollective properties of the considered materials. PACS numbers: 73.21.-b, 73.63.-b, 71.45.Gm, 73.20.Mf
I. INTRODUCTION
Despite the fact that the microscopic properties of various low-dimensional materials have been meticulously exam-ined over a fairly long period of time, only successful fabrication of graphene in 2004 stimulated an intriguinglynew research effort devoted to the study of atomically thin two-dimensional (2D) materials. In particular, it was byvirtue of its unique, yet unexpected massless Dirac electronic properties that led to high mobility (200,000 cm /V · s)and ballistic transport properties. In the corners of the Brillouin zone, referred to as K and K (cid:48) points, there is noenergy band gap and the dispersions represent a linear Dirac cone structure. Due to the existence of such an energyspectrum, opening a sufficiently large and tunable energy gap in graphene has become an important issue in orderto enable electron confinement. Researchers tried to achieve this by adjoining a variety of insulating substrates or even expos graphene to circularly-polarized radiation. In finite-width nanoribbons, their energy band structureand gap are modified by the type of insulating ‘cousin’ that is introduced.
In order to create a truly tunable band gap, one must use a material with large spin-orbit coupling or a buckledstructure. In this regard, silicene, a 2D silicon structure was deemed a good candidate. Single-monolayer Si possessesa buckled structure simply because of the larger ionic size of silicon compared to carbon. This results in a largespinorbit band gap 1 . meV and the possibility to modify its energy spectrum by applying an external perpendicularelectric field . These properties make it display an experimentally realizable KaneMele type of quantum spinHall effect, or a topological insulator state, because of the existence of time-reversal symmetry. Unlike graphene,the band structure of silicene and its nanoribbons directly depends on spin and valley indices which give lead toplenty of nanoelectronic, valleytronic and spintronic applications. a r X i v : . [ c ond - m a t . m e s - h a ll ] J u l Germanene, the most recently discovered and fabricated member of atomically thin buckled 2D honeycomb lat-tices, demonstrates substantially larger Fermi velocities and a band gap of 20 − meV . Grown by molecularbeam epitaxy and investigated with x-ray absorption spectroscopy, Ge layers demonstrated satisfactory agreementbetween the experimentally obtained and theoretically predicted results for its inter-atomic distance.Another important class of innovative 2D materials is represented by direct-gap transition metal dichalchogenides,or TMDC’s. Its chemical makeup consists of a transition metal atom M , such as molybdenum or tungsten, andtwo identical chalcogens C , i.e., sulfur, selenium or tellurium. Schematically, TMDC’s are described as MC . In ourconsideration, we are mostly going to focus on M oS , as their most studied representative. This material exhibits asemiconductor energy band structure with a very large direct gap 1 . eV , in contrast to its bulk states with indirectgap 1 . eV , and substantial spin-orbit coupling. Strictly speaking,
M oS is not a Dirac material since the massterms play a crucial role in its energy dispersions, however its low-energy Hamiltonian contains a t a Σ · k term,corresponding to the linear Dirac cone dispersion.An effective two-band continuum model and lattice Hamiltonian, based on the tight-binding model, accounts forthe hybridization of the d orbitals in M o and the p orbitals of sulfur atoms. It gives an adequate description of itslow-energy band structure and predicts large spin splitting. Due to the breaking of inversion symmetry and spin-orbit coupling, spin and valley physics is observed in all group-IV dichalcogenides, including
M oS . The low-energystates of such systems are no longer massive Dirac fermions since there is a a difference between electron and holemasses as well as trigonal warping effects. Strain engineering, used to tune optical and electronic properties ofconventional semiconductors, has also been applied to molybdenum disulfide, and its modified band structure hasbeen theoretically calculated. These unique electronic properties of a single-layer
M oS were later used to createhigh-performance transistors operating at room temperature. These electronic models and effective Hamiltonianhave been widely used to investigate the collective properties of TMDC’s and their influence on the gap transition. In optoelectronics, the band structure, spin and valley properties of molybdenum disulfide could be successfullycontrolled by off-resonant dressing field. Current many-body and quantum field theory methods in condensed matter physics have provided helpful waysto understand the electronic and transport properties of low-dimensional solids, including diverse bucked honeycombmaterials.
In most of these theories, we find the dynamical polarization function , or polarizability, to be themainstay, fundamental quantity, describing the screening of an external potential by interacting electrons.
Also, the dynamic polarization function plays a key role in calculating the plasmon excitations, due to the chargedensity oscillations, which occur in metals and semiconductors. Specifically, the plasmon dispersion relation alongwith their lifetimes have been theoretically investigated for a wide range of 2D Dirac systems.
The interestin graphene plasmons is due in part to the fact that these excitations have no classical counterpart. Therehas also been a considerable experimental effort for investigating graphene plasmons, gate-tuning, infrared nano-imaging and confinement
Graphene plasmonic resonances and instability at various wavelengths could be usedin photodetectors in the Terahertz range. All these techniques could be successfully applied to the recentlyfabricated materials, discussed in the present work.Plasmonic applications have been widely based on nanoscale hybrid systems , in which graphene plasmons arecoupled to a surface plasmon excitations in metals. Technology has now gone a long way in combining graphene withprefabricated plasmonic nanoarrays and metamaterials in order to produce plasmonics-based tunable hybrid opticaldevices. Therefore, accurate knowledge of plasmon mode dispersions in graphene interfaced with metallic substratesis crucial. Graphene-metal contacts are important components for all such devices. Consequently, exploration ofplasmon modes at these metallic interfaces is a mandatory step toward fabricating the devices. High-resolutionelectron energy loss spectroscopy (EELS) has been employed to investigate those excitations at the surface of Bi Se to disclose the interplay between surface and Dirac plasmons in topological insulators . Plasmons, their behavior,dispersions, quenching and environmental effects have been thoroughly studied in epitaxial graphene, in air-exposedgraphene- Ru contacts, Graphene on P t N i (111), graphene grown on Cu (111) foils. In all cases under consideration, we need to distinguish between extrinsic , or a sample initially doped at T = 0,and intrinsic materials, with zero Fermi energy and completely empty conduction band. In the latter case, both theplasmon excitations and the electrical conductivity are completely suppressed at zero temperature due to the absenceof free charge carriers. However, at finite temperature, the conduction band would receive thermally induced dopingin both cases. In graphene, with zero energy band gap, this density is enhanced as n (cid:119) T and the plasmondispersion behaves like Ω p (cid:119) qT .The properties of intrinsic finite-temperature plasmon excitations have been systematically examined for variousmaterials including silicene. In contrast, extrinsic or doped structure at finite temperature is associated with adifficulty to obtain a reliable and accurate value for temperature-dependent chemical potential µ ( T ). Generally, it isknown that µ ( T ) is decreased as the temperature is increased, and its value could be found based on carrier densityconservation In this work, our main objective is to obtain a set of non-integral, trancendental equations for awide class of Dirac gapped materials with linear density-of-states (DOS), i.e., gapped graphene, silicene, germaneneand transition metal dichalcogenides at arbitrary temperature.However, the range of our considered temperatures is limited by validity of linear or gapped Dirac cone approxima-tion for the energies, which recieve noticeable doping at those temperatures. Certain deviations start to build up atabout 0 . eV , leading to various effect on the plasmons, such as anisotropy, splitting and existence of additionalacoustic plasmon branch. Such energies are extremely far away from our range k B T (cid:118) E (0) F . In all our calcu-lations, the energy is measured in the units of a typical Fermi energy E (0) F = 5 . meV , corresponding to electrondensity n (0) = 1 . · m − .Once the chemical potential is known, one can obtain the finite-temperature dynamical polarization function byEq. (18), which is a key component for all relevant many-body calculations. These include optical absorption,electronic transport, plasmon excitations as well as electron exchange and correlation energies. Here, we pay closeattention to finite-temperature plasmons, demonstrating how much initial doping contributes to each branch locationat intermediate temperatures. Once the temperature becomes very high, k B T (cid:29) E F , thermally-induced dopingdominates and the contribution from the initial Fermi energy fades away.The rest of the paper is organized as follows. First, we derive the implicit analytic equations for the finite-temperature chemical potential for all types of Dirac structures with linear DOS in Sec. II. In Sec. III, we calculatethe dynamic polarization function which includes the single particle excitation mode frequencies. The single-particlemodes combine with a charge cloud to produce weakly interacting quasiparticles that vibrate collectively at thecharacteristic plasma frequencies. Emphasis has been placed on rather simple cases of gapped graphene and siliceneat small, but finite temperatures, when the zero temperature carrier doping plays a crucial role. We also brieflyexamine non-local, hybrid plasmons in an open system of a semi-infinite conductor, Coulomb-coupled to a 2D layerin the presence of finite doping doping of the electrons in the layer. Our concluding remarks and a concise discussionare presented in Sec. IV. We also provide two appendices with detailed derivations of the DOS for silicene and M oS - in Appendix A, and of the temperature-dependent chemical potential in Appendix B. II. CHEMICAL POTENTIAL
In this section, we discuss our analytical derivations for the electron chemical potential as a function of temperature.Being equivalent to the Fermi energy at T = 0, the chemical potential normally decreases with increasing temperature.Its specific value depends on multiple material parameters such as energy band gaps, Fermi velocities and the DOS ofthe of the electrons as well as the holes below the zero energy level. Thus, for a conventional 2D electron gas (2DEG)with no holes, the chemical potential could become negative at a certain temperature, which is not possible for a Diracsystem with symmetry between the electron and hole states. Here, we are going to provide closed-form trancendentalequations for the finite-temperature µ ( T ) for a number of Dirac systems: graphene, buckled honeycomb lattices andtransition metal dischacogenides. A. Buckled honeycomb lattices
One of the most outstanding features of silicene and other buckled honeycomb lattices is the existence of twogenerally double degenerate pairs of energy subbands and two inequivalent band gaps. These gaps There is a fixedintrinsic spin-orbit gap 2∆ SO and a tunable sublattice-asymmetry gap ∆ z which is induced by and proportional toan applied perpendicular electric field E z . For small fields, ∆ z = E z d ⊥ , where d ⊥ is the out-of-plane displacement ofa buckled lattice.The low-energy model Hamiltonian of a buckled honeycomb lattice has been found to be ˆ H ξ,σ = (cid:126) v F ( ξk x ˆ τ x + k y ˆ τ y ) ⊗ ˆ I × − ξ ∆ SO ˆΣ z ⊗ ˆ τ z + ∆ z ˆ τ z ⊗ ˆ I × , (1)where the Fermi velocity v F = 0 . · m/s is half that for graphene, ξ = ± K/K (cid:48) valley index, τ x,y,z andΣ x,y,z are Pauli matrices in two different spaces, attributed to pseudospin and real spin of the considered structure.Introducing spin index σ = ±
1, we can rewrite Eq. (1) in a block-diagonal matrix form ( ) b ( ) a ( ) c ( ) -5.0 -2.5 5.0 d FIG. 1: (Color online) Energy dispersions and density-of-states (DOS) for silicene (upper panels ( a ) and ( b )) and molybdenumdisulfide (plots ( c ) and ( d )). Left panels ( a ) and ( c ) represent low-energy dispersions in K valley ( τ = 1) for both materials.Linear dispersions, corresponding to zero band gap, are also provided for comparison. The range of the wave vector is extendedto ± k (0) F for MoS , in plot ( c ), in order to display finite curvature of the dispersion curves which is markedly suppressed due toa large band gap parameter ∆ = 1 . eV . The DOS for TMDC’s represented in panel ( d ), is calculated for the gapped grapheneapproximation, given by Eq. (13), and for a general model (12), which displays substantial difference between the two results. ˆ H ξ,σ = (cid:18) − ξσ ∆ SO + ∆ z (cid:126) v F ( k x − ik y ) (cid:126) v F ( k x + ik y ) ξσ ∆ SO − ∆ z (cid:19) . (2)The energy dispersions are ε γξ,σ ( k ) = γ (cid:113) ( ξσ ∆ z − ∆ SO ) + ( (cid:126) v F k ) , (3)where γ = ± | ∆ SO − ξσ ∆ z | , which will be later referred to as ∆ <,> = | ∆ SO ∓ ∆ z | . Clearly, bothenergy gaps depend on the perpendicular electrostatic field and the two subbands corresponding to the ξσ = ± E z is related to a topological insulator (TI)states with ∆ z < ∆ SO . Once the two gaps become equal, we observe a metallic gapless state with ∆ < = 0 and afinite ∆ > , defined as valley-spin polarized metal (VSPM). For larger fields, ∆ z (cid:118) E z would always exceed the constantintrinsic spin-orbit gap ∆ SO , which corresponds to the standard band insulator (BI) state.The DOS which is in general defined by ρ d ( E ) = (cid:90) d k (2 π ) (cid:88) γ = ± (cid:88) ξ,σ = ± δ (cid:104) E − ε γξ,σ ( k ) (cid:105) , (4)is immediately obtained for silicene (see Appendix A) as ρ d ( E ) = 1 π (cid:88) γ = ± E(cid:126) v F (cid:88) i = <,> Θ (cid:20) E γ − ∆ i (cid:21) , (5)in terms of the unit step function Θ( x ). We note that for systems sharing the same Dirac cone characteristics witharbitrary energy gap, the DOS is linear analogous to graphene. Experimentally obtained linear V-shaped DOS wasused to verify the Dirac cone dispersion for germanene. However, ρ d ( E ) has a finite value only above the energy gapsince only for this energy range electronic states exist and tha is how the band gap plays an important role.Finite-temperature chemical potential for an electronic system is obtained using conservation of carrier density n of electrons ( n ( e ) ) and holes n ( h ) concentrations at all temperatures, including T = 0. In this regard, we have n = n ( e ) + ( − n ( h ) = ∞ (cid:90) d E ρ d ( E ) f γ =1 ( E , T ) − (cid:90) −∞ d E ρ d ( E ) { − f γ =1 ( E , T ) } . (6)At zero temperature, the density n is related to the Fermi energy E F in a straightforward way. If only one subbandis occupied, we obtain n = 12 π E F − ∆ < (cid:126) v F . (7)Alternatively, if the doping density is such that both subbands are populate, then E F is obtained from n = 1 π (cid:126) v F (cid:20) E F − (cid:0) ∆ < + ∆ > (cid:1)(cid:21) , (8)and the critical density required to start filling the upper subband is n c = 2∆ SO ∆ z /π (cid:126) v F .We prove in Appendix B that the finite-temperature chemical potentail could be obtained from the followingequation n (cid:18) (cid:126) v F k B T (cid:19) = (cid:88) γ = ± γπ (cid:88) i = <,> − Li (cid:26) − exp (cid:20) γµ ( T ) − ∆ i k B T (cid:21)(cid:27) + ∆ i k B T ln (cid:26) exp (cid:20) γµ ( T ) − ∆ i k B T (cid:21) (cid:27) , (9)where Li ( x ) is a polylogarithm function. Connecting the doping density n with the Fermi energy through eitherEq. (7) or Eq. (8) depending on how many subbands are doped at zero temperature, we derive the chemical potentialwith a value equal to E F at T = 0, for all accessible temperatures. Although this equation is transcendental andcannot be resolved algebraically, a quasi-analytic or one-step numerical solutions could be easily provided for anyfinite temperature without having to perform an integration.Clearly, the chemical potential for silicene depends on the energy band gaps ∆ i , i = <, > . Our approach, discussedin Appendix B, is valid for a variety of materials with linear energy dependence for the DOS, including M oS .Specifically, Eq. (9) also describes µ ( T ) for gapped graphene with two degenerate subbands, or ∆ < = ∆ > = ∆ . Forgapless pristine graphene ∆ = 0 and πn = [ E F / ( (cid:126) v F )] , we have12 ( k B T ) E F = − (cid:88) γ = ± γ Li (cid:26) − exp (cid:20) γ µ ( T ) k B T (cid:21)(cid:27) . (10)If the temperature is low with k B T (cid:28) E F , Eq. (10) is reduced to the expressions, derived in Refs. [70,73] and [80].All our considered materials could be effectively classified by the existing or broken symmetries of certain kinds, and,consequently, by the degeneracy of their energy subbands, which may be generally different for electrons and holes.In this respect, graphene represents the simplest case with a fourfold spin and valley degeneracy of ± (cid:112) ( (cid:126) v F k ) + ∆ states. Silicene and germanene dispersions, yet showing complete electron/hole symmetry, exhibit spin- and valley-dependent pairs of subbands, each being double degenerate. Finally, M oS demonstrates broken symmetry betweenits electrons and holes, and a finite energy separation between two non-equivalent holes subbands.Typical energy dispersions and DOS for silicene and molybdenum disulfide are presented in Fig. 1. For bothmaterials, we consider the K valley with ξ = 1 so that the upper electron and lower hole subbands correspond to σ = 1 spin. Every time a new subband, or their degenerate manifold, begins to be doped, we see an immediateincrease. or discontinuity, of the DOS, as schematically shown for silicene in the insets of Fig. 1 ( b ). It is importantto observe that for both silicene and gapped graphene, ρ d ( E ) is directly proportional to energy E , i.e., the DOS forDirac materials with finite and zero gap are exactly the same as long as we measure it above the enery band gap.With no electronic states inside the gap region, we have ρ d ( E <
0) = 0. ( ) e silicene ( ) f silicene ( ) c ( ) b ( ) d ( ) a D E G g r a p h e n e FIG. 2: (Color online). Finite temperature chemical potential µ ( T ) for graphene and silicene. Panel ( a ) shows µ ( T ) dependencefor a number of representative cases - 2DEG with ε (cid:118) k , ρ d ( E ) = const and no holes (black solid line), gapless graphene with ρ d ( E ) (cid:118) E (solid red curve) and three others with broken electron-hole symmetry. Dash-dotted blue line describes graphenewithout holes ρ ( h ) d ( E <
0) = 0, green dashed curve corresponds to graphene with reduced hole DOS ρ ( h ) d ( E <
0) = 1 / ρ ( e ) d ( E > ρ ( e ) d ( E <
0) = 1 / ρ ( d ) d ( E > b ) and ( c ) present µ ( T ) for gapped graphene for electron and hole doping, respectively, for chosen Fermi energy. Eachcurve is matched to a specific energy gap, as described schematically in the insets. Chemical potential for silicene with bothsubbands filled ∆ < = 0 . E F and different ∆ > is presented in panel ( d ). Black curve describes ∆ > = ∆ < = 0 . E (0) F (gappedgraphene), red and short-dashed line - ∆ < = 0 . E (0) F , dotted blue line - ∆ < = 0 . E (0) F and dashed green curve correspondsto ∆ < = 0 . E (0) F . Plot ( e ) combines two degenerate-subband situations - graphene (∆ SO = ∆ z = 0 . E F = E (0) F ), blacksolid curve and gapped graphene with ∆ SO = 0 . E (0) F , ∆ z = 0 . E F = 1 . E (0) F (red dashed curve) with the two casesof silicene with split subbands - ∆ SO = ∆ z = 0 . E (0) F and ∆ z = 0 . E F , E F = 1 . E (0) F (cid:39) ∆ > , and , finally, ∆ SO = 0 . E (0) F ,∆ z = 1 . E (0) F and E F = 1 . E (0) F . The subband schematics and the Fermi level for each case are given in the insets. Panel( f ) shows µ ( T ) at low temperatures satisfying k B T ≤ . E (0) F for a number of cases of doping level located close to the uppersubband ∆ > ≈ E F . Here, ∆ SO = 0 . E (0) F for each curve, ∆ z = 0 .
70, 0 .
75, 0 .
80 and 0 . E (0) F , while the Fermi levels at T = 0are 1 . . . .
422 and 1 . E (0) F . In Fig. 2, we display our results for the finite-temperature chemical potential for graphene and buckled honeycomblattices. First, we show how dissimilar this temperature dependence might be for various electronic systems. In Fig. 2( a ), this situation is described for a 2DEG with a parabolic energy band and constant DOS, graphene with ρ s ( E ) (cid:118) E ,as well as two model structures - graphene with no hole states, with doubly prevailing electrons ρ d ( E >
0) = 2 ρ d (( E <
0) and doubly prevailing hole DOS. At low temperatures, all the curves, except the one for the 2DEG, are nearlyidentical since the holes do net play any important role. The hole distribution function is complimentary to that forelectrons, i.e., 1 − f e ( (cid:15) ) →
0, and, therefore, inconsequential. However, when the temperature becomes comparablewith the Fermi energy, the hole thermal excitations become crucial, causing an opposite effect compared with that forelectrons. They mitigate the reduction of the chemical potential and eventually prevent µ ( T ) from crossing the zeroenergy level. This is seen particularly well for a hole-dominating system ( > h (+) ), for which the chemical potentialstarts to increase and ultimately exceeds the initial E F value. We conclude that only total electron/hole symmetry,but not the energy gap, keeps the chemical potential positive for arbitrary high temperatures.For the remaining plots, we consider the behavior of µ ( T ) for graphene and silicene with different gaps. At T = 0,we may keep the Fermi energy fixed so that the actual electron density n differs, in which case the states with a largergap receive much smaller amount. Alternatively, we can dope the sample and µ ( T ) shows much stronger reduction,or we can fix the carrier density n so that the Fermi energy increases in the case with larger gap (see Eq. (7)). Theformer case is shown in panels ( b ), ( c ) and ( d ), whereas the latter at ( e ) and ( f ). In general, the carrier density n isaccepted to be the most meaningful parameter, determining the Fermi energy for each specific system.Plots ( b ) and ( c ) are designed to demonstrate complete reflection symmetry of the chemical potential betweenelectron and hole doping for all energy band gaps. This is a manifestation of γ = ± T = 0, the Fermi level is chosen so that only thelower subband is filled. In contrast, the upper one ∆ > is located so close to the Fermi level E F , that it starts gettingpopulated at all, even very small temperatures. Thermally-induced doping, received by the ∆ > -subband, results in amuch stronger reduction rate of µ ( T ), as shown in Fig. 2 ( e ). Such situation suggests very special thermal propertiesof silicene with closely located energy subbands, and shares this behavior with graphene having an additional spinand valley degeneracy. We investigate this phenomenon even further by adjusting the upper subband in the vicinityof the Fermi level, as shown in Fig. 2 ( f ). Now, each µ ( T ) curve demonstrates a significantly pronounced decreasewhenever the upper subband filling becomes essential. These two different types of µ ( T ) behavior in silicene could beused to achieve its additional tunability by introducing the upper energy subband ∆ > . B. Transition metal dichalcogenides
The low-energy electronic states in monolayer molybdenum disulfide (
M oS , ML-MDS), a prototype transitionmetal dichalcogenide, could be effectively described by a two-band model Hamiltonian ˆ H ξ,σd = (cid:18) ξσ λ + (cid:126) k m e α (cid:19) ˆ I × + (cid:18) ∆2 − ξσ λ + (cid:126) k m e β (cid:19) ˆΣ z + t a ˆ Σ ξ · k , (11)whose important feature is a major gap parameter ∆ = 1 . eV which results in the actual band gap (cid:119) . eV , largecompared to that for silicene. The spin-orbit coupling parameter λ = 0 .
042 ∆ represents a smaller, but essentialcorrection, to the single-particle excitation spectrum and the band gap. The energy subbands are now spin- andvalley- dependent since the corresponding degeneracy is lifted. The electron hopping parameter t = 0 .
884 ∆ and a = 1 .
843 ˚A shape the Dirac cone term in the Hamiltonian Eq .(11) as t a = 4 . × − J · m , counting up to ≈ .
47 of (cid:126) v F value in graphene.Next, we include the (cid:118) k mass terms with α = 2 .
21 = 5 . β in which m e is the free electron mass. Ourconsidered values of Fermi momentum at zero temperature are determined by the experimentally allowed electronand hole doping densities n = 10 ÷ m − as k F = √ πn (cid:119) ÷ m − . Anisotropic trigonal warpingterm t a ( ˆ Σ ξ · k ) ˆ σ x ( ˆ Σ ξ · k ) is clearly beyond out consideration since (cid:118) t = 0 . eV = 0 .
053 ∆ term does not causeany effect on the electron dispersions. The energy dispersion relations, corresponding to the Hamiltonian Eq. (11) ε ξ,σγ ( k ) = (cid:15) ξ,σ ( k ) + γ (cid:114)(cid:104) ∆ ξ,σ ( k ) (cid:105) + ( t a k ) , formally represent gapped graphene with a k − dependent gap term∆ ξ,σ ( k ) = (cid:126) k β/ (4 m e ) + ∆ / − ξσ λ /
2, and a band shift (cid:15) ξ,σ ( k ) = (cid:126) k α/ (4 m e ) + ξσ λ / Neglecting only the O ( k ) terms leads us to ε ξ,σγ ( k ) (cid:119) ξσλ + α (cid:126) m e k + γ (cid:110)(cid:2) (2 t a ) + (∆ − ξσ λ ) β (cid:126) /m e (cid:3) k + (∆ − ξσ λ ) (cid:111) / . (12)This is the principal model which we will use to describe the energy dispersions of M oS in our work. The (cid:119) k terms, trigonal warping and anisotropy are considered to be non-essential, even though as we will later see, causecertain discrepancies in the DOS. We also show (see Appendix A) that the curvature of the energy subbands inTMDC’s is so small that even the highest possible doping density 10 m − results in the Fermi energies (cid:118) λ . Thus,at zero or low temperatures, we do not need to consider any high-energy corrections to Eq. (12) On the other hand,inclusion of higher order terms into the dispersions, would enormously complicate the DOS calculation.At high temperatures, the electronic states far from the Dirac point would also receive substantial temperature-induced doping density due to the so-called Fermi tail. In that case, our model Hamiltonian in Eq. (11) andespecially, simplified dispersion relation (12) would no longer provide a satisfactory approximation. Consequently, ourprimary focus is on small but finite temperatures for which the initial doping density and E F still play an importantrole, beyond the O (cid:0) T /T F (cid:1) approximation discussed in Ref. [70]. ( ) a )( i ( ) b -8.0-11.0 10.0 30.0 )( i FIG. 3: (Color online). Temperature-dependent chemical potential for MoS . Panel ( a ) presents the µ ( T ) dependencefor various initial electron doping densities - n = 2 . · m − (black solid line), n = 3 . · m − (red and dashed), n = 5 . · m − (dotted blue curve) and n = 1 . · m − (green short-dashed line). Plot ( b ) shows the corresponsingdependence for hole doping with the same densities. Insets ( i
1) and ( i
2) illustrates how the Fermi energy depends on theelectron and hole doping densities at zero temperature.
Here we note that for
M oS , similar to the buckled honeycomb lattices, spin and valley indices always appeartogether as a product, so that taking into account a × ν = ξσ could beeffectively introduced. We are going to use only the index ν for the rest of our consideration. A valley- and spin-resolved gapped graphene approximation of dispersions (12) arises once we also neglect the mass terms in Eq. (12) ε νγ ( k ) (cid:119) νλ / γ (cid:112) ( t a ) k + (∆ − νλ ) / . (13)This approximation has a few important advantages including simplicity and its formal resemblance with gappedgraphene so that all the crucial quantities such as the DOS, wave function, polarizability and many others are alreadyknown. Furthermore, it gives a quite a suitable description of the energy band structure of M oS , taking into accounta large gap parameter and ν -dependent splitting of the two hole subbands. Nevertheless, the mass terms must betaken into account for a proper evaluation of the DOS and most of the temperature-dependent properties of TMDC’s.As we demonstrated in Appendix A, even in the simplest possible parabolic subbands approximation, valid only forsmall wave vectors k →
0, the mass terms make a contribution, comparable with the Dirac cone and band gap parts ofHamiltonian (11). Very small curvature of the energy subbands, mentioned above results in a tremendous ρ d (cid:118) m eff ,so that even (cid:119) α/ (4 m e ) correction, hardly noticeable on the electron band structure, becomes significant. It basicallydiscards the DOS results obtained from the gapped graphene model, even at zero or small wave vectors. Importanceof the mass and even higher order terms for the plasmon calculation was discussed in Ref. [38].Taking into account all the required terms in (11), rigorous numerical calculations give accurate results for theelectron DOS for M oS . In Fig. 1 ( d ), we present all three possible outcomes. Based on the gapped graphene model,the DOS is nearly twice as large as its numerical values . In contrast, our (cid:118) k model (12) demonstrates quitea good match, especially in the low-energy range. We also note that the numerically obtained dependence is clearlylinear for a wide range of energies, much exceeding our considered diapason.In summary, we consider a piecewise linear approximation ρ d ( E ) relatively close ( δ E (cid:119) λ ) for each of the threenon-degenerate subbands ρ d ( E ) = c ( i )0 + c ( i )1 E and ρ d ( E ) = 0 in the gap region − ∆ / λ < E < ∆ /
2. The expansioncoefficients are obtained as c (1)0 = 2 . E (0) F / ( (cid:126) v F ) = 0 . t − a − , c (1)1 = − .
397 ( (cid:126) v F ) − = − .
308 ( t a ) − for γ = − ν = − ≤ − ∆ / − λ ≥ E ≤ − ∆ / − λ . Finally, when γ = −
1, but ν = 1 and E ≤ − ∆ / λ ,the hole DOS coefficients are c (2)0 = 1 . E (0) F / ( (cid:126) v F ) = 0 . t − a − c (2)1 = − .
767 ( (cid:126) v F ) − = − .
169 ( t a ) − .For electrons with γ = +1 at E = ∆ / δ(cid:15) , the two quasi-degenerate subbands lead to the DOS equal to c (3)0 =5 . E (0) F / ( (cid:126) v F ) = 0 . t − a − and c (3)1 = 0 .
815 ( (cid:126) v F ) − = 0 .
179 ( t a ) − .We adopt these values for ρ d ( E ), which arise from the numerical calculations in order to achieve the highest possibleprecision and credibility for our finite-temperature derivations. However, our effective model, presented in AppendixA, gives the DOS results which show good agreement with these numerical values and could be used for decisiveestimates for various collective calculations for M oS .Once the DOS is known, we are in a position to calculate the Fermi energy for a given doping density n e (0) . Thenew point here is that ρ d ( E ) is not directly proportional to the energy so that E F is determined by n e (0) = (cid:18) E F − ∆2 (cid:19) (cid:34) c (3)0 + c (3)1 (cid:18) ∆2 + E F (cid:19)(cid:35) , (14)or E ( e ) F = 12 c (3)1 (cid:40) − c (3)0 + (cid:20)(cid:16) c (3)0 + c (3)1 ∆ (cid:17) + 8 n e (0) c (3)1 (cid:21) / (cid:41) . (15)For hole doping, the Fermi energy differs from the previously considered electron doping case, i.e., E ( h ) F = 1 c (2)1 − c (2)0 + (cid:40) − c (2)1 n h (0) + (cid:20) c (2)0 − c (2)1 (cid:18) ∆2 − λ (cid:19)(cid:21) (cid:41) / (16)Numerically, our results for the Fermi energy for electrons and holes are presented in insets ( i
1) and ( i
2) of Fig. 3. Here,both linear and quadratic terms in the doping density equations are present (see Eq. (14)), and, most importantly, there is no symmetry between the electron and hole states. Unlike graphene, the linear dependence here dominatesfor both cases due to the large energy band gap. Each curve starts from the corresponding band gap - ∆ / . E (0) F for electrons and − ∆ / λ = − . E F . The corresponding well-known result for gapped graphene E F − ∆ =2 π n ( (cid:126) v F ) are verified by putting c ( i )1 → λ → M oS is obtained in a similar way, as we have done for the buckledhoneycomb lattices, except that we need to evaluate four different terms related to the two separate hole subbands(see Eq. (B26)). The corresponding numerical results are described in Fig. 3. As discussed above, the most specialproperty of TMDC’s is the broken electron/hole symmetry. Consequently, the chemical potential for the electrondoping becomes negative at T ≈ . E F , while µ ( T ) for hole doping does not ever change its sign. Broken electron/holesymmetry leads to two substantially different types of temperature dependence of the chemical potential. Thus, MoS represents a special material with unique symmetry properties and chemical potential dependence, so far encounteredonly in certain types of semiconductors but not in Dirac materials. III. PLASMON EXCITATIONS AT FINITE TEMPERATURE
As one of the most relevant applications of our finite-temperature chemical potential formalism, we briefly considerplasmons for an extrinsic, substantially doped at T = 0, free-standing gapped Dirac cone material. The plasmondispersion relation is calculated from the zeros of the dielectric function (cid:15) ( q, ω ), expressed in the random phaseapproximation (RPA) as (cid:15) ( q, ω ) = 1 − v ( q ) Π T ( q, ω | µ ( T ) , T, ∆ i ) = 0 , (17)where v ( q ) = 2 πe / ( (cid:15) s q ) is the Fourier-transformed two-dimensional Coulomb potential, and (cid:15) s = 4 π(cid:15) (cid:15) r ith (cid:15) b isthe background dielectric constant in which the 2D material is embedded. At finite temperature, the dynamicalpolarization function Π T is given as an integral transformation of its T = 0 counterpart Π , i.e.,Π T ( q, ω | µ ( T ) , T, ∆ i ) = 12 k B T ∞ (cid:90) dξ Π ( q, ω | ξ, ∆ i )1 + cosh { [ µ ( T, E F ) − ξ ] / ( k B T ) } (18)The evaluation of the zero-temperature polarizability is quite similar for buckled honeycomb lattices and M oS , sincein both cases their low-energy bandstructure is represented by two generally inequivalent double-degenerate pairsof subbands which depend on the composite index ν = σ ξ . For any such pair, Π ( ν )0 is obtained in the one-loopapproximation ( g = 2) asΠ ( ν )0 ( q, ω | E F , ∆ ν ) = g π (cid:90) d k (cid:88) γ,γ (cid:48) = ± (cid:18) γγ (cid:48) k · ( k + q ) + ∆ ν ε ν ( k ) ε ν ( | k + q | ) (cid:19) f [ E − ε νγ ( k )] − f [ E − ε νγ (cid:48) ( | k + q | )] ε νγ ( k ) − ε νγ (cid:48) ( | k + q | ) , (19)where f [ E − γ E ν ( k )] is the Fermi-Dirac distribution function, showing electron and hole occupation numbers forchosen energy E . At T = 0, it is a Heaviside unit step function Θ (cid:2) E − ε νγ ( k ) (cid:3) . We note that the extra 1/2 comes from0 -0.43 -0.22 0.0 ( ) a ( ) b -0.61 -0.31 0.0-0.83 -0.42 0.0-0.66 -0.33 0.0 ( ) c ( ) d FIG. 4: (Color online) Signle-particle excitation regions or particle-hole modes, outlined by non-zero Im Π T ( q, ω ) at an arbitrarytemperature for silicene with ∆ SO = 0 . E (0) F and ∆ z = 0 . E (0) F . The upper panels ( a ) and ( b ) describe the situation for zerotemperature, while the lower plots ( c ) and ( d ) were obtained for T = 0 . E (0) F . Left panels ( a ) and ( c ) are for a system with E F = 1 . E (0) F whereas the right ones - with E F = 1 . E (0) F . Alternatively, the regions of Im Π T ( q, ω ) = 0 specify plasmonexcitations with no Landau damping. the form factor 1 / γγ (cid:48) ... which also depends on the absolute value of each energy ε ν ( k ) = 1 /γ ε νγ ( k ). We alsonote that since any valley or spin transitions are inadmissible and only one summation over the index ν is incuded,compared to the electron/hole indices γ and γ (cid:48) . Thus, for both types of materials, the dynamical polarization functionis obtained as a sum of terms obtained from Eq. (19) over ν Π ( q, ω | E F , ∆ i ) = (cid:88) ν Π ( ν )0 ( q, ω | E F , ∆ ν ) . (20)Our results for the polarization functions and plasmon excitations are presented in Figs. 4 and 5.First, we need to address the imaginary part of the polarizability since it specifies the regions and intensity ofthe plasmon damping. We see that at a finite temperature the plasmon dissipation generally increases and thedamping-free regions are nearly absent. Once the temperature becomes high, the imaginary part of the polarizationfunction is reduced as 1 /T , so there is no uniform temperature dependence of the plasmon damping. Dopingand proportional increase of the energy band gaps, in contrast, increase the regions free of single-particle excitationspectrum (see Fig. 4 ( b )), so the effect of a finite temperature and E F on the plasmon dissipation coul be opposite.The real part of Π ( q, ω | E F , ∆ i ), shown in Fig. 5 ( a )-( d ), designates the location and slope of the correspondingplasmon branches which are presented in panels ( e ) and ( f ). Here, the temperature and the doping produce a similareffects, shifting the location of the peaks to the right. The real part of the polarization function must be positive inorder to enable a real solution of Eq. (17) Plasmon branches are located at higher energies for a given wave vectordue to either an initial increase of the Fermi energy or thermally-induced doping. A. Non-local plasmons in an open system
Concluding our investigation of extrinsic 2D materials, we now turn to the plasmon excitations in so-called 2Dopen systems (2DMOS). Such a nanoscale hybrid arrangement being a part of graphene-based nanoscale devices, ( ) a -1.00.0-0.5 ( ) b ( ) c ( ) d -1.0-0.5 0.50.0 1.0 2.51.5 3.00.0 0.5 1.0 2.02.01.00.50.0 ( ) e ( ) f FIG. 5: (Color online) Dynamic polarization function and plasmon excitations for silicene with ∆ SO = 0 . E (0) F and ∆ z =0 . E (0) F . Panels ( a ) - ( d ) present constant frequency cuts of Re Π T ( q, ω ) as a function of q for Ω = 0 .
0, 0 .
5, 1 .
0, 1 . . E (0) F / (cid:126) , described by solid black, dashed red, short-dashed blue, long-dashed green and dash-dotted orange lines, respectively,for each plot. Panels ( a ) and ( b ) show the zero temperature case, ( c ) and ( d ) - T = 0 . E (0) F . Initial ( T = 0) doping is chosen toyield E F = 1 . E F T (0) for plots ( a ) and ( c ) while panels ( b ) and ( d ) present the case of E F = 1 . E F T (0) . Plasmon dispersionsat T = 0 . E (0) F /k B for E F = 1 . E (0) F is represented in plont ( e ), and for E F = 1 . E (0) F - in panel ( f ). consists of a 2D layer (it could be 2DEG, graphene, a buckled honeycomb layer, or M oS ), which is Coulomb coupled(not chemically bound) with a semi-infinite conductor and its surface plasmon. While a plasmon excitation in aclosed system is determined by two-particle Greens functions, in 2DMOS it more involved, depending on the Coulombinteraction with the environment.
The important feature of such a system is the screened Coulomb couplingbetween the electrons in graphene and the conducting substrate. Such a screened potential could be obtained usingthe nonlocal frequency-dependent inverse dielectric function
Consequently, two or more linear plasmonbranched have been obtained, confirming some previous experimental claims. For finite temperature, theircoupling to an external reservoir is reflected in existence of extra plasmon dissipation channels. The plasmon branches in such system are obtained as zeros of a so-called structure factor S ( q, ω | µ ( T ) , ∆ i ), playingthe role of the dielectric function (cid:15) ( q, ω ) in an isolated layer. It is obtained as S ( q, ω | µ ( T ) , ∆ i ) = 1 − v ( q ) Π (0) ( q, ω ; µ ) (cid:26) − (cid:15) B ( ω )1 + (cid:15) B ( ω ) Exp ( − a q ) + 1 (cid:27) , (21)where a is the distance between the 2D layer and the surface and the bulk dielectric function is given in the locallimit as (cid:15) B = 1 − Ω p /ω . The bulk-plasma frequency, defined as Ω p = ( n m e ) / ( (cid:15) o (cid:15) S m ∗ ), depends on on the electron2 ( ) a ( ) b FIG. 6: (Color online) Temperature-dependent non-local plasmon excitations for a silicene-based hybrid system with ∆ SO =0 . E (0) F and ∆ z = 0 . E (0) F . Density plots of the real part of S c ( T, ω + 0 + ) for a constant wave vector q = 0 . k (0) F are presented,so that the sought low-damped plasmons correspond to their peaks. Panel ( a ) corresponds to an intrinsic system with nozero-temperature doping, while plot ( b ) is related to E F = 2 . E (0) F . concentration n m , its effective mass m ∗ and the substrate dielectric constant (cid:15) S . This approximation stays valid fora large range of wave vectors q (cid:28) × m − since the Fermi wavelength in metals is comparable with the inverselattice constant. As a result the frequency of the upper plasmon branch in our system, equal to Ω p / √ q (cid:107) → in case of spin- and valey-dependent single-particle excitations in a 2D layer (whichis true for buckled honeyomb lattices and M oS ), such a hybrid structure could be effectively used to directlymeasure the dielectric properties or spin-orbit coupling parameters of such a layered material because the location ofeach plasmon branch, its damping rate and the signatures of the particle-hole modes are independently determined bythe material parameters of the 2D layer. Each of these properties is unique (for example, the two plasmon branches insilicene-based TDMOS depends on the energy band gaps as ∆ / i and ∆ / i , while the outermost PHM boundaries aredetermined solely by the lower gap ∆ < ), so that an additional linear plasmon branch provides us with the requiredearlier unknown piece of information about the specific material. This is not possible in the case of a single plasmonbranch in an isolated layer.In the present work, we additionally introduce finite doping and temperature into this consideration. Our numericalresults for the non-local plasmons for a extrinsic systems are presented in Fig. 6 We see that the location and dampingof both branches substantially changes in the presence of an initial Fermi energy. At low temperature, this increaseis especially apparent, similar to plasmons in an isolated layer. For a sample with stronger doping, E F = 1 . E (0) F atzero temperature, the upper plasmon branch is always located above the surface plasmon level Ω p / √ IV. CONCLUDING REMARKS
We have carried out an extensive investigation of extrinsic, or doped, Dirac gapped materials at arbitrary finitetemperatures and obtained a set of algebraic analytic equations determining the chemical potential. Our consideredsystems include graphene, with or without an energy band gap, buckled honeycomb lattices with spin- and valley-dependent energy subbands and reduced degeneracy, as well as the transition metal dichalcogenides having brokensymmetry between the electron and hole states. Our results could also be used to predict the finite-temperaturechemical potential for parabolic or quasi-parabolic eigenstates in semiconductors with light or heavy holes.
Ingeneral, our model is limited only by the linear dependence of the DOS which stays valid over a wide energy rangefor all the above mentioned materials.We have demonstrated that the chemical potential depends substantially on the energy band gap(s) of the consideredsystem since the DOS depends on the curvature of each subband. Specifically, we investigated structures with twonon-degenerate, separated spin- and valley- resolved energy subbands in both valence and conduction bands, such assilicene. The upper subband would receive thermally-induced doping even if it is undoped at zero temperature. Thisis always reflected in a higher reduction rate of the chemical potential whenever the second subband doping startsplaying a role. Consequently, one can tune the µ ( T ) dependence around any required temperature by bringing initialdoping close to the higher-energy subband. The number of such separated subbands contributing to the DOS could3be arbitrarily large for an electron in a quantum well or quasi-one-dimensional nanoribbons, and all these casescould be effectively treated by our model.The behavior of the chemical potential depends on whether there is symmetry between the electron and hole statesin the system. If the DOS for the electrons and holes is equal or symmetric around the Dirac point, then the chemicalpotentail does not change its sign even for high temperatures, i.e., it remains positive for electron doping or staysnegative if E F <
0. In fact, these two types of doping result in symmetric behavior, decreasing | µ ( T ) | , so that theelectron-hole symmetry persists at any finite temperature. Let us discuss the electron doping in more detail. Oncethe temperature is sufficiently high, the hole states become thermally excited and it has an opposite effect to that ofelectron, decreasing the reduction of the chemical potential. At extremely high temperatures, the two processes havealmost equivalent effects and µ ( T ) asymptotically tends to zero and never crosses the Dirac point. This situationchanges if the DOS of the electron and holes differs and the two hole subbands are energetically not equivalent, aswe observe for TMDC’s. We have shown that for MoS the chemical potential becomes negative, changing its signat T (cid:119) . E F . Alternatively, µ ( T ) could never reach the zero energy line if there is hole doping at zero temperatureand starts decreasing at sufficiently high temperatures.As a necessary intermediate step in our derivations, we obtained a piecewise-linear model for the DOS for transitionmetal dichalcogenides, directly from the Hamiltonian parameters of the considered system. This model gives exactresults for ρ d ( E ) at the band edges (next to each energy gap) and a fairly good approximation at higher energies. Thismodel significantly improves the results obtained from the spin- and valley-resolved gapped graphene approximation.Finally, we considered the way in which the initial doping affects the plasmon dispersions ( (cid:126) Ω p /E F ) (cid:119) Λ q , andthe effective length Λ depends on the doping. There could be an initial carrier density at T = 0, and the thermally-induced one at finite temperature. The latter doping type results in a finite polarization function and √ q T plasmondispersions even for intrinsic systems. We demonstrate how much each type of doping contributes this effectivelength.Several many-body calculations or collective electronic models requires reliable knowledge of the chemical potentialat finite temperature. In the absence of such information, much attention has been directed towards intrinsic, orundoped systems, or low temperatures, so that (cid:119) ( T /T F ) series expansions could have been applied. Our results aregoing to provide considerable assistance in transport studies, optical, thermally modulated conductivity for the mate-rials, discussed in our paper. Consequently, we expect our work to provide an important contribution to electronics,transport and plasmonics of these recently discovered structures, both theory and experiment. Acknowledgments
D.H. would like to thank the support from the Air Force Office of Scientific Research (AFOSR).
Appendix A: Density-of-states
The DOS for electrons and holes with energy dispersion ε ξ,σγ ( k ) is defined as ρ d ( E ) = (cid:90) d k (2 π ) (cid:88) γ = ± (cid:88) ξ,σ = ± δ (cid:2) E − ε ξ,σγ ( k ) (cid:3) , (A1)where ξ = ± σ = ± σξ , sothat a single composite index ν = σξ = ± (cid:80) ξ,σ = ± = ⇒ · (cid:80) ν = ± , which we will use throughout this work.For a number of cases of 2D structures, Eq.(4) could be immediately evaluated using the following property of adelta function δ ( f ( x )) = (cid:88) i δ ( x − x i ) | df ( x ) /dx x = x i | , (A2)where x i are the roots of f ( x ), formally given as f ( x ) = E − (cid:15) γν ( k ) for various energy dispersions (cid:15) γν ( k ). In the case ofsilicene and germanene, the result is straightforward with4 ρ d ( E ) = 1 π (cid:88) γ = ± γ E(cid:126) v F (cid:88) i = <,> Θ (cid:20) E γ − ∆ i (cid:21) , (A3)i.e., the DOS for Dirac gapped systems is linear similar to graphene. However, it is finite only above the energy gaps.This result also covers the case for gapped graphene if the two gaps are equal ∆ < = ∆ > = ∆ . Furthermore, wearrive at well-known V-shaped (cid:118) E DOS for pristine gapless graphene if ∆ <,> → Molybdenum disulfide
In our work, we discuss several effetive models of different complexity and accuracy, describing the energy dispersion ε νγ ( k ) for M oS . First, we consider a spin- and valley- resolved gapped graphene model given by Eq. (13), in which weleave out all the mass (cid:118) α and (cid:118) β terms. It gives quite accurate results for the energy eigenstates next to the cornersof Brillouin zone, as shown in Fig. 1. However, it is straightforward to see that the corresponding DOS obtained as ρ d ( E ) = 1 π ( t a ) (cid:88) γ = ± γ (cid:88) ν = ± (cid:16) E − ν λ (cid:17) Θ (cid:20) γ (cid:18) E − νλ (cid:19) −
12 (∆ − νλ ) (cid:21) , (A4)is V-shaped and does not match the numerical results even near the band edges. Once we get into the ”allowed“energy ranges outside the band gap ε ν = − , γ =1 ( k = 0) = ∆ / ε ν = − γ = − ( k = 0) = − ∆ / λ and ε ν =1 γ = − ( k =0) = − ∆ / − λ for the holes, the DOS experiences three different giant discontinuities due to each new contributingenergy subband. These giant leaps could be calculated using the parabolic subbands approximation, obtained for k (cid:28) k F ε νγ ( k ) = 12 [ νλ (1 − γ ) + γ ∆] + (cid:20) (cid:126) m e ( α + γβ ) + γ ( t a ) ∆ − νλ (cid:21) k . (A5)This result leads to the DOS given by ρ d ( E ) = 12 π (cid:126) (cid:88) γ, ν = ± (cid:12)(cid:12)(cid:12) α + γβ m e + γ ( t a ) (cid:126) (∆ − νλ ) (cid:12)(cid:12)(cid:12) − Θ (cid:20) γ (cid:18) E − νλ (cid:19) −
12 (∆ − νλ ) (cid:21) . (A6)The two terms in Eq. (A6) are of the same order of magnitude, consequently, each of them must be retained in ourcalculation. Physically it means that due to large gap ∆ the curvature of each subband at k = 0 is so small that eventhe (cid:119) α/ (4 m e ) correction is significant. It basically discards the DOS obtained for the gapped graphene model, evenat k →
0. The fact that the mass and even higher order terms must be taken into account for the plasmon calculationwas mentioned in Ref. [38].The principal model which yields quite reliable results for the DOS is derived by neglecting only non-essential O ( k )terms with ε νγ ( k ) (cid:119) ν λ + (cid:126) α m e k + γ (cid:115) (∆ − ν λ ) + (cid:20) (2 t a ) + (∆ − ν λ ) (cid:126) βm e (cid:21) k . (A7)Here, all (cid:118) k terms are retained and the final expression for the DOS appears to be quite complicated. We use ageneral equation from our previous work to come up with a linear approximation which, according to the mostprecise and generalized numerical results, is valid for all experimentally allowable electron and hole dopingdensities. For the dispersions Eq. (A7), we use Eq. (A1) to obtain ρ d ( E ) = 12 π (cid:88) j (cid:88) γ, ν = ± (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ α + γ ˜ A ν (∆ ± λ , β | a t )2 (cid:110) E − ˜ (cid:15) ν − ˜ α ξ ( j ) ν ( E ) (cid:111) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − Θ (cid:20) γ (cid:18) E − µλ (cid:19) −
12 (∆ − νλ ) (cid:21) , (A8)5where ˜ (cid:15) ν = νλ /
2, ˜∆ ν = (∆ − νλ ) /
2, ˜ A ν (∆ ± λ , β | a t ) = (∆ − νλ ) (cid:126) β/ (4 m e ) + ( t a ) , and ˜ α = (cid:126) α/ (4 m e ). ξ ( j ) ν ( E ) are the roots of f ( E , ξ ) = E − ˜ (cid:15) ν − αξ − γ (cid:113) ˜ A ν ξ + ˜∆ ν = 0 . (A9)This euqation could be solved by expressing it in quadratic form. However, one must bear in mind that if both partsof Eq. (A9) are squared, there might be additional non-physical solutions which must be disregarded. If this equationis written as ( αξ ) + B ξ + C = 0, where B = ˜ A ν + 2 α and C = ( E − ˜ (cid:15) ν ) − ˜∆ ν , the only appropriate solution is ξ (1) ( E ) = 1 / (2 α ) (cid:0) B + √ B − α C (cid:1) . The other solution corresponds to the E − ˜ (cid:15) ν − αξ = − γ (cid:113) ˜ A ν ξ + ˜∆ ν and isobviosly incorrect for α →
0. We also note that the electron/hole index γ is no longer present in this equation, sothat the two ± solutions are not associated with electron or hole states.In order to illustrate the physics behind selecting the only appropriate solution, let us consider a simple example ofgapless graphene with ε γ ( k ) = (cid:126) v F k with additional small, not depending on γ , the mass term (cid:118) αk , α (cid:28) (cid:126) v F /k F .The actual dispersion relation is now (cid:15) γ ( k ) = γ (cid:126) v F | k | + αk , and the DOS is calculated as ρ d ( E ) = 2 π (cid:88) γ = ± (cid:88) j k ( j ) | γ (cid:126) v F + 2 αk ( j ) | . (A10)The roots k ( j ) are the solutions of γ (cid:126) v F | k | + αk = 0. Even though such a quadratic equation generally has twoinequivalent solutions, only one of them ± k (1) (cid:119) E / ( γ (cid:126) v F ) − α E / ( (cid:126) v F ) satisfies the | k | -type equation. Consequently,we obtain the following expression ρ d ( E ) (cid:119) π (cid:88) γ = ± E γ ( (cid:126) v F ) − α (2 + γ ) E ( (cid:126) v F ) (A11)becoming equivalent to graphene DOS ρ d ( E ) = 2 /π E / (cid:2) γ ( (cid:126) v F ) (cid:3) for α →
0. For the holes with γ = −
1, however, thelinear and quadratic mass terms are competing, so that for k (cid:29) k F , another solution is present. However such wavevectors are beyond the Dirac cone model, and therefore, it is not physically acceptable. The obtained correction tothe DOS is a small decrease, as it is expected to be for an energetically elevated location of the subband. However,in our model for M oS , the mass terms (cid:118) α and (cid:118) β are not small and represent a finite correction to the DOS. Thesmall parameter which we used in our series expansions is the energy δ(cid:15) above each band gap.Now, we return to Eq. (A9) and present its solution as ξ (1) = 12 ˜ α (cid:26) ˜ A ν + 2 ˜ α (cid:16) E − ˜ ε (0) ν (cid:17) − (cid:104) ˜ A ν + 4 ˜ α ˜∆ ν + 4 ˜ α ˜ A ν (cid:16) E − ˜ ε (0) ν (cid:17)(cid:105) / (cid:27) . (A12)This solution is exact in the sense that no approximations have been made so far excpet for the (cid:118) k dispersions(A7). Substituted into Eq. (A8), it gives the DOS for an arbitrary energy, for both electrons and holes.As the next step, we substitute this result for ξ (1) into Eq. (A8). We are interested in obtaining a linear approxi-mation of the DOS next to each subband edge. Let us first consider electrons with ε ν ( k ) = ∆ / δ(cid:15) , δ(cid:15) (cid:28) E . In thiscase, ξ (1) ν (cid:119) m e α (cid:126) (cid:26) − ( a t ) + (cid:126) β/ (4 m e ) (∆ − νλ )( a t ) + (cid:126) / (4 m e ) [( α + β ) (∆ − νλ )] (cid:27) δ(cid:15) , (A13)and the DOS is now approximately given by ρ d ( E ) = 12 π (cid:88) ν = ± (A14)∆ − νλ ( a t ) + (cid:126) / (4 m e ) ( α + β )(∆ − νλ ) + 2 δ(cid:15) (cid:2) ( a t ) + (cid:126) β/ (4 m e ) (∆ − νλ ) (cid:3) { ( a t ) + (cid:126) / (4 m e ) [( α + β ) (∆ − νλ )] } . ρ d ( E ) = c (3)0 + c (3)1 (cid:18) E − ∆2 (cid:19) ,c (3)0 = 0 .
180 1 t a = 11 . E (0) F ( (cid:126) v F ) ,c (3)1 = 0 .
268 1( t a ) = 1 .
218 1( (cid:126) v F ) . (A15)In the valence band, we consider two separate hole subbands with ν = ±
1. If γ = − ν = 1 and E ≤ − ∆ / λ , theDOS is ρ d ( E ) = c (2)0 + c (2)1 (cid:2) E − (cid:0) ∆2 − λ (cid:1)(cid:3) , and the expansion coefficients are c (2)0 = 12 π ∆ − λ ( a t ) + ( β − α )(∆ − λ ) , (A16) c (2)1 = 1 π δ(cid:15) (cid:2) ( a t ) + (cid:126) β/ (4 m e ) (∆ − λ ) (cid:3) { ( a t ) + (cid:126) / (4 m e ) [( β − α ) (∆ − λ )] } < , or c (2)0 = 0 .
105 1 t a = 6 . E (0) F ( (cid:126) v F ) ,c (2)1 = − .
232 1( t a ) = − .
051 1( (cid:126) v F ) . (A17)Finally, for the lower hole subband with E (cid:47) − ∆ / − λ , we obtain ρ d ( E ) = c (2)0 + c (2)1 (cid:20) E − (cid:18) ∆2 + λ (cid:19)(cid:21) , (A18) c (1)0 = 12 π (cid:88) ν = ± ∆ − ν λ ( a t ) + ( β − α )(∆ − ν λ ) ,c (1)1 = 1 π (cid:88) ν = ± δ(cid:15) (cid:2) ( a t ) + (cid:126) β/ (4 m e ) (∆ − ν λ ) (cid:3) { ( a t ) + (cid:126) / (4 m e ) [( β − α ) (∆ − ν λ )] } < , and c (1)0 = 0 .
233 1 t a = 15 . E (0) F ( (cid:126) v F ) ,c (1)1 = − .
458 1( t a ) = 2 .
077 1( (cid:126) v F ) . (A19)It is straightforward to obtain our previous results for the gapped graphene model (A4) DOS if α = β →
0. We alsonote that the slope of the DOS in the conduction band is negative, as it should be accoring to Fig. 1 ( d ), and thesummation over the ν index is present in all cases excpet the upper hole subband in Eqs. (A16).Our results (A15) - (A19) (here, we move from the conduction electrons to the valence band, i.e., from the right toleft) represent a fairly good match with the previously obtained numerical values, specified in Sec. II and later usedfor all our finite temperature calculations. The coefficients c ( i )0 , i = 1 , , δρ d ( E ) at each subband edge or k = 0, except c (1)0 (cid:119) δρ (2) d ( − ∆ / λ ) + δρ (1) d ( − ∆ / − λ )), and, therefore,are accurate. The linear coefficients c ( i )1 , in fair agreement, are 20 −
25% larger compared with the numerical resultssince all the (cid:118) k terms of our energy dispersions are neglected. Inclusion of these terms leads to the higher energies7for chosen wave vector and a decrease of the DOS. This discrepancy is increased for higher energies, which is well seenfor c (1)1 for holes with E < − ∆ / − λ . However, the subbands for such energy ranges do not receive any substantialdoping unless the temperature becomes very high. For most considered situations, we are limited for δ(cid:15) ≈ λ withinthe band edges. In such a small range, c ( i )1 δ(cid:15) (cid:28) c ( i )0 , so that the actual DOS values remain almost unaffected and ourmodel yields accurate results. Appendix B: Chemical potential µ = µ ( T ) at a finite temperature We now derive a set of algebraic equations for the finite-temperature chemical potential µ ( T ). At zero temperature,it is equal to the Fermi energy E F = µ ( T ) (cid:12)(cid:12) T =0 . Our derivation is based on the total carrier density conservation,which includes both electrons and holes, for zero and any finite temperatures n = n ( e ) + ( − n ( h ) = ∞ (cid:90) d E ρ d ( E ) f γ =1 ( E , T ) − (cid:90) −∞ d E ρ d ( E ) { − f γ =1 ( E , T ) } . (B1)The electron and hole occupation probabilities are complimentary and for electron doping at T = 0 the hole statesterm has no effect on Eq. (B1) (for details see Ref. [100]).We begin with the relatively simple case for silicene with dispersions (3). The DOS ρ d ( E ) for which buckledhoneycomb lattices is given by Eq. (4). The expression for the Fermi energy E F for fixed electron doping density n at zero temperature depends on whether either one or both electron subbands are doped. The former case occurs fordoping densities n ≤ n c = 12 π ∆ > − ∆ < (cid:126) v F = 2 π (cid:126) v F ∆ SO ∆ z , (B2)and the Fermi energy is obtained form n = 12 π E F − ∆ < (cid:126) v F . (B3)Alternatively, if the doping density is sufficient to populate both subbands, E F is determined by n = 1 π (cid:126) v F (cid:20) E F − (cid:0) ∆ < + ∆ > (cid:1)(cid:21) . (B4)Once the temperature is set finite, Eq. (B1) leads us to n = I ( e ) (∆ i , T ) − I ( h ) (∆ i , T ), with the two terms corre-sponding to electron and hole components of the total carrier density. These integrals are presented as I ( e ) (∆ i , T ) = ∞ (cid:90) ∆ < d E A ( E , T ) + ∞ (cid:90) ∆ > d E A ( E , T ) , (B5)where A ( E , T ) = 1 π E(cid:126) v F (cid:26) exp (cid:20) E − µ ( T ) k B T (cid:21)(cid:27) − . (B6)Each of the these integrals could be easily evaluated. Using a variable substitution ξ = ( E − ∆ < ) /k B T , we obtain1 π (cid:126) v F ) ∞ (cid:90) ∆ < d E E (cid:26) exp (cid:20) E − µ ( T ) k B T (cid:21)(cid:27) − = k B Tπ ( (cid:126) v F ) ∞ (cid:90) dξ (∆ < + ξ k B T ) (cid:26) exp (cid:20) ξ − µ ( T ) − ∆ < k B T (cid:21)(cid:27) − (B7)8With the help of the following notation R ( p ) ( T, X ) = ∞ (cid:90) dξ ξ p / { exp [ ξ − X/ ( k B T )] } , (B8)we obtain the final result of the integration as k B T ∆ < π ( (cid:126) v F ) R (0) [ T, µ ( T ) − ∆ < ] + 1 π (cid:18) k B T (cid:126) v F (cid:19) R (1) [ T, µ ( T ) − ∆ < ] . (B9)For p = 0 and 1, corresponding to the 2DEG and gapless graphene, Eq. (B8) leads to R (0) ( T, X ) = ln (cid:26) exp (cid:20) Xk B T (cid:21)(cid:27) , R (1) ( T, X ) = − Li (cid:26) − exp (cid:20) Xk B T (cid:21)(cid:27) , (B10)where Li ( z ) is the second-order polylogarithm function or dilogarithm defined asLi p ( z ) = ∞ (cid:88) k =1 z k k p , Li ( z ) = − z (cid:90) ln(1 − t ) t dt . (B11)The second term of Eq.(B5), which only differs from the first one by its integration limits, is2 k B T ∆ > π ( (cid:126) v F ) ln (cid:26) exp (cid:20) µ ( T ) − ∆ > k B T (cid:21)(cid:27) − π (cid:18) k B T (cid:126) v F (cid:19) Li (cid:26) − exp (cid:20) µ ( T ) − ∆ > k B T (cid:21)(cid:27) . (B12)The remaining term I ( h ) (∆ i , T ), i = { <, > } , which describes the contribution from the holes, is also easily obtained I ( h ) (∆ i , T ) = I ( h )1 (∆ i , T ) + I ( h )2 (∆ i , T ) where , I ( h )1 (∆ i , T ) = 1 π k B T ( (cid:126) v F ) (cid:88) i = <,> ∆ i R (0) { T, − [ µ ( T ) + ∆ i ] } and I ( h )2 (∆ i , T ) = 1 π (cid:18) k B T (cid:126) v F (cid:19) (cid:88) i = <,> R (1) { T, − [ µ ( T ) + ∆ i ] } . (B13)Now, the total carrier denisty from Eq. (B1) could be written as n = (cid:18) k B T (cid:126) v F (cid:19) (cid:88) γ = ± γπ (cid:88) i = <,> R (1) [ T, γµ ( T ) − ∆ i ] + ∆ i k B T R (0) [ T, γµ ( T ) − ∆ i ] (B14)or, explicitly expressing the polylogarithm functions, we write n (cid:18) (cid:126) v F k B T (cid:19) = (cid:88) γ = ± γπ (cid:88) i = <,> − Li (cid:26) − exp (cid:20) γµ ( T ) − ∆ i k B T (cid:21)(cid:27) + ∆ i k B T ln (cid:26) exp (cid:20) γµ ( T ) − ∆ i k B T (cid:21) (cid:27) . (B15)9Using Eqs. (B3) or Eq. (B4) depending on whether only one or both subbands are filled at zero temperature, weobtain the equation which relates the finite-temperature chemical potential with its T = 0 value E F . Energy bandgap(s) obviously affects this result.The µ ( T ) for gapped graphene with two fourfold degenerate energy subbands is obtained if we substitute ∆ < =∆ > = ∆ and (cid:80) i = <,> = ⇒ ×
2. For gapless graphene, ∆ = 0 and πn = [ E F / ( (cid:126) v F )] , so that we write12 ( k B T ) E F = (cid:88) γ = ± γ R (1) [ T, γµ ( T )] = − (cid:88) γ = ± γ Li (cid:26) − exp (cid:20) γ µ ( T ) k B T (cid:21)(cid:27) . (B16)If the temperature is kept low with k B T (cid:28) E F , this result is simplified as R ( x ) (cid:119) (cid:18) x π (cid:19) Θ( x ) + x ln (cid:16) e −| x | (cid:17) , (B17)Finalizing our derivations for silicene, we briefly address the case of hole doping with E F <
0. The left part ofEq. (B1) is now modified as − n ( h ) = − π (cid:126) v F ) (cid:88) i = <,> − ∆ i (cid:90) −∞ d E | E | Θ( − E + E F ) . (B18)In analogy with electron doing, the Fermi energy depends on whether only the ∆ < -subband (which is now the higherone) or they are both doped. The equations determining the Fermi energy for given hole doping density n are exactlysimilar to Eqs. (B3) and (B4), which confirms complete symmetry between the electron and hole states in silicene. Theright part of Eq. (B1) remains unchanged except the chemical potential is negative µ < µ ( T ) for transition metal dichalcogenides In cotrast to the previously considered buckled honeycomb lattices and graphene, the electron/hole symmetry inTMDC’s (such as
M oS ) is clearly broken. Even in the simplest gapped graphene model given by Eq. (13), the twohole subbands are not degenerate and separated by λ at k = 0. For all reasonable doping densities n < m − , theFermi energy is such that the lower ε γ = − ν =1 ( k = 0) = − ∆ / − λ subband could not be populated at zero temperature.This could be verified by rewriting Eq. (B2) as n c = 2 π λ ∆( t a ) = 1 . · m − . (B19)From here on, we are going to use a picewise-linear model for the DOS with the empirical coefficients, providedin Sec. II. Our analytical model for the DOS, develop in the preceding Appendix in Eqs. (A15) - (A19) could alsobe employed here without losing much of precision. Let’s us first consider electron doping with density n e (0) at zerotemperature. The corresponding Fermi energy is determined by n e (0) = c (3)1 / (cid:0) E F − ∆ / (cid:1) + c (3)0 ( E F − ∆ / , (B20)or E eF = 1 c (3)1 − c (3)0 + (cid:34)(cid:18) c (3)0 + c (3)1 ∆2 (cid:19) + 2 n e (0) c (3)1 (cid:35) / . (B21)For hole doping the result is quite similar, except c (2)1 < − ∆ / λ :0 E ( h ) F = 1 c (2)1 − c (2)0 + (cid:40) − c (2)1 n h (0) + (cid:20) c (2)0 − c (2)1 (cid:18) ∆2 − λ (cid:19)(cid:21) (cid:41) / . (B22)Using this expression, we can further improve the result in Eq. (B19) for the critical hole doping density n c neededto reach the lower subband at zero temperature. The required Fermi energy is E ( h ) F ≤ − ∆ / − λ , so that thecorresponding hole density is n hc = − λ ∆ c (2)1 + 2 c (2)0 λ = 3 . · m − .This critical density value is still far above the experimentally realizable values (cid:119) . · m − , and for all ourcalculations, the lower hole subband is never populated at T = 0.Finally, we address the finite-temperature chemical potential for M oS . Once again, we start with the carrierdensity conservation equation (B1).The electron component is easily evaluated as k B T (cid:18) c (3)0 + ∆2 (cid:19) R (0) (cid:20) T, µ ( T ) − ∆2 (cid:21) + c (3)1 ( k B T ) R (1) (cid:20) T, µ ( T ) − ∆2 (cid:21) . (B23)The hole state integrals are written in the following form (cid:90) −∞ d E ρ d ( | E | ) { − f ( E , T ) } = − ∆ / λ (cid:90) −∞ d E (cid:104) − c (2)1 E + c (2)0 (cid:105) (cid:26) exp (cid:20) µ ( T ) − E k B T (cid:21)(cid:27) − + (B24)+ − ∆ / − λ (cid:90) −∞ d E (cid:104) − δc (1)1 E + δc (1)0 (cid:105) (cid:26) exp (cid:20) µ ( T ) − E k B T (cid:21)(cid:27) − , where δc ( i )1 = c ( i )1 − c ( i )2 , i = 1 ,
2. They are evaluated as I ( h ) (∆ , λ | T ) = (cid:88) j =1 I ( h ) j where , I ( h )1 (∆ , λ | T ) = k B T (cid:18) ∆2 − λ + c (0)2 (cid:19) R (0) (cid:26) T, − (cid:20) µ ( T ) + ∆2 − λ (cid:21)(cid:27) , I ( h )2 (∆ , λ | T ) = c (1)2 ( k B T ) R (1) (cid:26) T, − (cid:20) µ ( T ) + ∆2 − λ (cid:21)(cid:27) , (B25) I ( h )3 (∆ , λ | T ) = k B T (cid:18) ∆2 + λ + δc (0)1 (cid:19) R (0) (cid:26) T, − (cid:20) µ ( T ) + ∆2 + λ (cid:21)(cid:27) , I ( h )2 (∆ , λ | T ) = δc (1)1 ( k B T ) R (1) (cid:26) T, − (cid:20) µ ( T ) + ∆2 + λ (cid:21)(cid:27) . Combined with the electron terms (B23), these hole integrals (B26) form the right side of the carrier density conser-vation (B1). Its left side, corresponding to zero temperature, is given by Eq. (B20) for electron doping and by n e (0) = (cid:18) ∆2 + E F − λ (cid:19) (cid:40) c (0)2 + c (1)2 (cid:20) ∆2 − ( E F + λ ) (cid:21)(cid:41) , (B26)for hole doping. The symmetry between the electron and hole states is no longer present, which strongly affects thefinite-temperature behavior of the chemical potential in TMDC’s. ∗ Electronic address: [email protected] T. Ando, A. B. Fowler, and F. Stern, Reviews of Modern Physics , 437 (1982). F. Stern, Phys. Rev. Lett. , 546 (1967). K. Novoselov, A. K. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. Grigorieva, S. Dubonos, and A. Firsov, nature ,197 (2005). A. K. Geim and K. S. Novoselov, Nature materials , 183 (2007). A. C. Neto, F. Guinea, N. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics , 109 (2009). A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics , 109 (2009). P. Avouris, Nano letters , 4285 (2010). H.-R. Chang, J. Zhou, H. Zhang, and Y. Yao, Phys. Rev. B , 201411 (2014). G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly, and J. van den Brink, Physical Review B , 073103 (2007). N. Kharche and S. K. Nayak, Nano letters , 5274 (2011). Z. H. Ni, T. Yu, Y. H. Lu, Y. Y. Wang, Y. P. Feng, and Z. X. Shen, ACS nano , 2301 (2008). O. Kibis, Physical Review B , 165433 (2010). X. Li, X. Wang, L. Zhang, S. Lee, and H. Dai, science , 1229 (2008). Y.-W. Son, M. L. Cohen, and S. G. Louie, Physical review letters , 216803 (2006). M. Y. Han, B. ¨Ozyilmaz, Y. Zhang, and P. Kim, Physical review letters , 206805 (2007). M. Ezawa, New Journal of Physics , 033003 (2012). C. J. Tabert and E. J. Nicol, Phys. Rev. B , 195410 (2014). B. Lalmi, H. Oughaddou, H. Enriquez, A. Kara, S. Vizzini, B. Ealet, and B. Aufray, Applied Physics Letters , 223109(2010). C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 226801 (2005). M. Ezawa, Phys. Rev. Lett. , 055502 (2012). N. D. Drummond, V. Z´olyomi, and V. I. Fal’ko, Phys. Rev. B , 075423 (2012). B. Aufray, A. Kara, S. Vizzini, H. Oughaddou, C. Leandri, B. Ealet, and G. Le Lay, Applied Physics Letters , 183102(2010). P. De Padova, C. Quaresima, C. Ottaviani, P. M. Sheverdyaeva, P. Moras, C. Carbone, D. Topwal, B. Olivieri, A. Kara,H. Oughaddou, et al., Applied Physics Letters , 261905 (2010). L. Zhang, P. Bampoulis, A. van Houselt, and H. Zandvliet, Applied Physics Letters , 111605 (2015). A. Acun, L. Zhang, P. Bampoulis, M. Farmanbar, A. van Houselt, A. Rudenko, M. Lingenfelder, G. Brocks, B. Poelsema,M. Katsnelson, et al., Journal of Physics: Condensed Matter , 443002 (2015). L. Li, S.-z. Lu, J. Pan, Z. Qin, Y.-q. Wang, Y. Wang, G.-y. Cao, S. Du, and H.-J. Gao, Advanced Materials , 4820 (2014). M. D´avila, L. Xian, S. Cahangirov, A. Rubio, and G. Le Lay, New Journal of Physics , 095002 (2014). P. Bampoulis, L. Zhang, A. Safaei, R. Van Gastel, B. Poelsema, and H. J. W. Zandvliet, Journal of physics: Condensedmatter , 442001 (2014). M. Derivaz, D. Dentel, R. Stephan, M.-C. Hanf, A. Mehdaoui, P. Sonnet, and C. Pirri, Nano letters , 2510 (2015). F. dAcapito, S. Torrengo, E. Xenogiannopoulou, P. Tsipas, J. M. Velasco, D. Tsoutsou, and A. Dimoulas, Journal ofPhysics: Condensed Matter , 045002 (2016). R. Ganatra and Q. Zhang, ACS nano , 4074 (2014). H. Rostami, A. G. Moghaddam, and R. Asgari, Physical Review B , 085440 (2013). K. Ko´smider, J. W. Gonz´alez, and J. Fern´andez-Rossier, Phys. Rev. B , 245436 (2013). D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. , 196802 (2012). A. Korm´anyos, V. Z´olyomi, N. D. Drummond, P. Rakyta, G. Burkard, and V. I. Fal’ko, Phys. Rev. B , 045416 (2013). H. Rostami, R. Rold´an, E. Cappelluti, R. Asgari, and F. Guinea, Physical Review B , 195402 (2015). B. Radisavljevic, A. Radenovic, J. Brivio, i. V. Giacometti, and A. Kis, Nature nanotechnology , 147 (2011). A. Scholz, T. Stauber, and J. Schliemann, Phys. Rev. B , 035135 (2013). L. Wang, A. Kutana, and B. I. Yakobson, Annalen der Physik (2014). O. Kibis, K. Dini, I. Iorsh, and I. Shelykh, Physical Review B , 125401 (2017). G. Giuliani and G. Vignale,
Quantum theory of the electron liquid (Cambridge university press, 2005). G. Gumbs and D. Huang,
Properties of Interacting Low-Dimensional Systems (John Wiley & Sons, 2013). C. J. Tabert and E. J. Nicol, Physical review letters , 197402 (2013). C. J. Tabert and E. J. Nicol, Physical Review B , 235426 (2013). C. J. Tabert and E. J. Nicol, Physical Review B , 085434 (2013). T. Ando, Journal of the Physical Society of Japan , 074716 (2006). P. Pyatkovskiy, Journal of Physics: Condensed Matter , 025506 (2008). S. D. Sarma, E. Hwang, and H. Min, Physical Review B , 035201 (2015). O. Sobol, P. Pyatkovskiy, E. Gorbar, and V. Gusynin, Physical Review B , 115409 (2016). B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics , 318 (2006). E. H. Hwang and S. Das Sarma, Phys. Rev. B , 205418 (2007). S. Das Sarma and Q. Li, Phys. Rev. B , 235418 (2013). T. Stauber, J. Schliemann, and N. M. R. Peres, Phys. Rev. B , 085409 (2010). A. Scholz, T. Stauber, and J. Schliemann, Phys. Rev. B , 195424 (2012). T. Stauber, G. G´omez-Santos, and L. Brey, arXiv preprint arXiv:1706.09403 (2017). S. Das Sarma and E. H. Hwang, Phys. Rev. Lett. , 206412 (2009). H. Yan, T. Low, W. Zhu, Y. Wu, M. Freitag, X. Li, F. Guinea, P. Avouris, and F. Xia, Nature Photonics , 394 (2013). J. Chen, M. Badioli, P. Alonso-Gonz´alez, S. Thongrattanasiri, F. Huth, J. Osmond, M. Spasenovi´c, A. Centeno, A. Pesquera,P. Godignon, et al., Nature , 77 (2012). X. Luo, T. Qiu, W. Lu, and Z. Ni, Materials Science and Engineering: R: Reports , 351 (2013). A. Grigorenko, M. Polini, and K. Novoselov, Nature photonics , 749 (2012). J. H. Strait, P. Nene, W.-M. Chan, C. Manolatou, S. Tiwari, F. Rana, J. W. Kevek, and P. L. McEuen, Physical ReviewB , 241410 (2013). G. Gumbs, A. Iurov, D. Huang, and W. Pan, Journal of Applied Physics , 054303 (2015). A. Politano and G. Chiarello, Nanoscale , 10927 (2014). A. Politano, V. M. Silkin, I. A. Nechaev, M. S. Vitiello, L. Viti, Z. S. Aliev, M. B. Babanly, G. Chiarello, P. M. Echenique,and E. V. Chulkov, Phys. Rev. Lett. , 216802 (2015). A. Politano, A. R. Marino, V. Formoso, D. Far´ıas, R. Miranda, and G. Chiarello, Phys. Rev. B , 033401 (2011). A. Politano, A. R. Marino, and G. Chiarello, Phys. Rev. B , 085420 (2012). A. Politano and G. Chiarello, Applied Physics Letters , 201608 (2013). A. Politano and G. Chiarello, 2D Materials , 035003 (2017). A. Politano, I. Radovi´c, D. Borka, Z. Miˇskovi´c, H. Yu, D. Farias, and G. Chiarello, Carbon , 70 (2017). S. Das Sarma and Q. Li, Phys. Rev. B , 235418 (2013). D. K. Patel, S. S. Ashraf, and A. C. Sharma, physica status solidi (b) , 1817 (2015). J. Wu, S. Chen, and M. Lin, New Journal of Physics , 125002 (2014). E. H. Hwang and S. Das Sarma, Phys. Rev. B , 165404 (2009). L. Meng, Y. Zhang, W. Yan, L. Feng, L. He, R.-F. Dou, and J.-C. Nie, Applied Physics Letters , 091601 (2012). P. Plochocka, C. Faugeras, M. Orlita, M. L. Sadowski, G. Martinez, M. Potemski, M. O. Goerbig, J.-N. Fuchs, C. Berger,and W. A. De Heer, Physical review letters , 087401 (2008). M. Pisarra, A. Sindona, P. Riccardi, V. Silkin, and J. Pitarke, New Journal of Physics , 083003 (2014). V. Despoja, D. Novko, K. Dekani´c, M. ˇSunji´c, and L. Maruˇsi´c, Physical Review B , 075447 (2013). A. Iurov, G. Gumbs, and D. Huang, Journal of Physics: Condensed Matter , 135602 (2017). C. Walhout, A. Acun, L. Zhang, M. Ezawa, and H. Zandvliet, Journal of physics: Condensed matter , 284006 (2016). A. Iurov, G. Gumbs, D. Huang, and V. Silkin, Physical Review B , 035404 (2016). A. Iurov, G. Gumbs, D. Huang, and L. Zhemchuzhna, Journal of Applied Physics , 084306 (2017). D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Physical Review Letters , 196802 (2012). A. Scholz, T. Stauber, and J. Schliemann, Physical Review B , 035135 (2013). P. F. Maldague, Surface Science , 296 (1978). Y. Gogotsi and V. Presser,
Carbon nanomaterials (CRC press, 2013). H. Raza,
Graphene nanoelectronics: metrology, synthesis, properties and applications (Springer Science & Business Media,2012). W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature , 824 (2003). G. Gumbs, A. Iurov, and N. J. M. Horing, Phys. Rev. B , 235416 (2015). M. Campisi, P. H¨anggi, and P. Talkner, Rev. Mod. Phys. , 771 (2011). G. Schaller,
Open quantum systems far from equilibrium , vol. 881 (Springer, 2014). N. J. Horing, A. Iurov, G. Gumbs, A. Politano, and G. Chiarello, in
Low-Dimensional and Nanostructured Materials andDevices (Springer, 2016), pp. 205–237. N. J. M. Horing, Physical Review B , 193401 (2009). N. J. M. Horing, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and EngineeringSciences , 5525 (2010). N. J. M. Horing, E. Kamen, and H.-L. Cui, Physical Review B , 2184 (1985). G. Gumbs, A. Iurov, J.-Y. Wu, M. Lin, and P. Fekete, Scientific reports (2016). C. Kramberger, R. Hambach, C. Giorgetti, M. H. R¨ummeli, M. Knupfer, J. Fink, B. B¨uchner, L. Reining, E. Einarsson,S. Maruyama, et al., Phys. Rev. Lett. , 196803 (2008). E. d. A. e Silva, G. La Rocca, and F. Bassani, Physical Review B , 16293 (1997). M. Burt, Semiconductor science and technology , 739 (1988). M. Ezawa, Applied Physics Letters , 172103 (2013).
T. Izraeli, Final Project in the Computational Physics Course (2013).
E. W. Weisstein,
CRC concise encyclopedia of mathematics (CRC press, 2002).
I. S. Gradshteyn and I. M. Ryzhik,
Table of integrals, series, and products (Academic press, 2014).
Expression k F = √ πn holds true only for systems with a fourfold spin and valley degeneracy, such as graphene, i.e., g = g s g v = 4. The general equation for the Fermi momentum in a 2D material reads (2 π ) n = g πk F . Due to a large gap parameter ∆ = 1 . eV , the electronic states corresponding to large wave vectors are often involved andat k ≈ . k (0) F correction (cid:119) ∆ β k4