Plasmon modes of bilayer molybdenum disulfide: A density functional study

TL;DR

This study uses density functional theory and the random phase approximation to investigate plasmon modes in bilayer MoS2, revealing three distinct collective electronic excitations sensitive to doping levels.

Contribution

It provides a detailed ab initio analysis of plasmon modes in bilayer MoS2, including the prediction of three specific plasmon excitations and their dependence on doping.

Findings

Identification of three plasmon modes in bilayer MoS2

Doping significantly affects the electron energy-loss spectra

Prediction of a $ oot q$ dispersion for the optical mode

Abstract

We explore the collective electronic excitations of bilayer molybdenum disulfide (MoS$_2$) using the density functional theory together with the random phase approximation. The many-body dielectric function and electron energy-loss spectra are calculated using an {\it ab initio} based model involving material-realistic physical properties. The electron energy-loss function of bilayer MoS$_2$ system is found to be sensitive to either electron or hole doping and it is owing to the fact that the Kohn-Sham band dispersions are not symmetric for energies above and below the zero Fermi level. Three plasmon modes are predicted. A damped high-energy mode, one optical mode (in-phase mode) for which the plasmon dispersion exhibits $\sqrt q $ in the long wavelength limit originating from low-energy electron scattering and finally a highly damped acoustic mode (out-of-phase mode).

Tables

Table 1. Table 1: Calculated the hole and electron effective masses ( m h ∗ superscript subscript 𝑚 ℎ m_{h}^{*} , m e ∗ superscript subscript 𝑚 𝑒 m_{e}^{*} ) in units of the electron bare mass at the high symmetry points. The effective masses of MoS 2 are derived from the band structures as shown in Fig. 1

		Effective Masses
Structure	point	$m_e^{\ast }$	$m_h^{\ast }$
Monolayer	K	0.47	0.61
	Q	0.58
	$\Gamma$		3.5
Bilayer	K	0.53	0.63
	Q	0.56
	$\Gamma$		1.07

Equations

$${\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}^0(q,\omega )=\frac{2}{\Omega }\sum _{\mathbf{k},n,m}\frac{f_n(\mathbf{k})-f_m(\mathbf{k}+\mathbf{q})}{\hslash \omega +i\eta +{\epsilon }_n(\mathbf{k})-{\epsilon }_m(\mathbf{k}+\mathbf{q})}$$

$$M_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\mathbf{G})M_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}^{\ast }({\mathbf{G}}^{\prime })$$

$${\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}={\chi }_{\mathbf{G}{\mathbf{G}}^{\prime }}^0+\sum _{{\mathbf{G}}_1{\mathbf{G}}_2}{\chi }_{\mathbf{G}{\mathbf{G}}_1}^0{\nu }_{{\mathbf{G}}_1{\mathbf{G}}_2}{\chi }_{{\mathbf{G}}_2{\mathbf{G}}^{\prime }}$$

$$M_{n\mathbf{k},m\mathbf{k}+\mathbf{q}}(\mathbf{G})=<{\Phi }_{n\mathbf{k}}|e^{-i(\mathbf{q}+\mathbf{G})\cdot \mathbf{r}}|{\Phi }_{m\mathbf{k}+\mathbf{q}}{>}_{\Omega }$$

$${\Phi }_{n\mathbf{k}}(\mathbf{r})=\frac{1}{\Omega }\sum _{\mathbf{G}}{C}_{n\mathbf{k}}(\mathbf{G})e^{i(\mathbf{k}+\mathbf{G})\cdot \mathbf{r}}$$

$${\nu }_{\mathbf{G}{\mathbf{G}}^{\prime }}^0=\frac{2\pi e^2{\delta }_{\mathbf{g}{\mathbf{g}}^{\prime }}}{|\mathbf{q}+\mathbf{g}|}{\int }_{-L/2}^{L/2}dz{\int }_{-L/2}^{L/2}d{z}^{\prime }{e}^{i(G_{z^{\prime }}{z}^{\prime }-G_{z}z)-|\mathbf{q}+\mathbf{g}||z-z^{\prime }|}$$

