Determining quantum phase diagrams of topological Kitaev-inspired models on NISQ quantum hardware

Topological protection is employed in fault-tolerant error correction and in developing quantum algorithms with topological qubits. But, topological protection intrinsic to models being simulated, also robustly protects calculations, even on NISQ hardware. We leverage it by simulating Kitaev-inspired models on IBM quantum computers and accurately determining their phase diagrams. This requires constructing conventional quantum circuits for Majorana braiding to prepare the ground states of Kitaev-inspired models. The entanglement entropy is then measured to calculate the quantum phase boundaries. We show how maintaining particle-hole symmetry when sampling through the Brillouin zone is critical to obtaining high accuracy. This work illustrates how topological protection intrinsic to a quantum model can be employed to perform robust calculations on NISQ hardware, when one measures the appropriate protected quantum properties. It opens the door for further simulation of topological quantum models on quantum hardware available today.

Topological protection is employed in faulttolerant error correction and in developing quantum algorithms with topological qubits. But, topological protection intrinsic to models being simulated, also robustly protects calculations, even on NISQ hardware. We illustrate how this works by simulating Kitaev-inspired models on IBM quantum computers and accurately determining their phase diagrams. This requires constructing conventional quantum circuits for Majorana braiding to prepare the ground states of Kitaev-inspired models. The entanglement entropy is then measured to calculate the quantum phase boundaries. We show how maintaining particle-hole symmetry when sampling through the Brillouin zone is critical to obtaining high accuracy. This work illustrates how topological protection intrinsic to a quantum model can be employed to perform robust calculations on NISQ hardware, when one measures the appropriate protected quantum properties. It opens the door for further simulation of topological quantum models on quantum hardware available today.

introduction
Quantum computers are fragile and susceptible to rapid decoherence [1], and the use of topological protection has been proposed as a potential remedy. Typically, there are two propositions to make quantum computation practical for deep circuits: (i) faulttolerant quantum computers, which use logical qubits and error correcting codes (like the surface code) to correct all errors that appear during computation [2][3][4][5][6][7][8][9][10][11] and (ii) topological quantum computers, which use topological qubits, naturally protected from the environment, to perform intrinsically error-free computations [12][13][14]. Here, we propose a third avenue restricted to the simulation of specific models: when the model under study has non-trivial topological protections, the model's topological properties (rather than A. F. Kemper: akemper@ncsu.edu the hardware's) can be leveraged to identify its sharp phase transitions as the parameters are tuned.
A paradigmatic example of a ground state with non-trivial topological properties, occurs in the Kitaev spin model on the honeycomb lattice [15], which has spawned a class of related models with similar properties . Their solution involves the braiding and fusing of anyons (in this case Majorana fermions), just as algorithms employing topological qubits do. Ideally, topological quantum computers could be used to solve models of this type. Unfortunately, it has proved to be extremely difficult to create topological qubits in quantum hardware [42].
Instead, we illustrate that these Kitaev-inspired models can be solved on conventional quantum computers with the braiding operations encoded in lowdepth conventional quantum gates. Moreover, we can utilize the distinct topological properties of ground states to clearly distinguish the phase transitions in this class of models for system sizes where noise makes this difficult if one uses a local quantity instead. These topological quantum simulations are carried out on IBM's transmon-based "conventional" quantum computing hardware.
Typically, determining the phase boundaries in these models is challenging because they lack a local order parameter. The standard approach is to calculate the entanglement entropy [43,44], but this requires a large system size, beyond what is currently available in NISQ machines. Instead, we develop and use particle-hole symmetry-enforced circuits encoding twisted boundary conditions, enabling determination of ground-state properties through the entire Brillouin zone with only a modest number of qubits. Due to the discrete topological nature of the quantum phase transition, these circuits can identify the different phases by revealing discontinuities in the observables, even if the jump at the discontinuity is diminished due to noise.
The paper is organized as follows. We begin by discussing the general procedure for diagonalizing Kitaev-inspired models in Sec. 2.1. This is followed by a step-by-step translation of the procedure to a sequence of quantum circuits in Sec. 2.2, with a particular focus on the braiding operations. We finish the section with a discussion on the generalization and scaling of the quantum circuits in Sec. 2.3. With the circuits from Sec. 2.2 in hand, we apply these to the 1D Kitaev spin chain and the Kitaev honeycomb model in Sec. 3 and measure the entanglement entropy. In Sec. 4, we propose symmetry-enforced circuits to sample the Brillouin zone while maintaining a fixed number of qubits, and demonstrate the application of these circuits by measuring the entanglement spectra of 1D Kitaev spin chain via correlation functions. This enables the determination of the phase diagrams for the 1D Kitaev spin chain and Kitaev honeycomb models, as shown in Sec. 5. We conclude in Sec. 6.
2 Kitaev-inspired models and the quantum circuit to prepare their ground states 2

