Classifying topological neural network quantum states via diffusion maps
We discuss and demonstrate an unsupervised machine-learning procedure to
detect topological order in quantum many-body systems. Using a restricted
Boltzmann machine to define a variational ansatz for the low-energy spectrum,
we sample wave functions with probability decaying exponentially with their
variational energy; this defines our training dataset that we use as input to a
diffusion map scheme. The diffusion map provides a low-dimensional embedding of
the wave functions, revealing the presence or absence of superselection sectors
and, thus, topological order. We show that for the diffusion map, the required
similarity measure of quantum states can be defined in terms of the network
parameters, allowing for an efficient evaluation within polynomial time.
However, possible ''gauge redundancies'' have to be carefully taken into
account. As an explicit example, we apply the method to the toric code.
Yanting Teng, Subir Sachdev, Mathias S. Scheurer
# Classifying topological neural network quantum states via diffusion maps

###### Abstract

We discuss and demonstrate an unsupervised machine-learning procedure to detect topological order in quantum many-body systems. Using a restricted Boltzmann machine to define a variational ansatz for the low-energy spectrum, we sample wave functions with probability decaying exponentially with their variational energy; this defines our training dataset that we use as input to a diffusion map scheme. The diffusion map provides a low-dimensional embedding of the wave functions, revealing the presence or absence of superselection sectors and, thus, topological order. We show that for the diffusion map, the required similarity measure of quantum states can be defined in terms of the network parameters, allowing for an efficient evaluation within polynomial time. However, possible "gauge redundancies" have to be carefully taken into account. As an explicit example, we apply the method to the toric code.

## I Introduction

In the last few years, machine learning (ML) techniques have been very actively studied as novel tools in many-body physics [1; 2; 3; 4; 5; 6; 7]. A variety of valuable applications of ML has been established, such as ML-based variational ans\(\hat{\text{a}}\)ze for many-body wave functions, application of ML to experimental data to extract information about the underlying physics, ML methods for more efficient Monte-Carlo sampling, and employment of ML to detect phase transitions, to name a few. Regarding the latter type of applications, a particular focus has recently been on topological phase transitions [8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31]. This is motivated by the challenges associated with capturing topological phase transitions: by definition, topological features are related to the global connectivity of the dataset rather than local similarity of samples. Therefore, unless the dataset is sufficiently simple such that topologically connected pairs of samples also happen to be locally similar or features are used as input data that are closely related to the underlying topological invariant, the topological structure is hard to capture reliably with many standard ML techniques [11; 12].

In this regard, the ML approach proposed in Ref. [12], which is based on diffusion maps (DM) [32; 33; 34; 35], is a particularly promising route to learn topological phase transitions; it allows to embed high-dimensional data in a low-dimensional subspace such that pairs of samples that are smoothly connected in the dataset will be mapped close to each other, while disconnected pairs will be mapped to distant points. As such, the method captures the central notion of topology. In combination with the fact that it is unsupervised and thus does not require _a priori_ knowledge of the underlying topological invariants, it is ideally suited for the task of topological phase classification. As a result, there have been many recent efforts applying this approach to a variety of problems, such as different symmetry-protected, including non-Hermitian, topological systems [36; 37; 38; 39; 40; 41], experimental data [42; 39], many-body localized states [43], and dynamics [44]; extensions based on combining DM with path finding [36] as well as with quantum computing schemes [45] for speed-up have also been studied.

As alluded to above, another very actively pursued application of ML in physics are neural network quantum states: as proposed in Ref. [46], neural networks can be used to efficiently parameterize and, in many cases, optimize variational descriptions of wave functions of quantum many-body systems [47; 48; 49; 50; 51; 52; 53; 54; 55; 56]. In particular, restricted Boltzmann machines (RBMs) [4] represent a very popular neural-network structure in this context. For instance, the ground states of the toric code model [57] can be exactly expressed with a _local_ RBM ansatz [58], i.e., where only neighboring spins are connected to the same hidden neurons. When additional non-local extensions to the RBM ansatz of Ref. [58] are added, this has been shown to also provide a very accurate variational description of the toric code in the presence of a magnetic field [59].

In this work, we combine the DM approach of Ref. [12] with neural network quantum states with the goal of capturing topological order in an unsupervised way in interacting quantum many-body systems. We use a local network ansatz, with parameters \(\Lambda\), as a variational description for the wave functions \(\ket{\Psi(\Lambda)}\) of the low-energy subspace of a system with Hamiltonian \(\hat{\mathcal{H}}\). While we also briefly mention other possible ways of generating ensembles of states, we primarily focus on an energetic principle: we sample wavefunctions such that the probability of \(\ket{\Psi(\lambda)}\) is proportional to \(\exp(-\bra{\hat{\mathcal{H}}}_{\Lambda}/T)\) where \(\bra{\hat{\mathcal{H}}}_{\Lambda}=\bra{\Psi(\Lambda)}\hat{\mathcal{H}}\ket{ \Psi(\Lambda)}\). As illustrated in Fig. 1(a), the presence of superselection sectors in the low-energy spectrum of \(\hat{\mathcal{H}}\) implies that the ensemble of states decays into disconnected subsets of states for sufficiently small \(T\) (at least at fixed finite system size); these can be extracted, without need of prior labels, with dimensional reduction via DM (and subsequent \(k\)-means clustering), and thus allow to identify topological order. For sufficiently large \(T\), more and more high-energy states are included and all sectors are connected, see Fig. 1(b), as can also be readily revealed via DM-based embedding of the states.

Importantly, DM is a kernel technique in the sense that the input data \(x_{l}\) (in our case the states \(\ket{\Psi(\Lambda_{l})}\))
does not directly enter as a high-dimensional vector but only via a similarity measure \(S(x_{l},x_{l^{\prime}})\), comparing how "similar" two samples \(l\) and \(l^{\prime}\) are. In the context of applying DM to the problem of topological classification, it defines what a smooth deformation ("homotopy") of samples is. We discuss two possible such measures. The first one is just the quantum mechanical overlap, \(S_{\text{q}}(\Lambda_{l},\Lambda_{l^{\prime}})=|\langle\Psi(\Lambda_{l})|\Psi (\Lambda_{l^{\prime}})\rangle|^{2}\), of the wave functions. Although conceptually straightforward, its evaluation is computationally costly on a classical computer as it requires importance sampling. The local nature of our network ansatz allows us to also construct an alternative similarity measure that is expressed as a simple function of the network parameters \(\Lambda_{l}\) and \(\Lambda_{l^{\prime}}\) describing the two states to be compared. This can, however, lead to subtleties associated with the fact that two states with different \(\Lambda\) can correspond to the same wave functions (modulo global phase). We discuss how these "gauge redundancies" can be efficiently circumvented for generic states.

We illustrate these aspects and explicitly demonstrate the success of this approach using the toric code [57], a prototype model for topological order which has also been previously studied with other ML techniques with different focus [58, 59, 60, 15, 16, 17, 18]. We show that the DM algorithm learns the underlying loop operators wrapping around the torus without prior knowledge; at low \(T\), this leads to four clusters corresponding to the four ground states. At larger \(T\), these clusters start to merge, as expected. Interestingly, the DM still uncovers the underlying structure of the dataset related to the expectation value of the loop operators. Finally, we also show that applying a magnetic field leads to the disappearance of clusters in the DM, capturing the transition from topological order to the confined phase.

The remainder of the paper is organized as follows. In Sec. II, we describe our ML approach in general terms, including the local network quantum state description we use, the ensemble generation, a brief review of the DM scheme of Ref. [12], and the similarity measure in terms of neural network parameters. Using the toric code model as an example, all of these general aspects are then discussed in detail and illustrated in Sec. III. Finally, explicit numerical results can be found in Sec. IV and a conclusion is provided in Sec. V.

## II General algorithm

Here, we first present and discuss our algorithm [see Fig. 2(a)] in general terms before illustrating it using the toric code as an example in the subsequent sections. Consider a system of \(N\) qubits or spins, with associated operators \(\{\hat{\mathbf{s}}\}=\{\hat{\mathbf{s}}_{i},i=1,\cdots,N\}\), \(\hat{\mathbf{s}}_{i}=(\hat{s}_{\text{r}}^{x},\hat{s}_{i}^{y},\hat{s}_{i}^{z})\), and interactions governed by a local, gapped Hamiltonian \(\hat{\mathcal{H}}=\mathcal{H}(\{\hat{\mathbf{s}}\})\). We represent the states \(|\Psi(\Lambda)\rangle\) of this system using neural network quantum states [46],

\[|\Psi(\Lambda)\rangle=\sum_{\mathbf{\sigma}}\psi(\mathbf{\sigma};\,\Lambda)\ket{\mathbf{ \sigma}}, \tag{1}\]

where \(\mathbf{\sigma}\!=\!\{\sigma_{1},\sigma_{2},...,\sigma_{N}|\sigma_{i}=\pm 1\}\) enumerates configurations of the physical spin variables in a local computational basis (e.g. \(s^{z}\)-basis) and \(\Lambda\) is the set of parameters that the network \(\psi\) depends on to output the wavefunction amplitude \(\psi(\mathbf{\sigma};\,\Lambda)=\langle\mathbf{\sigma}|\Psi(\Lambda)\rangle\) for configuration \(\ket{\mathbf{\sigma}}\). Because the physical Hilbert space scales exponentially with the system size, there is a trade-off between the expressivity versus efficiency when choosing a network architecture (or ansatz) \(\psi\), so that the weights \(\Lambda\) can approximate the state \(|\Psi(\Lambda)\rangle\) to a reasonable degree and

Figure 1: (a) An illustration of a “low-energy” ensemble. Two (or more) initial states, \(|\Psi(\Lambda^{0})\rangle\) and \(|\Psi(\Lambda^{1})\rangle\), from two distinct topological sectors are chosen as “seeds” (green dots). The dots denote the dataset (later fed into the DM), which are a set of quantum states labeled by network parameters \(\Lambda\). This dataset is generated using the procedure outlined in Sec. II.1 and Algorithm. 1, where the next state \(\Lambda^{\prime}\) (blue dots at each arrow) is proposed by a random local perturbation and accepted with probability based on the energy expectation \(\langle H\rangle_{\Lambda^{\prime}}\). In the small-\(T\) regime, the full dataset is not inter-connected by such local perturbations and cluster among each topological sectors (at left and right valley). (b) An illustration of a “high-energy” ensemble. The states are generated using the same algorithm as before, however with a large hyperparameter \(T\) (compared to the energy gap \(\Delta\)). In this regime, the dataset include some of the low-energy states (blue dots), but also some high-energy states (red dots). Because the high-energy states are agnostic of the low-energy topological sectors, there exist paths (denoted by arrows among dots in the elliptical blob) such that the two initial seeds from distinct topological sectors effectively “diffuse” and form one connected cluster.
can at the same time be an efficient representation (with minimal number of parameters \(\Lambda\) that scale as a polynomial in \(N\)). To reach the ground state or, more generally, the relevant low-energy sector of the Hamiltonian \(\hat{\mathcal{H}}\) for the low-temperature physics, we minimize the energy in the variational subspace defined by Eq. (1) using gradient descent with a learning rate \(\lambda\),

