Broken adiabaticity induced by Lifshitz transition in MoS$_2$ and WS$_2$ single layers

TL;DR

This paper reveals how Lifshitz transitions in doped single-layer MoS2 and WS2 cause breakdowns in adiabaticity, affecting phonon modes and electron-phonon interactions, with implications for Raman spectroscopy in 2D materials.

Contribution

It demonstrates the link between Lifshitz transitions and adiabaticity breakdown in doped 2D transition metal dichalcogenides, highlighting the roles of electron-electron and electron-phonon interactions.

Findings

Nonadiabatic coupling affects Raman-active phonon modes.

Lifshitz transition triggers adiabaticity breakdown.

Doping modifies phonon linewidths via electron-hole scattering.

Abstract

The breakdown of the adiabatic Born-Oppenheimer approximation is striking dynamical phenomenon, however, it occurs only in a handful of layered materials. Here, I show that adiabaticity breaks down in doped single-layer transition metal dichalcogenides in a quite intriguing manner. Namely, significant nonadiabatic coupling, which acts on frequencies of the Raman-active modes, is prompted by a Lifshitz transition due to depopulation and population of multiple valence and conduction valleys, respectively. The outset of the latter event is shown to be dictated by the interplay of highly non-local electron-electron interaction and spin-orbit coupling. In addition, intense electron-hole pair scatterings due to electron-phonon coupling are inducing phonon linewidth modifications as a function of doping. Comprehending these intricate dynamical effects turns out to be a key for mastering characterization of electron doping in two-dimensional nano-devices by means of Raman spectroscopy.

Tables

Table 1. Table 1: Topology details of valence and conduction valleys. Energy differences between the top of the K and Γ Γ \Gamma valleys in the valence band Δ ​ ε K ​ Γ v Δ subscript superscript 𝜀 𝑣 K Γ \Delta\varepsilon^{v}_{{\rm K}\Gamma} and between the bottom of the Σ Σ \Sigma and K valleys in the conduction band Δ ​ ε Σ ​ K c Δ subscript superscript 𝜀 𝑐 Σ K \Delta\varepsilon^{c}_{\Sigma{\rm K}} for MoS 2 and WS 2 as obtained by means of several electron exchange and correlation approximations without and with spin-orbit coupling (SOC). In addition, the experimental results obtained in Refs.  (Miwa  et al. , 2015 ; Bussolotti  et al. , 2019 ; Eknapakul  et al. , 2014 ; Dendzik  et al. , 2015 ; Kastl  et al. , 2018 ) are shown. All the values are given in meV.

		xc approximation
		PBE	PBE+SOC	GW0	GW0+SOC	PBE0	PBE0+SOC	DF-cx0	DF-cx0+SOC	exp.
MoS2	$\Delta {\epsilon }_{K\Gamma }^v$	$-$19	51	194	264	92	171	203	273	140-310111Refs. Miwa et al. (2015); Bussolotti et al. (2019)
	$\Delta {\epsilon }_{\Sigma K}^c$	275	244	54	23	431	427	265	234	$\gtrsim 60$222Ref. Eknapakul et al. (2014)
WS2	$\Delta {\epsilon }_{K\Gamma }^v$	46	231	164	349	103	356	232	417	240-510333Refs. Dendzik et al. (2015); Kastl et al. (2018)
	$\Delta {\epsilon }_{\Sigma K}^c$	190	84	112	6	338	209	137	31	–

Equations

$$B_{\nu }(\mathbf{q},\omega )=-\frac{1}{\pi }\mathrm{Im}\left[\frac{2{\omega }_{\mathbf{q}\nu }^{\mathrm{A}}}{{\omega }^2-({\omega }_{\mathbf{q}\nu }^{\mathrm{A}}{)}^2-2{\omega }_{\mathbf{q}\nu }^{\mathrm{A}}{\tilde{\pi }}_{\nu }(\mathbf{q},\omega )}\right],$$

