$\lambda$ and $g$ interchange their roles. To derive it, we perform the
following unitary transformations; (1): $\pi$ rotation about the $\sigma_{y}$
axis, i.e., $\sigma^{x}_{ij}\rightarrow\sigma^{z}_{ij}$ ,
$\sigma^{z}_{ij}\rightarrow-\sigma^{x}_{ij}$; (2): $2\pi$ rotation about the
$\sigma_{z}$-axis, i.e., $\sigma^{x}_{ij}\rightarrow-\sigma^{x}_{ij}$,
$\sigma^{y}_{ij}\rightarrow-\sigma^{y}_{ij}$. In Eq.(82),
$A_{s}=\prod_{\langle ij\rangle\in+}\sigma^{x}_{ij}$ and $B_{p}=\prod_{\langle
ij\rangle\in\square}\sigma^{z}_{ij}$ are centered at the original ($r$) and
dual ($\tilde{r}=r+\frac{\hat{x}}{2}+\frac{\hat{y}}{2}$) lattice sites
respectively. Now, an elementary plaquette ($\square$) of a square lattice can
be thought of as a “star” ($+$)-like structure of the dual square lattice,
centered at the mid-point of the original lattice plaquette. Therefore,
$A_{s}(r)\rightarrow\tilde{B}_{p}(r)$ and
$B_{p}(\tilde{r})\rightarrow\tilde{A}_{s}(\tilde{r})$. Eq.(82) now becomes

$\displaystyle H^{dual}_{TC}$ $\displaystyle=-g\sum_{\langle
ij\rangle}\sigma^{z}_{ij}-\frac{1}{g}\sum_{+}\prod_{\langle
ij\rangle\in+}\sigma^{x}_{ij}$
$\displaystyle-\frac{1}{\lambda}\sum_{\square}\bigg{(}\prod_{\langle
ij\rangle\in+}\Theta_{ij}\bigg{)}\prod_{\langle
ij\rangle\in\square}\sigma^{z}_{ij}-\lambda\sum_{\langle
ij\rangle}\sigma^{x}_{ij}\Theta_{ij}$ (83)

Here, the $\sigma^{\alpha}_{ij}$ ($\Theta_{ij}$) are defined on the dual
(original) lattice links. We now rescale $\sigma^{\alpha}_{ij}$ by
$\Theta_{ij}$, by
$\sigma^{\alpha}_{ij}\rightarrow\Theta_{ij}\theta^{\alpha}_{ij}$ (where
$\theta^{\alpha}_{ij}$ are another set of Pauli spins). This scaling
transformation, keeps the spin commutation relations intact. The Hamiltonian
transforms accordingly,

$\displaystyle H^{dual}_{TC}=-\frac{1}{g}\sum_{+}\bigg{(}\prod_{\langle
ij\rangle\in+}$ $\displaystyle\Theta_{ij}\bigg{)}\prod_{\langle
ij\rangle\in+}\theta^{x}_{ij}-\frac{1}{\lambda}\sum_{\square}\prod_{\langle
ij\rangle\in\square}\theta^{z}_{ij}$ $\displaystyle-\lambda\sum_{\langle
ij\rangle}\theta^{x}_{ij}-g\sum_{\langle ij\rangle}\Theta_{ij}\theta^{z}_{ij}$
(84)

Notice that the $\Theta_{ij}$ fields were defined on the links (or the center
of the links) of original lattice. We can also think of these as sitting at
the center of the dual lattice links, since both are at the same location.
Hence, for the dual Toric code, $\Theta$ fields are now redefined on the dual
lattice links. Therefore, $H^{dual}_{TC}(g,\lambda)=H_{TC}(\lambda,g)$.

$\ast$Transforming the toric code model with external Zeeman fields (Eq.(82))
to a $Z_{2}$ gauge Hamiltonian coupled with $Z_{2}$ matter fields (Eq.(43)):

Although the local $Z_{2}$ symmetry ($[A_{s},H_{TC}]=0$) of the pure toric
code is explicitly broken in the presence of Zeeman fields, a modified local
invariance now appears. Notice that $\sigma^{z}_{ij}$ creates a $Z_{2}$
charge-pair ($A_{s}=-1$) at the lattice sites $i$ and $j$, beside flipping
$\sigma^{x}_{ij}$ (i.e., creating an electric field line connecting the two
$Z_{2}$ charges). We define the charge (number) operator $s^{x}_{i}$ which has
the eigenvalue $-1$, when there is a $Z_{2}$ charge sitting at site $i$. So,
in the presence of external (longitudinal) magnetic fields, we have the
following local $Z_{2}$ constraints on the physical Hilbert space [40],

$\displaystyle s^{x}_{i}A_{i}\ket{\text{Phys}}=\bigg{(}\prod_{\langle
ij\rangle\in+}\Theta_{ij}\bigg{)}\ket{\text{Phys}}$ (85)

The role of $\Theta_{ij}$ fields are now clear; it measures the static
background charge configuration of the model. To rewrite Eq.(82) in terms of
$s_{i}^{\alpha}$ and $\sigma^{\alpha}_{ij}$, we consider the following
transformations (under the constraint Eq.(85))[40]

$\displaystyle A_{i}\rightarrow s^{x}_{i}\bigg{(}\prod_{\langle
ij\rangle\in+}\Theta_{ij}\bigg{)}\ \ \ \ \sigma^{z}_{ij}\rightarrow
s^{z}_{i}\sigma^{z}_{ij}s^{z}_{j}$ (86)

The last transformation reflects the fact that $\sigma^{z}_{ij}$ creates a
charge-pair at lattice sites $i$ and $j$ ($s^{x}_{i,j}=-1$), and it keeps
$B_{p}$ unchanged. So, in the enlarged Hilbert space of gauge fields
($\sigma^{\alpha}_{ij}$) and matter ($s^{\alpha}_{i}$), the Hamiltonian
(Eq.(82)) becomes

$\displaystyle H_{g+m}$
$\displaystyle=-\frac{1}{\lambda}\sum_{i}s^{x}_{i}-\lambda\sum_{\langle
ij\rangle}s^{z}_{i}\Theta_{ij}\sigma^{z}_{ij}s^{z}_{j}$
$\displaystyle-g\sum_{\langle
ij\rangle}\sigma^{x}_{ij}-\frac{1}{g}\sum_{\square}\prod_{\langle
ij\rangle\in\square}\sigma^{z}_{ij}$ (87)

with the local constraint Eq.(85). Following a similar set of steps, we can
convert the dual toric code, Eq.(84) to the following Hamiltonian,

$\displaystyle H^{dual}_{g+m}=-\frac{1}{g}\sum_{i}v^{x}_{i}$
$\displaystyle-g\sum_{\langle
ij\rangle}v^{z}_{i}\Theta_{ij}\theta^{z}_{ij}v^{z}_{j}$
$\displaystyle-\lambda\sum_{\langle
ij\rangle}\theta^{x}_{ij}-\frac{1}{\lambda}\sum_{\square}\prod_{\langle
ij\rangle\in\square}\theta^{z}_{ij}$ (88)

now with the following local gauge constraint,

$\displaystyle
v^{x}_{i}\prod_{(ij)\in+}\theta^{x}_{ij}=\prod_{(ij)\in+}\Theta_{ij}\equiv\prod^{\prime}_{(ij)\in\square}\Theta_{ij}$
(89)

The prime on the product denotes that this product is performed on the dual
lattice links. In Eq.(88), the $v^{\alpha}_{i}$ denotes vison (or $Z_{2}$
flux) degrees of freedom of the original gauge theory Eq.(87).
Now, consider the limiting case, $g\rightarrow 0$ in Eq.(87) (pure Ising model
limit), The gauge fields are completely frozen, i.e. $\prod_{\langle
ij\rangle\in\square}\sigma^{z}_{ij}=+1$ for all plaquettes. One can choose the
gauge, $\sigma^{z}_{ij}=1$ for all the links.

$\displaystyle H_{g+m}(g\rightarrow
0)=H_{I}=-\frac{1}{\lambda}\sum_{i}s^{x}_{i}-\lambda\sum_{\langle
ij\rangle}s^{z}_{i}\Theta_{ij}s^{z}_{j}$ (90)

On the other hand,

$\displaystyle H^{dual}_{g+m}(g\rightarrow 0)=H_{gauge}=-\lambda\sum_{\langle
ij\rangle}\theta^{x}_{ij}-\frac{1}{\lambda}\sum_{\square}\prod_{\langle
ij\rangle\in\square}\theta^{z}_{ij}$ (91)

with the constraint, $\prod_{\langle
ij\rangle\in+}\theta^{x}_{ij}=\prod^{\prime}_{\langle
ij\rangle\in\square}\Theta_{ij}$. This is exactly the duality between the TFIM
with frustrated plaquettes and the Ising gauge theory with background charges.
Frustrated plaquettes of the Ising model map to the (static) background
charges of the $Z_{2}$ gauge theory.

##  Appendix C Brief description of TPSC equations

The equations for charge, spin susceptibilities, and single-particle self-
energy used in the main text are derived in this section. With the exception
of the fact that in our model pseudo-spin rotational symmetry is not preserved
(because of the directional hopping of orbitals), the derivation is almost
identical to that presented in Ref.[30].

