# Extracting Wilson loop operators and fractional statistics from a
single bulk ground state

Ze-Pei Cian  Department of Physics, University of Maryland, College Park,
Maryland 20742, USA  Joint Quantum Institute, College Park, 20742 MD, USA
Mohammad Hafezi  Department of Physics, University of Maryland, College Park,
Maryland 20742, USA  Joint Quantum Institute, College Park, 20742 MD, USA
Department of Electrical and Computer Engineering, University of Maryland,
College Park, Maryland 20742, USA  Maissam Barkeshli  Department of Physics,
University of Maryland, College Park, Maryland 20742, USA  Joint Quantum
Institute, College Park, 20742 MD, USA  Condensed Matter Theory Center,
Department of Physics, University of Maryland, College Park, 20742 MD, USA

###### Abstract

An essential aspect of topological phases of matter is the existence of Wilson
loop operators that keep the ground state subspace invariant. Here we present
and implement an unbiased numerical optimization scheme to systematically find
the Wilson loop operators given a single ground state wave function of a
gapped Hamiltonian on a disk. We then show how these Wilson loop operators can
be cut and glued through further optimization to give operators that can
create, move, and annihilate anyon excitations. We subsequently use these
operators to determine the braiding statistics and topological twists of the
anyons, yielding a way to fully extract topological order from a single wave
function. We apply our method to the ground state of the perturbed toric code
and doubled semion models with a magnetic field that is up to a half of the
critical value. From a contemporary perspective, this can be thought of as a
machine learning approach to discover emergent 1-form symmetries of a ground
state wave function. From an application perspective, our approach can be
relevant to find Wilson loop operators in current quantum simulators.

##  I Introduction

Topologically ordered phases are gapped quantum phases of matter that cannot
be characterized by local order parameters, but rather by long-range
entanglement and fractional statistics of quasiparticle excitations. For
decades, a major question has been how to properly diagnose and characterize
topological order in a quantum many-body system. While much progress has been
made wen04 ; nayak2008 ; wang2008 ; senthil2015 ; zeng2019 ; barkeshli2019 ;
kawagoe2020microscopic ; bulmashSymmFrac ; aasen2021characterization , an
outstanding question remains: Can we fully extract topological order from a
single bulk ground state wave function, with no access to the Hamiltonian?

Figure 1:  (a) WLO parameterized with a matrix product operator. (b) Numerical
procedure to optimize a WLO. (c) The expectation value of the WLO $W_{\rm FP}$
for the exactly solvable fixed point follows perimeter law when evaluated in
the perturbed ground state: $\langle\psi|W_{\rm FP}|\psi\rangle\sim
e^{-L/\xi}$, where $L$ and $\xi$ are the perimeter of $W_{\rm FP}$ and
correlation length, respectively. However, the expectation value of the
optimized Wilson loop operator does not decrease exponentially with the
perimeter. (d) The expectation value of the Wilson loop operator
$\langle\psi|W|\psi\rangle$ during the optimization iteration described in
(b). The typical total number of iterations is around $400$. $|\psi\rangle$ is
a ground state of the perturbed toric code model with $h_{x}=0.15$,
$h_{z}=0.05$. The Wilson loop operator is a rectangle with side length
$L_{x}=36$ and $L_{y}=6$ and thickness 1.

Apart from the fundamental interest in the above question, there is a growing
body of experimental effort in creating topologically ordered matter in
quantum simulators. Recent examples include the implementation of toric code
in superconducting qubit systems Googletoric2021 , and dimer models in Rydberg
arrays HarvardQSL1 ; bluvstein2021quantum . Since various kinds of
perturbation are present in any experimental implementation, the precise
Hamiltonian may not be known and may depart significantly from that of the
pristine, idealized models. It is thus important to find a systematic approach
to characterize topological order given a wave function, with minimal
knowledge of the Hamiltonian.

From a modern perspective, one key aspect of topological order is the
existence of an emergent, higher symmetry. To each curve in space, there
exists a set of Wilson line operators (WLOs), which correspond to
adiabatically transporting topologically non-trivial quasiparticles along
$\gamma$. If $\gamma$ is a contractible loop, the corresponding Wilson loop
operators, or closed WLOs, keep a particular ground state invariant
hastings2005quasiadiabatic , while a WLO with open ends creates quasiparticle
excitations near the two endpoints of $\gamma$. In this sense, the closed WLOs
on contractible loops can be thought of as an emergent symmetry of the ground
state.111Closed WLOs on non-contractible loops can be thought of as a
spontaneously broken emergent symmetry gaiotto2014 , because while they keep
the ground state subspace invariant, implying an emergent symmetry, they act
non-trivially on ground states, implying ‘spontaneous symmetry breaking.’  The
symmetry is emergent in general because, aside from certain exactly solvable
models kitaev2003 ; levin2005string , the Hamiltonian need not commute with
these WLOs. In contrast to ordinary symmetries, which are implemented by
operators with support over the entire space, the closed WLOs have support
only on loops; in the case where all topological quasiparticles are Abelian,
the WLOs can be thought of, in modern terminology, as emergent 1-form
symmetries of the system gaiotto2014 . These WLOs should in principle contain
all of the data that characterizes the topological order, however it is not
well-understood how to tease it out in practice.

In this work, we propose a numerical method to systematically search for a
complete set of Wilson loop operators for the case of Abelian topological
orders, using only the ground state of a gapped Hamiltonian defined on a disk-
like region of space. We do this by considering a variational ansatz for
Wilson loop operators in terms of a matrix product operator with support on a
ribbon along $\gamma$, as schematically shown in Fig. 1(a). We then set up a
cost function in terms of the variational parameters of the WLOs. The minima
of the cost function, which we numerically optimize for, gives the WLOs as
diagrammatically shown in Fig. 1(c). The obtained Wilson loop operator
expectation value can reach close to unity after a few hundred iterations
(Fig. 1(d)). We emphasize that our procedure is unbiased and assumes no prior
knowledge of the form of the WLOs.

Once the WLOs are obtained, we show how one can perform further optimization-
based schemes to find operators that can create, move, and annihilate the
anyons. Finally, we show that these operators can be utilized to extract the
modular $S$ and $T$ matrices of an Abelian topological order, which gives a
complete characterization of the topological order. In particular the $S$ and
$T$ matrices encode all of the information about the fractional statistics.

We successfully demonstrate our numerical protocol in models with non-zero
correlation lengths. For example, we show how one can extract the modular $S$
and $T$ matrices from only the ground states of the perturbed toric code and
doubled semion models, with a Zeeman field that is up to half of the critical
value.

To date, several invariants of two-dimensional topologically ordered states
have been shown to be obtainable from the ground state wave function through a
variety of methods. This includes the total quantum dimension measured through
topological entanglement entropy kitaev2006topological ; levin2006detecting ,
the many-body Chern number and Hall conductance dehghani2021extraction ;
cian2021many ; fan2022extracting , various invariants of symmetry-protected
topological (SPT) states Shiozaki2017 ; Shiozaki2018many ; elben2020many , and
the chiral central charge tu2013 ; zaletel2013 ; li2008 ; kim2021modular ;
kim2021chiral . The modular $S$ and $T$ matrices, which encode details of the
fractional statistics of the quasiparticles, can, under certain conditions, be
extracted from the full set of ground states on a torus zhang2012quasiparticle
; zhang2015 ; wen1990naberry  or in the presence of twist defects zhu2020 ,
but not to date from a single ground state on a disk.

We note that Ref. bridgeman2016 ; wahl2020local  also proposed to find WLOs
through an optimization approach, by searching for WLOs that commute with the
Hamiltonian. However generic systems are not expected to have WLOs that
commute with the Hamiltonian; instead as discussed above WLOs only appear as
emergent symmetries that keep the ground state subspace invariant. Our work,
in contrast to Ref.bridgeman2016 ; wahl2020local , uses only the ground state
without requiring knowledge of the Hamiltonian.

The paper is organized as follows. In Sec. II, we review the basic properties
of Wilson line operators and the algebraic theory of anyon. In Sec. III, we
provide the optimization scheme for probing closed Wilson loop operators. In
Sec. IV, we propose the scheme to create, move, annihilate anyons and measure
the topological twist. We present numerical simulations for abelian
topological order models in Sec. V. Finally, we provide an outlook for future
works in Sec. VI.

##  II Wilson loop operators and anyon data

In this section, we briefly review the basic properties of WLOs and the
algebraic theory of anyons. Since our goal is to extract topological
invariants from the bulk of the wave function, in this section we consider a
two dimensional system on an infinite plane.

The anyon theory consists of a collection of algebraic data that characterize
the universal topological properties, namely the fusion and braiding
properties, of the anyonic excitations of a many-body system. The precise
mathematical framework is that of a unitary modular tensor category (UMTC).
For reviews of UMTCs in the context of topological phases of matter, see for
example Ref. kitaev2006 ; Bonderson07b ; wang2008 ; barkeshli2019 ;
kawagoe2020microscopic . For a detailed discussion of how to relate the
algebraic data of the UMTC to the microscopic properties of a quantum many-
body system, see Ref. kawagoe2020microscopic .

The description that we provide below can be made exact in the context of
exactly solvable models, such as the toric code and its generalizations, the
quantum double and Levin-Wen models kitaev2003 ; levin2005string . The effect
of perturbations to these exactly solvable models can also be studied
systematically, using quasi-adiabatic continuation hastings2005quasiadiabatic
. For chiral topological phases, such as fractional quantum Hall states or
fractional Chern insulators, which have no description in terms of an exactly
solvable model, it is expected that the same discussion applies, although it
has not been explicitly studied outside of the context of field theory.

###  II.1 Anyons, Wilson line and loop operators

Since the system has a finite correlation length, we can define states with
quasiparticle excitations that are localized on the scale of the correlation
length. We can then group the quasiparticles into topological equivalence
classes: two quasi-particle excitations are equivalent if and only if there is
a local operator that can convert one into the other. The different
equivalence classes define a finite set of distinct anyon types, sometimes
also referred to as superselection sectors or topological charges,
$\\{I,a,b,c,...\\}$. The set of anyons contains the identity sector $I$, which
corresponds to excitations that can be created by local operators.

