Variational solutions to fermion-to-qubit mappings in two spatial dimensions

Through the introduction of auxiliary fermions, or an enlarged spin space, one can map local fermion Hamiltonians onto local spin Hamiltonians, at the expense of introducing a set of additional constraints. We present a variational Monte-Carlo framework to study fermionic systems through higher-dimensional (>1D) Jordan-Wigner transformations. We provide exact solutions to the parity and Gauss-law constraints that are encountered in bosonization procedures. We study the $t$-$V$ model in 2D and demonstrate how both the ground state and the low-energy excitation spectra can be retrieved in combination with neural network quantum state ansatze.


Introduction
Studying the mapping between fermionic operators and bosonic operators (and vice versa) is interesting both from a theory perspective, as well as for computational studies. As examples of the former, the 1D transverse-field Ising model and the Kitaev honeycomb model can be diagonalized after reformulating the Hamiltonian in terms of fermionic degrees of freedom [1]. Furthermore, transforming fermionic Hamiltonians into a set of spin operators is necessary to compute properties of fermionic systems using digital quantum devices. Especially in the NISQ era [2] of intermediate scale devices, it is highly advantageous to study efficient mappings that require fewer qubit resources.
The most natural and commonly used mapping between fermionic and spin degrees of freedom is the Jordan-Wigner transformation (JWT), which Jannes Nys: jannes.nys@epfl.ch follows as a natural consequence of the second quantization formalism of fermions [3]. After we have chosen a fermion ordering in the second quantization formalism, the JWT maps each fermion operator f † i onto spin operators as fol- are Pauli matrices applied to spin i. The operator chain S i = ⊗ j<i Z j is necessary to maintain the anti-commutation relations on the fermionic side using a set of Pauli matrix operators that themselves follow commutation relations, and are commonly referred to as Jordan-Wigner strings. Physical fermionic Hamiltonians that describe closed quantum systems consist of bilinear and quadratic terms in the creation/annihilation operators. Such Hamiltonians conserve the fermion parity P f = (−1) N f (with N f the number of fermions), and are referred to as even parity operators. When the original fermionic operators are both spatially local and have even parity (i.e. conserve P f ), the locality is trivially preserved in the resulting spin Hamiltonian in 1D, since the Jordan-Wigner strings S i of local fermionic operator pairs cancel each other. In higher dimensions (>1D), however, the chosen ordering of the fermions in the JWT becomes increasingly important. Local evenparity fermionic operators are no longer mapped onto a set of local products of Pauli matrices since the Jordan-Wigner strings S i of spatially local fermion operator pairs no longer cancel each other. When the dimensionality of the system increases, the Jordan-Wigner strings in the spin Hamiltonian become increasingly non-local and it therefore increasingly difficult to study these systems numerically [4].
The mapping of fermion operators onto quantum spin operators is not unique. One can use this freedom in order to generalizations to fermion-spin mappings in higher dimensions with the main aim to maintain locality in the operators, and thereby reducing the size of the Jordan-Wigner strings. One of the first studies that derived higher-dimensional generalizations to the Jordan-Wigner transformation dates back to the work of Wosiek [5]. In this work, Wosiek described a mapping of fermions moving on a 2D and 3D square lattice onto a set of local generalized Euclidean Dirac matrices. Thereby, he found the need to impose additional constraints on the system in order to remove redundant and unphysical sectors of the new Hilbert space. The constraints generated in this bosonization procedure were studied numerically only recently in Ref. [6]. Similar ideas were later explored by Bravyi and Kitaev [7] in the simulation of fermionic systems through local qubit gates. As in Ref. [8], they explored methods that increase the Hilbert space, while afterwards restricting the reachable quantum states to a physical sector of Hilbert space through a set of gauge conditions. Ball [9] demonstrated how these gauge conditions can be made local as well. Ball [9] and Verstraete-Cirac [8] both suggested the more explicit introduction of auxiliary fermionic modes to counteract the Jordan-Wigner strings. These auxiliary modes effectively store the parity nearby the interaction terms, which is otherwise captured by the Jordan-Wigner string [4]. The auxiliary Majorana fermions are subject to local interaction terms that commute with the physical Hamiltonian, which is necessary in order to keep the eigenspectrum of the original problem identifiable in the spectrum of the transformed Hamiltonian.
In recent years, we have witnessed a renewed interest in methods for simulating fermionic systems through a set of local qubit gates [10,11,12,13,14]. This more recent theoretical activity is again driven by two main motives. On one hand, local mappings are of practical interest in the implementation of quantum algorithms to simulate fermionic matter on increasingly available qubit-based digital quantum computers. On the other hand, a recent theoretical advance has been made in connecting bosonization in 2D and Z 2 lattice gauge theories with Chern-Simons-like Gauss laws [13] (later generalized to three and arbitrary spatial dimensions in Ref. [15,16]). Various followup works have built on this connection to derive new fermion mappings in higher dimensions [17,10,18,19]. Compared to earlier methods (such as Refs. [9,8]) where JWT were carried out explicitly, recent techniques [13,18,6] take a different approach where one first defines bosonic operators from the fermionic ones, which can then be mapped directly onto quantum spin operators without the need to order the fermions. The equivalence of these bosonization procedures in 2D was proven by Chen et al. [14].
Despite this recent theoretical progress, the practical application of these techniques remains elusive, both in the context of classical computational methods, and as a basis for quantum algorithms. The main difficulty lies in the fact that the auxiliary degrees of freedom must satisfy stringent constraints in order to correctly represent fermionic degrees of freedom. Efficiently satisfying these constraints is a particularly important task especially in applications involving variational searches of many-body fermionic states. This is for example relevant for both classical variational methods based on spin/qubit degrees of freedom and for variational quantum algorithms tailored to qubits. In this work, we specifically focus on the variational simulation of fermionic systems on classical computers, with suitable many-body quantum states of spins degrees of freedom. Specifically, we demonstrate that we can factorize the wave function into an exact solution to the constraint, and a physical wave function that can be determined variationally. Furthermore, maintaining spatial locality in the transformation can allow us to map the fermionic Hamiltonian onto a spin Hamiltonian which features the same symmetries [17,10], and therefore gives us access to the low-energy excitation spectrum. We then show that there is substantial variational freedom in parameterizing the resulting many-body state. In this context, we concentrate on on neural-network-based parameterizations of the many-body wave function, known as neural-network quantum states (NQS) [20]. We show how this approach can be used to approximate the energy eigenstates of the fermionic system, using methods that are available to approximate quantum spin states [20, 21].