Expressions for the charge and longitudinal spin susceptibilities:

For three-point susceptibility in the particle-hole channel, the Bethe-
Salpeter equation is given by [30],

$\displaystyle\chi_{\sigma\sigma^{\prime}}(1,3;2)=$ $\displaystyle-
G_{\sigma}(1,2)G_{\sigma}(2,3)\delta_{\sigma\sigma^{\prime}}$
$\displaystyle+G_{\sigma}(1,\bar{2})\Gamma_{\sigma\bar{\sigma}}^{ir}(\bar{2},\bar{3};\bar{4},\bar{5})\chi_{\bar{\sigma}\sigma^{\prime}}(\bar{4},\bar{5};2)G_{\sigma}(\bar{3},3)$
(92)

One of the main assumptions of TPSC theory is the replacement of the momentum
and frequency-dependent vertex functions by constant parameters, that are
determined self-consistently from the sum-rules.

$\displaystyle\Gamma_{\sigma\bar{\sigma}}^{ir}(\bar{2},\bar{3};\bar{4},\bar{5})=\Gamma_{\sigma\bar{\sigma}}\delta(\bar{2}-\bar{5})\delta(\bar{3}-\bar{4})\delta(\bar{4}^{+}-\bar{5})$
(93)

The two-point susceptibility is defined as
$\chi_{\sigma\sigma^{\prime}}(1,3\rightarrow
1^{+};2)\equiv\chi_{\sigma\sigma^{\prime}}(1,2)$, that satisfies the following
equation,

$\displaystyle\chi_{\sigma\sigma^{\prime}}(1,2)$
$\displaystyle=-G_{\sigma}(1,2)G_{\sigma}(2,1^{+})\delta_{\sigma\sigma^{\prime}}$
$\displaystyle+G_{\sigma}(1,\bar{2})\Gamma_{\sigma\bar{\sigma}}\chi_{\bar{\sigma}\sigma^{\prime}}(\bar{2},2)G_{\sigma}(\bar{2},1^{+})$
(94)

Hence,

$\displaystyle\chi_{dd}(1,2)=-G_{d}(1,2)G_{d}(2,1^{+})+G_{d}(1,\bar{2})\Gamma_{dd}\chi_{dd}(\bar{2},2)$
$\displaystyle\times
G_{d}(\bar{2},1^{+})+G_{d}(1,\bar{2})\Gamma_{df}\chi_{fd}(\bar{2},2)G_{d}(\bar{2},1^{+})$
(95)

and

$\displaystyle\chi_{fd}(1,2)=G_{f}(1,\bar{2})\Gamma_{ff}\chi_{fd}(\bar{2},2)G_{f}(\bar{2},1^{+})$
$\displaystyle+G_{f}(1,\bar{2})\Gamma_{fd}\chi_{dd}(\bar{2},2)G_{f}(\bar{2},1^{+})$
(96)

In the Hamiltonian Eq.(57), the anisotropy arises due to the different spatial
dependence of single particle Green’s functions. If we consider any local
quantity (like $\Gamma_{\sigma\bar{\sigma}}$ in TPSC), that should have the
following symmetry: $\Gamma_{dd}=\Gamma_{ff}$, $\Gamma_{df}=\Gamma_{fd}$. This
happens as the local Green’s functions are identical
$G^{loc}_{d}(i\omega)=G^{loc}_{f}(i\omega)$, although $k$-dependent Green’s
functions are (numerically) different. The pseudo-spin rotational symmetry is
preserved at the local level as there is no magnetic field term in the
Hamiltonian (or any kind of spontaneously broken symmetry).
From Eqs. (95), (96), we arrive (in the Fourier space) at the following
expression

$\displaystyle\chi_{\sigma\sigma}(q)=$ $\displaystyle\
\frac{[1+\Gamma_{\sigma\sigma}\chi^{0}_{\bar{\sigma}}(q)]\chi^{0}_{\sigma}(q)}{1+\Gamma_{\sigma\sigma}(\chi^{0}_{\sigma}(q)+\chi^{0}_{\bar{\sigma}}(q))+(\Gamma^{2}_{\sigma\sigma}-\Gamma^{2}_{\sigma\bar{\sigma}})\chi^{0}_{\sigma}(q)\chi^{0}_{\bar{\sigma}}(q)}$
(97)

Here, we have defined
$\chi^{0}_{\sigma}(1,2)=-G_{\sigma}(1,2)G_{\sigma}(2,1^{+})$, the so-called
“particle-hole bubble” or Lindhard function. Similarly,

$\displaystyle\chi_{\sigma\bar{\sigma}}(q)=$
$\displaystyle\frac{-\Gamma_{\sigma\bar{\sigma}}\chi^{0}_{\sigma}(q)\chi^{0}_{\bar{\sigma}}(q)}{1+\Gamma_{\sigma\sigma}(\chi^{0}_{\sigma}(q)+\chi^{0}_{\bar{\sigma}}(q))+(\Gamma^{2}_{\sigma\sigma}-\Gamma^{2}_{\sigma\bar{\sigma}})\chi^{0}_{\sigma}(q)\chi^{0}_{\bar{\sigma}}(q)}$
(98)

The charge and (longitudinal) spin susceptibilities are defined as
$\chi_{\alpha}(q)=\chi_{dd}(q)+\chi_{ff}(q)\pm 2\chi_{df}(q)$, where
$\alpha=[\text{ch},\ \text{sp}]$,

$\displaystyle\chi_{\alpha}(q)=$
$\displaystyle\frac{\big{[}\chi^{0}_{d}(q)+\chi^{0}_{f}(q)\mp
2\Gamma_{df}\chi^{0}_{d}(q)\chi^{0}_{f}(q)\big{]}+2\Gamma_{dd}\chi^{0}_{f}(q)\chi^{0}_{d}(q)}{\big{[}1+\Gamma_{dd}(\chi^{0}_{d}(q)+\chi^{0}_{f}(q))+(\Gamma^{2}_{dd}-\Gamma^{2}_{df})\chi^{0}_{d}(q)\chi^{0}_{f}(q)\big{]}}$
(99)

Now, using the canonical definitions of the charge and spin vertices,
$\Gamma_{ch}=\Gamma_{df}+\Gamma_{dd}$, $\Gamma_{sp}=\Gamma_{df}-\Gamma_{dd}$,
we get the following expressions for charge and (longitudinal) spin
susceptibilities

$\displaystyle\chi_{\alpha}(q)=\frac{\big{[}\chi^{0}_{d}(q)+\chi^{0}_{f}(q)\mp(\Gamma_{ch}+\Gamma_{sp})\chi^{0}_{d}(q)\chi^{0}_{f}(q)\big{]}+(\Gamma_{ch}-\Gamma_{sp})\chi^{0}_{f}(q)\chi^{0}_{d}(q)}{\bigg{[}1+\frac{(\Gamma_{ch}-\Gamma_{sp})}{2}[\chi^{0}_{d}(q)+\chi^{0}_{f}(q)]-\Gamma_{ch}\Gamma_{sp}\chi^{0}_{d}(q)\chi^{0}_{f}(q)\bigg{]}}$
(100)

($-$ve sign for charge and $+$ve sign for spin ($z$ component)
susceptibility). In the absence of pseudospin rotational symmetry,
$\Gamma_{ch}$ and $\Gamma_{sp}$ get mixed in $\chi_{ch}$, $\chi_{sp}$. We can
check the correctness of the above expression in the RPA limit, where
$\Gamma_{ch}=\Gamma_{sp}=U$, the previous equation reduces to,

$\displaystyle\chi^{\text{RPA}}_{\alpha}(q)=\frac{\chi_{d}^{0}(q)+\chi_{f}^{0}(q)\mp
2U\chi_{d}^{0}(q)\chi_{f}^{0}(q)}{1-U^{2}\chi_{d}^{0}(q)\chi_{f}^{0}(q)}$

which matches exactly with the RPA formula, which could be easily derived from
the Feynman diagrammatic calculations.
In TPSC theory, the unknown charge and spin vertex $\Gamma_{ch}\equiv U_{ch}$
and $\Gamma_{sp}\equiv U_{sp}$ are determined from the following sum rules,

$\displaystyle\frac{T}{L^{2}}\sum_{\mathbf{q}}\sum_{iq_{n}}[\chi_{ch}(\mathbf{q},iq_{n})+\chi_{sp}(\mathbf{q},iq_{n})]=2n-n^{2}$
(101)

here $n=\langle n_{d}\rangle+\langle n_{f}\rangle$. To determine
$U_{sp},U_{ch}$, we solve the following two integral equations,

$\displaystyle\frac{T}{L^{2}}\sum_{\mathbf{q}}\sum_{iq_{n}}\chi_{sp}(\mathbf{q},iq_{n},U_{ch},U_{sp})=n-2\langle
n_{d}n_{f}\rangle$ (102)
$\displaystyle\frac{T}{L^{2}}\sum_{\mathbf{q}}\sum_{iq_{n}}\chi_{ch}(\mathbf{q},iq_{n},U_{ch},U_{sp})=n+2\langle
n_{d}n_{f}\rangle-n^{2}$ (103)

