You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
abstract: 'We study vortex zero-energy bound states in presence of pairing between low-energy Dirac fermions on the surface of a topological insulator. The pairing symmetries considered include the $s$-wave, $p$-wave, and, in particular, the mixed-parity symmetry, which arises in absence of the inversion symmetry on the surface. The zero-mode is analyzed within the generalized Jackiw-Rossi-Dirac Hamiltonian that contains a momentum-dependent mass-term, and includes the effects of the electromagnetic gauge field and the Zeeman coupling as well. At a finite chemical potential, as long as the spectrum without the vortex is fully gapped, the presence of a single Fermi surface with a definite helicity always leads to one Majorana zero-mode, in which both electron’s spin projections participate. In particular, the critical effects of the Zeeman coupling on the zero-mode are discussed.'
author:
- 'Chi-Ken Lu'
- 'Igor F. Herbut'
title: 'Pairing symmetry and vortex zero-mode for superconducting Dirac fermions'
---

Introduction
============

The zero-energy Majorana modes have recently become a topic of many investigations in condensed matter physics, perhaps due to the concomitant non-Abelian statistics[@ReadGreen] and potential applications to fault-tolerant quantum computation[@nayak]. In a weakly-coupled triplet $p$-wave superconductor, for example, one such zero-mode was found in the core of a half-vortex, and two pairs of such vortices are required to perform the non-Abelian braiding operations[@Ivanov]. The vortex bound states were, of course, already studied a long ago by Caroli, de Gennes, and Matricon in a $s$-wave weak-coupling superconductor,[@deGennes] but a pair of true zero-modes was discovered only much later, by Kopnin and Salomaa in a triplet superfluid[@KopninSalomma]. Both studies were performed within the framework of Bogoliubov-de Gennes (BdG) equations which describe the low-energy electron-hole excitations around the *degenerate* Fermi surfaces in the presence of superconducting instability. The degenerate Fermi surfaces can be split, and the superconducting order parameter can contain both the spin singlet and triplet components when the inversion symmetry is broken[@Gorkov]. For example, the corresponding vortex bound states in a noncentrosymmetric superconductor,[@Lu2] or the Andreev bound states on the domain wall[@Tanaka] have also been studied. The existence of such zero-modes has been shown to depend on a Z$_2$ topological invariant associated the mixed order parameter[@Fujimoto; @ZhangTopology; @YipTopology].

The metallic surface of a topological insulator (TI) resembles the noncentrosymmetric material, as in both cases the Fermi surfaces are chiral, i.e. the projection of the spin is locked to the direction of the momentum. However, an important difference is that the surface of the topological insulator has a only one such Fermi surface [@TIrev]. Recently, an interesting system consisting of this peculiar surface subjected to the superconducting proximity effect has been proposed to bind one Majorana zero-mode inside a $s$-wave vortex core of unit vorticity[@FuKane]. This zero-mode is nothing but the mid-gap state of the Jackiw-Rossi-Dirac Hamiltonian[@Jackiw1] (JRD), studied before in the context of charge fractionalization in polyacetylene[@SSH] and graphene,[@HCM; @Jackiw2] with its existence assured by an index theorem[@IndexW]. In addition, the zero-modes of the JRD Hamiltonian have also been studied recently in the context of the insulating [@HCM; @IgorInsulating; @IgorDirac] and superconducting states on the graphene’s honeycomb lattice[@WilczekSC; @Kopnin; @Babak; @Bergerman; @Igor2010].

In this paper we study the existence and the form of the zero-modes of a class of the general JRD Hamiltonians describing paired Dirac fermions in presence of a superconducting vortex in the spin-triplet $p$-wave state, and in the mixed singlet-triplet state, when the inversion symmetry is absent. We pay special attention to the effects of a finite chemical potential[@Kopnin; @Bergerman] and the magnetic fields on the zero mode, which are all present in real systems. The problem under study is also relevant to graphene, where different symmetries of the superconducting phases have also been considered[@CastroNeto; @Annica; @Roy]. In graphene in particular, other exotic manifestations of the physics of zero-modes become possible due to their rich internal structure, consisting of spin, valley, sublattice, and the Nambu degrees of freedom[@HCM; @IgorInsulating; @Igor2010; @ghaemi]. It is therefore of importance to understand the general conditions for the appearance of the zero-energy states. Quite generally, we find that the zero-modes in the superconducting vortex may exist only for weak enough Zeeman coupling to the magnetic field. This is in accord with our previous study limited to the the vortex in the s-wave superconducting order parameter[@HerbutLu].

This paper has the following organization. We first investigate the spectrum associated with the spatially homogeneous Hamiltonian without vortex for the singlet $s$-wave, triplet $p$-wave and mixed-symmetry cases. The conditions for opening of a gap in all three cases are discussed. We then proceed to the equations for the zero-mode associated with different pairing symmetries, with the Zeeman term and the chemical potential included. In the sec. IV, the effects of the Zeeman field and the relation of the zero-mode with the antilinear operator that provides the reflection symmetry of the spectrum[@HerbutLu] in all symmetry cases are discussed. Finally, a brief summary is given.

Translation-invariant massive Dirac Hamiltonian
===============================================

![A single two-dimensional Fermi surface associated with the Dirac Hamiltonian, $H_0(\vec{k})=v_F(\hat{\rm{z}}\cdot{\vec{k}\times{\vec{\sigma}}})-\mu$, on the metallic surface of a TI. $\hat{\rm{z}}$, the out-of-plane direction, is the spin quantization axis. The dark (dashes) arrows indicate the directions of the spin associated with the momentum state $\vec{k}$ on the Fermi surface of positive (negative) helicity when $\mu$ is positive (negative).[]{data-label="SingleFS"}](SingleFS.eps)

