Investigating Topological Order using Recurrent Neural Networks
Recurrent neural networks (RNNs), originally developed for natural language
processing, hold great promise for accurately describing strongly correlated
quantum many-body systems. Here, we employ 2D RNNs to investigate two
prototypical quantum many-body Hamiltonians exhibiting topological order.
Specifically, we demonstrate that RNN wave functions can effectively capture
the topological order of the toric code and a Bose-Hubbard spin liquid on the
kagome lattice by estimating their topological entanglement entropies. We also
find that RNNs favor coherent superpositions of minimally-entangled states over
minimally-entangled states themselves. Overall, our findings demonstrate that
RNN wave functions constitute a powerful tool to study phases of matter beyond
Landau's symmetry-breaking paradigm.
Mohamed Hibat-Allah, Roger G. Melko, Juan Carrasquilla
# Investigating Topological Order using Recurrent Neural Networks

###### Abstract

Recurrent neural networks (RNNs), originally developed for natural language processing, hold great promise for accurately describing strongly correlated quantum many-body systems. Here, we employ 2D RNNs to investigate two prototypical quantum many-body Hamiltonians exhibiting topological order. Specifically, we demonstrate that RNN wave functions can effectively capture the topological order of the toric code and a Bose-Hubbard spin liquid on the kagome lattice by estimating their topological entanglement entropies. We also find that RNNs favor coherent superpositions of minimally-entangled states over minimally-entangled states themselves. Overall, our findings demonstrate that RNN wave functions constitute a powerful tool to study phases of matter beyond Landau's symmetry-breaking paradigm.

## I Introduction

Landau symmetry breaking theory provides a fundamental description of a wide range phases of matter and their phase transitions through the use of local order parameters [1]. Despite the fact that a great deal of our theoretical and experimental investigations of interacting quantum many-body systems have been developed with the aim of studying local order parameters, it is well-known that the most intriguing strongly correlated phases of matter may not be easily characterized through these observables. Instead, several states of matter seen in modern theoretical and experimental studies are characterized using non-local order parameters that rely on the phases' topological properties [2; 3; 4]. Topological order, in particular, refers to a type of order characterized by the emergence of quasi-particle anyonic excitations, topological invariants, and long-range entanglement, which typically do not appear in traditional forms of order. As a result of these properties, topologically ordered phases have been suggested as an important building block for the development of a protected qubit resistant to perturbations and errors [5; 6; 7]. Interestingly, such qubits have been devised recently at the experimental level [8].

While most manifestations of topological order are dynamical in nature--e.g. anyon statistics, ground state degeneracy, and edge excitations [9]-topological order can also be characterized directly in terms of the ground state wave function and its entanglement. In particular, a probe for topological order is the topological entanglement entropy (TEE) [9; 10], which offers a characterization of the global entanglement pattern of topological ground states not present in conventionally ordered systems. Notably, the TEE is readily accessible for large classes of topological orders [9; 11], in numerical simulations based on quantum Monte Carlo (QMC) [12; 13; 14] and density matrix renormalization group (DMRG) [15; 16], as well as in experimental realizations of topological order based on gate-based quantum computers [4]. Machine learning (ML) techniques offer an alternative approach study quantum many-body systems and have proved useful for a wide array of tasks including the classification of phases of matter [17; 18; 19; 20], quantum state tomography [21; 22], finding ground states of quantum systems [23; 24; 25; 26; 27; 28; 29; 30; 31; 32], studying open quantum systems [33; 34], and simulating quantum circuits [35; 36; 37], among many others [38; 39; 40; 41]. In particular, neural networks representations of quantum many-body states have been shown to be able of expressing topological order using, e.g., restricted Boltzmann machines [42; 43; 44; 45], convolutional neural networks [17] and autoregressive neural networks [46]. Here we use recurrent neural networks (RNN) [47; 48; 49] as an ansatz wave function [31; 32] to investigate topological order in 2D through the estimation of the TEE. We focus on two model Hamiltonians exhibiting topological order, namely Kitaev's toric code [6; 7] and a Bose-Hubbard model on the kagome lattice previously shown to host a gapped quantum spin liquid with non-trivial emergent \(\mathbb{Z}_{2}\) gauge symmetry [12; 14; 50]. In our study, we use Kitaev-Preskill constructions [10], and finite size-scaling analysis of the entanglement entropy to extract the TEE. We find convincing evidence that RNNs are capable of expressing ground states of Hamiltonians displaying topological order. We also find evidence that the RNN wave function is naturally biased toward finding superpositions of minimally entangled states, as reflected in the calculations of entanglement entropy and Wilson loop operators for the toric code. Overall, our results indicate that RNNs can represent phases of matter beyond the conventional Landau symmetry-breaking paradigm.
## II 2D RNNs

Our main aim is to study topological properties of Hamiltonians using an RNN wave function ansatz. Since the quantum systems we study are stoquastic [51], we consider an ansatz with positive amplitudes to model the ground state wave function [31]. Complex extensions of RNN wave functions for non-stoquastic Hamiltonians have been explored in Refs. [31; 32]. To model a positive RNN wave function, we write our ansatz in the computational basis as:

\[\Psi_{\mathbf{\theta}}(\mathbf{\sigma})=\sqrt{p_{\mathbf{\theta}}(\mathbf{\sigma})},\]

where \(\mathbf{\theta}\) denotes the variational parameters of the ansatz \(|\Psi_{\mathbf{\theta}}\rangle\), and \(\mathbf{\sigma}=(\sigma_{1},\sigma_{2},\ldots,\sigma_{N})\) is a basis state configuration. A key characteristic of the RNN wave function is its ability to estimate observables with uncorrelated samples through autoregressive sampling [31; 52]. This is achieved by parameterizing the joint probability \(p_{\mathbf{\theta}}(\mathbf{\sigma})\) with its conditionals \(p_{\mathbf{\theta}}(\sigma_{i}|\sigma_{<i})\), through the probability chain rule

\[p_{\mathbf{\theta}}(\mathbf{\sigma})=p_{\mathbf{\theta}}(\sigma_{1})p_{\mathbf{\theta}}( \sigma_{2}|\sigma_{1})\cdots p_{\mathbf{\theta}}(\sigma_{N}|\sigma_{N-1},\ldots, \sigma_{2},\sigma_{1}).\]

The conditionals are given by

\[p_{\mathbf{\theta}}(\sigma_{i}|\sigma_{<i})=\mathbf{y}_{i}\cdot\mathbf{\sigma}_{i},\]

