Superconducting states and Majorana modes in transition-metal dichalcogenides under inhomogeneous strain

TL;DR

This paper investigates how inhomogeneous strain influences superconductivity and Majorana modes in transition-metal dichalcogenides, revealing topological properties and conditions for Majorana mode emergence.

Contribution

It demonstrates the preservation of superconductivity within strain-induced pseudo Landau levels and predicts Majorana modes in strained dichalcogenides under specific conditions.

Findings

Superconductivity persists in strain-induced pseudo Landau levels.

Quantum critical behavior at integer fillings depends on LL spacing.

Majorana modes can emerge at edges under Zeeman field and valley symmetry breaking.

Abstract

We study the effect of inhomogeneous strain on transition-metal dichalcogenides with a large intrinsic gap in their spectrum. It is found that, by tuning the chemical potential, superconductivity can preserve within the strain-induced discrete pseudo Landau levels (LLs), which introduce interesting topological properties to these systems. As we show, the superconductivity for integer fillings is quantum critical, and the quantum critical coupling strength is determined by the spacing between the two LLs closest to the Fermi level. For partial fillings, the superconducting gap is scaled linearly with the coupling strength, and decreases rapidly when the chemical potential shifts away from the middle of each LL. In the presence of a Zeeman field, a pair of Majorana modes emerge simultaneously in the two valleys of strained dichalcogenides. When valley symmetry is further destroyed, a single Majorana mode can be expected to emerge at the edges of the strained monolayer dichalcogenides.

Equations

$$H=\int d\mathbf{x}\sum _{\xi }{\psi }_{\xi }^{†}(\mathbf{x})H_{\mathbf{p},\xi }{\psi }_{\xi }(\mathbf{x}),$$

$$H_{\mathbf{p},\xi }={\upsilon }_{\mathrm{F}}{\vec{\Pi }}^{\xi }\cdot {\vec{\sigma }}_{\xi }+\frac{{\Delta }_0}{2}{\sigma }_z-({\lambda }_{\mathrm{so}}\xi \frac{{\sigma }_z-1}{2}-m_{\mathrm{z}})s_z-\mu$$

$$\mathcal{H}_{\mathrm{BdG}}=\left(\begin{array}[c]{cc}H_{\mathbf{p},\xi}+m_{\mathrm{z}}s_{z}&\hat{\Delta}_{4}\\ \hat{\Delta}_{4}^{{\dagger}}&-\mathcal{T}H_{-\mathbf{p},-\xi}\mathcal{T}^{-1}-m_{\mathrm{z}}s_{z}\end{array}\right)$$

$$\widetilde{\mathcal{H}}_{\mathrm{BdG}}=\left(\begin{array}[c]{cccc}h_{\mathbf{p},\xi\uparrow}&\hat{\Delta}_{2}&0&0\\ \hat{\Delta}_{2}^{{\dagger}}&-h_{-\mathbf{p},-\xi\downarrow}^{\ast}&0&0\\ 0&0&h_{\mathbf{p},\xi\downarrow}&\hat{\Delta}_{2}\\ 0&0&\hat{\Delta}_{2}^{{\dagger}}&-h_{-\mathbf{p},-\xi\uparrow}^{\ast}\end{array}\right)$$

$$\psi_{k,\xi\sigma}^{(n)}=\frac{e^{ikx}}{\sqrt{2}}\left(\begin{array}[c]{c}s_{n}\alpha_{\xi\sigma,+}^{(n)}\phi_{|n|-1}(\zeta)\\ \alpha_{\xi\sigma,-}^{(n)}\phi_{|n|}(\zeta)\end{array}\right),$$

$${\alpha }_{\xi \sigma ,\pm }^{(n)}=\sqrt{1\pm {\lambda }_{\xi \sigma }/{\Omega }_{\xi \sigma }^n}$$

$$\widetilde{\mathcal{H}}_{\mathrm{BdG}}^{(m,n)}=\langle\Phi^{(m)}|\widetilde{\mathcal{H}}_{\mathrm{BdG}}|\Phi^{(n)}\rangle=\left(\begin{array}[c]{cc}\mathcal{H}_{+}&0\\ 0&\mathcal{H}_{-}\end{array}\right)\delta_{m,n}\ ,$$

$${\mathcal{N}}_{\chi }=-g{N}_{\varphi }\sum _n{\int }_{-\infty }^{\infty }d\omega \mathrm{Im}[G_{\chi ,11}^r(\omega )/\pi ]f(\omega ),$$

$$G_{\chi }^r(\omega )$$

$$\displaystyle=\frac{1}{2}\sum_{\eta=\pm}\frac{1}{\omega^{+}-\eta E_{n}^{\chi}}\left[\tau_{0}+\frac{1}{\eta E_{n}^{\chi}}\left(\begin{array}[c]{cc}\varepsilon_{n,\chi}&\Delta\\ \Delta&-\varepsilon_{n,\chi}\end{array}\right)\right]$$

$${\mathcal{N}}_{\chi }=\frac{g{N}_{\varphi }}{2}\sum _{n=-n_D}^{n_{\mathrm{F}}}\left(1-\frac{{\epsilon }_{n,\chi }\tanh [E_n^{\chi }/(2k_{B}T)]}{E_n^{\chi }}\right),$$

$$2({\nu }_{\chi }/g-n_{\mathrm{F}})=-\frac{{\epsilon }_{n_{\mathrm{F}},\chi }\tanh [E_{n_{\mathrm{F}}}^{\chi }/(2k_{B}T)]}{E_{n_{\mathrm{F}}}^{\chi }}$$

$${\mu }_{\chi }(T,{\nu }_{\chi })={\Omega }_{\chi }^{n_{\mathrm{F}}}+\frac{2\Delta (T,{\nu }_{\chi })({\nu }_{\chi }/g-n_{\mathrm{F}})}{\sqrt{1-4({\nu }_{\chi }/g-n_{\mathrm{F}}{)}^2}}.$$

