Topological Fulde-Ferrell Superfluids in Triangular Lattices

TL;DR

This paper theoretically explores topological Fulde-Ferrell superfluids in a triangular lattice with spin-orbit coupling and in-plane magnetic field, revealing various topological phases with distinct Chern numbers and chiral edge states.

Contribution

It introduces a model for topological FF superfluids on a triangular lattice with SOC, identifying multiple topological phases and their edge states, advancing understanding of topological superfluidity.

Findings

Identification of topological FF phases with Chern numbers ±1 and -2.

Phase boundaries determined by Pfaffian sign change.

Chiral edge states characterized for different topological phases.

Abstract

Fulde-Ferrell (FF) Larkin-Ovchinnikov (LO) phases were proposed for superconductors or superfluids in strong magnetic field. With the experimental progresses in ultracold atomic systems, topological FFLO phases has also been put forward, since it is a natural consequence of realizable spin-orbital coupling (SOC).In this work, we theoretically investigate a triangular lattice model with SOC and in-plane field. By constructing the phase diagram, we show that it can produce topological FF states with Chern numbers, $C=\pm1$ and $C=-2$. We get the phase boundaries by the change of the sign of Pfaffian. The chiral edge states for different topological FF phases are also elucidated.

Tables

Table 1. Table 1: Calculation of Chern number in different parameter region for the single particle system H s subscript 𝐻 𝑠 H_{s} .

parameter region	$o\left({𝐤}_i\right)p\left({𝐤}_i\right)$				$C$
	${𝐤}_1$	${𝐤}_2$	${𝐤}_3^{\pm }$	${𝐤}_4^{\pm }$
upper ellipse	1	-1	-1	-1	-2
lower ellipse	1	1	-1	1	1

Equations

$$H=\sum _{⟨i,j⟩}{\varphi }_i^{†}{H}_s{\varphi }_j+H_{\text{int}},$$

$$H_s=t^{i,j}{\sigma }_z+i{t}_{\mathrm{so}}^{i,j}({\mathbf{v}}_{i,j}\times \mathbf{\sigma }{)}_z+(\mathbf{h}\cdot \mathbf{\sigma }-\mu ){\delta }_{i,j},$$

$$H_{\text{int}}=-U\sum _i{c}_{i,\uparrow }^{†}{c}_{i,\uparrow }{c}_{i,\downarrow }^{†}{c}_{i,\downarrow },$$

$${\mathbf{\Delta }}_i=U⟨c_{i,\downarrow }{c}_{i,\uparrow }⟩=\Delta e^{i\mathbf{Q}\cdot {\mathbf{R}}_i}$$

$$H=\frac{1}{2}\sum _{\mathbf{k}}{\Psi }_{\mathbf{k}}^{†}{H}_{\text{BdG}}{\Psi }_{\mathbf{k}}+N\left(\frac{{\left|{\mathbf{\Delta }}_i\right|}^2}{U}-\mu \right),$$

$$H_{\text{BdG}}=$$

$$+(d_{\mathbf{k}}+h_x){\tau }_0\otimes {\sigma }_x+e_{\mathbf{k}}{\tau }_z\otimes {\sigma }_y+f_{\mathbf{k}}{\tau }_0\otimes {\sigma }_y$$

$$+\Delta {\tau }_x\otimes {\sigma }_0-\mu {\tau }_z\otimes {\sigma }_0$$

$$\Omega =$$

$$-k_{B}T\sum _{\alpha ,\mathbf{k}}\ln \left(1+e^{E_{\alpha ,\mathbf{k}}^{-}/(k_{B}T)}\right)\text{ ,}$$

$$C=\sum _{\alpha }\frac{1}{2\pi }\int d{k}_{x}d{k}_y\mathbf{z}\cdot {\nabla }_{\mathbf{k}}\times {\mathbf{A}}_{\alpha }^{-}(\mathbf{k})\text{,}$$

$$H_s=\sum _{\mathbf{k}}(d_x{\sigma }_x+d_y{\sigma }_y+d_z{\sigma }_z),$$

$$|\mathbf{d}(\mathbf{k})|=0,$$

$$\frac{h_x^2}{3t_{\mathrm{so}}^2}+{\left(\frac{h_z}{t}-\frac{5}{2}\right)}^2=\frac{1}{4}$$

$$\frac{4h_x^2}{3t_{\mathrm{so}}^2}+{\left(\frac{h_z}{t}+2\right)}^2=16,$$

