\section{introduction}
Recent progress in our understanding of non-interacting Bloch electrons~\cite{classification1, classification2,  mb, fkm, hm, qi1} reveals a large class of gapped topological phases, the so-called topological insulators and superconductors~\cite{hm, QiRev, KaneRev}. For example, a time-reversal-symmetric topological insulator is a band insulator that cannot be continuously tuned into a trivial atomic insulator, as long as time-reversal symmetry is respected. A topological insulator is featured by a single Dirac cone in its surface state spectrum. Typically, these topological phases are realized in systems with strong spin-orbit coupling. It is known that when topological insulators are combined with superconductivity via the proximity effect~\cite{fu1} or via phonon-mediated attractive interaction~\cite{fu2}, the interesting interplay between electron pairing and the spin-orbit coupling results in exotic superconductivity. For instance, when a topological insulator is doped and turned into a superconductor, an odd-parity topological superconductor is obtained~\cite{fu2} with a single gapless Majorana surface state, which is protected by time-reversal symmetry.

Another type of gapless `topological matter', the Weyl semimetal, is currently being studied intensively and is proposed to be realized in experiments~\cite{Ashvin, Burkov, Yang, cho, latticeWeyl}. Its electronic structure has an {\it even} number of Weyl nodes -- two cylyndrical 3D cones that touch at their apex, the Weyl point --  which carry non-trivial winding numbers ensuring their stability. These Weyl nodes can be thought of as 3D analogs of the two component Dirac fermions in graphene and at the surface of a 3D topological insulators.  They exhibit spin-momentum locking and thus require strong spin-orbit coupling to be realized. In analogy with superconductivity in topological insulators, it is natural to expect interesting superconducting states to emerge in these systems upon doping, resulting from their non-trivial topological winding numbers. To realize a Weyl semimetal phase requires either time-reversal symmetry~\cite{Ashvin} or inversion symmetry~\cite{murakamiweyl} to be broken.  In this paper we concentrate on the inversion-symmetric case, in which the two nodes connected by the inversion symmetry carry opposite chirality. Upon slight doping there are at least two disconnected components to the Fermi surface around the nodes, shifted in momentum space from the inversion-symmetric high-symmetry points such as the $\Gamma$-point. We will show that the interplay between the finite-momentum displacement and the non-trivial winding numbers around each Weyl node leads to interesting superconducting states. 

The finite momentum shift of the Fermi surface motivates the study of finite-momentum pairing states or Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states~\cite{FF,LO}. FFLO states break translational symmetry and have interesting physical properties~\cite{FF,LO, disorder3}. In the Weyl semimetals, the center of momentum of the FFLO pairs is fixed by the momentum of the Weyl nodes. Similarly, the non-trivial winding around the nodes and the broken time-reversal symmetry suggests the possibility of realizing even/odd-parity BCS states that are electronic analogues of the $^{3}$He-A phase~\cite{He3a1,He3a2}. Since the $^{3}$He-A phase has nodes with non-trivial winding number which guarantees the existence of a dispersionless surface states~\cite{He3a2}, the Weyl semimetal in these phases is also expected to support zero-energy surface flat bands, similar to a Weyl semimetal in proximity to a superconductor~\cite{He3a2,WeylSc,Fa,nodal}.

Surprisingly, when the attractive interaction is completely local in real space and represents a phonon-mediated interaction, we find from a self-consistent mean-field calculation that the {\it fully-gapped} finite-momentum pairing is energetically favored over the even-parity BCS state (both pairing states can be thought of as spin-`singlet' pairings though `singlet' is not a very exact terminology since spin-rotational symmetry is broken) and is stable against weak disorder. Hence, there is a good chance of experimentally observing these exotic phases. To be concrete, we concentrate on a specific lattice model realizing an inversion-symmetric Weyl semimetal and solve the gap equation of the lattice model in the BCS approximation. We also discuss the applicability of our result to other models realizing Weyl semimetals.

The proximity effect of an s-wave superconductor on an undoped Weyl semimetal has been studied by Meng and Balents in Ref.~\onlinecite{WeylSc}. In contrast to this work, where the superconductivity is {\it extrinsic}, we are interested in the {\it intrinsic} superconductivity of the doped Weyl semimetal. Note that if the Weyl semimetal is undoped, the intrinsic superconducting gap and the critical temperature are expected to be vanishingly small since the density of states goes to zero at the Weyl point.

\section{Model}
The model we consider in this work is given by the Hamiltonian 
\begin{equation}
	H = H_0 + V_{\rm ee}
	\label{eq:Model}
