Thermal Collective Excitations in Novel Two-Dimensional Dirac-Cone Materials

Abstract

The purpose of this chapter is to review some important, recent theoretical discoveries regarding the effect of temperature on the property of plasmons. These include their dispersion relations and Landau damping rates, and the explicit dependence of plasmon frequency on chemical potential at finite temperatures for a diverse group of recently discovered Dirac-cone materials. These novel materials cover gapped graphene, buckled howycomb lattices (such as silicene and germanene), molybdenum disulfide and other transition-metal dichalcogenides, especially the newest dice and α-T3 materials. The most crucial part of this review is a set of implicit analytical expressions about the exact chemical potential for each of considered materials, which greatly affects the plasmon dispersions and a lot of many-body quantum-statistical properties. We have also obtained the nonlocal plasmon modes of graphene which are further Coulomb-coupled to the surface of a thick conducting substrate, while the whole system is kept at a finite temperature. An especially rich physics feature is found for α-T3 materials, where each of the above-mentioned properties depends on both the hopping parameter α and temperature as well.

1. Introduction

Graphene, a two-dimensional (2D) carbon layer with a hexagonal atomic structure [1, 2, 3], has recently attracted outstanding attention from both academic scientists doing fundamental researches and engineers working on its technical applications [4]. Now, the scientific community is actively investigating the innovative semiconductors beyond graphene, with intrinsic spin-orbit interaction and tunable bandgap [5].

A remarkable feature of graphene is the absence of the bandgap in its energy dispersions. In spite of the obvious advantage of such bandstructure for novel electronic devices, electrons in graphene could not be confined due to the well-known Klein paradox [6]. To resolve this issue, graphene may be replaced with a material with a buckled structure and substantial spin-orbit interaction, such as silicene and germanene.

A new quasi-two-dimensional structure which has recently gained popularity among device scientists, is molybdenum disulfide monolayer, a honeycomb lattice which consists of two different molybdenum and sulfur atoms. It reveals a large direct band gap, absence of inversion symmetry and a substantial spin-orbit coupling. A summary of all recently fabricated materials beyond graphene is given in Figure 1. The last relevant example is black phosphorous (phosphorene) with a strong anisotropy of its composition and electron energy dispersion. Even though we do not study plasmons in phosphorene in the present chapter, there have been some crucial publications on that subject [7, 8].

Figure 1.

Recently discovered two-dimensional materials: graphene and beyond graphene.

Plasmons, or self-sustained collective excitations of interacting electrons in such low-dimensional materials, are especially important, since they serve as the basics for a number of novel devices and their applications [9, 10] in almost all fields of modern science, emerging nanofabrication and nanotechnology. Propagation and detection of plasmonic excitation in hybrid nanoscale devices can convert to or modify existing electromagnetic field or radiation [11, 12, 13, 14]. Localized surface plasmons are particularly of special interest considering their interactions with other plasmon modes in closely-located optoelectronic device as well as with imposed electromagnetic radiation [15].

Finite-temperature plasmons are of special interest for possible device applications. Among them is the possibility to increase the frequency (or energy) of a plasmon by an order of magnitude or even more, specifically, as a consequence of the raised temperature. As it was shown in Ref. [16], the dispersion of a thermal plasmon is given as $\sim \sqrt{q{k}_{B}T}$, where $q$ is the wave vector. This dispersion relation reveals the fact that the plasmon energy is monotonically increased with temperature and could be moved to the terahertz range and even above, which is crucial for imaging and spectroscopy.

At the same time, the damping rate, or broadening of the frequency, of such thermal plasmons varies as $\simeq 1/\sqrt{T}$, which means that we are dealing with a long-lived plasmon with low or even nearly zero Landau damping. Both plasmon frequency and the corresponding damping rate at finite temperature could be adjusted by electron doping, and could also be determined for an intrinsic material, where the chemical potential at $T=0K$ is located exactly at the Dirac point, while for zero temperature, intrinsic plasmons in graphene do not exist.

In this chapter, we will consider thermal behavior of plasmons, their dispersions and damping rates. By equipping with this information, it is possible to predict in advance the thermal properties of an electronic device designed for a particular temperature range. In spite of a number of reported theoretical studies on this subject [16, 17, 18], there is still a gap on demonstrating experimentally these unique thermal collective features of 2D materials. Therefore, our review can serve as an incentive to address this issue.

2. Novel two-dimensional materials beyond graphene

All the novel 2D materials considered here could be effectively assigned to an individual category based on their existing (or broken) symmetries and degeneracy in their low-energy band structure. We started with graphene having a bandgap ${\mathrm{\Delta }}_0$ and single-particle energy bands ${\epsilon }_{\pm }\left(k\right)=\pm \sqrt{{\left(ħv_{F}k\right)}^2+{\mathrm{\Delta }}_0^2}$, which are symmetric with respect to the Dirac point. Moreover, there is a total spin-valley degeneracy $g=g_s{g}_v=4$ for electrons and holes in each band.

Silicene and germanene, which represent buckled honeycomb lattices, possess subbands depending on valley and spin indices, and therefore are only doubly-degenerate. The electron-hole symmetry is broken for molybdenum disulfide and other transition-metal dichalcogenides (TMDC’s). For these situations, even though there exists a single electron-hole index $\gamma =\pm 1$, the energy of corresponding states does not have opposite values for each wave number, even at the valley point. In contrast to the electron states, the hole subbands reveal a splitting, as shown in Figure 2. All these partially broken symmetries strongly affect the chemical potential of 2D materials as well as its finite-temperature many-body properties. Black phosphorous, apart from all previously discussed broken symmetries, further acquires a preferred spatial direction in its atomic structure which leads to a strong anisotropy of its electronic states and band structures.

Figure 2.

Energy dispersions and density of states (DOS) ${\rho }_d\left(\mathbb{E}\right)$ of silicene [(a) and (b)] and molybdenum disulfide MoS₂ [(c) and (d)], where $E_F^{\left(0\right)}$ and $k_F^{\left(0\right)}$ are Fermi energy and wave number, respectively. For MoS₂, its dispersions and DOS, corresponding to parabolic band approximation in Eq. (23), are also shown for comparison.

2.1 Buckled honeycomb lattices

The energy dispersions of buckled honeycomb lattices, obtained from a Kane-Mele type Hamiltonian, appear as two inequivalent doubly-degenerate pairs of subbands with the same Fermi velocity $v_F=0.5\times {10}^6$ m/s and are given by

$${\epsilon }_{\xi ,\sigma }^{\gamma }\left(k\right)=\gamma \sqrt{{\left(\mathit{\xi \sigma }{\mathrm{\Delta }}_z-{\mathrm{\Delta }}_{\mathit{SO}}\right)}^2+{\left(ħv_{F}k\right)}^2},$$E1

where $\gamma =\pm 1$ labels symmetric electron and hole states. Here, two bandgaps [19, 20] ${\mathrm{\Delta }}_{<}=\mid {\mathrm{\Delta }}_{\mathit{SO}}-{\mathrm{\Delta }}_z\mid$ and ${\mathrm{\Delta }}_{>}={\mathrm{\Delta }}_{\mathit{SO}}+{\mathrm{\Delta }}_z$ are attributed to an intrinsic spin-orbit gap ${\mathrm{\Delta }}_{\mathit{SO}}=0.5-3.5$ meV [21, 22, 23, 24] and a tunable asymmetry bandgap ${\mathrm{\Delta }}_z$ proportional to applied electric field ${\mathcal{E}}_z$. The band structure, however, depends only on one composite index $\nu =\mathit{\sigma \xi }$, a product of spin $\sigma$ and valley $\xi$ index. At ${\mathcal{E}}_z=0$, two gaps become the same. As ${\mathcal{E}}_z\ne 0$, ${\mathrm{\Delta }}_{<}$ and ${\mathrm{\Delta }}_{>}$ change in opposite ways, and electrons stay in a topological insulator (TI) state. Additionally, ${\mathrm{\Delta }}_{<}$ decreases with ${\mathcal{E}}_z$ until reaching zero, corresponding to a new valley-spin polarized metal. On the other hand, if ${\mathcal{E}}_z$ further increases, both ${\mathrm{\Delta }}_{<}$ and ${\mathrm{\Delta }}_{>}$ will be enhanced, leading to a standard band insulator (BI) state for electrons. As we will see below, all of the single electronic and collective properties of buckled honeycomb lattices depend on both bandgaps ${\mathrm{\Delta }}_{<,>}$, and therefore could be tuned by varying a perpendicular electric field to create various types of functional electronic devices.

The wave function of silicene, corresponding to eigenvalue equation in Eq. (1), takes the form [25]

$$\begin{array}{l}{\mathrm{\Psi }}^{\gamma }\left(k\right)=\left[\begin{array}{c}{\mathrm{\Psi }}_{\xi =1,\sigma =+1}^{\gamma }\left(k\right)\\ {\mathrm{\Psi }}_{\xi =1,\sigma =-1}^{\gamma }\left(k\right)\end{array}\right],\\ {\mathrm{\Psi }}_{\xi ,\sigma }^{\gamma }\left(k\right)=\sqrt{\frac{\gamma }{2{\epsilon }_{\xi ,\sigma }^{\gamma }\left(k\right)}}\left[\begin{array}{c}\sqrt{\mid {\epsilon }_{\xi ,\sigma }^{\gamma }\left(k\right)+{\mathrm{\Delta }}_0^{\xi ,\sigma }\mid }\\ \\ \gamma \sqrt{\mid {\epsilon }_{\xi ,\sigma }^{\gamma }\left(k\right)-{\mathrm{\Delta }}_0^{\xi ,\sigma }\mid }{e}^{i{\theta }_k}\end{array}\right],\end{array}$$E2