$$C=\frac{1}{2}\sum _{{\mathbf{k}}_i}o({\mathbf{k}}_i)p({\mathbf{k}}_i),$$

$$P({\mathbf{k}}^{\prime })=\mathrm{sgn}[\mathrm{Pf}[W({\mathbf{k}}^{\prime })]],$$

$$[d({\mathbf{k}}^{\prime })+h_x{]}^2+[a({\mathbf{k}}^{\prime })+h_z{]}^2-\Delta (h_x,h_z{)}^2={\mu }^2-f({\mathbf{k}}^{\prime }{)}^2.$$

$$a_{\mathbf{k}}$$

$$b_{\mathbf{k}}$$

$$c_{\mathbf{k}}$$

$$d_{\mathbf{k}}$$

$$e_{\mathbf{k}}$$

$$f_{\mathbf{k}}$$

$${\mathbf{k}}_1$$

$${\mathbf{k}}_2$$

$${\mathbf{k}}_3^{\pm }$$

$${\mathbf{k}}_4^{\pm }$$

Full text

I Introduction

Cooper pair was first proposed in 1956 to explain superconductors Cooper . It describes a pair of fermions bound together due to attractive interaction. The fermions have opposite momentum so that the pair has zero momentum totally. However, Cooper pair with finite center-of-mass momentum may also exist in the presence of strong magnetic field. This consideration led to an exotic superconductor with inhomogeneous order parameter in real space, known as Fulde-Ferrell Fulde and Larkin-Ovchinnikov (FFLO) phases Larkin . There are two types of FFLO phases: phase modulated FF state and spatial modulated LO state. In the past two decades, FFLO phases attract tremendous interests in both experiment and theory Casalbuoni ; Kenzelmann ; Liao ; Li ; Hu1 ; Radzihovsky . But only ambiguous experiment evidences from heavy-fermion superconductors and organic superconductors are available Matsuda ; Shimahara ; Koutroulakis .

On the other hand, topology is also a hot topic in condensed matter field for several decades Xiao ; Hasan ; Qi . Recently, the spin-orbit coupling (SOC) has been realized by ultracold atoms as condensates or on an optical lattice Lin ; Wang ; Cheuk ; Wu , which paves the way to the topological FFLO states Liu ; Qu ; Zhang ; Chen ; Cao ; Hu2 . Theoretical researches show that topological FFLO states can be induced in one and two dimensional Fermi gas. And according to the bulk-edge correspondence, edge modes are supported when there are boundaries Hatsugai ; Bernvig .

In two dimensions, topological FF state with Chern number $C=1$ can be produced for cold atoms in a square lattice Qu ; Chen . To realize topological state with higher Chern number, one can resort to complicated hoppings Sarma or lattices Huang . Nonetheless, the simple triangular lattice favors some topologically nontrivial states, it can produce topological state with higher Chern number by merely the nearest-neighbor hoppings Kai2 ; Iskin . In this work, we investigate a system with SOC and in-plane Zeeman field on the triangular lattice to achieve topological FF states. We found, in the noninteracting case, our system with nearest-neighbor hoppings supports the gapped Chern insulators with Chern number $C=1$ as well as $C=-2$. The boundaries of both of the topological phases are ellipses in a two-parameter plane constituted by the in-plane and out-plane fields. In the presence of on-site attractive s-wave interactions, we solve the system by self-consistent equations. The non-uniform FF superfluid states are obtained in a large area. And we confirm that all the FF states are topologically nontrivial. Their Chern numbers are $C=1,-1$ or $-2$. We use a set of signs of Pfaffians at high symmetry points in the first Brillouin zone (1st BZ) to characterize the different topological phases. Each one of these signs of Pfaffians changes with the energy gap closing and reopening at the corresponding point. We calculate the chiral edge states, whose wave functions are spatially localized at the edges in open boundary situation, which can help us to confirm that the bulk is in a topologically nontrivial state. Different FF phases exhibit different pairs of chiral edge states. The edge current is directly determined by Chern number, i.e. the summation of chiralities of the edge modes.

This paper is organized as follows. In Sec. II, we introduce a system in the triangular lattice with SOC and in-plane field and solve it by self-consistent method. In Sec. III, we construct the phase diagram in the noninteracting and interacting case. The topological FF phases with $C=-1,1$ and $-2$ are obtained in the presence of interaction. The chiral edge states in distinct topological FF phases are also illustrated. At last, we give a brief summary.

II Model Hamiltonian