$$\Delta (T,{\nu }_{\chi })=-Ug{\overline{N}}_{\varphi }\sum _n{\int }_{-\infty }^{\infty }d\omega \mathrm{Im}[G_{\chi ,12}^r(\omega )/\pi ]f(\omega )$$

$$1=-(U/2)g{\overline{N}}_{\varphi }\sum _n\tanh [E_n^{\chi }/(2k_{B}T)]/E_n^{\chi }.$$

$$\Delta (0,{\nu }_{\chi })=\frac{\hslash {\upsilon }_{\mathrm{F}}\sqrt{1-4({\nu }_{\chi }/g-n_{\mathrm{F}}{)}^2}}{2l_B[1-{\gamma }_{n_{\mathrm{F}}}^{(1)}x]}x,$$

$${\gamma }_{n_{\mathrm{F}}}^{(k)}=\sum _{n<n_{\mathrm{F}}}\frac{1}{2|({\epsilon }_{\chi ,n_{\mathrm{F}}}-{\epsilon }_{\chi ,n})/\hslash {\omega }_c{|}^k}$$

$${\mu }_{\chi }(0,{\nu }_{\chi })={\Omega }_{\chi }^{n_{\mathrm{F}}}+\frac{\hslash {\upsilon }_{\mathrm{F}}({\nu }_{\chi }/g-n_{\mathrm{F}})}{l_B[1-{\gamma }_{n_{\mathrm{F}}}^{(1)}x]}x.$$

$$\Delta (0,{\nu }_{\chi }^{\mathrm{I}})=\hslash ({\upsilon }_{\mathrm{F}}/l_B)|{\Gamma }_{\chi }(n_{\mathrm{F}})|\sqrt{(x/x_c{)}^2-1}$$

$$x_c=\frac{1}{{\gamma }_{n_{\mathrm{F}}-1}^{(1)}+|{\Gamma }_{\chi }(n_{\mathrm{F}}){|}^{-1}}$$

$${\Gamma }_{\chi }(n_{\mathrm{F}})=\frac{{\epsilon }_{n_{\mathrm{F}},\chi }-{\epsilon }_{n_{\mathrm{F}}-1,\chi }}{2\hslash {\omega }_c}.$$

$$|{\epsilon }_{n,\chi }-{\epsilon }_{n-1,\chi }|=\frac{2(\hslash {\omega }_c{)}^2}{|{\Omega }_{\chi }^n+{\Omega }_{\chi }^{n-1}|}$$

$$f(\xi )=k_{B}T\sum _{m=-\infty }^{\infty }\frac{e^{i{\omega }_m\delta }}{i{\omega }_m-\xi }$$

$$1=x{k}_{B}T\sum _n\sum _{m=-\infty }^{\infty }\frac{\hslash {\omega }_c}{{\omega }_m^2+{\epsilon }_{n,\chi }^2+{\Delta }^2(T,{\nu }_{\chi })}.$$

$$1$$

$$-x\hslash {\omega }_c{k}_{B}T\sum _n\sum _{m=-\infty }^{\infty }\frac{{\Delta }^2(T,{\nu }_{\chi })}{({\omega }_m^2+{\epsilon }_{n,\chi }^2{)}^2}.$$

$$\Delta (T,{\nu }_{\chi }^{\mathrm{H}})=(\hslash {\upsilon }_{\mathrm{F}}/l_B{)}^{3/2}\sqrt{1-T/T_c}/[2k_B{T}_c{\gamma }_{n_{\mathrm{F}}}^{(3)}{]}^{1/2}$$

$$C_1^{\xi \sigma }$$

$$\times \frac{⟨{\Phi }_{m,\lambda }^{\xi \sigma }|{\upsilon }_{x,\xi \sigma }|{\Phi }_{n,{\lambda }^{\prime }}^{\xi \sigma }⟩⟨{\Phi }_{n,{\lambda }^{\prime }}^{\xi \sigma }|{\upsilon }_{y,\xi \sigma }|{\Phi }_{m,\lambda }^{\xi \sigma }⟩}{(E_{m,\lambda }^{\xi \sigma }-E_{n,{\lambda }^{\prime }}^{\xi \sigma }{)}^2},$$

$$E_{n,\lambda }^{\xi \sigma }=\lambda \sqrt{({\Omega }_{\xi \sigma }^n-{\mu }_{\xi \sigma }{)}^2+{\Delta }^2}+\sigma m_{\mathrm{z}},$$

$$\Phi_{n,\lambda}^{\xi\sigma}=\frac{1}{\sqrt{2}}\left(\begin{array}[c]{c}\beta_{\xi\sigma,+}^{(n,\lambda)}\\ \lambda\beta_{\xi\sigma,-}^{(n,\lambda)}\end{array}\right)\otimes\psi_{k,\xi\sigma}^{(n)}\ ,$$

$${\beta }_{\xi \sigma ,\pm }^{(n,\lambda )}=\sqrt{1\pm ({\Omega }_{\xi \sigma }^n-{\mu }_{\xi \sigma })/(E_{n,\lambda }^{\xi \sigma }-\sigma m_{\mathrm{z}})}.$$

Full text

I introduction

Since the remarkable discovery of grapheneNovoselov et al. (2004); Novoselov et al. (2005a); Zhang et al. (2005), the study of physics in atomically thin two dimensional (2D) crystals, which could be of great potential applications in next-generation nanoelectronic devicesNovoselov et al. (2005b); Lee et al. (2010), has attracted much attention on both theoretical and experimental sidesGoerbig (2011); Kotov et al. (2012). In graphene, the conduction and valence band touch at the corners, referred to as $K$ and $K^{\prime}$ points, of the 2D hexagonal Brillouin zone. The two inequivalent points constitute a binary index, termed as the valley index, for the low energy carriers. In the vicinity of the $K$ ($K^{\prime}$) points, the low-energy electronic excitations behave as massless Dirac quasiparticlesGomes et al. (2012), and the dispersions form a 2D Dirac cone, whose vertex is called Dirac point. The two valleys are separated far from each other in the momentum space. For electronic states closed to the Dirac points, the valley index is expected to be robust against scattering by perturbations. Therefore, the valley index can be served as a potential information carrier, the use of which leads to a new concept, valleytronicsGunawan et al. (2006); Yao et al. (2008); Xiao et al. (2007); Rycerz et al. (2007); Zhang et al. (2011). When the inversion symmetry is broken, valley Hall effect can emergeXiao et al. (2007), where carriers in different valleys flow to opposite transverse directions upon application of an electric field. Moreover, when including the spin-orbit interaction, one can explore spin physics and spintronics in grapheneZhang et al. (2011); Min et al. (2006); Yao et al. (2007); Avsar et al. (2014); Cummings et al. (2017).

