Restricted Boltzmann machine representation for the groundstate and excited states of Kitaev Honeycomb model

In this work, the capability of restricted Boltzmann machines (RBMs) to find solutions for the Kitaev honeycomb model with periodic boundary conditions is investigated. The measured groundstate (GS) energy of the system is compared and, for small lattice sizes (e.g. $3 \times 3$ with $18$ spinors), shown to agree with the analytically derived value of the energy up to a deviation of $0.09\%$. Moreover, the wave-functions we find have $99.89\%$ overlap with the exact ground state wave-functions. Furthermore, the possibility of realizing anyons in the RBM is discussed and an algorithm is given to build these anyonic excitations and braid them for possible future applications in quantum computation. Using the correspondence between topological field theories in (2+1)d and 2d CFTs, we propose an identification between our RBM states with the Moore-Read state and conformal blocks of the $2$d Ising model.


I. INTRODUCTION
The Honeycomb model we consider in this work is a 2-dimensional lattice spin system (see Fig. 1) which was studied in detail by Kitaev in 2006 [1] and is famous because of the topological quantum order due to a degenerate gapped groundstate which is persistent to local and finite-sized perturbations. The system also supports both Abelian and non-Abelian topological phases as demonstrated in the original proposal [1]. There is a wide range of applications for this model, from fault-tolerant quantum computation [2] to analytical study of strongly correlated systems [3] and quantum spin liquids [4].
The honeycomb lattice is not a Bravias lattice in its original structure, but it can be considered as a triangular Bravias lattice with a two-spin basis. The direct Bravias lattice and the primitive cells are illustrated in Fig. 1a with the dashed lines. Each primitive cell has a pair of odd and even indexed spins. Just for the sake of simplicity, we will denote each cell with the even indexed (empty circles) spins of the system. The primitive vectors of the lattice, a 1 and a 2 are also shown in Fig. 1a (see also Eq. (1)). The entire lattice can be tiled and covered with primitive cells using translations composed of different linear combinations of primitive vectors where a is the lattice constant. Basis vectors b 1 and b 2 for the reciprocal lattice can be obtained by solving Eq.
(2) * Corresponding author's email: babakhaghighat@tsinghua.edu.cn The solution is as below The reciprocal lattice is depicted in Fig. 1b. The Hamiltonian of the Kitaev model is a nearestneighbor interaction of Pauli matrices on a honeycomb lattice as written in Eq. (4), where r&r are indices for the nearest neighbour spins. The physics of the system is symmetric under permutation of coupling strengths J α with α = x, y, z and due to an-isotropic interaction, the model is a frustrated spin system, because a spin cannot satisfy conflicting demands of orientation from its three neighboring sites [1] Of central interest in this work is the property of the Kitaev model to support non-Abelian chiral spin liquid (CSL) states. A very symmetric way to realize this phase of the model is by setting J x = J y = J z = 1 and turning on a small magnetic field which gives rise to a further interaction term in the Hamiltonian (see Eq. (47)). In the present work, we will restrict ourselves to this non-Abelian phase but set -with the exception of Section V Cthe magnetic field to zero. When CSL systems are placed on topologically non-trivial Riemann surfaces with nonzero genus, one obtains multiple degenerate groundstates that are separated from the rest of the spectrum by an energy gap. This allows to describe the low-energy physics of such models using topological quantum field theories (TQFTs) (see for example [5] for a review). The non-Abelian phase of the Kitaev model is in the universality class of the Ising TQFT. When placed on the torus, which corresponds to periodic boundary conditions of the lattice and is the choice employed in this work, this results in a triply-degenerate groundstate. These 3 groundstates correspond to the 3 different ways to obtain a torus by  The generators of the lattice in position space are named e 1 and e 2 , as shown in the image above. Hence, the lattice points are L = n 1 e 1 + n 2 e 2 with n 1 , n 2 ∈ Z and the above lattice is drawn for N 1 = N 2 = 4 with toric boundary conditions. gluing the two ends of a cylinder in the Ising TQFT, namely by inserting the operators 1, σ and ψ at the two ends. For the topological phase discussed above, the Kitaev model has recently been numerically solved in the spin basis to a high degree of approximation using tensor networks [6,7]. In this paper we employ a different method based on the approach developed in [8] using RBMs to find groundstates and excited states of the Honeycomb model. For a different application of RBMs to the Kitaev model, see [9]. The RBM as a variational ansatz to find groundstates of spin lattice models was first proposed in [10], where the formalism is applied to the 1d transverse-field Ising model and the 2d anti-ferromagnetic Heisenberg model. The authors of [11] use RBMs and feed-forward neural networks to study excited states of the 1d Heisenberg model and the 1d Bose-Hubbard model. Further applications of RBMs in the literature can be found in [12][13][14][15][16]. An application of RBMs to systems with non-Abelian symmetries was presented in [17], where the authors sample in the basis of irreducible representations instead of individual spins.
In this work, the central focus is on the RBM as a variational ansatz to find accurate groundstate wave-functions for the 3 × 3 Kitaev Honeycomb lattice as this is the smallest lattice size where the model admits topological order resulting in a triply-degenerate groundstate. Along the way, we also use the techniques developed to solve accurately for groundstates of smaller lattice sizes and explore the accuracy for larger lattices as well as the possibility to realize excited states. It is, however, not the goal of this paper to develop a machine learning architecture which solves efficiently and accurately for groundstates of large lattice systems. As explained in the main body of the paper, such an endeavor would need the simultaneous deployment of several distinct techniques employed and is left for future work.
Apart from a variational ansatz to find groundstates, another deep insight into chiral topological order is provided by the work of Moore and Read [18], where the authors find that many wave-functions can be expressed as conformal blocks. Conformal blocks are chiral correlation functions in certain conformal field theories (CFTs) which are also known to arise as states in the Hilbert space of the 2 + 1-dimensional topological Chern-Simons theory [19]. This connection has been generalized to describe chiral topological states in lattice systems [20]. In our work, we explore the possibility to realize the 3 groundstates of the Kitaev model on the torus in terms of conformal blocks of the SU (2) 2 WZW CFT and show that our ansatz respects the symmetries of the Kitaev Hamiltonian (4). The consequences of the identification of the RBM states with the conformal black ansatz is then explored in the conclusions.
The outline of this paper is as follows. In Section II, we briefly review the analytical solution of the model. In Section III the mapping of the spin model to an RBM architecture is explained. In particular, we demonstrate how to calculate the groundstate and excited states of the Honeycomb model using machine learning techniques. To this end we employ both the NetKet environment [21] as well our own self-built package using PyTorch to train a restricted Boltzmann machine (RBM) in order to find a groundstate via gradient descent. In Section IV realization of vortices is discussed from both analytical and RBM points of view. In brief, a vortex pair can be created by adjusting the parameters of the RBM and subsequent braiding can be achieved by further adjustments of parameters. Section V is dedicated to results and relevant discussions and the reconstruction of a given quantum state from simple measurements known as quantum state tomography [22]. In Section VI we propose an identification of our RBM states with the Moore-Read states on the torus which can be computed using conformal blocks of the 2d Ising model. Finally we conclude in Section VII.