.1 Kitaev-inspired models and exact solvability
Based on the seminal work by Kitaev [15], significant effort was made to generalize the exact solution procedures to other strongly correlated spin models . These models can be solved exactly by using the same strategy proposed by Kitaev [15] and as such we refer to them as "Kitaev-inspired" models in this work.
The exact solution of these models begins with assigning a proper Jordan-Wigner (JW) string ordering for all the sites, along which the JW transformation is applied. After the JW transformation, the coupling terms involving transverse (x and y) components of spins can be transformed to be quadratic in terms of Majorana fermions (denoted here by γ operators).
The coupling terms between the z-component of spins are transformed into density-density interactions between the two sites. In terms of Majorana fermions, the occupation number at a generic site j is proportional to γ j η j , so that the z-component coupling term between sites i and j is ∝ γ i η i γ j η j . Here, the fermion at site i is composed of Majorana fermions η i and γ i . According to Kitaev [15,43], and in agreement with subsequent numerical verification [39] for the ground state in the thermodynamic limit, the product of the two η Majorana fermions reduces to a site-independent gauge field due to local conservations [43,45]. Consequently, as long as we are interested in ground states of the model, the z-component coupling terms can be treated as if they are only quadratic (in terms of γ Majorana fermions) and the quadratic in the η fermions can be replaced by the constant gauge field value.
To complete the exact solution, new fermion operators are defined from the γ Majorana fermions, and the resulting models can be further diagonalized by Fourier transformations followed with Bogoliubov transformations. In the following sections, we will denote the composite fermion constructed from two η Majorana fermions as the gauge fermion, and the composite fermion constructed from the two γ Majorana fermions as the bond fermion.
To explicitly demonstrate how to study these models on quantum hardware, we focus on the original Kitaev honeycomb model and its 1D simplification, the Kitaev spin chain (see Appendix A for more details), in the main text. Because of the generality of the exact solution procedures for the Kitaev-inspired models, the methods presented here can be readily applied to other Kitaev-inspired models. These include, for example, the generalizations of the original Kitaev honeycomb model to different lattices [16-19, 27, 34], different dimensions [20,23,24,[38][39][40], and higher spins [26,29]. In the Appendix D, one additional example of applying our approach to the 1D BCS-Hubbard model is provided.

Quantum circuit constructions
To illustrate how to construct the quantum circuit that prepares the ground states of Kitaev-inspired models, we focus on the original Kitaev model constructed on the honeycomb lattice with 8 sites (resulting in 2 × 2 unit cells) as shown in Fig. 1(a). To prepare the ground-state wave function in real space, we begin in momentum space (after all required transformations), where the Hamiltonian is diagonalized (and then go backwards through all of the steps required for solving the problem to end up back in real space). The ground state (in momentum space) is then written as a simple product state. In particular, for the model under consideration, the ground state is just the vacuum state |0000 , with the qubits representing the four momentum points of the elementary excitations with (from left to right) (k x , k y ) = (−π/2, −π/2), (k x , k y ) = (π/2, π/2), (k x , k y ) = (π/2, −π/2) and (k x , k y ) = (−π/2, π/2), respectively (We work with periodic boundary conditions). In the quantum circuit shown in Fig. 1(b), these momentum points correspond to the qubits indicated by the orange dots on the far left, ordered from bottom to top. The circuit construction for the Bogoliubov and fermionic Fourier transform (shown in Fig. 1(c)) are detailed in Refs. [46,47], which we reproduce in the Appendix B for completeness.
Here we divide the 8 sites into four unit cells as shown by the blue shaded regions in Fig. 1(a) so that each unit cell contains one gauge fermion and one bond fermion. In particular, we define the unit cell 1 containing sites 1 and 3, the unit cell 2 containing sites 2 and 5, the unit cell 3 containing sites 4 and 7, and the unit cell 4 containing sites 6 and 8. Then the 4 qubits after the Fourier transformations from bottom to top in Fig. 1(b) correspond to the bond fermions in the unit cells from 1 to 4. The gauge fermions in the unit cells from 1 to 4 are represented by the 4 red dots from bottom to top in Fig. 1  Majorana braiding inter-fermion braiding intra-fermion braiding Corresponding to the exact solution procedure that simplifies the z-component coupling terms, only the braiding of the Majorana fermions needs to be performed in the quantum circuit. To facilitate the implementation of the braiding operations, the bond fermion and the gauge fermion belonging to the same unit cell must be neighbors. In the quantum circuit, this is achieved by a series of fermion swap operations after the Fourier transformations as shown in Fig. 1(b). Then, the braiding operations are implemented on the 4 sets of neighboring qubits as shown by the purple boxes in Fig. 1(b). The braiding circuit is discussed in detail below. The final step is to rearrange the fermions according to the JW string ordering, and this is done by the 3 swap operations shown in Fig. 1(b).

Quantum Circuits for Braiding Operations
The strategy we use to create quantum circuits is to recognize that the braiding operations can be thought of as an intermediary that maps the conventional fermions to themselves-after identifying this mapping, we can then construct circuits within the conventional fermion language. We describe how this is done next.
The simplest braiding is the braiding of two Majorana fermions to construct a conventional fermion.
A conventional fermion, which can be written as f = (γ 2 +iγ 1 )/2, is connected to another conventional fermion defined byf = (γ 1 + iγ 2 )/2 via a braiding operation on the Majorana fermions. Recall that the braiding operator (in terms of Majorana fermions) is given by [45,49] In this first example, the two Majorana fermions come from the same conventional fermion. After expressing the Majorana fermions in terms of the conventional fermions, we find that where f f and f † f † vanish due to the Pauli principle.
To express this operator as a matrix in the qubit representation, we map the computational basis (|0 , |1 ) onto the conventional fermion number basis yielding In our second example, we consider the case of the braiding of two Majorana fermions that belong to two different conventional fermions. The conventional fermions are when expressed in terms of the different Majorana fermions.
There are four different ways to braid these Majorana fermions. Because braiding operations do not commute, we have to carefully specify the ordering scheme for each set of fermions. We choose an ascending order for the conventional fermions and also for the Majorana fermions within the same conventional fermion. With this specification, we find that the different inter-fermion braidings can be expressed as the product of the braidings within the same fermions and one more braiding, denoted by B ± 2,3 . This is the fundamental braiding between two conventional fermions and is denoted B ex in the later discussions. Replacing the Majorana operators by fermion operators, we find that Again, we using a conventional number basis , yields a matrix representation for the inter-fermion braiding Next, we decompose this braiding operation in terms of conventional quantum gates. Since the clockwise braiding is the Hermitian conjugate of the anticlockwise braiding, we focus on the clockwise braiding. We first consider the braiding within the same fermion, and we notice that so B in can be realized by the phase gate U 1 (π/2) as shown in the left panel in Fig. 1(c). Similarly, we find: So, the quantum circuit can be immediately constructed via the general scheme proposed by Vidal et al. [50]. The detailed circuit for B + ex is shown in the left panel in Fig. 1(c).
For the Kitaev-inspired models, the braiding operations are used to change the positions of Majorana fermions, allowing Majorana fermions belonging to two different conventional fermions in one unit cell to recombine into new fermions. This implies that all braiding processes can be constructed with one inter-fermion braiding and a few intra-fermion braidings. The inter-fermion braiding changes the original two-qubit states into a superposition of their particle-hole partners with a ±π/2 phase difference. Intra-fermion braidings behave differently depending on whether they are performed before or after the inter-fermion braiding. If the intra-fermion braiding is performed before the inter-fermion braiding, it introduces a ±π/4 phase to the fermion. However, if the intra-fermion braiding is performed after the interfermion braiding, it introduces a superposition of the original state (before inter-fermion braiding) and its particle-hole partner (with a ±π/2 phase difference). These basic properties of braiding help us construct braiding circuits based on local constraints.
The first constraint imposed on the braiding operations comes from the local z-bonds. For Kitaevinspired models, the local fermionic Hamiltonian coupling between two sites by a z-bond is where c 1 and c 2 are the annihilation operators for the fermions on the site 1 and site 2. The two-qubit states are denoted by |n 1 n 2 . By the braiding operation U B , this Hamiltonian can be re-expressed in terms of bond and gauge fermions as where D is the local gauge field and f stands for the bond fermion. Correspondingly, after the braiding transformation, the two-qubit states here are |n f n g with n g (n f ) denoting the occupation of the gauge (bond) fermion. This implies that the states |n f n g are transformed to the states |n 1 n 2 as follows: where α n f ng n1n2 denotes the transformation coefficient. The two Hamiltonians and their corresponding states are connected by braiding operations, so the matrix elements of the Hamiltonian (before and after transformation) needs to be consistent. A direct calculation yields and Therefore, the consistency imposes the constraint D = 1 for n g = 1 and D = −1 for n g = 0. This local constraint imposed by z-bonds is general for Kitaev-inspired models. But, the local constraints imposed by other coupling terms (i.e. x-, and y-bonds for the models considered in this work) depend on the JW strings, because those strings determine the ordering of bond fermions emerging from different braiding operations. By using the same analysis as shown here for the z-bond, the details of the braiding operations can be further determined from the constraints set by other coupling terms. For the Kitaev honeycomb model on the honeycomb lattice with the JW string in the ordering given in Fig. 1(a), the braiding operation can be implemented by the circuit shown in the left panel in Fig. 1(c). The details of constructing the braiding circuit for the 1D Kitaev spin chain can be found in Appendix C.

Generalization and scaling
We have illustrated how to construct the quantum circuit to prepare the ground state of the Kitaev honeycomb model with 2 × 2 unit cells in real space by reversing the methodology for the exact solution of the model. Because the exact solution procedures can be generally applied to a large class of Kitaev inspired models , the circuit construction method illustrated here can be readily extended for other Kitaevinspired models. The quantum circuits represented by the orange and green boxes in Fig. 1 are elementary components for the implementation of generic Bogoliubov and Fourier transformations. Though the braiding circuits shown in Fig. 1 would change according to different definitions of JW strings for the different models, the methodology to construct circuit will be the same. For example, in Appendix D we show that the circuit construction method can be applied to the 1D BCS-Hubbard model [39], which is another example of a Kitaev-inspired model.
Before we close this section, we would like to discuss the scaling of the circuit for the Kitaev honeycomb model. Because of the Fourier transformation, the system size must contain N = 2 n × 2 m unit cells with 2 n along the x direction and 2 m along the y direction. To prepare the ground state of such a system, the required number of qubits is 2N . The Bogliubov transformation is implemented on opposite momentum point pairs, so we require N /2 twoqubit Bogoliubov transformations (the corresponding circuit shown in Fig. 1(c)) in building the full circuit. We know that for N momentum points along one direction, the Fourier transformation would require N log N two-qubit Fourier transformations (the circuit as shown in Fig. 1(c)) [46][47][48]. To complete the Fourier transformation along both the x and y directions, the total number of two-qubit Fourier transformations is still N log N . The braiding operations are performed within each unit cell, so N two-qubit braiding operations shown in Fig. 1(c) are required. Finally, we count the number of swap boxes (shown in Fig. 1(c)) in the quantum circuit, and they would appear in the following four places according to the generic quantum circuit shown in Fig. 1(b): (i) after the Bogoliubov transformation to facilitate the Fourier transformation along the y direction; (ii) after the transformation in the y direction to facilitate the Fourier transformation along the x direction; (iii) when pairing the bond fermion with the gauge fermion in the same unit cell; and (iv) when reordering the fermions to map to the original spin representation. The number of swap boxes is then In the implementation of quantum circuit the swap boxes in (iv) can be saved by relabeling the qubits and counting the sign changes, which can be done classically. In summary, in the quantum circuit to prepare the ground state of the Kitaev honeycomb model, the number of two-qubit Bogoliubov transformation circuits and the number of two-qubit braiding circuits are each proportional to the lattice size as ∝ O(N ). The number of two-qubit Fourier transformation circuits is O(N log N ), and the number of swap boxes is proportional to the square of the lattice size as ∝ O(N 2 ). A similar analysis can be done for other Kitaev-inspired models as well, but we do not do so here.
3 Measuring the ground state properties of Kitaev-inspired models on real quantum hardware By using the procedures described in the section 2, we can construct the quantum circuits to prepare the ground states for the 1D Kitaev spin chain containing two unit cells (including 4 sites overall) and the Kitaev honeycomb model containing 2 × 2 unit cells (including 8 sites overall). The lattice configuration of the two models are shown in the the insets in Fig. 2(a) and (b) respectively. Indeed, these circuits can correctly prepare the ground states of the two models as shown in Appendix E. Given that the ground states of the models are topological, a useful quantity to characterize them is the entanglement entropy. To measure the entanglement entropy, the lattice of the model is divided into two subsystems; the entanglement entropy of the subsystem A (the green colored regions in the insets of Fig. 2) is found by projectively measuring the density matrix with the help of the maximum likelihood estimator [51]. While this approach clearly does not scale, it is feasible for these proof-of-principle size systems.
In the upper panel in Fig. 2(a) we plot the entanglement entropy for the 1D Kitaev spin chain as a function of the coupling strength J x by setting the other coupling strength J z = 1. We found that the results obtained from the quantum simulator provided by QISKIT [52] differs from the exact results only by statistical errors. This indicates that the quantum circuit can indeed correctly prepare the ground state of the model. Both the exact and simulator results show that the entanglement entropy for the 1D Kitaev spin chain increases with J x , and the derivative is largest at intermediate values of J x (around J x = 0.5). As a comparison, the entanglement entropy obtained from the IBMQ-Almaden quantum computer is shown in the lower panel in Fig. 2(a). The overall features of the entanglement entropy versus J x are captured in the real machine data, i. e. the entanglement entropy increases with J x with a larger derivative at intermediate values of J x . But the quantitative values deviate significantly due to the noise in the data.
In the upper panel in Fig. 2(b), the entanglement entropy for the Kitaev honeycomb model is plotted as a function of coupling strength J y by setting the other two coupling strengths to J x = 0.1 and J z = 1. The difference between the simulator results and the exact results is due to statistical errors. Because the Hilbert space of the ground state is larger for this model, the statistical errors appear more significant. We found that the entanglement entropy for the model on a 2×2 lattice is almost a constant and slightly decreases with increasing J y . The entanglement entropy measured on the machine IBMQ-Almaden is shown in the lower panel of Fig. 2(b). The real machine data deviates quite a bit from the exact results. Similar to the exact results, the entanglement entropy obtained from the quantum hardware is almost constant, but it increases slightly with the increasing J y (unlike the exact results). This larger error is due to the fact that the circuit depth is about three times longer here.

Scanning the Brillouin zone via symmetry-enforced circuits
The ranges of exchange parameters in Fig. 2 cover a regime where both models have quantum phase transitions (in the thermodynamic limit). We observe no signature of these transitions due to finite-size effects (see Appendix F for details). To go beyond this limitation of NISQ machines, we employ a grid of shifted momenta to capture the topological phase transition features in the bond fermion sector. With the help of this method, and given the fact that the contribution from the gauge fermion portion is featureless [43], the entanglement spectra as well as the entanglement entropy can be calculated from the correlation functions in the bond fermion sector only [53,54].

Symmetry-enforced circuit
The Kitaev-inspired models reduce to a blockdiagonal form in momentum space after they have been converted to the Majorana Fermion representation. But, because a Bogoliubov transformation is still needed to fully diagonalize the problem, one must be careful to preserve particle-hole (PH) symmetry when constructing the quantum circuit. Consider a system containing N unit cells. The discrete Fourier transformation to an arbitrary set of N momentum points in the BZ, K = {2πn/N | n ∈ [0, N )} + δk can be written as: where c n,K denotes real-space operators obtained by Fourier transformation of the set K. Here, δk is introduced to shift the momentum points. Physically, it can be achieved by introducing a twist boundary condition as is shown in Fig. 3(a). In the Fourier transformation circuits, the phase shift requires a phase correction [48]. PH symmetry requires that the information from the set with opposition momentum points −K must also be included. Hence, the symmetryenforced Fourier transformation is given by: This indicates that a more complex circuit with twice the number of qubits is required. The demonstration of the symmetry-enforced processes under the particle-hole symmetry for the 1D Kitaev spin chain with N = 2 unit cells is shown schematically in Fig. 3(a). The corresponding quantum circuit is  shown in Fig. 3(b). Since the Bogoliubov transformation mixes ±K, the correct momentum set is swapped before the transformation, and swapped back after.
In general, the symmetry-enforced circuit can not be simplified, but for special momentum shifts fulfilling the conditions stated in Appendix G, a simplification of the symmetry-enforced circuit is possible.

Measuring the entanglement of Kitaevinspired models from correlation functions
For Kitaev-inspired models, ground states have uniform gauge configurations, which contribute trivially to the entanglement entropy. The non-trivial quantum information is stored in the Brillouin zone of the bond fermions, whose entanglement properties can be efficiently obtained from the correlation functions.
To begin, the original spin model is mapped to a conventional fermionic model by the Jordan-Wigner transformation followed by a Majorana braiding. In general, the fermionic Hamiltonian can be written as where A i,j is the hopping matrix and B i,j is the pairing matrix. The explicit form of the two matrices depends on the details of models. The entanglement entropy of this Hamiltonian is easily obtained by introducing Majorana fermions f i = γ 2i−1 +iγ 2i for each physical fermion. In terms of the Majorana fermions, the entanglement entropy of a subsystem of the system is calculated from the eigenvalues of the correlation matrix C [53,54], whose matrix elements satisfy C m,n = 1 2 Tr [ργ m γ n ] = 1 2 γ m γ n with ρ denoting the density matrix of the system. Then the entanglement entropy is given by: with λ n the eigenvalues of C.
To measure the correlation matrix C on a quantum computer, we note that: where σ + n = 1 2 (σ x n + iσ y n ) and σ − n = 1 2 (σ x n − iσ y n ).

Then the four different expectation values become
with α, β = x, y.
With these general layouts, we show how to explicitly calculate the entanglement entropy of bond fermions by using the 1D Kitaev spin chain with two unit cells as an example. The lattice structure of the 1D Kitaev spin chain is shown in Fig. 3(a), where the two unit cells labeled by 0 and 1 in purple form the whole lattice. Here the unit cell 1 (labeled in green) of the left neighbor and the unit cell 0 (labeled in green) of the right neighbor are included to show the boundaries of the system explicitly. We begin the discussion with the Hamiltonian of the model in momentum space in terms of bond fermions (the derivation of this form can be found in Appendix A): Here D is the local gauge field defined in Eq. (10). The Bogoliubov transformation is where e iϕ k = i for the model with two unit cells, independent of k. This means that ϕ k = π/2, independent of k.
and f n f m are then determined from the ground-state expectation values of Bogoliubov fermions where N k is the number of momentum points in the Brillouin zone. Note that C 2m−1,2n−1 and C 2m,2n are the correlations of the real and imaginary parts of the fermions, respectively. As expected, we find C 2m−1,2n−1 = δ m,n and C 2m,2n = δ m,n . We also find that and C 2m,2n−1 = −C 2n−1,2m .
The correlation between the two bond fermions at the unit cells m and n is expressed by the following block matrix: where: This means that for a subsystem containing N sites, the correlation matrix C is given by (31) For the case with two unit cells, the correlation matrix becomes Particle-hole symmetry requires sin θ k = − sin θ −k . This implies that Using this expression, we obtain the exact curves in black shown in Fig. 3(c) and (d). To measure the entanglement properties on quantum hardware as a function of momentum, we need to use the symmetryenforced circuit.

Measuring the entanglement spectrum on quantum hardware by the symmetry-enforce circuit
We use this symmetry-enforced circuit to obtain the entanglement spectra of the 1D Kitaev spin chain containing two unit cells by sweeping across the Brillouin zone and measuring the fermionic correlation functions as described in the subsection 4.2. The results from the quantum hardware are qualitatively consistent with the exact results. When J x < J z , the entanglement spectra is always gapped [see Fig. 3(c)] indicating no edge Majorana modes, while when J x > J z , the entanglement spectrum is gapless corresponding to the Majorana fermion locating at the virtual boundary between the subsystems [see Fig. 3(d)]. We highlight that by using this symmetry-enforced circuit, the quantum hardware can even correctly determine the gap closing point at the high-symmetry point δk = π, as shown in Fig. 3(d).

Determining the quantum phase diagrams on real quantum hardware
The high symmetry points with δk = π in the BZ are the PH symmetric points, and we expect the quantum phase transition to be dominated by the behavior near these points. Hence, we must ensure the simulation includes the effects of these high-symmetry points, in order to efficiently determine the phase diagram on a small cluster. With the help of the symmetry-enforced circuit, this can be done straightforwardly. Using these specific momenta for the systems discussed in Sec. 3 and Fig. 2, we obtain the bond fermion entanglement entropy (by measuring correlation functions) as a function of the exchange parameters. In Fig. 4(a), the bond fermion entanglement entropy is shown as a function of exchange strength J x for the 1D Kitaev spin chain with J z = 1. We can find that the results from the exact solution, quantum simulators, and the IBMQ-Almaden quantum computer all show a discrete jump, and the entanglement entropy S A = √ 2/2, when J x > J z , indicates a Majorana fermion at the virtual boundary of the subsystem. Therefore, the sharp change at J x = J z indicates a topological quantum phase transition. Similar situations are found in the bond fermion entanglement entropy for the Kitaev honeycomb model as shown in Fig. 4(c). Here we set J z = 1 and J x = 0.1 in the calculations. The sharp jumps in S A are found to happen at J z = J x + J y and J y = J x + J z and identify two different topological phase transitions of the Kitaev honeycomb model. This efficient identification of the phase boundary relies on the difference of the topology of the single-particle wave functions in the two different phases. When |J i | < |J j + J k | (with i = j = k ∈ {x, y, z}), the single-particle spec-trum is gapped, and the ground state of the phase is known as the gapped Z 2 spin liquid in the sense that it can be mapped to the toric code model [15]. On the other hand, when |J i | ≥ |J j + J k |, the singleparticle spectrum is gapless, and the ground state of the phase is a spinon chiral topological superconductor [55]. The difference in the single-particle spectrum is revealed by the entanglement entropy of edge modes at the virtual boundaries of the subsystem: if virtual boundaries are introduced by dividing the system into subsystems, then the spinon chiral topological superconducting phase would support edge modes, while the gapped Z 2 spin liquid phase would not. Then using this step as the identifying feature, we are able to reconstruct the correct phase diagram of the two models as a function of exchange parameters as shown in Fig. 4(b) and (d). We have to note that the two gapped phases in the Kitaev honeycomb model have the same properties in their single-particle spectra, but the corresponding quantum states at the high symmetry points are quite different and thus expressed by different occupations in the qubit representation. Because the different qubits have different noise levels, the resultant states have quantitatively different entanglement entropy and result in the asymmetry seen in Fig. 4(d).

Conclusion
In this work, we illustrate how the topological properties intrinsic to a quantum model being simulated provides topological protection, resulting in robust quantum calculations, even on NISQ hardware. In particular, we demonstrate how to prepare the ground states of these Kitaev-inspired models through lowdepth braiding circuits of anyons. By using a particlehole symmetry-preserving methodology that includes high-symmetry points in the Brillouin zone, the simulation identifies the quantum phase diagram of both the Kitaev spin chain and the original Kitaev honeycomb model for systems larger than are otherwise possible on NISQ hardware (when only local quantities are measured). These results directly show how topological protection results in much more robust calculations on NISQ hardware when preparing and determining different quantum phases. The quantum phase diagram was determined by finding discontinuities in the entanglement entropy, which are reduced in magnitude, but still remained easy to identify in this noisy environment. Hence, the benefits of topological-inspired quantum computation can already be realized on NISQ-era conventional quantum computers and this provides a third methodology for topological quantum computation that explores a different realm than quantum error correction or topological quantum computing with topological qubits.

Acknowledgements
This work was supported by the Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Grant No. DE-SC0019469. J.K.F. was also supported by the McDevitt bequest at Georgetown. We acknowledge use of the IBM Q for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Q team. Access to the IBM Q Network was obtained through the IBM Q Hub at NC State. We acknowledge the use of the QISKIT software package [52] for performing the quantum simulations. The authors also thank Sonika Johri for insightful comments and Akhil Francis for some technical expertise. The codes and data are available at https://doi.org/10.5061/dryad.0zpc866vk.

A The exact solution of the 1D Kitaev spin chain
The Hamiltonian of the 1D Kitaev spin chain is given by Here we assume that the odd sites belong to the • sublattice while the even sites belong to the • sublattice (the configuration as shown by the inset in Fig. 2(a)). The factor of 4 is added for convenience, as we will see below. To solve this model, we first map the spin model into a fermionic model by the Jordan-Wigner transformation: where we use the standard mapping |0 = | ↑ and |1 = | ↓ . After fermionization, the Hamiltonian becomes which is identical to the fermionic expression of the Kitaev honeycomb model in the limit J y = 0 [15]. To solve this, we must transform the four-fermion term.
Here i labels unit cells, which have two sublattices • and •. Since the Majorana fermions do not have welldefined occupation numbers, we construct quantum circuits after mapping the Majorana fermions back to conventional fermions. We introduce bond fermions for the two sites connected by the z-bonds as follows: Finally, we re-express the Hamiltonian again and find Because of the translational symmetry, Lieb's theorem indicates that the fluxes D i must be uniform [15].
Since D i commutes with the HamiltonianH 1D , it is a constant of motion. Note that η 2 To diagonalize the Hamiltonian, we start with a Fourier transformation where the Hamiltonian becomes after dropping some irrelevant constants. Here, we introduce the "spinor" This Hamiltonian in momentum space can be easily diagonalized through a Bogoliubov transformation. Doing so gives us the final diagonal form where E k = (J z D + J x cos k) 2 + J 2 x sin 2 k, and b k (b † k ) denotes the annihilation (creation) operator for Bogoliubov fermions with momentum k.

B The Bogoliubov and Fourier circuit for the Kitaev honeycomb model B.1 The implementation of Bogoliubov transformation
The first step is to reverse the Bogoliubov transformation to transform the representation formed by bond fermions in momentum space. In general, the Bogoliubov transformation is a two-qubit operation [46,47]. By assuming that the order of qubits is from bottom to top, the quantum circuit to realize this two-qubit Bogoliubov transformation is drawn in the right middle panel in Fig. 1(c). Performing the pairing between opposite momentum points, the Bogoliubov transformation for the Kitaev honeycomb model is achieved by the two two-qubit Bogoliubov transformations (denoted by orange boxes in Fig. 1(b)) operating on neighboring qubits.

B.2 The implementation of Fourier transformation
We next convert from the bond fermion state in momentum space to real space by implementing two Fourier transformations sequentially along the y and x directions. Based on the fast Fourier transform method, the Fourier transformation for one of the directions containing N = 2 m qubits (m an integer) can be realized with a log-depth circuit and using at most two-qubit quantum gates [47,48]. For the Kitaev honeycomb model on the lattice containing a 2 × 2 = 4 unit cells, the Fourier transformation along the y-direction consists of two two-qubit Fourier transformations for k x = π/2 and k x = −π/2, and the transformation along the x-direction also contains two two-qubit transformations for k y = π/2 and k y = −π/2 respectively. In particular, the two-qubit Fourier transformation consists of a phase delay factor implemented by the single-qubit gate U 1 (φ) shown in the right upper panel in Fig. 1(c) and a two-qubit 'beam-splitter' operation implemented by the rest of the right upper panel in Fig. 1(c) except for U 1 (φ). Beginning from the states obtained after the Bogoliubov transformation, we need to reorder the qubits according to the momentum points as (−π/2, −π/2), (−π/2, π/2), (π/2, −π/2), and (π/2, π/2), which is achieved by the two swap operations after the Bogoliubov boxes in Fig. 1(b), so that the Fourier transformation in the y-direction can be implemented for the two sets of neighboring qubits, respectively, shown by the two green boxes following the two swap operations in Fig. 1(b). To facilitate the Fourier transformation in the x-direction, we need to switch the order of the two qubits for (−π/2, π/2) and (π/2, −π/2), which is achieved by the swap operation following the two green boxes in Fig. 1(b). In this way, the Fourier transformation along the x-direction can be implemented by the two green boxes following the swap box as shown in Fig. 1(b).
C The braiding circuit of the 1D Kitaev spin chain As we have discussed in the main text, the details of the braiding circuit are determined by the constraints set by local coupling terms of the Hamiltonian. The constraint set by the z-bond is general for the Kitaevinspired models, so here we will focus on the constraints from other coupling terms in the 1D Kitaev spin chain.
To make the discussion explicit, we focus on the x-bond connecting sites 2 and 3. The original Hamiltonian connecting across the x-bond is and it is re-expressed in terms of bond and gauge fermions asH where the braiding operations on c 1 and c 2 give the bond fermion f 1 and the gauge fermion g 1 , and similarly the braiding operations on c 3 and c 4 give the bond fermion f 2 and the gauge fermion g 2 . Hence, the states considered here contain four qubits. Start with the n g1 = n g2 = 1 case. The relevant states forH x are |n f1 n f2 n g1 n g2 = |n f1 n f2 11 f . They are transformed underH x as From the general form of braiding provided by Eq. (11), the states |n f1 n f2 11 f are transformed to |n 1 n 2 n 3 n 4 c (in the representation of fermions denoted by c operators) as follows: 11 α 01 01 |0001 c + α 00 11 α 10 01 |0010 c + α 11 11 α 01 01 |1101 c + α 11 11 α 10 01 |1110 c . (48) The action of H x on these states gives The consistency between n f1 n f2 11|H x |n f1 n f2 11 f and n f1 n f2 11|U † B H x U B |n f1 n f2 11 f implies (α 01 01 ) 2 = (α 10 01 ) 2 = −α 00 11 α 11 11 , (α 00 11 ) 2 = (α 11 11 ) 2 = −α 01 01 α 10 01 , α 01 01 = −α 10 01 , α 00 11 = −α 11 11 . (50) This yields α 01 Since α 01 01 and α 11 11 are imaginary, there are two intra-fermion braidings acting on the same fermion to adjust the global phase by ±i. The interfermion braiding introduces a π/2 phase difference. To increase the phase difference to π, one of the intrafermion braidings is performed after the inter-fermion braiding. We still need to determine on which fermion the intra-fermion braiding acts. By using Eq. (50), the transformation between the two representation can be written as: But, the phase difference introduced by the interfermion braiding operations are opposite for |01 and |11 . To keep the form of Eq. (51) up to a global minus sign, the only choice is that the intra-fermion braiding operation is applied to the first qubit, the occupation of which is different for the two states. Figure 5: The braiding circuit used for the 1D Kitaev spin chain.
From the analysis provided above, the braiding circuit must satisfy the following: 1. the circuit is constructed by one inter-fermion braiding and two intra-fermion braidings; 2. one of the intra-fermion braidings is performed after the inter-fermion braiding; 3. the two intra-fermion braiding operations act on the first qubit.
A circuit fulfilling these requirements is shown in Fig. 5. One can apply a similar analysis to the case with n g = 0. One finds the same requirements in that case.
which is constructed on a ladder lattice as shown in the inset of Fig. 6(b) with the subscript l denoting the two legs of the ladder. The JW string is defined so that the odd sites belongs to the • sublattice and the even ones are the • sublattice. Obviously, this model is a spin 1/2 model fulfilling the requirements stated in [40] and thus another example of Kitaev-inspired model. By the Jordan-Wigner transformation, we have that where we have identified |0 = | ↑ and |1 = | ↓ .
Hence, we find: Inserting the Jordan-Wigner transformation for σ z , we find the fermionic Hamiltonian: where we used the fact that all the odd sites are • sites, and all the even sites are • sites. To simplify the four-fermion interaction term, we assume that the fermions at the • sites are decomposed into Majorana fermions by following: while the fermions at the • sites are decomposed as: the sub-lattice symmetry make the k = 0 points always gapped, so Pf[H(k = 0)] does not change sign. This suggests that the quantum phase transition in the 1D BCS-Hubbard model depends on the sign of Pf[H B (k = π)], which is equal to the sign of J z . Therefore, the quantum phase transition in the 1D BCS-Hubbard model happens when J z changes sign, which is consistent with the numerical results shown in [39]. A straightforward check is that when J z changes sign, the filling at k = π changes as we expect it would. In this section, we prove that the quantum circuit proposed in the main text ( Fig. 1) correctly prepares the ground state of Kitaev-inspired models. The system sizes of the models considered here are the same as those shown in the main text. The lowest energy of the finite-size cluster is determined by exact diagonalization, and then compared to the results obtained by exact solution and by quantum circuit simulation. Results for the 1D Kitaev spin chain and the Kitaev honeycomb model are summarized in Fig. 7. It turns out that the lowest energies obtained from exact diagonalization are identical to those obtained from the exact solution. Results from the quantum simulators have small deviations due to finite-number sampling.

F Finite size effect on entanglement entropy
The advantage of calculating entanglement entropy with correlation functions is that it allows us to study the finite-size effects on the entanglement entropy. In Fig. 8, we show how the entanglement entropy of the bond fermions changes with increasing lattice size. We note that when the bond fermion lattice size N is smaller than 10, the entanglement entropy increases with the coupling strength J x monotonically, and no obvious signature appears when J x passes through the quantum phase transition point. If we further increase the lattice size (for example the N = 20 case), the non-monotonic behavior of the entanglement entropy with increasing coupling strength starts to appear. However, an accurate determination of the quantum phase transition requires the bond fermion lattice size N to be larger than about 80.
Similar procedures can be used to find the entanglement entropy of the bond fermions of the Kitaev honeycomb model. How the entanglement entropy changes with the increasing system size is shown in Fig. 9. It turns out that an accurate determination of the quantum phase transition requires a system size larger than L x × L y = 64 × 64.  Though in general, the symmetry-enforced circuit cannot be simplified, we show here that the circuit can be simplified for special momentum shifts. These special momentum shifts include the high symmetry points.
For a system containg N unit cells, the discrete Fourier transformation and the particle-hole symmetry constraint requires the momentum points to be in the following set: where N = 2 n , with n a positive integer. As demonstrated in the main text, if we shift the momentum points by δk, their particle-hole partners must shift by −δk to maintain particle-hole symmetry. Based on this, we note that δkN/2π = 1/2 is a special point, which makes the set of the momentum points by shifting δk to be identical to those that shift by −δk. Therefore, if the gap is vanishing at the high-symmetry points, the symmetry-enforced circuit can be decomposed into two identical copies. The case with N = 4 is shown in Fig. 10. The simplification requires three steps, which are shown in Fig. 10 (b) from left to right. In the first step, we identify the momentum points labeled by 3 + and 4 − , which is obvious from Fig. 10 (a); in the second step, due to the gap vanishing at the high-symmetry points, the two-qubit Bogoliubov transformation becomes two separate single-qubit identity operations, so we can change the position of the momentum points 4 + and 3 − at the cost of a trivial global phase; in the last step, we just renumber the −δk-shifted momentum points by using the numbers in red shown below the original labels (in purple) in Fig. 10 (a). Then one can easily find that the result yields two identical copies. We emphasize that the vanishing gap at the high-symmetry points is essential in making this simplification possible.
It turns out that the Kitaev honeycomb model fulfills this requirement, so the simplified circuit can be used to determine the different phases of the model. The details of the simplified circuits are shown in Fig. 11.