Global phase diagram, possible chiral spin liquid, and topological superconductivity in the triangular Kitaev–Heisenberg model

Recently, there has been enormous interest in the physics of the spin-1/2 Kitaev model on a honeycomb lattice [1], which has an exact Z₂ spin-liquid (SL) ground state (GS) supporting fractionalized excitations. One possible route to realize this highly anisotropic spin model is to include a strong relativistic spin-orbit coupling (SOC) in Mott insulators [2, 3]. Indeed, the interplay of SOC and electron interactions [4–8] gives rise to many novel phases [9–13], especially for the so-called relativistic Mott insulators (RMIs) whose physics may drastically differ from that of Mott insulators with weak SOC (e.g., cuprates) [3, 4, 14, 15]. Of particular interest are $5d$ transition metal oxides, such as iridates A₂IrO₃ (A = Na, Li) [16–18], where Na₂IrO₃ is interpreted as a novel RMI [19] and may also host the quantum spin Hall effect [20, 21]. The Kitaev–Heisenberg (KH) model on a honeycomb lattice, which has a rich phase diagram containing unconventional magnetic and Kitaev SL phases [22, 23], has been proposed to capture the low-energy properties of A₂IrO₃ [24, 25]. Meanwhile, experiments confirm a long-range zigzag spin order in Na₂IrO₃ [26–29], which is a natural GS of the KH model [22]. In addition, there are also studies on the $4d$ compound Li₂RhO₃, suggesting Li₂RhO₃ as a possible RMI with a spin-glass GS [30]. Theoretical studies also show that carrier doping into RMIs can induce unconventional superconducting pairings as well as the topological superconductivity (SC) [31–36].

In fact, the KH model can be generalized to the triangular lattice [37, 38]. Similar to the microscopic origin of the honeycomb KH model for A₂IrO₃, the triangular KH model can emerge from a class of ABO₂-type (where A and B are alkali and transition metal ions, respectively) layered compounds [2, 3], due to the joint effect of strong SOC, Coulomb interaction, orbital degeneracy, $t_{2g}^{5}$ configuration, and 90°-bonding geometry [22, 24]. However, to the best of our knowledge, no exact solution of the spin-1/2 Kitaev model on a triangular lattice has been obtained so far. Therefore, it remains conceptually interesting to investigate whether the SL could exist as a possible GS of the triangular KH model.

In this paper, by combining the numerical exact diagonalization (ED) calculation with the four-sublattice transformation (FST) [2, 22, 24] analysis, we demonstrate one possible exotic phase near the antiferromagnetic (AF) Kitaev point and four magnetically ordered phases, including a collinear stripe pattern and a noncollinear spiral pattern in the global phase diagram. For the exotic phase, resorting to the Schwinger-fermion mean-field (MF) method [39, 40], we find two local minimum solutions with an s-wave and a $d+{\rm i}d$-wave pairings, respectively, where the latter has a lower MF energy and is further identified as a Z₂ chiral SL state with a Chern number ±2. The effect of finite hole-doping is analyzed using the slave-boson MF theory, and a time-reversal (TR) invariant topological p-wave SC, an s-wave SC, and a ${{d}_{{{x}^{2}}-{{y}^{2}}}}+{\rm i}{{d}_{xy}}$-wave SC are found in the phase diagrams.

Let us begin with the spin-1/2 KH model defined on a triangular lattice

$$\begin{eqnarray}&&{{H}_{KH}}=\mathop{\sum }\limits_{\langle i,j\rangle }(-{{J}_{K}}S_{i}^{{{\alpha }_{ij}}}S_{j}^{{{\alpha }_{ij}}}+{{J}_{H}}{{{\boldsymbol{S}} }_{i}}\cdot {{{\boldsymbol{S}} }_{j}}),\end{eqnarray} \tag{ 1 }$$

where the index ${{\alpha }_{ij}}$ takes values $x,y,$ or z depending on the direction of the nearest-neighbor (NN) bond $\langle i,j\rangle$ (see figure 1(a)). This model consists of spin-anisotropic Kitaev interactions (the first term) and spin-isotropic Heisenberg interactions (the second term).

Zoom In Zoom Out Reset image size

For convenience, we parametrize the coupling constants in equation (1) by introducing the energy scale $A=\sqrt{J_{K}^{2}+J_{H}^{2}}$ and the angle ϕ via ${{J}_{K}}=A{\rm sin} \phi$ and ${{J}_{H}}=A{\rm cos} \phi$, and let ϕ vary from 0 to $2\pi$ to cover the global phase diagram.

To detect the quantum phase transitions, we perform a Lanczos ED calculation of the GS energy of the Hamiltonian (1) on a 24-site cluster with periodic boundary conditions, and the results are presented in figure 2. As indicated by the dashed lines in figure 2, the second derivative of the GS energy with respect to ϕ reveals five distinct phases separated by five transition points $\phi \approx 0.14\pi ,0.5\pi ,1.31\pi ,1.4\pi$, and $1.73\pi$, best visualized using the ϕ circle as in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At $\phi =0$, we are left with the AF Heisenberg model (${{J}_{K}}=0,{{J}_{H}}\gt 0$) exhibiting the 120° Neel order [41–44]; at $\phi =\pi$, equation (1) corresponds to the ferromagnetic (FM) Heisenberg model (${{J}_{K}}=0,{{J}_{H}}\lt 0$) with an FM GS. In addition to these well-known phases, we observe two more magnetic orders by using the so-called FST approach, which is in fact a spin-rotation transformation. Specifically, it is instructive to divide the triangular lattice into four sublattices (see figure 1(a)) and introduce the rotated spin operators $\tilde{{\boldsymbol{S}} }$: ${{\tilde{{\boldsymbol{S}} }}_{0}}={{{\boldsymbol{S}} }_{0}}$ in sublattice 0, ${{\tilde{{\boldsymbol{S}} }}_{1}}=(S_{1}^{x},-S_{1}^{y},-S_{1}^{z})$ in sublattice 1, ${{\tilde{{\boldsymbol{S}} }}_{2}}=(-S_{2}^{x},S_{2}^{y},-S_{2}^{z})$ in sublattice 2, and ${{\tilde{{\boldsymbol{S}} }}_{3}}=(-S_{3}^{x},-S_{3}^{y},S_{3}^{z})$ in sublattice 3. This transformation results in the new spin-$\tilde{{\boldsymbol{S}} }$ Hamiltonian of the same form as equation (1), but with the effective couplings ${{\tilde{J}}_{K}}={{J}_{K}}-2{{J}_{H}}$ and ${{\tilde{J}}_{H}}=-{{J}_{H}}$. For the angles, the mapping reads as ${\rm tan} \tilde{\phi }=-{\rm tan} \phi +2$.