II. SOLUTION OF THE MODEL
In this section an exact solution of the system is provided and an explicit form of the groundstate energy is derived where we will be following references [1,23,24]. First of all, we fermionize the Hamiltonian by performing a one-dimensional Jordan-Wigner transformation [25] defined in Eq. (5). A deformation of the hexagonal lattice to a brick-wall lattice (see Fig. 2) clarifies the mechanism of the Jordan-Wigner transformation and why it is onedimensional [26,27]. In the brick-wall lattice, each lattice site r is denoted by the coordinates (i, j).
This transformation maps the Hilbert space of spins to the Hilbert space of spinless complex fermions. Using the fact that σ ± = σ x ± iσ y one can expand the Hamiltonian in Eq. (4) into three terms and rewrite them as below [24] Therefore, the Hamiltonian transforms to The J x and J y terms are quadratic interactions in spinless fermions and are easy to solve, however, the J z term is a product of number density operators and can be further simplified by introducing the Majorana operators in Eq. (8) [24]. As we will see, this simplification can be done due to the presence of the conserved quantity called plaquette operator B p [24].
Then, the J z term can be rewritten using the Majorana operators Finally, using the circle indices for the odd and even lattice sites as defined in Eq. (8), the Hamiltonian transforms to the expression in Eq. (11).
If we write the Hamiltonian as a sum over unit cells we have where r is the position vector of z-bonds or unit cells (see Fig. 2). It makes no difference in the physics if we set r to be any point along the z-bond, whether we choose an even site or an odd site or some point in the middle, they are all related through translation. What will be important in what follows, is that now since we have grouped a pair of even and odd sites as an unit cell, to evaluate the summation over the unit cells, it is sufficient to just run over the even or odd sites. Furthermore, the α r = (id •,r d •,r ) operators are defined on each z-bond of the lattice (labeled by r) and they commute with the Hamiltonian and are good quantum numbers. Moreover, it can easily be shown that each plaquette operator can be written as On the other hand, from Lieb's theorem [28] we know that the groundstate manifold is obtained by setting B p = 1, ∀p. Thus the uniform choice of α r = 1, ∀r corresponds to a vortex-free sector, nevertheless all configurations leading to the same sector are equivalent. We also introduce a Dirac fermion on each z-link using the Majorana operators Using the inverse transformation we can rewrite the Hamiltonian as The Hamiltonian is translation invariant and can be transformed to momentum space in order to be diagonalized.
We define the Fourier transform for the Dirac fermion through Setting k → -k in d † r simplifies the calculations. We then have Using all previous transformations we obtain a Hamiltonian which is quadratic in the Dirac fermions living on each of the unit cells along the z-links where γ † k and γ k are quasiparticle creation and annihilation operators. The groundstate energy can easily be read off from the above Hamiltonian and is given by

