Investigating the Exchange of Ising Chains on a Digital Quantum Computer

The ferromagnetic state of an Ising chain can represent a two-fold degenerate subspace or equivalently a logical qubit which is protected from excitations by an energy gap. We study a a braiding-like exchange operation through the movement of the state in the qubit subspace which resembles that of the localized edge modes in a Kitaev chain. The system consists of two Ising chains in a 1D geometry where the operation is simulated through the adiabatic time evolution of the ground state. The time evolution is implemented via the Suzuki-Trotter expansion on basic single- and two-qubit quantum gates using IBM's Aer QASM simulator. The fidelity of the system is investigated as a function of the evolution and system parameters to obtain optimum efficiency and accuracy for different system sizes. Various aspects of the implementation including the circuit depth, Trotterization error, and quantum gate errors pertaining to the Noisy Intermediate-Scale Quantum (NISQ) hardware are discussed as well. We show that the quantum gate errors, i.e. bit-flip, phase errors, are the dominating factor in determining the fidelity of the system as the Trotter error and the adiabatic condition are less restrictive even for modest values of Trotter time steps. We reach an optimum fidelity $>99\%$ on systems of up to $11$ sites per Ising chain and find that the most efficient implementation of a single braiding-like operation for a fidelity above $90\%$ requires a circuit depth of the order of $\sim 10^{3}$ restricting the individual gate errors to be less than $\sim 10^{-6}$ which is prohibited in current NISQ hardware.


I. INTRODUCTION
The race for developing fault-tolerant quantum computing has only intensified in the last few years [1,2]. Recent experiments show quantum hardware are beginning to show their advantage, although in very strict cases, over classical computers [3]. While the main issue remains to be decoherence and local perturbations [4,5], there has been significant progress in two directions. First, regarding qubits that are engineered to autonomously correct errors with methods such as bosonic codes [6], and second, with proposed qubits based on non-abelian anyons that exhibit non-trivial topological statistics, which are immune to decoherence and local perturbations, such as Majorana qubits [7]. While these novel qubits are being developed, it is instructive can study their dynamics by simulating them on available noisy intermediate-scale quantum (NISQ) hardware.
Kitaev recognized that a 1D chain of spin-less fermions with superconducting pairing, known as a Kitaev chain, can host non-abelian anyons [8,9]. Recently, signatures of non-abelian anyons have been claimed to be experimentally observed [10,11]. The braiding, or series of exchanges, of Majorana fermions is the basis of topological quantum computing [12][13][14]. While in Kitaev chains, the ground states' degeneracy is not lifted under the influence of local perturbations, these chains are difficult to realize. A more accessible model is the Ising chain. The mapping between the Kitaev chain and the Ising chain can be done using the Jordan-Wigner transformation. However, under such mapping topological protection is lost. Nonetheless, one might be able to study Majorana-like systems implemented on Ising chain based models [15] as the mapping between the Ising Hamiltonian and the Kitaev Hamiltonian is exact and therefore the excitation gap is the same. While topological protection is lost, it has also been shown that some non-topological immunity to local noise and error can be achieved depending on the physical implementation, e.g., in photonic systems [16], and trapped ions [17]. Moreover, due to the ferromagnetic order of the ground state of the Ising chain, one can see that even in the presence of excitations such as bit flips, the ground state can be recovered through a majority rule [18].
Simulations of braiding Majorana fermions in topological superconductors are discussed in [19]. Further, dynamics of Majorana fermion braiding on a T-junction geometry using scaled two-qubit quantum gates on IBM quantum computers are investigated in [20]. Numerical simulations, as done in [21], even when performed with braiding approximations such as weaving, proved to be inefficient. The simulation of Majorana fermions has also been studied as part of systems such as a cavity QED lattices in [22] but without performing braiding or operations on the modes. Recently, Backens et al. [15] has shown how Majorana fermion braiding can be studied on Ising chains through fermionic mapping with a 1D geometry which is debatable given the non-topological nature of the system or particles. Here we study the system in terms of quantum gates. To simulate the quantum system in an appropriate environment practically and efficiently, we simulate the system on a quantum simulator running on a classical computer with the 1D geometry described in Backens et al. [15] using basic quantum logic gates to investigate the implementation of such a braiding protocol.
We provide a detailed analysis of the implementation of a braiding-like exchange protocol as well as the optimization and estimation of the resources required to achieve high fidelity on NISQ hardware. Even though [15] and [20] claim to preform "braiding of Majorana zero modes" using similar Ising chain-based models, we use the term "braiding" loosely in this work since true braiding involves an operation in 2D physical space and involves an intrinsic notion of chirality. Chirality plays a major role in the braiding of Majorana zero modes and depends on local characteristics of the system and the exact procedure used to move the modes [23,24]. In our model, braiding could be thought as the usual 2D braiding mapped out onto a local 1D system, an exchange of two particles.
For the digital quantum simulations, we map our physical system onto a set of qubits and implement our simulation using IBM's open-source software package, Qiskit, and run the simulation on IBM's Aer QASM simulator. However, our approach can be applied to any digital quantum computer. Current NISQ era quantum hardware have limited circuit depth which can be seen as the number of non-simultaneous gates between the input and the output. We also attempt to minimize the circuit depth and estimate the minimal required resources for a given system.
For Hamiltonians that involve multiple noncommuting parts, there is no simple circuit representation of the system's evolution operator.
To circumvent this issue, we will approximate the evolution using the Suzuki-Trotter decomposition, known as the Trotter expansion, described by Hatano and Suzuki [25]. This expansion approximates the total evolution operator using the product of the evolution operator for each component over a short time ∆t. Since this description is only an approximation, we have to choose the order of the expansion and the steps' duration to optimize the fidelity and reduce the circuit depth.
In this paper, we start by introducing the model and the braiding protocol in section II followed by the digital Hamiltonian simulation in Section III. Then, in section IV, we discuss the system's expected Suzuki-Trotter error before presenting the results and analysis of the digital simulation and the contributing errors.