Interestingly, recent theoretical and experimental studies showed that intrinsic superconductivity could be induced in graphene under the application of strain fields Levy et al. (2010); Gomes et al. (2012); Abanin and Pesin (2012); Ghaemi et al. (2012); Uchoa and Barlas (2013). The strain introduces pseudo Landau levels (LLs) into graphene, while the time-reversal (TR) symmetry remains. In the weak coupling regime, the critical temperature is found to scale linearly with the coupling strength, which is quite different from the conventional weak-coupling superconductors where the critical temperature decreases exponentially with the effective coupling Uchoa and Barlas (2013). By modulating the filling factor and magnitude of strain, one can control the superconducting transition temperature experimentally, which has profound significance for the manipulation of quantum states in solid states. Moreover, in the presence of superconductivity, the system can exhibit exotic topological properties, such as the emergence of Majorana modes, when the inversion and TR symmetries are broken spontaneously Wang and Wu (2016).

Although graphene has many extraordinary physical properties, the inversion symmetry is preserved and the spin-orbit coupling (SOC) is rather weak in graphene, which challenges some of its applications in valleytronics and spintronics. Instead, layered transition-metal dichalcogenides Novoselov et al. (2005b); Lee et al. (2010); Splendiani et al. (2010); Mak et al. (2010); Radisavljevic et al. (2011); Korn et al. (2011); Wakamura et al. (2018); Zihlmann et al. (2018), with broken inversion symmetry and strong SOC, represent an alternative class of 2D materials Novoselov et al. (2005b), which can provide excellent platforms towards the integration of valleytronics and spintronics Wang et al. (2012); Xiao et al. (2012). For example, monolayer MoS2 has similar hexagonal lattice structure as graphene Xiao et al. (2012), but the inversion symmetry is broken explicitly and the SOC is much stronger in MoS2, which makes it promising candidate materials for valleytronics and spintronics Garcia et al. (2017); Schmidt et al. (2016); Zhang et al. (2017). In fact, the physics in monolayers of group-VI dichalcogenides $MX_{2}$, with $M=$Mo and $X=$S, Se, is essentially the same, all of which are identified as direct-band-gap semiconductors Zhu et al. (2011). It is of importance to understand theoretically how the appearance of the intrinsic band gap and SOC influence the superconductivity and topological properties of these emergent 2D materials under strain fields.

In this paper, we investigate the superconductivity and topological properties of the electrons in strained transition-metal dichalcogenides. Taking MoS2 as an example, we generalize the superconductivity theory from gapless graphene Uchoa and Barlas (2013) to strained dichalcogenides with intrinsic band gaps in their spectrum. We find that the superconductivity can preserve within the discrete pseudo LLs for these gapped systems, and the resulting topological phenomena are very interesting. In the presence of a finite energy gap, the chemical potential plays a very important role in the occurrence of superconductivity, even for the $n=0$ LL with $n$ as the LL index. At half fillings, the chemical potential sitting at the middle of the LLs, the superconductivity gap is maximized. However, at integer fillings, the emergence of superconductivity requires a minimal quantum critical coupling, whose strength is determined by the spacing between the LLs closest to the Fermi level. Below the quantum critical coupling strength, the superconductivity is fully suppressed. Interestingly, in the presence of a Zeeman field, a pair of Majorana modes emerge simultaneously in the two valleys of strained dichalcogenides. A single Majorana mode can be expected to emerge at the edges of the sample, if the valley symmetry is further destroyed.

The rest of this paper is organized as follows. In the next section, we introduce the model Hamiltonian and method. The superconductivity in strained MoS2 is discussed in Sec. III, and its topological properties are analyzed in Sec. IV. The final section contains a summary.

II Model Hamiltonian and Method

As demonstrated in Refs. Xiao et al. (2012); Zhu et al. (2011), the underlying physics is the same for monolayers of group-VI dichalcogenides, such that we can take one of them, i.e., MoS2, as an example. The low-energy electronic excitation in strained monolayer MoS2 can be described by the Hamiltonian Ghaemi et al. (2012); Uchoa and Barlas (2013); Xiao et al. (2012)

[TABLE]

where $\psi_{\xi}(\mathbf{x})=(c_{A,\xi\uparrow},c_{B,\xi\uparrow},c_{A,\xi\downarrow},c_{B,\xi\downarrow})^{T}$, $c_{A(B),\xi\sigma}$ are electron annihilation operators, and

[TABLE]

with $\sigma_{i}$ and $s_{i}$ being the Pauli matrices for sublattice and spin, respectively, and $\xi=\pm$ representing the valley index. Here, $\upsilon_{\mathrm{F}}=at/\hbar$ is the Fermi velocity with $a$ and $t$ the lattice constant and electron hopping integral. $\Delta_{0}$ is the intrinsic energy gap due to the broken inversion symmetry, $2\lambda_{\mathrm{so}}$ is the spin splitting at the valence band top due to the spin-orbit coupling, $m_{\mathrm{z}}$ is the Zeeman field, and $\mu$ is the chemical potential. The valley-dependent gauge covariant momentum operator $\vec{\Pi}_{x,y}^{\xi}=\hat{p}_{x,y}+\xi eA_{x,y}$ is modified by the pseudo-vector-potential $\mathbf{A}=(\delta t_{x},\delta t_{y})/e\upsilon_{\mathrm{F}}$ generated by the strain and $\vec{\sigma}_{\xi}=(\xi\sigma_{x},\sigma_{y})$ is a vector of the Pauli matrices. In the presence of an effective attractive potential $U$, which stabilizes the superconducting state, the Bogoliubov-de Gennes (BdG) Hamiltonian is given by $H_{\mathrm{BdG}}=\frac{1}{2}\int d\mathbf{x}\sum_{\xi}\Psi_{\xi}^{{\dagger}}(\mathbf{x})\mathcal{H}_{\mathrm{BdG}}\Psi_{\xi}(\mathbf{x})$, where