where ${\theta }_k={tan}^{-1}\left(k_y/k_x\right)$ and ${\mathrm{\Delta }}_0^{\xi ,\sigma }=\mid \mathit{\xi \sigma }{\mathrm{\Delta }}_z-{\mathrm{\Delta }}_{\mathit{SO}}\mid$.

Germanene, another representative of buckled honeycomb lattices [26, 27, 28, 29, 30], demonstrates substantially higher Fermi velocities and an enhanced intrinsic bandgap $\sim 23$ meV. For a free-standing germanene, first-principles studies have revealed a buckling distances $\sim \mathrm{0.640.74}$ Å [31, 32].

2.2 Molybdenum disulfide and transition-metal dichalcogenides

MoS₂ is a typical representative of transition-metal dichalcogenide (TMDC) monolayers. TMDC’s are semiconductors with the composition of TC₂ type, where T refers to a transition-metal atom, such as Mo or W, while C corresponds to a chalcogen atom (S, Se or Te).

MoS₂ displays broken inversion symmetry and direct bandgaps. Its most crucial distinction from the discussed buckled honeycomb lattices is its broken symmetry between the electrons and holes so that the corresponding energy bands are no longer symmetric with respect to the Dirac point, but could still be classified by a single index $\gamma =\pm 1$. The absence of this particle-hole symmetry is expected to have a considerable effect on the plasmon branches at both zero and finite temperatures through the thermal convolution of the corresponding quantum states.

Specifically, the energy bands of MoS₂ can be described by a two-band model, i.e.,

$${\epsilon }_{\gamma }^{\xi ,\sigma }\left(k\right)\simeq \frac{1}{2}{\mathit{\xi \sigma \lambda }}_0+\frac{\alpha {ħ}^2}{4m_e}{k}^2+\frac{\gamma }{2}\sqrt{\left[{\left(2t_0{a}_0\right)}^2+\left(\mathrm{\Delta }-{\mathit{\xi \sigma \lambda }}_0\right)\beta {ħ}^2/m_e\right]k^2+{\left(\mathrm{\Delta }-{\mathit{\xi \sigma \lambda }}_0\right)}^2},$$E3

where $\mathrm{\Delta }=1.9$ eV is a gap parameter leading to an extremely large direct bandgap $\simeq 1.7$ eV. There is also substantial internal spin-orbit coupling ${\lambda }_0=0.042\mathrm{\Delta }$, $t_0{a}_0=4.95\times {10}^{-29}$ J m is a Dirac-cone term, where $t_0=0.884\mathrm{\Delta }$ is the electron hopping parameter while $a_0=1.843$ Å is the lattice constant. The $t_0{a}_0$ term is $\approx 0.47$ compared to $ħv_F$ in graphene. Similarly to silicene and germanene, the energy dispersions of TMDC’s depend on one composite valley-spin index $\nu =\mathit{\sigma \xi }$. There are also other less important but still non-negligible $\sim k^2$ mass terms with $\alpha =2.21=5.140\beta$, and $m_e$ represents the mass of a free electron.

In practical, we will neglect the $\simeq k^4$ terms, $\simeq t_1{a}_0^2$ trigonal warping term and anisotropy, which indeed have tiny or no effect on the density of states of the considered material, but would make our model much more complicated. The above dispersions could be presented in a form similarly to those for gapped graphene, i.e., ${\epsilon }_{\gamma }^{\nu }\left(k\right)={\mathbb{E}}_0^{\nu }\left(k\right)+\gamma \sqrt{{\left[{\mathrm{\Delta }}_0^{\nu }\left(k\right)\right]}^2+{\left(t_0{a}_{0}k\right)}^2}$, where ${\mathrm{\Delta }}_0^{\nu }\left(k\right)={ħ}^2{k}^2\beta /4m_e+\mathrm{\Delta }/2-{\mathit{\nu \lambda }}_0/2$ is a $k$-dependent “gap term” and the band shift is ${\mathbb{E}}_0^{\nu }\left(k\right)={ħ}^2{k}^2\alpha /4m_e+{\mathit{\nu \lambda }}_0/2$. A set of somewhat cumbersome analytical expressions for the components of the wave functions corresponding to dispersions in Eq. (3) can be found from Ref. [25].

Using Eq. (3), we can verify that the degeneracy of two hole subbands ($\gamma =-1$), corresponding to $\nu =\pm 1$, will be lifted and two subband will be separated by ${\lambda }_0\simeq 79.8$ meV. However, this is not the case for two corresponding electron states ($\gamma =1$). Consequently, the electron-hole asymmetry exists even at $k=0$ and becomes even more pronounced at finite $k$ values. One can clearly see this difference by comparing Figure 2(b) with Figure 2(d).

3. Thermal plasmons in graphene and other materials

One of the most important features in connection with plasmons at zero and finite temperatures is its dispersion relations, i.e., dependence of the plasmon frequency $\omega$ on wave number $q$. Physically, these complex relations can be determined from the zeros of a dielectric function ${\mathrm{ϵ}}_T\left(q,\omega \right)$, [19, 33] given by

$${\mathrm{ϵ}}_T\left(q,\omega \right)=1-v\left(q\right){\mathrm{\Pi }}_T\left(q,\omega |\mu \left(T\right)\right)=0,$$E4

where $v\left(q\right)=2{\mathit{\pi \alpha }}_r/q\equiv e^2/2{\mathrm{ϵ}}_0{\mathrm{ϵ}}_{r}q$ is the 2D Fourier-transformed Coulomb potential, ${\alpha }_r=e^2/4\pi {\mathrm{ϵ}}_0{\mathrm{ϵ}}_r$, and ${\mathrm{ϵ}}_r$ represents the dielectric constant of the host material.

The dielectric function introduced in Eq. (4) is determined directly by the finite-temperature polarization function, or polarizability, ${\mathrm{\Pi }}_T\left(q,\omega |\mu \left(T\right)\right)$, which is, in turn, related to its zero-temperature counterpart, ${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)$, by an integral convolution with respect to different Fermi energies [34], given by

$${\mathrm{\Pi }}_T\left(q,\omega |\mu \left(T\right)\right)=\frac{1}{2k_{B}T}\underset{0}{\overset{\infty }{\int }}\mathit{d\eta }\frac{{\mathrm{\Pi }}_0\left(q,\omega |\eta \right)}{1+cosh\left[\left(\mu -\eta \right)/k_{B}T\right]},$$E5

where the integration variable $\eta$ stands for the electron Fermi energy at $T=0$. This equation is derived for electron doping with $\eta =E_F>0$. We note that, in order to evaluate this integral, one needs to know in advance how the chemical potential $\mu \left(T\right)$ varies with temperature $T$. Such a unique $T$ dependence reflects a specific selection of a convolution path for a particular material band structure, which we will discuss in Section 4.

The zero-temperature polarizability, which is employed in Eq. (5), is quite similar for all 2D materials considered here. The only difference originates from the degeneracy level of the low-energy band structure, such as $g=g_v{g}_s=4$ for graphene with either finite or zero bandgap. We begin with the expression of the partial polarization function with two inequivalent doubly-degenerate pairs of subbands labeled by a composite index $\nu$

$${\mathrm{\Pi }}_0^{\left(\nu \right)}\left(q,\omega |E_F\right)=\frac{1}{4{\pi }^2}\int d^2\mathit{k}\sum _{\gamma ,{\gamma }^{\prime }=\pm 1}\left\{1+\gamma {\gamma }^{\prime }\frac{\mathit{k}\cdot \left(\mathit{k}+\mathit{q}\right)+{\mathrm{\Delta }}_{\nu }^2}{{\epsilon }_{\gamma }^{\nu }\left(k\right){\epsilon }_{{\gamma }^{\prime }}^{\nu }\left(|\mathit{k}+\mathit{q}|\right)}\right\}\times \frac{{\mathrm{\Theta }}_0\left[E_F-{\epsilon }_{\gamma }^{\nu }\left(k\right)\right]-{\mathrm{\Theta }}_0\left[E_F-{\epsilon }_{{\gamma }^{\prime }}^{\nu }\left(|\mathit{k}+\mathit{q}|\right)\right]}{ħ\left(\omega +i{0}^{+}\right)+{\epsilon }_{\gamma }^{\nu }\left(k\right)-{\epsilon }_{{\gamma }^{\prime }}^{\nu }\left(|\mathit{k}+\mathit{q}|\right)},$$E6

where ${\mathrm{\Theta }}_0\left(x\right)$ stands for a unit-step function, and $\gamma =\pm 1$ stands for the electron or hole state with energy dispersions above or below the Dirac point. Moreover, the index $\nu$, which equals to $\mathit{\sigma \xi }=\pm 1$ for buckled honeycomb lattices or molybdenum disulfide, specifies two different pairs of degenerate subbands from Eq. (1) or Eq. (3).

Finally, the full polarization function at zero temperature is obtained as

$${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)=\sum _{\nu =\pm 1}{\mathrm{\Pi }}_0^{\left(\nu \right)}\left(q,\omega |E_F\right).$$E7

If the dispersions of low-energy subbands do not depend on the valley or spin indices $\xi$ and $\sigma$, summation in Eq. (7) simply gives rise to a factor of two, as we have obtained for graphene.

Integral transformation in Eq. (5), which is used to obtain the finite-temperature polarization function from its zero-temperature counterparts with different Fermi energies, was first introduced in Ref. [34]. It could be derived in a straightforward way by noting that the only quantity which substantially depends on temperature in Eq. (6) is the Fermi-Dirac distribution function $n_F\left[{\epsilon }_{\gamma }^{\nu }\left(k\right)-\mu \left(T\right)\right]$. It changes to the unit-step functions ${\mathrm{\Theta }}_0\left[E_F-{\epsilon }_{\gamma }^{\nu }\left(k\right)\right]$ at $T=0$, as used in Eq. (6). As the temperature $T$ increases from zero, the distribution function in Eq. (6) evolves into [35, 36]

