Broken adiabaticity induced by Lifshitz transition in MoS2 and WS2 single layers

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. The Born-Oppenheimer approximation is a fundamental principle of quantum mechanics and describes electronic and nuclear motion, but there are systems where it doesn’t hold true. Here, the author theoretically describes the mechanisms behind the breakdown of the Born-Oppenheimer approximation in monolayer transition metal dichalcogenides, and demonstrates that their vibrational spectra are determined by strong non-adiabatic effects when multiple valleys are occupied.

T he 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 [1][2][3] , graphite intercalation compounds 4 , and magnesium diboride 4,5 , as well as in doped semiconductors, such as boron-doped diamond 6 . 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 7 . However, to make this characterization complete, the precise understanding of the microscopic processes underlying the NA modifications of phonons is needed 1,4,8,9 . 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 [10][11][12] . In other words, the gate-induced phonon renormalization can be qualitatively described by means of adiabatic density functional theory (DFT) calculations 11 . Despite this common belief, there are few worthwhile indications that point to the importance of the NA effects in doped TMDs [12][13][14][15] . For instance, recent study on carrierinduced 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 12 . Moreover, the theoretical considerations have predicted the importance of the NA coupling in charge transfer dynamics of TMDs heterostructures 13,14 as well as in exciton-mediated Raman scattering of MoS 2 15 . Additionally, the effect beyond adiabatic approximation, namely, the phonon-mediated superconductive state was shown to exist in single-and few-layer MoS 2 when multiple conduction valleys are partially occupied [16][17][18] . Considering the plethora of exceptional charge-induced optical properties [19][20][21][22] and good carrier mobility in doped single-layer TMDs 23 , 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 MoS 2 and WS 2 , 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) 18 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 9 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 1,4,6 .

