## Theoretical Framework
### Density-Density Response Function
\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline
**2D Structures** & **Lattice** & **a (Å)** & **b (Å)** & \multicolumn{3}{c}{**Gap (eV)**} & \multicolumn{3}{c}{**MP**} \\ \cline{4-9}  & & & & **Monolayer** & **Bilayer** & **Monolayer** & **Bilayer** \\ \hline graphene & Hexagonal & 2.46 & \(-\) & 0.00 & 0.00 & \(101\times 101\times 1\) & \(201\times 201\times 1\) \\ MoS\({}_{2}\) & Hexagonal & 3.15 & \(-\) & 1.80 & 1.05 & \(101\times 101\times 1\) & \(151\times 151\times 1\) \\ phosphorene & Rectangular & 4.62 & 3.30 & 0.98 & 0.63 & \(60\times\) & \(80\times 1\) & \(120\times\) & \(160\times 1\) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: The lattice constants (\(a,b\)) and band gap (in eV) of different 2D crystalline materials. Sampling of the reciprocal space Brillouin zone (BZ) is done by a Monkhorst-Pack (MP) grid for the considered 2D materials. 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 RPA follows from a microscopic approach whose main assumption is that the electrons respond not to the bare Coulomb potential, but to an effective potential resulting from the dynamical rearrangement of charges in response to the Coulomb potential. The long-range behavior of the Coulomb interaction allows non-negligible interactions between repeated planar arrays even at a large distance. This unphysical phenomenon can be removed by replacing \(v_{\mathbf{G}\mathbf{G}^{\prime}}\) with the truncated Fourier integral over the cutoff plane axis \((z)\)[33, 34], and thus, we have:  
\[v_{\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}dz^{\prime}e^{i(G_{z}z-G_{z}^{\prime}z^{\prime})-|\mathbf{q}+\mathbf{g}||z+z^{ \prime}|} \tag{9}\]  
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 zero owing to a uniform background of positive charge. In the framework of the 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:  
\[\varepsilon_{\mathbf{G}\mathbf{G}^{\prime}}=\delta_{\mathbf{G}\mathbf{G}^{ \prime}}-\sum_{\mathbf{G}_{1}}v_{\mathbf{G}\mathbf{G}_{1}}^{0}\chi_{\mathbf{G }_{1}\mathbf{G}^{\prime}}^{0} \tag{10}\]  
and the plasmon modes are established by the zero in the real part of the macroscopic dielectric function given by:  
\[\varepsilon(q,\omega)=\frac{1}{(\varepsilon^{-1})_{\mathbf{G}\mathbf{G}^{ \prime}}}|_{\mathbf{G}=\mathbf{G}^{\prime}=0} \tag{11}\]  
as long as there is no damping. The electron-energy loss (EEL) is proportional to the imaginary part of the inverse dielectric function, which is given by:  
\[E_{EEL}(q,\omega)=-\Im m[1/\varepsilon(q,\omega)] \tag{12}\]  
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}}\)[35] function. In addition, one may use the non-local dynamical conductivity to describe the electronic processes and light-matter interactions. This can be simply achieved by writing the longitudinal conductivity in terms of the density-density response function through the relation:  
\[\sigma(q,\omega)=\frac{ie^{2}\omega}{q^{2}}\chi(q,\omega) \tag{13}\]  
We note that Equations (11)-(13) are the fundamental physics describing the optoelectronic interactions in 2D crystalline materials.