# Detecting symmetry fractionalization in gapped quantum spin liquids by
magnetic impurities

Shuangyuan Lu  Department of Physics, The Ohio State University, Columbus OH
43210, USA  Yuan-Ming Lu  Department of Physics, The Ohio State University,
Columbus OH 43210, USA

###### Abstract

We study the Kondo effect of spin-$1/2$ magnetic impurities in gapped $Z_{2}$
spin liquids on two-dimensional lattices. We find that if the impurity is
placed at a high-symmetry location, a nontrivial spinon fractionalization
class of the impurity site symmetry group will necessarily lead to a non-
Kramers doublet in the Kondo screening regime, protected by associated
crystalline symmetries. This is in sharp contrast to a featureless screening
phase in the usual Kondo effect. We demonstrate this symmetry-protected
topological degeneracy by an exactly solvable model and by the large-$N$
theory. Based on this effect, we discuss how thermodynamic measurements in the
limit of dilute magnetic impurities can be used to detect symmetry
fractionalization in gapped $Z_{2}$ spin liquids.

##  I Introduction

Quantum spin liquids (QSLs) Balents (2010); Savary and Balents (2016); Zhou
_et al._ (2017) have attracted much interests in the past few decades due to
its exotic properties transcending the Landau scheme of symmetry breaking. In
particular, the presence of anyons which obey fractional statisticsWilczek
(1990), is among the most exciting manifestations of the topological order and
long-range entanglement in QSLsWen (2002, 2007), with potential applications
in topological quantum computations Nayak  _et al._ (2008). A number of QSL
candidate materials with various crystalline symmetries have been discovered
experimentally Kanoda and Kato (2022); Zhou  _et al._ (2017); Wen  _et al._
(2019); Takagi  _et al._ (2019).

Meanwhile, there is a gap between theoretical diagnosis and experimental
measurements to identify QSLs. On one hand, various theoretically computable
quantities has been proposed to sharply characterize topological orders, such
as the topological entanglement entropyKitaev and Preskill (2006); Levin and
Wen (2006) and modular matricesKESKI-VAKKURI and XIAO-GANG (2022); Zhang  _et
al._ (2012). On the other hand, most existing experiments aim at ruling out
long-range orders in the low temperature, deterred by the difficulty of
directly probing unique features of QSLs. In particular, compared to gapless
$U(1)$ spin liquids with clear signatures in inelastic neutron scattering
(INS) Benton  _et al._ (2012) or thermal transport Kanoda and Kato (2022);
Zhou  _et al._ (2017), a gapped symmetric QSL is more featureless and harder
to detect experimentally. While the long-range entanglement and fractional
statistics, as a definitive character of topological orders, is difficult to
access experimentally, the fractional symmetry quantum numberLaughlin (1999)
(formally known as symmetry fractionalizationWen (2002); Essin and Hermele
(2013); Tarantino  _et al._ (2016); Barkeshli  _et al._ (2019); Chen (2017))
of anyons provide extra features to characterize and identify the topological
order, which is usually easier to probe experimentally. In the well known
example of fractional quantum Hall effects (FQHEs), indeed the fractional
charge is experimentally observed in the nineties, much earlier than the
recently confirmed fractional statisticsFeldman and Halperin (2021). One
question arises naturally: can symmetry fractionalization be experimentally
detected as a direct evidence for a QSL state?

For gapped QSLs, which is the focus of this work, there are two major
challenges to experimental detection of fractionalization. First of all, the
fractionalized excitations such as spinons are charge neutral, therefore
insensitive to charge transport probes which played a crucial role in
identifying fractionalization in FQHEsFeldman and Halperin (2021). Secondly, a
gapped symmetric QSL usually has no features both in the bulk and on the
boundary, making it very hard to access experimentally. This is unlike the
$U(1)$ spin liquids, whose gapless excitations can be probed by INS in the
case of emergent photonsBenton  _et al._ (2012), or thermal transport in the
case of spinon Fermi surfacesKanoda and Kato (2022); Zhou  _et al._ (2017). Is
it possible to experimentally identify a gapped QSL? Previously, INS
spectroscopy has been proposed to exhibit features of fractional
statisticsMorampudi  _et al._ (2017) and spinon symmetry fractionalizationWen
(2002); Essin and Hermele (2014). In this work, we look into magnetic
impurities and Kondo effects in gapped QSLs, and show that they can provide
distinct thermodynamic signatures of symmetry frationalization in QSLs, in the
Kondo screening regime.

The Kondo effect in QSLs has previously been studied both in theories
Khaliullin and Fulde (1995); Dhochak  _et al._ (2010); Doretto and Vojta
(2009); Ribeiro and Lee (2011); Vojta  _et al._ (2016); Das  _et al._ (2016);
Kolezhuk  _et al._ (2006); He and Lee (2022) and in experiments Yamamoto  _et
al._ (2018); Gomilšek  _et al._ (2019); Chen  _et al._ (2022), focusing on
gapless QSLs. In this paper, we explore the Kondo effect in gapped $Z_{2}$
QSLs. Similar to the distinctions between Kondo effects in metals and in
insulators (with a vanishing density of states), the Kondo effect in gapped
$Z_{2}$ QSLs differs qualitatively from gapless QSLs. In particular, due to
the energy gap for spinon excitations, there is a finite threshold of Kondo
coupling strength to screen the magnetic impurityWithoff and Fradkin (1990);
Satori  _et al._ (1992); Saso (1992); Itoh (1993); Takegahara  _et al._
(1993); Chen and Jayaprakash (1998). Most remarkably, we find that when a
half-integer-spin impurity is placed at a high-symmetry location in the
crystal hosting a gapped $Z_{2}$ QSL, the Kondo screening phase will feature a
non-Kramers doublet localized at the impurity site, protected by
fractionalized crystalline symmetries in the $Z_{2}$ QSL. This symmetry
protected degeneracy lead to distinct signatures in the thermodynamics, such
as specific heat, which can serve as “smoking gun” evidence for symmetry
fractionalization in a gapped QSL. This phenomena is demonstrated by an
exactly solvable model and large-$N$ parton mean-field theory, which agree
with each other.

##  II Main results

Figure 1: High-symmetry impurity sites that can be used to detect the symmetry
fractionalization of spinons in a gapped symmetric $Z_{2}$ spin liquid, on the
(a) square, (b) kagome and (c) triangular lattices in two dimensions.

We first present the major results in Fig. 1 and TABLEs 1-3. Consider a gapped
symmetric $Z_{2}$ spin liquid on a two-dimensional lattice, whose Hilbert
space consists of a spin-$1/2$ (or a Kramers doublet) on each lattice site. A
spin-$1/2$ magnetic impurity located at certain high-symmetry position of the
lattice can be used to diagnose the symmetry fractionalization classEssin and
Hermele (2013); Barkeshli  _et al._ (2019); Tarantino  _et al._ (2016) of the
$Z_{2}$ spin liquid phase. Specifically, when such a magnetic impurity is
coupled symmetrically to the $Z_{2}$ spin liquid, in the Kondo screening
regime, there may or may not be a two-fold degeneracy (a non-Kramers doublet)
protected by the crystalline symmetry of the impurity site, depending on the
fractionalization class of spinons in the $Z_{2}$ spin liquid.

Fig. 1 illustrates the three lattices enumerated in this work, i.e. the
square, kagome and triangular lattices. In the presence of $SO(3)$ spin
rotational symmetry, the classification of symmetric $Z_{2}$ spin liquids on
these latticesLu (2018); Qi and Cheng (2018) are summarized in TABLE 1-3. A
part of the symmetry fractionalization data can be detected by
presence/absence of non-Kramers doublets for Kondo-screened magnetic
impurities located at different high-symmetry sites, such as a plaquette
center, a nearest-neighbor link center, or on a mirror plane.

Algebraic identity | $\omega\in\mathcal{H}^{2}(G,\mathcal{A})$ |  $\omega^{e}$Yang and Wang (2016) |  $\omega^{\epsilon}$Wen (2002) | Impurity site
---|---|---|---|---
$(R_{x})^{2}$ | $\omega_{R_{x},R_{x}}$ | $(-1)^{p_{4}}$ | $\eta_{\sigma}$ | -
$(M_{y})^{2}$ | $\omega_{M_{y},M_{y}}$ | $(-1)^{p_{3}+p_{4}}$ | $\eta_{\sigma}\eta_{xpx}$ | -
$(C_{4}R_{x})^{2}$ | $\omega_{C_{4}R_{x},C_{4}R_{x}}$ | $(-1)^{p_{4}+p_{7}}$ | $\eta_{\sigma}\eta_{\sigma C_{4}}$ | -
$M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}$ | $\frac{\omega_{M_{x},M_{y}}}{\omega_{M_{y},M_{x}}}$ | $(-1)^{p_{1}}$ | $\eta_{xy}$ | A
$(M_{y}\mathcal{T})^{2}$ | $\omega_{M_{y}\mathcal{T},M_{y}\mathcal{T}}$ | $(-1)^{p_{3}+p_{8}+1}$ | $-\eta_{t}\eta_{xpx}$ | B
$R_{x}M_{y}R_{x}^{-1}M_{y}^{-1}$ | $\frac{\omega_{R_{x},M_{y}}}{\omega_{M_{y},R_{x}}}$ | $(-1)^{p_{2}}$ | $\eta_{xpy}$ | C
Table 1: All gapped $Z_{2}$ spin liquids of spin-$1/2$’s on the square lattice ($G=p4m\times Z_{2}^{\mathcal{T}}$), characterized by $2^{6}$ fractionalization classesLu (2018); Qi and Cheng (2018), and their realizations in the Schwinger bosonYang and Wang (2016) and Abrikosov fermionWen (2002) representations. 3 of the 6 independent $Z_{2}$ invariants can be detected by magnetic impurities located at A, B and C sites in Fig. 1(a).  Algebraic identity | $\omega\in\mathcal{H}^{2}(G,\mathcal{A})$ |  $\omega^{e}$Yang and Wang (2016) |  $\omega^{\epsilon}$Wen (2002) | Impurity site
---|---|---|---|---
$(R_{x})^{2}$ | $\omega_{R_{x},R_{x}}$ | $(-1)^{p_{2}+p_{3}}$ | $\eta_{\sigma}$ | -
$(R_{y})^{2}$ | $\omega_{R_{y},R_{y}}$ | $(-1)^{p_{2}}$ | $\eta_{\sigma}\eta_{\sigma C_{6}}$ | -
$R_{x}R_{y}R_{x}^{-1}R_{y}^{-1}$ | $\frac{\omega_{R_{x},R_{y}}}{\omega_{R_{y},R_{x}}}$ | $(-1)^{p_{1}}$ | $\eta_{12}$ | A
Table 2: All gapped $Z_{2}$ spin liquids of spin-$1/2$’s on the kagome lattice ($G=p6mm\times Z_{2}^{\mathcal{T}}$), characterized by $2^{3}$ fractionalization classesLu (2018); Qi and Cheng (2018), and their realizations in the Schwinger bosonWang and Vishwanath (2006) and Abrikosov fermionLu  _et al._ (2011) representations. 1 of the 3 independent $Z_{2}$ invariants can be detected by magnetic impurities located at A site in Fig. 1(b).  Algebraic identity | $\omega\in\mathcal{H}^{2}(G,\mathcal{A})$ |  $\omega^{e}$Yang and Wang (2016) |  $\omega^{\epsilon}$Wen (2002) | Impurity site
---|---|---|---|---
$\sigma^{2}$ | $\omega_{\sigma,\sigma}$ | $(-1)^{p_{2}}$ | $\eta_{\sigma}$ | -
$(M_{y})^{2}$ | $\omega_{M_{y},M_{y}}$ | $(-1)^{p_{2}+p_{3}}$ | $\eta_{\sigma C_{6}}$ | -
$\sigma M_{y}\sigma^{-1}M_{y}^{-1}$ | $\frac{\omega_{\sigma,M_{y}}}{\omega_{M_{y},\sigma}}$ | $(-1)^{p_{1}}$ | $\eta_{12}$ | A
Table 3: All gapped $Z_{2}$ spin liquids of spin-$1/2$’s on the triangular
lattice ($G=p6mm\times Z_{2}^{\mathcal{T}}$), characterized by $2^{3}$
fractionalization classesLu (2018); Qi and Cheng (2018), and their
realizations in the Schwinger bosonWang and Vishwanath (2006) and Abrikosov
fermionLu (2016) representations. 1 of the 3 independent $Z_{2}$ invariants
can be detected by magnetic impurities located at A site in Fig. 1(c).

In a gapped system with a vanishing density of states, a finite Kondo coupling
strength is required to enter the Kondo screening regimeWithoff and Fradkin
(1990); Satori  _et al._ (1992); Saso (1992); Itoh (1993); Takegahara  _et
al._ (1993); Chen and Jayaprakash (1998), where the system typically reaches a
featureless paramagnetic ground state. The proposed non-Kramers doublet
protected by crystalline symmetries at the impurity site of the $Z_{2}$ spin
liquid is therefore a striking and unusual phenomenon. Below we describe the
physical picture behind this observation.

