††thanks<EMAIL_ADDRESS>
# Phase diagram for Hole-Doped Kitaev System on the Honeycomb Lattice

Su-Ming Zhang  Renmin University of China, Haidian District, Beijing, China,
100872  Zheng-Xin Liu  Renmin University of China, Haidian District, Beijing,
China, 100872

###### Abstract

In extended Kitaev models on the honeycomb lattice, off-diagonal interactions
(e.g. the $\Gamma,\Gamma^{{}^{\prime}}$ terms) can give rise to non-Kitaev
quantum spin liquids and several magnetically ordered phases. In the present
work, we dope holes to the system and study the resultant
$t$-$K$-$\Gamma$-$\Gamma^{{}^{\prime}}$ model using mean field theory. The
interplay between the charge and spin degrees of freedom results in a rich
phase diagram. Similar to doped cuprates, superconductors, pseudogap phases,
fermi liquid, strange metal and paramagnetic phase are generated. What is
different is that, we obtain more than one superconducting phase (including a
topological one) and more than one pseudogap phase no matter what the original
spin state is. The Chern number of the topological superconductor is either
$\nu=\pm 2$ or $\nu=\pm 1$, depending on the ratio $\Gamma/|K|$ in the spin
channel. We further find that an intermediate in-plane magnetic field can
slightly enlarge the size of the topological superconducting phase.

##  I INTRODUCTION

Quantum spin liquids (QSLs) are exotic phases of matter which exhibit long-
range entanglements but no traditional symmetry-breaking orders even at
extremely low temperatures.Balents (2010); Kitaev (2006) The strong quantum
fluctuations in systems with geometric frustrations or frustrating
interactions causes the melting of long-range magnetic orders, resulting in
disordered quantum state without breaking of spin and lattice symmetries. The
elementary excitations in a gapped QSL obey fractional statistical and have
potential applications in quantum computationsZhou et al. (2017); Kitaev
(2006). Therefore the realization of QSLs has attracted lots of interest in
condensed matter physics. Theoretical studies indicate that QSLs generally
result from the competition of different magnetic ordersChandra and Doucot
(1988); He et al. (2014); Yan et al. (2011). Experimentally, candidate QSL
materials have been reported on various lattice structures, including
triangular lattice, kagome lattice and even three dimensional magnetsShimizu
et al. (2003); Okamoto et al. (2007); Han et al. (2012); Chillal et al.
(2020).

QSLs also exist in magnets with strong spin-orbit couplings. In 2006, Kitaev
proposed a honeycomb lattice model which has an exactly solvable QSL ground
stateKitaev (2006). Later, candidate materials containing the Kitaev type
interactions, called the Kitaev materials, are proposed, including the well
studied $\alpha$-RuCl3 and Na2IrO3Jackeli and Khaliullin (2009); Chaloupka et
al. (2010). However, most of these materials exhibit zigzag type long-range
order at low temperatures, indicating the existence of non-Kitaev interactions
such as the Heisenberg exchangesJanssen et al. (2017); Winter et al. (2017);
Kimchi and You (2011); Winter et al. (2016); Singh et al. (2012) and off-
diagonal $\Gamma,\Gamma^{\prime}$ termsWinter et al. (2016); Rau and Kee
(2014); Chaloupka and Khaliullin (2015); Gordon et al. (2019); Lee et al.
(2020). Interestingly, external forces, such as magnetic fieldsBaek et al.
(2017); Zheng et al. (2017); Yadav et al. (2016) or mechanical
pressureMirebeau et al. (2002); Kozlenko et al. (2008); Shimizu et al. (2003);
Wang et al. (2018); Cui et al. (2017), can suppress the order and drive the
system into a spin-liquid like or a dimer-like disordered state.

An intuitional picture for a QSL with continuous spin rotation symmetry is the
Resonating Valence Bond(RVB) state raised by Anderson Anderson (1987). In a
RVB picture, the spins form singlet pairs and the ground state is a
superposition of all possible pairing configurations, resultantly the lattice
rotation and translation symmetries are restored. It was also proposed that
when doping charge carriers into the system, the previously paired electrons
behave like Cooper pairs and move coherently. Therefore, the QSL state is
driven into a superconductor Anderson (1987). This provides a scheme to
explain the mechanism of high-temperature superconductivity in cuprate
oxidesLee and Nagaosa (1992). In cuprates, the Neel order is easily destroyed
by doped charge carriers. When the doping concentration falls in the optimal
region, the cuprate becomes a superconductor with a high critical
temperatureBednorz and Müller (1986); Wu et al. (1987); Jiang et al. (2018);
Ruan et al. (2016); Mei et al. (2012); Weng (2011).

An interesting question is what happens if charge carriers are doped into the
Kitaev materials which have anisotropic spin-spin interactions with strong
spin-orbit couplings. Anderson’s theory indicates that superconductors should
appear. This was verified in an earlier study of the $t$-$K$-$J$ modelYou et
al. (2012). In the present work, we dope holes into Kitaev systems with off-
diagonal $\Gamma$ and $\Gamma^{\prime}$ exchange interactions, and study the
resultant $t$-$K$-$\Gamma$-$\Gamma^{\prime}$ model using mean field theory.
The phase diagram is very rich. No matter what the original spin state is,
more than one superconductor and more than one pseudogap phase emerge at
different doping concentration. Similar to cuprates, Fermi liquid phase or the
strange metal phase appear at heavy doping. Interestingly, several topological
superconductors are obtained, whose Chern number ($\nu=\pm 2,\pm 1$) are
dependent on the ratio $\Gamma/|K|$ in the spin-spin interactions. We further
find that intermediate in-plane magnetic fields can enlarge the range of the
SC phases.

The rest part of the present work is organized as follows. In section II, we
introduce the effective lattice model (1) as the starting point. In section
III, we introduce the slave particle representation and the related mean field
theory in details, we emphasize the existence of internal gauge symmetry in
the slave particle representation and specially introduce the symmetry group
of the mean-field description of QSLs — the projective symmetry groups. In
section IV, we present the mean-field phase diagrams at zero doping and finite
doping, respectively. Section V is devoted to the conclusions and discussions.

##  II The Model

In most Kitaev materials, the magnetic cations interact with each other via
super-exchanges mediated through nonmagnetic anions. Therefore, both the $d$
(or $f$) electrons in the cations and the $p$ electrons in the anions are
involved in the electronic structures, thus more than one electronic band
enters the low-energy physics. However, we still adopt the single band
assumption to describe the low energy properties of the system, just like in
the cuprate oxides.

To study the effect of charge carriers, we introduce the electron hopping term
(namely the $t$ term) in the $K$-$\Gamma$-$\Gamma^{\prime}$ effective spin
model, in analogy to the $t$-$J$ model for cuprates.

The Heisenberg exchange interactions are believed to be small and are hence
neglected in our discussion. Hence, the starting point of our study is the
following effective model,

$\displaystyle H=$ $\displaystyle-t\sum_{\sigma\langle
ij\rangle}(c^{\dagger}_{i\sigma}c_{j\sigma}+{\rm
H.c.})-\mu\sum_{i}c^{\dagger}_{i\sigma}c_{i\sigma}$
$\displaystyle+\sum_{\langle
ij\rangle\in\alpha\beta(\gamma)}\Big{[}KS^{\gamma}_{i}S^{\gamma}_{j}+\Gamma(S_{i}^{\alpha}S_{j}^{\beta}+S_{i}^{\beta}S_{j}^{\alpha})$
$\displaystyle+\Gamma^{{}^{\prime}}(S_{i}^{\alpha}S_{j}^{\gamma}+S_{i}^{\gamma}S_{j}^{\alpha}+S_{i}^{\beta}S_{j}^{\gamma}+S_{i}^{\gamma}S_{j}^{\beta})\Big{]},$
(1)

where $t$ is the electron hopping constant, $\mu$ is the chemical potential of
the electrons used to control the density of the doped holes, $K$ stands for
the Kitaev interaction, $\Gamma$ and $\Gamma^{\prime}$ are off-diagonal spin-
spin interactions(where $\alpha,\beta,\gamma=x,y,z$). For a singly occupied
site $i$, the spin operators $S_{i}^{\alpha}$ are related to the electron
operators as $S_{i}^{\alpha}=C_{i}^{\dagger}{\sigma_{\alpha}\over 2}C_{i}$
with $C_{i}=(c_{i\uparrow},c_{i\downarrow})$. In later discussions, we assume
$K<0,\Gamma>0$ Janssen et al. (2017); Liu and Normand (2018); Banerjee et al.
(2017); Wang et al. (2017); Cookmeyer and Moore (2018), and fix
$\Gamma^{\prime}=-0.02|K|$ (the $\Gamma^{\prime}$ interaction is introduced to
stabilize the zigzag order observed in Kitaev materialsWinter et al. (2016);
Rau and Kee (2014)).

##  III Mean Field Theory

###  III.1 Slave Particle Representation

In order to study the above Hamiltonian (1), we introduce the slave particle
(called the slave boson) representation to represent the electron creation
operators in forms of fermionic spinons $f$ and bosonic holons $b^{\dagger}$,
namely

$\displaystyle c_{i\sigma}=f_{i\sigma}b_{i}^{\dagger},$

where double occupancy is forbidden (otherwise another species of boson called
doublon need to be introduced). The number of the holons
$\sum_{i}b_{i}^{\dagger}b_{i}$ equals to the number of doped holes (i.e. the
number of removed electrons in the magnetic layer), and the ratio
$\delta={1\over N}\sum_{i}b_{i}^{\dagger}b_{i}$ is the doping concentration
relative to half-filling, where $N$ is the total number of sites.