$$n_F\left[{\epsilon }_{\gamma }^{\nu }\left(k\right)-\mu \left(T\right)\right]=\frac{1}{2}\left[1-tanh\left(\frac{{\epsilon }_{\gamma }^{\nu }\left(k\right)-\mu \left(T\right)}{2k_{B}T}\right)\right]=\underset{0}{\overset{\infty }{\int }}\mathit{d\eta }\frac{{\mathrm{\Theta }}_0\left[\mu \left(T\right)-{\epsilon }_{\gamma }^{\nu }\left(k\right)\right]}{4k_{B}T{cosh}^2\left\{\left[\mu \left(T\right)-\eta \right]/2k_{B}T\right\}}.$$E8

For accessible temperatures, the energy dispersions ${\epsilon }_{\gamma }^{\nu }\left(k\right)$, corresponding wave functions and their overlap factors are all temperature independent. As a result, the polarization function is expected to be modified by the same integral transformation, or a convolution, as each of the Fermi-Dirac distribution function in the numerator of Eq. (6).

We first look at intrinsic plasmons with $E_F=0$ at $T=0$. In this case, $\mu \left(T\right)$ also remains at the Dirac point for any temperature $T$. As $T$ increases to $k_{B}T\gg E_F>0$, on the other hand, ${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)$ for gapless graphene gives rise to a plasmon dispersion ${\omega }_p\simeq qT$ in the long-wave limit and the damping rate is $\sim q^{3/2}/\sqrt{T}$. As a result, the plasmon mode becomes well defined [6] for $q<16{\mathrm{ϵ}}_0{\mathrm{ϵ}}_r{k}_{B}T/\pi e^2$.

Additionally, finite-$T$ polarization function ${\mathrm{\Pi }}_T\left(q,\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }\right)$ of a 2D material is directly related to its optical conductivity ${\sigma }_O^{\left(T\right)}\left(\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }\right)$ through [19]

$${\sigma }_O^{\left(T\right)}\left(\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }\right)=\mathit{i\omega }{e}^2\underset{q\to 0}{lim}\frac{{\mathrm{\Pi }}_T\left(q,\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }\right)}{q^2},$$E9

where we introduce the notation ${\mu }_T\equiv \mu \left(T\right)$, ${\mathrm{\Pi }}_T\left(q,\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }\right)\sim q^2$ as $q\to 0$ for each of our considered 2D materials, regardless of their band structure, as given by Eq. (47) for $T=0$. This conclusion holds true even for finite $T$ and makes the optical conductivity independent of $q$, and therefore the $q\to 0$ limit in Eq. (9) becomes finite.

From Eq. (9), we find explicitly that

$$\begin{array}{cc}Im\left[{\sigma }_O^{\left(0\right)}\left(\omega |E_F,{\mathrm{\Delta }}_{\beta }\right)\right]& =\frac{e^2}{4\pi ħ}\sum _{\beta =\pm 1}\left\{\frac{4E_F}{ħ\omega }\left[1-{\left(\frac{{\mathrm{\Delta }}_{\beta }}{E_F}\right)}^2\right]\right.+\left.\left[1+{\left(\frac{2{\mathrm{\Delta }}_{\beta }}{ħ\omega }\right)}^2\right]ln\left|\frac{2E_F-ħ\omega }{2E_F+ħ\omega }\right|\right\},\\ Re\left[{\sigma }_O^{\left(0\right)}\left(\omega |E_F,{\mathrm{\Delta }}_{\beta }\right)\right]& =\frac{e^2}{4ħ}\mathrm{\Theta }\left(ħ\omega -2E_F\right)\sum _{\beta =\pm 1}\left[1+{\left(\frac{2{\mathrm{\Delta }}_{\beta }}{ħ\omega }\right)}^2\right].\end{array}$$E10

Here, the state-blocking effect due to Pauli exclusion principle directly results in the diminishing of the real part of the optical conductivity at zero temperature for $ħ\omega <2E_F$. However, if $T>0$, such state-blocking effect will not exist [37, 38, 39, 40] due to

$$\mathrm{\Theta }\left(E_F-{\epsilon }_{\beta }^{\gamma }\left(k\right)\right)\Rightarrow \frac{1}{2}\left\{1-tanh\left[\frac{{\epsilon }_{\beta }^{\gamma }\left(k\right)-{\mu }_T}{2k_{B}T}\right]\right\}.$$E11

Furthermore, for gapless (${\mathrm{\Delta }}_{\beta }=0$) but doped ($E_F>0$) graphene in the high-$T$ limit we obtain from Eq. (10)

$${\sigma }_O^{\left(T\right)}\left(\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }=0\right)\simeq \frac{e^2}{ħ}\left\{\frac{ħ\omega }{16k_{B}T}\left[1-\frac{1}{3}{\left(\frac{ħ\omega }{4k_{B}T}\right)}^2\right]\right.+i\frac{2ln2k_{B}T}{\pi ħ\omega }\left.\left[1+2ln2{\left(\frac{E_F}{4ln2k_{B}T}\right)}^4\right]\right\},$$E12

where we have used the high-$T$ limit [17] for the chemical potential ${\mu }_T\approx \left(E_F^2/4ln2k_{B}T\right)$. In either case above, we have to present analytical expression for $\mu \left(T\right)$ as a function of $T$ so as to gain the explicit $T$ dependence of optical conductivity. From Eq. (12) we conclude that $Im\left[{\sigma }_O^{\left(T\right)}\left(\omega |{\mu }_T,{\mathrm{\Delta }}_{\beta }=0\right)\right]$ depends weakly on $E_F$.

On the other hand, for gapped (${\mathrm{\Delta }}_{\beta }={\mathrm{\Delta }}_0$) but undoped ($E_F=0$) graphene at high $T$ ($k_{B}T\gg {\mathrm{\Delta }}_0$ and $ħ\omega$), we get its optical conductivity [41]

$$\begin{array}{l}Re{\left[{\sigma }_O^{\left(T\right)}\left(\omega |{\mu }_T,{\mathrm{\Delta }}_0\right)\right]}_{E_F=0}=\frac{e^2}{16ħ}\left(\frac{ħ\omega }{k_{B}T}\right)\left(1-\frac{{\mathrm{\Delta }}_0}{ħ\omega }\right),\\ Im{\left[{\sigma }_O^{\left(T\right)}\left(\omega |{\mu }_T,{\mathrm{\Delta }}_0\right)\right]}_{E_F=0}=\frac{4e^2}{\pi ħ}\left(\frac{k_{B}T}{ħ\omega }\right)\left\{2ln2-{\left(\frac{{\mathrm{\Delta }}_0}{k_{B}T}\right)}^2\left[{\mathbb{C}}_0-ln\left(\frac{{\mathrm{\Delta }}_0}{2k_{B}T}\right)\right]\right\},\end{array}$$E13

where the constant ${\mathbb{C}}_0\simeq 0.79$ appears due to a finite bandgap. [36].

4. Chemical potential at finite temperatures

As we have seen from Section 3, we need know $\mu \left(T\right)$ as a function of $T$ explicitly so as to gain $T$ dependence of polarization function, plasmon, transport and optical conductivities, or any other quantities related to collective behaviors of 2D materials at finite temperatures [42].

The density of states (DOS), which plays an important tool in calculating electron (or hole) Fermi energy $E_F$ and chemical potential $\mu \left(T\right)$, is defined as

$${\rho }_d\left(\mathbb{E}\right)=\sum _{\gamma =\pm 1}\sum _{\xi ,\sigma =\pm 1}\int \frac{d^2\mathit{k}}{{\left(2\pi \right)}^2}\delta \left(\mathbb{E}-{\epsilon }_{\xi ,\sigma }^{\gamma }\left(k\right)\right),$$E14

where $\delta \left(x\right)$ is Dirac delta function. Using Eq. (14), we immediately obtain a piece-wise linear function for silicene [43]

$${\rho }_d\left(\mathbb{E}\right)=\frac{1}{\pi }\sum _{\gamma =\pm 1}\frac{\gamma \mathbb{E}}{{ħ}^2{v}_F^2}\sum _{i=<,>}{\mathrm{\Theta }}_0\left(\frac{\mathbb{E}}{\gamma }-{\mathrm{\Delta }}_i\right).$$E15

This result is equivalent to the DOS of graphene except that there are no states within the bandgap region, as demonstrated by two unit-step functions ${\mathrm{\Theta }}_0\left(|\mathbb{E}|-{\mathrm{\Delta }}_{<}\right)$ and ${\mathrm{\Theta }}_0\left(|\mathbb{E}|-{\mathrm{\Delta }}_{>}\right)$.

Finally, the chemical potential $\mu \left(T\right)$ can be calculated using the conservation of the difference of electron and hole densities, [17] $n_e\left(T\right)$ and $n_h\left(T\right)$, for all temperatures including $T=0$, i.e.,

$$n=n_e\left(T\right)-n_h\left(T\right)=\underset{0}{\overset{\infty }{\int }}d\mathbb{E}{\rho }_d\left(\mathbb{E}\right)f_{\gamma =1}\left(\mathbb{E},T\right)-\underset{-\infty }{\overset{0}{\int }}d\mathbb{E}{\rho }_d\left(\mathbb{E}\right)\left[1-f_{\gamma =1}\left(\mathbb{E},\mathbb{T}\right)\right],$$E16