Therefore, at $\phi ={\rm arctan} 2$ and ${\rm arctan} 2+\pi$, FST maps equation (1) to a FM and an AF Heisenberg Hamiltonians for $\tilde{{\boldsymbol{S}} }$, respectively. Thus, transforming back to the original spin basis, we obtain a collinear stripe order at $\phi ={\rm arctan} 2$ and a noncollinear spiral order at $\phi ={\rm arctan} 2+\pi$, as shown in figure 3, where the spiral pattern has an enlarged magnetic unit cell containing 12 lattice sites.

Here, we would like to point out that the ED results for both the two pairs of transition points $(0.14\pi ,1.31\pi )$ and $(1.4\pi ,1.73\pi )$ match the FST mapping ${\rm tan} \tilde{\phi }=-{\rm tan} \phi +2$ very well (as indicated by the blue lines inside the ϕ circle, see figure 3), and the isolated transition point $\phi =0.5\pi$ is also consistent with FST, i.e., it is mapped to itself under FST.

Thus far, the numerical ED together with the analytical FST has revealed four magnetically ordered phases in the global phase diagram. However, we are still left with one possible exotic phase (i.e., the phase corresponds to the $1.4\pi \lt \phi \lt 1.73\pi$ arc of the ϕ circle) near the AF Kitaev point (i.e., the point at $\phi =3\pi /2$).

Under FST, it can be seen that each point inside the region $1.4\pi \lt \phi \lt 1.73\pi$ is mapped to another point still inside this region, where the AF Kitaev point ($\phi =3\pi /2$) is exactly invariant. Therefore, if the exotic phase ($1.4\pi \lt \phi \lt 1.73\pi$) is magnetically ordered, it would be somewhat subtle in the sense that the corresponding spin configuration must be invariant under FST. An alternative is a magnetically disordered state that may be favored by the frustration and quantum fluctuations embedded in equation (2).

As a representative of the phase in the region $1.4\pi \lt \phi \lt 1.73\pi$ and for simplicity, let us start from the AF Kitaev point of the Hamiltonian (1) to explore the nature of the exotic phase, say

$$\begin{eqnarray}&&{{H}_{K}}=-{{J}_{K}}\mathop{\sum }\limits_{\langle i,j\rangle}S_{i}^{{{\alpha }_{ij}}}S_{j}^{{{\alpha }_{ij}}},\;{\tt with}\;{{J}_{K}}\lt 0.\end{eqnarray} \tag{ 2 }$$

The model (2) is TR invariant. It is also invariant under the lattice-translation and -inversion and under the spin-rotation by π about the spin x, y, or z axes.

Unlike the honeycomb Kitaev model [1], the triangular Kitaev model (2) can not be solved exactly through the Majorana fermionization of spin-$1/2$ operators [45]. Instead, let us take the standard Schwinger-fermion representation of spin-1/2 operators: $S_{i}^{\alpha }=\frac{1}{2}{{\sum }_{\sigma ,{{\sigma }^{{\prime} }}=\uparrow ,\downarrow }}f_{i\sigma }^{\dagger }\tau _{\sigma \sigma ^{\prime} }^{\alpha }{{f}_{i{{\sigma }^{{\prime} }}}}$, with fermionic spinons ${{f}_{i\sigma }}$ and Pauli matrices ${{\tau }^{\alpha }}$. The physical-spin Hilbert space is then recovered by imposing the on-site constraint ${{\sum }_{\sigma }}f_{i\sigma }^{\dagger }{{f}_{i\sigma }}=1$, which can be enforced by a site-independent Lagrange multiplier λ at the MF level. Now, the spin quadratic terms in equation (2) can be decoupled into several different channels as

$$\begin{eqnarray}\begin{array}{rcl} S_{i}^{x}S_{j}^{x} & = & -\frac{1}{8}\mathop{\sum }\limits_{\sigma }t_{\sigma ,ij}^{\dagger }{{t}_{-\sigma ,ij}}+\frac{1}{8}\left( -\Delta _{0,ij}^{\dagger }{{\Delta }_{0,ij}} \right. \\ {} & {} & -\left. 2\Delta _{1,ij}^{\dagger }{{\Delta }_{1,ij}}+2\Delta _{2,ij}^{\dagger }{{\Delta }_{2,ij}}+\Delta _{3,ij}^{\dagger }{{\Delta }_{3,ij}} \right), \\ S_{i}^{y}S_{j}^{y} & = & -\frac{1}{8}\mathop{\sum }\limits_{\sigma }t_{\sigma ,ij}^{\dagger }{{t}_{-\sigma ,ij}}+\frac{1}{8}\left( -\Delta _{0,ij}^{\dagger }{{\Delta }_{0,ij}} \right. \\ {} & {} & +\left. 2\Delta _{1,ij}^{\dagger }{{\Delta }_{1,ij}}-2\Delta _{2,ij}^{\dagger }{{\Delta }_{2,ij}}+\Delta _{3,ij}^{\dagger }{{\Delta }_{3,ij}} \right), \\ S_{i}^{z}S_{j}^{z} & = & -\frac{1}{8}\mathop{\sum }\limits_{\sigma }t_{\sigma ,ij}^{\dagger }{{t}_{\sigma ,ij}}+\frac{1}{8}\left( -\Delta _{0,ij}^{\dagger }{{\Delta }_{0,ij}} \right. \\ {} & {} & +\left. 2\Delta _{1,ij}^{\dagger }{{\Delta }_{1,ij}}+2\Delta _{2,ij}^{\dagger }{{\Delta }_{2,ij}}-\Delta _{3,ij}^{\dagger }{{\Delta }_{3,ij}} \right), \\ \end{array}\end{eqnarray} \tag{ 3 }$$