A. RBM as a Neural Network Quantum State
With the ever growing applications of neural networks in sciences and the emergent new technologies to deploy and build physical neural networks, this is the right time to investigate potential applications of them in condensed matter systems. Recently, a new approach has been proposed for simulating topological quantum states using neural networks [8]. In this section we use the same approach and use a restricted Boltzmann machine (RBM) as a variational ansatz for the Honeycomb Kitaev model. There are many reasons why an RBM is chosen. This particular architecture has proved to be effective in many tasks such as dimensional reduction, classification, regression, collaborative filtering, feature learning and topic modeling and moreover a general theorem shows that RBMs are universal approximators of discrete distributions [29].
Having a set of spins on a lattice, Ξ = (σ 1 , σ 2 , · · · σ N ) (in our case the Jordan-Wigner chain of spins), our goal is to use an RBM to reduce the dimensionality of the Hilbert space and estimate the energy of the system in both groundstate and excited states to be able to classify different topological phases of the model. An RBM consists of two layers. The first layer (visible layer) has N nodes which are representing the physical spins in the Hamiltonian and the second layer (hidden layer) has M binary valued nodes (architecture of the network is shown in Fig. 4). The quantum state of the Honeycomb model, up to an irrelevant normalization factor, can be written as where To obtain the second equality we used the values for {h k } = {−1, 1} M which is the set of possible configurations of hidden layer nodes and Ω = (a k , b k , W kk ) is the set of weights and biases of the RBM which should be trained in such a way that the final RBM state represents the desired quantum state of the model (i.e. groundstate or the excited states), where we refer to [21], sections 2.2.3 and 2.2.4, for further details of the variational principle. The combined number of weight and bias parameters and nodes is polynomial in system size and computationally feasible.
The variational calculations for the model with RBM as ansatz is done twofold. One approach we follow is to apply the machinery already developed in the NetKet software package [21], while the second approach we employ is to develop our own package using PyTorch. In the case of NetKet, in order to train the network, the parameter space of the network is sampled using the Metropolis algorithm and the optimization iterations are based on stochastic gradient descent. It turns out that although the sampling algorithm is effective, the ground state energy cannot be reached to a satisfactory precision. The reasons will be explained as we proceed. To circumvent this problem, we develop our own software package using PyTorch which is designed tightly to suit the purposes of the Honeycomb model and avoid possible pitfalls.

B. Quantum State Tomography
Reconstruction of complex synthetic quantum states from experimental measurement data is computationally exponentially expensive. Therefore, the capability of neural network quantum states such as RBM to accurately and efficiently represent high-dimensional quantum states is beneficial to quantum state tomography. Apart from the numerous applications of QST, in our study case, the reconstructed states of the honeycomb model can later on be used for validation of particular quantum processes, e.g. braiding of vortices in the system, phase transitions, etc. at the experimental level and after implementation of the framework developed in Section IV. For instance, there are promising proposals for implementation of topological systems in cold atoms. A review has recently been published here [30].
In Section V D we present QST results for the groundstate of a small system size of the honeycomb model, as proof of concept. In a normal QST task, the data sets come from experimental measurements, but since we don't have experimental data, we need to first find the exact wave-function for the system through exact diagonalization methods and then sample that wave-function by performing single shot measurements in different bases to build the data set.