Since the anyonic excitations can be localized to within a correlation length
of a particular point in space, we can consider a state with anyon type $a$ at
position $x$, and denote it as $|a_{x}\rangle$. In defining $|a_{x}\rangle$,
we assume that far away from $x$, on the scale of the correlation length
$\xi$, $|a_{x}\rangle$ locally looks like the ground state. Furthermore, we
assume any other non-trivial topological charges are infinitely far away and
do not include them in labeling the state $|a_{x}\rangle$. Note that since $a$
refers to an equivalence class of excitations, there are many states that can
be labeled as $|a_{x}\rangle$, so our choice is not unique.

By construction, the expectation value of any local observable $O(x^{\prime})$
satisfies $\langle a_{x}|O(x^{\prime})|a_{x}\rangle=\langle
GS|O(x^{\prime})|GS\rangle$, as long as $x$ and $x^{\prime}$ are far away from
each other, $|x-x^{\prime}|\gg\xi$, where $|GS\rangle$ is the ground state of
the system and $\xi$ is the correlation length. The above equality holds up to
$O(e^{-|x-x^{\prime}|/\xi})$ corrections. Physically, this corresponds to the
fact that the state has short-range correlations, so that a disturbance in the
vicinity of $x$ has exponentially decreasing effects in the ground state
beyond a correlation length.

Figure 2:  (a) When the movement operator $W_{a}(\gamma)$ is applied to the
state $|a_{x}\rangle$, the anyon excitation is moved to $x^{\prime}$. (b) When
the splitting operator $W_{c\rightarrow a,b}(\gamma)$ is applied to the state
$|c_{x}\rangle$, the anyon splits into two anyons $a$ and $b$ at position $x$
and $x^{\prime}$, respectively.

The anyons can be moved from one place to another by applying an operator
along an arbitrary path. In particular, we define a Wilson line operator
$W_{a}(\gamma)$ along a path $\gamma$ starting at $x_{1}$ and ending at
$x_{2}$. $W_{a}(\gamma)$ moves the anyon excitation from position $x_{1}$ to
$x_{2}$ along the path $\gamma$ as shown in Fig. 2(a). Such operators are
referred to as “movement operators” in kawagoe2020microscopic . Specifically,

$\displaystyle|a_{x_{2}}\rangle=W_{a}(\gamma)|a_{x_{1}}\rangle.$ (1)

In general, $W_{a}(\gamma)$ has support, up to exponentially small
corrections, on a ribbon of thickness on the scale of $\xi$, centered on
$\gamma$. 222Specifically, $W_{a}(\gamma)$ can be approximated by a ribbon
operator with finite thickness $t$ hastings2005quasiadiabatic . The error of
the approximate Wilson line operator $W$ is of order
$\epsilon=|\langle\psi|W-W_{\rm exact}|\psi\rangle|\sim O(N_{s}e^{-t/\xi})$,
where $N_{s}$ is the number of sites in the support of $W$, and $W_{\rm
exact}$ is the exact, and presumably non-local WLO, for the ground state wave
function $|\psi\rangle$. The precise choice of the operator $W_{a}(\gamma)$ is
not unique in general, and a precise definition is also non-universal and
depends on the microscopic details of the system. Nevertheless,
$W_{a}(\gamma)$ encodes certain universal topological data that we wish to
extract.

Note that $W_{a}(\gamma)$ in general need not be a unitary or even invertible
operator, although in many simple examples, particularly when $a$ is an
Abelian anyon, $W_{a}(\gamma)$ can be chosen to be unitary. Moreover, when $a$
is Abelian, we can take

$\displaystyle W_{a}^{\dagger}(\gamma)=W_{\bar{a}}(\gamma)=W_{a}(-\gamma),$
(2)

where $-\gamma$ refers to the path $\gamma$ traversed in the opposite
direction and $\bar{a}$ is the anti-particle of $a$. That is,
$W_{a}^{\dagger}(\gamma)$ effectively takes $\bar{a}$ along $\gamma$ or,
equivalently, takes $a$ along the path $-\gamma$.

We can also define loop operators by picking $\gamma$ in $W_{a}(\gamma)$ to be
a closed loop. Physically this can be understood as creating $a$ and its dual
$\bar{a}$ out of the ground state, moving one around the loop $\gamma$, and
reannihilating. If $\gamma$ is a contractible loop in the space, such
operators should keep the ground state invariant. Therefore, for each loop
$\gamma$ we have a loop operator $W_{a}(\gamma)$, which keeps the ground state
invariant:

$\displaystyle W_{a}(\gamma)|GS\rangle=d_{a}|GS\rangle.$ (3)

Here $d_{a}\geq 1$ is referred to as the quantum dimension of the anyon $a$,
and is part of the universal data of the UMTC. We have chosen a convention
where $d_{a}\geq 1$ appears on the RHS; we could in principle absorb $d_{a}$
into the definition of $W_{a}(\gamma)$. The choice above allows us to make
contact with the fusion algebra of the UMTC description.

As with the line operators, the Wilson loop operators $W_{a}(\gamma)$ are in
general not unitary operators, unless $a$ is an Abelian anyon, in which case
we also have $d_{a}=1$.

Eq. 3 makes explicit that the ground state of a topologically ordered state
has emergent symmetries, as there are loop operators that keep the ground
state invariant. Importantly, the operators are supported, up to exponentially
small corrections, on a codimension-1 region, and therefore do not correspond
to ordinary global symmetries, which have support over the entire space. When
the anyons are Abelian, the Wilson loop operators form a group structure and
are referred to as 1-form symmetries gaiotto2014 ; more generally they are
referred to as categorical or non-invertible symmetries.

###  II.2 Fusion rules and splitting operators

The anyons define a fusion algebra

$\displaystyle a\times b=\sum_{c}N_{ab}^{c}c,$ (4)

where the fusion multiplicities $N_{ab}^{c}$ are non-negative integers, which
indicate the number of different ways the anyons $a$ and $b$ can be fused to
produce the anyon type $c$. Each anyon type $a$ has a unique anti-particle
$\bar{a}$, where $\bar{a}\in\mathcal{C}$ is such that $N_{a{\bar{a}}}^{I}\neq
0$. Note that we can define a fusion matrix $N_{a}$, with entries
$(N_{a})_{bc}=N_{ab}^{c}$; the quantum dimension $d_{a}$ is then the largest
eigenvalue of $N_{a}$.

An anyon $a$ is Abelian if and only if it gives a unique fusion outcome upon
fusing with another anyon $b$. That is, given $b$, $N_{ab}^{c}=1$ for a unique
$c$ and $N_{ab}^{c}=0$ otherwise.

The fusion rule leads to the following relation for the Wilson loop operator:

$W_{a}(\gamma)W_{b}(\gamma)|GS\rangle=\sum_{c}N_{ab}^{c}W_{c}(\gamma)|GS\rangle,$
(5)

whenever $\gamma$ is a loop. Note that with the conventions chosen above, this
implies $d_{a}d_{b}=\sum_{c}N_{ab}^{c}d_{c}$.

In addition to movement operators and loop operators, we can define splitting
operators. For simplicity, here we only introduce the spliting operator in the
case $N_{ab}^{c}\leq 1$; the generalization can be found in Ref.
kawagoe2020microscopic . Suppose that $c$ is contained in the fusion outcome
of $a$ and $b$, that is, $N_{ab}^{c}=1$. We can define a splitting operator,

$\displaystyle W_{c\rightarrow
a,b}(\gamma)|c_{x_{1}}\rangle=|a_{x1},b_{x_{2}}\rangle,$ (6)

where $|a_{x1},b_{x_{2}}\rangle$ denotes a state with two excitations: anyon
$a$ at position $x_{2}$ and anyon $b$ at position $x_{1}$ and
$|x_{2}-x_{1}|\gg\xi$ as shown in Fig. 2(b).

One can also create the anyon $a$ and its antiparticle $\bar{a}$ by applying

$W_{I\rightarrow\bar{a},a}(\gamma)|I_{x_{1}}\rangle=|\bar{a}_{x_{1}},a_{x_{2}}\rangle,$
(7)

where $|I\rangle$ is a state in the identity superselection sector.

Observe that if we start with a loop operator $W_{a}(\gamma)$, and we project
part of the operator along some segment of $\gamma$ to the identity, then we
obtain a cut operator that effective creates an anyon and its anti-particle
out of the vacuum. Therefore we can obtain a choice of $W_{I\rightarrow
a,\bar{a}}(\gamma_{\text{cut}})$ by starting with $W_{a}(\gamma)$ for a loop
$\gamma$ and implementing the above cutting procedure.

###  II.3 Modular $S$ matrix and twist product

A large portion of the universal data of a topological phase of matter is
encoded in the modular $S$ and $T$ matrices. In fact for almost all
topological phases of interest in physics, the $S$ and $T$ matrices provide a
complete set of invariants.

The modular $S$ matrix contains information about the mutual braiding
statistics between far separated anyon excitations, and also completely
defines the fusion coefficients $N_{ab}^{c}$.

In particular, $S_{ab}$ is the quantum mechanical amplitude of the process
where a particle of type $a$ and another particle of type $b$ are created and
separated, the particle $a$ is moved around the particle $b$, and then the
particle-anti-particle pairs are annihilated.

Given a set of closed Wilson loop operators $W_{a}$, where $a\in\mathcal{C}$,
we can define a matrix $\tilde{S}_{ab}$, as

$\displaystyle\tilde{S}_{a,b}=\frac{\langle GS|W_{a}\infty
W_{b}|GS\rangle}{\langle GS|W_{a}W_{b}|GS\rangle},$ (8)

where $\infty$ is the twist product (see e.g. Ref. haah2016invariant ) and
shown in Fig. 3. For arbitrary WLOs $P$ and $Q$ defined in regions $A$ and $B$
as shown in Fig. 3(a), $P=\sum_{k}P^{A}_{k}\otimes P^{B}_{k}$,
$Q=\sum_{l}Q^{A}_{l}\otimes Q^{B}_{l}$, the twist product $P\infty Q$ is
defined as

$\displaystyle P\infty Q=\sum_{kl}P_{k}^{A}Q_{l}^{A}\otimes
Q_{l}^{B}P_{k}^{B},$ (9)

where the product order is reversed in the region $B$.

Note that to define the above twist product, we need each operator $W_{a}$ to
have support strictly on a ribbon of finite thickness.

