Intervalley plasmons in crystals

TL;DR

This paper develops a theory for shortwave intervalley plasmons in multivalley 2D materials, incorporating local-field effects, and analyzes their dispersion and interaction with external charges.

Contribution

It introduces a comprehensive theory of intervalley plasmons in monolayer transition-metal dichalcogenides, including their dispersion and Coulomb interactions, which was previously less explored.

Findings

Derived the energy dispersion of intervalley plasmons.

Formulated the interaction between plasmons and external test charges.

Provided a single-plasmon pole approximation for the Coulomb potential.

Abstract

Collective charge excitations in solids have been the subject of intense research ever since the pioneering works of Bohm and Pines in the 1950s. Most of these studies focused on long-wavelength plasmons that involve charge excitations with a small crystal-momentum transfer, $q \ll G$, where $G$ is the wavenumber of a reciprocal lattice vector. Less emphasis was given to collective charge excitations that lead to shortwave plasmons in multivalley electronic systems (i.e., when $q \sim G$). We present a theory of intervalley plasmons, taking into account local-field effects in the dynamical dielectric function. Focusing on monolayer transition-metal dichalcogenides where each of the valleys is further spin-split, we derive the energy dispersion of these plasmons and their interaction with external charges. Emphasis in this work is given to sum rules from which we derive the interaction between intervalley plasmons and a test charge, as well as a compact single-plasmon pole expression for the dynamical Coulomb potential.

Figures

Tables

Table 1. Table 1: Characteristic values of β 𝛽 \beta and α 0 subscript 𝛼 0 \alpha_{0} in electron- and hole-doped ML-TMDs. We have assumed that η c = 0.2 subscript 𝜂 𝑐 0.2 \eta_{c}=0.2 and η v = 0.47 subscript 𝜂 𝑣 0.47 \eta_{v}=0.47 (see text).

	${\beta }_c$	${\beta }_v$	${\alpha }_{0,c}$	${\alpha }_{0,v}$
WSe2	0.379	-0.333	1.35	0.52
WS2	0.333	-0.280	1.11	0.47
MoSe2	-0.138	-0.143	1.80	0.92
MoS2	-0.064	-0.115	1.46	0.76
MoTe2	-0.164	-0.195	2.14	1.11

Table 2. Table 2: Values of the bare effective mass at the edges of the top or bottom valleys in the conduction ( m c ​ t , 0 subscript 𝑚 𝑐 𝑡 0 m_{ct,0} & \& m c ​ b , 0 subscript 𝑚 𝑐 𝑏 0 m_{cb,0} ) and valence bands ( m v ​ t , 0 subscript 𝑚 𝑣 𝑡 0 m_{vt,0} & \& m v ​ b , 0 subscript 𝑚 𝑣 𝑏 0 m_{vb,0} ). These results are taken from Ref. [ Kormanyos_2DMater15, ]. The units are the free electron mass. Also provided is our estimated value of the polaron parameter ( δ P subscript 𝛿 𝑃 \delta_{P} ).

	$m_{ct,0}$	$m_{cb,0}$	$m_{vt,0}$	$m_{vb,0}$	${\delta }_P$
WSe2	0.29	0.4	0.36	0.54	0.17
WS2	0.27	0.36	0.36	0.5	0.1
MoSe2	0.58	0.5	0.6	0.7	0.25
MoS2	0.47	0.44	0.54	0.61	0.15
MoTe2	0.61	0.51	0.62	0.77	0.46

Equations

$$V(\mathbf{q},\omega )=\frac{V_{\mathbf{q}}}{|\overline{\stackrel{ˉ}{ϵ}}(\mathbf{q},\omega )|},$$

$$V_{\mathbf{q}}^{\text{2D}}=\frac{2\pi e^2}{A{ϵ}_d(\mathbf{q})q}\text{or}{V}_{\mathbf{q}}^{\text{3D}}=\frac{4\pi e^2}{V{ϵ}_d(\mathbf{q})q^2},$$

$${ϵ}_{\mathbf{G},{\mathbf{G}}^{\prime }}(\mathbf{Q},z)={\delta }_{\mathbf{G},{\mathbf{G}}^{\prime }}-V_{\mathbf{Q}+\mathbf{G}}\cdot \sum _{\mathbf{k},\nu }\frac{f({\epsilon }_{b,\mathbf{k}})-f({\epsilon }_{t,\mathbf{k}+\bar{\mathbf{q}}}+\Delta )}{z_{\nu }-(\Delta +{\epsilon }_{t,\mathbf{k}+\bar{\mathbf{q}}}-{\epsilon }_{b,\mathbf{k}})}\times ⟨\mathbf{k}+\mathbf{Q}|e^{i(\mathbf{Q}+{\mathbf{G}}^{\prime })\mathbf{r}}|\mathbf{k}⟩⟨\mathbf{k}|e^{-i(\mathbf{Q}+\mathbf{G})\mathbf{r}}|\mathbf{k}+\mathbf{Q}⟩.$$

$$\overline{\stackrel{ˉ}{ϵ}}(\bar{\mathbf{q}},z)=\mathcal{I}-V_{{\mathbf{K}}_0}\chi (\bar{\mathbf{q}},z)W{U}^T.$$

$$\chi (\bar{\mathbf{q}},z)\equiv \sum _{\mathbf{k},\nu }\frac{f({\epsilon }_{b,\mathbf{k}})-f({\epsilon }_{t,\mathbf{k}+\bar{\mathbf{q}}}+\Delta )}{z_{\nu }-(\Delta +{\epsilon }_{t,\mathbf{k}+\bar{\mathbf{q}}}-{\epsilon }_{b,\mathbf{k}})},$$

$$W_{\mathbf{G}}$$

$$U_{{\mathbf{G}}^{\prime }}$$

$$\left|\overline{\stackrel{ˉ}{ϵ}}(\bar{\mathbf{q}},z)\right|=1-\frac{V_{{\mathbf{K}}_0}}{\eta }\sum _{\mathbf{k},\nu }\frac{f({\epsilon }_{b,\mathbf{k}})-f({\epsilon }_{t,\mathbf{k}+\bar{\mathbf{q}}}+\Delta )}{z_{\nu }-(\Delta +{\epsilon }_{t,\mathbf{k}+\bar{\mathbf{q}}}-{\epsilon }_{b,\mathbf{k}})},$$

$$\frac{1}{\eta }=\sum _{\mathbf{G}}\frac{V_{{\mathbf{K}}_0+\mathbf{G}}}{V_{{\mathbf{K}}_0}}|\mathcal{F}({\mathbf{K}}_0+\mathbf{G}){|}^2.$$

$$\mathcal{F}({\mathbf{K}}_0+\mathbf{G})=⟨{\mathbf{K}}_f|e^{i({\mathbf{K}}_0+\mathbf{G})\mathbf{r}}|{\mathbf{K}}_i⟩.$$

$${\alpha }_0=\frac{1}{\eta K_0}\frac{m_b{e}^2}{{ϵ}_d(K_0){\hslash }^2}.$$

$$\frac{\beta }{{\alpha }_0}=\sum _{\nu ,j}{S}_{\nu ,j}\ln \left(\frac{\sqrt{a_{\nu ,j}^2-2c_j{b}_{\nu ,j}+c_j^2}+c_j-b_{\nu ,j}}{|a_{\nu ,j}|-b_{\nu ,j}}\right),$$

$$a_{\nu ,j}$$

$$b_{\nu ,j}$$

$$c_j$$