\end{equation}
where $V_{\rm ee}$ is an electron-electron interaction term to be specified below. For the kinetic term, $H_0$, we take the minimal two-band lattice model~\cite{Yang} 
\begin{align}
H_0 = t& (\sigma^{x} \sin k_{x} + \sigma^{y} \sin k_{y})+t_{z}(\cos k_{z} - \cos Q)\sigma^{z}  \nonumber \\  &+m(2-\cos k_{x}-\cos k_{y})\sigma^{z} -\mu.
\label{lattice}
\end{align}
This model realizes a Weyl semimetal with two Weyl points at momenta ${\vec P}_\pm = (0,0,\pm Q)$. $\sigma^{x,y,z}$ are the Pauli sigma matrices (for later use we define $\sigma^0$ to be the $2\times2$ unit matrix), $t$ and $t_z \sin Q$ are the Fermi velocities at the Weyl points in the $x,y$ and $z$ directions respectively. Without a loss of generality, we assume $t = t_{z}\sin(Q)$ such that the Fermi velocity around the Weyl points is isotropic. We have explicitly included the chemical potential $\mu$ in the kinetic term. We are primarily interested in the parameter range $0 < |\mu/t| \ll Q$, when the Fermi surface consists of two disconnected spherical components around the Weyl points (see Fig.~\ref{Fig1}). In this case, the states on the Fermi surface have spin-momentum locking similar to the surface states of a strong topological insulator. This property will play an important role in our discussion of pairing states below.

The electron-electron interaction is short ranged and takes the form
\begin{equation}
V_{\rm ee} = V_{0} \sum_{i} n_{i}n_{i} + V_{1} \sum_{\langle ij\rangle}n_{i}n_{j} =  \sum_{{\vec k}} V({\vec k})n_{{\vec k}}n_{-{\vec k}},
\label{interaction}
\end{equation}
where $n_{i} = \sum_{\sigma} c^{\dagger}_{i, \sigma}c_{i,\sigma}$ is the number of electrons on site $i$, and the second sum is over nearest neighbors only. $V({\vec k}) = V_{0}+V_{1}(\cos k_{x} + \cos k_{y} + \cos k_{z})$ is the Fourier transform of the real-space interaction. $V_{0}$ represents a phonon-mediated attractive on-site interaction and $V_{1}$ the nearest-neighbor interaction. We are mainly interested in the case $V_{0}<0$ and $|V_0|\gg |V_{1}|$ when electrons form Cooper pairs and condense. The phenomenological interaction term~\eqref{interaction} captures the tendency of d-wave pairing for $V_{1}<0$ and $V_{0}>0$ in the context of high-Tc superconductors such as cuprates~\cite{MacDonald}. 

The point group symmetry of the model Hamiltonian, $H$, is $C_{4h} =\{I^{\eta_I}C_4^{\eta_4}|\eta_I=0,1;\eta_4=0,1,2,3\}$, where
\begin{align}
I&: \sigma^{z} H(-{\vec k})\sigma^{z} = H({\vec k}), \nonumber \\
C_{4}&: S^{\dagger} H[R_{\pi/2}({\vec k)}]S = H({\vec k}),\label{symmetry}
\end{align}
with $S = \frac{1}{\sqrt{2}}(\sigma^{0}+i\sigma^{z})$ and $R_{\pi/2}$ a rotation by an angle $\pi/2$ around the $z$-axis. $I$ is the inversion symmetry that takes ${\vec r} \rightarrow -{\vec r}$. Each spatial rotation is accompanied by an equal spin rotation, manifesting the spin-momentum locking due to the presence of strong spin-orbit coupling. The $C_4$-rotation symmetry in the $xy$-plane is not required to realize a Weyl semimetal, but is present in this model and similar lattice rotations are present in other models we discuss later. 

To make a connection to previous work~\cite{Ashvin, Yang, cho, Burkov} and to obtain a general understanding of Weyl semimetal phases, we derive the low energy effective theory corresponding to the lattice model~\eqref{eq:Model}. Expanding $H(\vec k)$ in the small momentum ${\vec q} = {\vec k} - {\vec P}_\pm$ around the two Weyl Points denoted by $\pm$, we obtain
\begin{equation}
H_0 = \sum_{{\vec k}}c^{\dagger}({\vec k}) H_0({\vec k}) c({\vec k}) \approx \sum_{a= \pm} \psi^{\dagger}_{a}({\vec q}) h_{a}({\vec q}) \psi_{a} ({\vec q}).
\label{weyl}
\end{equation}
The effective kinetic Hamiltonian $h_{\pm}({\vec q})$ is given by
\begin{equation}
h_{\pm}({\vec q}) = t (q_{x}\sigma^{x}+q_{y}\sigma^{y} \mp q_{z}\sigma^{z}) - \mu.
\label{low}
\end{equation}
Similarly, we obtain for the interaction term
\begin{equation}
V_{\rm ee} = \sum_{{\vec k},{\vec p}, {\vec q}}V^{ab;cd}({\vec q})\psi^{\dagger}_{a, \sigma}({\vec k}+{\vec q})\psi^{\dagger}_{b,\tau}({\vec p}-{\vec q})\psi_{c,\tau}({\vec p})\psi_{d,\sigma} ({\vec k}),
\end{equation}
where roman letters denote the nodal indices $\pm$ and $\sigma,\tau$ are spin indices. Here and henceforth, repeated indices are summed over. In the BCS channel (see App.~\ref{app:interactions} for details) 
\begin{equation}
V_{\rm ee} =\sum_{{\vec k},{\vec l}}V^{ab;cd}\psi^{\dagger}_{a, \sigma}({\vec k})\psi^{\dagger}_{b,\tau}(-{\vec k})\psi_{c,\tau}(-{\vec l})\psi_{d,\sigma} ({\vec l}),
	\label{eq:VeeBCS}
