## Theoretical Framework
### Computational Method  
There are several DFT packages that are available in the world, and each one mainly uses different basis sets. PWscf, a core component of the Quantum ESPRESSO distribution [24], performs many different kinds of self-consistent calculations of electronic structure properties such as ground-state energy and one-electron Kohn-Sham orbitals; within density functional theory, using a plane-wave basis set and pseudopotentials. The expression of the Kohn-Sham orbitals in the plane waves basis has the form:  
\[\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}}\]  
where \(\Omega\) represents the crystal volume, \(G\) is the reciprocal lattice vectors and \(k\) is the quasi-wavevectors of the first Brillouin zone (BZ). The coefficients \(C_{n\mathbf{k}}(\mathbf{G})\) are obtained by solving the LDA-KS equations self-consistently. Since the core electrons of an atom are highly localized, it would be difficult to implement them using the plane waves basis sets. Actually, a very large number of plane waves is required to expand their wave functions. Furthermore, the contribution of the core electrons to bonding compared to those of the valence electrons is usually negligible. Therefore, it is practically desirable to replace the atomic potential owing to the core electrons with a pseudopotential that has the same effect on the valence electrons [25]. There are mainly two kinds of pseudopotentials, norm-conserving soft pseudopotentials [26] and Vanderbilt ultra-soft pseudopotentials [27]. The first part of our calculations includes determining the KS ground-state of pristine 2D crystalline materials and the corresponding wave functions and energies. We carry out the first-principles simulations based on the DFT simulations implemented in the QUANTUM ESPRESSO code. The calculation of the density-density response functions \((\chi(q,\omega))\) is performed employing our own code. The computation of this quantity will be discussed in the next section. In this study, the electronic structures of 2D materials are computed using the Perdew-Zunger local-density approximation [28], unless otherwise stated. Furthermore, we use this throughout the norm-conserving pseudopotentials and the plane wave basis. The energy convergence criteria for electronic and ionic iterations are set to be \(10^{-5}\) eV and \(10^{-4}\) eV, respectively. Careful testing of the effect of the cutoff energy on the total energy can be implemented to determine a suitable cutoff energy. The cutoff energy is required to obtain a particular convergence precision. Well-converged results are found with a kinetic energy cutoff equal to 50 Ry. The \(k\)-point grid is another calculated parameter that must be considered, and it is used to approximate integrals of some property, e.g., the electron density over the entire unit cell. Notice that the integration is performed in reciprocal space (in the first Brillouin zone) for convenience and efficiency. We thus use a Monkhorst-Pack [29]\(k\)-point grid, which is essentially a uniformly-spaced grid in the Brillouin zone. Geometry optimization and ground-state calculations are carried out on the irreducible part of the first BZ, using a \(\Gamma\)-centered and unshifted Monkhorst-Pack grid of \(60\times 60\times 1\)\(k\)-points for graphene and MoS\({}_{2}\) and \(30\times 40\times 1\)\(k\)-points for phosphorene. The converged electron density is then used to calculate the KS electronic structure, i.e., the single-particle energies and orbitals on a denser \(k\)-point mesh. We perform calculations to check the convergence of the plasmon spectra with respect to the \(k\)-point sampling to obtain reliable spectra, and these are listed in Table 1.  
The equilibrium distance between two layers in bilayer materials is also determined by varying the interlayer distance, while keeping the in-plane lattice constant fixed at the monolayer value.