for a given value of the double occupancy $\langle n_{d}n_{f}\rangle$. We need
one more equation for $\langle n_{d}n_{f}\rangle$ to complete the above set.
In TPSC, that equation turns out to be a local equation which does not change
in the presence of pseudo-spin dependent hopping,

$\displaystyle U_{sp}=U\frac{\langle n_{d}n_{f}\rangle}{\langle
n_{d}\rangle\langle n_{f}\rangle}$ (104)

Solving (102), (103), and (104), we determine the self-consistent solutions
for $U_{sp},U_{ch}$, and $\langle n_{d}n_{f}\rangle$.

Expression for the transverse spin susceptibility:

Next, we derive the expression for transverse spin susceptibility. To derive
it, we at first consider it’s RPA-level expression. The RPA longitudinal
susceptibilities are found by summing the infinite set of “bubble” diagrams
and for the transverse ones ($\chi_{\pm\mp}\sim\langle
S^{\pm}_{i}(\tau)S^{\mp}_{j}(0)\rangle$), one has to sum a specific set of
“ladder”-like diagrams in the particle-hole channel upto infinite order.

$\displaystyle\chi^{+-}(q)=$
$\displaystyle\frac{1}{L^{2}}\sum_{\mathbf{k},\mathbf{k}^{\prime}}\int_{0}^{\beta}d\tau
e^{iq_{n}\tau}\langle
T_{\tau}d^{{\dagger}}_{\mathbf{k}}(\tau)f_{\mathbf{k}+\mathbf{q}}(\tau)f^{{\dagger}}_{\mathbf{k}^{\prime}}(0)d_{\mathbf{k}^{\prime}-\mathbf{q}}(0)\rangle$
(105)

After summing the infinite set of particle-hole ladder diagrams, one finds

$\displaystyle\chi^{+-}_{RPA}(q)=\frac{\Phi^{0}_{df}(q)}{1-U\Phi^{0}_{df}(q)}$
(106)

Here

$\displaystyle\Phi^{0}_{df}(\mathbf{q},iq_{n})$ $\displaystyle=-\frac{1}{\beta
L^{2}}\sum_{\mathbf{k},ik_{n}}G^{0}_{f}(\mathbf{k},ik_{n})G^{0}_{d}(\mathbf{k}+\mathbf{q},ik_{n}+iq_{n})$
$\displaystyle=-\frac{1}{\beta}\sum_{ik_{n}}G^{0,\text{loc}}_{f}(ik_{n})G^{0,\text{loc}}_{d}(ik_{n}+iq_{n})$
(107)

The non-interacting transverse spin susceptibility turns out to be a local
quantity. The reason behind this locality lies in the phase space decoupling
of the momentum summation, which arises from the directional nature of the
fermionic band structure. From Eq.(105), it becomes clear that a particle-hole
pair of opposite pseudo-spins created at site $i$ and $\tau=0$ will never meet
again at any later time $\tau>0$ and for any $j\neq i$ due to orthogonal band
motion of the fermions, when they are non-interacting.
In TPSC, the bare $U$ in (106) is replaced by a renormalized $U_{t}$. The
transverse susceptibility is expected to satisfy the following sum rule,

$\displaystyle\frac{1}{\beta}\sum_{iq_{n}}\frac{\Phi^{0}_{df}(iq_{n})}{1-U_{t}\Phi^{0}_{df}(iq_{n})}=\langle
d^{{\dagger}}_{i}f_{i}f^{{\dagger}}_{i}d_{i}\rangle=\langle
n_{d}\rangle-\langle n_{d}n_{f}\rangle$ (108)

We find that $U_{t}$ is negligible compared to $U_{ch}$ and $U_{sp}$ (in the
fermionic QCM). Due to momentum independence of $\chi^{x(y)}(iq_{n})$, we
don’t observe any singularity in this particular particle-hole channel. The
reason lies in the absence of pseudo-spin rotational symmetry. The
$\chi^{x(y)}(iq_{n})$ has some quantitative effects on the single-particle
self-energy (will give a non-singular, frequency dependent contribution), but
not in the description of ordering instabilities.

Figure 9: Transverse spin susceptibility, $\chi_{trans}(i\omega_{n})$ vs.
$\omega_{n}=2\pi nT$ for $T=0.01$. Compared to $\chi_{sp}$ (see Fig.10),
$\chi_{trans}$ does not show any singularity even at lower temperatures.

Expression for the self-energy:

The longitudinal part of the single-particle self-energy could be formally
expressed by the following relation,

$\displaystyle\Sigma^{l}_{\sigma}(1,2)=$ $\displaystyle
Un_{-\sigma}\delta(1-2)$
$\displaystyle+UG_{\sigma}(1-\bar{2})\Gamma^{ir}_{\sigma\bar{\sigma}}(\bar{2},2;\bar{4},\bar{5})\chi_{\bar{\sigma},-\sigma}(\bar{4},\bar{5},1)$
(109)

Using the local (real space) approximation ansatz for the vertex functions
(93), we find after some straightforward algebra the following,

$\displaystyle\Sigma^{l}_{\sigma}(k)=$ $\displaystyle
Un_{\bar{\sigma}}+\frac{U}{4\beta
L^{2}}\sum_{q}\big{[}U_{ch}\chi_{ch}(q)+U_{sp}\chi_{sp}(q)\big{]}G_{\sigma}(k+q)$
$\displaystyle-\frac{U}{4\beta
L^{2}}(U_{ch}+U_{sp})\sum_{q}[\chi_{\sigma\sigma}(q)-\chi_{\bar{\sigma}\bar{\sigma}}(q)]G_{\sigma}(k+q)$
(110)

In the above, the last term in $\Sigma^{l}$ arises due to absence of spin
rotational invariance. We find that this term has no significant effect in the
renormalized-classical regime where only $\chi_{sp}$ controls the leading
(singular) behavior.
The transverse component of the self-energy is local due to the momentum
independence of transverse susceptibility. Its expression is given as follows

$\displaystyle\Sigma^{(t)}_{\sigma}(ik_{n})=\frac{U}{\beta}\sum_{iq_{n}}\frac{U_{t}\Phi^{0}_{df}(iq_{n})}{1-U_{t}\Phi^{0}_{df}(iq_{n})}G^{0,loc}_{-\sigma}(ik_{n}+iq_{n})$
(111)

##  Appendix D Supplementary results from TPSC calculations

The non-interacting susceptibility of the individual orbitals $(d,f)$ show
logarithmic $T$-divergence due to the completely nested one dimensional Fermi
surface, i.e., $\chi^{\sigma}_{0}(Q)\sim N(0)\ln{(E_{0}/T)}$; here,
$N(0)=1/2\pi t$ is equal to the density of states at the Fermi surface (here
Fermi points) and $E_{0}\sim 2t$, the half-bandwidth. As a result, one expects
instabilities to occur at the corresponding wave-vectors. Here, we show the
interacting longitudinal spin susceptibility, $\chi_{sp}(\mathbf{q},i\nu=0)$
(see Eq.(100)) for certain $(T,U)$ values (one slightly inside the
$\xi>>\xi_{th}$ regime and the other one for the opposite regime).

Figure 10: Interacting $\chi_{sp}(q_{x},q_{y},i\nu=0)$ calculated by solving
TPSC self-consistency equations (Eqs. (102)-(104)). (a) $T=0.2$ and $U=1.0$,
in the $\xi<\xi_{th}$ regime (upper figure), (b) $T=0.2$ and $U=4.0$, in the
$\xi>\xi_{th}$ regime (lower figure). In the high-$T$ regime, a broad diffused
maximum in the scattering structure factor indicates a disordered state with
short correlation length (upper). At low-$T$, this maximum sharpens up at
$\mathbf{Q}=(\pi,\pi)$, testifying to onset of orbital order (lower).

.

As one can see from Fig.10 (lower one) that $\chi_{sp}$ obtain very large
numerical values near $(\pi,\pi)$ inside the renormalized classical regime
($\xi>>\xi_{th}$). This means the system develops a strong nearest-neighbour
$2d$ anti-ferro orbital correlations in this regime.
Next, We observe that $U_{ch}$ and $U_{sp}$ are almost same for very small
$U$. Then, as $U$ increases, the charge vertex rises rapidly, whereas the spin
vertex saturates to a finite value (see Fig.11). The rise of $U_{ch}$
corresponds suppression of charge collective modes with the increase in $U$.
The saturation of $U_{sp}$ is in accord with Kanamori-Bruckner renormalization
phenomenon.

Figure 11: $U_{ch}$, $U_{sp}$ vs. $U$ calculated from self-consistency
equations (Eqs.(101)-(103)) for different $T$ values.

We calculate the charge and spin correlation lengths, defined by
$\xi_{ch}\sim\sqrt{\chi_{ch}(\mathbf{Q},i\nu=0)/\chi^{0}_{ch}(\mathbf{Q},i\nu=0)}$
and
$\xi_{sp}\sim\sqrt{\chi_{sp}(\mathbf{Q},i\nu=0)/\chi^{0}_{sp}(\mathbf{Q},i\nu=0)}$,
where $\mathbf{Q}=(\pi,\pi)$ is the vector where $\chi(q,\nu=0)$ is maximum,
at half-filling. We see the suppression of charge fluctuations as $U$ is
increased, and also as temperature is lowered (see Fig.12). Although, TPSC
cannot capture metal-insulator phase transition (which involves frequency
dependent vertex corrections), it correctly produces the low-$T$ collective
behaviors.

