Introduction {#Sec1}
============

Plasmons are quantum collective motions of electrons in solids arising from the long-range Coulomb interaction. Since these propagating collective modes are strongly affected by boundary/geometric effects, plasmon modes can be tuned by assembling nanostructured materials with repeated patterns, or by selectively exciting modes that only exists on the surface of metals, such as the surface plasmon polaritons (SPPs)^[@CR1]--[@CR4]^. In addition, there has been great interest in plasmons in atomically thin (one to a few atomic layer) crystalline metals such as graphene^[@CR5]--[@CR9]^. Graphene plasmons are highly tunable and display the typical dispersion relation^[@CR10]^ $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _p \sim \sqrt q$$\end{document}$ of an ideal 2D electron gas for *q* → 0, and can be observed with techniques ranging from electron-energy loss spectroscopy (EELS) to nano-imaging of interference patterns generated with atomic-force microscopy (AFM) tips^[@CR1],[@CR2],[@CR4],[@CR11]^.

Besides graphene, plasmons in other atomically thin quasi-2D metals have been reported^[@CR12]--[@CR14]^, some displaying intriguing dispersion relations. Previous ab initio studies have shown that plasmons on monolayer metallic transition metal dichalcogenides (TMDs) disperse as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _p \sim \sqrt q$$\end{document}$ at small *q* and can have high energy (\~1 eV) that are nearly dispersionless for a large range of wavevectors (*q* \~ 0.1 Å^−1^ to 0.3 Å^−1^)^[@CR15]--[@CR18]^. Similar findings were found for other quasi-2D metals, such as borophene^[@CR19]^. The previous studies however have basically overlooked the origin and significance of this dispersionless behavior. As will become clear in the present paper, dispersionless quasi-2D plasmons could find a variety of exciting applications, as flat dispersion relations can translate into real-space localization of plasmon wave packets using practical excitations. Moreover, it is important to unravel what is the physical origin of these dispersionless plasmons and to assess whether they are (or can be controllably made to be) sufficiently long lived for technological applications^[@CR20]^.

Here, we demonstrate rigorously why the plasmon dispersion relation in real, quasi-2D metals is qualitatively different from that in idealized 2D metals (i.e., a homogeneous electron gas), in metallic slabs, and of the SPPs, tending to flatten at relatively large wavevectors for *any* quasi-2D crystals regardless of their chemical composition or nature of the energy bands crossing the Fermi energy. The local-fields screening arising from interband transitions in real quasi-2D crystalline metals^[@CR21]--[@CR28]^, which is not present in the 2D electron gas model, is the physics behind these dispersionless intrinsic plasmons. This behavior is universal for any atomically thin quasi-2D crystalline metal as long as there is a well-separated band (or set of bands) either completely occupied or unoccupied near the Fermi energy (as shown in Fig. [1a](#Fig1){ref-type="fig"}), a condition that is not fulfilled in doped graphene.Fig. 1Plasmon dispersion relation of an atomically thin quasi-2D crystalline metal.**a** Calculated GW quasiparticle band structure of a prototypical quasi-2D metal, monolayer TaS~2~. The arrows highlight the intraband and interband electronic transitions that contribute to a plasmon excitation. The Fermi energy is set at zero. **b** Calculations of the plasmon dispersion relation of monolayer TaS~2~. The green curve is the result of the full ab initio calculation, which displays a flattening of the plasmon dispersion relation; the blue curve corresponds to the result of a calculation including only the intraband transitions, which displays a dispersion relation much closer to \~$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sqrt q$$\end{document}$ (red curve). The inset shows the small-*q* behavior of the plasmon dispersion relation, where both calculations follow a $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sqrt q$$\end{document}$ dependence. **c** Effective 2D dielectric function *ε*~eff~(**q**) for two point charges located at the innermost portion of three different quasi-2D materials. Since *ε*~eff~(**q**) does not include intraband (but only interband) transitions, its behavior is similar to quasi-2D semiconductors. For very small wavevectors, there is no screening (*ε*~eff~(**q** = **0**) = 1)) since very long wavelength perturbations correspond to charges being so far apart that they only experience the screening of vacuum. For increasing wavevector **q**, the charges get closer, so more of the field lines connecting them are confined to the interior of the quasi-2D material. *ε*~eff~(**q**) increases linearly with *q* until it reaches near a maximum, and subsequently decreases with increasing **q** after the maximum due to the intrinsic inability of the medium to screen very short-wavelength excitations.

Results {#Sec2}
=======

Universal dispersionless plasmons in real quasi-2D metals {#Sec3}
---------------------------------------------------------

To shed light on this universal flat plasmon dispersion relation, we depict in Fig. [1](#Fig1){ref-type="fig"} the computed plasmon dispersion relation from first principles for monolayer TaS~2~. To perform the calculation, with atomistic local-field effects, we have to develop an algorithm to compute plasmons for arbitrarily small wavevectors, which heretofore was a bottleneck for ab initio calculations (see Methods section). As expected from a homogeneous 2D electron gas picture, we find that plasmons in monolayer TaS~2~ in Fig. [1b](#Fig1){ref-type="fig"} follows a $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sqrt q$$\end{document}$ dispersion relation up to a wavevector of *q* \~0.05 Å^−1^; however, the dispersion relation deviates significantly from this $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sqrt q$$\end{document}$ relation afterwards (it becomes flattened), with an onset that is large compared to the wavevector of visible light but still very small compared to the Brillouin zone size of \~1.1 Å^−1^. This deviation is related to the broken continuous translational symmetry in real crystals (which only have discrete crystal translational symmetry). A homogeneous 2D electron gas with continuous in-plane translation symmetry can only have a single partially occupied electronic band that extends to infinite *q*, and only the intraband electronic transitions contribute to the system's dielectric screening. On the other hand, both intraband and interband electronic transitions contribute to screening in a real quasi-2D metal such as monolayer TaS~2~, and, as we argue here, the interband transitions are the reason behind the appearance of virtually dispersionless plasmons. Indeed, the plasmon dispersion relation of monolayer TaS~2~ obtained by only considering intraband transitions (blue curve in Fig. [1b](#Fig1){ref-type="fig"}) follows a $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sqrt q$$\end{document}$ dispersion relation up to much larger wavevectors^[@CR16]^.

To gain some fundamental insights, we derive a general microscopic analytic expression for the plasmon dispersion relation for a quasi-2D metal, going beyond the 2D homogenous electron gas model. Our approach is to recast rigorously the plasmon excitations of a real quasi-2D material into contributions from two components: (1) an effective 2D electron gas (which is defined to consist of only those states of the system in the bands crossing the Fermi level); and (2) a nonlocal polarizable medium that incorporates the effects of the interband transitions in the real material (see Fig. [1a](#Fig1){ref-type="fig"}). Under the condition of long plasmon wavelengths, we derive (see SI section), by solving for the zeros of the total dielectric function, the following expression for the plasmon frequency,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _p^2 = \frac{{2\pi e^2}}{{\varepsilon _{{\mathrm{eff}}}\left( {\mathbf{q}} \right)}}\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}q,$$\end{document}$$where the symbolic notation $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}$$\end{document}$ is a function that is a **q**-dependent generalization of the ratio of the carrier density *n* by the effective mass $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${m^ \ast }$$\end{document}$ of an anisotropic 2D system (see Methods section), which only depends on the nature of the band (or set of bands) crossing the Fermi energy. The plasmon frequency is further dependent on an effective 2D dielectric constant, *ε*~eff~(**q**), which includes contributions to the screening due to *all* the intrinsic interband transitions (indicated by the green arrows in Fig. [1a](#Fig1){ref-type="fig"}) of the quasi-2D material, excluding the intraband contributions from the set of bands crossing the Fermi energy, as well as due to any external environment (such as a substrate). As demonstrated below, *ε*~eff~(**q**) is responsible for the flat dispersion relation discussed above.

The fundamental idea here is that the effective dielectric function *ε*~eff~(**q**) from the interband transitions of the material essentially has the same analytical expression and behavior as the dielectric function derived in the context of quasi-2D semiconductor screening (see, e.g., Ref. ^[@CR28]^), and it allows us to understand and determine, from Eq. ([1](#Equ1){ref-type=""}), plasmon dispersion relations in quasi-2D metals using established models for quasi-2D semiconductors^[@CR21]^ together with information from ab initio calculations^[@CR23],[@CR26]--[@CR30]^. In particular, the semiconductor-like interband effective 2D dielectric screening *ε*~eff~(**q**) of a monolayer material is not a constant, but increases linearly from unity with wavevector at small *q*, since very long-wavelength excitations effectively only experience screening from the vacuum or substrate (see Fig. [1c](#Fig1){ref-type="fig"}).

We show that, for long wavelength excitations, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}$$\end{document}$ does not depend on *q*, while $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{{\mathrm{eff}}}\left( {\mathbf{q}} \right)$$\end{document}$ can be approximated as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{{\mathrm{eff}}}\left( q \right) \approx \frac{{1 + \varepsilon _s}}{2} + \rho _0q$$\end{document}$, where *ε*~*s*~ is the dielectric constant of the insulating substrate and *ρ*~0~ is a screening length intrinsic to the quasi-2D material under consideration arising from the interband transitions (see Methods section). In this limit,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _{\boldsymbol{p}} \approx \sqrt {2\pi {\mathrm{e}}^2\left[ {\frac{{\it{n}}}{{m^ \ast }}} \right]_{{\mathrm{eff}}}} \sqrt {\frac{{\it{q}}}{{\frac{{1 + \varepsilon _{\boldsymbol{s}}}}{2} + \rho _0{\boldsymbol{q}}}}} .$$\end{document}$$Equation ([2](#Equ2){ref-type=""}) is one of the key results of this work. It shows that plasmons in real quasi-2D metals display a universal tendency to flatten for wavevectors above a critical wavevector *q*~*c*~ of the order of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\sim \frac{{\varepsilon _s + 1}}{{2\rho _0}}$$\end{document}$, resulting in a plasmon frequency of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _p^{{\mathrm{flat}}} = \sqrt {\frac{{2\pi {\mathrm{e}}^2}}{{\rho _0}}\left[ {\frac{{n}}{{{m}^ \ast }}} \right]_{{\mathrm{eff}}}}$$\end{document}$ that is independent of the substrate dielectric constant *ε*~s~, although importantly the substrate dielectric constant can dramatically change the value of *q*~*c*~ and the dispersion relation for *q* \< *q*~*c*~. We emphasize here that this phenomenon arises without considering interactions with photons, and hence is not a plasmon-polariton effect. Moreover, Eq. ([2](#Equ2){ref-type=""}) holds in practice even for larger plasmon wavevectors, since the effective ratio $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}} \right]_{{\mathrm{eff}}}\left( q \right)$$\end{document}$ decreases for larger values of *q* due to a typical reduction of the quasiparticle band velocity away from the Fermi surface, while the effective interband screening *ε*~eff~(*q*) also saturates for larger *q*^[@CR28]^. Since these two effects partially cancel each other, the main features predicted by Eq. ([2](#Equ2){ref-type=""})−that the plasmon dispersion relation flattens at larger *q* and has an asymptotic energy that does not depend on the substrate--holds from our ab initio calculations for all quasi-2D materials studied here up to *q* \~0.3 Å^−1^.

The plasmon dispersion relation in atomically thin quasi-2D metals is in fact conceptually different from that of metallic slabs with thickness large compared to interatomic distance, which support on each surface (within the classical dielectric response framework) a surface plasmon-polariton dispersion relation that flattens at a wavevector $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$q_c \sim \frac{{\omega _p^{{\mathrm{bulk}}}}}{{\sqrt 2 c}}$$\end{document}$, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _p^{{\mathrm{bulk}}}$$\end{document}$ is the bulk plasmon energy^[@CR31],[@CR32]^ (Fig. [2a](#Fig2){ref-type="fig"}). In contrast, as we demonstrated above, systems with atomically thin thickness support dispersionless quasi-2D plasmons with a very different *q*~*c*~ \~ 1/*ρ*~0~: the 2D screening length *ρ*~0~ in these systems (*ρ*~0~\~25 Å for monolayer metallic TMDs) is not directly related to the thickness of the material nor the bulk plasmon energy^[@CR28]^; it is a new length scale that emerges due to the quantum nature of these systems, which depends primarily on the energy for the interband transitions and the in-plane polarizability of the orbitals involved in these transitions. In addition, as mentioned above, neither retardation effects in the Coulomb interaction nor the hybridization of electronic transitions with the transverse electromagnetic fields are included in our calculations, so conventional hybrid plasmon-polariton modes cannot be responsible for the dispersionless plasmons in atomically thin metals. Further, the dispersionless quasi-2D plasmons here are also qualitatively and quantitatively different from plasmons in conventional thin metallic slabs, which support surface plasmons with an energy relation $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{{\hbar \omega _p^{{\mathrm{bulk}}}}}{{\sqrt 2 }}\sqrt {1 \pm e^{ - qd}}$$\end{document}$, where the positive (negative) sign corresponds to a mode of antisymmetric (symmetric) relative oscillations of the electron density on the opposite surfaces, and *d* is the thickness of the slab^[@CR31],[@CR33],[@CR34]^ (Fig. [2b](#Fig2){ref-type="fig"}). For materials such as bulk metallic TMDs, the bulk plasmon energy is of the order of 20 eV; a classical (and inappropriate) dielectric description of a metallic film of atomically thin thickness (*d* \~ 5Å) would yield a surface plasmon energy and a critical wavevector \~1/*d* for dispersionless surface plasmons that are both over one order of magnitude larger than the values found in our ab initio calculations for monolayer TaS~2~ (Fig. [2c](#Fig2){ref-type="fig"}). These differences highlight that the dispersionless plasmons in quasi-2D materials given by Eq. ([1](#Equ1){ref-type=""}) are conceptually different from classical surface plasmons or hybrid plasmon-polariton modes, and they have instead their origin in the broken translational symmetry in atomically thin metals.Fig. 2Comparison of plasmon excitations on the surface of a bulk metal, on a metallic thin film, and in atomically thin metals, which are different in nature.**a** Dispersion relation for a surface plasmon and the lower-energy surface plasmon-polariton (SPP) branch on the surface of a three-dimensional metal with a bulk plasmon energy $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _p^{{\mathrm{bulk}}}$$\end{document}$ \~20 eV, similar to that of bulk TaS~2~. Retardation effects of the Coulomb interaction, which give rise to hybrid SPPs modes, are noticeable around the intersect of the light cone with the surface plasmon dispersion (inset). The colored intensity sketch represents the electric potential associated with the surface plasmon. **b** Same as **a**, but for a thin metallic film with thickness *d* = 100 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\AA$$\end{document}$. The two distinct surface plasmon modes (red and blue solid lines), as well as the two lower-energy SPP modes (red and blue dashed lines), correspond to antisymmetric and symmetric relative oscillations of the electron density on the opposite surfaces, with the symmetric mode being lower in energy. **c** Dispersion relation for monolayer TaS~2~. The energy of the intrinsic plasmons in the dispersionless region is one order of magnitude smaller than that of the surface plasmons shown in **b**, and the critical wavevector 1/*ρ*~0 ~≈ 0.04 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\AA}^{ - 1}$$\end{document}$ that defines the flat dispersion in monolayer TaS~2~ is about an order of magnitude smaller than that of a classical, atomically thin metallic film of the same thickness, which is 1/*d*~1L ~= 0.2 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\AA}^{ - 1}$$\end{document}$ (*ρ*~0 ~= 25 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\AA$$\end{document}$ is the screening length and *d*~1L~ \~ 5 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\AA$$\end{document}$ the thickness of monolayer TaS~2~). Moreover, when interaction with the light field is considered, the intersection of the light cone with the intrinsic plasmon dispersion in monolayer TaS~2~ occurs at a wavevector that is orders of magnitude smaller than the critical wavevector shown in Fig. [1b](#Fig1){ref-type="fig"}. Thus, the dispersionless plasmons in real quasi-2D materials that is the focus of this work are not the result of hybridizing a plasmon with a photon, i.e., it is not the same as the traditional plasmon polariton (inset).

We illustrate the accuracy of our expression in Eq. ([1](#Equ1){ref-type=""}) by computing, from first principles^[@CR35],[@CR36]^, the plasmon dispersion relation of monolayer TaS~2~ on a number of substrates spanning a range of *ε*~s~ \~ 2 to 10 (Fig. [3b](#Fig3){ref-type="fig"}). Our novel algorithm for fine sampling of transitions near E~F~ enables the ab initio calculation of the polarizability at arbitrarily small wave vectors **q** and without requiring mappings to model Hamiltonians, solving an important bottleneck in previous ab initio calculations (see Methods section). We compare our direct ab initio results to those from Eq. ([1](#Equ1){ref-type=""})--that is, by computing $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}({\mathbf{q}})} \right]_{{\mathrm{eff}}}$$\end{document}$ and *ε*~eff~(**q**) for a suspended monolayer TaS~2~, while generalizing *ε*~eff~(**q**) for the system on different substrates using the Keldysh model but based on ab initio parameters without any additional fitting (see Methods section). The excellent agreement in Fig. [3b](#Fig3){ref-type="fig"} confirms the idea that plasmons in quasi-2D metals can be understood as an effective 2D electron gas that experiences an additional screening (due to interband transitions and substrates) of a form of that is characteristic of quasi-2D semiconductors. Figure [3b](#Fig3){ref-type="fig"} also highlights that the plasmon group velocity can be changed by roughly an order of magnitude with a suitable choice of substrate, in addition to the tunable plasmon energies.Fig. 3Ab initio results and model for plasmon dispersion relation of monolayer TaS~2~.**a** A real quasi-2D metal can be formulated as: a 2D electron gas with renormalized properties given by the energy band (or bands) that crosses the Fermi energy (thick band in Fig. [1a](#Fig1){ref-type="fig"}) embedded in an environment with an effective screening, which accounts for the interband transitions and environmental effects (Eq. ([1](#Equ1){ref-type=""}) and Fig. [1c](#Fig1){ref-type="fig"}). **b** Calculated plasmon dispersion relation of monolayer TaS~2~ on a variety of substrates (Polytetrafluoroethylene (PTFE), BN, WSe~2~)--calculated from first principles (solid colored lines) and computed with the expression given in Eq. ([1](#Equ1){ref-type=""}) (dash colored lines). For the latter, we fit the ab initio computed intrinsic effective dielectric screening *ε*~eff~(**q**) of suspended TaS~2~ to the form of a Keldysh model^[@CR21]^. The resulting curves of the model on different substrates are then obtained without any additional fitting. The Landau damping regions are also shown (gray shaded regions), which is associated with the electron-hole continuum of either intraband or interband transitions (Fig. [1a](#Fig1){ref-type="fig"}).

Lifetime of quasi-2D plasmons {#Sec4}
-----------------------------

We now address the fundamental, and technologically important, question of what the lifetime of these collective excitations is. Plasmons are strongly damped when they are inside the Landau damping region, wherein they can decay into free electron-hole pairs^[@CR37]^. Previous calculations^[@CR16]^ on monolayer metallic TMDs employed density-functional theory (DFT) to compute the Landau damping regions and showed that plasmons in such systems coexist with single-particle excitations for *q* ≳ 0.15 Å^−1^ (that is, they are inside Landau damping regions for such large *q* plasmons), leading to a conclusion of short lifetimes. However, as static DFT is a ground-state theory, it does not yield the needed accurate quasiparticle energies to determine the correct, relevant Landau damping regions^[@CR35],[@CR36]^ and hence plasmon lifetimes. By means of first-principles GW calculations for the quasiparticle excitations, we can accurately evaluate the Landau damping regimes in quasi-2D materials. Here, we focus on the monolayer metallic TMDs, as these materials are predicted (see Fig. [3](#Fig3){ref-type="fig"}) to have large plasmon energies^[@CR16],[@CR24]^ and are intrinsically metallic and stable above their charge-density-wave^[@CR38]^ and/or ferromagnetic^[@CR39]^ transition temperature, if any. We find that, in these systems, DFT underestimates the separation between the half-filled *d* shell bands and the rest of the occupied bands by almost 200 meV and does not yield the correct threshold for single-particle interband excitations. Our quasiparticle GW calculations, on the contrary, show that the plasmon dispersion relation is well outside of the Landau damping regions, and coupling to finite-*q* interband transitions (electron-hole pair creations) are *not* a limiting factor for the plasmon lifetime for the materials studied here. The plasmons only reach the continuum of intraband transitions (thus having a decay channel to electron-hole pairs) after *q* \~ 0.4 Å^−1^ (Fig. [3b](#Fig3){ref-type="fig"}). Hence, monolayer metallic TMDs represent potential systems to implement the concept of lossless metals to support long-lived, dispersionless plasmons^[@CR40]^. Other layered metals displaying a favorable electronic structure with interband transitions sufficiently separated in energy may also be good candidates^[@CR14]^.

Plasmons not in the Landau damping regions may also decay to free electron-hole pairs, through a higher order process, by additionally emitting or absorbing a phonon. While a first-principles calculation of this higher-order decay channel is quite involved, it can be estimated from a density of final states for decay, *D*(**q**), and an effective plasmon-phonon coupling matrix element $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g$$\end{document}$ (see Methods section). We compute *D*(**q**) and estimate $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g$$\end{document}$ for graphene and monolayer TaS~2~, both quantities from first principles. Even though *D*(**q**) is about two orders of magnitude larger for monolayer TaS~2~ than for graphene due to the larger band velocity in graphene, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {\bar g} \right|^2$$\end{document}$ is about 20 times larger for graphene due to the stronger chemical bonds. The product $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$D({\boldsymbol{q}})\overline g$$\end{document}$, which determines the plasmon decay rate via phonons, is roughly the same for both systems. We plot in Fig. [4](#Fig4){ref-type="fig"} the calculated **q**-resolved plasmon spectral function for monolayer TaS~2~. We obtain a plasmon lifetime *τ* \~ 2 ps from phonon-mediated decays, comparable to the lifetime of graphene plasmons measured in recent experiments^[@CR41]^. More importantly, because plasmons in monolayer metallic TMDs originate mostly from electron orbitals of the innermost transition metal layer, we expect that the plasmon wavefunction does not extend too much outside the structure, and hence plasmons in monolayer metallic TMDs are more robust against substrate phonons than plasmons in graphene. While our estimate of the lifetime is an upper bound for an intrinsic sample, extrinsic effects such as ripples and defects, which needs to be considered in applications, can also be partially reduced with substrate engineering and encapsulation^[@CR42]^.Fig. 4Plasmon spectral function with amplitude given in color scale for an intrinsic monolayer of TaS~2~ computed from first principles.We consider only decay processes solely from phonon-assisted scatterings, since the plasmon is outside of the Landau damping regions for the plotted range of wavevectors. The plasmon linewidth Γ becomes slightly sharper for modes with higher frequencies due to the smaller density of final states that the plasmon can decay to (see Methods section). The intrinsic plasmon linewidth decreases from 1 meV at \~0.15 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\AA^{{\mathrm{ - 1}}}$$\end{document}$ down to 0.3 meV at *q* \~ 0.3 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\AA^{{\mathrm{ - 1}}}$$\end{document}$.

Slow plasmon wave packets {#Sec5}
-------------------------

We now highlight how these finite-q, dispersionless quasi-2D plasmons can be used in experiments to create novel localized plasmon wave packets or use them for novel spectroscopic applications. We first consider a setup where plasmons are excited with an ultrafast laser pulse of duration *T* which is electromagnetically coupled to the sample through an AFM tip, commonly used to image plasmons^[@CR1],[@CR2],[@CR4],[@CR11]^. We model this external excitation from the AFM tip as an oscillating point charge with time dependence $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\rho _{{\mathrm{ext}}}\left( {{\mathrm{\Delta }}t} \right) \sim \sin \omega_0t \; e^{ - \frac{{{\mathrm{\Delta }}t^2}}{{2T^2}}}$$\end{document}$ at \~4 Å above the topmost S atoms and compute the induced charge density as a function of space and time on a monolayer TaS~2~ for $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _0\sim 1$$\end{document}$ eV and *T* \~ 80 fs. We find the resulting induced charge density moves very slowly through the material: after Δ*t* = 1 ps, the wave packet is still localized around a narrow (25 nm thick) ring with a radius of \~120 nm relative to the AFM tip (see inset of Fig. [5a](#Fig5){ref-type="fig"}). In addition, we can deduce that these plasmons are still well-defined excitations by computing the figure of merit *L*/*λ* \~ *v*~*g*~*τ*/*λ* \~ 60, where *L* is the plasmon propagation length and *λ* its wavelength, using the above discussed plasmon lifetime *τ* of \~2 ps and plasmon group velocity *v*~*g*~ \~ 0.6 Å/fs (see Methods section).Fig. 5Slowly moving plasmon wave packet excited on monolayer TaS~2~.**a** Plasmon modes are excited in monolayer TaS~2~ with an ultrafast laser pulse (energy $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _0 = 1.02$$\end{document}$ eV, modulated with a Gaussian profile of width of T = 80 fs), which is coupled to the sample through an AFM tip at \~4 Å above the topmost S atoms. The red and blue regions are schematic representations (not to scale) of the positively and negatively induced charge densities, respectively, that make up the plasmon wave packet at a time Δt = 1 ps after the external perturbation and are obtained from first principles. Inset shows a cross-sectional plot of the computed induced charge density which highlights that, 1 ps after the excitation, the plasmon wave packet only traveled \~100 nm and is still localized on a disk \~20 nm thick. **b** Field enhancement due to slow plasmons in monolayer TaS~2~. An external field is produced by an oscillating charge density, in the steady-state regime, with frequency *ω*~0~ distributed in a thin disk placed 10 Å below the bottom S layer. The field intensity enhancement, computed from first principles, is evaluated as the ratio of the intensity of the total electric field $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {E_{{\mathrm{tot}}}\left( {\mathbf{r}} \right)} \right|^2$$\end{document}$ to the maximum intensity of the external field $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {E_{{\mathrm{ext}}}^{{\mathrm{max}}}} \right|^2$$\end{document}$. The field intensity enhancement is two orders of magnitude larger and much more confined along the *z* direction perpendicular to the quasi-2D material when an external field with *ħω*~0~ = 0.86 eV excites slower and shorter-wavelength plasmons (right panel), compared to an external field with *ħω*~0~ = 0.16 eV which excites faster and longer-wavelength plasmons (left panel).

Slow plasmon wave packets excited such a way should allow one to use pump-probe timing information as a way to probe the electronic structure of complex heterostructures, since nano-imaging techniques employing AFM tips that image plasmon standing waves^[@CR1],[@CR2],[@CR4],[@CR11]^ can be adapted to measure plasmon wave packets on quasi-2D metals. For instance, a slow plasmon wave packet created at time t~1~ can bounce back from the edge of the qusi-2D metal and eventually reach the AFM tip at time t~2~; the time difference t~2 ~-- t~1~ is a function not only of the distance the plasmon traveled, but also of the plasmon group velocity, which is very sensitive to the local screening environment (see Fig. [3b](#Fig3){ref-type="fig"}). Hence, compared to typical plasmon nano-imaging techniques, time-resolved plasmon microscopy should have a very high sensitivity to changes in the local and non-local electronic structure and can therefore detect changes in the substrate supporting a quasi-2D metal and the presence of point defects and interfaces even if these features are not directly on the monolayer material.

Electric field enhancement {#Sec6}
--------------------------

Nearly flat plasmon dispersion relations also allow for a high degree of confinement of the electromagnetic radiation: our first-principles calculations show that the resulting field intensity enhancement on both sides of a monolayer TaS~2~ is extraordinarily large−exceeding 10^7^ for an external field generated by a thin disk placed 10 Å under the bottom S layer (modeling a defect or nanostructured pattern) and oscillating with a time-dependent charge density at a frequency of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _0 =$$\end{document}$ 0.86 eV (see Fig. [5b](#Fig5){ref-type="fig"} and Methods section for details). This giant field intensity enhancement is indeed originating from the strongly localized and slowly moving plasmons. The enhancement obtained at a lower frequency of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _0 =$$\end{document}$ 0.16 eV, wherein the plasmon group velocity is larger by a factor of \~20×, is only of the order of 10^4^. But more importantly, as illustrated in Fig. [5](#Fig5){ref-type="fig"}, the field enhancement can be experienced in a wide spatial region on top of the monolayer material−only possible as the enhancement originates from the intrinsic quasi-2D plasmons of the 2D materials as opposed to from Mie scattering. This would allow for much easier accessing of these plasmon-induced field enhancements in experiments for biosensing, single-molecular spectroscopy, and catalytic applications, among others. Moreover, we envision that excitations of these slowly moving plasmons with light in larger metamaterials would open the possibility of creating arrays of coupled localized bosonic excitations in extended systems, enabling the monitoring of plasmon-plasmon interactions and scattering processes largely evasive to detailed and systematic studies up to now.

Discussion {#Sec7}
==========

In summary, our ab initio calculations and effective model analysis show that plasmons in atomically thin quasi-2D metals follow a unique dispersion relation that flattens for larger wavevectors, owing to the characteristics of quasi-2D screening and broken translational symmetry in real atomically thin crystalline materials. We moreover show that monolayer metallic TMDs are especially suited to explore these nearly dispersionless plasmons, as they can support such collective excitations with long lifetime. The high tunability of these plasmons with varying substrates, as well as the flexibility to create plasmon wave packets with different group velocities, opens the possibility of exciting applications such as tracking plasmon wave packets in real time for time-resolved plasmonic imaging and supporting giant field enhancement (exceeding 10^7^) in extended quasi-2D metals over wide spatial regions.

Methods {#Sec8}
=======

New computational method for plasmon dispersion calculations {#Sec9}
------------------------------------------------------------

The key quantity required to obtain the plasmon dispersion is the proper, time-ordered polarizability matrix, which is given within the ring approximation by$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathbf{GG}}^{\prime} }^0\left( {{\mathbf{q}},\omega } \right) =  \frac{2}{{V_{{\mathrm{xtal}}}}}\mathop {\sum}\limits_{mn{\mathbf{k}}} \langle{n{\mathbf{k}}\left| {e^{ - i\left( {{\mathbf{q}} \,+\, {\mathbf{G}}} \right) \cdot {\mathbf{r}}}} \right|m{\mathbf{k}} + {\mathbf{q}}\rangle \langle m{\mathbf{k}} + {\mathbf{q}}\left| {e^{i\left( {{\mathbf{q}} \,+\, {\mathbf{G}}} \right) \cdot {\mathbf{r}}}} \right|n{\mathbf{k}}\rangle\,f_{m{\mathbf{k}} \,+\, {\mathbf{q}}}\left( {1-f_{n{\mathbf{k}}}} \right)} \\ \,   \times \left[ {\frac{1}{{\omega - \left( {E_{n,{\mathbf{k}}} - E_{m,{\mathbf{k}} \,+\, {\mathbf{q}}}} \right) + i0^ + }} - \frac{1}{{\omega + \left( {E_{n,{\mathbf{k}}} - E_{m,{\mathbf{k}} \,+\, {\mathbf{q}}}} \right) - i0^ + }}} \right],$$\end{document}$$where *V*~xtal~ is the crystal volume, *m* and *n* denote band indices, **k** and **q** are wavevectors in the Brillouin zone, and *f* is a Fermi-Dirac occupation factor. The dielectric matrix *ε* is related, within the RPA, to the proper polarizability by $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{{\mathbf{GG}}^{\prime} }\left( {{\mathbf{q}},\omega } \right) = \delta _{{\mathbf{GG}}^{\prime} } - v\left( {{\mathbf{q}} + {\mathbf{G}}} \right)\chi _{{\mathbf{GG}}^{\prime} }^0\left( {{\mathbf{q}},\omega } \right)$$\end{document}$, where *v*(**q** + **G**) is the Coulomb potential in reciprocal space and **G** is a reciprocal lattice vector.

In most ab initio codes, one typically needs to employ a uniform Monkhorst-Pack grid to sample all the **k** points for the virtual transitions involved in building *χ*^0^. Within this approach, the plasmon wavevector **q** is also typically made commensurate with the Monkhorst-Pack grid used to sample the **k** points. This leads to a serious computational challenge, as the number of **k** points *N*~*k*~ increases as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left( {\frac{1}{q}} \right)^2$$\end{document}$ for quasi-2D systems. Note that calculations on semiconductors do not exhibit this difficulty due to the presence of a gap. Here, we use an algorithm based on importance sampling to compute the polarizability matrix at arbitrarily small wavevectors, at a constant computational cost and with negligible and controllable error. First, we separate the proper polarizability matrix into intraband and interband contributions. The interband contribution, while may involve a large number of bands, can be easily sampled on a coarse **k**-point grid, since there is no Fermi surface to be resolved. The challenge is to efficiently sample the intraband contribution to *χ*^0^, which we denote by $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathrm{intra}}}^0$$\end{document}$.

In our new scheme to efficiently compute $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathrm{intra}}}^0$$\end{document}$, we sample only the set of **k** points that are involved with an intraband transition. We first distribute a set of special anchor points {**k**^*a*^} along the different pieces of the Fermi surface, in a way that: (i) **v**(**k**^*a*^)·**q** \< 0, where **v**(**k**) is the band velocity at **k** (this will later ensure that the product of the occupation factors *f*~*m***k** + **q**~(1−*f*~*n***k**~) is nonzero); and (ii) their spacing along the direction perpendicular to **q** is fixed as 1/*N*~*a*~. Note that *N*~*a*~ is not the total number of anchor points, but the maximum number of anchor points that a closed and concave Fermi surface could have. We choose *N*~*a*~ to be large enough to finely sample transitions along the Fermi energy, around *N*~*a*~ \~ 60; however, note that *N*~*a*~ does not depend on the magnitude of *q*, but on the size of the Fermi surface. We then distribute the actual **k** points that will go into the calculation of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathrm{intra}}}^0$$\end{document}$ in pairs around the anchor points, i.e., as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left\{ {{\mathbf{k}}_i = {\mathbf{k}}_i^a \pm \frac{{\mathbf{q}}}{2}} \right\}_i$$\end{document}$, so that each pair of **k** points will sample a transition along the various pieces of the Fermi surface.

We then perform the calculation of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathrm{intra}}}^0$$\end{document}$ with these special set of **k** points and by including only intraband transitions. This is enforced by keeping only transitions between states $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {|m{\mathbf{k}} + {\mathbf{q}}} \right\rangle$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {\left| {n{\mathbf{k}}} \right.} \right\rangle$$\end{document}$ with a large overlap, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {\langle {u_{m{\mathbf{k}} + {\mathbf{q}}}{\mathrm{|}}u_{n{\mathbf{k}}}}\rangle } \right|^2 \ge \frac{1}{2}$$\end{document}$. We then compute interband transitions keeping only transitions such that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {\langle {u_{m{\mathbf{k}} + {\mathbf{q}}}{\mathrm{|}}u_{n{\mathbf{k}}}}\rangle } \right|^2 < \frac{1}{2}$$\end{document}$. Since these interband transitions do not sample the Fermi surface, they can be computed with two independent set of 6×6×1 **k**-point grids shifted by the plasmon wave vector **q**. Note that, both for the interband and intraband calculations, new set of wavefunctions need to be calculated depending on the wavevector **q**. The **k**-point grids can be chosen so that the **q**-dependent wavefunctions only depend on the valence states, so they can be efficiently recomputed. However, the approximation of only sampling the intraband transitions with set of pairs of **k**-points around the Fermi surface is only valid for small values of *q*. For values of *q* larger than 1/60 of the reciprocal lattice vector, we compute the polarizability with an uniform 60×60×1 **k**-point grid for transition involving well-separated occupied and unoccupied bands which are both \~2 eV from the Fermi energy, and employing a uniform 6×6×1 **k**-point for the remaining transitions. Altogether, we can accurately and efficiently compute the polarizability matrix from first principles for all relevant wavevectors.

Finally, we obtain the plasmon excitations by finding the peaks of the loss function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$L\left( {{\mathbf{q}},\omega } \right) = - {\mathrm{Im}}\left[ {\frac{1}{{\varepsilon _m\left( {{\mathbf{q}},\omega } \right)}}} \right] = - {\mathrm{Im}}\, \varepsilon _{00}^{ - 1}\left( {{\mathbf{q}},\omega } \right)$$\end{document}$, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _m\left( {{\mathbf{q}},\omega } \right) = 1/\left[ {\varepsilon _{00}^{ - 1}({\mathbf{q}},\omega )} \right]$$\end{document}$ is the macroscopic dielectric function. Finding the poles of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{Im}}\, \varepsilon _{00}^{ - 1}({\mathbf{q}},\omega )$$\end{document}$ typically requires a very dense sampling of different frequencies *ω* around the plasmon. We avoid this by first decomposing $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$u({\mathbf{q}}) + iv({\mathbf{q}}) = 1/\varepsilon _{00}^{ - 1}({\mathbf{q}},\omega )$$\end{document}$. Both *u***(q)** and *v***(q)** are smooth, real-valued functions near the peak of the loss function, so they can be accurately interpolated to give the loss function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$L\left( {{\mathbf{q}},\omega } \right) = \frac{{v({\mathbf{q}})}}{{u({\mathbf{q}})^2 + v({\mathbf{q}})^2}}$$\end{document}$, the maximum of which we associate with a plasmon energy *ω*~*p*~. The real and imaginary parts of the **G**=**G**'=0 components of the inverse dielectric matrix, as well as the real part of the macroscopic dielectric matrix, are shown in Fig. [6](#Fig6){ref-type="fig"}. While the absolute values of dielectric matrices reported are not directly observables (i.e., they depend on the stacking spacing used in our supercell calculations, as well as broadening parameters), the poles of the dielectric matrices are not sensitive to these details, and are associated with plasmon excitations.Fig. 6Dielectric response matrices computed from first principles for monolayer TaS~2~ in a supercell arrangement with a cell length of *L* = 100 Å along the direction normal to the layer and with numerical broadening of 0.1 meV.**a** Real and (**b**) imaginary parts of the **G** = **G**' = 0 components of the inverse dielectric matrices of such dilutely stacked monolayers for three different in-plane wavevectors **q**, where **G** is a reciprocal lattice vector of the quasi-2D crystal. The imaginary part plotted in panel (**b**) is proportional to the loss function, the sharp peaks of which are associated with plasmon excitations. **c** Real part of the macroscopic dielectric matrix, the zeros of which are associated with plasmon collective excitations. While the absolute values of dielectric matrices depend on the stacking spacing and broadening parameter used in our supercell calculations, the poles in panels (**a**) and (**b**), as well as the zeros in panel (**c**), are not sensitive to such details.

GW quasiparticle band structure {#Sec10}
-------------------------------

Our mean-field DFT calculations are performed with the Quantum ESPRESSO package^[@CR43]^. We compute the quasiparticle band structure for monolayer TaS~2~ in the 2H phase, using the ab initio GW method^[@CR35]^ with the BerkeleyGW package^[@CR36]^, and using similar convergence parameters as we employed in previous studies on semiconducting monolayer TMDs^[@CR15]^. A few notable differences are: (i) we perform the calculations without resorting to any plasmon-pole model by computing the fully frequency-dependent dielectric matrix and self-energy, (ii) we include all unoccupied bands in the calculation of the polarizability matrix and self-energy, and (iii) in order to sample a very dense **k**-point grid, we use the nonuniform neck subsampling (NNS) method^[@CR44]^.

Microscopic analysis and dispersion relation of quasi-2D plasmons {#Sec11}
-----------------------------------------------------------------

We identify plasmon collective excitations with peaks in the loss function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$L\left( {{\mathbf{q}},\omega } \right) = - {\mathrm{Im}}\frac{1}{{\varepsilon _m\left( {{\mathbf{q}},\omega } \right)}} = - {\mathrm{Im}}\, \varepsilon _{00}^{ - 1}\left( {{\mathbf{q}},\omega } \right)$$\end{document}$, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _m({\mathbf{q}},\omega ) = \frac{1}{{\varepsilon _{00}^{ - 1}({\mathbf{q}},\omega )}}$$\end{document}$ is the macroscopic dielectric function, and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{00}^{ - 1}$$\end{document}$ is the **G** = **G**' = 0 component of the inverse dielectric matrix. The condition for a plasmon peak is that Re\[det *ε*\] = 0. Following the previous discussion, we can separate the proper polarizability into an interband contribution and an intraband contribution and assume for simplicity here that there is a single partially occupied band that is responsible for the intraband transitions.

The real part of the intraband contribution to the proper RPA polarizability, in the long wavelength limit and at zero temperature, is given by$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{Re}} \, \chi _{{\mathrm{intra}}}^0\left( {{\mathbf{r}},{\mathbf{r}}^{\prime} ,{\mathbf{q}},\omega } \right) =     \frac{2}{{V_{{\mathrm{xtal}}}}}\mathop {\sum }\limits_k u_{\mathbf{k}}^ \ast \left( {\mathbf{r}} \right)u_{{\mathbf{k}} \,+\, {\mathbf{q}}}\left( {\mathbf{r}} \right)u_{{\mathbf{k}} \,+\, {\mathbf{q}}}^ \ast \left( {{\mathbf{r}}^{\prime} } \right)u_{\mathbf{k}}\left( {{\mathbf{r}}^\prime } \right)\theta ( {E_{{\mathbf{k}} \,+\, {\mathbf{q}}} - E_F})\theta \left( {E_F - E_{\mathbf{k}}} \right)\\  \times \left[ {\frac{1}{{\omega - ( {E_{{\mathbf{k}} + {\mathbf{q}}} - E_{\mathbf{k}}})}} - \frac{1}{{\omega + ( {E_{{\mathbf{k}} + {\mathbf{q}}} - E_{\mathbf{k}}})}}} \right],$$\end{document}$$where *u*~**k**~(**r**) is the periodic part of the Bloch function for the band crossing the Fermi energy with wavevector **k**, and *E*~*F*~ is the Fermi energy.

Since we are interested in plasmon solutions, we know that we evaluate $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathrm{intra}}}^0$$\end{document}$ at *ω* = *ω*~*p*~ ≫ *E*~**k** + **q**~ − *E*~**k**~, the reason being that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _p \sim \sqrt q$$\end{document}$ but *E*~**k** + **q**~ − *E*~**k**~ scales like **q** in the long wavelength limit. If we define $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\rho _{\mathbf{k}}\left( {\mathbf{r}} \right) \equiv u_{\mathbf{k}}^ \ast \left( {\mathbf{r}} \right)u_{{\mathbf{k}} \,+\, {\mathbf{q}}}\left( {\mathbf{r}} \right)$$\end{document}$ and assume that we are in the long wavelength limit, we can rewrite $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathrm{intra}}}^0$$\end{document}$ in a basis-set independent way for *ω* near *ω*~*p*~ as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{Re}}\,\chi _{{\mathrm{intra}}}^0\left( {{\mathbf{q}},\omega \to \omega _p} \right) \approx \frac{{q^2\rho \rho ^T}}{{\omega ^2L_z}}\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}},$$\end{document}$$where we define$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}} \equiv \frac{1}{{2\pi ^2}}\mathop {\int}\nolimits_{{\mathrm{FS}}} {{\mathrm{d}}k\,v_k\left( {\hat q \cdot \hat v_k} \right)^2} ,$$\end{document}$$*L*~z~ is the length of the simulation supercell along the normal direction of the plane of the quasi-2D material, FS denotes an integration along the Fermi surface, **v**~*k*~ is the band velocity for a wavevector at the Fermi surface, and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}$$\end{document}$ is a **q**-dependent function that is a generalization of the ratio of the carrier density *n* to the effective mass $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${m^ \ast }$$\end{document}$.

The real part of dielectric matrix for frequencies near *ω*~*p*~ can now be written as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{Re}}\,\varepsilon = I - v\left( {\mathbf{q}} \right)\chi ^0\left( {\mathbf{q}} \right) = \varepsilon _{{\mathrm{inter}}}\left( {\mathbf{q}} \right) - v\left( {\mathbf{q}} \right)\frac{{q^2\rho \rho ^T}}{{\omega ^2L_z}}\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}},$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{{\mathrm{inter}}}\left( {\mathbf{q}} \right): = I - v\left( {\mathbf{q}} \right)\chi _{{\mathrm{inter}}}^0\left( {\mathbf{q}} \right)$$\end{document}$ is the interband contribution to the dielectric matrix, which we may safely assume is purely real in the frequency range of interest. Using Sylvester's determinant identity, the condition for a plasmon peak is equivalent to$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{Re}}[ {\det \varepsilon ( {{\mathbf{q}},\omega _p})}]    = 1 - \frac{{q^2}}{{\omega_p ^2L_z}}\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}\rho ^TW_{{\mathrm{inter}}}\left( {\mathbf{q}} \right)\rho = 0\\ W_{{\mathrm{inter}}}\left( {\mathbf{q}} \right)   \equiv \, \varepsilon _{{\mathrm{inter}}}^{ - 1}\left( {\mathbf{q}} \right)v\left( {\mathbf{q}} \right),$$\end{document}$$and *W*~inter~ is the Coulomb potential screened by the interband transitions.

In the long wavelength limit, if the charge density *ρ* is localized in the monolayer, we can approximate$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rho} ^TW_{{\mathrm{inter}}}\left( {\mathbf{q}} \right)\rho \approx {\mathop {\sum }\limits_{G_z,G^{\prime}_z}} \left[ {W_{{\mathrm{inter}}}\left( {\mathbf{q}} \right)} \right]_{G_z,G^{\prime}_z},$$\end{document}$$where *G*~*z*~ is a reciprocal-lattice vector along the confined direction of the material.

As done in the context of quasi-2D semiconductors^[@CR28]^, we define an effective quasi-2D dielectric function as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{1}{{\varepsilon _{{\mathrm{eff}}}({\mathbf{q}})}} = \frac{q}{{2\pi e^2L_z}}\mathop {\sum }\limits_{G_z,G^{\prime}_z} \left[ {W_{{\mathrm{inter}}}\left( {\mathbf{q}} \right)} \right]_{G_z,G^{\prime}_z},$$\end{document}$$and finally arrive at the expression for the plasmon frequency,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\omega _p^2 = \frac{{2\pi e^2}}{{\varepsilon _{{\mathrm{eff}}}\left( {\mathbf{q}} \right)}}\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}q.$$\end{document}$$

Comparison of first-principles calculations with analytical model {#Sec12}
-----------------------------------------------------------------

In order to assess the accuracy of the closed-form plasmon dispersion relation used in the model in Eq. ([1](#Equ1){ref-type=""}), we need to compute separately both the interband effective dielectric function, *ε*~eff~(**q**), and the effective ratio $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}$$\end{document}$ for a suspended monolayer TaS~2~. We evaluate *ε*~eff~(**q**) from first principles using the expression in the previous Methods section, while the ratio $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left[ {\frac{n}{{m^ \ast }}\left( {\mathbf{q}} \right)} \right]_{{\mathrm{eff}}}$$\end{document}$ is obtained by a numerical calculation of the plasmon dispersion relation including only intraband transitions also from first-principles. For *q* → 0, this is equivalent to the integral written in the previous Methods section, but it is easier to evaluate in practice.

Next, we obtain the substrate-dependent plasmon dispersion for monolayer TaS~2~. To do so, we first fit *ε*~eff~(**q**) obtained for a suspended monolayer TaS~2~ to the fully *q*-dependent expression for the Keldysh model^[@CR21]^,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{{\mathrm{eff}}}\left( q \right) = \frac{\varepsilon }{2}{\mathrm{sech}}\left( {\frac{{qd}}{2} + \eta _2} \right){\mathrm{sech}}\left( {\frac{{qd}}{2} + \eta _1} \right)\sinh \left( {qd + \eta _1 + \eta _2} \right),$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\eta _1 = \frac{1}{2}\log (\varepsilon + \varepsilon _s) - \frac{1}{2}\log (\varepsilon - \varepsilon _s)$$\end{document}$, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\eta _s = \frac{1}{2}\log (\varepsilon + 1) - \frac{1}{2}\log (\varepsilon - 1)$$\end{document}$, *ε*~*s*~ is the substrate dielectric constant, and *ε* and *d* are effectively fitting parameters intrinsic to the quasi-2D material that were obtained from the suspended calculation (for instance, *d* does not in general correspond to the thickness of the quasi-2D material). Note that, for small *q* and large ε, this expression is commonly approximated as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _{{\mathrm{eff}}}\left( q \right) \approx \frac{{1 \,+\, \varepsilon _s}}{2} + \rho _0q$$\end{document}$, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\rho _0 \approx \frac{{d\varepsilon }}{2}$$\end{document}$ whose value requires explicit first-principles calculation.

For TaS~2~, we obtain that a fit with *ε* = 5.52 and *d* = 10.3 Å reproduces the ab initio curve for *ε*~eff~(*q*) up to 0.4 Å^−1^. We then obtain the values of *ε*~eff~(*q*) for different substrates by varying *ε*~*s*~, which highlights the versatility of the model in Eq. ([1](#Equ1){ref-type=""}). We take values of *ε*~*s*~ from experiment ($\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _s^{{\mathrm{PTFE}}} \sim 2$$\end{document}$, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _s^{{\mathrm{hBN}}} \sim 4.5$$\end{document}$), with the exception of WSe~2~, for which we compute $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\varepsilon _s^{{\mathrm{WSe}}_2} \sim 15$$\end{document}$ from first principles. The accuracy of the resulting close-form dispersion relations is depicted in Fig. [3b](#Fig3){ref-type="fig"} of the main text.

Framework for estimating plasmon lifetime {#Sec13}
-----------------------------------------

We consider the decay rate of a plasmon with momentum **q** to an arbitrary final state *f* using Fermi's golden rule,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\Gamma \left( {{\mathbf{q}} \to f} \right) = \pi \mathop {\sum}\limits_f {\left| {M\left( {{\mathbf{q}} \to f} \right)} \right|^2\delta \left( {E_f - \omega _p\left( {\mathbf{q}} \right)} \right)} ,$$\end{document}$$where *M* is a coupling matrix element, and here we use atomic Rydberg units $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$(\hbar = 2m_e = e^2/2 = 1)$$\end{document}$.

We consider the case of the decay of a plasmon outside the Landau damping region, such as plasmons in monolayer TMDs (see Fig. [3b](#Fig3){ref-type="fig"}), so that a direct decay of a plasmon to an electron-hole pair cannot conserve both energy and momentum. Under these conditions, we need to consider decay processes that involve the emission of a phonon, i.e., in the low-temperature limit. There are three kinds of decay processes to consider that are first-order in the electron-phonon coupling matrix: (1) the direct decay of a plasmon to the electronic ground state by emission of a phonon, (2) the decay of a plasmon to another plasmon by emission of a phonon, and (3) the decay of a plasmon to an electron-hole pair plus a phonon. Process (1) is not relevant for plasmons close to the dispersionless regime, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega _p\left( {\mathbf{q}} \right)\sim 1\,{\mathrm{eV}}$$\end{document}$, since the phonon energy in monolayer TMDs is smaller by at least one order of magnitude. Process (2), while possible, has a small amplitude, since there is a tight constraint on momentum conservation because of the plasmon dispersion. Therefore, we focus on process (3), which could in principle be quite different between monolayer TMDs and graphene.

We label the final states as follows: the emitted phonon has a wavevector **Q**, branch index *λ* and frequency $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{\Omega }}_\lambda ^{{\mathrm{ph}}}({\mathbf{Q}})$$\end{document}$, while the final electron-hole pair consists of an electron with band index *c* and wavevector **k** + **q** − **Q**, and a hole with band index *v* and wavevector **k**, any combinations which satisfy wavevector conservation. We write the decay rate associated with this process as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\Gamma ({\mathbf{q}}) = \pi \mathop {\sum }\limits_{{\mathbf{Q}}\lambda } \mathop {\sum }\limits_{vc{\mathbf{k}}} \left| {M\left( {{\mathbf{q}} \to {\mathbf{Q}},\lambda ,v,c,{\mathbf{k}}} \right)} \right|^2\delta \left( {{\mathrm{\Omega }}_\lambda ^{{\mathrm{ph}}}\left( {\mathbf{Q}} \right) + {\it{\epsilon }}_{c{\mathbf{k}} \,+\, {\mathbf{q}} \,-\, {\mathbf{Q}}} - {\it{\epsilon }}_{v{\mathbf{k}}} - \omega _p\left( {\mathbf{q}} \right)} \right),$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\it{\epsilon }}$$\end{document}$~*n***k**~ is the quasiparticle energy for an electron or hole in band *n* and wavevector **k**.

The coupling matrix element *M* between an initial plasmon state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {|{\mathbf{q}}} \right\rangle$$\end{document}$ and a final state $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {|f} \right\rangle$$\end{document}$ mediated by the electron-phonon interaction Hamiltonian *H*~*ep*~ can be expressed as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$M\left( {{\mathbf{q}} \to {\mathbf{Q}},\lambda ,v,c,{\mathbf{k}}} \right) = \langle {f| {H_{ep}}|{\mathbf{q}}} \rangle$$\end{document}$$with$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$H_{ep} = N_k^{ - \frac{1}{2}}\mathop {\sum }\limits_{n\,m\,{\mathbf{k}}\,{\mathbf{Q}}\,\lambda } g_{nm\lambda }\left( {{\mathbf{k}}, - {\mathbf{Q}}} \right)c_{m{\mathbf{k}} - {\mathbf{Q}}}^\dagger c_{n{\mathbf{k}}}\left( {a_{ - {\mathbf{Q}}\lambda } + a_{{\mathbf{Q}}\lambda }^\dagger } \right),$$\end{document}$$where *N*~*k*~ is the number of **k** points in the Brillouin zone, g~*nmλ*~(**k**, **Q**) is the electron-phonon matrix element, and *a*^†^ and *c*^†^ denote the creation operator for a phonon and an electron, respectively. We write the initial plasmon state within the Tamm-Dancoff approximation in the basis of free electron-hole pairs, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {|{\mathbf{q}}} \right\rangle = \mathop {\sum }\limits_{v^{\prime} c^{\prime} {\mathbf{k}}^{\prime} } A_{v^{\prime} c^{\prime} {\mathbf{k}}^{\prime} }^{\mathbf{q}}|\left. {c^{\prime} {\mathbf{k}}^{\prime} + {\mathbf{q}},v^{\prime} {\mathbf{k}}^{\prime} } \right\rangle$$\end{document}$, where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$A_{v^{\prime} c^{\prime} {\mathbf{k}}^{\prime} }^{\mathbf{Q}}$$\end{document}$ is an expansion coefficient and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {\left| {c^{\prime} {\mathbf{k}}^\prime + {\mathbf{q}},v^\prime {\mathbf{k}}^\prime } \right.} \right\rangle$$\end{document}$ denotes a free electron-hole pair of Bloch states excited above the ground state, and we write the final state as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left. {\left| f \right.} \right\rangle = |\left. {c{\mathbf{k}} + {\mathbf{q}} - {\mathbf{Q}},v{\mathbf{k}},{\mathbf{Q}}\lambda } \right\rangle$$\end{document}$ as a *free* electron-hole pair plus a phonon. The normalization of the initial plasmon state requires that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$A^{\mathbf{q}} \sim N_k^{ - \frac{1}{2}}$$\end{document}$, so that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$M\left( {{\mathbf{q}} \to {\mathbf{Q}},\lambda ,v,c,{\mathbf{k}}} \right) \sim N_k^{ - 1}\left[ {g_{cc\lambda }\left( {{\mathbf{k}} + {\mathbf{q}}, - {\mathbf{Q}}} \right) - g_{vv\lambda }\left( {{\mathbf{k}}, - {\mathbf{Q}}} \right)} \right]$$\end{document}$.

A fully first-principle evaluation of Γ for a real quasi-2D material with multiple bands is challenging, as it involves the plasmon wavefunction and the fully band-dependent and wavevector-dependent electron-phonon coupling matrix elements. At this point, to make the calculation trackable, we approximate all this dependence into an effective electron-phonon coupling matrix element $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g$$\end{document}$, which is estimated from typical values of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {g_{mn}\left( {{\mathbf{k}},{\mathbf{Q}}} \right)} \right|$$\end{document}$ found for each material. This give us an order-of-magnitude estimate of the plasmon-phonon coupling matrix element. Note that this approximation is partially justified by the fact that we are dealing with metallic systems, as the electron-phonon coupling matrix elements are reasonably smooth and do not diverge for vanishing phonon wavevectors.

This also motivates us to define a momentum-integrated joint-density of states (MI-JDOS) as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\text{MI}}{\hbox{-}}{\text{JDOS}}\left( {\upomega} \right) = \frac{1}{{N_k^2}}\mathop {\sum }\limits_{vc{\mathbf{kk}}^\prime } \delta \left( {{\it{\epsilon }}_{c{\mathbf{k}}^\prime } - {\it{\epsilon }}_{v{\mathbf{k}}} - \omega } \right),$$\end{document}$$which is a measure of how many electron-hole pairs can be created with an energy $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\upomega}$$\end{document}$, regardless of the momentum.

We then arrive at the following simplified expression to estimate the plasmon decay rate,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathrm{\Gamma }} \sim \pi \,\bar g^2\,{\text{MI}}{\hbox{-}}{\text{JDOS}}\left( {\omega \left( {\mathbf{q}} \right) - {\mathrm{\Omega }}^{{\mathrm{ph}}}} \right),$$\end{document}$$where Ω^ph^ is a typical phonon frequency associated with $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g^2$$\end{document}$. For monolayer metallic TMDs, we simply set Ω^ph^ ≈ 0.

We compute from first principles the momentum-integrated JDOS using DFT and GW calculations and compute the electron-phonon matrix elements using density-functional perturbation-theory calculations for graphene and monolayer TaS~2~. The DFT and GW calculations are performed with the Quantum ESPRESSO^[@CR43]^ package and the BerkeleyGW^[@CR36]^ package, respectively. We find that the MI-JDOS decreases as a function of energy for monolayer materials for $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\hbar \omega \,\gtrsim\, 0.3\,{\mathrm{eV}}$$\end{document}$. At the dispersionless region, we find that the monolayer TMD $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\text{MI}}{\hbox{-}}{\text{JDOS}}\left( {\hbar {\upomega} = 1\,{\mathrm{eV}}} \right) \approx 0.1\,\left( {{\mathrm{eV}}} \right)^{ - 1}$$\end{document}$ per unit of TaS~2~. For graphene doped at *E*~*F*~ = 0.5eV, on the other hand, we find that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\text{MI}}{\hbox{-}}{\text{JDOS}}\left( {\hbar {\upomega} = 1\,{\mathrm{eV}}} \right) \approx 0.001\left( {{\mathrm{eV}}} \right)^{ - 1}$$\end{document}$ per two carbon atoms.

We estimate $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g^2$$\end{document}$ by computing the average of the electron-phonon coupling matrix elements over all possible momenta for states at the half-filled band. For monolayer TaS~2~, the values of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {g_\lambda ({\mathbf{k}},{\mathbf{q}})} \right|^2$$\end{document}$ range up to 0.035 eV^2^ with an average of \~0.0005 eV^2^, the latter of which we assign to $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g^2$$\end{document}$. For doped graphene, the values of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {g_\lambda ({\mathbf{k}},{\mathbf{q}})} \right|^2$$\end{document}$ range up to 0.26 eV^2^ with an average of \~0.01 eV^2^. Using these results for $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\bar g^2$$\end{document}$ together with the above calculated MI-JDOS, we obtain the plasmon linewidth given in Fig. [4](#Fig4){ref-type="fig"}.

The value of the plasmon lifetime $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\tau = (2\Gamma {\mathrm{/}}\hbar )^{ - 1} \sim 2\,{\mathrm{ps}}$$\end{document}$ allows us to estimate the figure of merit *L*/*λ* = *v*~*g*~*τ*/*λ*, where *L* is the plasmon propagation length, *λ* is its wavelength, and *v*~*g*~ the group velocity. For the wave packets used in the simulations to produce Fig. [5](#Fig5){ref-type="fig"}, we have that *q* \~ 0.3 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\AA^{ - 1}$$\end{document}$ and *v*~*g* ~\~ 0.6 Å/fs, so *L*/*λ* \~ 60. In fact, the plasmons should be well-defined excitations (i.e., *L*/*λ* ≥ 1) provided that their lifetimes are longer than \~35 fs.

Electric field enhancement {#Sec14}
--------------------------

We compute the electric field enhancement from first principles by computing the total electric field $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$E_{{\mathrm{tot}}}\left( {{\mathbf{q}} + {\mathbf{G}},\omega } \right) = \mathop {\sum }\nolimits_{{\mathbf{G}}^\prime } \varepsilon _{{\mathbf{G}},{\mathbf{G}}^\prime }^{ - 1}\left( {{\mathbf{q}},\omega } \right)E_{{\mathrm{ext}}}({\mathbf{q}} \,+ {\mathbf{G}}^\prime ,\omega )$$\end{document}$ given an external longitudinal field *E*~ext~ . One challenge to compute the total electric field in real space and from first principles is that one needs to perform a Fourier transform of a very sharp function, which depends on the matrix *ε*^−1^(**q**, *ω*). However, the dielectric matrices are only evaluated for a discrete number of wavevectors. We address this by writing the interacting RPA polarizability matrix close to a plasmon state and for a fixed frequency *ω*~0~ in terms of a spectral representation,$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\chi _{{\mathbf{G}},{\mathbf{G}}^\prime }\left( {{\mathbf{q}},\omega _0} \right) \approx \frac{{C_{{\mathbf{G}},{\mathbf{G}}^\prime }\left( {{\mathbf{q}}_0} \right)}}{{v_g}}\frac{1}{{q - q_0 - i\eta }},$$\end{document}$$where *q*~0~ is the plasmon wavevector corresponding to the frequency *ω*~0~, i.e., *ω*~*p*~(q~0~) = *ω*~0~, *v*~*g*~ is the plasmon group at the wavevector *q*~0~, *η* ≡ *τ*^−1^/*v*~*g*~, where *τ* is the plasmon lifetime (the results do not change for *τ* ranging from 1 to 10 ps), and *C*~**G,G′**~(**q**) is an Hermitian matrix that can be fit from the ab initio inverse dielectric matrix. Note that the field enhancement in real space, which depends on the Fourier transform of *χ*~**G,G′**~(**q**, *ω*~0~) over its spatial directions, gets enhanced for slower plasmons.

Because of the mismatch in their wavelengths, plasmons in quasi-2D extended materials will not directly couple to light, and thus require a nanostructure with a length scale comparable to that of the wavelength of the plasmon to be excited. This is modeled here as an infinitesimally thin disk with a radius of 10 Å. Because the total field still depends considerably on the geometry of this nanostructure, we define the external field as the longitudinal field created by this nanostructure, which is assumed to be generated by an oscillating charge density uniformly distributed over the disk. We compute the field intensity enhancement as the ratio of the intensity of the total electric field $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {E_{{\mathrm{tot}}}\left( {\mathbf{r}} \right)} \right|^2$$\end{document}$ to the maximum intensity of the external field $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {E_{{\mathrm{ext}}}^{{\mathrm{max}}}} \right|^2$$\end{document}$ generated by this nanostructure. This definition is likely a lower bound to the field enhancement $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {E_{{\mathrm{tot}}}\left( {\mathbf{r}} \right)} \right|^2/I_{{\mathrm{light}}}$$\end{document}$ defined with respect to the intensity *I*~light~ of an external planewave light source, since the nanostructure can also be engineered to enhance the intensity $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\left| {E_{{\mathrm{ext}}}^{{\mathrm{max}}}} \right|^2/I_{{\mathrm{light}}}$$\end{document}$.

**Peer review information** *Nature Communications* thanks Jacob Khurgin and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This work was supported by the Center for Computational Study of Excited-State Phenomena in Energy Materials (C2SEPEM) funded by the U.S. Department of Energy, Office of Basic Energy Sciences under Contract No. DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory, as part of the Computational Materials Sciences Program, which provided for theory development, code implementation, and calculations. Computational resources were provided by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation under Grant No. ACI-1053575. We acknowledge financial support from the European Research Council (ERC-2015-AdG-694097). The Flatiron Institute is a division of the Simons Foundation. The authors thank D. Basov, D.Y. Qiu, and H.S. Sen for helpful discussions.

S.G.L. and A.R. conceived the research. S.G.L. supervised the research and the new formalisms developed. F.H.J. developed the methods/formalisms and performed the DFT and GW studies. L.X. performed DFT calculations. All authors contributed to the analyses of results and writing of the paper.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Ab initio calculations were performed with the BerkeleyGW software package^[@CR36]^. The new algorithms proposed here and the post processing scripts that support the findings of this study are available from the corresponding authors upon reasonable request.

The authors declare no competing interests.