II. ISING CHAIN IMPLEMENTATION OF BRAIDING-LIKE EXCHANGE PROTOCOL
Majorana fermions have been predicted to exist in the 1D Kitaev model [9], which consists of spinless fermions with an attractive interaction. This system can be mapped to an Ising chain through the Jordan-Wigner transformation [26]. This mapping is exact, and the topological implications have equivalent counterparts in the spin space. However, the topological protection is lost in the transformation and the Ising chain is sensitive to local perturbations.
The Ising model describes a lattice of spins exhibiting nearest-neighbor interactions and interactions with an external transverse magnetic field. In the case of a homogeneous nearest-neighbor interaction, the Hamiltonian of the system can be written as [27]: where J represents the strength of the interaction between neighboring spins and h n represents the strength of the magnetic field at site n.
Given a weak magnetic field compared to the strength of the pair interaction, h n < J, the ground state of the system is ferromagnetic; when the opposite is true, h n > J, the system is paramagnetic. To create welldefined ferromagnetic and paramagnetic domains on an Ising chain, we set the magnetic field, h n , for each site n as h f erro J and h para J, respectively, as shown in Fig. 1.
1. An Ising chain with homogeneous nearest-neighbor interaction (J) and local transverse magnetic field (h). By controlling the transverse field's strength, the chain is split into two domains: left side is ferromagnetic (h J), and right side is paramagnetic (h J).
In the case of a three-site chain, the two states, |↑↑↑ and |↓↓↓ , of the ferromagnetic domain of an Ising chain, are degenerate and could span a qubit space. Consequently, the Ising chain's ferromagnetic domain can be used to mathematically represent the topological domain that hosts Majorana fermions [15]. For the ferromagnetic states of the Ising chain the mapped states, which can be either be empty (|0 ) or occupied (|1 ), would then follow: where |↑↑↑ represents a ferromagnetic domain of length 3 with all its spins in an |↑ state. Because of this mapping, we will now refer describe the topological domain as ferromagnetic domain.
Braiding Majorana fermions requires preforming exchange in real space as seen in Fig. 2a. Such an operation is strictly forbidden in 1D. When the Jordan-Wigner transformation is applied to multiple chains, issues arise that include non-local interactions and operator commutativity of different chains. To fix these issues, one can introduce an extra spin or coupler, S α , which ensures that the fermions of different chains properly anti-commute [15]. For the particular case of Ising chains, Backens et al. have shown [15] that with two Ising chains connected by a coupler and through a rotation of the coupler midway during the braiding procedure, it is possible to operate on a 1D chain. In the following, we will focus on this case, presented in Fig. 3, which is described by the following Hamiltonian: where σ z α (−1) and σ z α (1) corresponds to the last and first site of a chain α, respectively, with H CI representing the coupler interaction.
In addition to the coupler manipulation, the braiding procedure involves moving the ferromagnetic region along the chain. This can be achieved through manipulation of local fields h n . It is important here that each step of the procedure is done slowly enough, in an adiabatic fashion, so that the system remains in its ground state. By adiabatically decreasing h n on a site adjacent to the ferromagnetic domain from h para to h f erro , the ferromagnetic domain can be elongated by one site from right (left) [28]. The reverse process can be used to shorten the domain from left (right). Consequently, the ferromagnetic region can be moved along the chain by successively elongating the ferromagnetic domain from one end and then shortening from the other end, as shown in Starting from configuration a, the field on the fourth site is ramped down adiabatically from hpara to h f erro . As a consequence in b, the ferromagnetic region spans 4 sites. Ramping the field up on the leftmost site, the ferromagnetic region is contracted by one site as shown in c.
The complete braiding process is described in Fig. 5, starting with the ferromagnetic domain on the far end of the left chain. The coupler is rotated around the xaxis by π/2 so that it is in a superposition of S z = +1 and S z = −1. The ferromagnetic domain is then moved adiabatically, as described above, towards the far end of the right chain. On the way, it interacts with the coupler. The coupler is then rotated around the y-axis with angle θ, and the ferromagnetic domain is moved adiabatically back to the far end of the left chain.
The overall braiding operation induces a rotation of −φ around the z-axis of the (|0 , |1 ) subspace, whose states are given in Eq. 2. In particular, if we initialize the ferromagnetic region in |↑↑↑ and perform of rotation of φ = π the final state of the system is expected to be |↓↓↓ .