A gapped $Z_{2}$ spin liquid hosts three types of anyons (or superselection
sectors): bosonic spinon $e$, vison $m$, and their bound state
$\epsilon=e\times m$ known as a fermionic spinonKitaev (2003). In a symmetric
$Z_{2}$ spin liquid on a lattice with an odd number of spin-$1/2$’s in each
unit cell, spinons $e$ and $\epsilon$ must carry spin-$1/2$ each, while vison
$m$ is spinlessZaletel and Vishwanath (2015). As a result, when a spin-$1/2$
magnetic impurity is coupled to such a $Z_{2}$ spin liquid, to reach a spin-
singlet ground state, it can only be screened by a spinon $e$ or $\epsilon$.
The same conclusion holds if we replace spin-$1/2$ by a Kramers doublet with
$\mathcal{T}^{2}=-1$ in the argument.

In the presence of crystalline and time reversal symmetries, different gapped
$Z_{2}$ spin liquids are distinguished by their symmetry fractionalization
classes, classified by 2nd group cohomology $\mathcal{H}^{2}(G,\mathcal{A})$,
where the symmetry group is $G=SG\times Z_{2}^{\mathcal{T}}$ ($SG$ being the
space group), and $\mathcal{A}=Z_{2}\times Z_{2}$ is the fusion group of
Abelian anyons in the $Z_{2}$ spin liquidEssin and Hermele (2013). Thanks to
the $SO(3)$ spin rotational symmetry, the vison fractionalization class is
uniquely fixed on the three latticesQi  _et al._ (2015); Lu  _et al._ (2017),
leading to the classification shown in Table 1-3Qi and Cheng (2018); Lu
(2018), characterized by projective representations of $G$ carried by spinons.
Consequently, the singlet bound state formed by the impurity spin and the
screening cloud of spinon can carry a projective representation of the
impurity site symmetry group. This leads to the non-Kramers degeneracy at the
impurity site, protected by both crystalline and time reversal symmetries.

There are two types of impurity sites of particular interests to this work:
(i) the impurity site lies at the intersection of two mirror planes, such as
site A and C in Fig.1; (ii) the impurity site lies on a mirror plane, such as
site B in Fig. 1. In case (i), since the two mirror symmetries $M_{x}$ and
$M_{y}$ commute in the Hilbert space of the impurity spin, if they anticommute
on the spinon screening the impurity spin (i.e.
$\omega_{M_{x},M_{y}}/\omega_{M_{y},M_{x}}=-1$), the bound state of impurity
spin-$1/2$ and spinon will carry a projective representation of the site
symmetry group $G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$, leading to a 2-fold
degeneracy protected by the mirror symmetries. On the other hand, if the two
mirror actions commute on the spinon (i.e.
$\omega_{M_{x},M_{y}}/\omega_{M_{y},M_{x}}=+1$), the bound state will carry a
linear representation of the site symmetry and hence no degeneracy in the
screening regime. In case (ii), the impurity spin-$1/2$ on a mirror ($M_{x}$)
plane carries a projective representation $(M_{x}\mathcal{T})^{2}=-1$. If the
spinon screening the impurity carry a linear representation with
$\omega_{M_{x}\mathcal{T},M_{x}\mathcal{T}}=+1$, their bound state forms a
projective representation of the site symmetry group
$G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{\mathcal{T}}$, leading to a 2-fold non-
Kramers degeneracy protected by both mirror $M_{x}$ and time reversal
symmetries. In contrary, if the spinon carries a projective representation
with $\omega_{M_{x}\mathcal{T},M_{x}\mathcal{T}}=-1$, the bound state would
instead form a linear representation $(M_{x}\mathcal{T})^{2}=+1$ with no
degeneracy in the screened phase. This shows how the response of a $Z_{2}$
spin liquid to impurity spin-$1/2$’s can diagnose a part of the spinon
fractionalization data in the spin liquid phase, as summarized in Table 1-3.

Figure 2: An exactly solvable model for a spin-$1/2$ impurity in a $Z_{2}$
spin liquid with $SU(2)$ symmetry. The squares on link centers are qubits in
the toric codeKitaev (2003). Each black dot on a star/vertex represents a
Hilbert space of spin $0\oplus 1/2$. The red dot denotes the spin-$1/2$
impurity with a site symmetry $G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$sup .

##  III Methods

We use two methods to demonstrate the symmetry-protected non-Kramers
degeneracy induced by a spin-$1/2$ impurity in $Z_{2}$ spin liquids, in the
Kondo screening regime. We focus on case (i) of impurity site symmetry
$G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$, where the impurity spin-$1/2$ lies
at the intersection of two mirror planes $M_{x,y}$. First we construct an
exactly solvable model by modifying the toric code, to show the exact
degeneracy protected by two mirror symmetries. Next we use the large-$N$
approach to solve the Kondo problem in a symmetric $Z_{2}$ spin liquid, and to
compute the temperature dependence of thermodynamic quantities.

First we present an exactly solvable model illustrated in Fig. 2. The bulk
$Z_{2}$ spin liquid is described by:

$\displaystyle\hat{H}_{bulk}=-\sum_{s}A_{s}-\sum_{p}B_{p}$ (1)
$\displaystyle-\sum_{s}\Delta\big{[}\frac{(A_{s}+1)}{2}P_{s}(S=0)+\frac{(1-A_{s})}{2}P_{s}(S=1/2)\big{]}$

where $A_{s}$ and $B_{p}$ are the star and plaquette operators in Kitaev’s
toric codeKitaev (2003). In addition to one qubits on each link, there is a
3-dimensional Hilbert space of spin $0\oplus\frac{1}{2}$ on each site/vertex
(see Fig. 2). In the limit $\Delta\gg 1$, each $e$ and $\epsilon$ particle
will carry spin-$1/2$ of the site/vertex Hilbert space, while $m$ particles
are spinless. The fractionalization class associated with the impurity site
symmetry group $G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$ is given byEssin and
Hermele (2013)

$\displaystyle\frac{\omega^{e}_{M_{x},M_{y}}}{\omega^{e}_{M_{y},M_{x}}}=\frac{\omega^{m}_{M_{x},M_{y}}}{\omega^{m}_{M_{y},M_{x}}}=+1,~{}~{}~{}\frac{\omega^{\epsilon}_{M_{x},M_{y}}}{\omega^{\epsilon}_{M_{y},M_{x}}}=-1.$
(2)

Next we introduce a spin-$1/2$ impurity located at the center of plaquette
$(1234)$ in Fig. 2, which is coupled to the bulk spin liquid as follows:

$H_{imp}=J\sum_{i=1}^{4}\vec{S}_{i}\cdot\vec{S}_{imp}+E_{c}(A_{1}+A_{2}+A_{3}+A_{4}-3)^{2}+\Delta_{\epsilon}B_{p}[1234]$
(3)

In addition to the usual Kondo coupling $J$, we also introduce (i) a Coulomb
repulsion $E_{c}$ for spinons, which makes sure the impurity is screened by
one spinon, and (ii) an energy $\Delta_{\epsilon}$ coupled to the plaquette
opeartor on plaquette $[1234]$, to control which type of spinons ($e$ vs.
$\epsilon$) will screen the impurity spin. Assuming $E_{c}\gg 1,J$, the Kondo
screening regime happens when $J>4/3$, leading to 4 degenerate states in the
low energy manifold where a single spinon is located at one neighboring site
(out of $1,2,3,4$). When $\Delta_{\epsilon}<1$, the bosonic spinons cost less
energy and will screen the impurity, and the 4-fold degenerate can be
completely lifted with a unique paramagnetic ground state without breaking any
symmetrysup . This is consistent with the trivial fractionalization class of
$e$ particles in (2). When $\Delta_{\epsilon}>1$, however, the fermionic
spinons cost lower energy and are responsible for the Kondo screening. As
detailed in supplemental materialssup , the 4-dimensional low energy space can
be split into two doublets, each of which form an irreducible projective
representation of the impurity site symmetry group $G_{s}=Z_{2}^{M_{x}}\times
Z_{2}^{M_{y}}$. As a result, a 2-fold degeneracy protected by two mirror
symmetries $M_{x,y}$ will emerge in the Kondo screening regime, as indicated
by the nontrivial fractionalization class (2) of $\epsilon$ particles.
Therefore we have demonstrated the correspondence between nontrivial
fractionalization class of spinons screening the impurity, and protected
2-fold degeneracy in the Kondo screening regime.

Next we use a large-$N$ mean field theory to solve the Kondo problem in
symmetric $Z_{2}$ spin liquids. Both the bulk and impurity spins are
represented by fermonic partons with $Sp(2N)$ symmetryRan and Wen (2006):

$\displaystyle
S^{ab+}=\frac{1}{2}(c^{a\dagger}_{\uparrow}c_{\downarrow}^{b}+c^{b\dagger}_{\uparrow}c_{\downarrow}^{a}),~{}~{}~{}S^{ab,z}=\frac{1}{2}(c^{a\dagger}_{\uparrow}c_{\uparrow}^{b}-c^{b\dagger}_{\downarrow}c_{\downarrow}^{a})$
(4)

with $1\leq a,b\leq N$. They reduce to the familiar $SU(2)$ spin-$1/2$ case
when $N=1$. The model consists of a $Z_{2}$ spin liquid in the bulk described
by parton mean-field ansatz of $Sp(2N)$ partons

$H_{bulk}=\sum_{a=1}^{N}\sum_{i,j}J_{i,j}\psi_{i}^{a,\dagger}u_{i,j}\psi^{a}_{j}+h.c.$
(5)

where we denote
$\psi_{i}^{a}=(c_{i,\uparrow}^{a},c^{a\dagger}_{i,\downarrow})^{T}$, and the
Kondo coupling between a $Sp(2N)$ impurity spin and its neighboring spins:

$H_{imp}=\sum_{\langle
j,imp\rangle}\frac{J}{N}\mathbf{S}_{j}^{ab}\cdot\mathbf{S}_{imp}^{ba}+\frac{J^{\prime}}{N^{3}}(\mathbf{S}_{j}^{ab}\cdot\mathbf{S}_{imp}^{ba})^{2}$
(6)

As detailed in supplemental materialssup , the bulk parton ansatz can be
exactly realized in solvable models in analogy to Kitaev’s honeycomb
modelKitaev (2006), and choosing different link parameters $\\{u_{ij}\\}$ can
lead to either trivial or nontrivial fractionalization classes for fermionic
spinons $\\{\psi_{i}^{a}\\}$, with $M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=\pm 1$. A
self-consistent mean-field calculation, which becomes exact in the large $N$
limit, reveals a Kondo screening phase separated from the unscreened phase by
a Kondo temperature $T_{K}(J)$, for Kondo couplings beyond a finite threshold
$J>J_{c}$sup . In the Kondo screening regime with $T<T_{K}(J)$, while the
trivial fractionalization class ($M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=+1$) shows a
unique paramagnetic ground state, the nontrivial fractionalization class
($M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=-1$) exhibits two degenerate ground states
which cannot be mixed by any local perturbations preserving mirror symmetries
$M_{x,y}$sup . This again demonstrated our conclusion, that a nontrivial
spinon fractionalization class will lead to symmetry protected zero modes
localized at high-symmetry impurity sites.

Figure 3:  The temperature dependence of (a) uniform magnetic susceptibility
$\chi(T)$ and (b) entropy $S(T)$ contributed by Kondo impurities in different
regimes: the unscreened regime of free moments at the impurity sites (green),
the Kondo screening regime in $Z_{2}$ spin liquids with a trivial (blue) vs.
nontrivial (red) spinon fractionalization class. The calculations are
performed for two distant impurities with site symmetry
$G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$ on a $20\times 20$ latticesup .

##  IV Experimental implications

The large-$N$ mean field theory also allows us to predict distinct
experimental signatures of the anomalous Kondo screening phase described
above. The temperature dependence of uniform magnetic susceptibility $\chi(T)$
and the entropy $S(T)=\int_{0}^{T}\frac{C_{v}(t)}{t}\text{d}t$ are shown in
Fig. 3. The shown impurity contribution shown in the figure susceptibility
$\chi(T)$ and specific heat $C_{v}(T)$, by subtracting the bulk contribution
of Hamiltonian $H_{bulk}$ from the total amount of $H_{bulk}+H_{imp}$. Three
different regimes can be differentiated from each other by inspecting the
susceptibility and entropy (by integrating the specific heat) at low
temperatures: (i) In the unscreened regime, the magnetic impurity behaves as a
free moment, leading to $\chi(T)\sim 1/T$ and a finite entropy of $k_{B}\ln 2$
per impurity, colored green in Fig. 3. (ii) In the Kondo screening regime, for
$Z_{2}$ spin liquid with trivial fractionalization class, Kondo screening
leaves a unique paramagnetic ground state below the Kondo temperature, and
therefore exponentially decaying thermodynamic responses $\chi(T),C_{v}(T)\sim
e^{-\Delta/k_{B}T}$, as colored blue in Fig. 3. (iii) A $Z_{2}$ spin liquid
with a nontrivial fractionalization class, on the other hand, features a
symmetry protected non-Kramers doublet (2-fold degeneracy) localized at each
impurity site in the Kondo screening regime. As a result, while the
susceptiblity vanishes exponentially at low temperatures $\chi(T)\sim
e^{-\Delta/k_{B}T}$, there is a low energy entropy of $k_{B}\ln 2$ per
impurity below the Kondo temperature, as colored red in Fig. 3. The sharp
differences between the three scenarios provide clear experimental features to
detect a gapped $Z_{2}$ spin liquid with a nontrivial symmetry
fractionalization class for spinons.

Theoretically we only discussed the case of isolated impurities in the
large-$N$ self-consistent mean-field theory described earlier. In real
materials, we expect our predictions in Fig. 3 to hold in the case of _dilute
magnetic impurities_ , where the average distance $r$ between neighboring
impurities is much larger than the bulk correlation length $\xi$ of the gapped
spin liquid. In this case, the splitting of the symmetry protected zero modes
localized at impurity sites will be small $\sim Je^{-C_{0}r/\xi}$, where
$C_{0}$ is a constant of order one. This will lead to a peak in specific heat
at low temperature $T\sim Je^{-C_{0}r/\xi}$.