Figure 12: charge correlation length, $\xi_{ch}$ vs. $U$ at half-filling for
different $T$s, showing almost no sensitivity to $T$.

On the other side, $\chi_{sp}$ shows opposite behavior (see Fig.13), it grows
exponentially as $T$ is lowered. The system enters into the so-called
renormalized classical regime where the wavelength of orbital pseudo-spin
fluctuations becomes larger than thermal de-Broglie wavelength of the
excitations. One expects around this temperature, at half-filling, some
crossover happens between a Fermi liquid at high-$T$ to a low-$T$ pseudo-gap
phase driven by anti-ferro “orbital” correlations.

Figure 13: longitudinal spin correlation length, $\xi_{sp}$ vs. $U$ at half-
filling, for different $T$s. Onset of orbital order within TPSC is signalled
by a finite $T$ divergence of $\xi_{sp}$. The $T-U$ phase diagram in the main
text is constructed from this observation.

We have used $128\times 128$ spatial lattice grid and $N_{\tau}=8192$ time
grid to perform the above numerical calculations. As $\xi_{sp}$ diverges
exponentially at low-$T$, it becomes increasingly difficult to correctly
capture some low-$T$ data points with the above lattice size. The quantitative
accuracy could be increased by choosing much larger lattice grid without much
change in the qualitative behavior.
Finally, we show the double occupancy, $\langle n_{d}n_{f}\rangle$ vs. $U$,
for different $T$ (see Fig.14). It decreases continuously from $n^{2}/4$ as
$U$ is raised. Suppression of $\langle n_{d}n_{f}\rangle$ at higher-$T$ values
($\sim T=0.4$) compared to low-$T$ is consistent with the self-energy results
(shown in the main text) which finds incoherent, insulator-like features
(suppression of charge fluctuations) at higher-$T$ values. Here, the ($T,U$)
points in Fig.14 lie on the $\xi<\xi_{th}$ part of the phase diagram.

Figure 14: Double occupancy $\langle n_{\uparrow}n_{\downarrow}\rangle$ vs.
$U$ for different $T=0.1,0.2,0.4$, at half-filling (here $\uparrow\equiv d$,
and $\downarrow\equiv f$)

##  Appendix E Self-energy in the $\xi>\xi_{th}$ regime

In this section, we will derive the analytical expression of the single-
particle self-energy in the renormalized-classical (RC) regime, i.e.
$\xi>>\xi_{th}$ following Ref. [30]. In the present context, the $G^{(0)}$ has
one dimensional $k$-dependence, although the susceptibilities are two
dimensional functions. Close to this RC regime, spin correlations dominate,
and the long wavelength and low-energy longitudinal spin susceptibility
obtains the following form [85],

$\displaystyle\chi_{sp}(\mathbf{q},i\omega)\approx\frac{2\xi^{2}}{U_{sp}\xi_{0}^{2}}\frac{1}{1+\mathbf{(q-\mathbf{Q})}^{2}\xi^{2}-i\omega/\omega_{sp}}$
(112)

Here, $\omega_{sp}=D/\xi^{2}$, the characteristic spin relaxation frequency
and $D=(\tau_{0}/\xi^{2}_{0})^{-1}$, the microscopic diffusion co-efficient.
The microscopic correlation length is
$\xi_{0}=\frac{-1}{2\chi_{0}(\mathbf{Q})}\frac{\partial^{2}\chi_{0}(\mathbf{q})}{\partial
q^{2}}\big{|}_{q=\mathbf{Q}}$ and microscopic correlation time is
$\tau_{0}=\frac{1}{\chi_{0}(\mathbf{Q})}\frac{\partial\chi_{0}(\mathbf{Q,\omega})}{\partial\omega}\big{|}_{\omega=0}$

Close to the RC regime,
$\xi\sim\exp\big{[}\pi\xi^{2}_{0}\sigma^{2}(T)U_{sp}/T\big{]}$, where
$\sigma(T)$ is the (weakly) temperature dependent spin-stiffness. This implies
$\omega_{sp}<<T$ and the susceptibility is dominated by the classical
($i\omega_{n}=0$) part, $\chi_{sp}(\mathbf{q},0)$. The self-energy is mostly
affected by $\chi_{sp}$ compared to $\chi_{ch}$, which could be expressed as
following,

$\displaystyle\Sigma_{\sigma}(\mathbf{k},ik_{n})\approx$ $\displaystyle
Un_{-\sigma}+\frac{U}{4}\frac{1}{\beta
N}\sum_{\mathbf{q}}U_{sp}\chi_{sp}(\mathbf{q},0)G^{(0)}_{\sigma}(k_{\sigma}+q_{\sigma},ik_{n})$
$\displaystyle+\frac{U}{4}\frac{1}{\beta N}\sum_{\mathbf{q}}\sum_{q_{n}\neq
0}U_{sp}\chi_{sp}(\mathbf{q},iq_{n})G^{(0)}_{\sigma}(k_{\sigma}+q_{\sigma},ik_{n}+iq_{n})$
(113)

As argued above, we neglect the $iq_{n}\neq 0$ terms, as $\omega_{sp}<<T$ and
the classical (thermal) spin correlations dominate.

$\displaystyle\Sigma_{\sigma}(\mathbf{k},ik_{n})\approx$ $\displaystyle\ \ \
\frac{UT}{2\xi^{2}_{0}}\int\frac{d^{2}q}{(2\pi)^{2}}\frac{1}{q^{2}+\xi^{-2}}\frac{1}{ik_{n}-\xi_{\mathbf{k}+\mathbf{Q}}-q_{\sigma}v^{\sigma}_{\mathbf{k}+\mathbf{Q}}}$

Close to the Fermi surface, $k\approx k_{F}$, $\mathbf{Q}=(\pi,\pi)$, so
$\xi_{\mathbf{k}_{F}+\mathbf{Q}}=0$ and
$v^{\sigma}_{\mathbf{k}_{F}+\mathbf{Q}}=v_{F}$,

$\displaystyle\Sigma_{\sigma}(\mathbf{k},ik_{n})$
$\displaystyle=\frac{UT}{2\xi^{2}_{0}v_{F}}\int_{-\infty}^{\infty}\frac{dq_{\sigma}}{2\pi}\frac{1}{i(k_{n}/v_{F})-q_{\sigma}}\int_{-\infty}^{\infty}\frac{dq_{\bar{\sigma}}}{2\pi}\frac{1}{q_{\bar{\sigma}}^{2}+q_{\sigma}^{2}+\xi^{-2}}$
$\displaystyle=\frac{UT}{2\xi^{2}_{0}v_{F}}\int_{-\infty}^{\infty}\frac{dq_{\sigma}}{2\pi}\frac{1}{\sqrt{q_{\sigma}^{2}+\xi^{-2}}[i(k_{n}/v_{F})-q_{\sigma}]}$

Now,
$\int_{-\infty}^{\infty}\frac{dx}{\sqrt{x^{2}+a^{2}}(ib-x)}=-\frac{i}{\sqrt{a^{2}-b^{2}}}\cos^{-1}(b/a)$,
$\cos^{-1}{x}=\tan^{-1}(\sqrt{1-x^{2}}/x)$, and
$\tan^{-1}(z)=-\frac{1}{2i}\ln{[(i+z)/(i-z)]}$, using these identities, we
obtain

$\displaystyle\Sigma_{\sigma}(\mathbf{k}_{F},ik_{n})=$ $\displaystyle\ \ \
-i\frac{UT}{8\pi\xi^{2}_{0}\sqrt{k_{n}^{2}-v^{2}_{F}\xi^{-2}}}\ln{\bigg{[}\frac{k_{n}+\sqrt{k_{n}^{2}-v^{2}_{F}\xi^{-2}}}{k_{n}-\sqrt{k_{n}^{2}-v^{2}_{F}\xi^{-2}}}\bigg{]}}$
(114)

Now, $k_{n}$ is linear in $T$, but $\xi^{-1}\sim e^{-C/T}$ ($C>0$), so
$k_{n}>>\xi^{-1}$ whenever $\xi>>\xi_{th}$. The self-energy in Matsubara
frequency can be also expressed as

$\displaystyle\Sigma_{\sigma}(\mathbf{k}_{F},ik_{n})=-i$
$\displaystyle\frac{U}{8(2n+1)\pi^{2}\xi^{2}_{0}}\bigg{[}1-\frac{\xi^{2}_{th}}{\xi^{2}}\frac{1}{(2n+1)^{2}}\bigg{]}^{-\frac{1}{2}}$
$\displaystyle\ \ \
\times\ln{\bigg{[}\frac{1+\sqrt{1-\frac{\xi^{2}_{th}}{\xi^{2}}\frac{1}{(2n+1)^{2}}}}{1-\sqrt{1-\frac{\xi^{2}_{th}}{\xi^{2}}\frac{1}{(2n+1)^{2}}}}\bigg{]}}$
(115)

