Enhanced electron-phonon interaction in multi-valley materials

TL;DR

This paper reveals a new mechanism that significantly enhances electron-phonon interactions in multi-valley materials, with implications for superconductivity, supported by combined theoretical and experimental evidence.

Contribution

It identifies a phonon-induced out-of-phase energy shift mechanism that strengthens electron-phonon coupling in multi-valley materials, a novel insight not previously understood.

Findings

Enhanced electron-phonon coupling observed in transition-metal dichalcogenides

Mechanism reduces electrostatic screening via inter-valley charge transfer

Potential relevance to superconductivity phenomena

Abstract

Through a combined theoretical and experimental effort, we uncover a yet unidentified mechanism that strengthens considerably electron-phonon coupling in materials where electron accumulation leads to population of multiple valleys. Taking atomically-thin transition-metal dichalcogenides as prototypical examples, we establish that the mechanism results from a phonon-induced out-of-phase energy shift of the different valleys, which causes inter-valley charge transfer and reduces the effectiveness of electrostatic screening, thus enhancing electron-phonon interactions. The effect is physically robust, it can play a role in many materials and phenomena, as we illustrate by discussing experimental evidence for its relevance in the occurrence of superconductivity. (short abstract due to size limitations - full abstract in the manuscript)

Figures

Tables

Table 1. Table 1: Energy separation Δ ​ E Γ ​ K Δ subscript 𝐸 Γ K \Delta E_{\rm\Gamma K} between the top of the valence band at K and Γ Γ \Gamma in the monolayer TMDs investigated in this work. Experimental results are taken from angle-resolved photo-emission spectroscopy (ARPES) measurements in the literature, while DFT data refer to calculations done in the present work.

	$\Delta E_{\Gamma K}$ (meV)
	MoS2	WS2	WSe2
Exp.	15059, 310 60; 54	51061, 24062,30063	56064, 50065, 89052
DFT (this work)	60	230	490

Equations

$$\Delta {\omega }_{\nu }\approx -N({\epsilon }_F)⟨g_{\nu }^2{⟩}_{FS}.$$

$$\Delta {\omega }_{\nu }^{(\tau )}\approx \frac{\Delta {\omega }_{\nu }}{1+({\omega }_{\nu }\tau {)}^2}$$

$$\Delta \omega \approx \frac{\Delta {\mathcal{D}}_{\Gamma }}{2M{\omega }_{\mathrm{neutral}}}.$$

$$\Delta {\mathcal{D}}_{\Gamma }$$

$${\mathcal{F}}_{\Gamma }$$

$$\displaystyle+\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k},n}|\langle\mathbf{k}n|\Delta V_{\Gamma,\nu}|\mathbf{k}n\rangle|^{2}\frac{\partial f}{\partial\varepsilon}\Big{|}_{\varepsilon_{\mathbf{k},n}}$$

$$\Delta \omega$$

$$\displaystyle=\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k},n}|g^{(n)}_{\nu,\mathbf{k},\mathbf{k}}|^{2}\frac{\partial f}{\partial\varepsilon}\Big{|}_{\varepsilon_{\mathbf{k},n}},$$

$$|g_{\nu ,\mathbf{k},\mathbf{k}}^{(n)}{|}^2=\frac{|⟨\mathbf{k}n|\Delta V_{\Gamma ,\nu }|\mathbf{k}n⟩{|}^2}{2M{\omega }_{\mathrm{neutral}}}\text{leavevmode}.$$

$$\Delta \omega$$

$$=-N({\epsilon }_F)⟨g_{\nu }^2{⟩}_{FS},$$

$$N({\epsilon }_F)=\frac{1}{N_{\mathbf{k}}}\sum _{\mathbf{k},n}\delta ({\epsilon }_{\mathbf{k},n}-{\epsilon }_F)\text{leavevmode},$$

$$⟨g_{\nu }^2{⟩}_{FS}=\frac{{\sum }_{\mathbf{k},n}|g_{\nu ,\mathbf{k},\mathbf{k}}^{(n)}{|}^2\delta ({\epsilon }_{\mathbf{k},n}-{\epsilon }_F)}{{\sum }_{\mathbf{k},n}\delta ({\epsilon }_{\mathbf{k},n}-{\epsilon }_F)}\text{leavevmode}.$$

Full text

I Introduction

The electronic and elastic properties of solids are determined by the interaction between charge carriers and the elementary excitations associated with the vibration of the crystalline lattice, namely the phonons. In view of its fundamental relevance, the nature of the interaction between electrons and phonons has been investigated in great depth in many different contexts Ziman ; Grimvall ; Mahan . In the simplest possible terms, electron-phonon interaction can be understood as originating from the electrostatic potential generated by the lattice distortion associated to the ionic displacement in the presence of a phonon. In other words, exciting a phonon in a crystal displaces the charged ions forming the lattice from their equilibrium positions, thereby generating a local charge imbalance and a corresponding electrostatic potential (known as deformation potential Bardeen1950 ; Schockley1950 ), which in turns directly acts on all mobile electrons.

Depending on the specific type of material considered and on the level of detail needed for a precise microscopic understanding, the situation can be more complex. For instance, in some systems (among which, graphene) the coupling between the lattice deformation associated to certain phonon modes and charge carriers is more appropriately described in terms of a vector potential, rather than an electrostatic one Pietronero1980 ; Sohier2014 . In these cases, charge carriers interact with phonons more similarly to the way they interact with a magnetic field Vozmediano2010 . Also, the influence of electron-phonon interaction on physical phenomena depends strongly on whether the energy and momentum relaxation of electrons is fast or slow relative to the phonon frequency, corresponding to the so-called adiabatic or anti-adiabatic limit of electron-phonon coupling Engelsberg1963 ; Maksimov1996 .

Irrespective of the specific situation, what is key to understand electron-phonon interaction (and its strength) in solids is electrostatic screening Mahan . Indeed, the deformation potential generated by the excitation of a long wavelength phonon tends to be screened by the spatial redistribution of the electrons present in the system, which self-consistently reduce the strength of the electron-phonon coupling Bogulsawski1977 ; Khan1984 ; Kartheuser1986 . Screening is expected to be much more effective in the adiabatic regime –since then electrons are able to equilibrate sufficiently fast in response to the ionic motion– and a much weaker electron-phonon coupling is accordingly expected in this regime as compared to the antiadiabatic one. Conversely, if the nature of the electron-phonon coupling is properly described in terms of a vector potential, as we mentioned it is the case for some of the phonons in graphene ferrari_raman_2007 ; ferrari_raman_2013 , screening cannot influence electron-phonon coupling Sohier2014 ; Pisana2007 (just simply because no spatial distribution of charge carriers can screen a magnetic field). In this case, electron-phonon coupling is expected not to be significantly affected by the presence of charge carriers.

