Local-field effects on the plasmon dispersion of two-dimensional transition metal dichalcogenides

The isolation of graphene in 2004 [1] gave a new twist to the field of two-dimensional (2D) crystals and boosted new avenues in condensed matter physics [2–6]. The unconventional electronic properties of graphene, that can be described in terms of massless 2D Dirac particles, have attracted the enthusiastic attention of the scientific community [7, 8] and promise a large variety of exciting technological applications [9]. However, graphene lacks a bandgap, preventing its use in electronic devices like transistors. Hence, finding the best way to induce a gap in graphene has been a very active research topic in the past few years. This can be obtained, for example, by chemical functionalization as for graphane and other graphene derivatives [10–13], although at the price of losing some of the unconventional properties of pristine graphene. To overcome these difficulties while keeping the advantages of the reduced dimensionality, the research activity is rapidly extending to other families of layered materials. In particular, transition-metal dichalcogenides (TMDs) [14] are increasingly gaining attention [4, 5, 15]. Their chemical formula is MX₂, where M is a transition metal and X a chalcogen (S, Se and Te). In contrast to graphene, these materials show versatile electronic properties that vary from metallic to insulating, depending on the metal M [16].

TMDs can be shaped into monolayers, displaying distinct physical properties from their bulk counterpart. MoS₂, for example, has an indirect bandgap in bulk form, but in 2D, the bandgap increases and becomes direct [17, 18], being most favorable for optoelectronic applications. A 2D MoS₂ transistor has also recently been realized [19], demonstrating that TMDs can be used for flexible electronic devices. MoS₂ has shown to be promising also for spintronics applications [20], absent in graphene due to its vanishing spin–orbit coupling, and in 'valleytronics', exploiting the possibility of controlling the population of the valleys (i.e. conduction band minima) [21, 22].

Metallic TMDs such as NbSe₂ and TaS₂ exhibit remarkable low-temperature phenomena, including the competition between superconductivity and charge-density wave order [23]. In its single-layer form, the hope is that it can be efficiently used for plasmonics applications with the aim of controlling light on the nanoscale [24]. This typically requires surface plasmons that can confine the electromagnetic field at scales smaller than the wavelength of light. In fact, in comparison to conventional metals like gold or silver, 2D materials have a significant advantage: their plasmonic properties can be easily controlled by external means (e.g. doping or electrical voltage) [25–27].

In view of these potentialities, it is very timely to address the study of the elementary electronic excitations of this novel class of 2D materials on the basis of predictive and reliable first-principles methods. Hence here we make use of cutting-edge approaches based on time-dependent density-functional theory (TDDFT) [28] to take into account the details of the electronic interactions at play. We calculate the dispersion (i.e. the dependence on the momentum transfer) of plasmons in two prototypical TMDs: NbSe₂ and TaS₂. To disclose the effects of the reduced dimensionality, we compare the plasmon properties in the bulk [29–31] and in the single-layer forms, showing the prominent role played by the induced microscopic electronic charges (the 'local fields') in 2D. This suggests the possibility to effectively tune the dielectric properties of these novel materials for plasmonics applications.

Plasmon excitations can be identified by the peaks in the loss function L(q,ω) = −Im⁻¹_M(q,ω), which can be written in terms of the real and imaginary parts of the macroscopic dielectric function _M = ₁ + i₂:

$$\begin{equation} L(\mathbf{q},\omega)=-\textrm{Im} \frac{1}{\epsilon_{\rm M}(\mathbf{q},\omega)}=\frac{\epsilon_2(\mathbf{q},\omega)}{[\epsilon_1(\mathbf{q},\omega)]^2+[\epsilon_2(\mathbf{q},\omega)]^2}, \end{equation} \tag{ 1 }$$

where q is the momentum transfer. Plasmons energies ω_p(q) correspond to zeroes of ₁(q,ω_p) where the damping given by ₂(q,ω_p) is not too large, so that _M ≃ 0. In fact, this represents the condition that must be satisfied by an electronic system to sustain longitudinal collective excitations (i.e. the induced electric field is finite even in the absence of external sources) [32].

The macroscopic dielectric function _M can be calculated from first-principles, using methods based on Green's functions or within TDDFT [28], in which the Dyson-like equation for the linear-response polarizability χ reads