We consider a two-dimensional (2D) spin-orbit coupled fermionic gas trapped in a triangular lattice subjected to an external magnetic field. In the position space, the model Hamiltonian reads

[TABLE]

where $\phi_{i}^{\dagger}=\left(c_{i,\uparrow},c_{i,\downarrow}\right)$ with $c_{i,\sigma}^{\dagger}(c_{i,\sigma})$ being the creation (annihilation) operator for the spin $\sigma$ particle at site $i$, $H_{s}$ and $H_{\rm int}$ represent the single-particle and interaction Hamiltonian, respectively, and the summation is over the nearest-neighbor pairs. More specifically, we assume that the single-particle Hamiltonian takes the form

[TABLE]

where $\sigma_{\alpha}$ ($\alpha=x,y,z$) are a Pauli matrices, $t^{i,j}$ is the hopping matrix elements which, when combined with $\sigma$, gives rise to opposite signs for the spin-up and -down particles Xiongjun , $t_{\rm so}^{i,j}$ are the Rashba SOC coefficients, $\mathbf{v}_{i,j}$ are the vectors connecting lattice sites $i$ and $j$, $\mathbf{h}=(h_{x},0,h_{z})$ is the magnetic field which includes an in-plane component $h_{x}$ along the $x$ direction and an out-of-plane one $h_{z}$, and $\mu$ is chemical potential. The $s$-wave interaction Hamiltonian can be expressed as

[TABLE]

where $U>0$ represents attractive interaction. Finally, in the triangular lattice, the lattice basis vectors are defined as $\mathbf{a}_{1}=(1,0)$, $\mathbf{a}_{2}=(-1/2,\sqrt{3}/2)$, and $\mathbf{a}_{3}=(1/2,\sqrt{3}/2)$ which, as shown in Fig. 1(a), represent the displacements from one site to the nearest sites. Correspondingly, we obtain a hexagonal first Brillouin zone (1BZ) with reciprocal lattice basis vectors $\mathbf{G}_{1}=(0,4\pi/\sqrt{3})$ and $\mathbf{G}_{2}=(2\pi,2\pi\sqrt{3})$.

The attractive interaction may lead to the Bardeen-Cooper-Schrieffer (BCS) superfluid. Particularly, in the presence of the in-plane Zeeman field along the $x$ direction, the Fermi surface becomes asymmetric along the $y$ axis. Consequently, the BCS pairs may carry a finite center-of-mass momentum $Q_{y}$ along the $y$ direction Zheng . Such state is described by the FF order parameter, which, in the position space, is defined as

[TABLE]

with $\mathbf{Q}=(0,Q_{y})$ momentqy . After transformed into the momentum space, the mean-field Hamiltonian in the Nambu representation reads

[TABLE]

where $\Psi_{\mathbf{k}}=\left(c_{\mathbf{k}+\mathbf{Q}/2,\uparrow},c_{\mathbf{k}+\mathbf{Q}/2,\downarrow},-c_{-\mathbf{k}+\mathbf{Q}/2,\downarrow}^{\dagger},c_{-\mathbf{k}+\mathbf{Q}/2,\uparrow}^{\dagger}\right)^{T}$ is the Nambu spinor with ${\mathbf{k}}=(k_{x},k_{y})$, $N$ is the total number of lattice sites, and the Bogoliubov-de Gennes (BdG) Hamiltonian is

[TABLE]

with $\otimes$ being the Kronecker product and $\tau_{\alpha}$ ($\alpha=x,y,z$) being the Pauli matrices acting on the particle-hole space. The explicit expressions for the elements of $H_{\text{BdG}}$ matrix, $a_{\mathbf{k}}$, $b_{\mathbf{k}}$, $c_{\mathbf{k}}$, $e_{\mathbf{k}}$, $f_{\mathbf{k}}$, and $g_{\mathbf{k}}$, can be found in the Appendix A. We note that the BdG Hamiltonian only possesses a particle-hole symmetry, $\Xi H_{\text{BdG}}(\mathbf{k})\Xi^{-1}=\Lambda H_{\text{BdG}}^{\ast}(\mathbf{k})\Lambda=-H_{\text{BdG}}(-\mathbf{k})$, where $\Xi=\Lambda K$ with $\Lambda=i\sigma_{y}\otimes\tau_{y}$ and $K$ being the complex conjugation operator. It can be verified that $\Xi^{2}=1$. Apparently, the system belongs to the class D according to Wigner-Dyson symmetry classification of random matrix Altland which has a topological invariant Z in two dimension Schnyder ; Ryu ; Kiteav1 .

