Enhanced electron-phonon interaction in multi-valley materials

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)

We report a new mechanism responsible for the enhancement of electron-phonon coupling in doped semiconductors in which multiple inequivalent valleys are simultaneously populated. We have identified this mechanism through a combined experimental and theoretical investigation of the vibrational properties of atomically thin group-VI semiconducting transition metal dichalcogenides as a function of carrier density. Specifically, we have performed Raman spectroscopy experiments to measure different phonon modes of mono and bilayer MoS2, WS2, and WSe2 embedded in ionic-liquid-gated field-effect transistors, and probed vibrations of these systems over a wide range of carrier densities (∼ 5 × 10 13 cm −2 , both for electron and hole accumulation). We find that in-plane modes, such as the E2g mode, are insensitive to the introduction of charge carriers, whereas phonons with a dominant out-of-plane character exhibit strong softening upon electron accumulation, while remaining unaffected upon hole doping. This unexpected -but very pronounced-electron-hole asymmetry is systematically observed in all mono and bilayers. To understand its origin, we have performed first-principles simulations within a framework that accounts for the reduced dimensionality of the systems and gate-induced charge accumulation. Through this analysis we establish that the phonon softening occurs when multiple inequivalent valleys are populated simultaneously and that the phenomenon is a manifestation of an aspect of electron-phonon interaction that had not been previously identified. Accordingly, the observed electron-hole asymmetry in the evolution of the Raman spectra originates from the the much larger energy separation between valleys in the valence bands -as compared to the conduction band-that prevents the population of multiple valleys upon hole accumulation. We discuss the origin of the enhancement of the electron-phonon coupling in terms of a simple physical picture, and conclude that the phenomenon occurs because the population of multiple valleys acts to strongly reduce the efficiency of electrostatic screening for those phonon modes that cause the energy of the inequivalent valleys to oscillate out of phase. The strong enhancement of electron-phonon coupling in the presence of multivalley populations is a robust mechanism that is likely to play an important role in accounting for different physical phenomena and, based on existing experiments, we argue for its relevance in the occurrence of superconductivity -either gate-induced or spontaneous-in different transition metal dichalcogenides.

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 [1][2][3]. 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 [4,5]), 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 [6,7]. In these cases, charge carriers interact with phonons more similarly to the way they interact with a magnetic field [8]. 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 [9,10].
Irrespective of the specific situation, what is key to understand electron-phonon interaction (and its strength) in solids is electrostatic screening [3]. 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 [11][12][13]. Screening is expected to be much more effective in the adiabatic regimesince 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 [14,15], screening cannot influence electronphonon coupling [7,16] (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 electronphonon interaction to decrease, as observed in a variety of semiconductors (e.g., SrTiO 3 [17][18][19] and other transition metal oxides [20]). 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 [7]. 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 weakeningof 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; MoS 2 , WS 2 , and WSe 2 ) as a function of density of accumulated charge carriers. By integrating all these atomically-thin layers in ionic-liquidgated 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×10 13 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 electronhole asymmetry we have performed density-functionaltheory (DFT) calculations in the most realistic framework currently available, accounting for the reduced di-mensionality of the materials and for the presence of a finite density of charge carriers [21]. 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 outof-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 [22][23][24] and spontaneously occurring [25][26][27][28][29]) 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 (WSe 2 , WS 2 , and MoS 2 ) as a function of carrier density. The density and polarity of charge carriers is varied continuously by employing these atomicallythin crystals as active parts of ionic liquid gated fieldeffect 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 10 14 cm −2 ) of both electrons and holes in a same device [22,23,[30][31][32][33][34][35][36][37][38]. 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 V G (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 V G -we also apply a small bias volt- age, V SD between the source and the drain electrode to measure the current, I SD , passing through the transistor channel. This allows us to determine the threshold voltage for electron and hole conduction -V e th and V h thfrom which we directly estimate the density of electrons and holes at any given value of V G as n e = C|V G −V e th | e and n h = (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. 30, 33, and 35) and to Appendix A.
Characteristic Raman spectra from bare monolayers of WSe 2 (blue line), WS 2 (orange line), and MoS 2 (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 [39][40][41][42][43]. Using for simplicity the same nomenclature as in the bulk, these are the in-plane (E 2g ) and the out-of-plane (A 1g ) modes, as indicated in Fig. 1(b), which are seen in all group VI TMDs. For WS 2 and MoS 2 these features are spectrally separated while for WSe 2 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 WSe 2 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 Γ-point. However one of the peaks, visible in the spectra of both WSe 2 and WS 2 and labeled 2LA(M) (see Fig. 1(b)), is due to a double resonant process [44,45], which occurs because the frequency of the laser used in our measurements matches an electronic transition (C absorption for WSe 2 and the B-exciton for WS 2 ) [46]. 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 WSe 2 monolayers. The source-drain current I SD as a function of V G (for V SD = 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 I SD at positive values of V G indicates accumulation of electrons in the FET channel while the increase of I SD at negative V G 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 V G . Fig. 2(b) shows the Raman spectra measured for selected values of V G indicated by the vertical dashed lines in Fig. 2(a). For positive values of V G electrons are accumulated (green and red lines, i.e. V G = +1.2 and +1.5 V, corresponding respectively to n e ∼ 1.5 × 10 13 and 2.5 × 10 13 cm −2 ), for negative V G values the gate accumulates holes (orange line, i.e. V G = -1.5 V, corresponding to n h ∼ 6 × 10 13 cm −2 ), whereas for V G = 0 V (blue line) the chemical potential is in the gap and the semiconductor is neutral [47]. Upon electron accumulation the Raman peak originating from the hybridized A 1g /E 2g mode exhibits a clear softening, shifting towards lower wavenumbers by an amount ∆ω 3 − 4 cm −1 , as well as a decrease in intensity as compared to V G = 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 V G are observed for the resonant Raman peak originat- ing from the 2LA(M) mode, shown in the right panel of Fig. 2(b) [48]. Since in WSe 2 monolayers the degeneracy of the E 2g and A 1g 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 WS 2 , in which these modes are well separated in frequency. The corresponding Raman peak positions at V G = 0 V are located respectively at ω = 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 V G , in a range of carrier density somewhat larger than for monolayer WSe 2 . For values of V G 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 A 1g and the 2LA(M) modes exhibit a clear softening similarly to what is observed in WSe 2 (the position of the E 2g mode remains unaffected even for electron accumulation). Fig. 2 reveals a number of clear trends of the Raman active modes in WSe 2 monolayer and WS 2 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 gatedependent Raman spectroscopy measurements on mono and bilayers of WSe 2 , WS 2 and MoS 2 following a procedure identical to the one described here above for WSe 2 monolayer and WS 2 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 ∆ω 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 ∆ω versus V G − V e th and V G − V h th ). The results of these measurements for the A 1g modes in all the investigated mono-and bilayers are summarized in Fig. 3(a) and 3(b). The data for the E 2g 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 [49,50]-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 A 1g and E 2g modes measured at the highest value of applied gate voltage applied (relative to the corresponding threshold voltage), |V G − V e,h th | = 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 outof-plane component of the atomic displacement are influenced by the presence of charge carriers, whereas the E 2g 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 soft- ening 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 A 1g 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 [47]).

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 WSe 2 , WS 2 , and MoS 2 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 spinorbit 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 Γ 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 densityfunctional 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 [51][52][53][54], 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 I). 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 Γ with respect to the neutral semiconducting case reads (see appendix B and Refs. [55][56][57]): Here we assume the limit of zero temperature, N (ε F ) is the density of states at the Fermi energy, ε F , and g 2 ν F S is the screened electron-phonon coupling (squared) for the ν-th phonon mode averaged over Fermi surface (see Eq. (B9) in appendix B). It follows that softening (∆ω ν < 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 E 2g 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 g 2 ν F S ≈ 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 ω ν to be much smaller than the carrier relaxation rate /τ (ω ν τ 1), where τ is the carrier momentum-relaxation lifetime. In ionic-liquid gated FETs based on monolayer semiconducting TMDs the value of τ 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 µ ≈ 15 cm 2 /Vs for both electrons and holes. When compared to the characteristic phonon frequencies probed in our experiments, the resulting value of τ ≈ 10 fs gives ω ν τ ≤ 1. The condition to treat electronphonon interaction in monolayer TMDs using the Born-Oppenheimer approximation is therefore fairly well satisfied at room temperature. Note also that, even when ω ν τ 1, theory allows the phonon softening to be estimated as [10,56,58]: where ∆ω ν 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 [66]. 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 [21]. 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 WS 2 . 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 A 1g mode softens at Γ 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 E 2g 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 + (ω ν τ ) 2 ] accounting for deviations from the purely adiabatic approximation of DFT. In addition, softening appears also upon hole-doping for the A 1g mode at Γ, in contrast with the experimental results.
To investigate the origin of the discrepancy between measurements and simulations, we focus our attention on the A 1g which is characterized by an out-of-plane displacement of the chalcogen atoms as illustrated in Fig. 1(b) [48]. We computed the A 1g frequency at Γ 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 WS 2 (Fig. 6(a)) a large softening occurs on the electrons side, but the A 1g mode frequency decreases also upon hole doping (see also Fig. 5). For WSe 2 ( Fig. 6(b)) instead a significant softening is present only for electrons and not for holes, whereas the opposite is true for MoS 2 (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 Γ (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 A 1g mode occurs only when two inequivalent valleys are simultaneously occupied, whereas the frequency of this mode remains constant whenever the accumulated carri-ers occupy exclusively one of the valleys. Specifically, the rapid decrease of the A 1g frequency upon electron doping in WS 2 and WSe 2 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 MoS 2 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 MoS 2 , for which the K and Γ valleys are simultaneously occupied, whereas no softening is visible in WSe 2 where holes populate only the K valley (see Fig. 6(c) and (b), respectively). The case of WS 2 monolayers (Fig. 6(a)) confirms once more the underlying notion. In this case, the frequency of the A 1g 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 Γ valley. In all cases a very pronounced softening -larger than 10 cm −1 at |n| = 5×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 populationswe have performed additional calculations exploiting the extreme sensitivity of TMDs' band structures to the inplane lattice parameter and the vertical positions of the chalcogen atoms [54]. 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 A 1g mode is affected significantly. We have analyzed the behavior of MoS 2 monolayers using the experimentallyreported crystal structure (rather than the fully relaxed one obtained by minimizing forces and stresses at the DFT level), for which the energy separation ∆E ΓK between the K and Γ valleys in the valence band increases from 60 meV to 220 meV (see Fig. 7(c)), and the energy difference ∆E 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 [67] 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 A 1g 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  Fig. 2(a)). For WSe2 and MoS2 we also plot the results for the E2g mode (orange empty circles), which is found not to shift, in agreement with experiments (for WS2 the E2g mode is degenerate with the A1g mode which experimentally prevents its separate determination upon varying carrier density). Lower panels: relative fractional occupation P of the different valleys upon hole and electron accumulation, as extracted from ab-initio calculations done to obtain the density dependence of the phonon frequencies shown in the corresponding top panels. Different symbols refer to the K (orange squares), Q (blue circles), and Γ (green triangles) valleys. Note how in all cases -i.e., for all the different monolayers and upon both electron and hole accumulation-there is a perfect correlation between the shift of the A1g mode and the multiple occupation of valleys (i.e., the K and Γ valley in the valence band upon hole accumulation, and the K and Q valley in the conduction band upon electron accumulation).
and a simultaneous population of the K and Γ valleys are present for calculations performed using the DFT-relaxed structure. In contrast, no shift in frequency occurs in the calculations performed with the experimentallyreported 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 ∆E ΓK for the experimentally-reported crystal structure. Consistently, the reduction in ∆E KQ (see Fig. 7(c)) in the experimentally-reported structure leads to a nonnegligible occupation of the Q valley in addition to the lowest-lying K valley, causing a very large enhancement of the A 1g 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 A 1g 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 A 1g 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 super-conductivity 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 ∆E ΓK is sufficiently large in TMD monolayers to avoid the partial occupation of the Γ valley (in addition to the dominant K valley), for hole concentrations as large as ∼ 5 × 10 13 cm −2 . According to the same logic, the phonon softening observed upon electron accumulation indicates that ∆E KQ in the conduction band is sufficiently small, so that -at room temperatureboth the K and Q valleys are populated in monolayers of all TMDs upon accumulation of ∼ 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 ∆E KQ is smaller for WS 2 and WSe 2 than for MoS 2 .
All these conclusions are fully consistent with what is known about monolayers of semiconducting TMDs. The energy difference ∆E ΓK reported in Tab. I from ARPES measurements is consistent with having only the K valley occupied upon accumulation of a density of holes reaching up to ∼ 5 × 10 13 cm −2 [64,68,69]. 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 [70] that it is relatively easy to populate both the K and Q valleys in the conduction band (i.e., ∆E KQ in monolayer TMDs is significantly smaller than ∆E ΓK ). Moreover, the observation of a larger softening in W-based monolayers, associated with the K and Q valleys being closer in WS 2 and WSe 2 than in MoS 2 , is compatible with the fact that ∆E 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 MoS 2 (see Fig. 6(c)) originates from the large under-estimation of ∆E ΓK in DFT as compared to experiments (see Tab. I), which leads to an erroneous occupation of both the K and Γ valleys even at very small doping concentrations at room temperature. Something similar, although not as severe, happens also for WS 2 , for which the DFT prediction for ∆E ΓK is not much smaller than the experimental results in Tab. I, 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 ∆E ΓK and ∆E 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 [4,5,12]. Fig. 8 shows the effect of the A 1g mode on the band structure of undoped monolayer MoS 2 (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 A 1g mode has been intentionally exaggerated with respect to typical room-temperature meansquare 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 Γ (q → 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 longwavelength A 1g mode. In the adiabatic limit, the pop- ulation 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 Γ 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 screen-ing. 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 A 1g 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 E 2g 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 E 2g phonon near Γ 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 E 2g 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 gateinduced superconductivity observed in some of the TMDs (e.g., WS 2 and MoS 2 ) 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 WS 2 starts occurring at a lower electron density than in MoS 2 [22,24,37], because, as we discussed here above, the energy difference between the K and Q valley is smaller in WS 2 than in MoS 2 . It also explains why so far gate-induced superconductivity in semiconducting TMDs has been observed only upon electron accumulation and not with holes [22,23]: this is because with the hole concentrations that can be reached by fieldeffect doping only the K or Γ valleys are filled, so that the strength of the electron-phonon coupling remains weak. Indeed, in NbSe 2 -a metal having a similar valence band structure as WS 2 and MoS 2 , but the Fermi level very deep relative to the top of the valence band [71,72]-ro-bust superconductivity is observed [27][28][29] and both the K and the Γ valleys are occupied. The multi-valley mechanism responsible for the enhancement of the electronphonon coupling is likely at play for other superconducting materials, with MgB 2 providing a notable example. Indeed, also in MgB 2 multiple inequivalent valleys are known to be simultaneously occupied and to move out-ofphase in energy when atoms are displaced according to a specific phonon pattern [73]. 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 A 1g 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 electronphonon 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 multivalley 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 MoS 2 and WS 2 , or intrinsically present in NbSe 2 -cor-relates systematically with having multiple valley populated. Similar considerations can be made for other materials, with the case of MgB 2 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.

ACKNOWLEDGMENTS
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 MoS 2 , WS 2 and WSe 2 are obtained by mechanical exfoliation of bulk crystals and transferred onto Si/SiO 2 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 electronbeam 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 WSe 2 and Au evaporated films for WS 2 and MoS 2 , and are annealed at 200 • C in a flow of Ar:H 2 gas (100:10 sccm). As a last step, the ionic liquid is drop-casted on top of the structure and confined under a 50 µ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 ≈ 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(pentafluoroethyl) trifluorophosphate [P14] [FAP]. All measure-ments 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 N 2 -cooled CCD-array (Princeton Instruments). In the experiments, the gate voltage is slowly swept from positive to negative values. At V G intervals of 0.2 V the sweep is interrupted and Raman spectra are recorded at fixed V G values. After reaching the largest negative gate voltage, V G 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 [74,75], including a cutoff of the Coulomb interaction to implement the correct boundary conditions and gates to simulate doping [21]. We use a symmetric double-gate geometry; for monolayers a singlegate geometry leads to minimal changes on the results. The exchange-correlation functional is approximated using the Perdew-Burke-Ernzerhof (PBE) form [76] of the generalized-gradient approximation. To compute full phonon dispersions, we use the SSSP Accuracy library [77,78] (version 0.7), a 60 Ry energy cutoff, and a 32 × 32 × 1 grid with 0.02 Ry Marzari-Vanderbilt broadening [67]. We then use density-functional perturbation theory [66] -including Coulomb cutoff and gates-to compute the phonons on a 12 × 12 × 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 × 120 × 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 [79]. Phonopy [80] 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. [55][56][57] 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 Γ (phonon momentum q → 0) with respect to its value in the neutral case ∆ω = ω doped − ω neutral . The frequency of a given phonon mode corresponding to an eigenvector |ε of the dynamical matrix at Γ, D(q = 0), can be obtained as ω 2 = ε| D(q = 0)|ε / ε|ε = D Γ /M , where M = ε|ε is an effective mass associated with the phonon mode and D Γ = ε| D(q = 0)|ε = u| C(q = 0)|u is the projection of the dynamical matrix on the corresponding eigenvector or equivalently of the force constants C on the physical displacements |u of the phonon mode (with u|u = 1). Assuming that the variation of the frequency is small compared to the neutral frequency ∆ω ω neutral , we can differentiate (dω ≈ dω 2 /2ω) to obtain the following simple relation: 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. [56], we write: where F Γ is the part of the dynamical matrix at the Γ point involving the electronic response to the phonon perturbation, and thus electron-phonon interactions. The explicit expression for F Γ depends on the adiabaticity of the system. In the adiabatic limit relevant for TMDs, it reads [55][56][57]: N k k,n =m | km|∆V Γ,ν |kn | 2 f k,m − f k,n ε k,m − ε k,n (B3) + 1 N k k,n | kn|∆V Γ,ν |kn | 2 ∂f ∂ε ε k,n where ∆V q,ν is the perturbation of the effective potential due to phonon (q, ν), f k,n is the occupation of the state |kn corresponding to the band n at crystal momentum k with energy ε k,n . In the first term, there are only interband contributions (n = m) between occupied and empty states. We have verified that the numerator is much smaller than the denominator |ε k,m − ε 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: ∆ω ≈ ∆F Γ 2M ω neutral ≈ 1 N k k,n | kn|∆V Γ,ν |kn | 2 2M ω neutral ∂f ∂ε ε k,n = 1 N k k,n |g (n) ν,k,k | 2 ∂f ∂ε ε k,n , where we have introduced the square electron-phonon matrix element due to intraband transitions: 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 ε F and we obtain ∆ω ≈ − 1 N k k,n |g (n) ν,k,k | 2 δ(ε k,n − ε F ) which is Eq. (1) of the main text. Here we have introduced the density of states at the Fermi energy N (ε F ) = 1 N k k,n δ(ε k,n − ε F ) , (B8) and the square electron-phonon coupling averaged over the Fermi surface: ν,k,k | 2 δ(ε k,n − ε F ) k,n δ(ε k,n − ε F ) . (B9)