$$\begin{eqnarray}&&\fl \chi(\mathbf{r}_1,\mathbf{r}_2,\omega) = \chi_{\rm KS}(\mathbf{r}_1,\mathbf{r}_2,\omega)\nonumber \\ &&+ \int {\rm d}\mathbf{r}_3\,{\rm d}\mathbf{r}_4 \,\, \chi_{\rm KS}(\mathbf{r}_1,\mathbf{r}_3,\omega) [v(\mathbf{r}_3,\mathbf{r}_4)+f_{\rm xc}(\mathbf{r}_3,\mathbf{r}_4,\omega)]\chi(\mathbf{r}_4,\mathbf{r}_2,\omega), \end{eqnarray} \tag{ 2 }$$

where χ_KS is the Kohn–Sham non-interacting polarizability that can be expressed in terms of the Kohn–Sham energies _kii and orbitals ϕ_kii(r) as

$$\begin{eqnarray} &&\fl \chi_{\rm KS}(\mathbf{r}_1,\mathbf{r}_2,\omega) = 2\sum_{\mathbf{k}_1\mathrm{i}\mathbf{k}_j{\rm j}}\phi_{\mathbf{k}_i{\rm i}}(\mathbf{r}_1)\phi^*_{\mathbf{k}_j{\rm j}}(\mathbf{r}_1)\phi^*_{\mathbf{k}_i{\rm i}}(\mathbf{r}_2)\phi_{\mathbf{k}_j{\rm j}}(\mathbf{r}_2) \nonumber \\ &&\times \left[\frac{\theta(\epsilon_{\mathbf{k}_i{\rm i}}-\mu)\theta(\mu-\epsilon_{\mathbf{k}_j{\rm j}})}{\omega-(\epsilon_{\mathbf{k}_i{\rm i}}-\epsilon_{\mathbf{k}_j{\rm j}})+{\rm i}\eta}-\frac{\theta(\mu-\epsilon_{\mathbf{k}_i{\rm i}})\theta(\epsilon_{\mathbf{k}_j{\rm j}}-\mu)}{\omega-(\epsilon_{\mathbf{k}_i{\rm i}}-\epsilon_{\mathbf{k}_j{\rm j}})-{\rm i}\eta}\right], \end{eqnarray} \tag{ 3 }$$

v the Coulomb interaction and f_xc the exchange-correlation TDDFT kernel. In the calculations, we adopted the local-density approximation (LDA) for χ_KS and the random-phase approximation (RPA) for χ, which corresponds to setting f_xc = 0 in equation (2). From χ, one obtains the microscopic dielectric function since

$$\begin{equation} \epsilon^{-1}(\mathbf{r}_1,\mathbf{r}_2,\omega)=\delta(\mathbf{r}_1-\mathbf{r}_2) + \int {\rm d}\mathbf{r}_3 v(\mathbf{r}_1,\mathbf{r}_3)\chi(\mathbf{r}_3,\mathbf{r}_2,\omega) \end{equation} \tag{ 4 }$$

and the macroscopic dielectric function _M through [33, 34]