Bosonization
We define an L x × L y 2D lattice with square cells and such that all edges point either along the lattice basis vectors x or y. The resulting set of edges reads E = {(r, r + x)|r ∈ V} ∪ {(r, r + y)|r ∈ V}. Here, x and y represent the lattice vectors and r = (r x , r y ), where r x ∈ {0, ..., L x −1} and r x = L x maps onto r x = 0 due to periodicity (similar for r y ). We study the t-V model (also called spinless Fermi-Hubbard model) where f † r are fermionic operators, and we have introduced the usual number operator n r = f † r f r . On each site, the physical Hilbert space is 2 dimensional {|0 , f † r |0 }, and the Hamiltonian in Eq. (1) contains only even parity fermionic operators. We carry out the bosonization procedure defined in Refs. [17,10], since it keeps the symmetries of the Hamiltonian manifest. For the sake of completeness, we describe the procedure in Appendix B, and summarize the results here. Bosonizing the Hamiltonian in Eq. (1) results in: On the resulting square lattice, each site hosts a physical (1) and auxiliary qubit/spin (2), as demonstrated in Fig. 1. As shown by Ref. [14], other bosonisation procedures are equivalent in 2D.

Local constraints
In order to reduce the number of degrees of freedom, the auxiliary system is subject to a Gausslaw constraint of the form In our notation, constraints such as the ones in Eqs. |Ψ . In terms of Pauli operators, G r in Eq. (3) takes the form We can separate the physical (1) and auxiliary (2) system and rewrite Eq.
The operator on the left-hand side of this constraint also appears in Wen's plaquette model [22], and hence, Eq. (5) can be interpreted as a dynamical version this model, due to the right-hand side which depends on the physical system. The constraint in Eq. (5) resembles a Chern-Simons Gauss law, or flux attachment [13]. The constraint is diagonal in the physical system, and therefore, for each configuration of the physical system, the auxiliary system is in an eigenstate of the exactly solvable Wen's plaquette model (which is known to have robust topologically degenerate ground states), with different signs for the terms in the Hamiltonian. It is important to emphasize that the constraint in Eq. (5) is 'kinematic', meaning it only depends on the lattice topology, not on the Hamiltonian under consideration.