$${ϵ}_{\mathbf{G}{\mathbf{G}}^{\prime }}={\delta }_{\mathbf{G}{\mathbf{G}}^{\prime }}-\sum _{{\mathbf{G}}_1}{\nu }_{\mathbf{G}{\mathbf{G}}_1}^0{\chi }_{{\mathbf{G}}_1{\mathbf{G}}^{\prime }}^0$$

$$ϵ(q,\omega )=\frac{1}{({ϵ}^{-1}{)}_{\mathbf{G}{\mathbf{G}}^{\prime }}}{|}_{\mathbf{G}={\mathbf{G}}^{\prime }=0}$$

$$E_{EEL}(q,\omega )=-\Im m[1/ϵ(q,\omega )]$$

Full text

I Introduction

Despite being the most promising two-dimensional (2D) material, gapless graphene has limitations of its applications in nanoelectronics and nanophotonics. This leads to finding of other 2D-materials with finite band gap such as transition metal dichalcogenides systems (TMDCs) wang ; TMD-DF2 . Belonging to the family of layered TMCD, molybdenum disulfide (MoS2) has been widely used in numerous areas, such as hydrodesulfurization catalyst, photovoltaic cell, photocatalyst, nanotribology, lithium battery, and dry lubrication, due to their distinctive electronic, optical, and catalytic properties mos2 . In monolayer MoS2, a strong photoluminescence peak at about $1.90$ eV, together with peaks at about $1.90$ and $2.05$ eV of the adsorption spectrum, indicates that MoS2 undergoes an indirect to direct band gap transition when its bulk or multilayers is replaced by a monolayer DFT-GW . Therefore, MoS2 becomes a very interesting material owing to unique electronic and optical properties vasp ; ref_4-vasp ; ref_5-vasp ; ref_6-vasp .

The collective density oscillations of doped two-dimensional crystalline systems have recently received attention plasmon ; asgari ; wang2 . Interestingly, the propagation of graphene plasmons has been directly imaged in real space by using scattering-type scanning near-field optical microscopy fei in which the wavelength of the plasmon is much smaller than the free-space excitation wavelength, allowing an extreme concentration of electromagnetic energy. This allows that graphene plasmon properties can easily be gate tuned.

Moreover, the energy of the metallic plasmon mode is restricted by the achievable carrier concentration and this limits the application of graphene plasmonics to the tetrahertz regime. Metallic TMDCs have much higher charge carrier densities stemming from a quite flat band structure leading to plasmon energies of around 1 eV in bulk TMDCs ref_14-An ; ref_15-An . Control and manipulation of the plasmon properties of these two-dimensional systems are expected to guide the design of next-generation nanophotonic and nanoelectronic devices, with the enhanced capability to operate from terahertz (THz) to infrared frequencies Si .

Recently the angle-resolved optical absorption and transmission in MoS2 are studied by using current-current response tensor within the random phase approximation (RPA) Despoja . Also, the optical properties of MoS2 which include quasiparticle GW corrections, solving the Bethe-Salpeter equation, have been calculated GW1 ; GW2 . Their results are in very good agreement with those measured data in the experiment. The plasmon dispersions of monolayer MoS2 and other TMDCs have been already calculated by Scholz et al., Scholz using a model Hamiltonian and Andersen et al., Andersen within ab initio simulations. While most of the aforementioned theoretical works have been focused on monolayer MoS2, the plamon dispersion of bilayer MoS2 has not been investigated yet.

In this paper, we have been using a recently proposed theoretical formulation MLG ; BLG based on ab initio density functional theory (DFT) together with the RPA to investigate the electronic excitation spectrum of bilayer MoS2. For this purpose, the electronic ground-state of a periodically repeated slab is first determined and then a Dyson-like equation is solved within the RPA to calculate the density-density response function. Furthermore, a two-dimensional correction is applied to eliminate the artificial interaction between the replicas.