##  V Concluding remarks

In this work, we show that when magnetic impurities with half integer spins
are coupled to symmetric $Z_{2}$ spin liquids in an isotropic magnet with
$SU(2)$ symmetry, a spinon will form a singlet bound state with the impurity
in the Kondo screening regime. This bound state will feature a symmetry
protected 2-fold degeneracy, if the spinon fractionalization class of the
impurity site symmetry is nontrivial, therefore leading to a non-Kramers
doublet localized at the impurity site. We further show that this local
degeneracy in the Kondo screening regime can be distinguished from other
scenarios by the low temperature behaviors of magnetic susceptibiltiy and
specific heat, hence unveiling a new way to detect symmetry fractionalization
in QSLs.

In the future, it will be insightful to apply the new angle proposed in this
work to examine the candidate materials of QSLs, where magnetic impurities are
known to exist at high symmetry sites e.g. in HerbertsmithiteFu  _et al._
(2015). Theoretically, this work also provides a new idea to probe
fractionalization in models of strongly correlated electrons, by studying the
impurity problem e.g. using numerical methods.

###### Acknowledgements.

We thank Biao Huang for discussions and related collaborations at an early
stage of this work. This work is supported by National Science Foundation
(NSF) under award number DMR 1653769 (YML), by Center for Emergent Materials
at The Ohio State University, a NSF MRSEC through NSF Award No. DMR-2011876
(SL).

## References

  * Balents (2010) L. Balents, Nature 464, 199 (2010).
  * Savary and Balents (2016) L. Savary and L. Balents, Rep. Prog. Phys. 80, 016502 (2016).
  * Zhou  _et al._ (2017) Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Wilczek (1990) F. Wilczek, ed.,  _Fractional Statistics and Anyon Superconductivity_ (World Scientific, 1990).
  * Wen (2002) X.-G. Wen, Phys. Rev. B 65, 165113 (2002).
  * Wen (2007) X.-G. Wen,  _Quantum Field Theory of Many-body Systems: From the Origin of Sound to an Origin of Light and Electrons_ (Oxford University Press, 2007).
  * Nayak  _et al._ (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Kanoda and Kato (2022) K. Kanoda and R. Kato, Annu. Rev. Condens. Matter Phys. 2, 167 (2022).
  * Wen  _et al._ (2019) J. Wen, S.-L. Yu, S. Li, W. Yu, and J.-X. Li, npj Quantum Materials 4, 12 (2019).
  * Takagi  _et al._ (2019) H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Nature Reviews Physics 1, 264 (2019).
  * Kitaev and Preskill (2006) A. Kitaev and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006).
  * Levin and Wen (2006) M. Levin and X.-G. Wen, Phys. Rev. Lett. 96, 110405 (2006).
  * KESKI-VAKKURI and XIAO-GANG (2022) E. S. K. O. KESKI-VAKKURI and W. E. N. XIAO-GANG, Int. J. Mod. Phys. B 07, 4227 (2022).
  * Zhang  _et al._ (2012) Y. Zhang, T. Grover, A. Turner, M. Oshikawa, and A. Vishwanath, Phys. Rev. B 85, 235151 (2012).
  * Benton  _et al._ (2012) O. Benton, O. Sikora, and N. Shannon, Phys Rev B 86, 075154 (2012).
  * Laughlin (1999) R. B. Laughlin, Rev. Mod. Phys. 71, 863 (1999).
  * Essin and Hermele (2013) A. M. Essin and M. Hermele, Phys. Rev. B 87, 104406 (2013).
  * Tarantino  _et al._ (2016) N. Tarantino, N. H. Lindner, and L. Fidkowski, New J Phys 18, 035006 (2016).
  * Barkeshli  _et al._ (2019) M. Barkeshli, P. Bonderson, M. Cheng, and Z. Wang, Phys Rev B 100, 115147 (2019).
  * Chen (2017) X. Chen, Rev. Phys. 2, 3 (2017).
  * Feldman and Halperin (2021) D. E. Feldman and B. I. Halperin, Rep. Progr. Phys. 84, 076501 (2021).
  * Morampudi  _et al._ (2017) S. C. Morampudi, A. M. Turner, F. Pollmann, and F. Wilczek, Phys. Rev. Lett. 118, 227201 (2017).
  * Essin and Hermele (2014) A. M. Essin and M. Hermele, Phys. Rev. B 90, 121102 (2014).
  * Khaliullin and Fulde (1995) G. Khaliullin and P. Fulde, Phys. Rev. B 52, 9514 (1995).
  * Dhochak  _et al._ (2010) K. Dhochak, R. Shankar, and V. Tripathi, Phys Rev Lett 105, 117201 (2010).
  * Doretto and Vojta (2009) R. L. Doretto and M. Vojta, Phys Rev B 80, 024411 (2009).
  * Ribeiro and Lee (2011) P. Ribeiro and P. A. Lee, Phys Rev B 83, 235119 (2011).
  * Vojta  _et al._ (2016) M. Vojta, A. K. Mitchell, and F. Zschocke, Phys Rev Lett 117, 037202 (2016).
  * Das  _et al._ (2016) S. D. Das, K. Dhochak, and V. Tripathi, Phys Rev B 94, 024411 (2016).
  * Kolezhuk  _et al._ (2006) A. Kolezhuk, S. Sachdev, R. R. Biswas, and P. Chen, Phys Rev B 74, 165114 (2006).
  * He and Lee (2022) W.-Y. He and P. A. Lee, Phys. Rev. B 105, 195156 (2022).
  * Yamamoto  _et al._ (2018) T. D. Yamamoto, H. Taniguchi, and I. Terasaki, J. Phys.: Condens. Matter 30, 355801 (2018).
  * Gomilšek  _et al._ (2019) M. Gomilšek, R. Žitko, M. Klanjšek, M. Pregelj, C. Baines, Y. Li, Q. M. Zhang, and A. Zorko, Nat Phys 15, 754 (2019).
  * Chen  _et al._ (2022) Y. Chen, W.-Y. He, W. Ruan, J. Hwang, S. Tang, R. L. Lee, M. Wu, T. Zhu, C. Zhang, H. Ryu, F. Wang, S. G. Louie, Z.-X. Shen, S.-K. Mo, P. A. Lee, and M. F. Crommie, arXiv e-prints , arXiv:2202.07224 (2022), arXiv:2202.07224 [cond-mat.str-el] .
  * Withoff and Fradkin (1990) D. Withoff and E. Fradkin, Phys. Rev. Lett. 64, 1835 (1990).
  * Satori  _et al._ (1992) K. Satori, H. Shiba, O. Sakai, and Y. Shimizu, J. Phys. Soc. Jpn. 61, 3239 (1992).
  * Saso (1992) T. Saso, J. Phys. Soc. Jpn. 61, 3439 (1992).
  * Itoh (1993) Y. Itoh, J. Phys. Soc. Jpn. 62, 2184 (1993).
  * Takegahara  _et al._ (1993) K. Takegahara, Y. Shimizu, N. Goto, and O. Sakai, Physica B 186-188, 381 (1993).
  * Chen and Jayaprakash (1998) K. Chen and C. Jayaprakash, Phys. Rev. B 57, 5225 (1998).
  * Lu (2018) Y.-M. Lu, Phys Rev B 97 (2018), 10.1103/physrevb.97.094422.
  * Qi and Cheng (2018) Y. Qi and M. Cheng, Phys. Rev. B 97, 115138 (2018).
  * Yang and Wang (2016) X. Yang and F. Wang, Phys. Rev. B 94, 035160 (2016).
  * Wang and Vishwanath (2006) F. Wang and A. Vishwanath, Phys. Rev. B 74, 174423 (2006).
  * Lu  _et al._ (2011) Y.-M. Lu, Y. Ran, and P. A. Lee, Phys. Rev. B 83, 224413 (2011).
  * Lu (2016) Y.-M. Lu, Phys. Rev. B 93, 165113 (2016).
  * Kitaev (2003) A. Kitaev, Ann Phys-new York 303, 2–30 (2003).
  * Zaletel and Vishwanath (2015) M. P. Zaletel and A. Vishwanath, Phys. Rev. Lett. 114, 077201 (2015).
  * Qi  _et al._ (2015) Y. Qi, M. Cheng, and C. Fang, arXiv e-prints , arXiv:1509.02927 (2015), arXiv:1509.02927 [cond-mat.str-el] .
  * Lu  _et al._ (2017) Y.-M. Lu, G. Y. Cho, and A. Vishwanath, Phys Rev B 96, 205150 (2017).
  * (51) The details can be found in the supplemental materials .
  * Ran and Wen (2006) Y. Ran and X.-G. Wen, (2006), arXiv:cond-mat/0609620 [cond-mat.str-el] .
  * Kitaev (2006) A. Kitaev, Ann. Phys. 321, 2 (2006).
  * Fu  _et al._ (2015) M. Fu, T. Imai, T.-H. Han, and Y. S. Lee, Science 350, 655 (2015).
  * Qi and Fu (2015) Y. Qi and L. Fu, Phys. Rev. B 91, 100401 (2015).
  * Marston and Affleck (1989) J. B. Marston and I. Affleck, Phys. Rev. B 39, 11538 (1989).

Supplemental materials

##  Appendix A A brief review of symmetry fractionalization in $Z_{2}$ spin
liquids

Below we provide a brief summary on the classification of gapped symmetric
$Z_{2}$ spin liquids on the square lattice, with a single spin-$1/2$ on each
site.

Most generally, given a topological order described by a unitary modular
tensor category $\mathcal{C}$, a $G$-symmetry enriched topological (SET) phase
is mathematically classified by symmetry action $\rho:~{}G\rightarrow
Aut(\mathcal{C})$, and twisted group cohomology
$\mathcal{H}^{2}_{\rho}(G,\mathcal{A})$ where $\mathcal{A}$ represents all
Abelian anyons in category $\mathcal{C}$Essin and Hermele (2013); Barkeshli
_et al._ (2019); Tarantino  _et al._ (2016). In the special case where
symmetries do not permute different anyons in $\mathcal{C}$, different SET
phases are classified by group cohomology $\mathcal{H}^{2}(G,\mathcal{A})$,
which characterizes the symmetry fractionalization in the SET phaseEssin and
Hermele (2013). Physically, the localized symmetry $\\{U_{g}|g\in G\\}$ on the
anyons form a projective representationWen (2002) of the symmetry group $G$:

$\displaystyle
U_{g}U_{h}=\omega_{g,h}U_{g,h},~{}~{}~{}\omega_{g,h}\in\mathcal{A}.$ (7)

Acting on anyon $a$, the product of $U_{g}$ and $U_{h}$ actions differ from
$U_{gh}$ action by a $U(1)$ phase given by the braiding phase between anyon
$a$ and Abelion anyon $\omega_{g,h}$:

$\displaystyle\omega_{g,h}^{a}=\frac{S_{a,\omega_{g,h}}}{S_{a,0}}\in U(1)$ (8)

where $S_{a,b}$ is the modular S matrix of topological order $\mathcal{C}$.

###  A.1 Square lattice

As illustrated in Fig. 4, the space group symmetry $p4m$ of the square lattice
is generated by translation $T_{1}$, site-centered 4-fold rotation $C_{4}$ and
mirror reflection $\sigma$. All other symmetry elements can be generated by
them, for example in Fig. 4

$\displaystyle T_{2}=C_{4}T_{1}C_{4}^{-1},$ (9) $\displaystyle
M_{y}=T_{2}\sigma,$ (10) $\displaystyle R_{xy}=\sigma C_{4},$ (11)
$\displaystyle M_{x}=C_{4}M_{y}C_{4}^{-1}.$ (12)

The full symmetry group is the product of space group symmetry $p4m$, time
reversal symmetry $Z_{2}^{\mathcal{T}}$, and the spin rotational symmetry
$SO(3)$. On a square lattice of spin-$1/2$’s, the symmetry action of spin
rotations is fixed by the symmetry and Hilbert space in a relatively simple
fashion: $e$ and $\epsilon$ particles each carry spin-$1/2$, while $m$
particles carry spin-$0$. We consider the rest of the symmetry group
$G=p4m\times Z_{2}^{\mathcal{T}}$, where the symmetry fractionalization
classes of $Z_{2}$ spin liquids are classified by 2nd group cohomology

$\displaystyle\mathcal{H}^{2}(G,\mathcal{A}=Z_{2}\times
Z_{2})=\big{[}\mathcal{H}^{2}(G,Z_{2})\big{]}^{2}$ (13)

Here $\mathcal{A}$ is the fusion of Abelian anyons in a $Z_{2}$ spin liquid.
Each one of the two $[\omega]\in\mathcal{H}^{2}(G,Z_{2})$ classes can be
viewed as the classification for symmetry fractionalization on a specific
anyon, for example $[\omega^{e}]$ for bosonic spinon $e$ and
$[\omega^{\epsilon}]$ for fermionic spinon $\epsilon$. The associated symmetry
fractionalization class for the vison $m=e\times\epsilon$ can be obtained by
the fusion rule:

$\displaystyle\omega^{m}_{g,h}=\omega^{e}_{g,h}\omega^{\epsilon}_{g,h}\omega_{t}(g,h),~{}~{}~{}\forall~{}g,h\in
G.$ (14)