Figure 3: (a) Mutual-braiding statistics. To measure the mutual braiding
statistics, we calculate twist product for two Wilson loop operators $W_{a}$
and $W_{b}$. The supports of the two Wilson loop operators are divided into
two regions. In region A, the order of the product is $W_{a}W_{b}$ and in the
region B, the order of the product is reversed. (b) The matrix product
operator that measures the self-braiding statistic $T_{a}$.

We expect that $\tilde{S}_{ab}$ is related to $S_{ab}$ as

$\displaystyle\tilde{S}_{ab}=\frac{d_{a}d_{b}}{\mathcal{D}}S_{ab},$ (10)

where $\mathcal{D}=\sqrt{\sum_{a}d_{a}^{2}}$ is the total quantum dimension.
In this paper we only work with Abelian anyons, in which case $d_{a}=1$ and
$\mathcal{D}^{2}$ is the total number of distinct anyon types.

For Abelian anyons, the braiding phase between anyon $i$ and $j$ can be
measured from the phase of the twist product

$\displaystyle\phi_{i,j}={\rm arg}[\tilde{S}_{i,j}].$ (11)

###  II.4 Modular $T$ matrix

The modular $T$-matrix is a diagonal matrix,

$\displaystyle T_{ab}=\theta_{a}\delta_{ab}$ (12)

where $\theta_{a}$ is the topological twist of the anyon $a$. Due to the spin-
statistics theorem, $\theta_{a}$ also corresponds to the exchange statistics
of $a$. In order to exchange a pair of identical anyons, we first create two
anyon and anti-anyon pairs from the ground state, and we then move the two
identical anyons and exchange them. Finally we fuse the anyon and anti-anyon
and return the the ground state as shown in Fig. 3(b). If we normalize the
process properly, the net effect of this procedure gives the exchange
statistics of the anyons.

One can create and exchange anyons and measure the exchange statistics using
the Wilson loop operators. The detailed implementation of the extraction of
exchange statistics using the Wilson loop operators is given in Sec. IV.5.

##  III Optimization Scheme

Figure 4: Illustration of the tensor contraction for calculating
$\langle\psi|W|\psi\rangle$. The gray tensors represents corner transfer
matrix (CTM). The blue tensors are the contraction of the bra and ket of the
PEPS wave function. The two are connected by contraction of the physical bond.
The red tensors are the contraction of the PEPS wave function and the WLOs.

In order to study the WLOs in the bulk of a ground state wave function, we
propose a numerical scheme to extract contractible closed WLOs $W_{a}(\gamma)$
as defined in Eq. (3). We parametrize the WLOs by an ansatz based on matrix
product operators (MPOs) verstraete2008matrix ; bridgeman2016 ;
bultinck2017anyons . The ansatz is defined by two parameters: $(R,\chi)$,
where $R$ is a set of sites corresponding to the support of the WLO and $\chi$
is the bond dimension, as shown in Fig. 1. For a certain class of analytically
solvable topologically ordered states, e.g. the Levin-Wen model and the Kitaev
quantum double model, it is known that the WLOs can be efficiently expressed
by MPOs with a support region $R$ with small thickness $t$ and bond dimension
$\chi$ kitaev2003 ; levin2005string ; bultinck2017anyons .

To extract the closed WLOs using the MPO ansatz, we numerically optimize the
MPO ansatz in order to obtain a $W_{a}(\gamma)$ which approximately satisfies
Eq. (3). In order to efficiently optimize a Wilson loop, we note that for
Abelian topological orders, an operator $W$ and a wave function $|\psi\rangle$
satisfies Eq. (3) if and only if

$\displaystyle\langle\psi|W|\psi\rangle=1,$ (13)
$\displaystyle\langle\psi|W^{\dagger}W|\psi\rangle=1.$ (14)

The forward proof is trivial. For the backward proof, we assume that $W$ and
$|\psi\rangle$ satisfy Eqs. (13) and (14) and without loss of generality,
$W|\psi\rangle=a|\psi\rangle+b|\psi_{\perp}\rangle$, where $a$ and $b$ are
complex numbers and $\langle\psi|\psi_{\perp}\rangle=0$. Solving Eqs. (13) and
(14) we can obtain $a=1$ and $b=0$. Note that Eq. (11) implies that $a=1$ and
Eq. (12) ensures that $W|\psi\rangle$ is normalized to $1$, which then
requires $b=0$.

Therefore, we can define the cost function for a wave function $|\psi\rangle$
as

$\displaystyle\mathcal{C}(W)=[\langle\psi|W|\psi\rangle-1]^{2}+[\langle\psi|W^{\dagger}W|\psi\rangle-1]^{2}.$
(15)

It reaches a global minimum $\mathcal{C}(W)=0$ only when $W$ is an exact
Wilson loop operator.

To variationally optimize the WLO, we start by initializing a random MPO with
fixed $(R,\chi)$. Each tensor in the MPO is initialized randomly and
independently from each other. For a translationally invariant system, one may
naively expect that a translation symmetric closed WLO is a better ansatz.
However, we found that the translation symmetric ansatz tends to be unstable
numerically, leading to diverging or vanishing gradients in the optimization
procedure.

After the initialization, we minimize the cost function defined in Eq. (15)
through gradient based optimization. In this work, we apply the Adam algorithm
kingma2014adam  to minimize the cost function. In this work, we fix the
hypermeters of the Adam algorithm as $\beta_{1}=0.9,\beta_{2}=0.999$ and
learning rate $10^{-3}$. We iterate the optimization procedure until the cost
function converges, which typically takes a few hundred to a few thousand
iterations.

We repeat the initialization and minimization $N_{\rm sample}$ times to obtain
$N_{\rm sample}$ optimized WLOs $W_{k}$, where $1\leq k\leq N_{\rm sample}$.
Throughout the manuscript, we fix $N_{\rm sample}=20$. We then measure the
braiding phases and topological twists to classify the WLOs through the
equivalence relation described as follows : we compute the mutual-braiding
phases between $W_{i}$ and $W_{j}$, $\phi_{i,j}$, where $1\leq i,j\leq N_{\rm
sample}$ and topological twist for $W_{i}$ $\theta_{i}$, where $1\leq i\leq
N_{\rm sample}$. We say the two WLOs $W_{i}$ and $W_{j}$ are equivalent when
$\phi_{i,k}=\phi_{j,k}$ for $1\leq k\leq N_{\rm sample}$, and
$\theta_{i}=\theta_{j}$. After grouping the WLOs into equivalence classes, we
randomly pick one representative WLO from each equivalent class.

This equivalence relation assumes that if two WLOs have the same braiding
phase with the rest of the WLOs and identical topological twists, the two WLOs
are equivalent. We note that this condition only holds when we obtain a
complete set of WLOs in our optimization procedure. Missing one could result
in a false classification. However, one can verify whether a complete set of
WLOs is found by checking if the resulting $S$ matrix is a unitary matrix. If
one or more WLOs are missing, we can vary the hyper-parameters such as
increasing the thickness of the WLOs or the bond dimension $\chi$.

The bottleneck of the numerical optimization is in the tensor contraction when
calculating the expectation values $\langle\psi|W|\psi\rangle$ and
$\langle\psi|W^{\dagger}W|\psi\rangle$. Here, we briefly describe the tensor
contraction scheme and its computational time complexity. Given a ground state
wave function $|\psi\rangle$, represented by an infinite projective entangled
pair state (iPEPS) orus2009simulation ; liao2019differentiable ;
crone2020detecting  as shown in Fig. 4 , we use the following procedure to
evaluate the expectation value $\langle\psi|W|\psi\rangle$. Consider a closed
WLO that takes the form of a rectangular loop with side lengths $L_{x}$ and
$L_{y}$ as shown in Fig. 4, for a system on a square lattice. We first
contract all tensors at $x=0$ to form a tensor $M$ with $L_{y}$ bonds as shown
in Fig. 4. We then contract the tensors at $x=1$ with $M$ one by one from
$y=0$ to $y=L_{y}+1$. We then repeat this procedure for all $x\leq L_{x}+1$.
In this contraction procedure, the computational cost for contracting the
tensors scales linearly with the number of sites along the $x$ direction
($L_{x}$) and exponentially with number of sites along the $y$ direction
($L_{y}$). Therefore, the total computational cost is bounded by
$O(L_{x}L_{y}\chi^{L_{y}+5})$. Moreover, since the thickness $t$ and the size
of the hole need to be much larger than the correlation length, we have that
$L_{x},L_{y}>3\xi$. Thus the total computational cost for calculating the
expectation value scales up as $O(\xi^{2}\chi^{\alpha\xi})$ for some constant
$\alpha\gg 1$. Therefore, this method is particularly suitable for models with
small correlation length. Importantly, the complexity of the computation
scales with the correlation length, not the total system size.

##  IV Manipulation of anyons

Once we have obtained the closed WLOs, we can extract many non-trivial
properties of the anyons. In particular, we can further obtain operators that
create, move, and annihilate anyons.

In the following, we discuss how to manipulate anyons by starting from the
closed WLOs. For simplicity in this section, we assume that the thickness of
WLOs is $t=1$ . The idea can be easily generalized to $t>1$ as is the case for
our numerical results in Section V.

###  IV.1 Removing and adding a virtual bond

Before we proceed to the manipulation of anyons, we first define two basic
operations of a tensor $A$, which is that of removing and adding a virtual
bond. These two operations will be extensively used throughout this section.

Removing a virtual bond is useful for cutting open closed WLOs. To remove a
virtual bond in a tensor $A_{\alpha,\beta,\dots,\gamma}$, we define an edge
tensor $E_{\alpha}$ that describes the boundary condition and contract the
edge tensor $E_{\alpha}$ with $A_{\alpha,\beta,\dots,\gamma}$; the resulting
tensor is

$\displaystyle
A^{\prime}_{\beta,\dots,\gamma}=\sum_{\alpha}E_{\alpha}A_{\alpha,\beta,\dots,\gamma},$
(16)

and it has rank $n-1$.

Adding a virtual bond is useful when extending an open WLO or connecting two
open WLOs. To add a trivial virtual bond to an arbitrary tensor
$A_{\alpha,\beta,\dots,\gamma}$ with rank $n$, we define a new tensor
$\tilde{A}$ with rank $n+1$ as

$\displaystyle\tilde{A}_{\alpha^{\prime},\alpha,\beta,\dots,\gamma}=A_{\alpha,\beta,\dots,\gamma},$
(17)