where ${{t}_{\sigma ,ij}}=f_{i\sigma }^{\dagger }{{f}_{j\sigma }}$, ${{\Delta }_{0,ij}}=\frac{1}{\sqrt{2}}{{\sum }_{\sigma ,{{\sigma }^{{\prime} }}}}{{f}_{i\sigma }}{\bf i}\tau _{\sigma \sigma ^{\prime} }^{2}{{f}_{j{{\sigma }^{{\prime} }}}}$ (singlet pairing), ${{\Delta }_{l,ij}}=\frac{1}{\sqrt{2}}{{\sum }_{\sigma ,{{\sigma }^{{\prime} }}}}{{f}_{i\sigma }}{{[{\bf i}{{\tau }^{l}}{{\tau }^{2}}]}_{\sigma \sigma ^{\prime} }}{{f}_{j{{\sigma }^{{\prime} }}}}$ (triplet pairing), with ${\bf i}=\sqrt{-1}$, $l=1,2,3$, and the Pauli matrices ${{\tau }^{l}}$.

To proceed, we perform a MF approximation to equation (3) with translationally invariant MF Ansätze: $t_{\sigma }^{{{\alpha }_{ij}}}=\langle {{t}_{\sigma ,ij}}\rangle$, $\Delta _{\mu }^{{{\alpha }_{ij}}}=\frac{1}{4\sqrt{2}}\langle {{\Delta }_{\mu ,ij}}\rangle ,\mu =0,1,2,3$. By solving these MF parameters self-consistently, we find two local minima of the MF solutions: (1) $\Delta _{0}^{\alpha }=0.043\equiv {{\Delta }_{0}},\Delta _{l}^{\alpha }=t_{\sigma }^{\alpha }=\lambda =0$, leading to an s-wave MF Hamiltonian

$$\begin{eqnarray}&&{{H}_{s}}=-\frac{{{J}_{K}}{{\Delta }_{0}}}{2}\mathop{\sum }\limits_{\langle i,j\rangle }[(f_{i\uparrow }^{\dagger }f_{j\downarrow }^{\dagger }-f_{i\downarrow }^{\dagger }f_{j\uparrow }^{\dagger })+H.c.]\end{eqnarray} \tag{ 4 }$$

with a pseudo Fermi-surface; and (2) $\Delta _{0}^{z}={{e}^{-{\bf i}\frac{2\pi }{3}}}\Delta _{0}^{x}={{e}^{-{\bf i}\frac{4\pi }{3}}}\Delta _{0}^{y}=0.041,t_{\sigma }^{\alpha }=0.116\equiv {{t}_{0}},\Delta _{l}^{\alpha }=\lambda =0$, leading to a $d+{\rm i}d$-wave MF Hamiltonian

$$\begin{eqnarray}&&{{H}_{d+{\rm i}d}}=-\frac{{{J}_{K}}}{8}\mathop{\sum }\limits_{\langle i,j\rangle }[-{{t}_{0}}(f_{i\uparrow }^{\dagger }{{f}_{j\uparrow }}+f_{i\downarrow }^{\dagger }{{f}_{j\downarrow }})+4\Delta _{0}^{{{\alpha }_{ij}}}(f_{i\uparrow }^{\dagger }f_{j\downarrow }^{\dagger }-f_{i\downarrow }^{\dagger }f_{j\uparrow }^{\dagger })+H.c.]\end{eqnarray} \tag{ 5 }$$

with a finite bulk energy gap $0.044|{{J}_{K}}|$.

We further find that ${{H}_{d+{\rm i}d}}$ has a lower MF GS energy ($-0.144|{{J}_{K}}|$ per site) than H_s ($-0.107|{{J}_{K}}|$ per site), indicating that the GS of equation (5) is energetically favorable at the MF level. Therefore, we will focus on the GS properties of ${{H}_{d+{\rm i}d}}$ hereafter. The MF Hamiltonian (5) preserves the lattice-translation and -inversion and the spin-rotation symmetries, and thus the projected physical spin state also preserves these symmetries, which describes a SL phase.

The gauge structure of the MF Hamiltonian ${{H}_{d+{\rm i}d}}$ is described by the so-called invariant gauge group (IGG) [46]. More precisely, to calculate IGG, it is instructive to introduce the notation ${{\psi }_{i}}={{({{f}_{i\uparrow }},f_{i\downarrow }^{\dagger })}^{T}}$ and rewrite equation (5) as

$$\begin{eqnarray}&&{{H}_{d+{\rm i}d}}=\mathop{\sum }\limits_{\langle i,j\rangle }(\psi _{i}^{\dagger }{{\chi }_{ij}}{{\psi }_{j}}+H.c.),\end{eqnarray} \tag{ 6 }$$