[TABLE]

and $\Psi_{\xi}^{{\dagger}}(\mathbf{x})=(\psi_{\xi}^{{\dagger}},is_{y}\psi_{-\xi})$, with $\hat{\Delta}_{4}=\Delta s_{z}\otimes\sigma_{0}$, $\mathcal{T}=is_{y}\mathcal{K}$, and $\mathcal{K}$ denoting complex conjugation. In a proper basis order, the BdG Hamiltonian can be rewritten in the block diagonal form as

[TABLE]

with $h_{\mathbf{p},\xi\sigma}=\upsilon_{\mathrm{F}}\vec{\Pi}^{\xi}\cdot\vec{\sigma}_{\xi}+\lambda_{\xi\sigma}\sigma_{z}-\mu_{\xi\sigma}+\sigma m_{\mathrm{z}}$, where $\lambda_{\xi\sigma}=\frac{\Delta_{0}+\sigma\xi\lambda_{\mathrm{so}}}{2}$ and $\mu_{\xi\sigma}=\mu+\frac{\sigma\xi\lambda_{\mathrm{so}}}{2}$.

In the absence of the Zeeman field, the spin index $\sigma=\uparrow,\downarrow$ is locked to the valley index $\xi$, such that $h_{\mathbf{p},\xi\sigma}=h_{\overline{\mathbf{p}},\overline{\xi}\overline{\sigma}}^{\ast}$, where $\overline{\sigma}\equiv-\sigma$. As it shows, the strain-induced pseudomagnetic field does not break the TR symmetry for the system, i.e., $H_{-\mathbf{p},-\xi}=\mathcal{T}H_{\mathbf{p},\xi}\mathcal{T}^{-1}$, which is essentially different from a conventional magnetic field. Therefore, TR-symmetric states can pair up by the effective attractive interaction, which favors the formation of Cooper pairs. The pairing matrix $\hat{\Delta}_{2}=\Delta\sigma_{0}$ in Eq. (4), describing the formation of Cooper pairs, can be determined self-consistently by $\Delta^{2}=$ $U\mathrm{tr}\langle\psi_{k,\xi\sigma}|\hat{\Delta}_{2}|\psi_{-k,\overline{\xi}\overline{\sigma}}\rangle$, where $\psi_{k,\xi\sigma}$ is the two-component spinor for $h_{\mathbf{p},\xi\sigma}$. In the Landau gauge $\mathbf{A}=(-By,0)$, with $B$ as the pseudomagnetic field, the spinor for $h_{\mathbf{p},\xi\sigma}$ takes the following form

[TABLE]

when $\mu_{\xi\sigma}=m_{\mathrm{z}}=0$. Here, $s_{n}\equiv\mathrm{sgn}(n)$, $\phi_{|n|}(\zeta)$ is the harmonic wavefunction, and

[TABLE]

with $\Omega_{\xi\sigma}^{n}=s_{n}\sqrt{2|n|(\hbar\omega_{c})^{2}+\lambda_{\xi\sigma}^{2}}-\lambda_{\xi\sigma}\delta_{n,0}$ and $\zeta=\xi kl_{B}-y/l_{B}$. The cyclotron frequency is defined as $\omega_{c}=\upsilon_{\mathrm{F}}/l_{B}$, in which $l_{B}=\sqrt{\hbar/eB}$ denotes the magnetic length.

III Superconducting States

Since the Zeeman field is irrelevant to the emergence of superconductivity, we would omit it in the discussions of the properties of superconducting states, for simplicity, and it will be considered later in the discussions of the topological properties in Sec. IV. By projecting $\widetilde{\mathcal{H}}_{\mathrm{BdG}}$ into the Hilbert space $\{\Phi^{(n)}\}$, with $\Phi^{(n)}=$ $\left(\psi_{k,\xi\uparrow}^{(n)},\psi_{-k,\overline{\xi}\downarrow}^{(n)\ast},\psi_{k,\xi\downarrow}^{(n)},\psi_{-k,\overline{\xi}\uparrow}^{(n)\ast}\right)$, the BdG Hamiltonian can be written in the diagonal form as

[TABLE]

where $\mathcal{H}_{\chi}=(\Omega_{\chi}^{n}-\mu_{\chi})\tau_{z}+\Delta\tau_{x}$, with $\chi=\pm$ labeling the product of the valley and spin indices ($\xi\times\sigma$) and $\tau_{i}$ being the Pauli matrix for Bogoliubov isospin in the Nambu space. Some details of derivation are given in the appendixA, e.g., Eq. (51). In the presence of finite energy gap $\Delta_{0}$ and spin-orbit coupling $\lambda_{\mathrm{so}}$, the $n=0$ LL $\Omega_{\xi\sigma}^{0}-\mu_{\xi\sigma}$ in strained MoS2 shifts away from zero energy, which is different from the case in gapless graphene Uchoa and Barlas (2013). As a result, for $\mu=0$, to preserve the superconductivity in the discrete spectrum of LLs, a finite coupling is needed to overcome the energy gap between the electron and hole LLs, even for the $n=0$ LL. With tuning the chemical potential, the electron and hole LLs could approach each other and encounter at the Fermi level, such that the gap between the electron and hole LLs would decrease, which reduces the critical coupling strength. Therefore, in the presence of an intrinsic gap, the chemical potential plays a very important role in the occurrence of superconductivity.