$$\hslash {\omega }_{\mathrm{s}}(\bar{q}=0)=\sqrt{{\Delta }^2+\frac{2(\gamma -c_0)\beta }{\gamma -1}\Delta {\epsilon }_F+\frac{(\gamma -c_0^2){\beta }^2}{\gamma -1}{\epsilon }_F^2},$$

$$\hslash {\omega }_{\mathrm{s}}(\bar{\mathbf{q}})\approx \Delta +\frac{\gamma \beta }{\gamma -1}{\epsilon }_{\mathrm{F}}+\left(\frac{\gamma -1}{\beta }+\gamma \right){\epsilon }_{t,\bar{\mathbf{q}}}.$$

$$\bar{q}\lesssim \frac{\beta }{\gamma -1}\frac{k_F}{\beta +1}.$$

$$\frac{V_{{\mathbf{K}}_0+\mathbf{G}}}{V_{{\mathbf{K}}_0}}=\frac{K_0}{|{\mathbf{K}}_0+\mathbf{G}|}\frac{{ϵ}_d(K_0)}{{ϵ}_d(|{\mathbf{K}}_0+\mathbf{G}|)},$$

$${ϵ}_d(|{\mathbf{K}}_0+\mathbf{G}|)\simeq 1+({ϵ}_d(K_0)-1){\left(\frac{K_0}{|{\mathbf{K}}_0+\mathbf{G}|}\right)}^P.$$

$${\Sigma }_i(\mathbf{k},z)=-k_{B}T\sum _{\mathbf{q},z^{\prime }}\frac{V_{\mathbf{q}}^{\text{2D}}}{|\overline{\stackrel{ˉ}{ϵ}}(\mathbf{q},z^{\prime }-z)|}{G}_j(\mathbf{k}-\mathbf{q},z^{\prime }).$$

$$G_j(\mathbf{k}-\mathbf{q},z^{\prime })=\frac{1}{z^{\prime }-{\epsilon }_{j,\mathbf{k}-{\mathbf{q}}_1}-\Delta \cdot {\delta }_{j,t}+\mu },$$

$$\frac{1}{\widetilde{ϵ}(\bar{\mathbf{q}},\omega )}=1+\frac{r(\bar{\mathbf{q}})}{{\omega }^2-{\omega }_{\bar{\mathbf{q}}}^2}.$$

$$1-V_{{\mathbf{K}}_0}\frac{\chi (\bar{\mathbf{q}},\omega \to \infty )}{\eta }=1-\frac{r(\bar{\mathbf{q}})}{{\omega }^2}$$

$$r(\bar{\mathbf{q}})=2{\int }_0^{\infty }d\omega \omega \mathrm{Im}\left\{ϵ(\mathbf{q},\omega )\right\}.$$

$$r(\bar{\mathbf{q}})$$

$$1-V_{{\mathbf{K}}_0}\frac{\chi (\bar{\mathbf{q}},\omega =0)}{\eta }=1+\frac{r(\bar{\mathbf{q}})}{{\omega }_{\bar{\mathbf{q}}}^2-r(\bar{\mathbf{q}})}.$$

$${\omega }_{\bar{\mathbf{q}}}^2=r(\bar{\mathbf{q}})\left[1+\frac{|\beta |}{2{\alpha }_0{\mathcal{G}}_{\bar{\mathbf{q}}}}\right],$$

$${\mathcal{G}}_{\bar{\mathbf{q}}}=\ln \frac{1+|\beta |\mathcal{R}(\Delta +{\epsilon }_{t,\bar{\mathbf{q}}},(1+\beta ){\epsilon }_{t,\bar{\mathbf{q}}},{\epsilon }_F)}{1+|\beta |\Theta ({\epsilon }_F-\Delta )\mathcal{R}(\Delta -{\epsilon }_{b,\bar{\mathbf{q}}},{\epsilon }_{b,\bar{\mathbf{q}}},c_0{\epsilon }_F)},$$

$$\mathcal{R}({\epsilon }_1,{\epsilon }_2,{\epsilon }_3)=\frac{\sqrt{({\epsilon }_1+\beta {\epsilon }_3{)}^2-4{\epsilon }_2{\epsilon }_3}-({\epsilon }_1-|\beta |{\epsilon }_3)}{(|\beta |+\beta ){\epsilon }_1-2{\epsilon }_2}.$$

$${\hslash }^2{\omega }_{\bar{\mathbf{q}}}^2={\Delta }^2+2(1-c_0){\alpha }_0\Delta {\epsilon }_F$$

$$+2(1+c_0)\left[{\alpha }_0{\epsilon }_F+\frac{{\Delta }^2}{(1-c_0)\Delta -(1+c_0){\epsilon }_{\bar{\mathbf{q}}}}\right]{\epsilon }_{\bar{\mathbf{q}}}.$$

Full text

I Introduction

Plasmon is the quantum of charge excitations in many-electron systems, manifested as an organized collective oscillation of the entire electron gas. From its outset, the theory of this phenomenon in crystals was developed for the long-range part of the Coulomb interaction,Pines_RMP56 ; Bohm_PR51 ; Pines_PR52 ; Bohm_PR53 ; Pines_PR53 ; Nozieres_PR58 and soon after, was used to explain the observed energy loss of electrons passing through metal foils.Ferrell_PR56 The theory was later extended to include surface plasmons, quasiparticles formed by the coupling between plasmons and electromagnetic fields at the surface of a metal-dielectric interface,Ritchie_PR57 finding applications in sensing,Homola_SABC99 cancer therapy,Lai_ACR08 lasing,Berini_NatPhys12 and plasmonic waveguides.Bozhevolnyi_Nat06 To date, the bulk of experimental and theoretical works remained focused on long-wavelength plasmons.

The recent intense research on two-dimensional (2D) multi-valley crystals such as graphene and monolayer transition-metal dichalcogenides (ML-TMDs) has drawn attention to collective shortwave charge excitations.Tudorovskiy_PRB10 ; Dery_PRB16 ; Groenewald_PRB16 ; VanTuan_PRX17 ; VanTuan_arXiv18 Plasmons in these materials are governed either by long-wavelength charge excitations within the same valley (intravalley plasmons) or shortwave excitations between the time-reversed valleys (intervalley plasmons). ML-TMDs in particular, are an ideal testbed for intervalley plasmons because of the valley spin splitting, as shown in Fig. 1 for the conduction-band edge of an electron-doped sample. The spin splitting gives rise to an energy gap in the dispersion of shortwave plasmons,Dery_PRB16 ; Groenewald_PRB16 allowing one to tell them apart from the gapless intravalley plasmons in 2D systems.Haug_Koch_Book Recently, we have demonstrated that the gapped energy dispersion of intervalley plasmons leads to unique features in the optical spectra of ML-TMDs.VanTuan_arXiv18 ; VanTuan_PRX17

In this work, we present a comprehensive analysis of shortwave plasmons. We introduce in Sec. II an efficient procedure to calculate intervalley plasmon modes when local-field effects are included. The procedure is general and can be applied in 2D or 3D systems. Continuing with 2D systems with emphasis on ML-TMDs, we show the dielectric loss function at low temperatures, followed by the derivation of concise expressions for the plasmon energy and damping-free propagation range at zero temperature. In Sec. III, we use sum rules to replace the cumbersome form of the dynamical Coulomb potential in the random-phase approximation (RPA) with a compact single-plasmon pole (SPP) expression that includes local-field effects. We then quantify the Coulomb exchange and correlation contributions to the self-energy of electrons (or holes). In Sec. IV, we derive the $f$-sum rule for intervalley plasmons in ML-TMDs through which we express the plasmon interaction with a test charge. The latter can be a remote electron that passes through the crystal. It can also be a core or valence-band electron excited to the Fermi surface, leading to shake-up of the surrounding electron system, similar to $X$-ray catastrophe in metals or Fermi-edge singularity in degenerate semiconductors.Mahan_PR67a ; Mahan_PR67 ; Nozieres_PR69 ; Schotte_PR69 ; Combescot_JdP71 ; Hawrylak_PRB91 ; Skolnick_PRL87 The appendices include technical details on the calculations of local-field effects and parameter choices.