In this section, we consider the Dirac electrons on the metallic surface of a TI and their pairing due to the proximity to a superconducting material. The metallic surface of a TI has a single Fermi surface as shown in Fig. \[SingleFS\], and the spin-dependent spectrum can be described by a simple 2$\times$2 Hamiltonian, $H_0(\vec{k})=v_F\hat{\rm{z}}\cdot\vec{k}\times\vec{\sigma}-\mu$, where $\mu$ is the chemical potential and $v_F$ is the Fermi velocity. The quantization axis for spin is parallel with the normal of the surface $\hat{\rm{z}}$. As usual, the identity matrix and the three Pauli matrices, $\{I_2,\vec{\sigma}\}$, define a basis in the two-dimensional matrix space. One should note that $H_0(\vec{k})$ is invariant under the point group C$_{4{v}}$ in which the inversion symmetry is absent[@YipGarg]. At the mean-field level, the general pairing potential $V_p$ that respects the absence of inversion symmetry and takes different coupling channels (spin singlet and triplet) into account needs to be be expressed as a matrix.[@Gorkov] To this purpose we introduce the general mass matrix $M(\vec{k})=\Delta{I_2}+v_{\Delta}\vec{d}(\vec{k})\cdot\vec{\sigma}$, so that the pairing potential becomes $V_p(\vec{k})=iM(\vec{k})\sigma_y$. Here $\Delta$ stands for the singlet pairing gap while the velocity $v_{\Delta}$ is associated with the triplet pairing.

In the 4-dimensional spin-Nambu space, the BdG Hamiltonian $\hat{H}=\sum_{\vec{k}}\Psi^{\dag}_{\vec{k}}\hat{H}(\vec{k})\Psi_{\vec{k}}$, with $\Psi_{\vec{k}}=\{a_{k,\uparrow},a_{k,\downarrow},a^{\dag}_{-k,\downarrow},-a^{\dag}_{-k,\uparrow}\}^T$, then assumes the form,

()=(

[cccc]{} &H\_0()&M()\
&M\^()&-H\_0()

).\[fullHamiltonian\] For later convenience, we note that $\hat{H}(\vec{k})$ can also be written in terms of the cross-product representation, namely $\hat{H}(\vec{k})=(v_F\vec{n}(\vec{k})\cdot\vec{\sigma}-\mu)\tau_3+M(\vec{k})\tau_{+}
+M^{\dag}(\vec{k})\tau_{-}$, where $\sigma_i$’s act on the spin and $\tau_i$’s on the Nambu indices. The vector $\vec{n}(\vec{k})=(-k_y,k_x)$ lies on the x-y plane, as in Fig. \[SingleFS\]. The chemical potential $\mu$ enters along with the number operator $\rm{N}$, which can be identified as $\tau_3$ in this representation. As we discuss shortly, while they of course never commute with the number operator $\tau_3$, the mass terms may commute with the diagonal, kinetic energy terms. These cases will lead to a peculiar spectrum that is not possible in standard superconductors which have the quadratic energy-momentum dispersion in its normal state.

First, let us discuss different pairing symmetries for the mass term, in the translation-invariant situation without the vortex, and the resulting spectrum $E(\vec{k})$.

The simplest case corresponds to the singlet s-wave pairing in which $\Delta$ is constant and $v_{\Delta}$=0. The spectrum may be obtained most easily by squaring the Hamiltonian and using the anticommuting properties of the Pauli matrices; this gives $[\hat{H}(\vec{k})]^2=(v_F\vec{n}(\vec{k})\cdot\vec{\sigma}-\mu)^2+M^2$, and therefore $E(\vec{k})=\pm\sqrt{(v_Fk\pm\mu)^2+\Delta^2}$.[@FuKane; @Santos] Next we consider the triplet pairing, with the vector $\vec{d}$ coinciding with $\vec{n}$, in which the mass term $M=v_{\Delta}\vec{n}\cdot\vec{\sigma}$ still anticommutes with the kinetic energy. The same trick can be applied again, to obtain the spectrum,

E()=,\[TripletSpec\] which is gapped as long as $\mu$ is finite. At $\mu=0$, the spectrum is gapless, and in fact reduces to the Dirac cone, with an increased velocity. A somewhat unusual situation arises when $\vec{d}=(k_x,k_y)$, for which the mass term [*commutes*]{} with the kinetic energy: when $\vec{n}\bot\vec{d}$, $[M\tau_{\pm},\vec{n}\cdot\vec{\sigma}\tau_3]$=0. It can be shown that the spectrum is then given by $E(\vec{k})=\pm{v_F}k\pm\sqrt{\mu^2+(v_{\Delta}k)^2}$. In this case the spectrum is still gapless, and it is the chemical potential that becomes modified and momentum-dependent, loosely speaking[@Linder].

In the absence of the inversion symmetry the general “mixed" mass term, with both $\Delta$ and $v_{\Delta}$ finite, is also allowed. Let us assume the triplet component to be given by $\vec{d}=\vec{n}$. By squaring the Hamiltonian, one can obtain, $[\hat{H}(\vec{k})]^2=(v_F\vec{n}\cdot\vec{\sigma}-\mu)^2+(\Delta_s+v_{\Delta}\vec{n}\cdot\vec{\sigma})^2$, and consequently

E()=,\[MixedSpec\] where a gapped spectrum exists only for the positive helicity branch when $\Delta=v_{\Delta}k_F=\Delta_{+}/2$, with the Fermi wavevector $k_F=\mu/v_F$.

Equations for zero-energy quasi-particle state
==============================================