\[\Lambda\to\Lambda-\lambda\,\partial_{\Lambda}\left<\hat{\mathcal{H}}\right>_{ \Lambda},\quad\left<\hat{\mathcal{H}}\right>_{\Lambda}=\left<\Psi(\Lambda) \right|\hat{\mathcal{H}}\left|\Psi(\Lambda)\right>. \tag{2}\]

Here, the quantum mechanical expectation value \(\left<\hat{\mathcal{H}}\right>_{\Lambda}\) is evaluated using importance sampling (see Appendix B).

While there are exponentially many states in the Hilbert space, the low-energy sector of a local Hamiltonian is expected to occupy a small subspace where states obey area law entanglement [61, 62] whereas a typical state obeys volume law [63, 64]. Motivated by these considerations, we consider a class of networks that naturally describe quantum states that obey area-law entanglement. Pictorially, in such networks, the connections from the hidden neurons (representing the weights \(\Lambda\)) to the physical spins are _quasi-local_[51, 53, 54, 55]. In that case, it holds

\[\psi(\boldsymbol{\sigma},\Lambda)=\phi_{1}(\boldsymbol{\sigma}_{1},\Lambda_{ 1})\times\phi_{2}(\boldsymbol{\sigma}_{2},\Lambda_{2})\times\cdots, \tag{3}\]

where \(\boldsymbol{\sigma}_{j}\!=\!\{\sigma_{k}\}_{k\in j}\) denote (overlapping) subsets of neighboring spins with \(\cup_{j}\boldsymbol{\sigma}_{j}=\boldsymbol{\sigma}\) and \(\Lambda_{j}\) are the subsets of the network parameters (weights and biases) that are connected to the physical spins in j.

```
procedure(\(\{\Lambda\}_{n=1}^{N}\))  init: optimized parameters \(\Lambda\) for\(k\) independent times do: for\(n\) sampling steps do:  Propose new parameter \(\Lambda_{p}=f(\Lambda_{t})\)  Accept with probability determined by energy \(\left<\hat{\mathcal{H}}\right>_{\Lambda}\) and parameter \(T\): \(\Lambda_{t+1}=\mathbb{P}_{\mathrm{accept}}(\Lambda^{\prime}|\Lambda;T)\) return the last \(m\) states for each \(k\): \(\{\Lambda_{i}|i=n-m,...,n\}_{k}\)
```

**Algorithm 1** Ensemble generation

### Dataset: network parameter ensembles

The dataset we use for unsupervised detection of topological order consists of an ensemble of wavefunctions \(\{|\Psi(\Lambda)\rangle\}_{l}\), parameterized by the set of network parameters \(\{\Lambda\}_{l}\). While, depending on the precise application, other choices are conceivable, we generate this ensemble such that the relative occurrence of a state \(|\Psi(\Lambda)\rangle\) is given by \(\rho_{T}(\Lambda)=\exp(-\left<\hat{\mathcal{H}}\right>_{\Lambda}/T)/Z\), with appropriate normalization factor \(Z\). As such, a small value of the "temperature-like" hyperparameter \(T\) corresponds to a "low-energy" ensemble while large \(T\) parametrize "high-energy" ensembles.

In practice, to generate this ensemble, we here first optimize the parameters \(\Lambda\) via Eq. (2) to obtain wavefunctions with lowest energy expectation values. As Eq. (1) does not contain all possible states, this will, in general, only yield approximations to the exact low-energy eigenstates of \(\hat{\mathcal{H}}\). However, as long as it is able to capture all superselection sectors of the system as well as (a subset of) higher energy states connecting these sectors, Eq. (1) will be sufficient for our purpose of detecting topological order or the absence thereof. We perform this optimization several times, \(\Lambda\to\Lambda_{l}^{0}\), with different initial conditions, to obtain several "seeds", \(\Lambda_{l}^{0}\); this is done to make sure we have a low-energy representative of all superselection sectors. Ideally the dataset is sampled directly from the the target probability distribution \(\rho_{T}\), if for instance, one has access to an experimental system at finite temperature. Here, we adopt a Markov-chain-inspired procedure for generating the ensemble based on \(\rho_{T}\) for each of these seeds. Specifically, starting from a state \(\Lambda\), we propose updates on a randomly chosen local block of parameters connected to the spins at sites j,

\[\Lambda\,\to\,\Lambda^{\prime}=\{\Lambda_{1},\Lambda_{2},\cdots,u(\Lambda_{ j}),\cdots,\Lambda_{N}\}, \tag{4}\]

where the update \(u\) only depends on \(\Lambda_{j}\). The proposed parameter \(\Lambda^{\prime}\) given the current parameter \(\Lambda\) is accepted with probability

\[\mathbb{P}_{\mathrm{accept}}(\Lambda^{\prime}|\Lambda;T)=\min\Bigl{(}1,\,e^{- \frac{\left<\hat{\mathcal{H}}\right>_{\Lambda^{\prime}}-\left<\hat{\mathcal{H }}\right>_{\Lambda}}{T}}\Bigr{)}. \tag{5}\]

This means that if the proposed state \(\Psi(\Lambda^{\prime})\) has a lower energy expectation value than \(\Psi(\Lambda)\), then the proposal will be accepted; otherwise, it will be accepted with a probability determined by the Boltzmann factor. The entire ensemble generation procedure is summarized in Algorithm 1.

### Diffusion map

As proposed in Ref. [12], DM is ideally suited as an unsupervised ML algorithm to identify the presence and number of superselection sectors in a collection of states, such as \(\{|\Psi(\Lambda)\rangle\}_{l}\) defined above. To briefly review the key idea of the DM algorithm [32, 33, 34, 35] and introduce notation, assume we are given a dataset \(X=\{x_{l}|l=1,2,...,M\}\), consisting of \(M\) samples \(x_{l}\). Below we will consider the cases \(x_{l}=\Lambda_{l}\) and \(x_{l}=|\Psi(\Lambda_{l})\rangle\); in the first case, the samples are the network parameters parametrizing the wavefunction and, in the second, the samples are the wavefunctions themselves.

To understand DM intuitively, let us define a diffusion process among states \(x_{l}\in X\). The probability of state \(x_{l}\) transitioning to \(x_{l^{\prime}}\) is defined by the Markov transition matrix element \(p_{l,l^{\prime}}\). To construct \(p_{l,l^{\prime}}\), we introduce a symmetric and positive-definite kernel \(k_{\epsilon}(x_{l},x_{l^{\prime}})\) between states \(x_{l}\) and \(x_{l^{\prime}}\). Then the transition probability matrix
\(p_{l,l^{\prime}}\) is defined as

\[p_{l,l^{\prime}}=\frac{k_{\epsilon}(x_{l},x_{l^{\prime}})}{z_{l}},\quad z_{l}=\sum_ {l^{\prime}}k_{\epsilon}(x_{l},x_{l^{\prime}}), \tag{6}\]

where the factor \(z_{l}\) ensures probability conservation, \(\sum_{l^{\prime}}p_{l,l^{\prime}}=\!1\,\forall l\). Then spectral analysis on the transition probability matrix leads to information on the _global_ connectivity of the dataset \(X\), which, in our context of \(X\) containing low-energy states, allows to identify superselection sectors and, thus, topological order [12]. To quantify how strongly two samples \(x_{l}\) and \(x_{l^{\prime}}\) are connected, one introduces the \(2t\)-step diffusion distance [32, 33, 34, 35],

\[D_{2t}(l,l^{\prime})=\sum_{l^{\prime\prime}}\frac{1}{z_{l^{\prime\prime}}}[(p ^{t})_{l,l^{\prime\prime}}-(p^{t})_{l^{\prime},l^{\prime\prime}}]^{2}, \tag{7}\]

where \(p^{t}\) denotes the \(t\)-th matrix power of the transition probability matrix \(p\). It was shown that \(D_{2t}\) can be computed from the eigenvalues \(\lambda_{n}\) and right eigenvectors \(\psi_{n}\) of the transition matrix \(p\): with \(\sum_{l^{\prime}}p_{l,l^{\prime}}\left(\psi_{n}\right)_{l^{\prime}}\!=\! \lambda_{n}\left(\psi_{n}\right)_{l}\), and in descending ordering \(\lambda_{n}>\lambda_{n+1}\), it follows

\[D_{2t}(l,l^{\prime})=\sum_{n=1}^{M-1}\lambda_{n}^{2t}[(\psi_{n})_{l}-(\psi_{n} )_{l^{\prime}}]^{2} \tag{8}\]

after straightforward algebra [35]. Geometrically, this means that the diffusion distance is represented as a Euclidean distance (weighted with \(\lambda_{n}\)) if we perform the non-linear coordinate transformation \(x_{l}\rightarrow\{(\psi_{n})_{l},n=0,\ldots M-1\}\). Furthermore, as the global connectivity is seen from the long-time limit, \(t\rightarrow\infty\), of the diffusion distance, the largest eigenvalues are most important to describe the connectivity. To be more precise, let us choose a kernel \(k_{\epsilon}\) of the form

\[k_{\epsilon}(x_{l},x_{l^{\prime}})=\exp\left(-\frac{1-S(x_{l},x_{l^{\prime}})} {\epsilon}\right), \tag{9}\]

where \(S\) is a _local similarity measure_ which obeys \(S\in[0,1]\), \(S(x_{l},x_{l^{\prime}})=S(x_{l^{\prime}},x_{l})\), and \(S(x,x)=1\). Here "local" means that \(S(x_{l},x_{l^{\prime}})=\sum_{i}\mathcal{S}_{i}(x_{l},x_{l^{\prime}})\) where \(\mathcal{S}_{i}(x_{l},x_{l^{\prime}})\) only depend on the configuration of \(x_{l}\) and

