# Variational Tensor Wavefunctions for the Interacting Quantum Spin Hall Phase

Yixin Ma  Shenghan Jiang<EMAIL_ADDRESS>Chao Xu<EMAIL_ADDRESS>Kavli
Institute for Theoretical Sciences, University of Chinese Academy of Sciences,
Beijing 100190, China

###### Abstract

The quantum spin hall (QSH) phase, also known as the 2D topological insulator,
is characterized by protected helical edge modes arising from time reversal
symmetry. While initially proposed for band insulators, this phase can also
manifest in strongly-correlated systems where conventional band theory fails.
To overcome the challenge of simulating this phase in realistic correlated
models, we propose a novel framework utilizing fermionic tensor network
states. Our approach involves constructing a tensor representation of the
fixed-point wavefunction based on an exact solvable model, enabling us to
derive a set of tensor equations governing the transformation rules of local
tensors under symmetry operations. These tensor equations lead to the
anomalous edge theory, which provides a comprehensive description of the QSH
phase. By solving these tensor equations, we obtain variational ansatz for the
QSH phase, which we subsequently verify through numerical calculations. This
method serves as an initial step towards employing tensor algorithms to
simulate the QSH phase in strongly-correlated systems, opening new avenues for
investigating and understanding topological phenomena in complex materials.

_Introduction.–_ The discovery of the quantum spin Hall (QSH) phaseKane and
Mele (2005) has sparked research interest in studying the interplay between
symmetry and topology in quantum materialsQi and Zhang (2011); Hasan and Kane
(2010). Initially proposed as a topological band insulator, the QSH phase is
characterized by stable properties such as anomalous helical edge modes and
topological response to electromagnetic fieldsQi  _et al._ (2008a). It has
been found that the QSH phase can also be realized as a Mott insulator in
strongly-correlated systems, representing an example of interacting fermionic
symmetry-protected topological (SPT) phasesFidkowski and Kitaev (2010, 2011);
Chen  _et al._ (2011); Gu and Wen (2014); Gu and Levin (2014); Wang and Gu
(2018, 2020); Wang  _et al._ (2017); Cheng  _et al._ (2018); Wang  _et al._
(2014); Wang and Senthil (2014); Senthil (2015); Witten (2016); Freed and
Hopkins (2021). Solvable models based on commuting-projector Hamiltonians have
been used to construct various interacting fermionic SPT phases, including the
QSH phaseWare  _et al._ (2016); Tarantino and Fidkowski (2016); Wang  _et al._
(2018); Metlitski (2019); Son and Alicea (2019); Wang  _et al._ (2021).
However, these models only provide fixed-point wavefunctions and are hardly
useful for numerical simulations.

To construct generic variational wavefunctions beyond the fixed point, we turn
to fermionic tensor networksBarthel  _et al._ (2009); Corboz  _et al._ (2010);
Kraus  _et al._ (2010); Gu  _et al._ (2010); Schuch  _et al._ (2011); Bultinck
_et al._ (2017a, b); Wille  _et al._ (2017); Cirac  _et al._ (2021). Our
strategy is presented as following. Motivated by the interacting edge theory,
we introduce the fixed-point wavefunction proposed in Ref. Wang  _et al._
(2021), and then translate it to the tensor network state. With such tensor
representation, we extract a set of tensor equations for symmetry actions on
tensors. From tensor equations, we obtain algebraic data characterizing the
anomalous edge theory of the QSH phase. Finally, we apply our method to a
spin-1/2 fermionic system on honeycomb and square lattice: by listing and
solving tensor equations, we get variational ansatz for the QSH phase on such
systems. We further calculate topological invariantsShapourian  _et al._
(2017); Shiozaki  _et al._ (2018) of the variational wavefunction based on our
ansatz to show a parameter region of QSH phase which can be used for numerical
simulations.

_Interacting edge theory.–_ The QSH phase hosts charge conservation symmetry
generated by $n_{f}$ and time reversal symmetry $\mathcal{T}$, where

$\displaystyle\mathcal{T}^{2}=\exp[\operatorname{i}\pi n_{f}]\equiv
F~{},\quad\mathcal{T}\cdot n_{f}\cdot\mathcal{T}^{-1}=n_{f}~{}.$ (1)

Here, $F$ is the fermion parity operator.

To get intuition about the interacting bulk wavefunction, we start from its
anomalous edge states, which is described by massless helical Dirac fermions:

$\displaystyle H_{edge}$
$\displaystyle=\int\operatorname{d}\\!x\,(-\operatorname{i}v_{F})\left[\psi_{R}^{\dagger}(x)\partial_{x}\psi_{R}(x)-\psi_{L}^{\dagger}(x)\partial_{x}\psi_{L}(x)\right]~{},$

where $\psi_{L/R}$ is the left/right moving fermion mode, and $v_{F}$ the
fermion velocity. $\mathcal{T}$ acts as
$\psi_{R/L}\to\pm\operatorname{i}\psi_{L/R}$, forbidding mass terms opening a
gap.

The interacting edge theory can be analyzed by the bosonization methodHaldane
(1981). Conjugate fields $\phi(x)$ and $\theta(x)$ are introduced, both with
periodicity $2\pi$, where
$[\partial_{x}\theta(x),\phi(x^{\prime})]=2\pi\operatorname{i}\delta(x-x^{\prime})$.
With these hydrodynamic variables,
$\psi_{R/L}(x)\sim\exp[-(\operatorname{i}\phi(x)\pm\operatorname{i}\theta(x)/2)]$,
charge density $\delta\rho(x)=-\partial_{x}\theta(x)/2\pi$, and current
density $j(x)=\partial_{t}\theta(x)/2\pi$. Symmetry actions on $\theta$ and
$\phi$ are derived from its action on $\psi_{R/L}$, where

$\displaystyle U(\varphi):{}$
$\displaystyle\phi\to\phi+\varphi~{},\quad\theta\to\theta~{};$
$\displaystyle\mathcal{T}:{}$
$\displaystyle\phi\to-\phi~{},\quad\theta\to\theta+\pi~{},\quad\operatorname{i}\to-\operatorname{i}~{}.$
(2)

Lagrangian density for the interacting edge theory isWu  _et al._ (2006); Xu
and Moore (2006)

$\displaystyle\mathcal{L}_{edge}={}$
$\displaystyle\frac{1}{2\pi}\partial_{x}\theta\partial_{t}\phi-\frac{v_{F}}{4\pi}\left(\frac{1}{K}(\partial_{x}\theta)^{2}+K(\partial_{x}\phi)^{2}\right)$
$\displaystyle+\alpha\cos(2\theta-2\theta_{0})$ (3)

where $K$ is the Luttinger parameter, and for the non-interacting case $K=2$.
Due to Eq. (2), the most relevant symmetric scattering term is
$\alpha\cos(2\theta-2\theta_{0})$ with scaling dimension $2K$. It becomes
relevant when $K<1$, driving edge to a gapped phase. For the classical limit
where $\alpha\ll 0$, ground states are doubly degenerate, characterized by
$\langle\theta\rangle=\theta_{0}$ and $\langle\theta\rangle=\theta_{0}+\pi$
respectively. Note that these two states are exchanged under $\mathcal{T}$,
and thus spontaneously break $\mathcal{T}$ symmetry.

Topological defects of such edge symmetry breaking phase host anomalous
properties. We consider a time reversal domain wall at $x_{0}$, with domains
$\langle\theta(x<x_{0}-\epsilon)\rangle=\theta_{0}$ and
$\langle\theta(x>x_{0}+\epsilon)\rangle=\theta_{0}+\pi$, as shown in Fig. 1.
For region $(x_{0}-\epsilon,x_{0}+\epsilon)$, $\theta$ rotate
clockwise/counter-clockwise. Such domain wall carries $\pm 1/2$
chargeGoldstone and Wilczek (1981), as

$\displaystyle\int_{x_{0}-\epsilon}^{x_{0}+\epsilon}\operatorname{d}\\!x\,\delta\rho(x)=\int_{x_{0}-\epsilon}^{x_{0}+\epsilon}\operatorname{d}\\!x\,\left(-\frac{\partial_{x}\theta(x)}{2\pi}\right)=\pm\frac{1}{2}~{}.$
(4)

_The fixed-point wavefunction.-_ We now extend $\theta$-field to bulk. The
clockwise/counter-clockwise domain wall at edge is identified as $\pm 1/2$
vortex, as shown in Fig. 1. The half-charge edge domain wall motivates a
decorated vortex pictureChen  _et al._ (2014); Liu  _et al._ (2014): each
vortex core carries fermions with $n_{f}=n_{v}$, where $n_{v}$ is the winding
number. Note that $\mathcal{T}$ symmetry can be recovered by proliferating
vortices, and $n_{f}$ is conserved during this process due to conservation of
total vorticity.

Figure 1: Edge domain walls and bulk vortex of $\theta$-field with fermion
decoration.

With the decorated vortex picture, we introduce the fixed-point
wavefunctionWang  _et al._ (2021). As in Fig. 2, we consider a system with
spin-$1/2$ fermions $f_{\sigma}$’s at a honeycomb lattice, and Ising spins
$\lvert\tau\rangle$’s at the dual triangular lattice, where
$\sigma,\tau=\uparrow/\downarrow$. $\mathcal{T}$ flips both spins:

$\displaystyle\mathcal{T}:\lvert\uparrow\rangle\leftrightarrow\lvert\downarrow\rangle~{},~{}~{}f_{\sigma}\to\sigma^{y}_{\sigma\sigma^{\prime}}f_{\sigma^{\prime}}~{},~{}~{}\operatorname{i}\to-\operatorname{i}~{}.$
(5)

Figure 2: Configuration of the QSH phase’s fixed-point wavefunction.
$f_{\sigma}$ occupy the honeycomb lattice, while $|\tau\rangle$ on the dual
lattice. Crossing an Ising domain wall along/against oriented bonds of the
dual lattice introduces a $\pm\pi$ phase shift. Spins of fermions adhere to
the majority rule.

Here, Ising spins represent $\theta$-field, which rotate $\pm\pi$ when
crossing an Ising domain wall along/against the arrow on the bond. For arrows
in Fig. 2, an Ising domain wall going through site $(\bm{\mathbf{r}},u/v)$
leads to $n_{v}=\pm 1$ at this site. To match $n_{v}$, fermions at site
$(\bm{\mathbf{r}},u/v)$ are holes/electrons:

$\displaystyle[n_{f},f_{\bm{\mathbf{r}}s,\sigma}]=-(-1)^{s}f_{\bm{\mathbf{r}}s,\sigma}~{},$
(6)

where $(-1)^{s}=\pm 1$ for $s=u/v$. The fermion spin is enforced to follow the
majority Ising spins at adjacent plaquettes, as shown in Fig. 2. With such
majority rule, one can check that for each domain wall loop, number of
$f_{\uparrow}$ differs from number of $f_{\downarrow}$ by $\pm 6$. Let $c$ be
an Ising spin configuration and $\lvert\psi_{c}\rangle$ the corresponding
decorated fermion state, we have

$\displaystyle\mathcal{T}\lvert\psi_{c}\rangle=(-1)^{N_{dw}(c)}\lvert\psi_{\mathcal{T}c}\rangle~{},$
(7)

where $N_{dw}$ is the number of domain wall loops in $c$. The fixed-point
wavefunction is expressed asWang  _et al._ (2021)

$\displaystyle\lvert\Psi\rangle=\sum_{c}\Psi(c)\lvert
c\rangle\otimes\lvert\psi_{c}\rangle~{},$ (8)

where $\Psi(c)=\pm 1$ satisfies
$\Psi(c)=(-1)^{N_{dw}(c)}\Psi(\mathcal{T}c)$111 We mention that to fully
determine $\lvert\Psi\rangle$, we should explicitly write down the fermion
order and set the $\pm 1$ phase for each $\Psi(c)$ However, we will not
present such information as it is quite complicated and is unnecessary for the
following discussion. .

_Tensor network representation.-_ Constructing variational wavefunctions
beyond Eq. (8) is highly desirable for practical purposes. In the following,
we present a comprehensive framework based on fermionic projected entangled-
pair states (fPEPS). FPEPS are constructed using fermionic tensors, which are
quantum states residing in the fermionic tensor product ($\otimes_{f}$) of
physical and internal legs. The legs with inward/outward arrows correspond to
fermionic Hilbert spaces of ket/bra states, respectively. Fermionic tensor
contraction $\mathrm{fTr}$ are implemented by connecting outward and inward
internal legs, defined as

$\displaystyle\mathrm{fTr}[\langle i\rvert\otimes_{f}\lvert
j\rangle]=(-)^{\left\lvert i\right\rvert\left\lvert
j\right\rvert}\mathrm{fTr}[\lvert j\rangle\otimes_{f}\langle
i\rvert]=\delta_{ij}$ (9)

where $(-1)^{\left\lvert i\right\rvert}$ ($\left\lvert i\right\rvert=0/1$) is
the fermion parity of $\lvert i\rangle$. Physical wavefunctions are obtained
by contracting all internal legs. Site and bond tensors for fPEPS on honeycomb
lattice are drawn in Fig. 3, where all tensors are set to be _parity even_ in
this paper. More details about fPEPS are represented in Sec. I of Supplemental
Materials (SM)414141See Supplemental Material..

Figure 3: Site tensors $\hat{T}_{u,v}$ and bond tensors $\hat{B}_{x,y,z}$.
Physical spin-1/2 fermions live at sites, while physical Ising spins live at
bonds. Internal legs are represented by triple-lines. Lines on two sides are
internal Ising spins $\left|\tau_{0}\tau_{1}\right)$, whose order follows the
dashed arcs. The middle line is a spinless fermion $c$, where filled/empty
circles label filled/empty states.