$${\tilde{\pi }}_{\nu }^{0,\mathrm{intra}}(\omega )=\sum _{\mu \mathbf{k}}{g}_{\mu \mu ,\nu }^b(\mathbf{k},0)g_{\mu \mu ,\nu }^{\ast }(\mathbf{k},0)\left[-\frac{\partial f({\epsilon }_{\mu \mathbf{k}})}{\partial {\epsilon }_{\mu \mathbf{k}}}\right].$$

$${\tilde{\pi }}_{\nu }^{\mathrm{intra}}(\omega )$$

$$1/{\tau }_{\mathrm{ep}}(\omega )$$

Full text

The Born-Oppenheimer approximation, which assumes that electrons adiabatically follow the motion of lattice degrees of freedom, is a fundamental starting point in quantum mechanics. Recent studies have, however, shown that striking nonadiabatic (NA) effects visible in vibrational Raman spectra can be found in metallic layered materials, such as graphene Lazzeri and Mauri (2006); Pisana et al. (2007); Ferrante et al. (2018), graphite intercalation compounds Saitta et al. (2008), and magnesium diboride Quilty et al. (2002); Saitta et al. (2008), as well as in doped semiconductors, such as boron-doped diamond Caruso et al. (2017). The strength of this dynamical electron-phonon coupling (EPC) can, in fact, be modified as a function of carrier concentration, which makes the Raman spectroscopy a quite powerful tool for characterization of doped two-dimensional and layered materials Das et al. (2008). However, to make this characterization complete, the precise understanding of the microscopic processes underlying the NA modifications of phonons is needed Lazzeri and Mauri (2006); Saitta et al. (2008); Giustino (2017); Novko (2018). Such renormalization of phonons is actually considered quite rare, and apart from these few cases it is not clear in what other two-dimensional materials should adiabaticity break down.

Contrary to the aforesaid examples, the vibrational spectra of the single-layer transition metal dichalcogenides (TMDs) under bias voltage is believed to be adiabatic Zhang et al. (2015); Chakraborty et al. (2012); Sohier et al. (2019). In other words, the gate-induced phonon renormalization can be qualitatively described by means of adiabatic density functional theory (DFT) calculations Chakraborty et al. (2012). Despite this common belief, there are few worthwhile indications that point to the importance of the NA effects in doped TMDs Sohier et al. (2019); Zheng et al. (2017, 2018); Reichardt and Wirtz (2019). For instance, recent study on carrier-induced phonon modifications in atomically-thin TMDs had shown that the adiabatic DFT results are largely overestimating the corresponding frequency shifts of the Raman-active modes, but NA corrections were not taken into account Sohier et al. (2019). Moreover, the theoretical considerations have predicted the importance of the NA coupling in charge transfer dynamics of TMDs heterostructures Zheng et al. (2017, 2018) as well as in exciton-mediated Raman scattering of MoS2 Reichardt and Wirtz (2019). Additionally, the effect beyond adiabatic approximation, namely, the phonon-mediated superconductive state was shown to exist in single- and few-layer MoS2 when multiple conduction valleys are partially occupied Ge and Liu (2013); Costanzo et al. (2016); Piatti et al. (2018). Considering the plethora of exceptional charge-induced optical properties Mak et al. (2012); Jo et al. (2014); Wang et al. (2015, 2017) and good carrier mobility in doped single-layer TMDs Radisavljevic et al. (2011), as well as the pivotal role of the EPC in these features, it is of paramount importance to decipher the role of the NA effects and the EPC in the corresponding Raman spectra.