Following the standard treatment, we diagonalize the BdG Hamiltonian, $H_{\rm BdG}|\psi_{\alpha}^{\nu}({\mathbf{k}})\rangle=E_{\alpha,{\mathbf{k}}}^{\nu}|\psi_{\alpha}^{\nu}({\mathbf{k}})\rangle$, which leads to the quasiparticle eigenenergy $E_{\alpha,{\mathbf{k}}}^{\nu}$ and the quasiparticle wave function $|\psi_{\alpha}^{\nu}({\mathbf{k}})\rangle$, here $\nu=\pm$ represent the particle ($+$) and hole ($-$) bands, $\alpha=1$ and $2$ denote the upper ($1$) and lower ($2$) helicity branches. Now, the thermodynamic potential at temperature $T$ can be calculated through

[TABLE]

where $k_{B}$ the Boltzmann constant and the summation is restricted to the hole bands ($\nu=-$) due to the inherent particle-hole symmetry in the Nambu spinor representation. The order parameter can be numerically determined by the mean-field saddle equations $\partial\Omega/\partial\Delta=0$, $\partial\Omega/\partial Q_{y}=0$, as well as the equation for the conservation of the total particle number $\partial\Omega/\partial\mu=-N$. Finally, the topological property of the system is characterized by the Chern number Thouless ; Xiao

[TABLE]

where $\mathbf{A}_{\alpha}^{-}(\mathbf{k})=i\left\langle\psi_{\alpha}^{-}(\mathbf{k})\right|{\nabla}_{\mathbf{k}}\left|\psi_{\alpha}^{-}(\mathbf{k})\right\rangle$ is the vector potential.

III Results

In this section, we present the results about the zero-temperature quantum phases of our model. For convenience, we assume that the hopping matrix element is site independent, i.e., $t=t^{i,j}$ and $t_{\rm so}=t_{\rm so}^{i,j}$. Moreover, we select $t$ as the energy unit such that the model Hamiltonian is completely specified by the parameters $t_{\rm so}/t$, $h_{x}/t$, $h_{z}/t$, $U/t$, and $\mu/t$. In below, the value of the chemical potential is fixed at $\mu/t=1$ for simplicity.

To start, let us briefly discuss the quantum phases of our model in the absence of interaction. In the momentum space, the single-particle Hamiltonian reduces to

[TABLE]

where $d_{x}=\sqrt{3}t_{\rm so}(\sin{k_{3}-\sin{k_{2}}})+h_{x}$, $d_{y}=t_{\rm so}(2\sin{k_{1}}-\sin{k_{2}}-\sin{k_{3}})$, and $d_{z}=t(\cos{k_{1}}+\cos{k_{2}}+\cos{k_{3}})+h_{z}$ with $k_{\alpha}=\mathbf{k}\cdot\mathbf{a}_{\alpha}$. Due to the breaking of the time-reversal and chirality symmetries, the single-particle Hamiltonian belongs to the symmetry class C. Therefore, the possible topologically nontrivial ground state in 2D system is characterized by Z invariant Schnyder ; Ryu ; Kiteav1 . To identify the topological state, we focus on the gapless points defined by

[TABLE]

where ${\mathbf{d}}\equiv(d_{x},d_{y},d_{z})$. It can be shown that Eq. (7) gives rise to two ellipses,

[TABLE]

and

[TABLE]

which divide the $h_{x}h_{z}$ parameter plane into three regions. The topological property of each region is determined by the Chern number $C=\frac{1}{4\pi}\int_{\rm 1BZ}dk_{x}dk_{y}\mathbf{n}\cdot\partial_{k_{x}}\mathbf{n}\times\partial_{k_{y}}\mathbf{n}$, where $\mathbf{n}=\mathbf{d}/|\mathbf{d}|$ is a unit vector.

In Fig. 2, we plot the phase diagram of the Hamiltonian (6) for $t_{\rm so}/t=1$. Without loss of generality, only the result for $h_{x}>0$ is shown. As can be seen, the parameter regions enclosed by two ellipses are topologically nontrivial. Particularly, the Chern number inside the upper ellipse is $-2$. These topologically nontrivial states can be understood by noting that the Chern number can be reexpressed as Tretiakov ; Clement

[TABLE]