Figure 2: (a) Overview of the ML algorithm applied in this work: the “seeds” \(\{\Lambda^{0}\}\) are computed using variational Monte Carlo (see Appendix B), a Markov-chain algorithm is used to generate the network parameter ensemble dataset (Sec. II.1), then a similarity metric is used for the definition of kernels in the DM method (Sec. II.2 and Sec. II.3), and finally \(k\)-means is applied to the low-dimensional embedding in the subspace provided by the dominant DM eigenvector components. (b) The square lattice geometry for the toric code model, where the qubits \(\hat{s}_{i}\) are defined on the links of the lattice (grey dots). The Hamiltonian [given in Eq. (16)] is written in terms of the operators \(\hat{\mathcal{P}}_{P}\) (supported by spins on plaquette \(P\) denoted by the red square) and star \(\hat{\mathcal{S}}_{S}\) (supported by spins on star \(S\) denoted by the blue links). The two blue lines along \(x(y)\) directions denote the Wilson loop operators \(\hat{W}_{1,x}(\hat{W}_{2,g})\) along the straight paths \(\bar{x}(\bar{y})\). (c) An illustration of the quasi-local ansatz in Eq. (17). The ansatz is a product over local function \(\phi\) of spins in plaquette (or star), which depends on parameters \(\{w_{Xj},b_{X}\}\) for \(X=P(S)\) being plaquette (or star).
\(x_{l^{\prime}}\) in the vicinity of site \(i\). While we will discuss possible explicit forms of \(S\) for our quantum mechanical \(N\) spin/qubit system in Sec. II.3 below, a natural choice for a classical system of \(N\) spins, \(x_{l}=\{\mathbf{S}^{l}_{j},(\mathbf{S}^{l}_{i})^{2}=1,i=1,2,\ldots,N\}\), is \(S_{\text{cl}}(x_{l},x_{l^{\prime}})=\sum_{i}\mathbf{S}^{l}_{i}\cdot\mathbf{S}^{l}_{i}/N\). In Eq. (9), \(\epsilon\) plays the role of a "coarse graining" parameter that is necessary as we only deal with finite datasets \(X\): for given \(X\), we generically expect \(k_{\epsilon}(x_{l},x_{l^{\prime}})=p_{l,l^{\prime}}=\delta_{l,l^{\prime}}\) as \(\epsilon\to 0\), i.e., all samples are dissimilar if \(\epsilon\) is sufficiently small and all eigenvalues \(\Lambda_{n}\) approach 1. In turn, for \(\epsilon\to\infty\) the coarse graining parameter is so large that all samples become connected, \(k_{\epsilon}(x_{l},x_{l^{\prime}})\to 1\); as \(p_{l,l^{\prime}}\to 1/M\), we will have \(\lambda_{n>0}\to 0\), while the largest eigenvalue \(\lambda_{0}\) is always 1 (as a consequence of probability conservation). For values of \(\epsilon\) in between these extreme limits, the DM spectrum contains information about \(X\), including its topological structure: as shown in Ref. [12], the presence of \(k\in\mathbb{N}\) distinct topological equivalence classes in \(X\) is manifested by a range of \(\epsilon\) where \(\lambda_{1},\ldots\lambda_{k-1}\) are all exponentially close (in \(\epsilon\)) to 1, with a clear gap to \(\lambda_{n\geq k}\). Furthermore, the different samples \(l\) will cluster--with respect to the normal Euclidean measure, e.g., as can be captured with \(k\)-means--according to their topological equivalence class when plotted in the mapped \(k-1\)-dimensional space \(\{(\psi_{1})_{l},(\psi_{2})_{l},\ldots,(\psi_{k-1})_{l}\}\). In the following, we will use this procedure to identify the superselection sectors in the ensemble of wave functions defined in Sec. II.1. To this end, however, we first need to introduce a suitable similarity measure \(S\), to be discussed next.

### Local similarity measure

A natural generalization of the abovementioned classical similarity measure \(S_{\text{cl}}=\sum_{i}\mathbf{S}^{l}_{i}\cdot\mathbf{S}^{l}_{i}/N\), which can be thought of as the (Euclidean) inner product in the classical configuration space, is to take the inner product in the Hilbert space of the quantum system,

\[S_{\text{q}}(\Lambda_{l},\Lambda_{l^{\prime}})=|\langle\Psi(\Lambda_{l})| \Psi(\Lambda_{l^{\prime}})\rangle|^{2}. \tag{10}\]

While this or other related fidelity measures for low-rank quantum states could be estimated efficiently with quantum simulation and computing setups [65; 66; 67; 68], estimating \(S_{\text{q}}\) is generally a computationally expensive task on a classical computer, as it requires sampling over spin configurations for our variation procedure. To make the evaluation of the similarity measure more efficient, we here propose an alternative route that takes advantage of the fact that we use a local ansatz for \(\psi(\mathbf{\sigma};\Lambda)\), see Eq. (3). Our goal is to express the similarity measure directly as

\[S_{\text{n}}(\Lambda_{l},\Lambda_{l^{\prime}})=\frac{1}{N_{\text{j}}}\sum_{ \text{l}}f((\Lambda_{l})_{\text{j}},(\Lambda_{l^{\prime}})_{\text{j}}), \tag{11}\]

where \(f\) only compares a local block of parameters denoted by \(\text{j}\) and is a function that can be quickly evaluated, without having to sample spin configurations. Furthermore, \(S(x_{l},x_{l^{\prime}})=S(x_{l^{\prime}},x_{l})\) can be ensured by choosing a function \(f\) that is symmetric in its arguments and \(S\in[0,1]\) is also readily implemented by setting \(N_{\text{j}}=\sum_{\text{j}}\) and appropriate rescaling of \(f\) such that \(f\in[0,1]\). The most subtle condition is

\[S_{\text{n}}(\Lambda_{l},\Lambda_{l^{\prime}})=1\quad\Longleftrightarrow\quad |\Psi(\Lambda_{l})\rangle\propto|\Psi(\Lambda^{\prime}_{l})\rangle\,, \tag{12}\]

since, depending on the precise network architecture used for \(\psi(\mathbf{\sigma};\Lambda)\), there are "gauge transformations" \(g\in\mathcal{G}\) of the weights, \(\Lambda_{l}\to g[\Lambda_{l}]\), with

\[|\Psi(\Lambda_{l})\rangle=e^{i\vartheta_{g}}\,|\Psi(g[\Lambda_{l}])\rangle \tag{13}\]

for some global phase \(\vartheta_{g}\). We want to ensure that

\[S_{\text{n}}(\Lambda_{l},\Lambda_{l^{\prime}})=S_{\text{n}}(\Lambda_{l},g[ \Lambda_{l^{\prime}}])=S_{\text{n}}(g[\Lambda_{l}],\Lambda_{l^{\prime}}) \tag{14}\]

for all such gauge transformations \(g\in\mathcal{G}\). A general way to guarantee Eq. (14) proceeds by replacing,

\[S_{\text{n}}(\Lambda_{l},\Lambda_{l^{\prime}})\quad\longrightarrow\quad\max_{g, g^{\prime}\in\mathcal{G}}S_{\text{n}}(g[\Lambda_{l}],g^{\prime}[\Lambda_{l^{ \prime}}]). \tag{15}\]

However, in practice, it might not be required to iterate over all possible gauge transformations in \(\mathcal{G}\) due to the locality of the similarity measure. In the following, we will use the toric code and a specific RBM variational ansatz as an example to illustrate these gauge transformations and how an appropriate function \(f\) in Eq. (11) and gauge invariance (14) can be implemented efficiently.

Finally, note that, while we focus on applying DM in this work, a similarity measure in terms of neural network parameters can also be used for other kernel techniques such as kernel PCA. Depending on the structure of the underlying dataset, DM has clear advantage over kernel PCA: the former really captures the global connectivity of the dataset rather than the subspace with most variance that is extracted by the latter. This is why kernel PCA fails when identifying, e.g., winding numbers, in general datasets where DM still works well [12]. Specifically for our case study of the toric code below, we find that kernel PCA can also identify topological sectors for small \(T\) and without magnetic field, \(h=0\), as a result of the simple data structure; however, only DM works well when \(h\) is turned on, as we discuss below.

## III Example: Toric Code

Now we illustrate our DM-based ML algorithm using the toric code model [57], defined on an \(L_{x}\times L_{y}\) square lattice with spin-\(1/2\) operators or qubits on every bond, see Fig. 2(b), leading to a total of \(N=2L_{x}L_{y}\) spins; throughout this work, we will assume periodic boundary conditions. Referring to all four spins on the edges of an elementary square (vertex) of the lattice as plaquette \(P\) (star \(S\)), the plaquette and star operators are defined
as \(\hat{\mathcal{P}}_{P}=\prod_{i\in P}\hat{s}_{i}^{x}\) and \(\hat{\mathcal{S}}_{S}=\prod_{i\in S}\hat{s}_{i}^{x}\), respectively. The toric code Hamiltonian then reads as

\[\tilde{H}_{\rm tc}=-J_{P}\sum_{P}\hat{\mathcal{P}}_{P}-J_{S}\sum_{S}\hat{ \mathcal{S}}_{S}, \tag{16}\]

where the sums are over all plaquettes and stars of the lattice. All "stabilizers" \(\hat{\mathcal{P}}_{P}\), \(\hat{\mathcal{S}}_{S}\) commute among each other and with the Hamiltonian. Focusing on \(J_{P},J_{S}>0\), the ground states are obtained as the eigenstates with eigenvalue \(+1\) under all stabilizers. A counting argument, taking into account the constraint \(\prod_{S}\hat{\mathcal{S}}_{S}\!=\!\prod_{P}\hat{\mathcal{P}}_{P}\!=\!\mathds{1}\), reveals that there are four, exactly degenerate ground states for periodic boundary conditions.

To describe the ground-states and low-energy subspace of the toric code model (16) variationally, we parameterize \(\psi(\mathbf{\sigma};\,\Lambda)\) in Eq. (1) using the ansatz

\[\psi_{\rm rbm}(\mathbf{\sigma};\,\Lambda)= \prod_{P}\cos(b_{P}+\sum_{j\in P}w_{Pj}\sigma_{j})\] \[\times \prod_{S}\cos(b_{S}+\sum_{j\in S}w_{Sj}\sigma_{j}), \tag{17}\]