where \(\mathbf{y}_{i}=\mathrm{Softmax}(U\mathbf{h}_{i}+\mathbf{c})\) and '\(\cdot\)' is the dot product operation. The hidden state \(\mathbf{h}_{i}\) is calculated recursively as [49]

\[\mathbf{h}_{i}=f\!(W[\mathbf{\sigma}_{i-1};\mathbf{h}_{i-1}]+\mathbf{b}), \tag{1}\]

where the input \(\mathbf{\sigma}_{i-1}\) is a one-hot encoding of \(\sigma_{i-1}\) and the symbol \([\cdot;\cdot]\) corresponds to the concatenation of two vectors. Furthermore, \(U,W,\mathbf{b}\), and \(\mathbf{c}\) are learnable weights and biases and \(f\) is an activation function. The sequential operation of the RNN is shown in Fig. 1(a), where the RNN cell, i.e., the recurrent relation in Eq. (1), is depicted as a blue square. As each of the conditionals \(p_{\mathbf{\theta}}(\sigma_{i}|\sigma_{<i})\)[31] is normalized, the distribution \(p_{\mathbf{\theta}}\), and thus the quantum state \(|\Psi_{\mathbf{\theta}}\rangle\), are normalized. Notably, by virtue of the sequential structure built into the RNN ansatz, it is possible to obtain exact samples from \(p_{\mathbf{\theta}}\) by sampling the conditionals \(p_{\mathbf{\theta}}(\sigma_{i}|\sigma_{<i})\) sequentially as illustrated in Fig. 1(b). The sampling scheme is parallelizable and can produce fully uncorrelated samples distributed according to \(p_{\mathbf{\theta}}\) without the use of potentially slow Markov chains [31; 32].

As our aim is to study 2D quantum systems with periodic boundary conditions, we use 2D RNNs [31; 53], through the modification of the 1D relation in Eq. (1) to a recursion that encodes the 2D geometry of the lattice, i.e.,

\[\mathbf{h}_{i,j}=f\!\Big{(}W[\mathbf{\sigma}_{i-1,j};\mathbf{\sigma}_{i,j-1}; \mathbf{\sigma}_{\mathrm{mod}(i+1,L_{x}),j};\mathbf{\sigma}_{i,\mathrm{mod}(j+1,L_{y} )};\] \[\mathbf{h}_{i-1,j};\mathbf{h}_{i,j-1};\mathbf{h}_{\mathrm{mod}(i+1,L_{x}),j}; \mathbf{h}_{i,\mathrm{mod}(j+1,L_{y})}]+\mathbf{b}\Big{)}. \tag{2}\]

\(L_{x},L_{y}\) are respectively the width and the length of the 2D lattice. In our study, we choose \(L_{x}=L_{y}=L\). Additionally, \(\mathbf{h}_{i,j}\) is a hidden state with two indices for each site in the 2D lattice, which is computed based on the inputs and the hidden states of the nearest neighboring sites. Since \(\mathbf{h}_{i,j}\) contains information about of the history of generated variables \(\sigma_{i,j}\), it can be used to compute the conditionals

\[p_{\mathbf{\theta}}(\sigma_{i,j}|\sigma_{<i,j})=\mathrm{Softmax}(U\mathbf{h}_{i,j}+ \mathbf{c})\cdot\mathbf{\sigma}_{i,j}. \tag{3}\]

The additional variables \(\mathbf{\sigma}_{\mathrm{mod}(i+1,L_{x}),j}\), \(\mathbf{\sigma}_{i,\mathrm{mod}(j+1,L_{y})}\) and hidden states \(\mathbf{h}_{\mathrm{mod}(i+1,L_{x}),j}\), \(\mathbf{h}_{i,\mathrm{mod}(j+1,L_{y})}\) allows to model systems with periodic boundary conditions such that the ansatz accounts for the correlations between physical degrees of freedom across the boundaries. This approach has been also suggested and implemented in Ref. [46]. We note that during the process of autoregressive sampling, if either of the input vectors, in Eq. (2), have not been encountered yet, we initialize them to a null vector so that we preserve the autoregressive nature of the RNN wave function, as illustrated in Fig. 1(b). Furthermore, Fig. 1(c) illustrates the autoregressive sampling path in 2D as well as how information is being transferred among RNN cells. Importantly, we use an advanced version of 2D RNNs which incorporates a gating mechanism as previously done in Refs. [31; 46; 54; 55]. Additional details can be found in App. A. We also note that implementing lattice symmetries in our RNN ansatz can be done to improve the variational accuracy as shown in Refs. [31; 55], however, we do not pursue this direction in our study.

### Supplementing RNNs optimization with annealing

To train the parameters of the RNN, we minimize the energy expectation value \(E_{\mathbf{\theta}}=\bra{\Psi_{\mathbf{\theta}}}\hat{H}\ket{\Psi_{\mathbf{\theta}}}\) using Variational Monte Carlo (VMC) [56], where \(\hat{H}\) is a Hamiltonian of interest. In the presence of local minima in the optimization landscape of \(E_{\mathbf{\theta}}\), the VMC optimization may get stuck in a poor local optimum [57; 58]. To ameliorate this limitation, we supplement the VMC scheme with a pseudo-entropy whose objective is to help the optimization escape local minima [57; 59; 55; 57; 59]. The new objective function is defined as

\[F_{\mathbf{\theta}}(n)=E_{\mathbf{\theta}}-T(n)S_{\mathrm{classical}}(p_{\mathbf{\theta}}), \tag{4}\]

where \(F_{\mathbf{\theta}}\) is a variational pseudo-free energy. The Shannon entropy \(S_{\mathrm{classical}}\) of \(p_{\mathbf{\theta}}\left(\mathbf{\sigma}\right)\) is given by

\[S_{\mathrm{classical}}(p_{\mathbf{\theta}})=-\sum_{\mathbf{\sigma}}p_{\mathbf{\theta}}( \mathbf{\sigma})\ln\left(p_{\mathbf{\theta}}(\mathbf{\sigma})\right), \tag{5}\]

where the sum goes over all possible configurations \(\{\mathbf{\sigma}\}\) in the computational basis. The pseudo-entropy \(S_{\mathrm{classical}}\) and its gradients are evaluated by sampling the
RNN wave function. Furthermore, \(T(n)\) is a pseudo-temperature that is annealed from some initial value \(T_{0}\) to zero as follows: \(T(n)=T_{0}(1-n/N_{\text{annealing}})\) where \(n\in[0,N_{\text{annealing}}]\) and \(N_{\text{annealing}}\) is the total number of annealing steps. This scheme is inspired by the regularized variational quantum annealing scheme in Refs. [32; 55; 57]. More details about our training scheme are given in App. B. We also provide the hyperparameters in App. C.