where $[\omega_{t}(g,h)=\pm 1]$ are a set of twist factors determined by the
symmetry group and the underlying topological order.

Figure 4: Space group symmetries of the square lattice.

As shown in Table 4, the different cohomology classes of
$\mathcal{H}^{2}(G,Z_{2})=(Z_{2})^{10}$ are characterized by 10 different
$Z_{2}$-valued invariantsEssin and Hermele (2013). Naively, there will be 20
independent $Z_{2}$ invariants, 9 for bosonic spinons $[\omega^{e}]$ and 10
for fermionic spinons $[\omega^{\epsilon}]$. Luckily, as we will see below,
the $SO(3)$ spin rotational symmetry and spin-$1/2$ Hilbert space fix 14 of
them and reduce everything to only 6 independent $Z_{2}$ invariantsLu  _et
al._ (2017), summarized in Table 4.

Algebraic identities | $[\omega]\in\mathcal{H}^{2}(G,\mathcal{A})$ |  $\omega^{e}$Yang and Wang (2016) |  $\omega^{\epsilon}$Wen (2002) |  $\omega^{m}$Qi  _et al._ (2015) | Twist factor $\omega_{t}$Essin and Hermele (2013); Lu  _et al._ (2017)
---|---|---|---|---|---
$\mathcal{T}^{2}$ | $\omega_{\mathcal{T},\mathcal{T}}$ | $-1$ | $-1$ | 1 | 1
$T_{1}T_{2}T_{1}^{-1}T_{2}^{-1}$ | $\frac{\omega_{T_{1},T_{2}}}{\omega_{T_{2},T_{1}}}$ | $(-1)^{p_{1}}$ | $\eta_{xy}$ | -1 | 1
$\sigma^{2}$ | $\omega_{\sigma,\sigma}$ | $(-1)^{p_{4}}$ | $\eta_{\sigma}$ | 1 | -1
$(\sigma\mathcal{T})^{2}$ | $\omega_{\sigma\mathcal{T},\sigma\mathcal{T}}$ | $-1$ | $-\eta_{\sigma}\eta_{\sigma\mathcal{T}}\equiv-1$ | 1 | 1
$(M_{y})^{2}$ | $\omega_{M_{y},M_{y}}$ | $(-1)^{p_{3}+p_{4}}$ | $\eta_{\sigma}\eta_{xpx}$ | 1 | -1
$(M_{y}\mathcal{T})^{2}$ | $\omega_{M_{y}\mathcal{T},M_{y}\mathcal{T}}$ | $(-1)^{p_{3}+p_{8}+1}$ | $-\eta_{t}\eta_{xpx}\eta_{\sigma}\eta_{\sigma\mathcal{T}}\equiv-\eta_{t}\eta_{xpx}$ | 1 | 1
$(C_{4}\sigma)^{2}$ | $\omega_{C_{4}\sigma,C_{4}\sigma}$ | $(-1)^{p_{4}+p_{7}}$ | $\eta_{\sigma}\eta_{\sigma C_{4}}$ | 1 | -1
$(C_{4}\sigma\mathcal{T})^{2}$ | $\omega_{C_{4}\sigma\mathcal{T},C_{4}\sigma\mathcal{T}}$ | $-1$ | $-\eta_{\sigma}\eta_{\sigma\mathcal{T}}\eta_{\sigma C_{4}}\eta_{C_{4}\mathcal{T}}\equiv-1$ | 1 | 1
$M_{x}\sigma M_{x}^{-1}\sigma^{-1}$ | $\frac{\omega_{M_{x},\sigma}}{\omega_{\sigma,M_{x}}}$ | $(-1)^{p_{2}}$ | $\eta_{xpy}$ | 1 | -1
$(C_{4}^{2}\mathcal{T})^{2}$ | $\omega_{(C_{4})^{2}\mathcal{T},(C_{4})^{2}\mathcal{T}}$ | $-1$ | $-\eta_{C_{4}}\equiv-1$ | -1 | -1
Table 4: The classification of gapped $Z_{2}$ spin liquids of spin-$1/2$’s on
the square lattice, characterized by 6 independent $Z_{2}$-valued
invariantsYang and Wang (2016); Lu (2018), and their realizations in the
Schwinger bosonYang and Wang (2016) and Abrikosov fermionWen (2002)
representations.

First, similar to the $SO(3)$ spin rotational symmetry, the time reversal
symmetry action is fixed by the Hilbert space: $e$ and $\epsilon$ are Kramers
doublets with $\mathcal{T}^{2}=-1$, while each $m$ is a spin-0 particle with
$\mathcal{T}^{2}=+1$Essin and Hermele (2013). Secondly, the $SO(3)$ spin
rotational symmetry completely fix the vison fractionalization class
$[\omega^{m}]$Qi  _et al._ (2015), as shown in Table 4. This means only the
fermionic spinon fractionalization class $[\omega^{\epsilon}]$ is labeled by
independent invariants, and it fully determines the bosonic spinon
fractionalization class $[\omega^{e}]$ by fusion rule (14). Finally, for a
gapped $Z_{2}$ spin liquid, the spin-$1/2$ Hilbert space per site further
fixes $(M\mathcal{T})^{2}=\omega_{M\mathcal{T},M\mathcal{T}}\equiv-1$ on
spinons $e$ and $\epsilon$, as long as the reflection plane $M$ crosses an odd
number of sites per unit cellLu (2018); Qi and Cheng (2018). A similar
conclusion can be made for a 2-fold rotation (or inversion) $I_{s}$ centered
on a lattice site (with a spin-$1/2$ on it):
$(I_{s}\mathcal{T})^{2}=\omega_{I_{s}\mathcal{T},I_{s}\mathcal{T}}\equiv-1$
for an arbitrary anyon carrying spin-$1/2$Lu (2018); Qi and Cheng (2018). For
fermionic spinons, this means $\eta_{\sigma}\equiv\eta_{\sigma\mathcal{T}}$,
$\eta_{\sigma C_{4}}\equiv\eta_{C_{4}\mathcal{T}}$, and $\eta_{C_{4}}\equiv
1$. Therefore we reach the 6 independent $Z_{2}$ invariantsLu (2018)
summarized in Table 1.

Below we make a few comments on the relation between different algebraic
identities in the cohomology class, to elaborate on how Table 1 is obtained
from Table 4. First, for elements $\omega\in\mathcal{H}^{2}(G,Z_{2})$, we have

$\displaystyle\omega_{(C_{4})^{2}\mathcal{T},(C_{4})^{2}\mathcal{T}}\equiv\omega_{(C_{4})^{2},(C_{4})^{2}}\omega_{\mathcal{T},\mathcal{T}}$
(15)

hence the identity $(C_{4})^{4}=e$ does not give rise to an extra invariant
beyond Table 4. For spinons, this also indicates that

$\displaystyle\omega^{e,\epsilon}_{(C_{4})^{2},(C_{4})^{2}}\equiv 1$ (16)

Another algebraic relation is that

$\displaystyle\frac{\omega_{M_{x},M_{y}}}{\omega_{M_{y},M_{x}}}=\frac{\omega_{T_{1},T_{2}}}{\omega_{T_{2},T_{1}}}\omega_{(C_{4})^{2},(C_{4})^{2}}$
(17)

and hence $M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}$ does not give rise to a new
invariant. For spinons with
$\omega^{e,\epsilon}_{(C_{4})^{2},(C_{4})^{2}}\equiv 1$, we therefore have

$\displaystyle\frac{\omega^{e,\epsilon}_{M_{x},M_{y}}}{\omega^{e,\epsilon}_{M_{y},M_{x}}}=\frac{\omega^{e,\epsilon}_{T_{1},T_{2}}}{\omega^{e,\epsilon}_{T_{2},T_{1}}}$
(18)

as reported in Table 1.

Figure 5: Space group symmetries of the kagome and triangular lattices.

###  A.2 Kaogme and triangular lattices

Algebraic identities | $[\omega]\in\mathcal{H}^{2}(G,\mathcal{A})$ |  $\omega^{e}$Wang and Vishwanath (2006) |  $\omega^{\epsilon}$Lu  _et al._ (2011) |  $\omega^{m}$Qi  _et al._ (2015) | Twist factor $\omega_{t}$Qi and Fu (2015); Lu  _et al._ (2017)
---|---|---|---|---|---
$\mathcal{T}^{2}$ | $\omega_{\mathcal{T},\mathcal{T}}$ | $-1$ | $-1$ | 1 | 1
$T_{1}T_{2}T_{1}^{-1}T_{2}^{-1}$ | $\frac{\omega_{T_{1},T_{2}}}{\omega_{T_{2},T_{1}}}$ | $(-1)^{p_{1}}$ | $\eta_{12}$ | -1 | 1
$(R_{y})^{2}$ | $\omega_{R_{y},R_{y}}$ | $(-1)^{p_{2}}$ | $\eta_{\sigma}\eta_{\sigma C_{6}}$ | 1 | -1
$(R_{y}\mathcal{T})^{2}$ | $\omega_{R_{y}\mathcal{T},R_{y}\mathcal{T}}$ | $-1$ | $-\eta_{\sigma}\eta_{\sigma\mathcal{T}}\eta_{\sigma C_{6}}\eta_{C_{6}\mathcal{T}}\equiv-1$ | 1 | 1
$(R_{x})^{2}$ | $\omega_{R_{x},R_{x}}$ | $(-1)^{p_{2}+p_{3}}$ | $\eta_{\sigma}$ | 1 | -1
$(R_{x}\mathcal{T})^{2}$ | $\omega_{R_{x}\mathcal{T},R_{x}\mathcal{T}}$ | $-1$ | $-\eta_{\sigma}\eta_{\sigma\mathcal{T}}\equiv-1$ | 1 | 1
$M_{x}R_{y}M_{x}^{-1}R_{y}^{-1}$ | $\frac{\omega_{M_{x},R_{y}}}{\omega_{R_{y},M_{x}}}$ | $1$ | $\eta_{12}\eta_{C_{6}}\eta_{\sigma C_{6}}\equiv 1$ | -1 | -1
Table 5: The classification of gapped $Z_{2}$ spin liquids of spin-$1/2$’s on
the kagome lattice, characterized by 3 independent $Z_{2}$-valued invariantsQi
and Cheng (2018); Lu (2018), and their realizations in the Schwinger bosonWang
and Vishwanath (2006) and Abrikosov fermionLu  _et al._ (2011)
representations.

Similar to the case of spin-$1/2$’s on the square lattice, one can also
classify symmetric $Z_{2}$ spin liquids of spin-$1/2$’s on the kagome and
triangular lattices. Here we briefly review the classification results
reported previously in Ref.Qi and Cheng (2018); Lu (2018).

In the case of spin-$1/2$’s on the kagome lattice, the space group is $p6mm$,
generated by translation $T_{1}$, hexagon-centered 6-fold rotation $C_{6}$ and
mirror reflection $R_{y}$. The other symmetries shown in Fig. 5 can be
generated as follows:

$\displaystyle T_{2}=C_{6}T_{1}C_{6}^{-1},$ (19) $\displaystyle
R_{x}=(C_{6})^{3}R_{y},$ (20) $\displaystyle M_{x}=T_{1}R_{x}.$ (21)

We consider the symmetry group $G=p6mm\times Z_{2}^{\mathcal{T}}$ in addition
to $SO(3)$ spin rotational symmetry. The 2nd group cohomology has
$\mathcal{H}^{2}(G,Z_{2})=(Z_{2})^{7}$ different classes (see the 7 rows in
Table 5), naively leading to $(Z_{2})^{14}$ different fractionalization
classes. For spin-$1/2$’s on the kagome lattice with the $SO(3)$ spin
rotational symmetries, again we can fix many of the cohomology classes and
reduce the final classification to only $2^{3}=8$ distinct symmetric $Z_{2}$
spin liquidsLu (2018); Qi and Cheng (2018). They are characterized by three
$Z_{2}$-valued invariants, as summarized in TABLE 5.

In particular, since both reflection planes $R_{x,y}$ (see Fig. 5) cross one
site per unit cell, we have

$\displaystyle\omega^{e,\epsilon}_{R_{x}\mathcal{T},R_{x}\mathcal{T}}=\omega^{e,\epsilon}_{R_{y}\mathcal{T},R_{y}\mathcal{T}}=-1$
(22)

for bosonic spinon $e$ and fermionic spinon $\epsilon$. Meanwhile, since
reflection planes $M_{x}=T_{1}R_{x}$ and $R_{y}$ intersect at one site, we
have

$\displaystyle\frac{\omega^{e,\epsilon}_{M_{x},R_{y}}}{\omega^{e,\epsilon}_{R_{y},M_{x}}}=+1$
(23)

for both $e$ and $m$. As a result, for the two reflections plane $R_{x}$ and
$R_{y}$, using the following identity

$\displaystyle\frac{\omega_{M_{x},R_{y}}}{\omega_{R_{y},M_{x}}}\equiv\frac{\omega_{M_{x},R_{y}}}{\omega_{R_{y},M_{x}}}\cdot\frac{\omega_{T_{1},T_{2}}}{\omega_{T_{2},T_{1}}}$
(24)

we have

$\displaystyle\frac{\omega^{e,\epsilon}_{M_{x},R_{y}}}{\omega^{e,\epsilon}_{R_{y},M_{x}}}\equiv\frac{\omega^{e,\epsilon}_{T_{1},T_{2}}}{\omega^{e,\epsilon}_{T_{2},T_{1}}}$
(25)