The dynamical density-density response function of monolayer and bilayer MoS2 are calculated within the DFT and RPA theory. Having known the density-density response function, we therefore can calculate the macroscopic dielectric function whose imaginary part gives the optical absorption spectrum and the collective modes are established by the zero in the real part of the macroscopic dielectric function. Besides, not only the dielectric function is relevant for plasmon modes, but also for transport and the phonon spectra phonon is a useful quantity. In addition, the dielectric function also provides the connection between theory and experiment. The theoretical dielectric function is related to the electron energy-loss function, and it provides useful information about the optical properties of the system. Here, we are just interested in the low-energy excitation investigating the collective modes of the intraband excitations.

This paper is organized as follows. In Sec. II, we present the methodology used to calculate the ground-state of the system and also describe a full ab initio approach to calculate the plasmon modes, properties of two-dimensional systems, based on density functional theory, within the random phase approximation. In Sec. III we present and describe collective mode results of bilayer MoS2. Finally, we conclude and summarize our main results in Sec. IV

II THEORY AND COMPUTATIONAL METHODS

II.1 DFT calculations

In order to calculate the electronic structure, we use plane-wave basis in the local density approximation (LDA) and norm conserving pseudopotential within the Quantum Espresso (QE) packageQE . The ground-state calculations are carried on unshifted $60\times 60\times 1$ Monkhorst-Pack (MP) K-point mesh of the first Brillouin zone (BZ) and with a $50$ Ry plane-wave cutoff. In the case of a layered material, we implement the van der Waals interaction. The nearest-neighbor in-plane and inter-plane distances are optimized, respectively, to the values $a=3.15$ and $d=5.89$Åas defined in Fig. 1, however, for the unit cell constant, we use the experimental value of $a=3.16$Å aexp . To avoid the effects of the interaction between images in the periodic unit cell calculations, the lattice parameter in the direction perpendicular to one single bilayer plane is chosen to be $L=20$Å. We find the $AA^{\prime}$ stacking with Mo over S, as the most stable stacking order of bilayer MoS2, which is the most studied in the previous calculation AA' .

II.2 Density-density response function

A central quantity in the theoretical formulation of the many-body effects in electron liquid is the noninteracting dynamical response function. The noninteractiong density-density response function of a three-dimensional periodic of electrons in the reciprocal space is given by

[TABLE]

which is a consequence of the Kubo formula for periodic the system kubo .

Here, the factor of 2 accounts for the spin degeneracy and $f_{n}(\mathbf{k})=\theta(\varepsilon_{\rm F}-\varepsilon_{n}(\mathbf{k}))$ is the Fermi-Dirac distribution of the charge carrier with energy $\varepsilon_{n}(k)$ at $T=0$. In the theory, the linear combination of plan-waves are used to determine the Kohn-Sham (KS) single-particle orbitals of the DFT. The KS wavefunctions are normalized to unity in the crystal volume $\Omega$. The broadening parameter $\eta$ used in this calculation is $0.02$ eV. In the summation over k, we use $151\times 151\times 1$ K-point mesh sampling in the BZ and the $n$ and $m$ sum over $20$ bands for monolayer and $30$ bands for bilayer MoS2.

The density-density response function can be obtained in the framework of DFT, as follows DFT :

[TABLE]

where $\nu_{\mathbf{G}\mathbf{G}^{\prime}}$ represent the Fourier coefficients of an effective electron-electron interaction. For the electron liquid the bare Coulomb interaction is given by $\nu_{\mathbf{G}\mathbf{G}^{\prime}}^{0}=4\pi e^{2}\delta_{\mathbf{G}\mathbf{G}^{\prime}}/|\mathbf{q}+\mathbf{G}|^{2}$. The RPA procedure, an approximation valid in the high-density limit, takes into account electron interaction only to the extent required to produce the screening field and thus the response to the screened field is measured by $\chi^{0}$.

The matrix elements of Eq. (1) have the form

[TABLE]

where $\mathbf{q}$ is the momentum transfer vector parallel to the $x-y$ plane. Wave functions $\Phi_{n\mathbf{k}}(\mathbf{r})$ are the KS electron wave functions expanded in the plane-wave basis have the form

[TABLE]

where the coefficients $C_{n\mathbf{k}}(\mathbf{G})$ are obtained by solving the LDA-KS equations self-consistently.