Due to Pauli blocking, the electronic states are incompressible, so that the chemical potential would exhibit a discontinuous behavior with changing the pseudomagnetic field and filling factor. As shown in Eq. (7), $\mathcal{H}_{\pm}$ are block diagonal and there is no particle exchange between the two subspaces. Consequently, we can concentrate on one subspace, say $\mathcal{H}_{\chi}$, and the conclusions can be generalized to the other subspace straightforwardly. Following Ref. Uchoa and Barlas (2013), we calculate the chemical potential by fixing the number of particles $\mathcal{N}_{\chi}$, which can be determined by the fluctuation-dissipation theorem Deng et al. (2016)

[TABLE]

where $g=2$ is the valley$\times$spin degeneracy and $N_{\phi}=A/(2\pi l_{B}^{2})$ is degeneracy of the LL originating from the summation over $k$, with $A=L_{x}L_{y}$ as area of the sample. $f(\omega)=1/(1+e^{\omega/k_{B}T})$ is the Fermi-Dirac distribution function and the retarded Green’s function

[TABLE]

is a $2\times 2$ matrix in the Nambu space, with $E_{n}^{\chi}=\sqrt{\varepsilon_{n,\chi}^{2}+\Delta^{2}}$ and $\varepsilon_{n,\chi}=\Omega_{\chi}^{n}-\mu_{\chi}$, where we adopted the abbreviation $\omega^{+}=\omega+i0^{+}$. The subscripts of the retarded Green’s function, e.g., $G_{\chi,11}^{r}(\omega)$ and $G_{\chi,12}^{r}(\omega)$, stand for the matrix elements of $G_{\chi}^{r}(\omega)$. If the deep-energy states $n<n_{\mathrm{F}}$ are fully occupied, with $n_{\mathrm{F}}$ labeling the highest occupied LL, we can derive

[TABLE]

where $n_{D}$ is an ultravialet cutoff related to the bandwidth $D$. More details are presented in the appendixA.

In the low temperature and weak coupling limit $k_{B}T,\Delta\ll\hbar\omega_{c}$, we can further reduce Eq. (12) to

[TABLE]

with $\nu_{\chi}=N_{\chi}/N_{\phi}-g(n_{D}+1/2)$ being the filling factor. For $k_{B}T\ll\Delta$, the chemical potential can be approximated to be

[TABLE]

As can be seen, at half fillings $\nu_{\chi}/g=n_{\mathrm{F}}$, the low temperature chemical potential, pinned to the $n_{\mathrm{F}}$-th LL for relative small $\Delta(T,\nu_{\chi})$, is robust against the strain. At integer fillings $\nu_{\chi}/g=n_{\mathrm{F}}\pm 1/2$, however, Eq. (14) predicts an unphysical diverging chemical potential if $\Delta(T,\nu_{\chi})\neq 0$. This implies that the superconductivity must be fully suppressed for integer filling factors, i.e., $\Delta(T,\nu_{\chi})=0$ for $\nu_{\chi}/g=n_{\mathrm{F}}\pm 1/2$.

As analyzed, the superconducting gap and chemical potential are interactive, so that, to determine the superconductivity, the chemical potential must be accounted self-consistently into the equation of the superconducting gap. gap equation, as defined in Eq. (59), is obtained, self-consistently, by

[TABLE]

with $\overline{N}_{\phi}=N_{\phi}/A$ as the number of flux quanta per unit area. Substituting the retarded Green’s function, Eq. (57), into Eq. (15) leads to

[TABLE]

Combining Eqs. (14) and (16), we obtain for the zero temperature gap, in the weak coupling regime, as

[TABLE]

where

[TABLE]

is a constant, and $x=|U|g\overline{N}_{\phi}/(\hbar\omega_{c})$ is a dimensionless coupling strength of the attractive interaction. As shown by Eq. (17), $\Delta(0,\nu_{\chi})=0$ if $\nu_{\chi}/g=n_{\mathrm{F}}\pm 1/2$, which confirms the inference that the superconductivity is fully suppressed for integer filling factors. By substituting Eq. (17) into Eq. (14), the zero temperature chemical potential is obtained as

[TABLE]

As it shows, the zero temperature chemical potential at half fillings is robust to the strain for relative small $x$. At half fillings, the electron and hole density of states (DOSs) at the Fermi level are maximal, which is optimum for the emergence of superconductivity. Consequently, the superconducting gap reaches its maximum for half filling factors, as shown by Eq. (17). Interestingly, the chemical potential, away from the half fillings $|\nu_{\chi}/g-n_{\mathrm{F}}|>0$, is linearly scaled with the coupling strength $x$. With $|\nu_{\chi}/g-n_{\mathrm{F}}|$ increasing from [math], $\mu_{\chi}(0,\nu_{\chi})$ will shift away from the $n_{\mathrm{F}}$-th LL and the DOSs at the Fermi level decrease rapidly, which is unfavorable for the formation of Cooper pairs. As a result, $\Delta(0,\nu_{\chi})$ diminishes rapidly as the chemical potential shifts away from half fillings. Exactly at integer fillings, the superconductivity, as indicated by Eq. (17), is fully suppressed, i.e., $\Delta(0,\nu_{\chi})=0$. As a consequence, the electronic states, in the weak interaction regime, can not form the superconducting condensation at integer fillings.