The spin operators only involve the spinons operators, and can be written as
the familiar quadratic form

$S^{\alpha}_{i}=\frac{1}{2}f^{\dagger}_{i}\sigma_{\alpha}f_{i},$

where $f^{\dagger}_{i}=(f^{\dagger}_{i\uparrow},f^{\dagger}_{i\downarrow})$
and $\sigma_{\alpha},\alpha=x,y,z$ are the Pauli matrix. For convenience, we
introduce a matrix operator
$\psi_{i}=\begin{pmatrix}f_{i\uparrow}&f^{\dagger}_{i\downarrow}\\\
f_{i\downarrow}&-f^{\dagger}_{i\uparrow}\end{pmatrix}$, and write the spin
operators as $S^{\alpha}_{i}={\rm
Tr}(\psi^{\dagger}_{i}\frac{\sigma_{\alpha}}{4}\psi_{i})$.

The condition of no double occupancy requires that an on-site particle number
constraint should be imposed,

$\displaystyle\sum_{i\sigma}f_{i\sigma}^{\dagger}f_{i\sigma}+b_{i}^{\dagger}b_{i}=1.$
(2)

Obviously, under the phase transformation

$\displaystyle f_{i\sigma}\to f_{i\sigma}e^{i\theta_{i}},\ \ \ \
b_{i}^{\dagger}\to b_{i}^{\dagger}e^{-i\theta_{i}},$

the electron operator $c_{i\sigma}$ remains unchanged $c_{i\sigma}\to
c_{i\sigma}$. This reveals an emergent local $U(1)$ gauge symmetry in the
slave particle representation of electrons (at half filling the gauge symmetry
is enlarged into $SU(2)$ Lee et al. (2006); Affleck et al. (1988); Wen and Lee
(1996); Lee et al. (1998)). The spinons and the holons interact with each
other via the emergent $U(1)$ gauge field. Depending on the electron-electron
interactions and the doping concentration, the $U(1)$ gauge field may get
confined, resulting in a magnetically ordered phase or a superconducting
phaseLee and Nagaosa (1992). To begin with, we need to translate the
interactions between the electrons into the interactions between the slave
particles.

In the slave-particle representation, the electronic hopping term on a single
bond [see (1)] can be written in terms of $f_{i}$ and $b_{i}$ as

$\displaystyle-t\sum_{\langle
ij\rangle\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}=-t\sum_{\langle
ij\rangle\sigma}f_{i\sigma}^{\dagger}f_{j\sigma}b_{j}^{\dagger}b_{i}.$ (3)

Similarly, the spin-spin interactions can be be written in biquadratic form of
spionons. For instance, the Kitaev interaction can be written as

$\displaystyle S^{m}_{i}S^{m}_{j}=-\frac{1}{16}[{\rm
Tr}(\psi^{\dagger}_{j}\psi_{i}\psi^{\dagger}_{i}\psi_{j})+{\rm
Tr}(\psi^{\dagger}_{j}\sigma^{m}\psi_{i}\psi^{\dagger}_{i}\sigma^{m}\psi_{j})].$

The fermionic representation of the $\Gamma,\Gamma^{\prime}$ interactions are
given in Appendix A.

An advantage of the slave particle representation is that the various
interactions can be treated using mean-field approximation.

This approach was applied to uncover the mechanism of high-temperature
superconductivity in cupratesLee and Nagaosa (1992); Lee et al. (2006). In the
rest part of this work, we will study the doped Kitaev systems using mean
field theory.

###  III.2  QSLs at Half-filling

In mean field theory, the interactions are decoupled into non-interacting
Hamiltonians which can be easily solved. Before investigating the effect of
hole doping, we firstly introduce the mean field theory of the pure spin-spin
interaction,

$\displaystyle H=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\langle
ij\rangle\in\alpha\beta(\gamma)}KS^{\gamma}_{i}S^{\gamma}_{j}+\Gamma(S_{i}^{\alpha}S_{j}^{\beta}+S_{i}^{\beta}S_{j}^{\alpha})$
(4)
$\displaystyle+\Gamma^{{}^{\prime}}(S_{i}^{\alpha}S_{j}^{\gamma}+S_{i}^{\gamma}S_{j}^{\alpha}+S_{i}^{\beta}S_{j}^{\gamma}+S_{i}^{\gamma}S_{j}^{\beta}).$

Above Hamiltonian has the $D_{3d}\times Z_{2}^{T}$ group symmetry Wang et al.
(2019) in the sense that each element is a combination of lattice operation
and the corresponding spin operation. Owing to spin-orbit coupling, the spin
axes are locked with the lattice axes. For instance, the $c$-direction of the
lattice plane is parallel to ${1\over\sqrt{3}}(\hat{x}+\hat{y}+\hat{z})$
direction in the spin-frame, and the three lattice-bond directions are
parallel to ${1\over\sqrt{2}}(\hat{x}-\hat{y})$,
${1\over\sqrt{2}}(\hat{y}-\hat{z})$ and ${1\over\sqrt{2}}(\hat{z}-\hat{x})$,
respectively.

The model (4) have been studied in literatureWang et al. (2020), and it was
known that there are several magnetically ordered phases in the ground state
phase diagram. In mean field theory, the magnetically ordered states are
essentially the same as the classical ground states whose energy can be easily
calculatedRau et al. (2014). Therefore, we will focus on possible QSL phases.
Noticing that the pure Kitaev model is a special point in the model (4), it is
expected that the phase diagram contains at least one QSL phase, the Kitaev
Spin Liquid (KSL) phase. In the following we will investigate if there exist
other QSL phases.

Since the model (4) only contains nearest neighbor interactions, naturally
only nearest neighbor coupling terms need to be considered in the mean field
Hamiltonian. Suppose that the QSL contains zero fluxes in each hexagon (it has
been shown that the flux-states are not energetically favoredWang et al.
(2019)), and that $C_{3}$ symmetry is preserved in the mean field Hamiltonian.
The most general mean field Hamiltonian takes the following form,

$\displaystyle H_{\rm mf}^{\rm SL}\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\langle
ij\rangle\in\alpha\beta(\gamma)}\\!\\!\Big{\\{}[f_{i}^{\dagger}(t_{1}^{\gamma}R_{\alpha\beta}+t_{0}^{\gamma}+t_{2}^{\gamma}\sigma^{\gamma}+t_{3}^{\gamma}\sigma^{\gamma}R_{\alpha\beta})f_{j}$
(5)
$\displaystyle\\!\\!\\!+f_{i}^{\dagger}(\Delta_{1}^{\gamma}R_{\alpha\beta}+\Delta_{0}^{\gamma}+\Delta_{2}^{\gamma}\sigma^{\gamma}+\Delta_{3}^{\gamma}\sigma^{\gamma}R_{\alpha\beta})\bar{f}_{j}]+{\rm
H.c.}\Big{\\}}$
$\displaystyle\\!\\!\\!+\lambda\sum_{i}f_{i}^{\dagger}f_{i}+{\rm const},$

where the superscript $\gamma$ labels the orientation of the bond,
$R_{\alpha\beta}^{\gamma}=-\frac{i}{\sqrt{2}}(\sigma^{\alpha}+\sigma^{\beta})$
is a bond-dependent rotation matrixLiu and Normand (2018),
$\bar{f_{i}}=(f_{i\downarrow}^{\dagger},-f_{i\uparrow}^{\dagger})^{T}.$
$\Delta_{0,1,2,3}^{\gamma}$ and $t_{0,1,2,3}^{\gamma}$ are the mean field
parameters (complex numbers) representing the fermion pairing and hopping
terms, respectively. Owing to the $C_{3}$ rotation symmetry the values of
$t^{\gamma}$ and $\Delta^{\gamma}$ are independent on the bond orientation, so
we will only consider the $z$-bonds and omit the bond index in later
discussion. The self consistent mean-field equations for these parameters read

$\displaystyle t_{1}$ $\displaystyle=$ $\displaystyle g_{1}{2\over
N}\sum_{\langle ij\rangle\in z}\langle f_{i}^{\dagger}R_{xy}f_{j}\rangle^{*},$
(6) $\displaystyle t_{0,2}$ $\displaystyle=$ $\displaystyle g_{0,2}{2\over
N}\sum_{\langle ij\rangle\in z}\left[\langle
f_{i}^{\dagger}\sigma^{x}R_{xy}f_{j}\rangle^{*}\pm\langle
f_{i}^{\dagger}\sigma^{y}R_{xy}f_{j}\rangle^{*}\right],$ $\displaystyle t_{3}$
$\displaystyle=$ $\displaystyle g_{3}{2\over N}\sum_{\langle ij\rangle\in
z}\langle f_{i}^{\dagger}\sigma^{z}R_{xy}f_{j}\rangle^{*},$
$\displaystyle\Delta_{1}$ $\displaystyle=$ $\displaystyle
g_{1}^{{}^{\prime}}{2\over N}\sum_{\langle ij\rangle\in z}\langle
f_{i}^{\dagger}R_{xy}\bar{f}_{j}\rangle^{*},$ $\displaystyle\Delta_{0,2}$
$\displaystyle=$ $\displaystyle g_{0,2}^{{}^{\prime}}{2\over N}\sum_{\langle
ij\rangle\in z}\left[\langle
f_{i}^{\dagger}\sigma^{x}R_{xy}\bar{f}_{j}\rangle^{*}\pm\langle
f_{i}^{\dagger}\sigma^{y}R_{xy}f_{j}\rangle^{*}\right],$
$\displaystyle\Delta_{3}$ $\displaystyle=$ $\displaystyle
g_{3}^{{}^{\prime}}{2\over N}\sum_{\langle ij\rangle\in z}\langle
f_{i}^{\dagger}\sigma^{z}R_{xy}\bar{f}_{j}\rangle^{*},$

