---
abstract: 'We study the possibility of topological superconductivity in the noncentrosymmetric monolayer and trilayer TaS$_2$ with out-of-plane mirror symmetry.  A gapless time-reversal invariant f+s-wave pairing state with even mirror parity is found to be a promising candidate.  This mixing state holds 12(36) nodes at the Fermi pockets around $\Gamma$ for monolayer(trilayer) case and its unconventional superconductivity is consistent with the STM experiments observed in 2$H$-TaS$_2$ thin flakes.  Furthermore, with doping or under uniaxial pressure for trilayer 2$H$-TaS$_2$, large-Chern-number time-reversal symmetry breaking mixing states between d+id- and p-ip-wave pairings can be realized in the phase diagram.' 
author:
- Weipeng Chen
- Jin An
bibliography:
- 'ref.bib'
title: 'Topological Ising pairing states in monolayer and trilayer TaS$_2$'
---

Introduction {#sec1}
============

Dimensionality plays an important role in superconducting layered materials, which may behave quite differently in two-dimensional(2D) limit.  While the layered transition-metal dichalcogenide(TMD) bulk materials exhibit s-wave superconductivity[@Clayman1971; @Hess1990; @Corcoran1994; @Boaknin2003; @Huang2007], always coexisting with a competing charge density wave(CDW)[@Berthier1976; @Mutka1983; @Neto2001; @Yokoya2001; @Guillamon2011; @Xi2015N], their thin flakes have been found to show many new phenomena including the spin-valley locking[@Xiao2012; @Lu2013; @Suzuki2014; @Bawden2016; @Wu2016; @Dey2017], the valley Hall effect[@Xiao2012; @Mak2014; @Lee2016; @Yu2016], and Ising pairing which is supported by a great enhancement of the in-plane upper critical magnetic field[@Lu2015; @Xi2015; @Saito2016; @Xing2017; @Lu2018; @Sohn2018; @Barrera2018].  Recent investigations have also revealed the possibility of unconventional superconductivity in 2D TMDs[@Yuan2014; @Zhou2016; @Sharma2016; @Hsu2017; @Mockli2018]. 

One special member of TMDs family is 2$H$-TaS$_2$, which has attracted much attentions in recent years.  In the detached flakes of 2$H$-TaS$_2$, a zero-bias conductance peak(ZBCP) was observed by STM at 0.15 K, indicating unconventional superconductivity[@Galvis2014].  This leads to a f-wave pairing scenario in atomic thin 2$H$-TaS$_2$ layers.  On the other hand, neither superconductivity nor CDW was observed in the monolayer $H$-TaS$_2$ at 4.7 K[@Sanders2016], which means a strong suppression of CDW in 2D limit of the material.  Critical temperature T$_c$ of 2$H$-TaS$_2$ is also found to be strongly enhanced from 0.54 K to about 2.1 K, as the material thickness is decreased from bulk to about five molecular layers[@Navarro2016].  Furthermore, T$_c$ of the monolayer could reach 3 K[@Barrera2018].  Although the enhancement of T$_{c}$ is believed to result from the suppression of the CDW via dimension reduction[@Navarro2016; @Yang2018], what pairing symmetry the superconducting state in 2$H$-TaS$_2$ thin flakes can be and whether it supports topological superconductivity remain unclear. 

From the viewpoint of the crystal structure of 2$H$-TaS$_2$, the thin flakes may behave differently for even or odd number of layers, since while 2$H$-TaS$_2$ bulk material always has an inversion center, only the even-layer thin films are still inversion symmetric.  Recent experiment has found the difference in gap structure between the trilayer and four-layer 2$H$-NbSe$_2$[@Dvir2018].  The odd-layer 2$H$-TaS$_2$ thin films are noncentrosymmetric but always hold out-of-plane mirror symmetry.  In this work, we will demonstrate that the Ising superconductivity in monolayer and trilayer TaS$_2$ should be unconventional, and even topologically nontrivial.  We find a time-reversal invariant gapless f+s-wave pairing state can be a candidate and furthermore, the time-reversal symmetry(TRS) breaking mixing states between d+id- and p-ip-wave pairings can be induced by doping or uniaxial pressure. 

The paper is organized as follows.  In Sec.\[sec2\], we introduce our model.  In Sec.\[sec3\], we analyze the pairing symmetry and classify the gap functions according to crystal symmetry.  We also include an appendix to show the technical details of how to make use of the irreducible representations of crystal symmetry to solve the multilayer linearized gap equations.  In Sec.\[sec4\], we first give the pairing phase diagram, and discuss the Ising pairing states.  Then we investigate the topological features of these superconducting states as well as their pairing symmetry transitions with doping or under uniaxial pressure.  Finally, we summarize our results in Sec.\[sec5\]. 

Model {#sec2}
=====

TMDs have weak van der Waals couplings among TX$_2$($T$=Ta,Nb and $X$=S,Se) layers.  An $H$-TX$_2$ monolayer contains three atomic ones each of which forms a triangular lattice, with one $T$ layer sandwiched between two $X$ ones. “ 2$H$" in 2$H$-TX$_2$ material means the neighboring TX$_2$ layers in the bulk take the $AbA-CbC$ stacking sequence, preserving the global $D_{6h}$ symmetry.  The thin film we study in this paper consists of one or three TaS$_2$ layers.  In this and next sections we focus the discussion of model and symmetry analysis on the trilayer case, as shown schematically in Fig.  \[fig01\](a), since that of the monolayer one is straightforward and parallel.  We model the system by simply replacing the three TaS$_2$ layers by three atomic Ta layers in $bbb$ sequence, because superconductivity is believed to occur only within Ta layers.  Generically, a thin layer of 2$H$-TaS$_2$ always has Fermi surface(FS) sheets centered at $\Gamma$ and K(K$^\prime$).  Although the FS stems respectively from $d_{z^2}$ and $d_{xy/x^2-y^2}$ orbitals[@Ge2012; @Noat2015; @Heil2017], we simply simulate the bands by a tight-binding model based only on $d_{z^2}$ orbital.  Here we ignore CDW as it is believed to be strongly suppressed in the thin 2$H$-TaS$_2$ layers[@Navarro2016; @Yang2018].  For a trilayer system, the tight-binding Hamiltonian can thus be given by, $$\begin{aligned}
    H_0(\bm{k})=
    \begin{bmatrix}
    h_1(\bm{k}) & t_{\perp} & 0 \\
    t_{\perp} & h_2(\bm{k}) & t_{\perp}  \\
    0 & t_{\perp} & h_3(\bm{k})
    \end{bmatrix}
    .\end{aligned}$$ Here $t_{\perp}$ and $h_l(\bm{k})$ represent the interlayer hopping integral and the intralayer Hamiltonian for layer $l$($l=1,2,3$).  The latter can be expressed as $$\begin{aligned}
    h_l(\bm{k})=\epsilon (\bm{k})+h^{SO}_{l}(\bm{k})+h^R_{l}(\bm{k})
    ,\end{aligned}$$ where $\epsilon (\bm{k})=-\sum\limits_{j=1}^{3} \{2t_1 \cos k_{j}
    +2t_2 \cos(k_{j}-k_{j+1})\}  -\mu$ with $k_{j}=\bm{k}\cdot\bm{R_{j}}$, and the unit lattice vectors $\bm{R}_1=\bm{y}$, $\bm{R}_2=-\bm{y}/2-\sqrt{3}\bm{x}/2$, $\bm{R}_3=-\bm{y}/2+\sqrt{3}\bm{x}/2$ and $\bm{R}_4\equiv\bm{R}_1$. Here $\mu$ is the chemical potential, $t_1$, $t_2$ denote the nearest-neighbor(NN) and next-nearest-neighbor(NNN) hoping integrals.  We set $(t_1,t_2,t_\perp,\mu)=(-60,-140,-40,0)$ meV to fit the electronic band structure from experiments and DFT calculations [@Navarro2016; @Sanders2016; @Zhao2017; @Barrera2018]. 

![ \[fig01\](color online).(a)The side and top views of atomic structure of trilayer 2$H$-TaS$_2$. Three basis vectors $\bm{R}_j$ are given by the red arrows.  The band structures of (b)monolayer $H$-TaS$_2$ and (c)trilayer 2$H$-TaS$_2$. The inset is a blowup of the dotted square near $M$. (d)Spin-valley locking in trilayer 2$H$-TaS$_2$, where the arrows denote the effective magnetic field $B_{so}$ acting on the two valleys $K$ and $K^{'}$ of each layer.](Fig1.eps){width="7.cm"}

As depicted in Fig.  \[fig01\](a), the stacking form of sulfur atoms breaks the inversion symmetry in each $H$-TaS$_2$ layer, resulting in the Ising spin-orbit coupling(SOC).  Moreover, the middle TaS$_2$ layer can be viewed as being rotated by $180^\circ$ with respect to the top or bottom molecular ones.  This leads to a layer-dependent intrinsic SOC given below[@Youn2012; @Xi2015], $$\begin{aligned}
    h^{SO}_{l}(\bm{k})=(-1)^l\beta(\bm{k})\sigma_z,
    \end{aligned}$$ where $\beta(\bm{k})=\beta_{so}(\sin{k_1}+\sin{k_2}+\sin{k_3})$. For electrons near $K$($K^{'}$), this SOC can be viewed as an effective valley-dependent out-of-plane magnetic field $B_{so}$, favoring spin-up(-down) electrons and alternating among layers, as schematically shown in Fig.  \[fig01\](d).  This gives rise to the well-known spin-valley locking effect[@Xiao2012; @Lu2013; @Suzuki2014; @Bawden2016; @Wu2016; @Dey2017].  The magnitude of $B_{so}$ is estimated to be quite large, about 3000T for TaS$_2$, since the spin splitting at $K$($K'$) is about 348 meV[@Sanders2016], from which $\beta_{so}$ is fixed to be $\beta_{so}=67$ meV. Ising pairing occurs when Cooper pairs are formed between electrons around $K$ and $K^{'}$, spin oriented oppositely but along the out-of-plane direction.  Ising superconductivity in TMDs has been supported experimentally by the significant enhancement of the in-plane upper critical magnetic field[@Lu2015; @Xi2015; @Saito2016; @Xing2017; @Liu2017; @Lu2018; @Sohn2018; @Barrera2018]. 

Rashba SOC is also taken into account.  It takes the following form, $$\begin{aligned}
    h^R_{l}(\bm{k})=\alpha^{l}\bm{g(k)}\cdot \bm{\sigma}
    ,\end{aligned}$$ where $\bm{g}(\bm{k})=(-Im [\lambda(\bm{k})],-Re[\lambda(\bm{k})],0)$, with $\lambda (\bm{k})=\sin{k_1}+e^{-i\frac{2\pi}{3}}\sin{k_3}+e^{-i\frac{4\pi}{3}}\sin{k_3}$. Since the Rashba SOC is strongly dependent on the surface normal, the coupling constants should be layer-dependent, which are assumed to be $(\alpha^{1},\alpha^{2},\alpha^{3})=(\alpha_R,0,-\alpha_R)$ for a free-standing trilayer 2$H$-TaS$_2$, with $\alpha_R$ the coupling strength.  This kind of coupling also appears in the layered system with superlattice structures[@Mizukami2011; @Goh2012; @Shimozawa2014].  The Rashba SOC favors electrons’ spin oriented along the in-plane directions, which means it will compete with the Ising SOC and weaken the effect of spin-valley locking and hence Ising pairing.  As long as $\alpha_R\ll\beta_{so}$, Ising pairing is expected to be still dominant. 

The band structures for monolayer $H$-TaS$_2$ and trilayer 2$H$-TaS$_2$ without Rashba SOC are shown in Fig.  \[fig01\](b) and (c), respectively.  For monolayer $H$-TaS$_2$, the band is doubly degenerate at $\Gamma$ but spin split at $K$($K'$), as expected.  For trilayer 2$H$-TaS$_2$, there are four bands in total but two of the four are doubly degenerate due to accidental degeneracy.  The degeneracy can be lifted when the middle layer takes a slightly different value of $\beta_{so}$ to the top/bottom layer, which still preserves the crystal symmetry of trilayer 2$H$-TaS$_2$.

Pairing Interaction and Symmetry Analysis {#sec3}
=========================================

The trilayer 2$H$-TaS$_2$ system belongs to crystal $D_{3h}$, while each TaS$_2$ monolayer(especially the top/bottom layer) only has $C_{3v}$ symmetry.  Since no pairing occurs between layers, the pairing basis gap functions $\Delta^{\Gamma,\alpha}_{i}(\bm{k})$ can thus be classified according to the irreducible representations(IRs) of $C_{3v}$[@Yuan2014], as shown in Table \[tab1\].  According to whether the total spin $S$ of the Cooper pair is $1$ or $0$, $\Delta^{\Gamma,\alpha}_{i}(\bm{k})$ takes the following two different kinds of forms: $$\begin{aligned}
\Delta^{\Gamma,\alpha}_{i}(\bm{k})=
    \begin{cases}
    \Psi^{\Gamma,\alpha}_{i} i\sigma_y & {singlet \ pairing}\\
    \bm{d}^{\Gamma,\alpha}_{i} \cdot \bm{\sigma} i \sigma_y & {triplet \ pairing}
    \end{cases}
    ,\end{aligned}$$ where $\Psi^{\Gamma,\alpha}_{i}$ is the singlet order parameter while $\bm{d}^{\Gamma,\alpha}_{i}$ denotes the director for the triplet pairing.  Here $\Gamma$ represents the IR of $C_{3v}$ symmetry with dimension $d^{\Gamma}$ and it can be $A_{1}(d^{A_{1}}=1)$, $A_{2}(d^{A_{2}}=1)$ or $E(d^{E}=2)$. The index $\alpha$ is used to distinguish different types of representations in the same IR $\Gamma$, and $i$ runs from $1$ to $d^{\Gamma}$. As an example, for $\Gamma=A_1$, $\alpha$ can be ‘on’, ‘nn’, ‘z’ and ‘xy’.  All these $\Delta^{\Gamma,\alpha}_{i}(\bm{k})$ are so normalized $\frac{2}{N}\sum_{\bm{k}}\mid\Psi^{\Gamma,\alpha}_{i}\mid^{2}=1$(or $\frac{2}{N}\sum_{\bm{k}}\mid\bm{d}^{\Gamma,\alpha}_{i}\mid^{2}=1$) as to satisfy the following orthogonal relations between IRs: $$\begin{aligned}
\label{eqn6}
    \frac{1}{N}\sum_{\bm{k}} Tr\{\Delta^{ \dagger \Gamma^\prime,\alpha}_{i}(\bm{k}) \Delta^{\Gamma,\beta}_{j}(\bm{k})\}=\delta_{\Gamma^\prime \Gamma}\delta_{\alpha \beta}\delta_{ij}
    ,\end{aligned}$$ where N is the total number of unit cells in the system.  See Appendix for the detail. 

Although all the basis gap functions in Table \[tab1\] can be realized in principle in materials with crystal $C_{3v}$ symmetry, the in-plane triplet pairing components such as $\bm{d}^{A_1,xy}$, $\bm{d}^{A_2,xy}$, $\bm{d}^{E,xy}$, $\bm{d}^{E,\widetilde{xy}}$ actually never occur as long as Ising pairing is dominant in 2$H$-TaS$_2$, which is guaranteed by the strong intrinsic SOC.  This is actually also confirmed by our numerical calculations in Sec.\[sec4\]. 

Without considering the pairing mechanism, the pairing term for a multilayer 2$H$-TaS$_2$ can be generically written as: $$\begin{aligned}
    H_{int}&=&
    \frac{1}{2}\sum V_{s_1s_2s_3s_4}(\bm{k},\bm{k}^\prime)
    c^{\dagger}_{l,\bm{k}s_1}c^{\dagger}_{l,-\bm{k}s_2}\nonumber\\
    &&\times c_{l,\bm{k}^{\prime}s_3}c_{l,-\bm{k}^{\prime}s_4}
    ,\end{aligned}$$ where the sum runs over the repeated indices.  $V_{s_1s_2s_3s_4}(\bm{k},\bm{k}^\prime)$ is the pairing interaction and $c^{\dagger}_{l,\bm{k}s}$ the creation operator for an electron with spin $s$ at layer $l$. According to $C_{3v}$ symmetry and taking the on-site and NN parings into account, $V_{s_1s_2s_3s_4}(\bm{k},\bm{k}^\prime)$ can be expanded as follows,

$$\begin{aligned}
    &&V_{s_1s_2s_3s_4}(\bm{k},\bm{k}^\prime) \nonumber \\
    =&&v_0\Psi^{A_{1},on}\Psi^{*A_{1},on}(i\sigma_y)_{s_1s_2}(i\sigma_y)_{s_3s_4}+  \nonumber \\
    &&v_1 \sum_{\Gamma,i}\Psi^{\Gamma,nn}_{i} \Psi^{*\Gamma,nn}_{i}(i\sigma_y)_{s_1s_2}(i\sigma_y)_{s_3s_4}+  \nonumber \\
    &&v_1 \sum_{\Gamma,\alpha,i} [\bm{d}^{\Gamma,\alpha}_{i}\cdot \bm{\sigma}i\sigma_y]_{s_1s_2}[\bm{d}^{\Gamma,\alpha}_{i}\cdot \bm{\sigma}i\sigma_y]^*_{s_3s_4}
    ,\end{aligned}$$

where $v_0(v_1)$ is the pairing strength for the on-site(NN) pairing interaction. 

[cccc]{} $\Gamma$ && Singlet & Triplet\
$A_1$ &&

  ---------------------------------------
   $\Psi^{A_1,on}_{}=\frac{1}{\sqrt{2}}$
       $\Psi^{A_1,nn}_{}=C(\bm{k})$
  ---------------------------------------

  : Classification of spin-singlet and spin-triplet basis gap functions based on $C_{3v}$ symmetry, where $C(\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+\cos{k_2}+\cos{k_3}), C_+ (\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+e^{i\frac{2\pi}{3}}\cos{k_2}+e^{i\frac{4\pi}{3}}\cos{k_3}), S(\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+\sin{k_2}+\sin{k_3}), S_+ (\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+e^{i\frac{2\pi}{3}}\sin{k_2}+e^{i\frac{4\pi}{3}}\sin{k_3}), C_- (\bm{k})=C_+^* (\bm{k}), S_- (\bm{k})=S_+^* (\bm{k})$, with $\bm{x}_\pm=(\bm{x}\pm i\bm{y})/2$.

&

  ------------------------------------------------------------------
                 $\bm{d}^{A_1,z}_{}=S(\bm{k})\bm{z}$
   $\bm{d}^{A_1,xy}_{}=i[S_-(\bm{k})\bm{x}_+-S_+(\bm{k})\bm{x}_-$\]
  ------------------------------------------------------------------

  : Classification of spin-singlet and spin-triplet basis gap functions based on $C_{3v}$ symmetry, where $C(\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+\cos{k_2}+\cos{k_3}), C_+ (\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+e^{i\frac{2\pi}{3}}\cos{k_2}+e^{i\frac{4\pi}{3}}\cos{k_3}), S(\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+\sin{k_2}+\sin{k_3}), S_+ (\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+e^{i\frac{2\pi}{3}}\sin{k_2}+e^{i\frac{4\pi}{3}}\sin{k_3}), C_- (\bm{k})=C_+^* (\bm{k}), S_- (\bm{k})=S_+^* (\bm{k})$, with $\bm{x}_\pm=(\bm{x}\pm i\bm{y})/2$.

\
$A_2$ && & $\bm{d}^{A_2,xy}_{}=S_-(\bm{k})\bm{x}_++S_+(\bm{k})\bm{x}_-$\
E & &

  ------------------------------------------------------------------------------------------
   $\begin{cases} \Psi^{E,nn}_{1}=C_+(\bm{k})    \\ \Psi^{E,nn}_{2}=C_-(\bm{k})\end{cases}$
  ------------------------------------------------------------------------------------------

  : Classification of spin-singlet and spin-triplet basis gap functions based on $C_{3v}$ symmetry, where $C(\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+\cos{k_2}+\cos{k_3}), C_+ (\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+e^{i\frac{2\pi}{3}}\cos{k_2}+e^{i\frac{4\pi}{3}}\cos{k_3}), S(\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+\sin{k_2}+\sin{k_3}), S_+ (\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+e^{i\frac{2\pi}{3}}\sin{k_2}+e^{i\frac{4\pi}{3}}\sin{k_3}), C_- (\bm{k})=C_+^* (\bm{k}), S_- (\bm{k})=S_+^* (\bm{k})$, with $\bm{x}_\pm=(\bm{x}\pm i\bm{y})/2$.

&

  -----------------------------------------------------------------------------------------------------
   $\begin{cases} \bm{d}^{E,z}_{1}=S_+(\bm{k})\bm{z}\\ \bm{d}^{E,z}_{2}=S_-(\bm{k})\bm{z}\end{cases} $

                      $\begin{cases}\bm{d}^{E,xy}_{1}=\sqrt{2}S(\bm{k})\bm{x}_+ \\
                                \bm{d}^{E,xy}_{2}=-\sqrt{2}S(\bm{k})\bm{x}_-\end{cases}$

               $\begin{cases}\bm{d}^{E,\widetilde {xy}}_{1}=\sqrt{2}S_-(\bm{k})\bm{x}_-\\
                         \bm{d}^{E,\widetilde {xy}}_{2}=-\sqrt{2}S_+(\bm{k})\bm{x}_+\end{cases}$
  -----------------------------------------------------------------------------------------------------

  : Classification of spin-singlet and spin-triplet basis gap functions based on $C_{3v}$ symmetry, where $C(\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+\cos{k_2}+\cos{k_3}), C_+ (\bm{k})=\frac{1}{\sqrt{3}}(\cos{k_1}+e^{i\frac{2\pi}{3}}\cos{k_2}+e^{i\frac{4\pi}{3}}\cos{k_3}), S(\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+\sin{k_2}+\sin{k_3}), S_+ (\bm{k})=\frac{1}{\sqrt{3}}(\sin{k_1}+e^{i\frac{2\pi}{3}}\sin{k_2}+e^{i\frac{4\pi}{3}}\sin{k_3}), C_- (\bm{k})=C_+^* (\bm{k}), S_- (\bm{k})=S_+^* (\bm{k})$, with $\bm{x}_\pm=(\bm{x}\pm i\bm{y})/2$.

\

\[tab1\]

Because the trilayer 2$H$-TaS$_2$ system lacks an inversion center, generically there exists mixing between singlet and triplet pairings.  The corresponding $6\times6$ pairing gap function $\Delta(\bm{k})$ for the superconducting trilayer system can thus be expanded in the basis of $\Delta^{\Gamma,\alpha}_{i}(\bm{k})$ as, $$\begin{aligned}
\label{eqn9}
    \Delta(\bm{k})/\Delta&\equiv&
    \begin{bmatrix}
    \Delta_{l=1}(\bm{k})&& \\
    &\Delta_{l=2}(\bm{k})& \\
    &&\Delta_{l=3}(\bm{k}) \\
    \end{bmatrix} \nonumber
    \\
    &=&\sum\limits_{\alpha,i}
    b^{\Gamma,\alpha}_{i} \begin{bmatrix}
    \chi^{\Gamma,\alpha}_{l=1,i} && \\
    &\chi^{\Gamma,\alpha}_{l=2,i}& \\
    &&\chi^{\Gamma,\alpha}_{l=3,i}\\
    \end{bmatrix}\Delta^{\Gamma,\alpha}_{i}(\bm{k})
    ,\end{aligned}$$ where $\Delta$ is the gap value, and the expansion is made only for a definite IR $\Gamma$, as the mixing among different IRs is forbidden[@Sigrist1991].  Here $\Delta_{l}(\bm{k})$ is the $2\times2$ gap matrix for layer $l$. The expansion coefficients $b^{\Gamma,\alpha}_{i}$ can be normalized as $\sum \limits_{\alpha,i} |b^{\Gamma,\alpha}_{i}|^{2}=1$, while the relative paring amplitudes $\chi_{l,i}^{\Gamma,\alpha}$ are always normalized as $\sum \limits_l |\chi_{l,i}^{\Gamma,\alpha}|^{2}=1$. When the most energetically favorable pairing gap function $\Delta(\bm{k})$ is found, the BdG Hamiltonian of this superconducting trilayer system takes the standard form, $$\begin{aligned}
    H_{BdG}(\bm{k})=
    \begin{bmatrix}
    H_0(\bm{k}) & \Delta (\bm{k}) \\
    \Delta^\dagger (\bm{k}) & -H_0^T(\bm{-k})
    \end{bmatrix}
    .\end{aligned}$$

This $D_{3h}$ trilayer system also has an out-of-plane mirror symmetry, with the mirror plane lying in the middle Ta layer.  For the normal state, we have $[M_{xy},H_0(\bm{k})]=0$, with $$\begin{aligned}
    M_{xy}=
    \begin{bmatrix}
    0  & 0 & 1\\
    0  & 1 & 0\\
    1  & 0 & 0\\
    \end{bmatrix}
    i\sigma_z
    ,\end{aligned}$$ which exchanges the top layer with bottom layer and reverses the in-plane spins.  For the superconducting state, the mirror symmetry requires $[M^{\pm},H_{BdG}(\bm{k})]=0$, with mirror operators $M^{\pm}=\left(
    \begin{array}{cc}
    M_{xy}& 0\\
    0& \pm M_{xy}^* \\
    \end{array}
    \right)$, leading to a requirement that $$\begin{aligned}
    M_{xy}\Delta(\bm{k})M_{xy}^T=\pm \Delta(\bm{k}),
    \end{aligned}$$ where the symbol “+(-)” represents even(odd) mirror parity[@Yoshida2015].  Therefore in the trilayer 2$H$-TaS$_2$ system, for an even-mirror-parity state, its singlet components or triplet ones with $\bm{d}^{\Gamma}\parallel\bm{z}$, and for an odd-mirror-parity state, its triplet components with $\bm{d}^{\Gamma}\perp\bm{z}$, satisfy, $$\begin{aligned}
    \chi^{\Gamma,\alpha}_{l=1,i}(\bm{k})=\chi^{\Gamma,\alpha}_{l=3,i}(\bm{k}). 
    \end{aligned}$$ For an odd-mirror-parity state, its singlet components or triplet ones with $\bm{d}^{\Gamma}\parallel\bm{z}$, and for an even-mirror-parity state, its triplet components with $\bm{d}^{\Gamma}\perp\bm{z}$, satisfy instead the requirements, $$\begin{aligned}
    \begin{cases}
    \chi^{\Gamma,\alpha}_{l=1,i}(\bm{k})=-\chi^{\Gamma,\alpha}_{l=3,i}(\bm{k})=\frac{1}{\sqrt{2}} \\
    \chi^{\Gamma,\alpha}_{l=2,i}(\bm{k})=0
    \end{cases}
    .\end{aligned}$$ These symmetry properties are summarized in Table \[tab2\]. 

[p[1cm]{}&lt;p[3cm]{}&lt; p[4cm]{}&lt;]{} Mirror Parity&Components&$\bm{\chi}^{\Gamma,\alpha}_{i}\equiv (\chi^{\Gamma,\alpha}_{1,i},\chi^{\Gamma,\alpha}_{2,i},\chi^{\Gamma,\alpha}_{3,i})$\
Even &

  ----------------------------------------------------------------------
   $\psi^{\Gamma,\alpha}_{i}$ or $\bm{d}_{i}^{\Gamma,\alpha} || \bm{z}$
                $\bm{d}_{i}^{\Gamma,\alpha} \perp \bm{z}$
  ----------------------------------------------------------------------

  : Symmetry of relative pairing amplitudes $\bm{\chi}^{\Gamma,\alpha}_{i}$ for a pairing gap function of the trilayer 2$H$-TaS$_2$, classified by the out-of-plane mirror symmetry. 

&

  ---------------------------------------------------------
   $\chi^{\Gamma,\alpha}_{1,i}=\chi^{\Gamma,\alpha}_{3,i}$
      ($\frac{1}{\sqrt{2}}$, 0, -$\frac{1}{\sqrt{2}})$
  ---------------------------------------------------------

  : Symmetry of relative pairing amplitudes $\bm{\chi}^{\Gamma,\alpha}_{i}$ for a pairing gap function of the trilayer 2$H$-TaS$_2$, classified by the out-of-plane mirror symmetry. 

\
Odd &

  ----------------------------------------------------------------------
   $\psi^{\Gamma,\alpha}_{i}$ or $\bm{d}_{i}^{\Gamma,\alpha} || \bm{z}$
                $\bm{d}_{i}^{\Gamma,\alpha} \perp \bm{z}$
  ----------------------------------------------------------------------

  : Symmetry of relative pairing amplitudes $\bm{\chi}^{\Gamma,\alpha}_{i}$ for a pairing gap function of the trilayer 2$H$-TaS$_2$, classified by the out-of-plane mirror symmetry. 

&

  ---------------------------------------------------------
      ($\frac{1}{\sqrt{2}}$, 0, -$\frac{1}{\sqrt{2}}$)
   $\chi^{\Gamma,\alpha}_{1,i}=\chi^{\Gamma,\alpha}_{3,i}$
  ---------------------------------------------------------

  : Symmetry of relative pairing amplitudes $\bm{\chi}^{\Gamma,\alpha}_{i}$ for a pairing gap function of the trilayer 2$H$-TaS$_2$, classified by the out-of-plane mirror symmetry. 

\

\[tab2\]

In order to determine the detailed paring symmetry $\Delta(\bm{k})$ and T$_c$ of the superconducting state, we solve the following coupled linearized gap equations, $$\begin{aligned}
\label{eqn15}
    [\Delta_{l}(\bm{k})]_{s_1s_2}&=& \frac{-T_c}{N}\sum\limits_{\omega_{n}\bm{k}^\prime s_3 s_4}V_{s_1s_2s_3s_4}(\bm{k},\bm{k}^\prime) \nonumber\\
    && \times [G(k)\Delta(\bm{k}^\prime)G^{\tau}(-k^\prime)]_{\substack{l,l \\s_3s_4}}
    ,\end{aligned}$$\[eq15\] where $G(k)=[i\omega_n-H_0(\bm{k})]^{-1}$ with $k\equiv(\bm{k},i\omega_{n})$ is the Matsubara Green’s function for the normal state.  By solving the above eigenequation for each IR $\Gamma$, one obtain its eigenstates and corresponding T$_c$. The most favorable pairing state just corresponds to the eigenstate with the highest T$_c$. See Appendix for the detail. 

![ \[fig02\](color online). ( a) $v_0$ versus $v_1$ pairing phase diagram of monolayer $H$-TaS$_2$, where the white(orange) dots show the equal T$_c$ line with T$_c$=10 K(1$\times 10^{-2}$ K).  The FS and nodes of (b) sector $+i$ and (c) sector $-i$ for the $\bm{d}^{A_1,z}+\Psi^{A_1,nn}$ state.  The FS are colored in red(blue) to represent the positive(negative) gap functions.  The black(yellow) dots denote the nodes with winding number +1(-1).  Parameters are chosen to be $(v_{0},v_{1})=(0.1,0.16)$ eV, $(\Delta^s,\Delta^f)=\Delta(-0.252,0.968)$ with $\Delta=1$ K. (d)The corresponding edge band for sector $+i$ with open-boundary conditions along y direction, where for better view, $\Delta$ has been enlarged 20 times without changing the topological nature of the system. ]( Fig2.eps){width="7.5cm"}

Result and discussion {#sec4}
=====================

Monolayer $H$-TaS$_2$ and trilayer 2$H$-TaS$_2$
-----------------------------------------------

As a comparison, we first study the superconducting state for the monolayer $H$-TaS$_2$ without Rashba SOC.  This monolayer system preserves the out-of-plane mirror symmetry $i\sigma_{z}$ which has a requirement on the directors $\bm{d}^\Gamma$ of the triplet components: $\bm{d}^\Gamma$ is either parallel or normal to the TaS$_2$ plane, as shown in Table \[tab1\].  When $\bm{d}^\Gamma\parallel\bm{z}$($\bm{d}^\Gamma\perp\bm{z}$), the mirror parity of the superconducting state is even(odd).  Thus in this noncentrosymmetric monolayer system, the singlet components can be mixing merely with the triplet components having $\bm{d}^\Gamma\parallel\bm{z}$ to form an even-mirror-parity state, while the odd-mirror-parity state is a pure triplet state only consisting of components with $\bm{d}^\Gamma\perp\bm{z}$. The odd-mirror-parity state can be ruled out as Ising pairing is dominant due to strong Ising SOC.  Fig.  \[fig02\](a) is the pairing phase diagram calculated for the monolayer $H$-TaS$_2$. One finds two even-mirror-parity states, both of which are mixing ones: the mixing state between $\bm{d}^{A_1,z}$ and $\Psi^{A_1,nn}$, and that between $\Psi^{A_1,on}$ and $\Psi^{A_1,nn}$. The latter is a topologically trivial s-wave state, which is fully gapped on all the FS, while the former could be topologically nontrivial, since $\bm{d}^{A_1,z}$ is a $S_z=0$ f-wave pairing, which holds three nodal lines. 

![( color online).\[fig03\](a)The same as Fig.  \[fig02\], except that it is for trilayer 2$H$-TaS$_2$. Symbol ‘e’ represents even mirror parity. ( b)The FS and nodes of of sector $+i$ for the representative mixing state denoted by ‘P’ in (a), shown respectively for the two subsectors.  Here the gap value is chosen to be $\Delta \approx 1$ K. (c) The corresponding edge band for sector $+i$ .](Fig3.eps){width="7.5cm"}

[|p[4cm]{}&lt;|p[4cm]{}&lt;|p[4cm]{}&lt;|p[4cm]{}&lt;|c|]{} $b^{A_1,on}$; $\bm{\chi}^{A_1,on}$ & $b^{A_1,nn}$; $\bm{\chi}^{A_1,nn}$& $b^{A_1,z}$; $\bm{\chi}^{A_1,z}$ &$b^{A_1,xy}$; $\bm{\chi}^{A_1,xy}$&State\
& 0.179; (-0.642, 0.419,-0.642)&0.983; (0.586, 0.560, 0.586)&&$P$\
0.990; (0.572, 0.584, 0.572)&-0.044; (0.570, 0.580, 0.570) &&&$Q$\
& 0.151; (-0.627, 0.455,-0.627)&0.892; (0.583, 0.566, 0.583)&-0.424(0.707, 0,-0.707)&$R$\
0.967; (0.595, 0.539, 0.595)&0.100; (-0.570,-0.580,-0.570) & &-0.225(0.707, 0,-0.707)&$T$\

\[tab3\]

To investigate the topological feature of the $\bm{d}^{A_1,z}+\Psi^{A_1,nn}$ state, which can be favorable when the NN paring is dominant, one can block-diagonalize the monolayer Hamiltonian into $\pm$i(mirror parity) sectors according to the mirror symmetry, $ H_{BdG}(\bm{k}) \rightarrow \begin{bmatrix}
    H_{+i}(\bm{k}) &0 \\
    0 &H_{-i}(\bm{k}) \\
    \end{bmatrix}$. The two Hamiltonians are connected with each other by TRS: $H_{-i}(\bm{k})=H_{+i}(-\bm{k})$. The excitation spectrum is: $E_{+i}(\bm{k})=\sqrt{(\epsilon(\bm{k})-\beta(\bm{k}))^2+(\Delta^f(\bm{k})+\Delta^s(\bm{k}))^2}$, $E_{-i}(\bm{k})=E_{+i}(-\bm{k})$, with $\Delta^f(\bm{k})=\Delta^f S(\bm{k})$ and $\Delta^s(\bm{k})=\Delta^s C(\bm{k})$. Here $(\Delta^s,\Delta^f)=\Delta(b^{A_1,nn},b^{A_1,z})$, which are found to take real values.  There exist 12 nodal points and all nodes are located at the Fermi pockets around $\Gamma$, while the superconducting state is fully gapped with its order parameter taking opposite signs at the Fermi pockets around $K$ and $K'$(see Fig.  \[fig02\](b)-(c)).  Because the superconducting state is chiral symmetric, each node is characterized by a winding number(WN).  It can be calculated via, $$\begin{aligned}
\label{eqn17}
    \nu_{\pm i}=\frac{1}{2 \pi i}\oint dk Tr\partial_{k} lnq_{\pm i}(\bm{k})
    ,\end{aligned}$$ where the integration is along any small loop around the node and $q_{\pm i}(\bm{k})=i[\Delta^{f}(\bm{k})\pm\Delta^{s}(\bm{k})]+[\epsilon(\bm{k})\mp\beta(\bm{k})]$ appears in the off-diagonal representation of $H_{\pm i}(\bm{\bm{k}})$: $H_{\pm i}(\bm{k})\rightarrow
    \begin{bmatrix}
    0 & q_{\pm i}(\bm{k}) \\
    q^\dagger_{\pm i}(\bm{k}) & 0
    \end{bmatrix}$. The edge states calculated for $H_{+i}$ with zigzag boundary shown in Fig.  \[fig02\](d) give the Majorana flat edge band lying between the projections of pairs of nodes having the opposite WNs. 

Now we turn to study the superconducting state in trilayer 2$H$-TaS$_2$. Its pairing phase diagram is illustrated in Fig.  \[fig03\](a), which is similar to that of the monolayer $H$-TaS$_2$: There are two mixing Ising paring states, including a trivial full-gap s-wave state and a topologically nontrivial mixing one between $\bm{d}^{A_{1},z}$ and $\Psi^{A_{1},nn}$. The $\bm{d}^{A_{1},z}+\Psi^{A_{1},nn}$ state is dominated by the f-wave triplet component $\bm{d}^{A_{1},z}$, and the expansion coefficients for a representative state in the phase diagram is given in Table \[tab3\].  As before, we block-diagonalize the mirror-symmetric $H_{BdG}(\bm{k})$ into two TRS connected sectors $H_{\pm i}(\bm{k})$. Due to $S_{z}$ conservation, $H_{+i}(\bm{k})$ can be further reduced by a basis change to two smaller sectors ${\expandafter\@slowromancap\romannumeral 1@}$ and ${\expandafter\@slowromancap\romannumeral 2@}$: $H_{+i}(\bm{k})\rightarrow\begin{bmatrix}
    H_{+i}^{{\expandafter\@slowromancap\romannumeral 1@}}(\bm{k})& 0 \\
    0 &H_{+i}^{{\expandafter\@slowromancap\romannumeral 2@}}(\bm{k})\\
    \end{bmatrix}$ with

$$\begin{aligned}
        H_{+i}^{{\expandafter\@slowromancap\romannumeral 1@}}(\bm{k})=
        \begin{bmatrix}
        \epsilon(\bm{k})+\beta(\bm{k}) & \Delta_1^f(\bm{k})-\Delta_1^s(\bm{k}) \\
        \Delta_1^f(\bm{k})-\Delta_1^s(\bm{k}) & -\epsilon(\bm{k})-\beta(\bm{k}) \\
        \end{bmatrix}
        ,\end{aligned}$$

$$\begin{aligned}
        H_{+i}^{{\expandafter\@slowromancap\romannumeral 2@}}(\bm{k})=
        \begin{bmatrix}
        \epsilon(\bm{k})-\beta(\bm{k})  & \sqrt{2} t_\perp& 0& \Delta_2^f(\bm{k})+\Delta_2^s(\bm{k}) \\
        \sqrt{2} t_\perp&\epsilon(\bm{k})+\beta(\bm{k}) &  \Delta_1^f(\bm{k})+\Delta_1^s(\bm{k}) & 0 \\
        0 & \Delta_1^f(\bm{k})+\Delta_1^s(\bm{k})&  -\epsilon(\bm{k})+\beta(\bm{k})& -\sqrt{2} t_\perp \\
        \Delta_2^f(\bm{k})+\Delta_2^s(\bm{k}) &0 & -\sqrt{2} t_\perp& -\epsilon(\bm{k})-\beta(\bm{k}) \\
        \end{bmatrix},
        \end{aligned}$$

Here $H_{\pm i}^{{\expandafter\@slowromancap\romannumeral 1@}({\expandafter\@slowromancap\romannumeral 2@})}(\bm{k})$ is written in the basis of ($c_{-,\bm{k}\downarrow},c^{\dagger}_{-,\bm{-k}\uparrow}$) and ($c_{+,\bm{k}\uparrow},c_{2,\bm{k}\uparrow},c^{\dagger}_{+,-\bm{k}\downarrow},c^{\dagger}_{2,-\bm{k}\downarrow}$), respectively, with $c_{\pm,\bm{k}\sigma}=\frac{1}{\sqrt{2}}(c_{1,\bm{k}\sigma}\pm c_{3,\bm{k}\sigma})$. For a definite band $n$, the superconducting state on its FS can be described by the effective gap function, or the projected gap, which is given by the diagonal entry $\widetilde{\Delta}_{n,n}(\bm{k})$, where the effective gap function $\widetilde{\Delta} (\bm{k})$ is given by $\widetilde{\Delta} (\bm{k})=U^\dagger(\bm{k})\Delta (\bm{k})U^*(-\bm{k})$. Here the unitary matrix $U(\bm{k})$ diagonalizes $H_{0}(\bm{k})$: $U^{\dagger}(\bm{k})H_{0}(\bm{k})U(\bm{k})=D(\bm{k})$, with $D(\bm{k})$ a diagonal matrix.  At the FS around $K$ or $K'$, the superconducting state is fully gapped and the projected gap takes positive and negative values alternatively, as shown in Fig.  \[fig03\](b).  The same reason leads to the existence of nodes at the FS around $\Gamma$. Around each FS, the sign of the projected gap changes 6 times and there are 6(12) nodes for $H^{{\expandafter\@slowromancap\romannumeral 1@}}_{+i}(\bm{k})$($H^{{\expandafter\@slowromancap\romannumeral 2@}}_{+i}(\bm{k}))$ and 36 nodes in total for $H_{BdG}(\bm{k})$. Since both $H^{{\expandafter\@slowromancap\romannumeral 1@}}_{+i}(\bm{k})$ and $H^{{\expandafter\@slowromancap\romannumeral 2@}}_{+i}(\bm{k})$ hold chiral symmetry, each node has a well-defined WN $\pm1$ which can still be calculated via Eq.(\[eqn17\]), as exhibited in Fig.  \[fig03\](b). 

![ \[fig04\](color online).(a) The pairing phase diagram with strong Rashba SOC $\alpha_R=50$ meV. The inset is the corresponding band structure.  The FS around $\Gamma$ of (b)sector $+i$ and (c)sector $-i$ for the representative new mixing state denoted by ‘R’ in (a).  The arrows represent the projected gap functions $\widetilde{\Delta}_{n,n}(\bm{k})$ on the innermost FS, with their lengths(angles) denoting the absolute values(phases) of $\widetilde{\Delta}_{n,n}(\bm{k})$. The black(yellow) dots represent the nodes with WN +1(-1). ]( Fig4.eps){width="7.5cm"}

Taking the empirical temperature-dependent upper critical magnetic field of thin layer 2$H$-TaS$_2$, which can not be explained by a pure singlet or triplet pairing, into account[@Barrera2018], the nodal f+s-wave pairing state mentioned above can be a promising candidate.  This mixing state has been proposed in the gated superconducting MoS$_2$[@Yuan2014; @Hsu2017].  It is also consistent with the STM experiment on 2$H$-TaS$_2$, where a ZBCP in the superconducting TaS$_2$ detached flakes was observed[@Galvis2014]. 

In the above discussion, the Rashba SOC $\alpha_R$ has been neglected for simplicity.  If a relatively small $\alpha_R$ is taken into account, the above main results are qualitatively unchanged.  However, if $\alpha_R$ is assumed to be sufficiently large, a $\bm{d}^{A_{1},xy}$ pairing component is expected to be induced, since the Rashba SOC favors the triplet pairing with $\bm{d}^\Gamma \perp \bm{z}$. A detailed calculation confirms this point and gives the pairing phase diagram as shown in Fig.  \[fig04\](a) for $\alpha_R$ up to $50$ meV, which is so large that it is competing with the intrinsic SOC and hence strongly suppressing Ising pairing.  Here $S_z$ conversation is violated, so neither $H_{+i}(\bm{k})$ nor $H_{-i}(\bm{k})$ could be block-diagonalized as before.  This system with strong Rashba SOC is found to be still gapless but with reduced number of nodes: All the FS sheets are fully gapped except the two innermost ones around $\Gamma$, each of which has 6 nodes, as depicted in Fig.  \[fig04\](b)-(c), respectively. 

![ image](Fig5.eps){width="15cm"}

Doping and pressure effects of trilayer 2$H$-TaS$_2$
----------------------------------------------------

The superconducting behavior can be significantly tuned by doping.  Experimentally, the bulk TMDs can be chemically doped with Na and Cu or electron doped by the substrate[@Fang2005; @Wagner2008; @Albertini2017; @Lian2017], while the thin films on a substrate can also be effectively doped by a gate voltage[@Lu2015; @Xi2016; @Lu2018].  Here in the trilayer 2$H$-TaS$_2$, we assume a rigid band and consider the effect of p-type doping by fixing the chemical potential $\mu=-100$ meV, which is near the Van Hove singularity at $M$. Fig.  \[fig05\](a) gives the paring phase diagram, where a new even-mirror-parity mixing phase $\Psi^{E,nn}+ \bm{d}^{E,z}$ of the IR $E$ appears.  Because the IR $E$ is 2D, in this mixing state any combination between $\Psi^{E,nn}_1$ and $\Psi^{E,nn}_2$ (or $\bm{d}^{E,z}_1$ and $\bm{d}^{E,z}_2$) is allowed and shares identical T$_{c}$ determined by Eq.(\[eqn15\]).  In order to determine which combination is the most energetically favorable, one can make an energy minimization, which leads to a TRS breaking mixing state between d+id(d-id)- and p-ip(p+ip)-wave pairings.  The detailed expansion coefficients for a representative state for this phase is given in Table \[tab4\].  This mixing Ising pairing phase is fully gapped.  As an example, Fig.  \[fig05\](c) show the projected gap functions on the FS.  A mirror Chern number(MCN ) $C_{\pm i}$ can be defined for each sector $H_{\pm i}(\bm{k})$. More conveniently, MCN $C^{{\expandafter\@slowromancap\romannumeral 1@}({\expandafter\@slowromancap\romannumeral 2@})}_{\pm i}$ is also meaningful for each subsector $H_{\pm i}^{{\expandafter\@slowromancap\romannumeral 1@}({\expandafter\@slowromancap\romannumeral 2@})}(\bm{k})$, as total spin $S_{z}$ is still conserved in the superconducting state.  In the weak-coupling limit, each MCN can actually be viewed as the sum over the phase WNs of the projected order parameter $\widetilde{\Delta}_{n,n}(\bm{k})$ on each FS, as depicted in Fig.  \[fig05\](c).  Here we have a new WN defined on the n-th Fermi pocket $\gamma$:

$$\begin{aligned}
        \nu_n=\frac{1}{2 \pi i}\oint_\gamma dk \partial_{k} ln\widetilde{\Delta}_{n,n}(\bm{k})
        ,\end{aligned}$$

with the circuit integration along the FS.  Therefore, the total Chern number of the trilayer 2$H$-TaS$_2$ is -6, namely $(C_{+i},C_{-i})=(-3,-3)$, which is consistent with the corresponding chiral edge states of $H_{+i_{}}$ shown in Fig.  \[fig05\](e). 

[|p[4cm]{}&lt;|p[4cm]{}&lt;|p[4cm]{}&lt;|p[4cm]{}&lt;|c|]{} $b^{E,nn}_{1}$; $\bm{\chi}^{E,nn}_{1}$ & $b^{E,nn}_{2};\bm{\chi}^{E,nn}_{2}$ & $b^{E,z}_{1}$; $\bm{\chi}^{E,z}_{1}$ & $b^{E,z}_{2};\bm{\chi}^{E,z}_{2}$&State\
0.285; (-0.637, 0.434,-0.637)&0;$\backslash$&0.958; (0.571, 0.589, 0.571)&0;$\backslash$&$B$\
0.182; (-0.665, 0.341,-0.665)& 0;$\backslash$&0.983; (0.564, 0.602, 0.564)&0;$\backslash$&$C$\

\[tab4\]

On the other hand, since the couplings between TaS$_2$ layers are weak van der Waals forces, an uniaxial pressure can also be applied to tune the features of the material.  The pressure dependences of T$_c$ in 2$H$-TMDs have been measured recently[@Suderow2005; @Tissen2013; @Freitas2016; @Lian2017; @Grasset2018].  Here we assume that the only effect of the uniaxial pressure along z axis is the enhancement of interlayer coupling $t_\perp$. We set $t_\perp=-90$ meV here and get the pairing diagram as shown in Fig.  \[fig05\](b).  There are three mixing phases.  Except the trivial one, both the gapless $\Psi^{A_{1},nn}+\bm{d}^{A_1,z}$ and fully gapped $\Psi^{E,nn}_{1}+\bm{d}^{E,z}_{1}$ phases appear.  Similar to the doping case, the $\Psi^{E,nn}_1+ \bm{d}^{E,z}_1$ state under pressure is gapful and TRS breaking, but takes different Chern number, indicating it is topologically different from that in Fig.  \[fig05\](a).  In detail, the MCN of the two subsectors here are $C^{{\expandafter\@slowromancap\romannumeral 1@}}_{+i}=-3$ and $C^{{\expandafter\@slowromancap\romannumeral 2@}}_{+i}=-6$(see Fig.  \[fig05\](d)).  The total Chern number is thus -18, as $C_{+i}=C_{-i}=-9$, consistent with the corresponding currents-carrying chiral edge states of $H_{+i}(\bm{k})$ shown in Fig.  \[fig05\](f).  This mixing pairing phase has a rather high T$_{c}$(about 200K) and is expected to be possibly realized by high-pressure experiments on 2H-TaS$_2$ thin flakes. 

Conclusion {#sec5}
==========

In summary, based on the linearized gap equation and symmetry analysis, we have obtained the pairing phase diagram of the superconducting monolayer and trilayer TaS$_2$, and found a nodal f+s-wave state.  We suggest that this nodal pairing could be responsible for the anomalous tunneling conductance observed in STM experiments.  The nodal structure is so robust that even a strong mirror-symmetric Rashba SOC(up to 50 meV) cannot fully gap the system.  Besides, both p-type doping and uniaxial pressure along z axis could induce a TRS breaking mixing state between d+id- and p-ip-wave pairings, which has a large Chern number.  Our result indicates that the superconducting trilayer 2$H$-TaS$_2$ could be a promising candidate for realization of topological superconductors.  This study will be helpful to understand the unconventional superconductivity in the thin layer 2$H$-TaS$_2$ and other 2D TMDs.  However, final determination of the pairing symmetry of this kind of 2D Ising superconductors requires more theoretical and experimental efforts. 

Acknowledgments
===============

We thank Jinpeng Xiao, Feng Xiong, Qinli Zhu and Yao Zhou for useful discussions.  This work is supported by NSFC Project NO.111774126 and 973 Projects No.2015CB921202.\

Appendix {#appendix .unnumbered}
========

In this appendix we first prove the orthogonal relations Eq.(\[eqn6\]) and then show in detail how to solve the eiqenequation (\[eqn15\]) for a definite IR $\Gamma$, obtain its T$_c$ and pairing gap function $\Delta(\bm{k}) $ for a multilayer superconductor. 

To prove Eq.(\[eqn6\]), we define $f_{ij}(\Gamma,\alpha,\Gamma^\prime,\beta)$ as follows,

$$\begin{aligned}
        f_{ij}(\Gamma,\alpha,\Gamma^\prime, \beta)&\equiv&\sum_{\bm{k}}Tr\{\Delta^{\Gamma,\alpha}_{i}(\bm{k})\Delta^{\dagger \Gamma^\prime,\beta}_{j}(\bm{k})\} \nonumber\\
        & =&\sum_{\bm{k}}Tr\{D^\dagger(R)\Delta^{\Gamma,\alpha}_{i}(\bm{k})D^*(R)D^\tau(R)\Delta^{\dagger \Gamma^\prime,\beta}_{j}(\bm{k})D(R)\}\nonumber\\
        &=&\frac{1}{g}\sum_{R,m,n,\bm{k}}D^{\Gamma}_{im}(R)D^{*\Gamma^\prime}_{jn}(R)Tr\{\Delta^{\Gamma,\alpha}_{m}(R^{-1}\bm{k})\Delta^{\dagger \Gamma^\prime,\beta}_{n}(R^{-1}\bm{k})\} \nonumber\\
        & =&\sum_{m,n}\frac{1}{d^\Gamma}\delta_{ij}\delta_{mn}\delta_{\Gamma\Gamma^\prime}f_{m,n}(\Gamma,\alpha,\Gamma^\prime, \beta)\nonumber\\
        & =&\delta_{ij}\delta_{\Gamma\Gamma^\prime}\frac{1}{d^\Gamma} \sum_mf_{m,m}(\Gamma,\alpha,\Gamma, \beta),
        \end{aligned}$$

where $g$ is the order of the symmetry group $G$ of the system, and $R$ a group element with $D(R)$ its $2\times2$ spin-rotation representation.  The orthogonal relation between different IRs has been used here in the derivation.  Remarkably, $f_{ii}(\Gamma,\alpha,\Gamma,\beta)$ is independent of $i$. Moreover, for $\alpha\neq\beta$ in the same $\Gamma$, $f_{ii}(\Gamma,\alpha,\Gamma,\beta)$ is generically zero, so one has: $$\begin{aligned}
    f_{ij}(\Gamma,\alpha,\Gamma^\prime, \beta)=\delta_{\Gamma\Gamma^\prime}\delta_{\alpha\beta}\delta_{ij}f(\Gamma,\alpha),
    \end{aligned}$$ where the positive number $f(\Gamma,\alpha)\equiv f_{ii}(\Gamma,\alpha,\Gamma,\alpha)$ can be renormalized to be $1$ if $\Delta^{\Gamma,\alpha}_{i}(\bm{k})$ has been properly normalized.  Thus we come to the orthogonal relations of Eq.(\[eqn6\]).  These orthogonal relations can also be rewritten as, $$\begin{aligned}
    \frac{2}{N}\sum_{\bm{k}}\Psi^{\Gamma,\alpha}_{i}\Psi^{*\Gamma^\prime,\beta}_{j}=\delta_{\Gamma\Gamma^\prime}\delta_{\alpha\beta}\delta_{ij}\\
    \frac{2}{N}\sum_{\bm{k}}\bm{d}^{\Gamma,\alpha}_{i}\cdot \bm{d}^{*\Gamma^\prime,\beta}_{j}=\delta_{\Gamma\Gamma^\prime}\delta_{\alpha\beta}\delta_{ij}
    \end{aligned}$$ which can be easily confirmed by checking Table \[tab1\] for crystal C$_{3v}$.

Now by making use of the above orthogonal relations we try to solve Eq.(\[eqn15\]) for a multi-layer system.  Multiplying the two sides of the equation by $[\Delta^{* \Gamma,\alpha}_{i}(\bm{k})]_{s_1,s_2}$, then taking trace over spin indices and making a sum over $\bm{k}$, one has:

$$\begin{aligned}
        \frac{1}{N}\sum\limits_{\bm{k}}Tr\{\Delta^{\dagger \Gamma,\alpha}_{i}(\bm{k})\Delta_l(\bm{k})\}&=& \frac{-T_c}{N^{2}} \sum\limits_{\omega_n,\bm{k}s_1s_2} [\Delta^{*\Gamma,\alpha}_{i}(\bm{k})]_{s_1s_2}\sum\limits_{\bm{k}^\prime s_3s_4}V_{s_1s_2s_3s_4}(\bm{k},\bm{k}^\prime)[G(k)\Delta(\bm{k}^\prime)G^{\tau}(-k^\prime)]_{\substack{l,l\\s_3s_4}} \nonumber\\
        &=&\frac{-T_c}{N^2} \sum\limits_{\omega_n,\bm{k}s_1s_2}[\Delta^{*\Gamma,\alpha}_{i}(\bm{k})]_{s_1s_2}\sum \limits_{\substack{\bm{k}^\prime s_3s_4 \\\Gamma^\prime,\beta,j}} v^{\Gamma^\prime,\beta} [\Delta^{\Gamma^\prime,\beta}_{j}(\bm{k})]_{s_1s_2}[\Delta^{*\Gamma^\prime,\beta}_{j}(\bm{k}^\prime)]_{s_3s_4} \nonumber\\
        &&\times [G(k)\Delta(\bm{k}^\prime)G^{\tau}(-k^\prime)]_{\substack{l,l\\s_3s_4}} \nonumber\\
        &=&\frac{-T_c}{N^2} \sum \limits_{\omega_n,\bm{k}^\prime,\Gamma^\prime,\beta} v^{\Gamma^\prime,\beta} Tr\{\Delta^{\dagger\Gamma^\prime,\beta}_{i}(\bm{k}^\prime)[G(k)\Delta(\bm{k}^\prime)G^{\tau}(-k^\prime)]_{\substack{l,l}}\}\sum\limits_{\bm{k}}Tr\{\Delta^{\dagger\Gamma,\alpha}_{i}(\bm{k})\Delta^{\Gamma^\prime,\beta}_{i}(\bm{k})\}  \nonumber\\
        &=&\frac{-T_c}{N}  \sum\limits_{\omega_n,\bm{k}^\prime }v^{\Gamma,\alpha}Tr\{ [G(k)\Delta(\bm{k}^\prime)G^{\tau}(-k^\prime)]_{l,l}\Delta^{\dagger \Gamma,\alpha}_{i}(\bm{k}^\prime)\} \nonumber\\
        &=&\frac{-T_c}{N} \sum\limits_{\omega_n,\bm{k} }v^{\Gamma,\alpha}Tr\{ [G(k)\Delta(\bm{k})G^{\tau}(-k)]_{l,l}\Delta^{\dagger \Gamma,\alpha}_{i}(\bm{k})\},
        \end{aligned}$$

where N is the total number of unit cells of the system.  Substitute the expansion of $\Delta(\bm{k})$ Eq.(\[eqn9\]) into the above equation, one has a simplified version of the eigenequation, $$\begin{aligned}
    a^{\Gamma,\alpha}_{l,i}=\sum_{l^{\prime}\beta}{Q^{\Gamma,\alpha \beta}_{l,l^{\prime}}a^{\Gamma,\beta}_{l^{\prime},i}}
    ,\end{aligned}$$ where $a^{\Gamma,\alpha}_{l,i}=b^{\Gamma,\alpha} \chi^{\Gamma,\alpha}_{l,i}$. $Q^{\Gamma,\alpha \beta}_{l, l^{\prime}}$ is independent of i and reads,

$$\begin{aligned}
        Q^{\Gamma,\alpha \beta}_{l, l^{\prime}}&=&\frac{-T_c v^{\Gamma,\alpha}}{N} \sum\limits_{\beta,\bm{k},\omega_n} Tr\{G_{l l^{\prime}}(k) \Delta^{\Gamma,\beta}_{i}(\bm{k}) G_{l l^{\prime}}^{\tau}(-k)\Delta^{\dagger \Gamma,\alpha}_{i}(\bm{k})\} \nonumber\\
        &=&\frac{-T_c v^{\Gamma,\alpha}}{N} \sum_{\bm{k}\sigma \sigma ^{\prime}} \sum_{l_1 l_2} [\widetilde{\Delta}_{l^{\prime},i}^{\Gamma,\beta}(\bm{k})]_{ \substack{l_1 l_2 \\ \sigma \sigma ^{\prime}}}[\widetilde{\Delta}_{l,i}^{\Gamma,\alpha }(\bm{k})]^*_{ \substack{l_1l_2 \\ \sigma \sigma ^{\prime}}}\frac{\tanh[\beta_{c}\xi_{l_1,\sigma} (\bm{k})/2]+\tanh[\beta_{c}\xi_{l_2,\sigma ^{\prime}} (-\bm{k})/2]}{2 [\xi_{l_1,\sigma} (\bm{k})+\xi_{l_2,\sigma ^{\prime}} (-\bm{k})]},
        \end{aligned}$$

where $l$, $l^{\prime}, l_1$ and $l_2$ are the layer indices, with $l,l^\prime,l_1,l_2=1,2,\cdots,n$, $\sigma,\sigma^\prime$ the spin indices, $G_{l l^{\prime}}(k)$ is the $2\times2$ Matsubara Green’s function for the normal state, which takes the form, $$\begin{aligned}
    G_{l l^{\prime}}(k)=\sum_{l_1}U_{l l_1}(\bm{k})[i \omega _n -D_{l_1}(\bm{k})]^{-1}[U_{l^{\prime}l_1}(\bm{k})]^{\dagger} \nonumber \\
    .\end{aligned}$$ Here $U_{ll'}(\bm{k})$ is the $2\times2$ block of the unitary matrix $U(\bm{k})$, which diagonalize $H_0(\bm{k})$: $U^{\dagger}(\bm{k})H_0(\bm{k})U(\bm{k})=D(\bm{k})$. The $2n\times2n$ diagonal matrix $D(\bm{k})$ has the eigenenergies $\xi_{j}(\bm{k})(j=1,2,...,2n)$ of $H_0(\bm{k})$ as its diagonal entries.  The $2\times2$ matrix $D_{l}(\bm{k})$ is given by $D_{l}(\bm{k})=\left(
    \begin{array}{cc}
    \xi_{2l-1}(\bm{k}) & 0 \\
    0 & \xi_{2l}(\bm{k})\\
    \end{array}
    \right)$ with $\xi_{l,\uparrow}(\bm{k})\equiv\xi_{2l-1}(\bm{k})$, $\xi_{l,\downarrow}(\bm{k})\equiv\xi_{2l}(\bm{k})$. The effective gap matrix $[\widetilde{\Delta}_{l^{\prime},i}^{\Gamma,\alpha}(\bm{k})]_{l_1 l_2}$ is defined as, $$\begin{aligned}
    [\widetilde{\Delta}_{l^{\prime},i}^{\Gamma,\alpha}(\bm{k})]_{l_1 l_2}=U_{l^{\prime}l_1}^{\dagger}(\bm{k}) \Delta^{\Gamma,\alpha}_{i}(\bm{k}) U_{l^{\prime}l_2}^{*}(-\bm{k})
    .\end{aligned}$$

For a definite IR $\Gamma$, if the number of different $\alpha$ is $m$, then $Q$ is an $nmd^{\Gamma}\times nmd^{\Gamma}$ matrix.  Since $Q$ is independent of index $i$, it is actually a block-diagonalized one, with each block identical to each other.  Therefore we only have to solve the eigenequation corresponding to the reduced $nm\times nm$ $Q$ matrix, and then determine its T$_c$ and eigenstate(the expansion coefficients of the gap function). 

Now we demonstrate how to use Mirror symmetry to further reduce the $Q$ matrix.  We take n=2N$_0$(N$_0$ is an integer) as an example.  Only $a^{\Gamma,\alpha}_{1},\cdots ,a^{\Gamma,\alpha}_{N_0}$ are independent because the Mirror symmetry ensures that $$\begin{aligned}
\label{a10}
    a^{\Gamma,\alpha}_{n+1-l}=\eta_M \eta^{\Gamma,\alpha}_{}a^{\Gamma,\alpha}_{l}
    ,\end{aligned}$$ where $\eta_M=\pm 1$ denotes the mirror parity, while $\eta^{\Gamma,\alpha}_{}=+1(-1)$ for $\psi_i^{\Gamma,\alpha}$ or $\bm{d}_i^{\Gamma,\alpha}||\bm{z}$($\bm{d}_i^{\Gamma,\alpha}\perp \bm{z}$).  Thus we have: $$\begin{aligned}
 \label{a32}
    a^{\Gamma,\alpha}_{l} &=& \sum\limits_{l^{\prime}=1}^{N_0} \sum\limits_{\beta}^{m}{(Q^{\Gamma ,\alpha \beta}_{ l,l^\prime}+\eta_M \eta^{\Gamma,\alpha}_{}Q^{\Gamma,\alpha \beta}_{ l, n+1-l^\prime})a^{\Gamma,\beta}_{l^{\prime}}} \nonumber\\
    &=&\sum_{l^{\prime}=1}^{N_0} \sum_{\beta}^{m}{\widetilde{Q}^{\Gamma ,\alpha \beta}_{ l,l^\prime}a^{\Gamma,\beta}_{l^{\prime}}}\\
    \widetilde{Q}^{\Gamma,\alpha \beta}_{l,l^{\prime}}
    &=&Q^{\Gamma,\alpha \beta}_{ l,l^\prime}+\eta_M \eta^{\Gamma,\beta}Q^{\Gamma,\alpha  \beta}_{l,n+1-l^{\prime}}
    .\end{aligned}$$

Thus the reduced $\widetilde{Q}$ is an $N_0m \times N_0m$ matrix.