In semiconducting materials hosting only a small density of mobile charges, electrostatic screening is poor and a strong electron-phonon interaction is expected. Addition of charge carriers (e.g., by doping) can drastically improve screening and cause the strength of electron-phonon interaction to decrease, as observed in a variety of semiconductors (e.g., SrTiO3 kirzhnits_description_1973 ; van_mechelen_electron-phonon_2008 ; Wang2016 and other transition metal oxides Verdi2017 ). For those cases in which phonons are not effectively screened by the presence of mobile charges (see example above), it can also happen that addition of charge carriers leaves the strength of electron-phonon interaction unaffected Sohier2014 . In virtually no case, however, the strength of electron-phonon interaction in a semiconductor is expected to increase significantly upon adding charge carriers.

In contrast with this established understanding, here we uncover a yet unidentified mechanism that causes a very significant strengthening –and not a weakening– of electron-phonon coupling in atomically-thin semiconductors upon increasing electron density. Our work relies on a joint experimental and theoretical investigation of Raman spectroscopy performed on mono and bilayers of different semiconducting transition metal dichalcogenides (TMDs; MoS2, WS2, and WSe2) as a function of density of accumulated charge carriers. By integrating all these atomically-thin layers in ionic-liquid-gated field-effect transistors (FETs) we have measured the evolution of the Raman spectrum as a function of electron and hole density up to approximately 5$\times$1013 cm*-2*. In all investigated mono and bilayers, the experiments reveal unambiguously that the electron-phonon interaction systematically softens out-of-plane vibrational modes only when electrons are accumulated in the FET channel, while leaving them unaffected upon hole doping. To identify the origin of this unexpected electron-hole asymmetry we have performed density-functional-theory (DFT) calculations in the most realistic framework currently available, accounting for the reduced dimensionality of the materials and for the presence of a finite density of charge carriers Sohier2017 . The calculations allow us to conclude that the observed phonon softening originates from a pronounced increase in the strength of the electron-phonon coupling that occurs whenever the charge carriers in the systems populate simultaneously two inequivalent valleys. Besides accounting for the observed electron-hole asymmetry (as a direct consequence of the different valley structure in the conduction and valence bands of atomically-thin semiconducting TMDs), this finding reveals an aspect of electron-phonon coupling that had not yet been appreciated. Specifically, the strengthening of the electron-phonon coupling is caused by charge transfer between the inequivalent valleys induced by the deformation potential associated to out-of-plane vibrational modes: being local in real space, this inter-valley charge transfer decreases the density of charge that can be displaced, thereby suppressing the effectiveness of electrostatic screening. The phenomenon –expected to be of rather general validity– had not been identified so far, most likely because the models commonly used to describe electron-phonon interaction theoretically do not consider the presence of multiple inequivalent valleys. Equipped with the understanding resulting from the work presented here, we discuss how existing experiments on superconductivity in TMDs (both gate-induced ye_superconducting_2012 ; jo_electrostatically_2015 ; lu_full_2018 and spontaneously occurring Hamaue1986 ; Guillamon2008 ; Xi2016 ; Ugeda2016 ; Tsen2016 ) provide evidence correlating the enhancement in electron-phonon interaction due to multivalley populations to the occurrence of the superconducting state.

II Experimental results

Our experiments consist in measurements of the Raman spectrum of mono and bilayers of different semiconducting TMDs (WSe2, WS2, and MoS2) as a function of carrier density. The density and polarity of charge carriers is varied continuously by employing these atomically-thin crystals as active parts of ionic liquid gated field-effect transistors (see Fig. 1(a) for a schematic illustration). It has been shown in a multitude of experiments over the last several years that ionic-liquid gating is extremely effective in combination with semiconducting TMDs, as it allows the accumulation of large densities (even in excess of 1014 cm*-2*) of both electrons and holes in a same device braga_quantitative_2012 ; ye_superconducting_2012 ; zhang_electrically_2014 ; lezama_surface_2014 ; jo_mono-_2014 ; ubrig_scanning_2014 ; ponomarev_ambipolar_2015 ; jo_electrostatically_2015 ; zhang_ambipolar_2012 ; costanzo_gate-induced_2016 ; gutierrez-lezama_electroluminescence_2016 . One of the characteristic manifestations of the possibility to vary the charge density and type over such a large interval is the occurrence of ambipolar transport, which can be observed upon sweeping the gate voltage VG (see Fig. 2(a) for an example). In the present context, the occurrence of ambipolar transport is useful because it allows us to obtain a fairly accurate estimate of the carrier density in the devices used for the Raman measurements in-situ. In practice, in the experiments we record the Raman spectrum at many different values of applied gate voltage and –as we vary VG– we also apply a small bias voltage, VSD between the source and the drain electrode to measure the current, ISD, passing through the transistor channel. This allows us to determine the threshold voltage for electron and hole conduction –$V_{th}^{\rm e}$ and $V_{th}^{\rm h}$– from which we directly estimate the density of electrons and holes at any given value of VG as $n_{e}=\frac{C|V_{G}-V_{th}^{\rm e}|}{e}$ and $n_{h}=\frac{C|V_{\rm G}-V_{th}^{\rm h}|}{e}$ (where $C$ is the capacitance per unit area obtained from previous experiments on analogous devices, in which Hall effect measurements were performed to determine the density of charge carriers). We estimate the uncertainty in the value of carrier density extracted in this way to be approximately 20-30 %, which is sufficient for the purpose of this work (a more precise determination would require measuring the Hall resistance using a superconducting magnet which cannot be done in the same set-up that we use to perform Raman spectroscopy). For more details concerning the fabrication and operation of ionic-liquid FETs based on atomically-thin TMDs we refer the reader to our earlier work (Refs. braga_quantitative_2012, ; jo_mono-_2014, ; ponomarev_ambipolar_2015, ) and to Appendix A.