This approach is correct for the purely three-dimensional periodic system. The long-range behavior of the Coulomb interaction allows non-negligible interactions between repeated planar arrays even at large distance. This unphysical phenomenon can be removed by replacing $\nu_{\mathbf{G}\mathbf{G}^{\prime}}$ by the truncated Fourier integral over the cut of plane axis $(z)$ MLG ; BLG ; NRG and thus we have

[TABLE]

where the $\mathbf{g}$ and $G_{z}$ denote the in-plane and out-plane components of $\mathbf{G}$ and we assume that $q$ is never be zero owing to a uniform background of positive charge.

In the framework of linear response theory, the inelastic cross section corresponding to a process where the external perturbation creates an excitation of energy $\hbar\omega$ and wavevector $\mathbf{q}+\mathbf{G}$ is related to the diagonal elements of the dielectric function in the level of the RPA

[TABLE]

and the plasmon modes are established by the zero in the real part of the macroscopic dielectric function given by

[TABLE]

as long as there is no damping.

Electron energy-loss spectroscopy (EELS), on the other hand, is an analytical technique which is based on inelastic scattering of fast electrons in a thin sample. EELS offers unique possibilities for advanced materials analysis owing to the broad range of inelastic interactions of the high energy electrons with the specimen atoms, ranging from phonon interactions to ionisation processes. The low-loss or valence region of an electron energy-loss function (EEL) spectrum provides similar information to that provided by optical spectroscopy, containing valuable information about the band structure and in particular about the dielectric properties of a material e.g., band gap, surface plasmon modes. Interestingly, the most prominent peak comes from a plasma resonance of the valence atoms. The EEL is proportional to the imaginary part of the inverse dielectric function which is given by

[TABLE]

It is worth mentioning that the nonlocal field effects are included in EEL through the off-diagonal elements of the general $\chi_{\mathbf{G}\mathbf{G}^{\prime}}$nlc function.

III RESULTS AND DISCUSSION

In this section, we present our main numerical results based on first-principles simulations. Our aim is to explore the density-density response function and plasmon modes of bilayer MoS2. All the first-principles calculations are performed at zero temperature without considering the spin-orbit interaction. Our study is devoted to a consideration of the homogenous electron liquid in both the high-density and the intermediate-density regime. We do not consider the low-density charge carrier, primarily because its behavior is not especially relevant to the properties of metals.

To begin with, we examine the electronic structure and plasmon modes of monolayer MoS2 within the aforementioned theory. Our numerical results of the plasmon modes are in very good agreement with those results obtained in welling . It would be worth mentioning that the plasmon mode obtained within DFT-RPA approach for monolayer MoS2 differs with that calculated using the low-energy model Hamiltonian in Scholz and reasons may explain by having the multi-orbital and multiband structures of the system which discussed in asgari2 ; louie , thought the plasmon modes calculated by DFT-RPA approach for monolayer and bilayer graphene are in good agreement with those results obtained using the low-energy model Hamiltonian. Thus it would pretty accurate to calculate the collective modes based on the DFT-RPA specially for complex systems which contain several bands and orbitals. Having investigated these examinations, we can now survey bilayer MoS2 and explore its ground-state electronic properties and charge collective excitations.

III.1 Electronic structures of bilayer MoS2

The band structure of bilayer MoS2 along the high symmetry $\Gamma-M-K-\Gamma$ directions of the BZ is shown in Fig. 1. The blue and red lines correspond to the valence and conduction bands, respectively.

In contrast to monolayer MoS2, bilayer MoS2 is an indirect gap semiconductor that the maximum of the valence band is located at the $\Gamma$ point of the BZ while the conduction band minima is found between the $\Gamma$ and $K$ high symmetry points (Q point) BL . The conduction band states at K point are mainly owing to the localized $d$ orbitals of the Mo atoms, located in the middle of the S-Mo-S layer sandwiches and relatively unaffected by interlayer coupling. However, the states near the $\Gamma$ point are due to the combinations of the antibonding $p_{z}$ orbitals of S atoms and the $d$ orbitals of Mo atoms and have a strong interlayer coupling effect. Therefore, as the layer number changes, the direct excitonic states near the K point are relatively unchanged. The transition at the $\Gamma$ point shifts significantly from indirect gap to a larger and direct one. Therefore, changing in the band structure of layer number is due to quantum confinement and the resulting change in hybridization between the $p_{z}$ orbitals of S atoms and $d$ orbitals of Mo atoms kadan .