Let us work out fPEPS representation of Eq. (8). Imposing translational
symmetry, we focus on tensors in one unit cell, including site tensors
$\hat{T}_{u,v}$ and bond tensors $\hat{B}_{x,y,z}$, as in Fig. 3. Physical
spin-$\frac{1}{2}$ fermions live at sites, while two physical Ising spins live
at two sides of bond centers. Ising spins within a plaquette are enforced to
be the same, and thus are effectively plaquette spins.

An internal leg $(s\alpha)$ is represented as a triple-line, pointing from
site $s$ to bond $\alpha$, where the middle line is a spinless fermion mode
$c_{(s\alpha)}$, while lines at sides are Ising spins. Basis states are
represented as
$(c_{(s\alpha)}^{\dagger})^{n}\left|\tau_{0}\tau_{1}\right)_{(s\alpha)}$,
where vacuum $\left|0\right)_{(s\alpha)}$ is omitted for brevity.
$\tau_{0}\tau_{1}$ are ordered counter-clockwise/clockwise around the site
$u/v$, as indicated by directed dashed arcs in Fig. 3 .

As all spins within a plaquette are the same, an internal spin state of a site
tensor reads
$\left|\tau_{1}\tau_{2}\right)_{(sx)}\left|\tau_{2}\tau_{0}\right)_{(sy)}\left|\tau_{0}\tau_{1}\right)_{(sz)}$,
which is succinctly expressed as $\left|\tau_{0}\tau_{1}\tau_{2}\right)$. Site
tensors for Eq. (8) are

$\displaystyle\hat{T}_{u}=\lvert
0\rangle\otimes_{f}\Big{[}\left|\uparrow\uparrow\uparrow\right)+\left|\downarrow\downarrow\downarrow\right)\Big{]}$
(10) $\displaystyle+f_{u,\uparrow}^{\dagger}\lvert
0\rangle\otimes_{f}\Big{[}c_{(ux)}^{\dagger}\left|\uparrow\uparrow\downarrow\right)+c_{(uz)}^{\dagger}\left|\uparrow\downarrow\uparrow\right)+c_{(uy)}^{\dagger}\left|\downarrow\uparrow\uparrow\right)\Big{]}$
$\displaystyle+f_{u,\downarrow}^{\dagger}\lvert
0\rangle\otimes_{f}\Big{[}-c_{(uy)}^{\dagger}\left|\downarrow\downarrow\uparrow\right)+c_{(ux)}^{\dagger}\left|\downarrow\uparrow\downarrow\right)-c_{(uz)}^{\dagger}\left|\uparrow\downarrow\downarrow\right)\Big{]}$
$\displaystyle\hat{T}_{v}=\lvert
0\rangle\otimes_{f}\Big{[}\left|\uparrow\uparrow\uparrow\right)+\left|\downarrow\downarrow\downarrow\right)\Big{]}$
(11) $\displaystyle+f_{v,\uparrow}^{\dagger}\lvert
0\rangle\otimes_{f}\Big{[}c_{(vy)}^{\dagger}\left|\uparrow\uparrow\downarrow\right)+c_{(vx)}^{\dagger}\left|\uparrow\downarrow\uparrow\right)+c_{(vz)}^{\dagger}\left|\downarrow\uparrow\uparrow\right)\Big{]}$
$\displaystyle+f_{v,\downarrow}^{\dagger}\lvert
0\rangle\otimes_{f}\Big{[}c_{(vx)}^{\dagger}\left|\downarrow\downarrow\uparrow\right)-c_{(vz)}^{\dagger}\left|\downarrow\uparrow\downarrow\right)+c_{(vy)}^{\dagger}\left|\uparrow\downarrow\downarrow\right)\Big{]}$

Similarly, $\langle\tau_{0}\tau_{1}\rvert$ is short for a bond spin state
$\langle\tau_{0}\tau_{1}\rvert_{\alpha}\otimes\left(\tau_{0}\tau_{1}\right|_{(u\alpha)}\left(\tau_{1}\tau_{0}\right|_{(v\alpha)}$.
Bond tensors are expressed as

$\displaystyle\hat{B}_{\alpha}=\langle\uparrow\uparrow\rvert_{\alpha}+\langle\downarrow\downarrow\rvert_{\alpha}+\langle\downarrow\uparrow\rvert_{\alpha}-\langle\uparrow\downarrow\rvert_{\alpha}c_{(u\alpha)}c_{(v\alpha)}$
(12)

As we show in Sec. II of SM4141footnotemark: 41, two $\mathcal{T}$-related
state from contracting fPEPS satisfy Eq. (7), so it indeed gives the fixed-
point wavefunction.

_Tensor equations.-_ Here, we extract symmetry action rules on internal legs
for Eq. (11) and (12), which pave the way for wavefunctions beyond Eq. (8). We
assume that symmetry on physical legs are pushed to gauge transformation on
internal legsPérez-García  _et al._ (2010), see also SM4141footnotemark: 41.

  * •

To impose charge conservation, we require all tensors to be charge neutral,
which can be realized by assigning $c_{(s\alpha)}$ to carry charge
$(-1)^{1-s}$. Note that $f_{s,\sigma}$ carries charge $(-1)^{s}$, and thus

$\displaystyle\left(n_{f;s}+\sum_{\alpha=x,y,z}n_{f;(s\alpha)}\right)\cdot\hat{T}_{s}=\hat{B}_{\alpha}\cdot\left(\sum_{s=u,v}n_{f;(s\alpha)}\right)=0$

  * •

$\mathcal{T}$ action on $(s\alpha)$ are set as

$\displaystyle
W_{(s\alpha)}(\mathcal{T})=\left|\uparrow\uparrow\right)_{(s\alpha)}\left(\downarrow\downarrow\right|_{(s\alpha)}+\left|\downarrow\downarrow\right)_{(s\alpha)}\left(\uparrow\uparrow\right|_{(s\alpha)}$
(13)
$\displaystyle+\operatorname{i}c_{(s\alpha)}^{\dagger}\left|\uparrow\downarrow\right)_{(s\alpha)}\left(\downarrow\uparrow\right|_{(s\alpha)}+\left|\downarrow\uparrow\right)_{(s\alpha)}\left(\uparrow\downarrow\right|_{(s\alpha)}c_{(s\alpha)}$

which gives the following symmetric condition:

$\displaystyle\hat{T}_{s}={}$ $\displaystyle
U_{s}(\mathcal{T})\otimes_{f}W_{(sx)}(\mathcal{T})\otimes_{f}W_{(sy)}(\mathcal{T})\otimes_{f}W_{(sz)}(\mathcal{T})\cdot\hat{T}^{*}_{s}$
$\displaystyle\hat{B}_{\alpha}={}$ $\displaystyle V_{(\alpha
0)}(\mathcal{T})\otimes_{f}V_{(\alpha
1)}(\mathcal{T})\cdot\hat{B}^{*}_{\alpha}\cdot
W^{-1}_{(v\alpha)}(\mathcal{T})\otimes_{f}W^{-1}_{(u\alpha)}(\mathcal{T})$

Here, $U(\mathcal{T})$ and $V(\mathcal{T})$ are $\mathcal{T}$-action on
physical legs defined in Eq. (5). However, as $W(\mathcal{T})$’s are not
parity even, it may not lead to a symmetric wavefunction. In
SM4141footnotemark: 41 Sec. III, we show that this equation contains a hidden
Kasteleyn orientationCimasoni and Reshetikhin (2007); Tarantino and Fidkowski
(2016); Ware  _et al._ (2016); Ellison and Fidkowski (2019), which gives a
$\mathcal{T}$-symmetric condition, .

  * •

Besides the above physical symmetry, such local tensors also host a “gauge
symmetry”:

$\displaystyle\Big{(}n_{\lambda;(s\alpha
0)}+n_{\lambda;(s\widetilde{\alpha}1)}\Big{)}\cdot\hat{T}_{s}$
$\displaystyle=0~{},$
$\displaystyle\hat{B}_{\alpha}\cdot\Big{(}n_{\lambda;(u\alpha
a)}+n_{\lambda;(v\alpha a)}\Big{)}$ $\displaystyle=0~{},$ (14)

where $n_{\lambda;(s\alpha
a)}=(-1)^{s+a}\lvert\downarrow\rangle\langle\downarrow\rvert$ with $a=0/1$
labelling two side lines of $(s\alpha)$, and
$\widetilde{\alpha}=\alpha-(-1)^{s}$. $n_{\lambda}$’s action on all internal
legs in a plaquette $p$ imposes the same spin constraint within $p$, and
generates a $U(1)$ symmetry. We thus get $[U(1)]^{N_{p}}$ symmetry, where
$N_{p}$ is the number of plaquettes. Note that such symmetry acts trivially on
physical legs, and is called “invariant gauge group” (IGG)Wen (7071); Jiang
and Ran (2015, 2017), which is related to topological properties of the
phaseGu  _et al._ (2008); Schuch  _et al._ (2010, 2011) (see also
SM4141footnotemark: 41).

We now extract group relations between $n_{f}$, $W(\mathcal{T})$ and
$n_{\lambda}$, which are coined as _tensor equations_ in this work. Roughly
speaking, IGG gives possible action of the identity element on internal legs,
and then symmetry on internal legs satisfy Eq. (1) up to some IGG elementWen
(7071); Jiang and Ran (2015, 2017). From Eq. (13), the commutator between
$n_{f}$ and $\mathcal{T}$ on internal legs reads

$\displaystyle W_{(s\alpha)}(\mathcal{T})\cdot n_{f;(s\alpha)}\cdot
W_{(s\alpha)}^{-1}(\mathcal{T})=n_{f;(s\alpha)}+n_{D;(s\alpha)}$ (15)

where

$\displaystyle n_{D;(s\alpha)}$
$\displaystyle=(-1)^{s}\left(\lvert\downarrow\uparrow\rangle\langle\downarrow\uparrow\rvert-\lvert\uparrow\downarrow\rangle\langle\uparrow\downarrow\rvert\right)$
$\displaystyle=n_{\lambda;(s\alpha 0)}+n_{\lambda;(s\alpha 1)}$ (16)

Physically, $n_{D}$ gives $U(1)$ gauge theory, but due to the decomposition to
$n_{\lambda}$’s, such gauge theory is killed, leading to short-range entangled
phase4141footnotemark: 41.

For group relation $\mathcal{T}^{2}=F$, a naïve insertion of an IGG element
does not give the correct result. Instead, from Eq. (13) and (16), we have

$\displaystyle\exp\left[\operatorname{i}\frac{\pi}{2}n_{D;(s\alpha)}^{2}\right]\cdot
W_{(s\alpha)}(\mathcal{T})\cdot W_{(s\alpha)}^{*}(\mathcal{T})=F_{(s\alpha)}$
(17)

In Sec. IV of SM4141footnotemark: 41, we show that Eq. (17) is indeed
consistent with $\mathcal{T}^{2}=F$.

The commutator between $W(\mathcal{T})$ and $n_{\lambda}$ completes tensor
equations:

$\displaystyle W_{(s\alpha)}(\mathcal{T})\cdot n^{*}_{\lambda;(s\alpha
a)}\cdot W_{(s\alpha)}^{-1}(\mathcal{T})=-n_{\lambda;(s\alpha a)}+(-)^{s+a}$
(18)

_Edge theories from tensor equations.-_ In the following, we show that edge
properties of the QSH phase can be extracted from tensor equations from Eq.
(15) to Eq. (18). The anomalous edge theory is characterized by fusion of two
$\mathcal{T}$-fluxElse and Nayak (2014). To see this, we turn to the edge
theory in Eq. (3). By rotating $\theta(x)$ by $2\pi$ angle counter-clockwise
within region $[x_{0},x_{1}]$, we get the current density

$\displaystyle\int\operatorname{d}\\!t\,j(x)=\int\operatorname{d}\\!t\,\frac{\partial_{t}\theta}{2\pi}=\begin{cases}0&x\leq
x_{0}\text{ or }x\geq x_{1}\\\ 1&x_{0}<x<x_{1}\end{cases}$ (19)

So, a unit charge is pumped from $x_{0}$ to $x_{1}$Qi  _et al._ (2008b); Fu
and Kane (2006). According to Eq. (2), rotating $\theta$ by $\pi$ on
$[x_{0},x_{1}]$ is equivalent to acting $\mathcal{T}$ on this segment, which
creates $\mathcal{T}$-flux at two endsChen and Vishwanath (2015). The unit
charge pumping due to $2\pi$ rotation of $\theta$ is interpreted as _two
$\mathcal{T}$-flux fuses to an electron/hole_.

We now extract such fusion rule from tensor equations. As in Fig. 4, to obtain
edge theory of system $A$, we cut tensors within $A$ from the infinite fPEPS,
contract all internal legs within $A$, and obtain a large tensor
$\hat{T}_{A}$. $\hat{T}_{A}$ has $L$ boundary legs labeled by index
$j\in\partial A=\\{1,2,\dots,L\\}$, forming Hilbert space
$\mathbb{H}_{\partial A}$. As shown in Sec. V of SM4141footnotemark: 41, the
edge Hilbert space $\mathbb{H}_{edge}$ are formed by states in
$\mathbb{H}_{\partial A}$ that are invariant under IGG action. Let $P_{edge}$
be the projector from $\mathbb{H}_{\partial A}$ to $\mathbb{H}_{edge}$. In our
case, $P_{edge}$ identifies Ising spins belonging to the same plaquette:
$\tau_{j+\frac{1}{4}}=\tau_{j+\frac{3}{4}}$, where $\tau_{j\pm\frac{1}{4}}$
are spins at boundary leg $j$.

Figure 4: Left: Tensor $\hat{T}_{A}$ cutting from the infinite fPEPS, whose
boundary legs are numbered from $1$ to $L$. Right: Details of boundary legs of
$\hat{T}_{A}$. Leg $j$ is a triple-line, representing two Ising spins (thick
blue line) $\tau_{j-\frac{1}{4}}$ and $\tau_{j+\frac{1}{4}}$, and one spinless
fermion $c_{j}$ (thin orange line). Plaquette between $j$ and $j+1$ are
labeled as $p=j+\frac{1}{2}$.

By projecting $W_{\partial A}(\mathcal{T})\equiv\bigotimes_{f;j\in\partial
A}W_{j}(\mathcal{T})$ to $\mathbb{H}_{edge}$, one gets $\mathcal{T}$ action on
edge:

$\displaystyle U_{edge}(\mathcal{T})\mathcal{K}=W_{\partial
A}(\mathcal{T})\mathcal{K}\cdot P_{edge}$ (20)

Let $M=\\{2,3,\dots,l\\}$ be a subregion of $\partial A$. $\mathcal{T}$-flux
at ends of $M$ are created by a charge-neutral string operator
$U_{M}(\mathcal{T})\mathcal{K}$, where

$\displaystyle U_{M}(\mathcal{T})\mathcal{K}=P_{edge}\cdot w_{l+1}\cdot
w_{1}\cdot W_{M}(\mathcal{T})\mathcal{K}\cdot P_{edge}$ (21)

Here, $w_{1/(l+1)}$ are local operators on leg $1/(l+1)$. The charge-neutral
requirement for $U_{M}(\mathcal{T})\mathcal{K}$ puts the following constraint
on $w_{1/(l+1)}$4141footnotemark: 41

$[w_{1},n_{{}_{f;1}}]=n^{(0)}_{\lambda,\frac{3}{2}}\cdot
w_{1}\,;~{}[w_{l+1},n_{{}_{f;l+1}}]=n^{(1)}_{\lambda;l+\frac{1}{2}}\cdot
w_{l+1}\,,$ (22)

where $n_{\lambda;p}^{(0/1)}$ are IGG elements acting on
$\left|\tau_{p\mp\frac{1}{4}}\right)$.

Let $j=1$ and $l+1$ be $v$ site, and then we can solve

$\displaystyle w_{1}={}$
$\displaystyle\sum_{\tau_{\frac{3}{4}}}c_{1}\left|\tau_{\frac{3}{4}},\downarrow_{\frac{5}{4}}\right)\left(\tau_{\frac{3}{4}},\uparrow_{\frac{5}{4}}\right|+\left|\tau_{\frac{3}{4}},\uparrow_{\frac{5}{4}}\right)\left(\tau_{\frac{3}{4}},\downarrow_{\frac{5}{4}}\right|\,;$
$\displaystyle w_{l+1}={}$
$\displaystyle\sum_{\tau_{l+\frac{5}{4}}}c^{\dagger}_{l+1}\left|\downarrow_{l+\frac{3}{4}},\tau_{l+\frac{5}{4}}\right)\left(\uparrow_{l+\frac{3}{4}},\tau_{l+\frac{5}{4}}\right|$
$\displaystyle+\left|\uparrow_{l+\frac{3}{4}},\tau_{l+\frac{5}{4}}\right)\left(\downarrow_{l+\frac{3}{4}},\tau_{l+\frac{5}{4}}\right|\,.$
(23)

It is then straightforward to verify

$\displaystyle[U_{M}(\mathcal{T})\mathcal{K}]^{2}=P_{edge}\cdot\Omega_{1}\cdot\Omega_{l+1}\cdot\prod_{j=2}^{l}F_{j}\cdot
P_{edge}$ (24)

where $\Omega_{1}=c_{1}\exp\left[\operatorname{i}\frac{\pi}{2}\cdot
n^{(1)}_{\lambda;\frac{3}{2}}\right]$ is a hole, and
$\Omega_{l+1}=c_{l+1}^{\dagger}\exp\left[\operatorname{i}\frac{\pi}{2}\cdot
n^{(0)}_{\lambda;l+\frac{1}{2}}\right]$ an electron (see details in SM Sec.
VI4141footnotemark: 41). They can be viewed as quasi-particles from fusing two
$\mathcal{T}$-flux. We thus recover the anomalous edge theory, which suggests
that any fPEPS that satisfies tensor equations belongs to the QSH phase222 It
is possible that additional IGG elements emerge in the thermodynamic
limitDreyer  _et al._ (2020). In such case, fPEPS wavefunction may give other
phases, such as spontaneously symmetry breaking phases. .

_Variational tensor wavefunctions.-_ The fixed-point wavefunction in Eq. (8)
is quite artificial, as one needs additional plaquette Ising spins. In the
following, let us try to construct variational wavefunctions for
spin-$\frac{1}{2}$ fermions on a honeycomb lattice by solving tensor
equations. We will further demonstrate these wavefunctions give the desired
many-body topological invariants in the next part.

We start from fPEPS with two types of site tensors. Each site tensor has one
physical and three internal legs, which can be expressed as
$\hat{T}_{s}=(T_{s})_{ijk,p}\left|i\right)_{(s1)}\left|j\right)_{(s2)}\left|k\right)_{(s3)}\lvert
p\rangle_{s}$, with $s=u/v$ labeling the sublattices, and subindices $1$ to
$3$ ordering internal legs clockwise. Physical spin-$\frac{1}{2}$ fermions
$f_{s\sigma}$’s carry opposite charges on site $u$ and $v$. For simplicity, we
choose basis states of an internal leg $(s\alpha)$ to be
$\left\\{\left|\uparrow\uparrow\right)_{(s\alpha)},~{}c_{(s\alpha)}^{\dagger}\left|\uparrow\downarrow\right)_{(s\alpha)},~{}\left|\downarrow\uparrow\right)_{(s\alpha)},~{}\left|\downarrow\downarrow\right)_{(s\alpha)}\right\\}$.