as reported in Table 2.

Algebraic identities | $[\omega]\in\mathcal{H}^{2}(G,\mathcal{A})$ |  $\omega^{e}$Wang and Vishwanath (2006) |  $\omega^{\epsilon}$Lu (2016) |  $\omega^{m}$Qi  _et al._ (2015) | Twist factor $\omega_{t}$Qi and Fu (2015); Lu  _et al._ (2017)
---|---|---|---|---|---
$\mathcal{T}^{2}$ | $\omega_{\mathcal{T},\mathcal{T}}$ | $-1$ | $-1$ | 1 | 1
$T_{1}T_{2}T_{1}^{-1}T_{2}^{-1}$ | $\frac{\omega_{T_{1},T_{2}}}{\omega_{T_{2},T_{1}}}$ | $(-1)^{p_{1}}$ | $\eta_{12}$ | -1 | 1
$\sigma^{2}$ | $\omega_{\sigma,\sigma}$ | $(-1)^{p_{2}}$ | $\eta_{\sigma}$ | 1 | -1
$(\sigma\mathcal{T})^{2}$ | $\omega_{\sigma\mathcal{T},\sigma\mathcal{T}}$ | $-1$ | $-\eta_{\sigma}\eta_{\sigma\mathcal{T}}\equiv-1$ | 1 | 1
$R^{2}$ | $\omega_{R,R}$ | $(-1)^{p_{2}+p_{3}}$ | $\eta_{\sigma C_{6}}$ | 1 | -1
$(R\mathcal{T})^{2}$ | $\omega_{R\mathcal{T},R\mathcal{T}}$ | $-1$ | $-\eta_{\sigma C_{6}}\eta_{\sigma\mathcal{T}}\eta_{C_{6}\mathcal{T}}\equiv-1$ | 1 | 1
$\sigma R\sigma^{-1}R^{-1}$ | $\frac{\omega_{\sigma,R}}{\omega_{R,\sigma}}$ | $1$ | $\eta_{\sigma}\eta_{C_{6}}\eta_{\sigma C_{6}}\equiv 1$ | -1 | -1
Table 6: The classification of gapped $Z_{2}$ spin liquids of spin-$1/2$’s on
the triangular lattice, characterized by 3 independent $Z_{2}$-valued
invariantsQi and Cheng (2018); Lu (2018), and their realizations in the
Schwinger bosonWang and Vishwanath (2006) and Abrikosov fermionLu (2016)
representations.

In the case of the triangular lattice, the space group is still $p6mm$ with
the same $G=p6mm\times Z_{2}^{\mathcal{T}}$. The spin-$1/2$’s however are
located at different Wyckoff sites compared to the kagome lattice. As shown in
Fig. 5, the space group is generated by translation $T_{1}$, site-centered
6-fold rotation $C_{6}$ and mirror reflection $R$. The other symmetries in
Fig. 5 can be generated as

$\displaystyle\sigma=(C_{6})^{3}R,$ (26) $\displaystyle T_{2}=\sigma
T_{1}\sigma^{-1},$ (27) $\displaystyle M_{y}=T_{1}T_{2}R.$ (28)

Completely in parallel to the kagome lattice case, one can show the
$\mathcal{H}^{2}(G,\mathcal{A})=(Z_{2})^{14}$ cohomology classes summarized in
Table 6 can be reduced to $2^{3}$ distinct symmetric $Z_{2}$ spin liquids, by
requiring a spectrum gap and $SO(3)$ spin rotational symmetryQi and Cheng
(2018); Lu (2018). They are listed in Table 3. In particular, following the
group cohomology identity:

$\displaystyle\frac{\omega_{\sigma,M_{y}}}{\omega_{M_{y},\sigma}}\equiv\frac{\omega_{\sigma,R}}{\omega_{R,\sigma}}\cdot\frac{\omega_{T_{1},T_{2}}}{\omega_{T_{2},T_{1}}}$
(29)

we have

$\displaystyle\frac{\omega^{e,\epsilon}_{\sigma,M_{y}}}{\omega^{e,\epsilon}_{M_{y},\sigma}}\equiv\frac{\omega^{e,\epsilon}_{T_{1},T_{2}}}{\omega^{e,\epsilon}_{T_{2},T_{1}}}$
(30)

as reported in Table 3.

##  Appendix B Exactly solvable model

We use an exactly solvable model with spin-$1/2$ spinon excitations to
illustrate the symmetry protected non-Kramers doublet localized at the
impurity, in the Kondo screening phase. The model is based on Kitaev’s toric
code on the square lattice, with four types of topological excitations,
including the trivial sector $1$, bosonic spinon $e$, vison $m$, and fermionic
spinon $\epsilon=e\times m$. To incorporate the $SU(2)$ symmetry, we revise
the toric code to construct a $Z_{2}$ spin liquid where $e,\epsilon$
excitations each carry spin-$1/2$.

Figure 6: An exactly solvable model for a spin-$1/2$ impurity in a $Z_{2}$
spin liquid with $SU(2)$ symmetry. The squares on link centers are qubits in
the toric codeKitaev (2003). Each black dot on a star/vertex represents a
Hilbert space of spin $0\oplus 1/2$. The red dot denotes the spin-$1/2$
impurity with a site symmetry $G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$, which
is coupled to four neighboring spin-$1/2$’s.

The Hilbert space of our model consists of two parts in the bulk, as shown in
Fig. 6. Each squares on link center represents a qubit interacting with each
other as in the toric codeKitaev (2003). Each black dots on the vertices
stands for a three dimensional linear space constructed as spin $0\oplus 1/2$.
This Hilbert space allows us to attach a spin $1/2$ to each $e$ particle on
the vertex. In the ground state, every vertex spin is in the spin-$0$ state,
but if there is a $e$ excitation on a vertex, the vertex state is switched
into spin-$1/2$. The red dot represents a spin-$1/2$ impurity interacting with
four nearest spin-$1/2$’s.

To be more specific, the model is

$\hat{H}_{bulk}=-\sum_{s}A_{s}-\sum_{p}B_{p}-\sum_{s}\Delta\big{(}\frac{(A_{s}+1)}{2}P_{s}(S=0)+\frac{(1-A_{s})}{2}P_{s}(S=1/2)\big{)}$
(31)

where $s$ means star (or vertex) and $p$ means plaquette on the square lattice
of solid lines. $A_{s}$ and $B_{p}$ are corresponding star term and plaquette
term in the toric codeKitaev (2003):

$A_{s}=\prod_{j\in s}\sigma_{j}^{x},\quad B_{p}=\prod_{j\in p}\sigma_{j}^{z}$
(32)

The second part of Hamiltonian describes the coupling between spins and the
toric code model. We assume $\Delta\gg 1$. $P_{s}(S=0)$ means projection
operator on vertex $s$ into $S=0$ subspace if $A_{s}=1$ or no $e$ excitation.
Otherwise, if $A_{s}=-1$, the $S=1/2$ subspace has a lower energy because of
the $\Delta$ term. Note the difference between capital $P,S$ and $p,s$. They
have very different meanings here.

The $SU(2)$ symmetry is implemented in the Hilbert space of spin $0\oplus 1/2$
on each vertex. The model shares the same ground state as toric code where all
vertex spins are in the spin-$0$ state. If we choose a very large $\Delta\gg
1$, the low energy excitations are spin-0 $m$ particles, and $e,\epsilon$
particles each carry spin-$1/2$. This model preserves two mirror symmetries
$M_{x},M_{y}$ shown as dotted line in Fig. 6. The symmetry fractionalization
class regarding the site symmetry $G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$ is
given byEssin and Hermele (2013)

$\displaystyle\frac{\omega^{e}_{M_{x},M_{y}}}{\omega^{e}_{M_{y},M_{x}}}=\frac{\omega^{m}_{M_{x},M_{y}}}{\omega^{m}_{M_{y},M_{x}}}=+1,~{}~{}~{}\frac{\omega^{\epsilon}_{M_{x},M_{y}}}{\omega^{\epsilon}_{M_{y},M_{x}}}=-1.$
(33)

In other words, only the fermionic spinon $\epsilon$ has a nontrivial
fractionalization class associated with two mirror symmetries $M_{x,y}$.

The coupling between impurity and bulk also preserves $SU(2)$ symmetry. To
start, we couple the spin-$1/2$ impurity to the four nearly spin-$1/2$ degrees
of freedom in the bulk with antiferromagnetic Heisenberg terms:

$H_{imp}=J\sum_{\langle i,imp\rangle}\vec{S}_{i}\cdot\vec{S}_{imp}$ (34)

where $i$ is summed over four spins closest to the impurity. When $J$ is
positive, impurity couples to bulk spins antiferromagnetically. The full
Hamiltonian is

$\hat{H}=\hat{H}_{bulk}+\hat{H}_{imp}$ (35)

Now we study the ground state with an impurity. When $J$ is small, the
impurity spin remains an unscreened spin-$1/2$, and the decoupled bulk stays
in the original ground state as the toric code. This is the unscreened phase.
When $J>8/9$, the system enters into a new phase where $H_{imp}$ is minimized.
Since the impurity attracts four spin-$1/2$ $e$ particles around it, this is
an overscreened phase.

$|G_{\text{overscreened}}\rangle=|e_{1},e_{2},e_{3},e_{4}\rangle\otimes|S_{tot}=3/2,S_{4}=2,S_{imp}=1/2\rangle$
(36)

Here, $e_{i},i=1\cdots 4$ means four e particles on vertices $1\cdots 4$.
There are two good quantum numbers to label the ground state: $S_{4}$ denotes
the total spin of four spin $1/2$, $S_{tot}$ means the total spin of the whole
system.

In order to achieve a Kondo screening phase with a paramagnetic ground state,
we add another coupling to the impurity Hamiltonian:

$H_{imp}^{\prime}=H_{imp}+E_{c}(A_{1}+A_{2}+A_{3}+A_{4}-3)^{2}$ (37)

where $\hat{A}_{i}$ denotes the vertex term (32) on vertex $i$. The second
term can be viewed as a Coulomb repulsion of strength $E_{c}>0$ for spinons
near the impurity site, since $(A_{1}+A_{2}+A_{3}+A_{4}-3)^{2}$ is minimized
when there is no $e$ particles or only one $e$ particle. Let’s assume $E\gg
1,J$. When $J>4/3$, in the ground state the impurity spin-$1/2$ favors to form
a singlet bound state with one $e$ particle from the bulk, which also brings a
spin-$1/2$ at one neighboring site of the impurity.

$|G_{i}\rangle=|e_{i}\rangle\otimes(|\uparrow_{i}\downarrow_{imp}\rangle-|\downarrow_{i}\uparrow_{imp}\rangle)$
(38)

Here $i$ can be $1\cdots 4$. $|\uparrow_{i}\rangle$ means the vertex $i$ is in
the spin-up state in the 3-dimensional space. The ground state is four fold
degenerate. In particular, if we choose the following gauge

$\displaystyle|e_{2}\rangle=\sigma^{z}_{12}|e_{1}\rangle,~{}~{}~{}|e_{3}\rangle=\sigma^{z}_{12}\sigma^{z}_{23}|e_{1}\rangle,~{}~{}~{}|e_{4}\rangle=\sigma^{z}_{14}|e_{1}\rangle$
(39)

the four degenerate states transform under crystal symmetries as follows:

$\displaystyle|G_{i}\rangle\overset{M_{x}}{\longrightarrow}\sum_{j}U(M_{x})_{i,j}|G_{j}\rangle,~{}~{}~{}U({M_{x}})=\begin{pmatrix}0&0&1&0\\\
0&1&0&0\\\ 1&0&0&0\\\
0&0&0&B_{p}[1234]\end{pmatrix}=U^{-1}(C_{4})U(M_{y})U(C_{4}),$ (40)
$\displaystyle|G_{i}\rangle\overset{M_{y}}{\longrightarrow}\sum_{j}U(M_{y})_{i,j}|G_{j}\rangle,~{}~{}~{}U({M_{y}})=\begin{pmatrix}1&0&0&0\\\
0&0&0&1\\\ 0&0&B_{p}[1234]&0\\\ 0&1&0&0\end{pmatrix},$ (41)
$\displaystyle|G_{i}\rangle\overset{C_{4}}{\longrightarrow}\sum_{j}U(C_{4})_{i,j}|G_{j}\rangle,~{}~{}~{}U({C_{4}})=\begin{pmatrix}0&1&0&0\\\
0&0&1&0\\\ 0&0&0&B_{p}[1234]\\\ 1&0&0&0\end{pmatrix}$ (42)

where $C_{4}$ is the 4-fold rotation around the plaquette center (i.e. the
impurity site). We have defined the $Z_{2}$-valued variable:

$\displaystyle B_{p}[1234]\equiv\langle\hat{B}_{p}[1234]\rangle=\langle
0|\sigma^{z}_{12}\sigma^{z}_{23}\sigma^{z}_{34}\sigma^{z}_{41}|0\rangle=\pm 1$
(43)

where $|0\rangle$ is the ground state of the bulk model (31), and
$\hat{B}_{p}[1234]$ is the plaquette term for the plaquette including sites
$1,2,3,4$. It is instructive to work in the basis of $C_{4}$eigenstates and we
shall follow this strategy in the following calculations.

In the case of $B_{p}[1234]=+1$, the $C_{4}$ eigenstate $|n\rangle$ with
eigenvalue $e^{in\pi/2}$ is given by

$\displaystyle|n=0,1,2,3\rangle=\frac{1}{2}\sum_{j=1}^{4}e^{in(j-1)\pi/2}|G_{j}\rangle$
(44)