where ${{\chi }_{ij}}=-\frac{{{J}_{K}}}{2}{\tt Re}(\Delta _{0}^{{{\alpha }_{ij}}}){{\tau }^{x}}+\frac{{{J}_{K}}}{2}{\tt Im}(\Delta _{0}^{{{\alpha }_{ij}}}){{\tau }^{y}}+\frac{{{J}_{K}}{{t}_{0}}}{8}{{\tau }^{z}}$. The IGG is now defined as the set formed by all the $SU(2)$ gauge transformations that leave ${{\chi }_{ij}}$ unchanged, i.e., IGG $=\{{{G}_{i}}|{{G}_{i}}{{\chi }_{ij}}G_{j}^{\dagger }={{\chi }_{ij}},{{G}_{i}}\in SU(2)\}$. For equation (6), we get IGG$\;=\;\{{{G}_{i}}=\pm {\bf 1}\}$ with the 2 × 2 identity matrix ${\bf 1}$, indicating that the SL state of equation (5) is in fact a Z₂ SL.

The $d+{\rm i}d$ MF Hamiltonian (5) breaks the TR symmetry by itself (i.e., TR: ${{H}_{d+{\rm i}d}}\to {{H}_{d-{\rm i}d}}$). However, from a TR breaking MF Hamiltonian alone, we can not infer the TR symmetry violation [46, 47]. This is because the Schwinger-fermion representation of the spin-1/2 operators enlarges the spin Hilbert space and introduces an $SU(2)$ gauge redundancy. As a result, the MF GS after projection $\mathcal{P}$ (where $\mathcal{P}$ removes the unphysical states containing empty or doubly occupied sites) may have higher symmetries than the original MF Hamiltonian [46].

In practice, the TR symmetry is difficult to verify by writing down the projected GS wavefunction of equation (5) directly, although this can be done in principle. Following reference [39], we introduce the $SU(2)$ gauge-invariant loop variable to diagnose whether or not the TR symmetry is broken³ . The loop variable W_l is defined on a closed, oriented loop $l={{j}_{1}}\to {{j}_{2}}\to \cdots \to {{j}_{n}}\to {{j}_{1}}$ as

$$\begin{eqnarray}&&{{W}_{l}}={\tt tr}({{\chi }_{{{j}_{1}}{{j}_{2}}}}{{\chi }_{{{j}_{2}}{{j}_{3}}}}\cdots \;{{\chi }_{{{j}_{n}}{{j}_{1}}}}),\end{eqnarray} \tag{ 7 }$$

where ${\tt tr}(\cdot )$ represents the trace of a matrix. If the W_l configuration is changed under the TR transformation, then the TR symmetry would break down. One can show that under TR

$$\begin{eqnarray}&&{{W}_{l}}\to {\tt tr}[(-{{\tau }^{y}}{{\chi }_{{{j}_{1}}{{j}_{2}}}}{{\tau }^{y}})(-{{\tau }^{y}}{{\chi }_{{{j}_{2}}{{j}_{3}}}}{{\tau }^{y}})\cdots (-{{\tau }^{y}}{{\chi }_{{{j}_{n}}{{j}_{1}}}}{{\tau }^{y}})]={{(-1)}^{n}}{{W}_{l}}.\end{eqnarray} \tag{ 8 }$$

Thus, it can be seen that for a bipartite lattice (e.g., the square or honeycomb lattice) where all loops l have even length (i.e., n are even), the W_l configuration is invariant under TR, and hence the TR symmetry should be maintained. On the contrary, for a non-bipartite lattice (e.g., the triangular or kagome lattice) containing odd loops l (i.e., loops with odd length), the corresponding loop variables W_l (if nonzero) are changed by a sign under TR and hence the W_l configuration would be changed by TR (see equation (8)), implying the TR symmetry breaking [48].

Here, for the $d+{\rm i}d$ MF Hamiltonian (5) or (6) on a non-bipartite triangular lattice, we find that the loop variable around each triangular plaquette reads

$$\begin{eqnarray}&&{{W}_{\vartriangle }}=\pm {\bf i}3\sqrt{3}J_{K}^{3}{{t}_{0}}|\Delta _{0}^{\alpha }{{|}^{2}}/32,\end{eqnarray} \tag{ 9 }$$

which is a nonzero number. In fact, these nonzero triangular loop variables reflect the intrinsic frustration of a non-bipartite triangular lattice. Thus, the Z₂ SL state described by equation (5) indeed breaks the TR symmetry. In addition, it also breaks the parity symmetry in two spatial dimensions (i.e., reflection about the axis along the z link). Therefore, our Z₂ SL state is a chiral SL state that breaks both the TR symmetry and the parity symmetry [47].

The $d+{\rm i}d$ MF Hamiltonian (5) also has nontrivial band topology as can be seen by rewriting it in momentum space (up to nonessential constant terms): ${{H}_{d+{\rm i}d}}={{\sum }_{{\boldsymbol{k}} }}(f_{{\boldsymbol{k}} \uparrow }^{\dagger },{{f}_{-{\boldsymbol{k}} \downarrow }})[{\boldsymbol{h}} ({\boldsymbol{k}} )\cdot {\boldsymbol{ \tau }} ]{{(f_{{\boldsymbol{k}} \uparrow }^{\dagger },{{f}_{-{\boldsymbol{k}} \downarrow }})}^{\dagger }}$, where ${\boldsymbol{h}} ({\boldsymbol{k}} )=({\tt Re}{{\Delta }_{{\boldsymbol{k}} }},-{{\tt Im}{\Delta }_{{\boldsymbol{k}} }},{{\varepsilon }_{{\boldsymbol{k}} }})$, ${{\varepsilon }_{{\boldsymbol{k}} }}=\frac{{{J}_{K}}{{t}_{0}}}{4}[{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})]$, ${{\Delta }_{{\boldsymbol{k}} }}=-{{J}_{K}}[\Delta _{0}^{z}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+\Delta _{0}^{x}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{0}^{y}{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})]$, and the lattice unit vectors ${{{\boldsymbol{a}} }_{1}}=(1,0)$ and ${{{\boldsymbol{a}} }_{2}}=(-1/2,\sqrt{3}/2)$ specify the bond directions z and x, respectively. We see that the 2 × 2 Hamiltonian ${\boldsymbol{h}} ({\boldsymbol{k}} )\cdot {\boldsymbol{ \tau }}$ has two no-crossing bands where each band has a Chern number ±2.