Action of $\mathcal{T}$ on physical legs follows Eq. (5). Symmetries impose
constraints on site tensor $\hat{T}_{s}$ as

$\displaystyle\bigotimes_{\alpha=1}^{3}\mathop{{}_{\mkern-4.0muf}}W_{(s\alpha)}(\mathcal{T})\otimes_{f}U_{s}(\mathcal{T})\cdot\hat{T}^{*}_{s}$
$\displaystyle=\hat{T}_{s}$
$\displaystyle\Big{[}\sum_{\alpha=1}^{3}n_{f;(s\alpha)}+n_{f;s}\Big{]}\cdot\hat{T}_{s}$
$\displaystyle=0\,.$ (25)

In addition, site tensors should also be invariant under plaquette IGG
$n_{\lambda}$’s, as in Eq. (14). $W(\mathcal{T})$, $n_{f}$ and $n_{\lambda}$
satisfy tensor equations from Eq. (15) to (18), and can simply take the same
form as in the honeycomb example. By solving these tensor constraints, we
obtain $7$ linearly independent solutions for $\hat{T}_{u/v}$, as listed in SM
Sec. VIII4141footnotemark: 41. Bond tensors are set to be Eq. (12), which
satisfy all tensor equations4141footnotemark: 41.We also calculated the
variational ansatz for system in square lattice, see SM4141footnotemark: 41.

_Many-body topological invariants.-_ One can also diagnose the interacting QSH
wavefunction by calculating the many-body topological invariant proposed in
Ref. Shiozaki  _et al._ , 2018. The system is put on a cylinder periodic in
$y-$direction and open in $x-$direction, which is divided as follows The
topological invariant is given by

$\displaystyle Z$ $\displaystyle=\mathrm{Tr}\left[\rho^{+}_{R_{1}\cup
R_{3}}C^{R_{1}}_{T}(\rho^{-}_{R_{1}\cup
R_{3}})^{\mathsf{T}_{1}}[C^{R_{1}}_{T}]^{\dagger}\right]$ (26)

with

$\displaystyle\rho^{\pm}_{R_{1}\cup R_{3}}$
$\displaystyle=\mathrm{Tr}_{\overline{R_{1}\cup R_{3}}}\left[\exp{\frac{\pm
2\pi\mathrm{i}y\sum_{\mathrm{r}\in
R_{2}}n(\mathrm{r})}{L_{y}}}|\Phi\rangle\langle\Phi|\right]$

Here, $\mathsf{T}_{1}$ is the fermionic partial transpose of region $R_{1}$,
and $C^{R_{1}}_{T}$ is defined by $C^{R_{1}}_{T}c_{j\in
R_{1}}(C^{R_{1}}_{T})^{{\dagger}}=c^{{\dagger}}_{k\in
R_{1}}U(\mathcal{T})_{kj}$. The phase of $Z$ equals $0$/$\pi$ when
$|\Phi\rangle$ is in trivial/topological phase. Such quantity can be
calculated using numerics, where we present some in Fig. 5.

Figure 5: Calculation of $\mathrm{sgn}(Z)$ with respect to $c_{u}$, where
$c_{u}$ labels some tensor entry, see Supplemental Material for more details.
$\mathrm{sgn}(Z)=-1$ is a signature of topologically nontrivial phase, see
Ref. Shiozaki  _et al._ , 2018. Sizes of $R_{1,2,3}$ are set to be equal, with
length in $x$-direction to be $L_{R}$, and in $y$-direction to be $2$.

Note that the appearance of $\mathrm{sgn}(Z)=1$ may be due to various reasons,
e.g. the finite size effect, see details in SM4141footnotemark: 41.

_Discussion.-_ In this work, from the fPEPS representation of the fixed-point
wavefunction of the QSH phase, we extract tensor equations from Eq. (15) to
(18) characterizing symmetry rules on local tensors. By solving these
equations, we can obtain general forms for symmetry actions on internal legs.
Variational ansatz for the QSH phase are solved by imposing such symmetry
constraints on local tensors.

This work leaves several interesting future directions. Firstly, to express
variational ansatz for topological phases in half-filled spin-1/2 _electronic_
models, it is necessary to generalize our framework to include tensors with
odd parity. Additionally, developing variational numerical algorithms for
symmetric fPEPS wavefunctions obtained in this study would be desirable to
simulate the QSH phase in strongly correlated models. On the analytical side,
we aim to explore other fermionic topological phases, such as topological
superconductors and topologically ordered phases. Of particular interest is
the investigation of chiral phases, such as the $p+\operatorname{i}p$
topological superconductor Read and Green (2000); Ivanov (2001); Stone and Roy
(2004). The question of whether fPEPS can accurately represent these chiral
phases with a finite bulk gap remains an intriguing puzzle Wahl  _et al._
(2013, 2014); Dubail and Read (2015). Furthermore, tensor networks readily
incorporate spatial symmetries Jiang and Ran (2015), enabling the construction
of variational tensor wavefunctions for gapped electronic liquid phases and
high-order topological insulators/superconductorsRasmussen and Lu (2020); Song
_et al._ (2020); Else and Thorngren (2019). Exploring these possibilities
holds significant potential for advancing our understanding of exotic
topological phases.

We would like to thank Qing-Rui Wang and Xie Chen for helpful discussions. The
work is supported by MOST NO. 2022YFA1403901, NSFC NO. 12104451, and funds
from Strategic Priority Research Program of CAS (No. XDB28000000).