Parity constraints
We restrict the system to square even-by-even tori (i.e. L = L x = L y is even) with N = L 2 sites. Imposing the boundary conditions in the fermionic system, we obtain the additional constraints introduced by non-contractable Wilson loops. After bosonization, we obtain the following spin operator identities that need to be satisfied Furthermore, we fix the number of fermions to be N f , which is enforced through the constraint We can now set the chemical potential µ = 0, since it only adds a constant energy.

Variational Monte Carlo approach
The goal of this paper is to obtain variational solutions to the bosonized fermionic system. Hereby, we rely on Variational Monte Carlo (VMC). We briefly recap the concepts of VMC within our notation. For a more elaborate and pedagogical introduction, we refer to e.g. Ref.
[23]. After a short recap, we will introduce our novel approaches to solving the gauge constraints described above. We will study the system in Eq. (1) using the bosonized Hamiltonian in Eq. (2). Our aim is to obtain the ground and low-lying excited states |Ψ using the decomposition In this notation |σ represent basis states in the S z basis. The systems consists of a physical sys- r N ) and auxiliary system We will take σ r ∈ {−1, +1} for simplicity. The probability amplitudes in Eq. (10) are given by a parametrized function Ψ θ with parameters θ. These parameters are determined by minimizing the variational energy We will evaluate the expectation value in Eq. (11) by sampling spin configurations using Markov-Chain Monte Carlo (MCMC). However, the constraints in Eqs. (3), and (7)-(9) restrict the allowed Hilbert space, and hence we must enforce these restrictions either by imposing constraints on Ψ θ , and/or by designing Markov-Chain samplers that only generate samples within the allowed Hilbert space.