### Topological entanglement entropy

A powerful tool to probe topologically ordered states of matter is through the so-called topological entanglement entropy (TEE) [61; 62; 63; 10; 11; 12; 13; 64]. The TEE can be extracted by computing the entanglement entropy of a spatial bipartition of the system into \(A\) and \(B\), which together comprise the full system. For many phases of 2D matter, the Renyi-\(n\) entropy \(S_{n}(A)\equiv\frac{1}{1-n}\ln(\text{Tr}(\rho_{A}^{n}))\) satisfies the area law \(S_{n}(A)=aL-\gamma+\mathcal{O}(L^{-1})\). Here \(L\) is the size of the boundary between \(A\) and \(B\), \(\rho_{A}=\text{Tr}_{B}|\Psi\rangle\langle\Psi|\) is the reduced density matrix of subsystem \(A\), \(|\Psi\rangle\) is the state of the system, and \(\gamma\) is the TEE. The latter detects non-local correlations in the ground state wave function and plays the role of an order parameter for topological phases similar to the notion of a local order parameter in phases displaying long-range order. Interestingly, a measure of specific non-zero values of \(\gamma\) can be a clear sign of the existence of a topological order in a system of interest. Additionally, since the TEE is shown to be independent of the choice of Renyi index \(n\) for a contractible region \(A\)[62], we can use the swap trick [65] with our RNN wave function ansatz [66; 31] to calculate the second Renyi entropy \(S_{2}\) and extract the TEE \(\gamma\).

To access the TEE \(\gamma\), we can approximate the ground state of the system using an RNN wave function ansatz, i.e. \(|\Psi_{\mathbf{\theta}}\rangle\approx|\Psi\rangle\) for different system sizes followed by a finite-size scaling analysis of the second Renyi entropy. We can also make use of a TEE construction, e.g., the Kitaev-Preskill construction [10].

The Kitaev-Preskill construction prescribes dividing the system into four subregions \(A\), \(B\), \(C\), and \(D\) as illustrated in Fig. 2. The TEE

Figure 1: (a) An illustration of the positive RNN wave function. Each RNN cell (blue squares) receives an input \(\mathbf{\sigma}_{n-1}\) and a hidden state \(\mathbf{h}_{n-1}\) and outputs a new hidden state \(\mathbf{h}_{n}\). The latter is fed to a Softmax layer (denoted by S, red circle) that outputs a conditional probability \(P_{i}\) (orange circle). The conditional probabilities are multiplied after taking a square root to obtain the wave function \(\Psi_{\mathbf{\theta}}(\mathbf{\sigma})\) (pink circle). (b) Illustration of the sampling scheme for an RNN wave function. After obtaining the probability vector \(\mathbf{y}_{i}\) from the Softmax layer (S) at step \(i\), we sample it to produce \(\mathbf{\sigma}_{i}\) which is fed again to the RNN with the hidden state \(\mathbf{h}_{i}\) to produce the next configuration \(\mathbf{\sigma}_{i+1}\).(c) A 2D RNN with periodic boundary conditions. A bulk RNN cell receives two hidden states \(\mathbf{h}_{i,j-1}\) and \(\mathbf{h}_{i-1,j}\), as well as two input vectors \(\mathbf{\sigma}_{i,j-1}\) and \(\mathbf{\sigma}_{i-1,j}\) (not shown) illustrated by the black solid arrows. To handle periodic boundary conditions, RNN cells at the boundary receive an additional \(\mathbf{h}_{i,\text{mod}(j+1,L_{y})}\) and \(\mathbf{h}_{\text{mod}(i+1,L_{x}),j}\), as well as two input vectors \(\mathbf{\sigma}_{i,\text{mod}(j+1,L_{y})}\) and \(\mathbf{\sigma}_{\text{mod}(i+1,L_{x}),j}\) (not shown) illustrated by green solid arrows. The sampling path is illustrated with red dashed arrows. The initial memory state \(\mathbf{h}_{0}\) of the 2D RNN and the initial inputs \(\mathbf{\sigma}_{0}\) (not shown) are taken as null vectors.

Figure 2: A sketch of the parts \(A\), \(B\) and \(C\) that we use for Kitaev-Preskill construction to compute the TEE in a system of interest.
computing

\[\gamma =-S_{2}(A)-S_{2}(B)-S_{2}(C)+S_{2}(AB)\] \[\quad+S_{2}(AC)+S_{2}(BC)-S_{2}(ABC),\]

where \(S_{2}(A)\) is the second Renyi entropy of the subsystem \(A\), and \(AB\) is the union of \(A\) and \(B\) and similarly for the other terms. Finite-size effects on \(\gamma\) can be alleviated by increasing the size of the subregions \(A,B\) and \(C\)[10, 67]. Finally, we highlight the ability of the RNN wave function to study systems with fully periodic boundary conditions as a strategy to mitigate boundary effects, as opposed to cylinders used in DMRG [68, 69], which may potentially introduce edge effects that can affect the values of the TEE [70].

## III Results

### The toric code

We now focus our attention on the toric code Hamiltonian which is the simplest model that hosts a \(Z_{2}\) topological order [6, 61] and has a non-zero TEE equal to \(\gamma=\ln(2)\). The Hamiltonian is defined in terms of spin-1/2 degrees of freedom located on the edges of a square lattice (see Fig. 3(a)) and is given by

\[\hat{H}=-\sum_{p}\prod_{i\in p}\tilde{\sigma}_{i}^{z}-\sum_{v}\prod_{i\in v} \hat{\sigma}_{i}^{x},\]

where \(\hat{\sigma}_{i}^{x,z}\) are Pauli matrices. Additionally, the first summation is on the plaquettes and the second summation is on the vertices of the lattice [61]. Note that the lattice in Fig. 3(a) can be seen as a square lattice with a unit cell containing two spins. In our simulations, we use an \(L\times L\times 2\) array of spins where \(L\) is the number of plaquettes on each side of the underlying square lattice. It is possible to study the toric code with a 2D RNN defined on a primitive square lattice by merging the two spin degrees of freedom of the unit cell of the toric code into a single "patch" followed by an enlargement of the local Hilbert space dimension in the RNN from 2 to 4. This idea is illustrated in Fig. 3(a) and is similar in spirit to how the local Hilbert space is enlarged in DMRG to study quasi-1D systems [71]. We provide additional details about the mapping in App. A.