Here, I show that the vibrational spectra of single-layer MoS2 and WS2, as prototypical examples of semiconducting TMDs, doped with electrons and holes are indeed governed by strong NA effects. Namely, when more than one valley of either valence or conduction electronic bands are simultaneously partially occupied or empty a Lifshitz transition occurs (i.e., an abrupt change of the Fermi surface) Piatti et al. (2018) which then activates the significant dynamical part of the EPC. The latter is manifested as a considerable redshift of the $A_{1g}$ optical phonon frequency along with the large corresponding NA correction. In addition, there are other intriguing coupling mechanisms ruling the NA renormalization of phonons in TMDs. Firstly, the energy difference between the top (bottom) of the valence (conduction) band valleys, and thus the amount of carriers needed in order to empty (fill) both of them, is highly sensitive to the non-local electron-electron interaction. Furthermore, the spin-orbit coupling (SOC) is important in order to correctly capture the EPC in TMDs. And lastly, a high degree of the electron-hole pair scattering Novko (2018) is present in doped TMDs, which in turn slightly washes out the NA renormalization effects and is responsible for the enlargement of the phonon linewidths as a function of doping. All these processes are essential for recreating the experimental observations and therefore the breakdown of the Born-Oppenheimer approximation in doped TMDs turns out to be even more intricate than in the well-know examples Lazzeri and Mauri (2006); Saitta et al. (2008); Caruso et al. (2017).

Results

Topology of valence and conduction bands. According to the experimental studies Miwa et al. (2015); Bussolotti et al. (2019); Eknapakul et al. (2014); Dendzik et al. (2015); Kastl et al. (2018), the electronic structure of the single-layer semiconducting TMDs is characterized with the two valleys in both the conduction and valence bands. For the former, the global minimum is located in the K point of the Brillouin zone, while the remnant minimum is at the $\Sigma$ point, i.e., halfway along the $\Gamma-{\rm K}$ path [see Fig 1(a)]. The valence band has band maxima in the $\Gamma$ and K points, where the latter is the global maximum. The corresponding measured energy differences between these minima and maxima are presented in Table 1. Since the exact carrier density required for the Fermi level to cross the second valleys centered at the $\Sigma$ and $\Gamma$ points depends on this energy difference, it is first essential to capture the proper electronic band structures of single-layer TMDs before studying in detail the doping-induced phonon frequency shifts.

Table 1 shows the results for the energy difference between the top of the K and $\Gamma$ valleys in the valence band $\Delta\varepsilon^{v}_{{\rm K}\Gamma}$ as well as between the bottom of the $\Sigma$ and K valleys in the conduction band $\Delta\varepsilon^{c}_{\Sigma{\rm K}}$ for MoS2 and WS2 and several electron exchange and correlation approximations. It turns out that the standard density functional theory (DFT) calculations with semi-local Perdew, Burke and Ernzerhof (PBE) approximation are not enough to capture the right topology of the valence band in MoS2, since the order of the two valleys is reversed in this case (see also Supplementary Figures 1-3). Therefore, I utilize higher order of approximation for the electronic ground state. Namely, the self-consistent GW0 Hüser et al. (2013) as well as hybrid functionals PBE0 Perdew et al. (1996a) and vdW-DF-cx0 Berland et al. (2017). The corresponding results are, on the other hand, much closer to the experiments, which shows how the topology of the valleys in single-layer TMDs, and not just the energy band gap, is highly dictated by the strong quasi-particle physics, i.e., by highly non-local electron-electron interactions Kormányos et al. (2015). The SOC as well plays an important role here. Unfortunately, the DFT phonon calculations are numerically unfeasible when performed on top of the GW or the hybrid-functional electronic structure. In order to overcome this hurdle and having in mind that the strain of the unit cell can induce the valley modifications Yuan et al. (2016); Sohier et al. (2019); Ortenzi et al. (2018), the remaining calculations are performed with the PBE functional and with the exact amount of the strain necessary to recreate the correct order and topology of the valleys in MoS2 and WS2 (i.e., as obtained with the vdW-DF-cx0 functional). Further computational details can be found in the Methods section.

The corresponding Fermi surface cuts as a function of electron and hole dopings are shown in Figs. 1(a)-(d) and 1(e)-(h) for MoS2 and WS2, respectively. These Fermi surfaces clearly depict the transformation from the one-valley to multi-valley electronic band structure as the electron and hole densities are elevated (see also Supplementary Figure 4), in close agreement with the experiments Yuan et al. (2016); Piatti et al. (2018). The abrupt appearance of the $\Gamma$ and $\Sigma$ valleys (e.g., for MoS2 when charge density is around $\pm 0.08\,e/$unit cell) represents the standard case of the Lifshitz transition, which is believed to be responsible for the superconductivity in MoS2 Piatti et al. (2018). Here I show that this sudden change in the electron density of states induces the significant dynamical effects in the EPC contributing to the frequency of the Raman-active phonon modes.