The above analysis can be easily verified by the numerical results displayed in Fig. 1, which is calculated by Eqs. (8) and (15) self-consistently. As shown by the dark curve in Fig. 1(a), for a fixed $x$, the zero temperature chemical potential depends linearly on factor $\nu_{\chi}/g-n\in[-1/2,1/2)$, which characterizes the width of peaks of the zero temperature gap around half fillings, as shown by the red curve in Fig. 1(a). The numerical results are consistent with the analytical ones given by Eq. (19). As shown in Fig. 1(b), for small coupling strength $x\rightarrow 0$, the zero temperature gap peaks around half fillings $\nu_{\chi}/g=n$. The peaks are separated by zero plateaus in the vicinities of integer fillings $\nu_{\chi}/g=n\pm 1/2$, where formation of Cooper pairs is suppressed. With increasing the coupling strength, the zero plateaus reduce in width, and meanwhile, the peaks increase in height rapidly. If the coupling is strong enough $\Delta\rightarrow\hbar\omega_{c}$, transitions between different LLs become allowable, such that the zero plateaus could be fully filled by the peaks, leading to the emergence of superconductivity for integer fillings. Therefore, there exists a quantum critical coupling strength for the emergence of superconductivity at integer fillings.

The quantum critical coupling strength can be derived from the expression (17) for the zero temperature gap. For integer fillings, at which the chemical potential sits halfway between the two LLs nearest to the Fermi level, i.e., $\varepsilon_{n_{\mathrm{F}}-1,\chi}=-\varepsilon_{n_{\mathrm{F}},\chi}$, the zero temperature energy gap is derived to be

[TABLE]

where $\nu_{\chi}^{\mathrm{I}}=g(n_{\mathrm{F}}-1/2)$,

[TABLE]

is the critical coupling strength, and

[TABLE]

As can be seen, the critical coupling strength is determined by the spacing between the LLs nearest to the Fermi level. For $\Delta_{0}=\lambda_{\mathrm{so}}=0$, corresponding to case of a graphene sheet, the spacing of the LLs

[TABLE]

distributes symmetrically with respect to the $n=0$ LL, and decreases monotonously with increasing $|n|$. However, in the presence of finite $\Delta_{0}$ and $\lambda_{\mathrm{so}}$, corresponding to the MoS2 monolayer, the [math]-th LL $\varepsilon_{0,\chi}$ shifts away from the symmetry point $\mu=0$, which results in the asymmetric spacings between the $n=0$ and $n=\pm 1$ LLs, i.e., $|\varepsilon_{1,\chi}-\varepsilon_{0,\chi}|>$ $|\varepsilon_{0,\chi}-\varepsilon_{-1,\chi}|$, as demonstrated in Fig. 1(a). As a result, the critical coupling strength for $\nu_{\chi}/g=1/2$ is the largest. In other words, as $x$ increases from [math], the zero plateaus in Fig. 1(b) will vanish first for the higher LLs, later for $\nu_{\chi}/g=-1/2$, and last for $\nu_{\chi}/g=1/2$.

For finite temperatures, there exists a critical temperature $T_{c}$, above which the superconductivity vanishes, i.e., $\Delta(T\geq T_{c},\nu_{\chi})=0$. In Fig. 2, we plot the phase diagram of superconducting transition for strained monolayer MoS2 in the $T-x$ parameter space. In the critical regime $T\rightarrow T_{c}\gg\Delta$, by using the Poisson sum formula

[TABLE]

with $\omega_{m}=(2m+1)\pi k_{B}T$ and $\delta\rightarrow 0^{+}$, we can convert the gap equation to be

[TABLE]

For $T\rightarrow T_{c}$, $\Delta(T_{c},\nu_{\chi})\rightarrow 0$, we can expand the gap equation with respect to $\Delta^{2}(T_{c},\nu_{\chi})$. To the first order in $\Delta^{2}(T_{c},\nu_{\chi})$, we arrive at

[TABLE]

By using Eq. (16) and replacing the summation in Eq. (26) with an integral $\sum_{m}\rightarrow\frac{1}{k_{B}T}\int\frac{d\epsilon}{2\pi}$, we finally derive the superconducting gap at half filling factors to be

[TABLE]

with $\nu_{\chi}^{\mathrm{H}}=$ $gn_{\mathrm{F}}$, where the critical temperature is $T_{c}=\Delta(0,\nu_{\chi}^{\mathrm{H}})/(2k_{B})$. Consequently, in the weak coupling limit, $T_{c}\sim\hbar\upsilon_{\mathrm{F}}x/(4k_{B}l_{B})=g|U|eB/(4hk_{B})$, is linearly scaled with the interaction strength $U$ and the amount of strain $B$, which is similar to the case in strained graphene Uchoa and Barlas (2013). In fact, at partial filling of the LLs, $T_{c}\propto x$ in the $x\rightarrow 0$ limit, while, at the integer fillings, the transition is quantum critical below the critical coupling $x_{c}$, as seen in Fig. 2(b). The behavior of the critical coupling here is different from that for strained graphene in Ref. Uchoa and Barlas (2013). For strained graphene, the critical coupling strength $x_{c}$ distributes symmetrically with respect to the $n=0$ LL, due to the symmetrically-distributed LLs. For strained MoS2 monolayer, however, the [math]-th LL no longer locates at the symmetry point $\mu=0$, and as a result, the critical coupling is asymmetric for $\nu_{\chi}=1/2$ and $\nu_{\chi}=-1/2$, as reflected by the critical temperature in Fig. 2(b). From Fig. 2(b), we can also find that the critical coupling decreases monotonously with the LLs’ spacing, which is consistent with the analytical result presented in Eq. (21). The linearly scaled property of the critical temperature with $x$ is distinct from that in conventional weak coupling superconductors, where $T_{c}\sim\exp(-1/x)$ decreases exponentially with the effective coupling. With the coupling $x$ further increased, the system would cross over to the strong coupling regime $T_{c}\gtrsim\hbar\omega_{c}$.

IV Topological Index and Phase Diagram

In the presence of the Zeeman field, the TR symmetry is broken for the system. The BdG Hamiltonian without TR symmetry belongs to class DWang and Wu (2016); Schnyder et al. (2008) and is characterized in 2D by a topological invariant $C_{1}$, termed as the first Chern number. The topological index $C_{1}^{\xi\sigma}$ here can be calculated using the Kubo formulaZhang et al. (2014)

[TABLE]

