Enhancement of the upper critical field in disordered transition metal dichalcogenide monolayers

TL;DR

This paper investigates how impurities affect the superconducting phase diagram of transition metal dichalcogenide monolayers under in-plane magnetic fields, revealing that strong spin-orbit coupling enhances the upper critical field beyond the Pauli limit, with impurity effects depending on scattering types.

Contribution

It provides a detailed analysis of impurity effects on the upper critical field in disordered transition metal dichalcogenide monolayers, highlighting the roles of intravalley and intervalley scattering.

Findings

Upper critical field exceeds Pauli limit at low temperatures due to spin-orbit coupling.

Insensitivity of critical field to intravalley scattering.

Critical field limited by intervalley scattering.

Abstract

We calculate the effect of impurities on the superconducting phase diagram of transition metal dichalcogenide monolayers in the presence of an in-plane magnetic field. Due to strong intrinsic spin-orbit coupling, the upper critical field greatly surpasses the Pauli limit at low temperatures. We find that it is insensitive to intravalley scattering and, ultimately, limited by intervalley scattering.

Figures

Equations

$${\mathcal{H}}_{\mathbf{q}}$$

$${\mathcal{H}}_{\eta \mathbf{q}}={\xi }_{\mathbf{q}}+\eta {\Delta }_{so}(\mathbf{q})s_z+h{s}_x,$$

$$H_{BCS}=\sum _{\eta \mathbf{q}}{c}_{\eta \mathbf{q}}^{†}{\mathcal{H}}_{\eta \mathbf{q}}{c}_{\eta \mathbf{q}}+\Delta \sum _{\eta \mathbf{q}}{c}_{\eta \mathbf{q}\uparrow }^{†}{c}_{\bar{\eta }\bar{\mathbf{q}}\downarrow }^{†}+\text{h.c.},$$

$$\Delta =\frac{\lambda T}{4}\sum _{\eta ,\mathbf{q},|\omega |<\Omega }\text{Tr}[i{s}_y{\mathcal{G}}_{\bar{\eta }\bar{\mathbf{q}}\omega }^{+}\Delta i{s}_y{\mathcal{G}}_{\eta \mathbf{q}\omega }^{-}],$$

$${\mathcal{H}}_D(\mathbf{q}-{\mathbf{q}}^{\prime })=V_1(\mathbf{q}-{\mathbf{q}}^{\prime })+\sum _{i=x,y}{V}_{2i}(\mathbf{q}-{\mathbf{q}}^{\prime }){\eta }_i.$$

$$H_D=\sum _{\eta \mathbf{q}{\mathbf{q}}^{\prime }}{\mathcal{V}}_{1\mathbf{q}{\mathbf{q}}^{\prime }}{c}_{\eta \mathbf{q}}^{†}{c}_{\eta {\mathbf{q}}^{\prime }}+{\mathcal{V}}_{2\mathbf{q}{\mathbf{q}}^{\prime }}{c}_{\eta \mathbf{q}}^{†}{c}_{\bar{\eta }{\mathbf{q}}^{\prime }}+\text{h.c.},$$

$$\ln \frac{T}{T_c}$$

$$=2\pi T\sum _{\omega >0}\left[\frac{\omega (\omega +\frac{1}{{\tau }_2})+{\Delta }_{so}^2}{(\omega +\frac{1}{{\tau }_2})({\omega }^2+h_{c2}^2+{\Delta }_{so}^2)-\frac{{\Delta }_{so}^2}{{\tau }_2}}-\frac{1}{\omega }\right],$$

$$\ln\frac{T_{c}}{T}=\frac{h_{c2}^{2}}{h_{c2}^{2}+\Delta_{so}^{2}}\Re\bigg{[}\psi\bigg{(}\frac{1}{2}+\frac{i\sqrt{h_{c2}^{2}+\Delta_{so}^{2}}}{2\pi T}\bigg{)}-\psi\bigg{(}\frac{1}{2}\bigg{)}\bigg{]},$$

$$h_{c2}\propto\Delta_{so}\sqrt{\ln\frac{T_{c}}{T}}\qquad\text{for}\qquad\frac{T}{T_{c}}\ll\exp\bigg{(}-c\frac{\Delta_{0}^{2}}{\Delta_{so}^{2}}\bigg{)},$$

$$h_{c2}\approx {\Delta }_{so}\sqrt{\ln \frac{T_c}{T}/\ln \frac{2{\Delta }_{so}}{{\Delta }_0}}\text{for}\frac{T}{T_c}\ll \frac{{\Delta }_0}{{\Delta }_{so}}.$$

$$h_{c2}\approx {\Delta }_{so}\frac{1}{\sqrt{\ln \frac{2{\Delta }_{so}}{{\Delta }_0}}}\sqrt{1-\frac{T}{T_c}}.$$

$$h_{c2}\approx {\Delta }_{so}\sqrt{\ln ({\Delta }_0{\tau }_2)/\ln \frac{2{\Delta }_{so}}{{\Delta }_0}}$$

$$1=\frac{\lambda}{4}T\int\frac{d^{2}\,\mathbf{q}}{(2\pi)^{2}}\sum_{\eta\omega ss^{\prime}s^{\prime\prime}}\ s\langle\mathcal{G}_{\eta\mathbf{q}\omega}^{+,ss^{\prime}}\rangle\langle\mathcal{G}_{\bar{\eta}\mathbf{\bar{q}}\omega}^{-,\bar{s}s^{\prime\prime}}\rangle\bigg{[}s^{\prime}\delta_{\bar{s}^{\prime}s^{\prime\prime}}+\Pi_{\eta}^{s^{\prime}s^{\prime\prime}}(\theta)\bigg{]}.$$