## References

  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi  _et al._ (2008a) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008a).
  * Fidkowski and Kitaev (2010) L. Fidkowski and A. Kitaev, Phys. Rev. B 81, 134509 (2010).
  * Fidkowski and Kitaev (2011) L. Fidkowski and A. Kitaev, Phys. Rev. B 83, 075103 (2011).
  * Chen  _et al._ (2011) X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B 84, 235128 (2011).
  * Gu and Wen (2014) Z.-C. Gu and X.-G. Wen, Phys. Rev. B 90, 115141 (2014).
  * Gu and Levin (2014) Z.-C. Gu and M. Levin, Phys. Rev. B 89, 201113 (2014).
  * Wang and Gu (2018) Q.-R. Wang and Z.-C. Gu, Phys. Rev. X 8, 011055 (2018).
  * Wang and Gu (2020) Q.-R. Wang and Z.-C. Gu, Phys. Rev. X 10, 031055 (2020).
  * Wang  _et al._ (2017) C. Wang, C.-H. Lin, and Z.-C. Gu, Phys. Rev. B 95, 195147 (2017).
  * Cheng  _et al._ (2018) M. Cheng, N. Tantivasadakarn, and C. Wang, Phys. Rev. X 8, 011054 (2018).
  * Wang  _et al._ (2014) C. Wang, A. C. Potter, and T. Senthil, Science 343, 629 (2014).
  * Wang and Senthil (2014) C. Wang and T. Senthil, Phys. Rev. B 89, 195124 (2014).
  * Senthil (2015) T. Senthil, Annu. Rev. Condens. Matter Phys. 6, 299 (2015).
  * Witten (2016) E. Witten, Rev. Mod. Phys. 88, 035001 (2016).
  * Freed and Hopkins (2021) D. S. Freed and M. J. Hopkins, Geometry & Topology 25, 1165 (2021).
  * Ware  _et al._ (2016) B. Ware, J. H. Son, M. Cheng, R. V. Mishmash, J. Alicea, and B. Bauer, Phys. Rev. B 94, 115127 (2016).
  * Tarantino and Fidkowski (2016) N. Tarantino and L. Fidkowski, Phys. Rev. B 94, 115115 (2016).
  * Wang  _et al._ (2018) Z. Wang, S.-Q. Ning, and X. Chen, Phys. Rev. B 98, 094502 (2018).
  * Metlitski (2019) M. A. Metlitski, “A 1d lattice model for the boundary of the quantum spin-hall insulator,” (2019), arXiv:1908.08958 [cond-mat.str-el] .
  * Son and Alicea (2019) J. H. Son and J. Alicea, Phys. Rev. B 100, 155107 (2019).
  * Wang  _et al._ (2021) Q.-R. Wang, Y. Qi, C. Fang, M. Cheng, and Z.-C. Gu, “Exactly solvable lattice models for interacting electronic insulators in two dimensions,” (2021), arXiv:2112.15533 [cond-mat.str-el] .
  * Barthel  _et al._ (2009) T. Barthel, C. Pineda, and J. Eisert, Phys. Rev. A 80, 042333 (2009).
  * Corboz  _et al._ (2010) P. Corboz, G. Evenbly, F. Verstraete, and G. Vidal, Phys. Rev. A 81, 010303 (2010).
  * Kraus  _et al._ (2010) C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A 81, 052338 (2010).
  * Gu  _et al._ (2010) Z.-C. Gu, F. Verstraete, and X.-G. Wen, “Grassmann tensor network states and its renormalization for strongly correlated fermionic and bosonic states,” (2010), arXiv:1004.2563 [cond-mat.str-el] .
  * Schuch  _et al._ (2011) N. Schuch, D. Pérez-García, and I. Cirac, Phys. Rev. B 84, 165139 (2011).
  * Bultinck  _et al._ (2017a) N. Bultinck, D. J. Williamson, J. Haegeman, and F. Verstraete, Journal of Physics A: Mathematical and Theoretical 51, 025202 (2017a).
  * Bultinck  _et al._ (2017b) N. Bultinck, D. J. Williamson, J. Haegeman, and F. Verstraete, Phys. Rev. B 95, 075108 (2017b).
  * Wille  _et al._ (2017) C. Wille, O. Buerschaper, and J. Eisert, Phys. Rev. B 95, 245127 (2017).
  * Cirac  _et al._ (2021) J. I. Cirac, D. Perez-Garcia, N. Schuch, and F. Verstraete, Reviews of Modern Physics 93, 045003 (2021).
  * Shapourian  _et al._ (2017) H. Shapourian, K. Shiozaki, and S. Ryu, Phys. Rev. Lett. 118, 216402 (2017).
  * Shiozaki  _et al._ (2018) K. Shiozaki, H. Shapourian, K. Gomi, and S. Ryu, Phys. Rev. B 98, 035151 (2018).
  * Haldane (1981) F. Haldane, Journal of Physics C: Solid State Physics 14, 2585 (1981).
  * Wu  _et al._ (2006) C. Wu, B. A. Bernevig, and S.-C. Zhang, Phys. Rev. Lett. 96, 106401 (2006).
  * Xu and Moore (2006) C. Xu and J. E. Moore, Phys. Rev. B 73, 045322 (2006).
  * Goldstone and Wilczek (1981) J. Goldstone and F. Wilczek, Phys. Rev. Lett. 47, 986 (1981).
  * Chen  _et al._ (2014) X. Chen, Y.-M. Lu, and A. Vishwanath, Nature communications 5, 3507 (2014).
  * Liu  _et al._ (2014) Z.-X. Liu, Z.-C. Gu, and X.-G. Wen, Phys. Rev. Lett. 113, 267206 (2014).
  * Note (1) We mention that to fully determine $\lvert\Psi\rangle$, we should explicitly write down the fermion order and set the $\pm 1$ phase for each $\Psi(c)$ However, we will not present such information as it is quite complicated and is unnecessary for the following discussion.
  * Note (41) See Supplemental Material.
  * Pérez-García  _et al._ (2010) D. Pérez-García, M. Sanz, C. Gonzalez-Guillen, M. M. Wolf, and J. I. Cirac, New Journal of Physics 12, 025010 (2010).
  * Cimasoni and Reshetikhin (2007) D. Cimasoni and N. Reshetikhin, Communications in Mathematical Physics 275, 187 (2007).
  * Ellison and Fidkowski (2019) T. D. Ellison and L. Fidkowski, Phys. Rev. X 9, 011016 (2019).
  * Wen (7071) X.-G. Wen, Phys. Rev. B 65, 165113 (2002, cond-mat/0107071).
  * Jiang and Ran (2015) S. Jiang and Y. Ran, Phys. Rev. B 92, 104414 (2015).
  * Jiang and Ran (2017) S. Jiang and Y. Ran, Phys. Rev. B 95, 125107 (2017).
  * Gu  _et al._ (2008) Z.-C. Gu, M. Levin, and X.-G. Wen, Phys. Rev. B 78, 205116 (2008).
  * Schuch  _et al._ (2010) N. Schuch, I. Cirac, and D. Pérez-García, Annals of Physics 325, 2153 (2010).
  * Else and Nayak (2014) D. V. Else and C. Nayak, Phys. Rev. B 90, 235137 (2014).
  * Qi  _et al._ (2008b) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Nature Physics 4, 273 (2008b).
  * Fu and Kane (2006) L. Fu and C. L. Kane, Phys. Rev. B 74, 195312 (2006).
  * Chen and Vishwanath (2015) X. Chen and A. Vishwanath, Phys. Rev. X 5, 041034 (2015).
  * Note (2) It is possible that additional IGG elements emerge in the thermodynamic limitDreyer  _et al._ (2020). In such case, fPEPS wavefunction may give other phases, such as spontaneously symmetry breaking phases.
  * Read and Green (2000) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * Ivanov (2001) D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
  * Stone and Roy (2004) M. Stone and R. Roy, Phys. Rev. B 69, 184511 (2004).
  * Wahl  _et al._ (2013) T. B. Wahl, H.-H. Tu, N. Schuch, and J. I. Cirac, Phys. Rev. letters 111, 236805 (2013).
  * Wahl  _et al._ (2014) T. B. Wahl, S. T. Haßler, H.-H. Tu, J. I. Cirac, and N. Schuch, Phys. Rev. B 90, 115133 (2014).
  * Dubail and Read (2015) J. Dubail and N. Read, Phys. Rev. B 92, 205307 (2015).
  * Rasmussen and Lu (2020) A. Rasmussen and Y.-M. Lu, Phys. Rev. B 101, 085137 (2020).
  * Song  _et al._ (2020) Z. Song, C. Fang, and Y. Qi, Nature communications 11, 1 (2020).
  * Else and Thorngren (2019) D. V. Else and R. Thorngren, Phys. Rev. B 99, 115116 (2019).
  * Dreyer  _et al._ (2020) H. Dreyer, L. Vanderstraeten, J.-Y. Chen, R. Verresen, and N. Schuch, “Robustness of critical $u(1)$ spin liquids and emergent symmetries in tensor networks,” (2020), arXiv:2008.04833 [cond-mat.str-el] .

Supplementary Materials: Variational Tensor Wavefunctions for the Interacting
Quantum Spin Hall Phase

In this supplemental material, we provide a brief review of fermionic tensor
network states (Sec. I), detailed information for the fixed-point tensor
network wavefunction of QSH phase (Sec. II), Kasteleyn orientation (Sec. III),
detailed derivation of $\mathcal{T}^{2}=F$ on internal legs (Sec. IV), edge
theory from infinite PEPS (Sec. V), derivation of fusion of $\mathcal{T}$-flux
(Sec. VI), the variational ansatz for the QSH phase on a square system
(Sec.VII), and numerical calculation of many-body topological
invariant(Sec.VIII).

##  I Symmetric fermionic tensor network states

In this section, we review the fundamentals of fermionic tensor network
statesKraus  _et al._ (2010); Bultinck  _et al._ (2017b, a), and fix our
notations used in the main text.

###  I.1 Fermionic tensors and tensor contraction

Building blocks of fermionic tensor network states are fermionic tensors,
which live in fermionic tensor product (labeled as $\otimes_{f}$) of legs.
Legs with inward/outward arrows are local fermion Hilbert spaces for ket/bra
states, where the parity of state $\lvert i\rangle/\langle i\rvert$ is
$(-1)^{\left\lvert i\right\rvert}$ with $\left\lvert
i\right\rvert\in\\{0,1\\}$. Exchanging states of two legs gives $-1$ if these
two states are both parity odd:

$\displaystyle\lvert i\rangle_{a}\otimes_{f}\lvert
j\rangle_{b}=(-1)^{\left\lvert i\right\rvert\left\lvert j\right\rvert}\lvert
j\rangle_{b}\otimes_{f}\lvert i\rangle_{a}$ (S1)

As an example of fermionic tensors, let us consider tensor $\hat{T}$ with
three legs, say $a,b,c$:

$\displaystyle\hat{T}=(T_{abc})_{ijk}\lvert i\rangle_{a}\otimes_{f}\lvert
j\rangle_{b}\otimes_{f}\lvert k\rangle_{c}$ (S2)

Leg indices $abc$ are sometimes ignored when there is no confusion. We may
also omit $\otimes_{f}$’s and use a more compact form $\hat{T}=T_{ijk}\lvert
ijk\rangle$.

Ket and bra states can be contracted using $\mathrm{fTr}$, defined as

$\displaystyle\mathrm{fTr}[\langle i\rvert\otimes_{f}\lvert
j\rangle]=(-)^{\left\lvert i\right\rvert\left\lvert
j\right\rvert}\mathrm{fTr}[\lvert j\rangle\otimes_{f}\langle
i\rvert]=\delta_{ij}$ (S3)

It is noteworthy that the order of contracted states matters as extra $-1$ may
be produced. Generalization to tensor contractions is straightforward. As
shown in Fig. 1, for two fermionic tensors $\hat{M}=M^{ijk}\langle
i\rvert_{a}\langle j\rvert_{b}\langle k\rvert_{c}$, $\hat{N}=N_{lm}\lvert
l\rangle_{b}\lvert m\rangle_{d}$,

$\displaystyle\mathrm{fTr}_{b}[\hat{M}\otimes_{f}\hat{N}]\equiv(-1)^{\left\lvert
j\right\rvert\left\lvert k\right\rvert}M^{ijk}N_{lm}\delta_{jl}\langle
i\rvert\langle k\rvert\lvert m\rangle~{},$ (S4)

which is represented graphically by connecting inward and outward leg. We may
omit $\mathrm{fTr}$’s and $\otimes_{f}$’s and use $\hat{M}\cdot\hat{N}$ to
represent tensor contraction.

Figure 1: Graphic representation of contraction between fermionic tensors
$\hat{M}$ and $\hat{N}$. The outward legs are bra spaces and inward legs are
ket spaces. The intersection of leg $b$ and $c$ indicates a possible extra
minus sign due to exchanging fermions.

