# Unconventional gate-induced superconductivity in transition-metal dichalcogenides

Thibault Sohier, Marco Gibertini, Ivar Martin, Alberto F Morpurgo

► To cite this version:

Thibault Sohier, Marco Gibertini, Ivar Martin, Alberto F Morpurgo. Unconventional gate-induced superconductivity in transition-metal dichalcogenides. Physical Review Research, 2025, 7, pp.013290. 10.1103/PhysRevResearch.7.013290. hal-05370132

HAL Id: hal-05370132

https://hal.science/hal-05370132v1

Submitted on 18 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Unconventional gate-induced superconductivity in transition-metal dichalcogenides

Thibault Sohier,¹,*, Marco Gibertini,²,³ Ivar Martin,⁴,⁵ and Alberto F. Morpurgo⁶,†

¹Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, Montpellier, France

²Dipartimento di Scienze Fisiche, Informatiche e Matematiche, University of Modena and Reggio Emilia, IT-41125 Modena, Italy

³Centro S3, CNR-NANO Istituto Nanoscienze, IT-41125, Modena, Italy

⁴Materials Science Division, Argonne National Laboratory, Lemont, Illinois 60439, USA

⁵Department of Physics, University of Chicago, Chicago, Illinois 60637, USA

⁶Department of Quantum Matter Physics, University of Geneva, CH-1211 Geneva, Switzerland

<figure><img src="image_2.png" /></figure>

(Received 19 September 2024; accepted 20 February 2025; published 19 March 2025)

Superconductivity in few-layer semiconducting transition-metal dichalcogenides (TMDs) can be induced by field-effect doping through ionic-liquid gating. While several experimental observations have been collected over the years, a fully consistent theoretical picture is still missing. Here we develop a realistic framework that combines the predictive power of first-principles simulations with the versatility and insight of Bardeen-Cooper-Schrieffer gap equations to rationalize such experiments. The multivalley nature of semiconducting TMDs is taken into account, together with the doping- and momentum-dependent electron-phonon and Coulomb interactions. Consistently with experiments, we find that superconductivity occurs when the electron density is large enough that the *Q* valleys get occupied, as a result of a large enhancement of electron-phonon interactions. Despite being phonon driven, the superconducting state is predicted to be sensitive to Coulomb interactions, which can lead to the appearance of a relative sign difference between valleys and thus to a $s_{+−}$ character. We discuss qualitatively how such scenario may account for many of the observed physical phenomena for which no microscopic explanation has been found so far, including in particular the presence of a large subgap density of states, and the sample-dependent dome-shaped dependence of $T_c$ on accumulated electron density. Our results provide a comprehensive analysis of gate-induced superconductivity in semiconducting TMDs, and introduce an approach that will likely be valuable for other multivalley electronic systems, in which superconductivity occurs at relatively low electron density.

DOI: 10.1103/PhysRevResearch.7.013290

## I. INTRODUCTION

Gate-induced superconductivity at the surface of semiconducting transition-metal dichalcogenides (TMDs) has been discovered 10 years ago, in pioneering experiments on ionic gated transistors based on exfoliated $MoS_2$ crystals [1,2]. Subsequent work [3–6] showed that superconductivity occurs in the majority of the commonly studied semiconducting TMDs ($WS_2$, $MoSe_2$), starting at different, but comparable, values of accumulated electron density. The maximum critical temperature in the different compounds is also comparable, and reaches up to $T_c \simeq 11$ K in $MoS_2$ at $10^{14}$ electrons/cm$^{-2}$ [1,7]. The occurrence of gate-induced superconductivity in semiconducting TMDs is therefore a robust phenomenon, that manifests itself in a similar way in the majority of compounds in this class [8], including atomically

thin samples [6,7,9–11]. Despite its experimental robustness, the nature of the gate-induced superconducting state remains currently unknown.

A number of facts have been established experimentally. The in-plane critical magnetic field, for instance, significantly exceeds the Pauli limit [5,12], a consequence of the very strong spin-orbit interaction present in semiconducting TMDs [13–15] that causes the spin of the electrons forming Cooper pairs to point in the direction normal to the layers. This implies that superconductivity is of the Ising type [16], and that electron pairs have to form in a superposition of singlet and triplet states. It has also been established that the onset of the gate-induced superconducting state coincides with the filling of electronic states in the *Q* valleys [17,18]. Interestingly, in thick layers of common semiconducting TMDs, the conduction band edge at the *Q* point is located at an energy lower than that of the band edge at the *K* point [19,20], but in transistor devices the large perpendicular electric field that accompanies electron accumulation deforms the conduction band. This deformation causes the band edge at the *K* point to be shifted to lower energy than that at the *Q* point already at rather low electron density (typically, lower than $\simeq 10^{13}$ cm$^{-2}$), and experiments [17] correlating the evolution of normal-state transport with the occurrence of superconductivity indicate that gate-induced superconductivity first

<footnote>*Contact author: thibault.sohier@umontpellier.fr</footnote>