Away from half filling, we now consider the effect of finite hole-doping by adding the NN hopping terms to the undoped triangular KH model (1), i.e.,

$$\begin{eqnarray}&&{{H}_{h}}=-t\mathop{\sum }\limits_{\langle i,j\rangle ,\sigma }c_{i\sigma }^{\dagger }{{c}_{j\sigma }}-{{J}_{K}}\mathop{\sum }\limits_{\langle i,j\rangle }S_{i}^{{{\alpha }_{ij}}}S_{j}^{{{\alpha }_{ij}}}+{{J}_{H}}\mathop{\sum }\limits_{\langle i,j\rangle }({{{\boldsymbol{S}} }_{i}}\cdot {{{\boldsymbol{S}} }_{j}}-\frac{1}{4}{{\hat{n}}_{i}}{{\hat{n}}_{j}})-\mu \mathop{\sum }\limits_{i}{{\hat{n}}_{i}},\end{eqnarray} \tag{ 10 }$$

where the chemical potential μ is adjusted such that $\langle {{\hat{n}}_{i}}\rangle =1-\delta$ with the doping level δ per site. As in the t-J model for high-T_c cuprates, here, we consider the case that the double occupancy is prohibited due to the strong on-site repulsive interactions and adopt the slave-boson approach $c_{i\sigma }^{\dagger }=f_{i\sigma }^{\dagger }{{b}_{i}}$, with additional bosonic holons b_i that are assumed to be condensed, i.e., ${{b}_{i}}\approx b_{i}^{\dagger }\approx \sqrt{\langle b_{i}^{\dagger }{{b}_{i}}\rangle }=\sqrt{\delta }$. Thus, the quadratic terms in equation (10) are reduced to ${{H}_{T}}=-\delta t{{\sum }_{\langle i,j\rangle ,\sigma }}f_{i\sigma }^{\dagger }{{f}_{j\sigma }}-\mu {{\sum }_{i,\sigma }}f_{i\sigma }^{\dagger }{{f}_{i\sigma }}+\mu (1-\delta )N$. The spin-exchange terms in equation (10) can now be decoupled into both hopping and pairing channels, as in equation (3). However, because of the presence of the kinetic term H_T, the effect of decoupling spin interactions into hopping channels is not expected to qualitatively change the SC phase diagram at reasonably large doping [32]. For simplicity, we thus consider only the pairing channels at finite doping [32], say $S_{i}^{l}S_{j}^{l}=-\frac{1}{4}(\Delta _{0,ij}^{\dagger }{{\Delta }_{0,ij}}+2\Delta _{l,ij}^{\dagger }{{\Delta }_{l,ij}}-\sum _{m=1}^{3}\Delta _{m,ij}^{\dagger }{{\Delta }_{m,ij}})$, where the corresponding MF parameters are defined as in the undoped case, $l=1,2,3$ respectively corresponds to the spin component $x,y,z$, and the summation over l gives a Heisenberg term ${{{\boldsymbol{S}} }_{i}}\cdot {{{\boldsymbol{S}} }_{j}}-\frac{1}{4}{{\hat{n}}_{i}}{{\hat{n}}_{j}}=-\Delta _{0,ij}^{\dagger }{{\Delta }_{0,ij}}$. This slave-boson MF approach has also been widely used to study the doped KH model on the honeycomb lattice very recently [32–35], and subsequent unbiased numerical methods [36] have confirmed those MF results at the qualitative level.

Here, we focus on the following MF Ansätze: (1) p-wave pairing, $\Delta _{l}^{x}=\Delta _{l}^{y}=\Delta _{l}^{z}\equiv {{\Delta }_{l}}$; (2) s-wave pairing, $\Delta _{0}^{x}=\Delta _{0}^{y}=\Delta _{0}^{z}$; and (3) ${{d}_{{{x}^{2}}-{{y}^{2}}}}+{\rm i}{{d}_{xy}}$-wave pairing, ${{e}^{-{\bf i}\frac{2\pi }{3}}}\Delta _{0}^{x}={{e}^{-{\bf i}\frac{4\pi }{3}}}\Delta _{0}^{y}=\Delta _{0}^{z}$. Based on these MF Ansätze, the resulting phase diagram as a function of δ and ${{J}_{H}}/{{J}_{K}}$ is shown in figure 4, which is qualitatively analogous to the SC phase diagram of the honeycomb KH model at finite doping [32, 35, 36].

Zoom In Zoom Out Reset image size

We find that (1) the p-wave solutions are symmetric ${{\Delta }_{1}}={{\Delta }_{2}}={{\Delta }_{3}}$ and pure imaginary, and (2) the s-wave solutions are real. Thus, both the p-wave SC and the s-wave SC are TR invariant. We further note that the p-wave SC is fully gapped in the doping interval $0.05\leqslant \delta \leqslant 0.4$ (see the inset of figure 4(a)), indicating the absence of a topological phase transition in the finite-doping regime. In fact, the topological property of a fully gapped spin-triplet SC with TR symmetry is intimately related to the Fermi-surface topology in the normal state [48, 49]. More precisely, the Z₂ invariant is determined by the parity of the number of TR invariant points⁴ below the Fermi level. Here, the Fermi level is determined by the energy dispersion ${{\epsilon }_{{\boldsymbol{k}} }}$ of H_T, where ${{\epsilon }_{-{\boldsymbol{k}} }}={{\epsilon }_{{\boldsymbol{k}} }}$ due to the lattice inversion and TR symmetries of H_T. For the p-wave SC in the doping interval $0.05\leqslant \delta \leqslant 0.4$, we find that there is only one TR invariant point Γ below the Fermi level. According to the preceding criterion, the value of the Z₂ invariant is odd and the p-wave SC is thus topologically nontrivial [48, 49].