For local tensors in tensor network states, there are two types of legs:
internal ones and physical ones. States in physical legs is denoted by
$\lvert\bullet\rangle/\langle\bullet\rvert$, while states in internal legs by
$\left|\bullet\right)/\left(\bullet\right|$. To get a physical wavefunction,
all internal legs are contracted. By fixing parity of all local tensors,
physical wavefunctions also have fixed parity. In this work, we focus on the
case where all local tensors are _parity even_.

###  I.2 Gauge transformation and symmetries of fPEPS

We consider a particular type of fermionic tensor network – fermionic
projected entangled-pair state(s) (fPEPS)Schuch  _et al._ (2011). As shown in
Fig. 2(a), we focus on fPEPS with both site tensors $\hat{T}_{s}$ and bond
tensors $\hat{B}_{ss^{\prime}}$, where $s,s^{\prime}$ label neighbouring site
coordinates. Without loss of generality, we assume that internal legs of site
tensors all point inward (ket spaces), while those of bond tensors point
outward (bra spaces). Physical wavefunction then reads
$\lvert\Psi\rangle=\mathrm{fTr}\left[\mathbb{B}\otimes_{f}\mathbb{T}\right]$,
where $\mathbb{T}=\bigotimes_{s}\mathop{{}_{\mkern-4.0muf}}\hat{T}_{s}$ and
$\mathbb{B}=\bigotimes_{\langle
ss^{\prime}\rangle}\mathop{{}_{\mkern-4.0muf}}\hat{B}_{ss^{\prime}}$ Note that
as all tensors are parity even, different orders of tensor contraction give
the same state.

Figure 2:  (a) A bond tensor $\hat{B}_{ss^{\prime}}$ and its neighbouring site
tensors. (b) A $3\times 3$ fPEPS on square lattice with boundary legs. (c)
$\mathrm{IGG}$ invariance of a single site tensor.

Representation of one physical wavefunction $\lvert\Psi\rangle$ by fPEPS is
far from unique. In particular, different fPEPS represent the same
wavefunction if they are related by some gauge transformation, as

$\displaystyle\lvert\psi\rangle=\mathrm{fTr}\left[\mathbb{B}\otimes_{f}\mathbb{T}\right]=\mathrm{fTr}\left[\mathbb{B}\otimes_{f}\mathbb{W}^{-1}\otimes_{f}\mathbb{W}\otimes_{f}\mathbb{T}\right]$
(S5)

Here, $\mathbb{W}$ and $\mathbb{W}^{-1}$ are tensor products of gauge
transformation $W$’s on internal legs:

$\displaystyle\mathbb{W}$
$\displaystyle=W_{(s_{1}\alpha_{1})}\otimes_{f}W_{(s_{1}\alpha_{2})}\otimes_{f}\cdots\otimes_{f}W_{(s_{1}\alpha_{m})}\otimes_{f}W_{(s_{2}\alpha_{1})}\otimes_{f}\cdots\otimes_{f}W_{(s_{n}\alpha_{m})}~{},$
$\displaystyle\mathbb{W}^{-1}$
$\displaystyle=W^{-1}_{(s_{n}\alpha_{m})}\otimes_{f}\cdots\otimes_{f}W^{-1}_{(s_{n}\alpha_{1})}\otimes_{f}\cdots
W^{-1}_{(s_{n-1}\alpha_{m})}\otimes_{f}\cdots\otimes_{f}W^{-1}_{(s_{1}\alpha_{1})}$
(S6)

where $(s\alpha)$ labels internal leg, and

$\displaystyle
W_{(s\alpha)}\left|i\right)_{(s\alpha)}=\sum_{b}[W_{(s\alpha)}]_{ji}\left|j\right)_{(s\alpha)}~{},\quad\left(i\right|_{(s\alpha)}W^{-1}_{(s\alpha)}=\sum_{j}\left(j\right|_{(s\alpha)}\left[W_{(s\alpha)}^{-1}\right]_{ij}$
(S7)

$W$’s in general do not have fixed parity, and thus permuting $W$’s and
$W^{-1}$’s may lead to fermion swapping gate.

For the case where $\lvert\Psi\rangle$ is invariant under symmetry $g$, we
assume that $g$-action on physical legs can be pushed to gauge transformation
on internal legs of local tensors:

$\displaystyle
U_{s}(g)\otimes_{f}\left(\bigotimes_{\alpha}\mathop{{}_{\mkern-4.0muf}}W_{(s\alpha)}(g)\right)\cdot\hat{T}_{s}=\hat{T}_{s}~{},\quad\hat{B}_{ss^{\prime}}\cdot
W^{-1}_{(s\alpha)}(g)\otimes_{f}W^{-1}_{(s^{\prime}\alpha^{\prime})}(g)=\hat{B}_{ss^{\prime}}$
(S8)

The above equations give symmetry constraints for local tensors. We remark
that to get a symmetric wavefunction, orders of $W(g)$’s in the above equation
are essential. In particular, as we will show in Sec. III, a valid order of
$W(g)$’s gives a Kasteleyn orientation on the lattice.

As shown in Fig. 2(c), there exists a special kind of gauge transformation
$\mathbb{D}$, which leaves every single tensor invariant:

$\displaystyle\left(\bigotimes_{\alpha}\mathop{{}_{\mkern-4.0muf}}D_{(s\alpha)}\right)\cdot\hat{T}_{s}=\hat{T}_{s}~{},\quad\hat{B}_{ss^{\prime}}\cdot
D^{-1}_{(s\alpha)}\otimes_{f}D^{-1}_{(s^{\prime}\gamma)}=\hat{B}_{ss^{\prime}}~{}.$
(S9)

Such gauge transformation form invariant gauge group ($\mathrm{IGG}$). In this
work, we focus on the case where $D$’s are parity even.

Note that the group always have a trivial center $H$ formed by phase factors
$\chi_{(s\alpha)}$ that satisfy $\prod_{(s\alpha)}\chi_{(s\alpha)}=1$. In
addition, if the $\mathrm{IGG}$ is a $U(1)$ group generated by
$n_{D;(s\alpha)}$, namely, $D_{(s\alpha)}(\theta)=\exp[\operatorname{i}\theta
n_{D;(s\alpha)}]$. We then have
$\left(\sum_{\alpha}n_{D;(s\alpha)}\right)\cdot\hat{T}^{s}=0$

###  I.3 Plaquette $\mathrm{IGG}$ and vanishing long-range entanglement

In this work, all internal legs can be further decomposed to tensor product of
local Hilbert space:

$\displaystyle\mathcal{H}_{(s\alpha)}=\bigotimes_{a}\mathop{{}_{\mkern-4.0muf}}\mathcal{H}_{(s\alpha
a)}$ (S10)

Graphically, an internal leg are represented by multiple lines, and we use
index $a$ to label these lines.

As shown in Fig. 3, we assume that all elements of $\mathrm{IGG}$ are parity
even and have a “plaquette decomposition”Jiang and Ran (2017)

$\displaystyle
D_{(s\alpha)}=\bigotimes_{a}\mathop{{}_{\mkern-4.0muf}}D_{(s\alpha a)}$ (S11)

Here, $a=0/1$ denote lines at two sides.

Figure 3: Any $\mathrm{IGG}$ element can be decomposed to plaquette IGG.

In particular, for site tensors on honeycomb lattice, the plaquette IGG action
gives

(S12)

Nontrivial $\mathrm{IGG}$ element often leads to topological ground state
degeneracy. To see this, let us consider a fPEPS with periodic boundary
condition on square lattice:

$\displaystyle\lvert\Psi\rangle=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{0.5}[0.5]{{\leavevmode\hbox{\set@color{}}}}}}}$}~{},$

where the bond tensors, physical legs and leg orientations are omitted for
brevity. If such tensor network has a nontrivial $\mathrm{IGG}$ element, we
can insert “$\mathrm{IGG}$ loops” in internal legs, which leaves the
wavefunction invariant:

$\displaystyle\lvert\Psi\rangle=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{0.5}[0.5]{{\leavevmode\hbox{\set<EMAIL_ADDRESS>
However, by inserting non-contractible loops of $\mathrm{IGG}$ action, we
obtain a different state as

<EMAIL_ADDRESS>(S13)

As one can move the non-contractible loop without energy cost,
$\lvert\Psi_{v}\rangle$ is indistinguishable from $\lvert\Psi\rangle$ by local
observables. Thus, if $\lvert\psi\rangle$ is a ground state of a local
Hamiltonian, a nontrivial $\mathrm{IGG}$ indicates that the Hamiltonian has
topological ground state degeneracy. In other words, $\lvert\psi\rangle$ is a
long-range entangled state.

However, the plaquette decomposition of the $\mathrm{IGG}$ element kills the
long-range entanglement. To see this, we apply Eq. (S12) for
$\lvert\Psi_{v}\rangle$ in Eq. (S13), and obtain

$\displaystyle\lvert\Psi_{v}\rangle=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{0.5}[0.5]{{\leavevmode\hbox{\set@color{}}}}}}}$}~{}=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{0.5}[0.5]{{\leavevmode\hbox{\set@color{}}}}}}}$}=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{0.5}[0.5]{{\leavevmode\hbox{\set@color{}}}}}}}$}~{}=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{0.5}[0.5]{{\leavevmode\hbox{\set@color{}}}}}}}$}~{}=\lvert\Psi\rangle$

The non-contractible loop of an $\mathrm{IGG}$ element acts trivially on the
tensor network state if the $\mathrm{IGG}$ element has plaquette
decomposition. Thus, topological degeneracy is removed and we obtain a short-
range entangled state.

##  II Tensor representation of the fixed-point wavefunction for the QSH phase

In this section, we will prove that the wavefunction constructed from tensors
in Eq. (10), (11), and (12) in the main text is $\mathcal{T}$ symmetric.

The fPEPS physical wavefunction is
$\lvert\Psi\rangle=\mathrm{fTr}\left[\mathbb{T}\otimes_{f}\mathbb{B}\right]$.
Here,

$\displaystyle\mathbb{T}=\bigotimes_{\bm{\mathbf{r}}s}\mathop{{}_{\mkern-4.0muf}}\hat{T}_{\bm{\mathbf{r}}s}~{},\quad\mathbb{B}=\bigotimes_{\bm{\mathbf{r}}\alpha}\mathop{{}_{\mkern-4.0muf}}\hat{B}_{\bm{\mathbf{r}}\alpha}~{},$
(S14)

where $\bm{\mathbf{r}}$ is coordinate for unit cell, $s=u/v$, and
$\alpha=x/y/z$. The order of tensors is not important since site and bond
tensors are all parity even.

Such wavefunction can be organized according to plaquette Ising spin
configurations $c$:

$\displaystyle\lvert\Psi\rangle=\sum_{c}\Psi(c)\lvert
c\rangle\otimes\lvert\psi_{c}\rangle$ (S15)

Here, $c$ represents a Ising spin configuration, and $\lvert\psi_{c}\rangle$
is the spin-1/2 fermion decoration for $c$. We choose some fixed order for
physical fermions, and thus there is no ambiguous sign for $\Psi_{c}$. Entries
of site and bond tensors in the main take $\pm 1$, so $\Psi(c)$ also takes
$\pm 1$. In the following, we will show that for any configuration $c$ and its
$\mathcal{T}$ counterpart $\mathcal{T}c$,
$\Psi(c)=(-1)^{N_{dw}(c)}\Psi(\mathcal{T}c)$, where $N_{dw}(c)$ is the number
of domain wall loops for $c$. As we argue in Eq. (7) and (8) in the main text,
such state is symmetric under $\mathcal{T}$.

For the trivial configuration $c$ where all Ising spins points up,
$\lvert\psi_{c}\rangle$ is vacuum state, and $\Psi(c)=\Psi(\mathcal{T}c)=1$.
Let us consider _configuration $c$ with a single domain wall loop_, and assume
this loop contains $2L$ sites and $2L$ bonds. Overlapping $\lvert\Psi\rangle$
with $\lvert c\rangle$, we obtain $\Psi(c)\lvert\psi(c)\rangle$, which is a
new tensor network with site tensors $\hat{T}_{\bm{\mathbf{r}}s}^{c}$ and bond
tensors $\hat{B}^{c}_{\bm{\mathbf{r}}\alpha}$. Any internal leg of this new
tensor network only contains one state, which can be either parity even or
odd. Site and bond tensors away from the loop contains no fermion, and is a
pure spin state with coefficient $1$. We focus on tensors in the loop. We
label site tensors along the loop as $\hat{T}_{j}^{c}$, where $j\in\\{1\cdots
2L\\}$ increases clockwise along this loop. $\hat{B}^{c}_{j,j+1}$ labels the
bond tensor connecting $j$ and $j+1$. $\Psi(c)\lvert\psi(c)\rangle$ is then
schematically drawn as following:

$\displaystyle\Psi(c)\lvert\psi(c)\rangle=\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{\scalebox{1.0}[1.0]{{\leavevmode\hbox{\set@color{}}}}}}}$}~{},$
(S16)

where we only present tensors at the domain wall loop. Here, blue legs are
physical fermions, and light gray legs are internal spins connecting tensors
within domains. Yellow dots are fermion swapping gates. Internal legs
traveling from right to left are colored red, and give an additional $-1$ when
contracting odd parity states according to Eq. (S3).
$\Psi(\mathcal{T}c)\lvert\psi(\mathcal{T}c)\rangle$ can be represented in a
similar way.

Coefficients for these configurations come from two contributions: one from
swapping fermions, and the other one from the $\pm 1$ entries of tensors. Let
$j=2k-1$ to be $u-$sites, and $j=2k$ to be $v-$sites, we calculate these two
contributions respectively in the following.

  * •

We assume that the loop encircles an $\uparrow$ domain. From Eq. (12) in the
main text, $\hat{B}_{2k-1,2k}^{c}$ contains fermion modes $c_{2k-1}$ and
$c_{2k}$, while $\hat{B}_{2k,2k+1}^{c}$ carries zero fermion charge. From Eq.
(S16), contraction between $\hat{B}_{2k-1,2k}^{c}$ and $\hat{T}_{2k-1}^{c}$,
and contraction between $B_{2k-1,2k}^{c}$ and $T_{2k}^{c}$ both give a $-1$.
Thus, no sign factor is produced.

Configuration $\mathcal{T}c$ hosts a $\downarrow$ domain inside the loop.
$\hat{B}_{2k,2k+1}^{\mathcal{T}c}$ contains fermion modes $c_{2k}$ and
$c_{2k+1}$, while $\hat{B}_{2k-1,2k}^{\mathcal{T}c}$ is a pure spin state.
From Eq. (S16), one concludes that for $k<L$, the contraction between
$\hat{B}_{2k,2k+1}^{\mathcal{T}c}$, $\hat{T}_{2k}^{\mathcal{T}c}$ and
$\hat{T}_{2k+1}^{\mathcal{T}c}$ contributes $-1$, and contraction between
$\hat{B}_{2L,1}^{\mathcal{T}c}$, $\hat{T}_{1}^{\mathcal{T}c}$ and
$\hat{T}_{2L}^{\mathcal{T}c}$ contributes no phase factor. So, the sign
difference between $c$ and $\mathcal{T}c$ from fermion contraction is
$(-)^{L-1}$.

  * •

We now calculate contribution from tensor entries. Given site and bond tensors
in the main text, one can check that for any loop with length $2L$, under
$\mathcal{T}$ action, there are always $L$ site tensors and $2L$ bond tensors
on the loop change signs. This can be seen by noticing that on a domain wall
loop, there are always same number of $u-$sites and $v-$sites that locate
between domain walls of an $x-$bond and a $y-$bond (or $x-$bond and $z-$bond,
or $y-$bond and $z-$bond), while sign factors of these $u-$ and $v-$sites are
opposite under $\mathcal{T}$ action. Thus, tensor entries contribute
$(-1)^{L}$ to the sign difference between $c$ and $\mathcal{T}c$.

Combining above contributions, we conclude that $\Psi(c)=-\Psi(\mathcal{T}c)$
for $c$ with a single domain wall.

For configurations with multiple domain wall loops, one can first moving
tensors belonging to a single loop together. Coefficients for each loop can
then be calculated one by one. So, for any configuration $c$ with $N_{dw}(c)$
domain wall loops, we have $\Psi(\mathcal{T}c)=(-)^{N_{dw}(c)}\Psi(c)$.

##  III Time reversal symmetry and Kasteleyn orientation

When acting $\mathcal{T}$ on the honeycomb fPEPS, we get $W(\mathcal{T})$’s
and $W^{-1}(\mathcal{T})$’s on internal legs according to Eq. (14) in the main
text. For any internal leg $(s\alpha)$, $W_{(s\alpha)}(\mathcal{T})$ in
general are not in the neighbourhood of $W^{-1}_{(s\alpha)}(\mathcal{T})$, and
one should move them together to cancel each other in the tensor contraction.
Since $W(\mathcal{T})$ does not have fixed parity, permuting
$W(\mathcal{T})$’s and $W^{-1}(\mathcal{T})$’s lead to fermion-swap gates. Eq.
(14) in the main text in general does not lead to a $\mathcal{T}$-symmetric
physical wavefunction due to these swap gates. In this section, we will show
that local $\mathcal{T}$ constraint gives a symmetric physical wavefunction if
and only if a hidden Kasteleyn orientation can be extracted from such local
constraint.