III. DIGITAL SIMULATION OF BRAIDING IN 1D
The system's time evolution is described by a circuit of gate operations in the qubit space [29,30]. In this case, we map each spin in the Ising chain to a qubit and use an extra qubit to describe the coupler required by the protocol since we are dealing with spin-1/2s, which are essentially qubits. Simulating the braiding protocol described in II requires simulating the Hamiltonian dynamics of the system and performing adiabatic changes to its parameter. Next, we discuss each aspect and introduce the parameters that can be used to optimize the procedure.
The time evolution corresponding to the Hamiltonian of the system during a time interval ∆t is written as U = exp(−iH∆t). Since the circuit implementation is limited to single-and two-qubit gates, the Hamiltonian is decomposed into a minimal set of four terms, i.e. H = X H X , that can be exponentiated efficiently, that is with a fixed number of gates irrespective of the system size with cost O(1). These summands are the nearest-neighbour coupling interaction for even sites H J,e = −J n=2k σ z n σ z n+1 as well as odd ones H J,o = −J n=2k+1 σ z n σ z n+1 , the Zeeman terms H Z = − n h n σ x n , and the coupler interaction term H CI = −J C S z σ z 1 (1)σ z 2 (1). Assuming that the time step is sufficiently small, i.e., h∆t 1, one can approximate the time evolution operator as X exp(−iH X ∆t) using the Suzuki-Trotter expansion. The adiabatic condition and the Trotter error bounds are discussed in the next section.
The nearest-neighbor interaction terms of the Ising Hamiltonian, Jσ z n σ z n+1 described in Eq. 1, can be realized as follow. First, a controlled-NOT (CNOT) gate is added between qubit q n (control) and q n+1 (target). We then use an R z (J ∆t) gate to rotate q n around the z-axis. Finally, we apply another CNOT gate between q n (control) and q n+1 (target) to complete the time evolution of the nearest-neighbor interaction, i.e. exp(−iJσ z n σ z n+1 ∆t). The on-site magnetic field, h n σ x n , the second term of Eq. 1, is implemented using a single-qubit rotation in the x direction [29,30] as exp(−ih n σ x n ∆t) ≡ R x (h n ∆t). The circuit describing the time evolution of H j,e , H j,o , and H Z , for two chains of size three, is presented in Fig. 6. Note that to reduce the depth of the circuit, we apply the pair interaction to the q n , q n+1 pair and to the q n+2 , q n+3 pair simultaneously, allowing to describe the complete interaction of the chain with only two layers of gates. Finally, To implement the time evolution of the coupler interaction seen in Eq. 4, i.e., exp(−iJ C S z σ z 1 (1)σ z 2 (1)∆t), we use a combination of CNOT, Hadamard and R z gates as shown in Fig. 7 as described in [29]. We present the simulation of 1D braiding with a rotating coupler representing a phase gate, R z (π), in the basis of the logical states defined in Eq. 2. The braiding process involves four stages:

Initialization of Chain and Coupler
We begin by creating a circuit of qubits with one qubit for each site and an extra qubit for the coupler. We initialize the ferromagnetic domain in either one of the defined two logical states: |0 , |1 , or in one of their superpositions. To start from a logical state, we use a series of Hadamard, CNOT, and R y (− π 2 ) gates applied to the qubits of the ferromagnetic domain. Regarding the superpositions, the |↑↑↑ is trivially accessible as the system's default state. Hadamard gates are applied to all the qubits in the paramagnetic domain, and the coupler is initialized by applying an R x ( π 2 ) gate on the qubit. The initialization sequence relies on the fact that the chosen value of h in the ferromagnetic domain, h f erro , is much smaller than J, making the initial state a good approximation of the ground state. A similar argument applies to h para used in the paramagnetic domain.

Moving Ferromagnetic Domain from Left to Right
Given the initial configuration where the ferromagnetic domain is on the far left and the final configuration is on the far right, the intermediate field configurations are determined. To move the domain one way or the other, the domain is extended by one site on the side closer to the direction of the movement and contracted by one site on the other side, as seen in Fig. 4. In contrast to what was discussed in II, we perform the extension and the contraction of the ferromagnetic domain simultaneously rather than sequentially. We have verified that this does not degrade the fidelity of the process and allows for shallower circuits. As a consequence, each intermediate configuration represents the ferromagnetic domain's movement by one site. The on-site fields are updated in small steps between each configuration, where the Trotter expansion is used to describe the evolution. We chose a simple linear update for the on-site field, which in some cases could suffer from the large values required for h para . The size and frequency of the field update affect the adiabaticity of the evolution.

Mid-Braiding Manipulation
After the ferromagnetic domain is moved to one side, the coupler is rotated around the y-axis using an R y gate in a series of steps. With each rotation step, the site's interaction with the slightly rotated coupler is represented using the Trotter evolution. This approach assumes that the coupler rotation is not much faster than the Hamiltonian evolution.

Moving Ferromagnetic Domain Back from Right to Left
Following the moving protocol described in III 2, the ferromagnetic domain is moved from right to left. This concludes the R z (π) braiding process for which an initial state of |↑↑↑ and a total rotation of the coupler by π should lead to |↓↓↓ .
Based on the description of the simulation provided above, we can identify the following relevant parameters that could be susceptible to influence the fidelity of the procedure: • Trotter step -∆t: the duration for which we evolve the system in each iteration of the Trotter procedure. Using shorter steps should increase the approximation accuracy but would increase the overall depth of the circuit.
• Phase separation: represents how well the ground states of the ferromagnetic and paramagnetic region, defined by the simple states used in the initialization step, depend on h f erro /J and h para /J which need to be, respectively, much smaller and much larger than unity. Since the procedure needs to remain in the system's degenerate ground state, it is important that this approximation is accurate.
• Field update rate -∆h/T : the update of the onsite magnetic field when converting a site needs to be adiabatic. Consequently, we expect slower updates to give better fidelity. Since we perform step updates, this parameter is controlled by both the magnitude of the update (∆h) and their temporal distance (T ). There may exist an optimal choice of those two for a given rate. In particular, large steps are unlikely to yield accurate results.
• Coupler interaction strength -J C : The coupling strength between the coupler and the two chains, which corresponds to the chemical potential in a Kitaev chain, may influence the braiding's fidelity.
In particular, the back-action of the coupler rotation on the adjacent ferromagnetic domain may be mitigated by lower coupling strength.
• Coupler rotation speed -Γ: how fast the coupler's rotation is performed during the mid-braiding manipulation, after moving the ferromagnetic domain from left to right, may influence the fidelity of the braiding process.