proposed in Ref. [58], where every plaquette \(P\) (star \(S\)) is associated with a "bias" \(b_{P}\) (\(b_{S}\)) and four weights \(w_{P,j}\) (\(w_{S,j}\)), all of which are chosen to be real here, i.e., \(\Lambda=\{b_{P},b_{S},w_{P,j},w_{S,j}\}\). This ansatz can be thought of as an RBM [46] (see Appendix A), as illustrated in Fig. 2(c), with the same geometric properties as the underlying toric code model. It is clear that Eq. (17) defines a quasi-local ansatz as it is of the form of Eq. (3), with j enumerating all plaquettes and stars (and thus \(N_{j}=2N\)). For this specific ansatz, the gauge transformations \(g\in\mathcal{G}\), as introduced in Sec. II.3 above, are generated by the following set of operations on the parameters \(b_{P}\), \(b_{S}\), \(w_{P,j}\), and \(w_{S,j}\):

1. For \(X\) being any plaquette or star, multiplying all biases and weights of that plaquette or star by \(-1\) [see Fig. 3(a)], \[g_{X,-}:\,b_{X}\to-b_{X},\ w_{Xj}\to-w_{Xj},\] (18a) leaves the wave function invariant [\(\vartheta_{g}=0\) in Eq. (13)].
2. Adding \(\pi\) to either the bias or any of the weights associated with the plaquette or star \(X\) [see Fig. 3(b)], \[g_{X,\pi,b}: \!b_{X}\to b_{X}+\pi,\] (18b) \[g_{X,\pi,j}: \!w_{Xj}\to w_{Xj}+\pi,\quad j\in X,\] (18c) leads to an overall minus sign [\(\vartheta_{g}=\pi\) in Eq. (13)].
3. For any closed loop \(\ell\) (or \(\bar{\ell}\)) on the direct (or dual lattice), adding \(\frac{\pi}{2}\) to all weights of the stars (plaquettes) that are connected to the spins crossed by the string [see Fig. 3(c-d)], \[g_{\ell}:\,w_{Sj}\to w_{Sj}+\frac{\pi}{2},\quad Sj\in\ell,\] (18d) \[g_{\bar{\ell}}:\,w_{Pj}\to w_{Pj}+\frac{\pi}{2},\quad Pj\in\bar{ \ell},\] (18e) leads to \(\vartheta_{g}=0\) or \(\pi\) in Eq. (13) depending on the length of the string. Note that any loop configuration \(\mathcal{L}\), which can contain an arbitrary number of loops, can be generated by the set \(\{g_{S},g_{P},g_{x,y},g_{\bar{x},\bar{y}}\}\), where \(g_{S}\) (\(g_{P}\)) creates an elementary loop on the dual (direct) lattice encircling the star \(S\) (plaquette \(P\)), see Fig. 3(c,d), and \(g_{x,y}\) (\(g_{\bar{x},\bar{y}}\)) creates a non-contractible loop on the direct (dual) lattice along the \(x,y\) direction. Since the length of any contractible loop is even, \(\vartheta_{g}=0\) for any string transformations generated by \(g_{S}\) and \(g_{P}\). Meanwhile, on an odd lattice, the gauge transformations \(g_{x,y}(g_{\bar{x},\bar{y}})\) involve an odd number of sites and thus lead to \(\vartheta_{g}=\pi\).

Figure 3: Gauge freedom of RBM ansatz in Eq. (17). The following transformations only lead to a global phase: (a) Multiplying all the parameters of a plaquette (or star, not shown) by a minus sign, see Eq. (18a); (b) A \(\pi\) shift of a single parameter, see Eqs. (18b) and (18c); (c) A \(\pi/2\) shift to the weights crossed by a string \(\bar{l}\), defined by \(g_{\bar{l}}\) in Eq. (18e). The straight pink line represents the transformation on a non-contractible loop denoted by \(g_{y}\); (d) Same as (c) but for loops on the direct lattice and \(g_{\bar{l}}\) and \(g_{\bar{y}}\), cf. Eq. (18d).
A highly inefficient way of dealing with this gauge redundancy would be to use a choice of \(S_{n}\) in Eq. (11) which is not invariant under any of the transformations in Eq. (18); this would, for instance, be the case by just taking the Euclidean distance of the weights,

\[\begin{split} S_{\text{eu}}&(\Lambda_{l},\Lambda_{l^ {\prime}})\propto||\Lambda_{l}-\Lambda_{l^{\prime}}||^{2}\\ &=\sum_{X}\Bigl{[}(b_{X}^{l}-b_{X}^{l^{\prime}})^{2}+\sum_{j \in X}(w_{Xj}^{l}-w_{Xj}^{l^{\prime}})^{2}\Bigr{]},\end{split}\]

where the sum over \(X\) involves all plaquettes and stars. Naively going through all possible gauge transformations to find the maximum in Eq. (15) would in principle rectify the lack of gauge invariance. However, since the number of gauge transformations scales exponentially with system size \(N\) (holds for each of the three classes, 1.-3., of transformations defined above), such an approach would become very expensive for large \(N\). Luckily, locality of the ansatz and of the similarity measure allows us to construct similarity measures that can be evaluated much faster: as an example, consider

\[\begin{split} S_{n}&(\Lambda_{l},\Lambda_{l^{ \prime}})=\frac{1}{2}+\frac{1}{10N}\sum_{X}\max_{\tau_{X}=\pm}\Bigl{[}\\ &\sum_{j\in X}\cos 2(\tau_{X}w_{Xj}^{l}-w_{Xj}^{l^{\prime}})+ \cos 2(\tau_{X}b_{X}^{l}-b_{X}^{l^{\prime}})\Bigr{]},\end{split} \tag{19}\]

which clearly obeys \(S_{n}(\Lambda_{l},\Lambda_{l^{\prime}})=S_{n}(\Lambda_{l^{\prime}},\Lambda_{l})\), \(S_{n}(\Lambda_{l},\Lambda_{l^{\prime}})\in[0,1]\), and locality [it is of the form of Eq. (11) with j enumerating all \(X\)]. Concerning gauge invariance, first note that the choice of \(\cos(\cdot)\) immediately leads to invariance under Eq. (18a). Second, for each \(X\) we only have to maximize over two values (\(\tau_{X}\)) to enforce invariance under Eqs. (18b) and (18c), i.e., the maximization only doubles the computational cost.

The "string" redundancy, see Eqs. (18d) and (18e), however, is not yet taken into account in Eq. (19). It can be formally taken care of by maximizing over all possible loop configurations, denoted by \(\mathcal{L}\),

\[\begin{split} S_{\text{str}}&(\Lambda_{l},\Lambda_{ l^{\prime}})=\frac{1}{2}+\frac{1}{10N}\max_{\mathcal{L}}\Bigl{\{}\sum_{X} \max_{\tau_{X}=\pm}\Bigl{[}\\ &\sum_{j\in X}\mu_{Xj}^{\mathcal{L}}\cos 2(\tau_{X}w_{Xj}^{l}-w_{ Xj}^{l^{\prime}})+\cos 2(\tau_{X}b_{X}^{l}-b_{X}^{l^{\prime}})\Bigr{]}\Bigr{\}}, \end{split} \tag{20}\]

where \(\mu_{Xj}^{\mathcal{L}}=-1\) if \(Xj\) lives on a loop contained in \(\mathcal{L}\) and \(\mu_{Xj}^{\mathcal{L}}=1\) otherwise. While there is an exponential number of such strings, Ref. [12] has proposed an algorithm to efficiently find an approximate maximum value. In our case, this algorithm amounts to randomly choosing a plaquette \(P\) or a star \(S\) or a direction \(d=x,y\) and then applying \(g_{S}\) or \(g_{P}\) or \(g_{d=x,y}\) to \(\Lambda_{l}\) in Eq. (19). If this does not decrease the similarity, keep that transformation; if it decreases the similarity, discard the gauge transformation. Repeat this procedure \(N_{g}\) times. In Ref. [12], \(N_{g}\) between \(10^{3}\) and \(10^{4}\) was found to be enough for a large system consisting \(18\times 18\) square-lattice sites (total of \(N=2\times 18^{2}\) qubits). On top of this, \(g_{S}\) and \(g_{P}\) are local and, hence, the evaluation of the change of the similarity with the gauge transformation only requires \(\mathcal{O}(N^{0})\) amount of work.

In the numerical simulations below, using Eq. (19) without sampling over loop configurations \(\mathcal{L}\) turned out to be sufficient. The reason is that, for our Markov-chain-inspired sampling procedure of \(\Lambda_{l}\) (see Appendix C), updates that correspond to these loop transformations happen very infrequently. Furthermore, even if a few pairs of samples are incorrectly classified as distinct due to the string redundancy, the DM will still correctly capture the global connectivity and, hence, absence or presence of topological sectors.

Figure 4: (a) DM spectrum for topological phase at \(h=0\) and \(T=0.1\) using the neutral network similarity measure in Eq. (19). Inset left: associated leading DM components; color represents the loop observable expectations values defined in (c-d). Inset right: DM spectrum in descending order at \(\epsilon=0.01\) indicated by the dashed line. (b) Same as (a), but using exact overlaps \(S_{\text{q}}\) in Eq. (10) as metric. (c) Color map for the non-local loop values \(\langle\overline{W}_{1}\rangle,\langle\overline{W}_{2}\rangle\) in the left insets of (a) and (b). (d) Different straight Wilson loops \(\hat{W}_{1,x_{i}}\) (\(\hat{W}_{2,g_{i}}\)) along \(x\) (\(y\)) direction, denoted by blue (red) lines. The loop values in the color map in (c) are spatial averages over all straight-loop expectation values (as in the equations for \(\langle\overline{W}_{1}\rangle,\langle\overline{W}_{2}\rangle\)).
Numerical results

We next demonstrate explicitly how the general procedure outlined above can be used to probe and analyze topological order in the toric code. We start from the pure toric code Hamiltonian defined in Eq. (16) using the variational RBM ansatz in Eq. (17). An ensemble of network parameters is generated by applying the procedure of Sec. II.1 (see also Algorithm 1) for a system size of \(N\!=\!18\) spins; the hyperparameters for ensemble generation and more details including the form of \(u\) in Eq. (4) are given in Appendix C. From now on, we measure all energies in units of \(J_{P}\) and set \(J_{S}=J_{P}=1\).