To extract the TEE from our ansatz, we variationally optimize the 2D RNN wave function targetting the ground state of this model for multiple system sizes on a square lattice with periodic boundary conditions. After the optimization, we compute the TEE using system size extrapolation and using the Kitaev-Preskill scheme provided in Sec. II.2. More details about the regions chosen for this construction are provided in App. D. To avoid local minima during the variational optimization, we perform an initial annealing phase as described in Sec. II.1 (see additional details in App. C).

The results shown in Fig. 4(a) suggest that our 2D RNN wave function can describe states with an area law scaling in 2D. Linearized versions of the RNN wave function have been recently shown to display an entanglement area law [72]. For \(L=10\) (not included in the extrapolations in Fig. 4(a)), it is challenging to evaluate \(S_{2}\) accurately as the expectation value of the swap operator is proportional to \(\exp\left(-S_{2}\right)\), which becomes very small and is hard resolve accurately via sampling the RNN wave function. The improved ratio trick is an interesting alternative for enhancing the accuracy of our estimates [73, 65]. The use of conditional sampling is also another possibility for enhancing the accuracy of our measurements [66].

The extrapolation confirms the existence of a non-zero TEE whose value is close to \(\gamma^{\prime}=\ln(2)\) within error bars. Note that the sub-region we use to compute the TEE is half of the torus, namely a cylinder with two disconnected boundaries [74]. As shown in Ref. [63], the use of this geometry means that the expected TEE becomes state-dependent and given by

\[\gamma^{\prime}=2\gamma+\ln\left(\sum_{i}\frac{p_{i}^{2}}{d_{i}^{2}}\right) \tag{6}\]

for the second Renyi entropy. Here \(d_{i}\geq 1\) is the quantum dimension of an \(i\)-th quasi-particle. For the toric code, we have abelian anyons with \(d_{i}=1\) for \(i=1,2,3,4\). Additionally \(p_{i}=|\alpha_{i}|^{2}\) is the overlap of the computed ground state \(|\Psi\rangle\) with the \(i\)-th minimally entangled state (MES) \(|\Xi_{i}\rangle\) where

\[|\Psi\rangle=\sum_{i}\alpha_{i}\left|\Xi_{i}\right\rangle.\]

The observations above and the numerical result \(\gamma_{\text{RNN}}\approx\ln(2)\) suggest that the RNN wave functions optimized via gradient descent and annealing find a superposition of MES, as opposed to DMRG which preferentially collapses to a single MES for relatively low bond dimensions. For relatively large bond dimensions a superposition of

Figure 3: Mapping of 2D toric code lattice and kagome lattice to a square lattice that can be handled by a 2D RNN wave function. For the toric code lattice in panel (a), every two sites inside the dashed green ellipses are merged. For the kagome lattice in panel (c), every three sites in a unit cell enclosed by the dashed green circles are combined.
MES can be recovered in a DMRG simulation [15; 16]. The analysis provided in App. E demonstrates that our optimized RNN ansatz finds a uniform superposition of two MES which increases the entanglement in the state with respect to a single MES. Thus using Eq. (6), we expect \(\gamma^{\prime}=2\ln(2)+\ln\left(\frac{1}{4}+\frac{1}{4}\right)=\ln(2)\), which is consistent with our numerical observations.

We note that the exact autoregressive sampling procedure plays a key role in the ability of our RNN ansatz to sample a superposition of different topological sectors when this superposition is encoded in our ansatz. For wave functions representing the ground state of the toric code used in combination with Markov-chain Monte Carlo methods, the probability of sampling different topological sectors of the state is exponentially suppressed even if the exact wave function ansatz encodes different topological sectors. This observation can be illustrated using an exact convolutional neural network construction of the toric code ground state which contains an equal superposition of different topological sectors [17]. Although in principle such representation contains all topological sectors, its form is not amenable to exact sampling and uses Markov chains so that upon sampling with local moves the system chooses a fixed topological sector.

To further verify that our 2D RNN wave function can extract the correct TEE of the 2D toric code, we compute the TEE using the Preskill-Kitaev construction, that has contractible surfaces, and for which the TEE does not depend on the topological sector superposition [63; 13] (see App. D for details about the construction). The results reported in Fig. 4(b) demonstrate an excellent agreement between the TEE extracted by our RNN and the expected theoretical value for the toric code.

### Bose-Hubbard model on kagome lattice

We now turn our attention to a hard-core Bose-Hubbard model on the Kagome lattice, which has been shown to host topological order [50; 12; 14]. The Hamiltonian of this model is given by

\[\hat{H}=-t\sum_{\langle i,j\rangle}\left(b_{i}^{\dagger}b_{j}+b_{i}b_{j}^{ \dagger}\right)+V\sum_{\hexagon{O}}n_{\hexagon{O}}^{2}, \tag{7}\]

where \(b_{i}\) (\(b_{i}^{\dagger}\)) is the annihilation (creation) operator. Furthermore, \(t\) is the kinetic strength, \(V\) is a tunable interaction strength and \(n_{\hexagon{O}}=\sum_{i\in\hexagon{O}}(n_{i}-1/2)\). The first term corresponds to a kinetic term that favors hopping between nearest neighbors, whereas the second term promotes an occupation of three hard-core bosons in each hexagon of the kagome lattice. In our setup, we choose \(V\) in units of the kinetic term strength \(t\).

The atom configurations of this model correspond to an \(L\times L\times 3\) array of binary degrees of freedom where \(L\) is the size of each side of the kagome lattice. Following an analogous approach to the toric code, we combine three sites of the unit cell of the kagome lattice as input to the 2D RNN cell, as illustrated in Fig. 3(b). This allows us to map our kagome lattice with a local Hilbert space of 2 to a square lattice with an enlarged Hilbert space of size \(2^{3}=8\).

The model is known to host a \(Z_{2}\) spin-liquid phase for \(V\gtrsim 7\)[12; 14; 75]. To confirm this finding, we estimate \(\gamma\) for the system sizes \(6\times 6\times 3\) and \(8\times 8\times 3\). We use the Kitaev-Preskill construction [10]. The details of the construction of the regions \(A,B\) and \(C\) are provided in App. D. As the Hamiltonian in Eq. 7 has a \(U(1)\) symmetry associated with the conservation of bosons in the system, we impose this symmetry on our RNN wave function [31]. We also supplement the VMC optimization with annealing to overcome local minima as previously done for the 2D toric code (see App. B). For the system size \(8\times 8\times 3\), the RNN ansatz parameters were initialized using the optimized parameters of the system