where $f_{\gamma =1}\left(\mathbb{E},T\right)={\left\{1+exp\left[\left(\mathbb{E}-\mu \left(T\right)\right)/k_{B}T\right]\right\}}^{-1}$ is the Fermi function for electrons in thermal equilibrium. The hole distribution function is just $f_{\gamma =-1}\left(\mathbb{E}<0,T\right)=1-f_{\gamma =1}\left(\mathbb{E},T\right)$.

At $T=0$, it is straightforward to get the Fermi energy $E_F$ from Eq. (16) for silicene

$$E_F^2-\frac{1}{2}\left({\mathrm{\Delta }}_{>}^2+{\mathrm{\Delta }}_{<}^2\right)={\left(ħv_F\right)}^2\pi n.$$E17

where we have assumed that both subbands are occupied for simplicity. The discussions of other cases can be found from Ref. [18]. Consequently, minimum electron density required to occupy the upper subband of silicene is $n_c=2{\mathrm{\Delta }}_{\mathit{SO}}{\mathrm{\Delta }}_z/\pi {ħ}^2{v}_F^2$.

On the other hand, by applying Eq. (16), in combination with DOS in Eq. (15), for silicene, a transcendental (non-algebraic) equation [43, 44] could be obtained for $\mu \left(T\right)$, that is

$${\left(\frac{ħv_F}{k_{B}T}\right)}^{2}n=\frac{1}{\pi }\sum _{\gamma =\pm 1}\gamma \sum _{i=<,>}\left\{-{\mathrm{Li}}_2\left(-exp\left[\frac{\gamma \mu \left(T\right)-{\mathrm{\Delta }}_i}{k_{B}T}\right]\right)\right.+\left.\frac{{\mathrm{\Delta }}_i}{k_{B}T}ln\left(1+exp\left[\frac{\gamma \mu \left(T\right)-{\mathrm{\Delta }}_i}{k_{B}T}\right]\right)\right\},$$E18

where ${\mathrm{Li}}_2\left(x\right)$ is a polylogarithm or dilogarithm function, defined mathematically by

$${\mathrm{Li}}_2\left(z\right)=-\underset{0}{\overset{z}{\int }}\mathit{dt}\left[\frac{ln\left(1-t\right)}{t}\right].$$E19

Interestingly, the right-hand side of Eq. (18) contains terms corresponding to both pristine and gapless graphene, using which we find from Ref. [17].

$$\frac{1}{2{\left(k_{B}T\right)}^2}{E}_F^2=-\sum _{\gamma =\pm 1}\gamma {\mathrm{Li}}_2\left\{-exp\left[\frac{\gamma \mu \left(T\right)}{k_{B}T}\right]\right\},$$E20

as well as a well-known analytical expression for ${\mu }_0\left(T\right)$ of 2D electron gas with Schrödinger-based electron dynamics

$${\mu }_0\left(T\right)=k_{B}Tlog\left(1+exp\left[\frac{\pi {ħ}^2{n}_0}{m^{\ast }{k}_{B}T}\right]\right).$$E21

An advantage of Eq. (18) is that it could be solved even without taking an actual integration. In fact, one can either readily solve it numerically using some standard computational algorithms, or introduce an analytical approximation to the sought solution near specific temperature assigned.

Numerical results for $\mu \left(T\right)$ of silicene are presented in Figure 3. In all cases, $\mu \left(T\right)$ decreases with increasing $T$ from zero. However, it is very important to notice that $\mu \left(T\right)$ never reaches zero or changes its sign in the systems with an electron–hole symmetry due to increasing contribution from holes in Eq. (16). All of displayed results in Figure 3 depend on individual bandgaps ${\mathrm{\Delta }}_{<}$ and ${\mathrm{\Delta }}_{>}$. The special case with ${\mathrm{\Delta }}_{<}={\mathrm{\Delta }}_{>}$ corresponds to gapped graphene, for which plasmon modes at $T=0$ were studied in Ref. [33]. The graphene Fermi wave number is $k_F=\sqrt{\mathit{\pi n}}$, irrelevant to its bandgap. The general relation between $k_F$ and $n$ in 2D materials is given by ${\left(2\pi \right)}^{2}n=g\pi k_F^2$. The experimentally allowable electron (or hole) doping is within the range of $n={10}^{10}-{10}^{12}$ cm${}^{-2}$, leading to $k_F={10}^6-{10}^7$ cm${}^{-1}$. For two pair of inequivalent subbands, such as in silicene or MoS₂, there are two different Fermi wave numbers for these subbands. Moreover, the numerically calculated $\mu \left(T\right)$ as functions of $T$ for electron and hole doping are presented in Figure 4. In this case, however, there exists no electron-hole symmetry, and therefore the resulting $\mu \left(T\right)$ can be zero and change its sign as $T$ increases, in contrast to the results in Figure 3.

Figure 3.

Temperature dependence of the chemical potential $\mu \left(T\right)$ for silicene with two inequivalent energy subbands with various bandgaps and a fixed doping density $n=1\times {10}^{11}$cm${}^{-2}$. Panel (a) highlights the situation close to $\mathrm{T}=0$, while panel (b) shows the whole temperature range. Here, ${\mathrm{E}}_0$ is the Fermi energy of graphene.

Figure 4.

Temperature dependence of the chemical potential $\mu \left(T\right)$ for molybdenum disulfide for cases of electron (a) and hole (b) doping with various doping densities. $\mu \left(T\right)$ might change its sign in contrast to the previously considered graphene and silicene. The two insets demonstrate how the Fermi energy depends on the electron and hole doping densities, correspondingly.

Eq. (18) could also be applied to a wide range of 2D materials if its DOS has a linear dependence on energy $\mathbb{E}$. Particularly, it is valid for calculating the finite-$T$ chemical potential of TMDC’s with an energy dispersion presented in Eq. (3). However, we are aware that some terms in Eq. (3) for TMDC’s, which might be insignificant for dispersions of other 2D materials, become essential in DOS because of very large bandgap and mass terms around $k=0$. As an estimation, for $k/k_0\approx 5.0$, the correction term $\simeq \beta \mathrm{\Delta }{k}^4$ must be taken into account. Meanwhile, the highest accessible doping $n={10}^{13}$ cm${}^{-2}$ only gives rise to a Fermi energy $E_F\sim {\lambda }_0$, comparable to spin-orbit coupling.

Now, we turn to calculate $\mu \left(T\right)$ as a function of $T$ for MoS₂ with a much more complicated band structure. After taking into account the $\sim k^2$ mass terms, we are able to write down [18]

$${\rho }_d\left(\mathbb{E}\right)=\frac{1}{2\pi {ħ}^2}\sum _{\gamma ,\nu =\pm 1}{\left|\frac{\alpha +\mathit{\gamma \beta }}{4m_e}+\frac{\gamma {\left(t_0{a}_0\right)}^2}{{ħ}^2\left(\mathrm{\Delta }-{\mathit{\nu \lambda }}_0\right)}\right|}^{-1}{\mathrm{\Theta }}_0\left(\gamma \left[\mathbb{E}-\frac{{\mathit{\nu \lambda }}_0}{2}\right]-\frac{1}{2}\left(\mathrm{\Delta }-{\mathit{\nu \lambda }}_0\right)\right),$$E22

where the calculation is based on a parabolic-subband approximation, i.e.,

$${\epsilon }_{\gamma }^{\nu }\left(k\right)=\frac{1}{2}\left[{\mathit{\nu \lambda }}_0\left(1-\gamma \right)+\gamma \mathrm{\Delta }\right]+\left[\frac{{ħ}^2}{4m_e}\left(\alpha +\mathit{\gamma \beta }\right)+\frac{\gamma {\left(t_0{a}_0\right)}^2}{\mathrm{\Delta }-{\mathit{\nu \lambda }}_0}\right]k^2.$$E23

From Eq. (22), we further seek an explicit expression for DOS in the form a piecewise-linear function of energy $\mathbb{E}$: ${\rho }_d\left(\mathbb{E}\right)=A_i+B_i\mathbb{E}$. A complete set of expressions for DOS of MoS₂ has been reported in Ref. [18]. Here, we merely provide and discuss these DOS expression around the lower hole subband with $\mathbb{E}\approx -\mathrm{\Delta }/2-{\lambda }_0$, yielding

$${\rho }_d\left(\mathbb{E}\right)=c_0^{\left(1\right)}+\left[\mathbb{E}-\left(\frac{\mathrm{\Delta }}{2}+{\lambda }_0\right)\right]c_1^{\left(1\right)},$$

$$c_0^{\left(1\right)}=\frac{1}{2\pi }\sum _{\nu =\pm 1}\frac{\mathrm{\Delta }-\nu {\lambda }_0}{{\left(a_0{t}_0\right)}^2+\left(\beta -\alpha \right)\left(\mathrm{\Delta }-\nu {\lambda }_0\right)},$$

$$c_1^{\left(1\right)}=\frac{1}{\pi }\sum _{\nu =\pm 1}\frac{{\left[{\left(a_0{t}_0\right)}^2+{ħ}^2\beta /\left(4m_e\right)\left(\mathrm{\Delta }-\nu {\lambda }_0\right)\right]}^2\mathit{\delta \epsilon }}{{\left\{{\left(a_0{t}_0\right)}^2+{ħ}^2/\left(4m_e\right)\left[\left(\beta -\alpha \right)\left(\mathrm{\Delta }-\nu {\lambda }_0\right)\right]\right\}}^3}<0,$$E24

or numerically,

$$c_0^{\left(1\right)}=0.233\frac{1}{t_0{a}_0^2}=15.17\frac{E_0}{{\left(ħv_F\right)}^2},$$

$$c_1^{\left(1\right)}=-0.458\frac{1}{{\left(t_0{a}_0\right)}^2}=-2.077\frac{1}{{\left(ħv_F\right)}^2}.$$E25