In this case both $|n=0\rangle$ and $|n=2\rangle$ states are invariant under
$M_{x,y}$ operations and hence can be the unique ground state of the system.
For example, to select $|n=0\rangle=\sum_{j}|G_{j}\rangle/4$ as the unique
ground state, we can add the following perturbation

$\displaystyle\Delta H_{imp}=-V\sum_{i}|G_{i}\rangle\langle G_{i+1}|+h.c.$
$\displaystyle=-V\sum_{i=1}^{4}\sigma^{z}_{i,i+1}\otimes\big{[}(|\uparrow_{i}\downarrow_{imp}0_{i+1}\rangle-|\downarrow_{i}\uparrow_{imp}0_{i+1}\rangle)(\langle\uparrow_{i+1}\downarrow_{imp}0_{i}|-\langle\downarrow_{i+1}\uparrow_{imp}0_{i}|)+h.c.\big{]}$
(45)

which is symmetric w.r.t. all site symmetries at the impurity site. In the
presence of this term, in the Kondo screened phase, the system will have a
unique paramagnetic ground state with no protected degeneracy at the impurity
site.

In the case of $B_{p}[1234]=-1$ which implies one $m$ particle present in the
plaquette, on the other hand, the anyon screening the spin-$1/2$ impurity is
not $e$, but instead a fermionic spinon $\epsilon=e\times m$. This can be
achieved by adding one extra term to (37)

$\displaystyle
H^{\prime\prime}_{imp}=H^{\prime}_{imp}+\Delta_{\epsilon}\hat{B}_{p}[1234]$
(46)

When $\Delta_{\epsilon}>1$, the ground state of the impurity Hamiltonian would
favor the low energy fermionic spinon $\epsilon$ over bosonic spinon $e$, and
therefore screen the impurity spin-$1/2$ by one $\epsilon$ particle. We can
similarly label the 4 ground states in the low-energy manifold of the full
Hamiltonian $H_{bulk}+H^{\prime\prime}_{imp}$ as

$|G_{i}\rangle=|\epsilon_{i}\rangle\otimes(|\uparrow_{i}\downarrow_{imp}\rangle-|\downarrow_{i}\uparrow_{imp}\rangle)$
(47)

Choosing the same gauge, their symmetry transformation laws still follow Eqs.
(40)-(42). If we label them in the basis of $C_{4}$ eigenstates

$\displaystyle C_{4}|\pm 1\rangle=e^{\pm i\pi/4}|\pm 1\rangle,~{}~{}~{}|\pm
1\rangle\equiv\frac{e^{\pm i\pi/4}|G_{2}\rangle\pm i|G_{2}\rangle+e^{\pm
i3\pi/4}|G_{3}\rangle+|G_{4}\rangle}{2};$ (48) $\displaystyle C_{4}|\pm
3\rangle=e^{\pm i3\pi/4}|\pm 3\rangle,~{}~{}~{}|\pm 3\rangle\equiv\frac{e^{\pm
i3\pi/4}|G_{2}\rangle\mp i|G_{2}\rangle+e^{\pm
i\pi/4}|G_{3}\rangle+|G_{4}\rangle}{2}.$ (49)

Each of the two doublets, i.e. $|\pm 1\rangle$ (or $|\pm 3\rangle$) forms a
2-dimensional irreducible projective representation of the site symmetry group
$G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$. In particular, in the doublet $|\pm
1\rangle$ with $C_{4}$ eigenvalues $e^{\pm i\pi/4}$, the two mirror actions
are represented by Pauli matrices:

$\displaystyle\langle\pm 1|U(M_{x})|\pm 1\rangle=-\sigma_{x},~{}~{}\langle\pm
1|U(M_{y})|\pm 1\rangle=\sigma_{y}.$ (50)

Similarly, in the doublet $|\pm 3\rangle$ with $C_{4}$ eigenvalues $e^{\pm
i3\pi/4}$, the two mirror actions are represented by Pauli matrices:

$\displaystyle\langle\pm 3|U(M_{x})|\pm 3\rangle=-\sigma_{x},~{}~{}\langle\pm
3|U(M_{y})|\pm 3\rangle=-\sigma_{y}.$ (51)

Indeed $M_{x}$ and $M_{y}$ operations anticommute with each other in each
doublet pair, consistent with the fractionalization class (33) of fermionic
spinon $\epsilon$ in the toric code. As a result, any symmetry-preserving
local Hamiltonian cannot lift this two-fold degeneracy, which is protected by
two mirror symmetries $M_{x,y}$.

##  Appendix C Large-$N$ theory and thermodynamics

To further illustrate the effect of nontrivial symmetry fractionalization
classes on Kondo impurity problem in $Z_{2}$ spin liquids, from the
perspective of parton constructionLu (2018), we construct a large-$N$ model
for the Kondo impurity problem in a symmetric $Z_{2}$ spin liquid. We then
carry out large-$N$ mean-field calculations to obtain the phase diagram, and
the temperature dependence of thermodynamic quantities including the specific
heat (and hence entropy) and uniform susceptibility.

We start from a solvable model similar to Kitaev’s honeycomb model Kitaev
(2006) for the bulk $Z_{2}$ spin liquid, where the exact spectrum of the
system is given by parton mean-field ansatz under different $Z_{2}$ flux
configurations. In these models, the fractionalization class
$[\omega^{\epsilon}]\in\mathcal{H}^{2}(G,z_{2})$ of fermionic spinons
$\epsilon$ manifest itself as a projective symmetry group (PSG)Wen (2002) of
the symmetry group $G$ carried by fermionic partons. This model is defined on
the square lattice and can support gapped $Z_{2}$ spin liquid ground states
with either nontrivial or trivial PSGs of fermionic partons (or spinons). The
model is not exactly solvable anymore after the impurity spin is added and
coupled to the bulk, but we can still obtain analytical results in large-$N$
limit.

###  C.1 Single impurity in the $SU(2)$ parton construction

The parton construction of quantum spin liquidsWen (2002) provides another
aspect to look at the Kondo effect in $Z_{2}$ spin liquids. We consider a
square lattice with a spin-$1/2$ on each site. In the $SU(2)$ parton (or slave
particle) construction, spin $1/2$ can be represented by a doublet of fermions
$\\{c_{\uparrow},c_{\downarrow}\\}$ on each site with the single occupancy
constraint
$c^{\dagger}_{\uparrow}c_{\uparrow}+c^{\dagger}_{\downarrow}c_{\downarrow}=1$Wen
(2002). In this way, the physical spin state can be obtained by projecting a
many-fermion state into the physical spin Hilbert space, by enforcing the
single-occupancy constraint on each site. Below, we use the $SU(2)$ parton
construction on the square lattice, as shown in Fig.7, to clarify the physical
reason behind the symmetry protected non-Kramers doublet localized at the
impurity site, in the Kondo screening phase.

Figure 7:  (Left) A spin-$1/2$ impurity at a plaquette center on a square
lattice, preserving two mirror symmetries $M_{x,y}$. (Right) The spectrum of
the parton BdG Hamiltonian with a nontrivial PSG (52), where two mirror
symmetries enforce a pair of zero modes on an even by even square lattice.
Among the 4 many-body states obtained by filling the two zero modes, only the
both filled and the both empty states are physical, corresponding to the
unscreened phase with a free spin at the impurity site.

In the parton construction, in order to obtain a symmetric spin state, the
fermion state only needs to be invariant under symmetry up to a gauge
transformation. In other words, the fermions preserve a projetive symmetry
group (PGS), which is an extension of the physical symmetry group by an
invariant gauge group (IGG)Wen (2002). The IGG is $Z_{2}$ in our case of
$Z_{2}$ spin liquids. Specifically, we consider two mirror symmetries
$M_{x},M_{y}$, which always commute with each other in the physical spin
system. However for fermionic spinons with a nontrivial fractionalization
class, they satisfy the following PSG:

$M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=(-)^{\hat{N}_{F}}$ (52)

where $\hat{N}_{F}$ is total fermion number operator. In contrast, for a
$Z_{2}$ spin liquid with a trivial fractionalization class, the fermionic
spinons satisfy $M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=1$.

Due to $SU(2)$ spin rotational symmetry, a quadratic mean-field ansatz for a
symmetric $Z_{2}$ spin liquid has the following formWen (2002)

$H=\sum_{i,j}J_{i,j}\psi_{i}^{\dagger}u_{i,j}\psi_{j}+h.c.$ (53)

where we defined 2-component fermions
$\psi_{i}=(c_{\uparrow},c^{\dagger}_{\downarrow})^{T}$. The particle-hole pair
$\gamma_{\pm E}$ with energy $E$ and $-E$ always appears together in the
spectrum of the parton BdG Hamiltonian in the Nambu basis, as guaranteed by
time reversal symmetry $\mathcal{T}$ or spin rotational symmetry $e^{i\pi
S^{y}}$. Both $\mathcal{T}$ and $e^{i\pi S_{y}}$ symmetries map $\psi$ to
$(c_{\downarrow},-c^{\dagger}_{\uparrow})^{T}=i\tau^{y}\psi_{i}^{\ast}$,
therefore mapping a mode $\gamma_{E}$ to $\gamma_{-E}^{\dagger}$.

Now we place a spin-$1/2$ impurity to one plaquette center (see Fig.7) and
couple it to neighboring spins. Similar to bulk spin-$1/2$’s, the impurity
spin-$1/2$ can also be represented by a doublet of fermions. We consider a
$N\times N$ square lattice with $N=$ even, and one impurity at the center the
torus. Together with the impurity site, overall there is an odd number of
spin-$1/2$s, and $2(N^{2}+1)$ energy levels in the spectrum of the BdG
Hamiltonian.

In the case of a nontrivial fractionalization class with
$M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=(-1)^{\hat{N}_{F}}$, each energy level in the
spectrum of the parton BdG Hamiltonian must be at least 2-fold degenerate,
forming a two dimensional projective representation of mirror symmetries
$Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$. Here we discuss two different scenarios,
corresponding to the unscreened phase and the Kondo screening phase in the
parton construction.

First, if the mean-field ansatz (including impurity) preserves mirror
symmetries $M_{x,y}$ mentioned above, there must be a pair of zero modes in
the parton BdG spectrum. This is because the each energy level must be 2-fold
degenerate due to PSG (52), and therefore any particle-hole symmetric spectrum
must have $N^{2}$ negative energy levels, $N^{2}$ positive energy levels, and
2 degenerate zero modes, as shown in Fig.7. Naively there are four degenerate
many-body ground states, because those two zero modes can be either filled or
empty, as shown in Fig. 7. Both filled or both empty in the Nambu basis
corresponds to unscreened free spin $\uparrow$ or $\downarrow$ state, while
filling only one zero mode leads to paramagnetic ground states. Here due to
the single occupancy constraint on each site, the physical Hilbert space must
have an odd number of fermions on a lattice with an odd number of sites, and
hence an even particle number in the Nambu basis. As a result, only the both
filled and the both empty states are physical, and adiabatically connected to
the free spin state where partons on the impurity site do not couple to bulk
partons in the BdG ansatz. This corresponds to an unscreened phase with a free
spin on the impurity site.

The second scenario is when the parton BdG spectrum is gapped because the
mean-field ansatz around the impurity site spontaneously breaks the mirror
symmetries, even though the bulk spin liquid preserves both mirrors. In this
case, due to lack of mirror symmetries in the full BdG ansatz, there should be
$N^{2}+1$ negative levels, and $N^{2}+1$ positive ones, and the BdG spectrum
is gapped with a unique many-body ground state. This corresponds to the Kondo
screening phase with paramagnetic ground states. Here we need to consider the
two symmetry-breaking mean-field ansatz $|A_{1}\rangle$ and $|A_{2}\rangle$
related by mirror operation $M_{x}$ (or $M_{y}$). As illustrated in Fig. 10.
One can show that due to nontrivial PSG (52), quantum tunneling between these
two states are forbidden by the mirror symmetries, leading to a 2-fold
degeneracy that cannot be lifted by any symmetry-preserving local Hamiltonian.
This exactly corresponds to the mirror symmetry protected non-Kramers doublet
at the impurity site, in the Kondo screening phase. More details will be
discussed later in the self-consistent calculations of the $Sp(2N)$ model.

For $Z_{2}$ spin liquid with a trivial spinon PSG, i.e. when
$M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=1$ for partons, the symmetry fractionalization
of $\epsilon$ particle is trivial. In this case, the BdG spectrum is
generically gapped with $N^{2}+1$ positive and negative levels. The unique
ground state preserves both mirror symmetries and is a paramagnet,
corresponding to a usual Kondo screening phase, with no symmetry protected
degeneracy. This physical picture based on $SU(2)$ parton construction will be
manifested in more detail soon in the large-$N$ $Sp(2N)$ model.

Up to this point, we have only discussed a single spin-$1/2$ impurity on a
lattice with an even number of spin-$1/2$s. In this situation, in fact, the
Kondo screening phase will also excite another delocalized spinon in the bulk,
when one spinon is localized at the impurity site to form a singlet bound
state with the impurity. In a realistic finite size calculation, to achieve a
paramagnetic ground state, we always consider two high-symmetry impurity
sites, separated from each other by a largest distance. In this case, there
will be one non-Kramers doublet at each impurity site, in the Kondo screening
phase.

###  C.2 The $Sp(2N)$ model

The full Hamiltonian of the large $N$ theory contains the bulk $Z_{2}$ spin
liquid and an impurity coupled to bulk.