where $\upsilon_{x(y),\xi\sigma}=\hbar^{-1}\partial\widetilde{\mathcal{H}}_{\mathrm{BdG}}/\partial k_{x(y)}$, $f_{n,\lambda}^{\xi\sigma}=f(E_{n,\lambda}^{\xi\sigma})$ represents the Fermi-Dirac distribution function, and

[TABLE]

[TABLE]

are eigenenergy and wavefunction of Hamiltonian (4), respectively, with $\lambda=\pm$ being the band index and

[TABLE]

By using the expressions above, we can derive the topological index to be

[TABLE]

where $f_{n}^{\xi\sigma}=f_{n,+}^{\xi\sigma}+f_{n,-}^{\xi\sigma}$. The first line of Eq. (32) vanishes after the summation over $n$ and, finally, $C_{1}^{\xi\sigma}$ can be rewritten as

[TABLE]

At zero or very low temperatures and for the case of $m_{\mathrm{z}}=0$, according to Eq. (29), $f_{n,+}^{\xi\sigma}=0$ and $f_{n,-}^{\xi\sigma}=1$, such that $C_{1}^{\xi\sigma}=0$ always satisfies. However, for a finite Zeeman field, $m_{\mathrm{z}}\neq 0$, by tuning the parameters, $E_{n,+}^{\xi\sigma}<0$ or $E_{n,-}^{\xi\sigma}>0$ could occur, which changes $C_{1}^{\xi\sigma}$ from zero to an integer, with the phase boundaries determined by $m_{\mathrm{z}}^{2}-(\mu_{\xi\sigma}-\Omega_{\xi\sigma}^{n})^{2}=\Delta^{2}$. For example, when $n_{\mathrm{F}}=0$, the topological index reduces to

[TABLE]

such that $C_{1}^{\xi\sigma}=\xi$ ($-\xi$) if $E_{0,+}^{\xi\sigma}<0$ ($E_{0,-}^{\xi\sigma}>0$). For $|m_{\mathrm{z}}|>\sqrt{(\mu+\frac{\Delta_{0}}{2}+\sigma\xi\lambda_{\mathrm{so}})^{2}+\Delta^{2}}$, $|C_{1}^{\xi\sigma}|=1$, which can also be seen from the numerical results displayed in Figs. 3 (a) and (b). As a result, the Chern number for $\xi$ valley $C_{1}^{\xi}=C_{1}^{\xi\uparrow}+C_{1}^{\xi\downarrow}$ will be $|C_{1}^{\xi}|=1$, when $(\mu+\frac{\Delta_{0}}{2}+\sigma\xi\lambda_{\mathrm{so}})^{2}<m_{\mathrm{z}}^{2}-\Delta^{2}<(\mu+\frac{\Delta_{0}}{2}-\sigma\xi\lambda_{\mathrm{so}})^{2}$, as shown in Figs. 3 (a) and (b). Odd values of the Chern number $|C_{1}^{\xi}|$ for the BdG Hamiltonian implies the emergence of Majorana modes for valley $\xi$. With tuning the chemical potential, the filling factors will change, and larger values of $C_{1}^{\xi\sigma}$ will emerge. However, the difference between $C_{1}^{\xi\uparrow}$ and $-C_{1}^{\xi\downarrow}$ only can be [math] or $\pm 1$, as shown in Fig. 3(c). Therefore, $|C_{1}^{\xi}|=0$ and $1$ occur alternately in the parameter space, as seen from Fig. 3 (c). Due to the valley symmetry $E_{n,\overline{\lambda}}^{\overline{\xi}\overline{\sigma}}=-E_{n,\lambda}^{\xi\sigma}$, $C_{1}^{\xi\sigma}=C_{1}^{\overline{\xi}\overline{\sigma}}$. As a result, the total Chern number for the system $C_{1}=C_{1}^{+}+C_{1}^{-}=2C_{1}^{\xi}$ is even, as illustrated in Fig. 3 (c). In other words, a pair of Majorana modes emerge simultaneously in two valleys, when $C_{1}=\pm 2$. The realization of odd values of the total Chern number for the present system requires to further break the valley symmetry. For example, when the Rashba spin-orbit interaction is taken into account, the intervalley coupling would destroy the valley symmetry Wang and Wu (2016), and then we can expect the emergence of a single Majorana mode for the edges of the present system.

We have assumed that the $s$-wave pairing for the convenience of calculation. As proposed in Ref. Uchoa and Barlas (2013), in the presence of substrates, superconductivity can be triggered by conventional electron-phonon coupling, meaning that the $s$-wave pairing can possibly be realized in such a system. Our conclusions, being insensitive to the phase information of the pairing potential, may not be limited to the s-wave pairing. There could possibly be other mechanismKotov et al. (2012); Uchoa and Castro Neto (2007); Black-Schaffer and Doniach (2007); Honerkamp (2008); Kopnin et al. (2011), such as density waveHonerkamp (2008), leading to the superconductivity, but it may not affect our discussions on the topological properties of the superconductivity phase.

Theoretically, a topological superconductor in two-dimensions with odd integer Chern numbers is predicted to host topologically protected gapless chiral Majorana edge modes. In experiments, while intensive efforts have been made to search for the chiral Majorana edge modes, the exclusive signature for such exotic fermions is still under debate. Very recently, following the theoretical proposal in Ref.Wang et al. (2015), He $et$ $al.$ He et al. (2017) have observed the characteristic half-integer longitudinal conductance. However, the half-integer longitudinal conductance, as argued by Ji $et$ $al.$ Ji and Wen (2018), is only a necessary condition for identifying the Majorana edge modes, but not a sufficient condition. Here, we propose that the Majorana edge modes can be realized by splitting of the pseudo LLs. The strained MoS2 monolayer can be realized by using the method proposed in Ref. Levy et al. (2010). By placing the grown sample on the setup proposed in Refs. Wang et al. (2015); He et al. (2017), we can expect to observe an integer to half-integer transition of the longitudinal conductance, when the Zeeman field is turned on gradually, for an appropriate chemical potential. The half-integer longitudinal conductance will be a valuable signature, but possibly not an exclusive evidence, of the Majorana edge modes. Searching for unambiguous experimental fingerprint of the Majorana edge modes is still a challenging task at the research front in the condensed-matter physics.