Solving constraints in VMC
We now describe our novel approach to satisfy all constraints. To the best of our knowledge, there have so far been no attempts to solve local fermion-to-qubit transformed Hamiltonians in the VMC framework. The constraints in Eq. (7), (8) and (9) can be fulfilled exactly through a suitable sampling procedure since the constraints are diagonal in the S z basis. More specifically, MCMC samples can be generated through a set of sample updates based on the (free-fermion) Hamiltonian. Hence, given a sample, we can make the following Markov transitions r+x ): flip two neighboring qubits of system (1) along the x-axis: X r+y ): flip four qubits (2 physical, 2 auxiliary) on two neighboring sites along the y-axis: X It is also straightforward to generate initial random samples that fulfill the parity and number constraints, which are necessary to initiate the Markov chains. First, a physical system (1) fulfilling the N f constraint in Eq. (9) can be constructed. Next, we generate (L x − 1) × (L y − 1) random auxiliary qubit states, and infer the remaining auxiliary qubits by imposing the periodicity constraints in Eqs. (7) and (8). The resulting scheme is summarized in Algorithm 1. Alternatively, these parity constraints can be captured by a Restricted Boltzmann Machine (RBM) quantum state by introducing a hidden neuron with specific weight function [24].
The Gauss law in Eq. (3) cannot be satisfied through a sampling procedure alone, since the operators are not diagonal in the computational basis, and they therefore impose stringent constraints on the wave function itself.
Within the basis that fulfills the constraints in Eq. (7) and (8) one can obtain the eigenspectrum of the original Hamiltonian through eigendecomposition of P −1 G HP G , where P G is the projection operator to the Gaussian constraint in Eq. (3) While the representation of the P G operator in the S z basis is block diagonal, it is not sparse, since the size of each block scales exponentially with the system size. Furthermore, since constraint Eq.
(3) must be satisfied for all r ∈ V, the projection generates non-local effects in the auxiliary system, even though the individual constraints are local. Therefore, we must find an analytical form for the probability amplitude of a general quantum state that lies within the manifold of the Gauss law. Another approach is to implement the constraint in the Hamiltonian by adding the terms When the coupling K is taken to be sufficiently large, G r = 1 can in principle be satisfied by minimizing the total energy of the augmented Hamiltonian H = H + H c (this procedure is also suggested in Ref. [8]). In practice, however, we find that a soft constraint does not result in quantum states lying on the manifold dictated by the Gauss law, and therefore the spectrum does not reliably represent the physics of the fermionic system.
As mentioned, exactly respecting Eq. (3) is essential in order to restrict our solution to the physical Hilbert space representing fermionic degrees of freedom. Since the r.h.s. operator in Eq. (5) is diagonal in the physical states, each Gauss constraint can be regarded as a dynamical constraint on the auxiliary system. It is important to obtain variational ansatzes which obey the Gauss law constraints by construction. Our (non-symmetrized) variational ansatz in general assumes a factorized form where θ represent a set of variational parameters, and ξ(σ) is purely a sign factor with |ξ| = 1. In this representation, we choose ξ(σ) to be the only factor with an explicit dependence on the auxiliary system. Notice that there are no constraints on Φ θ with respect to anti-symmetry, since this property is entirely covered by the ξ parity factor. Hence, ξ must be chosen in such a way that Eq. (3) is exactly fulfilled. However, we point out that many solutions are feasible due to the freedom of extracting additional sign structure from Φ θ , and absorbing it into ξ. Indeed, since G r is diagonal in the physical system (1), we obtain independent Gauss-law constraints for each physical configuration σ (1) . After defining the sign-generating function ξ, the sign structure of Φ(σ (1) ) is determined by the local spin Hamiltonian in Eq. (2). Furthermore, the main challenge is to find an exact solution that satisfies all constraints, including the parity constraints in Eqs. (7)-(8).

Canonical sample reduction
As a first solution to the Gauss law constraint, we define a consistent approach which uses a repeated application of G r to a quantum state to transform each auxiliary configuration σ (2) to a reduced auxiliary sample α (2) . One can verify that the following ansatz is an eigenstate of the constraint operators G r : Here, m i ∈ {0, 1} is determined by the iterative procedure followed to obtain the reduced sample. Therefore, we define an ordering for the lattice sites (r 1 , ..., r N ) and iteratively set m i = 1 if the auxiliary qubit at position r i is in the |1 state, and m i = 0 otherwise. When m i = 1, we apply operator C r i to the auxiliary plaquette attached to r i and continue with the next site in the sequence r i+1 . Algorithm 2 describes the sequential steps to obtain either directly ξ(σ), or (m 1 , ..., m N ) to evaluate Eq. (14). Also note that to each physical system configuration σ (1) , there corresponds only a single α (2) . The time complexity of this approach to satisfy the Gauss constraint is O (N ) for each sample σ.