Let us first introduce Kasteleyn orientation. For a given oriented graph,
Kasteleyn orientation means that for any face in such graph, the number of
clockwise-oriented edges bounding it must be odd. It has been proven that
Kasteleyn orientation exists for any planar graph with an even number of
verticesCimasoni and Reshetikhin (2007). The choice of Kasteleyn orientation
is far from unique: given an Kasteleyn orientation, one can obtain another one
by flipping arrows on all edges connecting a given vertex $v$. For each
vertex, one can choose to flip or not to flip arrows on edges connecting this
vertex, and it exhausts all possible choices of Kasteleyn orientation. Thus,
there are total $2^{N_{v}}$ choices, where $N_{v}$ is the total number of
vertices.

We now turn to the honeycomb example in the main text. Let us explain rules to
extract orientation from action of $W(\mathcal{T})$’s on local tensors. We
first expand all sites of the honeycomb lattice to triangles as

$\displaystyle\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}\rightarrow\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}$
(S17)

Vertices in the new lattice live on internal legs of the honeycomb fPEPS, and
thus can be labeled as $(s\alpha)$. Arrows on blue triangles in Eq. (S17) are
extracted from $\mathcal{T}$-action on site tensors, while arrows on magenta
lines in Eq. (S17) of the original hexagons from $\mathcal{T}$-action on bond
tensors. More specifically, from $\mathcal{T}$-symmetric condition in the main
text, for $\mathcal{T}$-action on tensor at site $s$

$\displaystyle
U_{s}(\mathcal{T})\cdot\hat{T}_{s}^{*}=\left[W_{(s\alpha)}(\mathcal{T})\otimes_{f}W_{(s\beta)}(\mathcal{T})\otimes_{f}W_{(s\gamma)}(\mathcal{T})\right]^{-1}\cdot\hat{T}_{s}\,,$
(S18)

the arrows on the corresponding blue triangle is
$(s\alpha)\rightarrow(s\beta)$, $(s\beta)\rightarrow(s\gamma)$, and
$(s\alpha)\rightarrow(s\gamma)$. For example,

$\displaystyle
U_{s}(\mathcal{T})\cdot\hat{T}^{*}_{s}=\left[W_{(sx)}(\mathcal{T})\otimes_{f}W_{(sy)}(\mathcal{T})\otimes_{f}W_{(sz)}(\mathcal{T})\right]^{-1}\cdot\hat{T}_{s}~{}~{}~{}\Rightarrow\scalebox{2.0}[2.0]{\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}}\,.$
(S19)

If $\mathcal{T}$-action on bond tensor $\hat{B}_{\alpha}$ reads

$\displaystyle V_{(\alpha 0)}(\mathcal{T})\otimes_{f}V_{(\alpha
1)}(\mathcal{T})\cdot\left(\hat{B}_{\alpha}\right)^{*}=\hat{B}_{\alpha}\cdot
W_{(s\alpha)}(\mathcal{T})\otimes_{f}W_{(s^{\prime}\alpha)}(\mathcal{T})\,,$
(S20)

the arrow is drawn from $(s\alpha)$ to $(s^{\prime}\alpha)$. Following this
rule, $\mathcal{T}$-symmetry constraint in the main text leads to the oriented
graph in Eq. (S17), and one can check that it indeed gives one Kasteleyn
orientation.

In the following, let us prove that the physical wavefunction is
$\mathcal{T}$-symmetric only when the right hand side of Eq. (S17) is a
Kasteleyn orientation. According to Eq. (16) in the main text, the commutator
between $W(\mathcal{T})$ and parity $F$ reads

$\displaystyle W^{*}(\mathcal{T})\cdot
F\cdot\left[W^{-1}(\mathcal{T})\right]^{*}=D\cdot F~{},\quad\text{where
}D=\exp[\operatorname{i}\pi n_{D}^{*}]$ (S21)

Here, $D$ can be expressed as $P_{e}-P_{o}$, where $P_{e/o}$ are projectors to
internal states without/with Ising domain wall. Using these projectors, we can
decompose site and bond tensors to orthogonal sectors, e.g.
$P_{e,x}\otimes_{f}P_{o,y}\otimes_{f}P_{o,z}\cdot\hat{T}$. By inserting
$P_{e}+P_{o}$ on all internal legs of the tensor network wavefunction, the
wavefunction equals summation of tensor contractions for different sectors.
However, as $D$ is an IGG element of the tensor network, local tensors vanish
if acted by odd number of $P_{o}$’s. So, for the whole tensor network,
$P_{o}$’s must form loops, which gives domain wall loops for physical Ising
spins after contraction.

Given a loop configuration $d$, we act $P_{o}$’s on internal legs along loops
and $P_{e}$’s on internal legs within domains, and then obtain site tensors
$\hat{T}^{d}$’s and bond tensors $\hat{B}^{d}$’s. By contracting internal legs
of $\hat{T}^{d}$’s and $\hat{B}^{d}$’s, we get physical state
$\lvert\psi_{d}\rangle$. Physical wavefunction is obtained by summation over
all loop configurations $\lvert\Psi\rangle=\sum_{d}\lvert\psi_{d}\rangle$.

In the following, we will prove that $\forall
d,~{}\mathcal{T}\lvert\psi_{d}\rangle=\lvert\psi_{d}\rangle$ if there is
hidden Kasteleyn orientation, and thus $\lvert\Psi\rangle$ is invariant under
$\mathcal{T}$.

Let us start from configuration $d$ without any loop. In this case, we only
get $W_{e}(\mathcal{T})$’s when acting $\mathcal{T}$, which are parity even
and free to permute. $\lvert\psi_{d}\rangle$ is apparently invariant under
$\mathcal{T}$: $\mathcal{T}\lvert\psi_{d}\rangle=\lvert\psi_{d}\rangle$.

For configuration $d$ with a single loop, let the number of internal ket legs
to be $2L$. We label internal legs along this loop counter-clockwise by $l$,
and assume site tensors sitting between $l=2k-1$ and $l=2k$. Site tensors
along this loop for $\lvert\psi_{d}\rangle$ can then be named as
$\hat{T}_{2k-1,2k}^{d}$, while bond tensors as $\hat{B}_{2k,2k+1}^{d}$. As all
tensors are parity even, we rearrange tensors along the loop together in the
following way:

$\displaystyle\lvert\psi_{d}\rangle=\mathrm{fTr}[\cdots\otimes_{f}\hat{B}_{2,3}^{d}\otimes_{f}\hat{B}_{3,4}^{d}\otimes_{f}\cdots\otimes_{f}\hat{B}_{2L,1}^{d}\otimes_{f}\hat{T}_{2L-1,2L}^{d}\otimes_{f}\hat{T}^{d}_{2L-3,2L-2}\otimes_{f}\cdots\otimes_{f}\hat{T}_{1,2}^{d}\otimes_{f}\cdots]$
(S22)

We now act $\mathcal{T}$ on $\lvert\psi_{d}\rangle$, and according to Eq.
(S18) and Eq. (S20), it gives $W(\mathcal{T})$’s on internal legs of
$\hat{T}^{d}$’s and $\hat{B}^{d}$’s. We define $W_{e/o}(\mathcal{T})\equiv
P_{e/o}\cdot W(\mathcal{T})$, which is parity even/odd sector of
$W(\mathcal{T})$. From definition of $\lvert\psi_{d}\rangle$,
$W(\mathcal{T})$’s act as $W_{o}(\mathcal{T})$’s on internal legs along the
loop, while acting as $W_{e}(\mathcal{T})$’s on internal legs away from the
loop. Fermion signs come from permuting $W_{o}(\mathcal{T})$’s and
$[W_{o}(\mathcal{T})]^{-1}$’s, and thus we focus on contraction of
$W(\mathcal{T})$’s along the loop. We arrange the order of
$W_{o}(\mathcal{T})$’s contraction according to Eq. (S22) as

$\displaystyle\mathrm{fTr}\bigg{\\{}$ $\displaystyle\Big{(}(-1)^{s_{2,3}}\cdot
W_{o,2}(\mathcal{T})\otimes_{f}W_{o,3}(\mathcal{T})\Big{)}\otimes_{f}\cdots\otimes_{f}\Big{(}(-1)^{s_{2L,1}}\cdot
W_{o,2L}(\mathcal{T})\otimes_{f}W_{o,1}(\mathcal{T})\Big{)}\bigotimes\mathop{{}_{\mkern-4.0muf}}$
$\displaystyle\bigg{[}\Big{(}(-1)^{s_{1,2}}\cdot
W_{o,1}(\mathcal{T})\otimes_{f}W_{o,2}(\mathcal{T})\Big{)}\otimes_{f}\cdots\otimes_{f}\Big{(}(-1)^{s_{2L-1,2L}}\cdot
W_{o,2L-1}(\mathcal{T})\otimes_{f}W_{o,2L}(\mathcal{T})\Big{)}\bigg{]}^{-1}\bigg{\\}}$
$\displaystyle={}$ $\displaystyle(-1)^{1+\sum_{l}s_{l,l+1}}=1$ (S23)

where $s_{l,l+1}=0/1$ if the arrow at $(l,l+1)$ is along/against the direction
of the loop (counter-clockwise/clockwise direction). The last equation is from
the definition of Kasteleyn orientation: there are always odd number of arrows
against direction of the loop. So, for configuration $d$ with a single loop,
$W(\mathcal{T})$’s and $[W(\mathcal{T})]^{-1}$’s cancels, and
$\lvert\psi_{d}\rangle$ is $\mathcal{T}$-symmetric. In contrast, if the
orientation is not Kasteleyn, one can always find a loop configuration $d$,
such that the last line of Eq. (S23) gives $-1$, making $\lvert\Psi\rangle$
break $\mathcal{T}$.

For configuration $d$ with multiple loops, we can arrange all tensors
belonging to one loop together, and repeat the above calculation for every
loop. Thus, such $\lvert\psi_{d}\rangle$ is also $\mathcal{T}$-symmetric. In
conclusion, $\lvert\Psi\rangle$ is $\mathcal{T}$-symmetric if and only if the
orientation extracted from $W(\mathcal{T})$’s is a Kasteleyn orientation.

As we mentioned in the beginning of this part, by flipping arrows on all edges
connecting to certain vertices, one gets another Kasteleyn orientation. In the
tensor language, flipping arrows for edges connecting to vertex $(s\alpha)$
corresponds to modifying $W_{(s\alpha)(\mathcal{T})}$ to $D_{(s\alpha)}\cdot
W_{(s\alpha)}(\mathcal{T})$. To see this, we consider the following
$\mathcal{T}$-action on site tensors:

$\displaystyle\left[D_{(sx)}\cdot
W_{(sx)}(\mathcal{T})\otimes_{f}W_{(sy)}(\mathcal{T})\otimes_{f}W_{(sz)}(\mathcal{T})\right]^{-1}\cdot\hat{T}_{s}=\left[W_{(sy)}(\mathcal{T})\otimes_{f}W_{(sz)}(\mathcal{T})\otimes_{f}W_{(sx)}\right]^{-1}\cdot\hat{T}_{s}~{}~{}~{}\Rightarrow\scalebox{2.0}[2.0]{\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}}$

Similar logic works for bond tensors. Following these rules for arrows, there
is one-to-one correspondence between gauge transformation $W(\mathcal{T})$ and
Kasteleyn orientation.

We focus on trivalent lattice in the above argument, where each vertex
connects three bonds, and there is at most one domain wall travelling through
a vertex. We now generalize the above argument to generic lattices, where more
than one domain walls may meet at sites.

Let us present rules for extracting orientations in generic lattices. Similar
as Eq. (S17), we first construct a new planar graph, where a site connecting
$n$-bonds in the original lattice splits to $n$ vertices in the new graph, and
each vertex is labeled by the internal leg index. Each pair of these $n$
vertices are connected by new edges. Given $\mathcal{T}$-action on internal
legs, arrows on edges of the new planar graph follows similar rules presented
in Eq. (S19) and Eq. (S20). For example, consider site tensor $s$ with four
internal legs, arrows can be read from $W(\mathcal{T})$ action as

$\displaystyle
U(\mathcal{T})\cdot\hat{T}^{*}_{s}=\left[W_{(sa)}(\mathcal{T})\otimes_{f}W_{(sb)}(\mathcal{T})\otimes_{f}W_{(sc)}\otimes_{f}W_{(sd)}\right]^{-1}\cdot\hat{T}_{s}~{}~{}~{}\Rightarrow\scalebox{2.0}[2.0]{\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}}$

To proceed, let us focus on a particular choice of Kasteleyn orientation. We
number the $n$ vertices from $1$ to $n$ clockwise, and let the arrow pointing
from $i$ to $j$ if $i<j$. It is easy to verify that any loop within these $n$
vertices matches the condition for Kasteleyn orientation. Arrows on bond
tensors are chosen to satisfy conditions for Kasteleyn orientation on larger
loops.

For configurations without loop crossing, using similar argument presented in
the honeycomb case, we conclude that $W(\mathcal{T})$’s cancels without
additional sign.

Figure 4: Two domain wall loops (red and blue) intersect at site $s_{1}$ and
$s_{2}$. Here, $s_{1,2}$ are not “true crossing points” of these two loops:
they are separate at these two points.

.

We consider configuration $d$ where two loops (colored blue and red) intersect
at site $s_{1}$ and $s_{2}$. As shown in Fig. 4, colors are chosen such that
loops can be separate at these two sites, and there are no “true crossing
points” between these two loops.