III.2 Collective excitation spectra in extrinsic bilayer MoS2

We investigate the dielectric properties of doped bilayer MoS2. The extrinsic systems are simulated by slightly changing the level populations in Eq. (1) and assuming that the doping has a negligible effect on the KS electronic structure.

The static limit of the response function $\chi^{0}(\mathbf{q},\omega=0)$, which contains a number of noteworthy features, is purely real as illustrated in Fig. 2(a) in units of the noninteracting density of states at the Fermi surface for electron doping with $E_{\rm F}=0.06$ eV. The $\lim_{q\to 0}\chi^{0}(\mathbf{q},\omega=0$) is finite and equal to the density of states at the Fermi energy, $N(0)$ a measure of the number of excited states. Apparently, the structure of the static noninteracting response function, in the vicinity of the $q=2k_{\rm F}$, is very similar to that found in three-dimensional electron gas. The static limit of the response function exhibits a Kohn anomaly at $q=2k_{\rm F}$ and it is responsible for several interesting phenomena such as Friedel oscillations. Basically, the static response function $\chi^{0}(\mathbf{q},\omega=0)$ can be Fourier transformed to obtain $\chi^{0}(\mathbf{r},\omega=0)$, which physically represents the response of the density of the noninteracting electron gas to a potential that is located at the origin of the coordinate system vignale . Furthermore, we plot the density-density response function $\chi(\mathbf{q},\omega=0)=\chi^{0}(\mathbf{q},0)/\epsilon(q,0)$ in the same figure. The interacting response function is different from the noninteracting one. The difference with the $\chi^{0}(\mathbf{q},\omega=0)$ is noticeable in the long wavelength limit, where the diverging dielectric function makes $\chi(\mathbf{q},\omega=0)$ vanish as $q$. Also, the real part (black) and the imaginary part (red) of the noninteracting response function $-\chi^{0}(\mathbf{q},\omega)/N(0)$ as a function of $\hbar\omega$ are plotted in Fig. 2(b) for $q=0.6k_{\rm F}$. As we expected the peak of the imaginary part of the response function, $\Im m\chi^{0}(\mathbf{q},\omega)$ is corresponding to a zero in the real part of the response function, $\Re e\chi^{0}(\mathbf{q},\omega)$. Notice that $\Re e\chi^{0}(q,\omega)$ changes sign from negative to positive as $\omega$ moves across the electron-hole continuum. In addition, the density-density response function originates from the inter- and intra-band contributions to this quantity. There are three different peaks in those quantities. One strong peak of the intra-band contribution refers to a main collective mode and two others which are small and may lead two different plasmon modes at small and high energy regions. Basically, the effect of those peaks would have appeared in the electron loss function and we will discuss this point later.

In Fig. 3, we plot the static dielectric function, $\epsilon(\mathbf{q},0)$ of bilayer MoS2 for electron and hole doped cases where the charge density is $6.9\times 10^{13}$ cm*-2*. Interestingly enough, the $\epsilon(\mathbf{q},0)$ for the hole doped case is smaller than the electron doped one and shows a stronger many-body screening in the electron doped case. There is a clear dependence of the static dielectric function with the type of doping, showing a different behavior for electron or hole doping.

Here, we study the collective modes of bilayer MoS2, which are defined by the zeroes of the dielectric function $\epsilon(\mathbf{q},\omega)$. The dispersion relation of the plasmon modes is defined from Re $\epsilon(\mathbf{q},\omega)=0$ which leads to poles in the EEL function $-\Im m\frac{1}{\epsilon(\mathbf{q},\omega)}$, that can be measured by means of the electron-loss spectroscopy.