Nonadiabatic frequency renormalization. The calculated doping-induced frequency shifts of the $E_{2g}$ and $A_{1g}$ optical phonon modes in MoS2 and WS2 single layers are shown in Figs. 1(i)-(l). In particular, I show the results obtained by means of the adiabatic (A) density functional perturbation theory (DFPT) as well as the frequencies corrected with the DFT-based NA method Novko (2018); Caruso et al. (2017); Giustino (2017) (see the Method section and Supplementary Note 1). The corresponding absolute values of the adiabatic frequencies for pristine systems are given in the Supplementary Table 1. From this the two distinct regimes emerge in the doped single-layer TMDs. Namely, the A regime where only the K point valleys intersect the Fermi level and the NA regime where the Lifshitz transition occurs and sizable dynamical contribution to the EPC is triggered. This intriguing phenomenon is present in both systems and for the both optical phonon modes. However, the NA effects are more pronounced for the $A_{1g}$ phonons, indicating a larger EPC in this case Chakraborty et al. (2012). For instance, in both TMDs when carrier density is $\sim 10^{14}$ cm*-2* the difference between the NA and A frequencies ($\Delta\omega_{\rm NA}$) of the $A_{1g}$ mode is $\sim 30$ cm*-1*, while only $\sim 5$ cm*-1* for the $E_{2g}$ mode. Note here that the relative NA frequency correction for the $A_{1g}$ mode ($\Delta\omega_{\rm NA}/\omega_{\rm A}\sim 8\%$) is even larger than for the well-studied case of the $E_{2g}$ mode in graphene ($\Delta\omega_{\rm NA}/\omega_{\rm A}\sim 3\%$) Lazzeri and Mauri (2006) for the similar amount of carrier density. The potential importance of the NA effects in pristine TMDs was in fact discussed in L. Ortenzi et al. Ortenzi et al. (2018), where a strong entanglement between the lattice degrees of freedom and the electronic excitations was found by using electronic band structure calculations. However, the NA corrections to the phonon spectrum were not discussed nor calculated. The results also show how the NA corrections are crucial for achieving a quantitative agreement with the experimental frequency shifts Sohier et al. (2019), particularly for WS2. A more detailed comparison with the experimental results is depicted in the Supplementary Figure 5. Moreover, the exact carrier density windows defining the two regimes differ in MoS2 and WS2, which is due to different topology of the valence and conduction bands in these two TMDs. Consequently, the latter demonstrates how the Raman spectroscopy could be utilized in order to extract the right amount of carrier density required for one- to multi-valley transformations in TMDs and similar semiconductors. Also, since both superconductive state and the breakdown of the adiabaticity occur at the outset of the Lifshitz transition, it is likely that the nonadiabatic channels are open in the superconductivity of the TMDs (Migdal’s theorem is not valid in the nonadiabatic regime Engelsberg and Schrieffer (1963)), and thus the framework beyond the standard Migdal-Eliashberg theory Grimaldi et al. (1995) might be necessary to explain the corresponding superconductive pairing as well as the discrepancy between the theoretically- and experimentally-obtained transition temperatures $T_{c}$ Piatti et al. (2018).

Note that the significant doping-induced redshift of the adiabatic phonon frequencies in TMDs was as well obtained in T. Sohier et al. Sohier et al. (2019), where such behaviour was explained in terms of reduced electrostatic screening caused by the population of the multiple valleys (see the Supplementary Figures 6 and 7 for comparison and further discussion). However, it should be emphasized that the NA effects were disregarded in that case and thus the frequency shifts were considered as a pure adiabatic phenomenon. Similar reasoning was as well delivered in the earlier work Chakraborty et al. (2012). On the other hand, here I show that the significant redshift of the adiabatic $A_{1g}$ phonon frequency is reduced due to strong dynamical effects [e.g., see Fig. 1(l)], which in turn improves the agreement with the experiment.