By inserting $P_{o}$’s on internal legs at these two loops, and $P_{e}$’s on
other internal legs, we obtain tensors $\hat{T}^{d}$’s and $\hat{B}^{d}$’s,
and wavefunction $\lvert\psi_{d}\rangle$ from contracting $\hat{T}^{d}$’s and
$\hat{B}^{d}$’s. Let the number of internal ket legs of the red/blue loop to
be $2L_{1}/2L_{2}$. We label internal legs along these two loops by index $l$,
where $1\leq l\leq 2L_{1}$ labels internal legs for the red loop and
$2L_{1}+1\leq l\leq 2L_{1}+2L_{2}$ for the blue loop. Due to the intersecting
sites $s_{1/2}$, Eq. (S22) do not directly apply to the case here. However, at
any intersection point, two $W_{o}^{-1}$’s belonging to one loop are always
neighbour to each other (or can be moved as neighbour without sign). Thus, we
can always move all $W_{o}^{-1}$’s belonging to one loop together without
extra sign, and Eq. (S23) still applies for every single loop.

The above argument can be easily generalized to any loop configurations, as we
can always choose loops such that there are no “true crossing points”. We
conclude that for a $\mathcal{T}$ symmetric fPEPS of the QSH phase in any
lattice, the planar graph extracted from $W(\mathcal{T})$’s satisfies
Kasteleyn orientation.

##  IV $\mathcal{T}^{2}=F$ on internal legs

In this part, we show that Eq. (18) in the main text is consistent with
$\mathcal{T}^{2}=F$ when acting on tensors. We consider a tensor $\hat{T}_{A}$
with $L$ internal bra legs and physical Hilbert space $\mathbb{H}_{A}$.
According to the main text, the $\mathcal{T}$-symmetric condition for
$\hat{T}_{A}$ reads

$\displaystyle
U_{A}(\mathcal{T})\cdot\hat{T}_{A}^{*}=\hat{T}_{A}\cdot\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W_{j}(\mathcal{T})$
(S24)

By acting $\mathcal{T}$ twice, we obtain

$\displaystyle F_{A}\cdot\hat{T}_{A}=U_{A}(\mathcal{T})\cdot
U_{A}^{*}(\mathcal{T})\cdot\hat{T}_{A}=\hat{T}_{A}\cdot\left(\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W_{j}(\mathcal{T})\right)\cdot\left(\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W^{*}_{j}(\mathcal{T})\right)$
(S25)

Due to the unfixed parity of $W(\mathcal{T})$’s, one cannot permute
$W(\mathcal{T})$’s directly. Instead, from Eq. (S21), we have

$\displaystyle
F_{A}\cdot\hat{T}_{A}=\hat{T}_{A}\cdot\left(\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W_{j}(\mathcal{T})\right)\cdot\left(\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W^{*}_{j}(\mathcal{T})\right)$
$\displaystyle={}$
$\displaystyle\hat{T}_{A}\cdot\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}(-)^{n_{D,j}\sum_{k<j}n_{D,k}}W_{j}(\mathcal{T})\cdot
W_{j}^{*}(\mathcal{T})=\hat{T}_{A}\cdot(-)^{\sum_{k<j}n_{{D,k}}n_{{D,j}}}\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W_{j}(\mathcal{T})\cdot
W_{j}^{*}(\mathcal{T})$ $\displaystyle={}$
$\displaystyle\hat{T}_{A}\cdot\mathrm{e}^{-\operatorname{i}\frac{\pi}{2}\sum_{j\neq
k}n_{{D,k}}n_{{D,j}}}\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}W_{j}(\mathcal{T})\cdot
W_{j}^{*}(\mathcal{T})=\hat{T}_{A}\cdot\bigotimes_{j=1}^{L}\mathop{{}_{\mkern-4.0muf}}\mathrm{e}^{\operatorname{i}\frac{\pi}{2}n_{{D,j}}^{2}}W_{j}(\mathcal{T})\cdot
W_{j}^{*}(\mathcal{T})\,,$ (S26)

where the last line comes from the condition that $\sum_{j}n_{{D,j}}=0$ when
acting on $\hat{T}_{A}$. Then on each internal leg we have

$\displaystyle\exp\left[\operatorname{i}\frac{\pi}{2}n_{D}^{2}\right]\cdot
W(\mathcal{T})\cdot W^{*}(\mathcal{T})=F$ (S27)

or at most differ from $F$ up to an IGG element. So, by acting
$\mathcal{T}^{2}=F$ on internal legs of tensors, we get consistent result as
Eq. (18) in the main text.

##  V Edge theories from infinite PEPS

In this section, we will identify Hilbert space and symmetry action of the
edge theory from infinite PEPS.

We cut a finite region $A$ from an infinite PEPS. By contracting all internal
legs within $A$, we obtain a linear map $\hat{T}_{A}$ from virtual legs at
boundary of $A$ – labeled as $\mathbb{H}_{\partial A}$ – to physical legs in
the bulk of $A$ – labeled as $\mathbb{H}_{A}$:

$\displaystyle\hat{T}_{A}=\sum(T_{A})_{i_{b}i_{e}}\lvert
i_{b}\rangle\left(i_{e}\right|~{},\quad\lvert
i_{b}\rangle\in\mathbb{H}_{A},~{}\left(i_{e}\right|\in\mathbb{H}_{\partial A}$
(S28)

Here, without loss of generality, we choose all boundary legs to be bra
spaces. For large enough region $A$,
$\dim\mathbb{H}_{A}\gg\dim\mathbb{H}_{\partial A}$, so the map can never be
surjective.

We can write down a symmetric Hamiltonian for a system on $A$, whose low-
energy space is image of $\hat{T}_{A}$, which is isomorphic to
$\mathbb{H}_{\partial A}/\ker{\hat{T}_{A}}$. As bulk excitations are gapped,
low energy states are identified as edge modes, and thus
$\mathbb{H}_{edge}\cong\mathbb{H}_{\partial A}/\ker{\hat{T}_{A}}$.

If $\hat{T}_{A}$ is injective, $\mathbb{H}_{edge}=\mathbb{H}_{\partial A}$,
and it naturally leads to a tensor product structure of the edge Hilbert
space. If IGG is nontrivial, $\hat{T}_{A}$ will no longer be injective. Given
an IGG element whose action on $\partial A$ is $J_{\partial A}$, according to
the definition of IGG, $T_{A}\cdot(\hat{1}_{\partial A}-J_{\partial A})=0$,
and we have $\ker{\hat{T}_{A}}\supset\operatorname{imag}(\hat{1}_{\partial
A}-J_{\partial A})\neq 0$.

In this work, we further assume that
$\ker{\hat{T}_{A}}=\left\\{\operatorname{imag}(\hat{1}_{\partial
A}-J_{\partial A})\middle|\forall J\in\mathrm{IGG}\right\\}$. In other words,
_IGG determines the edge Hilbert space_ :

$\displaystyle\mathbb{H}_{edge}=\left\\{\lvert\psi_{\partial
A}\rangle~{}\middle|~{}J_{\partial A}\lvert\psi_{\partial
A}\rangle=\lvert\psi_{\partial A}\rangle,~{}\forall J\in\mathrm{IGG}\right\\}$
(S29)

We define $\hat{T}_{A}^{-1}:\mathbb{H}_{A}\to\mathbb{H}_{\partial A}$ as
pseudo-inverse of $\hat{T}_{A}$, which satisfies

$\displaystyle\hat{T}_{A}\cdot\hat{T}_{A}^{-1}=P_{l}~{},\quad\hat{T}_{A}^{-1}\cdot\hat{T}_{A}=P_{edge}~{},$
(S30)

where $P_{l}$ is the projector from $\mathbb{H}_{A}$ to
$\operatorname{imag}\hat{T}_{A}$, while $P_{edge}$ is the projector from
$\mathbb{H}_{\partial A}$ to $\mathbb{H}_{edge}$.

In the following, let us work out how symmetries act on edge. Here, we focus
on onsite symmetry group $G$. For $g\in G$, we have

$\displaystyle U_{A}(g)\mathcal{K}^{s(g)}\cdot T_{A}=T_{A}\cdot W_{\partial
A}(g)\mathcal{K}^{s(g)}$ (S31)

So, it is natural to identify $U_{edge}(g)\mathcal{K}^{s(g)}\equiv
P_{edge}\cdot W_{\partial A}(g)\mathcal{K}^{s(g)}\cdot P_{edge}$ as symmetry
action on $\mathbb{H}_{edge}$.

Note that $\forall J\in\mathrm{IGG}$, we have

$\displaystyle U_{A}(g)\mathcal{K}^{s(g)}\cdot\hat{T}_{A}=\hat{T}_{A}\cdot
W_{\partial
A}\mathcal{K}^{s(g)}=U_{A}(g)\mathcal{K}^{s(g)}\cdot\hat{T}_{A}\cdot
J_{\partial A}=\hat{T}_{A}\cdot W_{\partial A}\mathcal{K}^{s(g)}\cdot
J_{\partial A}$ (S32)

Together with Eq. (S29), we conclude

$\displaystyle U_{edge}(g)\mathcal{K}^{s(g)}=P_{edge}\cdot W_{\partial
A}(g)\mathcal{K}^{s(g)}=W_{\partial A}(g)\mathcal{K}^{s(g)}\cdot P_{edge}$
(S33)

##  VI Fusion of $\mathcal{T}-$flux

As in Sec. V, let us consider a region $A$ described by a large tensor
$\hat{T}_{A}$ with physical Hilbert space $\mathbb{H}_{A}$ and boundary legs
$\mathbb{H}_{\partial A}$, where legs at $\partial A$ are labeled by
$j=\\{1,2,\dots,L\\}$.

We create $\mathcal{T}$-flux at $j=1$ and $l+1$ by inserting a charge-neutral
operator $U_{M}(\mathcal{T})\mathcal{K}\equiv P_{edge}\cdot w_{l+1}\cdot
w_{1}\cdot W_{M}(\mathcal{T})\mathcal{K}\cdot P_{edge}$, where
$M=\\{2,3,\dots,l\\}\in\partial A$. Here, $w_{1/(l+1)}$ are local operators at
the ends and

$W_{M}(\mathcal{T})=\bigotimes_{j=2}^{l}\mathop{{}_{\mkern-4.0muf}}W_{j}(\mathcal{T})$

We choose $w_{1/(l+1)}$ such that $w_{l+1}\cdot w_{1}\cdot
W_{M}(\mathcal{T})\mathcal{K}$ commute with $P_{edge}$.

To make $U_{M}(\mathcal{T})K$ charge neutral, we require commutator between
$w_{1/(l+1)}$ and $n_{f}$ to be

$\displaystyle[w_{1},n_{{}_{f;1}}]={n_{\lambda;\frac{3}{2}}^{(0)}}\cdot
w_{1}~{};\quad[w_{l+1},n_{{}_{f;l+1}}]={n_{\lambda;l+\frac{1}{2}}^{(1)}}\cdot
w_{l+1}$ (S34)

The charge neutral condition is satisfied as
$[W_{j}(\mathcal{T}),n_{{}_{f;j}}]=n_{D;j}\cdot W_{j}(\mathcal{T})$ and
$\left(n_{\lambda;\frac{3}{2}}^{(0)}+n_{\lambda;l+\frac{1}{2}}^{(1)}+\sum_{j=2}^{l}n_{D,j}\right)\cdot
P_{edge}=0$.

We now calculate fusion of two $\mathcal{T}$-flux. Without loss of generality,
we assume that $j=1$ and $j=l+1$ legs belong to $v$-sublattice. By performing
similar derivation presented in Eq. (S26), we have

$\displaystyle U_{M}(\mathcal{T})\cdot U_{M}^{*}(\mathcal{T})$ (S35)
$\displaystyle=$ $\displaystyle P_{edge}\cdot\mathrm{e}^{\operatorname{i}\pi
n_{\lambda;\frac{3}{2}}^{(0)}}\mathrm{e}^{\operatorname{i}\frac{\pi}{2}\left(n_{\lambda;l+\frac{1}{2}}^{(1)}\right)^{2}}w_{l+1}\cdot
w_{l+1}^{*}\otimes_{f}\mathrm{e}^{\operatorname{i}\frac{\pi}{2}\left(n_{\lambda;\frac{3}{2}}^{(0)}\right)^{2}}w_{1}\cdot
w_{1}^{*}\bigotimes_{j=2}^{l}\mathop{{}_{\mkern-4.0muf}}\mathrm{e}^{\operatorname{i}\frac{\pi}{2}n_{D,j}^{2}}W_{j}(\mathcal{T})\cdot
W_{j}^{*}(\mathcal{T})\cdot P_{edge}\,.$

As
$\left(n_{\lambda;l+\frac{1}{2}}^{(1)}\right)^{2}=n_{\lambda;l+\frac{1}{2}}^{(1)}$
and
$\left(n_{\lambda;\frac{3}{2}}^{(0)}\right)^{2}=-n_{\lambda;\frac{3}{2}}^{(0)}$
for $v-$sites, we obtain

$\displaystyle(U_{M}(\mathcal{T})\mathcal{K})^{2}=P_{edge}\cdot\Omega_{l+1}\Omega_{1}\cdot\prod_{j=2}^{l}F_{j}\cdot
P_{edge}\,,$ (S36)

Here,

$\displaystyle\Omega_{l+1}=$
$\displaystyle{}P_{edge}\cdot\mathrm{e}^{\operatorname{i}\frac{\pi}{2}n_{\lambda;l+\frac{1}{2}}^{(1)}}w_{l+1}\cdot
w_{l+1}^{*}\cdot
P_{edge}=\mathrm{e}^{\operatorname{i}\frac{\pi}{2}n_{\lambda;l+\frac{1}{2}}^{(1)}}c^{\dagger}_{l+1}\,,$
(S37) $\displaystyle\Omega_{1}=$
$\displaystyle{}P_{edge}\cdot\mathrm{e}^{\operatorname{i}\frac{\pi}{2}n_{\lambda;\frac{3}{2}}^{(0)}}w_{1}\cdot
w_{1}^{*}\cdot
P_{edge}=\mathrm{e}^{\operatorname{i}\frac{\pi}{2}n_{\lambda;\frac{3}{2}}^{(0)}}c_{1}\,.$
(S38)

where we use Eq. (24) in the main text to fix the final form of
$\Omega_{1/(l+1)}$.

##  VII Variational ansatz for the QSH phase on square lattice spin-1/2
electronic system

In this section, we present detailed derivation for solving tensor equations
on a spin-$\frac{1}{2}$ fermionic system on the bipartite square lattice. Each
site tensor has four internal legs and one physical spin-$\frac{1}{2}$
fermion, which is represented as

$\displaystyle\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}=T_{ijkl,p}\left|i\right)_{1}\left|j\right)_{2}\left|k\right)_{3}\left|l\right)_{4}\lvert
p\rangle\,.$