Figure 4: Entanglement properties of the 2D toric code. (a) Second Renyi entropy scaling is computed using our RNN wave function on the 2D toric code for different lengths \(L\) where the total system size is given as \(L\times L\times 2\). (b) TEE computed with the Kitaev-Preskill construction (see App. D). The values found by the RNN are very close to \(\ln(2)\). Error bars correspond to one standard deviation and are smaller than the symbol size.
size \(6\times 6\times 3\) (see details about the hyperparameters in App. C). This pre-training technique has been motivated in Refs. [32, 46, 55].

The results are provided in Fig. 5. The computed TEEs for \(L=6,8\) show a saturation of \(\gamma_{\text{RNN}}\) for large values of the interaction strength \(V\). We observe that the saturation values of \(\gamma_{\text{RNN}}\) are in good agreement with the expected TEE \(\gamma=\ln(2)\) of a \(Z_{2}\) spin-liquid [12]. Additionally, the negative values of \(\gamma_{\text{RNN}}\) observed for \(V\leq 6\) in the superfluid phase [12] may be related to the presence of Goldstone modes that manifest themselves as corrections to the area law in the entanglement entropy and can be seen as a negative contribution to the TEE [76]. We note that the QMC methods are capable of obtaining a consistent value with the exact TEE for this model at \(V=8\) for very large system sizes [14] using finite-size extrapolation. This observation suggests that our RNN ansatz is still limited by finite-size effects at \(V=8\) (see Fig. 5) for which the TEE is not yet saturated to \(\ln 2\). Other sources of error in our calculation may be due to inaccuracies in the variational calculations and statistical errors due to the sampling. However, we note that our variational calculation is performed at zero temperature, which makes our calculations insensitive to temperature effects as opposed to QMC [12].

## IV Conclusions and Outlooks

We have demonstrated a successful application of neural network wave functions to the task of detecting topological order in quantum systems. In particular we use RNNs borrowed from natural language processing as ansatz wave functions. RNNs enjoy the autoregressive property which allows to sample uncorrelated configurations. They are also capable of estimating second Renyi entropies using the swap trick [31] with which we computed TEEs using finite-size scaling and Kitaev-Preskill constructions [10]. Furthermore, the structural flexibility of the RNN offers the possibility to handle a wide variety of geometries including periodic boundary conditions in any spatial dimension which alleviate boundary effects on the TEE.

We have empirically demonstrated that 2D RNN wave functions support the 2D area law and can find a non-zero TEE for the toric code and for the hard-core Bose-Hubbard model on the Kagome lattice. We also find that RNNs favor coherent superpositions of minimally-entangled states over minimally-entangled states themselves. The success of our numerical experiments hinges on the combination of the exact sampling strategy used to compute observables, the structural properties of the RNN wave function, and the use of annealing as a strategy to overcome local minima during the optimization procedure.

The accuracy improvement of our findings can be achieved through the use of more advanced versions of RNNs and autoregressive models in general [72, 55], or even a hybrid approach that combines QMC and RNNs [77, 78]. Similarly, the incorporation of lattice symmetries provides a strategy to enhance the accuracy of our calculations [31, 55, 79]. Although our results match the anticipated behaviour of the toric code and Bose-Hubbard spin liquid models, we highlight that the RNN wave function may be susceptible to spurious contributions to the TEE [64] and we have not addressed this issue in our work.

Finally, our methods can be applied to study other systems displaying topological order, such as the Rydberg atoms array [3, 80, 81], either through variational methods or in combination with experimental data. To experimentally study topological order, it is possible to use quantum state tomography with RNNs [22]. This involves using experimental data to reconstruct the state seen in the experiment followed by an estimation of the TEE using the methods outlined in our work. Overall, our findings suggest that RNN wave function ansatzes have promising potential for discovering new phases of matter with topological order.

###### Acknowledgements.

 We would like to thank Giacomo Torlai, Jeremy Cote, Schuyler Moss, Roeland Wiersema, Ejaaz Merali, Isaac De Vlugt and Arun Paramekanti for helpful discussions. Our RNN implementation is based on Tensorflow [82] and NumPy [83]. Computer simulations were made possible thanks to the Vector Institute computing cluster and Compute Canada. M.H acknowledges support from Mitacs Atecelerate. J.C. acknowledges support from Natural Sciences and Engineering Research Council of Canada (NSERC), the Shared Hierarchical

Figure 5: A plot of the topological entanglement entropy against the interaction strength \(V\) (in units of \(t\)) of Hard-core Bose-Hubbard model on kagome lattice for system sizes \(N=L\times L\times 3\) where \(L=6,8\). The calculations were performed using the Kitaev-Preskill construction (see App. D). The continuous black horizontal line corresponds to a zero TEE, and the dashed blue horizontal line for a \(\ln(2)\) TEE.
Academic Research Computing Network (SHARCNET), Compute Canada, and the Canadian Institute for Advanced Research (CIFAR) AI chair program. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade.

## Appendix A 2D periodic gated RNNs

In this appendix, we describe our implementation of a 2D gated RNN wave function for periodic systems, that can be used to approximate the ground states of the Hamiltonians considered in this paper. If we define

\[\mathbf{h}^{\prime}_{i,j} =[\mathbf{h}_{i-1,j};\mathbf{h}_{i,j-1};\mathbf{h}_{\mathrm{mod}(i+1,L_{x}), j};\mathbf{h}_{i,\mathrm{mod}(j+1,L_{y})}],\] \[\mathbf{\sigma}^{\prime}_{i,j} =[\mathbf{\sigma}_{i-1,j};\mathbf{\sigma}_{i,j-1};\mathbf{\sigma}_{\mathrm{ mod}(i+1,L_{x}),j};\mathbf{\sigma}_{i,\mathrm{mod}(j+1,L_{y})}],\]

then our gated 2D RNN wave function ansatz is based on the following recursion relations:

\[\tilde{\mathbf{h}}_{i,j} =\tanh\Bigl{(}W[\mathbf{\sigma}^{\prime}_{i,j};\mathbf{h}^{\prime}_{i,j}] +\mathbf{b}\Bigr{)},\] \[\mathbf{u}_{i,j} =\mathrm{sigmoid}\Bigl{(}W_{g}[\mathbf{\sigma}^{\prime}_{i,j};\mathbf{h}^ {\prime}_{i,j}]+\mathbf{b}_{g}\Bigr{)},\] \[\mathbf{h}_{i,j} =\mathbf{u}_{i,j}\odot\tilde{\mathbf{h}}_{i,j}+(1-\mathbf{u}_{i,j})\odot(U \mathbf{h}^{\prime}_{i,j}).\]