where $N$ is the total number of sites, the ‘coupling constants’ $g_{0,1,2,3}$
are related to specific form of the spin-spin interactions and will be
specified later. Finally, the Lagrangian multiplier $\lambda$ is introduced to
tune the averaged particle number,

$\displaystyle 1={1\over N}\sum_{i}\langle f_{i}^{\dagger}f_{i}\rangle.$ (7)

Strictly speaking, the ‘particle number constraint’ contains three components
since the SU(2) gauge group have three generators, here we only consider the
third one. At half-filling, the self-consistent solution of $\lambda$ is
usually zero owing to particle-hole symmetry.

In above mean field theory, we only considered the $C_{3}$ rotation symmetry
and the translation symmetry. Actually, a spin liquid should preserve all the
point group symmetries of the spin Hamiltonian. However, in the mean field
description of a QSL phase, the symmetry group of the mean-field Hamiltonian
is the projective symmetry group (PSG) whose group elements are space group
operations followed by $SU(2)$ gauge transformationsWen (2002). We will
restrict to the PSG of the exactly solvable Kitaev point (called the Kitaev
PSG) You et al. (2012); Wang et al. (2019) in the following discussion.

To analyze the PSG symmetry, it is helpful to rewrite the mean-field
Hamiltonian (5) in a matrix form,

$\displaystyle H^{\rm SL}_{\rm mf}\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\sum_{\langle
i,j\rangle\in\alpha\beta(\gamma)}\\!\\!\\!{\rm
Tr}[U^{(0)}_{ji}\psi^{\dagger}_{i}\psi_{j}]+{\rm
Tr}[U^{(1)}_{ji}\psi^{\dagger}_{i}(iR^{\gamma}_{\alpha\beta})\psi_{j}]$ (8)
$\displaystyle+{\rm
Tr}[U^{(2)}_{ji}\psi^{\dagger}_{i}\sigma^{\gamma}\psi_{j}]+{\rm
Tr}[U^{(3)}_{ji}\psi^{\dagger}_{i}\sigma^{\gamma}R^{\gamma}_{\alpha\beta}\psi_{j}]+{\rm
H.c.}$ $\displaystyle+\lambda\sum_{i}{\rm Tr}(\psi_{i}{\tau^{z}\over
4}\psi_{i}^{\dagger})+{\rm const},$

where $U_{ji}^{(0,1,2,3)}$ are the matrix form of the mean-field parameters
$\Delta_{0,1,2,3},t_{0,1,2,3}$. To satisfy the Kitaev PSG, only eight real
mean field parameters remain, namely, $\eta_{0,1,2}$, $\rho_{a,b,c,d}$ and the
Lagrangian multiplier $\lambda$. These eight parameters are real because they
are expected values of quadratic majorana fermion operators, for details see
Appendix B. The matrices $U^{(0,1,2,3)}_{ji}$ can be expanded by these
parameters as

$\displaystyle U_{ji}^{(0)}=i\eta_{0}+i(\rho_{a}+\rho_{b}),$ $\displaystyle
U_{ji}^{(1)}=i\eta_{1}(\tau^{x}\\!+\\!\tau^{y}\\!+\\!\tau^{z})+i(\rho_{a}\\!-\\!\rho_{b}\\!+\\!\rho_{c}\\!+\\!2\rho_{d})(\tau^{\alpha}\\!+\\!\tau^{\beta}),$
$\displaystyle
U_{ji}^{(2)}=i(\rho_{a}\\!+\\!\rho_{b})\tau^{\gamma}+i\rho_{d}(\tau^{\alpha}\\!+\\!\tau^{\beta})+i\eta_{2}(\tau^{x}\\!+\\!\tau^{y}\\!+\\!\tau^{z}),$
$\displaystyle
U_{ji}^{(3)}=i(\rho_{b}\\!-\\!\rho_{a}\\!+\\!\rho_{c})(\tau^{\alpha}-\tau^{\beta}),$
(9)

where $\tau^{x,y,z}$ are Pauli matrices generating the $SU(2)$ gauge group.

Comparing (5) and (6) with (8) and (III.2), we can express $\eta_{0,1,2}$ and
$\rho_{a,b,c,d}$ by the parameters $t_{0,1,2,3},\Delta_{0,1,2,3}$, namely,

$\displaystyle\eta_{0}=-i(t_{0}-\Delta_{2})-it_{2}(1-i),$
$\displaystyle\eta_{1}=i(t_{1}-\Delta_{1}),$ $\displaystyle\eta_{2}=-it_{2},$

and

$\displaystyle\rho_{a}=\frac{i}{2}(1-i)t_{2}+\frac{i}{2}(-t_{1}-t_{3}+\Delta_{1}-\Delta_{2}),$
$\displaystyle\rho_{b}=\frac{i}{2}(1-i)t_{2}+\frac{i}{2}(t_{1}+t_{3}-\Delta_{1}-\Delta_{2}),$
$\displaystyle\rho_{c}=i(-t_{1}+t_{3}+\Delta_{1}),$
$\displaystyle\rho_{d}=\frac{i}{\sqrt{2}}\\{[e^{-i{\pi\over
4}}t_{0}+e^{i{\pi\over 4}}t_{2}+e^{i{\pi\over 4}}\Delta_{1}-e^{-i{\pi\over
4}}\Delta_{3}]-{\rm h.c.}\\}.$

Substituting (6) and (7) into above relations, we obtain the self-consistent
mean-field equations for the set of eight real parameters
$\eta_{0,1,2},\rho_{a,b,c,d}$, $\lambda$, where the coupling constants
$g_{0,1,2,3}$ are given by

$\displaystyle g_{0}=g^{\prime}_{0}=-\frac{1}{8}\Gamma,$ $\displaystyle
g_{1}=-\frac{1}{8}(2\Gamma+|K|+4\Gamma^{{}^{\prime}}),\qquad
g_{1}^{{}^{\prime}}=-\frac{1}{8}(2\Gamma+|K|),$ $\displaystyle
g_{2}=-\frac{1}{8}\Gamma,\qquad
g_{2}^{{}^{\prime}}=-\frac{1}{8}(\Gamma+4\Gamma^{{}^{\prime}}),$
$\displaystyle g_{3}=g^{\prime}_{3}=-\frac{1}{8}|K|.$

Solving the above equations, we obtain a family of QSL solutions, two of which
appear in the phase diagram (see Fig.1 in Sec. IV).

###  III.3 Charge-doped QSLs

In this section, we consider the situations in which holes are doped into the
QSLs. When holes are present, the electrons can lower their energy by hoping
to the vacant sites. From the slave particle representation (3), we can
perform the following mean field decoupling to the electron hoping term,

$\displaystyle H_{\rm mf}^{\rm hop}\\!=\\!-t\sum_{\langle
ij\rangle,\sigma}\\!\Big{(}\langle f_{i\sigma}^{\dagger}f_{j\sigma}\rangle
b_{j}^{\dagger}b_{i}+\langle b_{j}^{\dagger}b_{i}\rangle
f_{i\sigma}^{\dagger}f_{j\sigma}+{\rm H.c.}\Big{)}.$ (10)

In above expression, the boson operators $b,b^{\dagger}$ enters the mean field
theory. We need to introduces one more parameter, namely,

$\displaystyle\chi_{b}^{\gamma}=-t{2\over N}\sum_{\langle
ij\rangle\in\gamma}\langle b_{i}^{\dagger}b_{j}\rangle.$

Another parameter $\chi^{\gamma}=-t{2\over N}\sum_{\langle
ij\rangle\in\gamma,\sigma}\langle f_{i\sigma}^{\dagger}f_{j\sigma}\rangle$ is
not new since it is proportional to $t_{0}$ with $\chi^{\gamma}=-{t\over
g_{0}}t_{0}$. Owing to the $C_{3}$ rotation symmetry, the parameters take the
same values on different bonds $\chi^{x}=\chi^{y}=\chi^{z}$ and
$\chi^{x}_{b}=\chi^{y}_{b}=\chi^{z}_{b}$, so the superscript $\gamma$ will be
omitted later.

Combining of spinon part and the chargon part, we obtain the total mean field
Hamiltonian,

$H_{\rm mf}^{\rm tot}=H_{\rm mf}^{\rm hop}+H_{\rm mf}^{\rm
SL}+\lambda\sum_{i}(f_{i}^{\dagger}f_{i}+b_{i}^{\dagger}b_{i})-\mu\sum_{i}b_{i}^{\dagger}b_{i}.$

Since we are more interested in the region where the magnetic orders are
destroyed by the charge fluctuations(see III.4 for discussion), we do not
consider the magnetic order in the mean field Hamiltonian. Many ansatz have
been tried in determining the phase diagram.

At light doping, namely when $\delta$ is small, the spin-spin interactions
play important roles in the low energy physics. The ansatz (8) preserving the
Kitaev PSG is lowest in free energy at small $\delta$. With the increasing of
$\delta$, the charge part is more important, we find that the system favors a
general QSL ansatz (5) without the restriction of the Kitaev PSG. Especially,
when $\delta$ is very large, the system becomes metalic and the free energy
favors the ansatz in which the paring parameters $\Delta_{0,1,2,3}$ vanish and
the hoping parameters $t_{0,1,2,3}$ are nonzero.

As an example, we write down the total mean field Hamiltonian in the light
doping region. From (8) and (10), we obtain,