so, when $\xi_{th}<<\xi$,
$\sqrt{1-\frac{\xi^{2}_{th}}{\xi^{2}}\frac{1}{(2n+1)^{2}}}\approx
1-\frac{\xi^{2}_{th}}{2\xi^{2}}\frac{1}{(2n+1)^{2}}$. Thus,

$\displaystyle\Sigma_{\sigma}(\mathbf{k}_{F},ik_{n})\approx-i\frac{U}{8(2n+1)\pi^{2}\xi^{2}_{0}}\ln{\big{[}4(2n+1)\xi/\xi_{th}\big{]}}$
(116)

Hence, for $k_{n}=\pi T$,
$\Sigma_{\sigma}(\mathbf{k}_{F},ik_{n=0})\sim-i\frac{U}{8\pi^{2}\xi^{2}_{0}}\ln{(4\xi/\xi_{th})}$
which diverges as $T\rightarrow 0$. When $\xi_{th}/\xi\sim O(1)$, the
matsubara self-energy at lowest Matsubara frequency $k_{n}=\pi T$, we may have
a different behavior. We want to see what happens if $k_{n=0}\rightarrow
v_{F}/\xi$ in (114) (which is same as $\xi\rightarrow\xi_{th}$), we find
$\Sigma_{\sigma}(k_{F},ik_{n=0})=-i\frac{UT}{2\pi\xi^{2}_{0}}+f(i\omega_{n})$,
where $f(\mathbf{k}_{F},i\omega_{n})$ is the regular part of the self-energy
which comes from the $iq_{n}\neq 0$ contributions of $\chi$, in Eq.(113). So,
expanding $f(k_{F},i\omega_{n})\approx-C$ near $\omega_{n}=0$,

$\displaystyle-\text{Im}\Sigma_{\sigma}(k_{F},i\omega_{n}\approx
0)\bigg{|}_{\xi/\xi_{th}\sim O(1)}\approx C+\frac{UT}{2\pi\xi^{2}_{0}}$ (117)

##  Appendix F Derivation of the strong-coupling effective spin Hamiltonian
from fermionic QCM

We consider the following generalized version of fermionic QCM which has the
$p$-wave Cooper pairing of strength $\lambda$ (here, $|\lambda|\leq 1$),

$\displaystyle H=$ $\displaystyle t\sum_{i}(\lambda
d_{i}d_{i+\hat{x}}-d^{{\dagger}}_{i}d_{i+\hat{x}}+h.c.)$
$\displaystyle+t\sum_{i}(\lambda
f_{i}f_{i+\hat{y}}-f^{{\dagger}}_{i}f_{i+\hat{y}}+h.c.)$
$\displaystyle+U\sum_{i}(2d^{{\dagger}}_{i}d_{i}-1)(2f^{{\dagger}}_{i}f_{i}-1)$
(118)

We write the complex fermions in the Majorana basis,
$f_{i}=(A^{b}_{i}+iA^{w}_{i})/2$, $d_{i}=(B^{w}_{i}+iB^{b}_{i})/2$. Here,
$(A^{b}_{i})^{2}=(A^{w}_{i})^{2}=(B^{b}_{i})^{2}=(B^{w}_{i})^{2}=1$. We get
the following Majorana Hamiltonian, $H=H_{0}+V_{x}+V_{y}$, where

$\displaystyle H_{0}=-UA^{w}_{i}A^{b}_{i}B^{b}_{i}B^{w}_{i}$ (119)
$\displaystyle
V_{x}=\frac{it}{2}(1+\lambda)(B^{b}_{i}B^{w}_{i+\hat{x}}+A^{w}_{i}A^{b}_{i+\hat{y}}),$
(120) $\displaystyle
V_{y}=-\frac{it}{2}(1-\lambda)(B^{w}_{i}B^{b}_{i+\hat{x}}+A^{b}_{i}A^{w}_{i+\hat{y}})$
(121)

We take the limit, $U>>t$ and develop a strong-coupling perturbative
expansion. The ground state subspace, in the limit, $U>>t$, corresponds to
either $\ket{n_{id}=0,n_{if}=1}$ or $\ket{n_{id}=1,n_{if}=0}$ at each lattice
site, $i$. It also satisfies, $A^{w}_{i}A^{b}_{i}B^{b}_{i}B^{w}_{i}=1$, i.e.
even Majorana parity sector. One can treat the local subspace as some
effective spin-1/2 degrees of freedom, which in the Majorana basis could be
represented by the following relations,
$\sigma^{z}_{i}=iB^{b}_{i}B^{w}_{i}=iA^{b}_{i}A^{w}_{i},\sigma^{x}_{i}=iB^{b}_{i}A^{w}_{i}=iB^{w}_{i}A^{b}_{i},\sigma^{y}_{i}=iA^{w}_{i}B^{w}_{i}=-iA^{b}_{i}B^{b}_{i}$.
The above representation is a faithful one under the constraint,
$A^{w}_{i}A^{b}_{i}B^{b}_{i}B^{w}_{i}=1$, which has to be satisfied in the
large $U$ limit. We define the projector onto the $2^{N}$ dimensional ground
state subspace,
$P=\frac{1}{2^{N}}\prod_{i=1}^{N}(1+A^{w}_{i}A^{b}_{i}B^{b}_{i}B^{w}_{i})$.
Under the action of $V_{x}$ and $V_{y}$, the Majorana fermion parities flip on
both the ends of the link $(i,i+\alpha)$. In order to determine the effective
low energy Hamiltonian, we have to figure out all possible ways of restoring
the GS parity, at various orders of the perturbation theory. One can verify
that the odd order terms don’t contribute because the GS parity restriction is
not maintained. In the second order, there are two contributions, (1) which
just gives a constant shift, this arise when we apply the same link hopping
twice to return to the GS parity sector. (2) Another contribution gives an
anti-ferromagnetic Ising coupling. For example, consider the following process
along the x-directions,

$\displaystyle\frac{1}{4}i(-i)(1-\lambda^{2})B^{b}_{i}B^{w}_{i+\hat{x}}B^{w}_{i}B^{b}_{i+\hat{x}}=-(1-\lambda^{2})\sigma^{z}_{i}\sigma^{z}_{i+x}$

The excitation energy to flip the parities at any two site is $4U$. So, in the
second order,

$\displaystyle
H_{eff}^{(2)}=\frac{(1-\lambda^{2})t^{2}}{8U}\sum_{i}\big{(}\sigma^{z}_{i}\sigma^{z}_{i+\hat{x}}+\sigma^{z}_{i}\sigma^{z}_{i+\hat{y}}\big{)}$
(122)

Similarly, the fourth-order effective Hamiltonian can be derived from

$\displaystyle H^{(4)}_{eff}=$ $\displaystyle
PV(1-P)\frac{1}{E_{0}-H_{0}}V(1-P)\frac{1}{E_{0}-H_{0}}$ $\displaystyle\times
V(1-P)\frac{1}{E_{0}-H_{0}}(1-P)VP$ (123)

Here, there are several contributions; these come from terms like (1):
$\tilde{V}_{x}^{4}$, (2): $\tilde{V}_{x}^{3}\tilde{V}_{y}$, (3):
$\tilde{V}_{x}\tilde{V}^{3}_{y}$, (4): $\tilde{V}_{x}^{2}\tilde{V}_{y}^{2}$,
(5): $\tilde{V}_{y}^{4}$, where the tilde denotes the projected operators. The
appropriate energy denominator will be multiplied later.

$\displaystyle(1)=\frac{t^{4}(1+\lambda)^{4}}{16}\sigma^{x}_{i}\sigma^{y}_{i+\hat{x}}\sigma^{x}_{i+\hat{x}+\hat{y}}\sigma^{y}_{i+\hat{y}}$
$\displaystyle(2)+(3)=-\frac{2}{16}t^{4}(1-\lambda^{4})[\sigma^{y}_{i}\sigma^{x}_{i+\hat{x}}\sigma^{x}_{i+\hat{x}+\hat{y}}\sigma^{y}_{i+\hat{y}}$
$\displaystyle\ \
+\sigma^{x}_{i}\sigma^{x}_{i+\hat{x}}\sigma^{y}_{i+\hat{x}+\hat{y}}\sigma^{y}_{i+\hat{y}}+\sigma^{y}_{i}\sigma^{y}_{i+\hat{x}}\sigma^{x}_{i+\hat{x}+\hat{y}}\sigma^{x}_{i+\hat{y}}$
$\displaystyle\ \
+\sigma^{x}_{i}\sigma^{y}_{i+\hat{x}}\sigma^{y}_{i+\hat{x}+\hat{y}}\sigma^{x}_{i+\hat{y}}]$
$\displaystyle(4)=-\frac{2}{16}t^{4}(1-\lambda^{4})[\sigma^{y}_{i}\sigma^{y}_{i+\hat{x}}\sigma^{y}_{i+\hat{x}+\hat{y}}\sigma^{y}_{i+\hat{y}}$
$\displaystyle\
+\sigma^{x}_{i}\sigma^{x}_{i+\hat{x}}\sigma^{x}_{i+\hat{x}+\hat{y}}\sigma^{x}_{i+\hat{y}}-\sigma^{y}_{i}\sigma^{x}_{i+\hat{x}}\sigma^{y}_{i+\hat{x}+\hat{y}}\sigma^{x}_{i+\hat{y}}]$
$\displaystyle(5)=\frac{t^{4}(1-\lambda)^{4}}{16}\sigma^{x}_{i}\sigma^{y}_{i+\hat{x}}\sigma^{x}_{i+\hat{x}+\hat{y}}\sigma^{y}_{i+\hat{y}}$