Doubly canonical sample reduction
In the abovementioned solution, the symmetry properties of Φ(σ (1) ) are elusive. Alternatively, we can choose to optimize the sign structure of Φ by incorporating knowledge from the relevant symmetry group. Therefore, the abovementioned procedure can also be carried out on the reduced samples obtained by listing the samples obtained by applying all symmetry elements g(σ) in the symmetry group g ∈ G. We then take the sample σ red with the smallest lexicographic encoding of the physical system and apply the above-mentioned sequence of reductions on σ red to obtain α (2) . The representation of symmetry operations will be discussed in more detail in Section 4.

Vacuum reduction
A final method approaches the problem differently by defining the correct sign structure ξ only explicitly on the vacuum state |0 (1) σ (2) , and inferring the signs of other configurations by consistently relating them to the vacuum. The latter is defined in terms of spin states as where c σ (2) = ξ(0 (1) , σ (2) ) = ±1 is obtained via the approach in Section 3.2, and thus satisfies Eq. (3). In the vacuum, the constraint in Eq. (4) reduces to C r = 1.
In order to obtain ξ(σ) from the definition of |Ω q , we take the set of occupation numbers (n 1 , ..., n N ) corresponding to σ (1) using n i = (1 − σ (1) i )/2), and rewrite these in terms of fermion operators in the usual way |n 1 , ..., can be transformed into a set of spin operator through the same bosonisation procedures that led to the spin Hamiltonian in Eq. (2). Our ansatz is now inspired by the above-mentioned mapping, and hence we assume where B represents the bosonisation procedure from Ref. [10], using the mappings in Eq. (31). Notice that to consistently bosonize non-local fermionic operators f † r N f ...f † r 1 in Eq. (17), this procedure requires an ordering of the sites, similarly to JWT. The above-mentioned solution also forms a solution to tackle the standard JWT Hamiltonian. The main conceptual difference that in the current formalism, the ansatz is varied to optimize a local Hamiltonian with manifest symmetries. However, as pointed out in Ref.
[18], despite the fact that the constraints in Eq. (3) are local, they can introduce long-range correlations when they are satisfied for all plaquettes (which shows up in the sign structure). It is important to point out that the variational part of our ansatz Φ θ (σ (1) ) is a function of the σ (1) , which is equivalent to an occupation configuration.
The resulting ansatz exactly fulfills the Gauss law in Eq. (3), and hence, the physical part of the wave function Φ(σ (1) ) is constraint free, meaning it does not require anti-symmetrization and therefore our method is determinant free. Hence, we may use a universal function approximator such as a Neural Network, to represent Φ θ (σ (1) ). In this work, we adopt a simple Restricted Boltzmann Machine (RBM) ansatz with complex weights and N hidden spins (thus with a hidden-spin density of α = 1) [20].

Quantum state symmetries
It is expected that imposing symmetries will result in more reliable and accurate predictions of the ground state. We therefore turn to the representation of symmetry transformations within the bosonisation framework. A representation T g of an element g of a symmetry group G, can be decomposed into two components: a 'bare' transformation T b , and an auxiliary-mode transformation V T , such that U T = V T T b . The V T are tensor products of local single-site operators: V T = r∈V V T,r . These effectively replace the (non-local) parity factor one would encounter in the Jordan-Wigner transformation. Once we are able to determine V T,r for all elements in the C 4v group, we can obtain a symmetric quantum state ansatz (lying withing a chosen irrep I of the group), using where χ g represents the character of the chosen irrep. For translations, we have More details on symmetry operations and the V T,r for translation, rotation and reflection symmetry are deferred to Appendix C. Notice that the factor ξ(g(σ)) inside Ψ is the alternative for a factor which determines the parity of the T g operator on σ. Hence, no sorting is required, and the parity is determined through local operators, since parity information is captured by the auxiliary modes. The above-mentioned methods can therefore be carried out in O (N ) time, contrary to other determinant-free methods Ref. [25]. The latter would also be necessary when carrying out 1D Jordan-Wigner procedures, in which one must compute the parity of the translation operator in terms of fermion operators.