The calculated numerical results for DOS in all regions are listed in Table 1. All introduced coefficients $c_0^{\left(i\right)},c_1^{\left(i\right)}$ for $i=1,2,3$ can be deduced from the calculated parameters $A_i$ and $B_i$ using a similar correspondence as in Eq. (24).

Range index	Energy range	$\gamma$	$\nu$	$A_i$[$1/\left(t_0{a}_0^2\right)$]	$B_i$[$1/{\left(t_0{a}_0\right)}^2$]
$i=1$	$\mathbb{E}<-\mathrm{\Delta }/2-{\lambda }_0$	$-1$	$+1$	$0.0174$	$-0.169$
$i=2$	$\mid \mathbb{E}+\mathrm{\Delta }/2\mid <{\lambda }_0$	$-1$	$-1$	$0.043$	$-0.308$
$i=3$	$\mathbb{E}>\mathrm{\Delta }/2$	$+1$		$0.078$	$+0.179$

Table 1.

Linearized density of states (DOS) ${\rho }_d\left(\mathbb{E}\right)=A_i+B_i\mathbb{E}$ of MoS₂ for all three energy regions. Here, the DOS within the gap region, $-\mathrm{\Delta }/2+{\lambda }_0<\mathbb{E}<\mathrm{\Delta }/2$, is zero.

The critical doping density which is required to populate the lower hole subband in MoS₂ is found to be

$$n_c=\frac{2}{\pi }\frac{{\lambda }_0\mathrm{\Delta }}{{\left(t_0{a}_0\right)}^2}=1.0\times {10}^{14}{\mathrm{cm}}^{-2}.$$E26

Therefore, for most experimentally accessible densities $n\le {10}^{13}$ cm${}^{-2}$, the lower hole subband still could not be populated at $T=0$.

Next, we would evaluate both sides of Eq. (16) for MoS₂. As an example, we consider electron doping with density $n_e>0$. The electron Fermi energy $E_F^e$ is determined by the following relation

$$n_e=\frac{c_1^{\left(3\right)}}{2}\left[{\left(E_F^e\right)}^2-\frac{{\mathrm{\Delta }}^2}{4}\right]+c_0^{\left(3\right)}\left(E_F^e-\frac{\mathrm{\Delta }}{2}\right).$$E27

From Eq. (27), we can easily find the electron Fermi energy $E_F>0$ as

$$E_F^e=\frac{1}{c_1^{\left(3\right)}}\left[-c_0^{\left(3\right)}+\sqrt{{\left(c_0^{\left(3\right)}+c_1^{\left(3\right)}\frac{\mathrm{\Delta }}{2}\right)}^2+2n_e{c}_1^{\left(3\right)}}\right].$$E28

In a similar way, for hole doping with density $n_h$ and the Fermi energy $E_F^h$ located between two hole subbands (region 2), we find

$$\mid E_F^h\mid =\frac{1}{c_1^{\left(2\right)}}\left\{\sqrt{{\left[c_0^{\left(2\right)}-c_1^{\left(2\right)}\left(\frac{\mathrm{\Delta }}{2}-{\lambda }_0\right)\right]}^2-2n_h{c}_1^{\left(2\right)}}-c_0^{\left(2\right)}\right\},$$E29

where $c_1^{\left(2\right)}<0$ and $\mathbb{E}<-\mathrm{\Delta }/2+{\lambda }_0$. From Eq. (29) we easily find the doping density

$$n_h=\left(\frac{\mathrm{\Delta }}{2}+E_F^h-{\lambda }_0\right)\left[c_0^{\left(2\right)}+\frac{c_1^{\left(2\right)}}{2}\left(\frac{\mathrm{\Delta }}{2}-\left(E_F^h+{\lambda }_0\right)\right)\right].$$E30

The right-hand side of Eq. (16) for TMDC’s could be expressed as a combination of electron and hole contributions ${\mathcal{I}}_e\left(\mathrm{\Delta }|T\right)-{\mathcal{I}}_h\left(\mathrm{\Delta },{\lambda }_0|T\right)$. Here, we will introduce two self-defined functions

$$\begin{array}{l}{\mathcal{A}}_0\left(\mathbb{E},T\right)={\left[1+exp\left(\frac{\mathbb{E}-\mu \left(T\right)}{k_{B}T}\right)\right]}^{-1},\\ {\mathcal{A}}_1\left(\mathbb{E},T\right)=\mathbb{E}{\mathcal{A}}_0\left(\mathbb{E},T\right)=\mathbb{E}{\left[1+exp\left(\frac{\mathbb{E}-\mu \left(T\right)}{k_{B}T}\right)\right]}^{-1},\end{array}$$E31

so that

$${\mathcal{I}}_e\left(\mathrm{\Delta }|T\right)=\sum _{j=0}^1{c}_j^{\left(3\right)}\underset{\mathrm{\Delta }/2}{\overset{\infty }{\int }}d\mathbb{E}{\mathcal{A}}_j\left(\mathbb{E},T\right).$$E32

For convenience, we introduce another function ${\mathcal{R}}_p\left(T,X\right)$

$${\mathcal{R}}_p\left(T,X\right)=\underset{0}{\overset{\infty }{\int }}\mathit{d\xi }\frac{{\xi }^p}{1+exp\left(\xi -X/k_{B}T\right)},$$E33

where $\xi =\left(\mathbb{E}-{\mathrm{\Delta }}_{<}\right)/k_{B}T$. Consequently, we are able to rewrite Eq. (32) as

$${\mathcal{I}}_e\left(\mathrm{\Delta }|T\right)=k_{B}T\left(c_0^{\left(3\right)}+\frac{\mathrm{\Delta }}{2}\right){\mathcal{R}}_0\left(T,\mu \left(T\right)-\frac{\mathrm{\Delta }}{2}\right)+c_1^{\left(3\right)}{\left(k_{B}T\right)}^2{\mathcal{R}}_1\left(T,\mu \left(T\right)-\frac{\mathrm{\Delta }}{2}\right),$$E34

where two terms with $p=0,1$ are physically related to a 2D electron gas. Explicitly, Eq. (33) leads to

$$\begin{array}{l}{\mathcal{R}}_0\left(T,X\right)=ln\left(1+exp\left(\frac{X}{k_{B}T}\right)\right),\\ {\mathcal{R}}_1\left(T,X\right)=-{\mathrm{Li}}_2\left(-exp\left(\frac{X}{k_{B}T}\right)\right),\end{array}$$E35

Using these self-defined functions and their notations, we finally arrive at the “hole term” ${\mathcal{I}}_h\left(\mathrm{\Delta },{\lambda }_0|T\right)$ in Eq. (16)

$$\begin{array}{ll}\underset{-\infty }{\overset{0}{\int }}d\mathbb{E}{\rho }_d\left(|\mathbb{E}|\right)\left[1-f_1\left(\mathbb{E},\mathbb{T}\right)\right]& =\underset{-\infty }{\overset{-\mathrm{\Delta }/2+{\lambda }_0}{\int }}d\mathbb{E}\left(c_0^{\left(2\right)}-c_1^{\left(2\right)}\mathbb{E}\right){\left[1+exp\left(\frac{\mu \left(T\right)-\mathbb{E}}{k_{B}T}\right)\right]}^{-1}\\ & +\underset{-\infty }{\overset{-\mathrm{\Delta }/2-{\lambda }_0}{\int }}d\mathbb{E}\left(\delta c_0^{\left(1\right)}-\delta c_1^{\left(1\right)}\mathbb{E}\right){\left[1+exp\left(\frac{\mu \left(T\right)-\mathbb{E}}{k_{B}T}\right)\right]}^{-1}\\ & \equiv \sum _{j=1}^4{\mathcal{I}}_h^{\left(j\right)}\left(\mathrm{\Delta },{\lambda }_0|T\right),\end{array}$$E36

where $\delta c_i^{\left(1\right)}=c_i^{\left(1\right)}-c_i^{\left(2\right)}$ for $i=0,1$, and

$$\begin{array}{l}{\mathcal{I}}_h^{\left(1\right)}\left(\mathrm{\Delta },{\lambda }_0|T\right)=k_{B}T\left(\frac{\mathrm{\Delta }}{2}-{\lambda }_0+c_2^{\left(0\right)}\right){\mathcal{R}}_0\left(T,-\left[\mu \left(T\right)+\mathrm{\Delta }/2-{\lambda }_0\right]\right),\\ {\mathcal{I}}_h^{\left(2\right)}\left(\mathrm{\Delta },{\lambda }_0|T\right)=c_2^{\left(1\right)}{\left(k_{B}T\right)}^2{\mathcal{R}}_1\left(T,-\left[\mu \left(T\right)+\mathrm{\Delta }/2-{\lambda }_0\right]\right),\\ {\mathcal{I}}_h^{\left(3\right)}\left(\mathrm{\Delta },{\lambda }_0|T\right)=k_{B}T\left(\frac{\mathrm{\Delta }}{2}+{\lambda }_0+\delta c_1^{\left(0\right)}\right){\mathcal{R}}_0\left(T,-\left[\mu \left(T\right)+\mathrm{\Delta }/2+{\lambda }_0\right]\right),\\ {\mathcal{I}}_h^{\left(4\right)}\left(\mathrm{\Delta },{\lambda }_0|T\right)=\delta c_1^{\left(1\right)}{\left(k_{B}T\right)}^2{\mathcal{R}}_1\left(T,-\left[\mu \left(T\right)+\mathrm{\Delta }/2+{\lambda }_0\right]\right).\end{array}$$E37