Let us first focus on the low-energy ensemble and choose \(T=0.1\) in Eq. (5). For the simple similarity measure in Eq. (19), that can be exactly evaluated at a time linear in system size \(N\), we find the DM spectrum shown in Fig. 4(a) as a function of \(\epsilon\) in Eq. (9). We observe the hallmark feature of four superselection sectors [12]: there is a finite range of \(\epsilon\) where there are four eigenvalues exponentially close to \(1\). The association of samples (in our case states) and these four sectors is thus expected to be visible in a scatter plot of a projected subspace spanned by the first three non-trivial eigenvectors \(\psi_{1,2,3}\)[12]; note the zeroth eigenvector \((\psi_{0})_{l}=C\) is always constant with eigenvalue \(\lambda=1\) from probability conservation. In fact, we can see these clusters already in the first two components, see left inset in Fig. 4(a). Then a standard \(k\)-means algorithm is applied onto this projected subspace to identify the cluster number for each data point. To verify that the ML algorithm has correctly clustered the states according to the four physical sectors, we compute the expectation value for each state of the string operators,

\[\hat{W}_{1,\bar{x}}=\prod_{i\in\bar{x}}\hat{s}_{i}^{x},\quad\hat{W}_{2,\bar{y} }=\prod_{i\in\bar{y}}\hat{s}_{i}^{x}, \tag{21}\]

where \(\bar{x}(\bar{y})\) are loops defined on the dual lattice winding along the \(x(y)\) direction, shown as blue lines in Fig. 2(b). We quantify the association of a state to physical sectors by the average of a set of straight loops \(\mathcal{X}(\mathcal{Y})\) winding around the \(x(y)\) direction, shown as blue (red) lines in Fig. 4(d). Indicating this averaged expectation value \(\langle\overline{W}_{1}\rangle,\langle\overline{W}_{2}\rangle\) in the inset of Fig. 4(a) using the color code defined in Fig. 4(c), we indeed see that the clustering is done correctly.

To demonstrate that this is not a special feature of the similarity measure in Eq. (19), we have done the same analysis, with result shown in Fig. 4(b), using the full quantum mechanical overlap measure in Eq. (10). Quantitative details change but, as expected, four superselection sectors are clearly identified and the clustering is done correctly. We reiterate that the evaluation of the neural-network similarity measure in Eq. (19) [exact evaluation \(\mathcal{O}(N)\)] is much fast than that in Eq. (10) [exact evaluation \(\mathcal{O}(2^{N})\), but we can compute it approximately with importance sampling] on a classical computer. Note, however, that once \(S_{n}\) is computed for all samples, the actual DM-based clustering takes the same amount of computational time for both approaches. Consequently, suppose there is a quantum simulator that can efficiently measure the quantum overlap in Eq. (10) or any other viable similarity measure for that matter, then we can equivalently use the "measured" similarity for an efficient clustering of the superselection sectors via the DM scheme. As a next step, we demonstrate that the superselection sectors are eventually connected if we take into account states with sufficiently high energy. To this end, we repeat the same analysis but for an ensemble with \(T=1\). As can be seen in the resulting DM spectrum in Fig. 5(a), there is no value of \(\epsilon\) where more than one eigenvalue is (exponentially) close to \(1\) and separated from the rest of the spectrum by a clear gap. Here we used again the simplified measure in Eq. (19), but have checked nothing changes qualitatively when using the overlap measure. To verify that this is the correct answer for the given dataset, we again computed the expectation value of the loop operators in Eq. (21) for each state in the ensemble. This is shown in Fig. 5(b), where we also use color to indicate the energy expectation value for each state. We can clearly see the four low-energy (blue) sectors (with \(|W_{1,2}|\simeq 1\)) are connected via high-energy (red) states (with \(|W_{1,2}|\ll 1\)). This agrees with the DM result that

Figure 5: (a) DM spectrum for the high-energy ensemble at \(h\!=\!0\) and \(T\!=\!1\). The inset is the spectrum at \(\epsilon=0.03\) indicated by the dashed line in the main panel; (b) Spatially averaged straight Wilson loops \(\langle\overline{W}_{1(2)}\rangle\) [see Fig. 4(c-d)] along two directions for the states in (a), where the color encodes energy density \(\langle H\rangle/N\); (c) Leading DM components where the color of the dots encodes \(\langle\overline{W}_{1(2)}\rangle\) using the color map in Fig. 4(d); (d) DM spectrum for the trivial phase at \(h\!=\!1.0\) and \(T\!=\!0.1\) using the quantum metric \(S_{\rm q}\).
all states are connected within the ensemble (topological order is lost). We can nonetheless investigate the clustering in the leading three non-trivial DM components \(\psi_{1,2,3}\). Focusing on a 2D projection in Fig. 5(c) for simplicity of the presentation, we can see that the DM reveals very interesting structure in the data: the four lobes roughly correspond to the four colors blue, red, orange, and green associated with the four superselection sectors and the states closer to \(|W_{1,2}|=1\) (darker color) appear closer to the tips. Finally, note that the colors are arranged such that the red and green [orange and blue] lobes are on opposite ends, as expected since they correspond to \((W_{1},W_{2})\simeq(1,-1)\) and \((-1,1)\) [\((-1,-1)\) and \((1,1)\)].

Another route to destroying topological order proceeds via application of a magnetic field. To study this, we extend the toric code Hamiltonian according to

\[\hat{H}^{\prime}_{\rm tc}=\hat{H}_{\rm tc}-h\sum_{i}\hat{s}_{i}^{z}. \tag{22}\]

Clearly, in the limit of \(h\to\infty\), the ground state is just a state where all spins are polarized along \(\hat{s}^{z}\) and topological order is lost. Starting from the pure toric model (\(h=0\)) and turning on \(h\) reduces the gap of the "charge excitations" defined by flipping \(\hat{\mathcal{S}}_{S}\) from \(+1\) in the toric code groundstate to \(-1\). Their condensation leads to a second-order quantum phase transition [69, 70, 71, 72].

Before addressing the transition, let us study the large-\(h\) limit. We first note that our ansatz in Eq. (17) does not need to be changed as it can capture the polarized phase as well. For instance, denoting the "northmost" (and "southmost") spin of the plaquette \(P\) (and star \(S\)) by \(j_{0}(P)\) (and \(j_{0}(S)\)), respectively, the spin polarized state is realized for [see also Fig. 8(a) in the Appendix]

\[b_{P}=b_{S}=-\frac{\pi}{4},\ \ w_{Xj}=\begin{cases}\frac{\pi}{4},&j=j_{0}(X), \\ 0,&\text{otherwise}.\end{cases} \tag{23}\]

In fact, the spin polarized state has many representations within our RBM ansatz in Eq. (17), including representations that are not just related by the gauge transformations in Eq. (18). For instance, the association \(j\to j_{0}(X)\) of a spin to a plaquette and star can be changed, e.g., by using the "easternmost" spin. As discussed in more detail in Appendix A.2, this redundancy is a consequence of the product from of \(\psi_{\rm rbm}(\mathbf{\sigma})\) in Eq. (17) and the fact that \(\psi_{\rm rbm}(\mathbf{\sigma})\) is _exactly_ zero if there is a single \(j\) with \(\sigma_{j}=-1\); consequently, it is a special feature of the simple product nature of the spin-polarized ground state. While in general there can still be additional redundancies besides the aforementioned gauge transformations, we do not expect such a structured set of redundancy to hold for generic states. There are various ways of resolving this issue. The most straightforward one is to replace the simple overlap measure \(S_{\rm n}\) in Eq. (11) by the direct overlap \(S_{\rm q}\) in Eq. (10) for a certain fraction of pairs of samples \(l\) and \(l^{\prime}\). If this fraction is large enough, the DM algorithm will be able recognize that clusters of network parameters that might be distinct according to \(S_{\rm n}\) actually correspond to identical wave functions. We refer to Appendix A.3 where this is explicitly demonstrated. We note, however, that kernel PCA will not work anymore in this case; it will incorrectly classify connected samples as distinct as it's based on the variance of the data rather than connectivity. For simplicity of the presentation, we use \(S_{\rm q}\) for all states in the main text and focus on DM.

Figure 6: DM spectra for low-energy ensembles with \(T\!=\!0.3\) at finite field \(h\). (a) First 10 eigenvalues for various field values \(h\!=\!0.475,0.55,0.575,0.6,0.7\) at \(\epsilon\!=\!0.05\). The dot marker (\(h\!=\!0.475\)) shows that the eigenvalue spectra have four-fold degeneracy, indicating signature for topological order. In comparison, for spectra marked by the the triangular markers (\(h\geq 0.55\)), such degeneracy is absent. A transition field value \(h_{t}\simeq 0.55\) is identified by observing that a gap opens in the degenerate eigenvalue spectra. This is consistent with what we have observed in the fidelity using the same dataset [see Appendix B.1]. (b) Projected eigenvectors onto the first two components for \(h\!=\!0.475\). The color encodes \(\langle\overline{W}_{1(2)}\rangle\) with the color scheme of Fig. 4(c). The black cross marks the \(k\)-means centers. (c) Same as (b) for \(h=0.7\). (d) Expectation for averaged straight Wilson loops \(\langle\overline{W}_{1(2)}\rangle\) along two directions for the states in (b). The color encodes the clustering results from \(k\)-means in the projected subspace of the eigenvectors shown in (b). (e) Same as (d) for ensemble shown in (c).
The DM spectrum for large magnetic field, \(h=1\), and low temperatures, \(T=0.1\), is shown in Fig. 5(d). Clearly, there is no value of \(\epsilon\) for which there is more than one eigenvalue close to \(1\) while exhibiting a gap to the rest of the spectrum. This shows that, as expected, the magnetic field \(h\) has lead to the loss of topological order.

To study with our DM algorithm the associated phase transition induced by \(h\), we repeat the same procedure for various different values of \(h\). The resulting spectra for selected \(h\) are shown in Fig. 6(a). We see that there are still four sectors for \(h=0.55\) in the data that are absent for \(h=0.575\) and larger values. While the associated critical value of \(h\) is larger than expected [69; 70; 71], this is not a shortcoming of the DM algorithm but rather a consequence of our simple local variational ansatz in Eq. (17). By computing the fidelity as well as loop-operator expectation values, we can see that a critical value around \(h=0.55\) is the expected answer for our dataset (see Appendix B.1). More sophisticated ansatze for the wavefunction are expected to yield better values, but this is not the main focus of this work. More importantly, we see in Fig. 6(b) that the DM clustering of the states correctly reproduces the clustering according to the averaged loop operator expectation values \(\langle\overline{W}_{j}\rangle\) (again indicated with color). Alternatively, this can be seen in Fig. 6(d) where \(\langle\overline{W}_{j}\rangle\) is indicated for the individual samples. Using four different colors for the four different clusters identified by the DM, we see that all states are clustered correctly. As expected based on the eigenvalues, there are no clear clusters anymore for larger \(h\), Fig. 6(c); nonetheless, naively applying \(k\)-means clustering in \(\psi_{1,2,3}\) manages to discover some residual structure of the wavefunctions related to \(\langle\overline{W}_{j}\rangle\) as demonstrated in Fig. 6(e).