where $\mathbf{k}_{i}$ are the roots of the equations $d_{x}(\mathbf{k})=d_{y}(\mathbf{k})=0$ (see Appendix B), $o(\mathbf{k}_{i})$ is the chirality of the vector field $(d_{x}(\mathbf{k}),d_{y}(\mathbf{k}))$ around $\mathbf{k}_{i}$, and $p(\mathbf{k}_{i})=\text{sgn}(d_{z}(\mathbf{k}_{i}))$ is the polarity at $\mathbf{k}_{i}$. In Tab. 1, we summarize the contribution of each root to the Chern number.

We now turn to study the superfluid phases of the system by taking into account the attractive $s$-wave interaction. Fig. 3 summarizes the quantum phases in the $h_{x}h_{z}$ parameter plane with $U/t=5$ and $\mu/t=1$. Here the quantum phases are characterized by the superfluid order parameter $\Delta$, the center-of-mass momentum $Q_{y}$, the Chern number $C$, and the chemical potential $\mu$. Specifically, when $\Delta=0$, we may either have an insulating (IN) phase or a normal gas (NG) phase. For the former state, the chemical potential locates in the band gap; while, for the latter state, it lies in the band such that the excitation is gapless Xu . Moreover, an IN phase is Chern insulator (CI) if it is topologically nontrivial ($C=1$ here). Next, a superfluid state ($\Delta\neq 0$) with nonzero $Q_{y}$ is denoted by ‘FF’. In our model, the appearance of FF states is due to the in-plane magnetic field $h_{x}$ that deforms the Fermi surface Zheng . In fact, the standard Bardeen-Cooper-Schrieffer superfluid ($Q_{y}=0$) only exists at $h_{x}=0$. Additionally, all superfluid phases are found to be topologically nontrivial and they are further specified by the letter ‘t’ and by Chern numbers in the square bracket. Finally, the dashed line in Fig. 3 marks the boundary between gapped and gapless superfluid phases. More specifically, for small $h_{x}$, the energy of the lower helicity particle branch $E_{2,{\mathbf{k}}}^{+}$ is always positive for the gapped superfluid phases. However, $E_{2,{\mathbf{k}}}^{+}$ may become less than zero as $h_{x}$ is increased, which leads to the gapless superfluid phases Fan ; Deng .

In Fig. 4, we plot the $h_{z}$ dependence of the quasiparticle band gap $E_{g}\equiv{\rm max}\{0,{\rm min}(E_{2,{\mathbf{k}}}^{+})\}$, the superfluid order parameter $\Delta$, and the pairing momentum $Q_{y}$ for a fixed $h_{x}/t=0.3$. As can be seen, for small $h_{z}$, the Fermi surface locates in the band gap such that the system remains in the insulating phase. Then as $h_{z}$ is increased, $E_{2,{\mathbf{k}}}^{+}$ moves downward, while $E_{2,{\mathbf{k}}}^{-}$ moves oppositely, which leads to a vanishing $E_{g}$, signaling the entering of the NG phase. As one further increases $h_{z}$, the tFF phases emerge. Within the superfluid phases, the energy gap closes and reopens whenever a topological phase transition is encountered. Eventually, at very large $h_{z}$, the large population difference between the spin-$\uparrow$ and -$\downarrow$ particles prohibits the formation of the superfluid pairs. As a result, the system falls into insulting phases again.

The phase boundaries between topological superfluid phases in Fig. 3 can be alternatively determined by considering the high symmetry points $\mathbf{k}^{\prime}$ that satisfy $\Xi H_{\text{BdG}}(\mathbf{k}^{\prime})\Xi^{-1}=-H_{\text{BdG}}(\mathbf{k}^{\prime})$. To this end, we introduce an auxiliary matrix $W(\mathbf{k}^{\prime})\equiv H_{\text{BdG}}(\mathbf{k}^{\prime})\Lambda$ Kiteav2 ; Ghosh that is antisymmetric, i.e., $W^{T}({\mathbf{k}}^{\prime})=-W({\mathbf{k}}^{\prime})$. A topological index can then be defined as

[TABLE]

where ${\rm Pf}[W(\mathbf{k}^{\prime})]$ denotes the Pfaffian of $W(\mathbf{k}^{\prime})$. Since ${\rm Pf}[W(\mathbf{k}^{\prime})]=\pm\sqrt{\det[H_{\text{BdG}}(\mathbf{k}^{\prime})]}$, $P(\mathbf{k}^{\prime})$ will never change its sign unless the energy gap at $\mathbf{k}^{\prime}$ is closed, indicating that $P(\mathbf{k}^{\prime})$ is indeed topologically protected. This also suggests that the phase boundaries are determined by the condition $P({\mathbf{k}}^{\prime})=0$ or, formally,