II General Formalism

Plasmons are studied through the dynamically-screened Coulomb potential,

[TABLE]

where $\omega$ is the angular frequency and $\mathbf{q}$ is the crystal momentum (wavevector). $V_{\mathbf{q}}$ is the bare potential and $|\bar{\bar{\epsilon}}(\mathbf{q},\omega)|$ is the determinant of the dynamical dielectric function, which comes into a matrix form when local-field effects are considered.Adler_PR1962 ; Wiser_PR1963 The bare potential reads

[TABLE]

where $A$ ($V$) is the sample area (volume) in the case of a 2D (3D) system. $\epsilon_{d}(\bm{q})$ is the non-local dielectric function whose role is to capture the $q$-dependence of the effective dielectric constant due to material parameters of the ML and its surrounding.Latini_PRB15 ; Qiu_PRB16 It is not related and should not be confused with the static limit $\omega\rightarrow 0$ of $|\bar{\bar{\epsilon}}(\mathbf{q},\omega)|$. The role of the dynamical dielectric matrix is to capture the response of delocalized electrons (or holes) in the ML to a test charge, and in the limit of zero charge density $\bar{\bar{\epsilon}}(\mathbf{q},\omega)$ becomes the identity matrix.

We study the dynamical dielectric matrix of the two-valley system in Fig. 1. The valleys are centered around distinct points in the Brillouin zone (BZ), marked by $\bm{K}_{i}$ and $\bm{K}_{f}$. $\bm{K}_{0}$ is the wavevector that connects the valley centers. Each of the valleys is spin-split, where $\Delta$ is the spin splitting energy at the valley center. We assume parabolic energy dispersion for electronic states, $\varepsilon_{b,\bm{k}}=\hbar^{2}k^{2}/2m_{b}$ and $\varepsilon_{t,\bm{k}}=\hbar^{2}k^{2}/2m_{t}$, where $m_{b}$ and $m_{t}$ are the effective masses in the bottom and top valleys, respectively. Using the RPA, the elements of the dynamical dielectric matrix of the two-valley system readDery_PRB16

[TABLE]

$\mathbf{G}$ and $\mathbf{G}^{\prime}$ are reciprocal lattice vectors, and the wavevector $\bm{Q}=\mathbf{K}_{0}+\bar{\mathbf{q}}$ is restricted to the first BZ. $\nu=\{0,1\}$ are the two spin configurations that contribute to intervalley excitations (Fig. 1), and $z_{\nu}=(-1)^{\nu}\hbar\omega$. The state $|\mathbf{k}\rangle$ resides in the $\bm{K}_{i}$-point valley, where the valley center is the origin for $\mathbf{k}$. $|\mathbf{k}+\mathbf{Q}\rangle$ resides in the $\bm{K}_{f}$-point valley, whose center is the origin for $\mathbf{k}+\bar{\mathbf{q}}$. $f(\varepsilon_{b,\mathbf{k}})$ and $f(\varepsilon_{t,\mathbf{k}+\bar{\mathbf{q}}}+\Delta)$ are Fermi-Dirac distributions in the bottom and top valleys, respectively.

Shortwave plasmon modes are found from values of $z$ as a function of $\bar{\mathbf{q}}$ for which the determinant vanishes, $\left|\bar{\bar{\epsilon}}(\bar{\mathbf{q}},z=\hbar{\omega_{s}(\bar{\mathbf{q}}}))\right|=0$. The determinant calculation is greatly simplified if the matrix elements in Eq. (3) can be evaluated by replacing $|\mathbf{k}\rangle$ and $|\mathbf{k}+\mathbf{Q}\rangle$ with the valley-center states, $|\mathbf{K}_{i}\rangle$ and $|\mathbf{K}_{f}\rangle$. This condition is met when $\bar{q},k_{F}\ll K_{0}$, where $k_{F}$ is the Fermi wavenumber, and it is typically the case in semiconductors and bad metals. We can then write the dynamical dielectric matrix in a compact form,

[TABLE]

$\mathcal{I}$ is the identity matrix,

[TABLE]

is the density response function, while $W$ and $U$ are column vectors with elements

[TABLE]

Applying Sylvester’s determinant theorem,Harville1997 $\left|\mathcal{I}-\chi(\bar{\mathbf{q}},z)WU^{T}\right|=\left|\mathcal{I}-\chi(\bar{\mathbf{q}},z)U^{T}W\right|$, it is far simpler to calculate the determinant through the scalar $U^{T}W$ instead of the square matrix $WU^{T}$. We get

[TABLE]

where $\eta^{-1}=U^{T}W$ is a scalar that lumps together local-field effects

[TABLE]

The sum is over reciprocal lattice vectors ($\bm{G}$), and

[TABLE]

II.1 Applications

The equation set (7)-(9) is general and can be applied to multi-valley 2D or 3D materials. For example, finding low-energy shortwave plasmons in bulk elemental semiconductors such as silicon or germanium is straightforward, wherein $\Delta=0$ due to the space inversion symmetry of diamond-structure crystals. In electron-doped germanium, one can choose any pair of the four distinct valleys in the bottom of the conduction band because the zone-edge $L$-points (valley centers) are equivalent.Pezzoli_PRB13 ; Pengke_PRL13 ; Jacoboni_PR81 ; Yu_Cardona_Book ; Streitwolf_PSS70 In electron-doped silicon, on the other hand, there can be two types of low-energy intervalley excitations depending on which two valleys are selected from the six distinct valleys in the bottom of the conduction band. Borrowing the nomenclature of intervalley scattering types in silicon,Yu_Cardona_Book ; Streitwolf_PSS70 ; Li_PRL11 ; Song_PRL14 $g$-type intervalley plasmons arise if the chosen pair resides on the same crystal axis and $f$-type arises when they reside on different axes. The two types differ in the values of $K_{0}$ and $\eta$. The $g$-type excitation is similar to the one in graphene or ML-TMDs in the sense that the two valleys are related through time-reversal symmetry.Tudorovskiy_PRB10 ; Dery_PRB16 ; Groenewald_PRB16 As a result, the orbital compositions in the valley centers are similar.

Hereafter, we focus on 2D semiconductor systems. In this case, the determinant in Eq. (7) depends on the parameter

[TABLE]

The role of intervalley plasmons is negligible if $\alpha_{0}\ll 1$. We note that the materials below and above the two-valley system do not affect the value of $\epsilon_{d}(K_{0})$ because of the shortwave nature of $K_{0}$.