$\displaystyle H_{\rm mf}^{\rm tot}$ $\displaystyle=\sum_{\langle
ij\rangle\in\alpha\beta(\gamma)}\big{(}\chi^{*}b_{i}^{\dagger}b_{j}+\chi_{b}^{*}f_{i}^{\dagger}f_{j}+{\rm
H.c.}\big{)}$ (11) $\displaystyle+\sum_{\langle
i,j\rangle\in\alpha\beta(\gamma)}\Big{\\{}{\rm
Tr}[U^{(0)}_{ji}\psi^{\dagger}_{i}\psi_{j}]+{\rm
Tr}[U^{(1)}_{ji}\psi^{\dagger}_{i}(iR^{\gamma}_{\alpha\beta})\psi_{j}]$
$\displaystyle\qquad\quad+{\rm
Tr}[U^{(2)}_{ji}\psi^{\dagger}_{i}\sigma^{\gamma}\psi_{j}]+{\rm
Tr}[U^{(3)}_{ji}\psi^{\dagger}_{i}\sigma^{\gamma}R^{\gamma}_{\alpha\beta}\psi_{j}]\Big{\\}}$
$\displaystyle+\quad\lambda\sum_{i}(f_{i}^{\dagger}f_{i}+b_{i}^{\dagger}b_{i})-\mu\sum_{i}b_{i}^{\dagger}b_{i}+{\rm
const}.$

Notice that the mean field parameter $\chi_{b}=-t\langle
b_{i}^{\dagger}b_{j}\rangle$ (average of boson kinetic term) is generally a
complex number which spontaneously breaks the time reversal symmetry.
$\lambda$ is a Lagrange multiplier for the particle number constraint (2), and
$\mu$ is the chemical potential to adjust the doping concentration
$\delta={1\over N}\sum_{i}b_{i}^{\dagger}b_{i}$.

At half-filling, namely as $\delta=0$, Eq. (11) reduces to the mean field
description of QSLs with the condensation of the spinon pairs. Each spinon
pair carries 2 units of gauge charge, thus the condensation of the spinon
pairs will gap out the $U(1)$ gauge field via Anderson-Higgs mechanism.
Resultantly, in the low-energy limit the gauge fields are described by $Z_{2}$
gauge theory. Upon finite doping, once the bosonic holons undergoes Bose-
Einstein condensation at low temperatures, the $Z_{2}$ gauge field is further
confined. The confinement of the $Z_{2}$ gauge charge indicating that unpaired
spinons cannot be separated, only paired spinons are allowed low-energy
excitations. On the other hand, the $U(1)$ gauge ‘symmetry’ can never be
broken, only the $U(1)$ gauge invariant quantities can have nonzero
expectation values Elitzur (1975). To obtain a $U(1)$ gauge invariant
quantity, one can combine a pair of spinons (which are confined with each
other) and two holons, which are essentially equivalent to two electrons, to
form a Cooper pair. The $U(1)$ gauge invariance of the Cooper pair can be seen
recalling that each spinon $f_{i}$ carries $+1$ gauge charge and each holon
$b_{i}^{\dagger}$ carries $-1$ gauge charge such that a Cooper pair carries
zero gauge charge. The coherent moving of Cooper pairs gives rise to a
superconductor. The nonzero expectation value of the Cooper pair in the ground
state can be considered as the order parameter of the superconducting phase.

Different from the $d$-wave superconductors in cuprates, here the spinon
pairing channel includes both singlet pairing and triplet pairing. Therefore,
the pairing symmetry is neither $s(d)$-wave or $p(f)$-wave, and the spatial
inversion symmetry is spontaneously broken.

###  III.4 Magnetically ordered phases

At half-filling. At zero temperature and at half-filling, the mean field
description of the magnetic ordered states is essentially the classical ground
state since the magnetic order are generally fully saturated. For this reason,
to estimate the ground state energy in the ordered phase, we just use the
classical method (such as single-Q approximationRau et al. (2014)) instead of
taking mean field theory.

To estimation the critical temperature, we adopt the mean-field Hamiltonian
(8) with a background field $\boldsymbol{M}_{i}$ which induces the magnetic
order. Namely, we add an extra term $H_{mf}=\frac{1}{2}\sum_{i}\bm{M}_{i}\cdot
f_{i}^{\dagger}\bm{\sigma}f_{i}$ into (8). The ordering pattern $\bm{M}_{i}$
is determined by the classical solution within the single-$Q$ approximation,
and the magnitude $M=|\boldsymbol{M}_{i}|$ is treated as a variational
parameter determined by minimizing the free energy. The threshold temperature
at which $M$ reduces to zero is the critical temperature $T_{c\rm M}$. In this
approach, the spin-spin interactions are decoupled in two different channels,
namely the spin-liquid channel and the magnetically ordered channel, both
being included in the mean field Hamiltonian. Therefore, when calculating the
free energy of the system from mean-field theory, we have counted the
contribution from the spin-spin interactions twice. This may result in an
overcounting of the free energy and may introduce a systematic error in
estimating the critical temperature $T_{c\rm M}$.

Finite doping at zero temperature $T=0$. In the absence of doping, the energy
per site in the classically ordered state is equal to half of the total energy
of the three bonds connected to it. When the doping concentration is nonzero
$\delta\neq 0$, we assume that the holes are uniformly distributed among the
lattice sites (namely, we assume that the holes are well separated). Then we
approximately count the total energy of the ordered phase by assuming that the
each hole increases the total energy by the amount which the interactions on
three neighboring bonds would have if there were no holes.

On the other hand, at finite doping the superconducting states are competing
in energy. When the energy of the ordered state meets that of a
superconducting state, then the critical concentration $\delta_{c\rm M}$ is
obtained.

(a) phase diagram

(b) FM (c) Zigzag (d) spinon dispersion

Figure 1: (a) Phase diagram at half-filling with fixed
$\Gamma^{{}^{\prime}}/|K|=-0.02$. KSL denotes the Kitaev spin liquid, FM the
ferromagnetic phase,QSL14 the 14-cone QSL phase, IS the incommensurate spiral
phase, and Zigzag the zigzag ordered phase. (b) and (c) show the spin
configurations of the FM and the Zigzag order, respectively. (d) The spinon
dispersion in the QSL14 phase.

##  IV The PHASE DIAGRAMs

The phase diagram is obtained by minimizing the free energy. We have tried
different self-consistent solutions of the mean field equations. For a given
set of interaction parameters, the one with the lowest free energy (or energy
at zero temperature) defines the phase in the phase diagram.

###  IV.1 Quantum Phase Diagram ($T=0$) at Half-filling

The spin model with $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ interactions have bee
profoundly studied. Since most phases are magnetically ordered, we will focus
on possible QSLs. In our mean field theory, we indeed obtain two QSL phases.
Besides the well known KSL, we obtain an additional gapless QSL phase with 14
majorana cones in the spinon excitation spectrum(a similar 14-cone QSL state
called PKSL14 was obtained in Ref.Wang et al. (2019)). The zero-temperature
mean field phase diagram is shown in Fig.1.

At the small $\Gamma$ side, the KSL locates in the region $0\leq\Gamma/|K|\leq
0.08$. The KSL with $K<0$ has a ferromagnetic instability. At larger $\Gamma$
with $0.08<\Gamma/|K|<0.25$, a ferromagnetic(FM) phase is obtained Wang et al.
(2019, 2020); Lee et al. (2020).

At the large $\Gamma$ side, the system favors the zigzag order for
$\Gamma/|K|>0.78$. The zigzag order has been observed in serval Kitaev
candidatesSears et al. (2015); Johnson et al. (2015); Cao et al. (2016); Ye et
al. (2012); Choi et al. (2012); Williams et al. (2016). It was shown that the
zigzag phase is absent in the $K$-$\Gamma$ model, for this reason, we
introduce a small $\Gamma^{\prime}=-0.02|K|$ to stabilize it. With the
decreasing $\Gamma$, the zigzag order is replaced by an incommensurate spiral
(IS) phase at $0.48<\Gamma/|K|<0.78$.

Finally, resulting from the competition between the FM phase and the IS phase,
there exist another interesting QSL phase (labeled as QSL14) living in the
region $0.25<\Gamma/|K|<0.48$. This QSL contains 14 majorana cones (see Fig.1
(d)) and shares the same PSG as the KSL phase. Although a similar state was
obtained from variational Monte Carlo method, the two 14-cone states are
different. When a small magnetic field is applied along
$\hat{c}={1\over\sqrt{3}}(\hat{x}+\hat{y}+\hat{z})$ direction, the 14-cone QSL
from self-consistent solution is turned into a gapped state with Chern number
$\nu=-1$, but the PKSL14 acquires a Chern number $\nu=5$ under the same
magnetic field. In this region, there is another spin liquid solution with
eight cones (labeled as QSL8), but with a higher energy.

Figure 2: Phase diagrams of doped Kitaev system (with fixed
$\Gamma^{\prime}/|K|=-0.02$). (a) doped zigzag phase with $\Gamma/|K|=1$; (b)
doped QSL14 phase with $\Gamma/|K|=0.4$; (c) doped FM phase with
$\Gamma/|K|=0.15$; (d) doped IS phase with $\Gamma/|K|=0.6$. SC denotes the
superconducting phase, SC0 means the Chern number is $\nu=0$, and SC2zz stands
for a SC phase with Chern number $\nu=2$ doped from the zigzag phase, so on
and so forth. SG denotes the spin gapped phase, FL the fermi liquid phase, SM
the strange metal phase and PM the paramagnetic phase.

