Finite-temperature Coulomb Excitations in Extrinsic Dirac Structures
Andrii Iurov Godfrey Gumbs Danhong Huang Ganesh Balakrishnan
cond-mat.mes-hall





aiurov@unm.edu

Center for High Technology Materials, University of New Mexico,
1313 Goddard SE, Albuquerque, NM, 87106, USA





Department of Physics and Astronomy, Hunter College of the City
University of New York, 695 Park Avenue, New York, NY 10065, USA


Donostia International Physics Center (DIPC),
P de Manuel Lardizabal, 4, 20018 San Sebastian, Basque Country, Spain





Air Force Research Laboratory, Space Vehicles Directorate,
Kirtland Air Force Base, NM 87117, USA


	Center for High Technology Materials, University of New Mexico,
	1313 Goddard SE, Albuquerque, NM, 87106, USA





        Center for High Technology Materials, University of New Mexico,
        1313 Goddard se, Albuquerque, NM, 87106, USA








We have derived algebraic, analytic expressions for the chemical potential without any restriction on temperature
for all types of doped, or extrinsic, gapped Dirac cone materials including gapped graphene, silicene, germanene 
and single-layer transition metal dichalcogenides. As an important intermediate step of our derivations, 
we have established a reliable piecewise-linear model for calculating the density-of-states in molybdenum disulfide, 
showing good agreement with previously obtained numerical results. For the spin- and valley-resolved 
band structure, we obtain an additional decrease of the chemical potential due to thermally induced doping 
of the upper subband at finite temperature. It has been demonstrated that since the symmetry between the electron and hole
states in MoS_2 is broken, the chemical potential could cross the zero-energy level at sufficiently high 
temperature. These results allow us to investigate the collective properties, polarizability, plasmons and their 
damping.  Emphasis is placed on low temperatures, when initial electron doping plays a crucial role.
We clearly demonstrated the contribution of the initial doping to the finite-temperature collective properties 
of the considered materials.


73.21.-b, 73.63.-b, 71.45.Gm, 73.20.Mf



Finite-temperature Coulomb Excitations in
Extrinsic Dirac Structures

    Ganesh Balakrishnan
    December 30, 2023
======================================================================








§ INTRODUCTION



Despite the fact that the microscopic properties of various low-dimensional materials have been 
meticulously examined over a fairly long period of time,  <cit.> only successful fabrication of graphene 
in 2004  <cit.> stimulated an intriguingly new  research effort devoted to the study of atomically thin 
two-dimensional (2D) materials. In particular, it was by virtue of its unique, yet unexpected massless Dirac
electronic properties that led to high mobility (200,000 cm^2/V·s) and ballistic transport properties. 
 <cit.> In the corners of the Brillouin zone, referred to as K and K' points, 
there is no energy band gap and the dispersions represent a linear Dirac cone structure. Due to the
existence of such an energy spectrum, opening a sufficiently large  and tunable energy gap in 
graphene has become an important issue in order to enable electron confinement. Researchers tried to 
achieve this by adjoining  a variety of  insulating 
substrates  <cit.> or even  expos graphene to circularly-polarized radiation.
 <cit.>  In finite-width nanoribbons, their energy band structure and gap are modified by the type