Results
Topology of valence and conduction bands. According to the experimental studies [24][25][26][27][28] , the electronic structure of the singlelayer 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 Σ point, i.e., halfway along the Γ − K path [see Fig 1a]. The valence band has band maxima in the Γ 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 Σ and Γ 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 Γ valleys in the valence band Δε v KΓ as well as between the bottom of the Σ and K valleys in the conduction band Δε c ΣK for MoS 2 and WS 2 and several electron exchange and correlation approximations. It turns out that the standard DFT calculations with semi-local Perdew, Burke and Ernzerhof (PBE) approximation are not enough to capture the right topology of the valence band in MoS 2 , since the order of the two valleys is reversed in this case (see also Supplementary Figs. 1-3). Therefore, I utilize higher order of approximation for the electronic ground state. Namely, the self-consistent GW0 29 as well as hybrid functionals PBE0 30 and vdW-DF-cx0 31 . 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 nonlocal electron-electron interactions 32 . 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 12,33,34 , 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 MoS 2 and WS 2 (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 Fig. 1a-d and e-h for MoS 2 and WS 2 , 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 Fig. 4), in close agreement with the experiments 18,33 . The abrupt appearance of the Γ and Σ valleys (e.g., for MoS 2 when charge density is around ±0.08 e/unit cell) represents the standard case of the Lifshitz transition, which is believed to be responsible for the superconductivity in MoS 2 18 . 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 MoS 2 and WS 2 single layers are shown in Fig. 1i-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 6,8,9 (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 11 . For instance, in both TMDs when carrier density is~10 14 cm −2 the difference between the NA and A frequencies (Δω NA ) of the A 1g mode is~30 cm −1 , while only~5 cm −1 for the E 2g mode. Note here that the relative NA frequency correction for the A 1g mode (Δω NA /ω A~8 %) is even larger than for the well-studied case of the E 2g mode in graphene (Δω NA /ω A~3 %) 1 for the similar amount of carrier density. The potential importance of the NA effects in pristine TMDs was in fact discussed in Ortenzi et al. 34 , 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 12 , particularly for WS 2 . A more detailed comparison with the experimental results is depicted in the Supplementary Fig. 5. Moreover, the exact carrier density windows defining the two regimes differ in MoS 2 and WS 2 , 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 35 ), and thus the framework beyond the standard Migdal-Eliashberg theory 36 might be necessary to explain the corresponding superconductive pairing as well as the discrepancy between the theoretically and experimentally obtained transition temperatures T c 18 . Energy differences between the top of the K and Γ valleys in the valence band Δε v KΓ and between the bottom of the Σ and K valleys in the conduction band Δε c ΣK for MoS2 and WS2 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. [24][25][26][27][28] are shown. All the values are given in meV. a Miwa et al. 24 and Bussolotti et al. 25 . b Eknapakul et al. 26 . c Dendzik 27 and Kastl et al. 28 . Note that the significant doping-induced redshift of the adiabatic phonon frequencies in TMDs was as well obtained in Sohier et al. 12 , where such behavior was explained in terms of reduced electrostatic screening caused by the population of the multiple valleys (see the Supplementary Figs. 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 11 . 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. 1l], 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 4,9 on phonons in TMDs. In Fig. 2a-d the results of the phonon dispersions for different dopings and for WS 2 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=τ / λ ω h i 9 (see the Methods section), the momentum-and mode-resolved EPC strengths weighted with frequencies λ qν ω qν 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., K → K with q % Γ) are present. On the other hand, when multiple valleys are partially occupied or emptied both intra-and inter-valley EHP scatterings occur (e.g., K ↔ Σ with q % M for electron-doped and K ↔ Γ with q % K for hole-doped cases) 16,18 . In Fig. 2e-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 λ ω h i as a function of doping (e.g., through sudden increase for 0.08 e ∕ unit cell in MoS 2 and for −0.08 e/unit cell in WS 2 ). 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 37,38 (see also the Supplementary Fig. 8).
Finally, in Fig. 2g-h and i-j, I show the effect of these EPCinduced EHP scatterings on the E 2g and A 1g phonon frequencies as a function of doping for MoS 2 and WS 2 , respectively. The EHP scattering processes slightly washe out the NA corrections 4,9 , improving, for example, the agreement with the experiment in the case of WS 2 (see also the Supplementary Fig. 5). However, a more striking and direct ramification of the EHP scattering events is the appearance of the phonon linewidth 4,9 (see Fig. 2k, l). In close agreement with the experiments 11,39 , 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 11,39 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 MgB 2 where nonadiabatic effects as well play a significant role 9 . 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 In addition, size of the orange circles shows the strength of the momentum-and mode-resolved electron-phonon coupling constants weighted with the corresponding frequencies λ qν ω qν . e The Eliashberg function α 2 F(ω) for 0.08 e/unit cell and for MoS 2 and WS 2 when spin-orbit coupling is included (purple and green) or not (light purple and light green). f The total amount of the electron-hole pair scattering rate due to electron-phonon coupling for several dopings and for MoS 2 (purple) and WS 2 (green). The light shades of purple and green are the results without the spin-orbit coupling. The dopings ±0.02e, ±0.04e, and ±0.08e correspond to carrier densities of ±2.34 × 10 13 cm −2 , ±4.69 × 10 13 cm −2 , and ±9.38 × 10 13 cm −2 , respectively. g, h Doping-induced frequency shifts of the E 2g and A 1g phonons in MoS 2 as obtained with the pure nonadiabatic (NA, orchid squares) theory and with the NA theory corrected with the electron-hole scattering processes (NA+τ, green squares). The experimental results from ref. 12 are shown with blue circles. i, j Same as g, h but for WS 2 . k, l show the accompanying full-width-at-half-maximum (i.e., phonon linewidth) due the electron-hole pair scattering events for MoS 2 and WS 2 , respectively. 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 [10][11][12] .
The adiabatic and nonadiabatic regimes of single-layer transition metal dichalcogenides in the case of electron doping are schematically depicted in Fig. 3a. 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 Fig. 3b, 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. 3a). On the other hand, at the outset of the Lifshitz transition (nonadiabatic regime), when the bottom of the Σ valley is only slightly occupied, the atom displacements along the A 1g mode give rise to charge modifications only for the adiabatic approximation (cf. Fig. 3d, e). In that case, part of the charge from the K valley is distributed into the Σ valley, where the electron-phonon coupling is much stronger (see the Supplementary Fig. 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 [1][2][3][4][5][6] .
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 twodimensional metallic (e.g., Nb 2 C 40 ) and semiconducting (e.g., black phosphorus 41 ) 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 6,8,9 , where q and ν are the phonon momentum and band index, respectively, ω A qν is the adiabatic phonon frequency, and e π ν is the nonadiabatic phonon self-energy due to the EPC (see the Supplementary Note 1). The real part of e π ν gives the renormalization of the phonon frequency due to the NA coupling, i.e., qν Re e π ν ðq; ωÞ, while the imaginary part is the NA phonon linewidth, i.e., γ qν ¼ À2Im e π ν ðq; ωÞ. For simulating the Raman spectra only the q ≈ 0 limits of e π ν and B ν 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 1,4,8,9 , The electron band index and momentum are μ and k, respectively, ε μk is the corresponding electron energy, f(ε μk ) is the Fermi-Dirac distribution function, while g b μμ;ν and g μμ,ν are the bare and screened intraband electron-phonon coupling function. The dynamical interband contribution at q ≈ 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 8 , while the usual assumption is that both vertex functions are screened, i.e., The results of Eq. (2) for WS 2 and different dopings obtained with g Ã μμ 0 ;ν g b μμ 0 ;ν and g μμ 0 ;ν 2 are shown in Supplementary Fig. 9. The corresponding difference is small, as it was already assumed in previous studies 4,42 .
When EPC-induced electron-hole pair scattering processes up to all orders are taken into account 4,9,43 , the intraband phonon self-energy acquires the following form 9 , 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 9,44 , where k B is the Boltzmann constant and α 2 F(Ω) is the Eliashberg function. The dynamical energy renormalization parameter λ ep (ω) is obtained by performing the Kramers-Kronig transformation of 1 ∕ τ ep (ω). For T = 0 K and for frequencies much higher than the characteristic frequencies of the system, the damping rate acquires the following form 1=τ ep ð1Þ ¼ π P qν λ qν ω qν ¼ πλ ω h i, where λ is the EPC constant and ω h i is the so-called first moment of the phonon spectrum 45 .
Computational details. The ground-state calculations were done by means of the QUANTUM ESPRESSO (QE) package 46 with a plane-wave cutoff energy of 80 Ry. Fully relativistic norm-conserving pseudopotentials from the PseudoDojo project 47 were used with the PBE exchange-correlation functional 48 . Spin-orbit coupling was also included. A (24 × 24 × 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 16 .
In order to investigate the effects of the non-local electron-electron interactions on the electronic structure of the TMDs, the hybrid functionals PBE0 30 and vdW-DF-cx0 31 were employed within QE. In addition, the corresponding electronic structures of MoS 2 and WS 2 were calculated within the self-consistent GW0 29 (i.e., where eigenvalues in the electronic Green's function are included self-consistently) implemented in GPAW 49 . Moderate differences between non-self-consistent G0W0 and self-consistent GW0 results are shown in the Supplementary Fig. 1. The GPAW calculations were done using the PAW pseudopotentials with energy cutoff of 600 eV, PBE functional, and a (18 × 18 × 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 Figs. 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 spinorbit 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 33,34 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 MoS 2 the unit cell strain of ϵ = −1.59% (a = 3.135 Å) and ϵ = 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 WS 2 the unit cell strain of ϵ = −1.64% (a = 3.14 Å) and ϵ = −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 Figs. 1-4. Adiabatic phonon frequencies at the Γ point were calculated on a dense (72 × 72 × 1) grid using density functional perturbation theory 50 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. 51 ). The Eliashberg functions needed for calculating the relaxation rate Eq. (4) were calculated with the EPW code 52 , while the corresponding input quantities (electron energies, phonon frequencies, and electron-phonon coupling elements) were interpolated using maximally localized Wannier functions 53 . The corresponding summations were done on (480 × 480 × 1) electron and (48 × 48 × 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 12 . The phonon calculations were performed without and with 2D Coulomb interaction truncation, which eliminates the spurious interaction between periodic images (i.e., layers) 54,55 . However, no difference between the two types of calculations were obtained (see the Supplementary Fig. 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.