\end{equation}
with
\begin{align}
&V^{-+;+-}=V^{+-;-+}= V_{0} + 3V_{1} - \frac{V_{1}}{2}({\vec k}-{\vec l})^{2},\nonumber \\
&V^{-+;-+}= V_{0}+2V_{1} +V_{\perp} + V^{+}_{\parallel},\label{ContinuumInteraction}\\
&V^{+-;+-} =V_{0}+2V_{1} +V_{\perp} + V^{-}_{\parallel}, \nonumber 
\end{align}
and 
\begin{align}
V_{\perp} &= -\frac{V_{1}}{2} ({\vec k}_{\perp}-{\vec l}_{\perp})^{2},\\
V^{+}_{\parallel}/V_{1} &= [1 - \frac{1}{2}(k_{z}-l_{z})^{2}]\cos2Q + (k_{z}-l_{z})\sin2Q,\nonumber \\
V^{-}_{\parallel}/V_{1} &= [1 - \frac{1}{2}(k_{z}-l_{z})^{2}]\cos2Q - (k_{z}-l_{z})\sin2Q,\nonumber
\end{align}
where ${\vec k}_{\perp} = (k_{x},k_{y},0)$. These expressions will be used in the next section.

\section{Mean field theory and pairing channels}
\label{MFEnergy}
We treat the interaction term $V_{\rm ee}$ in a mean-field approximation and solve the resulting gap equations self-consistently. In addition to the more standard BCS paring, we also study finite-momentum or FFLO pairing.