$H=H_{bulk}+H_{imp}$ (54)

As shown in Fig. 8, on each site of the two-dimensional square lattice, we
define a physical “spin” (more precisely, a qudit) composed of $2N$ fermions
and $4$ Majorana fermions with a total fermion parity constraint shown below.
We have $4$ Majorana operators $\gamma^{\alpha}_{i}$, where $i$ labels lattice
sites and $\alpha=1\cdots 4$ labels four species of Majorana fermions on the
same site. We also have $N$ flavors of spin-$1/2$ complex fermions
$c_{i\alpha}^{a}$, where $i$ labels lattice sites,
$\alpha=\uparrow,\downarrow$ labels spin and $a=1\cdots N$ labels flavor. The
aforementioned constraint is that the total fermion parity on each site is
even:

$\prod_{\alpha=1}^{4}\gamma^{\alpha}(-1)^{\sum_{a=1}^{N}(f_{i\uparrow}^{\dagger}f_{i\uparrow}+f_{i\downarrow}^{\dagger}f_{i\downarrow})}=1$
(55)

Figure 8: An illustration of the solvable $Sp(2N)$ model of a $Z_{2}$ spin
liquid. Each blue circle represents a physical spin composed of four Majorana
operators (small black dots) and $2N$ complex fermions (big black dot). They
live on a square lattice, with the two mirror symmetries shown by dashed
lines. Red dot denotes the impurity and it only couples to the four nearest
spins.

The solvable model is illustrated in Fig. 8. On each site blue circle
represents spin, four small dots represent four Majorana operators and the big
dot at the center represents $2N$ fermions. Coupling between Majorana
operators happens between neighbour sites as shown by solid lines in the
Figure. We label two Majorana fermions which forms a dimer between neighboring
sites $i,j$ as $\gamma_{i}^{\alpha(i,j)}$ and $\gamma_{j}^{\beta(i,j)}$, where
we defined

$\displaystyle\alpha(i,i+\hat{x})$ $\displaystyle=1\quad\beta(i,i+\hat{x})$
$\displaystyle=3$ (56) $\displaystyle\alpha(i,i+\hat{y})$
$\displaystyle=2\quad\beta(i,i+\hat{y})$ $\displaystyle=4$
$\displaystyle\alpha(i,i-\hat{x})$ $\displaystyle=3\quad\beta(i,i-\hat{x})$
$\displaystyle=1$ $\displaystyle\alpha(i,i-\hat{y})$
$\displaystyle=4\quad\beta(i,i-\hat{y})$ $\displaystyle=2$

The solvable Hamiltonian for the bulk $Z_{2}$ spin liquid is

$\displaystyle H_{bulk}=t$ $\displaystyle\sum_{\langle
i,j\rangle}\hat{u}_{i,j}(ic_{i\uparrow}^{a\dagger}c^{a}_{j\uparrow}+ic_{i\downarrow}^{a\dagger}c^{a}_{j\downarrow}+h.c.)$
(57) $\displaystyle+t^{\prime}$ $\displaystyle\sum_{j=i\pm
2\hat{x}+2\hat{y}}\hat{u}_{i,j}(c_{i\uparrow}^{a\dagger}c^{a}_{j\uparrow}+c_{i\downarrow}^{a\dagger}c^{a}_{j\downarrow}+h.c.)$
$\displaystyle+$
$\displaystyle\Delta\sum_{i}c^{a}_{i\uparrow}c^{a}_{i\downarrow}+g\sum_{i}\hat{u}_{i,i+\hat{x}}\hat{u}_{i+\hat{x},i+\hat{x}+\hat{y}}\hat{u}_{i+\hat{x}+\hat{y},i+\hat{y}}\hat{u}_{i+\hat{y},i}$

The first term is the nearest neighbour coupling of physical spins, where
$\hat{u}_{i,j}$ is a conserved quantity (vector potential of the emergent
$Z_{2}$ gauge field) of this Hamiltonian, defined asKitaev (2006)

$u_{i,j}=i\gamma^{\alpha(i,j)}_{i}\gamma^{\beta(i,j)}_{j}$ (58)

Here $\alpha(i,j)$ and $\beta(i,j)$ are two functions of two nearby sites
$i,j$ defined previously. The second term with strength $t^{\prime}$ is
interaction between $4$th nearest neighbours. To construct this interaction
terms we also use $\hat{u}_{i,j}$ for two sites $i,j$ that are not nearest
neighbours. We construct it by multiplying a series of connected nearest
neighbor $u_{\langle i,j\rangle}$s. For example, we define

$u_{i,i+2\hat{x}+2\hat{y}}=u_{i,i+\hat{x}}u_{i+\hat{x},i+2\hat{x}}u_{i+2\hat{x},i+2\hat{x}+\hat{y}}u_{i+2\hat{x}+\hat{y},i+2\hat{x}+2\hat{y}}$
(59)

We can also change all $+$ sign to $-$ sign before $\hat{x}$ or $\hat{y}$ in
the equation above to define $u_{i,i-2\hat{x}+2\hat{y}}$,
$u_{i,i+2\hat{x}-2\hat{y}}$ and $u_{i,i-2\hat{x}-2\hat{y}}$. The third term
represents for the onsite pairing. The fourth term with a positive $g>0$
selects the flux around each plaquette in the ground state. It is assumed to
be much larger than the other 3 terms to make sure that $\pi$ flux has the
lowest energy.

This Hamiltonian preserves time reversal symmetry $\mathcal{T}$. Time reversal
symmetry does not change $\gamma_{i}^{\alpha}$ and it is represented as
$\mathcal{T}=i\sigma^{2}K$ on the fermions, where $\sigma^{2}$ is Pauli matrix
of spin index and $K$ is the complex conjugation operation. Translational
symmetries $T_{x},T_{y}$ and link centered mirror symmetry $M_{x},M_{y}$
(labeled in Fig. 8) act conventionally on the lattice site indices and do not
act on the spin indices.

Since the link variables $\\{\hat{u}_{i,j}\\}$ commute with each other and
Hamiltonian, they can be diagonalized simultanenously and fixed to
$u_{i,j}=\pm 1$. In the absence of vison excitations, in the ground state with
a $\pi$ flux per plaquette, the effective Hamiltonian of the bulk for the
complex fermions is described by

$H_{bulk}=\sum_{i,j}J_{i,j}\psi_{i}^{\dagger}u_{i,j}\psi_{j}+h.c.$ (60)

where

$\displaystyle u_{i,i+\hat{x}}=it\tau^{0}$ (61) $\displaystyle
u_{i,i+\hat{y}}=i(-)^{i_{x}}t\tau^{0}$ $\displaystyle u_{i,i}=\Delta\tau^{1}f$
$\displaystyle u_{i,i\pm 2\hat{x}\pm 2\hat{y}}=t^{\prime}\tau^{3}$

Here, we defined $\psi_{i}=(c_{i,\uparrow},c^{\dagger}_{i,\downarrow})^{T}$Wen
(2002) to write $H_{bulk}$ in a more compact form. Parameters are chosen as
$t=0.5,\Delta=2,t^{\prime}=0.5$ to obtain a gapped $Z_{2}$ spin liquid.

For the quadratic mean-field ansatz of fermionic partons, symmetries are
realized projectively by a PSGWen (2002). More precisely, each symmetry is
implemented together with a gauge transformation on the partons

$\displaystyle G_{\mathcal{T}}(x,y)=(-)^{x+y}$ (62) $\displaystyle
G_{T_{x}}(x,y)=(-)^{y},\quad G_{T_{y}}(x,y)=1$ $\displaystyle
G_{M_{x}}(x,y)=(-)^{x+y},\quad G_{M_{y}}(x,y)=(-)^{y}$

where $(x,y)\in\mathbb{Z}^{2}$ labels a site on the square lattice. It is easy
to verify that

$M_{x}M_{y}M_{x}^{-1}M_{y}^{-1}=-1$ (63)

This implies a nontrivial fractionalization class (or PSG) of the mirror
symmetries $M_{x,y}$ for fermionic spinons $\epsilon$ in the $Z_{2}$ spin
liquid, which plays an important role in the Kondo effect.

For comparison, we can construct another model with a trivial PSG
(fractionalization class) of mirror symmetries for fermionic spinons. Its
Hamiltonian has the same form as Eq. 60 and 66, but with a different bulk mean
field ansatz

$\displaystyle u_{i,i+\hat{x}}=t_{1}\tau^{3}$ (64) $\displaystyle
u_{i,i+\hat{y}}=t_{2}\tau^{2}+t_{1}\tau^{3}$ $\displaystyle
u_{i,i}=\Delta\tau^{1}$ $\displaystyle
u_{i,i\pm\hat{x}\pm\hat{y}}=t^{\prime}\tau^{3}$

Its PSG is trivial for all space group symmetries. The following parameter set
$t_{1}=t_{2}=0.5,t^{\prime}=0.2,\Delta=0.5$ yields a gapped $Z_{2}$ spin
liquid.

Next, we consider adding the impurity spin and how it is coupled to the bulk
spin liquid. We make use of the $Sp(2N)$ model discussed in Ref.Ran and Wen
(2006) with $SU(2)$ gauge structure. This non-Abelian gauge structure makes it
possible for the partons on the impurity site to satisfy PSG relation 63. The
$Sp(2N)$ impurity spin is represented by $N$ flavor spin $1/2$ fermions
$f_{\alpha}^{a}$ with the following constraint

$\displaystyle f^{a\dagger}_{\alpha}f^{a}_{\alpha}=N$ (65) $\displaystyle
f^{a}_{\alpha}f^{a}_{\beta}\epsilon_{\alpha,\beta}=0$ $\displaystyle
f^{a\dagger}_{\alpha}f^{a\dagger}_{\beta}\epsilon_{\alpha,\beta}=0$

where $a=1\cdots N$ labels the flavor index and $\alpha=\uparrow,\downarrow$
labels the spin index. Repeated indices are summed over in Eq. (65) following
the Einstein convention. $\epsilon$ is antisymmetric tensor.

The coupling between the impurity and bulk $Z_{2}$ spin liquid is

$H_{imp}=\sum_{\langle
j,imp\rangle}\frac{J}{N}\mathbf{S}_{j}^{ab}\cdot\mathbf{S}_{imp}^{ba}+\frac{J^{\prime}}{N^{3}}(\mathbf{S}_{j}^{ab}\cdot\mathbf{S}_{imp}^{ba})^{2}$
(66)

where $sp(2N)$ spin $\mathbf{S}$ ion the bulk is defined asRan and Wen (2006)

$\displaystyle
S^{ab+}=\frac{1}{2}(c^{a\dagger}_{\uparrow}c_{\downarrow}^{b}+c^{b\dagger}_{\uparrow}c_{\downarrow}^{a}),~{}~{}~{}S^{ab,z}=\frac{1}{2}(c^{a\dagger}_{\uparrow}c_{\uparrow}^{b}-c^{b\dagger}_{\downarrow}c_{\downarrow}^{a})$
(67)

and the impurity spin ${\bf S}_{imp}$ is defined as

$\displaystyle
S_{imp}^{ab+}=\frac{1}{2}(f^{a\dagger}_{\uparrow}f_{\downarrow}^{b}+f^{b\dagger}_{\uparrow}f_{\downarrow}^{a}),~{}~{}~{}S_{imp}^{ab,z}=\frac{1}{2}(f^{a\dagger}_{\uparrow}f_{\uparrow}^{b}-f^{b\dagger}_{\downarrow}f_{\downarrow}^{a})$
(68)

$\langle j,imp\rangle$ denotes the 4 nearest neighbors $j$ to the impurity
site, i.e. impurity is on the center of plaquette with those 4 spins on the
corner, as shown in Fig. 8.

Here, we introduce the $J^{\prime}$ term to address a technical issue. With
only the $J$ term, the mean field saddle points will be highly degenerate,
which means infinitely many different mean field ansatz share the same ground
state energy, an artifact of the $J$ only mean-field theory. The introduction
of $J^{\prime}$ term reduces this infinite degeneracy to a finite fold. This
technique is introduced in Ref.Marston and Affleck (1989). Following the
standard procedure, we first perform a Hubbard-Stratonovich transformation on
the $J^{\prime}$ term and write Hamiltonian in the basis of fermions asRan and
Wen (2006)

$\displaystyle H_{imp}=-\sum_{\langle
j,imp\rangle}\frac{J(1+2\Phi_{j}J^{\prime}/J)}{4N}\left(\hat{\eta}_{j,imp}^{aa\dagger}\hat{\eta}_{j,imp}^{bb}+\hat{\chi}_{j,imp}^{aa\dagger}\hat{\chi}_{j,imp}^{bb}\right)-J^{\prime}\Phi_{j}^{2}$
(69)

Where,

$\displaystyle\hat{\chi}_{j,imp}^{aa}=c_{j\uparrow}^{a\dagger}f_{\uparrow}^{a}+c_{j\downarrow}^{a\dagger}f_{\downarrow}^{a}$
(70)
$\displaystyle\hat{\eta}_{j,imp}^{aa}=c_{j\uparrow}^{a}f_{\downarrow}^{a}-c_{j\downarrow}^{a}f_{\uparrow}^{a}$
(71)

Then, we do Hubbard-Stratonovich transformation to get mean field quadratic
Hamitonian and ignore fluctuation of mean field, the Hamiltonian can be
written as