Characteristic Raman spectra from bare monolayers of WSe2 (blue line), WS2 (orange line), and MoS2 (green line) are shown in Fig. 1(b). We focus on the most prominent features in these spectra, whose assignment in terms of Raman active vibrational modes has been discussed extensively in the literature lee_anomalous_2010 ; Molina-Sanchez2011 ; berkdemir_identification_2013 ; sahin_anomalous_2013 ; saito_raman_2016 . Using for simplicity the same nomenclature as in the bulk, these are the in-plane (E2g) and the out-of-plane (A1g) modes, as indicated in Fig. 1(b), which are seen in all group VI TMDs. For WS2 and MoS2 these features are spectrally separated while for WSe2 the modes hybridize due to an accidental degeneracy in the phonon dispersion relation, and appear as a single peak at about 250 cm*-1*. The degeneracy is also present in the calculated phonon dispersion relation of WSe2 in Fig. 1(b), on which we point to the modes observed in the Raman spectroscopy data. For most peaks the observed frequency corresponds to the value of the corresponding mode at the $\Gamma$-point. However one of the peaks, visible in the spectra of both WSe2 and WS2 and labeled 2LA(M) (see Fig. 1(b)), is due to a double resonant process terrones_new_2014 ; mitioglu_second-order_2014 , which occurs because the frequency of the laser used in our measurements matches an electronic transition (C absorption for WSe2 and the B-exciton for WS2) carvalho_band_2013 . Although it is not the main point of interest here, the possibility to observe this mode is relevant because it allows the investigation of the coupling between electrons and phonons at finite momentum, i.e. away from center of the Brillouin zone, as indicated by the corresponding label LA(M) in Fig. 1(c).

To illustrate the results of our Raman spectroscopy measurements as a function of accumulated charge density we start by discussing the case of WSe2 monolayers. The source-drain current ISD as a function of VG (for VSD = 50 mV) measured in the device used for the Raman measurements, is presented in Fig. 2(a) and exhibits clear ambipolar behavior. The increase of ISD at positive values of VG indicates accumulation of electrons in the FET channel while the increase of ISD at negative VG corresponds to the accumulation of holes. As explained above, we use these measurements to estimate the density of electrons and holes as a function of VG. Fig. 2(b) shows the Raman spectra measured for selected values of VG indicated by the vertical dashed lines in Fig. 2(a). For positive values of VG electrons are accumulated (green and red lines, i.e. VG = +1.2 and +1.5 V, corresponding respectively to $n_{e}\sim 1.5\times 10^{13}$ and $2.5\times 10^{13}$ cm*-2*), for negative VG values the gate accumulates holes (orange line, i.e. VG = -1.5 V, corresponding to $n_{h}\sim 6\times 10^{13}$ cm*-2*), whereas for VG = 0 V (blue line) the chemical potential is in the gap and the semiconductor is neutral capa . Upon electron accumulation the Raman peak originating from the hybridized A1g/E2g mode exhibits a clear softening, shifting towards lower wavenumbers by an amount $\Delta\omega\simeq 3-4$ cm*-1*, as well as a decrease in intensity as compared to VG = 0 V case. In contrast, accumulation of holes leaves the Raman spectrum virtually unchanged: the peak does not exhibit any appreciable shift, and only its height decreases slightly. The same trends upon varying VG are observed for the resonant Raman peak originating from the 2LA(M) mode, shown in the right panel of Fig. 2(b) 2LAm .

Since in WSe2 monolayers the degeneracy of the E2g and A1g modes prevents us to study separately the evolution of the two individual Raman peaks upon accumulation of charge carriers, in Fig. 2(c) we show data from bilayer WS2, in which these modes are well separated in frequency. The corresponding Raman peak positions at VG = 0 V are located respectively at $\omega$ = 355 cm*-1* –Fig. 2(c), left panel– and 417 cm*-1* –Fig. 2(c), right panel. The same figures also show the Raman spectra measured at selected values of VG, in a range of carrier density somewhat larger than for monolayer WSe2. For values of VG corresponding to hole accumulation only a negligible shift in position, and at most a small change in amplitude, are seen for all peaks. In contrast, upon electron accumulation, the A1g and the 2LA(M) modes exhibit a clear softening similarly to what is observed in WSe2 (the position of the E2g mode remains unaffected even for electron accumulation).

Fig. 2 reveals a number of clear trends of the Raman active modes in WSe2 monolayer and WS2 bilayer, such as the insensitivity of all phonon modes to hole accumulation, and the softening of out-of-plane modes upon electron accumulation. We have performed gate-dependent Raman spectroscopy measurements on mono and bilayers of WSe2, WS2 and MoS2 following a procedure identical to the one described here above for WSe2 monolayer and WS2 bilayer, and found that these trends are generically present in all the investigated atomically thin crystals of semiconducting TMDs. To illustrate and compare quantitatively the results of measurements performed on different systems we plot for each one of them the change $\Delta\omega$ in the frequency of the different phonon modes relative to the value measured at threshold for electron or hole accumulation (i.e., for each system we plot $\Delta\omega$ versus $V_{\rm G}-V_{\rm th}^{\rm e}$ and $V_{\rm G}-V_{\rm th}^{\rm h}$). The results of these measurements for the A1g modes in all the investigated mono- and bilayers are summarized in Fig. 3(a) and 3(b). The data for the E2g modes in monolayers are shown in Fig. 3(c). It is clear from this systematic experimental analysis –which goes well beyond what had been done until now  chakraborty_symmetry-dependent_2012 ; lu_gate-tunable_2017 – that the same behavior is common to all atomically-thin TMDs.

To quantify the strength of the influence of electron and hole accumulation on the frequency of the different phonons, we plot for all systems investigated the frequency shift of the A1g and E2g modes measured at the highest value of applied gate voltage applied (relative to the corresponding threshold voltage), $|V_{\rm G}-V_{\rm th}^{\rm e,h}|=0.6$ V, for which we have obtained data for both mono and bilayer, as well as for electron and hole accumulation. Fig. 3(d) demonstrates clearly that, for both mono and bilayers, only Raman active phonon modes with an out-of-plane component of the atomic displacement are influenced by the presence of charge carriers, whereas the E2g mode –for which the atomic displacement is in-plane– is left unaffected by the presence of free charge. The data also make the strong electron-hole asymmetry in the softening of the out-of-plane modes very apparent for all the different TMDs investigated. The order of magnitude of the effect –which is what we discuss here– is the same in all cases, corresponding to a shift of 2-4 cm*-1* for the A1g optical Raman mode (for a more detailed quantitative analysis it should be recalled that –since the value of the capacitance used is different for mono and bilayers and for electron and hole accumulation– the same value of gate voltage relative to threshold does not correspond to the same carrier density capa ).

III Theoretical analysis

The experiments discussed in the previous section show unambiguously that a strong electron-hole asymmetry in the softening of the phonon frequency is a generic property of atomically thin semiconducting TMDs. The microscopic mechanism at the origin of this asymmetry, however, is not known and, especially for monolayers, finding such a pronounced asymmetry is unexpected. Indeed, the low-energy electronic states in WSe2, WS2, and MoS2 monolayers at the bottom of the conduction band and at the top of the valence bands –centered around the K/K’ points– are often described in terms of massive Dirac fermions, a conceptual framework that may lead one to expect a high degree of symmetry in the physical response upon accumulation of electrons and holes. Even though the very different strength of spin-orbit interaction in the conduction and valence band limits the validity of these considerations, whether and how spin-orbit interaction can cause such a strong asymmetry in the phonon properties upon electron and hole accumulation is very far from obvious.