<footnote>†Contact author: alberto.morpurgo@unige.ch</footnote>
occurs when both the *K* and *Q* valleys are populated (note that intercalated MoS₂ bilayer devices show superconductivity of different nature already at lower density, with only the *K*/*K'* and valley occupied). It has also been concluded from gate-dependent Raman spectroscopy experiments [21] that a very large enhancement of the electron-phonon coupling strength takes place in transistors of all the investigated semiconducting TMDs when the *Q* valley starts to be populated by electrons. Experiments therefore strongly suggest that the transition to the gate-induced superconducting state is driven by electron-phonon coupling.

Additional experiments providing useful indications, albeit with less clear-cut microscopic implications, have also been reported. For instance, the critical temperature of the gate-induced superconducting state appears to exhibit a domelike shape as a function of accumulated electron density *n*, reaching a maximum at an optimal value of *n*, before decreasing for larger electron accumulation. Details, however, seem to depend on the material and on the specific experiment reported. In some cases the “dome” is unambiguously present (e.g., in MoS₂ [1]) whereas in others it is less pronounced, i.e., it is less clear whether *T*<sub>c</sub> actually decreases or just tends to saturate (e.g., for the highest mobility WS₂ monolayers [6,11]). Systematic tunneling conductance measurements have been performed on MoS₂ transistors and have shown a large subgap density of states, with the zero-bias tunneling conductance scaling linearly with temperature down to *T* ≪ *T*<sub>c</sub> [22]. Such a behavior is incompatible with the presence of a fully gapped superconducting state, and suggests that gate-induced superconductivity has unconventional character. Whether this behavior is common to all TMDs is, however, not yet established, as controlled tunneling conductance measurements have so far been performed exclusively on MoS₂.

Our goal here is to investigate theoretically the nature of the superconducting state by combining first-principles calculations with Bardeen-Cooper-Schrieffer (BCS) theory [23], to rationalize the experimental observations. Semiconducting TMDs are particularly well suited to implement such a strategy because a large body of work has already been done to analyze their normal-state properties in terms of *ab initio* calculations [24–27], and to compare results to experiments [28,29]. This past work has established that all important aspects of the low-energy band structure of common semiconducting TMDs, its dependence on perpendicular electric field, and details of the electron-phonon and Coulomb interactions in the presence of accumulated electron density are correctly captured by first-principles calculations, with an overall excellent quantitative agreement. It should therefore be expected that *ab initio* simulations also allow capturing the key ingredients needed to describe the gate-induced superconducting state, namely, the multivalley structure of these materials, the strong and strongly momentum-dependent electron-phonon interaction, and Coulomb repulsion.

Taking into account these ingredients at a sufficient level of detail is important because most earlier attempts to describe gate-induced superconductivity in semiconducting TMDs were based on too drastic approximations. Examples include early work that largely overestimated $T_c$ by only considering the effect of electron-phonon interaction [30,31], or phenomenological models that included exclusively Coulomb

repulsion to analyze the possible occurrence of unconventional superconductivity [32,33]. Typically, this early work only considered states in the *K* and *K'* valleys, which, as we now understand, is not appropriate to explain experimental results. Over the years, more refined models, even generalized to hole doping [34,35], have been introduced, for instance including spin and density fluctuations [36], employing quantum Monte Carlo [37] or using a valley-dependent $\mu$ parameter to include Coulomb repulsion effects [38], but still with limited agreement with experiments. Very recently, a first-principles approach [39] based on a momentum-resolved Migdal-Eliashberg theory [40,41] and including anharmonic effects [42] has been applied to the case of TMDs [43], but with a treatment of Coulomb repulsion that misses the crucial role of nonlocal interactions in two dimensions (2D) [44] and multivalley effects [45].

Focusing on monolayers, we formulate a realistic model of multivalley effects within the simplest BCS approach, by employing *ab initio* techniques to analyze the momentum-resolved electron-phonon and Coulomb interaction as a function of accumulated electron density. The results of this analysis are used to determine the interaction coupling constants for intravalley and intervalley processes responsible for electron pairing, which are subsequently included in the BCS equation for the critical temperature. Such an approach represents a sensible approximation at the electron densities reached by electrostatic gating, which are much smaller than those of common metallic superconductors (and which therefore result in relatively small Fermi surfaces). It is particularly effective, as it allows important physical effects to be captured while keeping the formalism simple, which greatly facilitates their interpretation.

We find that the superconducting transition is driven by the large enhancement of the electron-phonon interaction strength that occurs when the *Q* valleys start to be populated, in agreement with both existing Raman spectroscopy experiments [21] (which directly probe electron-phonon interaction) and *T*<sub>c</sub> measurements [1–7,17] (the value of *T*<sub>c</sub> that we obtain are comparable to the experimental ones). We further find that, even though superconductivity is driven by electron-phonon interaction, Coulomb repulsion can profoundly affect the nature of the superconducting state, as it can cause a sign change in some of the intervalley interaction coupling constants. When that happens, the relative sign of the order parameter in different valleys becomes opposite and the superconducting state acquires an *s*<sub>+</sub> character. Our model therefore supports an unusual scenario, in which *s*<sub>+</sub> superconductivity occurs despite the dominant interaction driving superconductivity being electron-phonon attraction. We discuss qualitatively how, in actual devices, such scenario may account for some of the observed physical phenomena for which no microscopic explanation has been found so far, such as the presence of a large subgap density of states, and the sample-dependent dome-shaped dependence of *T*<sub>c</sub> on accumulated electron density.

The paper is structured as follows. In Sec. II we first introduce the key quantities and ingredients entering the theoretical model and their evaluation from first principles. The model predictions concerning gate-induced multivalley superconductivity in semiconducting TMD monolayers are then
<figure><img src="image_3.png" /><figcaption>FIG. 1. Spin and valley structure of WS₂. Opposite spins are projected on opposite sides of the valleys. We assume independent spin channels and consider only one spin (blue) in the model. We group the 3 *Q* (or *Q*') valleys that are degenerate in energy. Note that only two nonequivalent *K* points are represented, such that all states of the first BZ are represented once and only once. The scattering of a Cooper pair is also represented. The Cooper pair initially consisting of states (*k*, ↑) and (−*k*, ↓) is scattered into the pair (*k* + *q*, ↑) and (−*k* − *q*, ↓).</figcaption></figure>

presented in Sec. III. Finally, in Sec. IV we summarize the main findings and discuss their implications towards the interpretation of existing experiments.

## II. MODEL

This work focuses on semiconducting TMD single layers. In this section, we first describe the rich electronic structure of those monolayers and then propose a BCS-type model for superconductivity adapted to their multivalley character. This model is able to capture how the interplay between electron-phonon attraction and Coulomb repulsion determines the pairing within and between valleys. It is then reported how those two competing interactions are computed from first-principles simulations, as a function of momentum and doping (i.e., Fermi level).

### A. Multivalley structure of TMDs

Semiconducting TMD monolayers display a rich electronic structure, with multiple valleys [27] and a distinct spin texture [14,46], both in the valence and in the conduction bands. In the following we shall consider mainly electron doping, and we thus focus on the lowest conduction states that are occupied by electrons upon gating. In Fig. 1 we show the corresponding bottoms of the conduction bands across the hexagonal Brillouin zone (BZ), as computed in the case of WS₂ by density-functional-theory (DFT) simulations (see the Appendix for more details). Multiple valleys are present, whose structure is further enriched by a sizable spin-orbit coupling that, combined with the lack of inversion symmetry, gives rise to a net spin splitting. The spin polarization is essentially orthogonal to the layers, with spin-up and spin-down states represented, respectively, in blue and red in Fig. 1. Although results are illustrated with the specific

case of WS₂, the valley and spin structure is quite general and applies to all the TMDs considered here (MoS₂, MoSe₂, WS₂, and WSe₂), with only the quantitative details for the energy separation between different valleys and different spins varying [19,26].

The lowest-energy conduction states are located around two high-symmetry points at the opposite BZ vertices, typically denoted as *K* and *K'*, giving rise to two inequivalent valleys (the other BZ vertices being equivalent by a reciprocal lattice vector translation). The spin splitting in these valleys is small and opposite at *K* and *K'* because of time-reversal symmetry. At *K*/*K'*, the splitting has different signs for W-based and Mo-based TMDs, but spin-split bands cross in Mo-TMDs relatively close to the band edge because they have different effective masses, resulting in the same spin ordering of the bands sufficiently far from the BZ vertices in the energy range relevant for this study. Additional valleys appear higher in energy, approximately located midway between the BZ vertices and the BZ center *Γ*, around six inequivalent *k* points typically denoted with the letter *Q* (or less frequently *Λ*). For these valleys the spin splitting is much larger and alternates in sign in a way that is consistent with the threefold symmetry and time-reversal invariance. Overall, we thus have four distinct band edges, corresponding to the *K* and *Q* points for the two spin polarizations, giving rise to a total of 16 pockets, as shown in Fig. 1: 4 = 2 × 2 arising from spin-up and -down states around *K* and *K'* and 12 = 2 × 6 associated with opposite spins at the six *Q* valleys.

Given that Coulomb interaction does not flip spin and spin-conserving electron-phonon interactions are much stronger than the spin-flipping ones even in the presence of spin-orbit coupling in TMDs [47], the spins of individual electrons comprising Cooper pairs are assumed to be conserved during scattering. The spin-orbit interaction in these materials is of Ising type, with an out-of-plane spin-quantization axis. Opposite spin pairing is expected to be favored energetically due to the time-reversal symmetry, which ensures that there is an equal energy partner with opposite spin and momentum available for each electron (note that this pairing is neither singlet nor triplet, but rather a linear combination of them [48]). In this context, spin-up and -down channels can be considered to be equivalent and independent. In the following we shall therefore focus on a single channel and we identify Cooper pairs by just one of the partners, that we take to be spin up, i.e., we consider the blue valleys in Fig. 1.

For what concerns the BZ corners, we thus have a valley at *K* and one at *K'*, with the same spin but different energies. The six *Q* valleys with spin up can be separated into two groups that lie at different energies, as shown in Fig. 1. The valleys within each group are related by the threefold rotation symmetry, which enforces their degeneracy. In the following we shall consider superconducting solutions where this rotation symmetry is preserved, so it is convenient to pool together valleys within a given group into a single set. The lowest three valleys, which form a triangle that points towards the *K'* point, will be called *Q*, while the upper three, pointing towards the *K* point, will be denoted *Q'* (the choice of prime notation will become clear when discussing electron-phonon interactions). The superconducting gap equations will be written assuming
that the electronic states belong to one of these four valleys: K, K', Q, and Q'.

The relevance of these four valleys depends on the doping level, which can be efficiently tuned in semiconducting TMDs using electrostatic gating [18,49–51]. Mimicking experiments, the doping density (and Fermi energy) will thus be the main external parameter that we tune in our simulations. At low doping, only the K/K' valleys are occupied. With increasing density also the Q valleys start to get populated. The precise value for the onset of Q occupation with density does not affect qualitatively the present analysis and we do not claim to be predictive as this onset depends on the energy difference between the K and Q valleys, which is not known precisely at a quantitative level [52] (large variations in the values are reported in both experiments and DFT simulations). At even larger density, the Q' valleys could also be occupied. However, the energy difference between Q and Q' is generally large and it is unlikely that electrons visit the Q' valleys in existing experiments, i.e., it is unlikely that the Q' valley can be populated at the maximum density accumulated by ionic gating. Q' states are nevertheless kept in the model for completeness.

## B. Basic aspects of multiband SC and microscopic interactions

We now describe our approach to modeling superconductivity. We neglect finite-momentum pairing, which is unlikely in the trigonal-warped conduction bands at high carrier density [34]. The Ising spin-orbit interaction guarantees that for every electron with momentum spin (k, σ) there is also an electron with the same energy at (−k, −σ). These electrons can form zero-momentum Cooper pairs via a weak coupling instability, in a way directly analogous to the standard BSC superconductivity. What makes the situation nonstandard in TMDs is that there are multiple Fermi surfaces, controlled by the doping level selected by the voltage applied to the ionic gate, such that the superconducting order parameter may vary both within each one of them and even more among them. Given that these Fermi surfaces are rather small relative to the Brillouin zone size, we can use a simplified description where each Fermi pocket has a unique value of the order parameter Δ<sub>v</sub> (v is the valley index).

Different Δ<sub>v</sub> are connected to each other by the self-consistent BCS gap equations. These equations generalize the standard BCS equations for a single band [53] and have the form [45,54]

$$ \Delta_v = M_{vw} \Delta_w. \quad (1) $$

In general, the matrix M depends both on temperature and the order parameters Δ<sub>v</sub> themselves, M = M(T, Δ<sub>v</sub>). Setting the temperature to zero leads to a nonlinear system of equations for the order parameter in the ground state; setting Δ<sub>v</sub> to zero in M yields a linearized equation in the order parameters that allows computing the superconducting transition temperature det[M(T<sub>c</sub>) − 1] = 0.

Graphically, the matrix M can be represented as depicted in Fig. 2.

Individual coupling matrix elements include, for a pair of valleys vw, the competition between phonon-mediated

<figure><img src="image_4.png" /><figcaption>FIG. 2. Multivalley superconductivity is modeled via a matrix containing the average interactions between electrons in different valleys.</figcaption></figure>

attraction (G<sub>vw</sub>) and Coulomb repulsion (C<sub>vw</sub>):

$$ \begin{align} M_{vw} &= \ln \left[ \frac{\omega_{vw}^{\text{eff}}}{k_b T_c} \right] (G_{vw} - C_{vw}), \nonumber \\ G_{vw} &= \left\langle \sum_{\mathbf{k}' \in w} \sum_{\lambda} \frac{g_{\mathbf{k}\mathbf{k}', \lambda}^2}{\omega_{\mathbf{q},\lambda}} \delta(\varepsilon' - \varepsilon_F) \right\rangle_{\mathbf{k} \in v}, \nonumber \\ C_{vw} &= \alpha \times \left\langle \sum_{\mathbf{k}' \in w} V(q = |\mathbf{k}' - \mathbf{k}|) \delta(\varepsilon' - \varepsilon_F) \right\rangle_{\mathbf{k} \in v} \tag{2} \end{align} $$

with v, w ∈ {K, K', Q, Q'}. M<sub>wv</sub> represents an average matrix element for a Cooper pair scattering between valleys v and w, all near the Fermi level ε<sub>F</sub>. A constituent electron of momentum k and energy ε scatters to state k' = k + q (ε') either via phonons or Coulomb interaction. ⟨f(k)⟩<sub>k∈v</sub> stands for the average of f over k momenta in v at ε<sub>F</sub>. However, very minimal dependence on the choice of k ∈ v was found, as shown in methods, due to the fairly small Fermi surfaces. In practice, one reference state k ∈ v is picked from the ε = ε<sub>F</sub> isoenergetic line. The matrix M is computed for each material at various carrier concentrations by evaluating from first principles the phonon and the Coulomb components of M as described in the Secs. IIC and IID. In particular, we compute the electron-phonon matrix element g<sub>kk',λ</sub>, which describes the coupling between a phonon of mode λ at momentum q and frequency ω<sub>q,λ</sub> with two electronic states at k and k' = k + q. To account for Coulomb scattering, we compute the screened Coulomb potential within the isotropic approximation

$$ V(q) = \frac{2\pi e^2}{q\epsilon(q)}, \quad (3) $$

where q = |q| and ε(q) is the dielectric function in the static approximation. It includes screening by itinerant electrons (see Sec. IID for details), and offsets the attractive interaction provided by phonons. Within the BCS treatment, it thus comes in with the same energy range of integration and therefore the same logarithmic prefactor as the attraction G.

There is a possibility of a more accurate Eliashberg-type treatment [41] that allows for fully frequency-dependent order parameter, up to the electronic bandwidth frequencies, which can be adapted to a first-principles approach [39,55]. In that case, the bare Coulomb interaction is taken as unscreened and the screening is generated automatically in the course of solving the frequency-dependent Eliashberg equations for the normal and anomalous Greens functions and self-energies [56]. Our approach is more qualitative and
simple: we assume that the dominant pairing is due to phonons, which defines the relevant frequency window around the Fermi level, while the Coulomb interaction is taken as electronically screened. For the latter we take advantage of the first-principles methods as well; this allows us to be material specific not only as far as the electron phonon interaction is concerned, but also for the Coulomb (pseudo)potential. The limitation of our approach is that is only justifiable if phonon-mediated attraction is the dominant pairing mechanism.

To uncover the role of Coulomb interaction, we introduce a scaling factor $\alpha$. It allows us to gradually turn on the Coulomb repulsion and explore how the nature of superconductivity changes with the interaction strength. In practice, the scaling factor $\alpha$ can approximately capture the dependence of Coulomb on the dielectric properties of the surrounding substrate or ionic liquid. Note that a more sophisticated dynamical treatment of electron-phonon and Coulomb interaction would essentially fix the value of $\alpha$. Nevertheless, given that this value would strongly depend on the dielectric environment, the use of a tunable parameter brings more generality to the model.