IV. BRAIDING
In this section we want to explore the possibility to create anyonic excitations in the RBM representation and describe anyon braiding. Such a representation is important for various reasons, in particular for performing quantum computation [5].
First of all, let us describe the creation of these quasiparticles by starting with a groundstate sector and then subsequently changing the eigenvalue of the B p operators from +1 to -1 locally and in the region where we want to realize the vortices in the system. In an arbitrary configuration of spins (not necessarily in the groundstate sector) it is easy to prove that B p = +1, ∀p. Hence, we can just build the vortices in pairs. For example, if you apply the operatorÔ 1 to the groundstate, it produces two vortices in plaquettes 1 and 2 (see Fig. 5) [31].
Another possibility is to apply the following operator which produces two vortices along theẑ direction in plaquettes 3 and 4. Let us explain how these operators create vortices from the viewpoint of Kitaev's original Majorana fermion representation [1]. We assume that there are two fermionic modes living on each lattice site corresponding to the four creation and annihilation operators a † m,i and a m,i , where m ∈ {1, 2} is the index for the modes and i is the index for the lattice sites. Here unlike before, we use one index for lattice site to avoid unnecessary complications. Then we will separate the imaginary and real parts of these operators, to define the Majorana fermions as below Notice that the spins have a two dimensional space (being up or down) and now that we are representing them with four Majorana modes (two complex fermionic modes), we need to project out the unphysical states. The Fock space of the complex fermionic modes can be represented as {|00 , |01 , |10 , |11 }. We make the following correspondence between the states of spins and fermions One can define the projector P i on site i [32] to do this job for us One can easily show that the relation between Majorana operators and the original Pauli operators is as follows In fact, by the above projection, we satisfied the extra condition coming from the algebra of Pauli matrices, Therefore, the eigenvalue of the projector P i for physical states is 1 and for unphysical states is 0.
Using Eq. (28) the Hamiltonian can be written in terms of Majorana operators as The u ij are antisymmetric Hermitian link operators with eigenvalues ±1 The link operators commute with the Hamiltonian, [H, u ij ] = 0, so they are local symmetries. In this form, the Hamiltonian is representing a tight binding model of free Majorana fermions hopping along the lattice, with tunneling couplings that depend on the eigenvalues of link operators which can be thought of as a classical Z 2 gauge field. One can assign a configuration {u} to all link operators and diagonalize the resulting quadratic Hamiltonian in Majorana operators directly and obtain the spectrum and groundstate energy for H{u}. But which configuration {u} corresponds to the global minimum of the energy? This question has been answered by Lieb (1994) [28]. The groundstate lies in the sector with B p = 1 ∀ p and since the B p operators are the only gauge invariant objects, every gauge fixing choice {u} which leads to this particular configuration for B p is equivalent. Therefore, our goal is to find a configuration {u} which leads to B p = 1 ∀ p. Next, we define the vortex operatorŴ p aŝ The vortex operator can be simplified as beloŵ where we used ic x j c y j = −σ z j D j and cyclic permutations.
Therefore, the uniform choice of u ij = 1 for every link will lead us to the groundstate in the extended (not yet projected) space with B p = 1. In order to find the physical groundstate, we apply the global projector D to project out the unphysical states Now we are at a point to explain in detail why the operatorsÔ 1 andÔ 2 can create vortices in the system. As mentioned earlier, if we flip the sign of the eigenvalue of the vortex operator at p, a vortex will be created in the system. Looking back at Eq. (31), we see that the vortex configuration {Ŵ p } is created by fixing the gauge at every link {u ij }. Therefore, by changing the link operators, we can move between vortex sectors. Starting from a vortexfree sector with all link operators being +1, applyingÔ 1 on site a will leave the σ z a untouched while flipping the σ x a and σ y a . Consequently, this will change the sign of the eigenvalue of the u ij for x-link and y-link attached to site a, leaving the system with B p1 = −1 and B p2 = −1. We can also move the vortices around in the system, by consecutively applying these operators on the particular sites in a path for which we want to move the vortices along. As an example, the red path shown in the Fig. 6 corresponds to the operator In general, we can define string operators to move the vortices anywhere in the system along the string, similar to the example above. To create such a string FIG. 6: Two vortices connected by a string operator passing through several links. The white dots denote even lattice sites while the black dots are odd ones. By convention, the string operator is built up from operators acting on even sites.
operator is straightforward now. If the string passes through a link, we need to apply an operator on one of the sites connected to that α − link which carries σ α in the exponent. Before continuing, we need to make a convention. A sharp reader might point out that applying the operatorÔ = exp (−i π 2 σ z l ) exp (−i π 2 σ x k ) instead of O creates the same pair of vortices in the system. For consistency reasons, we therefore use the convention that if a string is passing through a link, the corresponding string operator is built up by picking the even site on the link. Having said all of this, a generic string operator is defined asŜ This string operator will create two vortices at the ending points of the string which passes through the α−links. The manipulation of link operators using these string operators is practically equivalent to changing the sign of the couplings J ij corresponding to the link u ij in the definition of the A ij (see Eq. (29)). Therefore, in the RBM language one can simply redefine the Hamiltonian with the desired signs of J ij for a particular vortex sector and then try to find the RBM representation for this new Hamiltonian through the approach explained in Section III. The only practical disadvantage of this vortex realization approach is that it is computationally expensive as the RBM should be trained again from scratch to find the representation for the excited state. An ideal approach is one where, given the representation for the groundstate, we would be able to directly and without further optimization steps find the representation for the excited state.
In principle, it is possible to directly manipulate the weight and bias parameters of the RBM, Ω = (a k , b k , W kk ) to realize the vortices and find the excited state representation. Looking at Eq. (22), we notice that there is a symmetry between the eigenvalues of physical spins and the value of parameters of the RBM. In other words, flipping a particular spin in site k is equivalent to flipping the sign of the RBM parameter sitting next to the σ α k in Eq. (22) and/or adding a phase to it. In general, for a string operatorŜ α l = exp (−i π 2σ α l ) we havê Next, one needs to evaluate how (−iσ α a ) acts on |Ξ . Here we subdivide the three cases α = z, x and y. If α = z, spin up states will get a phase −i and spin down states a phase i. The overall effect of this can be absorbed in a shift of a l by −i π 2 . If α = x, then spins up and down are just exchanged and the −i is an overall phase which can be ignored. If α = y, we get a combination of the previous two cases. In equation form, the above explanations take the following form V. RESULTS AND DISCUSSION