Here, both ${\mathcal{I}}_e\left(\mathrm{\Delta }|T\right)$ in Eq. (34) and ${\mathcal{I}}_h\left(\mathrm{\Delta },{\lambda }_0|T\right)$ in Eq. (36) comprise a finite-temperature part for the right-hand side of Eq. (16). Its left-hand side has already been given by Eq. (27). From these results, it is clear that there exists no symmetry between the electron and hole states at either zero or finite $T$. Finally, $\mu \left(T\right)$ of TMDC’s could be computed from a transcendental equation in Eq. (18), similarly to finding $\mu \left(T\right)$ for silicene.

By using the calculated $\mu \left(T\right)$, the plasmon dispersions and their Landau damping, determined from Eqs. (4) and (5), are displayed in Figure 5 for silicene at different $T$. Comparison of panels (a) and (b) indicates that the $T$ dependence of plasmon damping is not uniform even on a fixed convolution path $\mu \left(T\right)$. The doping density, on the other hand, widens the plasmon damping-free regions. Therefore, both the thermal and doping effects are found to compete with each other in dominating the plasmon dampings through selecting different convolution paths $\mu \left(T\right)$ with various doping densities or Fermi energies. Furthermore, the plasmon energy in (c) is pushed up slightly by increasing doping density at finite $T$.

Figure 5.

Particle-hole modes and plasmon branch for extrinsic (doped) silicene layer at a finite temperature. Panels (a) and (b) show two comparative graphs for $Im\left[{\mathrm{\Pi }}_0(q,\omega |\mu \left(T\right))\right]$ at zero and finite $T$, respectively, while plot (c) presents the finite-$T$ plasmon branch with $\mu \left(T\right)$ calculated from Eq. (16).

5. Dice lattice and $\alpha$-${\mathcal{T}}_3$ materials

In addition to graphene and silicene, another type of Dirac-cone materials is the one with fermionic states in which multiple Dirac points evolve into a middle flat band. One of the first fabricated materials with such a flat band is a dice or a ${\mathcal{T}}_3$ lattice, for which its atomic composition consists of hexagons similarly to graphene, but with an additional atom at the center of each hexagon. In a dice lattice, the bond coupling between a central site and three nearest neighbors is the same as that between atoms on corners, while for an $\alpha$-${\mathcal{T}}_3$ model the ratio $\alpha$ between hub-rim and rim-rim hopping coefficients can vary [45, 46] within the range of $0<\alpha <1$.

The low-energy electronic states of $\alpha$-${\mathcal{T}}_3$ materials are specified by a $3\times 3$ pseudospin-1 Dirac Hamiltonian, which results in three solutions for the energy dispersions and includes one completely flat and dispersionless band with ${\epsilon }_0\left(k\right)\equiv 0$. The other two bands are equivalent to Dirac cone ${\epsilon }_{\pm }\left(k\right)=\pm ħv_{F}k$ in graphene with the same Fermi velocity $v_F={10}^8$ cm/s. All of three bands touch at the corners of the first Brillouin zone, and therefore the band structure becomes metallic. In addition, the flat band has been shown to be stable against external perturbations, magnetic fields and structure disorders [47].

The $\alpha$-${\mathcal{T}}_3$ model was initially considered only as a theoretical contraption, an interpolation between graphene and a dice lattice. As parameter $\alpha \to 0$, this structure approaches graphene and a completely decoupled system of the hub atoms at the centers of each hexagon. A bit later, first evidence of really existing or fabricated materials with $\alpha$-${\mathcal{T}}_3$ electronic structure began mounting up. This includes Josephson arrays, optical arrangement based on the laser beams, Kagome and Lieb lattices with optical waveguides, Hg_1−xCd_xTe for a specific electron doping density, dielectric photonic crystals having zero-refractive index and a few others [48, 49]. So far, $\alpha$-${\mathcal{T}}_3$ model is believed to be the most promising innovative low-dimensional systems, and is one of the mostly investigated material in modern condensed matter physics. The most important technological application of $\alpha$-${\mathcal{T}}_3$ rests on the availability of materials with various $\alpha$ values, i.e., with small and large rim-hub hopping coefficients, ranging from $\alpha =0$ for graphene up to $\alpha =1$ for a dice lattice.

The low-energy Dirac-Weyl Hamiltonian for the $\alpha$-${\mathcal{T}}_3$ model is [45]

$${\widehat{\mathbb{H}}}_{\xi }^{\varphi }\left(\mathit{k}\right)=ħv_F\left[\begin{array}{ccc}0& k_{-}^{\xi }cos\varphi & 0\\ k_{+}^{\xi }cos\varphi & 0& k_{-}^{\xi }sin\varphi \\ 0& k_{+}^{\xi }sin\varphi & 0\end{array}\right],$$E38

where $\mathit{k}=\left(k_x,k_y\right)$ is the electron wave vector, $k_{\pm }^{\xi }=\xi k_x\pm {\mathit{ik}}_y$, $\xi =\pm 1$ corresponds to two different valleys, and $v_F$ is the Fermi velocity. Here, the parameter $\alpha$ is related to the geometry phase $\varphi ={tan}^{-1}\alpha$ which directly enters into the Hamiltonian in Eq. (38). The phase $\varphi$ possesses a fixed, one-to-one correspondence to the Berry phase of electrons in $\alpha$-${\mathcal{T}}_3$ model. In particular, for $\alpha =1$ or $\varphi =\pi /4$ we get a dice lattice with its Hamiltonian given by [50].

$${\widehat{\mathbb{H}}}_{\xi }^d\left(\mathit{k}\right)=\frac{ħv_F}{\sqrt{2}}\left[\begin{array}{ccc}0& k_{-}^{\xi }& 0\\ k_{+}^{\xi }& 0& k_{-}^{\xi }\\ 0& k_{+}^{\xi }& 0\end{array}\right].$$E39

Three energy bands from Hamiltonian in Eq. (38) or Eq. (39) are ${\epsilon }_0^{\gamma }\left(k\right)=\gamma ħv_{F}k$ for valence ($\gamma =-1$), conduction ($\gamma =+1$) and flat ($\gamma =0$) bands. These energy bands are degenerate with respect to $\xi$ and phase $\varphi$. The corresponding wave functions for the valence and conduction bands take the form

$${\mathrm{\Psi }}_0^{\gamma =\pm 1}\left(\mathit{k}|\xi ,\varphi \right)=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\xi cos\varphi {\mathrm{e}}^{-i{\mathit{\xi \theta }}_{\mathbf{k}}}\\ \gamma \\ \xi sin\varphi {\mathrm{e}}^{i{\mathit{\xi \theta }}_{\mathbf{k}}}\end{array}\right],$$E40

where $tan{\theta }_{\mathbf{k}}=k_y/k_x$. Meanwhile, for the flat band, we find

$${\mathrm{\Psi }}_0^{\gamma =0}\left(\mathit{k}|\xi ,\varphi \right)=\left[\begin{array}{c}\xi sin\varphi {\mathrm{e}}^{-i{\mathit{\xi \theta }}_{\mathbf{k}}}\\ 0\\ -\xi cos\varphi {\mathrm{e}}^{i{\mathit{\xi \theta }}_{\mathbf{k}}}\end{array}\right].$$E41

Here, the components of wave functions in Eqs. (40) and (41) depend on valley index $\xi$ and phase $\varphi$, which leads to the same dependence on all collective properties of an $\alpha$-${\mathcal{T}}_3$ materials, including plasmon dispersion.

Now we turn to deriving plasmon branches and their damping rates at finite $T$ in $\alpha$-${\mathcal{T}}_3$ model. The computation procedure is quite similar to that in the case of two non-equivalent doubly degenerate subband pairs, including silicene, germanene and MoS₂ discussed in Section 4.

For $\alpha$-${\mathcal{T}}_3$ model, the finite-$T$ polarization function ${\mathrm{\Pi }}_T\left(q,\omega |\mu \left(T\right)\right)$ can be obtained by an integral transformation of its zero-temperature counterpart ${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)$, as presented in Eq. (5). In this case, the zero-$T$ counterpart ${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)$ is calculated as

$${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)=\frac{1}{{\pi }^2}\sum _{\gamma ,{\gamma }^{\prime }=0,\pm 1}\int d^2\mathit{k}{\mathbb{O}}_{\gamma ,{\gamma }^{\prime }}\left(\mathit{k},\mathit{k}+\mathit{q}|\varphi \right)\times \frac{{\mathrm{\Theta }}_0\left(E_F-{\epsilon }_{\gamma }\left(k\right)\right)-{\mathrm{\Theta }}_0\left(E_F-{\epsilon }_{{\gamma }^{\prime }}\left(|\mathit{k}+\mathit{q}|\right)\right)}{ħ\left(\omega +i{0}^{+}\right)+{\epsilon }_{\gamma }\left(k\right)-{\epsilon }_{{\gamma }^{\prime }}\left(|\mathit{k}+\mathit{q}|\right)}.$$E42

Structurally, Eq. (42) looks quite similarly to Eq. (6) for buckled honeycomb lattices and TMDC’s. The most significant difference comes as the existence of an additional flat band with $\gamma =0$ so that the summation index runs over $\pm 1$ and $0$ instead of two. On the other hand, the overall expression for ${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)$ in Eq. (42) is simplified because the $4$-fold degeneracy of each energy band independent of valley and spin index.

Here, we would limit our consideration to the case of electron doping with $n>0$ and apply the random-phase approximation theory only for that case. For electron doping with $n>0$, we can neglect the transitions within the valence band and also the transitions between the flat and valences bands due to full occupations of these electronic states. On the other hand, the overlap of initial and final electron transition states is defined by [51] ${\mathbb{O}}_{\gamma ,{\gamma }^{\prime }}^{\xi }\left(\mathit{k},{\mathit{k}}^{\prime }|\varphi ,{\lambda }_0\right)$ with respect to the initial ${\mathrm{\Psi }}_{\gamma }^{\xi }\left(\mathit{k},{\lambda }_0\right)$ and the final ${\mathrm{\Psi }}_{{\gamma }^{\prime }}^{\xi }\left({\mathit{k}}^{\prime },{\lambda }_0\right)$ states with a momentum transfer $\mathit{q}={\mathit{k}}^{\prime }-\mathit{k}$, i.e.,

