## Theoretical Framework
### Density-Density Response Function  
A central quantity in the theoretical formulation of the many-body effects in electronic systems is the non-interacting dynamical polarizability function [30]\(\chi^{(0)}(\mathbf{q},\omega)\) for a finite chemical potential, \(\mu\). Here, we would like to emphasize that we have calculated the non-interacting polarization function for doped graphene, MoS\({}_{2}\) and phosphorene, both monolayer and bilayer structures based on DFT simulations. The expression of the non-interacting density-density response function of a three-dimensional periodic electrons in the reciprocal space is:  
\[\chi^{0}_{\mathbf{G}\mathbf{G}^{\prime}}(\mathbf{q},\omega)=\frac{2}{\Omega} \sum_{\mathbf{k},\
u,\
u,\sigma}\frac{f_{\
u}(\mathbf{k})-f_{\
u^{\prime}}( \mathbf{k}+\mathbf{q})}{\hbar\omega+\varepsilon_{n}(\mathbf{k})-\varepsilon_ {m}(\mathbf{k}+\mathbf{q})+i\eta}\rho^{\mathbf{k}\mathbf{q}}_{\
u\
u^{\prime}} (\mathbf{G})\rho^{\mathbf{k}\mathbf{q}}_{\
u\
u^{\prime}}(\mathbf{G}^{\prime}) ^{*} \tag{5}\]  
which is obtained from the Adler-Wiser periodic system [31]. Here, \(\varepsilon_{n}(\mathbf{k})\) and \(\varepsilon_{m}(\mathbf{k}+\mathbf{q})\) denote the empty and filled bands and \(f_{n}(\mathbf{k})=\theta(E_{\mathrm{F}}-\varepsilon_{n}(\mathbf{k}))\) is the Fermi-Dirac distribution of the charge carrier with energy \(\varepsilon_{n}(k)\) at temperature zero. Furthermore, \(G\) and \(G^{\prime}\) are the three-dimensional (3D) reciprocal lattice vectors, and \(\omega\) is the frequency. In this theory, the linear combination of plane-waves is used to determine the KS single-particle orbitals of the DFT. The KS wave functions are normalized to unity in the crystal volume \(\Omega\). The sum is over a special set of \(k\) vectors and energy bands (\(\
u\) and \(\
u^{\prime}\)). In Table 1, the sampling \(k\) point in the BZ is indicated for different 2D materials in order to fully converge the results. The factor of two accounts for the spin degeneracy, and \(\eta\) is a small (positive) lifetime broadening parameter. The matrix elements of Equation (5) have the form:  
\[\rho^{\mathbf{k}\mathbf{q}}_{\
u\
u^{\prime}}(\mathbf{G})=<\Phi_{\
u\mathbf{k }}|\epsilon^{-i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}}|\Phi_{\
u^{\prime} \mathbf{k}+\mathbf{q}}>\Omega \tag{6}\]  
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 and when expanded in the plane-wave basis have 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}} \tag{7}\]  
where the coefficients \(C_{n\mathbf{k}}(\mathbf{G})\) are obtained by solving the LDA-KS equations self-consistently. The exact interacting density-response function can be obtained in the framework of the DFT, as follows [32]:  
\[\chi_{\mathbf{G}\mathbf{G}^{\prime}}=\chi^{0}_{\mathbf{G}\mathbf{G}^{\prime} }+\sum_{\mathbf{G_{1}G_{2}}}\chi^{0}_{\mathbf{G_{1}}}{}^{\upsilon}\mathbf{G_{ 1}}{}^{\upsilon}\mathbf{G_{2}}{}^{\chi}\mathbf{G_{2}}{}^{\prime} \tag{8}\]  
where \(\upsilon_{\mathbf{G}\mathbf{G}^{\prime}}\) represents the Fourier coefficients of an effective electron-electron interaction. In the electron liquid, the bare Coulomb interaction is given by \(\upsilon^{0}_{\mathbf{G}\mathbf{G}^{\prime}}=4\pi e^{2}\delta_{\mathbf{G} \mathbf{G}^{\prime}}/\epsilon|\mathbf{q}+\mathbf{G}|^{2}\) where \(\epsilon\) is the average dielectric constant of the environment. In all our numerical results, we consider \(\epsilon=1\). The RPA procedure, an approximation valid in the high-density limit, takes into account electron