Figure 2 shows an example of the shortwave dielectric loss function at $T=10$ K and $\alpha_{0}=1.45$. The loss function is the imaginary part of $\left|(\bar{\bar{\epsilon}}(\bar{\bm{q}},\hbar\omega+i\delta))\right|^{-1}$, where we have assigned $\delta=k_{B}T$ to represent broadening. The effective mass parameters are $m_{b}=0.5m_{0}$ and $m_{t}=0.4m_{0}$, the spin-splitting energy is $\Delta=35$ meV, and the charge density is $n\simeq 4.2\times 10^{12}$ cm*-2* corresponding to a Fermi energy $\varepsilon_{F}=\pi\hbar^{2}n/m_{b}=20$ meV. The inset shows the calculation while neglecting the contribution from the off-resonance term [$\nu=1$ in Eq. (7)]. Comparing the inset and the main figure, the added effect of the off-resonance term is to somewhat shrink the amplitude (see color-bar scales), range of damping-free plasmon propagation (extension of the plasmon branch along the $x$-axis), and the plasmon energy (onset of the branch on the $y$-axis).

II.2 Zero-temperature ($T=0$)

The plasmon modes at $T=0$ are easier to calculate because the sum over $\bm{k}$ in Eq. (7) is handled analytically. The plasmon modes are found from the solution of

[TABLE]

where $\beta=m_{b}/m_{t}-1$ is the valley mass asymmetry, and the sum has four terms, $\nu=\{0,1\}$ and $j=\{0,1\}$. $S_{\nu,j}=(-1)^{j+1}\operatorname{{\mathrm{s}gn}}(a_{\nu,j})$ is a sign function, and

[TABLE]

$\Theta(\varepsilon_{F}-\Delta)$ is the step function and $\delta_{j,i}$ is the Kronecker delta. The contribution of terms with $j=0$ in Eq. (11) vanishes when $\varepsilon_{F}<\Delta$ because $c_{0}=0$ and the argument of the logarithm becomes one. The damping-free propagation range, $\bar{q}\leq q_{\text{max}}$, is defined by the values of $\bar{q}$ for which the solution of Eq. (11), $\hbar\omega=\hbar\omega_{\mathrm{s}}(\bar{q})$, is a real-value plasmon energy.

The solution of the transcendental Eq. (11) becomes analytical when $\bar{q}=0$, corresponding to intervalley transitions with wavenumber $K_{0}$. The energy of this plasmon mode reads

[TABLE]

where $\gamma=\exp(\beta/\alpha_{0})$. In mass symmetric valleys in which the electron masses in the top and bottom valleys are the same, $\beta=0$ and $\gamma=1$, we get that $\beta/(\gamma-1)\rightarrow\alpha_{0}$, and accordingly $\hbar\omega_{\mathrm{s}}(\bar{q}\!=\!0)\rightarrow\sqrt{\Delta^{2}+2(1\!-\!c_{0})\alpha_{0}\Delta\varepsilon_{F}}$. Note that $a_{\nu,j}$ and $b_{\nu,j}$ in Eq. (12) are ill-defined when $\beta=0$, so that a general mode ($\bar{q}\neq 0$) in mass symmetric systems can be studied by applying the limit $\beta\rightarrow 0^{+}$ or $\beta\rightarrow 0^{-}$.

II.3 $\Delta>\varepsilon_{\mathrm{F}}\gg k_{B}T$

The plasmon energy-dispersion relation, $\hbar\omega_{\mathrm{s}}(\bar{q})$, and the damping-free propagation range, $\bar{q}\leq q_{\text{max}}$, can be approximated by compact analytic expressions in the regime $\Delta>\varepsilon_{\mathrm{F}}\gg k_{B}T$. Here, the condition $\Delta,\varepsilon_{\mathrm{F}}\gg k_{B}T$ implies that we can still use the zero-temperature approximation, while $\Delta>\varepsilon_{\mathrm{F}}$ implies that only terms with $j=1$ contribute to the sum in Eq. (11). The solution becomes analytical if we further neglect the contribution from the term ($j=1$ and $\nu=1$) because of its relatively small contribution, as shown by Fig. 2. Keeping only the resonance term ($j=1$ and $\nu=0$), one finds after some algebra that the dispersion relation can be compactly expressed by

[TABLE]

The damping-free propagation range of intervalley plasmons is limited to

[TABLE]

In mass symmetric valleys, $\beta=0$ and $\gamma=1$, we get that $\bar{q}\lesssim\alpha_{0}k_{F}$. The numerical solution of Eq. (11) when one includes the off-resonance term, $j=1$ and $\nu=1$, is very close to the compact expressions in Eqs. (14)-(15) when $\Delta\gg\varepsilon_{\mathrm{F}}$. The validity of these expressions degrades as $\varepsilon_{\mathrm{F}}$ continues to grow because the effect of the spin-splitting energy is slowly washed-out.

II.4 ML-TMDs

Intervalley plasmons in ML-TMDs have been studied both through DFT calculations and analytically.Groenewald_PRB16 ; Dery_PRB16 In Ref. [Dery_PRB16, ], the approach was to neglect the mass asymmetry and local-field effects (i.e., assign $\eta=1$ in Eq. (7) and use $\beta=0$). However, the inclusion of local-field effects in these materials is not cosmetic but rather crucial because of the shortwave nature of intervalley plasmons. For example, the $C_{3}$ symmetry of the honeycomb crystal dictates that at least three reciprocal lattice vectors provide similar amplitudes for $\mathbf{K}_{0}+\mathbf{G}$. These include in the lowest order $|\mathbf{K}_{0}|=|\mathbf{K}_{0}+\bm{G}_{-}|=|\mathbf{K}_{0}-\bm{G}_{+}|=4\pi/3a$ where $a\simeq 3.2~{}\AA$ is the lattice constant, $\bm{K}_{0}=(2\pi/a)\left[0,\tfrac{2}{3}\right]$, and $\bm{G}_{\pm}=(2\pi/a)\left[\sqrt{\tfrac{1}{3}},\pm 1\right]$ are the basis vectors of the reciprocal lattice. Therefore, the approach taken in Ref. [Dery_PRB16, ], where umklapp processes were dispensed altogether (considering only the term $\bm{G}=\bm{G}^{\prime}=0$) evidently underestimates the amplitude of intervalley plasmons in the dynamical dielectric function and the damping-free propagation range.

The two-valley system in Fig. 1 and the formalism we have presented so far directly apply to electron-doped ML-TMDs. The spin-splitting energy in the conduction band is $\Delta=\Delta_{c}$, and the effective masses at the top and bottom conduction-band valleys are $m_{t}=m_{\text{ct}}$ and $m_{b}=m_{\text{cb}}$, respectively. Local-field effects can be approximated by employing the orbital $d_{z^{2}}$ of the transition-metal atom to represent the atomic compositions in the valley centers [Eq. (9)]. The formalism for hole-doped samples is similar but with three changes. Conduction-band parameters are replaced by valence-band ones ($c\rightarrow v$). The index of the top and bottom valleys is exchanged ($b\leftrightarrow t$) because electrons first populate the bottom valleys in the conduction band (when $\Delta>\varepsilon_{\mathrm{F}}$), while holes populate the top valleys in the valence band. Finally, we use the orbital $d_{(x\pm iy)^{2}}$ instead of $d_{z^{2}}$ when dealing with valence-band states.Zhu_PRB11

There are two main differences between intervalley plasmons in electron- and hole-doped ML-TMDs. The first one is the value of the energy gap in the dispersion relation, $\Delta$ in Eqs. (12) and (14). Its value is of the order of a few hundreds meV in the valence band, coming from the relatively strong spin-orbit interaction involved with the orbital $d_{(x\pm iy)^{2}}$ of transition-metal atoms.Zhu_PRB11 ; Xiao_PRL12 Accordingly, the condition $\Delta<\varepsilon_{\mathrm{F}}$ cannot be readily met in hole-doped ML-TMDs. On the other, $\Delta$ is about ten-fold smaller in the conduction band due to the vanishing spin-orbit coupling involved with $d_{z^{2}}$. Yet, the value of $\Delta$ is finite because the crystal field slightly distorts the atomic orbitals.Song_PRL13 In low-temperature ML-WSe2, for example, electron population of the top valley starts at charge densities around $6\times 10^{12}$ cm*-2*.Wang_NatNano2017 ; Wang_PRL18