for all $1\leq\alpha^{\prime}\leq\chi$, where $\chi$ is the bond dimension of
the new virtual bond.

Using the above two basic operations on tensors, we can create, move and
annihilate anyons, as described in the following sections.

###  IV.2 Creation of anyon and anti-anyon pairs

To create an anyon and anti-anyon pair, we simply discard a segment of the
closed WLO and apply an arbitrary edge tensor to terminate its boundary as
shown in Fig. 5. Specifically, we start with a closed WLO obtained using the
procedure described in Sec. III for a give ground state wave function
$|\psi\rangle$, which is of the form

$\displaystyle
W=\sum_{\\{\alpha\\}}A^{1}_{\alpha_{0},\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{L-1}_{\alpha_{L-2},\alpha_{L-1}}A^{L}_{\alpha_{L-1},\alpha_{0}},$ (18)

where $L$ is the length of the closed WLO.

We then keep the tensors from sites $1,\cdots,l$, as shown in step I of Fig.
5, so that the new MPO is of the form

$\displaystyle\tilde{W}_{\beta,\beta^{\prime}}=\sum_{\\{\alpha\\}}A^{1}_{\beta,\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{l-1}_{\alpha_{l-2},\alpha_{l-1}}A^{l}_{\alpha_{l-1},\beta^{\prime}}.$ (19)

Figure 5: Creation of an anyon and anti-anyon pair. Step I : We truncate a
closed WLO by discarding a segment of it. Step II : We contract two two
boundary tensor with edge tensor to terminate the open WLO.

We then remove the virtual indices $\beta$ and $\beta^{\prime}$
$\tilde{W}_{\beta,\beta^{\prime}}$ by applying the edge tensor $E$. The edge
tensors can be chosen arbitrarily, since different choices are related to each
other by a local operator. Throughout out this article, we use

$\displaystyle E_{\gamma}=\begin{cases}1,&\text{if}~{}\gamma=1\\\
0,&\text{otherwise}\end{cases}.$ (20)

After removing the edge virtual bond shown in step II in Fig. 5, the open WLO
is of the form

$\displaystyle W^{\rm
open}=\sum_{\beta,\beta^{\prime}}\tilde{W}_{\beta,\beta^{\prime}}E_{\beta}E_{\beta^{\prime}}.$
(21)

When another WLO passes through the open WLO, the system acquires an anyonic
braiding phase as long as the crossing point is away from the boundary of the
open WLO by an $O(\xi)$.

Finally, in order to preserve the norm of the wave function, we normalize the
open WLO defined above by a factor $\sqrt{\langle\psi|W^{\rm
open\dagger}W^{\rm open}|\psi\rangle}$.

###  IV.3 Moving anyons

Let us imagine we have a WLO $W^{\text{open}}_{x_{1},x_{2}}$ that creates an
anyon $a$ at one endpoint $x_{1}$ and its conjugate $\overline{a}$ at the
other endpoint $x_{2}$. We can use this to construct a different operator
$W^{\text{open}}_{x_{1}^{\prime},x_{2}}$, effectively moving $a$ from $x_{1}$
to $x_{1}^{\prime}$. We can do this as follows.

We start with an open WLO as shown in step I of Fig. 6. The open WLO with
length $L$ is of the form

$\displaystyle W^{\rm
open}=\sum_{\\{\alpha\\}}A^{1}_{\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{L-1}_{\alpha_{L-2},\alpha_{L-1}}A^{L}_{\alpha_{L-1}}.$ (22)

The open WLO can be obtained using the procedure described in Sec. IV.2.

We initialize another random MPO with length $L^{\prime}$ that is a loop
complement of the WLO $W^{\rm open}$ as shown in step II of Fig. 6 and is of
the form

$\displaystyle W^{\rm
open}_{C}=\sum_{\\{\beta\\}}C^{1}_{\beta_{1}}C^{2}_{\beta_{1},\beta_{2}}\cdots
C^{L^{\prime}-1}_{\beta_{L^{\prime}-2},\beta_{L^{\prime}-1}}C^{L^{\prime}}_{\beta_{L^{\prime}-1}},$
(23)

Figure 6: Moving anyons. Step I : We start with an open WLO that can be
generated using the method described in Sec. IV.2. Step II : We randomly
initialize an open MPO with virtual bond at the boundary tensors. The
initialized open MPO and the open WLO form a closed MPO. Step III: We vary the
closed MPO to minimize the cost function defined in Eq. (15). Step IV and V
Finally, we cut the closed WLO to create two excitations that can be located
at different sites from the anyons in Step I.

To connect $W^{\rm open}$ and $W^{\rm open}_{P}$ at the boundary, we add a
trivial virtual bond on each boundary tensor. By adding trivial virtual bond
on all the boundary tensors $A^{1},A^{L},C^{1},C^{L^{\prime}}$, we have

$\displaystyle\tilde{W}^{\rm open}_{\gamma_{0},\gamma_{L}}$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha\\}}\tilde{A}^{1}_{\gamma_{0},\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{L-1}_{\alpha_{L-2},\alpha_{L-1}}\tilde{A}^{L}_{\alpha_{L-1},\gamma_{L}},$
$\displaystyle\tilde{W}^{\rm open}_{C,\mu_{0},\mu_{L^{\prime}}}$
$\displaystyle=$
$\displaystyle\sum_{\\{\beta\\}}\tilde{C}^{1}_{\mu_{0},\beta_{1}}C^{2}_{\beta_{1},\beta_{2}}\cdots
C^{L^{\prime}-1}_{\beta_{L^{\prime}-2},\beta_{L^{\prime}-1}}\tilde{C}^{L^{\prime}}_{\beta_{L^{\prime}-1},\mu_{L^{\prime}}}.$

We can therefore have a closed WLO of the form

$\displaystyle
W=\sum_{\gamma_{0},\gamma_{L},\mu_{0},\mu_{L^{\prime}}}\tilde{W}^{\rm
open}_{\gamma_{0},\gamma_{L}}\tilde{W}^{\rm
open}_{C,\mu_{0},\mu_{L^{\prime}}}\delta_{\gamma_{0},\mu_{L^{\prime}}}\delta_{\gamma_{L},\mu_{0}},$
(25)

where $\delta_{i,j}$ is a Kronecker delta function.

In step III of Fig. 6, we minimize the cost function of Eq. (15) for the
closed WLO defined above. We note that in addition to vary the tensors in the
$\tilde{W}^{\rm open}_{C,\mu_{0},\mu_{L^{\prime}}}$, we also have to vary
tensors at the boundary of $\tilde{W}^{\rm open}_{\gamma_{0},\gamma_{L}}$ i.e.
$\tilde{A}^{1}$, $\tilde{A}^{L}$,…etc. in order to eliminate the anyon
excitation at the boundary.

We define a length parameter $b$, which is of order $O(\xi)$. In the
optimization procedure, we vary all the tensors $C$s and $\tilde{C}$s in
$W_{C}^{\rm open}$ and the boundary tensor of $W^{\rm open}$, $\tilde{A}^{1}$,
$A^{2}$, $\dots$, $A^{b}$ and $A^{L-b+1}$, $\cdots$ $A^{L-1}$,
$\tilde{A}^{L}$, while fixing the tensors $A^{i}$ for $b+1\leq i\leq L-b+1$.

Once we find the optimal solution that approximately satisfies Eq. (3), we can
cut the closed WLO to create two ends as described in Sec. IV.2 and
effectively move the anyon as shown in step IV and V in Fig. 6.

###  IV.4 Annihilation of anyon and anti-anyon pairs

In this section, we describe a procedure to fuse an anyon and anti-anyon pair
to identity. Given two open WLOs $W_{a}(\gamma)$ and $W_{a}(\gamma^{\prime})$,
we show how to join them into a single WLO, as shown in Fig. 7, by effectively
bringing together two endpoints of $\gamma$ and $\gamma^{\prime}$ and
annihilating the anyons.

We consider two open WLOs of the form

$\displaystyle W^{\rm open}_{1}$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha\\}}A^{1}_{\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{L_{1}-1}_{\alpha_{L_{1}-2},\alpha_{L_{1}-1}}A^{L_{1}}_{\alpha_{L_{1}-1}},$
$\displaystyle W^{\rm open}_{2}$ $\displaystyle=$
$\displaystyle\sum_{\\{\beta\\}}B^{1}_{\beta_{1}}B^{2}_{\beta_{1},\beta_{2}}\cdots
B^{L_{2}-1}_{\beta_{L_{2}-2},\beta_{L_{2}-1}}B^{L_{2}}_{\beta_{L_{2}-1}}.$
(26)

For simplicity we assume that the tensors $A^{L_{1}-1}$ and $B^{1}$ are
located at nearest-neighbor sites as shown in step I of Fig. 6. If this is not
the case, we can move the end of $W^{\rm open}_{2}$ using the procedure
described in Sec. IV.3.

Figure 7: Annihilation of anyon and anti-anyon pairs Step I: We start with two
open WLOs The boundary of the two open WLOs are at nearest neighbor sites.
Step II: We randomly initialize an open MPO with virtual bonds at the boundary
tensors. The two open WLOs and the MPO form a closed MPO. Step III: We vary
the closed MPO to minimize the cost function defined in Eq. (15). Step IV and
V : Finally, we cut the closed WLO and keep the tensors on supports of the two
open WLOs in Step I.

We then initialize a random MPO with length $L^{\prime}$ that is a loop
complement of $W^{\rm open}_{1}$ and $W^{\rm open}_{2}$, which takes the form

$\displaystyle W^{\rm
open}_{C}=\sum_{\\{\gamma\\}}C^{1}_{\gamma_{1}}C^{2}_{\gamma_{1},\gamma_{2}}\cdots
C^{L^{\prime}-1}_{\gamma_{L^{\prime}-2},\gamma_{L^{\prime}-1}}C^{L^{\prime}}_{\gamma_{L^{\prime}-1}}.$
(27)

We can connect the three open WLOs by adding a trivial virtual bond on each
boundary tensor as in step III of Fig. 7. After adding a virtual bond, the
WLOs become

$\displaystyle\tilde{W}^{\rm open}_{1,\mu_{0},\mu_{L}}$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha\\}}\tilde{A}^{1}_{\mu_{0},\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{L_{1}-1}_{\alpha_{L_{1}-2},\alpha_{L_{1}-1}}\tilde{A}^{L_{1}}_{\alpha_{L_{1}-1},\mu_{L_{1}}},$
$\displaystyle\tilde{W}^{\rm open}_{2,\nu_{0},\nu_{L}}$ $\displaystyle=$
$\displaystyle\sum_{\\{\beta\\}}\tilde{B}^{1}_{\nu_{0},\beta_{1}}B^{2}_{\beta_{1},\beta_{2}}\cdots
B^{L_{2}-1}_{\beta_{L_{2}-2},\beta_{L_{2}-1}}\tilde{B}^{L_{2}}_{\beta_{L_{2}-1},\nu_{L_{2}}},$
$\displaystyle\tilde{W}^{\rm open}_{C,\sigma_{0},\sigma_{L}}$ $\displaystyle=$
$\displaystyle\sum_{\\{\gamma\\}}\tilde{C}^{1}_{\sigma_{0},\gamma_{1}}C^{2}_{\gamma_{1},\gamma_{2}}\cdots
C^{L^{\prime}-1}_{\gamma_{L^{\prime}-2},\gamma_{L^{\prime}-1}}\tilde{C}^{L^{\prime}}_{\gamma_{L^{\prime}-1},\sigma_{L^{\prime}}}.$
(28)

We can therefore connect the three WLO and have a closed MPO of the form

$\displaystyle
W=\sum_{\mu_{0},\mu_{L_{1}},\nu_{0},\nu_{L_{2}},\sigma_{0},\sigma_{L^{\prime}}}\tilde{W}^{\rm
open}_{1,\mu_{0},\mu_{L_{1}}}\tilde{W}^{\rm
open}_{2,\nu_{0},\nu_{L_{2}}}\tilde{W}^{\rm open}_{C,\sigma_{0},\sigma_{L}}$
$\displaystyle\times\delta_{\mu_{L_{1}},\nu_{0}}\delta_{\nu_{L_{2}},\sigma_{0}}\delta_{\sigma_{L^{\prime}},\mu_{0}}.$
(29)

To minimize the cost in Eq. (15) for the WLO above, we vary the MPO above
except for the tensors $A^{i}_{\alpha_{i-1},\alpha_{i}}$,
$B^{i}_{\beta_{i-1},\beta_{i}}$, where $b\leq i\leq L_{1}-b$ . After obtaining
a closed WLO as shown in step III of Fig. 7, we cut the WLO at the support of
$W_{C}^{\rm open}$ using the procedure described in Sec. IV.2. The resulting
open WLO is the fusion of $W_{1}^{\rm open}$ and $W_{2}^{\rm open}$.

###  IV.5 Topological twist (exchange statistics)

Figure 8: Topological twist. Step I. We cut a closed WLO and keep two open
WLO. Step II. We flip the right open WLO to reverse the direction of the anyon
and anti-anyon. The boundary sites are denoted as $i$, $j$, $m$ and $n$. Step
III. We move the anyons at $i$ and $n$ to $m$ and $j$ respectively. Step IV.
We connect the WLO to form a closed loop. Step V.  Finally, we contract the
physical indices in the self-intersecting region.  Figure 9: (a) Procedure for
creating the support of WLO for process $P_{2}$ defined in Sec. IV.5 from
process $P_{1}$. We start with the support of process $P_{1}$. We rotate the
right half of the support and stack it on top of the left half of the support.
This WLO is equivalent to a trivial loop. (b) Transformation from the WLO for
process $P_{1}$ to process $P_{2}$, using the procedure described in (a).

Here we present the scheme that we use to extract the topological twists of
the anyons. To calculate the topological twists, we calculate the ratio of the
amplitude for the following two processes. In the first process ($P_{1}$), we
create an anyon and anti-anyon pair from the ground state wave function. We
then exchange them, and finally, we annihilate the pair of anyon and anti-
anyon. In the second process ($P_{2}$), we create and annihilate the pair
directly without any exchange

For the first process, we start with a closed WLO with length $L_{o}$ of the
form

$\displaystyle
W=\sum_{\\{\alpha\\}}A^{1}_{\alpha_{0},\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{L_{o}-1}_{\alpha_{L_{o}-2},\alpha_{L_{o}-1}}A^{L_{o}}_{\alpha_{L_{o}-1},\alpha_{0}}.$
(30)

We then cut the closed WLO and keep two segments with length $l$ as shown in
step I of Fig. 8. These segments represent a creation of an anyon and anti-
anyon. The end of the two segments of the WLOs are away from each other by
distance $b$ which is an integer larger than $O(\xi)$.

The WLOs are of the form

$\displaystyle W^{\rm open}_{1}$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha\\}}A^{1}_{\alpha_{1}}A^{2}_{\alpha_{1},\alpha_{2}}\cdots
A^{l-1}_{\alpha_{l-2},\alpha_{l-1}}A^{l}_{\alpha_{l-1}},$ $\displaystyle
W^{\rm open}_{2}$ $\displaystyle=$
$\displaystyle\sum_{\\{\alpha\\}}A^{b+1}_{\alpha_{b+1}}A^{b+2}_{\alpha_{b+1},\alpha_{b+2}}\cdots
A^{b+l-1}_{\alpha_{b+l-2},\alpha_{b+l-1}}A^{b+l}_{\alpha_{b+l-1}}.$

