Simulating the exchange of Majorana zero modes with a photonic system

The realization of Majorana zero modes is in the centre of intense theoretical and experimental investigations. Unfortunately, their exchange that can reveal their exotic statistics needs manipulations that are still beyond our experimental capabilities. Here we take an alternative approach. Through the Jordan–Wigner transformation, the Kitaev's chain supporting two Majorana zero modes is mapped to the spin-1/2 chain. We experimentally simulated the spin system and its evolution with a photonic quantum simulator. This allows us to probe the geometric phase, which corresponds to the exchange of two Majorana zero modes positioned at the ends of a three-site chain. Finally, we demonstrate the immunity of quantum information encoded in the Majorana zero modes against local errors through the simulator. Our photonic simulator opens the way for the efficient realization and manipulation of Majorana zero modes in complex architectures.

and e. The corresponding experimental final states (after the second dissipative evolution DE0) with the dark green dots labeled as 1 , cyan dots labeled 2 , magenta dots labeled as 3 , dark yellow dots labeled as 4 , violet dots labeled as 5 , and navy dots labeled as 6 in the Bloch spheres for the cases of flip-error protection and phase-error protection, respectively.
The black dots in the poles of the Bloch spheres represent the corresponding theoretical predictions with the states |xxx (+Z respectively. b. The corresponding fidelities of the states in a. c. Theoretical predictions (black dots) and experimental results (magenta dots) for the nigh-ground states of H 1 (after dissipative evolution DE1) in a Bloch sphere. The corresponding theoretical states (corresponding numbers from 1 to 9) are represented as |zxx , 3 The left-and right-hand columns in each pane represent the real (Re) and imaginary (Im) parts of the corresponding density matrices. j. The fidelities of the four-mode output states. The numbers from 1 to 9 represent the cases from a to i.
Prepared initial states Fixed fermion parity (|0 3f or |1 3f ) Any superposition of |0 3s and |1 3s Flip error where ω represents an experimental parameter and c j (c † j ) represents the fermionic annihilation (creation) operator. If we introduce the Majorana operators, (ω is set to 1 in the main text). It can be seen that the Majorana fermions γ 2j and γ 2j+1 are paired by interactions, but γ 1 and γ 2L , which correspond to two Majoranas, remain free. The basis of the two-fold-degenerate ground-state space of this system can be defined as the eigenstates of the operatorc The parameter λ adiabatically changes from 0 to 1, the Hamiltonian will change from with Majorana mode A located at the site 2m − 1 to the Hamiltonian with Majorana mode A located at the site 2m + 1. As a result, the Majorana mode A moves a step from left to right. With K such processes, Majorana mode A moves from the left boundary to the site 2K − 1.
As shown in Supplementary Figure 1, the box in the right denotes a small system to demonstrate the braiding. A four-fermion system has been proposed with a superconducting wire [2]. After that, more theoretical schemes for an implementation or simulation of the Kitaev model are proposed [3][4][5]. Here, without considering a concrete physical realization, we theoretically proposed that the three-fermion system is also possible to complete the braiding (details of the braiding processes in the box will be introduced in Supplementary Note 5).
When the braiding process in the box is finished, Majorana mode B will be located at 2K − 1 and Majorana mode A is located at 2L. Finally, Majorana mode B at site 2K − 1 can be driven back to site 1 to complete the entire braiding process by a set of the inverse adiabatic evolutions introduced before.
It has been shown that the braiding transformations will introduce different geometric phases in the basis of the ground-space of the initial system. Measuring the geometric phases in this basis is a key part in demonstrating the statistical property of the MZMs [6]. Here, we propose to measure the geometric phases by a set of projective measurements P 1 , P 2 , · · · , P n [7]. The projective measurement P i is uniquely determined by the Hamiltonian H i as e −H i t with a large time t. The operator e −H i t preserves the fermion parity of the initial state which guarantees the adiabatic condition during the braiding. The projective measurements can be realized by an imaginary time evolution algorithm, such as, the cooling algorithm in [8].
Unfortunately, in the one-chain fermionic system, the braiding effect (geometric phase as an overall phase) in this process can not be directly measured, because of the fixed fermion parity. However, through the Jordan-Wigner (JW) transformation, a general KCM can be transformed into a 1D transverse-field Ising model (TFIM) [9, 10] whose ground states can be prepared in any superposition state and the relative geometric phases can be finally measured. The method to probe the geometric phases in the KCM is presented in Supplementary Figure 2.

SUPPLEMENTARY NOTE 2. PRINCIPLE OF THE IMAGINARY TIME EVOLUTION
Here we explain in detail the basic idea of the imaginary time evolution (ITE) algorithm based on the cooling algorithm [8]. Each ITE operator can be realized by this method. A ITE operator can be successfully realized by a set of basic steps with some probability which is dependent on the character of the Hamiltonian and the initial state. Fortunately, in the braiding process for the KCM, the ITE operator can be efficiently realized.
The circuit for one step of the imaginary time evolution is shown in Supplementary Figure   3. The key components of the quantum circuit consists of four gates [8]: (a) two Hadamard applying to the ancilla qubit (|0 and |1 are the corresponding two levels) at the beginning and at the end of the quantum circuit; (b) a local phase gate, where the parameter α is chosen to optimize the efficiency of the algorithm; (c) a two-qubit controlled unitary operation, where 1 1 is the identity operator, and is the real time evolution operator for the system. H s is the Hamiltonian of the considered system (L-site KCM here), and t is the time of evolution which is another parameter we can use to optimize the efficiency of the algorithm. For a many-body Hamiltonian, we can well approximate the unitary evolution operator U(t) by the product of a set of local unitary operators through the Trotter-Suzuki expansion [11].
For any given initial state of the system, |ψ in = Ns k=1 √ p k n k l=1 β k,l |e k,l , where N s is the number of the eigenvector subspace of Hamiltonian H s , n k ( Ns k=1 n k = 2 L ) is the degeneracy of the k-th eigenvector subspace of the Hamiltonian H s with eigenvalue E k , and |e k,l is the l-th eigenvector in this subspace. The probability to find the state in the k-th eigenvector subspace is denoted by p k . For each k, n k l=1 β k,l |e k,l is normalized, i.e. n k l=1 |β k,l | 2 = 1. Generally, we do not need to know the exact form of |e k,l and the phase between different eigenvector subspaces are not important in the algorithm.
The quantum circuit then produces the following output state: A measurement on the ancilla qubit in the computational basis {|0 , |1 } yields the states and respectively, where A 0 (A 1 ) is the normalization factor. It is clear that the coefficient of the eigenvector subspace is modified by the factor on the results of the ancilla qubit. As a result, the weight of the eigenvector subspace, especially the weight of the ground-state subspace, will be modified. To clarify this point, let us consider the module of the factor (1 ± ie −i(E k t−α) ), If the parameters α and t are properly chosen, such that, for all the eigenvalue E k , the function sin(E k t − α) will increase with the energy. Therefore, As a result, the weight of ground state of the system will increase along with the measurement result |0 on the ancilla qubit, and the weight of the ground state will decrease along with the the measurement result |1 on the ancilla qubit.
To make the parameters satisfy the condition introduced before, we normalize the Hamil- By setting α = 0 and supposing the parameter t satisfies the condition: E k t 1, the weight of the eigenvectors of the Hamiltonian change as which is the imaginary time evolution operator for small time t.
To realize a ITE with a large t, we divide it into k steps, where each step satisfies the condition introduced before. Thus, a total of k steps of ITE introduced before are applied on the initial state (we make no measurement during the cooling; we measure the qubits at the end of the manipulation). The final state we obtain is where |ϕ ji is a properly normalized state of the ancilla. a i = 1 2 sin(E i t − α) and α = 0. After the manipulations, we make measurement on the ancilla. The probability to get j |0 (the number of the success cooling manipulation) is, This is a mixture of several binomial distributions with the first one corresponding to the ground state. For a standard binomial distribution: C j k p j (1 − p) k−j , the expected value is kp and the variance is np(1 − p). Thus, the concentration interval of the binomial distribution In order to prepare the ground state of the system, the intersection between different binomial distribution should be very small. In other words which is equivalent to √ k(a 2 − a 1 ) It should be noted that the other binomial distribution corresponding to higher energy have much less intersection with the distribution corresponding to ground states. Under the condition E k t 1 (k = 1, 2, · · · , 2 L ), we have and the intersection condition can be simplify to where ∆ is the gap of the system. If this condition is satisfied, the binomial distribution of the ground state is sufficiently separated from the others, and the number of |0 outcomes during the measurement on the ancilla will be concentrated at k(1 − 1 2 a 1 ) with probability p 1 . Thus, the number of measurements that successfully obtain the ground state of the system scale as 1 The gap of the Hamiltonian during the braiding decreases polynomially with L. In addition, the overlap between the ground state of H MF i and H MF i+1 is independent on L. Therefore, the ITE operator e −H i t with large t is polynomial scaled with L in the KCM Braiding situation.

SUPPLEMENTARY NOTE 3. THE RELATION BETWEEN KITAEV CHAIN MODEL AND TRANSVERSE-FIELD ISING MODEL
The Jordan-Wigner transformation gives a one-to-one mapping between the Kitaev chain model (KCM) and a subspace of the transverse-field Ising model (TFIM) which consists of all the operators with Z 2 symmetry. As mentioned in the main text, the geometric phase due to the braiding in the KCM can be directly mapped into this subspace. In addition, the information encoded on the ground space of the TFIM is robust against the local noise in this subspace.
The Z 2 -invariant Ising model has two-degenerate ground states. There are two zero modes which are topologically protected by the boundary condition [9,10]. As an example, the without Z direction magnetic field has two basis in its ground state space: | →→ · · · → and | ←← · · · ← . If a local small Z direction magnetic field is introduced, the local flip of the basis of σ x is occurred. These two basis can be connected in the following manner: | →→ · · · → | → · · · →← | → · · · →←← · · · | →→← · · · ← | →← · · · ←← | ←← · · · ← . Thus the transformation between these two basis states can be induced and the tunneling rate is about e −L/ζ where L is the length of the chain and ζ is the correlation length. When the length of the chain is very large, the degeneracy of the spin will be topologically protected.
Although the KCM and TFIM share the same mathematical structure they are governed by different physics. The major differences between them include: the initial state in the TFIM can be prepared in the superposition of |0 Ls and |1 Ls ; the degeneracy in the TFIM can be lifted by noises without Z 2 symmetry (such as j σ x j ) and the excitations in the TFIM that corresponds to the MZMs in the KCM are non-local. The relation between KCM and TFIM is given schematically in Supplementary Figure 4. We further summarize the similarities and differences between them in Supplementary Table 1.
With the well controlled spin model, the whole system can always be in the degenerate subspace. Generally, the local noise can induce the tunneling between the MZMs whose amplitude is exponentially decaying with the distance. In our simulated spin system, the distance is very short. The tunneling between the MZMs will be induced if all the local noises are present. Fortunately, the noises in the optical system can be well controlled and the tunneling can be well suppressed by reducing some local noises. The robustness of the information encoded in the degenerate subspace against single local noises have been demonstrated in our experiment.

GEOMETRIC PHASE
Generally, the angle (ϕ B ) obtained from the Bargmann invariant [12] associated with the three states ρ 1 , ρ 2 and ρ 3 in the ray space is arg(tr(ρ 1 ρ 2 ρ 3 )) [13][14][15]. ϕ B is gauge invariant and is equal to the negative geometric phase of ϕ g which is associated with the geodesic triangle of the three states [13]. Actually, following the line of Pancharatnam [16], it can be extended to the situation with many states.
In KCM system, the geometric phase, associated with the basis |i gf of the ground-state space of H MF 0 , obtained during the braiding can be determined by [7] where ϕ g0 (ϕ g1 ) is the geometric phase of the basis |0 (|1 ) in the ground space of where γ i represents the corresponding Majorana operator. The geometric phases become It can be expressed in the spin model as where |0 3s (|1 3s ) is the mapping of |0 3f (|1 3f ) under JW transformation, and it is the basis of the ground space of H 0 . All the physical quantities have been transformed into a spin system.
To determine the geometric phases, we prepare the initial state in |φ 0 = α|0 3s + β|1 3s with α and β representing the complex numbers (not normalized). We then apply the operator U ex = e −H 2 t e −H 1 t to this state. The final state is determined and the relative geometric phase are obtained. If we prepare the initial state as |0 3s or |1 3s , the state after the operator U ex is still the same which can be used to determine the off-diagonal elements of the braiding matrix [7]. In this experiment, we perform a tomography to measure the effect of the operator U ex in the ground space of H 0 , and the information of the braiding matrix is thus determined [7]. A benefit of this method is that we can also prepare the initial state by the imaginary evolution method which needs only the Hamiltonian H 0 .

SUPPLEMENTARY NOTE 5. THE ADIABATIC EXCHANGE PROCESS IN THE SIMPLEST KITAEV CHAIN MODEL
As mentioned in the main text and above, the braiding of two MZMs can be completed in the three-fermion KCM [10]. Now, we explicitly illustrate the adiabatic braiding process in this system (see Figure 1 in the main text). Throughout the entire braiding process, the Hamiltonian can be described in the general form where c + i and c i represent the creation and annihilation fermionic operators at position i. t ij represents the interaction amplification between the particles in positions i and j. µ 1 (≥ 0) represents the chemical potential at site 1.
The Majorana operators can be introduced as follows: The inverse relations between the Fermi operators and the Majorana operators are The general Hamiltonian can be rewritten in terms of Majorana operators as (overall constants are ignored) The initial Hamiltonian is defined by setting t 13 = 0, t 12 = t 23 = 1 and µ 1 = 0 as follows: The ground state of the Hamiltonian H MF 0 is two-fold degenerate. The basis denoted by |0 3f and |1 3f consist of the eigenstates of the operatorc † fc f (c f = 1 2 (γ 1 + iγ 6 )). Explicitly, |Vac where |Vac denotes the vacuum of c i s [2,17]. The Majorana zero modes (MZMs) γ 1 and γ 6 are located at the endpoints of the chain. Now we will show that the two isolated MZMs denoted by A and B, located at sites 1 and 6, in the KCM can be driven to complete the braiding transformation. The process is accomplished by adiabatically tuning the parameters in equation (28), where different geometric phases will be added to the states |0 3f and |1 3f . The entire braiding process can be divided into the following three components: (1) Let the parameters in the Hamiltonians described by equations (28) be The Hamiltonian then becomes When the parameter λ adiabatically changes from 0 to 1, the Hamiltonian will be trans- At the same time, the Majorana mode A will move from site 1 to site 3.
The Majorana mode B, initially located at site 6, will be driven to site 1. This adiabatic process is the key element of the braiding.
(3) Let The Hamiltonians described by equations (28) becomes When λ adiabatically changes from 0 to 1, the Hamiltonian will be transformed from H MF 2 back to H MF 0 , and the Majorana mode A will be driven from site 3 to site 6. The braiding of the MZMs A and B is thus completed. The basis of the ground-state space will develop various geometric phases, which correspond to the non-Abelian characteristics of MZMs [2,17].
Currently, the direct implementation of braiding in fermionic systems remains a considerable challenge. In addition, the geometric phases obtained in one chain model can not be directly measured in fermionic system. By employing the Jordan-Wigner (JW) transformation, we can transform the fermionic Hamiltonian into a spin-1/2 system. In the three-fermionic KCM, the JW transformation can be written as follows: where σ + i (= σ x i + iσ y i ) and σ − i (= σ x i − iσ y i ) represent the raising and lowering operators of the spin at site i. After the JW transformation, a general KCM will be transformed into a transverse-field Ising model (TFIM). The ground-state degeneracy in the ferromagnetic phase of the TFIM correspond to that of MZMs at both ends of the chain [10]. However, because of the non-local nature of the JW transformation, the MZMs in the KCM correspond to nonlocal excitations in the spin model, that is, The braiding of non-local excitations is not well defined. However, the geometric phases and It is noted that all the interactions in the KCM are local (three-body next-nearestneighbor interaction at most). The ground state of the Hamiltonian H 0 is also two-fold degenerate and the basis in the degenerate space is denoted by |0 3s and |1 3s . Explicitly, |0 3s = 1 2 (| ↑↓↓ + | ↓↑↓ + | ↓↓↑ + | ↑↑↑ ) and |1 3s = 1 2 (| ↓↓↓ + | ↑↑↓ + | ↑↓↑ + | ↓↑↑), where | ↓ and | ↑ represent spin up and spin down, respectively.

STATES ENCODED IN OPTICAL SPATIAL MODES
For a given Hamiltonian H with a complete set of eigenstates |e k and the corresponding eigenvalues E k , any arbitrary pure state |φ can be expressed as with q k representing the corresponding complex amplitude. Here, we focus on pure states, but the argument is also valid for mixed states. The corresponding imaginary-time evolution (ITE) operator (U ) [18] on the state becomes The evolution time t is chosen to be 5, which is long enough to drive the input state to the ground state of H in our analysis. After the ITE, the amplitude q k is changed to be q k exp(−5E k ). We can see that the decay of the amplitude is strongly (exponentially) dependent on the energy: the higher energy, the faster the decay of the amplitude. Therefore, only the ground states (with lowest energy) survives this evolution. Furthermore, due to the fact that the Hamiltonian H 0 , H 1 and H 2 can be diveded into two commuted parts, we can separate the ITE operator to two ITE operators whose eigenvectors are very easy to obtain (see equation (2) in the main text).
We then show the detailed processes to probe the geometric phases and braiding matrix, and demonstrate the local noises immunity of MZMs.

Geometric phases and braiding matrix
(1) In our experiment, the state information is encoded in the optical spatial modes. The basis of a two-level system can be expressed in the eigenstates of σ z (denoted by |z , with an eigenvalue of 1, and |z , with an eigenvalue of -1), σ y (denoted by |y , with an eigenvalue of 1, and |ȳ , with an eigenvalue of -1) or σ x (denoted by |x , with an eigenvalue of 1, and |x , with an eigenvalue of -1). The eigenstates of H 0 (= −(σ x 1 σ x 2 + σ x 2 σ x 3 )) can be expressed as {|xxx , |xxx , |xxx , |xxx , |xxx , |xxx , |xxx , |xxx }. For the term of −σ x 1 σ x 2 , the largest factor e 5 is obtained when the products of eigenvalues for particles 1 and 2 are equal to 1, according to the ITE (equation (42)). For the other cases the products of eigenvalues for particles 1 and 2 are equal to −1. Hence, the added amplitude factor becomes e −5 , which is negligible compared to e 5 (for any state with nonzero ground state probability, we can always increase its amplitude to arbitrary high level by increasing the efficiency of the ITE, i.e., by increasing the evolution time from 5 to ∞). As a result, only the terms involving the basis of {|xxx , |xxx , |xxx , |xxx } are preserved after the projection of e −σ x 1 σ x 2 * 5 . Similarly, after the projector e −σ x 2 σ x 3 * 5 , the largest factor e 5 is obtained when the products of eigenvalues for particles 2 and 3 are equal to 1 and the preserved state would ivolve in the basis of {|xxx , |xxx }. As a result, the ITE for H 0 can be written as which is normalized with complex coefficients α and β (|α| 2 + |β| 2 = 1). For simplicity, we would omit below the normalization of the output states.
Robustness against perturbation operation of 1 4 (iσ y 1 σ x 2 + σ y 1 σ y 2 + σ x 1 σ x 2 − iσ x 1 σ y 2 ) (1) For any input state with nonzero amplitude on the ground states of H 0 , when it is subjected to the ITE U 0 , the output state is the ground state of H 0 with the form of (2) After the perturbation operation of 1 4 (iσ y 1 σ x 2 +σ y 1 σ y 2 +σ x 1 σ x 2 −iσ x 1 σ y 2 ), the state becomes The simulation of the perturbation operation is realized by the rotation of the half-wave plates in the corresponding spatial modes.
(3) The same ITE U 0 is applied to the state |φ 0 , where only the terms |xxx and |xxx are preserved. As a result, the final state becomes which is the same as the initial state and the protection of local site flip is shown.
Robustness against perturbation operation of (σ z + 1)/2 (1) For any input state with nonzero amplitude of ground states of H 0 , when it is subjected to the ITE U 0 , the output state is the ground state of H 0 with the form of |φ 0 = α|xxx + β|xxx .
(2) The perturbation operation 1 2 (σ z 1 +1) means that there is a probability of 0.5 to rotate the particle 1 along σ z direction and a probability of 0.5 to do nothing. As a result, the state would becomes where the subscripts np and pe represent the cases without and with perturbation, respectively. The setup is similar to the case of transforming the basis of particle 1 from σ x to σ z , as expressed in equation (45). However, the four output modes are denoted by {|xxx , |xxx , |xxx , |xxx } for the perturbation operation and {|zxx , |zxx , |zxx , |zxx } for the basis rotation.
(3) The same ITE U 0 is applied to the state |φ 0 , where only the terms |xxx and |xxx are preserved. As a result, the final state becomes which is the same as the initial state and the property of local phase-error immunity is shown.
To illustrate the encoding process more clearly, the cross sections of the output spatial modes of each process are presented in Supplementary Figure  ference between different spatial modes [19]. To achieve high visibility of the interference between different spatial modes, the gradients of the beam displacers are carefully adjusted, and phase compensation using additional wave plates (not shown in Figure 2 in the main text or Supplementary Figures 12 and 13) between different spatial modes is applied. In our experiment, the fidelity for the investigation of the geometric phases is approximately 94.13 ± 0.04%, and the fidelities for the investigation of flip-error immunity and phase-error immunity are approximately 97.91 ± 0.03% and 96.99 ± 0.04%, respectively.

SUPPLEMENTARY NOTE 10. EXPERIMENTAL SETUP FOR SIMULATING THE EXCHANGE OF MZMS
As shown in the Pre pane in Figure 2 in the main text, the polarization of a single photon is rotated using a half-wave plate (HWP), and the photon is then sent to a BD30 with a vertical beam displacement of 3.0 mm. Both output beams are further rotated using HWPs and sent to another BD30 with a horizontal beam displacement of 3.0 mm. There are now four output beams, which are further rotated using HWPs and pass through a BD60 with a horizontal beam displacement of 6.0 mm. We then obtain eight output beams, and the distances between neighboring beams are all equal to each other, which facilitates the measurement of interferences between different optical modes. The relative amplitudes of the eight beams are controlled using HWPs in the corresponding modes. The eight optical beams represent the basis of the eight-dimensional space of the systems, which are initially expressed in the basis of σ x as {|xxx , |xxx , |xxx , |xxx , |xxx , |xxx , |xxx , |xxx }.
The precision of the dissipative evolution is theoretically dependent on the scale of the evolution time t. In our protocol using a polarization beam splitter (PBS), it is dependent on the ratio between the reflected and transmitted parts of the vertical polarization after the PBS, which can be higher than 500:1. As a result, the optical modes with vertical polarization are discarded directly after each dissipative evolution process (the reflection of photons with vertical polarizations is not shown in Figure 2 in the main text). For the first dissipative evoultion DE0 (the projection of exp(σ x 2 σ x 3 * t) exp(σ x 1 σ x 2 * t)), the normalized output state would have the form of |φ 0 = α|xxx + β|xxx . To perform DE1, we must express the system located at site 1 in the basis of σ z , which is implemented by BR1, which contains two HWPs and a BD60. Consider the dissipative evolution of H 2 : two HWPs and a BD30 in BR2 are used to transform the state of the system located at site 1 back to σ x . The following setup contains HWPs and a BD60 to transform the state of the system located at site 2 to σ z . To transform the state of the system located at site 3 to σ y , two BD30s and two HWPs with a quarter-wave plate (QWP) are used. An inverse base rotation is implemented in BR3 to transform the entire state back to σ x to implement DE0 once again.
For the case of two-mode measurement (TM), the two spatial modes are denoted by H and V , respectively. The amplitudes of the two modes, which correspond to the amplitudes of the two polarizations, can be measured directly. For the measurement of R and D, the interference between these two modes is needed. By applying a HWP in the mode denoted by V , the spatial modes can be combined again using a BD30. The interference of the optical modes is transformed to the interference of polarizations, which can be measured using the standard polarization-analysis setup. For the four-mode measurement (FM) case, the four modes are denoted by HH, V H, HV and V V . Measurements of interference between the optical modes are carried out using a step-by-step process. For example, the measurement of RR is obtained by first combining HH and V H to obtain RH and then combining HV and V V to obtain RV . RH and RV are then further combined to obtain RR. As a result, the requirement of the 16 measurement configurations for the standard two-qubit-state tomography can be achieved.
To fully characterize the behavior of the entire setup, we experimentally performed full quantum process tomography [20,21]. By expanding the output state E(ρ) with a complete set of basisÊ m of the Pauli operators {I(identity), σ x (X), σ y (Y ), σ z (Z)}, the operation of the quantum process can be expressed as E(ρ) = mn χ mnÊm ρÊ † n . The physical process E is completely and uniquely characterized by the 4-by-4 matrix χ [22]. For simulating the exchange of MZMs, the final state can be written as 1 √ 2 (α−iβ)|xxx + 1 √ 2 (β −iα)|xxx when the initial state is α|xxx +β|xxx . The experimental χ-matrix for demonstrating non-trivial statistics is denoted by χ e n , while χ t n represents the theoretical prediction. The fidelity of the experimental result χ e n is approximately 94.13 ± 0.04%, which can be calculated from (Tr √ χ e n χ t n √ χ e n ) 2 [22]. The experimental χ-matrix for demonstrating the immunity against flip errors is denoted by χ e f . χ t f represents the theoretical prediction, which corresponds to