The second important difference between intervalley plasmons in electron- and hole-doped ML-TMDs deals with local-field effects, expressed through the parameter $\eta$. A strong local-field effect amounts to $\eta\rightarrow 0$ and as a result to $\gamma\rightarrow 1$ and $\alpha_{0}\gtrsim 1$. An accurate but demanding calculation of $\eta$ from Eq. (8) requires first-principles simulations of the electronic states and non-local dielectric function. However, one can achieve a good estimate for $\eta$ by using two approximations. The first one is to calculate $\mathcal{F}({q})$ from Eq. (9) by assuming hydrogen-like wavefunctions for the $4d$ ($5d$) orbitals in Mo atoms (W atoms). Appendix A provides details on the straightforward calculation. This calculation is analytical and depends only on the parameter $r_{0}=4a_{0}/Z_{\text{Mo}}$ for Mo atoms and $r_{0}=5a_{0}/Z_{\text{W}}$ for W atoms, defined by the Bohr radius of hydrogen $a_{0}=0.529$ Å, and the effective nuclear charges, $Z_{\text{Mo}}$ and $Z_{\text{W}}$. The angular dependence of the orbitals follows spherical harmonics, $Y_{\ell,m}(\theta,\phi)$, with $\{\ell=2,m=0\}$ for electrons ($d_{z^{2}}$) and $\{\ell=2,m=\pm 2\}$ for holes ($d_{(x\pm iy)^{2}}$). Figure 3 shows plots of $\mathcal{F}({q})$ as a function of $qr_{0}$.

The second approximation we use to estimate $\eta$ is to express the potential ratio in Eq. (8) by

[TABLE]

where the non-local dielectric function follows

[TABLE]

The power-law parameter $P$ denotes how fast the non-local dielectric function decays to unity because of the vanishing induced potential when $G\rightarrow\infty$. The form in Eq. (17) is extracted from the results of DFT calculations in ML-MoS2,Latini_PRB15 ; Qiu_PRB16 where $\epsilon_{d}(K_{0}\simeq 1.3\,\AA^{-1})\sim 2.5$ and the power-law is about quadratic.

Table  1 lists values for the valley mass asymmetry ($\beta$), and amplitude measure ($\alpha_{0}$) in various ML-TMDs. We have assumed that $P=2$, $\epsilon_{d}(K_{0})=2.5$, $Z_{\text{Mo}}=14$ and $Z_{\text{W}}=23$. These parameters yield $\eta_{c}\simeq 0.2$ and $\eta_{v}\simeq 0.47$ for all ML-TMDs. Appendix A includes further details on the parameter choices for $Z_{\text{Mo}}$ and $Z_{\text{W}}$. The values of $\beta_{e}=m_{cb}/m_{ct}-1$ and $\beta_{h}=m_{vt}/m_{vb}-1$ are based on DFT calculations of the effective masses.Kormanyos_2DMater15 The values of $\alpha_{0}$, evaluated from Eq. (10), also include the polaron mass effect.VanTuan_PRB18 Namely, $m_{b(t)}=(1+\delta_{P})m_{b(t),0}$ where $m_{b(t),0}$ is the DFT result for the band-edge effective mass and $\delta_{P}$ is the polaron parameter, ranging from $\sim$0.1 in WS2 to as high as $\sim$0.5 in MoTe2 (see Appendix B for details). Table  1 shows that $\alpha_{0}$ is systematically larger in the electron-doped ML-TMDs, stemming from the differences in the orbital composition of electrons and holes that lead to $\eta_{c}^{-1}>\eta_{v}^{-1}$. This behavior can be understood from the faster decay of $\mathcal{F}({q})$ in the hole-doped case (Fig. 3), indicating that higher-order umklapp processes in Eq. (8), $|\mathbf{K}_{0}+\bm{G}|>K_{0}$, are more effective in enhancing the damping-free propagation range of intervalley plasmons in electron-doped conditions.

All in all, the effects of valley mass asymmetry and especially of local-field effects are evident in ML-TMDs. For example, the values of $\alpha_{0}$ in Table  1 decrease by an order of magnitude when employing $\eta=1$. The latter amounts to dispensing local-field effects altogether, which is the approach taken in Ref. [Dery_PRB16, ]. We will show in Secs. III and IV that the enhancement of $\alpha_{0}$ due to local-field effects has important implications on the amplitudes of the shortwave SPP Coulomb potential and the interaction of a test charge with intervalley plasmons.

III Single plasmon pole approximation and energy renormalization

The dynamical dielectric matrix is used in this section to evaluate the self-energies of electrons in the 2D Fermi sea. Employing a finite-temperature Green’s function formalism,MahanBook the lowest-order self-energy term follows

[TABLE]

$z$ and $z^{\prime}$ are odd Matsubara energies, and the index $i=\{b,t\}$ denotes bottom or top valley. The unperturbed Green’s function follows

[TABLE]

where $\mu$ is the chemical potential. In the long-wavelength regime ($q\rightarrow 0$), we assign $\bm{q}_{1}=\bm{q}$ and $j=i$ to denote intravalley excitations between electronic states below and above the Fermi surface of the same valley. On the other hand, in the shortwave regime where $q=|\bm{K}_{0}+\bar{\bm{q}}|$, we assign $\bm{q}_{1}=\bar{\bm{q}}$ and $j\neq i$ because the self-energy of an electron in the bottom valley ($i=b$) is affected by intervalley excitations with the top valley ($j=t$) and vice versa, as shown in Fig. 1.

The self-energy in the long-wavelength limit has been well studied and it corresponds to band-gap renormalization.Haug_Koch_Book ; VanTuan_arXiv18 ; Haug_SchmittRink_PqE84 ; SchmittRink_PRB86 ; Scharf_arXiv19 The dominant contribution to the self-energy comes from zeros of the determinant in Eq. (18), representing plasmon-induced intravalley virtual transitions. In the shortwave limit, on the other hand, the plasmon induces an intervalley virtual transition. In addition to a slight band-gap renormalization, we will show that the main resulting effect is distinct resonance features in the electron’s self-energy, attributed to the gapped energy-dispersion relation of intervalley plasmons.

III.1 Single-plasmon pole (SPP) approximation

The SPP approximation is a compact way to replace the relatively cumbersome RPA excitation spectrum by a single collective mode, $\omega_{\bm{q}}$. The dielectric function under the SPP approximation reads,

[TABLE]

The plasmon-pole frequency is denoted by $\omega_{\bar{\bm{q}}}$ and $r(\bar{\bm{q}})$ is the residue that represents the weight carried by the summation over $\bm{k}$ in the density response function. The residue can be found from the asymptotic behavior of the RPA dielectric function at high-frequencies ($\mathrm{Re}\left\{\epsilon(\bm{q},\omega\rightarrow\infty)\right\}$),

[TABLE]

Alternatively, employing the Kramers-Kronig relation, the residue can also be extracted from the conductivity sum-rule,

[TABLE]

Using Eq. (5) for the density response function, we get that

[TABLE]