A. Groundstate
In this section, we present the results for the groundstate energy estimate of different system sizes for the proposed mapping discussed in Section III. These results are also compared to the exact analytical expression given in Eq. (20). The absolute value of the energies is plotted in Fig. 11 in Appendix C. All results for this and the following sections are for the phase J x = J y = J z = 1.
We used both NetKet and PyTorch to build our RBM. While NetKet provides build-in functions to evaluate physical observables, we build our own version of these functions in the PyTorch environment. When training with NetKet, the configuration space of the network is sampled for hundreds of times, depending on the lattice size, using the Metropolis algorithm. While when training with Py-Torch, the configuration space of the network is sampled from all of the Hilbert space. This is easier to implement and enables us to focus on the explicit RBM representation. The loss function of the network has been minimized using Stochastic Gradient Descent and the learning rates are set to decay using the Adaptive Moment (ADAM) estimation algorithm. The groundstate energies are obtained by fine-tuning the optimisation algorithm parameters.
The numerical results obtained are shown in Fig. 11 in Appendix C. The different methods for computing the groundstate energy listed in the Table I  Based on the results obtained for the energy of the groundstate of the system, training of the network is tested and optimized by changing the learning rate for different training phases and iteration numbers. The most effective initial learning rate is found to be lr init ≈ 0.01. This value is small enough to avoid over-fitting and large enough to train the network effectively. This learning rate is rather large compared to standard machine learning implementations, because the number of parameters in our RBM is relatively small compared to hundreds of millions of parameters used in typical deep learning applications.
A central ingredient of our approach in PyTorch is to use projection onto translation invariant states in order to restrict RBM parameters to the physical parameter space. This dramatically improves the results we can achieve. The theoretical reason for this projection is that ground states of our Hamiltonian reside in such translation invariant subspaces of the Hilbert space as explained in Section VI. Here we provide the technical details of the implementation. Suppose |Φ is the RBM state mentioned in Eq. (21). Then the state used in calculating the energy is as follows a Notice that the energies obtained via the analytical formula do not match the exact diagonalization results for 2 × 2 and 2 × 4 lattices. This is not a problem because the energy obtained via the analytical formula (method a.i) is still in the spectrum, which means this is a case where Lieb's theorem fails -the vortex free sector no longer contains the ground state. This is an edge case and will not appear when m, n > 2. whereT ma1+na2 are operators translating m steps along the a 1 direction and n steps along the a 2 direction, with a 1 , a 2 defined in Eq. (1).T ma1+na2 can be realized either as a permutation matrix acting on quantum states or a permutation of parameters in RBM. The hyperbolic trigonometric functions appearing in the weights of the RBM are numerically unstable. They sometimes lead to divergence and we have to use an extremely small learning rate to avoid the divergence. To cure this instability, we use the following pre-train algorithm. We first perform the substitution as the polynomials x and x 2 are more stable under numerical changes than the exponential e x which appears in the hyperbolic trigonometric functions. Then this modified model is trained for several thousand epochs. After we get close enough to the minimum, we change back to the standard RBM representation using cosh and sinh to refine our result. We find that this two-step training shortens the time to reach a given accuracy. Remember that RBM parameters are complex numbers. However, the meaning of the real part and the imaginary part of parameters are different. The real part contributes by cosh / sinh while imaginary parts show up as arguments of cos / sin, i.e. they are phases. So they should be initialized with different standard deviations in the sampling distribution and trained with different learning rates. In other words they are "renormalized" differently along the flow.
On the other hand, when working with NetKet, we find that the ground state energy cannot be reached beyond a certain accuracy. Thus the results obtained with our PyTorch algorithm described above are always far more accurate than the best results NetKet can achieve. We believe this is due to the following identity used by NetKet, in order to accelerate sequential multiplication. This is theoretically correct but numerically unstable as in an RBM representation there are always extremely large or small values of cosh / sinh (this is the means by which an RBM is capable of approximating the correct physical state!) leading to truncation errors when using the functions log and exp. So it is not surprising that NetKet cannot reach a high enough accuracy and gets stuck at a too high energy. We find that when using the identity (40) in our PyTorch implementation, our ground state energy becomes similarly inaccurate as the one obtained by NetKet.
In order to check the effect of the relative number of hidden and visible nodes on the outcome of the network, different values for α were tested, keeping other parameters fixed. The α parameter was changed just by changing the number of hidden nodes for a fixed system size (and hence fixed visible layer size). The results (for a 5 × 5 lattice) are summarized in the Table IV in Appendix C.

-The performance of FFNN vs. no-visible-bias RBM
In Table V in Appendix C, we give a brief summary of the comparison results between an FFNN with one layer and a no-visible-bias RBM (referred to as novb-RBM in the table). The two architectures are similar and for the same computational resources (such as CPU and memory) and system size these results have been obtained. The conclusions we reach are: 1. The FFNN trains faster than the RBM with relatively fewer free parameters and performs better than the RBM with no bias.

But the FFNN's result deviates by a greater margin from the analytic result.
-RBM's parameters: an example To get a better feeling for how the RBM approximates quantum states so well, let us have a closer look at an RBM for a small system size, as shown in Eq. (41). These are the parameters for an RBM with 5 hidden nodes for a 2×2 system which has 8 spinors. So the parameter matrix W kk is 5 × 8 in size. The indexing rule for the spinors is shown in Fig. 9. The biases a k and b k in RBM are set to 0 manually. For illustration purposes, this RBM is trained without the projection method. After training, the energy of the resulting parameter set is −6.928201 while the true ground state energy is −6.92820323. Its wave-function has 99.999982% overlap with the groundstate.
To begin with, notice that, for example, for a state with all spins up, the RBM amplitude is given by where we are using the convention that spin up corresponds to the value σ z = 1 and spin down to σ z = −1.
If we focus on the factor corresponding to k = 5 and use the identification π/4 ≈ 0.7854, we see from Eq. (41) that its amplitude is given by On the other hand, if we flip one of the spins, the corresponding amplitude becomes That means the second and fifth line of Im(W ) impose strong restrictions on spinor configurations and in this particular case imply that up and down spins are paired in even numbers. This restriction also happens in the first, third and fourth line of Im(W ), although this is not apparent from first sight. This can be seen from observing that Furthermore, there are many numerical coincidences for parameters between spinor pairs (0, 7) (1, 6) (2, 5) and (3,4). This reflects the "parity" symmetry of the 2 × 2 lattice. More specifically, the 2 × 2 lattice can be drawn as FIG. 7: 2 × 2 honeycomb lattice can be drawn as a cubic.
The indexing rule for the spinors is shown in Fig. 9.
a cubic as shown in Fig. 7. Then it is apparent that there are strong correlations between the pairs (0, 7) (1, 6) (2, 5) and (3,4). This small example shows that an RBM is able to find restrictions and symmetries of our quantum system efficiently. We can imagine that this will also happen for larger systems, although it is hard to observe such symmetries explicitly since the number of parameters becomes too large. c) transforming into auxiliary fermion representation and fixing u ij s [2,23].
We applied all three methods to a 3x3 lattice and computed the corresponding energy levels for excited states.
The results we get are shown in Table II.  We find that method (a) is more stable than method (b), which might be due to the random initialization of RBM parameters in method (b) for the training process. Both methods, (a) and (b), deviate by a margin from the analytic values, namely method (c). This can be explained by the fact that the state generated by the spin flip operation is not the groundstate of the vortex sector but rather an excited state carrying quasi-particles. We will comment on this further in the next sub-sections about the expectation values of plaquette operators.