$\displaystyle H_{imp}=$ $\displaystyle\sum_{\langle
j,imp\rangle}-J^{\prime}N\Phi_{j}^{2}+\frac{J}{4}(N|{\eta}_{j,imp}|^{2}+N|{\chi}_{j,imp}|^{2})$
(72)
$\displaystyle-\frac{J}{4}\sqrt{1+2\Phi_{j}J^{\prime}/J}(\eta_{j,imp}\hat{\eta}^{aa\dagger}_{j,imp}+\chi_{j,imp}\hat{\chi}^{aa\dagger}_{j,imp}+h.c.)$

Note that our convention is a bit different from Ref.Marston and Affleck
(1989). For simplicity, we define $\tilde{J}=J\sqrt{1+2\Phi_{j}J^{\prime}/J}$.
Because of the constraint (65) on the impurity site, we add Lagrangian
multipliers $a^{\pm}=a^{1}\pm ia^{2}$ and $a_{3}$. Since the fluctuation can
be ignored in the large $N$ limit, we have

$\displaystyle H_{imp}=$ $\displaystyle\sum_{\langle
j,imp\rangle}-J^{\prime}N\Phi_{j}^{2}+\frac{J}{4}(N|{\eta}_{j,imp}|^{2}+N|{\chi}_{j,imp}|^{2})$
(73)
$\displaystyle-\frac{\tilde{J}}{4}(\eta_{j,imp}(c_{j\uparrow}^{a\dagger}f_{\downarrow}^{a\dagger}-c_{j\downarrow}^{a\dagger}f_{\uparrow}^{a\dagger})+\chi_{j,imp}(c_{j\uparrow}^{a\dagger}f_{\uparrow}^{a}+c_{j\downarrow}^{a\dagger}f_{\downarrow}^{a})+h.c.)$
$\displaystyle+a^{+}f^{a}_{\uparrow}f^{a}_{\downarrow}+a^{-}f^{a\dagger_{\uparrow}}f^{a\dagger}_{\downarrow}+\frac{1}{2}a^{3}(f^{a\dagger}_{\uparrow}f^{a}_{\uparrow}+f^{a\dagger}_{\downarrow}f^{a}_{\downarrow})$
$\displaystyle=$ $\displaystyle\sum_{\langle
j,imp\rangle}-J^{\prime}N\Phi_{j}^{2}+\frac{J}{8}N\
\mathbf{tr}(u_{j,imp}^{\dagger}u_{j,imp})$
$\displaystyle+\frac{\tilde{J}}{4}(\psi_{j}^{\dagger}u_{j,imp}\psi_{imp}+h.c)+\psi_{imp}^{a\dagger}a_{i}\tau^{i}\psi_{imp}^{a}.$

Here, we use the notation in Ref.Wen (2002)

$u_{j,imp}=\left(\begin{array}[]{cc}\chi_{j,imp}^{\dagger}&\eta_{j,imp}\\\
\eta_{j,imp}^{\dagger}&-\chi_{j,imp}\end{array}\right)$ (74)

Combined with the bulk Hamiltonian (60), we have a quadratic Hamiltonian of
fermionic partons. We need to find the mean field saddle point solution for
$\chi,\eta$ and $\vec{a}$ that extremizes the free energy. In this case, a
saddle point is found where the energy is the minimum w.r.t. $\chi,\eta$ and
the maximum w.r.t. $\vec{a}$ and $\Phi$. We obtained the saddle point solution
by the gradient descent method.

###  C.3 Large-$N$ mean field solutions

Now we discuss the results of our large $N$ mean-field calculations. The
result is not sensitive to parameter $J^{\prime}$ and we fixed it to be
$J^{\prime}=1.2J$. The realistic situation of a spin-$1/2$ impurity added to a
square lattice of spin-$1/2$’s is just the $N=1$ case of the $Sp(2N)$ model
solved here. The main result of interest here, i.e. a symmetry protected non-
Kramers doublet localized at the impurity site, is qualitatively the same for
different $N\in 2\mathbb{Z}+1$. Therefore we expect the same observable
features of large $N$ results to also apply to the more realistic situation.

Figure 9: The phase diagram of isolated Kondo impurities in a $Z_{2}$ spin
liquid with a nontrivial PSG. There is a Kondo temperature $T_{K}(J)$, denoted
by the blue line, which separates the unscreened phase at high temperatures
and the Kondo screening phase at low temperatues.

As shown in Fig. 9, the $J-T$ phase diagram features two phases separated by a
Kondo temperature $T_{K}(J)$: (i) the unscreened phase where the magnetic
impurity behaves as a free magnetic moment, at a high temperature
$T>T_{K}(J)$, and (ii) the Kondo screening phase where the magnetic impurity
is screend by spinons, at a low temperature $T<T_{K}(J)$. Note that the Kondo
screening phase only happens when the Kondo coupling $J$ is larger than a
finite threshold $J_{c}\approx 5$, as expected for Kondo effects in a gapped
bulkWithoff and Fradkin (1990); Satori  _et al._ (1992); Saso (1992); Itoh
(1993); Takegahara  _et al._ (1993); Chen and Jayaprakash (1998). The criteria
for the Kondo screening is the magnitude of mean-field parameters $|\chi|$ and
$|\eta|$. When both parameters are zero, there is a free spin at the impurity
site, pointing to the unscreened phase. When the mean-field parameters become
nonzero, the impurity becomes a part of the bulk spin liquid, with a gapped
spectrum for the full Hamiltonian. This corresponds to the Kondo screening
phase.

In the Kondo screening phase, we obtain two saddle point solutions that
minimize the free energy:

$\displaystyle u_{(1,1),imp}=u_{1}i\tau^{0}$ (75) $\displaystyle
u_{(2,1),imp}=\frac{u_{1}}{\sqrt{2}}i\tau^{0}-\frac{u_{1}}{\sqrt{2}}\tau^{1}$
$\displaystyle
u_{(1,2),imp}=-\frac{u_{1}}{\sqrt{2}}i\tau^{0}-\frac{u_{1}}{\sqrt{2}}\tau^{1}$
$\displaystyle u_{(2,2),imp}=u_{1}\tau^{1}$ $\displaystyle a=a_{1}\tau^{1}$

and

$\displaystyle
u_{(1,1),imp}=-\frac{u_{1}}{\sqrt{2}}\tau^{2}+\frac{u_{1}}{\sqrt{2}}\tau^{3}$
(76) $\displaystyle u_{(2,1),imp}=u_{1}\tau^{2}$ $\displaystyle
u_{(1,2),imp}=u_{1}\tau^{3}$ $\displaystyle
u_{(2,2),imp}=\frac{u_{1}}{\sqrt{2}}\tau^{2}+\frac{u_{1}}{\sqrt{2}}\tau^{3}$
$\displaystyle a=-a_{1}\tau^{1}$

where $(x,y)$ with $x,y=1,2$ labels the four neighboring spins closest to the
impurity. $u_{1}$ and $a_{1}$ are two positive numbers that are functions of
$J$ and $T$. We label there two mean-field ansatz $A_{1}$ and $A_{2}$. They
are related by mirror symmetries $M_{x}$ and $M_{y}$, as illustrated in Fig.
10. The gauge transformation associated with $M_{x}$ and $M_{y}$ in the bulk
are shown in Eq. (62). On the impurity site, the gauge transformations
associated with the two mirrors are $G_{M_{x}}=i\tau^{2}$ and
$G_{M_{y}}=-i\tau^{3}$.

On a square lattice of $N^{2}$ sites ($N$= even) and 1 impurity spin, with an
odd number of fermions due to on-site constraint (65), due to the nontrivial
PSG in (52), $M_{x}$ and $M_{y}$ operators anticommute with each other in the
two mean-field ground states $|A_{1}\rangle$ and $|A_{2}\rangle$Lu (2018).
Since the two mean-field states are related by mirror symmetry $M_{x}$ (or
$M_{y}$), as shown in Fig. 10, the representation of the two mirror symmetries
in this two-dimensional space must have the following form in the basis of
$(|A_{1}\rangle,|A_{2}\rangle)^{T}$:

$\displaystyle R(M_{x})=\sigma_{x},~{}~{}~{}R(M_{y})=e^{i\phi}\sigma_{y}$ (77)

by a proper gauge choice, where $\sigma_{x,y,z}$ are the Pauli matrices. As a
result, any local perturbation $\hat{H}$ preserving $M_{x,y}$ symmetries must
have the following form:

$\displaystyle\langle A_{i}|\hat{H}|A_{j}\rangle=h\delta_{i,j}$ (78)

In other words, the two mean-field states form a projective representation of
the impurity site symmetry group $G_{s}=Z_{2}^{M_{x}}\times Z_{2}^{M_{y}}$,
analogous to the exactly solvable model discussed previously in section B.
This means the 2-fold degeneracy between saddle point solutions
$|A_{1}\rangle$ and $|A_{2}\rangle$ cannot be split by any symmetry-preserving
local perturbations. For example, we have explicitly checked that in our case

$\langle A_{1}|\sum_{\langle
j,imp\rangle}\vec{S}_{j}\cdot\vec{S}_{imp}|A_{2}\rangle\equiv 0$ (79)

As a result, for a $Z_{2}$ spin liquid with a nontrivial symmetry
fractionalization class, the Kondo screening phase is characterized by a
symmetry protected non-Kramers doublet at the impurity site.

Figure 10: Illustration of symmetry acting on ground states.
$|A_{1}\rangle,|A_{2}\rangle$ mean two states corresponding to two mean-field
ansatz.

To compare, we also consider the case of a bulk $Z_{2}$ spin liquid with a
trivial PSG, shown in Eq. (64), where $M_{x}$ and $M_{y}$ commute with each
other when acting on a single parton operator. The unscreened phase is still
featured by $\chi=\eta=0$, leaving a free spin at the impurity site.
Meanwhile, in the Kondo screening phase, there is a unique mean-field ansatz
as a gapped symmetric ground state. The associated mean-field ansatz is

$\displaystyle u_{(1,1),imp}=u_{1}i\tau^{0}$ (80) $\displaystyle
u_{(2,1),imp}=u_{1}i\tau^{0}$ $\displaystyle u_{(1,2),imp}=-u_{1}i\tau^{0}$
$\displaystyle u_{(2,2),imp}=-u_{1}i\tau^{0}$ $\displaystyle
a=a_{1}\tau^{1}+a_{2}\tau^{2}+a_{3}\tau^{3}$

where parameters $u_{1},a_{i},i=1\cdots 3$ are functions of $T,J$. This mean-
field ansatz is invariant under both mirror symmetries as we choose gauge
transformations $G_{M_{x}}=1,G_{M_{y}}=-1$ on the impurity site. This unique
paramagnetic ground state distinguishes itself from the 2-fold degeneracy in
the case of a nontrivial PSG, because here there is no local degeneracy
protected by mirror symmetries.

###  C.4 Thermodynamic quantities

We perform a self-consistent mean-field calculation in the large $N$ limit, on
a square lattice of $20\times 20$ sites. To avoid an extra spinon lurking in
the bulk, we consider two magnetic impurities separated by a largest distance
($d=10\sqrt{2}$) between each other.

In Fig. 11, we show the thermodynamic quantities from the mean-field
calculation to illustrate the experimental implications of the anomalous Kondo
effect discussed above.

The uniform magnetic susceptibility $\chi(T)$ behaves as $1/T$ for free spins,
and therefore $\chi\cdot T$ approaches a constant at low temperatures in the
unscreened phase. In the Kondo screening regime, however, $\chi(T)$ decays
exponentially as $\sim e^{-\Delta/k_{B}T}$ below the Kondo temperature
$T_{K}(J)$, where $\Delta$ is the spin gap of the system. Since the ground
states always feature a finite spin gap for both the case of a trivial and
that of a nontrivial PSG, the susceptibility can only differentiate unscreend
phase from the Kondo screening phase, but cannot differentiate a trivial PSG
from a nontrivial one.

On the other hand, the specific heat $C_{v}(T)$, or the entropy

$\displaystyle S(T)=\int_{0}^{T}\frac{C_{v}(t)}{t}\text{d}t$ (81)

at low temperatures can measure the density of states at a low energy, and can
be used to differentiate a trivial PSG from a nontrivial one. To be precise,
in the unscreened phase, each impurity site with a free moment contribute a
$k_{B}\ln 2$ at a low temperature. Meanwhile, in the Kondo screening regime,
for a trivial PSG, the system is gapped everywhere, leading to an
exponentially vanishing specific heat and entropy. This is in sharp contrast
to a nontrivial impurity, where each impurity site features a symmetry-
protected 2-fold non-Kramers doublet, and contributes a $k_{B}\ln 2$ entropy
at low temperatures.

Figure 11:  The temperature dependence of (a) uniform magnetic susceptibility
$\chi(T)$ and (b) entropy $S(T)$ contributions from Kondo impurities in
different regimes: the unscreened regime of free moments at the impurity sites
(green), the Kondo screening regime in $Z_{2}$ spin liquids with a trivial
(blue) vs. nontrivial (red) spinon fractionalization class. The calculations
are preformed on a $20\times 20$ square lattice with 2 spatially separated
magnetic impurities.

In reality, in a quasi-2d material with a dilute concentration of impurities
satisfying $\rho\xi^{2}\ll 1$, where $\rho$ is the impurity density and $\xi$
is the correlation length, we expect the observable signatures discussed above
to remain valid. In particular, there will be a low temperature peak in the
specific heat at $T\sim Je^{-C_{0}/\xi\sqrt{\rho}}$, where $C_{0}$ is a
constant of order one. The entropy plateau will persist between this low
temperature scale and the Kondo temperature $T_{K}(J)$.