Each $M_{vw}$ has its own effective phonon ("Debye") frequency under the logarithm, which corresponds to the energy of the typical phonons at the momenta involved in the intervalley or intravalley transitions. In practice we average over the phonon frequencies weighted by their coupling to the electrons:

$$ \omega_{vw}^{\text{eff}} = \frac{\sum_{k' \in w, k \in v} \sum_{\lambda} \int \frac{g_{kk',\lambda}^2}{\omega} \delta(\varepsilon' - \varepsilon_F)}{\sum_{k' \in w, k \in v} \sum_{\lambda} \int \frac{g_{kk',\lambda}^2}{\omega} \delta(\varepsilon' - \varepsilon_F)}. \quad (4) $$

The superconducting transition is signaled by the appearance of a nonzero solution in the linearized (i.e., $\Delta \to 0$) system of equations; it corresponds to the temperature for which the matrix $M$ develops a unit eigenvalue. The corresponding eigenvector represents the relative signs and sizes of the order parameter in the different valleys. For a purely attractive interaction, by Perron-Frobenius theorem, all components $\Delta_v$ of the eigenvector with the largest eigenvalue have the same sign. This corresponds to the standard s-wave state, adjusted for the fact that there is Ising spin-orbit coupling. However, the elements of $M$ can change sign for some valley pairs if the Coulomb repulsion exceeds the phonon-mediated attraction. In this case, the largest eigenvalue (the highest transition temperature) can have an eigenvector with both positive and negative elements, implying an order parameter that has different signs on different Fermi sheets. This would indicate the formation of a $s_\pm$ state.

The added value of the approach described above is that it provides a good compromise between complexity and accuracy. Electron-phonon interactions are highly dependent on both the magnitudes and the relative angles between momenta. In relatively crude estimations of superconducting properties, these interactions are averaged over all Fermi surfaces and combined into a single parameter [9,30,57] (usually called $\lambda$). The absence of valley resolution then limits the physical insight into the multivalley aspect of the coupling mechanisms. In more sophisticated estimations, the full momentum dependency is retained within the Eliashberg framework [43] (see also [58,59] in other families of TMDs).

While quantitatively more complete, the size and complexity of the resulting momentum-dependent order parameter makes it difficult to draw qualitative conclusions. Also, the role of different phonons can become obscured. Here, we retain only the principal momentum structure of the couplings in the superconducting problem by averaging over valleys. This level of coarse graining is found to be simple enough to provide physical insight, and detailed enough to understand the underlying multivalley mechanisms. Furthermore, we have the full momentum and energy dependence of the couplings and can identify the phonons leading to dominant couplings within a given valley pair. Concerning the Coulomb repulsion, most works simply use a constant parameter [31,57] $\mu$. This is not satisfactory, as already demonstrated by earlier works that do account for the variation of $\mu$ between intravalley and intervalley processes [38]. Recognizing the importance of the momentum dependence of interactions, here we treat the Coulomb repulsion on the same footing as phonon-mediated attraction, defining a momentum-dependent interaction and summing it on the Fermi surface.

### C. Phonon-mediated attraction