Results
We carry out numerical studies of the bosonization procedures by relying on the abovementioned solutions to all constraints. Hereby, we first investigate their ability to represent the ground state and the effect of different choices of ξ(σ) on the complexity of learning Φ θ (σ (1) ).
Through these studies, we mainly probe the sensitivity of numerical and variational approaches to the dynamical character of the Gauss law constraint in Eq. (3), which was interpreted as a dynamical Wen plaquette model. Consider first an unsymmetrized variational ansatz. In Fig. 2, we show how the different methods to satisfy the Gauss law perform for a 2 × L lattices. We compare the results to the ones obtained through a Jordan-Wigner transformation mapping the two dimensional lattice onto a one-dimensional one through "snaking" along the L direction. By increasing L, we increase the degree of non-locality in the Jordan-Wignertransformed Hamiltonian, where Jordan-Wigner strings appear explicitly in the kinetic terms. Indeed, hopping operators transform under JWT as f † r f r+y → Q − r r∈P i→j Z r Q + r+y , where P i→j represents all sites on the path along the snake direction connecting site r and r + y. In the extreme case of r = (0, 0), the length of this path is 2L − 2. Although the canonical and doubly canonical methods are valid solutions to the Gauss law constraints, the results demonstrate that the corresponding physical wave function factor Φ θ (σ (1) ) are challenging to optimize. The underlying reason is that the ξ factors impose a non-trivial sign structure on Φ θ . Although the doubly canonical method should simplify the sign structure to be learned for samples that are connected though the considered symmetry operations, the ground state does not necessarily correspond to trivial irrep characters. Furthermore, the variational wave function Φ θ must still capture the signs of samples that are not connected through symmetry operators, an effect that is reflected in the inaccurate ground-state representations with this method. The vacuum reduction, however, only uses the non-local re- (0, π 2 ) (0,π) ( π 2 , π 2 ) ( π 2 ,π) (π,π) 10 −2 10 −1 10 0 10 1 (0, π 2 ) (0,π) ( π 2 , π 2 ) ( π 2 ,π) (π,π) duction method in Section 3.2 on the vacuum state (by solving the corresponding Wen plaquette model [22]). This results in energies closely matching the Exact Diagonalization (ED) results, and closely match those obtained with the snaked Jordan-Wigner approach. Similar conclusions can be drawn from Fig. 2b, where we increased the short dimension of the lattice, as well as the interaction strength. A difference between the canonical and vacuum method becomes increasingly visible for the 4×L lattice at large L, where the vacuum reduction method tends to result in lower energies. Notice that the canonical and vacuum methods result in similar energies for 2 × L lattices since the canonicalization procedure operates on a single chain of plaquettes. On the other hand, for 4 × L lattices, we follow a snaking procedure to canonicalize the plaquettes.
Next, we consider a 4 × 4 system, and investigate the effect of embedding symmetries in the variational ansatz. We focus on the vacuum reduction method from Section 3.4, which generated superior results in the previous experiment. We compare the results to the ones obtainable through ED in Fig. 3. Note that the Hamiltonian is diagonalized on the manifold determined by the total projector P G P x P y P N f , where the projectors P x , P y , P N f enforce the parity constraints in Eqs. (7), (8) and (9) respectively. The unsymmetrized ansatz generates accurate groundstate energies, except for N f = 6. In the latter case, the ground state has a challenging sign structure, since the corresponding translationsymmetry sectors are ( π 2 , π 2 ) and (0, π). On the other hand, N f = 4 has an isolated ground state, and the ground-state for N f = 8 contains the simpler (0, 0) sector. Therefore, the ground states of N f = 4, 8 are more accessible during the VMC optimization. Hence, while the variational ansatz can accurately represent the ground state, and is only hindered by the optimization procedure. Symmetrizing the ansatz removes the local minimum encountered in the optimization, and therefore returns significantly better ground state estimates. When the energy gap is small, the unsymmetrized ansatz does not accurately represent the ground state. The latter is not necessarily due to a limited representational power, but is directly related to the optimization procedure via Stochastic Reconfiguration (SR), which heavily depends on the size of the gap. Symmetrizing the ansatz and restricting to a given irrep opens the gap and therefore generates more accurate results.