The preceding results focus on the effects of hole-doping. In the case of electron doping, if one wants to study the effects still in the context of t-J-type model like equation (10), technically one needs to perform a particle-hole transformation ${{c}_{i,\sigma }}\to c_{i,-\sigma }^{\dagger }$ such that the constraint now means no double occupancy of holes under electron doping (or equivalently no empty sites of electrons). The resulting electron-doped Hamiltonian then takes the form

$$\begin{eqnarray}&&{{H}_{e}}=t\mathop{\sum }\limits_{\langle i,j\rangle ,\sigma }c_{i\sigma }^{\dagger }{{c}_{j\sigma }}-{{J}_{K}}\mathop{\sum }\limits_{\langle i,j\rangle }S_{i}^{{{\alpha }_{ij}}}S_{j}^{{{\alpha }_{ij}}}+{{J}_{H}}\mathop{\sum }\limits_{\langle i,j\rangle }({{{\boldsymbol{S}} }_{i}}\cdot {{{\boldsymbol{S}} }_{j}}-\frac{1}{4}{{\hat{n}}_{i}}{{\hat{n}}_{j}})-\mu \mathop{\sum }\limits_{i}{{\hat{n}}_{i}}.\end{eqnarray} \tag{ 11 }$$

Note that as compared with equation (10), the minus sign before the hopping t disappears and the parameter t is positive here. The operator $c_{i\sigma }^{\dagger }$ in equation (11) now creates a hole instead of an electron, and the corresponding number operator ${{\hat{n}}_{i}}$ now denotes the number of holes. And the chemical potential μ is adjusted such that $\langle {{\hat{n}}_{i}}\rangle =1-\delta$ with the electron doping level $\delta \gt 0$ from half filling. Thus, the subsequent mathematical procedure to solve the electron-doped Hamiltonian (11) is the same as that for the hole-doped Hamiltonian (10) by using the slave-boson formalism. Now the p-wave SC, the s-wave SC, and the ${{d}_{{{x}^{2}}-{{y}^{2}}}}+{\rm i}{{d}_{xy}}$-wave SC MF Ansätze under electron doping are defined as the same form as those under hole-doping. Based on those MF Ansätze, a self-consistent calculation yields the SC phase diagram of electron doping, as shown in figure 5. We find that both the p-wave SC and the s-wave SC are TR invariant, and the p-wave SC is fully gapped in the electron doping interval $0.05\leqslant \delta \leqslant 0.4$. We further note that there are three TR invariant points (say ${{M}_{1}},{{M}_{2}},$ and M₃ as shown in figure 1(b)) below the Fermi level for the p-wave SC, indicating that it is also topologically nontrivial [48, 49] under electron doping.

Zoom In Zoom Out Reset image size

It is noteworthy that the electron-doped phase diagram (figure 5) is qualitatively similar to but quantitatively different from the hole-doped phase diagram (figure 4(a)). As shown, there also exist three SC phases that are the same as those under hole-doping but the region of the triplet SC phase is strongly suppressed and the region of the singlet SC phase is considerably enlarged. This asymmetry between the electron-doped and the hole-doped phase diagrams is due to the absence of particle-hole symmetry of model (10) on the triangular lattice.

Based on the Schwinger-fermion MF theory, we show that the possible candidate phase near the AF Kitaev point for the undoped triangular KH model is a gapped Z₂ chiral SL. We note that [37] suggests a nematic order near the AF Kitaev point. Since the investigation in [37] is based on a classical triangular KH model, this difference may be due to the different treatments of the spin dynamics, and whether the quantum fluctuations would distort the classical nematic order requires further study. On the other hand, the 120° Neel order has been shown to be destabilized by an infinitesimal Kitaev term in [37] with the Monte Carlo simulations. The resulting phase consists of non-coplanar spin patterns remaining locally close to the 120° pattern but distorted at larger distances, which hosts Z₂ vortices. For the local feature of the spin correlations, our result is consistent with that in [37]. Besides, the magnetically ordered phases such as the collinear stripe and the noncollinear spiral patterns obtained in our work are consistent with those in [37].

It is also instructive to compare our results for the undoped triangular KH model with the honeycomb Kitaev model that has an exact solution. The broken time-reversal symmetry of our Z₂ chiral SL stems from the existence of triangular loops △ with odd length and the corresponding nonzero loop variables equation (9). A similar reason also leads to the emergence of an exact chiral SL state in the Kitaev model defined on a decorated honeycomb lattice that also contains odd triangular loops [50]. By contrast, the exact Kitaev SL in a bipartite honeycomb lattice is TR invariant since all loops there have even length and hence the corresponding loop variables (which are called Wilson loops in [1]) are invariant under the TR transformation [48]. It is also known that the Kitaev SL in the honeycomb lattice is insensitive to the signs of coupling constants, which is again due to the bipartite nature of the honeycomb lattice. Reversing those signs is simply equivalent to a gauge transformation [1], which will not affect the physical SL state. This is in sharp contrast to the Kitaev model on a non-bipartite triangular lattice, where we find that the FM Kitaev point (i.e., the point at $\phi =\pi /2$) itself is a critical point separating the stripe and the FM phases while the AF Kitaev point represents the exotic phase which is suggested as a Z₂ chiral SL.

In conclusion, we investigate the GS phase diagrams of the undoped and doped KH model on a triangular lattice. In the undoped case, a numerical ED calculation shows five transition points separating five distinct phases in the global phase diagram, where two phases are known to be the FM order and the 120° Neel order. The FST sheds further light on the nature of the phase diagram and reveals other two magnetic orders, namely, a collinear stripe pattern and a noncollinear spiral pattern. Lastly, based on the Schwinger-fermion MF theory, an energetically favorable Z₂ chiral SL with a Chern number ±2 is proposed as a potential candidate for the exotic phase near the AF Kitaev point. Yet, finite hole-doping induces an s-wave, a ${{d}_{{{x}^{2}}-{{y}^{2}}}}+{\rm i}{{d}_{xy}}$-wave, and a TR invariant topological p-wave superconducting states.