$$\begin{array}{l}{\mathbb{O}}_{\gamma ,{\gamma }^{\prime }}^{\xi }\left(\mathit{k},\mathit{k}+\mathit{q}|\varphi ,{\lambda }_0\right)=\mid {\mathbb{S}}_{\gamma ,{\gamma }^{\prime }}^{\xi }{\left(\mathit{k},\mathit{k}+\mathit{q}\left|\varphi ,{\lambda }_0\right|\right)}^2,\\ {\mathbb{S}}_{\gamma ,{\gamma }^{\prime }}^{\xi }\left(\mathit{k},\mathit{k}+\mathit{q}|\varphi ,{\lambda }_0\right)=⟨{\mathrm{\Psi }}_{\gamma }^{\xi }\left(\mathit{k},{\lambda }_0\right)|{\mathrm{\Psi }}_{{\gamma }^{\prime }}^{\xi }(\mathit{k}+\mathit{q},{\lambda }_0)⟩,\end{array}$$E43

where ${\beta }_{\mathbf{k},{\mathbf{k}}^{\prime }}={\theta }_{{\mathbf{k}}^{\prime }}-{\theta }_{\mathbf{k}}$ is the scattering angle between two electronic states and $k^{\prime }=\sqrt{k^2+q^2+2\mathit{kq}cos{\beta }_{\mathbf{k},{\mathbf{k}}^{\prime }}}$. Moreover, we find from Eq. (43) [52]

$${\mathbb{O}}_{\gamma ,{\gamma }^{\prime }}^{\xi }\left(\mathit{k},\mathit{k}+\mathit{q}|\varphi ,{\lambda }_0\right)=\frac{1}{4}\left[{\left(1+cos{\beta }_{\mathbf{k},{\mathbf{k}}^{\prime }}\right)}^2+{cos}^2\left(2\varphi \right){sin}^2{\beta }_{\mathbf{k},{\mathbf{k}}^{\prime }}\right]$$E44

for an arbitrary value of $\varphi$ or $\alpha$. It is easy to verify the known results $\left(1+cos{\beta }_{\mathbf{k},{\mathbf{k}}^{\prime }}\right)/2$ for graphene and ${\left(1+cos{\beta }_{\mathbf{k},{\mathbf{k}}^{\prime }}\right)}^2/4$ for a dice lattice as two limiting cases of our general result in Eq. (44) as $\alpha \to 0$ or $\alpha \to 1$, respectively. Furthermore, we find from Eq. (44) that the overlap does not depend on valley index $\xi$, even though individual wave function does, and then this index can be dropped. However, the valley-dependence in ${\mathbb{O}}_{\gamma ,{\gamma }^{\prime }}^{\xi }\left(\mathit{k},\mathit{k}+\mathit{q}|\varphi ,{\lambda }_0\right)$ persists if $\alpha$-${\mathcal{T}}_3$ material is irradiated by circularly- or elliptically-polarized light. This incident radiation permits creating an valleytronic filter or any other types of valleytronic electron device.

Density plots for Landau damping with $Im\left[{\mathrm{\Pi }}_O(q,\omega |\mu \left(T\right))\right]\ne 0$ is presented in Figure 6, where we find plasmon branch will be completely free from damping within the region determined by $ħ\omega /E_0\le 1$ and $q/k_0\le 1$, independent of geometry phase $\varphi$. On the other hand, another region with $ħ\omega \le ħv_{F}q$ (below the diagonal) becomes always Landau damped. Increasing $T$ is able to increase greatly the damping in the region below the diagonal, as seen in Figure 6(c).

Figure 6.

Particle-hole modes, determined by non-zero $Im\left[{\mathrm{\Pi }}_O(q,\omega |\mu \left(T\right))\right]$ within the $q$-$\omega$ plane, for an $\alpha$-${\mathcal{T}}_3$ layer with $\varphi =\pi /10$ (in (a), (c)) and $\varphi =\pi /7$ in (b). Panels (a) and (b) correspond to $T=0$, while plot (c) is for $k_{B}T=1.0E_F$.

In a correspondence to the damping of plasmons presented in Figure 6, we show in Figure 7 the density plots for plasmon dispersions at $T=0$ in (a), (b) and $k_{B}T=E_F$ in (c), (d). Comparing Figure 7(a) with Figure 7(c) we have clearly seen the thermal suppression of Landau damping for plasmon mode entering into a high-frequency region beyond $ħ\omega =E_F$. To visualize a full plasmon dispersion clearly, we also include damped counterpart in Figure 7(b) and (d) at $T=0$ and $k_{B}T=E_F$, respectively, where a significant enhancement of plasmon energy appears for large $q$ values, moving upwards from the diagonal.

Figure 7.

Plasmon branches for an isolated $\alpha$-${\mathcal{T}}_3$ layer with $\varphi =\pi /10$. Panels (a) and (c) only show undamped plasmons, while (b) and (d) display full plasmon branches including damped ones. Left panels (a) and (b) corresponds to $T=0$, while $k_{B}T=1.0E_F$ for right panels (c) and (d).

6. Plasmons in $\alpha$-${\mathcal{T}}_3$ layer coupled to conducting substrate

In the last part of THIS CHAPTER, WE WOULD LIKE TO FOCUS ON finite-$T$ plasmons in a so-called nanoscale-hybrid structure consisting of a 2D layer, such as, graphene, silicene or a dice lattice, which is Coulomb-coupled to a large, conducting material. Physically, the Coulomb coupling between the 2D layer and the conductor results in a strong hybridization of graphene plasmon and localized surface-plasmon modes. This structure, which is referred to as an open system, could be realized experimentally or even by a device fabrication.

Our schematics for an open system is shown in Figure 8. The dynamical screening to the Coulomb interaction between electrons in a 2D layer and in metallic substrate is taken into account by a nonlocal and dynamical inverse dielectric function $\mathcal{K}\left(\mathit{r},{\mathit{r}}^{\prime };\omega \right)$, as demonstrated in Refs. [53, 54]. This nonlocal inverse dielectric function is connected to a dielectric function $\mathrm{ϵ}\left(\mathit{r},{\mathit{r}}^{\prime };\omega \right)$ in Eq. (4) by

Figure 8.

Schematics for a silicene-based open system and numerical results for the two plasmon branches and their damping in this system with ${\mathrm{\Delta }}_{\mathit{SO}}=0.3E_0$ and $0.1E_0$, where $E_0=ħv_F\sqrt{\pi n_0}=54.6$ meV.

$$\int d^3{\mathit{r}}^{\prime }\mathcal{K}\left(\mathit{r},{\mathit{r}}^{\prime };\omega \right)\mathrm{ϵ}\left({\mathit{r}}^{\prime },{\mathit{r}}^{″};\omega \right)=\delta \left(\mathit{r}-{\mathit{r}}^{″}\right),$$E45

and the resonances in $\mathcal{K}\left(\mathit{r},{\mathit{r}}^{\prime };\omega \right)$ reveal the nonlocal hybridized plasmon modes supported by both 2D layer and the conducting surface as a single quantum system.

By using the Drude model for metallic substrate, the dielectric function can be written as ${\mathrm{ϵ}}_B\left(\omega \right)=1-{\mathrm{\Omega }}_p^2/{\omega }^2$, where ${\mathrm{\Omega }}_p=\sqrt{n_0{e}^2/{\mathrm{ϵ}}_0{\mathrm{ϵ}}_b{m}^{\ast }}$ is the bulk-plasma frequency for the conductor, $n_0$ electron concentration and $m^{\ast }$ is the effective mass of electrons. Drude model describes electron screening in the long-wavelength limit. Based on the previously developed mean-field theory [53, 55, 56], we are able to calculate plasmon dispersions in this 2D open system. For this, the plasmon dispersions are obtained from the zeros of the so-called dispersion factor${\mathbb{S}}_C\left(q,\omega \right)$, instead of the dielectric function in Eq. (4). ${\mathbb{S}}_C\left(q,\omega \right)$ for this open system is given by [25, 54, 57]

$${\mathbb{S}}_C\left(q,\omega |E_F\right)=1-\frac{2{\mathit{\pi \alpha }}_r}{q}{\mathrm{\Pi }}_0\left(q,\omega |E_F\right)\left[1+\frac{{\mathrm{\Omega }}_p^2}{2{\omega }^2-{\mathrm{\Omega }}_p^2}exp\left(-2\mathit{qa}\right)\right],$$E46

where $a$ is the separation between the 2D layer and the conducting surface. Most important, we should emphasize that the second term in Eq. (46) does not have a full analogy with polarization function of an isolated layer, and the resulting plasmon dispersions in open system represents a hybridized plasmon mode with the environment. Therefore, these plasmon dispersions are expected to be sensitive to Coulomb coupling to electrons in the conducting substrate through a factor $\sim exp\left(-2\mathit{qa}\right)$ in Eq. (46), similarly to what we have found for coupled double graphene layers [16]. The strong Coulomb coupling leads to a linear dispersion of plasmon in this open system [54, 58], which is in contrast with well-known $\sim \sqrt{q}$ dependence in all 2D materials.

As a special example, let us consider a silicene 2D layer with two bandgaps ${\mathrm{\Delta }}_{<,>}$ and an electron doping density $n$. We start with seeking for a non-interacting polarization function in the long-wave limit $q\ll k_F^{<,>}$ for doping density $n$ and assume a high-enough $n$ to keep the Fermi level $E_F=\sqrt{{\left(ħv_F{k}_F^{\beta }\right)}^2+{\mathrm{\Delta }}_{\beta }^2}>{\mathrm{\Delta }}_{>}$ above the large bandgap. Under this assumption, we get