In Step II, we flip the direction of the anyon transport for the WLO $W^{\rm
open}_{2}$ by applying Hermitian conjugation.

After cutting and flipping the WLO, we move the anyon located at the open ends
$i$ and $n$ to sites $m$ and $j$ respectively, as shown in step III of Fig. 8.
We then connect the WLOs by annihilating anyon and anti-anyon pairs in step IV
of Fig. 8. The WLO becomes a self-intersecting closed loop of the form

$\displaystyle W^{\rm
twist}_{P_{1}}=\sum_{\\{\beta\\}}B^{1}_{\beta_{0},\beta_{1}}B^{2}_{\beta_{1},\beta_{2}}\cdots
B^{L-1}_{\beta_{L-2},\beta_{L-1}}B^{L}_{\beta_{L},\beta_{0}},$ (32)

where $L$ is the path length of the closed WLO $W^{\rm twist}_{P_{1}}$. The
labels for the supports are shown in Fig. 8. Finally, we contract the physical
indices in the self-intersecting region as shown in step V of Fig. 8.

In step I of process $P_{1}$, we cut and discard two segments of a closed WLO
as described in Eq. (LABEL:eq_twist_two_segs). This step explicitly breaks the
gauge symmetry of a matrix product operator, which introduces non-universal
complex phases to $W^{\rm open}_{1}$ and $W^{\rm open}_{2}$. The non-universal
complex phase depends on the details of the implementation. In order to cancel
the non-universal complex phase, we calculate the amplitude of the second
process $P_{2}$ using the same tensors in Eq. (32), but we contract the
tensors without exchanging anyons. We conjecture and numerically verify that
the process $P_{2}$ has the same non-universal complex phase as the process
$P_{1}$. The calculation of process $P_{2}$ can be achieved by rotating the
right half of the support of $W^{\rm twist}_{P_{1}}$ around the center of
self-intersecting region and stacking it on top of the left half of the
support as shown in Fig. 9(a). Instead of creating, exchanging and
annihilating the anyon and anti-anyon pair, this process creates an anyon and
anti-anyon pair and the anyon travels around the left half of the support
twice. And subsequently it is annihilated with the anti-anyon. We note that in
order to perform the rotation, presumably the ground state wave function must
have rotation symmetry around the center of the self-intersecting region and
translation symmetry. We have not systematically studied how the above
procedure would need to be modified if the translation and rotation symmetries
of the system are broken.

The WLO after the rotation is of the following form

$\displaystyle W^{\rm
twist}_{P2}=\sum_{\\{\beta\\}}C_{\beta_{1},\beta_{\frac{L}{2}-1},\beta_{\frac{L}{2}+1},\beta_{L-1}}\prod_{i=2}^{\frac{L}{2}}D^{i}_{\beta_{i-1},\beta_{L-i+1},\beta_{i},\beta_{L-i}}$

where

$\displaystyle
C_{\beta_{1},\beta_{\frac{L}{2}-1},\beta_{\frac{L}{2}+1},\beta_{L-1}}$
$\displaystyle=$
$\displaystyle\sum_{\beta_{0},\beta_{\frac{L}{2}}}B^{1}_{\beta_{0},\beta_{1}}*B^{\frac{L}{2}}_{\beta_{\frac{L}{2}-1},\beta_{\frac{L}{2}}}$
$\displaystyle*$ $\displaystyle
B^{\frac{L}{2}+1}_{\beta_{\frac{L}{2}},\beta_{\frac{L}{2}+1}}*B^{L}_{\beta_{L},\beta_{0}},$
$\displaystyle D^{i}_{\beta_{i-1},\beta_{L-i+1},\beta_{i},\beta_{L-i}}$
$\displaystyle=$ $\displaystyle
B^{i}_{\beta_{i-1},\beta_{i}}*B^{L-i+1}_{\beta_{L-i},\beta_{L-i+1}}$ (34)

and the tensor multiplication $*$ denotes the contraction of physical indices
as shown in Fig. 9(b).

Therefore, the topological twist that represents the exchange statistics is
calculated from the ratio of the expectation value of $W^{\rm twist}_{P_{1}}$
and $W^{\rm twist}_{P_{2}}$,

$\displaystyle\tilde{T}_{a}=\frac{\langle\psi|W^{\rm
twist}_{P_{1};a}|\psi\rangle}{\langle\psi|W^{\rm
twist}_{P_{2};a}|\psi\rangle},$ (35)

where the label $a$ denotes the anyon type.

##  V Numerical results