Electron-hole pair scattering processes. In what follows, I investigate the impact of the EPC-induced electron-hole pair (EHP) scattering processes Saitta et al. (2008); Novko (2018) on phonons in TMDs. In Figs. 2(a)-(d) the results of the phonon dispersions for different dopings and for WS2 are shown along high-symmetry points of the Brillouin zone. Since the EHP rate is proportional to the EPC strength weighted with the first moment of the phonon spectrum, i.e., $1/\tau\propto\lambda\left\langle\omega\right\rangle$ Novko (2018) (see the Methods section), the momentum- and mode-resolved EPC strengths weighted with frequencies $\lambda_{\mathbf{q}\nu}\omega_{\mathbf{q}\nu}$ are shown as well (orange circles). The doping induces changes in the EHP scattering rates in accordance with the transformations of the Fermi surface cuts. In particular, when solely K valleys intersect the Fermi level only the intra-valley EHP scatterings (i.e., $\mathrm{K}\rightarrow\mathrm{K}$ with $\mathbf{q}\approx\overline{\Gamma}$) are present. On the other hand, when multiple valleys are partially occupied or emptied both intra- and inter-valley EHP scatterings occur (e.g., $\mathrm{K}\leftrightarrow\Sigma$ with $\mathbf{q}\approx\overline{\rm M}$ for electron-doped and $\mathrm{K}\leftrightarrow\Gamma$ with $\mathbf{q}\approx\overline{\rm K}$ for hole-doped cases) Ge and Liu (2013); Piatti et al. (2018). In Figs. 2(e)-(f) the results show how the EPC-induced EHP scattering rate as a function of doping depends greatly on the topology of the bands as well as on the SOC. For instance, the fingerprints of the Lifshitz transition are distinctly visible in the change of the $\lambda\left\langle\omega\right\rangle$ as a function of doping (e.g., through sudden increase for $0.08\,e/$unit cell in MoS2 and for $-0.08\,e/$unit cell in WS2). Furthermore, the SOC reduces the EHP scattering rates by more than factor of two in some cases. This is due to spin-orbit splitting of the valleys, where, e.g., the lower valence valley is more coupled to the phonons than the upper valence valley Hinsche et al. (2017); Mahatha et al. (2019) (see also the Supplementary Figure 8).

Finally, in Figs. 2(g)-(h) and 2(i)-(j), I show the effect of these EPC-induced EHP scatterings on the $E_{2g}$ and $A_{1g}$ phonon frequencies as a function of doping for MoS2 and WS2, respectively. The EHP scattering processes slightly washe out the NA corrections Saitta et al. (2008); Novko (2018), improving, for example, the agreement with the experiment in the case of WS2 (see also the Supplementary Figure 5). However, a more striking and direct ramification of the EHP scattering events is the appearance of the phonon linewidth Saitta et al. (2008); Novko (2018) [see Figs. 2(k) and 2(l)]. In close agreement with the experiments Chakraborty et al. (2012); Lu et al. (2017), the linewidth of the $A_{1g}$ mode shows a steeper increase as a function of doping, and is generally larger than the $E_{2g}$ linewidth. This proves that the major part of the experimentally observed doping-induced phonon linewidth modifications in TMDs Chakraborty et al. (2012); Lu et al. (2017) is underlain by the EHP scattering processes due to the EPC (other contributions might be anharmonic coupling and scattering on inhomogeneities), as it is, for example, the case in MgB2 where nonadiabatic effects as well play a significant role Novko (2018). The linewidth, similarly to the frequency shift, shows a great sensitivity on the topology of the valleys, and therefore can also be used in Raman spectroscopy as an indicator of the Lifshitz transition.

Discussion

