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

Two-dimensional transition-metal dichalcogenides (TMDs) are gaining increasing attention as alternative to graphene for their very high potential in optoelectronics applications. Here we consider two prototypical metallic 2D TMDs, NbSe$_2$ and TaS$_2$. Using a first-principles approach, we investigate the properties of the localised intraband $d$ plasmon that cannot be modelled on the basis of the homogeneous electron gas. Finally, we discuss the effects of the reduced dimensionality on the plasmon dispersion through the interplay between interband transitions and local-field effects. This result can be exploited to tune the plasmonic properties of these novel 2D materials.


Introduction
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,3,4,5,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 chemically functionalisation as for graphane and other graphene derivatives [10,11,12,13], although at the price of loosing 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 [15,4,5]. Their chemical formula is MX 2 , where M is a transition metal and X a chalcogen (S, Se,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 2 , 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 2 transistor has also recently realized [19], demonstrating that TMDs can be used for flexible electronic devices. MoS 2 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 2 and TaS 2 exhibit remarkable low-temperature phenomena including the competition between superconductivity and charge-density wave (CDW) order [23]. In the single-layer form the hope is that they 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,26,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 2 and TaS 2 . To disclose the effects of the reduced dimensionality, we compare the plasmon properties in the bulk [29, 30,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.

Methods
Plasmon excitations can be identified by the peaks in the loss function L(q, ω) = −Imǫ −1 M (q, ω), which can be written in terms of the real and imaginary parts of the macroscopic dielectric function ǫ M = ǫ 1 + iǫ 2 : where q is the momentum transfer. Plasmons energies ω p (q) correspond to zeroes of ǫ 1 (q, ω p ) where the damping given by ǫ 2 (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 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: where χ KS is the Kohn-Sham non-interacting polarizability that can be expressed in terms of the Kohn-Sham energies ǫ k i i and orbitals φ k i i (r) as: 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 randomphase approximation (RPA) for χ, which corresponds to setting f xc = 0 in Eq. (2). From χ one obtains the microscopic dielectric function ǫ since: and the macroscopic dielectric function ǫ M through [33,34]: where G and G ′ are reciprocal-lattice vectors and ǫ −1 G,G ′ (q, ω) is the Fourier transform to reciprocal space of ǫ −1 (r 1 , r 2 , ω) (with the momentum transfer q belonging to the first Brillouin zone). The loss function L(q, ω) is thus: (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 ǫ −1 , reflecting the so-called crystal local-field effects (LFE). Alternatively, the macroscopic dielectric function can be directly calculated from [28]: where the modified polarizabilityχ satifies a Dyson equation like (2) with a modified Coulomb interactionv that in the reciprocal space is equal to v for all the G components but G = 0, for whichv(G = 0) = 0: In the RPA with f xc = 0 setting alsov = 0 in the Dyson equation (8) means neglecting the LFE in the calculation of ǫ M , since in this caseχ = χ KS . Instead, with the inclusion of the LFE, χ KS is screened by the short-range charge fluctuations associated to the induced Hartree potential. Plasmon were first studied in the homogeneous electron gas [35], where the plasmon energy evaluated in RPA has a positive dispersion, which is parabolic as a function of q in 3D and behaves as √ q in 2D [36]. The HEG model, which is characterized by a single parabolic band, is often taken as 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: different electronic structure (i.e. non-parabolicity), effect of oscillator-strenght matrixelements, non-homogeneity, presence of interband transitions, correlation effects beyond RPA [i.e. due to f xc in Eq.
(2)]. In fact, results at variance with HEG, including negative dispersions, have been reported in many cases, as for example in heavy alkali metal like Cs [37,38], doped molecular crystals [39] and 3D TMDs [29, 30,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.

Computational details
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 have used Abinit [40] for the ground-state calculations and Yambo [41] for the RPA spectra. All the systems have been simulated using a two-dimensional triangular lattice with a basis of one MX 2 unit. We have used the theoretical lattice constant (3.37Å for TaS 2 and 3.38Å for NbSe 2 ) 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 has been set to 40Å. For the calculation of the dielectric function we have used a 30 × 30 × 1 grid of k point with 80 band. Finally, local-field effects have been included inverting a matrix of rank 300 G vectors. In Fig. 1 we show the calculated low-energy RPA loss function of the single layers of TaS 2 and NbSe 2 as a function of the momentum transfer q along the ΓM direction ‡. In the following we focus on TaS 2 , but the same analysis holds also for NbSe 2 , which has an electronic band structure that is similar to TaS 2 (see Fig. 2). The main peak in ‡ For the high-energy part of the spectrum at q = 0, comprising the π and π + σ plasmons, we refer to Ref. [42]. However, those calculations do not take into account LFE, which are negligible in the bulk, but are the dominant contribution in 2D, as we discuss in the following. On the other hand at q = 0 the intra-band plasmon is expected to have zero energy for dimentionality reasons. the spectrum, which for TaS 2 at q = 0.07Å −1 is located at ω = 0.94 eV is related to a zero of the real part of the dielectric function ǫ 1 (see dashed green curve in Fig. 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 Fig. 2). In fact, these intraband transitions give rise to a peak of ǫ 2 (at ω = 0.70 eV for this small q, see solid green curve in Fig. 5) just below to 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 the 3D bulk systems [30,31].

Results and discussion
By increasing q we observe a positive dispersion of the plasmon peak, which looses intensity and becomes broader, entering the particle-hole continuum at q = 0.35Å −1 (i.e. where the loss function gets to coincide to ǫ 2 ). The plasmon bandwidth is rather small, only ∼ 0.1 eV, implying that the plasmon is quite localised in real space. This positive dispersion, that is common to all the metallic single layers of the 2H family, is opposite to the bulk TMDs, where it has recently shown that the plasmon has a negative dispersion [29, 30,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 strenghts 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 solid black line in Fig. 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 more in details for the 2D case. In the HEG the positive dispersion of the plasmon frequency can be simply 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 (JODS) with strong positive dispersion with q (see dashed lines Fig. 4). Since in the 2D HEG ǫ 2 is proportional to the JDOS, ǫ 2 (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 √ 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 solid lines in Fig. 4). Since the numerator of χ intra KS [see Eq. (3)] does not play an essential role in these materials, this results in a negative bare plasmon dispersion (see solid black line in Fig. 3), as it has been shown for the 3D TMDs in Ref. [30]. Therefore we conclude that, even though the HEG is often taken as 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 behaviour of the 2D HEG in the range of q studied.
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 intrinsic strong inhomogeneity. LFE are the result of the mixing of different electron-hole transitions. We analyse them in two steps: first taking into account only the intraband contribution and then including the effect of the interband transitions. Therefore in Fig. 5 we compare the imaginary part of the dielectric function ǫ 2 (see solid lines) in three different approximations: (i) including only intraband transitions without LFE (black lines); (ii) including only intraband transitions with LFE (red lines); (iii) the full calculation adding also interband transitions (green lines). We see that the main effect of local fields on ǫ 2 derive from the intraband transitions (compare black and red lines in Fig. 5). It consists into a blueshift (by about 0.8 eV) and a strong damping (by more than 80%) of the peak of ǫ 2 . This depolarization effect on ǫ 2 is due to the repulsive nature ofv in Eq. (8). From the expression ofχ in Eq. (8) we can also understand that LFE are important when χ KS is large (i.e. the system is strongly polarizable) and when the electronic wavefunctions are strongly localized and the system is inhomogeneous. In fact, a similar depolarization effect has been observed for example also 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 ǫ 2 reflect also on ǫ 1 through the Kramers-Kronig relations (see dashed lines in Fig. 5). In fact, due to the suppression of ǫ 2 , the zero of ǫ 1 that sets the bare plasmon energy gets closer to the continuum of electron-hole transitions (i.e. the peak in ǫ 2 ).
Thus when LFE are taken into account the plasmon is blueshifted (see Fig. 3) and damped by the independent-particle motion. LFE cause a strong positive dispersion of the intraband peak in ǫ 2 at small q (see red dashed line in Fig 3). Moreover, due to the strong depolarization, the plasmon follows the behavior of the electron-hole continuum (compare dashed and solid red lines in Fig. 3). In fact, in a localised system, where LFE are very large, the long-range part of the Coulomb interaction, which is the difference between v andv, becomes negligible andχ and χ becomes very similar. In turn, since χ gives the loss function −Imǫ  becomes completely damped. We now consider also 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 Fig. 5 (green line), the effect of interband transitions is a redshift of the continuum of electronhole transitions. As a consequence, also the zero of ǫ 1 and hence the plasmon frequency results redshifted by about 0.7 eV. Moreover, the positive plasmon dispersion, found taking into account only intraband transitions and LFE, is strongly reduced (see green line in Fig. 3) with an increase of the damping.
To better understand the interplay between interband transitions and 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: where W inter (ω) = ǫ −1 inter (ω)v is the Coulomb interaction that is dynamically screened by the interband transitions. Then the total polarizability χ can be written as: 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 Eq. (10) simplifies to: In this energy range, where ǫ inter is real, we thus obtain for the loss function the following expression [see Eq. (6)]: By retaining only the zero-order contribution of the off-diagonal terms of χ intra from Eq. (9), we finally arrive at the following expresion for the loss functions: From Eq. (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 Eq. (13). The second term in Eq. (13) is instead related to the LFE: it is responsible for summing the offdiagonal 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 Eq. (13) become more and more important, ehnancing 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 interband transitions and LFE induces two competing effects on the dispersion of the intraband plasmon that manifest differently in 2D and 3D structures.
In the bulk 3D TMDs the LFE are negligible and the first term of Eq. (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 Eq. (13) becomes the most prominent contribution. In this case, interband transitions reduce the positive dispersion of the intraband plasmon (compare red and green lines in Fig. 3).

Conclusions
In summary, we have calculated the charge-carrier plasmon dispersion in two prototypical single-layer metallic transition metal dichalcogenides (TMD): TaS 2 and NbSe 2 . As in the bulk, the nondispersive metallic band prevents using the homogeneous electron gas as a model to describe the plasmon excitations in these materials. Our analysis has shown that the plasmon dispersion of the 2D TMDs can be even more easily tuned than their 3D counterparts. In the bulk the plasmon dispersion is mainly determined by the joint density of the states 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 local-field effects and interband transitions, which gives rise to a small bandwidth corresponding to a localised plasmon. By acting also 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.