IV. SIMULATION PERFORMANCE AND FIDELITY
Here we present the results of implementing the phase gate, R z (π) in the {|0 , |1 } subspace, using the braiding protocol discussed in Section III. To investigate the performance of the Hamiltonian simulation and the simulation's fidelity, F , we start by looking at five different situations and measure the fidelity in each case on a chain of 6 sites, N s = 6, with a ferromagnetic domain of 3 sites: • Two separate Ising chains with no coupler spin and simply performing the movement of the ferromagnetic region from the left to the right and back. For this operation, we start with the ferromagnetic domain initialized either in |↑↑↑ or |0 .
• Two Ising chains connected through a coupler and performing the same operation as above with the coupler initialized along the +x-axis, as in the braiding protocol. We consider again two initial states: |↑↑↑ or |0 .
• The full braiding applied on |↑↑↑ which involves the coupler's rotation after the first translation.
In each case above, except for the braiding, we expect to recover the original state at the end of the procedure. All the following quantum circuit simulations were run on IBM's Qiskit Aer QASM Simulator with 10000 shots. The QASM simulator acts as an ideal quantum computer exhibiting only sampling errors. Another main source of error in the simulation is attributed to the termination of the Trotter expansion as will be discussed. Other types of errors that arise from NISQ hardware are discussed also in the last section. Through a process of optimization starting with reasonable estimates, the optimized parameters, reported in Table I, are obtained and used for all simulations where we consider units such that = 1. The nearest-neighbour interaction strength is always set at J = 1 throughout all the simulations. Fidelity ∆t h f erro hpara ∆h T JC Γ > 99% 0.2 0.01 5.0 0.05 2 0.3 π/3 > 90% 0.7 0.01 1.5 0.1 2 0.3 π/2