The purpose of this section is to present a theoretical analysis showing that the strength of the electron-phonon interaction in atomically-thin semiconducting TMDs is governed by whether or not multiple valleys (the $\Gamma$ and K valleys in the valence band and the K and Q valleys in the conduction band, see Fig. 4) are simultaneously populated by the accumulated charge carriers. In other words, it is the rather different valley structure of the conduction and the valence bands that is ultimately responsible for the difference in the coupling to phonons of electrons and holes, and that causes the observed electron-hole asymmetry in the softening of the phonon frequencies. In this regard, spin-orbit interaction does play a prime role, because its strongly asymmetric strength between valence and conduction bands governs the energetic alignment of the different valleys, and therefore determines whether –at any given value of carrier density– one or multiple valleys are populated.

This conclusion is based on DFT and density-functional perturbation theory (DFPT) calculations, which provide an extremely powerful tool to identify all trends in the observed dependence of the phonon frequency on the density of accumulated charge. DFT, however, cannot reproduce the experiments at a detailed quantitative level. This is because a full quantitative agreement would require a precise prediction of the relative energetic alignment of the valleys present in the valence and conduction band, which is well-known to be extremely sensitive to multiple parameters MolinaSanchez2015 ; Le2015 ; Kormanyos2015 ; Yuan2016 , from the choice of exchange-correlation functional to details of the crystal structure. Indeed, even experimental attempts to determine precisely the energetics of the different valleys seem to give scattered results, depending –for instance– on the substrate employed in the experiments (see Table 1). The situation is even more complex for bilayers, where the electric field associated with the FET doping influences the relative energy of the valleys in a way that cannot be precisely predicted at this stage. As such, our analysis aims exclusively at exploiting the possibility to tune parameters and conditions in the controlled environment of DFT simulations to identify the physical processes that explain our experimental results, and not at reproducing in quantitative detail the observed carrier density dependence of the phonon frequencies. That is also why –in view of the subtle effects caused by the gate electric field in bilayers– we focus our theoretical analysis exclusively on TMD monolayers.

To start discussing the effect of the accumulated carrier density on the phonons, we note that the softening observed in the experiments is expected when describing electron-phonon interactions within the Born-Oppenheimer approximation. In such an approximation scheme, electrons are assumed to follow the ionic motion, remaining at all times in the ground state corresponding to the instantaneous lattice configuration. This implies that phonons can be considered as static perturbations acting on the electrons, which in turn affect the vibrational properties of the system. In this regime, the expression for the change in frequency of a phonon at $\Gamma$ with respect to the neutral semiconducting case reads (see appendix B and Refs. Lazzeri2006 ; Saitta2008 ; Calandra2010 ):

[TABLE]

Here we assume the limit of zero temperature, $N(\varepsilon_{F})$ is the density of states at the Fermi energy, $\varepsilon_{F}$, and $\langle g_{\nu}^{2}\rangle_{FS}$ is the screened electron-phonon coupling (squared) for the $\nu$-th phonon mode averaged over Fermi surface (see Eq. (B9) in appendix B). It follows that softening ($\Delta\omega_{\nu}<0$) is expected for phonon modes that are strongly coupled with the electronic states on the Fermi surface, whereas no frequency shift is expected for modes with a negligible electron-phonon coupling. For instance, this is the case of the pure E2g longitudinal optical mode, which affects electrons by generating a long-range scalar potential (known as Fröhlich interaction) that in the presence of free carriers (either electrons or holes) is perfectly screened in the long-wavelength limit, so that $\langle g_{\nu}^{2}\rangle_{FS}\approx 0$.

The assumption that carriers follow the ionic motion remaining in the instantaneous ground state implies that the Born-Oppenheimer approximation is a good description only for systems that are in the adiabatic limit. This requires the phonon energy $\hbar\omega_{\nu}$ to be much smaller than the carrier relaxation rate $\hbar/\tau$ ($\omega_{\nu}\tau\ll 1$), where $\tau$ is the carrier momentum-relaxation lifetime. In ionic-liquid gated FETs based on monolayer semiconducting TMDs the value of $\tau$ can be estimated from the experimentally determined carrier mobility that –under the conditions of the experiments (i.e., with ionic-liquid gated devices)– is typically $\mu\approx$ 15 cm2/Vs for both electrons and holes. When compared to the characteristic phonon frequencies probed in our experiments, the resulting value of $\tau\approx$ 10 fs gives $\omega_{\nu}\tau\leq 1$. The condition to treat electron-phonon interaction in monolayer TMDs using the Born-Oppenheimer approximation is therefore fairly well satisfied at room temperature. Note also that, even when $\omega_{\nu}\tau\simeq 1$, theory allows the phonon softening to be estimated as Maksimov1996 ; Marsiglio1992 ; Saitta2008 :

[TABLE]

where $\Delta\omega_{\nu}$ is given by Eq. (1). It follows directly from this expression that in the regime of our experiments the shift of the phonon frequency upon accumulation of charge carriers is correctly captured by the value in the fully adiabatic limit within a factor of two.

The validity of the adiabatic Born-Oppenheimer approximation makes DFT ideally suited to calculate the phonon frequency as a function of carrier density, because in DFT calculations the ground-state density is obtained by keeping the ionic positions fixed. We thus first investigate which phonon modes are sensitive to accumulation of charge carriers by computing phonon dispersions of neutral, electron-doped and hole-doped monolayer TMDs using DFPT Baroni2001 . Since simulations require overall charge neutrality, the accumulated charge is compensated by the presence of gate electrodes, and the correct boundary conditions for 2D materials are provided by cutting off the Coulomb interactions in the direction perpendicular to the layers Sohier2017 . In these initial calculations, spin-orbit coupling is not included and a relatively large electronic temperature is used to smear the Fermi surface (see methods in appendix A).

Representative DFPT results are shown in Fig. 5 for WS2. When a finite density of carriers is included, we observe a softening in various branches at several phonon wave vectors. In particular, we see that the A1g mode softens at $\Gamma$ both for electron and hole accumulation, and at M for electron accumulation only. Similarly, the frequencies of the in-plane acoustic modes (LA and TA) decrease at K both upon electron and hole accumulation, and at M exclusively for electron accumulation. On the contrary, the E2g mode remains unaffected. For all the phonon modes detected in the experiments, therefore, the simulations succeed in reproducing which modes are affected by the accumulation of charge carriers, and correctly captures the fact that carrier accumulation causes softening of the phonon frequencies. However, the overall agreement between calculations and experiments is poor. The magnitude of the softening is largely overestimated, well beyond the factor $[1+(\omega_{\nu}\tau)^{2}]$ accounting for deviations from the purely adiabatic approximation of DFT. In addition, softening appears also upon hole-doping for the A1g mode at $\Gamma$, in contrast with the experimental results.