All in all, the phonon dynamics in doped single-layer transition metal dichalcogenides visible in Raman spectra was shown to be quite complex and governed by the interplay between electron-electron, spin-orbit, and dynamical electron-phonon interactions. In particular, at the outset of the Lifshitz transition, i.e., when the Fermi level crosses multiple conduction or valence valleys, the adiabatic Born-Oppenheminer approximation breaks down and significant nonadiabatic effects are prompted. Further, the outset of the Lifshitz transition was shown to be dictated by the non-local electron-electron interaction, which determines the right topology of the valleys. In addition to that, the strong electron-hole pair scatterings due to electron-phonon coupling are present and responsible for the phonon linewidth increase as a function of doping. Note that these results are at odds with the previous studies, where the phonon spectrum of transition metal dichalcogenides was considered to be adiabatic Zhang et al. (2015); Chakraborty et al. (2012); Sohier et al. (2019).

The adiabatic and nonadiabatic regimes of single-layer transition metal dichalcogenides in the case of electron doping are schematically depicted in Fig. 3(a). When only the K valley is partially occupied (the adiabatic regime), the modifications of the electronic structure induced by the $A_{1g}$ phonon does not generate the changes in the charge density, as shown in Figs. 3(b) and 3(c). Consequently, in this adiabatic regime, both adiabatic and nonadibatic approximations give the same results, i.e., zero frequency shifts [see the inset in Fig. 3(a)]. On the other hand, at the outset of the Lifshitz transition (nonadiabatic regime), when the bottom of the $\Sigma$ valley is only slightly occupied, the atom displacements along the $A_{1g}$ mode give rise to charge modifications only for the adiabatic approximation [cf. Figs. 3(d) and 3(e)]. In that case, part of the charge from the K valley is distributed into the $\Sigma$ valley, where the electron-phonon coupling is much stronger (see the Supplementary Figure 8). For the nonadiabatic approximation, carriers remain in the same valley. As a result, the corresponding phonon frequency shifts turn out to be large for the adiabatic, while small for the nonadiabatic approximation. It is important to stress that such multi-valley nonadiabatic mechanism was not discussed thus far Lazzeri and Mauri (2006); Pisana et al. (2007); Ferrante et al. (2018); Saitta et al. (2008); Quilty et al. (2002); Caruso et al. (2017).

The insights given here might help elucidate the superconductive pairing mechanism in transition metal dichalcogenides, as well as puzzling discrepancy between the theoretically- and experimentally-obtained transition temperatures $T_{c}$. Even more, it is possible that the nonadiabatic channels investigated here could be universal and important in other novel multi-valley two-dimensional metallic (e.g., Nb2C Huang et al. (2019)) and semiconducting (e.g., black phosphorus Chakraborty et al. (2016)) materials in the presence of the external excess charge.

Methods

Phonon self-energy. Using the many-body perturbation theory, the phonon spectral function can be obtained with the following formula Novko (2018); Giustino (2017); Caruso et al. (2017),

[TABLE]

where $\mathbf{q}$ and $\nu$ are the phonon momentum and band index, respectively, $\omega_{\mathbf{q}\nu}^{\mathrm{A}}$ is the adiabatic phonon frequency, and $\widetilde{\pi}_{\nu}$ is the nonadiabatic phonon self-energy due to the EPC (see the Supplementary Note 1). The real part of $\widetilde{\pi}_{\nu}$ gives the renormalization of the phonon frequency due to the NA coupling, i.e., $\omega^{2}=(\omega_{\mathbf{q}\nu}^{\mathrm{A}})^{2}+2\omega_{\mathbf{q}\nu}^{\mathrm{A}}\mathrm{Re}\,\widetilde{\pi}_{\nu}(\mathbf{q},\omega)$, while the imaginary part is the NA phonon linewidth, i.e., $\gamma_{\mathbf{q}\nu}=-2\mathrm{Im}\,\widetilde{\pi}_{\nu}(\mathbf{q},\omega)$. For simulating the Raman spectra only the $\mathbf{q}\approx 0$ limits of $\widetilde{\pi}_{\nu}$ and $B_{\nu}$ are relevant.

In the absence of the electron-hole pair scatterings the dynamical correction over the adiabatic phonon spectral function comes from the dynamical bare intraband phonon self-energy Lazzeri and Mauri (2006); Saitta et al. (2008); Giustino (2017); Novko (2018),

[TABLE]