Figure 3: First order phase transitions with increasing $\delta$ at $T=0$
illustrated by the Energy curves. (a) Doped zigzag phase with
$\Gamma/|K|=1,\Gamma^{{}^{\prime}}/|K|=-0.02$. The level crossings indicate a
series of phase transitions from the Zigzag phase to the SC2zz, from the SC2zz
to the SC0 and from the SC0 phase to the FL. (b) Doped QSL14 with
$\Gamma/|K|=0.4,\Gamma^{{}^{\prime}}/|K|=-0.02$. The level crossing indicates
the transition from the SC1SL phase to the SC0SL phase.

Figure 4:  Phase transitions of doped QSL14 (with
$\Gamma/|K|=0.4,\Gamma^{{}^{\prime}}/|K|=-0.02$) at finite temperature
illustrated by the free energy curves. (a) $\delta=0.1$, the level crossing
indicates a first order transition from the SC0SL to the SC0. (b)
$\delta=0.02$, the converging of the two curves indicates a second order
continuous phase transition from the SG0 phase to the PM phase. The insert
illustrates the first order derivative of the free energy.

###  IV.2 Finite Doping

In this section, we study the physical consequence when holes are doped into
the system. Many interesting phases have been obtained, including
(topological) superconductors, spin-gapped phases, the fermi liquid phase, the
strange metal phase and the paramagnetic phase. Most of these phases appeared
in the $t$-$J$ model as discussed in Ref.Lee and Nagaosa (1992); Ubbens and
Lee (1992), but our phase diagram is even richer [see Fig.2(a)$\sim$(d)] owing
to spin-orbit coupling in the spin-spin interactions.

Different mean field phases. In mean field theory, the fermionic spinons and
the bosonic holons are decoupled, so the total free energy is a sum of that of
the spinons and that of the holons. We tried several self-consist solutions,
the one with the lowest free energy is adopt to describe the state of the
system.

Here we provide the mean-field description of the phases listed above. The
superconducting (SC) phase is characterized by the Bose-Einstein condensation
of the holons and the condensation of the spinon pairs. If the spinon pairing
parameter is nonzero but the holons are not condensed due to thermal
fluctuations, then the resultant state belongs to the spin-gapped phase (SG)
(also called the pseudogapped phase). On the other hand, if the holons
condense but the spinon pairing parameter is zero, then the resultant phase is
the fermi liquid (FL) phase. If holons are not condensed, and if the spinons
do not form pairs but can coherently hop, then the resultant state is the
strange metal (SM) phase. When both the spinons and the holons lose their
coherence (namely, the all of the parameters $t_{0,1,2,3}$ and
$\Delta_{0,1,2,3}$ vanish), and the system enters the high temperature
paramagnetic (PM) phase.

In the following we illustrate the steps of obtaining the phase diagram.

Firstly we investigate the phase boundary of the magnetically ordered phases.
In section III.4, we have estimated the critical temperature $T_{c\rm M}$ and
the critical concentration $\delta_{c\rm M}$ (see Fig.3(a) for example) for
the magnetic transition. From $T_{c\rm M}$ and $\delta_{c\rm M}$, an
approximate phase boundary of the magnetic ordered phase can be figured out by
drawing a smooth curve linking $T_{c\rm M}$ and $\delta_{c\rm M}$.

Then we estimate the critical temperature $T_{c\rm HC}$ for the holon
condensation. In two-dimensions bosonic systems cannot spontaneously break
continuous symmetries at any finite temperature. In our study, the $T_{c\rm
HC}$ stands for the Kosterlitz-Thouless transition temperature and is
evaluated approximatelyKosterlitz and Thouless (1973).

With the increasing of doping concentration $\delta$, another important
feature is the critical concentration $\delta_{c\rm SP}$ at which the spinons
pairing parameter drops to zero. To this end, we compare the ansatz (8) with
nonzero spinon-pairing (the SC states at lower $T$ or the SG states at higher
$T$) and another ansatz (5) with vanishing spinon pairing terms (the FL state
at lower $T$ or the SM state at higher $T$). It turns out that the latter
ansatz is favored by the free energy when the doping concentration exceeds the
critical $\delta_{c\rm SP}$. An example of such transition is shown in
Fig.3(a).

As mentioned in section III.3, different types of spinon ansatz have been
considered to describe the SC phases. It turns out that at small $\delta$
(after the magnetic order being suppressed) and at low $T$, the spinon ansatz
preserving the Kitaev PSG is lower in free energy. Above a critical $\delta$
or a critical $T$, the spinons prefer another pairing state outside the Kitaev
PSG class. Fig.4 (a) shows a first order transition between two SC phases with
increasing $T$, where the low temperature SC0SL preserves the Kitaev PSG and
the high temperature SC0 does not. Another example can be found in Fig.3(a),
where a first order transition between SC2zz (preserving the Kitaev PSG) and
SC0 (violating the Kitaev PSG) appears at a critical $\delta$ at zero
temperature. Owing to this subtle difference, doped Kitaev system contains
more than one SC phase and more than one pseudogap phase.

Finally, the paramagnetic phase is characterized by the vanishing of all
parameters, in which the spinons and holons lose their coherence. There is a
continuous phase transition between SM and PM (see Fig.4 (b)).

From the above procedure, the global phase diagram is obtained. In the
following, we analyze the doped ordered phases and doped QSLs one-by-one. The
doped KSL has been studied in Ref. You et al. (2012) and will be skipped here.

Doped zigzag phase. The doped zigzag state is most interesting because it may
be relevant to several Kitaev materials. Fig. 2(a) shows the phase diagram of
doped zigzag phase with $\Gamma/|K|=1,\Gamma^{{}^{\prime}}/|K|=-0.02$. The
solid lines denote phase transitions obtained by minimizing the free energy of
self-consistent mean-field solutions, while the dotted lines are obtained by
estimation using the method given previously.

Interestingly, if the doping concentration is close to $\delta=0.1$, we obtain
a topological superconductor SC2zz with Chern number $\nu={2}$ or $\nu=-2$. In
this region, the spinon states preserves the Kitaev PSG, and the nonzero Chern
number results from the spontaneous breaking of time reversal symmetry in the
complex parameter $\chi_{b}$ and $\chi_{f}$. The critical temperature of the
superconductivity is around $0.08|K|$ (which is of the order of $10$ Kelvin
for $\alpha$-RuCl3Ran et al. (2017); Li et al. (2021)). The SG2zz phase
appears when the temperature exceeds the $T_{c\rm HC}$, in which the holons
are not condensed but the spinons stay in the same state as in the SC2zz
phase. With increasing $\delta\sim 0.2$, another superconductor SC0 with Chern
number $\nu=0$ and the corresponding SG0 are obtained, in which the spinons
still form pairs but are no longer preserving the Kitaev PSG. As
$\delta>0.32$, the system becomes metalic and enters the FL or SM phase.

When the doping concentration reaches some value, the resultant physics is not
sensitive to the spin-spin interactions. Therefore, the phase diagrams
Fig.2(a)$\sim$(d) in the region $\delta>0.14$ are very similar.

Doped QSL14. Fig. 2(b) shows the phase diagram of doped QSL14 phase with
$\Gamma/|K|=0.4,\Gamma^{{}^{\prime}}/|K|=-0.02$. Since in the QSL14 phase the
spinons have already formed pairs, the system becomes superconducting once
holes are doped with $0<\delta\leq 0.32$.

Three different SC phases are obtained, one has Chern number $\nu=\pm{1}$ and
the other two have $\nu=0$. The first two SC phases, SC1SL and SC0SL respect
the Kiatev PSG and are separated at $\delta=0.06$ with a first-order phase
transition (see Fig.3 (b)). In the third SC phase SC0, the spinon state falls
outside the Kitaev PSG class and the pairing gap decreases with increasing
$\delta$.

Doped FM and IS phases. The phase diagrams of the doped FM and IS states are
shown in Fig. 2(c) and (d), respectively, which are very similar to the case
with doped zigzag state shown in Fig.2(a). The main difference is that the
Chern number of the topological SC phase doped from the FM and IS phases are
$\nu=\pm{1}$.

###  IV.3 Effect of In-plane Magnetic Fields

Below the magnetic transition temperature, the zigzag order in Kitaev
materials (such as $\alpha$-RuCl3) can be easily suppressed by in-plane
magnetic field, resulting in a QSL-like disordered phase. In this section, we
study the effect of in-plane magnetic field in the doped zigzag phase, which
may be relevant for future experimental study.

We consider the in-plane magnetic field, say, along the
${1\over\sqrt{2}}(\hat{\boldsymbol{x}}-\hat{\boldsymbol{y}})$ direction. In
this case, the spinons gain Zeeman energy

$H^{\rm zeeman}_{\rm mf}={\mu_{B}B\over
2\sqrt{2}}\sum_{i}f_{i}^{\dagger}(\sigma^{x}-\sigma^{y})f_{i},$

but the cyclotron motion of the charge carriers is avoided. Thus the holons
are not affected by the in-plane magnetic field. Here we only discuss the
situations at zero temperature.

Since the in-plane magnetic field frustrates the zigzag order, it will need
less doped holes to destroyed the order and to drive the system to the
superconducting phase. Namely, the in-plane magnetic field can reduce the
critical doping concentration $\delta_{c\rm M}$ to the SC phase. Similarly,
doped holes also decreases the lower critical magnetic field strength in
suppressing the magnetic order.

Figure 5: Zero temperature phase diagram with in-plane magnetic fields
$\boldsymbol{B}\parallel{1\over\sqrt{2}}(\hat{x}-\hat{y})$ applied to the
hole-doped zigzag state.