V summary

In summary, we have investigated the superconductivity and topological properties in strained dichalcogenides. We generalized the superconducting theory for gapless graphene to dichalcogenides with an intrinsic band gap. It is found that superconductivity can emerge in the pseudo LLs induced by strain. In the weak coupling limit, the superconducting gap is linearly-scaled with the coupling strength for the partial fillings, in contrast to conventional weak coupling superconductors. The superconductivity gap is maximized when the LLs are half-filled, but for integer fillings the superconductivity is fully suppressed, when the coupling strength is below a quantum critical value. We find the quantum critical coupling strength is determined by the spacing between the two LLs closest to the Fermi level. Interestingly, in the presence of a Zeeman field, a pair of Majorana modes emerge simultaneously in the two valleys of strained dichalcogenides. A single Majorana mode can be realized, if the valley symmetry is further lifted.

VI acknowledgements

We thank Prof. Tao Zhou for helpful discussions. This work was supported by the State Key Program for Basic Researches of China under Grants No. 2015CB921202, and No. 2017YFA0303203 (D.Y.X), the National Natural Science Foundation of China under Grants No. 11674160 (L.S.), No. 11574155 (R.M.), No. 11474149 (R.S.), No. 11804130 (W.L.), No. 11474106 (R.-Q.W) and the Key Program for Guangdong NSF of China under Grant No. 2017B030311003 (R.-Q.W) and GDUPS (2017).

Appendix A Derivation for the number of particles and superconducting gap

equation

In the Hilbert space $\{\Phi^{(n)}\}$, with $\Phi^{(n)}=$ $\left(\psi_{k,\xi\uparrow}^{(n)},\psi_{-k,\overline{\xi}\downarrow}^{(n)\ast},\psi_{k,\xi\downarrow}^{(n)},\psi_{-k,\overline{\xi}\uparrow}^{(n)\ast}\right)$, the matrix elements of the BdG Hamiltonian are calculated by

[TABLE]

where $\psi_{-k,\xi\sigma}^{(n)\ast}$ is the wavefunction corresponding to the hole Hamiltonian $h_{-\mathbf{p},\xi\sigma}^{\ast}$. By using the relations below

[TABLE]

we can derive the matrix elements of the BdG Hamiltonian to be

[TABLE]

with $\Lambda_{n,\xi\sigma}=\Omega_{\xi\sigma}^{n}-\mu_{\xi\sigma}$. Since the spin index $\sigma$ is locked to the valley index $\xi$, $\Lambda_{n,\xi\sigma}=\Lambda_{n,\overline{\xi}\overline{\sigma}}$ and, by labeling $\chi=\xi\sigma$ for brevity, we can express the BdG Hamiltonian as presented in Eq. (7) of the paper.

To calculate the number of particles, we express the BdG Hamiltonian in the form of the second quantization as

[TABLE]

where $c_{n,k,\xi\sigma}^{{\dagger}}$ ($c_{n,k,\xi\sigma}$) creates (annihilates) an electron in the state $|\psi_{k,\xi\sigma}^{(n)}\rangle$. The number of particles now can be defined as $\mathcal{N}_{\xi\sigma}=\sum_{nk}\langle c_{n,k,\xi\sigma}^{{\dagger}}c_{n,k,\xi\sigma}\rangle$. According to the fluctuation-dissipation theorem Deng et al. (2016), i.e.,

[TABLE]

where $\langle\langle c_{n,k,\xi\sigma}|c_{n,k,\xi\sigma}^{{\dagger}}\rangle\rangle_{\omega}^{r}$ represents the retarded Green’s function and $f(\omega)=1/(1+e^{\omega/k_{B}T})$ is the Fermi-Dirac distribution function, the number of particles can be expressed as

[TABLE]

where $g=2$ is the valley$\times$spin degeneracy and $N_{\phi}=A/(2\pi l_{B}^{2})$ is degeneracy of the LL originating from the summation over $k$, with $A=L_{x}L_{y}$ as area of the sample. The retarded Green’s function, with respect to $\mathcal{H}_{\chi}$, is defined as

[TABLE]

with $E_{n}^{\chi}=\sqrt{\varepsilon_{n,\chi}^{2}+\Delta^{2}}$ and $\varepsilon_{n,\chi}=\Omega_{\chi}^{n}-\mu_{\chi}$, where we adopted the abbreviation $\omega^{+}=\omega+i0^{+}$. Therefore, we can obtain for

[TABLE]

If the deep-energy states $n<n_{\mathrm{F}}$ are fully occupied, with $n_{\mathrm{F}}$ labeling the highest occupied LL, we can reduce Eq. (58) to Eq. (12) in the paper. In the low temperature and weak interaction limit $T,\Delta\ll\hbar\omega_{c}$, we can approximate $\tanh[\sqrt{\varepsilon_{n,\chi}^{2}+\Delta^{2}}/(2k_{B}T)]\simeq 1$ for $n<n_{\mathrm{F}}$ and $\sqrt{\varepsilon_{n_{\mathrm{F}},\chi}^{2}+\Delta^{2}}\simeq|\varepsilon_{n_{\mathrm{F}},\chi}|$. Therefore, we can further reduce Eq. (12) to Eq. (13). For $k_{B}T\ll\Delta$, $\tanh[E_{n_{\mathrm{F}}}^{\chi}/(2k_{B}T)]\simeq 1$ and from Eq. (13), it is easy to solve the chemical potential, as given by Eq. (14).

The gap equation, similar to Eq. (8), is obtained self-consistently by the mean-field approximation, i.e.,

[TABLE]

with $\overline{N}_{\phi}=N_{\phi}/A$ as the number of flux quanta per unit area. By substituting the retarded Green’s function, Eq. (57), into the above equation, we can derive Eq. (16).