\subsection{BCS pairing}
\begin{figure}
\includegraphics[width=1\columnwidth]{SC_Weyl3}
\caption{A schematic diagram of the spin texture around the Weyl nodes in momentum space and the pairing states. (a) The spin direction of the eigenstates is given by thick arrows. The double-headed arrows labeled with (1) and (2) indicate the partner states in the BCS pairing. The spin state is maximally anti-parallel for (1) and parallel for (2), indicating that there will be nodes in the latter case if the pairing is in the singlet channel. Contrary to the BCS pairing, the FFLO pairing (3) connects two states within the same node (`intra-nodal' pairing). The two states connected by the FFLO pairing have the opposite spin directions. (b) Position of the nodes for the even-parity state. The nodes of the same chirality are on the same component of the Fermi surface, with their partner nodes of opposite chirality on the other. The filled circle represents a node of chirality $+1$ and the crossed circle represents a node of chirality $-1$. Hence, there are four nodal points on the Fermi surface of the even-parity paired BCS state.}
\label{Fig1}
\end{figure}

The symmetry classification of different BCS pairing order parameters in a doped Weyl semimetal, according to the lattice symmetry~\eqref{symmetry}, is summarized in Table \ref{table1} (see App.~\ref{app:symm} for more details). There are three fully-gapped BCS pairing order parameters ($\Gamma^1$ and $\Gamma^{3,\pm}$) and one that has nodal lines ($\Gamma^2$). 
\begin{table}
\begin{tabular} {l|l|l|l}
\hline
{}&IRR&$C_{4h}$&pairing function\\ \hline
$\Gamma^{1}$&$A$&$1$& $Z^2(i\sigma^{y}),(X-iY)(\sigma^0+\sigma^z),(X+iY)(\sigma^0-\sigma^z)$\\ \hline
$\Gamma^{2}$&$B$&$-1$&$XY(i\sigma^y),(X+iY)(\sigma^0+\sigma^z),(X-iY)(\sigma^0-\sigma^z)$\\ \hline
$\Gamma^{3,+}$&$E_{+}$&$i$&$(X+iY)Z(i\sigma^y),Z(\sigma^0+\sigma^z),XYZ(\sigma^0-\sigma^z)$\\ \hline
$\Gamma^{3,-}$&$E_{-}$&$-i$&$(X-iY)Z(i\sigma^y),XYZ(\sigma^0+\sigma^z),Z(\sigma^0-\sigma^z)$\\
\hline
\end{tabular}
\caption{Symmetry classification of the BCS pairing order parameters for the model (\ref{lattice}) according to the different irreducible representations (IRRs) of the group $C_{4h}$. $X$, $Y$, and $Z$ are basis functions for the momentum-space pairing function which we take to be $\sin p_{x}$, $\sin p_{y}$, and $\sin p_{z}$ respectively, and can be realized by nearest-neighbor pairings. Among the pairing functions, $i\sigma^y$ means singlet pairing, $\sigma^x$ is the spinfull triplet pairing, $\sigma^0+\sigma^z$ the triplet pairing for polarized $\uparrow$ spins and $\sigma^0-\sigma^z$ the triplet pairing for polarized $\downarrow$ spins. The paired state $\Gamma^2$ has nodal lines, while the other three states are gapless.}
\label{table1}
\end{table}

In the continuum theory, the pairing terms of Table~\ref{table1} take the form
\begin{equation}
\sum_{{\vec k}}\Delta_{\sigma\tau}({\vec k})c^{\dagger}_{\sigma}({\vec k})c^{\dagger}_{\tau}(-{\vec k})\approx \sum_{{\vec q}}\Delta^{a,b}_{\sigma\tau} ({\vec q}) \psi^{\dagger}_{a,\sigma}({\vec q})\psi^{\dagger}_{b,\tau} (-{\vec q}),
\end{equation}
The standard BCS pairing term connects two Weyl nodes in the effective theory. The explicit form of $\Delta^{a,b}_{\sigma\tau}$ and $\Delta_{\sigma\tau}$ can be found in Table \ref{table1} and in Eq.~\eqref{Pair} in App.~\ref{app:interactions}. The self-consistent gap equation takes the form
\begin{equation}
\Delta^{ab}_{\sigma\tau}({\vec p}) = \sum_{{\vec k}}V^{ab;cd}({\vec p} - {\vec k}) \langle\psi_{c,\tau}(-{\vec k})\psi_{d,\sigma}({\vec k})\rangle,
\end{equation}
where the expectation value is taken with respect to the mean-field superconducting state (see App.~\ref{app:interactions} for more explicit expressions for the gap equations).

\subsection{FFLO pairing}
In the doped Weyl semimetal, the Fermi surface is formed around the Weyl points ${\vec P}_\pm$, and it is natural to expect a finite-momentum pairing to compete with the standard BCS-paired states. We therefore introduce a FFLO state with a center of momentum at $2{\vec P}_\pm$, which paring function satisfies 
\begin{equation}
\Delta^{\pm}_{\rm FFLO}({\vec r}) \propto \exp(2i{\vec P_+}\cdot {\vec r}) \pm \exp(2i{\vec P}_-\cdot {\vec r})
\label{fflo}
\end{equation}
The self-consistent equations for these pairing order parameters take the form
\begin{equation}
\Delta_{\sigma\tau}({\vec p}; \pm {\vec P}) = \sum_{{\vec k}}V({\vec p} -{\vec k}) \langle \psi_{\pm,\tau}(-{\vec k})\psi_{\pm,\sigma}({\vec k})\rangle,
\end{equation}
where the two nodes $\pm$ are decoupled. These FFLO states correspond to the {\it intra-node} pairing, in contrast to the BCS case which is {\it inter-node} pairing (see Fig.~\ref{Fig1}). The two states of the pairings $\Delta^{\pm}_{\rm FFLO}$ in Eq.~\eqref{fflo} with a relative phase of $\pm1$ between the two components of the Fermi surface have the {\it same} mean-field energy since the two nodes are decoupled in the mean-field theory.
\begin{table}
\begin{tabular} {l|l|l}
\hline
{}&$C_{4}$&Pairing function\\ \hline
$\Gamma^{1}$&$1$& $Z^2(i\sigma^{y}),(X-iY)(\sigma^0+\sigma^z),(X+iY)(\sigma^0-\sigma^z)$\\ \hline
$\Gamma^{2}$&$-1$&$XY(i\sigma^y),(X+iY)(\sigma^0+\sigma^z),(X-iY)(\sigma^0-\sigma^z)$\\ \hline
$\Gamma^{3,+}$&$i$&$(X+iY)Z(i\sigma^y),Z(\sigma^0+\sigma^z),XYZ(\sigma^0-\sigma^z)$\\ \hline
$\Gamma^{3,-}$&$-i$&$(X-iY)Z(i\sigma^y),XYZ(\sigma^0+\sigma^z),Z(\sigma^0-\sigma^z)$\\
\hline
\end{tabular}
\caption{Classification of the FFLO states of superconducting Weyl fermions based on the lattice symmetry $C_{4}$. The notation is the same as in Table~\ref{table1}. We assume that the center of momentum for the pairing is at ${\vec P}_\pm$. Note that the symmetry is only $C_{4}$ on the $xy$-plane without the inversion because `inversion' is already encoded by the ansatz Eq.~\eqref{fflo}. This classification is essentially the same as that of BCS-type pairing order parameters.}
\label{table2}
\end{table}


\subsection{Mean field energy}
Having identified the possible superconducting states, we compute their free energy by solving the self-consistent gap equations numerically. We are interested in the case $|V_0|\gg |V_1|$ with $V_0 < 0$ where the spin singlet is preferred. We denote the pairing term $\propto i\sigma^y$ in Table~\ref{table1} as 'singlet' and the other terms $\propto i{\vec \sigma}\sigma^{y}$ as 'triplet'. The singlet and triplet components have a different dependence on the interaction parameters $V_{0}$ and $V_{1}$. The gap of the triplet components depends only on the value of $V_{1}$, while the singlet component depends only on $V_{0}+3V_{1} \approx V_{0}$ or $V_{0}+ V_{1}(5 + \cos 2Q)/2 \approx V_{0}$. We therefore consider in the following only the singlet component $\propto i\sigma^y$ of $\Gamma^1$ for the BCS and FFLO states. 

For these two states the BCS mean-field approximation is
\begin{equation}
	H = H_0 + V_{\rm ee}^{\rm pair},
	\label{eq:Hmf}
\end{equation}
where $H_0$ is given by Eq.~\eqref{lattice} and $V_{\rm ee}^{\rm pair}$ is the effective {\it projected} pair potential derived from the lattice interaction in Eq.~\eqref{interaction}. For the $\Gamma^{1}$-BCS state we have 
\begin{align}
V_{\rm ee}^{\rm pair} &= -U_{\rm BCS}\sum_{{\vec k}, {\vec p}}P^{\dagger}_{{\vec k}}P_{-{\vec p}}, \notag \\		
U_{\rm BCS} &= V_{0}+V_{1}\frac{5 + \cos(2Q)}{2},\\
P^{\dagger}_{{\vec k}} &= \psi^{\dagger}({\vec k}) \tau^{x}i\sigma^{y} \psi^{*}(-{\vec k}), \notag
\end{align}
and the gap equation is 
\begin{equation}
\Delta = -\frac{U_{\rm BCS}}{4} \int_{{\vec k}} \langle\psi_{a,\alpha}({\vec k}) (\tau^{x})^{ab}(-i\sigma^{y})^{\alpha\beta}\psi_{b,\beta}(-{\vec k})\rangle.
\end{equation}
For the $\Gamma^{1}$-FFLO state we obtain
\begin{align}
V_{\rm ee}^{\rm pair} &= -U_{\rm FFLO}\sum_{{\vec k}, {\vec p}}P^{\dagger}_{{\vec k}}P_{{\vec p}}  \notag \\
U_{\rm FFLO} &= V_{0}+3V_{1}, \\
P^{\dagger}_{{\vec k}} &= \psi^{\dagger}({\vec k}) i\sigma^{y} \psi^{*}(-{\vec k}), \notag 
\end{align}
with a gap equation 
\begin{equation}
\Delta =  -\int_{{\vec k}} \frac{U_{\rm FFLO}}{2} \langle\psi_{a,\alpha}({\vec k})(-i\sigma^{y})^{\alpha\beta}\psi_{a,\beta}(-{\vec k})\rangle
\end{equation}

In this standard BCS-type approximation, we can evaluate the energy $E$ of the pairing states with the pairing amplitude $\Delta_{\vec k} = -U\langle P_{-{\vec k}}\rangle$ (with the effective pairing interaction strength $U$) by computing
\begin{align}
E &= E_{\rm el} + E_{\rm sc}, \nonumber\\
E_{\rm el} &= \sum_{{\vec k} \in d\Omega, \epsilon_{\rm sc} <0} \epsilon_{\rm sc}({\vec k})n_{e}[\epsilon_{\rm sc}({\vec k})] \nonumber \\
& -\sum_{{\vec k} \in d\Omega, \epsilon_{\rm fs}<0} \epsilon_{\rm fs}({\vec k})n_{e}[\epsilon_{\rm fs}({\vec k})], \nonumber\\
E_{\rm sc} &= - \sum_{{\vec k} \in d\Omega}\frac{\Delta_{{\vec k}}\Delta^{*}_{-{\vec k}}}{2U} + h.c.
\label{Energy}
\end{align}
Here $n_{e}$ is the filling of the electron for the state at ${\vec k}$ of energy $\epsilon({\vec k})$, $\epsilon_{\rm sc}$ is the energy of the filled band of the BdG quasiparticle with mean-field gap $\Delta$, and $\epsilon_{\rm fs}$ is the energy of the filled bands of the free Weyl electrons without pairing, i.e.\ the energy of the normal state. Thus, the second line of Eq.~\eqref{Energy} represents the energy  gain of the superconducting state relative to the normal state by opening up a gap near the Fermi surface. The last line of Eq.~\eqref{Energy} represents the contribution from the pairing interaction labeled by the momentum ${\vec k}$. The range of the summation is restricted to a shell $d\Omega$ around the Fermi surface, which width is determined by the strength of the attractive interaction. 

In Fig.~\ref{Fig2} we plot the mean-field energy for the two parings, as obtained from Eq.~\eqref{Energy}, as a function of the interaction strength $V_0$. The $\Gamma^{1}$-FFLO state has a larger gap than the $\Gamma^{1}$-BCS state and is energetically favored. This result can be understood by considering the spin-momentum locking around the Fermi surface. For the even-parity pairing state, the state $|{\vec k}, \alpha\rangle$ ($\alpha$ is the spin state) is paired with the inversion partner state $|-{\vec k}, \sigma^{z}\alpha\rangle$. The pairing amplitude is of the form $\sim \langle c^{\dagger}({\vec k})i\sigma^{y}c^{*}(-{\vec k}) \rangle$ which takes the maximum value if the two states at ${\vec k}$ and $-{\vec k}$ have opposite spins. However, the spins at ${\vec k}$ and $-{\vec k}$ are not anti-parallel and even become parallel at the poles (which is the origin of the nodes, see Fig.~\ref{Fig1}) which tends to reduce the superconducting gap. In contrast, the FFLO state connects the states $|{\vec k}+{\vec Q}, \alpha \rangle$ and $|-{\vec k}+{\vec Q}, \beta \rangle$ via the spin singlet channel with $\beta = -\alpha$ (anti-parallel spins). A gap opens up everywhere at the Fermi surface with a larger gap than the even-parity BCS state. This is very similar to the surface of topological insulator which we discuess Appendix B (See also Ref. [\onlinecite{FFLOsurface}]). The similar finite-momentum pairing (``intra-valley" pairing or ``Kekule" pairing) can happen in a graphene in the presence of a nearest-neighbor attractive interaction~\cite{Gr1, Gr2}.      

\section{Discussion}

In this section, we discuss the nodal structure of the $\Gamma^{1}$-BCS state and the effect of the disorder on the $\Gamma^{1}$-FFLO state. Like in the last section, we consider only the singlet components of these states, assuming $|V_{0}|\gg |V_{1}|$. 

\subsection{$\Gamma^{1}$-BCS state}
\label{BCS state}

The pair potential term in the mean field Hamiltonian of the singlet component of the $\Gamma^1$-BCS state is
\begin{align}
H^{\rm BCS}_{\rm pair} &= \sum_{{\vec k}}\Delta c^{\dagger}_{\alpha}({\vec k}) (i\sigma^{y})^{\alpha\beta} c^{\dagger}_{\beta}(-{\vec k}) + h.c. \\
&= \sum_{\vec q}\Delta \psi^{\dagger}_{a,\alpha}({\vec q})(\tau^{x})^{ab}(i\sigma^{y})^{\alpha\beta} \psi^{\dagger}_{b,\beta}(-{\vec q})+h.c. \notag
\label{singlet1}
\end{align}
The second form is obtained in the low energy theory. This superconducting state is an even-parity state and has four point nodes on the Fermi surface at $q_{x}=q_{y}=0$ and $q_{z} = \pm \sqrt{\Delta^{2}+\mu^{2}}$ (see Fig.~\ref{Fig1}). The gap remains closed even when the triplet pairings of the $\Gamma^{1}$-BCS state in Table~\ref{table1} are included.

\begin{figure}
\includegraphics[width=1\columnwidth]{SC_Energy5}
\caption{Mean-field energy $E$ of the even-parity $\Gamma^{1}$-BCS and the $\Gamma^{1}$- FFLO states as a function of interaction strength $V_0$. Other model parameters used to obtain this plot are $\mu/t = 0.3$, $Q = 0.7$, $V_1 = 0$, and $d\Omega = 0.2$}
\label{Fig2}
\end{figure}

The two nodal points near Weyl node ${\vec P}_+$ (${\vec P}_-$) carry a winding number of +1 (-1). To demonstrate this we write down the Bogoliubov-de Gennes (BdG) Hamiltonian $H = \sum_{{\vec k}} \Phi_{{\vec k}} {\tilde H}_{{\vec k}} \Phi_{{\vec k}}$ for $\Phi_{{\vec k}} = (c_{{\vec k}}, i\sigma^{y}c^{*}_{-{\vec k}})^{T}$. In the continuum limit at ${\vec P}_+$ (similar expressions are obtained for ${\vec P}_-$) 
\begin{equation}
{\tilde H} =
\begin{pmatrix}
h_{+}({\vec q})& \Delta \sigma^{0} \\
\Delta\sigma^{0}& -h_{-}({\vec q})
\end{pmatrix},
\label{nodes}
\end{equation}
with $h_{\pm}({\vec q})$ defined in Eq.~\eqref{low}. The quasiparticle spectrum corresponding to this BdG Hamiltonian is 
\begin{equation}
E({\vec q}) = \pm [q^{2} + \Delta^{2} + \mu^{2} \pm 2(\Delta^{2}q^{2}_{z} + \mu^{2}q^{2})^{1/2}]^{1/2}
\label{spec}
\end{equation}
which has nodes at $q_x = q_y = 0$, $q_{z} =\pm \sqrt{\Delta^{2}+\mu^{2}}$, both with chirality of $+1$. Near the nodes $|q_x|,|q_y| \ll |q_z|,|\mu|$, we obtain the anisotropic Weyl spectrum
\begin{equation}
E({\vec q}) \approx \pm \left[ (q_{z} \pm \sqrt{\Delta^{2}+\mu^{2}})^{2} + q^{2}_{\perp} (1+ \frac{\mu^{2}}{\mu^{2}+\Delta^{2}})\right]^{1/2},
\end{equation}
with $q_\perp = (q_x,q_y)$. At zero chemical potential, this is similar to the results of Meng and Balents~\cite{WeylSc} who considered the proximity effect of undoped Weyl semimetals. The effect of nonzero chemical potential is to simply shift the Weyl nodes located at $q_{z} = \pm \Delta$ at $\mu =0$, to $q_{z} = \pm \sqrt{\Delta^{2}+\mu^{2}}$. 

Because of the non-trivial winding number carried by the nodes, the nodal points are robust against small perturbations.  The only way to gap out the nodes is to undergo a pair-annihilation of nodes with the opposite winding numbers, and the nodal points are {\it topologically} stable as long as they are separated enough in momentum space. Strikingly, this nodal structure implies that there will be a zero-energy state on the surface which should be detectable in experiment. This is similar to $^{3}$He-A which is an odd-parity pairing state, while our superconducting phase is realized by the even-parity pairing.

\subsection{$\Gamma^{1}$-FFLO state}
\label{FFLO state}

The singlet component of the $\Gamma^{1}$- FFLO state is fully gapped with a mean-field pair potential 
\begin{equation}
H^{\rm FFLO}_{\rm pair} = \Delta c^{\dagger}_{\alpha}({\vec k}+{\vec P}_+)(i\sigma^{y})^{\alpha\beta}c^{\dagger}_{\beta}(-{\vec k}+{\vec P}_+) \pm ({\vec P}_+\leftrightarrow {\vec P}_-)
\label{FFLO}
\end{equation}
with center-of-momentum of $2{\vec P}_{\pm}$. In the low-energy theory, it can be represented by the {\it intra-node} pairing $\sim \Delta \sum_{{\vec q}}\psi^{\dagger}_{a,\alpha}({\vec q}) (i\sigma^{y})^{\alpha\beta} \psi^{\dagger}_{a,\beta}(-{\vec q})$. 

It is known that some two-dimensional FFLO states with strong spin-orbit coupling and parallel magnetic field are unstable against weak disorder~\cite{disorder1,disorder3}. In contrast, the FFLO state discussed in this paper is found to be robust against weak disorder. In fact, the structure of the FFLO state Eq.~\eqref{fflo} and Eq.~\eqref{FFLO} is more similar to the even/odd-parity state of the doped topological insulators studied in Ref.~\onlinecite{disorder2} than usual FFLO states in the two spatial dimension. This similarity is manifested if we write down the pairing for the Weyl fermions in the continuum limit in the helicity eigenstates
\begin{equation}
\Delta_{\pm} \propto e^{i\phi} [\langle \psi_{+}({\vec q})\psi_{+}(-{\vec q}) \rangle \pm \langle \psi_{-}({\vec q})\psi_{-}(-{\vec q}) \rangle],
\label{FFLO2}
\end{equation}
which corresponds to Eq.~(5) of Ref.~\onlinecite{disorder2}. Within this Cooper channel, we add a scalar disorder term to the Hamiltonian
\begin{equation}
H_{\rm imp} = V_{\rm imp} \sum_{{\vec k}, {\vec p} \in FS}c^{\dagger}_{{\vec k},\sigma}c_{{\vec p},\sigma} = \sum_{{\vec q}, {\vec l} \in FS} V^{ab}_{{\vec q}, {\vec l}} \psi^{\dagger}_{a,{\vec q}}\psi_{b,{\vec l}}
\end{equation}
The matrix element $V^{ab}_{{\vec q}, {\vec l}}$ is given by
\begin{equation}
V^{ab}_{{\vec q}, {\vec l}} = V_{imp} 
\begin{pmatrix}
\langle{\hat q}|{\hat l}\rangle& \langle{\hat q} | {\bar l}\rangle  \\
\langle {\bar q}|{\hat l}\rangle & \langle{\bar q}|{\bar l}\rangle
\end{pmatrix},
\label{disorders}
\end{equation}
where we have used the standard normalized spin state ${\hat q} \cdot {\vec \sigma}|{\hat q}\rangle = |{\hat q}\rangle$ and ${\bar q} \cdot {\vec \sigma} |{\bar q}\rangle = |{\bar q}\rangle$ with ${\hat q} = {\vec q}/|{\vec q}|$ and ${\bar q} = ({\vec q}_{\perp}, -q_{z})/ |{\vec q}|$. With this impurity scattering, the self-energy can be worked out in the self-consistent Born approximation, and we find that the correction to the self-energy and the Cooperon diagram due to disorder are exactly of the same form as obtained by Michaeli and Fu~\cite{disorder2}. In fact, the only difference between the FFLO state $\Delta_{\pm}$ in Eq.~\eqref{FFLO} and the even/odd-parity paired states of Ref.~\onlinecite{disorder2} is phase factors in the matrix elements of $V^{ab}_{{\vec k}, {\vec l}}$ in Eq.~\eqref{disorders} which does not show up in the corrections to the self-energy, the Cooperon diagram, and the pairing susceptibilities. Thus we conclude that the critical temperature of $\Gamma^{1}$- FFLO states is not affected by the disorder and thus $\Gamma^{1}$- FFLO state is robust. 


\section{conclusion}

In conclusion, we have studied the possible superconducting states of doped inversion-symmetric Weyl semimetals. We considered a concrete lattice model realizing a Weyl semimetal and found that the FFLO state has a lower energy than the even-parity state if the interaction is phonon-mediated, and the phase is argued to be stable against disorder. Though the even-parity state is less favored in energy than the FFLO state, it interestingly provides an electronic analogue of $^{3}$He-A phase. 

We remark briefly on the implication of our work for superconducting states of Weyl semimetal models other than the one studied in this paper. Among the many proposals for the Weyl semimetal phase, we restrict ourselves to the models based on the topological insulators~\cite{cho, Burkov} with the time-reversal breaking perturbation.
\begin{equation}
H = v \tau^{z}{\vec \sigma}\cdot {\vec k}_{\perp} + \tau^{x}k_{z} + m\sigma^{z}
\label{real}
\end{equation}
The typical symmetry of the model is $I \times C_{n}$ ($\times M$, Mirror symmetry) where $I$ is the inversion symmetry and $C_{n}$ is the $n$-fold lattice rotation symmetry along a certain axis (for the model based on Bi$_{2}$Se$_{3}$~\cite{cho}, we have $n=3$). Due to the strong spin-orbit interaction, the spatial symmetry operation involves the spin/orbital operations, e.g., $I: {\vec k}\rightarrow -{\vec k}$ should involve $\tau^{y}$, and $\tau^{y}H(-{\vec k}) \tau^{y} = H({\vec k})$ (the lattice rotation will involve a spin rotation). Note that these symmetry considerations already manifest the similarity between the realistic model and the simplified model Eq.~\eqref{lattice}, and this similarity becomes much clearer if we go to the low-energy theory of Eq.~\eqref{real}. It is not difficult to confirm that the low-energy theory is identical to Eq.~\eqref{low}, and hence we will have similar superconducting states, FFLO and electronic analogues of $^{3}$He-A, in the more realistic model. Hence, we predict that the superconducting states we found should show up in other proposals for Weyl semimetals.

Note that FFLO state Eq.~\eqref{FFLO} shows a density modulation pinned by the momentum of the Weyl nodes (which is reminiscent of the field-induced charge density wave~\cite{Yang} of the Weyl semimetals). Many experimentally available Weyl semimetals have a large number of Weyl nodes, 
for example the irridates which have 24 nodes~\cite{Ashvin} (or an inversion-symmetry broken Weyl semimetal has {\it at least} four Weyl nodes~\cite{murakamiweyl}). While our minimal model calculation here does not guarantee that the FFLO state will be the lowest energy state in such systems, 
at minimum it suggests that it will be a competing state. In this case the FFLO state can have multiple centers of momenta. This directly implies that there will be interesting density modulation patterns which are fully determined by the position of the Weyl nodes. (This is true at least at the level of mean-field theory which ignores the effect of $O(\Delta^{4})$ terms in the Landau-Ginzburg theory. $O(\Delta^{4})$ terms can potentially {\it melt} this pattern). 



We also note that FFLO state can host interesting half-quantum vorticies discussed in Ref. [\onlinecite{VortexFFLO}]. In the FFLO state, we have {\it two} indepedent superconducting order parameters $\Delta(\pm {\vec P}) \propto \exp(\pm 2i{\vec P}\cdot {\vec r})$,i.e., the order parameter space is $S^{1}\times S^{1}$. The half-quantum vortex corresponds to a unit ``winding" of the phase of the $\Delta({\vec P})$ while the phase of the $\Delta(-{\vec P})$ does not wind. On the other hand, the Fermi surface around the Weyl node at $\vec P$ encloses the $\pi$-Berry phase~\cite{Hosur_Vishwanath} which signals that there will be a gapless ``chiral'' Majorana mode at the core of the half-quantum vortex. Furthermore, this implies that a full quantum vortex will be a composite of the two half-quantum vortices and each half-quantum vortex will have a chiral mode. Thus the full quantum vortex will host a helical Majorana mode. In contrast to the related case~\cite{Hosur_Vishwanath}, this helical Majorana mode is not symmetry protected and is therefore generally gapped out. Furthermore, the helical Majorana mode can be understood as the critical point between a weak pairing state and a strong pairing state in a $1$D p-wave superconductor~\cite{Kitaev}. There are two possible phases for the full quantum vortex depending on the sign of ``mass gap'' for the helical mode~\cite{Kitaev} and in a nontrivial phase there will be a Majorana fermion at the end of the vortex.    


\acknowledgements
The authors thank Pavan Hosur, Ashvin Vishwanath, Sid Parameswaran, Eun Gook Moon, Yong Baek Kim, and Tarun Grover for helpful discussion and Daniel Agterberg for useful comments from which we learn about the half-quantum vortex in FFLO states. The authors acknowledge support from NSF DMR-1206515 (G. Y. C. and J. E. M.), Office of BES, Materials Sciences Division of the U.S. DOE under contract No.DE-AC02-05CH1123 (Y. M. L.), and the LBNL Thermoelectrics Program (J. H. B.) of DOE BES.