Figure 10: (a) Schematic of toric code model on an infinite plane. The red
square is the region supporting the WLO. (b) The red and blue square are
regions supporting the WLOs for calculating braiding statistic as described in
II.3. The circled region is region B defined in Eq. (9). (c) The curve is the
region supporting the WLOs for calculating topological twist as described in
IV.5. In the twist region, the blue curve is on top of the green curve. (d)
The minimum costs as function of uniform magnetic field $h_{x}$ for three
anyon types in the toric code model with $h_{z}=0.05$ and $\chi=1$. For $t=1$,
the support of the Wilson loops is the perimeter of a $(L_{x},Ly)=(6,4)$
square. For $t=2$, the support is depicted in (a). (e) The expectation value
of the closed WLOs in renormalization group fixed points ($h_{x}=h_{z}=0$),
$W_{\rm FP}$ and the optimized closed WLOs as function of the magnetic field
$h_{x}$. In this figure, $h_{z}=0.05$, $L_{x}=36$, $L_{y}=6$, $\chi=1$ and
$t=1$. (f) The error of the Wilson loop operators $\epsilon=|\langle W-W_{\rm
exact}\rangle|$ as function of the inverse of correlation length ($1/\xi$) of
the ground state wave function of toric code in uniform magnetic field with
$h_{z}=0.05$, $0.04\leq h_{x}\leq 0.1$, $\chi=1$. (g) The error of the Wilson
loop operators $\epsilon$ as function of the number of sites $N_{s}$ with
$h_{x}=0.1$, $h_{z}=0.05$, $t=1$, $\chi=1$, $L_{y}=4$. We vary $N_{s}$ by
increasing the side length $L_{x}$ from $L_{x}=4$ to $L_{x}=8$.

In this section, we present the numerical results for extracting WLOs and
braiding statistics for various systems.

###  V.1 $\mathbb{Z}_{2}$ toric code model in a magnetic field

We first consider the $Z_{2}$ toric code model in a magnetic (Zeeman) field
wu2012phase ; vidal2009low ; tupitsyn2010topological ; halasz2012probing ;
ritz2021wegner . The Hamiltonian is of the form

$H_{{\rm TC}}=-\sum_{p\in{\rm plaquette}}B_{p}-\sum_{v\in{\rm
vertex}}A_{v}-\sum_{i}(h_{x}X_{i}+h_{z}Z_{i}),$ (36)

where the plaquette operators $B_{p}=\prod_{i\in p}Z_{i}$, the vertex
operators $A_{v}=\prod_{i\in v}X_{i}$ and $h_{x}$ and $h_{z}$ are the magnetic
fields along the $x$ and $z$ directions respectively.

In the following, we consider a dual lattice, so that the spins are on the
lattice sites instead of the bonds. The plaquette operators $A_{p}$ and the
vertice operators $B_{v}$ of the toric code model are then on alternating
plaquettes, as depicted in gray and white respectively in Fig. 10 (a).

When the magnetic fields $h_{x}=h_{z}=0$, the ground state of the toric code
model can be solved analytically. The ground state has zero correlation length
and the Wilson loop operators can be solved exactly (cost $\mathcal{C}=0$)
with bond dimension $\chi=1$ and thickness $t=1$ for any size of the closed
WLOs $L_{x}$ and $L_{y}$.

There exist four distinct Wilson loop operators for the anyon types : $I$,
$e$, $m$ and $\psi$, where $I$ is the identity sector, $e$ and $m$ are bosons
with trivial self-braiding phase and a $\pi$ mutual-braiding phase and
$\psi=e\times m$ is the fermion. The modular $S$ matrix characterizing the
mutual braiding statistics, listing the anyons in the order $I$, $e$, $m$,
$\psi$ is of the form,

$\displaystyle S=\frac{1}{2}\begin{pmatrix}1&1&1&1\\\ 1&1&-1&-1\\\
1&-1&1&-1\\\ 1&-1&-1&1\end{pmatrix}.$ (37)

The topological twists, which characterize the exchange statistics of the
anyons, is given by

$\displaystyle T=diag(1,1,1,-1).$ (38)

Figure 10 (b) and (c) illustrate the measurement of modular $S$ and $T$
matrices with thickness $t=1$ on 2D square lattice.

With non-zero magnetic fields, the ground state wave function has a finite
correlation length $\xi$. The Wilson loop operators with finite thickness $t$
can no longer represent the exact WLOs. However, the optimization scheme can
still find approximate WLOs with error of order
$\mathcal{O}(N_{s}e^{-t/\xi})$, where $N_{s}$ is the number of sites.

We first obtain the ground state wave function by minimizing the infinite
projective entangled pair state (iPEPS) using a recently proposed differential
programming approach liao2019differentiable ; crone2020detecting  with corner
transfer matrix renormalization group (CTMRG) orus2009simulation .

Next, we follow the protocol described in Sec. III, minimizing the cost
function defined in Eq. (15) starting from random MPO initializations, to find
the Wilson loop operators. Figure 10(d) shows the minimum cost as a function
of $h_{x}$ with a fixed $h_{z}=0.05$ for thickness $t=1$ and $t=2$. The cost
decreases when the thickess $t$ increases, demonstrating that a larger
thickness gives a better approximation to the true Wilson loop operator.

In Fig. 10(e), we compare the optimized closed WLOs to the WLOs known for the
fixed point toric code Hamiltonian ($h_{x}=h_{z}=0$). The expectation values
of the optimized closed WLOs stay close to $1$ as the magnetic field is
increased. However, the expectation values of the fixed point WLOs decrease
with increasing $h_{x}$. We note that this is remarkable given that we do not
use prior knowldge of the WLOs of the fixed point Hamiltonian; our scheme is
completely unbiased and uses no prior information aside from the ground state
of the perturbed Hamiltonian whose WLOs we are trying to find.

While we do not know the form of the exact WLOs in the presence of non-zero
$h_{x}$, $h_{z}$, we do know that their expectation value in the ground state
is 1. Therefore, we can compute how close our optimized WLOs are to the exact
WLOs. We define the error

$\displaystyle\epsilon=|\langle\psi|W-W_{\rm exact}|\psi\rangle|,$ (39)

where we use the fact that the exact WLO, $W_{\rm exact}$, satisfies $\langle
W_{\rm exact}\rangle=1$. Here $W$ represents a WLO found through our
optimization protocol.

We present the error $\epsilon$ of the WLOs as a function of the inverse of
the correlation length $1/\xi$ and number of sites $N_{s}$ in Fig. 10(f) and
(g), respectively The correlation length is computed by through the
exponential decay of the correlation funcion, $\langle
X_{x_{0},y_{0}}X_{x_{0}+d,y_{0}}\rangle-\langle X_{x_{0},y_{0}}\rangle\langle
X_{x_{0}+d,y_{0}}\rangle=Ae^{-d/\xi}$. The largest correlation length we reach
in our simulation is $\xi=1.21$.

In Fig. 10(f), we vary the correlation length by varying the magnetic field
$h_{x}$. We show that the error drops exponentially as function of $1/\xi$.
Fig. 10 (g) shows how the error $\epsilon$ scales with the number of sites
$N_{s}$ in the support of the WLO with $t=1$, indicating that the error scales
up linearly with $N_{s}$. Therefore the error scaling is consistent with the
error bound $\mathcal{O}(N_{s}e^{-t/\xi})$ hastings2005quasiadiabatic .

Finally, we numerically evaluate the twist product matrix $\tilde{S}_{i,j}$
and the topological twist $\tilde{T}_{i}$. The twist product matrix is
consistent with Eq. (37) up to $10^{-5}$ error and the topological twist has
error up to $10^{-2}$. For example, the twist product matrix for $h_{x}=0.1$,
$h_{z}=0.05$, $t=2$, $L_{x}=6$, $L_{y}=4$, $\chi=1$ is

$\displaystyle\tilde{S}=\begin{pmatrix}1.0&1.0&1.0&1.0\\\
1.0&1.0&-1.00001&-0.99996\\\ 1.0&-0.99999&1.0&-0.99977\\\
1.0&-1.00004&-1.00024&1.0\\\ \end{pmatrix},$ (40)

and

$\displaystyle\tilde{T}=diag(1.0,0.963,0.912,-0.974).$ (41)

###  V.2 Double semion model

Figure 11: (a) Schematic of the plaquette and vertex for double semion model.
The “legs” of a plaquette are the edges that stick out of the plaquette. (b)
Schematic of the L-vertex, R-vertex and R-leg for double semion model. (c) and
(d) The support of the Wilson loop operators with thickness $t=1$ and $t=2$
respectively. (e) The minimum cost as a function of uniform magnetic field
$h_{x}$ for two semions in the double semion model. The support of the Wilson
loops are depicted in (c) and (d). (f) The minimum cost as a function of the
bond dimension $\chi$ for the two semions for $h_{x}=0.1$ and $t=1$

In this section, we present our numerical results for the double semion model.

The Hamiltonian of the double semion model levin2005string , which we take to
be on the honeycomb lattice, is of the form

$\displaystyle H_{DS}$ $\displaystyle=$ $\displaystyle-\sum_{v\in
vert.}A_{v}-\sum_{i}h_{x}X_{i},$ (42) $\displaystyle+$
$\displaystyle(\sum_{p\in
plaq.}B_{p}\prod_{j\in\rm{legs}(p)}i^{\frac{1-X_{j}}{2}}+h.c.)$

where the plaquette operators $B_{p}=\prod_{i\in p}Z_{i}$, the vertex
operators $A_{v}=\prod_{i\in v}X_{i}$ and legs$(p)$ is the legs of plaquette
$p$ as shown in Fig. 11(a).

For $h_{x}=0$, the ground state and the Wilson loop operators can be obtained
analytically levin2005string . There exists four distinct WLOs denoted by $I$,
$s$, $s^{\prime}$, $b=s\times s^{\prime}$, which represent identity, right and
left-handed semions and a boson. The two semions have $\pi$ self-braiding
statistics and trivial mutual braiding statistics. The modular $S$ matrix,
listing the anyons in the order $I$, $s$, $s^{\prime}$, $b$, is of the form

$\displaystyle S=\frac{1}{2}\begin{pmatrix}1&1&1&1\\\ 1&-1&1&-1\\\
1&1&-1&-1\\\ 1&-1&-1&1\end{pmatrix}.$ (43)

The topological twists that characterize the exchnage statistics of the double
semion model is given by

$\displaystyle T=diag(1,i,-i,1).$ (44)

The exact closed WLOs for $h_{x}=0$ are of the form