A. Parameter Dependence and Optimization
Results for the impact of different system and evolution parameters, i.e. trotter time step size(∆t), value of the magnetic field on sites in the paramagnetic domain(h para ), field update rate(∆h) and the strength of the coupler interaction (J C ), on the simulation's fidelity are presented in Figs. [8][9][10][11]. Different curves correspond to each of the five situations described above to show the impact of the initial state and the presence of the coupler on the fidelity of the operations.  (2J)T /∆h. As seen from Table I, the optimized parameters are sufficiently far from being non-adiabatic and, therefore, the drop in the fidelity is mostly related to the truncation error of the Trotter expansion. As seen from Eq. 10 the error is proportional to ∆t. The sudden drop in the fidelity can therefore be attributed to the non-commutativity of the summands of the Hamiltonian which becomes more pronounced as the time step increases. Similarly, the dependence of the fidelity on the initial state can be understood through the Trotter expansion. The leading Trotter's error is proportional to the commutators of the summands of the Hamiltonian which, in turn, have different expectation values depending on the initial states.  Table I.
The dependence of the fidelity on the Zeeman field of the paramagnetic domain, h para , is depicted in Fig.  9. At small values of the field, i.e. h para < 2, the fidelity drops drastically, especially for braiding. This is due to the fact that a successful translation and braiding operation requires a clear distinction between the ferromagnetic and paramagnetic domains which is achieved through a high Zeeman field, i.e. h para J. This can also be understood through the fact that the Zeeman field in the Ising model is equivalent to the chemical potential in the Kitaev model [15] where it is used to define trivial/non-trivial regions locally. Interestingly, even though the energy difference between the states |↑↑↑←←← and |↑↑←←←← is |J − h|, which if the middle spin is in a ferromagnetic region is about J, stabilizing the paramagnetic region requires fields leading to a much larger energy separation. Conversely, for a very large value of the field, h para > 50, the fidelity starts to drop (not shown in the figure). The reason is that the Trotter error is quadratic in h para since h para appears in  Table I. field update increments ∆h, which is proportional to the update rate ∆h/T , in Fig. 10. Optimum fidelities are achieved only for update rates ∆h/T smaller than 0.25. The drop in the fidelity at ∆h > 0.25 suggests the onset of the violation of the adiabatic condition as the update rate becomes comparable to the excitation gap of the system. The Trotter error is inversely proportional to ∆h because smaller update rate corresponds to a larger number of steps. However, this asymptotic behaviour is not observed for the values shown in the Fig. 10.  Table I.
Lastly, we investigate the effect of the coupler interaction strength J C on the fidelity. Fig. 11 shows the fidelity dropping at lower values, i.e., J C < 0.75, which can be understood as the breaking of the ferromagnetic/antiferromagnetic interaction mediated by the coupler which drives the braiding process. As the coupling between the chains is weakened, it gets more difficult for the ground state to rotate in the degenerate subspace and therefore degrades the braiding process. For large values of J C where it becomes considerable compared to J, the fidelity drops as well. As discussed in the previous section, this may be related to the back-action of the coupler rotation. This is also consistent with the Trotter error containing the commutator [H Z , H CI ] which is proportional to J C . The optimal value of the coupling strength is therefore chosen according to Fig. 11 from the interval 0.25 < J C < 0.75.
The coupler rotation speed Γ appeared to only minimally affect the fidelity of the braiding operation in our simulations. This is interesting since Ref. [15] suggests only to manipulate the coupler when the ferromagnetic domain is far from the coupler, and here we find that even when the domain is adjacent to the coupler, in a chain of 6 sites and a ferromagnetic domain of 3 sites, relatively fast manipulations remain possible without losing fidelity.
FIG. 12. The circuit depth, defined as the number of nonsimultaneous gates executed, for optimum braiding fidelity, F , as a function of the size of the system. The optimum fidelity achieved is noted above the bars, and the ferromagnetic domain's size is half the system's size. The depth of the circuit can be estimated in terms of the simulation parameters. According to the circuits in Figs. 6 and 7, the total circuit depth of each Trotter time step that includes the full Hamiltonian interaction is D T = 12, as seen from the number of non-simultaneous gates and calculated by the Qiskit simulator. Shifting the ferromagnetic domain by one site requires h para /∆h field updates each with T /∆t Trotter steps. The rotation of the coupler requires π/Γ updates for each Trotter step. Therefore, the braiding operation, i.e., shifting the ferromagnetic domain by N s sites along with the coupler rotation, results in the circuit depth Here N s is the total number of Ising spins in the system. We note that this is an upper bound to the depth, and a shallower circuit can be obtained in practice by combining commuting gates. For a given fidelity F and system size N s , we minimize the circuit depth by examining T /∆t, h para /∆h, and π/Γ. A braiding fidelity of F > 99% can be achieved through optimization with the depth being as low as 56518 with parameter values shown in Table I. Using our results, we also found parameters, detailed in Table  I, yielding above 90% fidelity with a depth reduced to 2679 which is significantly more accessible to the current NISQ hardware.
For larger systems, we investigated the depth while maintaining F ≈ 99% for a system of size up to 22 sites, corresponding to a ferromagnetic domain of 11 sites. The circuit depth is seen to increase linearly with the size of the system reaching a circuit depth above 200k for a system of 22 sites while the fidelity slightly drops as expected, as seen in Fig. 12.