Here '\(\odot\)' is the Hadamard (element-wise) product. A hidden state \(\mathbf{h}_{i,j}\) can be obtained by combining a candidate state \(\tilde{\mathbf{h}}_{i,j}\) and the neighbouring hidden states \(\mathbf{h}_{i-1,j},\mathbf{h}_{i,j-1},\mathbf{h}_{\mathrm{mod}(i+1,L_{x}),j},\mathbf{h}_{i, \mathrm{mod}(j+1,L_{y})}\). The update \(\mathbf{u}_{i,j}\) determines how much of the candidate hidden state \(\tilde{\mathbf{h}}_{i,j}\) will be taken into account and how much of the neighboring states will be considered. With this combination, it is possible to circumvent some limitations of the vanishing gradient problems [84; 85]. The weight matrices \(W,W_{g},U\) and the biases \(b,b_{g}\) are variational parameters of our RNN ansatz. The size of the hidden state \(\mathbf{h}_{i,j}\) is a hyperparameter called the number of memory units or the hidden dimension and it is denoted as \(d_{h}\). Finally, to motivate the use of the gating mechanism, we highlight that a gated 2D RNN was found to be better than a non-gated 2D RNN for the task of finding the ground state of a 2D Heisenberg model [55].

Since we use an enlarged Hilbert space in our 2D RNN, we take \(\mathbf{\sigma}_{ij}\) to be the concatenation of the one-hot encodings of the binary physical variables in each unit cell. In this case, if \(m\) is the number of physical variables per unit cell, then \(\mathbf{\sigma}_{ij}\) has size \(2m\). We also note that the size of the Softmax layer is taken as \(2^{m}\) so that we can autoregressively sample each unit cell variables at once.

## Appendix B Variational monte carlo and variance reduction

To optimize the energy expectation value of our RNN wave function \(\ket{\Psi_{\mathbf{\theta}}}\), we use the Variational Monte Carlo (VMC) scheme, which consists of using importance sampling to estimate the expectation value of the energy \(E_{\mathbf{\theta}}=\bra{\Psi_{\mathbf{\theta}}}\hat{H}\ket{\Psi_{\mathbf{\theta}}}\) as follows [56; 31]:

\[E_{\mathbf{\theta}}=\frac{1}{M}\sum_{i=1}^{M}E_{\mathrm{loc}}(\mathbf{\sigma^{(i)}}),\]

where the local energies \(E_{\mathrm{loc}}\) are defined as

\[E_{\mathrm{loc}}(\mathbf{\sigma})=\sum_{\mathbf{\sigma}^{\prime}}H_{\mathbf{\sigma}\mathbf{ \sigma}^{\prime}}\frac{\Psi_{\mathbf{\theta}}(\mathbf{\sigma}^{\prime})}{\Psi_{\mathbf{ \theta}}(\mathbf{\sigma})}.\]

Here the configurations \(\{\mathbf{\sigma^{(i)}}\}_{i=1}^{M}\) are sampled from our ansatz using autoregressive sampling. The choice of \(M\) is a hyperparameter that can be tuned. Furthermore, \(E_{\mathrm{loc}}(\mathbf{\sigma})\) can be efficiently computed for local Hamiltonians. Furthermore, the gradients can be estimated as [31]

\[\partial_{\mathbf{\theta}}E_{\mathbf{\theta}}=\frac{2}{M}\mathfrak{Re}\left(\sum_{i=1} ^{M}\partial_{\mathbf{\theta}}\log\left(\Psi_{\mathbf{\theta}}^{*}(\mathbf{\sigma^{(i)}}) \right)\left(E_{\mathrm{loc}}(\mathbf{\sigma^{(i)}})-E_{\mathbf{\theta}}\right)\right).\]

For a stoquastic Hamiltonian \(\hat{H}\)[51], we can use a positive RNN wave function where the use of the real part \(\mathfrak{Re}\) is not necessary. Importantly, subtracting the variational energy \(E_{\mathbf{\theta}}\) is helpful to achieve convergence as it reduces the variance of the gradients near convergence without biasing its expectation value as shown in Ref. [31]. The subtracted term is referred to as a baseline, which is typically used for the same purpose in the context of Reinforcement learning [86]. To demonstrate the noise reduction more rigorously compared to the intuition provided in Ref. [31], let us focus on the variance of the gradient with respect to a parameter \(\theta\) in the set of the variational parameters \(\mathbf{\theta}\), after subtracting the baseline. Here we focus on the case of a positive ansatz wave function \(\Psi_{\mathbf{\theta}}(\mathbf{\sigma})=\sqrt{P_{\mathbf{\theta}}(\mathbf{\sigma})}\) that we used in our study. First of all, we define:

\[O_{\theta}(\mathbf{\sigma})\equiv\partial_{\theta}\log\left(\Psi_{\mathbf{\theta}}^{*} (\mathbf{\sigma})\right)=\frac{1}{2}\partial_{\theta}\log(P_{\mathbf{\theta}}(\mathbf{ \sigma})).\]

Thus, the gradient with a baseline can be written as:

\[\partial_{\theta}E_{\mathbf{\theta}}=2\bra{O_{\theta}(\mathbf{\sigma})\overline{E}_{ \mathrm{loc}}(\mathbf{\sigma})},\]

where \(\overline{E}_{\mathrm{loc}}(\mathbf{\sigma})\equiv E_{\mathrm{loc}}(\mathbf{\sigma})- E_{\mathbf{\theta}}\) and \(\langle.\rangle\) denotes an expectation value over the Born distribution \(\ket{\Psi_{\mathbf{\theta}}(\mathbf{\sigma})}^{2}\). To estimate the gradients' noise, we look at the variance of the gradient estimator, which can be decomposed as follows:

\[\mathrm{Var}(O_{\theta}\overline{E}_{\mathrm{loc}}) =\mathrm{Var}(O_{\theta}E_{\mathrm{loc}})\] \[-2\mathrm{Cov}(O_{\theta}E_{\mathrm{loc}},O_{\theta}E_{\mathbf{ \theta}})+E_{\mathbf{\theta}}^{2}\mathrm{Var}(O_{\theta}).\]

Thus the variance reduction \(R\), after subtracting the baseline, is given as:

\[R \equiv\mathrm{Var}(O_{\theta}\overline{E}_{\mathrm{loc}})-\mathrm{ Var}(O_{\theta}E_{\mathrm{loc}}),\] \[=-2E_{\mathbf{\theta}}\mathrm{Cov}(O_{\theta}E_{\mathrm{loc}},O_{ \theta})+E_{\mathbf{\theta}}^{2}\mathrm{Var}(O_{\theta}).\]
Since the gradients' magnitude tends to near-zero values close to convergence, statistical errors are more likely to make the VMC optimization more challenging. We focus on this regime for this derivation to show the importance of the baseline in reducing noise. Thus, we assume that \(E_{\text{loc}}(\mathbf{\sigma})=E_{\mathbf{\theta}}+\xi(\mathbf{\sigma})\), where the supremum of the local energies fluctuations is much smaller compared to the variational energy, i.e., \((\sup_{\mathbf{\sigma}}|\xi(\mathbf{\sigma})|)\ll E_{\mathbf{\theta}}\). From this assumption, we can deduce that:

\[R =-2E_{\mathbf{\theta}}^{2}\text{Cov}(O_{\theta},O_{\theta}) \tag{11}\] \[-2E_{\mathbf{\theta}}\text{Cov}(O_{\theta}\xi,O_{\theta})+E_{\mathbf{ \theta}}^{2}\text{Var}(O_{\theta}),\] (12) \[=-E_{\mathbf{\theta}}^{2}\text{Var}(O_{\theta})-2E_{\mathbf{\theta}} \text{Cov}(O_{\theta}\xi,O_{\theta}). \tag{13}\]

The second term can be decomposed as follows:

\[\text{Cov}(O_{\theta}\xi,O_{\theta})=\langle O_{\theta}^{2}\xi \rangle-\langle O_{\theta}\xi\rangle\langle O_{\theta}\rangle. \tag{14}\]

Since \(\langle O_{\theta}\rangle=\frac{1}{2}\langle\partial_{\theta}\log(P_{\mathbf{ \theta}})\rangle=0\)[86; 87; 31], then we can bound the covariance term from above as:

\[\text{Cov}(O_{\theta}\xi,O_{\theta}) \leq\left(\sup_{\mathbf{\sigma}}|\xi(\mathbf{\sigma})|\right)\langle O_{ \theta}^{2}\rangle,\] \[=\left(\sup_{\mathbf{\sigma}}|\xi(\mathbf{\sigma})|\right)\text{Var}(O_{ \theta}),\] \[\ll E_{\mathbf{\theta}}\text{Var}(O_{\theta}).\]

Thus, we can conclude that the variance reduction \(R\) in Eq. (13) is negative. This observation highlights the importance of the baseline in reducing the statistical noise of the energy gradients near convergence.

For a complex ansatz wave function \(\Psi_{\mathbf{\theta}}(\mathbf{\sigma})=\sqrt{P_{\mathbf{\theta}}(\mathbf{\sigma})}\exp(\text {i}\phi_{\mathbf{\theta}}(\mathbf{\sigma}))\), the expectation value \(\langle O_{\theta}\rangle=\frac{1}{2}\langle\partial_{\theta}\log(P_{\mathbf{ \theta}})\rangle-\text{i}\langle\partial_{\theta}\phi_{\mathbf{\theta}}\rangle=- \text{i}\langle\partial_{\theta}\phi_{\mathbf{\theta}}\rangle\) is no longer equal zero in general. We leave the investigation of this case for future studies.

Similarly to the stochastic estimation of the variational energy using our ansatz, we can do the same for the estimation of the variational pseudo free energy \(F_{\mathbf{\theta}}\) in Eq. (4). More details can be found in the supplementary information of Ref. [57]. Finally, we note that the gradient steps in our numerical simulations are performed using Adam optimizer [88].

## Appendix C Hyperparameters

For all models studied in this paper, we note that for each annealing step, we perform \(N_{\text{train}}=5\) gradient steps. Concerning the learning rate \(\eta\), we choose \(\eta=10^{-3}\) during the warmup phase and the annealing phase and we switch to a learning rate \(\eta=10^{-4}\) in the convergence phase. We finally note that we set the number of convergence steps as \(N_{\text{convergence}}=10000\). In Tab. 1, we provide further details about the hyperparameters we choose in our study for the different models. The meaning of each hyperparameter related to annealing is discussed in detail in Refs. [55; 57].

Additionally, we use \(M_{e}=2\times 10^{6}\) samples for the estimation of the entanglement entropy along with their error bars for the toric code. For the Bose-Hubbard model we use \(M_{e}=10^{7}\) samples to reduce the error bars on the TEE in Fig. 5. To estimate the TEE uncertainty from the Kitaev-Preskill construction, we use the standard deviation expression of the sum of independent random variables [89].

Finally, we note that to avoid fine-tuning the learning rate for each value of \(V\) (between 4 and 13) in the Bose-Hubbard model, we target the normalized Hamiltonian

\[\hat{H}=-\frac{1}{V}\sum_{\langle i,j\rangle}\left(b_{i}^{\dagger}b_{j}+b_{i} b_{j}^{\dagger}\right)+\sum_{\bigcirc}n_{\bigcirc}^{2} \tag{15}\]

in our numerical experiments.

## Appendix D Kitaev-Preskill constructions

In this appendix, we provide details about the subregions used to calculate the TEE using the Kitaev-Preskill construction (see Sec. II.2). For the 2D toric code, we use three spins for each subregion and for the Bose-Hubbard model in the Kagome lattice we increase the subregions sizes to mitigate finite size effects [10] as opposed to the 2D toric code that does not suffer from this limitation [11]. The illustrations of these subregions are provided in Fig. 6.

## Appendix E RNNs and MES

The results in Sec. III.1 indicate that the RNN wave function encodes a superposition of minimally entangled states (MES). Here we further investigate this statement by analyzing the expectation values of the average Wilson loop operators and the average 't Hooft loop operators.

We define the average Wilson loop operators as

\[\hat{W}_{d}^{z}=\frac{1}{L}\left(\sum_{\mathcal{C}_{d}}\prod_{\sigma_{j}\in \mathcal{C}_{d}}\hat{\sigma}_{j}^{z}\right). \tag{16}\]

Here \(d=h,v\) and \(\mathcal{C}_{h},\mathcal{C}_{v}\) are closed non-contractible loops illustrated in Fig. 7. A set of degenerate ground states of the toric code are eigenstates of the operators \(\hat{W}_{h}^{z}\), \(\hat{W}_{v}^{z}\) with eigenvalues \(\pm 1\). Additionally, the two eigenvalues uniquely determine the topological sector of the ground state. In this case, the topological ground states can be labeled as \(|\xi_{ab}\rangle\) with \(a,b=0,1\)[63; 90].

We can also define the average 't Hooft loop operators on non-contractible closed loops [90], such that