After adding all the above contributions (with the appropriate energy
denominators) and multiplying a combinatorial factor, which counts the
different possible orders of the link operators ($V_{x},V_{y}$) in the above
(1)-(5) terms, we finally arrive at the effective low-energy Hamiltonian (in
the large $U/t$) given by (58) in the main text. Notice that the link
operators all commute with each other, so their actual order of operation is
unimportant. This is in contrast to the Kitaev model, where such
simplifications do not happen.
When $\lambda\rightarrow 1$ (equal hopping and pairing limit), we obtain the
Wen’s plaquette model (which is the same as toric code).
To transform Wen’s plaquette model to the canonical toric code form (see Eqn.
(10)), we apply some local spin basis rotations. To perform it, we divide the
square lattice into the union of alternate plaquettes, which are labeled as
$a$ and $b$. On the “$a$” plaquettes , we rotate the spins on the sites
$i+\hat{x}$ and $i+\hat{z}$ by $R_{z}(\pi/2)$ and leave the two others, i.e.
$\sigma^{x}_{i+\hat{\alpha}}\rightarrow\sigma^{y}_{i+\hat{\alpha}}$ and
$\sigma^{y}_{i+\hat{\alpha}}\rightarrow-\sigma^{x}_{i+\hat{\alpha}}$, where
$\alpha=x,y$. Similarly, on the “$b$” plaquettes, we apply $R_{z}(\pi/2)$ on
the sites $i$ and $i+\hat{d}$. Basically, the rotations have to be performed
on two diagonally opposite points. Therefore,

$\displaystyle H_{TC}=-J_{w}\sum_{a}\prod_{i\in
a}\sigma^{x}_{i}-J_{w}\sum_{b}\prod_{i\in b}\sigma^{y}_{i}$ (124)

Where $J_{w}=\frac{3t^{4}}{32U^{3}}(3\lambda^{2}+\lambda^{4})$.
When, $\lambda\rightarrow 0$ (zero pairing), The $U(1)$ symmetry is restored,

$\displaystyle H_{\lambda=0}=J_{I}\sum_{\langle
ij\rangle}\sigma^{z}_{i}\sigma^{z}_{j}-J_{p}\sum_{i}\big{(}\sigma^{+}_{i}\sigma^{-}_{i+\hat{x}}\sigma^{+}_{i+\hat{x}+\hat{y}}\sigma^{-}_{i+\hat{y}}+h.c.\big{)}$

Here, $J_{I}=\frac{t^{2}}{8U}$, and $J_{p}=\frac{3t^{4}}{32U^{3}}$. The sign
of the ring-exchange coupling term has been flipped without affecting the
signature of the Ising coupling. We have rotated the spins on alternate sub-
lattices by $\pi/2$ with respect to $z$ axis. So,
$\sigma^{+}_{i}\rightarrow-i\sigma^{+}_{i}$ and
$\sigma^{+}_{i+\hat{x}+\hat{y}}\rightarrow-i\sigma^{+}_{i+\hat{x}+\hat{y}}$.

##  Appendix G HQCM with staggered Zeeman field

Consider the HQCM Hamiltonian with a staggered Zeeman magnetic field coupling,

$\displaystyle H=J_{x}\sum_{i\in b}(S^{x}_{i,b}S^{x}_{i+\hat{x},w}+$
$\displaystyle S^{x}_{i,b}S^{x}_{i-\hat{x},w})+J_{z}\sum_{i\in
b}S^{z}_{i,b}S^{z}_{i+\hat{z},w}$
$\displaystyle+h_{z}\sum_{i,\alpha=b,w}(-1)^{i_{\alpha}}S^{z}_{i,\alpha}$
(125)

Here, we choose the negative sign for the Zeeman term on the $w$ sites, i.e.
$i_{w}=1,i_{b}=0$. On the dual square lattice (joined by zz bonds of the brick
wall lattice), the Hamiltonian could be re-written as follows,

$\displaystyle
H=J_{x}\sum_{i}(S^{x}_{i,b}S^{x}_{i+\hat{e}_{x},w}+S^{x}_{i,w}S^{x}_{i+\hat{e}_{y},b})$
$\displaystyle\ \ \
+J_{z}\sum_{i}S^{z}_{i,b}S^{z}_{i,w}+h_{z}\sum_{i}(S^{z}_{i,b}-S^{z}_{i,w})$
(126)

Under the JW transformation and the subsequent Majorana fermion decompositions
(as described at the beginning of the main text, see Eqs.(2)-(4)),
$S^{z}_{i,b}=-(i/2)B^{b}_{i}A^{b}_{i}$ and
$S^{z}_{i,w}=-(i/2)A^{w}_{i}B^{w}_{i}$. So, the Zeeman term transforms to the
following,

$\displaystyle
H_{z}=\frac{h_{z}}{2}\sum_{i}(iA^{w}_{i}B^{w}_{i}-iB^{b}_{i}A^{b}_{i})$ (127)

It commutes with the $J_{z}\rightarrow\infty$ ground state subspace projectors
$P\sim\prod_{i=1}^{N}(1+A^{w}_{i}A^{b}_{i}B^{b}_{i}B^{w}_{i})$, i.e.
$[P,H_{z}]=0$. Therefore, it should not affect the structure of the strong-
coupling ($J_{z}>>J_{x}$) Hamiltonian in the absence of $h_{z}$. The effective
strong-coupling Hamiltonian after including $h_{z}$ now becomes,

$\displaystyle
H^{s}_{TC}=-J_{w}\sum_{P}\sigma^{x}_{i}\sigma^{y}_{i+x}\sigma^{x}_{i+x+y}\sigma^{y}_{i+y}-h\sum_{i}(-1)^{i}\sigma^{y}_{i}$
(128)

Here, the effective spins are defined as earlier (see Appendix F), and the
Zeeman term act as a staggered field along $y$. We now rotate the spins of one
sub-lattice (by $\pi$), $\sigma^{y}\rightarrow-\sigma^{y}$ and
$\sigma^{z}\rightarrow-\sigma^{z}$; this transforms the staggered coupling to
a translationally symmetric term,

$\displaystyle
H^{s}_{TC}=-J_{w}\sum_{P}\sigma^{x}_{i}\sigma^{y}_{i+x}\sigma^{x}_{i+x+y}\sigma^{y}_{i+y}-h\sum_{i}\sigma^{y}_{i}$
(129)

This model belongs to the critical universality class of $1d$ TFIM.
Instead of using the JW fermionization in the intermediate step, we can
perform the strong-coupling expansion directly on the spin Hamiltonian (125).
The conclusion should not differ. The ground state (GS) subspace, when
$J_{z}>>J_{x},h_{z}$ belongs to
$\\{\ket{\uparrow\downarrow}_{i},\ket{\downarrow\uparrow}_{i}\\}$. Since the
model has local symmetries $S^{x}_{i,w}S^{x}_{i,b}$ (when $h_{z}=0$), we
choose the symmetric and anti-symmetric combinations,
$\ket{\phi^{\pm}_{i}}\sim\ket{\uparrow\downarrow}_{i}\pm\ket{\downarrow\uparrow}_{i}$
as our local basis states. The projector onto this subspace is
$P=2^{-N_{b}}\prod_{i\in N_{b}}(1-\sigma^{z}_{i,b}\sigma^{z}_{i,w})$. Notice
that $[H_{z},P]=0$ (so it will remain unchanged) but $[H_{x},P]\neq 0$. Here,
we call the Hamiltonian involving $J_{x}$ couplings as $H_{x}$. When $H_{x}$
acts on the GS subspace, it creates transitions between
$\\{\ket{\phi^{\pm}_{i}}\\}$ and
$\\{\ket{\psi^{\pm}_{i}}\sim(\ket{\uparrow\uparrow}_{i}\pm\ket{\downarrow\downarrow}_{i})\\}$
on the two neighbouring zz links, $(i,i+\alpha)$, ($\alpha=x,y$). All the odd-
order perturbative corrections are zero. The second-order correction is just a
constant, because one has to unflip the flipped spins to return to the same GS
space.
The fourth order process has the following structure (we show here one
particular case): Take the following initial states on a plaquette,
$\ket{\phi^{+}}_{i}$, $\ket{\phi^{-}}_{i+x}$, $\ket{\phi^{+}}_{i+x+y}$, and
$\ket{\phi^{-}}_{i+y}$. Then,