In this section we shall consider the inhomogeneity in the mass term $M$ generated by an isolated vortex in which the magnetic field applied along $\hat{\rm{z}}$ is threading through its core. The superconducting phase changes by $2\pi{Q}$ on circling once around the core. Here $Q$ labels the winding number of the vortex. In the context of superconductivity, the mass term here is in general a function of both the center-of-mass position $\vec{r}$ and relative momentum $\vec{k}$ of the Cooper pair, namely, $M=M(\vec{r},\vec{k})$. More explicitly, in the polar coordinate $\vec{r}=(r\cos\phi,r\sin\phi)$, the vortex enters the mass via $M(\vec{r},\vec{k})=e^{{\it{i}}Q\phi}M(\vec{k})$. In the following subsections, we shall employ the BdG equation to investigate the zero-energy bound state in the vortex core for various pairing symmetries. In particular, we also consider the effects of an external Zeeman field $H_{\rm{Z}}=h\sum_{\vec{k}}\Psi^{\dag}_{\vec{k}}(\sigma_3\tau_0)\Psi_{\vec{k}}$, in which $\tau_0$ stands for the identity matrix acting on the Nambu space.

Singlet $s$-wave pairing {#singlet}
------------------------

The mass generated by the vortex of singlet $s$-wave pairing potential with the winding number of unity can be written as $M=e^{i\phi}|\Delta(r)|$, in which the magnitude increases from zero at the core and reaches a positive value of $\Delta_{\infty}$ when the distance $r$ is sufficiently large than the superconducting coherence length. In terms of the Dirac fermion $\Psi$ we defined, the quasi-particle states which diagonalize the BdG Hamiltonian and give exactly zero-energy have the corresponding coefficients ${\bf{u}}=\{u_1(\vec{r}),u_2(\vec{r})\}$ for the electron sector and the corresponding components ${\bf{v}}$ for the hole sector satisfy the following equations,

+e\^[i]{}|(r)|[**[v]{}**]{}=0,\
e\^[-i]{}|(r)|[**[u]{}**]{}+=0, where $\hbar$=$c$=1 is taken for simplicity. The momentum operators, $k_i=-i\partial_{i}$, are coupled to the gauge field $\vec{A}=(-y,x)B(r)/2$, where $B$ varies on a length scale of penetration length $l_p$ inside a superconductor. Using the identity $(\partial_x\pm{i}\partial_y)$=$e^{\pm{i}\phi}(\partial_r\pm (i/r)\partial_{\phi})$, and the angular decomposition, $[{\bf{u}},{\bf{v}}]^T$=$[u_1(r),e^{i\phi}u_2(r),e^{-i\phi}v_1(r),v_2(r)]^T$, one can decouple the set of four first-order differential equations into two blocks, namely,

f\_a +(-h)f\_b=0,\[Forb1\]\
f\_b -(+h)f\_a=0,\[Forb2\]\
g\_a +(-h)g\_b=0,\[decay1\]\
g\_b -(+h)g\_a=0,\[decay2\] where the definitions $u_2\pm{v_1}$=$f_a$($g_a$) and $u_1\mp{v_2}$=$f_b$($g_b$) are made. The explicit dependence of magnetic field $B$ and pairing potential $\Delta$ on the radial coordinate $r$ is not shown to simplify the notation. In fact, it will be seen in the following asymptotic analysis that the effects of coupling to the gauge field on the zero-mode are negligible.

Far from the core the gap function $\Delta(r)\rightarrow\Delta_{\infty}$ and moreover, only the derivative terms are relevant in the kinetic energy. The terms related to $B(r)$ can be dropped since $B(r)$ decay faster than $1/r$ in the case of a isolated vortex[@deGennesB]. Therefore, by defining the new variable $x\equiv\sqrt{|\mu^2-h^2|}r/v_F$, Eqs. (\[Forb1\]) and (\[Forb2\]) are reduced to the following:

(1)=0, where the plus(negative) signs are indicative of that $|\mu|$ is larger(smaller) than $|h|$, respectively. In both cases, as long as $h^2<\mu^2+\Delta^2_{\infty}$, the two solutions for $f_a$ are both exploding at infinity, which, from Eq.(\[Forb1\]), is also true for $f_b$. Thus one concludes that $u_1=v_2$ and $u_2=-v_1$, to ensure the normalizability. The remaining components can be determined then from Eq. (\[decay1\]) and (\[decay2\]), which can be reduced to,

(1)=0, where two decaying solutions are available under the same condition $h^2<\mu^2+\Delta^2_{\infty}$.

Near the core where the argument $x\ll1$, it is legitimate to keep only the singular and the derivative terms. This is because $\Delta(r)$ is vanishingly small there, and the inherent logarithmic divergence from $B(r)$[@deGennesB] is suppressed in the vector potential. Eqs. (\[decay1\]) and (\[decay2\]) can also be written as,

(+-)g\_a(x)=0,\[CoreEq\] where the plus sign still means that $|\mu|$ is larger than $|h|$. Now one can see that $u_2(x)$=$J_1(x)$, the first-order Bessel function, for small argument. Note that the other solution, $N_1(x)$, is forbidden because of its singularity of $1/x$ for small argument. For $|\mu|<|h|$, $u_2(x)$=$I_1(x)$, the modified Bessel function which is regular at small $x$. The other component $u_1(x)$ is not independent and is given by,

u\_1(x)=-(+)u\_2(x) , for small argument $x$. If we think of the pair potential $\Delta(r)$ as a step function which vanishes within some certain range and is a constant outside, the obtained solution from the second-order differential equation for $x\ll1$ can be smoothly joined with the solution for $x\gg1$ at the boundary by matching the values of the function and its first derivative. This can be accomplished when $h^2<\mu^2+\Delta^2_{\infty}$ since, given one function on the left-hand side of the boundary, there are two available decaying solutions on the right-hand side. For an arbitrary $\Delta(r)$ the solutions can be obtained by adiabatic deformation of the step-like pair potential that leaves the asymptotic behaviors intact.

We conclude that, as long as $h^2<\mu^2+\Delta^2_{\infty}$, there exists one zero-energy solution for the $s$-wave pairing vortex of winding number one, with the following radial wave function,

(

[cccc]{} u\_1(r)\
u\_2(r)\
-u\_2(r)\
u\_1(r)

),\[S\_wave\_function\] where the electron components $u_1$ and $u_2$ for small argument are given by,

(

[cccc]{} u\_1(x)\
u\_2(x)

)= (

[cccc]{} J\_0(x)\
-J\_1(x)

),\[s\_wf\_1\] when $h^2<\mu^2$. As for $\mu^2<h^2<\mu^2+\Delta^2_{\infty}$, $J_0(x)$ and $J_1(x)$ are simply replaced with the corresponding modified Bessel functions $I_0(x)$ and $I_1(x)$, respectively. For a large argument,

(

[cccc]{} u\_1(x)\
u\_2(x)

)= (

[cccc]{} (A-B)\
-(A+B)

)e\^[-x]{},\[s\_wf\_2\] which is valid only for $h^2<\mu^2$. For $\mu^2<h^2<\mu^2+\Delta^2_{\infty}$, the trigonometric functions are replaced with the corresponding hyperbolic ones. The coefficients $A$ and $B$ are to match $u_2(x)$ at the boundary and an additional overall coefficient for normalization is not shown explicitly.

A special case of above is when $\mu=h=0$ with the solution given by,

(

[cccc]{} u\_1(r)\
u\_2(r)

)= (

[cccc]{} 1\
0

)e\^[-\^r|(r’)|dr’]{},\[s\_wf\_FuKane\] which is exact for all range of $r$. Comparing with Eqs. (\[s\_wf\_1\]) and (\[s\_wf\_2\]), one should note that the ambiguity arising at the zero chemical potential is removed by the fact that $J_1$ is vanishing for zero argument. Besides, the zero-mode in presence of a finite chemical potential appears with both spin projections in contrast to the zero-mode in Eq. (\[s\_wf\_FuKane\]) where only the up-spin components show up[@FuKane]. For the case with an anti-vortex, $M\sim|\Delta(r)|e^{-i\phi}$, different angular decomposition, $[{\bf{u}},{\bf{v}}]^T$=$[e^{-i\phi}u_1(r),u_2(r),v_1(r),e^{i\phi}v_2(r)]^T$, has to be employed, and, consequently, $u_1\sim{J_1}$ and $u_2\sim{J_0}$ near the core. Moreover, in the special case with vanishing $\mu$ and $h$, the Majorana zero-mode have similar form as in Eq. (\[s\_wf\_FuKane\]) but with down-spin components only.

Triplet $p$-wave pairing {#triplet}
------------------------

Now we consider the triplet order parameter specified by $\vec{d}=(-k_y,k_x)$, which is parallel to the spin direction $\vec{n}$ in the kinetic energy. With this Rashba type $p$-wave pairing, an isotropic superconducting gap is generated for nonzero $\mu$ according to Eq. (\[TripletSpec\]) in contrast to the case with $\vec d\perp\vec n$. The quasi-particle states with zero energy in the presence of vortex in the triplet order parameter with winding number of unity have to satisfy the following equations,

v\_F(\_r+)f\_a+(-h-v\_\_r)f\_b=0\[f1\],\
v\_F\_rf\_b-f\_a=0\[f2\],\
v\_F(\_r+)g\_a+(-h+v\_\_r)g\_b=0\[g1\],\
v\_F\_rg\_b-g\_a=0\[g2\], where identical angular dependence has been employed, and $f_{a(b)}$ and $g_{a(b)}$ have the same definitions as before. $v_{\Delta}(r)$ is also a radial function. When $\mu=h=0$, it can readily be deduced that both $f_a$ and $g_a$ must be identically zero; otherwise the operator $(\partial_r+1/r)$ will lead to the divergence of $\sim 1/r$ in the solution near the origin. Consequently, the zero-energy state must satisfy $\partial_rf_b=\partial_rg_b=0$ for all $r$, which leads to the conclusion that the zero-energy bound state does not exist when $\mu=0$. For general $\mu$ and $h$, on the other hand, we may proceed in the same manner as in previous subsection. Near the vortex core the analysis is identical since the gap velocity $v_{\Delta}$ is vanishingly small there. Far from the core, it is useful to rewrite Eq. (\[f1\]) and (\[f2\]) in terms of $f_a\pm{i}f_b$, as,

(f\_a+if\_b)-ih(f-if\_b)=0,\[faa\]\
(f\_a-if\_b)+ih(f+if\_b)=0,\[fbb\] which also can be recast as a second-order differential equation. For $\mu^2-h^2>0$, the desired equation is $(d^2/dy^2-2\lambda{d}/dy+1)(f_a+if_b)=0$, where the variable $y=\sqrt{|\mu^2-h^2|/(v_F^2+v_{\Delta}^2)}r$ and the constant $\lambda=\mu{v_{\Delta}}/\sqrt{(v_F^2+v_{\Delta}^2)|\mu^2-h^2|}$. The solution then has a plane wave form, $e^{iky}$, with $k=-i\lambda\pm\sqrt{1-\lambda^2}$. Without loss of generality $v_{\Delta}$ may be assumed to be positive. We then first consider $\mu$ to be positive: there are then two exploding solutions for $f_a+if_b$ at infinity when $0<\lambda<1$, or equivalently $\mu^2>(1+v^2_{\Delta}/v^2_F)h^2$. For the case of $\lambda>1$, or $h^2$ is more close to $\mu^2$, the two solutions are still exploding but with $k=-i(\lambda\mp\sqrt{\lambda^2-1})$. On the other hand, for $\mu^2<h^2$, the differential equation becomes $(d^2/dy^2-2\lambda{d}/dy-1)(f_a+if_b)=0$, with the asymptotic solution of wave vector $k=-i(\lambda\pm\sqrt{\lambda^2+1})$. In this case one solution is always decaying, while the other is exploding. Thus we conclude that $f_a=f_b=0$ as long as $\mu^2>h^2$, and we take $u_1=v_2$ and $u_2=-v_1$. The remaining task then is to obtain the $u_i$’s from the equations for $g_a$ and $g_b$, which are identical to Eqs. (\[faa\]) and (\[fbb\]), except that $v_{\Delta}\rightarrow-v_{\Delta}$ and hence two decaying solutions at infinity are available as long as $\mu^2>h^2$.

For $\mu$ negative, the above conclusion has to be reversed, so that $g_a=g_b=0$ when $\mu^2>h^2$, and therefore $u_1=-v_2$ and $u_2=v_1$. The rest of analysis is the same as with positive $\mu$. The complete wave function for zero-energy state has the radial part summarized as following:

(

[cccc]{} u\_1(r)\
u\_2(r)\
-[[sgn]{}]{}()u\_2(r)\
[[sgn]{}]{}()u\_1(r)

),\[p\_wf\_1\] and the asymptotic behaviors are given by

(

[cccc]{} u\_1(y)\
u\_2(y)

)= (

[cccc]{} (A-B)\
(A+B)

)e\^[-||y]{},\[p\_wf\_2\] which has an oscillation on top of a decaying form when $0<\lambda<1$. It will be monotonically decaying at infinity when $\lambda>1$, and the trigonometric functions need to be replaced with the corresponding hypergeometric ones. The coefficients are determined by matching the function $u_2+iu_1$ to the solution near the core.

Mixed $s$-wave and $p$-wave pairing {#mixed}
-----------------------------------

Finally, we consider the combination of the problems in the previous two subsections. In absence of the inversion symmetry the kinetic energy in the normal state can have its chiral form with the pairing potential as the mixture in two order parameters of different parities. We will thus assume it to be a mixture of the singlet and triplet components, with the latter specified by $\vec{d}=(-k_y,k_x)$. For simplicity, we also neglect the Zeeman effect. Far from the vortex core, the corresponding equations for $f$’s are, similarly as before,

(v\_F\_r-\_)f\_a+(||-v\_\_r)f\_b=0,\
(v\_F\_r-\_)f\_b-(||-v\_\_r)f\_a=0, in which $\alpha$ specifies the sign of chemical potential, and we take $\alpha=1$ first. The above equations are identical to those for $g_a$ and $g_b$, except that $\Delta_{\infty}\rightarrow-\Delta_{\infty}$ and $v_{\Delta}\rightarrow-v_{\Delta}$. It follows that at infinity both $f_a\pm{if_b}$ are in the form of $e^{(\kappa\pm{iq}){r}}$ with $\kappa=(\alpha|\mu|{v_{\Delta}+\Delta_{\infty}v_F})/(v^2_{\Delta}+v^2_F)$ positive, and $q=(\alpha|\mu|{v_F-\Delta_{\infty}v_{\Delta}})/(v^2_{\Delta}+v^2_F)$. It is useful to define the parameters

\_=\_\_[k]{}\_F,\[DefDeltaPM\] where $v_Fk_F=|\mu|$ is the Fermi momentum for nonzero chemical potential and the renormalized Fermi velocity $\bar{v}_F=v_F(1+v^2_{\Delta}/v^2_F)$. The exploding solutions can then be written as $\sim e^{\Delta_{+}r/{\bar v}_F}$, which is similar to the case of $s$-wave pairing. Of course, the exploding solutions are forbidden, and therefore $u_1=v_2$ and $u_2=-v_1$. It can be seen that the set of equations for $g_a$ and $g_b$ will lead to two asymptotic decaying solutions in the form $e^{-\Delta_{+}r/{\bar v}_F}e^{\pm{iqr}}$, which in turn determines the zero-mode wavefunctions to be,

(

[cccc]{} u\_1(r)\
u\_2(r)

)= (

[cccc]{} Ae\^[i(qr+)]{}-Be\^[-i(qr-)]{}\
Ae\^[iqr]{}+Be\^[-iqr]{}

)e\^[-\_[+]{}[r]{}/|[v]{}\_F]{},\[mixed\_wf\] with the coefficients determined by the similar boundary conditions for matching $u_2$ to its piece near the core.

For a negative chemical potential, by similar reasoning, the set of equations associated with $f$’s gives two asymptotic solutions as $e^{\Delta_{-}r/{\bar v}_F}e^{\pm{iqr}}$, which will be forbidden if $\Delta_{-}>0$. Then, the zero-mode is determined from the set of equations associated with the $g$’s, which give similar solution as in Eq. (\[mixed\_wf\]), except that $\Delta_{+}\rightarrow\Delta_{-}$. However, the zero-mode will have $u_1=-v_2$ and $u_2=v_1$ in contrast to the previous case when $\Delta_{-}<0$, in which the solutions from the $f$’s were decaying and allowed, while those from $g$’s were forbidden. The wavefunction for $\mu<0$ therefore has a general form,

(

[cccc]{} u\_1(r)\
u\_2(r)\
-[[sgn]{}]{}(\_[-]{})u\_2(r)\
[[sgn]{}]{}(\_[-]{})u\_1(r)

),\[mix\_wf\_1\] where one can see that a vanishing $\Delta_{-}$ will remove the zero-mode. The corresponding $u_1$ and $u_2$ are determined in a similar manner.

Discussions and Conclusions
===========================

In the preceding sections, the existence of zero-energy bound states in the presence of various symmetries of the mass-gap is investigated in the context of superconducting Dirac fermions on the surface of a TI. It may also be useful to write the BdG Hamiltonian in presence of the vortex in terms of the 4-dimensional Dirac $\Gamma$-matrices. The $s$-wave vortex, for example, is described by the Hamiltonian $\hat{H}=\Gamma_1p_1+\Gamma_2p_2+|\Delta(r)|(\Gamma_3\cos\phi+\Gamma_5\sin\phi)+ih\Gamma_1\Gamma_2-i\mu\Gamma_3\Gamma_5$. The last two terms are the Zeeman field and the chemical potential, respectively. Choosing the representation for the matrices as $\Gamma_1=\sigma_z\otimes\sigma_y$, $\Gamma_2=-\sigma_z\otimes\sigma_x$, $\Gamma_3=\sigma_x\otimes{I_2}$, and $\Gamma_5=-\sigma_y\otimes{I_2}$ will then yield the original differential equations. The four matrices, together with $\Gamma_0
= \sigma_z\otimes\sigma_z$, are all Hermitian, and anticommuting among themselves, $\{\Gamma_i,\Gamma_j\}=2\delta_{ij}$ for $i,j=0,1,2,3,5$. It is easy to see that in the absence of the last two terms the matrix $\Gamma_0$ anticommutes with the Hamiltonian, and thus ensures that the spectrum has the reflection symmetry around zero. In this limit there is a total of $|n|$ zero-modes for the vortex of vorticity $n$[@Jackiw1; @IndexW]. It can be shown from the general properties of the Clifford algebra satisfied by the Dirac matrices that there also exists one antilinear operator that anticommutes with $\hat{H}$ even in the presence of the last two, $\sim \Gamma_1\Gamma_2$ and $\sim
\Gamma_3\Gamma_5$ terms[@footnote]. In the representation defined above, in which $\Gamma_1$ and $\Gamma_5$ are the only two matrices among five which are imaginary, the antilinear operator in question is $A=i\Gamma_1\Gamma_3K$, where $K$ denotes the operation of complex conjugation. Explicitly,

A=(

[cccc]{} 0&i\_y\
-i\_y&0

)K. We recognize that the zero-mode $\Psi_0$ in Eq.(\[S\_wave\_function\]) is its eigenvector with eigenvalue $1$: $A\Psi_0=\Psi_0$. Here, it is worthwhile noting that $\Gamma_0$ and $A$ serve as the sublattice (chiral) and particle-hole symmetry operators, respectively. In the presence of finite chemical potential and/or Zeeman field, only the particle-hole symmetry is present and the JRD Hamiltonian belongs to the class D in the symmetry classifications.[@Altland] In the case of vanishing chemical potential and zero magnetic field, the chiral symmetry is restored and the JRD Hamiltonian belongs to the class DIII.

![(Color online) The phases of zero-mode for $h^2$ versus $\mu^2$ in the (a) $s$-wave pairing and (b) $p$-wave pairing cases. The left-handed sides with respect to the empty circles denote the existence of the zero-mode, while the right-handed sides, including the empty circle, indicate that the zero-mode is not normalizable and hence does not exist. The blue regions stand for the fact that the wavefunction has, in addition to the decaying, an asymptotic oscillation part, which is gone when entering the green regions.[]{data-label="PurePhase"}](PhaseDecay.eps)

The effects of the Zeeman field on the zero-mode may be summarized as in Fig. \[PurePhase\](a). For $|h|>\sqrt{\mu^2+\Delta^2_{\infty}}$, specified by the empty circle, the zero-mode is not normalizable, and hence disappears from the spectrum. The full circle, on the other hand, denotes the end of oscillations in the zero-mode in Eq. (\[s\_wf\_1\]) and (\[s\_wf\_2\]). Lastly, the presence of a single Fermi surface is the key difference in comparison with the same problem of the vortex in a noncentrosymmetric superconductor, in which the spin-orbit coupling is just the chiral kinetic energy, while the ordinary kinetic energy enter the Hamiltonian as $i(k^2/2m-\mu)\Gamma_3\Gamma_5$. It is easy to see that while the operator $A$ still anticommutes with the Hamiltonian, no zero-mode exists in the case for the $s$-wave symmetry of the order parameter[@Lu2].

In the triplet case, the corresponding BdG Hamiltonian has a momentum-dependent mass term. In the Dirac notation, the Hamiltonian is now $\hat{H}=H_D+iMH_D+ih\Gamma_1\Gamma_2-i\mu\Gamma_3\Gamma_5$, where $H_D$ and $M$ are of the same form as, respectively, the kinetic energy and $s$-wave mass in the previous paragraph. It can be shown that the operator $A$ still facilitates the reflection symmetry of the energy spectrum, due to the factor of $i$ in the order-parameter term, which on the other hand is necessary to make it Hermitian. We showed here that the zero-mode does not exist when $\mu=0$, which can also be understood in terms of the symmetry of the zero-mode $\Psi_0$ in Eq. (\[p\_wf\_1\]) under $A$: for positive and negative chemical potential $\mu$, we find that $A\Psi_0=\pm\Psi_0$, respectively. In fact, the disappearance of the zero-mode at $\mu=0$ suggests a topological transition[@DasSarma], because the spectrum in Eq.(\[TripletSpec\]) in the translationally invariant case becomes gapless right at $\mu=0$. In Fig. \[PurePhase\](b) we summarize the effects of Zeeman field on the zero-mode. For $h^2>\mu^2$, there is no zero-mode, as again denoted by the empty circle. Similarly, the full circle at $h^2=\mu^2/(1+v^2_{\Delta}/v^2_F)$, represents the disappearance of oscillations under the overall exponential decaying in the zero-mode.

In a real type-II superconductor, the Zeeman field $h$ is not uniform and is proportional to the distribution of magnetic field $B(r)$, a function riding on the scale of penetration length $l_p$. The obtained zero-mode $\Psi_0$ possesses two length scales $k_F^{-1}=v_F/\sqrt{|\mu^2-h^2|}$ and $\xi=v_F/\Delta$, which are assumed to be much smaller than the penetration length in the present problem. Therefore, it is valid to regard $h(r)$ as a *quasi*-constant throughout the space. In fact, the distribution $B(r)=\frac{\phi_0}{2\pi{l}_p^2}K_0(r/\lambda)$ for an isolated vortex, where $\phi_0$ is the magnetic flux. Hence, it could happen that near the core $|h|>|\mu|$ as a result of $|B|\sim\ln\frac{{l}_p}{r}$ when $r\ll{l}_p$[@deGennesB]. Correspondingly, the zero-mode will behave monotonically near the core, which is just the same as the the modified Bessel function with a length scale of $v_F/|h|$. Nevertheless, the zero-mode will still be sustainable with the asymptotic $h(r\rightarrow\infty)$ within the range in Fig.\[PurePhase\].

![(Color online) Different phases for the zero-mode in the mixed pairing symmetry, as a function of $\mu$ and the parameter $\Delta_{-}$, assuming $\Delta_{+}$ is positive. The dotted axis is to exclude the special cases with $\mu=0$. The grey and blue regions denote different eigenvalues of $A$ for the zero-mode $\Psi_0$ as specified in text. At the dashed line $\mu<0$, $\Delta_- =0$, in between the two phases where $\Delta_-$ changes sign (as denoted by the green arrows) the localized zero-mode does not exist.[]{data-label="MixedPhase"}](PhaseMixed.eps)

The mixed $s$-wave and $p$-wave symmetry is of relevance when it comes to the 2D surface considered here. In this case, the corresponding mass term is $M+iMH_D$, which still anticommutes with $A$ and the principles outlined above can be applied. The amplitudes of the respective pairings are determined by their asymptotic values $\Delta_{\infty}$ and $v_{\Delta}$. First, note that the Fermi surfaces for $\mu>0$ and $\mu<0$ have opposite helicities, and the corresponding zero-mode wave functions in Eq.(\[mixed\_wf\]) have the decaying length as ${\bar
v}_F/\Delta_{+}$ and ${\bar v}_F/|\Delta_{-}|$, respectively. As shown in Fig. \[MixedPhase\], the eigenvalue $\lambda$ of the zero-mode, $A\Psi_0=\lambda\Psi_0$, may be different for opposite signs of the product $\Delta_{+}\Delta_{-}$. Assuming $\mu$ is nonzero and $\Delta_{+}$ is positive, $\lambda$ remains the same when one changes $\mu=|\mu|$ to $\mu=-|\mu|$ if $\Delta_{-}>0$, which is similar to the $s$-wave case. However, $\lambda$ changes sign if $\Delta_{-} < 0$, and this resembles the situation in the $p$-wave case. Besides, the zero-mode stops being normalizable when $\mu<0$ and the parameters $\Delta_{\infty}$ and $k_Fv_{\Delta}$ have $\Delta_{-}$ vanishing, which corresponds to a gapless condition in Eq. (\[MixedSpec\]). This is specified by the dashed line in Fig. \[MixedPhase\]. The meaning of the parameters $\Delta_{\pm}$ is clear when one considers the helicity basis, $\{a^{\dag}_{\rm{k}+},a^{\dag}_{\rm{k}-}\}$, which can be transformed back to the ordinary basis, $\{a^{\dag}_{\rm{k}\uparrow},a^{\dag}_{\rm{k}\downarrow}\}$, with the following,

(

[cccc]{} a\^\_[[k]{}+]{}\
a\^\_[[k]{}-]{}

)= (

[cccc]{} 1 &ie\^[i\_[**k**]{}]{}\
ie\^[-i\_[**k**]{}]{}& 1

) (

[cccc]{} a\^\_[[k]{}]{}\
a\^\_[[k]{}]{}

), where the angle $\phi_{\bf k}$ denotes the direction associated with the momentum on the plane.[@YipTopology; @Santos] Therefore, the time-reversal invariant pairing Hamiltonian can also be written as,[@YipTopology]

H\_=\_[**[k]{}**]{}(i\_[+]{}e\^[-i\_[[k]{}]{}]{}a\^\_[[k]{}+]{}a\^\_[-[k]{}+]{}- i\_[-]{}e\^[i\_[[k]{}]{}]{}a\^\_[[k]{}-]{}a\^\_[-[k]{}-]{})   + h.c., where the two independent order parameters $\Delta_{+,-}$ are real, and corresponding to the pairing order parameters on the two helicity bands. Transforming back to the usual spin basis, the explicit singlet/triplet order parameters are then given by, $\Delta_{s(p)}=(\Delta_{+}\pm\Delta_{-})/2$, respectively, which then yields the relation in Eq. (\[DefDeltaPM\]). It is worthwhile comparing with the noncentrosymmetric superconductors in which the presence of quadratic dispersion term $k^2/2m$ in the normal state leads to appearance of *both* helicity bands at the Fermi level. It follows that those with same (opposite) signs of $\Delta_{\pm}$ resemble the singlet $s$-wave (triplet $p$-wave) superconductors.[@YipTopology] Only the triplet superconductors are associated with topological states such as Andreev bound or zero-energy vortex bound states, which can exist near the topological defects. For superconducting surface Dirac fermions considered here, the appearance of a *single* helicity band near the Fermi level is the key difference.[@Alicea] All the pairing symmetries that we consider possess topological states as long as the superconducting gap is not closed by the Zeeman field or the decreasing chemical potential. Nevertheless, the mixed-parity pairing symmetry cases with same (opposite) signs of $\Delta_{\pm}$ still resemble the pure $s$-wave ($p$-wave) case, manifested by the eigenvalue of antilinear spectrum-reflecting operator $A$ as in Fig.\[MixedPhase\].

In conclusion, we determined the zero-modes for a variety of the Jackiw-Rossi-Dirac Hamiltonians representing pairing of electrons on the surface of a topological insulator in presence of superconducting gaps of different symmetries.

[**Acknowledgment**]{}

The work has been supported by NSC Taiwan (C.K.L.) and NSERC Canada (I.F.H.). C.K.L. also thanks Dr. <PRESIDIO_ANONYMIZED_PERSON> for drawing Ref.[@FuKane] to his attention and the relevant discussions.

[plain]{}

N. Read and D. Green, Phys. Rev. B [**[61]{}**]{}, 10267 (2000).

C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. [**80**]{}, 1083 (2008).

D. A. Ivanov, Phys. Rev. Lett. [**[86]{}**]{}, 268 (2001).

C. Caroli, P. G. de Gennes, and J. Matricon, Phys. Lett. [**[9]{}**]{}, 307 (1964).

N. B. Kopnin and M. M. Salomaa, Phys. Rev. B [**[44]{}**]{}, 9667 (1991).

L. P. Gor’kov and E. I. Rashba, Phys. Rev. Let. [**[87]{}**]{}, 037004 (2001).

C.-K. Lu and S.-K. Yip, Phys. Rev. B [**[78]{}**]{}, 132502 (2008).

Y. Tanaka, T. Yokoyama, A. V. Balatsky, and N. Nagaosa, Phys. Rev. B [**[79]{}**]{}, 060505(R) (2009).

M. Sato and S. Fujimoto, Phys. Rev. B [**79**]{}, 094504 (2009).

X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B [**[81]{}**]{}, 134508 (2010).

S.-K. Yip, J. Low Temp. Phys. [**[160]{}**]{}, 12 (2010).

For a recent review on topological insulator, see M. Z. Hasan and C. L. Kane, arXiv:1002.3895.

L. Fu and C. L. Kane, Phys. Rev. Lett. [**[100]{}**]{}, 096407 (2008).

R. Jackiw and P. Rossi, Nucl. Phys. B [**[190]{}**]{}, 681 (1981).

W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. [**42**]{}, 1698 (1979).

C.-Y. Hou, C. Chamon, and C. Mudry, Phys. Rev. Lett. [**[98]{}**]{}, 186809 (2007).

R. Jackiw and S.-Y. Pi, Phys. Rev. Lett. [**[98]{}**]{}, 266402 (2007).

E. J. Weinberg, Phys. Rev. D [**[24]{}**]{}, 2669 (1981).

I. F. Herbut, Phys. Rev. Lett. [**99**]{}, 206404 (2007).

I. F. Herbut, Phys. Rev. B [**[81]{}**]{}, 205429 (2010); C-Y. Hou, C. Chamon, and C. Mudry, Phys. Rev. B [**81**]{}, 075427 (2010); G. W. Semenoff, arXiv:1005.0572.

B. Seradjeh, H. Weber, and M. Franz, Phys. Rev. Lett. [**[101]{}**]{}, 246404 (2008).

P. Ghaemi and F. Wilczek, arXiv:0709.2626 (unpublished).

I. M. Khaymovich, N. B. Kopnin, A. S. Mel’nikov, and I. A. Shereshevskii, Phys. Rev. B [**[79]{}**]{}, 224506 (2009).

D. L. Bergman and K. Le Hur, Phys. Rev. B [**[79]{}**]{}, 184520 (2009).

I. F. Herbut, Phys. Rev. Lett. [**[104]{}**]{}, 066404 (2010).

B. Uchoa and A. H. Castro Neto, Phys. Rev. Lett. [**[98]{}**]{}, 146801 (2007).

A. M. Black-Schaffer and S. Doniach, Phys. Rev. B [**75**]{}, 134512 (2007).

B. Roy and I. F. Herbut, Phys. Rev. B [**82**]{}, 035429 (2010).

P. Ghaemi, S. Ryu, and D.-H. Lee, Phys. Rev. B [**81**]{}, 081403 (2010).

I. F. Herbut and Chi-Ken Lu, Phys. Rev. B [**82**]{}, 125402 (2010).

S.-K. Yip and A. Garg, Phys. Rev. B [**[48]{}**]{}, 3304 (1993).

L. Santos, T. Neupert, C. Chamon, and C. Mudry, Phys. Rev. B [**81**]{}, 184502 (2010).

J. Linder, Y. Tanaka, T. Yokoyama, A. Sudbo, and N. Nagaosa, Phys. Rev. Lett. [**[104]{}**]{}, 067001 (2010).

P. G. de Gennes, *Superconductivity of metals and alloys* (Benjamin, New York, 1966)

The most general Hamiltonian with the antilinear reflection symmetry of its spectrum can also contain the covariant momentum coupled to gauge field by $p_i\rightarrow{p_i-iA_i\Gamma_3\Gamma_5}$, besides the terms of chemical potential and Zeeman field [@HerbutLu].

A. Altland and M. R. Zirnbauer, Phys. Rev. B [**55**]{}, 1142 (1997); A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B [**78**]{}, 195125 (2008).

J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma, Phys. Rev. Lett. [**104**]{}, 040502 (2010).

J. Alicea, Phys. Rev. B [**81**]{}, 125318 (2010).