$\displaystyle W_{I}$ $\displaystyle=$ $\displaystyle I,$ $\displaystyle
W_{s}$ $\displaystyle=$
$\displaystyle\prod_{l\in\gamma}Z_{l}\prod_{k\in\rm{L-vertex}}(-1)^{\frac{1}{4}(1-X_{i})(1+X_{j})}\prod_{l\in\rm{R-leg}}i^{\frac{1}{2}(1-X_{l})}$
$\displaystyle W_{s^{\prime}}$ $\displaystyle=$
$\displaystyle\prod_{l\in\gamma}Z_{l}\prod_{k\in\rm{L-vertex}}(-1)^{\frac{1}{4}(1-X_{i})(1+X_{j})}\prod_{l\in\rm{R-leg}}(-i)^{\frac{1}{2}(1-X_{l})}$
$\displaystyle W_{b}$ $\displaystyle=$
$\displaystyle\prod_{l\in\mathcal{\gamma}^{\prime}}X_{l},$ (45)

where $\mathcal{\gamma}$ is a path of a closed WLO, $\gamma^{\prime}$ is a
path of a closed WLO on the duel lattice, $i$ and $j$ are two legs attached to
L-vertex $k$ as shown in Fig. 11(b).

Note that the $W_{s}$ and $W_{s^{\prime}}$ operators above have thickness
$t=2$ and bond dimension $\chi=1$. Interestingly, using our unbiased numerical
optimization, we found an equivalent way to represent the exact WLOs at the
fixed point with $t=1$ and $\chi=2$. We present its analytical form in
Appendix A.

In order to numerically optimize the WLOs, we numerically optimize the iPEPS
ground state wave function corboz2012simplex  and consider WLOs with thickness
$t=1$ and $t=2$ as shown in Fig. 11(c) and 11(d) respectively. Fig. 11(e)
shows the minimum cost achieved as a function of the magnetic field $h_{x}$.
The optimizer converges to these minimum costs after roughly 1000 iterations.
Fig. 11(f) shows the minimum cost as a function of the bond dimension $\chi$
for the two semions for $h_{x}=0.1$ and $t=1$.

Figure 12: (a) The red and blue curves are regions supporting the WLOs for
calculating braiding statistic for double semion model. (b) The twisted curve
is the regions supporting the WLOs for calculating topological twist as
described in IV.5. In the overlaping region, the blue curve is on top of the
green curve.

We then numerically evaluate the twist product matrix $\tilde{S}_{i,j}$ and
the topological twist $\tilde{T}_{i}$. The support of $\tilde{S}$ and
$\tilde{T}$ are shown in Fig. 12(a) and (b) respectively. The twist product
matrix is consistent with Eq. (43) up to $10^{-4}$ error and the topological
twist has error up to $10^{-2}$. For example, the twist product matrix for
$h_{x}=0.1$, $t=1$, $\chi=2$ is

$\displaystyle\tilde{S}=\begin{pmatrix}1.0&1.0&1.0&1.0\\\
1.0&-1.0&1.001&-0.994\\\ 1.0&0.999&-1.0&-1.002\\\
1.0&-1.006&-0.998&1.0\end{pmatrix},$ (46)

and

$\displaystyle\tilde{T}={\rm diag}(1.0,0.9351e^{i\pi 0.502},0.9412e^{-i\pi
0.489},0.993).$ (47)

The correlation length of the ground state, calculated through the exponential
decay of the correlation funcion, $\langle
X_{x_{0},y_{0}}X_{x_{0}+d,y_{0}}\rangle-\langle X_{x_{0},y_{0}}\rangle\langle
X_{x_{0}+d,y_{0}}\rangle$, with $h_{x}=0.1$ is $\xi=0.87$.

##  VI Summary and outlook

In conclusion, we propose a numerical optimization-based scheme to extract
Wilson loop operators of anyons from a single ground state wave function
defined on a simply connected region of space. We show how, after extracting
closed loop operators, one can then modify them to obtain Wilson line
operators that create, move, and annihilate anyons. This allows us to
ultimately extract the braiding statistics and topological twists of the
anyons from a single bulk ground state wave function. While our protocol for
extracting the modular $S$ matrix is expected to be general, our protocol for
extracting the topological twists may be benefiting from the lattice
symmetries of our models; we leave a systematic investigation of this for
future work.

Our algorithm fully succeeds only when all distinct equivalence classes of
Wilson loop operators have been found. We expect that in general, for a large
enough bond dimension, all Wilson loop operators can be captured by our matrix
product operator ansatz. In practice it may be the case that some Wilson loop
operators might be more complicated than others, in the sense of having higher
operator entanglement or requiring larger bond dimension. In this case, if
there is an implicit bias of the optimization procedure towards the simpler
loop operators, then the algorithm may never succeed in discovering a complete
set of loop operators starting from random initialization. In this case,
further work needs to be done on either finding improved initializations or
modifying the optimization procedure to remove the implicit biases.

Our work also raises an interesting question of whether all Wilson loop
operators for anyons can always be described by an MPO with finite bond
dimension. This is particularly intriguing to study for chiral topological
orders, such as fractional quantum Hall states, where loop operators have
never been explicitly computed in MPO form and the ground state wave functions
cannot be described by a PEPS with finite bond dimension.

So far, our optimization scheme is tailored to Abelian topological orders. It
would be interesting to generalize it to the case of non-abelian topologically
ordered phases, which can have anyons with quantum dimension greater than one.
One possible direction is to design an optimization scheme to extract both the
WLOs and the quantum dimension of anyons simultaneously.

As we discussed, the Wilson line operators can be generalized to more generic
movement and splitting operators. If we discover such generic movement and
splitting operators through a similar optimization approach to what we have
described, it may also be possible to extract the $F$ and $R$ symbols of the
underlying modular tensor category using the results of Ref.
kawagoe2020microscopic .

Looking further, similar optimization procedures applied to symmetry defect
creation and movement operators may eventually allow us to extract the full
$G$-crossed braided tensor category barkeshli2019  that describes a given
symmetry-enriched topological ground state. This would allow the extraction of
all possible topological invariants from a single bulk ground state wave
function.

More broadly, recent developments in quantum simulators allow the realization
of topologically ordered states that might not occur in conventional
electronic matter Googletoric2021 ; HarvardQSL1 . Given this opportunity, it
is intriguing to develop measurement protocols for probing topological
properties of a ground state wave function associated with a prior unknown
gapped Hamiltonian. Our optimization protocol may be particularly relevant in
this context. Since our scheme requires measuring several observables for a
given wave function, it may be useful when combined with shadow tomography.
huang2020predicting .

##  VII Acknowledgements

This work is supported by NSF CAREER DMR- 1753240 (MB), AFOSR-MURI FA9550-
19-1-0399, ARO W911NF2010232, ONR N00014-20-1-2325, National Science
Foundation QLCI grant OMA-2120757, U.S. Department of Energy, Quantum Systems
Accelerator program and the Simons and Minta Martin Foundations (MH,ZC).

##  Appendix A Alternative definition of Wilson loop operator for double
semion model

In Eq. (45), the operators applied on the R-leg $l$ are $(\pm
i)^{\frac{1}{2}(1-X_{l})}$. The operator in the exponent
$\frac{1}{2}(1-X_{l})$ is equivalent to $\frac{1}{2}(1+X_{a}X_{b})$ in the
space spanned by close string configurations as shown in fig. 13(a), where $a$
and $b$ are the two sites attached to the same vertex of the R-leg $l$.
Therefore, the Wilson loop operators can be rewritten as

$\displaystyle W_{s}$ $\displaystyle=$
$\displaystyle\prod_{i\in\gamma}Z_{i}\prod_{k\in
L}(-1)^{\frac{1}{4}(1-X_{i})(1+X_{j})}\prod_{l\in
R}i^{\frac{1}{2}(1-X_{a}X_{b})}$ $\displaystyle W_{s^{\prime}}$
$\displaystyle=$ $\displaystyle\prod_{i\in\gamma}Z_{i}\prod_{k\in
L}(-1)^{\frac{1}{4}(1-X_{i})(1+X_{j})}\prod_{l\in
R}(-i)^{\frac{1}{2}(1-X_{a}X_{b})},$

where $\mathcal{\gamma}$ is a path of a closed WLO, $R$ denotes R-vertex and
$a$ and $b$ are two legs attached to R-vertex $l$ as shown in fig. 11(b).

Figure 13: Alternative representation of $W_{s}$ and $W_{s^{\prime}}$.

To show that the WLOs in Eq. (LABEL:Wilson_DS_alternative) can be represented
by MPOs with $t=1$, $\chi=2$, we consider an open Wilson operator shown in
fig. 13(b) which creates semions on the both ends of the Wilson operator for
illustration. A closed WLO can be defined in a similar fashion. To simplify
the notation, we define $O_{R}=(\pm i)^{\frac{1}{2}(1-X_{a}X_{b}})$ for right
and left-handed semion respectively and
$O_{L}=(-1)^{\frac{1}{4}(1-X_{i})(1+X_{j})}$. The Wilson loop operators can be
represented in the tensor network notation as shown in 13(b). The operators
$O_{R}$ and $O_{L}$ can be decomposed through singular value decomposition
(SVD) as

$\displaystyle O_{R}=\sum_{\alpha=1}^{2}O_{R1,\alpha}O_{R_{2},\alpha}$
$\displaystyle O_{L}=\sum_{\alpha=1}^{2}O_{L1,\alpha}O_{L_{2},\alpha},$ (49)

where $\alpha$ is the index auxiliary bond as shown in fig. 13(c) and

$\displaystyle O_{R_{1}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{i}{2}}\begin{pmatrix}I&-iX\end{pmatrix},~{}O_{R_{2}}=\begin{pmatrix}I\\\
X\end{pmatrix}$ $\displaystyle O_{L_{1}}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}I&-2P_{-}\end{pmatrix},~{}O_{L_{2}}=\begin{pmatrix}I\\\
P_{+}\end{pmatrix},$ (50)

where $P_{\pm}=\frac{1}{2}(1\pm X)$. After SVD, the Wilson operators can be
expressed by a two-site periodic MPO terminated by edge tensors :

$\displaystyle
W_{s}=\sum_{\\{\alpha\\}}E_{1\alpha_{1}}A_{\alpha_{1}\alpha_{2}}B_{\alpha_{2}\alpha_{3}}\dots
A_{\alpha_{L-2}\alpha_{L-1}}B_{\alpha_{L-1}\alpha_{L}}E_{2\alpha_{L}},$