$(1)\ \sigma^{x}_{i,b}\sigma^{x}_{i+x,w}\
\\{\ket{\phi^{+}}_{i},\ket{\phi^{-}}_{i+x}\\}\ \
\rightarrow\\{\ket{\psi^{+}}_{i},\ket{\psi^{-}}_{i+x}\\}$
$(2)\ \sigma^{x}_{i+x,w}\sigma^{x}_{i+x+y,b}\
\\{\ket{\psi^{-}}_{i+x},\ket{\phi^{+}}_{i+x+y}\\}\rightarrow\\{\ket{\phi^{-}}_{i+x},\ket{\psi^{+}}_{i+x+y}\\}$
$(3)\ \sigma^{x}_{i+x+y,w}\sigma^{x}_{i+y,b}\
\\{\ket{\psi^{+}}_{i+x+y},\ket{\phi^{-}}_{i+y}\\}\rightarrow\\{\ket{\phi^{+}}_{i+x+y},\ket{\psi^{-}}_{i+y}\\}$
$(4)\ \sigma^{x}_{i,w}\sigma^{x}_{i+y,b}\
\\{\ket{\psi^{+}}_{i},\ket{\psi^{-}}_{i+y}\\}\rightarrow\\{\ket{\phi^{+}}_{i},\ket{\phi^{-}}_{i+y}\\}$

The staggered field term ($\sim h_{z}$) causes transition from
$\ket{\phi^{+}}_{i}$ to $\ket{\phi^{-}}_{i}$ (and vice-versa). Define the
effective (bond) spins,
$\tau^{x}_{i}=\ket{\uparrow\downarrow}\bra{\downarrow\uparrow}+\ket{\downarrow\uparrow}\bra{\uparrow\downarrow}$,
and
$\tau^{z}_{i}=\ket{\uparrow\downarrow}\bra{\uparrow\downarrow}-\ket{\downarrow\uparrow}\bra{\downarrow\uparrow}$,
the fourth-order effective Hamiltonian obtains,

$\displaystyle
H_{eff}=-J_{w}\sum_{i}\tau^{z}_{i}\tau^{z}_{i+x+y}-h_{z}\sum_{i}\tau^{x}_{i}$
(130)

Here $J_{w}\sim O(J^{4}_{x}/J_{z}^{3})$. This is just a collection of
decoupled TFIM chains along the diagonal of the dual square lattice. The
domain wall excitations of these Ising chains corresponds to the anyons of the
Wen’s model [73] which are dispersive due to finite $h_{z}$. Both (129) and
(130) fall in the same $1d$ quantum Ising universality class.