$$\begin{equation} \epsilon_{\rm M}(\mathbf{q},\omega) = \frac{1}{\epsilon^{-1}_{\mathbf{G}=0,\mathbf{G}'=0}(\mathbf{q},\omega)}, \end{equation} \tag{ 5 }$$

where G and G' are reciprocal-lattice vectors and ⁻¹_G,G'(q,ω) is the Fourier transform to reciprocal space of ⁻¹(r₁,r₂,ω) (with the momentum transfer q belonging to the first Brillouin zone). The loss function L(q,ω) is thus

$$\begin{equation} -\textrm{Im} \epsilon^{-1}_{\rm M}(\mathbf{q},\omega) = -v_0(\mathbf{q})\chi_{0,0}(\mathbf{q},\omega), \end{equation} \tag{ 6 }$$

where the G and G' indexes are understood. In the homogeneous electron gas (HEG), the microscopic dielectric function is diagonal in G and G' and _M(q,ω) = _G=0,G'=0(q,ω). In a system that is inhomogeneous and polarizable on the microscopic scale, instead, off-diagonal elements of are not zero and they all contribute to the head element (G = 0,G' = 0) of ⁻¹, reflecting the so-called crystal local-field effects (LFE).

Alternatively, the macroscopic dielectric function can be directly calculated from [28]

$$\begin{equation} \epsilon_{\rm M}(\mathbf{q},\omega) = 1- v_0(\mathbf{q})\bar{\chi}_{0,0}(\mathbf{q},\omega), \end{equation} \tag{ 7 }$$

where the modified polarizability $\bar \chi$ satisfies a Dyson equation like (2) with a modified Coulomb interaction $\bar v$ that in the reciprocal space is equal to v for all the G components but G = 0, for which $\bar v(\mathbf {G}=0)=0$

$$\begin{eqnarray} &&\fl \bar \chi(\mathbf{r}_1,\mathbf{r}_2,\omega) = \chi_{\rm KS}(\mathbf{r}_1,\mathbf{r}_2,\omega) \nonumber \\ &&+ \int {\rm d}\mathbf{r}_3\,{\rm d}\mathbf{r}_4 \,\, \chi_{\rm KS}(\mathbf{r}_1,\mathbf{r}_3,\omega) [\bar v(\mathbf{r}_3,\mathbf{r}_4)+f_{\rm xc}(\mathbf{r}_3,\mathbf{r}_4,\omega)]\bar \chi(\mathbf{r}_4,\mathbf{r}_2,\omega). \end{eqnarray} \tag{ 8 }$$

In the RPA with f_xc = 0, setting also $\bar v=0$ in the Dyson equation (8) means neglecting the LFE in the calculation of _M, since in this case $\bar \chi = \chi _{\mathrm { KS}}$. Instead, with the inclusion of the LFE, χ_KS is screened by the short-range charge fluctuations associated with the induced Hartree potential.

Plasmon were first studied in the HEG [35], where the plasmon energy evaluated in RPA has a positive dispersion, which is parabolic as a function of q in three-dimensional (3D) and behaves as $\sqrt {q}$ in 2D [36]. The HEG model, which is characterized by a single parabolic band, is often taken as a reference system also to study and interpret the collective excitations in real metallic materials. However, in realistic materials, deviations from the HEG RPA behavior can have different origins: a different electronic structure (i.e. non-parabolicity), the effect of oscillator-strength matrix elements, non-homogeneity, the presence of interband transitions and correlation effects beyond RPA (i.e. due to f_xc in equation (2)). In fact, results at variance with HEG, including negative dispersions, have been reported in many cases, as for example in heavy alkali metals like Cs [37, 38], doped molecular crystals [39] and 3D TMDs [29–31].

Therefore, especially in view of the negative plasmon dispersion found in the bulk materials, it makes sense to address here in detail the investigation of plasmons in the 2D TMDs using a microscopic approach based on a first-principles calculation of the realistic band structure beyond the HEG model.

In our approach, the dielectric function has been evaluated on top of LDA calculations implemented in a plane-wave-basis framework using norm-conserving Troullier–Martins pseudopotentials. We used Abinit [40] for the ground-state calculations and Yambo [41] for the RPA spectra. All the systems were simulated using a 2D triangular lattice with a basis of one MX₂ unit. We used a theoretical lattice constant (3.37 Å for TaS₂ and 3.38 Å for NbSe₂), calculated using a 10 × 10 × 1 k point mesh and an energy cutoff for the plane-wave expansion between 30 and 40 Ha. Moreover, in the supercell approach, the distance between adjacent layers was set to 40 Å. For the calculation of the dielectric function, we used a 30 × 30 × 1 grid of k point with 80 band. Finally, LFE have been included inverting a matrix of rank 300 G vectors.

In figure 1, we show the calculated low-energy RPA loss function of the single layers of TaS₂ and NbSe₂ as a function of the momentum transfer q along the ΓM direction⁶. In the following, we focus on TaS₂, but the same analysis holds also for NbSe₂, which has an electronic band structure that is similar to TaS₂ (see figure 2). The main peak in the spectrum, which for TaS₂ at q = 0.07 Å⁻¹ is located at ω = 0.94 eV, is related to a zero of the real part of the dielectric function ₁ (see the dashed green curve in figure 5). We thus identify this main feature as a plasmon resonance that derives from electron–hole pairs belonging to the d band crossing the Fermi energy (see figure 2). In fact, these intraband transitions give rise to a peak of ₂ (at ω = 0.70 eV for this small q, see the solid green curve in figure 5) just below the plasmon frequency. This implies that the plasmon is a collective excitation involving the charge density related to this metallic band. Therefore, it has the same nature as the charge–carrier plasmon of 3D bulk systems [30, 31].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By increasing q, we observe a positive dispersion of the plasmon peak, which loses intensity and becomes broader, entering the particle–hole continuum at q = 0.35 Å⁻¹ (i.e. where the loss function gets to coincide to ₂). The plasmon bandwidth is rather small, only ∼0.1 eV, implying that the plasmon is quite localized in real space. This positive dispersion, which is common to all the metallic single layers of the 2D family, is opposite to the bulk TMDs, where it has recently shown that the plasmon has a negative dispersion [29–31]. The similarity between the electronic band structure in 2D and 3D TMDs suggests that the opposite slope of the plasmon dispersion in the bulk and in the monolayer must be a direct consequence of the different dimensionality, as we are now going to show.

The non-interacting polarizability χ_KS can be split into the sum of an intraband contribution χ^intra_KS from the band crossing the Fermi level and an interband term χ^inter_KS from all the other electron–hole transitions: χ_KS = χ^intra_KS + χ^inter_KS. In the first stage, we calculate the 'bare plasmon', considering only intraband transitions and neglecting the LFE. The plasmon is thus determined only by the details of the band structure and the matrix elements of the oscillator strengths entering χ^intra_KS. As in the bulk, and in contrast to the HEG, we find that the intrinsic dispersion of the bare plasmon is negative (see the solid black line in figure 3). This is a consequence of the presence of non-dispersive d bands at the Fermi level both in 3D and 2D [30], which makes these materials very different from the HEG, as we now discuss in more detail for the 2D case.

Zoom In Zoom Out Reset image size

In the HEG, the positive dispersion of the plasmon frequency can simply be understood in terms of its band structure. The presence of a wide parabolic energy band gives rise to a sharp peak in the joint density of states (JDOS) with strong positive dispersion with q (see the dashed lines figure 4). Since in the 2D HEG ₂ is proportional to the JDOS, ₂(q,ω)∝JDOS(q,ω)/|q|, this behavior is reflected, through the Kramers–Kronig relations, on the real part of and thus on the plasmon causing the $\sqrt {q}$ dependence. In the real 2D TMDs, in contrast to the 2D HEG, the presence of a narrow d band, gives rise to an intense and non-dispersive peak in the JDOS (see the solid lines in figure 4). Since the numerator of χ^intra_KS (see equation (3)) does not play an essential role in these materials, this results in a negative bare plasmon dispersion (see the solid black line in figure 3), as has been shown for 3D TMDs in [30]. Therefore, we conclude that, even though the HEG is often taken as a reference system to describe collective excitations in metals, a picture based on the HEG model is not valid for the TMDs either in 3D or in 2D. In fact, also in the full calculation, the plasmon dispersion does not follow the behavior of the 2D HEG in the range of q studied.

Zoom In Zoom Out Reset image size

The main difference between the bulk and the single layer is given by the effect of the crystal local fields. In the bulk TMDs, LFE do not affect the loss function for in-plane momentum transfers: the plasmon dispersion remains negative also when LFE are taken into account. On the contrary, LFE play a key role in the single layers that are characterized by an intrinsically strong inhomogeneity. LFE are the result of the mixing of different electron–hole transitions. We analyze them in two steps: first taking into account only the intraband contribution and then including the effect of the interband transitions. Therefore, in figure 5, we compare the imaginary part of the dielectric function ₂ (see the solid lines) in three different approximations: (i) including only intraband transitions without LFE (black lines); (ii) including only intraband transitions with LFE (red lines); and (iii) the full calculation adding also interband transitions (green lines). We see that the main effect of local fields on ₂ derive from the intraband transitions (compare the black and red lines in figure 5). It consists of a blueshift (by about 0.8 eV) and a strong damping (by more than 80%) of the peak of ₂. This depolarization effect on ₂ is due to the repulsive nature of $\bar v$ in equation (8). From the expression of $\bar {\chi }$ in equation (8), we can also understand that LFE are important when χ_KS is large (i.e. the system is strongly polarizable), when the electronic wavefunctions are strongly localized and the system is inhomogeneous. In fact, a similar depolarization effect has been observed, for example, for the π plasmon in graphene [43]. In 2D TMDs, we find an even larger effect of local fields. While graphene is characterized by a zero bandgap and π states, the metallic band in 2D TMD is related to d states that are more strongly localized. These depolarization effects on ₂ reflect also on ₁ through the Kramers–Kronig relations (see the dashed lines in figure 5). In fact, due to the suppression of ₂, the zero of ₁, which sets the bare plasmon energy, gets closer to the continuum of electron–hole transitions (i.e. the peak in ₂).

Zoom In Zoom Out Reset image size

Thus, when LFE are taken into account, the plasmon is blueshifted (see figure 3) and damped by the independent-particle motion. LFE cause a strong positive dispersion of the intraband peak in ₂ at small q (see the red dashed line in figure 3). Moreover, due to the strong depolarization, the plasmon follows the behavior of the electron–hole continuum (compare the dashed and solid red lines in figure 3). In fact, in a localized system, where LFE are very large, the long-range part of the Coulomb interaction, which is the difference between v and $\bar v$, becomes negligible and $\bar \chi$ and χ become very similar. In turn, since χ gives the loss function −Im⁻¹_M (see equation (6)) and $\bar \chi$ gives Im_M (see equation (7)), this explains why at small q the plasmon follows the electron–hole continuum and then becomes completely damped.

We now consider the contribution of interband transitions by adding χ^inter_KS to χ^intra_KS and solving the Dyson equation for χ to include LFE. As we can see in figure 5 (green line), the effect of interband transitions is a redshift of the continuum of electron–hole transitions. As a consequence, the zero of ₁ and hence the plasmon frequency results also redshifted by about 0.7 eV. Moreover, the positive plasmon dispersion, found by only taking into account the intraband transitions and LFE, is strongly reduced (see the green line in figure 3) with an increase of the damping.

To better understand the interplay between the interband transitions and the LFE on the plasmon dispersion, it is useful to introduce the intraband polarizability χ^intra, which, dropping for the moment the spatial indexes for simplicity, is defined as

$$\begin{equation} \chi^{\rm intra}(\omega)=\chi^{\rm intra}_{\rm KS}(\omega)+\chi^{\rm intra}_{\rm KS}(\omega)W_{\rm inter}(\omega)\chi^{\rm intra}(\omega), \end{equation} \tag{ 9 }$$

where W_inter(ω) = ⁻¹_inter(ω)v is the Coulomb interaction that is dynamically screened by the interband transitions. Then the total polarizability χ can be written as

$$\begin{equation} \chi(\omega)=\epsilon_{\rm inter}^{-1}(\omega) \{\chi^{\rm intra}(\omega)+[1-\chi^{\rm intra}_{\rm KS}(\omega)W_{\rm inter}(\omega)]^{-1}\chi^{\rm inter}_{\rm KS}(\omega)\}. \end{equation} \tag{ 10 }$$

In absence of interband transitions, χ^inter_KS = 0 and _inter = 1. Thus χ reduces to the intraband response function χ^intra evaluated with the unscreened Coulomb potential, since W_inter = v. At frequencies around the plasmon energy, χ^inter_KS is still very small and equation (10) simplifies to

$$\begin{equation} \chi(\omega) \simeq \epsilon_{\rm inter}^{-1}(\omega) \chi^{\rm intra}(\omega). \end{equation} \tag{ 11 }$$

In this energy range, where _inter is real, we thus obtain for the loss function the following expression (see equation (6)):

$$\begin{eqnarray} &&\fl -\textrm{Im} \epsilon^{-1}_{\rm M}(\mathbf{q},\omega) = -v_0(\mathbf{q}) [\epsilon_{\rm inter}(\mathbf{q},\omega)]_{00}^{-1}\, \textrm{Im}\chi^{\rm intra}_{00}(\mathbf{q},\omega) -v_0(\mathbf{q}) \sum_{\mathbf{G}\neq 0}[\epsilon_{\rm inter}(\mathbf{q},\omega)]_{0\mathbf{G}}^{-1}\,\textrm{Im}\chi^{\rm intra}_{\mathbf{G}0}(\mathbf{q},\omega).\nonumber\\ \end{eqnarray} \tag{ 12 }$$

By retaining only the zero-order contribution of the off-diagonal terms of χ^intra from equation (9), we finally arrive at the following expression for the loss functions:

$$\begin{eqnarray}&&\fl -\textrm{Im} \epsilon^{-1}_{\rm M}(\mathbf{q},\omega) = -v_0(\mathbf{q})\, [\epsilon_{\rm inter}(\mathbf{q},\omega)]_{00}^{-1}\, \textrm{Im}\chi^{\rm intra}_{00}(\mathbf{q},\omega) \nonumber\\ &&-v_0(\mathbf{q}) \sum_{\mathbf{G}\neq 0}[\epsilon_{\rm inter}(\mathbf{q},\omega)]_{0\mathbf{G}}^{-1}\,[\textrm{Im}\chi^{\rm intra}_{\rm KS}]_{\mathbf{G}0}(\mathbf{q},\omega) + \cdots \end{eqnarray} \tag{ 13 }$$

From equation (13), it is thus clear that interband transitions have two effects. First interband transitions screen χ^intra causing a redshift of the plasmon and the continuum of the electron–hole transitions. This is the effect of the first term in equation (13). The second term in equation (13) is instead related to the LFE: it is responsible for summing the off-diagonal part of χ^intra_KS, dominated by single-particle excitations, to the diagonal part of χ^intra, dominated by collective excitations. This term causes a redshift of the spectrum and makes the plasmon damping stronger.

Moreover, as q increases, the screening from interband transitions becomes weaker, reducing its redshift effect on the plasmon energy. As a result, the positive dispersion of the plasmon becomes stronger. However, at the same time, as q increases, LFE contained in the second term of equation (13) become more and more important, enhancing the coupling between the plasmon and the independent electron–hole transitions at lower energy. In turn, this tends to reduce the positive dispersion of the plasmon. We thus see that the interplay between the interband transitions and the LFE induces two competing effects on the dispersion of the intraband plasmon, which manifest differently in 2D and 3D structures.

In the bulk 3D TMDs, the LFE are negligible and the first term of equation (13) is always dominant. As a consequence, the resulting effect of the interband transitions is to temper the negative dispersion of the bulk intraband plasmon. On the other hand, in 2D systems, due to the large LFE, the second term in equation (13) becomes the most prominent contribution. In this case, interband transitions reduce the positive dispersion of the intraband plasmon (compare the red and green lines in figure 3).

In summary, we have calculated the charge–carrier plasmon dispersion in two prototypical single-layer metallic TMD: TaS₂ and NbSe₂. As in the bulk, the non-dispersive metallic band prevents using the HEG as a model to describe the plasmon excitations in these materials. Our analysis has shown that the plasmon dispersion of 2D TMDs can be even more easily tuned than their 3D counterparts. In the bulk, the plasmon dispersion is mainly determined by the JDOS crossing the Fermi level. The negative dispersion in the bulk can be switched to positive by doping with electrons or holes [30]. In the monolayers, instead, the plasmon character is also the result of the interplay between the LFE and the interband transitions, which gives rise to a small bandwidth corresponding to a localized plasmon. By also acting on the interband transitions, it is, in principle, possible to tune the plasmon nature through the microscopic charge response that constitute the strong 'local fields' acting on these materials. Therefore, metallic 2D TMDs seem a promising platform to investigate possible applications in nanoplasmonics and demand a deeper consideration.

We acknowledge financial support from the European Research Council Advanced Grant DYNamo (ERC-2010-AdG-267374), Spanish Grants (2010-21282-C02-01 and PIB2010US-00652), Grupos Consolidados UPV/EHU del Gobierno Vasco (IT578-13) and European Commission projects CRONOS (grant number 280879-2 CRONOS CP-FP7). Computational time was granted by BSC Red Espanola de Supercomputacion and GENCI (project number 544).