To investigate the origin of the discrepancy between measurements and simulations, we focus our attention on the A1g  which is characterized by an out-of-plane displacement of the chalcogen atoms as illustrated in Fig. 1(b) 2LAm . We computed the A1g frequency at $\Gamma$ as a function of carrier density using a finite-difference DFT scheme (see appendix A) for the three monolayer TMDs investigated here. At this stage, we made calculations more realistic by including spin-orbit coupling through fully-relativistic pseudopotentials, and by describing the occupation of electronic states according to the room-temperature Fermi-Dirac distribution (for which we made a much denser sampling of the Brillouin zone). Again, all calculations of total energy and forces include a cutoff of the Coulomb interaction to implement the correct periodic boundary conditions, and gates to maintain overall charge neutrality and properly simulate a FET configuration.

Results for monolayers of the three different TMDs investigated experimentally are shown in Fig. 6, both for electron and hole accumulation, in a range of carrier densities comparable to that explored in the experiments. In contrast to the very systematic results of the Raman measurements (see Fig. 3), the evolution of the calculated phonon frequency upon electron and hole accumulation is different for the individual TMDs, and appears at first sight to offer no systematics. In particular, we find that for WS2 (Fig. 6(a)) a large softening occurs on the electrons side, but the A1g mode frequency decreases also upon hole doping (see also Fig. 5). For WSe2 (Fig. 6(b)) instead a significant softening is present only for electrons and not for holes, whereas the opposite is true for MoS2 (Fig. 6(c)). Despite this seemingly erratic behavior, we succeeded to identify common trends enabling a rationalization of the results of the calculations. To this end, it is essential to realize that the nature of the electronic states involved plays an important role, because the softening of the phonon frequencies arises from the coupling of phonons with electrons on the Fermi surface (see Eq. (1)). That is why in the bottom panels of Fig. 6 we also plot how charge carriers in the different TMD monolayers are distributed among the relevant valleys at K, Q and $\Gamma$ (see Fig. 4), by showing their fractional occupation $P$.

A comparison of the evolution of the phonon frequency with the population of the different valleys does indeed reveal a systematic behavior common to monolayers of all different TMDs, which holds true for both electron and hole accumulation. In all cases a strong softening of the A1g mode occurs only when two inequivalent valleys are simultaneously occupied, whereas the frequency of this mode remains constant whenever the accumulated carriers occupy exclusively one of the valleys. Specifically, the rapid decrease of the A1g frequency upon electron doping in WS2 and WSe2 originates from the simultaneous occupation of the K and Q valleys (see orange squares and blue circles in the bottom panels of Fig. 6 (a) and (b), respectively). Consistently with this idea, in MoS2 only the K-valley is populated, and no significant softening is observed. Analogous trends can be seen upon hole accumulation: A large softening is present in MoS2, for which the K and $\Gamma$ valleys are simultaneously occupied, whereas no softening is visible in WSe2 where holes populate only the K valley (see Fig. 6(c) and (b), respectively). The case of WS2 monolayers (Fig. 6(a)) confirms once more the underlying notion. In this case, the frequency of the A1g mode is insensitive to hole doping for small concentrations –when only the K valley is occupied– and starts to soften past 0.02 holes/unit cell, when holes begin to populate both the K and the $\Gamma$ valley. In all cases a very pronounced softening –larger than 10 cm*-1* at $|n|=5\times 10^{13}$ cm*-2*– is always seen when two distinct valleys are populated, whereas whenever only one valley is populated the change in phonon frequency is typically at least one order of magnitude smaller.

It follows from these considerations that the energy separation between valleys, which governs their relative occupation, plays an essential role in determining the magnitude of the shift in phonon frequencies (i.e., the strength of the electron-phonon coupling) upon carrier accumulation. To confirm this conclusion –and to provide further support for the effect of multivalley populations– we have performed additional calculations exploiting the extreme sensitivity of TMDs’ band structures to the in-plane lattice parameter and the vertical positions of the chalcogen atoms Yuan2016 . The goal of these calculations is to artificially tune the energy separation between valleys in order to check whether, as we expect, the softening of the A1g mode is affected significantly. We have analyzed the behavior of MoS2 monolayers using the experimentally-reported crystal structure (rather than the fully relaxed one obtained by minimizing forces and stresses at the DFT level), for which the energy separation $\Delta E_{\rm\Gamma K}$ between the K and $\Gamma$ valleys in the valence band increases from 60 meV to 220 meV (see Fig. 7(c)), and the energy difference $\Delta E_{\rm KQ}$ between the K and Q valleys in the conduction band decreases from 230 meV to 30 meV. To reduce the spurious effects of thermal broadening in these calculations, we intentionally determine band occupations by using a cold-smeared distribution Marzari1999 with $T=150$ K, instead of a broad Fermi-Dirac distribution at $T=300$ K as above.

The outcome of the calculations is illustrated in Fig. 7, where we compare the carrier-density dependence of (a) the frequency of the A1g mode and (b) the fractional valley occupations obtained with either the DFT-relaxed structure (empty symbols, see also Fig. 6(c)) or the experimentally-reported one (full symbols). It is apparent that the results are distinctly different in the two cases. Upon hole accumulation a large frequency shift and a simultaneous population of the K and $\Gamma$ valleys are present for calculations performed using the DFT-relaxed structure. In contrast, no shift in frequency occurs in the calculations performed with the experimentally-reported crystal structure, where the fractional occupation of the K valley remains equal to one, i.e. only the K valley is populated. This is exactly the behavior that we would have anticipated as a result of the increased $\Delta E_{\rm\Gamma K}$ for the experimentally-reported crystal structure. Consistently, the reduction in $\Delta E_{\rm KQ}$ (see Fig. 7(c)) in the experimentally-reported structure leads to a non-negligible occupation of the Q valley in addition to the lowest-lying K valley, causing a very large enhancement of the A1g softening upon electron accumulation.

We therefore conclude that the scenario hypothesized above is correct: for all monolayers investigated here, and for both electron and hole accumulation, a very sizable softening of the A1g mode occurs only when charge carriers populate multiple inequivalent valleys. Through Eq. (1), this directly implies that also the strength of the electron-phonon coupling for the A1g mode in semiconducting TMDs is large only when two valleys are occupied. This conclusion is particularly interesting because such a strong dependence of the electron-phonon coupling strength on the population of multiple inequivalent valleys in an individual band had not been identified earlier. It is nevertheless extremely likely that this mechanism is at play in many other material systems besides semiconducting TMDs, and relevant to explain a variety of physical phenomena (see, e.g., the discussion of superconductivity in section IV).