$${\mathrm{\Pi }}_0\left(q,\omega |E_F\right)=\frac{q^2}{\pi {ħ}^2{\omega }^2}\sum _{\beta =>,<}{k}_F^{\beta }{\left|\frac{\partial {\mathbb{E}}_{\beta }\left(k\right)}{\partial k}\right|}_{k=k_F^{\beta }}=\frac{E_F}{\pi }\left(2-\frac{{\mathrm{\Delta }}_{<}^2}{E_F^2}-\frac{{\mathrm{\Delta }}_{>}^2}{E_F^2}\right)\frac{q^2}{{ħ}^2{\omega }^2},$$E47

where $k_F^{\beta }=\sqrt{2\pi n_{\beta }}$ are two different Fermi wave numbers associated with a single Fermi energy $E_F$, and $n_{\beta }$ is the electron density for each subband satisfying $n=n_{<}+n_{>}$.

In the limit of $a\to \infty$, the plasmon branch of an isolated silicene layer can be recovered from Eq. (46), yielding

$${\omega }_p^2\left(q\right)=\frac{4{\alpha }_r}{{ħ}^2{E}_F}\left(E_F^2-\frac{{\mathrm{\Delta }}_{>}^2+{\mathrm{\Delta }}_{<}^2}{2}\right)q\equiv \mathrm{\Xi }q,$$E48

where for convenience we introduced a coefficient

$$\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)=\frac{4{\alpha }_r}{{ħ}^2{E}_F}\left(E_F^2-\frac{{\mathrm{\Delta }}_{>}^2+{\mathrm{\Delta }}_{<}^2}{2}\right).$$E49

We notice from Eq. (48) that ${\omega }_p\left(q\right)\sim \sqrt{q}$, disregarding of the energy bandgaps ${\mathrm{\Delta }}_{<,>}$ or doping density $n$. On the other hand, the Fermi energy $E_F$ for silicene is given by Eq. (17).

Furthermore, using the notation defined by Eq. (49), we get from Eqs. (46) and (47) that

$$1-\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)\frac{q}{{\omega }^2}\left[1+\frac{{\mathrm{\Omega }}_p^2}{2{\omega }^2-{\mathrm{\Omega }}_p^2}exp\left(-2\mathit{qa}\right)\right]=0,$$E50

which leads to a bi-quadratic equation

$$2{\left(\frac{{\omega }^2}{{\mathrm{\Omega }}_p^2}\right)}^2-\left[1+\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)\frac{2q}{{\mathrm{\Omega }}_p^2}\right]\left(\frac{{\omega }^2}{{\mathrm{\Omega }}_p^2}\right)+\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)\frac{q}{{\mathrm{\Omega }}_p^2}\left[1-exp\left(-2\mathit{qa}\right)\right]=0.$$E51

Eq. (51) gives rise to two solutions

$$\frac{4}{{\mathrm{\Omega }}_p^2}{\omega }_{p,\pm }^2=\left(1+\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)\frac{2q}{{\mathrm{\Omega }}_p^2}\right)\pm \sqrt{{\left(1+\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)\frac{2q}{{\mathrm{\Omega }}_p^2}\right)}^2-\frac{8q\mathrm{\Xi }\left(E_F,{\mathrm{\Delta }}_{\beta }\right)}{{\mathrm{\Omega }}_p^2}\left[1-exp\left(-2\mathit{qa}\right)\right]},$$E52

where $\pm$ terms correspond to in-phase and out-of-phase plasmon modes, respectively. Two hybrid plasmon modes in Eq. (52) become

$$\begin{array}{l}{\omega }_{p,+}\left(q\right)\simeq \frac{{\mathrm{\Omega }}_p}{\sqrt{2}}+\frac{\mathrm{\Xi }}{\sqrt{2}{\mathrm{\Omega }}_p}q-\frac{\mathrm{\Xi }\left(\mathrm{\Xi }+4a{\mathrm{\Omega }}_p^2\right)}{2\sqrt{2}{\mathrm{\Omega }}_p^3}{q}^2+\mathcal{O}\left(q^3\right),\\ {\omega }_{p,-}\left(q\right)\simeq q\sqrt{2a\mathrm{\Xi }}-\frac{\mathrm{\Xi }\sqrt{2a\mathrm{\Xi }}}{{\mathrm{\Omega }}_p^2}{q}^2+\mathcal{O}\left(q^3\right).\end{array}$$E53

In Eq. (53), both plasmon branches contain a linear $\sim q$ term, and ${\omega }_{p,+}\left(q\right)$ approaches a constant as $q\to 0$, i.e., an optical mode for plasmons. Two independent bandgaps, ${\mathrm{\Delta }}_{\mathit{SO}}$ and ${\mathrm{\Delta }}_z$, together with doping density $n$, play a crucial role on shaping the plasmon dispersions, as well as the particle-hole mode damping regions. The outer boundaries of a particle-hole mode region specify an area within the $q$-$\omega$ plane in which the plasmon modes become damping free and are solely determined by ${\mathrm{\Delta }}_{<}$, while the group velocity of plasmon mode depends on both ${\mathrm{\Delta }}_{<}$ and ${\mathrm{\Delta }}_{>}$. Since each bandgap could be experimentally tuned by applying a perpendicular electric field, we acquire a full control of both plasmon dispersions and their damping-free regions at the same time.

Numerical results for thermal plasmons in open system are presented in Figure 9. Similarly to what we have found for graphene and silicene, there are two plasmon branches, both of which depend linearly on $q$ with a finite slope as $q\ll k_0$. The acoustic-plasmon branch starts from the origin, while the optical-plasmon branch from ${\mathrm{\Omega }}_p/\sqrt{2}$. The dispersion of each branch also varies with parameter $\alpha$ (or $\varphi ={tan}^{-1}\alpha$), which is observed for the upper branch, as shown in (a) and (c) of Figure 9. In addition, we see a much smaller slope for the lower plasmon branch in Figure 9(b) and (d) due to enhanced Coulomb coupling with a reduced separation $a$. The finite-temperature upper plasmon branches in Figure 9(e) and (f) are shifted up greatly, as it is expected to be true for all finite-temperature plasmons, which is further accompanied by enhanced damping below the diagonal as seen in Figure 6(c). Meanwhile, the lower plasmon branch seems much less affected by finite temperatures, as demonstrated by both upper and lower rows of Figure 9 for different separations, except for enhanced damping in Figure 9(f) below the surface-plasmon energy $ħ{\mathrm{\Omega }}_p/\sqrt{2}$.

Figure 9.

Nonlocal hybridized plasmon dispersions for $\alpha$-${\mathcal{T}}_3$ layer coupled to a closely-located surface of a semi-infinite conductor. Panels (a)–(d) are for $T=0$, while plots (e) and (f) for $k_{B}T=E_0$. All the upper-row plots correspond to the separation $a=1.0k_F^{-1}$, and the lower-row ones to $a=0.5k_F^{-1}$. Additionally, middle-column plots, (b) and (d), correspond to $\varphi =\pi /7$, and all other columns to $\varphi =\pi /10$.

7. Summary and remarks

In conclusion, we have developed a general theory for finite-temperature polarization function, plasmon dispersions and their damping for all known innovation 2D Dirac-cone materials with various types of symmetries and bandgaps. We have also derived a set of explicit transcendental equations determining the chemical potential as a function of temperature, which serves as a key part in calculating finite-temperature polarization function through the so-called thermal convolution path. The selection of a particular path with a specific $\mu \left(T\right)$ could be employed for studying the temperature dependence of plasmon modes in each of the considered 2D materials. The fact that a chemical potential keeps its sign is true only for materials with symmetric energy bands of electrons and holes, but can cross the zero line for TMDC’s with asymmetric electron and hole bands.

Using the calculated finite-temperature polarization function, we have further found the dispersions of hybrid plasmon-modes in various types of open systems including a 2D material coupled to a conducting substrate. The obtained plasmon dispersions in these 2D-layer systems are crucial for measuring spin-orbit interaction strength and dynamical screening to Coulomb interaction between electrons in 2D materials, as well as for designing novel surface-plasmon based multi-functional near-field opto-electronic devices.

We have generalized our developed theory for 2D materials further to most recently proposed $\alpha$-${\mathcal{T}}_3$ lattices, in which the characteristic parameter $\alpha$ is the ratio of hub-rim to hub-hub hopping coefficients and can vary from $0$ to $1$ continuously corresponding to different material properties. For $\alpha$-${\mathcal{T}}_3$ materials, we have obtained the hybrid plasmon modes for different $\alpha$ values at both zero and finite temperatures and demonstrated that the resulting hybridized plasmon dispersions could be tuned sensitively by geometry phase, temperature, and separation between $\alpha$-${\mathcal{T}}_3$ layer and conducting surface. Such tunability has a profound influence on performance of $\alpha$-${\mathcal{T}}_3$ material based quantum electronic devices.

Acknowledgments

A.I. thanks Liubov Zhemchuzhna for helpful and fruitful discussions, and Drs. Armando Howard, Leon Johnson and Ms. Beverly Tarver for proofreading the manuscript and providing very useful suggestions on the style and language. G.G. would like to acknowledge the financial support from the Air Force Research Laboratory (AFRL) through grant FA9453-18-1-0100 and award FA2386-18-1-0120. D.H. thanks the supports from the Laboratory University Collaboration Initiative (LUCI) program and from the Air Force Office of Scientific Research (AFOSR).

Conflict of interest

All the authors declare that they have no conflict of interest.