## 3 Result and Discussion
### Bilayer MoS\({}_{2}\)  
The AA' stacking, in which the layers are exactly aligned, with Mo over S is the stablest stacking for the MoS\({}_{2}\) bilayer, and it is shown in Figure 11a,b. From the band structure plotted in Figure 11c, based on our DFT simulations, along the high symmetry \(\Gamma-M-K-\Gamma\) directions, it can be seen that bilayer MoS\({}_{2}\) has an indirect band gap in contrast to the direct band gap of the corresponding monolayer. In fact, the conduction band minima are located between the \(\Gamma\) and \(K\) high symmetry points (\(Q\) point), while the valence band maxima are located at the \(\Gamma\) point of the BZ, revealing the indirect band gap, and this is in good agreement with previous calculations [77; 78; 79]. More analysis of the band structure of MoS\({}_{2}\) can show why this material experiences an indirect to direct band gap transition when its bulk or multilayers are replaced by a monolayer. In fact, the states originating from mixing of Mo(\(d_{z^{2}}\)) orbitals and the S(\(p_{z}\)) orbitals at \(\Gamma\) are fairly delocalized and have an antibonding nature. With increasing the separation between consecutive MoS\({}_{2}\) layers, the layer-layer interaction decreases and lowers the energy of the antibonding states, and consequently, the valence band maximum shifts downward. The states at \(K\), which are of the \(d_{xy}-d_{x^{2}-y^{2}}\) character, are mostly unaffected by interlayer spacing. Thus, in the limit of widely-separated planes, i.e., monolayer MoS\({}_{2}\), the material becomes a direct gap semiconductor [80]. The non-interacting response function of bilayer MoS\({}_{2}\) for \(q=0.084\) A\({}^{-1}\) with the charge carrier concentration of \(n=6.9\times 10^{13}\) cm\({}^{-2}\) shown in Figure 12. We can see that in bilayer MoS\({}_{2}\), both interband and intraband contributions of non-interacting response function are considerable. The intraband contribution of this quantity is quite strong and refers to the main collective mode, although its interband term is weak, and this can lead to new modes in the system; thus, in the following, we will explain it further. In Figure 13, we display the EEL function and plasmon spectrum of bilayer MoS\({}_{2}\) for the electron doping with \(n=6.9\times 10^{13}\) cm\({}^{-2}\) and for energies below 1 eV and momenta \(\mathbf{q}\) along the high-symmetry \(\Gamma-K\) direction. We observe that plasmonic features in electron-doped bilayer MoS\({}_{2}\) are mainly characterized by a square root mode in small \(q\). In the paper by [PERSON], plasmons in bilayer NbS\({}_{2}\) were studied [70], and they obtained very similar plasmon dispersions. Figure 11: **(a,b)** Top and side view of bilayer molybdenum disulfide (BL_MoS\({}_{2}\)). **(c)** The band structure of BL_MoS\({}_{2}\) along the high symmetry \(\Gamma-K-M-\Gamma\) directions. The Fermi level is set at 0 eV, and interlayer separation \(d\) is 5.89 Å.  
Figure 12: The non-interacting response function in units of the Fermi-level density of states as a function of \(\hbar\omega\) of bilayer MoS\({}_{2}\) for \(q=0.084\) Å\({}^{-1}\). As previously mentioned, in the two-layer systems, we can find two plasmon modes, which can be regarded as symmetric and antisymmetric combinations of two unperturbed monolayer plasmons. One of them has a nearly linear dispersion, referred to as the acoustic plasmon mode. In this case, we find that the low-energy dynamical dielectric function in bilayer MoS\({}_{2}\) is controlled by both interband and intraband contributions, leading to an extra collective mode. Also, there are a damped high energy mode and a highly damped acoustic mode, which originates from the interband transitions. More details of these calculations can be found in [81]. In Figure 14, we display the damping parameter of plasmon modes in bilayer MoS\({}_{2}\) for the electron doping value \(n=6.9\times 10^{13}\) cm\({}^{-2}\) for both the acoustic and higher plasmon modes.