The insight gained on the importance of the simultaneous occupation of multiple valleys can now be exploited to discuss in more detail the results of our Raman spectroscopy experiments (see Fig. 3). To start with, the complete absence of phonon softening upon hole accumulation indicates that $\Delta E_{\rm\Gamma K}$ is sufficiently large in TMD monolayers to avoid the partial occupation of the $\Gamma$ valley (in addition to the dominant K valley), for hole concentrations as large as $\sim 5\times 10^{13}$ cm*-2*. According to the same logic, the phonon softening observed upon electron accumulation indicates that $\Delta E_{\rm KQ}$ in the conduction band is sufficiently small, so that –at room temperature– both the K and Q valleys are populated in monolayers of all TMDs upon accumulation of $\sim 10^{-2}$ electrons per unit cell. Additionally, it is clear from Fig. 3(a) that phonon softening upon electron accumulation is much more pronounced in W-based TMD monolayers, indicating that $\Delta E_{\rm KQ}$ is smaller for WS2 and WSe2 than for MoS2.

All these conclusions are fully consistent with what is known about monolayers of semiconducting TMDs. The energy difference $\Delta E_{\rm\Gamma K}$ reported in Tab. 1 from ARPES measurements is consistent with having only the K valley occupied upon accumulation of a density of holes reaching up to $\sim 5\times 10^{13}$ cm*-2*  Zhang2016 ; Gao2017 ; Katoch2018 . Similarly, although no sufficiently systematic ARPES study of the conduction band of monolayer TMDs has been reported yet, analogous measurements on the doped surface of bulk TMDs (where doping should be limited to the first few layers) have shown Kang2017 that it is relatively easy to populate both the K and Q valleys in the conduction band (i.e., $\Delta E_{\rm KQ}$ in monolayer TMDs is significantly smaller than $\Delta E_{\rm\Gamma K}$). Moreover, the observation of a larger softening in W-based monolayers, associated with the K and Q valleys being closer in WS2 and WSe2 than in MoS2, is compatible with the fact that $\Delta E_{\rm KQ}$ is largely controlled by spin-orbit coupling (stronger in W-based compounds).

These considerations also make clear why the strong sensitivity of the phonon frequencies to the precise energy separation between valleys severely affects the ability of DFT to make reliable quantitative predictions. In particular, the spurious softening predicted on the hole side of MoS2 (see Fig. 6(c)) originates from the large underestimation of $\Delta E_{\rm\Gamma K}$ in DFT as compared to experiments (see Tab. 1), which leads to an erroneous occupation of both the K and $\Gamma$ valleys even at very small doping concentrations at room temperature. Something similar, although not as severe, happens also for WS2, for which the DFT prediction for $\Delta E_{\rm\Gamma K}$ is not much smaller than the experimental results in Tab. 1, but it is still sufficient to lead to a spurious softening at large hole doping. Finally, the discrepancy with experiments might arise also from a partial failure of DFT in accounting for band-structure renormalization effects that influence the doping dependence of $\Delta E_{\rm\Gamma K}$ and $\Delta E_{\rm KQ}$.

IV Discussion

Finding that the strength of the coupling of electrons to phonons is drastically amplified when multiple inequivalent valleys in a given band are populated is a phenomenon that had not been identified in the past, and it is important to build a robust physical understanding. To this end, recall that electron-phonon coupling can be qualitatively understood by considering how the band structure is perturbed by the displacement pattern of the ions associated to the phonon under consideration Bardeen1950 ; Schockley1950 ; Khan1984 . Fig. 8 shows the effect of the A1g mode on the band structure of undoped monolayer MoS2 (for better clarity, spin-orbit coupling is not included; we checked that this has no qualitative consequences on the arguments here below). Although the amplitude (up to 0.12 Å) of the A1g mode has been intentionally exaggerated with respect to typical room-temperature mean-square displacements to illustrate this effect, the large variation of both the valence and conduction bands signals the presence of a very strong electron-phonon coupling. This strength can be quantified by taking the first derivative of the band energy reported with respect to the atomic displacement, which is proportional to the bare electron phonon coupling for that band (see the right panel of Fig. 8).

The introduction of free carriers –in our experiments by means of field-effect doping– allows for the possibility of electrostatic screening, and generically tends to decrease the strength of the electron-phonon coupling (as discussed in Section I). In a simple Thomas-Fermi picture, whether or not the bare electron-phonon coupling is screened effectively by the free carriers depends on whether the evolution of the band induced by the phonon displacement pattern leads to a spatial variation of the density of charge carriers. In this regard, the situation can be very different depending on whether one or more valleys are populated.

To understand why, consider a phonon near $\Gamma$ ($\mathbf{q}\to 0$) as a frozen atomic displacement with very long, but finite, wavelength. The static displacement associated to such a phonon shifts the valleys up and down in energy locally in space, following the local amplitude of the long-wavelength A1g mode. In the adiabatic limit, the population of the electronic states will change to keep the system in the lowest energy state. Now, if only a single valley is occupied (say the K valley in the conduction band), the only possibility for the free carriers is to follow the spatially varying position of the bottom of the valley to keep the Fermi level constant (or more appropriately the electrochemical potential, since we are considering the atomic displacement as frozen). This implies a spatial variation of the charge of the free carriers that self-consistently screens the potential of the lattice displacement, thereby strongly suppressing electron-phonon coupling.

When at least two valleys are occupied, the situation can be very different if the atomic displacements associated with the phonon mode lead to a relative out-of-phase energy shift. This is indeed so in monolayer TMDs for the Q and K valleys in the conduction band and for the K and $\Gamma$ valley in the valence band (see Fig. 8). In this case, in the presence of the atomic displacement (and when electrons/holes have the time to relax their energy and momentum during the phonon perturbation, i.e. in the adiabatic limit) charge transfer between the valleys must also occur locally in real space for the system to stay in the ground state. Since part of the charge is now involved in this “local” process, the charge that can be displaced to screen the spatially varying potential generated by the phonon is reduced, hindering the effectiveness of screening. In other words, in the presence of two valleys that shift out-of-phase in response to long-wavelength phonon excitations, screening becomes much less effective in reducing the strength of the electron-phonon coupling.