of insulating `cousin' that is introduced.<cit.>




In order to create a truly tunable band gap, one must use a material with large spin-orbit coupling
or a buckled structure. In this regard, silicene, a 2D silicon structure was deemed a good 
candidate. Single-monolayer Si possesses a buckled structure  simply because of the larger ionic size of silicon
compared to carbon. This results in a large spin–orbit band gap 1.55 meV and the possibility to modify
its energy spectrum by applying an external perpendicular electric field  <cit.>. These properties 
make it  display an experimentally realizable Kane–Mele type of quantum spin Hall effect, or a topological insulator state, 
because of the existence of time-reversal symmetry.  <cit.> Unlike graphene, the band structure of 
silicene and its nanoribbons <cit.>  directly depends on spin and valley indices
which give lead to plenty of nanoelectronic, valleytronic and spintronic applications.




Germanene, the most recently discovered and fabricated member of atomically thin buckled 2D honeycomb lattices, 
 <cit.> demonstrates substantially larger
Fermi velocities and a band gap of 20 - 90   meV. Grown by molecular beam epitaxy  <cit.> and 
investigated with x-ray absorption spectroscopy, Ge layers demonstrated satisfactory agreement between
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, and two identical chalcogens C, i.e., sulfur, selenium or tellurium. Schematically, TMDC's are 
described as MC_2. In our consideration, we are mostly going to focus on MoS_2, as
their most studied representative. This material exhibits a semiconductor energy band structure with a very large 
direct gap 1.78 eV, in contrast to its bulk states with indirect gap 1.3 eV, and substantial spin-orbit
coupling.  <cit.> Strictly speaking, MoS_2 is not a Dirac material since the mass terms play a crucial 
role in its energy dispersions, however its low-energy Hamiltonian contains a t_0 a_0  Σ· k term,
corresponding to the linear Dirac cone dispersion.  



An effective two-band continuum model and lattice Hamiltonian, <cit.> based on the tight-binding model,
accounts for the hybridization of the d orbitals in Mo and the p orbitals of sulfur atoms.
It gives an adequate description of its low-energy band structure and predicts large spin splitting. <cit.>
Due to the breaking of inversion symmetry and spin-orbit coupling, spin and valley physics is observed
in all group-IV dichalcogenides, including MoS_2. <cit.>  The low-energy states of such systems
are no longer massive Dirac fermions since there is a a difference between electron and hole masses as well as 
trigonal warping effects.  <cit.> Strain engineering, used to tune optical and electronic properties
of conventional semiconductors, has also been  applied to molybdenum disulfide, and its modified band structure
has been theoretically calculated.  <cit.> These unique electronic properties of a single-layer 
MoS_2  were later used to create high-performance transistors operating at room temperature.  <cit.> 
These electronic models and effective Hamiltonian have been widely used to investigate the collective properties 
of TMDC's  <cit.> and their influence on the gap transition.  <cit.> In optoelectronics, the 
band structure, spin and valley properties of molybdenum disulfide could be successfully controlled by 
off-resonant dressing field.  <cit.>




Current many-body and quantum field theory methods in condensed matter physics <cit.> have provided 
helpful ways to understand the electronic and transport properties of low-dimensional solids, including 
diverse bucked honeycomb materials.  <cit.>  In most 
of these theories, we find the dynamical polarization function, or polarizability, to be the mainstay, 
fundamental quantity, describing the screening of an external potential by interacting electrons.  <cit.>    
Also, the dynamic polarization function plays a key role in calculating the plasmon excitations, due to 
the charge  density oscillations, which occur in metals and semiconductors. Specifically, 
the plasmon dispersion relation  along with their lifetimes  have been theoretically investigated for a wide range of 2D 
Dirac systems.  <cit.> 
The interest in graphene plasmons is due in part to the fact that  these excitations
have no classical counterpart.  <cit.> There has also been a considerable experimental effort 
for investigating graphene plasmons,  gate-tuning, infrared nano-imaging and 
confinement  <cit.>  
Graphene plasmonic resonances and instability at various wavelengths could be used in photodetectors 
in the Terahertz range.  <cit.> All these techniques could be successfully applied to the recently fabricated
materials, discussed in the present work. 
 

 Plasmonic applications have been widely based on nanoscale hybrid systems, in which graphene 
plasmons are coupled to a surface plasmon excitations in metals. Technology has now gone a long way in combining 
graphene with prefabricated plasmonic nanoarrays and metamaterials in order to produce plasmonics-based tunable
hybrid optical devices.  <cit.> Therefore, accurate knowledge of plasmon mode dispersions in graphene interfaced 
with metallic substrates is crucial. Graphene-metal contacts are important components for all such devices. 
Consequently, exploration of plasmon modes at these metallic interfaces is a mandatory step toward fabricating the devices. 
High-resolution electron energy loss spectroscopy (EELS)  has been employed to investigate those excitations at the surface of 
Bi_2 Se_3 to disclose the interplay between surface and Dirac plasmons in topological insulators  <cit.>. Plasmons, 
their behavior, dispersions, quenching and environmental effects have been thoroughly studied in epitaxial
graphene, in air-exposed graphene-Ru contacts, Graphene on Pt_3Ni (1 1 1), graphene grown on Cu (1 1 1) 
foils.  <cit.>



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 the plasmon excitations and the electrical conductivity are completely suppressed 
at zero temperature due to the absence of free charge carriers. However, at  finite temperature, the
conduction band would receive thermally induced doping in both cases.  <cit.> In graphene, 
with zero energy band gap, this density is enhanced as n ⋍ T^2 and the plasmon dispersion 
behaves like Ω_p^2 ⋍ q T.


The properties of intrinsic finite-temperature plasmon excitations have been systematically examined 
for various materials including silicene.  <cit.> In contrast, extrinsic or 
doped structure at finite temperature is associated with a difficulty to obtain a reliable 
and accurate  value for temperature-dependent chemical potential μ(T). 
Generally, it is known that μ(T) is decreased as the temperature is increased, and its value could be found 
based on carrier density conservation  <cit.> In this work, our main objective is to obtain a set 
of non-integral, trancendental equations for a wide class of Dirac gapped materials with linear 
density-of-states (DOS), i.e., gapped graphene, silicene, germanene and transition metal dichalcogenides at arbitrary temperature.

However, the range of our considered temperatures is limited by validity of linear or gapped Dirac cone 
approximation for the energies, which recieve noticeable doping at those temperatures. Certain deviations
start to build up at about 0.5   eV, <cit.> leading to various effect on the plasmons, such 
as anisotropy, splitting and existence of additional acoustic plasmon branch.  <cit.> Such energies
are extremely far away from our range k_B T ∽ E_F^(0). In all our calculations, the energy is measured in
the units of a typical Fermi energy  E_F^(0) = 5.22  meV, corresponding to electron density n^(0) = 1.0 · 10^15 m^-2.

Once the chemical potential is known, one can obtain the finite-temperature dynamical polarization 
function by Eq. (<ref>), 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.  <cit.> Here, we pay close attention to 
finite-temperature plasmons, demonstrating how much initial doping contributes to each branch location at intermediate
temperatures. Once the temperature becomes very high, k_B T ≫ E_F, thermally-induced doping dominates 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.  
<ref>. In Sec. <ref>, we calculate the dynamic polarization function which includes the single
particle excitation mode frequencies. The single-particle modes combine with a charge cloud to 
produce weakly interacting quasiparticles that vibrate collectively at the characteristic plasma frequencies.   
Emphasis has been placed on rather simple cases of gapped graphene and silicene at small, but finite 
temperatures, when the zero temperature carrier doping plays a crucial role. We also briefly examine 
non-local, hybrid plasmons in an open system of a semi-infinite conductor, Coulomb-coupled to a 
2D layer in the presence of finite doping doping of the electrons in the layer. Our concluding 
remarks and a concise discussion are presented in Sec. <ref>. We also provide two appendices with detailed 
derivations of the DOS for silicene and MoS_2 - in Appendix <ref>, and of the 
temperature-dependent chemical potential in Appendix <ref>.







§ 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 of the 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 Dirac system with symmetry between
the electron and hole states. Here, we are going to provide closed-form trancendental equations for the 
finite-temperature μ(T) for a number of Dirac systems: graphene, buckled honeycomb lattices and transition 
metal dischacogenides.
  
 



 §.§ Buckled honeycomb lattices


One of the most outstanding features of silicene and other buckled honeycomb lattices is the existence of two generally
double degenerate pairs of energy subbands and two inequivalent band gaps. These gaps There is a fixed intrinsic 
spin-orbit gap 2 Δ_SO and a tunable sublattice-asymmetry gap Δ_z  which is induced by and proportional 
to an applied perpendicular electric field E_z. For small fields, Δ_z = E_z   d_, where 
d_ is the out-of-plane displacement of a buckled lattice.


The low-energy model Hamiltonian of a buckled honeycomb lattice has been found to be  <cit.> 


    Ĥ_ξ,σ = ħ v_F ( 
     ξ k_x τ̂_x + k_y τ̂_y
     ) ⊗Î_2 × 2 - ξΔ_SOΣ̂_z ⊗τ̂_z +
    Δ_z τ̂_z ⊗Î_2 × 2  ,

where the Fermi velocity v_F = 0.5 · 10^6 m/s is half that for graphene, ξ =  ± 1
is the K/K' 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. (<ref>) in a block-diagonal matrix form


    Ĥ_ξ,σ = ( [      - ξσΔ_SO + Δ_z ħ v_F (k_x - i k_y); ħ v_F (k_x + i k_y)        ξσΔ_SO - Δ_z ])
                                     .

The energy dispersions are


    ε_ξ,σ^γ(k) = γ√((
     ξσΔ_z - Δ_SO)^2 +
     (ħ v_F k )^2 
     )  ,

where γ = ± 1 determines the electron or hole state similar to graphene, with a finite or zero band gap.
These dispersions, given by Eq. (<ref>), represent two pairs of spin-dependent energy subbands in
a chosen valley with the two generally different band gaps |Δ_SO - ξσΔ_z |,
which will be later referred to as Δ_<,> = |Δ_SO∓Δ_z |. 
Clearly, both energy gaps depend on the perpendicular electrostatic field and the two subbands corresponding to the
ξσ = ± 1 indices switch their locations with increasing field strength. Small or zero 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 a finite Δ_>, defined as valley-spin polarized metal (VSPM).
For larger fields, Δ_z ∽E_z would always exceed the constant intrinsic spin-orbit gap Δ_SO, 
which corresponds to the standard band insulator (BI) state. 

 
 








The DOS which is in general defined by


    ρ_d(E)= ∫d^2 k/(2 π)^2∑_γ = 
    ± 1 ∑_ξ,σ= ± 1 δ[ E -
    ε_ξ,σ^γ(k) ]   ,

is immediately obtained for silicene (see Appendix <ref>) as


    ρ_d(E) = 1/π ∑_γ = ± 1E/ħ^2 v_F^2∑_i = <,>Θ[ 
    E/γ - Δ_i
     ]   ,

in terms of the unit step function Θ(x). 
We note that for systems sharing the same Dirac cone characteristics with arbitrary energy gap, 
the DOS is linear analogous  to  graphene. Experimentally obtained linear V-shaped DOS was used 
to verify the Dirac cone dispersion for germanene. <cit.> However, ρ_d (E) has a 
finite value only above the energy gap since 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) + (-1) n^(h) = ∫_0^∞ d E ρ_d(E) f_γ = 1(E,T) - 
     ∫_-∞^0 d E ρ_d( E )  { 1 - f_γ = 1(E,T) }  .


At zero temperature, the density n is related to the Fermi energy E_F in a straightforward way.
If only one subband is occupied, we obtain

    n = 1/2 π E_F^2 - Δ_<^2/ħ^2 v_F^2  .

Alternatively, if the doping density is such that both subbands are populate, then E_F is obtained from


    n = 1/π 1/ħ^2 v_F^2 [ E_F^2 - 1/2 ( Δ_<^2 + Δ_>^2) ]   ,


and the critical density required to start filling the upper subband is n_c = 2 Δ_SOΔ_z / πħ^2 v_F^2.



We prove in Appendix <ref> that  the finite-temperature chemical potentail could be obtained from the following
equation 


    n  ( ħ v_F/k_B T)^2 =  ∑_γ = ± 1 γ/π ∑_i = <,>
    - Li_ 2{ - exp[γμ(T) - Δ_i/k_B T] }  +
    Δ_i/k_B T ln{ 1 + exp[γμ(T) - Δ_i/k_B T] 
     }  ,

where Li_ 2 (x) is a polylogarithm function. Connecting the doping density n with the Fermi energy through either 
Eq. (<ref>) or Eq. (<ref>) depending on how many subbands are doped at zero temperature, we derive the 
chemical potential with a value equal to E_F at T = 0, for all accessible  temperatures.
Although this equation is transcendental and cannot be resolved algebraically, a quasi-analytic or 
one-step numerical solutions could be easily provided for any finite temperature without having to perform an integration. 



Clearly, the chemical potential for silicene depends on the energy band gaps Δ_i, i = < ,>. Our approach, discussed 
in Appendix <ref>, is valid for a variety of materials with linear energy dependence for the DOS, including MoS_2.  
Specifically, Eq. (<ref>) also describes μ(T) for gapped graphene with two degenerate subbands, or Δ_< = 
Δ_> = Δ_0. For gapless pristine graphene Δ_0 = 0 and π n = [ E_F/ (ħ v_F) ]^2, we have 


    1/2 ( k_B T )^2  E_F^2  
     = - ∑_γ = ± 1γ Li_ 2{ - exp[γ μ(T)/k_B T] }  .

If the temperature is low with k_B T ≪ E_F, Eq. (<ref>) is reduced to the expressions,  
derived in Refs. [SDSS, SDSLi] and [myT].

  






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 
±√((ħ v_F k)^2 + Δ_0^2) states. Silicene and germanene dispersions, yet showing complete electron/hole
symmetry, exhibit spin- and valley-dependent pairs of subbands, each being double degenerate. Finally, MoS_2 demonstrates
broken symmetry between its 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. <ref>.
For both materials, 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 immediate increase.
or discontinuity, of the DOS, as schematically shown for silicene in the insets of Fig. <ref> (b). It is important 
to observe that for both silicene and gapped graphene, ρ_d (E) is directly proportional to energy E, i.e.,
the DOS for Dirac 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. 






In Fig. <ref>, we display our results for the finite-temperature chemical potential for graphene and buckled 
honeycomb lattices. First, we show how dissimilar this temperature dependence might be for various electronic 
systems. In Fig. <ref> (a), this situation is described for a 2DEG with a parabolic energy band and constant
DOS, graphene with ρ_s(E) ∽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 nearly identical since the 
holes do net play any important role. The hole distribution function is complimentary to that for electrons, i.e., 
1-f_e(ϵ) → 0, and, therefore, inconsequential. 
However, when the temperature becomes comparable with the Fermi energy, the hole thermal excitations become crucial, 
causing an opposite effect compared with that for electrons. They mitigate the reduction of the chemical potential 
and eventually prevent μ(T)  from  crossing the zero energy level. This is seen particularly well for a 
hole-dominating system (> h^(+)), for which the chemical potential starts 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 larger 
gap 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. (<ref>)). The former case is shown in panels (b), (c) and (d),
whereas the latter at (e) and (f). In general, the carrier density n  is accepted 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 between electron 
and hole doping for all energy band gaps. This is a manifestation of γ = ± 1 symmetry of these states persisting
at arbitrary temperatures. Silicene with both filled subbands exhibits qualitatively similar finite-temperature behavior as 
gapped graphene. More generally, only the actual gaps are relevant, but not the exact type of states (such as topological 
insulator or regular band insulator), or the way in which such a state has been achieved.  



We now pay attention to the following situation in silicene. At T=0, the Fermi level is chosen so that only the lower
subband is filled. In contrast, the upper one Δ_> is located so close  to the Fermi level E_F, that it starts 
getting populated at all, even very small temperatures. Thermally-induced doping, received by the Δ_>-subband, results
in a much stronger reduction rate of μ(T), as shown in Fig. <ref> (e). Such situation suggests very special 
thermal properties of silicene with closely located energy subbands, and shares this behavior with graphene having an additional spin and
valley degeneracy. We investigate this phenomenon even further by adjusting the upper subband in the vicinity of the Fermi
level, as shown in Fig. <ref> (f). Now, each μ(T) curve demonstrates a significantly pronounced decrease whenever the 
upper subband filling becomes essential. These two different types of μ(T) behavior in silicene could be used to achieve 
its additional tunability by introducing the upper energy subband Δ_>. 






 §.§ Transition metal dichalcogenides


The low-energy electronic states in monolayer molybdenum disulfide (MoS_2 , ML-MDS), a prototype transition metal 
dichalcogenide, could be effectively described by a two-band model Hamiltonian  <cit.>


    Ĥ_d^ξ,σ = ( 
    1/2 ξσ λ_0 + ħ^2 k^2/4 m_eα) Î_2× 2
    +
    (
    Δ/2 - 1/2 ξσ λ_0 + ħ^2 k^2/4 m_eβ) Σ̂_z + 
    t_0 a_0 Σ̂_ξ· k  ,

whose important feature is a major  gap parameter Δ = 1.9 eV which results in the actual 
band gap ⋍ 1.7   eV,  large compared to that for silicene. The spin-orbit coupling 
parameter λ_0=0.042 Δ represents a  smaller,
but essential correction, to the single-particle excitation spectrum and the band gap. The energy subbands are now spin- and valley-
dependent since the corresponding degeneracy is lifted. The electron hopping parameter t_0 = 0.884 Δ and a_0=1.843Å shape 
the Dirac cone term in the Hamiltonian Eq .(<ref>) as t_0 a_0= 4.95× 10^-29J· m, 
counting up to ≈ 0.47  of ħ v_F value in graphene. 


Next, we include  the ∽ k^2 mass terms with α =2.21 = 5.140 β in which m_e is the
 free electron mass. Our considered values of Fermi momentum at zero temperature are determined 
by the experimentally allowed electron and hole doping densities n = 
10^14÷ 10^16 m^-2 as k_F = √(π n)⋍ 10^8÷ 10^9 m^-1.
 [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 π)^2   n = g  π k_F.]
Anisotropic trigonal warping term t_1 a_0^2 (Σ̂_ξ· k) σ̂_x (Σ̂_ξ· k) is clearly 
beyond out consideration since 
∽ t_1 = 0.1 eV = 0.053 Δ term does not cause any effect on the electron dispersions.   
The energy dispersion relations, corresponding to the Hamiltonian 
Eq. (<ref>) ε^ξ,σ_γ(k) = ϵ_0^ξ,
σ(k) + γ√([ Δ_0^ξ,σ(k) ]^2 + (t_0 a_0 k)^2) , 
formally represent   gapped graphene with 
a k-dependent gap term  Δ_0^ξ, σ(k) = ħ^2 k^2 β / (4 m_e)
 + Δ / 2 - ξσ λ_0 / 2, and a band
shift ϵ_0^ξ,σ(k) = ħ^2 k^2 α /(4 m_e) + ξσ λ_0 / 2.  <cit.>





Neglecting only the O(k^4) terms leads us to 


    ε^ξ,σ_γ(k)⋍1/2 ξσλ_0 + αħ^2/4 m_e k^2 + 
    γ/2{[
    (2 t_0 a_0)^2+( Δ - ξσ λ_0  )  βħ^2 / m_e
    ] k^2 +
    ( Δ - ξσ λ_0 )^2
    }^1/2  .

This is the principal model which we will use to describe the energy dispersions of MoS_2 in our work. 
The ⋍ k^4 terms, trigonal warping and anisotropy are considered to be non-essential, even though
 as we will later see, cause certain discrepancies in the DOS. 
[Due to a large  gap parameter Δ = 1.9 eV, the electronic states corresponding to large 
wave vectors are often involved and at k ≈ 5.0 k_F^(0) correction 
 ⋍Δ β   k^4 is no longer negligible.]
We also show (see Appendix <ref>) that the curvature of the energy subbands in TMDC's is so small that
 even the highest  possible doping density 10^17 m^-2 results in the Fermi energies 
∽λ_0. Thus, at zero or low temperatures, we do not need to consider any high-energy 
corrections to Eq. (<ref>)  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.  <cit.> In
that case, our model Hamiltonian in Eq. (<ref>) and especially, simplified dispersion 
relation (<ref>) would no longer provide a satisfactory   approximation. Consequently,
our primary focus is on small but finite temperatures for which the initial doping density and E_F still play an 
important role, beyond the O(T^2/T_F^2 ) approximation discussed in Ref. [SDSLi]. 



Here we note that for MoS_2, similar to the buckled honeycomb lattices, spin and valley indices always 
appear together as a product, so that taking into account a × 2 degeneracy, a single composite index
ν = ξσ could be effectively 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 (<ref>) arises 
once we also neglect the  mass terms in Eq. (<ref>)  <cit.>


    ε^ν_γ(k)⋍νλ_0/2+γ√((t_0 a_0)^2 k^2+(Δ-νλ_0)^2/4)  .

This approximation has a few important advantages  including simplicity and its formal resemblance 
with gapped graphene so that all the crucial quantities such as the DOS, wave function,  polarizability 
and many others are already known. Furthermore, it gives a quite a suitable description of the energy 
band structure of MoS_2, taking into account a large gap parameter and ν-dependent splitting of the two hole subbands.  
Nevertheless, the mass terms must be taken into account for a proper evaluation of the DOS and most of the 
temperature-dependent properties of TMDC's. As we demonstrated in Appendix <ref>, even in the simplest possible parabolic
subbands approximation, valid only for small wave vectors k → 0, the mass terms make a contribution,
comparable with the Dirac cone and band gap parts of Hamiltonian (<ref>). Very small curvature of the energy subbands, 
mentioned above results in a tremendous ρ_d ∽ m_eff, so that even ⋍α/(4 m_e) correction, hardly
noticeable on the electron band structure,   becomes significant. It basically discards the DOS results obtained 
from the gapped graphene model, even at zero or small wave vectors. Importance of the mass and even higher order terms for the
plasmon calculation was discussed in Ref. [Scholz1].  



Taking into account all the required terms in (<ref>),   rigorous numerical calculations give   
accurate results for the electron DOS for MoS_2. In Fig. <ref> (d), we present all three possible
outcomes.  Based on the gapped graphene model, the DOS is nearly twice as large as its numerical values
.  <cit.> In contrast, our ∽ k^2 model (<ref>) demonstrates quite a good match,
 especially in the low-energy range.   We also note that the numerically obtained dependence is clearly linear 
for a wide range of energies, much exceeding our considered diapason.


In summary, we consider a piecewise linear approximation 
ρ_d (E) relatively close (δE⋍λ_0 ) for each of
the three non-degenerate subbands ρ_d(E) = c_0^(i) + c_1^(i) E and ρ_d(E) = 
0 in the gap region -Δ/2 + λ_0 < E < Δ/2. The expansion coefficients 
are obtained as c_0^(1) =  2.837   E_F^(0)/(ħ v_F)^2 = 0.043  t_0^-1 a_0^-2, c_1^(1) = - 
1.397   (ħ v_F)^-2 = -0.308  (t_0 a_0)^-2 
for γ = -1, ν = -1 and ≤ -Δ/2 - λ_0 ≥E≤  -Δ/2 - λ_0.
Finally, when γ = -1, but ν = 1 and E≤ -Δ/2 + λ_0,
the hole DOS coefficients are c_0^(2) = 1.132   E_F^(0)/(ħ v_F)^2 = 0.0174  t_0^-1 a_0^-2
c_1^(2) = -0.767   (ħ v_F)^-2 = -0.169  (t_0 a_0)^-2. For electrons with γ=+1 at 
E= Δ/2 + δϵ, the two quasi-degenerate subbands lead to the DOS 
equal to c_0^(3) = 5.110   E_F^(0)/(ħ v_F)^2 = 0.078   t_0^-1 a_0^-2 and c_1^(3) = 
0.815   (ħ v_F)^-2 = 0.179  (t_0 a_0)^-2. 


We adopt these values for ρ_d(E), which arise from the numerical calculations in order 
to achieve the highest possible precision and credibility for our finite-temperature derivations. However, 
 our effective model, presented in Appendix <ref>, gives the DOS results which show good agreement with 
these numerical values and could be used for decisive estimates for various collective calculations for MoS_2. 

 





Once the DOS is known, we are in a  position to calculate the Fermi energy for a given doping density n_(0)^e.  
The new point here is that ρ_d (E) is not directly proportional to the energy so that E_F is determined by


    n_(0)^e = (
     E_F - Δ/2)  [ c_0^(3) +
    c_1^(3)/2 (
    Δ/2 + E_F
    )
    ]   ,

or 


    E_F^(e) = 1/2   c_1^(3) {-2  c_0^(3) + [ ( 2   c_0^(3) + c_1^(3) Δ)^2 + 8 n_(0)^e 
     c_1^(3)]^1/2}  .

For hole doping, the Fermi energy differs from the previously considered electron doping case, i.e., 


    E_F^(h) =  1/c_1^(2) { 
     - c_0^(2) + { 
     -2 c_1^(2)  n_(0)^h + [
     c_0^(2) - c_1^(2) (
     Δ/2 - λ_0
     )
     ]^2 
     }^1/2 }

Numerically, our results for the  Fermi energy for electrons and holes are presented in insets 
(i1) and (i2)  of Fig. <ref>. Here, both linear and quadratic terms in the doping density 
equations are present (see Eq. (<ref>)), and, most importantly, there is no symmetry
between the electron and hole states. Unlike graphene, the linear dependence here dominates for both 
cases due to the large energy band gap. Each curve starts from the corresponding band 
gap - Δ/2 = 8.13 E_F^(0) for electrons and -Δ/2 + λ_0 = - 7.45 E_F^0. 
The corresponding well-known result for gapped graphene E_F^2 - Δ_0^2 = 
2 π  n  (ħ v_F )^2 are verified by putting c_1^(i)→ 0 and λ_0 → 0. 


The finite-temperature chemical potential for MoS_2 is obtained in a similar way, as we have 
done for the buckled honeycomb lattices, except that we need to evaluate four different 
terms related to the two separate hole subbands (see Eq. (<ref>)).  The corresponding 
numerical results are described in Fig. <ref>. As discussed above, the most special 
property of TMDC's is the broken electron/hole symmetry. Consequently, the chemical potential for 
the electron doping becomes negative at T ≈ 2.5 E_F, while μ(T) for hole  
doping does not ever change its sign. Broken electron/hole  symmetry leads to two substantially 
different types of temperature dependence of the chemical potential. Thus,   MoS_2 
represents a special material with unique symmetry properties and chemical potential dependence, 
so far encountered only in certain types of semiconductors but not in Dirac materials. 








§ PLASMON EXCITATIONS AT FINITE TEMPERATURE



As one of the most relevant applications of our finite-temperature chemical potential formalism, 
we briefly consider plasmons for an extrinsic, substantially doped at T=0, free-standing gapped Dirac cone material.
The plasmon dispersion relation is calculated from the zeros of the dielectric function ϵ(q, ω), 
expressed in the random phase approximation (RPA) as


    ϵ(q, ω) = 1 - v(q)  Π_T(q, ω | μ(T), T, Δ_i) = 0   ,

where v(q) = 2 π e^2 / (ϵ_s q) is the Fourier-transformed two-dimensional Coulomb potential, 
and ϵ_s = 4 πϵ_0 ϵ_r ith  ϵ_b is the background dielectric constant  
in which the 2D material is embedded.
At finite temperature, the dynamical polarization function Π_T is given as an integral transformation <cit.>
of its T=0 counterpart Π_0, i.e., 


    Π_T(q, ω | μ(T), T, Δ_i) = 1/2 k_B T ∫_0^∞  d ξ Π_0 (q,ω | ξ, Δ_i)/1 + cosh{ [ μ(T,E_F) - ξ] / (k_B T) }

The evaluation of the zero-temperature polarizability is quite similar for buckled honeycomb lattices and MoS_2,
since in both cases their low-energy bandstructure is represented by two generally inequivalent double-degenerate
pairs of subbands which depend on the composite index ν = σ ξ. For any such pair,   Π_0^(ν) is
obtained in the one-loop approximation (g_0 = 2) as


    Π_0^(ν) (q,ω |  E_F, Δ_ν) = g_0/8 π^2∫ d^2 k ∑_γ,γ' = ± 1( 
      1 + γγ'   k · ( k +  q) + Δ_ν^2 /ε^ν (k)   ε^ν(| k +  q|) )
      f[ E - ε^ν_γ (k)] - f[E - ε^ν_γ'(| k +  q|)]/ε^ν_γ (k) - ε^ν_γ' (| k +  q|)  ,

where f[E - γ E_ν(k)] is the Fermi-Dirac distribution function, showing electron and hole 
occupation numbers for chosen energy E. At T=0, it is a Heaviside unit step function 
Θ[ E - ε^ν_γ
(k) ]. We note that the extra 1/2   comes from the form factor 1/2   [ 1 + γγ' ... 
 which also depends on the absolute value of each energy 
ε^ν (k) = 1 / γ  ε^ν_γ (k).  We also note 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 γ'.
Thus, for both types of materials, the dynamical polarization function is obtained as a sum of terms obtained from 
Eq.  (<ref>) over ν   <cit.> 


    Π_0 (q,ω |  E_F, Δ_i) = ∑_ν Π_0^(ν) (q,ω |  E_F, Δ_ν)   .

Our results for the polarization functions and plasmon excitations are presented in Figs. <ref> and <ref>. 

   



  



First, we need to address the imaginary part of the polarizability  since it specifies the regions and intensity
of the plasmon damping. We see that at a finite temperature the plasmon dissipation generally increases and the 
damping-free regions are nearly absent. Once the temperature becomes   high, the imaginary part of the polarization
function is reduced as 1/T,  <cit.> so there is no uniform temperature dependence of the plasmon damping. 
Doping and proportional increase of the energy band gaps, in contrast, increase the regions free of single-particle 
excitation spectrum (see Fig. <ref> (b)), so the effect of a finite temperature and E_F on the plasmon dissipation
coul be opposite.
The real part of Π_0 (q,ω |  E_F, Δ_i), shown in Fig. <ref> (a)-(d), designates the 
location and slope of the corresponding plasmon branches which are presented in panels (e) and (f). Here,
the temperature and the doping produce a similar effects, shifting the location of the peaks to the right.
The real part of the polarization function must be positive in order to enable a real solution of Eq. (<ref>) 
Plasmon branches are located at higher energies for a given wave vector due to either an initial increase of the Fermi 
energy or thermally-induced doping. 



 





 §.§ Non-local plasmons in an open system




Concluding our investigation of extrinsic 2D materials, we now turn to the     
plasmon excitations in so-called 2D open systems (2DMOS).
Such a nanoscale hybrid arrangement being a part of graphene-based nanoscale devices, <cit.> consists
of a 2D layer (it could be 2DEG, graphene, a buckled honeycomb layer, or MoS_2), which is Coulomb coupled
(not chemically bound) with a semi-infinite conductor and its surface plasmon.  <cit.>
While a plasmon excitation in a closed system is determined by two-particle Green’s functions, 
in 2DMOS it more involved, depending on the  Coulomb interaction with the environment.  <cit.>
The important feature of such a system is the screened Coulomb coupling between the electrons in graphene and 
the conducting substrate.<cit.> Such a screened potential could be obtained using the nonlocal frequency-dependent 
inverse dielectric function  <cit.>  Consequently, two <cit.> or 
more <cit.> linear plasmon branched have been obtained, confirming some previous experimental  
claims.  <cit.>   For finite temperature, their coupling to an external reservoir is
 reflected in existence of extra  plasmon dissipation channels. <cit.>


The plasmon branches in such system are obtained as zeros of a so-called structure factor 
S (q,ω | μ(T), Δ_i ), playing the role of the dielectric function 
ϵ (q, ω) in an isolated layer. It is obtained as  <cit.>


    S (q,ω | μ(T), Δ_i) = 1 - v(q)  Π^(0)(q,ω;μ) {1 - ϵ_B(ω)/1 + ϵ_B(ω) Exp(-2 a   q) + 1
     }  ,

where a is the distance between the 2D layer and the surface and the bulk dielectric function
is given in the local limit as ϵ_B = 1- Ω_p^2/ω^2. The bulk-plasma frequency,
defined as Ω_p^2 = (n_m e^2)/(ϵ_o ϵ_S m^∗), depends on on the electron concentration
n_m, its effective mass m^∗ and the substrate dielectric constant ϵ_S. This approximation 
stays valid for a large range of wave vectors q ≪ 2 × 10^9 m^-1 since the Fermi 
wavelength in metals is comparable with the inverse lattice constant. As a result the frequency of the 
upper plasmon branch in our system, equal to Ω_p/√(2) at q_∥→ 0, might range from 
ultra violet down to infrared or even terahertz, depending on the substate material and must stay commensurate
with the energy band gap in the 2D layer. 

  

 
 
 


Previously, we reported that  <cit.> in case of spin- and valey-dependent single-particle excitations
in a 2D layer (which is true for buckled honeyomb lattices and MoS_2 ), such a hybrid structure
could be effectively used to directly measure the dielectric properties or spin-orbit coupling parameters 
of such a layered material because the location of each plasmon branch, its damping rate and the signatures of the 
particle-hole modes are independently determined by the material parameters of the 2D layer. 
Each of these properties is unique (for example, the two plasmon branches in silicene-based TDMOS depends on 
the energy band gaps as Δ_i^1/2 and Δ_i^1/4, while the outermost PHM  boundaries are 
determined solely by the lower gap Δ_<), so that an additional linear plasmon branch provides us 
with the required earlier unknown piece of information about the specific material. This is not possible 
in the case of a single plasmon branch in an isolated layer.


In the present work, we additionally introduce finite doping and temperature into this consideration. Our numerical
results for the non-local plasmons for a extrinsic systems are presented in Fig. <ref> We see that the location 
and damping of both branches substantially changes in the presence of  an initial Fermi energy. At low temperature, this
increase is especially apparent, similar to plasmons in an isolated layer. For a sample with stronger doping, 
E_F = 1.5 E_F^(0) at zero temperature, the upper plasmon branch is always located above the surface plasmon level
Ω_p/√(2) and is never damped. 




§ CONCLUDING REMARKS




We have carried out an extensive investigation of extrinsic, or doped, Dirac gapped  materials
at arbitrary finite temperatures and obtained a set of  algebraic analytic equations determining the chemical 
potential. Our considered systems 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 broken symmetry between the electron and hole states. Our results could also be used 
to predict the finite-temperature chemical potential for parabolic or quasi-parabolic eigenstates in 
semiconductors with light or heavy holes.  <cit.> In general, our model is limited
only by the linear dependence of the DOS which stays valid over a wide energy range for all the 
above mentioned materials.  




We have demonstrated that the chemical potential depends substantially on the energy band gap(s) of the considered 
system since the DOS depends on the curvature of each subband. Specifically, we investigated structures
with two non-degenerate, separated spin- and valley- resolved energy subbands in both valence and conduction bands, 
such as silicene. The upper subband would receive thermally-induced doping even if it is undoped at zero temperature. This is
always reflected in a higher reduction rate of the chemical potential whenever the second subband doping starts playing 
a role. Consequently, one can tune the μ(T) dependence around any required temperature by bringing initial doping close
to the higher-energy subband. The number of such separated subbands contributing to the DOS could be arbitrarily large for 
an electron in a quantum well or quasi-one-dimensional nanoribbons,  <cit.> and all these cases could be effectively
treated by our model. 


The behavior of the chemical potential  depends on whether there is symmetry between the electron and hole states in the system. If 
the DOS for the electrons and holes is equal or symmetric around the Dirac point, then the chemical potentail does not 
change its sign even for high temperatures, i.e., it remains positive for electron doping or stays negative if E_F < 0. 
In fact, these two types of doping result in symmetric behavior, decreasing |μ(T) |, so that the 
electron-hole symmetry persists at any finite temperature. Let us discuss the electron doping in more detail. Once the 
temperature is sufficiently high, the hole states become thermally excited and it has an opposite effect to that of electron,
decreasing the reduction of the chemical potential. At extremely high temperatures, the two processes have almost equivalent
effects and μ(T) asymptotically tends to zero and never crosses the Dirac point.
This situation changes if the DOS of the electron and holes differs and the two hole subbands are energetically 
not equivalent, as we observe for TMDC's. We have shown that for MoS_2 the chemical potential becomes negative, changing its sign
at T ⋍ 2.5 E_F. Alternatively, μ(T) could never reach the zero energy line if there is hole doping
at zero temperature and starts decreasing at sufficiently high temperatures. 


As a necessary intermediate step in our derivations, we obtained a piecewise-linear model for the DOS 
for transition metal dichalcogenides, directly from the Hamiltonian parameters of the considered system. 
This model gives exact results for ρ_d(E) at the band edges (next to each energy gap) and a 
fairly good approximation at higher energies. This model 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 
(ħΩ_p/ E_F)^2  ⋍Λ  q, and the 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) plasmon dispersions even for
intrinsic systems.  <cit.> We demonstrate how much each type of doping contributes this effective length.     


Several many-body calculations or collective electronic models requires reliable knowledge of 
the chemical potential at finite temperature. In the absence of such information, much attention has
been directed towards intrinsic, or undoped systems, or low temperatures, so that ⋍ (T/T_F)^2 
series expansions could have been applied. Our results are going to provide considerable assistance   in
transport studies, optical, thermally modulated conductivity for the  materials, 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.    


D.H. would like to thank the support from the Air Force 
Office of Scientific Research (AFOSR).
















§ DENSITY-OF-STATES



The DOS for electrons and holes with energy dispersion ε^ ξ,σ_γ(k) 
is defined as


    ρ_d(E)= ∫d^2 k/(2 π)^2∑_γ = ± 1 ∑_ξ,σ = ± 1 δ[ E -ε^ ξ,σ_γ(k) ]   ,

where ξ = ± 1 and σ = ± 1 are valley and spin indices which in our considerations
always appear as a product σξ, so that a single composite index ν = σξ = ± 1 
could be introduced. This leads to a double degeneracy of all considered dispersion relations and 
a transformation  ∑_ξ,σ = ± 1⟹ 2
·∑_ν = ± 1  , which we will use throughout this work.



For a number of cases of 2D structures, Eq.(<ref>) could be immediately evaluated 
using the following property of a delta function

    δ(f(x)) = ∑_iδ(x-x_i)/| df(x)/dx_ x=x_i|  ,

where x_i are the roots of f(x), formally given as f(x) = E - ϵ_ν^γ(k) for various 
energy dispersions ϵ_ν^γ(k). In the case of silicene and germanene, the result is 
straightforward with


    ρ_d(E) = 1/π ∑_γ = ± 1γ E/ħ^2 v_F^2∑_i = <,>Θ[ E/γ - Δ_i
     ]   ,

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 Δ_< = Δ_> = Δ_0. Furthermore, we arrive at well-known V-shaped 
∽E DOS for pristine gapless graphene if Δ_<,>→ 0.  







 §.§ Molybdenum disulfide


In our work, we discuss several effetive models of different complexity and accuracy, describing the 
energy dispersion ε^ν_γ(k) for MoS_2. 
First, we consider a spin- and valley- resolved gapped graphene model given by Eq. (<ref>), in which
we leave out all the mass ∽α and ∽β terms. It gives quite accurate results 
for the energy eigenstates next to the corners of Brillouin zone, as shown in Fig. <ref>.  
However, it is straightforward to see that the corresponding DOS obtained as 


    ρ_d(E) = 1/π  (t_0 a_0)^2∑_γ = ± 1 1/γ ∑_ν= ± 1  (
     E - ν/2λ_0 
     )  Θ[
    γ(E -νλ_0/2) - 1/2( Δ - νλ_0 )
    ]   ,
     
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,1_γ = 1(k = 0) = Δ/2 for electrons,
ε^ν=-1_γ = -1(k = 0) = -Δ/2 + λ_0 and ε^ν=1_γ 
= -1(k = 0) = -Δ/2 - λ_0 for the holes, the DOS experiences three different
giant discontinuities due to each new contributing energy subband. These giant leaps could be   
calculated using the parabolic subbands approximation, obtained for k ≪ k_F 


    ε_γ^ν(k)  = 1/2[ 
    νλ_0 (1 - γ) + γ Δ] + [
    ħ^2/4 m_e (α + γβ ) + γ  (t_0 a_0)^2/Δ -νλ_0] k^2   .

This result leads to the DOS given by  <cit.>


    ρ_d(E) = 1/2πħ^2∑_γ, ν = ± 1| α+γβ/4 m_e + γ (t_0a_0)^2/ħ^2(Δ - νλ_0)|^-1Θ[
    γ(E -νλ_0/2) - 1/2( Δ -νλ_0 )
    ]  .

The two terms in Eq. (<ref>) are of the same order of magnitude, 
consequently, each of them must be retained in our calculation. Physically it means that due to
large gap Δ the curvature of each subband at k = 0 is so small that even the ⋍α/(4 m_e) 
correction is significant. It basically discards the DOS obtained for the gapped graphene model,
even at k → 0. The fact that the mass and even higher order terms must be taken into account 
for the plasmon calculation was mentioned in Ref. [Scholz1]. 




The principal model which yields quite reliable results for the DOS   
is derived by neglecting only non-essential O(k^4) terms with


    ε^ν_γ(k)⋍1/2 ν λ_0 + ħ^2 α/4 m_e k^2 + 
    γ/2√((Δ - ν λ_0 )^2
    +
    [
    (2 t_0 a_0)^2+( Δ - ν λ_0  ) ħ^2  β/m_e] k^2 
    )    .

Here, all ∽ k^2 terms are retained and the final expression for the DOS appears to be quite complicated. 
We use a general equation from our previous work  <cit.> to come up with a linear approximation 
which, according to the most precise and generalized numerical results,  <cit.> is valid 
for all experimentally allowable electron and hole doping densities. For the dispersions Eq. (<ref>), 
we use Eq. (<ref>) to obtain


    ρ_d(E) = 1/2π ∑_j ∑_γ, ν=± 1|
    α̃ + γ Ã_ν(Δ±λ_0, β |  a_0 t_0) / 2 {E-ϵ̃_ν -α̃ ξ_ν^ (j) (E)
    }|^-1Θ[
    γ(E -μλ_0/2) - 1/2( Δ - νλ_0 )
    ]  ,

where ϵ̃_ν=νλ_0/2, Δ̃_ν=(Δ -νλ_0)/2,
Ã_ν (Δ±λ_0, β |  a_0 t_0)=(Δ -νλ_0)ħ^2 β /(4 m_e) +
(t_0 a_0)^2, and α̃=ħ^2α/(4 m_e). ξ_ν^ (j)(E) are the roots of 


    f(E,ξ) = E - ϵ̃_ν - αξ - γ √(Ã_νξ + 
    Δ̃_ν^2) = 0   .

This euqation could be solved by expressing it  in quadratic form. However, one must bear in mind that if both 
parts of Eq. (<ref>) are squared, there might be additional non-physical solutions which must be disregarded.  
If this equation is written as (αξ)^2 + B  ξ + C = 0, where B = Ã_ν + 2  α and C = (E - ϵ̃_ν)^2 - Δ̃_ν^2, the only appropriate 
solution is ξ^ (1)(E) = 1/(2 α)^2 ( B + √(B^2 - 4 α^2 C)). 
The other solution
corresponds to the  E - ϵ̃_ν - αξ  = - γ √(Ã_νξ + Δ̃_ν^2) and is obviosly incorrect
for α→ 0.  We also note that the electron/hole index γ is no longer present 
in this equation, so that 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 of gapless graphene with ε_γ(k) = ħ v_F k with additional small, not depending on
γ, the mass term ∽α k^2, α≪ħ v_F/ k_F. The actual dispersion relation is now 
ϵ^γ(k) = γħ v_F | k | + α k^2, and the DOS is calculated as


    ρ_d(E) = 2/π ∑_γ = ± 1 ∑_j
     k^ (j)/|γ ħ v_F + 2 α k^ (j)|  .

The roots k^ (j) are the solutions of γħ v_F | k | + α k^2 = 0. Even though
 such a quadratic equation generally has two inequivalent solutions, only one of them 
± k^(1)⋍E/(γħ v_F) - α E^2/(ħ^3 v_F^3)
satisfies the | k |-type equation. Consequently, we obtain the following expression


    ρ_d(E) ⋍2/π ∑_γ = ± 1 E/γ (ħ v_F)^2 - 
      α (2 + γ)  E^2/
       (ħ v_F)^4

becoming equivalent to graphene DOS ρ_d(E) = 2/π E/ [γ  (ħ v_F)^2 ] 
for α→ 0. For the holes with γ = -1, however, the linear and quadratic mass terms are 
competing, so that for k ≫ k_F, another solution is present. However such wave vectors are beyond the Dirac 
cone model, and therefore, it is not physically acceptable. The obtained correction to the DOS is a small decrease, 
as it is expected to be for an energetically elevated location of the subband.
However, in our model for MoS_2, the mass terms ∽α and ∽β are not small and represent a 
finite correction to the DOS. The small parameter which we used in our series expansions is the energy 
δϵ above each band gap.



Now, we return to Eq. (<ref>) and present its solution as


    ξ^ (1) = 1/2  α̃^2 {Ã_ν + 2  α̃( E - ε̃_ν^(0))
     - [
     Ã_ν^2 + 4  α̃^2 Δ̃_ν^2 + 4  α̃ Ã_ν ( E -
     ε̃_ν^(0))
     ]^1/2}  .

This solution is exact in the sense that no approximations have been made so far excpet for 
the ∽ k^2 dispersions (<ref>). Substituted into Eq. (<ref>), 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. (<ref>).
We are interested in obtaining a linear approximation of the DOS next to 
each subband edge. Let us first consider electrons with  ε^ν_1(k) = 
Δ/2 + δϵ, δϵ≪E. In this case,


    ξ_ν^ (1)⋍4 m_e/α ħ^2 { 
     1 - 
     (a_0   t_0)^2 + ħ^2 β/(4 m_e)  ( Δ - νλ_0 
     )
     /
     (a_0   t_0)^2 + ħ^2 /(4 m_e)  [
     (α + β) ( Δ - νλ_0 
     )
     ]
     } δϵ  ,

and the DOS is now approximately given by  


    ρ_d (E) = 1/2 π ∑_ν=± 1 
          Δ - νλ_0
     /
     (a_0 t_0)^2 + ħ^2/(4 m_e)  (α + β)(Δ - νλ_0)
      + 
      2  δϵ  [
     (a_0   t_0)^2 + ħ^2 β/(4 m_e)  ( Δ - νλ_0 
     )  ]^2
     /{ 
     (a_0   t_0)^2 + ħ^2 /(4 m_e)  [
     (α + β) ( Δ - νλ_0 
     )
     ]  }^3
       .

The actual numerical resuts are  determined from
  

    ρ_d (E) = c_0^(3) + c_1^(3)( E - Δ/2)   , 
           c_0^(3) = 0.180 1/t_0 a_0^2 = 11.74  E_F^(0)/(ħ v_F)^2  , 
           c_1^(3) = 0.268  1/(t_0 a_0)^2 = 1.218  1/(ħ v_F)^2  .

In the valence band, we consider two separate hole subbands with ν = ± 1. 
If γ = -1, ν = 1 and E≤ -Δ/2 + λ_0,
the DOS is ρ_d(E) = c_0^(2) + c_1^(2) [ E -
( Δ/2 - λ_0  ) ], and the expansion coefficients are



    c_0^(2) = 1/2 π Δ - λ_0/(a_0 t_0)^2 + (β - α) (Δ - λ_0)  , 
           c_1^(2) = 1/π δϵ  [
     (a_0   t_0)^2 + ħ^2 β/(4 m_e)  ( Δ - λ_0 
     )  ]^2
     /{ 
     (a_0   t_0)^2 + ħ^2 /(4 m_e)  [
     (β - α) ( Δ - λ_0 
     )
     ]  }^3
      < 0   ,

or


    c_0^(2) = 0.105 1/t_0 a_0^2 = 6.847  E_F^(0)/(ħ v_F)^2  , 
           c_1^(2) = -0.232  1/(t_0 a_0)^2 = - 1.051  1/(ħ v_F)^2  .




Finally, for the lower hole subband with E⪅ - Δ/2 - λ_0, 
we obtain


    ρ_d(E) = c_0^(2) + c_1^(2) [ E - ( Δ/2 + λ_0  ) ]   , 
           c_0^(1) = 1/2 π ∑_ν = ± 1Δ -ν λ_0/(a_0 t_0)^2 + (β -
     α) (Δ - ν λ_0)  , 
           c_1^(1) = 1/π  ∑_ν = ± 1δϵ  [
     (a_0   t_0)^2 + ħ^2 β/(4 m_e)  ( Δ -ν λ_0 
     )  ]^2
     /{ 
     (a_0   t_0)^2 + ħ^2 /(4 m_e)  [
     (β - α) ( Δ - ν λ_0 
     )
     ]  }^3
      < 0   ,

and 


    c_0^(1) = 0.233 1/t_0 a_0^2 = 15.17  E_F^(0)/(ħ v_F)^2  , 
           c_1^(1) = -0.458  1/(t_0 a_0)^2 = 2.077  1/(ħ v_F)^2  .

 It is straightforward to obtain our previous results for the gapped graphene model (<ref>) DOS 
if α = β→ 0. We also note that the slope of the DOS in the conduction 
band is negative, as it should be accoring to Fig. <ref> (d), and the summation over the
ν index is present in all cases excpet the upper hole subband in Eqs. (<ref>).



Our results  (<ref>) - (<ref>) (here, we move from the conduction electrons to the valence band,
i.e., from the right to left) represent a fairly good match with the previously obtained numerical
values, specified in Sec. <ref> and later used for all our finite temperature calculations. The 
coefficients c_0^(i), i = 1, 2, 3 are equal to the giant discontinuities of the DOS
δρ_d (E) at each subband edge or k = 0, except c_0^(1)⋍δρ_d^(2)
(-Δ/2 + λ_0) + δρ_d^(1) (-Δ/2 - λ_0)), and, therefore, are accurate.
The linear coefficients c_1^(i) , in fair agreement, are 20-25 % larger compared with the 
numerical results since all the ∽ k^4 terms of our energy dispersions are neglected. Inclusion 
of these terms leads to the higher energies for chosen wave vector and a decrease of the     
DOS. This discrepancy is increased for higher energies, which is well seen for c_1^(1) for holes with 
E < -Δ/2 - λ_0. However,  the subbands for such energy ranges do not receive any substantial 
doping unless the temperature becomes very high. For most considered situations, we are limited for δϵ≈λ_0 within the band edges. In such a small range, c_1^(i) δϵ≪ c_0^(i), so 
that the actual DOS values remain almost unaffected and our model yields accurate results. 











§ 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)|_T=0. Our derivation
is based on the total carrier density conservation, which includes both electrons and holes, for zero and any
finite temperatures  <cit.>


    n = n^(e) + (-1) n^(h) = ∫_0^∞ d E ρ_d(E) f_γ = 1(E,T) - 
     ∫_-∞^0 d E ρ_d( E )  { 1 - f_γ = 1(E,T) }  .

The electron and hole occupation probabilities are complimentary and for electron
doping at T = 0 the hole states term has no effect on Eq. (<ref>) (for details see 
Ref. [proj]).



We begin  with the relatively simple case for silicene with dispersions (<ref>). The DOS 
ρ_d(E) for which buckled honeycomb lattices is given by Eq. (<ref>). 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 for doping densities 


    n ≤ n_c = 1/2 π Δ_>^2 - Δ_<^2/ħ^2 v_F^2 = 2/πħ^2 v_F^2 Δ_SOΔ_z   ,

and the Fermi energy is obtained form 


    n = 1/2 π E_F^2 - Δ_<^2/ħ^2 v_F^2  .

Alternatively, if the doping density is sufficient to populate both subbands, E_F is determined by

    n = 1/π 1/ħ^2 v_F^2 [ E_F^2 - 1/2 ( Δ_<^2 + Δ_>^2) ]   .




Once the temperature is set finite, Eq. (<ref>) leads us to  n = I^(e)(Δ_i,T)
-I^(h)(Δ_i,T), with the two terms corresponding to electron and hole components of 
the total carrier density. These integrals are presented as 


    I^(e)(Δ_i,T)  =  ∫_Δ_<^∞ d E A (E,T) + 
      ∫_Δ_>^∞ d E A (E,T)   ,

where


    A(E,T) = 1/π E/ħ^2 v_F^2 {1  +  exp[ E - μ(T)/k_B T] }^-1  .

Each of the these integrals could be easily evaluated. Using a variable substitution
 ξ = (E - Δ_< )/k_B T,   we obtain 


    1/π 1/(ħ v_F )^2∫_Δ_<^∞ d E E {1  +  exp[ E - μ(T)/k_B T] }^-1
    = k_B T/π(ħ v_F )^2∫_0^∞ d ξ  (Δ_< + ξ  k_B T)  { 1 + exp[ξ - μ(T) - Δ_</k_B T] }^-1

With the help of the following notation


    R^(p) (T,X) = ∫_0^∞ d ξ  ξ^p/ {1 + exp[ξ - X / (k_B T) ] }  ,

we obtain the final result of the integration as


    k_B T  Δ_</π( ħ v_F)^2 R^(0) [T,  μ(T) - Δ_<] +
      1/π (k_B T/ħ v_F)^2  R^(1) [T,  μ(T)-Δ_< ]   .

For p=0 and 1, corresponding to the  2DEG and gapless graphene, Eq. (<ref>) leads to  <cit.>


    R^(0) (T, X) = ln{ 1 + exp[X/k_B T] }  , 
          R^(1) (T, X) = - Li_ 2{ - exp[X/k_B T] }  ,

where Li_ 2 (z) is the second-order polylogarithm function or 
dilogarithm defined as


    Li_ p (z) = ∑_k=1^∞z^k/k^p  , 
          Li_ 2 (z) = - ∫_0^zln (1-t)/t  dt   .

The second term of Eq.(<ref>), which only differs from the first one by its integration limits, is 


    2 k_B T  Δ_>/π( ħ v_F )^2 ln{ 1 + 
    exp[μ(T) - Δ_>/k_B T] } - 1/π ( k_B T/ħ v_F)^2  Li_ 2{ - exp[μ(T)-Δ_>/k_B T] }  .

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) 0.1in where , 
          I^(h)_1 (Δ_i,T) = 1/π k_B T /( ħ v_F )^2 ∑_i = <,>Δ_i  R^(0){ T,   - [μ(T) + Δ_i ] }0.1in and
          I^(h)_2 (Δ_i,T) = 1/π  ( k_B T/ħ v_F)^2  ∑_i = <,>R^(1){ T,   - [μ(T) + Δ_i ] }  .

Now, the total carrier denisty from Eq. (<ref>) could be written as


    n  = ( k_B T/ħ v_F)^2 ∑_γ = ± 1 γ/π ∑_i = <,>R^(1)[ T,  γμ(T) - Δ_i  ] +
      Δ_i/k_B T R^(0)[ T,  γμ(T) - Δ_i ]

or, explicitly expressing the polylogarithm functions, we write


    n  ( ħ v_F/k_B T)^2 =  ∑_γ = ± 1 γ/π ∑_i = <,>
    - Li_ 2{ - exp[γμ(T) - Δ_i/k_B T] }  +
    Δ_i/k_B T ln{ 1 + exp[γμ(T) - Δ_i/k_B T] 
     }  .


Using Eqs. (<ref>) or Eq. (<ref>) depending on whether only one or both  subbands are filled at zero temperature, 
we obtain the equation which relates the finite-temperature chemical potential with its T = 0 value E_F.   
Energy band gap(s) obviously affects this result. 



The μ(T) for gapped graphene with two fourfold degenerate energy subbands is obtained if we substitute  
Δ_< = Δ_> = Δ_0 and ∑_i = <,>⟹× 2.
For gapless graphene, Δ_0 = 0 and π n = [ E_F/ (ħ v_F) ]^2, so that we write


    1/2 ( k_B T )^2  E_F^2  =  ∑_γ = ± 1γ R^(1)[T,  γμ(T) ]
     = - ∑_γ = ± 1γ Li_ 2{ - exp[γ μ(T)/k_B T] }  .

If the temperature is kept low with k_B T ≪ E_F, this result is simplified as  <cit.>


    R^1 (x) ⋍( x^2/2+ π^2/6)  Θ(x) + x  ln(1 + e^- | x |)   ,






Finalizing our derivations for silicene, we briefly address the case of hole doping with E_F < 0. The left 
part of Eq. (<ref>) is now modified as


    - n^(h) = - 2/π1/(ħ v_F)^2∑_i = < ,> ∫_- ∞^-Δ_i
     d  E |E|  Θ(-E + E_F)   .

In analogy with electron doing, the Fermi energy depends on whether only the Δ_<-subband (which 
is now the higher one) or they are both  doped. The equations determining the Fermi energy for given hole 
doping density n are exactly similar to Eqs. (<ref>) and (<ref>), which confirms complete symmetry 
between the electron and hole states in silicene. The right part of Eq. (<ref>) remains unchanged 
except the chemical potential is negative μ < 0 for any finite temperature. 





 §.§ μ(T) for transition metal dichalcogenides
 

In cotrast to the previously considered buckled honeycomb lattices and graphene, the electron/hole symmetry 
in TMDC's (such as MoS_2) is clearly broken. Even in the simplest gapped graphene model given by 
Eq. (<ref>), the two hole subbands are not degenerate and separated by λ_0 at k=0. For all 
reasonable doping densities n < 10^17 m^-2, the Fermi energy is such that the lower ε_ν 
= 1^γ = -1(k=0)= -Δ/2 - λ_0 subband could not be populated at zero temperature. 
This could be verified by rewriting Eq. (<ref>) as


    n_c = 2/πλ_0  Δ/(t_0   a_0)^2  = 1.0 · 10^18 m^-2  .




From here on, we are going to use a picewise-linear model for the DOS with the 
empirical coefficients, provided in Sec. <ref>. Our analytical model for the DOS, 
develop in the preceding Appendix in Eqs. (<ref>) - (<ref>) could also be employed here without
losing much of precision.   Let's us first consider electron doping with density n_(0)^e at zero temperature.  
The corresponding Fermi energy is determined by


    n_(0)^e = c_1^(3)/2  ( E_F^2 - Δ^2/4 ) + c_0^(3) ( E_F - Δ/2 )   ,

or 


    E_F^ e = 1/c_1^(3) {- c_0^(3) + [ ( c_0^(3) + c_1^(3) Δ/2)^2 + 2 n_(0)^e 
     c_1^(3)]^1/2}  .

For hole doping the result is quite similar, except c_1^(2) < 0  and the upper valence 
band gap is -Δ/2 + λ_0:


    E_F^(h) =  1/c_1^(2) { 
     - c_0^(2) + { 
     -2 c_1^(2)  n_(0)^h + [
     c_0^(2) - c_1^(2) (
     Δ/2 - λ_0
     )
     ]^2 
     }^1/2 }  .

Using this expression, we can further improve the result in Eq. (<ref>) for the critical hole doping
density n_c needed to reach the lower subband at zero temperature. The required Fermi energy is E_F^(h)≤ - Δ/2 - λ_0, so that the corresponding hole density is  
n_c^h = - λ_0  Δ   c_1^(2) + 2 c_0^(2) λ_0 = 3.3 · 10^17  m^-2.

This critical density value is still far above the experimentally realizable values ⋍ 1.0 · 10^17  
m^-2, and for all our calculations, the lower hole subband is never populated at T=0.



Finally, we address the finite-temperature chemical potential for MoS_2. Once again, we start with the carrier 
density conservation equation (<ref>).









The electron component is easily evaluated as









    k_B T  (c_0^(3) + Δ/2)  R^(0)[T,  μ(T) - Δ/2] +
    c_1^(3)(k_B T)^2  R^(1)[T,  μ(T)-Δ/2]   .

 

The hole state integrals are written in the following form


    ∫_-∞^0 d E ρ_d( |E|)  {1 - f_1(E,T) }  = 
    ∫_-∞^-Δ/2 + λ_0 d E [ - c_1^(2)E + c_0^(2)]
    {1  +  exp[ μ(T) - E/k_B T] }^-1 + 
           + ∫_-∞^-Δ/2 - λ_0 d E [ - δ c_1^(1)E + δ c_0^(1)] 
    {1  +  exp[ μ(T) - E/k_B T] }^-1  ,

where δ c_1^(i) = c_1^(i) - c_2^(i),   i = 1,2. They are evaluated as


    I^(h)(Δ,λ_0  |  T) = ∑_j = 1 ^4I^(h)_j0.1in where , 
          I^(h)_1 (Δ,λ_0  |  T) = k_B T (
      Δ/2 - λ_0 + c_2^(0))  R^(0){ T,   - [μ(T) + Δ/2 - λ_0 ] }  ,
          I^(h)_2 (Δ,λ_0  |  T) = c_2^(1)( k_B T )^2  R^(1){ T,   - [μ(T) + Δ/2 - λ_0 ] }  ,
          I^(h)_3 (Δ,λ_0  |  T) = k_B T (
      Δ/2 + λ_0 + δ c_1^(0))  R^(0){ T,   - [μ(T) + Δ/2 + λ_0 ] }  , 
          I^(h)_2 (Δ,λ_0  |  T) = δ c_1^(1)( k_B T )^2  R^(1){ T,   - [μ(T) + Δ/2 + λ_0 ] }  .

Combined with the electron terms (<ref>), these hole integrals (<ref>) form the right side 
of the carrier density conservation (<ref>). Its left side, corresponding to zero temperature,
is given by Eq. (<ref>) for electron doping and by 


    n_(0)^e = ( 
     Δ/2 + E_F - λ_0
     ) 
     {
     c_2^(0) + c_2^(1)/2[
     Δ/2 - (
     E_F + λ_0
     )
     ]
     }  ,

for hole doping. The symmetry between the electron and hole states is no longer present, which strongly
affects the finite-temperature behavior of the chemical potential in TMDC's. 