The electron band index and momentum are $\mu$ and $\mathbf{k}$, respectively, $\varepsilon_{\mu\mathbf{k}}$ is the corresponding electron energy, $f(\varepsilon_{\mu\mathbf{k}})$ is the Fermi-Dirac distribution function, while $g^{b}_{\mu\mu,\nu}$ and $g_{\mu\mu,\nu}$ are the bare and screened intraband electron-phonon coupling function. The dynamical interband contribution at $\mathbf{q}\approx 0$ is negligible for the TMDs, since the interband gap is much larger than the relevant phonon frequencies (see also the Supplementary Table 2). Equation 2 only contributes to the renormalization of the phonon frequency. Note that the correct treatment of the Coulomb screening in the derivation of the phonon self-energy is employed here, i.e., with one screened and one bare electron-phonon coupling function Giustino (2017), while the usual assumption is that both vertex functions are screened, i.e., $g_{\mu\mu^{\prime},\nu}^{*}g_{\mu\mu^{\prime},\nu}^{b}\rightarrow\left|g_{\mu\mu^{\prime},\nu}\right|^{2}$. The results of Eq. (2) for WS2 and different dopings obtained with $g_{\mu\mu^{\prime},\nu}^{*}g_{\mu\mu^{\prime},\nu}^{b}$ and $\left|g_{\mu\mu^{\prime},\nu}\right|^{2}$ are shown in Supplementary Figure 9. The corresponding difference is small, as it was already assumed in previous studies Saitta et al. (2008); Calandra et al. (2010).

When EPC-induced electron-hole pair scattering processes up to all orders are taken into account Novko (2018); Cappelluti (2006); Saitta et al. (2008), the intraband phonon self-energy acquires the following form Novko (2018),

[TABLE]

The scattering parameter describes the damping rate of the excited electron-hole pairs and for the case of electron-phonon scattering can be written as Novko (2018); Shulga et al. (1991),

[TABLE]

where $k_{B}$ is the Boltzmann constant and $\alpha^{2}F(\Omega)$ is the Eliashberg function. The dynamical energy renormalization parameter $\lambda_{\rm ep}(\omega)$ is obtained by performing the Kramers-Kronig transformation of $1/\tau_{\rm ep}(\omega)$. For $T=0$ K and for frequencies much higher than the characteristic frequencies of the system, the damping rate acquires the following form $1/\tau_{\rm ep}(\infty)=\pi\sum_{\mathbf{q}\nu}\lambda_{\mathbf{q}\nu}\omega_{\mathbf{q}\nu}=\pi\lambda\left\langle\omega\right\rangle$, where $\lambda$ is the EPC constant and $\left\langle\omega\right\rangle$ is the so-called first moment of the phonon spectrum Novko (2017).

Computational details. The ground-state calculations were done by means of the quantum espresso (QE) package Giannozzi et al. (2017) with a plane-wave cutoff energy of 80 Ry. Fully relativistic norm-conserving pseudopotentials from the PseudoDojo project van Setten et al. (2018) were used with the PBE exchange-correlation functional Perdew et al. (1996b). Spin-orbit coupling was also included. A ($24\times 24\times 1$) Monkhorst-Pack grid was used for sampling the Brillouin zone (with Gaussian smearing of 0.02 Ry). Electron and hole dopings were simulated by adding and removing, respectively, the electrons and introducing the compensating homogeneous charged background Ge and Liu (2013).