The quantities needed to estimate the phonon-mediated attraction $G_{vw}$ from first principles are the values of the squared electron-phonon coupling divided by the phonon frequency $g_{kk',\lambda}^2/\omega_{q,\lambda}$ at the Fermi level and between states in valleys $v$ and $w$. As we have recently shown [21,52], and will confirm in this section, the important feature of electron-phonon interactions in TMDs is the activation of strong intervalley $KQ$ and intravalley $KK$ couplings when carrier density is increased enough to occupy the $Q$ states.

Here we consider the case of high doping, when the Fermi level crosses all valleys such that all scattering transitions are present. The values of $g_{kk',\lambda}^2/\omega_{q,\lambda}$ between states at the Fermi level are shown in Fig. 3 for MoS$_2$, although similar results are obtained for all other TMDs. Differently from what is typically done in the literature, this is computed by including doping explicitly. This is achieved by performing density functional theory simulations (DFPT) in a field-effect setup that allows for doping [60] and constructing a model that takes into account the unusual dependency of screening on the valley occupations [21,52], as detailed in Appendix 2.

To identify the most relevant couplings, Fig. 3 shows $\sum_{\lambda} g_{k_i, k_f}^2 / \omega_{k_f - k_i, \lambda}$ for initial states $k_f$ in the different valleys as a function of the possible final states $k_f$ on the Fermi surface. For an initial state in $K$ (main panel), intravalley $KK$ and intervalley $KQ$ couplings are stronger than $KK'$ or $KQ'$ by nearly an order of magnitude. Moreover, the couplings are rather uniform when considering final states within a given valley, which justifies the approximation of dealing with valley-averaged quantities on the Fermi surface.

The other smaller panels in Fig. 3 show the same quantity but for initial states in different valleys. This highlights the symmetry in the pairs of strongly coupled valleys. While the main panel showed that states in $K$ mostly couple to the $K$ and $Q$ valleys, we can now see that states in $K'$ mostly couple to the $K'$ and $Q'$ valleys. Inversely, states in $Q$ ($Q'$) are mainly coupled to $K$ ($K'$). We can now rationalize the prime notation as a pairing of strongly phonon-coupled valleys ($Q$ with $K$ and
<figure><img src="image_5.png" /><figcaption>FIG. 3. Electron-phonon interactions between states at the Fermi level for MoS$_2$ at a high electron doping of $n = 2 \times 10^{14}$ cm$^{-2}$. $a_0$ is the Bohr radius. An initial state $\mathbf{k}_i$ (black star) is picked in one of the four valleys and the strength of the phonon-mediated interaction $\sum_\lambda g^2_{\mathbf{k}_i, \mathbf{k}_f / \omega_{\mathbf{k}_f - \mathbf{k}_i, \lambda}}$ is computed for final states $\mathbf{k}_f = \mathbf{k}_i + \mathbf{q}$ at the Fermi level. Here the Fermi level crosses all four valleys, but at lower doping that might not be the case.</figcaption></figure>

$Q'$ with $K'$. $QQ'$ (and $Q'Q$) couplings are also strong in the case of MoS$_2$. However, the carrier density necessary to reach the $Q'$ valley is not likely to be achieved experimentally. The most relevant result is thus the predominance of intravalley ($KK$ and $K'K'$) and intervalley ($KQ$ and $K'Q'$) couplings with respect to intervalley $KK'$.

The origin of the strong intravalley $KK$ and $K'K'$ couplings is the $A_1$ optical phonon mode, corresponding to a breathing mode involving the out-of-plane movement of the chalcogens, Fig. 4. This can be described as a variation of the layer's thickness, which directly implies a variation of the electronic density and a strong electron-phonon interaction akin to a deformation potential. In addition, as shown in [21,52], the strong electron-phonon coupling of this mode is not screened by free carriers in the conduction band when the Fermi level is above the $Q$ valley.

The origin of the strong $KQ$ (and $K'Q'$) coupling is a particularly large interaction with the zone border longitudinal acoustic phonon around the $q = M$ point. This applies to all the TMDs studied here. Similarly to the zone center $A_1$ mode, its strong coupling can be linked to a variation of the layer's thickness, but this time at a momentum that links $K$ and $Q$ valleys. More precisely, this mode is characterized by in-plane displacements of the transition metal combined with out-of-plane breathing displacements of the chalcogens

(see, e.g., phonon visualizer [61], and Fig. 4). In addition to a variation of the electronic density, this pattern seems to maximize the variations in the angle of the bonds between transition metal (Mo, W) and chalcogens (S, Se). The K and Q valleys correspond to out-of-plane and in-plane orbitals [62], respectively. Strong variations of the bond angle are consistent with strong coupling between in- and out-of-plane orbitals, and thus transitions between K and Q valleys.

This strong coupling comes with a doping-dependent softening of the LA phonon at $q = M$ that is not included in the above results. More precisely, the softening starts when $Q$ is occupied, then gradually increases with doping [43]. The softening is difficult to evaluate accurately for different temperatures and dopings, especially due to anharmonic effects [43]. Here, this softening will be approximately included by applying a uniform 30% reduction to the corresponding phonon frequency. This corresponds to the result of the anharmonic calculations in Ref. [43] for MoS$_2$ at $n \sim 10^{14}$ cm$^{-2}$ where superconductivity is typically found. This only impacts our results when $Q$ is occupied, and we note that the softening could very well be more pronounced at higher dopings.

Figure 3 shows electron-phonon couplings at a high electron doping of $n = 2 \times 10^{14}$ cm$^{-2}$. For completeness, it is useful to describe the low doping case, when the $Q$ and $Q'$ valleys are not yet crossed by the Fermi level. In that situation, the only relevant couplings involve the $K$ and $K'$ valleys. Comparing with their respective high doping values, the $KK'$ coupling is similar whereas the $KK$ coupling is significantly weaker. In fact, $KK$ and $KK'$ couplings are comparable at low doping. A more quantitative study of doping dependence is performed in Sec. III A.

The data shown in Fig. 3 can be collected to reconstruct the matrix elements $G_{vw}$ upon averaging over the initial states in valley $v$ and summing over final states in valley $w$, Eq. (2). In practice, as already mentioned, the limited variation with initial state within a valley (see also Fig. 9) makes the averaging procedure not necessary, in favor of taking just one representative initial state. On the other hand, the sum over

<figure><img src="image_6.png" /><figcaption>FIG. 4. Displacement patterns associated with the strongest coupling phonon modes, the $A_1$ mode at $\Gamma$ and the zone border LA mode at $M$. $M$ and $X$ stand for transition metal and chalcogen, respectively.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 5. Momentum and doping dependency of the static, screened Coulomb interaction in electron-doped MoS$_2$. $a_0$ is the Bohr radius. The color bar shows the Fermi level, with the bottom of the K' valleys as a reference. The position of the valleys are indicated.</figcaption></figure>

final states at the Fermi level is indeed performed. We use a triangular integration technique [63] rather than replacing the Dirac deltas $\delta(\epsilon' - \epsilon_F)$ by some finite-width distribution which would imply the need for convergence studies.

### D. Coulomb repulsion

The main ingredient to obtain the Coulomb repulsion contribution $C_{vw}$ from first principles is the screened Coulomb interaction $V(q)$ between states at the Fermi level in valleys $v$ and $w$. We compute this quantity as a function of the momentum transfer between an initial and a final state $q = |\mathbf{k}_f - \mathbf{k}_i|$, as in Eq. (3), that is as the bare Coulomb potential divided by the macroscopic static dielectric function $\epsilon(q)$ of the doped TMD. The latter is obtained within the random-phase approximation (RPA), assuming the independent particle susceptibility to be the sum of contributions from the neutral material and from the extra free carriers [52]. Concerning the latter contribution, we are mostly interested in the metallic screening effects when $q$ is small enough to correspond to intravalley transitions. At large momenta, we assume the dielectric function to be dominated by interband transitions, i.e., those that make up the dielectric function of the neutral material. Those assumptions are reflected in the choice of a simple complementary error function to model the wave-function overlaps in the calculation of the free carriers contribution (see Appendix 2 b). Furthermore, we assume an isotropic dielectric function. Since the Fermi surfaces of all valleys are roughly isotropic, this is a good approximation for small momenta corresponding to intravalley transitions. It is also a good approximation at larger momenta where the interband contributions are assumed to be dominant.

Figure 5 shows results for the static screened Coulomb interaction $V(q)$ for different Fermi levels, i.e., different doping concentrations. Intravalley and intervalley processes are characterized by a very different momentum dependence. At small $q$, relevant for intravalley scattering, and in particular for $q < 2k_F$, with $k_F$ being the Fermi wave vector, there is a strong momentum dependence that arises from the metallic screening associated with free carriers. On the contrary, in the regime of large momentum typical of intervalley processes,

and in particular for $q > 2k_F$, the screened interaction is slowly varying with $q$, with a behavior that is reminiscent of three-dimensional (3D) dielectric screening, as indeed at large momentum $1/q$ is smaller than the thickness of the 2D layer.

A distinct behavior between intravalley and intervalley processes is found also for what concerns the doping dependence. Intravalley coupling is very sensitive to doping as in the long-wavelength limit the screened Coulomb potential $V(q \to 0)$ is inversely proportional to the density of states (DOS), as captured already by the Thomas-Fermi approximation. At low doping, only the K/K' valleys are occupied, with a relatively small DOS. On the contrary, at high doping also the Q valleys get occupied, with a much larger DOS (both because of the higher degeneracy and a heavier effective mass), with a significant suppression of the screened interaction. The suppression associated with metallic screening extends to finite $q$ but still within the Fermi surface, with a flattening of $V(q)$ and a reduction of Coulomb repulsion. Intervalley interactions characterized by large momentum transfer $q$ are not affected by metallic screening, and the Coulomb repulsion is doping independent.

The screened Coulomb interaction is then summed over possible final states on the Fermi surface for a given valley to obtain the Coulomb matrix elements $C_{vw}$ through Eq. (2).

### III. RESULTS

The procedure outlined above to compute $G_{vw}$ and $C_{vw}$ is repeated for different doping concentrations and in different materials. The linearized form of Eq. (1) is then solved to find the corresponding $T_c$, and from the eigenstates we gain information on the nature of the superconducting state. However, before illustrating the main results, it is important to stress that the gap equations depend on the value of the $\alpha$ parameter, which tunes the relative strength of Coulomb repulsion and phonon-mediated attraction. The model in Sec. II is intended to describe phonon-driven superconductivity, with Coulomb repulsion as a perturbation. The starting point is thus $\alpha = 0$, corresponding to the complete neglect of Coulomb interactions. In that case the gap equations are fully determined by phonon attraction and any superconducting solution found in this regime is obviously purely phonon driven. This is still true for finite but small values of $\alpha$. In the opposite limit of large $\alpha$, all $C_{vw}$ matrix elements eventually become larger than the corresponding $G_{vw}$, meaning that the valley couplings are all repulsive, and the gap equations, as written in Eq. (1) assuming a dominating phonon attraction, fail to be valid.

We thus restrict our investigations to a maximum value of $\alpha = 0.3$. Although there are small variations among TMDs, this roughly corresponds to the point where the strongest phonon-mediated coupling ($G_{KQ}$) is half compensated by the Coulomb repulsion ($M_{KQ} = G_{KQ} - \alpha C_{KQ} \sim G_{KQ}/2$), while the second strongest ($G_{KK}$) is fully compensated ($M_{KK} \sim 0$), all other elements of $M$ being already negative. Note that those criteria only make sense when $Q$ is occupied, but as we will see, that is a necessary condition for superconductivity anyway, even for $\alpha = 0$. The Coulomb repulsion was simply computed from the screened Coulomb interaction for single-layer TMDs in vacuum. One may consider its reduction by $0 \le \alpha \le 0.3$ as including the contribution from a dielectric
<figure><img src="image_8.png" /></figure>

FIG. 6. Evolution of the intervalley phonon-mediated attraction, Coulomb repulsion, and their difference with electron doping in MoS$_2$. The horizontal axis shows when the Fermi level crosses the bottom of each valley. The values of the Fermi levels and carrier densities corresponding to each point can be found in Fig. 7. Here, the tuning prefactor $\alpha$ for the Coulomb repulsion is set at 0.2. Only active couplings, i.e., between valleys that are populated at that Fermi energy, are reported.

environment. In particular, in the common case of an ionic gate, the liquid's highly polarizable ligands are very close to the layer, closer than the average distance between the added carriers in the conduction band.

In the following we set $\alpha = 0.2$ when studying the doping dependence of the interaction matrix elements. The observed trends are actually robust for a wide range of $\alpha > 0.05$, and $\alpha = 0.2$ is simply chosen to make the effects of Coulomb repulsion quantitatively obvious, without reaching the maximum considered as the threshold for validity ($\alpha = 0.3$). The solutions to the gap equations are then investigated for different doping levels and at five representative values in the entire range of $\alpha$ (= 0, 0.05, 0.1, 0.2, 0.3).

A. Doping dependence of interaction matrix elements

The matrix elements of Eq. (2) describe phonon-mediated attraction and Coulomb repulsion. Figure 6 reports representative results for the case of MoS$_2$, for $\alpha = 0.2$, showing the relative importance of those two contributions as a function of doping. For simplicity, we plot only the matrix elements with $K$ as a first index. Those, along with the symmetry of the electron-phonon and Coulomb interactions [$\frac{g_{KK'+q}^2}{\omega_{q,\lambda}}$ and $V(q)$]

with respect to valley indices, are enough to describe the trends for the full matrix.

The top panel of Fig. 6 shows the evolution of $G_{K,K}$, $G_{K,K'}$, $G_{K,Q}$, and $G_{K,Q'}$ as a function of the Fermi energy. At low doping, the Fermi level does not cross the $Q$ valley, such that only matrix elements involving $K$ and $K'$ can be nonzero. Those are finite but very small ($<0.01$) on the scale of Fig. 6.

Interestingly, however, intravalley $KK$ (and $K'K'$) coupling becomes strong as $Q$ is occupied, with an enhancement by more than an order of magnitude. This is the consequence of the peculiar multivalley screening effect discussed earlier by some of us [21,52]. In particular, electron-phonon coupling is dominated by the $A_1$ phonon mode. The corresponding coupling mechanism can be thought of as a shift in energy of the $K$ and $Q$ valleys following the periodic phonon displacement pattern and out-of-phase with one another. The relative amplitude of the $K$ and $Q$ shifts are such that the charge density perturbations associated to each valley, opposite in sign, are equal in amplitude. Thus, they compensate each other. With a net charge density perturbation that vanishes, the screening from the conduction band electrons also vanishes when $Q$ is occupied, leading to a large electron-phonon coupling.

Intervalley $KQ$ scattering is obviously present only when doping is high enough for the Fermi level to cross the $Q$ valley. Once active, it is an order of magnitude stronger than the others. We have verified numerically that, as expected from symmetry, primed valleys $K'$ and $Q'$ also couple strongly to each other. The appearance of the large phonon-mediated couplings discussed earlier, which follow the prime notation, thus coincides with the occupation of the $Q$ valley.

The middle panel of Fig. 6 reports analogous results for Coulomb coupling matrix elements $C_{K,K}$, $C_{K,K'}$, $C_{K,Q}$, and $C_{K,Q'}$, with the aforementioned choice of $\alpha = 0.2$. The intervalley ($KK'$) component is only slightly dependent of the Fermi level location. This is consistent with the constant $V(q)$ in Fig. 5 associated to a slight increase of the density of states as the Fermi level increases. The intravalley component $C_{KK}$ instead gradually decreases with doping, from a factor of 2 above the intervalley $C_{KK'}$ to slightly below it. When the doping level is increased so that also the $Q$ valley is occupied, the suppression of $C_{KK}$ is slightly accelerated. More importantly, intervalley couplings to the $Q$ valley become active, and they are significantly larger than the other matrix elements. This is because the main ingredient determining the relative magnitude of the matrix elements from valley to valley is the DOS in the final valley $w$. In particular, the larger DOS of $Q$ and $Q'$ makes $C_{v,Q}$ and $C_{v,Q'}$ matrix elements significantly larger. Importantly, the Coulomb matrix elements are not sensitive to the primed or unprimed distinction between the valleys. Not only $KK$ and $KK'$ matrix elements have comparable orders of magnitude, but also $KQ$ and $KQ'$ components have similar values.

Plotting $G_{vw} - C_{vw}$ in the bottom panel of Fig. 6 highlights some interesting features of the competition between phonon-mediated attraction and Coulomb repulsion. Indeed, some elements of $M$ change sign as a function of doping. At small carrier concentration, the matrix elements are either null when the corresponding intervalley transition is not on the Fermi surface, or small and negative when Coulomb repulsion dominates. As $Q$ is occupied, both $G_{KQ}$ and $G_{KK}$ increase
<figure><img src="image_9.png" /><figcaption>FIG. 7. Critical temperature as a function of Fermi level (with respect to the bottom of K valley) for the four TMDs, for different values of the prefactor to Coulomb repulsion, $\alpha$. The top axis indicates the electrostatic doping level corresponding to each point. Markers indicate the type of solution found. The *Q* and *Q'* valleys of MoS₂ have been lowered by 100 meV with respect to the DFT results.</figcaption></figure>

strongly, overcoming the Coulomb repulsion and leading to positive $M_{KQ}$ and $M_{KK}$ elements. As we will see in the next section, superconductivity emerges around those conditions.

Phonon-mediated attraction does not overcome Coulomb repulsion for all valley pairs: $M_{KK'}$ and $M_{KQ'}$ remain negative. Strong phonon-mediated couplings follow the prime notation, while the Coulomb repulsion is similar among and between primed and unprimed valleys, and varies mostly according to the density of states associated with the second valley index. As a result, we observe different outcomes of their interplay in different valley pairs. In particular, once *Q* is occupied, phonon-mediated attraction systematically prevails among primed or unprimed valleys, while the *KK'* coupling remains repulsive. The *KQ'* coupling is also negative when active, although it is less relevant in practice since the *Q'* valley is unlikely to be occupied in typical experiments. The main effect of introducing Coulomb repulsion is thus to invert the sign of the valley couplings where phonon-mediated attraction is weak, while keeping the couplings attractive between the most strongly phonon-coupled valleys.

The pivotal sign inversion of the *KK'* coupling, illustrated here with $\alpha = 0.2$, actually happens over a broad range of values of $\alpha$, starting already from $\alpha \sim 0.05$ (depending on the TMD and doping level). This is ultimately due to the *KK'* phonon-mediated attraction being quite weak for any doping level. For values of $\alpha$ above the chosen threshold of 0.3, the dominance of phonon-mediated attraction is challenged even for *KK* and *KQ* couplings, a regime which is beyond the scope of this work. As we will see in the next section, the dominance of *KK* and *KQ* phonon-mediated attraction as a function of

doping, i.e., when *Q* is occupied, and the possible coexistence of attractive and repulsive couplings between different valley pairs are the key features to understand the onset and nature of superconductivity in TMDs.

## B. Superconducting solutions

With the *G* and *C* matrices computed, the linearized form of Eq. (1) is solved for all four TMDs considered here (MoS₂, MoSe₂, WS₂, and WSe₂), at different carrier concentrations, and for four representative values of $\alpha$ (0.0, 0.05, 0.1, 0.2, 0.3) within the range discussed above. The critical temperature $T_c$ is found as the highest temperature for which Eq. (1) allows for a nontrivial solution. We only look for solutions with a finite and experimentally relevant $T_c > 10^{-4}$ K. The behavior of $T_c$ as a function of doping and the nature of the solution is determined by the interplay between attractive phonon-mediated and repulsive Coulomb interactions. This interplay is different for each pair of valleys within the multivalley Fermi surface, and strongly depends on doping. A summary of all corresponding results is reported in Fig. 7.

### 1. Onset of superconductivity

There are several robust features common to all TMDs studied here. The onset of a critical temperature above the chosen threshold for experimental significance ($>0.1$ mK) systematically coincides with the Fermi level crossing the bottom of the *Q* valley. This is true whatever the value of the tuning parameter for the Coulomb repulsion $\alpha$. This holds in spite of the variations in the carrier density at which *Q* is
reached, i.e., it does not depend on the energy of the bottom of the *Q* valley with respect to the *K* and *K'* valleys. Clearly the emergence of phonon-driven superconductivity is triggered by the activation of *KQ* and the enhancement of *KK* couplings as the *Q* valley is occupied.

Once the *Q* valley is occupied, we find a fairly constant $T_c$ of a few Kelvins for all TMDs. This value is maximal in the absence of Coulomb repulsion ($\alpha = 0$) and gradually decreases as $\alpha$ is turned on, in a similar way for all TMDs. For the maximal value of $\alpha = 0.3$, $T_c$ is decreased by an order of magnitude.

The precise value of $T_c$ and its doping dependence are quite sensitive to the softening of the LA phonons at $q = M$, which is delicate to evaluate *ab initio*, notably due to anharmonic effects [43]. Indeed, the associated phonon-mediated coupling between the *K* and *Q* valleys is very large and drives the superconductivity. As already mentioned, for this work we simply apply a 30% reduction of the LA(*M*) frequency, based on the results of Ref. [43]. Realistically, one expects a doping-dependent softening, starting from the occupation of the *Q* valley and steadily increasing, likely past the 30% used here, with a corresponding increase in $T_c$. Nevertheless, we expect the current treatment to give a good estimation in the relevant doping range $n = 10^{13}-10^{14} \text{ cm}^{-2}$. Most importantly, it is sufficient to understand the underlying physics.

Thus, despite variations in their chemical composition, electronic bands and phonon properties, all the semiconducting TMDs considered here consistently display phonon-driven superconductivity with a $T_c$ on the order of a few Kelvins, the onset of which systematically coincides with the occupation of the *Q* valley.

## 2. Emergence of exotic solutions

While superconductivity is driven by electron-phonon interactions, the Coulomb repulsion plays an important role. Interestingly, even at small tuning parameter $\alpha$, it modifies the coupling between valleys significantly enough to change the nature of the superconducting solutions, as characterized by the eigenvector $\Delta_v$ corresponding to a unit eigenvalue in the linearized gap equation (1). Its components provide information on the relative sign and magnitude of the superconducting order parameter in each valley. In Fig. 7, this is illustrated by a switch from filled to empty markers, used to denote qualitatively different solutions. This consistently happens for all finite values of $\alpha$ considered.

To investigate further the change in the nature of the solution, Fig. 8 shows the components $\Delta_v$ as a function of $\alpha$. Doping is chosen to correspond to the standard experimental conditions in which superconductivity is observed, when the Fermi level is about ~10 meV above *Q*. For $\alpha = 0$, all components of $\Delta_v$ are positive. The *K* valley dominates, followed by *Q*, while *K'* has a small but positive contribution. For fairly small values of $\alpha \lesssim 0.05$, the *K'* component goes through zero then changes sign. This coincides with the change of sign for the $G_{KK'} - C_{KK'}$ matrix elements mentioned before. The *K* and *Q* components stay positive. Conventional phonon-driven superconductivity corresponds to a positive order parameter in all valleys, as obtained for $\alpha = 0$. The introduction of Coulomb repulsion thus systematically leads to the appear-

<figure><img src="image_10.png" /><figcaption>FIG. 8. Components of the eigenvector $\Delta_v$ corresponding to the superconducting solution as a function of the Coulomb prefactor $\alpha$, in $\text{MoS}_2$. The electron doping is fixed and corresponds to the second data point above *Q*, with $\epsilon_F \approx 0.14$ eV and $n = 8 \times 10^{13} \text{ cm}^{-2}$.</figcaption></figure>

ance of exotic solutions, i.e., an eigenvector with components of different signs, which corresponds to the order parameter changing sign for electrons in different valleys. Note, however, that the solutions we find do not break the point-group symmetry (*s* wave), which in conventional superconductors, such as lead or aluminum, ensures robustness with respect to disorder [64]. In particular, both solutions belong to the same $A'_1$ irreducible representation [65]. Although the exotic solution is sometimes referred to as *f* wave in the literature [32,34,54,65], we prefer to call it $s_{+-}$ to better convey the message that there are two different pockets of the Fermi surface (the *K* and *K'* valleys) with opposite sign, but also different magnitude, of the order parameter.

A necessary (but not sufficient) condition for the existence of those exotic solutions is to have coexisting attractive and repulsive elements in the matrix *M*. The strength of Coulomb relative to the phonon-mediated interaction determines the sign of the coupling between each pair of valleys. Depending on momentum, doping, and $\alpha$, the repulsive Coulomb term $C_{vw}$ overcomes the phonon-mediated term $G_{vw}$ in some but not all of the matrix elements. Due to their strong phonon-mediated coupling, *K* and *Q* valleys remain attractively coupled. The opposite-sign channels are made up of (*KQ*) and (*K'*$Q'$). Those channels are attractive within themselves, but repulsive between each other. This corresponds to a pairing parameter that has opposite signs in the *K* and *K'* valleys. The *Q* (*Q'* valley shares the sign of the *K* (*K'* valley, although with a smaller magnitude.

Going back to Fig. 7, note that in the peculiar case of $\text{WSe}_2$ where the *Q* valley is reached before the *K'* valley, the standard solutions persist for all values of $\alpha$ when *Q* is occupied but *K'* is not. This is consistent with *K'* (and $Q'$ when applicable) being the negative channel. Exotic solutions only emerge once *K*, *K'*, and *Q* are occupied.

## 3. Summary of the results

The main findings of this work with respect to the current understanding of superconductivity are twofold. First, while the onset of superconductivity at *Q* valley occupation was already anticipated [11,17,18,21], the underlying physical mechanism is unambiguously clarified and shown to be robust. Second, we propose a paradigm for the nature of
superconductivity, with a phonon-driven mechanism leading to exotic solutions in the presence of Coulomb repulsion.

The emergence of exotic solutions is an important result because their properties differ from the standard solutions, leading to different implications for experimental observations. The sign change in the order parameter while going between different valleys allows the system to take advantage of the intervalley Coulomb repulsion, which can dominate the phonon-mediated attraction for large enough $\alpha$. The overall superconducting state is still fully gapped; however, such a state with sign-changing order parameter is more fragile against atomic-scale disorder. That is because disorder can scatter electrons between regions with different gap signs, thereby reducing $T_c$ and introducing a finite density of states near the Fermi energy [66,67].

## IV. COMPARISON AND IMPLICATIONS FOR EXPERIMENTS

The results of the theoretical analysis presented above are in overall agreement with the main experimental results reported on different semiconducting TMDs, and provide clear indications to explain aspects of the experiments that are currently not understood. In this section we discuss experimental observations in the context of the framework provided by the theoretical model discussed in the paper.

A first essential point of our model is that it predicts that robust superconductivity occurs when the *Q* valleys start to be populated with electrons, and that electron-phonon interaction is the main driving force for the superconducting transition. Both the role of *Q*-valley filling and of electron-phonon interaction have been established by experiments over the last few years. Piatti et al. [17] first and Ding et al. [11] later have performed gate-dependent transport measurements (respectively on MoS₂ multilayers and WS₂ monolayers), showing that the onset for superconductivity indeed corresponds to filling of the *Q* valley. In independent Raman spectroscopy measurements on monolayers of different semiconducting TMDs performed as a function of electron density [21], it was further shown that populating the *Q* valley with electrons leads to a drastic increase in the strength of the electron-phonon coupling. This drastic increase, reproduced by our simulations and discussed here in Sec. III A, is indeed the mechanism driving superconductivity. We conclude that both the assumptions underlying our model, and the result of our calculations are fully consistent with experiments, and provide a microscopic explanation for the key aspects of the gate-induced superconducting transition that are common to the different TMDs investigated. Indeed, not only do our calculations show that the superconducting transition is brought into an experimentally accessible range of critical temperatures when the *Q* valley is filled, but also they give a range of critical temperatures for the different TMDs that matches the order of magnitude observed experimentally.

The second important conclusion of our model is that superconductivity, despite being driven by electron-phonon interaction, has a $s_{+-}$ character, i.e., an order parameter with different sign in the different populated bands and valleys, due to the non-negligible role of electron-electron interaction. We argue that the $s_{+-}$ nature of the superconducting state is essen-

tial to interpret different aspects of the experimental results that are currently not understood. That is because in a $s_{+-}$ state disorder has very different physical effects as compared to conventional superconducting states (even multiband ones), in which the order parameter has the same sign everywhere on the Fermi surface. The situation has been summarized in detail by Mazin and Schmalian [68], in their review of pnictide superconductors (many of which are also considered to host a $s_{+-}$ superconducting state) that we follow here to guide our discussion.

How disorder affects $s_{+-}$ superconductors depends strongly on the nature and the density of defects present in the system. Strong and short-range disorder potentials (such as those caused by atomic vacancies) introduce subgap states, without strongly affecting the superconducting transition temperature. This is theoretically predicted to happen up to rather large densities of defects, sufficient to suppress the “hard” gap of the pure material by filling it with states at all energies (all the way to $E=0$, i.e., at the Fermi energy). The situation is different for weak, long-range disorder potentials, such as those caused by inhomogeneity in a background charge distribution (e.g., as the one generated by spatial fluctuations in the density of dopants, in a common doped semiconductor). In $s_{+-}$ superconductors, such weak and long-range disorder potentials suppress the critical temperature and eventually kill superconductivity.

For gate-induced superconductivity in TMDs, both types of disorder are normally present. Chalcogen vacancies, which are unavoidably present in semiconducting TMDs with densities that depend on the specific compound, naturally lead to strong, short-range potentials. In MoS₂, for instance, it is common to have up to $10^{13} \text{ cm}^{-2}$ sulfur vacancies; in other TMDs the density of these vacancies is likely smaller, but still sizable. Based on the general behavior of $s_{+-}$ superconductors discussed above, these chalcogen vacancies create subgap states and can explain why a large density of states (DOS) at subgap energies has been observed in experiments on MoS₂, which would not be possible to explain in terms of a usual fully gapped s-wave superconducting state. To substantiate more this conclusion, it will be important in the future to measure the subgap DOS in other TMDs with lower density of chalcogen vacancies, or to find ways to correlate the subgap DOS with the density of these vacancies.

Weak long-range disorder, instead, originates from potential fluctuations associated to variations in the density of ions in the ionic liquid used for electrostatic gating. In a conventional s-wave superconductor, these potential fluctuations would have no important physical effect, but in a $s_{+-}$ they can suppress $T_c$ and weaken the superconducting state. We believe that this may be the reason for the superconducting dome that is seen in some experiments. Indeed, in some case (e.g., MoS₂), a dome with consistent properties has been reported in experiments by different groups, on identical kinds of devices. In other experiments, however, whether a dome is present or not seems to depend on the quality of the devices investigated (as measured, for instance, by the value of the electron mobility). In particular, in WS₂ monolayers, devices realized on hBN substrates exhibit larger $T_c$ and no dome (data are compatible with a saturation of $T_c$), whereas in devices with the WS₂ monolayer in direct con-
tact with a SiO₂ substrate the critical temperature is lower, and a clear dome-shaped dependence on electron density is observed. The latter experiments, as well as a comparison of reported data (see Fig. 1 in Ref. [43], the latter experiments [6,11]), suggest that the domelike dependence of $T_c$ on electron density is an extrinsic phenomenon. This is consistent with the predictions of our model, namely, that in the absence of disorder $T_c$ approximately saturates upon increasing electron density, once the Q valley starts to be populated, and that the dome may be understood as an intrinsic effect due to weak, long-range disorder due to the $s_{+-}$ nature of the superconducting state.

In summary, known experimental facts are consistent with all basic aspects of the behavior predicted by our model and the $s_{+-}$ nature of superconductivity provides a scenario that allows us to rationalize phenomena that were so far difficult to understand theoretically.

## V. CONCLUSION

We develop a model of superconductivity in semiconducting TMDs that is sophisticated enough to capture the momentum dependency of phonon-mediated and Coulomb interactions within a complex multivalley structure, yet simple enough to draw clear physical insight. We find systematic, robust results that apply to all the TMDs studied here, i.e., electron-doped MoS₂, WS₂, MoSe₂, and WSe₂. We find that superconductivity is phonon driven, with a doping-dependent onset corresponding to strong electron-phonon interactions being activated as the Q valley is occupied. In particular, zone border acoustic phonons mediate a strong coupling between K and Q valleys, while intravalley couplings are enhanced by a vanishing screening of the zone center $A_1$ mode. While phonons drive the emergence of superconductivity, Coulomb repulsion plays a crucial role as it is shown to change the nature of the superconducting solution. In particular, the introduction of Coulomb repulsion, even relatively weak, leads to exotic $s_{+-}$ solutions where the order parameter changes sign between valleys. All robust experimental facts currently reported in the literature are consistent with the results of our model, and the emergence of $s_{+-}$ solutions explains some aspects that were puzzling up to now. The physical insight gained on the nature and mechanism of superconductivity in TMDs can be used to better control it. For example, one could imagine tuning the Coulomb repulsion via dielectric engineering of the environment to cross over from standard to exotic superconducting solutions. The model is also general enough to draw conclusions on different systems. First, for hole doping of the same TMDs, we can expect similar qualitative conclusions with the Γ valley playing the role of the Q valleys. The fact that the Γ valley is much more difficult to reach with hole doping would then explain the lack of experimental evidence for superconductivity on the hole side. Many conclusions could apply to superconductivity in metallic TMDs such as NbSe₂ [69–71], where the Γ valley is instead typically populated and the importance of the multivalley and multiorbital character for superconductivity has already been reported [45,54,65,72]. Finally, a similar modeling procedure can be applied to other multivalley 2D materials, leveraging

the general strength of intervalley electron-phonon coupling to unveil new superconducting systems.

## ACKNOWLEDGMENTS

A.F.M. acknowledges financial support from the Swiss National Science Foundation under Project No. 200020_178891 and the EU Graphene Flagship project for support. M.G. acknowledges support from Ministero Italiano dell'Università della Ricerca through the PNRR project ECS_00000033_ECOSISTER and the PRIN2022 project SECSY. I.M. acknowledges support by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

## APPENDIX

### 1. Multivalley gap equation

In this Appendix we elaborate on our treatment of superconductivity. It is a simple generalization of the standard BCS treatment (see for instance [73]). The standard BCS gap equation is

$$ \Delta(k) = \sum_q W(q) \frac{\Delta(k-q)}{E_{k-q}} \tanh \frac{E_{k-q}}{2k_B T}. \quad (A1) $$

Here $E_k = \sqrt{\epsilon_k^2 + \Delta_k^2}$ is the Bogoliubov quasiparticle energy. The gap function $\Delta(k)$ corresponds to the pairing between spin-up electrons with momentum $k$, and spin-down electrons with momentum $-k$. The interaction $W(q)$ is assumed to be operational only within a window of energies $\pm\omega_D$ around the Fermi level. (In BCS, this is because the phonon-mediated interaction is attractive below the phonon characteristic phonon frequency, and decays rapidly above that frequency.)

If the gap function can be assumed to be momentum independent, $\Delta$, the equation for the superconducting transition temperature $T_c$ becomes

$$ \begin{aligned} \Delta &= \sum_q W(q) \frac{\Delta}{E_{k-q}} \tanh \frac{E_{k-q}}{2k_B T_c} \\ &\approx \Delta \ln \left[ \frac{\omega_D}{k_B T_c} \right] \left\langle \sum_q W(q) \delta(\epsilon_{k+q} - \epsilon_F) \right\rangle_{k \in \text{FS}} . \end{aligned} \quad (A2) $$

In our case, there are several symmetry-unrelated Fermi sheets, which in general will have different gap functions. For spin-up electrons, we label them K, K', Q, Q' (see Fig. 2). Thanks to the time-reversal symmetry, for spin-down electrons, the pockets have the same shape but are reflected as $k \to -k$. Note also that there are three symmetry-related Q and Q' Fermi surfaces, but only one K and K'. Given that all these Fermi sheets are much smaller than the size of the Brillouin zone, we can safely assume that within each Fermi pocket the gap function can be well approximated by just a single value $\Delta_v$.
These gap functions combine into a matrix gap equation that describes their couplings,

$$ \Delta_v \approx \sum_w \Delta_w \ln \left[ \frac{\omega_{vw}}{k_B T_c} \right] \left\langle \sum_{q,k+q \in w} W(q) \delta(\epsilon_{k+q} - \epsilon_F) \right\rangle_{k \in v} . \quad (A3) $$

Since there are three distinct symmetry-related Q and Q' pockets, it is conceivable that the order parameter could have a nonzero angular momentum. However, that would imply that the K and K' pockets remain nonsuperconducting. This is not energetically favorable; thus, we assumed that all Q (and Q') pockets have the same value of gap function. Consequently, there are four distinct values of the gap, which leads to the system of four equations, Eq. (1) and Fig. 2. Both the interacting strength and the attraction frequency cutoff for every pair of valleys are chosen based on the first-principles simulations, which we describe next.

## 2. First-principles simulations

### a. Electron-phonon interactions

The electron-phonon coupling matrix elements are obtained using the method described in Ref. [52], starting from the same ab initio simulations (the results of which were provided there). We report some aspects of this method here for completeness. Some repetition of the content of Ref. [52] is inevitable.

The ground-state and phonon properties of each semiconducting TMD (MoS₂, MoSe₂, WS₂, WSe₂) are simulated within density functional (perturbation) theory using the QUANTUM ESPRESSO package [74,75]. The full simulation is only performed at one carrier concentration of $n = 5 \times 10^{12} \text{ cm}^{-2}$, on the relaxed structures. We use norm-conserving, fully relativistic pseudopotentials with Perdew-Burke-Ernzerhof (PBE) functionals [76] from the pseudodoJO library [77], with kinetic energy cutoffs of 50 (WS₂), 70 (MoS₂), or 80 (WSe₂, MoSe₂). We use 2D boundary conditions and symmetric gates to induce doping [60]. The electronic momentum grid is nonuniform, with a sampling ranging from 12 × 12 to 96 × 96 around the Fermi surface. The electronic occupations follow the Fermi-Dirac distribution taken at room temperature, with a Fermi level such that the carrier density in the conduction band is $n = 5 \times 10^{12} \text{ cm}^{-2}$. The electron-phonon coupling matrix elements involve initial and final electronic momenta, assumed to be confined within a certain energy window, from the bottom of the conduction band up to around 0.3 eV above, as represented in Fig. 1. Final states are sampled on a 120 × 120 Monkhorst-Pack grid, while the initial states are sampled on a grid twice as coarse. The phonon momenta to compute are obtained by linking initial states to final states. The electron-phonon coupling matrix elements at other carrier densities are evaluated using the model of Ref. [52]. Notably, this accounts for the variations in free-carrier screening of the intravalley electron-phonon interactions. Indeed, as shown experimentally and explained qualitatively in Ref. [21], the occupation of both K and Q valleys leads to a suppressed screening and thus an enhancement of the intravalley electron-phonon couplings in TMDs. We then modeled and quantified this peculiar

screening mechanism as a function of Fermi level and valley occupations in Ref. [52]. Standard free-carrier screening is also accounted for. We effectively end up with a model of electron-phonon interactions as a function of both the Fermi level and the relative position of the Q (Q') valleys with respect to the K valleys.

Note that the position of the Q valley for MoS₂ in the calculations of Ref. [52] was fairly high (>0.2 eV). For convenience, and since we do not claim nor aim for accuracy in the relative positions of the valleys, the Q and Q' valleys of MoS₂ were lowered by 0.1 eV.

The phonons that significantly couple at the Fermi level are intravalley LA, TA, LO, and A₁, as well as intervalley zone border acoustic phonons. As the coupling between valleys of opposite spins is at least an order of magnitude smaller than same-spin couplings, it is neglected. On the one hand, standard free-carrier screening increases with doping, diminishing the strength of field-mediated, small-momenta (intravalley) couplings like Fröhlich and piezoelectric interactions. On the other hand, the screening of intravalley coupling to the A₁ and LA modes can counterintuitively decrease as a function of doping when both K and Q types of valleys are occupied. The couplings are then enhanced.

We do not compute the phonon softening coming from the occupation of the Q valley. The phonon frequencies are those obtained with only the K valleys occupied. It has been shown that the LA(M) phonon softens significantly when the corresponding electronic transitions are allowed, i.e., when the Fermi level crosses the Q valley. According to anharmonic phonon calculations in multilayer MoS₂ [43], we can expect the LA mode to soften at M. The frequency decreases as a function of carrier concentration and potentially vanishes, leading to a charge density wave. The computation of a realistic softening value is challenging and out of scope here. Our objective here is more the description of the nature of the solutions and the understanding of the mechanisms behind. We simply apply a 30% softening of the LA(M) mode when Q is occupied.

Equation (2) describes the couplings (phonon-mediated attraction or Coulomb repulsion) between different valleys. The final states are summed on the arrival valley w. In principle, one would average on the initial states at the Fermi level in v. Figure 9 shows the matrix elements $G_{vw}$ and $C_{vw}$ for all initial states sampled in the irreducible wedge of the Brillouin zone. Plotted as a function of the Fermi level, we obtain well-defined lines of minimal widths. This indicates that those quantities are mostly dependent on the energy of the initial state and not so much on its momentum. This justifies the procedure of simply picking an initial state at the Fermi level rather than averaging over all states at the Fermi level.

### b. Screened Coulomb interaction

To evaluate the dielectric function of the doped system, we simply sum the irreducible response functions associated to the neutral material $\chi_{\text{neut}}^0$ and the added electrons in the conduction band $\chi_{\text{cond}}^0$. The former is extracted from DFPT on the neutral material as in Refs. [78,79]. The latter is computed
<figure><img src="image_11.png" /><figcaption>FIG. 9. Matrix elements of phonon-mediated attraction and Coulomb repulsion for all initial states $k_i$ of energy $\epsilon_i$ in the irreducible wedge of the Brillouin zone, showing negligible angular variations of the matrix elements.</figcaption></figure>

from the band structure as

$$ \chi_{\text{cond}}^{0}(q) = \sum_{n} \int \frac{d^2\mathbf{k}}{(2\pi)^2} \frac{n_{\epsilon_{nk}}^{\text{FD}} - n_{\epsilon_{nk+q}}^{\text{FD}}}{\epsilon_{nk} - \epsilon_{nk+q}} |(u_{nk}|u_{nk+q})|^2, \quad (\text{A4}) $$

where $n_{\epsilon_{nk}}^{\text{FD}}$ is the Fermi-Dirac occupation for electronic state $\mathbf{k}$ of energy $\epsilon_{nk}$. The wave-function overlaps are not computed. In general, it is 1 for $\mathbf{q} \to 0$ (same electronic state), and

decreases as $|\mathbf{q}|$ increases. It is expected to be much smaller than 1 for intervalley transitions, e.g., when $\mathbf{q} \sim \mathbf{M}$, such that $\chi_{\text{cond}}^{0}(q)$ becomes smaller for $|\mathbf{q}|$ larger than the size of the valleys. We assume the contribution from $\chi_{\text{cond}}^{0}(q)$ to be negligible with respect to $\chi_{\text{neut}}^{0}$, and use a simple complementary error function to model the wave-function overlaps, with a threshold roughly corresponding to the size of the valleys: 0.3 bohrs$^{-1}$.

[1] J. T. Ye, Y. J. Zhang, R. Akashi, M. S. Bahramy, R. Arita, and Y. Iwasa, Superconducting dome in a gate-tuned band insulator, Science 338, 1193 (2012).

[2] K. Taniguchi, A. Matsumoto, H. Shimotani, and H. Takagi, Electric-field-induced superconductivity at 9.4 K in a layered transition metal disulphide MoS$_2$, Appl. Phys. Lett. 101, 042603 (2012).

[3] W. Shi, J. Ye, Y. Zhang, R. Suzuki, M. Yoshida, J. Miyazaki, N. Inoue, Y. Saito, and Y. Iwasa, Superconductivity series in transition metal dichalcogenides by ionic gating, Sci. Rep. 5, 12534 (2015).

[4] S. Jo, D. Costanzo, H. Berger, and A. F. Morpurgo, Electrostatically induced superconductivity at the surface of WS$_2$, Nano Lett. 15, 1197 (2015).

[5] J. M. Lu, O. Zheliuk, I. Leermakers, N. F. Q. Yuan, U. Zeitler, K. T. Law, and J. T. Ye, Evidence for two-dimensional Ising superconductivity in gated MoS$_2$, Science 350, 1353 (2015).

[6] J. Lu, O. Zheliuk, Q. Chen, I. Leermakers, N. E. Hussey, U. Zeitler, and J. Ye, Full superconducting dome of strong Ising protection in gated monolayer WS$_2$, Proc. Natl. Acad. Sci. USA 115, 3551 (2018).

[7] D. Costanzo, S. Jo, H. Berger, and A. F. Morpurgo, Gate-induced superconductivity in atomically thin MoS$_2$ crystals, Nat. Nanotechnol. 11, 339 (2016).

[8] D. Qiu, C. Gong, S. Wang, M. Zhang, C. Yang, X. Wang, and J. Xiong, Recent advances in 2D superconductors, Adv. Mater. 33, 2006124 (2021).

[9] Y. Fu, E. Liu, H. Yuan, P. Tang, B. Lian, G. Xu, J. Zeng, Z. Chen, Y. Wang, W. Zhou, K. Xu, A. Gao, C. Pan, M. Wang, B. Wang, S.-C. Zhang, Y. Cui, H. Y. Hwang, and F. Miao, Gated tuned superconductivity and phonon softening in monolayer and bilayer MoS$_2$, npj Quantum Mater. 2, 52 (2017).

[10] O. Zheliuk, J. M. Lu, Q. H. Chen, A. A. E. Yumin, S. Golightly, and J. T. Ye, Josephson coupled Ising pairing induced in suspended MoS$_2$ bilayers by double-side ionic gating, Nat. Nanotechnol. 14, 1123 (2019).

[11] D. Ding, Z. Qu, X. Han, C. Han, Q. Zhuang, X.-L. Yu, R. Niu, Z. Wang, Z. Li, Z. Gan, J. Wu, and J. Lu, Multivalley superconductivity in monolayer transition metal dichalcogenides, Nano Lett. 22, 7919 (2022).

[12] Y. Saito, Y. Nakamura, M. S. Bahramy, Y. Kohama, J. Ye, Y. Kasahara, Y. Nakagawa, M. Onga, M. Tokunaga, T. Nojima, Y. Yanase, and Y. Iwasa, Superconductivity protected by spin-valley locking in ion-gated MoS$_2$, Nat. Phys. 12, 144 (2016).

[13] Z. Y. Zhu, Y. C. Cheng, and U. Schwingenschlögl, Giant spin-orbit-induced spin splitting in two-dimensional transition-metal dichalcogenide semiconductors, Phys. Rev. B 84, 153402 (2011).

[14] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of MoS$_2$ and other group-VI dichalcogenides, Phys. Rev. Lett. 108, 196802 (2012).

[15] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, 2D transition metal dichalcogenides, Nat. Rev. Mater. 2, 17033 (2017).

[16] B. T. Zhou, N. F. Q. Yuan, H.-L. Jiang, and K. T. Law, Ising superconductivity and Majorana fermions in transition-metal dichalcogenides, Phys. Rev. B 93, 180501(R) (2016).

[17] E. Piatti, D. De Fazio, D. Daghero, S. R. Tamalampudi, D. Yoon, A. C. Ferrari, and R. S. Gonnelli, Multi-valley superconductivity in ion-gated MoS$_2$ layers, Nano Lett. 18, 4821 (2018).

[18] E. Piatti, D. Romanin, and R. S. Gonnelli, Mapping multi-valley lifshitz transitions induced by field-effect doping in strained MoS$_2$ nanolayers, J. Phys.: Condens. Matter 31, 114002 (2019).

[19] R. Roldán, J. A. Silva-Guillén, M. P. López-Sancho, F. Guinea, E. Cappelluti, and P. Ordejón, Electronic properties of single-
layer and multilayer transition metal dichalcogenides MX$_2$ (M = Mo, W and X = S, Se), Ann. Phys. **526**, 347 (2014).

[20] H. Yuan, Z. Liu, G. Xu, B. Zhou, S. Wu, D. Dumcenco, K. Yan, Y. Zhang, S.-K. Mo, P. Dudin, V. Kandyba, M. Yablonskikh, A. Barinov, Z. Shen, S. Zhang, Y. Huang, X. Xu, Z. Hussain, H. Y. Hwang, Y. Cui et al., Evolution of the valley position in bulk transition-metal chalcogenides and their monolayer limit, Nano Lett. **16**, 4738 (2016).

[21] T. Sohier, E. Ponomarev, M. Gibertini, H. Berger, N. Marzari, N. Ubrig, and A. F. Morpurgo, Enhanced electron-phonon interaction in multivalley materials, Phys. Rev. X **9**, 031019 (2019).

[22] D. Costanzo, H. Zhang, B. A. Reddy, H. Berger, and A. F. Morpurgo, Tunnelling spectroscopy of gate-induced superconductivity in MoS$_2$, Nat. Nanotechnol. **13**, 483 (2018).

[23] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of superconductivity, Phys. Rev. **108**, 1175 (1957).

[24] T. Li and G. Galli, Electronic properties of MoS$_2$ nanoparticles, J. Phys. Chem. C **111**, 16192 (2007).

[25] E. S. Kadantsev and P. Hawrylak, Electronic structure of a single MoS$_2$ monolayer, Solid State Commun. **152**, 909 (2012).

[26] A. Kormányos, G. Burkard, M. Gmitra, J. Fabian, V. Zólyomi, N. D. Drummond, and V. Fal'ko, **k·p** theory for two-dimensional transition metal dichalcogenide semiconductors, 2D Mater. **2**, 022001 (2015).

[27] T. Brumme, M. Calandra, and F. Mauri, First-principles theory of field-effect doping in transition-metal dichalcogenides: Structural properties, electronic structure, Hall coefficient, and electrical conductivity, Phys. Rev. B **91**, 155436 (2015).

[28] A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. Chim, G. Galli, and F. Wang, Emerging photoluminescence in monolayer MoS$_2$, Nano Lett. **10**, 1271 (2010).

[29] Y. Wang, T. Sohier, K. Watanabe, T. Taniguchi, M. J. Verstraete, and E. Tutuc, Electron mobility in monolayer WS$_2$ encapsulated in hexagonal boron-nitride, Appl. Phys. Lett. **118**, 102105 (2021).

[30] Y. Ge and A. Y. Liu, Phonon-mediated superconductivity in electron-doped single-layer MoS$_2$: A first-principles prediction, Phys. Rev. B **87**, 241408(R) (2013).

[31] M. Rösner, S. Haas, and T. O. Wehling, Phase diagram of electron-doped dichalcogenides, Phys. Rev. B **90**, 245105 (2014).

[32] R. Roldán, E. Cappelluti, and F. Guinea, Interactions and superconductivity in heavily doped MoS$_2$, Phys. Rev. B **88**, 054515 (2013).

[33] N. F. Q. Yuan, K. F. Mak, and K. T. Law, Possible topological superconducting phases of MoS$_2$, Phys. Rev. Lett. **113**, 097001 (2014).

[34] Y.-T. Hsu, A. Vaezi, M. H. Fischer, and E.-A. Kim, Topological superconductivity in monolayer transition metal dichalcogenides, Nat. Commun. **8**, 14985 (2017).

[35] R. Oiwa, Y. Yanagi, and H. Kusunose, Theory of superconductivity in hole-doped monolayer MoS$_2$, Phys. Rev. B **98**, 064509 (2018).

[36] T. Das and K. Dolui, Superconducting dome in MoS$_2$ and TiSe$_2$ generated by quasiparticle-phonon coupling, Phys. Rev. B **91**, 094510 (2015).

[37] J. Wang, X. Zhang, R. Ma, G. Yang, E. V. Castro, and T. Ma, Spin-triplet superconducting pairing in doped MoS$_2$, Phys. Rev. B **106**, 134513 (2022).

[38] G. Schönhoff, M. Rösner, R. E. Groenewald, S. Haas, and T. O. Wehling, Interplay of screening and superconductivity in low-dimensional materials, Phys. Rev. B **94**, 134504 (2016).

[39] E. R. Margine and F. Giustino, Anisotropic Migdal-Eliashberg theory using Wannier functions, Phys. Rev. B **87**, 024505 (2013).

[40] A. B. Migdal, Interaction between electrons and lattice vibrations in a normal metal, Zh. Eksp. Teor. Fiz., **34**, 1438 (1958) [Sov. Phys.-JETP **7**, 996 (1958)].

[41] G. M. Eliashberg, Interactions between electrons and lattice vibrations in a superconductor, Zh. Eksp. Teor. Fiz. **34**, 1438 (1958) [Sov. Phys.-JETP **11**, 696 (1960)].

[42] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, The stochastic self-consistent harmonic approximation: Calculating vibrational properties of materials with full quantum and anharmonic effects, J. Phys.: Condens. Matter **33**, 363001 (2021).

[43] G. Marini and M. Calandra, Phonon mediated superconductivity in field-effect doped molybdenum dichalcogenides, 2D Mater. **10**, 015013 (2023).

[44] M. Simonato, M. I. Katsnelson, and M. Rösner, Revised Tolmachev-Morel-Anderson pseudopotential for layered conventional superconductors with nonlocal coulomb interaction, Phys. Rev. B **108**, 064513 (2023).

[45] S. Das, H. Paudyal, E. R. Margine, D. F. Agterberg, and I. I. Mazin, Electron-phonon coupling and spin fluctuations in the Ising superconductor NbSe$_2$, npj Comput. Mater. **9**, 66 (2023).

[46] X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Spin and pseudospins in layered transition metal dichalcogenides, Nat. Phys. **10**, 343 (2014).

[47] Z. Wang, A. Molina-Sánchez, P. Altmann, D. Sangalli, D. De Fazio, G. Soavi, U. Sassi, F. Bottegoni, F. Ciccacci, M. Finazzi, L. Wirtz, A. C. Ferrari, A. Marini, G. Cerullo, and S. Dal Conte, Intravalley spin-flip relaxation dynamics in single-layer WS$_2$, Nano Lett. **18**, 6882 (2018).

[48] L. P. Gor'kov and E. I. Rashba, Superconducting 2D system with Lifted spin degeneracy: Mixed singlet-triplet state, Phys. Rev. Lett. **87**, 037004 (2001).

[49] H. Zhang, C. Berthod, H. Berger, T. Giamarchi, and A. F. Morpurgo, Band filling and cross quantum capacitance in ion-gated semiconducting transition metal dichalcogenide monolayers, Nano Lett. **19**, 8836 (2019).

[50] C. H. Ahn, A. Bhattacharya, M. Di Ventra, J. N. Eckstein, C. D. Frisbie, M. E. Gershenson, A. M. Goldman, I. H. Inoue, J. Mannhart, A. J. Millis, A. F. Morpurgo, D. Natelson, and J.-M. Triscone, Electrostatic modification of novel materials, Rev. Mod. Phys. **78**, 1185 (2006).

[51] C. Cao, M. Melegari, M. Philippi, D. Domaretskiy, N. Ubrig, I. Gutiérrez-Lezama, and A. F. Morpurgo, Full control of solid-state electrolytes for electrostatic gating, Adv. Mater. **35**, 2211993 (2023).

[52] T. Sohier, P. M. M. C. de Melo, Z. Zanolli, and M. J. Verstraete, The impact of valley profile on the mobility and Kerr rotation of transition metal dichalcogenides, 2D Mater. **10**, 025006 (2023).

[53] M. Tinkham, *Introduction to Superconductivity*, Dover Books on Physics Series (Dover, New York, 2004).
[54] S. Hörhold, J. Graf, M. Marganska, and M. Grifoni, Two-bands Ising superconductivity from Coulomb interactions in monolayer NbSe₂, *2D Mater.* **10**, 025008 (2023).

[55] R. Lucrezi, P. P. Ferreira, S. Hajinazar, H. Mori, H. Paudyal, E. R. Margine, and C. Heil, Full-bandwidth anisotropic Migdal-Eliashberg theory and its application to superhydrides, *Commun. Phys.* **7**, 33 (2024).

[56] P. Morel and P. W. Anderson, Calculation of the superconducting state parameters with retarded electron-phonon interaction, *Phys. Rev.* **125**, 1263 (1962).

[57] S. Zeng, Y. Zhao, G. Li, and J. Ni, Strongly enhanced superconductivity in doped monolayer MoS₂ by strain, *Phys. Rev. B* **94**, 024501 (2016).

[58] C. Heil, S. Poncé, H. Lambert, M. Schlipf, E. R. Margine, and F. Giustino, Origin of superconductivity and latent charge density wave in NbS₂, *Phys. Rev. Lett.* **119**, 087003 (2017).

[59] H. Paudyal, S. Poncé, F. Giustino, and E. R. Margine, Superconducting properties of MoTe₂ from Ab Initio anisotropic Migdal-Eliashberg theory, *Phys. Rev. B* **101**, 214515 (2020).

[60] T. Sohier, M. Calandra, and F. Mauri, Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene, *Phys. Rev. B* **96**, 075448 (2017).

[61] D. Campi, N. Mounet, M. Gibertini, G. Pizzi, and N. Marzari, The Materials Cloud 2D database (MC2D), www материcloud.org/discover/mc2d/details/MoS2-MoS2.

[62] E. Cappelluti, R. Roldán, J. A. Silva-Guillén, P. Ordejón, and F. Guinea, Tight-binding model and direct-gap/indirect-gap transition in single-layer and multilayer MoS₂, *Phys. Rev. B* **88**, 075409 (2013).

[63] J. A. Ashraff and P. D. Loly, The triangular linear analytic method for two-dimensional spectral functions, *J. Phys. C: Solid State Phys.* **20**, 4823 (1987).

[64] P. W. Anderson, Theory of dirty superconductors, *J. Phys. Chem. Solids* **11**, 26 (1959).

[65] D. Möckli and M. Khodas, Magnetic-field induced s + if pairing in Ising superconductors, *Phys. Rev. B* **99**, 180505(R) (2019).

[66] I. Mazin, D. J. Singh, s. M. Johannes, and M.-H. Du, Unconventional superconductivity with a sign reversal in the order parameter of LaFeAsO₁₋ₓFₓ, *Phys. Rev. Lett.* **101**, 057003 (2008).

[67] R. M. Fernandes, A. I. Coldea, H. Ding, I. R. Fisher, P. J. Hirschfeld, and G. Kotliar, Iron pnictides and chalcogenides: A new paradigm for superconductivity, *Nature (London)* **601**, 35 (2022).

[68] I. Mazin and J. Schmalian, Pairing symmetry and pairing state in ferropnictides: Theoretical overview, *Phys. C (Amsterdam)* **469**, 614 (2009).

[69] X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Ising pairing in superconducting NbSe₂ atomic layers, *Nat. Phys.* **12**, 139 (2016).

[70] A. W. Tsen, B. Hunt, Y. D. Kim, Z. J. Yuan, S. Jia, R. J. Cava, J. Hone, P. Kim, C. R. Dean, and A. N. Pasupathy, Nature of the quantum metal in a two-dimensional crystalline superconductor, *Nat. Phys.* **12**, 208 (2016).

[71] M. M. Ugeda, A. J. Bradley, Y. Zhang, S. Onishi, Y. Chen, W. Ruan, C. Ojeda-Aristizabal, H. Ryu, M. T. Edmonds, H.-Z. Tsai, A. Riss, S.-K. Mo, D. Lee, A. Zettl, Z. Hussain, Z.-X. Shen, and M. F. Crommie, Characterization of collective ground states in single-layer NbSe₂, *Nat. Phys.* **12**, 92 (2016).

[72] D. Möckli and M. Khodas, Robust parity-mixed superconductivity in disordered monolayer transition metal dichalcogenides, *Phys. Rev. B* **98**, 144518 (2018).

[73] G. D. Mahan, *Many-Particle Physics* (Springer, New York, 2013).

[74] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. De Gironcoli, P. Delugas, R. A. Distasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo *et al.*, Advanced capabilities for materials modelling with QUANTUM ESPRESSO, *J. Phys.: Condens. Matter* **29**, 465901 (2017).

[75] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, Davide Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, Anton Kokalj, M. Lazzeri, L. Martin-Samos *et al.*, QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials, *J. Phys.: Condens. Matter* **21**, 395502 (2009).

[76] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, *Phys. Rev. B* **88**, 085117 (2013).

[77] M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, The PSEUDODOJO: Training and grading a 85 element optimized norm-conserving pseudopotential table, *Comput. Phys. Commun.* **226**, 39 (2018).

[78] T. Sohier, M. Gibertini, and M. J. Verstraete, Remote free-carrier screening to boost the mobility of Fröhlich-limited two-dimensional semiconductors, *Phys. Rev. Mater.* **5**, 024004 (2021).

[79] F. Macheda, F. Mauri, and T. Sohier, Ab Initio Van der waals electrodynamics: Polaritons and electron scattering from plasmons and phonons in BN-capped Graphene, *Phys. Rev. B* **110**, 115407 (2024).