On the other hand, the magnetic field tends to polarize the spinons and breaks
the spinon pairing, so if the magnetic field reaches the upper critical value,
the spin of the electrons will be polarized and the system enters the trivial
polarized phase. For the same reason, the in-plane magnetic field results in a
smaller critical $\delta_{c\rm SP}$ for the system to enter the fermi liquid
phase. The phase diagram is shown in Fig.5. It can be seen that the size of
the topological SC2zz phase is slightly enlarged while the nontopological SC0
phase is shortened at the intermediate region of intensity of the magnetic
field.

It should be noted that vertical line with $\delta=0$ is special, where the
intermediate phase between the Zigzag phase and the polarized phase is a
gapless quantum spin liquid phase with two Majorana cones (similar results are
obtained previouslyWang et al. (2020) using variational Monte Carlo method).
Once holes are doped into the system such that $\delta\neq 0$, the holons will
condense and at the mean while the cones will open a gap, resulting in the
topological superconducting phase SC2zz.

Comparing Fig.5 and Fig.2(a), it can be seen that the effect of in-plane
magnetic field is very similar to that of the temperature. Similar results can
be obtain in the doped FM phase or the doped IS phase and will not be shown
here.

##  V Conclusions and Discussions

In summary, to explore the effect of doped holes in Kitaev systems, we studied
the $t$-$K$-$\Gamma$-$\Gamma^{{}^{\prime}}$ model via mean field theory.
Firstly we (almost) reproduce the zero temperature phase diagram of the pure
spin model with $K$-$\Gamma$-$\Gamma^{{}^{\prime}}$ interactions at zero
doping, which contains three magnetically ordered phases and two QSL phases
(the KSL and the QSL14 which contains 14 majorana cones). By doping holes to
these phases, we obtain the superconducting phases, the pseudogap phases, the
fermi liquid phase, the strange metal phase and the paramagnetic phase.
Topological superconductors are obtained, where the Chern number is dependent
on the original spin state. By doping the zigzag phase one has Chern number
$\nu=\pm{2}$, by doping the FM phase or the IS phase one get $\nu=\pm{1}$.
Interestingly, no matter what is the original spin state, more than one
superconductor and more than one spin gapped state are obtained after doping.
We studied the effect of applied in-plane magnetic field and find that it can
slightly enlarge the size of the topological superconducting phase.

It should be clarified that the conclusions obtained from our mean field
theory is very preliminary. Except for the magnetically ordered phases, we
only consider the ansatz which preserving the $C_{3}$ and translation
symmetry. Other symmetry breaking orders, such as the charge/spin density
wave, nematicity, and valence bond solid order may appear at low temperatures
but are not considered. Furthermore, quantum fluctuations around the mean
field solution may quantitatively or qualitatively change the results. For
instance, in the variational Monte Carlo phase diagram using Gutzwiller
projected states as trial states, a 14-cone QSL called PKSL14 (proximate
Kitaev spin liquid with 14 cones) is obtained. But this PKSL14 is different
from the QSL14 obtained from the self-consistent mean field solution, since
the two states have different Chern number in a small magnetic field applied
along $\boldsymbol{B}\parallel{1\over\sqrt{3}}(\hat{x}+\hat{y}+\hat{z})$
direction. Therefore, it will be interesting to study the present model using
partially projected wave functions as trial states and determine the
parameters variationally. We leave this for future study.

Finally, the $t$-$K$-$\Gamma$-$\Gamma^{\prime}$ model is based on the single-
band approximation. It is uncertain if this approximation is good to describes
the low energy physics of doped Kitaev materials. If not, it deserves further
study to find the true effective model at finite doping.

Owing to the obtained rich phase diagram, our work may intrigue interest of
future experimental studies.

Acknowledgement – We thank J.-C. Wang, Q.R. Zhao and F. Yang for valuable
discussions. This work is supported by the Ministry of Science and Technology
of China (Grant No. 2016YFA0300504), the NSF of China (Grants No. 11574392 and
No. 11974421), and the Fundamental Research Funds for the Central Universities
and the Research Funds of Renmin University of China (No. 19XNLG11).

##  Appendix A General form of Mean field Hamiltonian and projective symmetry
groups

In the main text, we introduce the fermionic spinon representation where the
spin operators can be written as the familiar quadratic form of fermions
$S^{\alpha}_{i}=\frac{1}{2}f^{\dagger}_{i}\sigma_{\alpha}f_{i}$ under the
constraint $f^{\dagger}_{i}f_{i}=1$, where
$f^{\dagger}_{i}=(f^{\dagger}_{i\uparrow},f^{\dagger}_{i\downarrow})$ and
$\sigma_{\alpha},\alpha=x,y,z$ are the Pauli matrix. The two spinon species
may further be expressed in terms of four Majorana fermions,

$\displaystyle f_{\uparrow}=\frac{1}{2}(b^{z}+ic),\qquad
f_{\downarrow}=\frac{1}{2}(b^{x}+ib^{y}),$

which satisfy the anti-commutation relations
$\\{b^{\alpha},b^{\beta}\\}=2\delta^{\alpha\beta}$
($\alpha,\beta=0,x,y,z;b^{0}\equiv c$). In this basis, the spin operator takes
the form $S^{m}=ib^{m}c$.

There is a SU(2) gauge symmetryCookmeyer and Moore (2018) in the spinon
representation. To see this, we introduce
$\bar{f}_{i}=(f^{\dagger}_{i\downarrow}-f^{\dagger}_{i\uparrow})^{T}$, the
time reversal partner of $f^{\dagger}_{i}$, and a matrix operator
$\psi_{i}=(f_{i},\bar{f}_{i})$.

The mixing between $f$ and $\bar{f}$, namely, multiplying $\psi$ by a SU(2)
matrix from the right hand side, does not change the spin operators
$S^{m}_{i}={\rm Tr}(\psi^{\dagger}_{i}\frac{\sigma_{m}}{4}\psi_{i})$, which
defines the SU(2) gauge symmetry.

For rotationally invariant interactions. The Heisenberg exchange interaction
can be written as

$\displaystyle\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{j}\\!=\\!-\frac{1}{8}{\rm
Tr}(\psi^{\dagger}_{j}\psi_{i}\psi^{\dagger}_{i}\psi_{j})\\!=\\!\frac{1}{8}{\rm
Tr}(\psi^{\dagger}_{j}\bm{\sigma}\psi_{i}\cdot\psi^{\dagger}_{i}\bm{\sigma}\psi_{j})$
(1)

up to constant terms. The anisotropic interactions, such as the AFM Ising
interaction, can be decoupled as

$\displaystyle S^{m}_{i}S^{m}_{j}=-\frac{1}{16}[{\rm
Tr}(\psi^{\dagger}_{j}\psi_{i}\psi^{\dagger}_{i}\psi_{j})+{\rm
Tr}(\psi^{\dagger}_{j}\sigma^{m}\psi_{i}\psi^{\dagger}_{i}\sigma^{m}\psi_{j})]$
(2)

up to a constant term, where $m=x,y,z$.

To make more transparent the connection to the spinon representation, we note
that the singlet matrix operator $\psi^{\dagger}_{i}\psi_{j}$, and the triplet
operator $\psi^{\dagger}_{i}\bm{\sigma}\psi_{j}$, can be expanded as

$\displaystyle\psi^{\dagger}_{i}\psi_{j}=\begin{pmatrix}f^{\dagger}_{i}f_{j}&f^{\dagger}_{i}\bar{f}_{j}\\\
\bar{f}^{\dagger}_{i}f_{j}&\bar{f}^{\dagger}_{i}\bar{f}_{j}\end{pmatrix},\psi^{\dagger}_{i}\bm{\sigma}\psi_{j}=\begin{pmatrix}f^{\dagger}_{i}\bm{\sigma}f_{j}&f^{\dagger}_{i}\bm{\sigma}\bar{f}_{j}\\\
\bar{f}^{\dagger}_{i}\bm{\sigma}f_{j}&\bar{f}^{\dagger}_{i}\bm{\sigma}\bar{f}_{j}\end{pmatrix},$

respectively.

Actually, from above expression of the spin-spin interactions, we obtain the
most general mean-field Hamiltonian (with only nearest neighbor couplings) for
a spin-liquid state with spin-orbit coupling,

$\displaystyle H^{\rm SL}_{\rm mf}$ $\displaystyle=$
$\displaystyle\sum_{ij}{\rm
Tr}\Big{(}U^{(0)}_{ji}\psi^{\dagger}_{i}\psi_{j}+\sum_{m=x,y,z}{U}^{(m)}_{ji}\psi^{\dagger}_{i}{\sigma^{m}}\psi_{j}+{\rm
H.c.}\Big{)}$ (3)
$\displaystyle+\sum_{i}\boldsymbol{\lambda}\cdot\boldsymbol{\varLambda_{i}},$

where $U_{ji}^{(0)}$ is the same as that in (8), and the matrices
$U_{ji}^{(x,y,z)}$ are linearly related to the ones $U_{ji}^{(1,2,3)}$. For
example, in the $z$-bond (namely for $\langle ij\rangle\in z$), one has