C. Perturbation with non-trivial magnetic field
As originally discussed by Kitaev [1], one can add the following perturbation to the Hamiltonian (4) which corresponds to an effective magnetic field with coupling K, and the sum is over all vertices with links (i, j, k). This term leads to an energy gap in the spectrum for asymptotically large lattices and renders the vortex sectors of the model topological. In this phase, the vortices of the honeycomb lattice bind localized Majorana modes and behave as Ising anyons. For the 3 × 3 lattice, we find that the gap in the spectrum separating the groundstates from excited states becomes larger for non-zero K. For K = 0.2, exact diagonalization gives a GS energy of −17.5260. We find that using our RBM approach, with the same number of parameters and training procedures discussed above, we get an energy value of −17.393 (0.75% deviation) which is quite close to the analytical value, as shown in Table I. This shows that the RBM can be an efficient and accurate representation for the non-trivial topological phases of the Kitaev honeycomb model. As in this work we are mainly interested in the groundstate, which already displays topological order for the 3 × 3 lattice, we leave further work on the phase with nontrivial K for the future.

D. Quantum State Tomography for Groundstate
In this section we apply the technique discussed in Section III B to reconstruct the wave function from mea-surements . In order to generate measurement data set, we sample from an exact groundstate by random measurements in different basis. To this end, we focus on the 2 × 2 lattice in order to deal with the computational resource requests and also being able to find the exact groundstate wave function through exact diagonalization. The last step is to train the RBM using the data set obtained. The result is shown in Fig. 8. There are many intricate parameters which one can adjust during training, such as learning rate, batch size, number of Metropolis samples, and number of single measurements. Only with a suitable choice of parameters and large enough training data set, can one achieve a satisfactory result. For the example depicted in Fig. 8, in order to build the training data set, the exact wave function was sampled 2000 times for 256 randomly chosen basis x, y, z and the hidden node density α was set equal to 2. With this combinations of parameters, QST leads to an overlap of 97 percent between the exact wave function and the learned one.

-Plaquette Operator expectation values and spin flipping
Starting from the state obtained from tomography as discussed above, one can then realize vortices in the system by modifying parameters of the RBM as described in Section IV. In order to verify that modification of RBM parameters takes us to the right excited state, we measure the relevant plaquette operators in the system. The convention for labeling plaquette operators and sites is as shown in Fig. 9.
For the trained RBM state after performing QST for groundstate of 2 × 2 lattice, the expectation values of its four plaquette operators as well as the measured energies are shown in Table III equivalently, applyingŜ x 5 ) will produce a pair of vortices at B and D plaquettes. Therefore, the expectation value of theŴ p operator in these plaquettes also flips sign. Similarly, applyingŜ y 7 , will kill the vortex in plaquette D and produce one in plaquette A. The changes in the sign of expectation values for these plaquettes is shown in the Table III. Finally, by applyingŜ z 3 we return back to the no-vortex sector. These results, prove the effectiveness of the framework developed in Section IV for realization of vortices in RBM states, in both QST and usual learning approaches.
Note that the energy after each flip changes by an amount ∆E ∼ 2. This can be explained by the fact that each vortex sector has its own groundstate but also excited states with higher energies created by fermionic quasi-particle excitations. For example, the no-vortex sector has a spectrum given by the Hamiltonian in Eq. (20) and while the groundstate is the lowest lying state in this sector, there will be excited quasi-particle states carrying momentum k which are created by the γ † k operators. Similarly, vortex sectors with a non-trivial vortex number can be in an excited state by applying quasi-particle creation operators. A rough estimate shows that the difference ∆E ∼ 2 is of the order of several such a quasi-particle excitations and thus the vortices we are creating carry such excitations. This also explains why when returning to the no-vortex sector we do not arrive at the original energy.