Conclusions
We showed that bosonization procedures which allow one to transform local fermionic operators into local qubit operators result in Gauss law constraints that must be fulfilled in order to reduce the Hilbert space. Fulfilling all constraints at the same time is challenging, and we introduced multiple approaches that solve these constraints ex-actly, without the need to apply projection operators. We applied our method to the t-V model on a square lattice in 2D. We argued that constraints related to parity and boundary conditions can be fulfilled through the sampling procedure. Our experiments demonstrate that satisfying the Gauss constraints directly affects our ability to reliably represent the ground state of a given Hamiltonian. We found that imposing a sign structure on the vacuum state simplifies this challenge. However, the gauge constraint introduces non-local effects, as would also be expected from a standard Jordan-Wigner transformation. By following the bosonization procedure from Ref. [10], we keep the symmetries of the fermionic system manifest, which allows us to restrict the variational ansatz to a chosen symmetry irrep using techniques that are commonly used to study quantum spin systems. We demonstrate that this allows us to study the low-energy spectrum using neuralnetwork quantum state ansatze. Since the Gauss constraint depends only on the lattice topology, our approach can in the future directly be applied to other Hamiltonians on square lattices. We foresee extensions of our work to other lattices and higher dimensions, as well as studies of the bosonic systems that are equivalent to fermionic systems. and auxiliary majorana fermions". Phys. Rev. B 106, 115109 (2022).
[20] Giuseppe Carleo and Matthias Troyer. "Solving the quantum many-body problem with artificial neural networks".   Figure 4: Conventions for the ordering of the auxiliary χ i r modes, and the directions of the edges that determine the sign of the Λ iµ r Λ jν r terms. Blue boxes correspond to a single location vector r.

B.2 Parton decomposition
In a last step, we again move to fermionic (parton) operators. For the physical modes, we have and for the auxiliary modes Notice that our parton Hilbert space is now 2 3 dimensional at each site. The on-site Hilbert space reads where n i=a,b,c ∈ {0, 1} and |Ω represents the vacuum.

B.3 On-site Jordan-Wigner transformation
The bosonic operators defined in Eq.
where we identified |0 = |↑ and |1 = |↓ and As we will show in the next section, using the abovementioned mapping, one can rewrite the Hamiltonian in Eq. (27) in terms of Pauli operators using the definitions in Eq. (22). However, we first demonstrate how to eliminate the effect of one of the auxiliary qubits in Eq. (31).

B.4 Constraints
The Hilbert space is enlarged due to the b † and c † auxiliary fermion modes. The on-site parity operator Γ r for any site commutes with the Hamiltonian in Eq. (2). To reduce the Hilbert space, we restrict the on-site parton parity operator (which commutes with the Hamiltonian). Hence, the third spin is redundant and we can remove it by absorbing its effect in the other spins Using Eq. (35) to eliminate the third qubit, we obtain the expression for the Hamiltonian in Eq. (2).

(36)
While the above reduces to the identity operator in the fermionic formalism, it introduces a (gauge) constraint on the bosonic side. We can rewrite the above in terms of bosonic operators using Eq. (24)-(25), and ultimately in terms of the gauge operators Υ using Eq. (22). After some algebra, we obtain the constraint Υ 24 r Υ 32 r+x Υ 13 r+x+y Υ 41 r+y c = −1 (37) Hence, the auxiliary system is subject to a Gauss-law constraint of the form Imposing the boundary conditions in the fermionic system, we obtain the additional constraints introduced by non-contractable Wilson loops W x,y : W x,y = (−1) Lx,y . To see this, we carry out the same procedure as in Eq. (37) on these loops (for even-by-even tori). After bosonisation, we obtain the following spin operator identities that need to be satisfied