$$\nu \int d{\xi }_{\mathbf{q}}⟨{\mathcal{G}}_{\eta \mathbf{q}\omega }^{+,s_1{s}_2}⟩⟨{\mathcal{G}}_{\bar{\eta }\bar{\mathbf{q}}\omega }^{-,s_3{s}_4}⟩$$

$$\displaystyle\Pi_{\eta}^{s_{1}s_{2}}(\theta)=\int\frac{d^{2}\,\mathbf{q^{\prime}}}{(2\pi)^{2}}\sum_{ss^{\prime}}\bigg{[}$$

$${\Pi }_{\eta }^{s_1{s}_2}(\theta )={\Pi }_{\eta 0}^{s_1{s}_2}+{\Pi }_{\eta 1}^{s_1{s}_2}\cos \theta +{\Pi }_{\eta 2}^{s_1{s}_2}\sin \theta .$$

$$\displaystyle\Pi_{\eta 0}^{s_{1}s_{2}}=\int\frac{d\,\xi_{\mathbf{q}^{\prime}}}{2\pi}\sum_{ss^{\prime}}\bigg{[}\frac{1}{\tau_{1}}\langle\mathcal{G}_{\eta\mathbf{q^{\prime}}\omega}^{+,s_{1}s}\rangle$$

$$\ln \frac{T_c}{T}=\frac{h_{c2}^2}{{\rho }^2}\ln \frac{4e^{\gamma }\rho }{2\pi T},$$

$$\ln \frac{2h_{c2}}{{\Delta }_0}\approx \frac{{\Delta }_{so}^2}{h_{c2}^2}\ln \frac{T_c}{T}.$$

$$\ln\frac{T_{c}}{T}=\frac{\Delta_{so}^{2}}{\rho^{2}}\bigg{[}\psi\bigg{(}\frac{1}{2}+\frac{h_{c2}^{2}}{2\pi\tau_{2}\rho^{2}T}\bigg{)}-\psi\bigg{(}\frac{1}{2}\bigg{)}\bigg{]}+\frac{h_{c2}^{2}}{\rho^{2}}\Re\bigg{[}\psi\bigg{(}\frac{1}{2}+\frac{i\rho}{2\pi T}\bigg{)}-\psi\bigg{(}\frac{1}{2}\bigg{)}\bigg{]}.$$

$$\ln \frac{2h_{c2}}{{\Delta }_0}\approx \frac{{\Delta }_{so}^2}{h_{c2}^2}\ln {\tau }_2{h}_{c2},$$

$$\ln\frac{T}{T_{c}}=2\pi T\sum_{\omega>0}\bigg{[}\frac{\omega+\tau_{2}\Delta_{so}^{2}}{\omega^{2}+h_{c2}^{2}}-\frac{1}{\omega}\bigg{]}=\Re\bigg{[}\psi\bigg{(}\frac{1}{2}\bigg{)}-\psi\bigg{(}\frac{1}{2}+\frac{ih_{c2}}{2\pi T}\bigg{)}\bigg{]}+\frac{\pi\Delta_{so}^{2}\tau_{2}}{2h_{c2}}\tanh\frac{h_{c2}}{2T},$$

$$\ln \frac{2h_{c2}}{{\Delta }_0}=\frac{\pi {\Delta }_{so}^2{\tau }_2}{2h_{c2}}.$$

$$\ln\frac{2h_{c2}(0)}{\Delta_{0}}=\int_{0}^{\infty}d\omega\bigg{[}\frac{\omega(\omega+\frac{1}{\tau_{2}})+\Delta_{so}^{2}}{(\omega+\frac{1}{\tau_{2}})(\omega^{2}+h_{c2}^{2}+\Delta_{so}^{2})-\frac{\Delta_{so}^{2}}{\tau_{2}}}-\frac{\omega}{\omega^{2}+h_{c2}^{2}}\bigg{]}.$$

$$h_{c2}\simeq \sqrt{\frac{T_c(T_c-T)}{C}}$$

$$C=2\pi T_c^3\sum _{\omega >0}\frac{\omega +\frac{1}{{\tau }_2}}{{\omega }^2[\omega (\omega +\frac{1}{{\tau }_2})+{\Delta }_{so}^2]}$$

$$H_{FFLO}=\sum _{\eta \mathbf{q}}{c}_{\eta \mathbf{q}}^{†}{\mathcal{H}}_{\eta \mathbf{q}}{c}_{\eta \mathbf{q}}+{\Delta }_{\mathbf{p}}\sum _{\tau \mathbf{q}}{c}_{\tau \mathbf{q}+\mathbf{p}/2\uparrow }^{†}{c}_{\bar{\tau }\bar{\mathbf{q}}+\mathbf{p}/2\downarrow }^{†}+h.c..$$

$${\Delta }_{\mathbf{p}}=\frac{\lambda T}{4}\sum _{\eta ,\mathbf{q},|\omega |<\Omega }\text{Tr}[i{s}_y{\mathcal{G}}_{\bar{\eta }\bar{\mathbf{q}}+\frac{\mathbf{p}}{2}\omega }^{+}{\Delta }_{\mathbf{p}}i{s}_y{\mathcal{G}}_{\eta \mathbf{q}+\frac{\mathbf{p}}{2}\omega }^{-}].$$

$$\ln\frac{T}{T_{c}}=\int\frac{d\theta}{2\pi}\Re\bigg{\{}\frac{\Delta_{so}^{2}}{\rho^{2}}\bigg{[}\psi\bigg{(}\frac{1}{2}\bigg{)}-\psi\bigg{(}\frac{1}{2}-\frac{iv_{F}\mathbf{p}.\hat{\mathbf{q}}}{4\pi T}\bigg{)}\bigg{]}+\frac{h_{c2}^{2}}{\rho^{2}}\bigg{[}\psi\bigg{(}\frac{1}{2}\bigg{)}-\psi\bigg{(}\frac{1}{2}+\frac{2i\rho-iv_{F}\mathbf{p}.\hat{\mathbf{q}}}{4\pi T}\bigg{)}\bigg{]}\bigg{\}}.$$

$$\ln\frac{T}{T_{c}}=-\frac{h_{c2}^{2}}{\rho^{2}}\Re\bigg{[}\psi\bigg{(}\frac{1}{2}+\frac{i\rho}{2\pi T}\bigg{)}-\psi\bigg{(}\frac{1}{2}\bigg{)}\bigg{]}+\frac{v_{F}^{2}|\mathbf{p}|^{2}T_{c}^{2}}{16\pi^{2}h_{c2}^{2}T^{2}}F_{1}(T,\Delta_{so}),$$

Full text

Due to the heavy constituent atoms, TMDC monolayers exhibit a very large intrinsic spin-orbit coupling (SOC), often called Ising SOC Zhou et al. (2016), which acts as an effective Zeeman field perpendicular to the plane of the material with opposite orientations in the two valleys Zhu et al. (2011); Kormányos et al. (2015); Xiao et al. (2012). As a consequence, a large valley-dependent spin-splitting occurs in the valence band as well as in the conduction band, though with a much smaller magnitude (due to the predominant $d_{x^{2}-y^{2}}\pm id_{xy}$ and $d_{z^{2}}$ orbital character of carriers, respectively). Recent optical investigations have confirmed that electrons from different valleys can be excited selectively with circularly polarized light Mak et al. (2012); Zeng et al. (2012). These properties open the door for novel applications in spintronics and so-called valleytronics Wang et al. (2012b); Xiao et al. (2012).

The coupling between the spin and valley degrees of freedom has remarkable repercussions for the superconducting properties that have been reported in these materials. In heavily $n$-doped ionic-gated $\text{MoS}_{2}$ flakes, a 2D superconducting phase with a critical temperature $T_{c}$ around 10 K has been observed Saito et al. (2016); Lu et al. (2015). Interestingly, the in-plane upper critical field {\color[rgb]{0,0,0}H_{c2}} reaches up to 60 T, thus greatly surpassing the Pauli limit, $H_{P}=\sqrt{2}{\Delta_{0}}/(g\mu_{B})$, where $\Delta_{0}\simeq 1.76k_{B}T_{c}$, $\mu_{B}$ is the Bohr magneton, and $g$ the $g$-factor. This is interpreted as a consequence of the interplay of Ising SOC and the 2D nature of the materials. Namely, in 2D materials the orbital pair-breaking effect is largely suppressed for an in-plane magnetic field Tinkham (1996). On the other hand, due to the Ising SOC, the in-plane magnetic field is also not efficient for breaking Cooper pairs by the paramagnetic effect, as they are formed of electrons in opposite valleys with strongly pinned out-of-plane spins. Related results have also been reported in superconducting $\text{NbSe}_{2}$ monolayers in the $p$-doped regime Xi et al. (2016). Existing theories only considered the clean case Saito et al. (2016); Frigeri et al. (2004) and do not describe all the experimental observations.

In this work, we establish the theory of the upper critical field for Ising superconductors at arbitrary disorder strength, assuming a conventional $s$-wave Cooper pairing. The enhancement of the in-plane upper critical field above the Pauli limit is a general feature of superconductivity in the presence of SOC. It has been predicted in disordered 2D superconductors with spin-orbit scattering Maki (1966); Klemm et al. (1975) or Rashba SOC Barzykin and Gor’kov (2002). The latter has been invoked to interpret recent experiments on oxide interfaces Reyren et al. (2009) and Pb monolayer films Sekihara et al. (2013). Ising SOC leads to qualitatively different effects. Using a simple model, we show that Ising superconductivity results in a much larger enhancement and exhibits several new and interesting properties. Namely, $H_{c2}$ diverges at low temperature – a phenomenon which is robust to intravalley scattering. On the other hand, intervalley scattering provides an effective spin-flip scattering mechanism and leads to the saturation of $H_{c2}$, consistent with the experimental findings. Furthermore, we show that, in contrast with Rashba SOC Edelstein (1989), Ising SOC does not stabilize a spatially non-uniform superconducting phase.

Understanding the role of impurities is important for future applications of TMDCs and their incorporation into van der Waals heterostructures Geim and Grigorieva (2013). The doping techniques used to prepare the superconducting samples, as well as the defects in the crystal lattice could be the source of significant disorder. Recent weak localization measurements in $\text{MoS}_{2}$ monolayers suggest substantial intervalley scattering, attributed to a high concentration of sulfur vacancies in the monolayer Schmidt et al. (2016). Our result will contribute to identifying the superconducting properties of these materials in the presence of disorder and assess their potential for applications such as in superconducting spintronics.

The model.- The Hamiltonian describing the normal state of TMDC monolayers in the vicinity of the $\pm\mathbf{K}$ points in the presence of a parallel magnetic field is Kormányos et al. (2015) (we use units, where $\hbar=k_{B}=1$):

[TABLE]

Here, $\mathbf{q}=(q_{x},q_{y})$ is a small momentum measured from $\pm\mathbf{K}$, $v$ is the velocity associated with the linearized kinetic dispersion, $E_{g}$ is the difference in on-site energy responsible for the opening of a band gap, and $\Delta_{A}$ and $\Delta_{B}$ are spin-splitting parameters on two different sublattices. Furthermore, $\sigma_{x,y,z}$, $s_{x,y,z}$, and $\eta_{x,y,z}$ are Pauli matrices acting in sublattice, spin, and valley spaces, respectively, and the Zeeman energy $h=\frac{1}{2}g\mu_{B}B$ is related with the amplitude of the magnetic field and the in-plane $g$-factor.

We proceed by projecting the Hamiltonian (1) to the conduction band ($n$-doped regime). The spin-independent part of the Hamiltonian  (1) gives the dominant contribution to the energy of the system and has a simple spectrum $\xi_{\mathbf{q}}=\sqrt{v^{2}|\mathbf{q}|^{2}+E_{g}^{2}}$. It is diagonalized by the unitary transformation $U_{\mathbf{q}}=\exp[i\alpha_{\mathbf{q}}\sigma_{y}\eta_{z}]\exp[i\beta_{\mathbf{q}}\sigma_{z}\eta_{z}]$, with $\tan(2\alpha_{\mathbf{q}})=v|\mathbf{q}|/E_{g}$ and $\tan(2\beta_{\mathbf{q}})=q_{y}/q_{x}$. By assuming $\xi_{\mathbf{q}}\gg\Delta_{A},\Delta_{B}$ and by projecting $U_{\mathbf{q}}\mathcal{H}_{\mathbf{q}}U_{\mathbf{q}}^{\dagger}$ onto the conduction band, we obtain

[TABLE]

where $\eta=\pm 1$ represents the valley index and \Delta_{so}(\mathbf{q})=1/2{\color[rgb]{0,0,0}[\Delta_{A}-\Delta_{B}+(\Delta_{A}+\Delta_{B})E_{g}/\xi_{\mathbf{q}}]} is an effective SOC parameter 111Similar considerations in graphene, where $\Delta_{B}=\Delta_{A}$ and $E_{g}=0$, lead to the suppression of Ising SOC in the projected basis. A Hamiltonian similar to Eq. (2) with $\xi_{\mathbf{q}}\to-\xi_{\mathbf{q}}$ holds in the valence band ($p$-doped regime). Below, we will assume that the chemical potential $\mu$ is sufficiently far above $E_{g}$ 222Ref. Sosenko et al. (2017) analyzes the effect of disorder on a distinct topological superconducting phase when $\mu$ lies in the vicinity of $E_{g}$. on the relevant energy scales determining the superconducting properties ($h,\Delta,T,\dots\ll|\mu-E_{g}|$), so that $\Delta_{so}$ can be taken as a constant. The energy spectrum of Hamiltonian (3) is illustrated in Fig. 1.

Superconductivity is described using a standard BCS Hamiltonian, where the singlet pairing of electrons into Cooper pairs is necessarily intervalley. The corresponding mean-field Hamiltonian reads

[TABLE]

where $c_{\eta\mathbf{q}}=(c_{\eta\mathbf{q}\uparrow},c_{\eta\mathbf{q}\downarrow})^{T}$ is an annihilation operator for spin-up and spin-down electrons, and $\Delta$ is the superconducting order parameter. For compactness, we use the abbreviations $\bar{\eta}=-\eta$ and $\bar{\mathbf{q}}=-\mathbf{q}$. In the vicinity of the second-order superconducting phase transition, $\Delta$ solves the linearized self-consistent gap equation

[TABLE]

where $\lambda$ is the BCS pairing amplitude and $\Omega$ is a cut-off frequency. The particle and hole Green’s functions are given as $\mathcal{G}^{-}_{\eta\mathbf{q}\omega}=(i\omega-\mathcal{H}_{\eta\mathbf{q}})^{-1}$ and $\mathcal{G}^{+}_{\eta\mathbf{q}\omega}=(-i\omega-\mathcal{H}_{\eta\mathbf{q}}^{T})^{-1}$, respectively, with Matsubara frequencies $\omega=(2n+1)\pi T$, $n\in\mathbbm{Z}$, at temperature $T$.

Scattering potential.- The effect of impurities in 2D hexagonal lattices was extensively studied in the case of graphene Aleiner and Efetov (2006); Altland (2006). Here we study the dominant scattering terms, namely spin-independent intra- and intervalley scattering, which may originate from long-range Coulomb interactions or defects in the lattice. Their contribution, expressed in the same basis as the normal state Hamiltonian (1), reads

[TABLE]

The random disorder potentials are characterized by the Gaussian correlators $\langle V_{\alpha}(\mathbf{q})V_{\beta}(\mathbf{\bar{q}^{\prime}})\rangle=w_{\alpha}^{2}\delta_{\alpha,\beta}\delta_{\mathbf{q},\mathbf{q^{\prime}}}$, where $\alpha,\beta=1,2x,2y$, and the brackets denote disorder averaging. The ratio $w_{2i}/w_{1}\sim 1/(|\mathbf{K}|^{2}R^{2})$, which is related with the range $R$ of the potential created by a single impurity, is typically small for remote impurities and of order 1 for lattice defects.

Then, projecting the Hamiltonian $U_{\mathbf{q}}\mathcal{H}_{D}(\mathbf{q-q^{\prime}})U^{\dagger}_{\mathbf{q^{\prime}}}$ onto the conduction band, we find that Eq. (3) has to be supplemented with the disorder term

[TABLE]

where $\mathcal{V}_{1\mathbf{qq^{\prime}}}=V_{1}(\mathbf{q-q^{\prime}})\cos\frac{\theta-\theta^{\prime}}{2}$ and $\mathcal{V}_{2\mathbf{qq^{\prime}}}=[iV_{2x}(\mathbf{q-q^{\prime}})+V_{2y}(\mathbf{q-q^{\prime}})]\sin\frac{\theta+\theta^{\prime}}{2}$ with ${\bf q}^{(^{\prime})}=|{\bf q}^{(^{\prime})}|(\cos\theta^{(^{\prime})},\sin\theta^{(^{\prime})})$.

The disorder-averaged Green’s function can be calculated from the Dyson equation represented diagrammatically in Fig. 2. That is, $\langle\mathcal{G}^{\pm}_{\eta\mathbf{q}\omega}\rangle=({\mathcal{G}^{\pm}_{\eta\mathbf{q}\omega}}^{-1}-\Sigma_{\eta}^{\pm})^{-1}$, where the self-energy $\Sigma_{\eta}^{\pm}$ is obtained using the self-consistent Born approximation. As a result, we find \Sigma_{\eta}^{\pm}=\mp{\color[rgb]{0,0,0}i}[1/(2\tau_{1})+1/(2\tau_{2})]\text{sgn}(\omega), where $1/\tau_{1}=\pi\nu w_{1}^{2}$ and $1/\tau_{2}=\pi\nu(w_{2x}^{2}+w_{2y}^{2})$ are the intra- and intervalley elastic scattering rates, respectively, $\nu=\mu/(2\pi v_{F}^{2})$ is the density of states per spin at the Fermi level, and $v_{F}$ is the Fermi velocity.

Upper critical line.- The upper critical field $h_{c2}(T)=\frac{1}{2}g\mu_{B}H_{c2}(T)$ can be calculated from the disorder-averaged gap equation (4). This involves finding the averages of the products of two Green’s functions, which we do in the standard ladder approximation, as shown in Fig. 3(a).

The first diagram represents the bare vertex, while the second one is a ladder diagram, expressed in terms of eight vertex functions $\Pi_{\eta}^{ss^{\prime}}$(4 combinations of spin indices and 2 values of the valley index), which solve coupled Bethe-Salpeter equations, see Fig. 3(b). Solving these equations yields a remarkably simple expression SM :

[TABLE]

where we used the standard BCS result $T_{c}\simeq\ 1.13\Omega e^{-1/(\lambda\nu)}$. Note that Eq. (7) holds for arbitrary values of the intravalley scattering rate. Below we analyze this equation, which is the main result of this Letter.

Without intervalley scattering.- In the absence of intervalley scattering, $1/\tau_{2}=0$, the critical line given by Eq. (7) does not depend on disorder (Anderson theorem). In that case, Eq. (7) can be alternatively expressed as

[TABLE]

where $\psi(z)$ is the digamma function. In this form, it ressembles – and generalizes to arbitrary disorder – an expression derived by Frigeri et al. Frigeri et al. (2004) in the clean case. It also reproduces the results of Ref. Saito et al. (2016), where the linearized gap equation was solved numerically in the disorder-free case, using a complex multi-band model.

The result can be understood as follows. The effective magnetic field in the two valleys is given by {\color[rgb]{0,0,0}\mathbf{h}}^{\rm eff}_{\eta}=h{\bf e}_{x}+\eta\Delta_{so}{\bf e}_{z}, where ${\bf e}_{i}$ is a unit vector in the $i$-direction. Electrons that are both aligned or anti-aligned with their respective “local” field have the same energy when their momenta are opposite and, thus, their contribution to pairing is not affected by the field. As at finite $h$ the local fields are not along the same axis, however, they enter the gap equation with a suppressed weight $\Delta_{so}^{2}/(h^{2}+\Delta_{so}^{2})$, determined by the overlap of their spin directions. If one electron is aligned whereas the other electron is anti-aligned with their respective “local” field, they have an energy difference of {\color[rgb]{0,0,0}2}\sqrt{h^{2}+\Delta_{so}^{2}} when their momenta are opposite and, thus, their contribution to pairing is suppressed by the field. Here the weight is given as $h^{2}/(h^{2}+\Delta_{so}^{2})$, which is what appears in Eq. (8). Since intravalley scattering does not allow for spin flips, it does not change the result.

As seen in Fig. 4, $h_{c2}$ is enhanced in the presence of spin-orbit coupling, especially at low temperatures. In fact, it diverges in the zero-temperature limit for finite $\Delta_{so}$. Physically, this can be understood as a consequence of the inability of the Zeeman field to completely align the electron spins in the in-plane orientation, due to the anti-parallel out-of-plane field provided by the Ising SOC.

For weak Ising SOC ($\Delta_{so}\ll\Delta_{0}$), the critical curve $h_{c2}(T)$ significantly deviates from the conventional one only at very low temperature, where it diverges (in logarithmic accuracy) as

[TABLE]

where $c$ is a constant of order 1. Close to $T_{c}$, one obtains the standard result $h_{c2}\approx 2.16T_{c}\sqrt{1-T/T_{c}}$.

In the more interesting case of large Ising SOC, $\Delta_{so}\gg\Delta_{0}$, Eq. (8) yields a logarithmic divergence starting at higher temperatures,

[TABLE]

Close to $T_{c}$, the critical field exhibits a standard square-root dependence on temperature, but with an enhanced prefactor,

[TABLE]

With intervalley scattering.- At finite magnetic field, intervalley scattering provides an effective spin-flip mechanism, since electrons scattered between two valleys “feel” opposite values of the Ising SOC field. This pair-breaking effect leads to a saturation of $h_{c2}$ at zero temperature, as illustrated in Fig. 5. For weak intervalley disorder, ${1}/{\tau_{2}}\ll\Delta_{0}\ll\Delta_{so}$, we estimate the zero-temperature critical field as

[TABLE]

in logarithmic accuracy. In the vicinity of $T_{c}$, the critical line is still described by Eq. (11) in that parameter regime.

On the other hand, the standard expression for the second-order paramagnetically-limited critical line, given by Eq. (8) at $\Delta_{so}=0$, is recovered at large disorder strength, $1/\tau_{2}\gg\Delta_{so}^{2}/\Delta_{0}$. In this regime, electrons are frequently scattered between two valleys and do not “feel” the effect of valley-dependent Ising SOC anymore.

Nature of the transition.- In the absence of SOC, the phase transition is a second-order transition into a uniform superconducting state only at sufficiently high temperature. A non-uniform (Fulde-Ferrell-Larkin-Ovchinnikov, or FFLO) phase could possibly contribute to the enhancement of $h_{c2}$, as it was recently discussed in clean bilayer TMDC superconductors Liu (2017). In order to study the nature of the transition in the clean case, we generalize the expression for the gap equation (4) by adding quadratic corrections in a finite modulation wavevector and cubic corrections in the gap amplitude $\Delta$ SM . We find that both do not affect the transition when $\Delta_{so}\gtrsim\Delta_{0}$. Moderate disorder is not expected to change these conclusions FN .

Discussion and conclusion.- Experiments Saito et al. (2016); Lu et al. (2015); Xi et al. (2016) revealed superconductivity in TMDC well above the Pauli limit. The measured fields remained, however, below the values expected in the clean case with Ising SOC only. In Ref. Saito et al. (2016), Rashba SOC was considered as a possible mechanism for the suppression of $h_{c2}$ at low temperature. However, the model required an irrealistically large amplitude for the Rashba SOC. Our work shows that moderate intervalley scattering, which has already been invoked in Ref. Schmidt et al. (2016) in the normal state, could provide an alternative scenario for the saturation. The curve corresponding to $1/(\tau_{2}T_{c})=1.5$ shown in Fig. 5(c) gives a good fit of the experimental data from Ref. Saito et al. (2016) using their estimate for $\Delta_{so}/T_{c}$.

In conclusion, we have studied the effect of disorder on TMDC monolayer superconductors. We have predicted that the large enhancement of the upper critical magnetic field is robust to intravalley scattering. Furthermore, we have identified intervalley scattering as a likely mechanism for the more moderate enhancement of $h_{c2}$ observed in experiment. Interestingly, TMDCs have been identified as a possible platform for topological superconductivity and Majorana fermions Zhou et al. (2016); Sosenko et al. (2017) provided that unconventional pairing takes place. The role of disorder within these scenarios can be investigated within the theory frame provided by our work.

I I. Evaluation of diagrams

In Fig. S1, we show the Feynman rules for the ladder diagrams used in the main text.

The upper and lower branch in all the ladder diagrams represent the Green’s functions $\mathcal{G}^{+}$ and $\mathcal{G}^{-}$, respectively, and all internal momenta (in-between scattering events) are integrated over.

The diagrams for the disorder-averaged self-consistent gap equation, Fig. 3(a) in the main text, translate to

[TABLE]

We evaluate the integrals over the momenta using the residue theorem with the substitution $\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\approx\frac{\nu}{2\pi}\int d\xi_{\mathbf{q}}\int d\theta$, valid in the regime of chemical potentials specified in the main text. In particular, we obtain

[TABLE]

where we used the notation $\tilde{\omega}=\omega+{1}/({2\tau_{1}})+{1}/({2\tau_{2}}$).

We note that, in general, the renormalized vertex function $\Pi_{\eta}$ is dependent on the polar angle $\theta$, due to the anisotropy of the projected disorder potential, Eq. (6). The $\Pi_{\eta}$ are determined by a system of Bethe-Salpeter equations [corresponding to diagrams in Fig. 3(b) in the main text]:

[TABLE]

We readily check that $\Pi_{\eta}(\theta)$ is fully determined by its first harmonics in $\theta$:

[TABLE]

Namely, by combining Eqs. (S4) and (S3), after integration over $\theta^{\prime}$, we verify that no higher harmonics are generated. Furthermore, the equation for the constant part of the vertex functions $\Pi_{\eta 0}$ is decoupled from the angle-dependent parts $\Pi_{\eta 1}$ and $\Pi_{\eta 2}$. Replacing the $\Pi_{\eta}(\theta)$ in Eq. (S1), we see that the angle-dependent contributions vanish after the integration over momenta. Therefore, it is sufficient to compute only $\Pi_{\eta 0}$.

The eight different $\Pi_{\eta 0}$ are determined from the linear system of equations obtained after integrating (S3) over angles:

[TABLE]

Inserting the vertex functions that solve Eq. (S5) into Eq. (S1) yields Eq. (7) in the main text.

II II. Limiting cases

In the following, starting from the full expression for the critical curve $h_{c2}(T)$, given by Eqs. (7) and (8) in the main text, we derive simple analytical expressions in the zero-temperature limit as well as close to the transition temperature $T_{c}$, for various strengths of the intervalley scattering rate. We compare them with the results of the numerical calculations of $h_{c2}$ for arbitrary disorder strength in the same temperature regimes.

II.1 A. Limit $T\rightarrow 0$

In the absence of intervalley scattering, $1/\tau_{2}=0$, the critical curve $h_{c2}(T)$ is described by Eq. (8) in the main text. At low temperatures, using the asymptotic behavior of the digamma function $\psi(z)\approx\ln|z|$ for $|z|\gg 1$, Eq. (8) in the main text thus yields

[TABLE]

where $\gamma\approx 0.577$ and we introduced $\rho=\sqrt{h_{c2}^{2}+\Delta_{so}^{2}}$. In the limit $T\to 0$, $h_{c2}$ diverges, as illustrated in Fig. 4 in the main text. Thus, we can approximate our results assuming $h_{c2}\gg\Delta_{so},\Delta_{0}$. Then, Eq. (S6) yields

[TABLE]

In the limit of strong Ising SOC, $\Delta_{so}\gg\Delta_{0}$, the condition $h_{c2}\gg\Delta_{so}$ allows one to simplify Eq. (S7) in the temperature regime $T/T_{c}<\Delta_{0}/\Delta_{so}$ to Eq. (10) in the main text, in logarithmic accuracy. In the limit of weak Ising SOC, $\Delta_{so}\ll\Delta_{0}$, we see that the condition $h_{c2}\gg\Delta_{0}$ is fulfilled for $(\Delta_{so}/\Delta_{0})^{2}\ln(T_{c}/T)\gg 1$. Then, we are able to give a rough estimate provided by Eq. (9) in the main text.

At weak disorder, ${1}/{\tau_{2}}\ll\Delta_{so},\Delta_{0}$, the poles in the $\omega$-dependent terms in Eq. (7) in the main text can be evaluated perturbatively to yield:

[TABLE]

By comparing with Eq. (8) in the main text, we see that the main effect of weak intervalley scattering is to provide an effective pair-breaking rate for pairs of electrons that are both aligned or anti-aligned with their respective local fields, which yields an additional mechanism for the suppression of $h_{c2}$. As a consequence, $h_{c2}$ now saturates at $T\to 0$. Assuming $1/\tau_{2}\ll\Delta_{0}\ll\Delta_{so}$ such that $h_{c2}\gg\Delta_{so}$ still holds, we find from Eq. (S8) that

[TABLE]

which evaluates to Eq. (12) in the main text, in logarithmic accuracy. For weak Ising SOC, $\frac{1}{\tau_{2}}\ll\Delta_{so}\ll\Delta_{0}$, similar estimates for $h_{c2}$ can be made up to a constant factor, in the same manner as in Eq. (9).

For large disorder strength, Eq. (7) in the main text evaluates to

[TABLE]

provided that only Matsubara frequencies $\omega\lesssim h_{c2}\ll 1/\tau_{2}$ contribute to the sum. At $T\to 0$, Eq. (S10) simplifies to

[TABLE]

Thus, at large Ising SOC, $\Delta_{so}\gg\Delta_{0}$, we find the standard paramagnetic limit $h_{c2}=\Delta_{0}/2$ for the second-order phase transition if $1/\tau_{2}\gg\Delta_{so}^{2}/\Delta_{0}$. At weaker disorder, $\Delta_{so}\ll 1/\tau_{2}\ll\Delta_{so}^{2}/\Delta_{0}$, we obtain $h_{c2}=\pi\tau_{2}\Delta^{2}_{so}/[2\ln(\tau_{2}\Delta_{so}^{2}/\Delta_{0})]$. Further decreasing disorder, the condition $h_{c2}\ll 1/\tau_{2}$ breaks down and Eq. (S10) cannot be used anymore.

To compare these results with the numerical evaluation of $h_{c2}(0)$, it is convenient to obtain the later from the following equation:

[TABLE]

We plot $h_{c2}(0)$ obtained this way in Fig. S2(a) as a function of $1/\tau_{2}$. We verify that Eq. (12) is in good agreement with these results in the relevant parameter regime and that $h_{c2}(0)$ reaches $\Delta_{0}/2$ for strong disorder.

II.2 B. Limit $T\rightarrow T_{c}$

Close to $T_{c}$, the critical field $h_{c2}$ is small, and the assumptions $\Delta_{so}\gg h_{c2}$ and $\ln\frac{T_{c}}{T}\approx\frac{T_{c}-T}{T_{c}}$ hold. Expanding Eq. (7) in the main text in this regime yields a square-root behavior

[TABLE]

with

[TABLE]

and $\omega=(2n+1)\pi T_{c}$ ($n$ integer). Analyzing Eq. (S14) in various regimes, we find $C=7\zeta(3)/4\pi^{2}$ at weak disorder and weak Ising SOC, $1/\tau_{2},\Delta_{so}\ll\Delta_{0}$, as well as at large disorder, $1/\tau_{2}\gg\Delta_{so}^{2}/\Delta_{0}\gg\Delta_{so}$. Outside those ranges of disorder, $C$ is suppressed – and, thus, $h_{c2}$ is enhanced. For instance, $C=(T_{c}/\Delta_{so})^{2}\ln(2\Delta_{so}/\Delta_{0})$ at weak disorder and large Ising SOC, $1/\tau_{2}\ll\Delta_{0}\ll\Delta_{so}$, and $C=\pi T_{c}/(4\tau_{2}\Delta_{so}^{2})$ in an intermediate disorder range and at large Ising SOC, $\Delta_{0}\ll 1/\tau_{2},\Delta_{so}\ll\Delta_{so}^{2}/\Delta_{0}$.

We explore a wider window of disorder strengths numerically using Eq. (S14), as shown in Fig. S2. The results match the approximate formulas close to $T_{c}$ given above at weak and strong disorder.

III III. FFLO phase and the First-order phase transition

In the main text, we only considered the second-order transition into a homogeneous superconducting state. In the following, we examine the conditions for realizing a FFLO phase and/or first-order phase transition in clean Ising superconductors. We show that both are absent for $\Delta_{so}\gtrsim 0.52\Delta_{0}$. We do not expect moderate disorder, $1/\tau_{2}\ll{\rm max}(\Delta_{0},\Delta_{so}^{2}/\Delta_{0})$, to lead to a reappearance of the FFLO phase and/or first-order phase transition.

III.1 A. FFLO phase

In the FFLO phase, the superconductor is spatially modulated. At the second-order transition into that state, we account for an exponentially modulated order parameter $\Delta(\mathbf{x})=\Delta_{\mathbf{p}}e^{i\mathbf{p}.\mathbf{x}}$, where $\mathbf{p}$ is the modulation wavevector, by modifying the pairing term in the BCS Hamiltonian,

[TABLE]

Then, in the vicinity of the phase transition, the amplitude of the order parameter should solve the linearized self-consistent gap equation

[TABLE]

In order to consider the instability toward an FFLO state along the upper ciritcal line $h_{c2}(T)$, we further assume that the modulation wavevector is small, $v_{F}|\mathbf{p}|\ll\Delta_{0}$. The small momentum shift $\mathbf{p}/2$ in the Green’s functions can be accounted for by shifting the energy $\xi_{\mathbf{q}+\frac{\mathbf{p}}{2}}\approx\xi_{\mathbf{q}}+v_{F}\mathbf{p}.\hat{\mathbf{q}}/2$, where we introduced $\hat{\mathbf{q}}=\mathbf{q}/|\mathbf{q}|$. Then, integration over the momenta can be carried out using the results given in Eq. (I) by taking $\tilde{\omega}\to\omega-iv_{F}\mathbf{p}.\hat{\mathbf{q}}/2$. After summing over Matsubara frequencies, we obtain

[TABLE]

Expanding the above expression in small $\mathbf{p}$ and integrating over angles yields

[TABLE]

where

[TABLE]

Here, $\psi^{(2)}(z)$ is the second derivative of the digamma function.

The last term in Eq. (S18) is the correction to the result in the uniform case, Eq. (8), due to the modulation. The instability toward the FFLO state is determined by the sign of $F_{1}$ along the line $h_{c2}(T)$ for the uniform state. Namely, if $F_{1}>0$ (resp. $F_{1}<0$), $h_{c2}$ decreases (resp. increases) when the order parameter is modulated.

We evaluate $F_{1}$ along the upper critical line derived for the uniform state in Fig. S3(a). At $\Delta_{so}=0$, $F_{1}$ changes sign at $T^{*}=0.56T_{c}$, signaling a transition into the FFLO state below that temperature. At small $\Delta_{so}$, we find that $F_{1}$ changes sign at two temperatures $T^{*}_{1}$ and $T^{*}_{2}$, with $T^{*}_{1}<T<T^{*}_{2}$. The range of temperatures $T^{*}_{1}<T<T^{*}_{2}$, where the FFLO state can be expected shrinks as $\Delta_{so}$ increases, and it eventually disappears at $\Delta_{so}\gtrsim 0.30\Delta_{0}$, thus excluding the possibility of an FFLO phase at larger $\Delta_{so}$.

III.2 B. First-order phase transition

In order to study the possibility of a first-order phase transition, the linearized self-consistency equation is not sufficient and we need to include higher order terms in $\Delta$. Thus, we write the self-consistent gap equation for the uniform phase [described by Eq. (3)] up to third order in $\Delta$,

[TABLE]

After evaluating the integrals over the products of four Green’s functions using the residue theorem [in a similar fashion as in Eq. (I)], and summing over Matsubara frequencies, we obtain

[TABLE]

where

[TABLE]

The last term in Eq. (S21) is the correction to the linearized gap equation (8) due to a finite amplitude of the order parameter in the vicinity of the transition. The order of the transition is determined by the sign of $F_{2}$ along the second-order transition line $h_{c2}(T)$. Namely, if $F_{2}>0$ (resp. $F_{2}<0$), the transition remains second-order (resp. a change of the order of the transition occurs).

At $\Delta_{so}=0$, we find that $F_{2}=(h_{c2}/T_{c})^{2}F_{1}$. Thus, the sign change occurs at the same temperature and, as a consequence, for $T<T^{*}=0.56T_{c}$ the transition into the FFLO state is in competition with a first-order transition. At finite $\Delta_{so}$, we evaluate $F_{2}$ along the upper critical line derived for the uniform state in Fig. S3(b). We find that its temperature dependence is qualitatively similar to, though quantitatively different from $F_{1}$. Thus, a change of the transition order may occur in a finite temperature range, if Ising SOC is weak. On the other hand, $F_{2}$ remains positive at all temperatures if $\Delta_{so}\gtrsim 0.52\Delta_{0}$, and therefore the transition remains a second-order transition at larger Ising SOC.