where $A_{\alpha,\beta}=ZO_{L_{1},\beta}O_{R_{2},\alpha}$,
$B_{\alpha,\beta}=ZO_{L_{2},\alpha}O_{R_{1},\beta}$, $E_{1(2)}=ZO_{R_{1(2)}}$
and $\alpha,\beta=1,2$. $W_{s^{\prime}}$ can be defined similarly. Therefore,
the WLO for creating semions can be expressed by MPO with $t=1$ and $\chi=2$.

## References

  * [1] Xiao-Gang Wen.  Quantum Field Theory of Many-Body Systems.  Oxford Univ. Press, Oxford, 2004.
  * [2] Chetan Nayak, Steven H. Simon, Ady Stern, Michael Freedman, and Sankar Das Sarma.  Non-abelian anyons and topological quantum computation.  Rev. Mod. Phys., 80:1083, 2008.
  * [3] Zhenghan Wang.  Topological Quantum Computation.  American Mathematical Society, 2008.
  * [4] T. Senthil.  Symmetry-protected topological phases of quantum matter.  Annual Review of Condensed Matter Physics, 6:299–324, 2015.
  * [5] Bei Zeng, Xie Chen, Duan-Lu Zhou, Xiao-Gang Wen, et al.  Quantum information meets quantum matter.  Springer, 2019.
  * [6] Maissam Barkeshli, Parsa Bonderson, Meng Cheng, and Zhenghan Wang.  Symmetry fractionalization, defects, and gauging of topological phases.  Phys. Rev. B, 100:115147, Sep 2019.
  * [7] Kyle Kawagoe and Michael Levin.  Microscopic definitions of anyon data.  Physical Review B, 101(11):115113, 2020.
  * [8] Daniel Bulmash and Maissam Barkeshli.  Fermionic symmetry fractionalization in 2+1 dimensions.  Physical Review B, 105(12), mar 2022.
  * [9] David Aasen, Parsa Bonderson, and Christina Knapp.  Characterization and classification of fermionic symmetry enriched topological phases.  arXiv preprint arXiv:2109.10911, 2021.
  * [10] K. J. Satzinger et al.  Realizing topologically ordered states on a quantum processor.  Science, 374(6572):1237–1241, 2021.
  * [11] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletic, and M. D. Lukin.  Probing topological spin liquids on a programmable quantum simulator.  Science, 374(6572):1242–1247, 2021.
  * [12] Dolev Bluvstein, Harry Levine, Giulia Semeghini, Tout T Wang, Sepehr Ebadi, Marcin Kalinowski, Alexander Keesling, Nishad Maskara, Hannes Pichler, Markus Greiner, et al.  A quantum processor based on coherent transport of entangled atom arrays.  arXiv preprint arXiv:2112.03923, 2021.
  * [13] Matthew B Hastings and Xiao-Gang Wen.  Quasiadiabatic continuation of quantum states: The stability of topological ground-state degeneracy and emergent gauge invariance.  Physical review b, 72(4):045141, 2005.
  * [14] Davide Gaiotto, Anton Kapustin, Nathan Seiberg, and Brian Willett.  Generalized global symmetries.  2014\.
  * [15] A.Yu. Kitaev.  Fault-tolerant quantum computation by anyons.  Annals Phys., 303:2–30, 2003.
  * [16] Michael A Levin and Xiao-Gang Wen.  String-net condensation: A physical mechanism for topological phases.  Physical Review B, 71(4):045110, 2005.
  * [17] Alexei Kitaev and John Preskill.  Topological entanglement entropy.  Physical review letters, 96(11):110404, 2006.
  * [18] Michael Levin and Xiao-Gang Wen.  Detecting topological order in a ground state wave function.  Physical review letters, 96(11):110405, 2006.
  * [19] Hossein Dehghani, Ze-Pei Cian, Mohammad Hafezi, and Maissam Barkeshli.  Extraction of the many-body chern number from a single wave function.  Physical Review B, 103(7):075102, 2021.
  * [20] Ze-Pei Cian, Hossein Dehghani, Andreas Elben, Benoît Vermersch, Guanyu Zhu, Maissam Barkeshli, Peter Zoller, and Mohammad Hafezi.  Many-body chern number from statistical correlations of randomized measurements.  Physical Review Letters, 126(5):050501, 2021.
  * [21] Ruihua Fan, Rahul Sahay, and Ashvin Vishwanath.  Extracting the quantum hall conductance from a single bulk wavefunction.  arXiv preprint arXiv:2208.11710, 2022.
  * [22] Ken Shiozaki and Shinsei Ryu.  Matrix product states and equivariant topological field theories for bosonic symmetry-protected topological phases in (1+1) dimensions.  Journal of High Energy Physics, 2017(4):100, Apr 2017.
  * [23] Ken Shiozaki, Hassan Shapourian, Kiyonori Gomi, and Shinsei Ryu.  Many-body topological invariants for fermionic short-range entangled topological phases protected by antiunitary symmetries.  Phys. Rev. B, 98:035151, Jul 2018.
  * [24] Andreas Elben, Jinlong Yu, Guanyu Zhu, Mohammad Hafezi, Frank Pollmann, Peter Zoller, and Benoît Vermersch.  Many-body topological invariants from randomized measurements in synthetic quantum matter.  Science advances, 6(15):eaaz3666, 2020.
  * [25] Hong-Hao Tu, Yi Zhang, and Xiao-Liang Qi.  Momentum polarization: An entanglement measure of topological spin and chiral central charge.  Physical Review B, 88(19), Nov 2013.
  * [26] Michael P. Zaletel, Roger S. K. Mong, and Frank Pollmann.  Topological characterization of fractional quantum hall ground states from microscopic hamiltonians.  Physical Review Letters, 110(23), Jun 2013.
  * [27] Hui Li and F. D. M. Haldane.  Entanglement spectrum as a generalization of entanglement entropy: Identification of topological order in non-abelian fractional quantum hall effect states.  Phys. Rev. Lett., 101:010504, Jul 2008.
  * [28] Isaac H Kim, Bowen Shi, Kohtaro Kato, and Victor V Albert.  Modular commutator in gapped quantum many-body systems.  arXiv preprint arXiv:2110.10400, 2021.
  * [29] Isaac H Kim, Bowen Shi, Kohtaro Kato, and Victor V Albert.  Chiral central charge from a single bulk wave function.  arXiv preprint arXiv:2110.06932, 2021.
  * [30] Yi Zhang, Tarun Grover, Ari Turner, Masaki Oshikawa, and Ashvin Vishwanath.  Quasiparticle statistics and braiding from ground-state entanglement.  Physical Review B, 85(23):235151, 2012.
  * [31] Yi Zhang, Tarun Grover, and Ashvin Vishwanath.  General procedure for determining braiding and statistics of anyons using entanglement interferometry.  Physical Review B, 91(3), Jan 2015.
  * [32] X.G. Wen.  Topological orders in rigid states.  Int. J. Mod. Phys., B:239, 1990.
  * [33] Guanyu Zhu, Mohammad Hafezi, and Maissam Barkeshli.  Quantum origami: Transversal gates for quantum computation and measurement of topological order.  Physical Review Research, 2(1), Mar 2020.
  * [34] Jacob C. Bridgeman, Steven T. Flammia, and David Poulin.  Detecting topological order with ribbon operators.  Physical Review B, 94(20), Nov 2016.
  * [35] Thorsten B Wahl and Benjamin Béri.  Local integrals of motion for topologically ordered many-body localized systems.  Physical Review Research, 2(3):033099, 2020.
  * [36] Alexei Kitaev.  Anyons in an exactly solved model and beyond.  Ann. Phys., 321(1):2 – 111, 2006.
  * [37] Parsa Hassan Bonderson.  Non-Abelian Anyons and Interferometry.  PhD thesis, California Institute of Technology, 2007.
  * [38] Jeongwan Haah.  An invariant of topologically ordered states under local unitary transformations.  Communications in Mathematical Physics, 342(3):771–801, 2016.
  * [39] Frank Verstraete, Valentin Murg, and J Ignacio Cirac.  Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems.  Advances in physics, 57(2):143–224, 2008.
  * [40] Nick Bultinck, Michael Mariën, Dominic J Williamson, Mehmet B Şahinoğlu, Jutho Haegeman, and Frank Verstraete.  Anyons and matrix product operator algebras.  Annals of physics, 378:183–233, 2017.
  * [41] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [42] Román Orús and Guifré Vidal.  Simulation of two-dimensional quantum systems on an infinite lattice revisited: Corner transfer matrix for tensor contraction.  Physical Review B, 80(9):094403, 2009.
  * [43] Hai-Jun Liao, Jin-Guo Liu, Lei Wang, and Tao Xiang.  Differentiable programming tensor networks.  Physical Review X, 9(3):031041, 2019.
  * [44] SPG Crone and P Corboz.  Detecting a z 2 topologically ordered phase from unbiased infinite projected entangled-pair state simulations.  Physical Review B, 101(11):115143, 2020.
  * [45] Fengcheng Wu, Youjin Deng, and Nikolay Prokof’ev.  Phase diagram of the toric code model in a parallel magnetic field.  Physical Review B, 85(19):195104, 2012.
  * [46] Julien Vidal, Sébastien Dusuel, and Kai Phillip Schmidt.  Low-energy effective theory of the toric code model in a parallel magnetic field.  Physical Review B, 79(3):033109, 2009.
  * [47] IS Tupitsyn, Alexei Kitaev, NV Prokof Ev, and PCE Stamp.  Topological multicritical point in the phase diagram of the toric code model and three-dimensional lattice gauge higgs model.  Physical Review B, 82(8):085114, 2010.
  * [48] Gábor B Halász and Alioscia Hamma.  Probing topological order with rényi entropy.  Physical Review A, 86(6):062330, 2012.
  * [49] Anna Ritz-Zwilling, Jean-Noël Fuchs, and Julien Vidal.  Wegner-wilson loops in string nets.  Physical Review B, 103(7):075128, 2021.
  * [50] Philippe Corboz, Karlo Penc, Frédéric Mila, and Andreas M Läuchli.  Simplex solids in su (n) heisenberg models on the kagome and checkerboard lattices.  Physical Review B, 86(4):041106, 2012.
  * [51] Hsin-Yuan Huang, Richard Kueng, and John Preskill.  Predicting many properties of a quantum system from very few measurements.  Nature Physics, 16(10):1050–1057, 2020.