\[\hat{W}_{d}^{x}=\frac{1}{L}\left(\sum_{\mathcal{C}_{d}}\prod_{\sigma_{j}\in \mathcal{C}_{d}}\hat{\sigma}_{j}^{x}\right), \tag{17}\]
where \(d=h,v\) and \(\tilde{\mathcal{C}}_{h}\), \(\tilde{\mathcal{C}}_{v}\) correspond to horizontal and vertical loops as illustrated in Fig. 7. These operators satisfy the anti-commutation relations \(\{\hat{W}_{h}^{z},\hat{W}_{v}^{x}\}=0\) and \(\{\hat{W}_{v}^{z},\hat{W}_{h}^{x}\}=0\).

From the optimized RNN wave function (\(L=8\)), we find \(\langle\hat{W}_{h}^{z}\rangle=0.0009(2)\) and \(\langle\hat{W}_{v}^{z}\rangle=-0.0039(2)\) which are consistent with vanishing expectation values. We also obtain \(\langle\hat{W}_{h}^{x}\rangle=0.999847(5)\) and \(\langle\hat{W}_{v}^{x}\rangle=0.999785(5)\) for the 't Hooft loop operators, which are consistent with \(+1\) expectation values. These results are in part due to the use of a positive RNN wave function which forces the expectation values \(\langle\hat{W}_{h}^{x}\rangle\) and \(\langle\hat{W}_{v}^{x}\rangle\) to strictly positive values and rules out the possibility to obtain, e.g., \(\langle\hat{W}_{h}^{x}\rangle=-1\).

By expanding the optimized RNN wave function in the \(\ket{\xi_{ab}}\) basis, where \(a,b\) are binary variables, we obtain

\[\ket{\Psi_{\text{RNN}}}\approx\sum_{ab}c_{ab}\ket{\xi_{ab}}.\]

Here \(\{\ket{\xi_{ab}}\}\) correspond to the four topological sectors and they are mutually orthogonal. \(c_{ab}\) are real numbers since we use a positive RNN wave function. Additionally, the basis states \(\ket{\xi_{ab}}\) satisfy

\[\hat{W}_{h}^{z}\ket{\xi_{ab}} =(-1)^{a}\ket{\xi_{ab}},\] \[\hat{W}_{v}^{z}\ket{\xi_{ab}} =(-1)^{b}\ket{\xi_{ab}}.\]

\begin{table}
\begin{tabular}{|c|c|c|} \hline Figures & Parameter & Value \\ \hline \multirow{3}{*}{2D toric code} & Number of memory units & \(d_{h}=60\) \\  & Number of samples & \(M=100\) \\  & Initial pseudo-temperature & \(T_{0}=2\) \\  & Number of annealing steps & \(N_{\text{anneasing}}=4000\) \\ \hline \multirow{3}{*}{Bose-Hubbard model (L = 6)} & Number of memory units & \(d_{h}=100\) \\  & Number of samples & \(M=500\) \\ \cline{1-1}  & Pseudo-temperature & \(T_{0}=0\) \\ \cline{1-1}  & Number of steps & \(10000\) \\ \hline \end{tabular}
\end{table}
Table 1: A summary of the hyperparameters used to obtain the results reported in this paper. Note that the number of samples \(M\) corresponds to the batch size used during the training phase. Additionally, the \(L=8\) RNN model is initialized with the \(L=6\) optimized RNN for the Bose-Hubbard model.

Figure 6: An illustration of the sub-regions \(A,B,C\) chosen for the Kitaev-Preskill constructions. (a) For the 2D toric code, each sub-region has 3 spins. We also do the same for the lattice sizes \(L=6,8,10\) in Fig. 4(b). For the hard-core Bose-Hubbard model on the Kagome lattice, we target two different system sizes. Panel (b) shows the construction for L = 6 and panel (c) provides the construction for \(L=8\). In panels (b) and (c), each site corresponds to a block of three bosons.

Figure 7: An illustration of the vertical and the horizontal loops used to compute the Wilson loop operators (see Eq. 1) and the ’t Hooft loop operators (see Eq. 2).
From the anti-commutation relations, we can show that:

\[\hat{W}_{h}^{x}\ket{\xi_{ab}} =\ket{\xi_{\text{a}\!\text{b}}},\] \[\hat{W}_{v}^{x}\ket{\xi_{ab}} =\ket{\xi_{\text{a}\!\text{b}}},\]

where \(\bar{a}=1-a\) and \(\bar{b}=1-b\). By plugging the last two equations in the \(\hat{W}_{h}^{x}\), \(\hat{W}_{v}^{x}\) expectation values of our optimized RNN wave function, we obtain:

\[2c_{00}c_{01}+2c_{10}c_{11} \approx 1,\] \[2c_{00}c_{10}+2c_{01}c_{11} \approx 1.\]

From the normalization constraint \(1=\sum_{ab}c_{ab}^{2}\), we deduce that:

\[(c_{00}-c_{01})^{2}+(c_{10}-c_{11})^{2} \approx 0,\] \[(c_{00}-c_{10})^{2}+(c_{01}-c_{11})^{2} \approx 0.\]

As a consequence, we conclude that \(c_{00}\approx c_{01}\approx c_{10}\approx c_{11}\), which means that the optimized RNN wave function is approximately a uniform superposition of the four topological ground states \(\ket{\xi_{ab}}\). This observation is also consistent with vanishing expectation values of the operators \(\hat{W}_{h}^{z}\), \(\hat{W}_{v}^{z}\).

Furthermore, from Ref. [63] the MES of the toric code are given as follows:

\[\ket{\Xi_{1}} =\frac{1}{\sqrt{2}}\left(\ket{\xi_{00}}+\ket{\xi_{01}}\right),\] \[\ket{\Xi_{2}} =\frac{1}{\sqrt{2}}\left(\ket{\xi_{00}}-\ket{\xi_{01}}\right),\] \[\ket{\Xi_{3}} =\frac{1}{\sqrt{2}}\left(\ket{\xi_{10}}+\ket{\xi_{11}}\right),\] \[\ket{\Xi_{4}} =\frac{1}{\sqrt{2}}\left(\ket{\xi_{10}}-\ket{\xi_{11}}\right).\]

Thus, our RNN wave function can be written approximately as a uniform superposition of the MES \(\ket{\Xi_{1}}\) and \(\ket{\Xi_{3}}\), i.e.

\[\ket{\Psi_{\text{RNN}}}\approx\frac{1}{\sqrt{2}}\left(\ket{\Xi_{1}}+\ket{\Xi _{3}}\right).\]