VI. RELATION TO CONFORMAL BLOCKS AND MOORE-READ STATE
As originally described in [1], the non-Abelian sector of the Kitaev Honeycomb model which we are addressing in the current paper is in the same universality class as the Ising conformal field theory. That is the topological excitations of the model correspond to the primary operators of the chiral sector of the Ising CFT, namely the vacuum corresponds to the identity operator 1, the vortex corresponds to the σ-field and the quasi-particle corresponds to the ψ-field. This relation is part of a broader framework connecting topological phases of 2 + 1-dimensional quantum field theories to chiral sectors of WZW conformal field theories that can be traced back to the work of Witten on Chern-Simons theory and WZW models [19,33], see [5] for a review. Within this correspondence wave-functions of the 2 + 1-dimensional topological theory are mapped to conformal blocks of the WZW model and vice versa. One very well-known application of this paradigm is the description of fractional Quantum Hall states with filling fraction using conformal blocks in an SU (p) k WZW model [34], where p and k are positive integers, and n is an arbitrary positive integer. For p = 1, these are simply the Laughlin states. The case of Ising anyons corresponds to the choice p = k = 2 and leads to the famous Moore-Read state [18].
Let us now describe this map in some degree of detail which is suitable for our purposes. The ground state wave-function of the topological theory on the 2-sphere is described by the identification [35] where in order to obtain a unique ground state, the con- is computed in the representation R = Sym k , namely the kth symmetric repre-sentation of SU (p), which has the property In the case of quantum hall states, l B = /eB and the z i are positions of "electrons". Specifying to the case p = k = 2, we obtain the Pfaffian state [35] Ψ Pf 0 ({z i }) = Pf (51) where we have taken O R to be the spin 1 field ψ of the Ising model and the Pfaffian of a 2m-by-2m antisymmetric matrix is defined as The state |ij 1 appearing in the Pfaffian is a spin singlet formed from two spin 1 states We will see the significance of this choice shortly. The non-Abelian sector of the Kitaev model corresponds to the value n = 0 as quasi-particles carry the same electric charge as an electron With these choices equation Eq. (51) is a spin singlet version of the Moore Read state [18]. Next, we would like to interpret Eq. (51) as a ground state wave-function of the Kitaev model. That such a correspondence is possible for spin lattice models has been previously shown in the case of Laughlin spin-liquid states on lattices [20,36]. Passing over to the torus with periodic boundary conditions, which is the situation examined in the present paper, the Pfaffian state lifts by the replacements to the following expression involving theta functions [37], where the modular parameter of the theta-functions is chosen to be the one of the lattice 1 , namely in our case τ = e πi/3 . The center of mass contribution has been modified in order to accommodate for the fact that ν = 1 and, using N s = L x L y /2πl 2 B , takes the form and we refer the reader to Appendix B for conventions and definitions of the theta functions. Note that on the torus, the ground state is triply degenerate with the different states corresponding to the values α = (1/2, 0), (0, 0), and (0, 1/2). The transition between these ground states can be understood as creation of vortex/anti-vortex pairs which are then successively transported along the cycles of the torus and then annihilated, as beautifully described in [39]. Moreover, for B = 0 we have l B → ∞ and N s = 0, resulting in periodic boundary conditions for the wave-function. Now let us make the connection between what we have said above about the Moore-Read state and the RBM ground state wave-functions constructed in this paper. First of all, note that the RBM wave-function is in the spin basis of the electrons at each lattice site. Recall that, as discussed in Section IV, the Hilbert space of our spin 1 2 electrons can be described by two spinless complex fermions a 1 and a 2 . In particular, spin up will correspond to both modes unoccupied, while spin down arises from both modes occupied: with a 1 |00 = a 2 |00 = 0 and |11 = a † 1 a † 2 |00 . However, we remark here that | ↑ and | ↓ are not the standard spin up and down basis vectors with respect to σ z , but rather rotated basis vectors which respect the symmetries of the Honeycomb Hamiltonian which for J x = J y = J z = 1 is invariant under permutation of x, y and z. In particular, we make the choice of basis satisfying where α = z, x, y. A concrete solution in the σ z -basis satisfying the above constrains is given by (60) The representation (58) is faithful if we restrict to the subspace of fermionic states |Ψ that satisfy This projects out the unphysical states |10 and |01 . As pointed out in [35], the unprojected system can also describe the degrees of freedom of a spin 1 particle Thinking reversely, we propose to identify our RBM state |Φ with the Pfaffian state by projecting onto physical states 2 where the positions z i of spin 1 operators ψ are identified with complex lattice point coordinates of lattice site i on the torus. Put another way, an electron at site i of the Honeycomb lattice is represented in the conformal block by ψ(z i ). Note that the number of lattice sites is twice the number of plaquettes and therefore always even such that the Pfaffian is well-defined. Moreover, note that the projection operator P physical projects our spin 1 singlet onto the following spin 1 2 tensor product The relation (63) should not be thought of as an exact equality but rather as an approximation. The reason is that due to the complexity of the Pfaffian state it is not possible for us, yet, to show that the projected state respects all the symmetries of the Honeycomb Hamiltonian. In fact, our Ansatz (56) is highly symmetric and is invariant under the 120 degree rotational symmetry of the Honeycomb lattice as well as the symmetry operator which permutes the Pauli matrices, as described in Appendix A. Together, these symmetries leave the Hamiltonian invariant and should thus also leave the groundstate sector invariant, which in Appendix A is used to conjecture a more precise version of the correspondence (63). Moreover, we note that our Ansatz satisfies three key criteria which the ground state wave-function has to respect. Namely, it is invariant under translations by torus lattice generators, it is triply degenerate, and it is lying in a 2 m -dimensional sub-space of the 2 2m -dimensional Hilbert space which also happens to be the dimension of the vortex-free sector.