Local field effects are lumped together in the parameter $\alpha_{0}\propto\eta^{-1}$. Recalling Eq. (8), one can also write that $r=\sum_{\bf G}r_{\bf G}$. Namely, the residue includes the contribution of different umklapp processes to the screening of the macroscopic field.

The single collective mode, $\omega_{\bar{\bm{q}}}$, is found from comparing the static limits of the RPA and SPP dielectric functions,

[TABLE]

A straightforward calculation yields

[TABLE]

where

[TABLE]

and

[TABLE]

In mass symmetric or nearly symmetric systems ($\beta\rightarrow 0$ and $\varepsilon_{\bar{\bm{q}}}=\varepsilon_{b,\bar{\bm{q}}}=\varepsilon_{t,\bar{\bm{q}}}$), and when $\Delta\gg\varepsilon_{F}\varepsilon_{\bar{\bm{q}}}$, we get that

[TABLE]

Similar to the long-wavelength case,Overhauser_PRB71 the values of $\omega_{\bar{\bm{q}}}$ in the SPP Coulomb potential approaches that of the plasmon mode when $\beta=0$ and $\varepsilon_{\bar{\bm{q}}}=0$, as can be seen by comparing Eqs. (13) and (28). Their values differ when $\bar{q}\neq 0$ because the single collective mode, $\omega_{\bar{\bm{q}}}$, is derived from the static limit of the dielectric function.

III.2 Self-energy

Substituting Eq. (20) in (18), we divide the self-energy into exchange- and correlation-related contributions,

[TABLE]

The exchange contribution comes from the bare potential

[TABLE]

and the correlation contribution comes from the dynamical part of the potential,

[TABLE]

We have used the fact that $V_{|\bm{K}_{0}+\bm{k}-\bar{\bm{q}}|}\simeq V_{\bm{K}_{0}}$ to take out the shortwave potential from the sum. Replacing the sum over Matsubara frequencies ($z^{\prime}$) with contour integration,

[TABLE]

the exchange part at low temperatures reads

[TABLE]

In ML-TMDs, the result is about $1$ meV redshift of the top valley per electron density of $n\sim 10^{12}$ cm*-2* in the bottom valley.

To evaluate the correlation part, we make use of the spectral representation of the SPP-based dielectric function,

[TABLE]

Using this definition in Eq. (31), replacing the sum over Matsubara frequencies with contour integration [Eq. (32)], and using Dirac identity,

[TABLE]

with $\mathcal{P}$ denoting the Cauchy principal value, we get

[TABLE]

The sum is limited to the damping-free propagation range ($\bar{q}\leq q_{\text{max}}$), and $g(x)$ denotes the Bose-Einstein distribution. As before, we either have $\{i=b,\,j=t\}$ or $\{i=t,\,j=b\}$. Figure 4 shows the correlation contribution to the self-energy at the band edge ($k=0$) for a system with the following parameters: $T=10$ K, $\varepsilon_{F}=20$ meV, $\Delta=30$ meV, $m_{b}=0.47m_{0}$, $m_{t}=0.34m_{0}$, $\alpha_{0}=1.35$ and $\eta=0.2$. In addition, we have used the thermal energy for broadening, $z=E+ik_{B}T$. The self-energy of the top valley (left panel) includes resonance features below the continuum. The singular region in the left panel of Fig. 4 lies at the interval $[-(\Delta-\varepsilon_{b,q_{\text{max}}}+\hbar\omega_{q_{\text{max}}}),-(\Delta+\hbar\omega_{\bar{q}=0})]$, where $q_{\text{max}}$ is the largest possible value for damping-free plasmon propagation. The singularity arises from the Fermi-distribution term in the first expression of the square brackets in Eq. (36). The second term in square brackets does not lead to a resonance feature because $q_{\text{max}}<k_{F}$ in this case, and hence $f(\varepsilon_{b,\bar{q}})+g(-\hbar\omega_{\bar{\bm{q}}})\simeq 0$ for the entire integration range. Figure 4 shows that the magnitude of the correlation resonance below the continuum reaches a few tens meV. This value is much larger than the one calculated in Ref. [Dery_PRB16, ], where local-field effects were neglected. Here, these effects are included through the residue $r(\bm{q})$, and the enhanced integration region ($q_{\text{max}}$).

III.3 Renormalization of $\Delta$

We have seen that the finite density leads to different energy shifts of the top and bottom valleys. As a result, the spin-splitting energy has charge-density-dependent contributions from exchange and correlation in addition to one from spin-orbit coupling,

[TABLE]

It is important to recognize that the exchange and correlation contributions can change the value of $\Delta$ only when $\Delta_{0}\neq 0$. In systems where the spin-orbit coupling does not lead to spin-split valleys (e.g., $\Delta_{0}=0$ in crystals with space inversion symmetry), the density in the spin-up and spin-down valleys is similar, and hence their energy shift is similar. That is, $\Delta_{x}(n)=\Delta_{c}(n)=0$ in systems where $\Delta_{0}=0$.

In systems with a finite $\Delta_{0}$, the spin-splitting energy has contributions from both intervalley and intravalley excitations,

[TABLE]

The spin splitting energy is evaluated from the valley edge states ($k=0$), denoted by the zero arguments of the self-energies. The first (second) term in square brackets on the right-hand sides of Eqs. (38)-(39) is the contribution from intervalley (intravalley) excitations. For the intravalley case, we have split the self-energy to contributions from screened-exchange and Coulomb-hole parts ($\Sigma_{sx}$ and $\Sigma_{ch}$).Haug_SchmittRink_PqE84 ; Scharf_arXiv19 ; VanTuan_arXiv18 The Coulomb-hole refers to the lack of charge next to a charged particle due to Pauli exclusion principle. The screened-exchange self-energy is calculated by using the SPP potential in the long-wavelength limit instead of the bare non-local potential.Scharf_arXiv19

The difference in the energy shifts of the top and bottom valleys mostly comes from the exchange contributions,

[TABLE]

The term associated with the factor $1/2$ stems from the intravalley screened-exchange interaction,Scharf_arXiv19 $\Sigma_{t,sx}(0)-\Sigma_{b,sx}(0)$, and the term $\eta\alpha_{0}$ stems from the intervalley exchange contribution, $\Sigma_{t,x}-\Sigma_{b,x}$ [see Eq. (33)]. The correlation term $\Delta_{c}(n)$, on the other hand, is very small for the following reasons. While the Coulomb-hole self-energies can be very large, they are similar for the top and bottom valleys regardless of the difference in their population. Furthermore, whether electrons or holes are present, both the conduction- and valence-bands shift by a similar magnitude with the only difference that the former (latter) shifts down (up). As a result, the band-gap energy shrinks when the charge-density increases.Scharf_arXiv19 The difference in the energy shifts of the top and bottom valleys due to the intervalley correlation terms, $\Sigma_{t,c}(0,E-\Delta-\mu)$ and $\Sigma_{b,c}(0,E-\mu)$, is also very small when $E\rightarrow 0$, as can be seen from Fig. 4. The main effect in this case is the resonance features that lie well below the continuum energy edge of the valleys.

IV The interaction of a test-charge with intervalley plasmon

In the previous section we have evaluated the self-energy correction of electrons in the Fermi sea due to intervalley plasmons. The plasmons are generated by the Fermi-sea electrons, and as such the electrons and plasmons do not damp the kinetic motion of one another (the self-energy renormalization comes from plasmon-induced virtual transitions). The scenario changes when a particle, external to the Fermi sea, shows up. Such a particle can be a remote electron that passes through the crystal. It can also be an exciton (photoexcited electron-hole pair) or a core electron that is excited to the Fermi surface. In this case, the scattering between the test-charge and the plasmon can induce a real transition with distinct initial and final states.