## References

  * [1]
  * Xiao-Gang Wen [2000] X. G. Wen,  _Quantum Field Theory of Many-Body Systems: From the Origin of Sound to an Origin of Light and Electrons,_ Oxford Scholarship Online (2010).
  * L. Savary, and Leon Balents [2016b] L. Savary, and L. Balents, Rep. Prog. Phys. 80, 016502 (2016b).
  * Xiao-Gang Wen [2003b] X-G. Wen, Phys. Rev. Lett. 90, 016803 (2003b).
  * Michael Levin, and Xiao-Gang Wen [2006b] M. Levin, and X-G. Wen, Phys. Rev. Lett. 96, 110405 (2006b).
  * Alexei Kitaev, and John Preskill [2006b] A. Kitaev, and J. Preskill, Phys. Rev. Lett. 96, 110404 (2006b).
  * A. Kitaev [2006b] A. Kitaev, Annals of Physics 321, 2–111 (2006b).
  * Yogesh Singh, Manni, S., ,Reuther, J., Berlijn, T., Thomale, R., Ku, W., Trebst, S., and Gegenwart, P. [2012b] Yogesh Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Phys. Rev. Lett. 108, 127203 (2012b).
  * Do, Seung-Hwan, Park, Sang-Youn, Yoshitake, Junki _et al._ [2017b] S-H. Do, S-Y. Park, and J. Yoshitake,  _et al._ Nature Physics 13, 1079–1084 (2017b).
  * Banerjee, A., Lampen-Kelley, P., Knolle, J. _et al._ [2018b] A. Banerjee, P. Lampen-Kelley, and J. Knolle,  _et al._ npj Quantum Materials 3, 8 (2018b).
  * You, J. Q., Shi, Xiao-Feng, Hu, Xuedong, and Nori, Franco [2010b] J. Q. You, Xiao-Feng Shi, X. Hu, and F. Nori, Phys. Rev. B 81, 014505 (2010b).
  * M. Barkeshli, and J. D. Sau [2015b] M. Barkeshli, and J. D. Sau, arXiv:1509.07135 (2015b).
  * Bacon, Dave [2006b] D. Bacon Phys. Rev. A 73, 012340 (2006b).
  * Armin Rahmani, Xiaoyu Zhu, Marcel Franz, and Ian Affleck [2015b] A. Rahmani, X. Zhu, M. Franz, and I. Affleck, Phys. Rev. Lett 115, 166401 (2015b).
  * Ian Affleck, Armin Rahmani, and Dmitry Pikulin [2017b] I. Affleck, A. Rahmani, and D. Pikulin, Phys. Rev. B 96, 125121 (2017b).
  * Y. Kamiya, A. Furusaki, J. C. Y. Teo, and G. W. Chern [2018b] Y. Kamiya, A. Furusaki, J. C. Y. Teo, and G. W. Chern, Phys. Rev. B 98, 161409(R) (2018b).
  * Jason Alicea [2012b] J. Alicea Rep. Prog. Phys. 75, 076501 (2012b).
  * Alexander Ziesen, Fabian Hassler, and Ananda Roy [2019b] A. Ziesen, F. Hassler, and A. Roy, Phys. Rev. B 100, 104508 (2019b).
  * Y You, Daniel Litinski, and Felix von Oppen [2019b] Y. You, D. Litinski, and F. von Oppen, Phys. Rev. B 100, 054513 (2019b).
  * J Fu [2018b] J. Fu Phys. Rev. B 98, 214432 (2018b).
  * Zohar Nussinov, Gerardo Ortiz, and Emilio Cobanera [2012b] Z. Nussinov, G. Ortiz, and E. Cobanera, Phys. Rev. B 86, 085415 (2012b).
  * Zhao, Erhai, and Liu, W. Vincent [2008b] E. Zhao, and W. V. Liu, Phys. Rev. Lett. 100, 160403 (2008b).
  * Nussinov, Zohar, and van den Brink, Jeroen [2015b] Z. Nussinov, and J. van den Brink, Rev. Mod. Phys. 87, 1–59 (2015b).
  * Leon Balents, Matthew P. A. Fisher, and Chetan Nayak [1999b] L. Balents, M. P. A. Fisher, and C. Nayak, Phys. Rev. B 60, 1654 (1999b).
  * A. Rüegg, S. D. Huber, and M. Sigrist [2010b] A. Rüegg, S. D. Huber, and M. Sigrist, Phys. Rev. B 81, 155118 (2010b).
  * L. de’Medici, A. Georges, and S. Biermann [2005b] L. de’Medici, A. Georges, and S. Biermann, Phys. Rev. B 72, 205124 (2005b).
  * Eduardo Fradkin, and Stephen H. Shenker [1979b] E. Fradkin, and S. H. Shenker, Phys. Rev. D 19, 3682 (1979b).
  * G. Wendin, and V. S. Shumeiko [2005b] G. Wendin, and V. S. Shumeiko, arXiv:cond-mat/0508729v1 (2005b).
  * Siewert, J., Fazio, R., Palma, G.M. _et al._ [2000b] J. Siewert, R. Fazio, and G. M. Palma,  _et al._ Journal of Low Temperature Physics 118, 795–804 (2000b).
  * Y. M. Vilk and A. MS. Tremblay [1997b] Y. M. Vilk, and A. MS. Tremblay, J. Phys. I France 7, 11 (1997b).
  * Mattis and Nam [1972b] D. C. Mattis and S. B. Nam, Journal of Mathematical Physics 13, 1185-1189 (1972b).
  * M. P. A. Fisher [2004b] M. P. A. Fisher, _Duality in low dimensional quantum field theories._, Strong interactions in low dimensions, Physics and Chemistry of Materials with Low-Dimens, eds. D. Baeriswyl, and D. Degiorgi, vol. 25, (Springer, Dordrecht, 2004).
  * P Pfeuty and R J Elliott [1971b] P. Pfeuty and R. J. Elliott, Journal of Physics C: Solid State Physics 4, 2370–2385 (1971b).
  * [34] To get the four robust corner modes, we require additional edge perturbations (that hybridizes the MZMs along the edges) to make the boundary edges fully gapped, leaving unpaired MZMs only at the high-symmetry corners of the lattice.
  * Vladimir Dobrosavljevic  _et. al._ [2012] _Conductor-Insulator Quantum Phase Transitions_, eds. V. Dobrosavljevic, N. Trivedi, and James M. Valles, (Oxford Scholarship Online, 2012).
  * [36] Since the JW fermion model Eq.(6) explicitly contains a BCS Cooper-pairing term, the global charge $U(1)$ symmetry is reduced to $Z_{2}$, which corresponds to total fermion number parity conservation of the system..
  * [37] Paramagnetic phase of Ising models can be described as off-diagonal long-range order of some appropriate disorder operators (observables), see Ref.[32]..
  * Steven A. Kivelson, Daniel S. Rokhsar, and James P. Sethna [1987b] Steven A. Kivelson, D. S. Rokhsar, and J. P. Sethna, Phys. Rev. B 35, 8865(R) (1987b).
  * D. S. Rokhsar and Steven A. Kivelson [1988b] D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. 61, 2376 (1988b).
  * I. S. Tupitsyn, A. Kitaev, N. V. Prokof’ev, and P. C. E. Stamp [2010b] I. S. Tupitsyn, A. Kitaev, N. V. Prokof’ev, and P. C. E. Stamp, Phys. Rev. B 82, 085114 (2010b).
  * R. Moessner, S. L. Sondhi, and Eduardo Fradkin [2001b] R. Moessner, S. L. Sondhi, and E. Fradkin, Phys. Rev. B 65, 024504 (2001b).
  * Rodolfo A. Jalabert and Subir Sachdev [1991b] R. A. Jalabert and S. Sachdev, Phys. Rev. B 44, 686 (1991b).
  * R. Moessner and S. L. Sondhi [2001b] R. Moessner and S. L. Sondhi, Phys. Rev. B 63, 224401 (2001b).
  * Thierry Giamarchi [2003] T. Giamarchi,  _Quantum Physics in One Dimension,_ Oxford University Press (2003).
  * T. Senthil and Matthew P. A. Fisher [2000b] T. Senthil and M. P. A. Fisher, Phys. Rev. B 62, 7850 (2000b).
  * Jorge V. José, L. P. Kadanoff, Scott Kirkpatrick, and David R. Nelson [1977b] J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Phys. Rev. B 16, 1217 (1977b).
  * C. Dasgupta and B. I. Halperin [1981b] C. Dasgupta and B. I. Halperin, Phys. Rev. Lett. 47, 1556 (1981b).
  * K. I. Kugel and D. I. Khomskii [1982b] K. I. Kugel and D. I. Khomskii, Soviet Physics Uspekhi 25, 231–256 (1982b).
  * Lacroix  _et. al._ [2011] _Introduction to Frustrated Magnetism. Materials, Experiments, Theory_, eds. C. Lacorix, P. Mendels, and F. Mila, Springer Series in Solid State Sciences, vol. 164, (Springer-Verlag Berlin Heidelberg, 2011).
  * Hajime Yamamoto, Sachiko Kamiyama, Ikuya Yamada, and Hiroyuki Kimura [2022b] H. Yamamoto, S. Kamiyama, I. Yamada,, and H. Kimura, J. Am. Chem. Soc. 144, 1082–1086 (2022b).
  * Park, J., Tan, TY., Adroja, D. _et al._ [2016b] J. Park, TY. Tan, and D. Adroja,  _et al._ Sci Rep 6, 25238 (2016b).
  * Maekawa, Sadamichi, Tohyama, Takami, Barnes, Stewart Edward, Ishihara, Sumio, Koshibae, Wataru, and Khaliullin, Giniyat [2004b] S. Maekawa, T. Tohyama, and S. E. Barnes  _et al._ ,  _Physics of transition metal oxides_, vol. 144, (Springer Science & Business Media, 2004).
  * Schlittler, Thiago M., Mosseri, Rémy, and Barthel, Thomas [2017b] T. M. Schlittler, R. Mosseri, and T. Barthel, Phys. Rev. B 96, 195142 (2017b), and references therein..
  * P W Anderson, P A Lee, M Randeria, T M Rice, N Trivedi, and F C Zhang [2004b] P. W. Anderson, P. A. Lee, M. Randeria, T. M. Rice, N. Trivedi, and F. C. Zhang, Journal of Physics: Condensed Matter, IOP Publishing 16, R755–R769 (2004b).
  * Nandkishore, R., Levitov, L., and Chubukov, A. [2012b] R. Nandkishore, L. Levitov, and A. Chubukov, Nature Phys 8, 158–163 (2012b).
  * Vojta, Matthias [2002b] M. Vojta, Phys. Rev. B 66, 104505 (2002b).
  * Pathak, Sandeep, Shenoy, Vijay B., and Baskaran, G. [2010b] S. Pathak, V. B. Shenoy, and G. Baskaran, Phys. Rev. B 81, 085431 (2010b).
  * G. Baskaran [2002b] G. Baskaran, Phys. Rev. B 65, 212505 (2002b).
  * Abanov, A. G. and Wiegmann, P. B. [1997b] A. G. Abanov and P. B. Wiegmann, Phys. Rev. Lett. 78, 4103–4106 (1997b).
  * Wang, Zhan, Zhou, Sen, Chen, Weiqiang, and Zhang, Fu-Chun [2020b] Z. Wang, S. Zhou, W. Chen, and F-C. Zhang, Phys. Rev. B 101, 180509(R) (2020b).
  * A. Yu. Kitaev [2003b] A. Yu. Kitaev, Annals of Physics 303, 2-30 (2003b).
  * P. W. Anderson [1987b] P. W. Anderson, Science 235, 1196-1198 (1987b).
  * T. Moriya [2012] T. Moriya,  _Spin Fluctuations in Itinerant Electron Magnetism,_ Springer Series in Solid-State sciences, 56 (2012).
  * N. D. Mermin and H. Wagner [1966b] N. D. Mermin and H. Wagner, Phys. Rev. Lett 17, 1133 (1966b).
  * Chakravarty, Sudip, Halperin, Bertrand I., and Nelson, David R. [1989b] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. B 39, 2344–2371 (1989b).
  * Subir Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions,_ 2nd edition, Cambridge University Press, Cambridge (2011).
  * S. Allen and A.-M. S. Tremblay [2001b] S. Allen and A.-M. S. Tremblay, Phys. Rev. B 64, 075115 (2001b).
  * B. Davoudi and A.-M. S. Tremblay [2006b] B. Davoudi and A.-M. S. Tremblay, Phys. Rev. B 74, 035113 (2006b).
  * Dominic Bergeron, Vasyl Hankevych, Bumsoo Kyung, and A.-M. S. Tremblay [2011b] D. Bergeron, V. Hankevych, B. Kyung, and A.-M. S. Tremblay, Phys. Rev. B 84, 085128 (2011b).
  * Wright, Katherine [2021b] Katherine Wright, Physics 14, 159 (2021b).
  * Bruun, G. M and Baym, Gordon [2006b] G. M. Bruun and G. Baym, Phys. Rev. A 74, 033623 (2006b).
  * Hart, R., Duarte, P., Yang, TL. _et al._ [2015b] R. Hart, P. Duarte, and TL. Yang,  _et al._ Nature 519, 211–214 (2015b).
  * Han-Dong Chen, and Nussinov, Zohar [2008b] H-D. Chen, and Z. Nussinov, J. Phys. A: Math. Theor. 41, 075001 (2008b).
  * Chen and Hu [2007b] H. D. Chen and J. Hu, Phys. Rev. B 76, 193101 (2007b).
  * Vidal  _et al._ [2005b] J. Vidal, R. Thomale, K. P. Schmidt, and S. Dusuel, Phys. Rev. B 80, 081104(R) (2009).
  * S. Sur, M. S. Laad, and S. R. Hassan [2021b] S. Sur, M. S. Laad, and S. R. Hassan, Phys. Rev. B 103, 144419 (2021b).
  * Brzezicki and Andrzej [2005b] W. Brzezicki and A. M. Oleś, Phys. Rev. B 82, 060401(R) (2010b).
  * Paramekanti, Arun, Balents, Leon, and Fisher, Matthew P. A. [2002b] A. Paramekanti, L. Balents, and M. P. A. Fisher, Phys. Rev. B 66, 054526 (2002b).
  * S. Doniach and D. Das [2003b] S. Doniach and D. Das, Braz. J. Phys. 33, 4 (2003b).
  * Paredes, Belén and Bloch, Immanuel [2008b] B. Paredes and I. Bloch, Phys. Rev. A 77, 023603 (2008b).
  * Dai, HN., Yang, B., Reingruber, A. _et al._ [2017b] HN. Dai, B. Yang, and A. Reingruber  _et al._ Nature Phys 13, 1195–1200 (2017b).
  * Sandvik, A. W., Daul, S., Singh, R. R. P., and Scalapino, D. J. [2002b] A. W. Sandvik, S. Daul, R. R. P. Singh, and D. J. Scalapino, Phys. Rev. Lett. 89, 247201 (2002b).
  * Tay, Tiamhock and Motrunich, Olexei I. [2011b] T. Tay and O. I. Motrunich, Phys. Rev. B 83, 205107 (2011b).
  * Ma, Han and Pretko, Michael [2018b] H. Ma and M. Pretko, Phys. Rev. B 98, 125105 (2018b).
  * Anne-Marie Daré,, Vilk, Y. M., and Tremblay, A. -M. S. [1996b] Anne-Marie Daré, Y. M. Vilk, and A. M. S. Tremblay, Phys. Rev. B 53, 14236–14251 (1996b).
  * Wen-Long You, Yimin Wang, Tian-Cheng Yi, Chengjie Zhang, and Andrzej M. Oleś [2018b] W-L. You, Y. Wang, T-C. Yi, C. Zhang, and A. M. Oleś, Phys. Rev. B 97, 224420 (2018b).
  * Jing You, Su-Peng Kou, and Xiao-Gang Wen [2008b] J. You, S-P. Kou, and X-G. Wen, EPL 84, 17004 (2008b).