VII. CONCLUSION
In this paper, we have constructed RBM realizations of several sectors of the Kitaev Honeycomb model. In particular, we have focused on the gap-less non-Abelian phase corresponding to the parameter choice J x = J y = J z and used machine learning algorithms to train restricted Boltzmann machines to find groundstates and excited states of such systems. We used two different approaches to find these states, namely NetKet and PyTorch. We find that using PyTorch and our own software package, we can achieve a very accurate result for the groundstate energy for small lattice sizes. On the other, the RBM trained with NetKet does not give an accurate enough result for the energy but the sampling algorithm employed is very efficient and allows NetKet to also proper larger lattices. An ideal approach would be to combine NetKet's sampling algorithm with our own groundstate search algorithm. We leave this task for the future. Since our RBM wave-functions have 99.8% overlap with the exact groundstate (this result is for the 3 × 3 lattice), we find that we can transition to vortex sectors by manipulating RBM parameters as described in Section IV. This transitioning into vortex sectors has been also demonstrated starting from a wave-function obtained from quantum state tomography. Here we are able to reconstruct an exact wave function from single-shot measurements and use this method to "detect" the vortex state of the system. This is useful for applications in topological quantum computation where one would need to "decode" input data, then perform a quantum gate operation, and finally create an output.
Of course, in the case of the Honeycomb model, we have an exact solution of the spectrum using Kitaev's Majorana fermion representation. However, such a representation, though giving exact values for the energies of given states in the Hilbert space, is at the same time losing information. For example, it is difficult and in some cases unknown how to compute expectation values of string operators composed of Pauli matrices. Moreover, one may lose the Majorana representation in cases where a complicated perturbation is added to the Hamiltonian which obscures such a simple representation. For all these reasons it is desirable to obtain realizations of the quantum state in terms of spin eigenstates and the RBM gave in Eq. (21) is such a representation.
One important application of our results which we present in Section VI, is to relate our RBM groundstate wave-function to the Moore-Read Pfaffian state of the Ising CFT. This relation is part of the broader framework of correspondence between 2+1 dimensional topological field theories and 2d CFTs which goes back to the work of Witten on WZW models and Chern-Simons theory [19,33]. Our relation is a proposal which passes first non-trivial tests but should be examined in more detail to find further evidence and/or possible proof for it. We leave this for future work. Here, we note that such a correspondence is of great significance as it opens a door for neural networks to represent quantum field theory correlation functions of a high degree of complexity. On the one hand, the parameters of the RBM wave-function grow only polynomially with the spin number, while the CFT correlation function to which it is identified is of exponential complexity in the spin number. We hope that in the future one can generalize this framework/recipe to other models. It would be also very interesting to see whether there is a connection between what we propose, namely using an RBM to compute CFT conformal blocks, and the proposal made in [40] regarding a novel relation between neural networks and quantum field theory.

ACKNOWLEDGMENTS
The work of MN is partially supported by Trinity-Henry Barlow Scholarship. MN would like to thank Yau Mathematical Science Center (YMSC) for hospitality during the early stages of this work. He also thanks Martin Duy Tat for helpful discussions. The work of BH is supported by the National Thousand-Young-Talents Program of China. BH would also like to thank the Bethe Center for Theoretical Physics in Bonn, where part of this work was completed, for hospitality. We would furthermore like to thank Jin Chen, Bartek Czech and Ce Shen for valuable discussions. Part of the numerical calculations were performed on the SARMAD computing cluster at Shahid Beheshti University but the majority of the work was done using the Lambda cluster at YMSC.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. where in the first line we have applied the SL(2, Z) transformation and in the last line we have made use of the identity −1 τ −1 = τ for our particular choice of τ corresponding to the honeycomb lattice. Similarly, one can derive the identity Combining equation (A6) and (A8), we then see that the ratio of the two transforms as follows under the action of thê C 6 operatorĈ 6 ϑ[α](τ, z i − z j ) ϑ 1 (τ, z i − z j ) = ϑ[α](τ, (τ − 1)(z i − z j )) ϑ 1 (τ, (τ − 1)(z i − z j )) = κ([α], γ) −1 κ 1/2 1/2 , γ ϑ[β](τ, z i − z j ) where the resulting spin structure β of the theta function after the transformation is given as follows Thus we see, that as the two factors of κ just contribute an overall phase which is the same for all summands of the Pfaffian in expression (56), the three groundstate wave functions just get permuted under the action ofĈ 6 . This is, however, a too large symmetry group as the groundstate space should be only invariant under the combined action ofĈ 6 andÛ C6 . To understand what has been happening and how to modify the Ansatz such that it is only invariant under the combined symmetry, it is useful to have another look at the projection (64). When performing such projections, there is a freedom for the choice of phase factors and it is a-priory not clear which ones to choose. In (64) we have made the simplest choice for such factors. Here, we want to modify our choice and present a more precise conjecture for the relation between the Pfaffian state and the groundstate wave-functions. To this end, note that the Honeycomb lattice can be divided into three conical regions as shown in Fig. 10. These regions are connected by 120 degree rotations and invariant under inversions. When applyingP physical , we introduce phases depending on in which region the complex value x i − z j lies: P physical |ij 1 = e 2πi/3 | ↓ i | ↓ j if z i − z j ∈ region III. (A13) Next, note the above three spin states are each eigenvectors ofÛ C6 with eigenvalues given as followŝ Now it is easy to see that a permutation which moves P physical |ij 1 from region I to region II is combined with an action ofÛ C6 with eigenvalue 1, a permutation which moves the state from region II to region III will be combined with an eigenvalue e 2πi/3 from the action ofÛ C6 resulting in a state in region III, and a permutation moving the state from region III to I will lead to a cancellation of phases upon the action ofÛ C6 resulting in a target state with the correct phase. Thus it is apparent that the Pfaffian state (56) is invariant (up to a permutation of theta function spin structures as described above) under the combined action ofĈ 6 andÛ C6 while it is not invariant under individual actions of these operators. This mirrors the symmetry properties of the RBM state under the action of these operators exactly. We therefore conjecture the Pfaffian state to be an exact groundstate under the refined action of P physical as described above.