Finally, a few remarks are in order concerning the possible experimental realization of our results and the relevant materials. The recently synthesized iridate Ba₃IrTi₂O₉ [51] consisting of a layered triangular arrangement of Ir⁴⁺ ions with an effective magnetic moment ${{J}_{\tt eff}}=\frac{1}{2}$ possesses almost all the necessary ingredients for Kitaev-type exchange couplings (as mentioned in the Introduction.), and thus it is a promising candidate for a microscopic realization of the triangular KH model. Experimental measurements on the Ba₃IrTi₂O₉ show no magnetic ordering down to 0.35 K and suggest that Ba₃IrTi₂O₉ probably hosts a SL GS [51], and hence the proposed Z₂ chiral SL phase may be observed in this material. In addition, a class of ABO₂-type (where A and B are alkali and transition metal ions, respectively) transition metal compounds [2, 3] and the possible material Na_xIrO₂ [38] may be potential candidates for realizing the triangular KH model. It was also found very recently that both Pt or Pd intercalations and substitutions of layered iridium ditelluride IrTe₂ could induce bulk SC with T_c up to ∼3 K [52], and the SOC in these compounds is expected to be strong due to the large atomic numbers of both Ir and Te [52, 53]. Thus, a layered IrTe₂ with Pt or Pd intercalations and substitutions may be a candidate for the unconventional SC proposed here.

Note added. Upon the completion of this manuscript, we became aware of an independent work [54] containing related results on the undoped case. The four magnetically ordered phases obtained here are consistent with those in [54]. However, [54] starts with a classical GS at the AF Kitaev point and further includes the effects of quantum fluctuations via a numerical analysis; consequently, it suggests a nematic phase. Therefore, the issue related to this difference is whether the quantum fluctuations are strong enough to spoil the classical nematic order into a disordered phase (e.g., a quantum SL).

This work was supported by the National Natural Science Foundation of China (91021001, 11190023, and 11204125) and the Ministry of Science and Technology of China (973 Project Grants No. 2011CB922101 and No. 2011CB605902).

In this appendix, we provide some details of the MF calculations in the momentum space. After the MF decoupling, in both the undoped and doped cases, the Fourier transformed MF Hamiltonians with translationally invariant Ansätze take the following form

$$\begin{eqnarray}\begin{array}{rcl} {{H}_{MF}} & = & \mathop{\sum }\limits_{{\boldsymbol{k}} ,\sigma }{{\varepsilon }_{{\boldsymbol{k}} \sigma }}f_{{\boldsymbol{k}} \sigma }^{\dagger }{{f}_{{\boldsymbol{k}} \sigma }}+\mathop{\sum }\limits_{{\boldsymbol{k}} }\left( {{\Delta }_{{\boldsymbol{k}} \uparrow }}f_{{\boldsymbol{k}} \uparrow }^{\dagger }f_{-{\boldsymbol{k}} \uparrow }^{\dagger }+{{\Delta }_{{\boldsymbol{k}} \downarrow }}f_{{\boldsymbol{k}} \downarrow }^{\dagger }f_{-{\boldsymbol{k}} \downarrow }^{\dagger } \right. \\ {} & {} & +\left. \Delta _{{\boldsymbol{k}} }^{t}f_{{\boldsymbol{k}} \uparrow }^{\dagger }f_{-{\boldsymbol{k}} \downarrow }^{\dagger }+\Delta _{{\boldsymbol{k}} }^{s}f_{{\boldsymbol{k}} \uparrow }^{\dagger }f_{-{\boldsymbol{k}} \downarrow }^{\dagger }+H.c. \right)+{{E}_{0}} \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

with the triplet pairing functions $({{\Delta }_{{\boldsymbol{k}} \uparrow }},{{\Delta }_{{\boldsymbol{k}} \downarrow }},\Delta _{{\boldsymbol{k}} }^{t})$, the singlet pairing function $\Delta _{{\boldsymbol{k}} }^{s}$, and a constant term E₀. Specifically, in the undoped case, we have