[TABLE]

In our model, the high symmetry points include $\mathbf{\Gamma}=(0,0)$, $\mathbf{M}_{1}=(0,2\pi/\sqrt{3})$, $\mathbf{M}_{2}=(\pm\pi,\pi/\sqrt{3})$, and $\mathbf{\Sigma}=(\pm 2\cos^{-1}[-\cos(\sqrt{3}Q_{y}/4)/2],0)$, as shown in Fig. 1(a). Here $P(\mathbf{k}^{\prime}=\Gamma)=-1$ remains unchanged and is irrelevant to the topological phase transitions. However, it can be verified that other high symmetry points $\mathbf{M}_{1}$, $\mathbf{M}_{2}$, and $\mathbf{\Sigma}$ leads to the three phase boundaries between the topological superfluid phases through Eq. (12).

To gain more insight into these topological phases, we further consider the edge modes of the system by imposing a hard-wall boundary condition along the $y$ direction. Fig. 5 shows the quasiparticle spectra of distinct tFF phases and the probability distribution of the corresponding edge states. Consider, for example, the tFF[-1] phase in Fig. 5(a), the lines of the edge modes cross the Fermi level three times. To distinguish each edge state clearly, we plot a horizontal dashed line slightly above the Fermi level to reveal the Hall current. It shows that the dashed line intersects with the lines of edge modes at A, B, C, D, E, and F successively. Correspondingly, the probability distributions over the lattice sites for these states are plotted in the right panel of Fig. 5(a), which shows that states A, D, and E locate at the lower edge and states B, C, and F at the upper edge. These edge states also give rise to edge currents, $I\propto-\partial E/\partial k_{x}$, at the upper or lower edges. To proceed further, let us focus on the edge states A and F, which locate at upper and lower edges, respectively. Due to the translational symmetry along the $x$ direction, the band is symmetrical about the $k_{x}$ axis, which suggests that the group velocities of the states A and F have the same magnitude but with opposite signs. As a result, states A and F form a clockwise current loop at the edges and contribute a Chern number, $C_{\rm AF}=-1$, to the system. Similar analysis can be carried out for other pairs of states and leads to the Chern numbers $C_{\rm CD}=1$ and $C_{\rm BE}=-1$. Consequently, the total Chern number of the system is $C=C_{AF}+C_{CD}+C_{BE}=-1$. Likewise, from Fig. 5(b) and (c), it can be verified that the the Chern numbers of the tFF[1] and tFF[-2] phases are indeed $C=1$ and $-2$, respectively.

IV Summary

In summary, we have investigated a triangular lattice with SOC and in-plane magnetic field. We construct the phase diagram in the noninteracting andinteracting cases respectively. In the noninteracting case, we distinguish different phases by Chern numbers. The Chern insulators with Chern number $C=1$ and $C=-2$ are limited within two ellipses. In the interacting case, we get a very rich phase diagram including gapped and gapless topological FF superfluid phases with Chern numbers $C=-1,1,-2$. The topological FF phases are characterized by phase modulating order parameters and non-zero Chern numbers. We disclose edge states in different FF states and also discuss the connection between the edge state and the Chern number.

ACKNOWLEDGMENTS

We acknowledge useful discussions with Yan He and Yuangang Deng. This work was supported by the NSFC (Grants No. 11074177, No. 11421063, and No. 11674334) and SRF for ROCS SEM (20111139-10-2).

Appendix A Components of the BdG Hamiltonian

The components of the BdG Hamiltonian (3) are

[TABLE]

where $Q=\sqrt{3}Q_{y}/4$. It can be seen that $a_{\mathbf{k}}$, $d_{\mathbf{k}}$, and $e_{\mathbf{k}}$ are even functions of $\mathbf{k}$; while $b_{\mathbf{k}}$, $c_{\mathbf{k}}$, and $f_{\mathbf{k}}$ are odd functions.

Appendix B Roots of the equations $d_{x}(\mathbf{k})=d_{y}(\mathbf{k})=0$

Explicitly, it can be easily evaluated that

[TABLE]

where $\eta={2h_{x}}/{\sqrt{3}t_{\rm so}}$. The condition for existence of roots $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ is $\left|\eta\right|/4\leq 1$ and $\left|\eta\right|\leq 1$ for $\mathbf{k}_{3}^{\pm}$ and $\mathbf{k}_{4}^{\pm}$.