This is the main mechanism responsible for the strong softening of the A1g mode that we observe in the Raman measurements in all cases when two valleys are simultaneously occupied. For this mechanism to have a strong effect, it is essential that the phonon mode considered causes an out-of-phase shift in the the energy of the two valleys, because this maximizes the amount of inter-valley charge transfer that occurs locally in real space. For the E2g mode detected in our Raman measurements on semiconducting TMDs, for instance, this is not the case: the long-range Fröhlich potential associated with the E2g phonon near $\Gamma$ shifts the energy of the two valleys in the same direction and by comparable amounts so that no (or minimal) inter-valley charge transfer occurs locally. As a result, despite the presence of multiple inequivalent valleys, the response of the system to E2g phonons is similar to the canonical case with one single valley populated, with a strong electrostatic screening that suppresses the electron-phonon coupling.

It follows from these considerations that the underlying physical processes responsible for the enhancement of electron-phonon coupling in the presence of simultaneously populated inequivalent valleys is of quite general validity. As such, it is likely to occur in a multitude of materials and to have different experimental manifestations. It is quite remarkable that such a mechanism had not been identified earlier. We believe that this is because the analytic theoretical models commonly employed to treat electron-phonon interaction consider only one family of carriers (i.e., one single valley) coupled to phonons, and in these models the mechanism is simply not present.

Superconductivity provides an example of a physical phenomenon for which the enhanced electron-phonon coupling mechanism that we have identified may be at play. It appears that much of the existing data on gate-induced superconductivity observed in some of the TMDs (e.g., WS2 and MoS2) is consistent with having the superconducting instability setting in only when multiple valleys are occupied. Having multiple valleys occupied would explain, for instance, why gate-induced superconductivity in WS2 starts occurring at a lower electron density than in MoS2 ye_superconducting_2012 ; lu_full_2018 ; costanzo_gate-induced_2016 , because, as we discussed here above, the energy difference between the K and Q valley is smaller in WS2 than in MoS2. It also explains why so far gate-induced superconductivity in semiconducting TMDs has been observed only upon electron accumulation and not with holes ye_superconducting_2012 ; jo_electrostatically_2015 : this is because with the hole concentrations that can be reached by field-effect doping only the K or $\Gamma$ valleys are filled, so that the strength of the electron-phonon coupling remains weak. Indeed, in NbSe2 –a metal having a similar valence band structure as WS2 and MoS2, but the Fermi level very deep relative to the top of the valence band Lebegue2009 ; Nakata2018 – robust superconductivity is observed Ugeda2016 ; Xi2016 ; Tsen2016 and both the K and the $\Gamma$ valleys are occupied. The multi-valley mechanism responsible for the enhancement of the electron-phonon coupling is likely at play for other superconducting materials, with MgB2 providing a notable example. Indeed, also in MgB2 multiple inequivalent valleys are known to be simultaneously occupied and to move out-of-phase in energy when atoms are displaced according to a specific phonon pattern An2001 . More work is needed to fully substantiate the relevance of multivalley population for the occurrence of superconductivity in all these systems, but the known experimental phenomenology does provide significant circumstantial evidence supporting such a scenario.

V conclusions

In conclusion, we have identified a strong and systematic electron-hole asymmetry in the vibrational properties of mono and bilayers of semiconducting TMDs, by performing Raman measurements on ionic-liquid-gated transistors. The experiments very systematically show that, in these systems, out-of-plane modes –such as the A1g mode near zone center– soften significantly as the concentration of free electrons in the conduction band is increased, while they are left unaffected upon hole accumulation. By performing a theoretical analysis based on first-principles simulations, we have revealed that the phonon softening originates from a previously unidentified mechanism that amplifies pronouncedly the strength of electron-phonon coupling when two inequivalent valleys are simultaneously populated. This mechanism explains the asymmetry in phonon softening observed experimentally in the Raman spectroscopy measurements as due to a much larger energy separation between valleys in the valence bands as compared to the conduction band. As a result, in the range of carrier densities accessed in the experiments, simultaneous occupation of multiple valleys occurs only upon electron accumulation and not for holes.

It is surprising that this mechanism had not been identified earlier. This is most likely due to the fact that analytic theoretical models used in the discussion of electron-phonon interactions commonly treat the case of a single family of carriers (i.e., a single valley) that is coupled the lattice vibrations. It is nevertheless clear from our discussion that the microscopic processes responsible for the enhancement of electron-phonon interaction in multi-valley systems are physically robust. Hence, the same mechanism should be at play in many different material systems (i.e., not only in semiconducting transition metal dichalcogenides), and manifests itself in a variety of physical phenomena. This is likely the case, for example, of superconductivity. Existing experiments on transition metal dichalcogenides appear to indicate that the occurrence of superconductivity –either gate-induced in MoS2 and WS2, or intrinsically present in NbSe2– correlates systematically with having multiple valley populated. Similar considerations can be made for other materials, with the case of MgB2 being particularly notable. These examples underscore how the identification of a previously unappreciated mechanism enhancing the strength of electron-phonon interaction can have an important impact in the interpretation of many different phenomena observed experimentally.

Acknowledgements.

We sincerely thank Francesco Mauri, Matteo Calandra, Haijing Zhang, and Jérémie Teyssier for fruitful discussions. We gratefully acknowledge Alexandre Ferreira for continuous and precious technical support. A.F.M. gratefully acknowledges financial support from the Swis National Science Foundation (Division II) and from the EU Graphene Flagship project. M.G. and N.U. acknowledge support from the Swiss National Science Foundation through the Ambizione program. Simulation time was provided by CSCS on Piz Daint (project ID s825) and by PRACE on Marconi at CINECA, Italy (project ID 2016163963).

Appendix A Methods

Te realization of the field-effect transistors used to perform the experiments described in the main text relies on a combination of conventional nanofabrication processes and of techniques that are commonly employed to manipulate atomically thin crystals, i.e. 2D materials. Monolayers and bilayers of MoS2, WS2 and WSe2 are obtained by mechanical exfoliation of bulk crystals and transferred onto Si/SiO2 substrates. The layer thickness is identified by looking at optical contrast under an optical microscope, and confirmed by photoluminescence measurements. To attach electrical contacts to the atomically thin layers, a conventional process based on electron-beam lithography, electron-beam evaporation and lift-off is used (together with the electrodes, also a large side gate electrode is deposited). The contacts consist of Pt/Au evaporated films for WSe2 and Au evaporated films for WS2 and MoS2, and are annealed at 200${}^{\circ}\,$C in a flow of Ar:H2 gas (100:10 sccm). As a last step, the ionic liquid is drop-casted on top of the structure and confined under a 50 $\mu$m thick glass plate to ensure a sharp optical image during the optical measurements. As rapidly as possible after depositing the liquid, the device is mounted in the vacuum chamber with optical access (Cryovac KONTI cryostat) where the measurements are performed (p $\approx$ 10*-6* mbar). The devices are left under vacuum at least overnight prior to the application of any gate bias to ensure that moisture or oxygen possibly present are pumped out of the system. The ionic liquid used in all experiments is 1-butyl-1-methylpyrrolidiniumtris(penta-fluoroethyl) trifluorophosphate [P14][FAP]. All measurements discussed here have been performed at room temperature; at lower temperature the same procedure cannot be applied because of freezing of the [P14][FAP] ionic liquid.