$$\begin{eqnarray} {{\varepsilon }_{{\boldsymbol{k}} \sigma }} & = & \frac{{{J}_{K}}}{4}[t_{\sigma }^{z}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+t_{\sigma }^{x}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+t_{\sigma }^{y}{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})]-\lambda , \\ {{\Delta }_{{\boldsymbol{k}} \uparrow }} & = & -{\bf i}{{J}_{K}}[\Delta _{1}^{z}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}-\Delta _{1}^{x}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}-\Delta _{1}^{y}{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})] \\ {} & {} & +{{J}_{K}}[\Delta _{2}^{z}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+\Delta _{2}^{x}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{2}^{y}{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ {{\Delta }_{{\boldsymbol{k}} \downarrow }} & = & {\bf i}{{J}_{K}}[\Delta _{1}^{z}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}-\Delta _{1}^{x}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}-\Delta _{1}^{y}{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})] \\ {} & {} & +{{J}_{K}}[\Delta _{2}^{z}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+\Delta _{2}^{x}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{2}^{y}{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \Delta _{{\boldsymbol{k}} }^{t} & = & -{\bf i}{{J}_{K}}[\Delta _{3}^{z}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}-\Delta _{3}^{x}{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{3}^{y}{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \Delta _{{\boldsymbol{k}} }^{s} & = & -{{J}_{K}}[\Delta _{0}^{z}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+\Delta _{0}^{x}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{0}^{y}{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \Delta _{{\boldsymbol{k}} }^{s} & = & -{{J}_{K}}[\Delta _{0}^{z}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+\Delta _{0}^{x}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{0}^{y}{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \end{eqnarray}$$

and the constant term

$$\begin{eqnarray*} \begin{array}{rcl} {{E}_{0}} & = & -\frac{N{{J}_{K}}}{8}(t_{\uparrow }^{x*}t_{\downarrow }^{x}+t_{\uparrow }^{y*}t_{\downarrow }^{y}+H.c.+|t_{\uparrow }^{z}{{|}^{2}}+|t_{\downarrow }^{z}{{|}^{2}})-4N{{J}_{K}}\mathop{\sum }\limits_{\alpha }|\Delta _{0}^{\alpha }{{|}^{2}} \\ {} & {} & +8N{{J}_{K}}\mathop{\sum }\limits_{\alpha }|\Delta _{1}^{\alpha }{{|}^{2}}+8N{{J}_{K}}\mathop{\sum }\limits_{\alpha }|\Delta _{2}^{\alpha }{{|}^{2}}+4N{{J}_{K}}\mathop{\sum }\limits_{\alpha }|\Delta _{3}^{\alpha }{{|}^{2}} \\ {} & {} & -16N{{J}_{K}}|\Delta _{1}^{x}{{|}^{2}}-16N{{J}_{K}}|\Delta _{2}^{y}{{|}^{2}}-8N{{J}_{K}}|\Delta _{3}^{z}{{|}^{2}}+N\lambda . \\ \end{array} \end{eqnarray*}$$

While in the finite doped case, we have

$$\begin{eqnarray*}\begin{array}{rcl} {{\varepsilon }_{{\boldsymbol{k}} \sigma }} & = & -2\delta t[{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})]-\mu , \\ {{\Delta }_{{\boldsymbol{k}} \uparrow }} & = & -{\bf i}{{J}_{K}}{{\Delta }_{1}}[{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}-{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}-{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})] \\ {} & {} & +{{J}_{K}}{{\Delta }_{2}}[{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ {{\Delta }_{{\boldsymbol{k}} \downarrow }} & = & {\bf i}{{J}_{K}}{{\Delta }_{1}}[{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}-{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}-{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})] \\ {} & {} & +{{J}_{K}}{{\Delta }_{2}}[{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \Delta _{{\boldsymbol{k}} }^{t} & = & -2{\bf i}{{J}_{K}}{{\Delta }_{3}}[{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}-{\rm sin} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+{\rm sin} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \Delta _{{\boldsymbol{k}} }^{s} & = & 2(4{{J}_{H}}-{{J}_{K}})[\Delta _{0}^{z}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{1}}+\Delta _{0}^{x}{\rm cos} {\boldsymbol{k}} \cdot {{{\boldsymbol{a}} }_{2}}+\Delta _{0}^{y}{\rm cos} {\boldsymbol{k}} \cdot ({{{\boldsymbol{a}} }_{1}}+{{{\boldsymbol{a}} }_{2}})], \\ \end{array}\end{eqnarray*}$$

and the constant term

$$\begin{eqnarray*}&&{{E}_{0}}=8N(4{{J}_{H}}-{{J}_{K}})\mathop{\sum }\limits_{\alpha }|\Delta _{0}^{\alpha }{{|}^{2}}+8N{{J}_{K}}\mathop{\sum }\limits_{l}|{{\Delta }_{l}}{{|}^{2}}+N\mu (1-\delta ).\end{eqnarray*}$$

The quadratic Hamiltonian (A.1) can be diagonalized via the Bogoliubov transformation and the resulting fermion spectrum is given by

$$\begin{eqnarray*}&&{{E}_{\pm }}({\boldsymbol{k}} )=\frac{1}{2}|{{\varepsilon }_{{\boldsymbol{k}} \uparrow }}-{{\varepsilon }_{{\boldsymbol{k}} \downarrow }}|+\sqrt{\frac{1}{4}{{({{\varepsilon }_{{\boldsymbol{k}} \uparrow }}+{{\varepsilon }_{{\boldsymbol{k}} \downarrow }})}^{2}}+\frac{1}{2}({{A}_{{\boldsymbol{k}} }}\pm {{B}_{{\boldsymbol{k}} }})}\end{eqnarray*}$$

with ${{A}_{{\boldsymbol{k}} }}$ $\;=\;4{{\sum }_{\sigma }}|{{\Delta }_{{\boldsymbol{k}} \sigma }}{{|}^{2}}+|{{\Delta }_{{\boldsymbol{k}} }}{{|}^{2}}+|{{\Delta }_{-{\boldsymbol{k}} }}{{|}^{2}}$, ${{B}_{{\boldsymbol{k}} }}$ $\;=\;\left[ A_{{\boldsymbol{k}} }^{2}-64 \right.$ $|{{\Delta }_{{\boldsymbol{k}} \uparrow }}{{|}^{2}}|{{\Delta }_{{\boldsymbol{k}} \downarrow }}{{|}^{2}}-4|{{\Delta }_{{\boldsymbol{k}} }}{{|}^{2}}|{{\Delta }_{-{\boldsymbol{k}} }}{{|}^{2}}$ ${{\left. -32{\tt Re}(\Delta _{{\boldsymbol{k}} \uparrow }^{*}\Delta _{{\boldsymbol{k}} \downarrow }^{*}{{\Delta }_{{\boldsymbol{k}} }}{{\Delta }_{-{\boldsymbol{k}} }}) \right]}^{\frac{1}{2}}}$,and ${{\Delta }_{{\boldsymbol{k}} }}=\Delta _{{\boldsymbol{k}} }^{t}+\Delta _{{\boldsymbol{k}} }^{s}$.

And the GS energy of equation (A.1) is given by

$$\begin{eqnarray*}&&{{E}_{g}}=\frac{1}{2}\mathop{\sum }\limits_{{\boldsymbol{k}} ,\sigma }{{\varepsilon }_{{\boldsymbol{k}} \sigma }}-\frac{1}{2}\mathop{\sum }\limits_{{\boldsymbol{k}} }[{{E}_{+}}({\boldsymbol{k}} )+{{E}_{-}}({\boldsymbol{k}} )]+{{E}_{0}}.\end{eqnarray*}$$