In this section, we derive the interaction Hamiltonian between a test-charge and intervalley plasmons. As before we will focus on electron-doped samples, and the case of hole-doped samples is similar but with the three changes mentioned in Sec. II.4. The starting point is Poisson’s Equation, from which we can write the interaction between charge excitations in the crystal and a test-charge,Overhauser_PRB71

[TABLE]

where $V_{\bm{q}}$ is the potential of the test-charge and $\rho_{\bm{q}}$ is the Fourier component of the charge density operator. The sum is not bound, and $\bm{q}$ can also take values outside the first BZ because the test-charge can be everywhere in the crystal and not only in lattice sites. Our goal is to express plasmons in terms of the charge density operator. We start by using second quantization and focusing on charge excitations due to spin-conserving intraband transitions (i.e., within the conduction band in electron-doped samples or within the valence band in hole-doped ones). We can then write that

[TABLE]

where $a^{\dagger}_{s,\bm{k}}$ and $a_{s,\bm{k}}$ are, respectively, the creation and annihilation operators of the charged particle in a state defined by spin and crystal momentum quantum numbers ($s$, $\bm{k}$). The summation over $\bm{k}$ is restricted to states in the first BZ, and so $\bm{G}_{\bm{k},\bm{q}}$ is the reciprocal lattice vector needed to bring back $\bm{k}+\bm{q}$ to the first BZ. Often times, one is interested in long-wavelength charge excitations, $\bm{q}\rightarrow 0$ and $\bm{G}_{\bm{k},\bm{q}}=0$, where the charge-density operator reduces to its familiar form $\rho_{\bm{q}\rightarrow 0}=\sum_{\bm{k},s}a^{\dagger}_{s,\bm{k}+\bm{q}}a_{s,\bm{k}}$. Here, on the other hand, we focus on spin-conserving shortwave excitations from low-energy states in the $\bm{K}_{i}$ valley to the $\bm{K}_{f}$ valley, as shown in Fig. 1. Accordingly, we use the notations,

[TABLE]

and rewrite the sum in Eq. (42) as

[TABLE]

where the sum runs over $\bm{q}_{G}$. Here $\bm{G}$ can be any reciprocal lattice vector, and we have that $\bar{q}\ll K_{0}$ due to the relatively small range of free-plasmon propagation. Equation (43) is then rewritten by adding the valley quantum number $\tau$ for states near $\bm{K}_{f}$ and $-\tau$ for states near $\bm{K}_{i}$,

[TABLE]

where $\bm{k}$ and $\bm{k}+\bar{\bm{q}}$ are measured from the valley centers. The approximation made on the right-hand side, where $\mathcal{F}(\mathbf{q})=\langle\mathbf{K}_{f}|e^{i\mathbf{q}\mathbf{r}}|\mathbf{K}_{i}\rangle$, is similar to the one we made when the dynamical dielectric matrix in Eq. (3) was replaced with (4). That is, the value of the matrix element calculated with the orbital composition of the valley-center states does not change appreciably for off-center states as long as $\bar{q},k\ll K_{0}$.

We now can repeat the approach taken by Nozières and Pines and later by Overhauser to find the interaction of plasmons with a test-charge.Nozieres_PR58 ; Overhauser_PRB71 In our case, the test-charge is an external perturbation to the two-valley electron system from which intervalley plasmons emerge. The derivation of the interaction relies on two ways from which one can calculate the double-commutator matrix element $\langle 0|[[H,\rho_{\bm{q}}],\rho^{\dagger}_{\bm{q}}]|0\rangle$. $H$ is the Hamiltonian of the unperturbed system in Fig. 1,

[TABLE]

The first way to calculate the double-commutator matrix element is a straightforward approach by using Eqs. (46) and (47). One gets after some algebra that

[TABLE]

where $A$ is the sample area and $r(\bar{\bm{q}})$ is the residue in Eq. (23). The second approach to calculate the double-commutator matrix element is by inserting a complete set of projection operators, $\mathcal{I}=\sum_{j}|j\rangle\langle j|$. Comparing results from both approaches, one gets

[TABLE]

where $\hbar\omega_{j0}$ is the energy difference between the excited and ground states ($|j\rangle$ and $|0\rangle$). Equation (49) is the $f$-sum rule for intervalley plasmons in the two-valley model (Fig. 1). $f$-sum rules like Eq. (49) take their name from an equivalent sum rule for dipole oscillator strengths in atomic physics, the celebrated Thomas–Reiche–Kuhn sum rule. They are a consequence of the fundamental laws of quantum mechanics and can be employed in various physical systems, like here, where we use Eq. (49) to calculate the interaction between a test charge and intervalley plasmons.

The next step in the derivation is the evaluation of the matrix element on the left-hand side of Eq. (49). We can do it by writing the charge-density operator in terms of plasmon creation and annihilation operators, $b^{\dagger}$ and $b$,

[TABLE]

Substituting the expression for $\rho_{\bm{q}}$ from Eq. (50) into the left-hand side of Eq. (49) and assuming a single collective excitation for a given $\bar{\mathbf{q}}$, one gets

[TABLE]

where $\omega_{\bar{\bm{q}}}$ is the single collective mode that we have found when deriving the Coulomb potential under the SPP approximation. The use of the plasmon mode from the numerical solution of Eq. (11) instead of the single collective mode from Eq. (25) is not suitable because the test-charge polarizes the crystal, and this effect is accounted for by the single collective mode that was derived under the static screening limit of the dielectric function. A similar scenario arises in the long wavelength limit.Overhauser_PRB71 ; Kato1999 ; KrstajicPRB2012 ; KrstajicPRB2013 Finally, using Eqs. (50)-(51) to rewrite Eq. (45), the interaction between plasmons and a test charge reads

[TABLE]

V Conclusions

Key aspects of intervalley plasmons in crystals were analyzed in this work. Using a two-band valley model, we have first studied the dynamical dielectric function matrix under the random-phase approximation. Unlike the case of long-wavelength (intravalley) plasmons, local-field effects become important due to the shortwave nature of intervalley plasmons, and as a result, umklapp processes contribute to intervalley charge excitations. We have introduced an effective method to incorporate local-field effects in the dynamical dielectric function, allowing one to extract the wavevector-energy dispersion relation of intervalley plasmons in two and three dimensional systems from the roots of a compact equation instead of a matrix determinant.

Focusing on two-dimensional problems, we have introduced the parameter $\alpha_{0}$, which is a measure for the magnitude of intervalley plasmons in a given material and for the extent of their damping-free propagation range. we have studied the dielectric loss function and derived the dispersion relation of intervalley plasmons with emphasis on monolayer transition-metal dichalcogenides in which the effective masses are different in the top and bottom spin-split valleys. Electrons in these materials generate intervalley collective excitations more effectively than holes because of the orbital composition of electronic states in the conduction and valence bands.

Finally, we have replaced the excitation spectrum of the dynamical Coulomb potential under the random-phase approximation by employing single-plasmon pole approximation. From the found single collective mode, we were able to analytically evaluate the self-energy of electrons (or holes) in the Fermi sea. Importantly, the single collective mode allows us to evaluate the $f$-sum rule from which we have derived the interaction between a test charge particle and intervalley plasmons. Such a particle can be a remote electron that passes through the crystal. It can also be an exciton (photoexcited electron-hole pair) or a core electron that is excited to the Fermi surface and shakes up the plasma during photon absorption.Mahan_PR67a ; Skolnick_PRL87 ; Hawrylak_PRB91 ; VanTuan_arXiv18 ; VanTuan_PRX17