## V Summary and discussion

In this work, we have described an unsupervised ML algorithm for quantum phases with topological order. We use neural network parameters to efficiently represent an ensemble of quantum states, which are sampled according to their energy expectation values. To uncover the structure of the superselection sectors in the quantum states, we used the dimensional reduction technique of diffusion map and provided a kernel defined in terms of network parameters. As opposed to a kernel based on the overlap of wavefunctions (or other quantum mechanical similarity measures of states for that matter), this metric can be evaluated efficiently (within polynomial time) on a classical computer.

We illustrated our general algorithm using a quasi-local restricted Boltzmann machine (RBM) and the toric code model in an external field; the choice of network ansatz was inspired by previous works [58; 59] showing the existence of efficient representations of the low-energy spectrum in terms of RBMs. Allowing for spatially inhomogeneous RBM networks, we identified the "gauge symmetries" of the ansatz, i.e., the set of changes in the network parameters that do not change the wavefunction, apart from trivial global phase factors. We carefully designed a similarity measure that is gauge invariant--a key property as, otherwise, identical wavefunctions represented in different gauges would be falsely identified as being distinct. We showed that the resultant unsupervised diffusion-map-based embedding of the wavefunctions is consistent with the expectation values of loop operators; it correctly captures the presence of superselection sectors and topological order at low energies and fields, as well as the lack thereof when higher-energy states are involved and/or the magnetic field is increased. We also verified our results using the full quantum mechanical overlap of wavefunctions as similarity measure.

On a more general level, our analysis highlights the importance of the following two key properties of diffusion maps: first, in the presence of different topological sectors, the leading eigenvectors of diffusion maps capture the connectivity rather than, e.g., the variance as is the case for PCA. For this reason, the clustering is still done correctly even if a fraction of pairs of wavefunctions are incorrectly classified as being distinct due to the usage of an approximate similarity measure. This is why complementing the neural-network similarity measure, which has additional, state-specific redundancies in the large-field limit, by direct quantum mechanical overlaps for a certain fraction of pairs of states is sufficient to yield the correct classification. The second key property is that diffusion map is a kernel technique. This means that the actual machine learning procedure does not require the full wavefunctions as input; instead, only (some measure of) the kernel of all pairs of wavefunctions in the dataset is required. We have used this to effectively remove the gauge redundancy in the RBM parametrization of the states by proper definition of the network similarity measure in Eq. (20). Since the evaluation of full quantum mechanical similarity measures, like the wavefunction overlap, are very expensive on classical computers, an interesting future direction would be to use the emerging quantum-computing resources to evaluate a similarity measure quantum mechanically. This could then be used as input for a diffusion-map-based clustering.

We finally point out that the ensemble of states we used in this work, which was based on sampling states according to their energy with respect to a Hamiltonian, is only one of many possibilities. The proposed technique of applying diffusion map clustering using a gauge-invariant kernel in terms of network parameters of a variational description of quantum many-body wavefunctions can be applied more generally, in principle, to any ensemble of interest. For instance, to consider arbitrary local perturbations, one could generate an ensemble using finite depth local unitary circuits. Alternatively, one could generate an ensemble based on (Lindbladian) time-evolution to probe the stability of topological order against time-dependent perturbations or the coupling to a bath. We leave the investigation of such possibilities for future works.
Code and data availability