$\displaystyle U_{ji}^{(x)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(U^{(1)}_{ji}-U^{(3)}_{ji}),$ $\displaystyle
U_{ji}^{(y)}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(U^{(1)}_{ji}+U^{(3)}_{ji}),$ $\displaystyle
U_{ji}^{(z)}$ $\displaystyle=$ $\displaystyle U^{(2)}_{ji}.$

When expanding the matrices $U_{ji}^{(0,1,2,3)}$ by the Pauli bases, the
coefficients correspond to the parameters $t_{0,1,2,3}$ and $\Delta_{0,1,2,3}$
in (5). The advantage of the nation $U_{ji}^{(0,1,2,3)}$ is that if the mean
field Hamiltonian preserves $C_{3}$ symmetry, then the values of
$U_{ji}^{(0,1,2,3)}$ are independent on the bond directions.

The operators $\boldsymbol{\Lambda}_{i}=\frac{1}{4}{\rm
Tr}(\psi_{i}\bm{\tau}\psi^{\dagger}_{i})$ are the generators of the SU(2)
gauge group, and the third component of $\boldsymbol{\lambda}$ is the
Lagrangian multiplier of the particle number constraint appeared in the main
text.

As shown in appendix B, above Hamiltonian (3) is the general expression of the
mean-field Hamiltonian even if the system contain off-diagonal interactions
(with only nearest neighbor interactions).

The mean-field Hamiltonian of Eq. (3) is not in general invariant under an
arbitrary SU(2) gauge transformation. The subgroup of the SU(2) gauge group
under which Eq. (3) remains invariant is called the invariant gauge group
(IGG) of the spin-liquid state. On the other hand, a QSL should respect all
the symmetries of the spin Hamiltonian, but at the mean-field level this
constraint can be relaxed in the following sense. Under a symmetry operation
$g$, $H_{mf}$ may be transformed to a different expression,
$\hat{g}H_{mf}\hat{g}^{-1}=H^{{}^{\prime}}_{mf}\neq H_{mf}$, but if it can be
transformed back to its original form by an SU(2) gauge transformation then
this mean-field Hamiltonian still describes a spin-liquid state. Specifically,

$\displaystyle\psi_{i}\to\hat{g}^{\dagger}\psi_{g(i)}W_{i}(g),$ (4a)
$\displaystyle{\rm Tr}[U_{ji}\psi^{\dagger}_{i}\psi_{j}]$
$\displaystyle\to{\rm
Tr}[W_{j}U_{ji}W^{\dagger}_{i}\psi^{\dagger}_{g(i)}\hat{g}\hat{g}^{\dagger}\psi_{g(j)}],$
$\displaystyle={\rm Tr}[U_{g(j)g(i)}\psi^{\dagger}_{g(i)}\psi_{g(j)}],$ (4b)
$\displaystyle{\rm Tr}[\bm{U}_{ji}\cdot\psi^{\dagger}_{i}\bm{\sigma}\psi_{j}]$
$\displaystyle\to{\rm
Tr}[W_{j}\bm{U}_{ji}W^{\dagger}_{i}\cdot\psi^{\dagger}_{g(i)}\hat{g}\bm{\sigma}\hat{g}^{\dagger}\psi_{g(j)}],$
$\displaystyle={\rm
Tr}[\bm{U}_{g(j)g(i)}\cdot\psi^{\dagger}_{g(i)}\bm{\sigma}\psi_{g(j)}]$ (4c)

where $\hat{g}$ is the double-valued representation of $g$,
$U^{n}_{g(j)g(i)}=\sum_{m}R_{nm}(g)W_{j}(g)U^{m}_{ji}W^{\dagger}_{i}(g)$,
$U_{g(j)g(i)}=W_{j}(g)U_{ji}W^{\dagger}_{i}(g)$ and $R(g)$ is a vector
representation of $g$. The new symmetry operations, each of which involves a
symmetry operation $g$ followed by a gauge transformation $W_{i}(g)$, form a
larger group that is known as the projective symmetry group (PSG)Wang et al.
(2019, 2020).

Now we provide the PSG of the Kitaev Spin Liquid (namely, the Kitaev PSG),
which can be read out in the Majorana representation. The symmetry group
$G=D_{3d}\times Z^{T}_{2}$ has three generators

$S_{6}=(C_{3})^{2}P,\quad M=C^{x-y}_{2}P,\quad T=i\sigma^{y}K,$

where $C_{3}$ is a threefold rotation around the direction
$\hat{c}\equiv\frac{1}{\sqrt{3}}(\hat{x}+\hat{y}+\hat{z})$, $C^{x-y}_{2}$ is a
twofold rotation around $\frac{1}{\sqrt{2}}(\hat{x}-\hat{y})$, and $P$ is
spatial inversion. From (4a), one obtain the Kitaev PSG,

$\displaystyle
W_{A}(S_{6})=-W_{B}(S_{6})=\exp[-i\frac{4\pi}{3}\frac{1}{2\sqrt{3}}(\tau^{x}+\tau^{y}+\tau^{z})],$
$\displaystyle
W_{A}(M)=-W_{B}(M)=\exp[-i\pi\frac{1}{2\sqrt{2}}(\tau^{x}-\tau^{y})],$
$\displaystyle W_{A}(T)=-W_{B}(T)=i\tau^{y},$

where A and B denote the two sublattices of the honeycomb lattice.

##  Appendix B Decoupling of the $K$-$\Gamma$-$\Gamma^{\prime}$ interactions

In the Majorana representation, the mean-field decouplings of the Kitaev
interactions which preserve the Kitaev PSG read

$\displaystyle H^{K}_{mf}$ $\displaystyle=\sum_{\langle
i,j\rangle\in\alpha\beta(\gamma)}\rho_{a}(ic_{i}c_{j})+\rho_{c}(ib^{\gamma}_{i}b^{\gamma}_{j})$
$\displaystyle=\sum_{\langle i,j\rangle\in\alpha\beta(\gamma)}i\rho_{a}{\rm
Tr}(\psi^{\dagger}_{i}\psi_{j}+\tau^{x}\psi^{\dagger}\sigma^{x}\psi_{j}+\tau^{y}\psi^{\dagger}_{i}\sigma^{y}\psi_{j}$
$\displaystyle\qquad\quad\quad+\tau^{z}\psi^{\dagger}\sigma^{z}\psi_{j})+i\rho_{c}{\rm
Tr}(\psi^{\dagger}_{i}\psi_{j}+\tau^{\gamma}\psi^{\dagger}\sigma^{\gamma}\psi_{j}$
$\displaystyle\qquad\quad\quad-\tau^{\alpha}\psi^{\dagger}_{i}\sigma^{\alpha}\psi_{j}-\tau^{\beta}\psi^{\dagger}\sigma^{\beta}\psi_{j})+{\rm
H.c.}$ (5)

Similarly, the $\Gamma$ interaction decouples as

$\displaystyle H^{\Gamma}_{mf}$ $\displaystyle=\sum_{\langle
i,j\rangle\in\alpha\beta(\gamma)}i\rho_{d}(b^{\alpha}_{i}b^{\beta}_{j}+b^{\beta}_{i}b^{\alpha}_{j})$
$\displaystyle=\sum_{\langle i,j\rangle\in\alpha\beta(\gamma)}i\rho_{d}{\rm
Tr}(\tau^{\alpha}\psi^{\dagger}_{i}\sigma^{\beta}\psi_{j}+\tau^{\beta}\psi^{\dagger}_{i}\sigma^{\alpha}\psi_{j})+{\rm
H.c.}$ (6)

and for the $\Gamma^{{}^{\prime}}$ interaction

$\displaystyle H^{\Gamma^{{}^{\prime}}}_{mf}$ $\displaystyle=\sum_{\langle
i,j\rangle\in\alpha\beta(\gamma)}i\rho_{f}(b^{\alpha}_{i}b^{\gamma}_{j}+b^{\gamma}_{i}b^{\alpha}_{j}+b^{\beta}_{i}b^{\gamma}_{j}+b^{\gamma}_{i}b^{\beta}_{j})$
$\displaystyle=\sum_{\langle i,j\rangle\in\alpha\beta(\gamma)}i\rho_{f}{\rm
Tr}(\tau^{\alpha}\psi^{\dagger}_{i}\sigma^{\gamma}\psi_{j}+\tau^{\gamma}\psi^{\dagger}_{i}\sigma^{\alpha}\psi_{j}$
$\displaystyle\qquad\qquad+\tau^{\beta}\psi^{\dagger}_{i}\sigma^{\gamma}\psi_{j}+\tau^{\gamma}\psi^{\dagger}_{i}\sigma^{\beta}\psi_{j})+{\rm
H.c.}$ (7)

However, the most general coefficients preserving the $C_{3}$ rotation
symmetry also contain multiples of the uniform ($\tau^{0}$) and
$\tau^{x}+\tau^{y}+\tau^{z}$ gauge components. Thus we need to add extra terms
in the mean-field Hamiltonian as long as they are symmetry allowed. The terms
we found that preserve all the Kitaev PSG symmetries contain three more
parameters $\eta_{0,1,2}$,

$\displaystyle H^{\rm ex}_{\rm mf}\\!$
$\displaystyle=\\!\\!\\!\\!\sum_{\langle
i,j\rangle\in\alpha\beta(\gamma)}\\!\\!\\!i\eta_{0}{\rm
Tr}(\psi^{\dagger}_{i}\psi_{j})\\!+\\!i\eta_{2}{\rm
Tr}[(\tau^{x}\\!+\\!\tau^{y}\\!+\\!\tau^{z})\psi^{\dagger}_{i}\sigma^{\gamma}\psi_{j}]$
$\displaystyle\qquad\qquad+i\eta_{1}{\rm
Tr}[(\tau^{x}+\tau^{y}+\tau^{z})\psi^{\dagger}_{i}(iR^{\gamma}_{\alpha\beta})\psi_{j}]+{\rm
H.c.}$ (8)

Combining the equations (5), (6), (7) and (8), we obtain the complete mean-
field Hamiltonian (8) and (III.2) which describe QSL states preserving the
Kitaev PSG.

## References

  * Balents (2010) L. Balents, Nature 464, 199 (2010).
  * Kitaev (2006) A. Kitaev, Annals of Physics 321, 2 (2006).
  * Zhou et al. (2017) Y. Zhou, K. Kanoda, and T.-K. Ng, Reviews of Modern Physics 89, 025003 (2017).
  * Chandra and Doucot (1988) P. Chandra and B. Doucot, Physical Review B 38, 9335 (1988).
  * He et al. (2014) Y.-C. He, D. Sheng, and Y. Chen, Physical review letters 112, 137202 (2014).
  * Yan et al. (2011) S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
  * Shimizu et al. (2003) Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, and G. Saito, Physical review letters 91, 107001 (2003).
  * Okamoto et al. (2007) Y. Okamoto, M. Nohara, H. Aruga-Katori, and H. Takagi, Physical review letters 99, 137207 (2007).
  * Han et al. (2012) T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Nature 492, 406 (2012).
  * Chillal et al. (2020) S. Chillal, Y. Iqbal, H. O. Jeschke, J. A. Rodriguez-Rivera, R. Bewley, P. Manuel, D. Khalyavin, P. Steffens, R. Thomale, A. N. Islam, et al., Nature communications 11, 1 (2020).
  * Jackeli and Khaliullin (2009) G. Jackeli and G. Khaliullin, Physical review letters 102, 017205 (2009).
  * Chaloupka et al. (2010) J. Chaloupka, G. Jackeli, and G. Khaliullin, Physical review letters 105, 027204 (2010).
  * Janssen et al. (2017) L. Janssen, E. C. Andrade, and M. Vojta, Physical Review B 96, 064430 (2017).
  * Winter et al. (2017) S. M. Winter, K. Riedl, P. A. Maksimov, A. L. Chernyshev, A. Honecker, and R. Valentí, Nature communications 8, 1 (2017).
  * Kimchi and You (2011) I. Kimchi and Y.-Z. You, Physical Review B 84, 180407 (2011).
  * Winter et al. (2016) S. M. Winter, Y. Li, H. O. Jeschke, and R. Valentí, Physical Review B 93, 214431 (2016).
  * Singh et al. (2012) Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Physical review letters 108, 127203 (2012).
  * Rau and Kee (2014) J. G. Rau and H.-Y. Kee, arXiv preprint arXiv:1408.4811 (2014).
  * Chaloupka and Khaliullin (2015) J. Chaloupka and G. Khaliullin, Physical Review B 92, 024413 (2015).
  * Gordon et al. (2019) J. S. Gordon, A. Catuneanu, E. S. Sørensen, and H.-Y. Kee, Nature communications 10, 1 (2019).
  * Lee et al. (2020) H.-Y. Lee, R. Kaneko, L. E. Chern, T. Okubo, Y. Yamaji, N. Kawashima, and Y. B. Kim, Nature communications 11, 1 (2020).
  * Baek et al. (2017) S.-H. Baek, S.-H. Do, K.-Y. Choi, Y. S. Kwon, A. Wolter, S. Nishimoto, J. Van Den Brink, and B. Büchner, Physical review letters 119, 037201 (2017).
  * Zheng et al. (2017) J. Zheng, K. Ran, T. Li, J. Wang, P. Wang, B. Liu, Z.-X. Liu, B. Normand, J. Wen, and W. Yu, Physical review letters 119, 227208 (2017).
  * Yadav et al. (2016) R. Yadav, N. A. Bogdanov, V. M. Katukuri, S. Nishimoto, J. Van Den Brink, and L. Hozoi, Scientific reports 6, 1 (2016).
  * Mirebeau et al. (2002) I. Mirebeau, I. Goncharenko, P. Cadavez-Peres, S. Bramwell, M. Gingras, and J. Gardner, Nature 420, 54 (2002).
  * Kozlenko et al. (2008) D. Kozlenko, I. Mirebeau, J.-G. Park, I. Goncharenko, S. Lee, J. Park, and B. Savenko, Physical Review B 78, 054401 (2008).
  * Wang et al. (2018) Z. Wang, J. Guo, F. Tafti, A. Hegg, S. Sen, V. A. Sidorov, L. Wang, S. Cai, W. Yi, Y. Zhou, et al., Physical Review B 97, 245149 (2018).
  * Cui et al. (2017) Y. Cui, J. Zheng, K. Ran, J. Wen, Z.-X. Liu, B. Liu, W. Guo, and W. Yu, Physical Review B 96, 205147 (2017).
  * Anderson (1987) P. W. Anderson, science 235, 1196 (1987).
  * Lee and Nagaosa (1992) P. A. Lee and N. Nagaosa, Physical Review B 46, 5621 (1992).
  * Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Wu et al. (1987) M.-K. Wu, J. R. Ashburn, C. Torng, P. H. Hor, R. L. Meng, L. Gao, Z. J. Huang, Y. Wang, and a. Chu, Physical review letters 58, 908 (1987).
  * Jiang et al. (2018) H.-C. Jiang, Z.-Y. Weng, and S. A. Kivelson, Physical Review B 98, 140505 (2018).
  * Ruan et al. (2016) W. Ruan, C. Hu, J. Zhao, P. Cai, Y. Peng, C. Ye, R. Yu, X. Li, Z. Hao, C. Jin, et al., Science bulletin 61, 1826 (2016).
  * Mei et al. (2012) J.-W. Mei, S. Kawasaki, G.-Q. Zheng, Z.-Y. Weng, and X.-G. Wen, Physical Review B 85, 134519 (2012).
  * Weng (2011) Z.-Y. Weng, New Journal of Physics 13, 103039 (2011).
  * You et al. (2012) Y.-Z. You, I. Kimchi, and A. Vishwanath, Physical Review B 86, 085145 (2012).
  * Liu and Normand (2018) Z.-X. Liu and B. Normand, Physical review letters 120, 187201 (2018).
  * Banerjee et al. (2017) A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Science 356, 1055 (2017).
  * Wang et al. (2017) W. Wang, Z.-Y. Dong, S.-L. Yu, and J.-X. Li, Physical Review B 96, 115103 (2017).
  * Cookmeyer and Moore (2018) J. Cookmeyer and J. E. Moore, Physical Review B 98, 060412 (2018).
  * Lee et al. (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Reviews of modern physics 78, 17 (2006).
  * Affleck et al. (1988) I. Affleck, Z. Zou, T. Hsu, and P. Anderson, Physical Review B 38, 745 (1988).
  * Wen and Lee (1996) X.-G. Wen and P. A. Lee, Physical Review Letters 76, 503 (1996).
  * Lee et al. (1998) P. A. Lee, N. Nagaosa, T.-K. Ng, and X.-G. Wen, Physical Review B 57, 6003 (1998).
  * Wang et al. (2019) J. Wang, B. Normand, and Z.-X. Liu, Physical review letters 123, 197201 (2019).
  * Wang et al. (2020) J. Wang, Q. Zhao, X. Wang, and Z.-X. Liu, Physical Review B 102, 144427 (2020).
  * Rau et al. (2014) J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Physical review letters 112, 077204 (2014).
  * Wen (2002) X.-G. Wen, Physical Review B 65, 165113 (2002).
  * Elitzur (1975) S. Elitzur, Phys. Rev. D 12, 3978 (1975), URL https://link.aps.org/doi/10.1103/PhysRevD.12.3978.
  * Sears et al. (2015) J. A. Sears, M. Songvilay, K. Plumb, J. Clancy, Y. Qiu, Y. Zhao, D. Parshall, and Y.-J. Kim, Physical Review B 91, 144420 (2015).
  * Johnson et al. (2015) R. D. Johnson, S. Williams, A. Haghighirad, J. Singleton, V. Zapf, P. Manuel, I. Mazin, Y. Li, H. O. Jeschke, R. Valentí, et al., Physical Review B 92, 235119 (2015).
  * Cao et al. (2016) H. B. Cao, A. Banerjee, J.-Q. Yan, C. Bridges, M. Lumsden, D. Mandrus, D. Tennant, B. Chakoumakos, and S. Nagler, Physical Review B 93, 134423 (2016).
  * Ye et al. (2012) F. Ye, S. Chi, H. Cao, B. C. Chakoumakos, J. A. Fernandez-Baca, R. Custelcean, T. Qi, O. Korneta, and G. Cao, Physical Review B 85, 180403 (2012).
  * Choi et al. (2012) S. Choi, R. Coldea, A. Kolmogorov, T. Lancaster, I. Mazin, S. Blundell, P. Radaelli, Y. Singh, P. Gegenwart, K. Choi, et al., Physical review letters 108, 127204 (2012).
  * Williams et al. (2016) S. Williams, R. Johnson, F. Freund, S. Choi, A. Jesche, I. Kimchi, S. Manni, A. Bombardi, P. Manuel, P. Gegenwart, et al., Physical Review B 93, 195158 (2016).
  * Ubbens and Lee (1992) M. U. Ubbens and P. A. Lee, Physical Review B 46, 8434 (1992).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Ran et al. (2017) K. Ran, J. Wang, W. Wang, Z.-Y. Dong, X. Ren, S. Bao, S. Li, Z. Ma, Y. Gan, Y. Zhang, et al., Phys. Rev. Lett. 118, 107203 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.107203.
  * Li et al. (2021) H. Li, H.-K. Zhang, J. Wang, H.-Q. Wu, Y. Gao, D.-W. Qu, Z.-X. Liu, S.-S. Gong, and W. Li, Nature Communications 12, 4007 (2021), URL https://doi.org/10.1038/s41467-021-24257-8.