The Raman scattering spectroscopy measurements are performed using a home tailored confocal microscope in a back-scattering geometry, i.e. collecting the emitted light with the same microscope used to couple the laser beam (at 514.5 nm) onto the device (see Fig. 1c for a sketch of the experimental setup). The Back-scattered light is sent to a Czerni-Turner spectrometer equipped with a 1800 groves/mm grating which resolves the optical spectra with a precision of 0.5 cm*-1*. The light is detected with the help of a N2-cooled CCD-array (Princeton Instruments). In the experiments, the gate voltage is slowly swept from positive to negative values. At VG intervals of 0.2 V the sweep is interrupted and Raman spectra are recorded at fixed VG values. After reaching the largest negative gate voltage, VG is swept back to positive values, and spectra are again recorded at fixed intervals to ensure that the experimental data are free from hysteresis. Finally, we have checked that the ionic liquid does not possess any characteristic Raman feature in the spectral range of interest.

DFT calculations are performed using Quantum ESPRESSO Giannozzi2009 ; Giannozzi2017 , including a cutoff of the Coulomb interaction to implement the correct boundary conditions and gates to simulate doping Sohier2017 . We use a symmetric double-gate geometry; for monolayers a single-gate geometry leads to minimal changes on the results. The exchange-correlation functional is approximated using the Perdew-Burke-Ernzerhof (PBE) form PBE of the generalized-gradient approximation. To compute full phonon dispersions, we use the SSSP Accuracy library prandini_standard_2018 ; Prandini2018 (version 0.7), a $60$ Ry energy cutoff, and a $32\times 32\times 1$ grid with $0.02$ Ry Marzari-Vanderbilt broadening Marzari1999 . We then use density-functional perturbation theory Baroni2001 –including Coulomb cutoff and gates– to compute the phonons on a $12\times 12\times 1$ grid in the Brillouin zone. For zone-center phonon calculations, we significantly reduce smearing and increase the sampling accordingly: Namely, we use a dense $120\times 120\times 1$ k-point grid and Fermi-Dirac smearing corresponding to room temperature (0.002 Ry). Spin-orbit interaction is included by using fully relativistic norm-conserving pseudopotentials from the Pseudo-Dojo library vanSetten2018 . Phonopy Togo2015 is then used to generate the atomic-displacement patterns and to post-process forces to obtain phonon frequencies.

Appendix B Justification of Eq. (1) in the main text

In the the following, we show how to use expressions derived in Refs. Lazzeri2006 ; Saitta2008 ; Calandra2010 to obtain Eq. (1) with the help of few assumptions valid in the case of semiconducting TMDs. We are interested in the doping-induced frequency variations of a phonon mode at $\Gamma$ (phonon momentum $\mathbf{q}\to 0$) with respect to its value in the neutral case $\Delta\omega=\omega_{\rm{doped}}-\omega_{\rm{neutral}}$. The frequency of a given phonon mode corresponding to an eigenvector $|\varepsilon\rangle\,$ of the dynamical matrix at $\Gamma$, $\mathcal{D}(\mathbf{q}=0)$, can be obtained as $\omega^{2}=\langle\varepsilon|\,\mathcal{D}(\mathbf{q}=0)|\varepsilon\rangle\,/\langle\varepsilon|\varepsilon\rangle=\mathcal{D}_{\Gamma}/M$, where $M=\langle\varepsilon|\varepsilon\rangle$ is an effective mass associated with the phonon mode and $\mathcal{D}_{\Gamma}=\langle\varepsilon|\,\mathcal{D}(\mathbf{q}=0)|\varepsilon\rangle\,=\langle u|\,\mathcal{C}(\mathbf{q}=0)|u\rangle\,$ is the projection of the dynamical matrix on the corresponding eigenvector or equivalently of the force constants $\mathcal{C}$ on the physical displacements $|u\rangle\,$ of the phonon mode (with $\langle u|u\rangle=1$). Assuming that the variation of the frequency is small compared to the neutral frequency $\Delta\omega\ll\omega_{\rm neutral}$, we can differentiate ($d\omega\approx d\omega^{2}/2\omega$) to obtain the following simple relation:

[TABLE]

Neglecting the variations with doping of the contributions from the second-order derivatives of the ionic potential and the exchange-correlation kernel with respect to the atomic displacements in Eq. (3) of Ref. Saitta2008 , we write:

[TABLE]

where $\mathcal{F}_{\Gamma}$ is the part of the dynamical matrix at the $\Gamma$ point involving the electronic response to the phonon perturbation, and thus electron-phonon interactions. The explicit expression for $\mathcal{F}_{\Gamma}$ depends on the adiabaticity of the system. In the adiabatic limit relevant for TMDs, it reads Saitta2008 ; Lazzeri2006 ; Calandra2010 :

[TABLE]

where $\Delta V_{\mathbf{q},\nu}$ is the perturbation of the effective potential due to phonon $(\mathbf{q},\nu)$, $f_{\mathbf{k},n}$ is the occupation of the state $|\mathbf{k}n\rangle$ corresponding to the band $n$ at crystal momentum $\mathbf{k}$ with energy $\varepsilon_{\mathbf{k},n}$. In the first term, there are only interband contributions ($n\neq m$) between occupied and empty states. We have verified that the numerator is much smaller than the denominator $|\varepsilon_{\mathbf{k},m}-\varepsilon_{\mathbf{k},n}|$, which is equal to the direct gap in the neutral case, so that this term can be neglected. The second term involves intraband transitions between states close to the Fermi surface and it is non-zero only in doped systems, while vanishing in neutral TMDs. This is the main contribution to doping-induced variations, so that we can write:

[TABLE]

where we have introduced the square electron-phonon matrix element due to intraband transitions:

[TABLE]

The derivative of the Fermi distribution is negative (the phonon frequency softens) and non-zero only in a small energy window of order $k_{B}T$ around the chemical potential, so that only electronic states belonging to the temperature-smeared Fermi surface contribute. In the zero temperature limit the derivative becomes a delta function at the Fermi energy $\varepsilon_{F}$ and we obtain

[TABLE]

which is Eq. (1) of the main text. Here we have introduced the density of states at the Fermi energy

[TABLE]

and the square electron-phonon coupling averaged over the Fermi surface:

[TABLE]