A. Suzuki-Trotter Error
The time evolution of the system during the time period T can be decomposed as U (0, T ) = U (T − ∆t, T ) · · · U (∆t, 2∆t)U (0, ∆t), where each shorter time evolution is written as a time-ordered exponential Here H(t) is the time-dependent Hamiltonian of the system. During a time step ∆t, the field at the ends of the chain changes by ∆h/(T /∆t) which should be sufficiently smaller than the excitation gap, i.e., 2J, to make sure one is in the adiabatic regime, that is ∆t (2J)T /∆h. When this condition holds, the integral above can be replaced with H(t j )∆t. Using the Baker-Campbell-Hausdorff formula along with the decomposition H(t j ) = X H X , one obtains where U (∆t) = X exp(−iH X ∆t) is the product formula implemented in the quantum circuit. The leading error is written as [31,32] U (∆t) − U (∆t) ≤ 1 2 2iJ(h j σ y j σ z j+1 + h j+1 σ z j σ y j+1 ), This error corresponds to a single Trotter step during which the field is changed by ∆h/(T /∆t). The total error of changing the field by ∆h during time period T adds a factor of T /∆t. Ramping up the field to h para on each site requires another factor of h para /∆h. Lastly, the movement of the ferromagnetic domain from one end of the chain to the other end and back adds a factor of N s to the error bound. Therefore, the total leading error due to the truncation of the Suzuki-Trotter expansion is This error bound demonstrates the asymptotic behaviour of the error in terms of the parameters of the simulation. Overall, the fidelity of the simulations seems to generally depend on both the adiabatic condition and the truncation error of the Trotter's expansion. However, plugging the default parameters into Eq. 10 results in a loose error bound which is several orders of magnitude greater than what is seen from Figs. 8 to 11. As discussed in Refs. [33,34], in a Hamiltonian with local interactions, the error may not be dominated by the second order term and therefore the true experimental error observed in the simulations are by orders of magnitude smaller than the second order error bounds.

B. NISQ Hardware Error
In addition to errors arising from the adiabatic conditions and the Trotter expansion, which is intrinsic to the algorithm and limits the ultimate fidelity of the simulation, there are extrinsic errors that arise from the nonidealities of the hardware implementation. Using Qiskit's predefined noise models to introduce bit-flip, phase, and measurement errors on all single-qubit and two-qubit gates for all qubits we are able to estimate the simulation's tolerance to the NISQ hardware's noise. We study the optimized case where N s = 6 where the errors are applied to all the qubits and gates in the system. Fig.  13 illustrates the fidelity as a function of the error for different types of errors. The value of the error being applied is the same for both single-qubit and two-qubit gates. As expected, the errors are found to come mostly from two-qubit gates rather than single-qubit gates. We further observe that measurement errors alone do not have a significant effect on the fidelity as seen in Fig. 13 given that the current achievable measurement errors are of the order of ∼ 10 −2 [35]. Phase and bit-flip errors, on the other hand, are seen to be the dominant sources of error in the system. Achieving fidelities above 90% requires quantum gates with an average phase and bit-flip errors of < 10 −6 . This agrees with the fact that the circuit depth is of the order of 10 5 and suggests that high fidelities cannot be achieved on current NISQ hardware as the errors are several orders of magnitude greater than what is needed for a circuit depth corresponding to a single braid protocol without taking into consideration other sources of errors.

VI. CONCLUSION
In a one-dimensional system of two Ising chains, we studied a braiding-like operation that moves the wavefunction inside the degenerate subspace of the ferromagnetic state representing a logical qubit. The braiding operation occurs in the parameter space of the Ising chain which resembles the real space braiding of zero modes in a Kitaev chain. The adiabatic time evolution of the system was implemented by using the Suzuki-Trotter expansion and we achieved a fidelity of > 99% for systems of up to 11 sites. We studied the trends of the fidelity of the simulations as a function of different system evolution parameters and obtained optimum sets of parameters for high fidelity as well as high efficiency implementations.
The circuit depth scales linearly as a function of system size N s , but one should note that the depth of those circuits is large. Our optimum parameters yield a depth of ∼ 57000 for N s = 6, far exceeding the capabilities of the current state of the art devices, which can handle the depth of a few 10s [3]. In an effort to bridge this gap, the complexity of the system could be significantly reduced and with the parameter values detailed in Table I, to achieve ∼ 90% fidelity with a depth reduced to about ∼ 2700. This is a considerable improvement but remains beyond present-day capabilities.
We showed that the Trotter error bounds are, by orders of magnitude, looser than the error observed in the simulations. This suggests that the error might not be dominated by the second order terms in this expansion and therefore a higher fidelity with a lower circuit depth is achievable. Further, for our simulations bit-flip and phase errors contribute the most error to the simulation and a NISQ hardware needs to have an error of the order ∼ 10 −6 per gate for a optimum fidelity or above which is orders of magnitude lower than errors of quantum hardware we have access to currently.