We begin with Fig. 4 (upper panel), which exhibits the real and imaginary parts of the dielectric function together with their electron energy-loss function of bilayer MoS2 and momentum transfer along $\Gamma K$ direction, namely $q=0.079$Å*-1* for $E_{\rm F}=0.05$ eV ($n_{e}=6.9\times 10^{13}$cm*-2*). Notice that value of the $\Re e\epsilon(q,\omega)$ or $\Im m\epsilon(q,\omega)$ depends strongly on the Fermi energy. Moreover, there is a value of $\omega$ for which both $\Re e\epsilon(q,\omega)$ and $\Im m\epsilon(q,\omega)$ are zero, which it corresponds to the plasmon mode with a long lived mode of the system. Furthermore, as shown in the Fig. 4 numerically (bottom panel), the peak in the electron energy-loss spectra corresponding to a zero in the real part of the dielectric function, at a frequency where the imaginary part, $\Im m\epsilon(q,\omega)$ is almost zero.

In Fig. 5(a), we display the EEL of bilayer MoS2 for the electron doping value $n_{e}=6.9\times 10^{13}$ cm*-2* and for energies $\hbar\omega$ below 1 eV and momenta $\mathbf{q}$ along the high-symmetry $\Gamma K$ direction. It is clear that the EEL is dominated by a narrow peak, which is placed at $\omega\thicksim 0.25-0.35$ eV. Furthermore, we explore other peaks in the system by studying the EEL at small and large energies as shown in Fig. 5(b). There are two peaks in the EEL, which originate from the inter-band contribution of the response function, leads to extra collective modes. Moreover, the $\Im m\epsilon(q,\omega)$ at position of those peaks are finite. Therefore, it turns out that those collective modes are damped and the damping parameter, $\gamma$ would be momentum dependent. In this case, we should verify the condition where $\epsilon(q,\omega_{p}-i\gamma)=0$ is satisfied. This equations leads to two separate equations in which the $\gamma$ is defined by $\gamma(q)=\Im m\chi^{0}(q,\omega)/\partial\Re e\chi^{0}(q,\omega)/\partial\omega$ at $\omega=\omega_{p}$ and the collective mode is obtained by $1-\nu_{q}\Re e\chi^{0}(q,\omega)-\gamma\nu_{q}\partial\Im m\chi^{0}(q,\omega)/\partial\omega=0$ again at $\omega=\omega_{p}$.

We plot three plasmon dispersions of the system obtained from the peaks in the EEL in Fig. 5(c). A damped high-energy mode, one optical mode (in phase) for which the plasmon dispersion exhibits $\sqrt{q}$ in the long wavelength limit originating from low-energy electron scattering and finally a highly damped acoustic mode (out-of-plane) are shown in this figure. The shaded areas show the electron-hole continuum, the region of the $q,\omega$ plane in which $\Im m\chi^{0}(q,\omega)$ differs from zero, that its boundary can be obtained at any specific momentum transfer $\mathbf{q}$ by the difference between system energy at $\mathbf{k_{\rm F}}$ and $\mathbf{k_{\rm F}+q}$. Besides, two damped modes lie on the electron-hole continuum regions. In the optical mode and for longer momentum transfer, the strength of the plasmon is significantly reduced owing to the screening by interband transitions Scholz ; Andersen . Notice that as one goes to larger values of $q$, the optical plasmon energy increases rather slowly; the maximum single-pair energy increases more rapidly. Furthermore, the optical plasmon oscillator strength vanishes as its dispersion approaches the edge of the electron-hole continuum at $q=q_{c}$. In a vicinity of $q\rightarrow q_{c}$, the quantity $|\partial\epsilon(q,\omega)/\partial\omega|$ tends to infinity too.

In Fig. 6 the optical mode and damped acoustic dispersions are shown at a given carrier concentration for (a) electron and (b) hole doping (we neglect to show the other high energy mode). The electron-hole continuum are shown by the solid and dashed lines. It is worth mentioning that the plasmon modes in electron and hole doping graphene, at a given carrier concentration, show the same dynamics. However, it is quite different in bilayer MoS2 Scholz . More importantly, the change of the plasmon modes in the hole doped case is much stronger than that in the electron doped one at given density. As Fig. 1 shows, the structure of the valence bands is quite different compared to the conduction bands. Therefore, there is no electron-hole symmetry in MoS2, and intriguingly the plasmon dispersion of electron and hole doped cases are completely different. The electron and hole effective band masses of monolayer and bilayer MoS2 given in Table 1 that coincide with the effective masses reported in Ref. Yun . The effective band masses of bilayer MoS2 for an electron in the conduction band minima at the Q point and for a hole at the valence band maximum at the $\Gamma$ point are determined to be $m_{e}^{*}=0.56$ and $m_{h}^{*}=1.07$, respectively, in units of the electron bare mass, that thus confirms the anisotropy of the valence and conduction bands. Furthermore, the plasmon modes depend on the Fermi energy values and increases by increasing the Fermi energy.