The Monte Carlo simulations in this work were implemented in JAX [73]. Python code and data will be available at [https://github.com/teng10/ml_toric_code/](https://github.com/teng10/ml_toric_code/).

## Acknowledgements

Y.T. acknowledges useful discussions with Dmitrii Kochkov, Juan Carrasquilla, Khadijeh Sona Najafi, Maine Christos and Rhine Samajdar. Y.T. and S.S. acknowledge funding by the U.S. Department of Energy under Grant DE-SC0019030. M.S.S. thanks Joaquin F. Rodriguez-Nieva for a previous collaboration on DM [12]. The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.

## Appendix A Variational Ansatz: Restricted Boltzmann Machine

The variational ansatz in Eq. (17) is a _further-restricted_ restricted Boltzmann machine (RBM), first introduced by Ref. [58]. RBM is a restricted class of Boltzmann machine with an "energy" function \(E_{\text{RBM}}(\boldsymbol{\sigma},\boldsymbol{h};\Lambda)\) dependent on the network parameters \(\Lambda\), where \(\boldsymbol{\sigma}\) are physical spins and \(\boldsymbol{h}=\{h_{1},h_{2},\cdots,h_{N}\mid h_{i}=\pm 1\}\) are hidden spins (or hidden neurons) that are Ising variables. The parameters \(\Lambda\) define the coupling strength among the physical and hidden spins. The restriction in RBM is that the couplings are only between the physical spin \(\sigma_{i}\) and hidden spin \(h_{j}\) with strength \(-w_{ij}\), so that the "energy" function takes the form \(E_{\text{RBM}}(\boldsymbol{\sigma},\boldsymbol{h};\Lambda)\!=\!-\sum_{i}a_{i} \sigma_{i}-\sum_{i}b_{i}h_{i}-\sum_{ij}w_{ij}\sigma_{i}h_{j}\). It is a generative neural network that aims to model a probability distribution \(\mathbb{P}\) based on the Boltzmann factor,

\[\mathbb{P}(\boldsymbol{\sigma};\Lambda) =\frac{1}{Z}\sum_{\boldsymbol{h}}e^{-E_{\text{RBM}}(\boldsymbol{ \sigma},\boldsymbol{h};\Lambda)}, \tag{18a}\] \[\text{normalization}\quad Z =\sum_{\boldsymbol{\sigma},\boldsymbol{h}}e^{-E_{\text{RBM}}( \boldsymbol{\sigma},\boldsymbol{h};\Lambda)}. \tag{18b}\]

For the task of modeling a quantum wavefunction amplitude \(\psi(\boldsymbol{\sigma};\Lambda)\), RBMs can be used as a variational ansatz by extending the parameters \(\Lambda\) to complex numbers.

Further restricting parameters to the interlayer connections to the plaquette and star geometry in the toric code model [cf. Fig. 2(c)] and taking all parameters \(\Lambda\) to be purely imaginary, we recover the ansatz in Eq. (17) (up to normalization factor \(\widetilde{Z}\)),

\[\psi(\boldsymbol{\sigma};\Lambda) =\frac{1}{\widetilde{Z}}\sum_{X=P,S}\sum_{h_{X}=\pm 1}e^{-i \sum_{X}(w_{Xj}\sigma_{j}+b_{X})h_{X}},\] \[=\frac{1}{\widetilde{Z}}\prod_{X=P,S}\cos(\sum_{j\in X}w_{Xj} \sigma_{j}+b_{X}). \tag{19}\]

The \(\cos(\cdot)\) factors come from summing over the hidden neurons and the ansatz factorizes into the product of individual plaquette (star) terms because of the restricted connections. The estimation of physical observables of a wave function based on the RBM ansatz requires Monte Carlo sampling procedure which we discuss in Appendix B.

### Ground states representation in different topological sectors

Placing the toric code model in Eq. (16) on the torus geometry, it is useful to define the loop operators,

\[\hat{W}_{1} =\prod_{i\in\bar{l}_{x}}\hat{s}_{i}^{x},\quad\hat{W}_{2}=\prod_{i \in\bar{l}_{y}}\hat{s}_{i}^{x}, \tag{20a}\] \[\hat{V}_{1} =\prod_{i\in l_{x}}\hat{s}_{i}^{z},\quad\hat{V}_{2}=\prod_{i\in l _{y}}\hat{s}_{i}^{z}, \tag{20b}\]

where \(l_{x,y}\) is a non-contractible loop along \(x\), \(y\) direction, and \(\bar{l}_{x,y}\) is similar on the dual lattice. Note the loop operators along two directions do not commute with each other as \(\left[\hat{W}_{1},\hat{V}_{2}\right]\neq 0\) and \(\left[\hat{W}_{2},\hat{V}_{1}\right]\neq 0\). However, since the hamiltonian commute with these loop operators \(\left[\hat{W}_{1,2},\hat{H}_{\text{tc}}\right]\!=\!\left[\hat{V}_{1,2},\hat{H }_{\text{tc}}\right]\!=\!0\), it follows that the ground state subspace is four-fold degenerate and spanned by the eigenvectors of the loop operators.

Suppose we work in the eigenbasis of \(\hat{W}_{1,2}\); we define the four orthogonal ground states \(\left|\psi_{i}\right\rangle(i=0,1,2,3)\) that

Figure 7: RBM representations of the four toric code ground states in the eigenbasis [Eq. (18)] of loop operators \(\hat{W}_{1},\hat{W}_{2}\) in Eq. (20a).
span \(\mathcal{L}\) as,

\[\hat{W}_{1}\ket{\psi_{0}}=\ket{\psi_{0}},\quad\hat{W}_{2}\ket{\psi_{0} }=\ket{\psi_{0}}, \tag{10a}\] \[\hat{W}_{1}\ket{\psi_{1}}=\ket{\psi_{1}},\quad\hat{W}_{2}\ket{\psi_{1 }}=-\ket{\psi_{1}},\] (10b) \[\hat{W}_{1}\ket{\psi_{2}}=\ket{\psi_{2}},\quad\hat{W}_{2}\ket{\psi_{2 }}=-\ket{\psi_{2}},\] (10c) \[\hat{W}_{1}\ket{\psi_{3}}=-\ket{\psi_{3}},\quad\hat{W}_{2}\ket{\psi_ {3}}=-\ket{\psi_{3}}. \tag{10d}\]

The RBM ansatz in Eq. (11) can represent eigenstates of \(\hat{W}_{1,2}\) with eigenvalues \((W_{1},W_{2})=(\pm 1,\pm 1)\). Ref. [58] gave an representation of \(\ket{\psi_{3}}\) with parameters,

\[w_{Pj}=\frac{\pi}{4},\quad b_{P}=0,\quad w_{Sj}=\frac{\pi}{2},\quad b_{S}=0.\] (11a) On a system with odd number of sites along \[x\] and \[y\] direction, the other three degenerate states can be realized analogously by fixing the weights associated to stars to be \[w_{Sj}\!=\!0,b_{S}\!=\!0\]. Then the four states can be chosen by changing the \[w_{Pj}\] and \[b_{P}\] as shown in Fig. 7.

### Network parameter redundancies in polarized phase

In Sec. III, we identified a set of gauge transformations Eq. (18) that leave a generic wavefunction parameterized by the RBM ansatz in Eq. (17) invariant up to a global phase [Eq. (13)]. Such gauge transformations should be taken into consideration when evaluating the similarity measure \(S_{n}\). Moreover, we have numerically verified that for states generated close to the exact toric code wave functions, \(S_{n}\) is a good proxy for the quantum measure \(S_{\rm q}\) after explicit removals of such redundancies via \(S_{n}\) in Eq. (19). However, as alluded to in the discussions of the large-\(h\) limit, there are state-specific redundancies that are generally not related by the gauge transformations in Eq. (18).

Let us illustrate such redundancies here for the polarized state \(\ket{\Psi}=\ket{1,\cdots,1}_{z}\) which has all spin pointing up in the \(z\)-basis. Notice that there is the same number of \(\cos(\cdot)\) factors in the wavefunction ansatz as the number of spins. As a result, we can define a "_covering_" by assigning each individual spin to a single factor, and choosing the weights to ensure all spins are pointing up. Any such "covering" is a valid representation of the polarized state. For example, one representation is given by,

\[b_{P}=b_{S}=-\frac{\pi}{4},\ \ w_{Sj}=\begin{cases}\frac{\pi}{4},&j=j_{s}(S),\\ 0,&\text{otherwise},\end{cases}\quad\text{and}\ \ \ w_{\rm Pj}=\begin{cases}\frac{\pi}{4},&j=j_{n}(P),\\ 0,&\text{otherwise}.\end{cases} \tag{11b}\]

where \(j_{s}(S)\) denotes the "southmost" spin in the star \(S\) and \(j_{n}(P)\) denotes the "northmost" spin in the plaquette \(P\) [see Fig. 8(a)]. Any such coverings of the spins will correspond to a polarized state. For example, performing a "rotation" leads to a different covering in Fig. 8(b). Actually, because most amplitudes in local-\(z\) basis are \(0\) so there are so few constraints in the wave function amplitudes, a continuous set of weights exist to represent the polarized state, so there are an infinite amount of redundancies for completely polarized state.

To illustrate this, let us consider the simplest example of just two spins [the boxed region in Fig. 8(c)] with the same RBM ansatz, which can be easily generalized to more spins. For two spins, such ansatz is given by,

\[\psi_{\Lambda}(\sigma_{A},\sigma_{B})=\cos(b_{S}+w_{SA}\sigma_{A}+w_{SB} \sigma_{B})\cos(b_{P}+w_{PA}\sigma_{A}+w_{PB}\sigma_{B}), \tag{12}\]

where the weights \(\Lambda=\{\Lambda_{S}=\{b_{S},w_{SA},w_{SB}\},\,\Lambda_{P}=\{b_{P},w_{PA},w_{PB }\}\}\) with \(\Lambda_{Xj}\in[0,\pi)\) for \(X=S\) or \(P\) fully determine the two-qubits physical state. For example, the following two choices of weights [\(\Lambda_{1}\) and \(\Lambda_{2}\) pictorially in

Figure 8: (a-b) Two RBM representations Eq. (10) of the polarized state. (c) A path that connects the presentation for two spins in (a-b), which is explicitly shown in Table. 1.
Fig. 8(c)] both parametrize the polarized state:

\[\Lambda_{1}=\{b_{S}=-\frac{\pi}{4},w_{SA}=0,w_{SB}=\frac{\pi}{4},b_{P} =-\frac{\pi}{4},w_{PA}=\frac{\pi}{4},w_{PB}=0\}, \tag{10a}\] \[\Lambda_{2}=\{b_{S}=-\frac{\pi}{4},w_{SA}=\frac{\pi}{4},w_{SB}=0,b _{P}=-\frac{\pi}{4},w_{PA}=0,w_{PB}=\frac{\pi}{4}\},\] (10b) \[\psi_{\Lambda_{1,2}}=\begin{cases}1,&\sigma_{A}=\sigma_{B}=1,\\ 0,&\text{otherwise}.\end{cases} \tag{10c}\]

Now to illustrate the continuous redundancies, we construct a path in the parameter space to go from \(\Lambda_{1}\) to \(\Lambda_{2}\). The path is composed of three steps [Fig. 8(c)],

\[\Lambda_{1}\xrightarrow{\text{path}\,1}\Lambda_{3}\xrightarrow{\text{path}\,2 }\Lambda_{4}\xrightarrow{\text{path}\,3}\Lambda_{2}, \tag{11}\]

where the intermediate parameters are given by,

\[\Lambda_{3}=\{b_{S}=0,w_{SA}=\frac{\pi}{4},w_{SB}=-\frac{\pi}{4}, b_{P}=-\frac{\pi}{4},w_{PA}=\frac{\pi}{4},w_{PB}=0\}, \tag{12}\] \[\Lambda_{4}=\{b_{S}=0,w_{SA}=\frac{\pi}{4},w_{SB}=-\frac{\pi}{4}, b_{P}=-\frac{\pi}{4},w_{PA}=0,w_{PB}=\frac{\pi}{4}\}. \tag{13}\]

Along each path component, referred to as path 1 through 3 in Table 1, the parameters of \(S\) (or \(P\)) are varied and the other held fixed, while remaining in the exactly polarized state. The path is continuous except at a singular point on path 1 where the wave function vanishes at \(\Lambda_{\text{singular}}=\{b_{S}=0,w_{SA}=\frac{\pi}{4},w_{SB}=-\frac{\pi}{4},b_{P}=-\frac{\pi}{4},w_{PA}=\frac{\pi}{4},w_{PB}=0\}\).

### Resolving the special redundancies

In Appendix A.2, we explicitly showed that there can be a large set of redundancies given a polarized state. Hence, for simplicity in the main text, we have used the direct overlap \(S_{q}\) in Eq. (10) as the relevant measure at finite field values. As discussed in the main text, a straightforward way to alleviate the redundancies in the similarity measure \(S_{n}\) in Eq. (19) of the network parameters is to complement it with the direct overlap. By using a combination of both measures, we are able to reduce the amount of computational cost of the direct overlap

\begin{table}
\begin{tabular}{|c||c|c|c|} \hline Path 1 & \(w_{SB}=b_{S}+w_{SA}-\frac{\pi}{2}\) & \(\Lambda_{P}\) fixed & product \(\psi=\psi_{S}\times\psi_{P}\) \\ \(\Lambda_{1}\rightarrow\Lambda_{3}\) & \(w_{SA}:[0,\frac{\pi}{4}),w_{SB}:[\frac{\pi}{4},-\frac{\pi}{4}),b_{S}:[-\frac{ \pi}{4},0)\) & \(w_{PA}=\frac{\pi}{4},w_{PB}=0,b_{P}=-\frac{\pi}{4}\) & \\ \hline \hline \(\cos(b_{X}+w_{XA}+w_{XB})\) & \(\neq 0\) if \(b_{S}+w_{SA}\neq\frac{n}{2}\pi,n\in\mathbb{Z}\to 0\to 1\) & 1 & \(\to 0\to 1\) \\ \hline \(\cos(b_{X}+w_{XA}-w_{XB})\) & & & \(0\) ✓ \\ \hline \(\cos(b_{X}-w_{XA}+w_{XB})\) & \(\cos(2b_{S}-\frac{\pi}{2})\to 0\) & & \(0\) ✓ \\ \hline \(\cos(b_{X}-w_{XA}-w_{XB})\) & & & \(0\) ✓ \\ \hline \hline Path 2 & \(\Lambda_{S}\) fixed & \(w_{PB}=b_{P}-w_{PA}+\frac{\pi}{2}\) & & \\ \(\Lambda_{3}\rightarrow\Lambda_{4}\) & \(w_{SA}=\frac{\pi}{4},w_{SB}=-\frac{\pi}{4},b_{S}=0\) & \(w_{PA}:[\frac{\pi}{4},0],w_{PB}:[0,\frac{\pi}{4}],b_{P}=-\frac{\pi}{4}\) & \\ \hline \(\cos(b_{X}+w_{XA}+w_{XB})\) & 1 & 1 & 1 \\ \hline \(\cos(b_{X}+w_{XA}-w_{XB})\) & 0 & \(\cos(2w_{PA}-\frac{\pi}{2})\to 0\) & 0 ✓ \\ \hline \(\cos(b_{X}-w_{XA}+w_{XB})\) & 0 & & \(0\) ✓ \\ \hline \(\cos(b_{X}-w_{XA}-w_{XB})\) & & & \(0\) ✓ \\ \hline \hline Path 3 & \(w_{SB}=-b_{S}+w_{SA}+\frac{\pi}{2}\) & \(\Lambda_{P}\) fixed & & \\ \(\Lambda_{4}\rightarrow\Lambda_{2}\) & \(w_{SA}=\frac{\pi}{4},w_{SB}:(-\frac{\pi}{4},0],b_{S}:(0,-\frac{\pi}{4})\) & \(w_{PA}=0,w_{PB}=\frac{\pi}{4},b_{P}=-\frac{\pi}{4}\) & \\ \hline \(\cos(b_{X}+w_{XA}+w_{XB})\) & 1 & 1 & 1 \\ \hline \(\cos(b_{X}+w_{XA}-w_{XB})\) & & 0 & 0 ✓ \\ \hline \(\cos(b_{X}-w_{XA}+w_{XB})\) & 0 & & 0 ✓ \\ \hline \(\cos(b_{X}-w_{XA}-w_{XB})\) & & & 0 & ✓ \\ \hline \end{tabular}
\end{table}
Table 1: A path going from \(\Lambda_{1}\) to \(\Lambda_{2}\) is composed of three steps. Path 1 (\(\Lambda_{1}\rightarrow\Lambda_{3}\)) is smooth except at the point \(w_{SA}=\frac{\pi}{4},w_{SB}=-\frac{\pi}{4},b_{S}=0\), where the wavefunction vanishes. This is further denoted by the red arrows first decreasing to 0 before increasing to 1 in the first row. Path 2 and 3 are both smooth. The last column illustrates that the wavefunction \(\psi\) remains in the polarized state along the path.
by a fraction as the similarity is easy to compute. More specifically, we define a mixed measure \(S_{m}\) by replacing a random fraction (given by \(f\)) of the similarity measure pairs \(\left\{l,l^{\prime}\right\}\) by a rescaled overlap measure \(\widetilde{S}_{q}\) such that,

\[S_{m}(l,l^{\prime})=\begin{cases}\widetilde{S}_{q}(l,l^{\prime})&\text{with probability }f,\\ S_{n}(l,l^{\prime})&\text{with probability }1-f.\end{cases} \tag{38}\]

The following rescaling of the overlap measure \(S_{q}\) is necessary as we want to include the two measures on an equal-footing given by,

\[\widetilde{S}_{q} =\frac{S_{q}-n_{q}}{m_{q}-n_{q}}\cdot(m_{n}-n_{n})+n_{n}, \tag{39a}\] \[m_{q} =\max(S_{q}),\quad n_{q}=\min(S_{q}),\] (39b) \[m_{n} =\max(S_{n}),\quad n_{n}=\min(S_{n}). \tag{39c}\]

For example, we see that the minimum of the rescaled overlap is the same as the minimum of the similarity \(\min(\widetilde{S}_{q})=\min(S_{n})\).

In Fig. 9, we demonstrate that by using a mixed measure with a fraction of \(f=0.4\) replacement, our algorithm with DM is able to identify the presence (indicated by the shaded blue region for smaller field values \(h=0.475\) and \(h=0.55\)) and absence (\(h=0.7\)) of superselection sectors across various field values, consistent with the predictions of the algorithm using direct overlap (shown in Fig. 6). We note that in the case with a mixed measure, DM is a natural technique as the algorithm looks for connectivity; whereas kernel PCA would fail to identify such transition (since a fraction of pairs of wave functions are incorrectly considered to be dissimilar by \(S_{n}\), the leading kernel PCA components still show four separated clusters up to the largest magnetic field, \(h=1\)).

## Appendix B Optimization with Variational Monte Carlo

To find the ground state \(\ket{\Psi(\Lambda^{0})}\propto\sum_{\mathbf{\sigma}}\psi(\mathbf{\sigma};\Lambda^{0}) \ket{\mathbf{\sigma}}\), we wish to minimize the energy expectation \(\bra{E}=\bra{\Psi}\hat{H}\ket{\Psi}/\bra{\Psi}\) (omitting the variational parameters \(\Lambda^{0}\) in this section), which is bounded by the ground state energy by the variational principle. An exact computation \(\bra{E}_{\text{exact}}\) is costly as the summation enumerates over exponentially many spin configurations \(\mathbf{\sigma}\) as the system size increases. Here we use variational Monte Carlo (VMC) importance sampling algorithm to estimate such expectation values. The idea is to compute relative probability between different configurations and sample from the true wavefunction probability density \(\ket{\psi(\mathbf{\sigma})}^{2}\), without having to compute \(\ket{\psi(\mathbf{\sigma})}^{2}\) for all \(\mathbf{\sigma}\). To perform this algorithm, we initialize \(M\) random configurations \(\{\mathbf{\sigma}_{i}\}_{i=1}^{M}\) and continue each with random walks based on previous configurations, hence forming \(M\) Markov chains.

In particular, the Metropolis-Rosenbluth algorithm [74] is used to propose the next configuration \(\mathbf{\sigma}^{\prime}_{i}\) that is locally connected to \(c_{i}\) according to function \(g(\mathbf{\sigma}^{\prime}|\mathbf{\sigma})\). For the toric code model, we use two types of proposals: spin flips and vertex flips. Here, we will assume a probability of \(p\) for proposing spin flips and analogously \(1-p\) for vertex flips that are equally likely at all sites:

\[g(\mathbf{\sigma}^{\prime}|\mathbf{\sigma})=\begin{cases}\frac{p}{n_{s}},&\text{for spin flips}\\ \frac{1-p}{n_{v}},&\text{for vertex flips}\end{cases} \tag{40}\]

where \(n_{s}\) and \(n_{v}\) are the number of all possible spin and vertex flips. The acceptance of \(\mathbf{\sigma}^{\prime}\) is determined by a probability,

\[\mathbb{P}_{\text{accept}}(\mathbf{\sigma}\rightarrow\mathbf{\sigma}^{\prime})=\min \left(|\frac{\psi(\mathbf{\sigma}^{\prime})}{\psi(\mathbf{\sigma})}|^{2},\,1\right). \tag{41}\]

Figure 9: DM spectra for different field values \(h=0.475,0.55,0.7\) at \(T=0.3\) using a mixed similarity measure \(S_{m}\) with a fraction \(f=0.4\) in Eq. (38). The blue shaded regions highlight the existence of a range of \(\epsilon\) with spectral gap between the degenerate eigenvalues and the decaying eigenvalues, indicating underlying superselection sectors. As the field value approaches the transition field \(h_{c}\), the range of such region shrinks and disappears at high field \(h=0.7\), indicating the absence of sectors.
The random walks will be repeated long enough so that the final configurations at the tail of the chains \(\Sigma_{\text{MC}}=\{\mathbf{\sigma}_{f}\}_{i=b}^{M}\) approximate samples drawn from the probability distribution \(|\psi(\mathbf{\sigma})|^{2}\). A certain number \(b\) of walkers in each chain are discarded to reduce the biases from initialization of the chains. Then the expectation of an observable \(\hat{O}\) is given by,

\[\left\langle\hat{O}\right\rangle_{\text{MC}} =\frac{\sum_{\mathbf{\sigma}}\psi(\mathbf{\sigma})^{*}\langle\mathbf{\sigma} |\hat{O}|\Psi\rangle}{\sum_{\mathbf{\sigma}}|\psi(\mathbf{\sigma})|^{2}}, \tag{10a}\] \[=\frac{\sum_{\mathbf{\sigma}}|\psi(\mathbf{\sigma})|^{2}\frac{\langle\mathbf{ \sigma}|\hat{O}|\Psi\rangle}{\psi(\mathbf{\sigma})}}{\sum_{\mathbf{\sigma}}|\psi(\mathbf{ \sigma})|^{2}},\] (10b) \[=\frac{1}{M}\sum_{\mathbf{\sigma}\in\Sigma_{\text{MC}}}\frac{\langle \mathbf{\sigma}|\hat{O}|\Psi\rangle}{\psi(\mathbf{\sigma})}. \tag{10c}\]

Defining a local value of the operator \(\hat{O}\) as,

\[O_{\text{loc}}=\frac{\langle\mathbf{\sigma}|\hat{O}|\Psi\rangle}{\psi(\mathbf{\sigma} )}, \tag{11}\]

then the Monte Carlo estimation is the average of the local values in the Markov chain: \(\left\langle\hat{O}\right\rangle_{\text{MC}}=\frac{1}{M}\sum_{\mathbf{\sigma}\in \Sigma_{\text{MC}}}O_{\text{loc}}\).

Next, to minimize \(\langle E\rangle\), we can compute its gradient with respect to the weights \(\Lambda^{0}\) in terms of the local energy \(E_{\text{loc}}\) and wavefunction amplitude derivative \(D_{i}\):

\[\partial_{\Lambda_{i}}\langle E\rangle =\langle E_{\text{loc}}D_{i}\rangle-\langle E_{\text{loc}} \rangle\langle D_{i}\rangle \tag{12a}\] \[E_{\text{loc}} =\frac{\langle\mathbf{\sigma}|\,H\,|\Psi\rangle}{\psi(\mathbf{\sigma})}, \quad D_{i}=\frac{\partial_{\Lambda_{i}}\psi(\mathbf{\sigma})}{\psi(\mathbf{\sigma})} \tag{12b}\]

Finally, we use gradient descent with learning rate \(\lambda\),

\[\Lambda_{i}\rightarrow\Lambda_{i}-\lambda\partial_{\Lambda_{i}}\langle E\rangle, \tag{13}\]

to minimize the energy expectation value. The gradient descent is performed by using an adaptive Adam optimizer [75]. We repeat this training step until empirical convergence.

Note that the RBM ansatz can get stuck in local minima. To find the toric code ground state, we initialize the network parameters close to the analytic solutions in Eq. (10).

### Fidelity

To find the approximate ground states at finite field values \(h\) with step size \(\Delta h\), we initialize the weights to be those from the previous field value \(h-\Delta h\), and then use the current optimized weights as the initialization for the next step \(h+\Delta h\). A good indication of a quantum phase transition is by inspecting the fidelity \(\mathcal{F}(h)\) defined as,

\[\mathcal{F}(h)=|\langle\psi(h)|\psi(h+\Delta h)\rangle|^{2}. \tag{14}\]

The critical field \(h_{c}\) is identified as a dip in the fidelity, indicating an abrupt change in the ground state wavefunction. A field value of \(h_{c}\simeq 0.57\) (at dashed line in Fig. 10) is found for the RBM ansatz. Note that one can get more accurate field value by including loop expectations in the ansatz as done in Ref. [59].

## Appendix C Ensemble generation

Using the algorithm outlined in Sec. 1, we can generate ensembles that deviate from the initial optimized parameters by setting hyper-parameter \(T=0.1,0.3,1\). The other choices of hyper-parameters for the ensembles are number of independent chains \(k=2\), length of each chain \(n=250\), and number of samples kept \(m=n\). The parameter proposal function we use consists of with probability \(p_{m}\) randomly apply minus sign or randomly adding local noise at a single spin site \(\jmath\). More precisely,

\[f(\Lambda,\xi) =\begin{cases}f_{-,\jmath},&\text{with probability}:p_{m},\\ f_{\text{local},\jmath},&\text{with probability}:1-p_{m},\end{cases} \tag{15a}\] \[f_{-,\jmath} =\begin{cases}-(\Lambda)_{i},&i\in\jmath\\ (\Lambda)_{i},&i\not\in\jmath\end{cases}\] (15b) \[f_{\text{local},\jmath} =\begin{cases}\text{uniform}(0,\xi)+(\Lambda)_{\text{i}},&\text{ i}\in\jmath\\ (\Lambda)_{i},&i\not\in\jmath\end{cases} \tag{15c}\]

In the exact toric code state, \(f_{-,\jmath}\) corresponds to act \(\sigma_{x}\) operator at site \(\jmath\) to create a pair of m-particles. In the trivial phase, depending on the parametrization of the state, \(f_{-,\jmath}\) could correspond to a single spin flip at site \(\jmath\). The hyperparameters are chosen to be \(p_{m}=0.3\) and \(\xi=0.2\). In Fig. 11, we visualize the ensembles by computing their loop expectations \(\langle\overline{W}_{j}\rangle\) at different field values.

Figure 10: Fidelity \(\mathcal{F}\) as a function of field \(h\). The red dashed line is drawn to guide the eye, where the dip in fidelity indicates the critical field value \(h_{c}\simeq 0.57\).