In order to investigate the effects of the non-local electron-electron interactions on the electronic structure of the TMDs, the hybrid functionals PBE0 Perdew et al. (1996a) and vdW-DF-cx0 Berland et al. (2017) were employed within QE. In addition, the corresponding electronic structures of MoS2 and WS2 were calculated within the self-consistent GW0 Hüser et al. (2013) (i.e., where eigenvalues in the electronic Green’s function are included self-consistently) implemented in gpaw Enkovaara et al. (2010). Moderate differences between non-self-consistent G0W0 and self-consistent GW0 results are shown in the Supplementary Figure 1. The gpaw calculations were done using the PAW pseudopotentials with energy cutoff of 600 eV, PBE functional, and a ($18\times 18\times 1$) electron-momentum grid. Fermi-Dirac smearing was used with $T=300$ K. The GW0 calculations were done on the highest valence and lowest conduction bands with the 100 eV plane wave cutoff. The corresponding results of the valence and conduction band topology are in Table 1. Note that the comparison between different xc approximations is made without spin-orbit coupling in Supplementary Figures 1-3. This is because the calculation with spin-orbit coupling along with hybrid DF-cx0 functional is not implemented in QE (the DF-cx0+SOC values presented in Table 1 are actually the DF-cx0 values corrected with the energy modifications due to spin-orbit coupling as obtained with the PBE). Also, the band structure modifications induced by spin-orbit coupling are very similar at the GW0 and PBE levels, and are in agreement with ARPES experiments. Therefore, it is justified to make the fitting procedure described below without spin-orbit coupling, and then include it a posteriori (at the PBE level).

Since the phonon calculations are unfeasible with the hybrid functionals as well as with the GW approaximation, the further calculations were done with the PBE functional and with the exact amount of unit cell strain Yuan et al. (2016); Ortenzi et al. (2018) required to reproduce the same topology of the valleys as obtained with the hybrid vdW-DF-cx0 functional (the conclusions of the present paper would not change at all if the reference topology of valleys is as obtained with GW0 or PBE0). In particular, in the case of MoS2 the unit cell strain of $\epsilon=-1.59\%$ ($a=3.135$ Å) and $\epsilon=0\%$ ($a=3.1856$ Å) were used, respectively, to simulate the right topology of the valence and conduction bands (which are then used for the hole and electron doping, respectively). For WS2 the unit cell strain of $\epsilon=-1.64\%$ ($a=3.14$ Å) and $\epsilon=-0.54\%$ ($a=3.175$ Å) were used, respectively, to simulate the right topology of the valence and conduction bands. The distances between the Mo/W and S atoms were relaxed for each doping until forces were less than $10^{-5}$ Ry/a.u. The spin-orbit coupling was then introduced for these strained structures and in the following phonon calculations. For further information see the Supplementary Figures 1-4.

Adiabatic phonon frequencies at the $\Gamma$ point were calculated on a dense ($72\times 72\times 1$) grid using density functional perturbation theory Baroni et al. (2001) as implemented in QE and with the Gaussian smearing corresponding to $T=300$ K. The intraband nonadiabatic phonon self-energies Eq. (2) were extracted from the DFPT calculations as described in Supplementary Note 1 (the corresponding method is based on work presented in Y. Nomura et al. Nomura and Arita (2015)). The Eliashberg functions needed for calculating the relaxation rate Eq. (4) were calculated with the epw code Poncé et al. (2016), while the corresponding input quantities (electron energies, phonon frequencies, and electron-phonon coupling elements) were interpolated using maximally localized Wannier functions Marzari et al. (2012). The corresponding summations were done on $(480\times 480\times 1)$ electron and $(48\times 48\times 1)$ phonon momentum grids. The temperature of the Fermi-Dirac distribution in Eqs. (2) and (3) as well as the temperature entering Eq. (4) was set to 300 K as in the experiment Sohier et al. (2019). The phonon calculations were performed without and with 2D Coulomb interaction truncation, which eliminates the spurious interaction between periodic images (i.e., layers) Sohier et al. (2016, 2017). However, no difference between the two types of calculations were obtained (see the Supplementary Figure 10).

Data Availability

The data that support the findings of this study are available from the author on reasonable request.

Code Availability

The codes used to produce the data presented in this study are available from the author on reasonable request.

Acknowledgments

The author gratefully acknowledges financial support from the European Regional Development Fund for the “Center of Excellence for Advanced Materials and Sensing Devices” (Grant No. KK.01.1.1.01.0001). Financial support by Donostia International Physics Center (DIPC) during various stages of this work is also highly acknowledged. Computational resources were provided by the DIPC computing center.

Author contributions

D.N. conceived the idea, developed the method, performed theoretical calculations and analysis, created visualizations, and wrote the manuscript.

Competing interests

The author declare no competing interests.

Additional information

Supplementary information is available for this paper at [URL].