To the best of our knowledge, a direct detection of intervalley plasmons in monolayer transition-metal dichalcogenides has not been demonstrated yet (i.e., not through the exciton optical transitions). Reflection electron energy loss spectroscopy, resonant Raman or THz spectroscopies are possible experiments to detect these plasmons. In Raman or THz spectroscopies one should recall that a photon can only couple to two counter-propagating shortwave plasmons in order to conserve momentum (the photon wavenumber is far smaller than the wavenumber that connects the valley centers). Electrostatic doping can be used to tell apart the signature of intervalley plasmons from that of optical phonons in the far-infrared spectrum. The gate voltage tunes the charge density and hence the plasmon energy and its amplitude (wider damping-free propagation range). Optical phonons, on the other hand, hardly change their energies.

Acknowledgements.

This work was mostly supported by the Department of Energy under Contract No. DE-SC0014349. The work at the University at Buffalo was supported by the Department of Energy, Basic Energy Sciences under Grant No. DESC0004890. The work in Würzburg was supported by the German Science Foundation (DFG) via Grant No. SFB 1170 “ToCoTronics” and by the ENB Graduate School on Topological Insulators.

Appendix A A simple calculation of $\eta_{c}$ and $\eta_{v}$

The value of $\eta_{c}$ for electron-doped samples and $\eta_{v}$ for hole-doped ones, is calculated from $j=\{c,v\}$

[TABLE]

$\bm{K}_{0}$ and the two-dimensional reciprocal lattice vectors in the sum, $\bm{G}$, are defined by

[TABLE]

where $a$ is the lattice constant, $m_{1}$ and $m_{2}$ take integer values, and $\bm{G}_{\pm}$ are the basis vectors of the reciprocal lattice.

The matrix elements in Eq. (53) are evaluated by considering a simple tight-binding model where the overlap between atomic orbitals of different lattice sites is neglected. Given that the conduction-band (valence-band) states near the $K$ and $K^{\prime}$ points are governed by the $d_{z^{2}}$ ($d_{(x\pm iy)^{2}}$) orbital of the transition-metal atom, we can write that

[TABLE]

where $\mathbf{q}$ is a two-dimensional wavevector ($q_{z}=0$), $R_{n}(r)$ is the radial part of the orbital, and $Y_{j}(\theta,\phi)=Y_{\ell,m}(\theta,\phi)$ is the spherical harmonics function where electrons (holes) are modeled by $\ell=2$ and $m=0$ ($\ell=2$ and $m=\pm 2$),

[TABLE]

Assuming hydrogen-like wavefunctions for the radial part, $R_{n}(r)$, where $n=4$ mimics the 4$d$ orbitals in MoSe2 and $n=5$ for the 5$d$ orbitals in WSe2, we get

[TABLE]

where $y=r/r_{0}$ and $r_{0}=na_{0}/Z_{\text{eff}}$ is an effective radius, defined by the Bohr radius in hydrogen $a_{0}=0.529$ Å, the energy level ($n=4$ for MoSe2 and $n=5$ for WSe2), and the effective nuclear charge seen by the $d$-orbital electrons, $Z_{\text{eff}}$. Substituting Eqs. (56)-(57) in (55), we get that

[TABLE]

where $x=qr_{0}$. Figure 3 in the main text shows plots of these expressions.

Using Eqs. (16), (17), (54) and (58), one can estimate the values of $\eta_{c}$ and $\eta_{v}$ from Eq. (53). Figure 5 shows the results for $\eta_{c}^{-1}$ (solid lines) and $\eta_{v}^{-1}$ (dashed lines) as a function of $Z_{\text{eff}}$. We have used $\epsilon_{d}(K_{0})=2.5$ and $P=2$ in Eq. (17).Latini_PRB15 ; Qiu_PRB16 In addition, the results are shown for both the 2D and 3D Coulomb potential forms. The 3D potential yields a smaller local-field effect due faster decay of umklapp processes ($q^{-2}$ in 3D vs $q^{-1}$ in 2D).

The most important parameter to evaluate the local-field effect in ML-TMDs is the effective nuclear charge $Z_{\text{eff}}$ from which we determine $r_{0}$. Choosing its value is somewhat subtle for the following reason. While the Slater’s rule for $5d$-orbitals in tungsten atoms yields $Z_{\text{eff}}\simeq 17$, it underestimates the lanthanide-contraction effect (poor screening of $4f$ electrons) and the strong spin-orbit coupling in tungsten. The lanthanide-contraction results in a smaller empirical radii of the outer electrons (i.e., larger value of $Z_{\text{eff}}$ is required to reproduce empirical values). In the main text, we have used $Z_{\text{eff}}=23$ because it yields very good agreement between theory and absorption-type experiments.VanTuan_arXiv18 The value we chose for the $4d$ orbitals in Mo, $Z_{\text{eff}}=14$, is $\sim$20-25% larger than the one estimated by employing the simple Slater’s rule ($Z_{\text{eff}}=11.392$). We use this value because it yields the same $\eta$ parameters in MoSe2 and WSe2, thereby reducing the number of parameters we use in the main text ($\eta_{c}\sim 0.2$ and $\eta_{v}\sim 0.47$ are assumed for all ML-TMDs; see Fig. 5).

Appendix B Effective masses and polaron effects in ML-TMDs

The electron-phonon Fröhlich interaction is relatively strong in ML-TMDs because of their polar nature.VanTuan_PRB18 ; Sohier_PRB16 One important consequence of the Fröhlich interaction is the polaron effect, manifested as a mass increase of charge particles,

[TABLE]

$m_{b(t),0}$ is the bare effective mass at the edge of the bottom (top) valley. $\delta_{P}$ is the polaron parameter, and we use this notation instead of the conventional $\alpha$-polaron parameter to prevent confusion with the $\alpha_{0}$-parameter of intervalley plasmons that we focus on in this work. Recently, we have used the polaron mass effect to fit the trion binding energies with empirical results, getting that $\delta_{P}=0.17$ in ML-WSe2 and $\delta_{P}=0.25$ in ML-MoSe2.VanTuan_PRB18 These values are commensurate with the Born effective charges of the metal atoms, $\delta_{P}\propto Z_{M}$, and the larger value in ML-MoSe2 compared with ML-WSe2 is inline with calculated DFT results of $Z_{M}$.Sohier_PRB16 Using this proportionality and the Born effective charges from Ref. [Sohier_PRB16, ], we can also estimate the polaron parameters of other ML-TMDs (Table 2).

Table 1 in the main text shows the estimated values for the mass asymmetry between the top and bottom valleys in the conduction band ($\beta_{c}=m_{cb}/m_{ct}-1$) and valence band ($\beta_{v}=m_{vt}/m_{vb}-1$). The Fröhlich interaction depends on the charge of the electron or hole, and therefore, we assume the polaron effect to be the same in the conduction and valence bands as well as for states in the top or bottom valleys. Accordingly, $\beta_{c}$ and $\beta_{v}$ are independent of the polaron effect and can be calculated from the bare effective masses. The values of the valley mass asymmetry in Table 1 of the main text are evaluated from DFT-based calculations of the bare effective masses at the valley edges.Kormanyos_2DMater15 Table 2 summarizes the bare mass values in ML-TMDs.