For the sake of completeness, we compare effective band masses between monolayer and bilayer Mo2 system. From monolayer to bilayer MoS2, the $m_{h}^{*}$ at the $K$ point dramatically increases at the $\Gamma$ point. While $m_{e}^{*}$ at the $K$ point changes slightly at the Q point (Table 1). The change of the conduction minimum band from the $K$ point to the $Q$ point may have no significant effect because of their similar effective masses. In contrast the change of the valence maximum band from the $K$ point in the $\Gamma$ point, induces the dramatic increase of the effective mass.

In order to explore the dependence of the collective modes on the symmetry direction in the BZ, we compare the plasmon dispersions (optical and acoustic modes) of electron-doped bilayer MoS2 for $n_{e}=9.2\times 10^{13}$ cm*-2* along the high-symmetry inequivalent directions $\Gamma M$ and $\Gamma K$ in Fig. 7. In particular, the plasmon dispersions along the $\Gamma M$ is energetically higher than the plasmon dispersions along the $\Gamma K$.

The distribution of the charge densities of the valence band maximum (VBM) and conduction band minimum (CBM) states is shown in Fig. 8 for bilayer MoS2. The CBM states at the $K$ point are distributed within the Mo sublattice while the charge on the VBM at the $\Gamma$ point accumulated on the inner S atoms. This indicates that the states at the $\Gamma$ point have strong interlayer coupling, while the states around the $K$ valley remain largely unaffected. This may also explain why there is quantum effect difference in the plasmon mode for the electron and hole doped systems as shown in Fig. 6.

It is important to mention that when the extrinsic Fermi level crosses two bands, two Fermi surfaces per inequivalent K point, two collective modes in the spectrum emerge. Therefore, two different oscillations of charge carriers occur, leading to the appearance of two collective modes: one in which the oscillation is in phase (the optical mode) and another one is the out of phase oscillation (acoustic mode) Eguiluz . In our numerical study, different charge-carrier concentrations are considered such a way that the extrinsic Fermi levels cross only one band, for instance, there is only one intersection around $Q$ point in the electron doped case.

IV CONCLUSION

We have used a full DFT simulations together with RPA analysis to calculate the energy loss functions and plasmon dispersions of extrinsic bilayer MoS2 for energies up to 1 eV. The interacting density-density response function as well as dielectric function are calculated using the theory in $q,\omega$ plane. In extrinsic bilayer MoS2, the electron energy-loss function is calculated and we have explored that the function is sensitive to doping is either positive or negative, as there is no electron-hole symmetry in MoS2. Furthermore, we have calculated the collective modes of the system for various charge carrier densities. We have found three different collective modes in the system. One mode occur at high-energy and for which the mode is damped. Other mode refers to the optical mode in which the plasmon dispersion exhibits a $\sqrt{q}$ dispersion as the conventional plasmon of a two-dimensional electron gas originating from low-momentum carrier scattering. An acoustic mode in the system is observed however this mode is heavily damped. In addition, we have observed that the plasmon modes of electron and hole doping of bilayer MoS2 are not equivalent and the discrepancy is owing to the fact that the Kohn-Sham band dispersions are not symmetric for energies above or below the zero Fermi level. Finally, we should emphasize that the material-specific dielectric function considering the multi-orbital and multiband structures are needed to obtain realistic plasmon dispersions in MoS2. Our numerical finding can be investigated by current experiments.

V acknowledgments

Z. T. would like to thank Iran Nanotechnology Initiative Council for their support. This work is also partially supported by Iran Science Elites Federation.