Sites on the bipartite square lattice can be grouped to two types, where we
use $\hat{T}_{u/v}$ to label site tensors on different sublattices. Physical
legs are spin-1/2 fermions $f_{\sigma}$, and charges carried by $f_{\sigma}$’s
are opposite in $u-$ and $v-$site. As in the honeycomb example, an internal
leg $(s\alpha)$ can also be represented by a triple-line, where the middle
line is a spinless fermion mode $c_{(s\alpha)}$, and two side lines are Ising
spins $\left|\tau\right)_{(s\alpha a)}$, where $s=u/v$,
$\alpha\in\\{1,2,3,4\\}$ and $a=0/1$. Orders of $a$ is counter-
clockwise/clockwise for site $u/v$. For simplicity, we focus on four of eight
basis states:
$\\{\left|\uparrow\uparrow\right)_{(s\alpha)},~{}c_{(s\alpha)}^{\dagger}\left|\uparrow\downarrow\right)_{(s\alpha)},~{}\left|\downarrow\uparrow\right)_{(s\alpha)},~{}\left|\downarrow\downarrow\right)_{(s\alpha)}\\}$.

Symmetry constraints on site tensor $\hat{T}_{s}$ are

$\displaystyle\bigotimes_{\alpha=1}^{4}\mathop{{}_{\mkern-4.0muf}}W_{(s\alpha)}(\mathcal{T})\otimes_{f}U_{s}(\mathcal{T})\cdot\hat{T}^{*}_{s}$
$\displaystyle=\hat{T}_{s}$
$\displaystyle\Big{[}\sum_{\alpha=1}^{4}n_{f;(s\alpha)}+n_{f;s}\Big{]}\cdot\hat{T}_{s}$
$\displaystyle=0$ (S39)

where $W(\mathcal{T})$, $n_{f}$, and $n_{\lambda}$ take the same form as those
in the honeycomb example. To take care of sign factors when acting
$W(\mathcal{T})$’s on site tensor $\hat{T}_{s}$, a Jordan-Wigner string
$J=\exp\\{\operatorname{i}\pi n_{f}\\}$ is introduced. Then, the “bosonized”
$W(\mathcal{T})$ is

$\displaystyle W^{b}_{(sj)}(\mathcal{T})=P_{e}\cdot
W_{(sj)}(\mathcal{T})+\bigotimes_{k=1}^{j-1}J_{k}~{}\otimes P_{o}\cdot
W_{(sj)}(\mathcal{T})\,,$ (S40)

$\mathcal{T}$ symmetry constraint on a site tensor gives

$\prod_{\alpha=1}^{4}W_{(s\alpha)}^{b}(\mathcal{T})\cdot U(\mathcal{T})\cdot
T^{*}_{s}=T_{s}$ (S41)

Here $T_{s}$ without hat is a “bosonic tensor” with entries
$(T_{s})_{\alpha\beta\gamma\delta,p}$.

Site tensor $\hat{T}_{s}$ should also satisfy the plaquette IGG condition:

$\displaystyle\Big{(}n_{\lambda;(s\alpha
0)}+n_{\lambda;(s\widetilde{\alpha}1)}\Big{)}\cdot\hat{T}_{s}=0,~{}\forall\alpha\,.$
(S42)

where $\widetilde{\alpha}=\alpha+(-1)^{s}$. This equation identifies internal
states within a plaquette:

$\displaystyle\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}\to\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}\,.$

The dimension of the tensor $T_{s}$ is 1024, and after solving the tensor
equations in Eq. (S39) and (S42) there are only 14 linearly independent
solutions. And the solution for $\hat{T}_{v}$ reads

$\displaystyle\hat{T}_{v}=\sum_{l}^{14}c_{l}\hat{t}_{l}\,,$ (S43)

where $c_{l}$’s are real numbers, and $\hat{t}_{l}$ can be represented
graphically as

Magenta arrow indicates a fermion on the internal leg, and the
$\uparrow/\downarrow$ at 4 corners are the identified internal states within
the same plaquette. $\hat{T}_{u}$’s solution is given by flipping all
plaquette Ising spins of $\hat{T}_{v}$.

Now, let us discuss the bond tensors. As shown in Sec. III, to make the tensor
network $\mathcal{T}$ symmetric, $\mathcal{T}$ action on bonds should be
chosen to satisfy Kasteleyn orientation:

$\displaystyle\hat{B}_{\alpha}\cdot
W_{(u\alpha)}(\mathcal{T})\otimes_{f}W_{(v\beta)}(\mathcal{T})$
$\displaystyle=\hat{B}_{\alpha}^{*}$ (S44)

Here, we impose rotational symmetry, and thus all four types of bond tensors
share the same form:

$\displaystyle\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}\,=\,\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}\,=\,\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}\,=\,\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}=B_{\alpha\beta}\left(\alpha\right|\left(\beta\right|$
(S45)

Bond tensors should also be invariant under plaquette IGG:

$\displaystyle\Big{(}n_{\lambda;(u\alpha a)}+n_{\lambda;(v\alpha
a)}\Big{)}\cdot\hat{B}_{\alpha}$
$\displaystyle=0~{},\quad\forall\alpha~{}\&~{}a~{},$ (S46)

which identifies the internal states within a plaquette:

$\displaystyle\scalebox{1.0}[1.0]{\leavevmode\hbox{$\vbox{\hbox{\resizebox{}{}{{\leavevmode\hbox{\set@color{}}}}}}$}}\,.$
(S47)

By further imposing charge neutral condition, we get solution for
$\hat{B}_{\alpha}$ as

$\displaystyle\hat{B}_{\alpha}=d_{e}\hat{b}_{e}+d_{o}\hat{b}_{o}\,,$ (S48)

where
$\hat{b}_{e}=\left(\uparrow\uparrow\right|_{u}\left(\uparrow\uparrow\right|_{v}+\left(\downarrow\downarrow\right|_{u}\left(\downarrow\downarrow\right|_{v}$,
and
$\hat{b}_{o}={\left(\uparrow\downarrow\right|_{u}\left(\uparrow\downarrow\right|_{v}}c_{u}c_{v}-\left(\downarrow\uparrow\right|_{u}\left(\downarrow\uparrow\right|_{v}$.
Here, $d_{e/o}$ are real parameters. By performing gauge transformation, we
can always absorb $d$’s to site tensors and simply set $d_{e}=d_{o}=1$.
Namely, bond tensors are maximal entangled states, which share the same form
as Eq. (12) in the main text.

##  VIII Numerical calculation of many-body topological invariants

In this appendix, we present the ansatz tensor for the honeycomb lattice with
a bond dimension $4$ and perform explicit numerical calculations of the many-
body topological invariant proposed by Shiozaki et alShiozaki  _et al._
(2018). This analysis demonstrates the existence of a parameter region
corresponding to the Quantum Spin Hall (QSH) phase.

By solving the symmetry and $\mathrm{IGG}$ restrictions (Eq. (S39) and (S42)),
we get ansatz tensors. The $u$-site variational tensor can be written as:

$\displaystyle\hat{T}_{u}=\sum_{l}c^{u}_{l}t^{u}_{l}$ (S49)

where $c^{u}_{l}$’s are real numbers, and $\hat{t}^{u}_{l}$’s are shown
graphically below:

For $v$-site variational tensor, we have

$\displaystyle\hat{T}_{v}=\sum_{l}c^{v}_{l}t^{v}_{l}$ (S50)

with $c^{v}_{l}$’s real numbers and $\hat{t}_{n}$’s below

And bond tensor is chosen to be
$\hat{B}=\left(\uparrow\uparrow\right|_{u}\left(\uparrow\uparrow\right|_{v}+\left(\downarrow\downarrow\right|_{u}\left(\downarrow\downarrow\right|_{v}-{\left(\uparrow\downarrow\right|_{u}\left(\uparrow\downarrow\right|_{v}}c_{u}c_{v}+\left(\downarrow\uparrow\right|_{u}\left(\downarrow\uparrow\right|_{v}$.
In contrast to the fixed-point wavefunction discussed in the main text, the
variational wavefunction remains a superposition of loop configurations with
decorated fermions. However, it now exhibits distinct coefficients and does
not involve any additional physical Ising spins. Notably, the spin degrees of
freedom become decoupled from the fermion degrees of freedom, existing solely
as internal degrees of freedom.

To calculate the topological invariant, we set the $y$ direction of our 2+1D
wavefunction to be periodic, and keep the $x$ direction open to get a
cylindrical geometry. Then we divide the cylinder as below

The topological invariant of state $|\phi\rangle$ is obtained from the
following formula:

$\displaystyle Z$ $\displaystyle=\mathrm{Tr}[\rho^{+}_{R_{1}\cup
R_{3}}C^{R_{1}}_{T}(\rho^{-}_{R_{1}\cup
R_{3}})^{\mathsf{T}_{1}}[C^{R_{1}}_{T}]^{\dagger}]$
$\displaystyle\rho^{\pm}_{R_{1}\cup R_{3}}$
$\displaystyle=\mathrm{Tr}_{\overline{R_{1}\cup R_{3}}}[\exp{\frac{\pm
2\pi\mathrm{i}y\sum_{\mathrm{r}\in
R_{2}}n(\mathrm{r})}{L_{y}}}|\phi\rangle\langle\phi|]$ (S51)

Where $C^{R_{1}}_{T}(\rho^{-}_{R_{1}\cup
R_{3}})^{\mathsf{T}_{1}}[C^{R_{1}}_{T}]^{\dagger}$ is the time-reversal
partial transpose of $\rho^{-}_{R_{1}\cup R_{3}}$. In Fock space, the time-
reversal partial transpose of $|\\{n_{j}\\}_{j\in R_{1}},\\{n_{j}\\}_{j\in
R_{2}}\rangle\langle\\{\bar{n}_{j}\\}_{j\in R_{1}},\\{\bar{n}_{j}\\}_{j\in
R_{2}}|$ is defined asShiozaki  _et al._ (2018)

$\displaystyle(-i)^{[\tau_{1}+\bar{\tau}_{1}]}(-1)^{(\tau_{1}+\bar{\tau}_{1})(\tau_{2}+\bar{\tau_{2}})}U^{R_{1}}_{T}|\\{\bar{n}_{j}\\}_{j\in
R_{1}},\\{n_{j}\\}_{j\in R_{2}}\rangle\langle\\{n_{j}\\}_{j\in
R_{1}},\\{\bar{n}_{j}\\}_{j\in R_{2}}|[U^{R_{1}}_{T}]^{\dagger}$ (S52)

with $[x]=0$ for even $x$ and $[x]=1$ for odd $x$ and $\tau_{1/2}=\sum_{j\in
R_{1/2}}n_{j}$. $U^{R_{1}}_{T}$ is the unitary part of time-reversal action on
region $R_{1}$. The angle of $Z$ is $\pi$ for QSH phase and $0$ for trivial
phase.

We set the length of the $y$-direction, $L_{y}=2$ (two unit cells with four
sites in the $y$-direction). Additionally, we set the length of three middle
regions in the $x$-direction to be $L_{x}$. By tuning $c^{u}_{2}$ while
keeping other coefficients to be $1$, for different $L_{x}$, we get the sign
and amplitude of $Z$ as shown in Fig. 5.

Figure 5: The sign (a) and amplitude (b) of $Z$ versus $c^{u}_{2}$ while other
coefficients are set to be 1.

We can also enter QSH phase by tuning other coefficients, like $c^{v}_{3}$,
while keeping other coefficients to be $1$ as shown in Fig. 6. Tuning other
coefficients by the same means may not give us QSH phase region for the
wavefunction are already out of QSH phase and those coefficients are
irrelevant. But if we allow two coefficients to vary simultaneously, it is
possible to make every coefficients relevant. To see this, we compare the sign
of $Z$ versus $c^{v}_{6}$ and different $c^{u}_{2}$ while keeping other
coefficients to be $1$. As shown in Fig.7 , we can’t find QSH phase region
when $c^{u}_{2}$ and other coefficients are $1$, but it is possible when
$c^{u}_{2}=2$. Consequently, there is a wide parameter region in QSH phase
that can be tuned for variational simulations.

Figure 6: The sign (a) and amplitude (b) of $Z$ versus $c^{v}_{3}$ while other
coefficients are set to be 1. Figure 7: The sign (a)/(c) and amplitude (b)/(d)
of $Z$ versus $c^{v}_{6}$ when $c^{u}_{2}=1$/$c^{u}_{2}=2$. Other coefficients
are set to be $1$.

The results show that there exists a stable parameter region of QSH phase when
the system size is large enough. The emergence of the trivial phase region can
be attributed to two factors. Firstly, it is influenced by the limitation of
system size, as finite-size effects become significant beyond the fixed point.
Secondly, our ansatz possesses the capability to represent long-range
entangled GHZ (Greenberger-Horne-Zeilinger) states because emergent
$\mathrm{IGG}$ cannot be ruled out solely based on local tensor equations.