Simulation of anyonic statistics and its topological path independence using a 7-qubit quantum simulator

Anyons, quasiparticles living in two-dimensional spaces with exotic exchange statistics, can serve as the fundamental units for fault-tolerant quantum computation. However, experimentally demonstrating anyonic statistics is a challenge due to the technical limitations of current experimental platforms. Here, we take a state perpetration approach to mimic anyons in the Kitaev lattice model using a 7-qubit nuclear magnetic resonance quantum simulator. Anyons are created by dynamically preparing the ground and excited states of the 7-qubit Kitaev lattice model, and are subsequently braided along two distinct, but topologically equivalent, paths. We observe that the phase acquired by the anyons is independent of the path, and coincides with the ideal theoretical predictions when decoherence and implementation errors are taken into account. As the first demonstration of the topological path independence of anyons, our experiment helps to study and exploit the anyonic properties towards the goal of building a topological quantum computer.


Introduction
It is a fundamental question to investigate the physical properties of exchanging two identical particles. In onedimensional space, the exchange is trivial as the two particles inevitably collide. In three-dimensional spaces, a wave function of a system acquires either +1 or −1 phase factor for bosons and fermions upon the exchange, respectively. Alternatively, in two-dimensional spaces, exotic statistical properties emerge when exchanging two identical particles, leading to the theoretical existence of anyons. For instance, the wave function can acquire an arbitrary phase factor q e i , ranging continuously from +1 to −1 for Abelian anyons, or undergo non-trivial unitary evolutions for non-Abelian anyons. This exotic feature of anyons called fractional statistics has attracted great interest over the past few decades [1][2][3][4][5][6].
As truly two-dimensional systems do not exist in nature, anyons appear in an effective two-dimensional system as quasi-particles, a collective behaviour of a group of fundamental particles behaving as a single entity. The experimental evidence of quasi-particle anyons was first discovered in the fractional quantum hall effect in the 1980s [2]. Since then, an experimental quest for anyons has interested researchers not only for their fundamentally intriguing feature, but also for their application in performing protected quantum information processing (QIP). As the goal of QIP is to exploit quantum mechanical properties for computation, manipulating quantum properties in a precise manner is critical [7]. Topological properties of anyons have potential to achieve such a goal, as these properties are resilient to small fluctuations. Therefore, the prospect of utilizing anyons' fractional statistics to achieve robust control has gained much attention. In the past two decades, many powerful quantum computing schemes using anyons have been proposed. Such schemes constitute topological quantum computing (TQC) [5,[8][9][10] and topological quantum error correction [11,12]. Of the many proposals, the toric code of spin lattice [12] is one of the most renowned. By artificially designing a 2. Theoretical model 2.1. Toric code model In this model, qubits are located on the edges of a two-dimensional lattice as shown in figure 1(a). The Hamiltonian of the system consists of two different types of four neighbouring-qubit interactions, XXXX and ZZZZ (X and Z are Pauli matrices) at a vertex v and plaquette p, respectively. Hence Applying an operator X i changes the eigenvalues of two of the B p operators nearby the ith qubit, thus, shifting the system to a higher energy. Similarly, applying Z i changes the eigenvalues of two of the A v operators nearby the ith qubit. The particular vertices or plaquettes that have been excited by the ith qubit are marked by red and blue dots respectively, and can be identified as having defects. These defects can be represented as localized quasi-particle anyons. The anyons on the vertices and plaquettes are called e and m particles, respectively. The excited qubit which creates the pair of 'defects' can be identified as a bisector of the string connecting the pair.
where ( ) , which are referred to as stabilizer operators. Here, star (v) is a set of four spins that share a link with the vertex v, and bond (p) is a set of four spins placed at the edges of the plaquette p. Since all the stabilizer operators commute with each other, the ground state |y ñ g of this Hamiltonian is a+1 eigenstate of the A v and B p operators (note the minus sign in the Hamiltonian in equation (1)): The case with a periodic boundary condition on the lattice exemplifies a toric code [12] where the ground states are four-fold degenerate. The degenerate ground states form a protected subspace from possible noise-induced excited states.
Excited states of this Hamiltonian are created by applying single-qubit operators X and/or Z to the ground state. These operations create two types of quasiparticles, e particle at a vertex or m particles on a plaquette, as described in figure 1(b). Subsequently, as shown in figures 2(a) and (b), one can move a m particle to a different plaquette by applying X to a qubit nearby. Particles created on the same site annihilate each other so the X operation effectively moves the m particle. Analogously, applying Z to a relevant nearby qubit moves an e particle.
One can demonstrate anyonic statistics between e and m particles by moving one around the other, making a closed loop as shown in figure 2(c). This braiding operation is equivalent to the two successive particle exchanges. Note that it is not possible to exchange their positions once, since one is located at a vertex and the other at a plaquette. It can be shown that the wave function acquires a −1 phase factor (corresponding to a π phase) after such braiding, indicating that a single exchange of Abelian anyon e and m particles would result in a p 2 phase. Therefore, the anyonic statistics of e and m particles is p 2.

The seven-qubit toric code
The planar seven-qubit model used to demonstrate the path independent property of anyonic braiding is shown in figure 3(a). For the case of a periodic lattice, the Hamiltonian consists of the four body interactions (equation (1)). However, for the seven-qubit model, we consider a lattice with a rough boudary, which results in two-body ZZ interactions at the boundary. Therefore, the Hamiltonian H 7 of this system is Moreover, due to the absence of periodic boundary conditions, the ground state of this model |y ñ g 7 is nondegenerate. Since the state | ñ Ä 0 7 in the computational basis is already a+1 eigenstate of  B 1, , 5 , the ground state is given by projecting | ñ Ä 0 7 on to the +1 eigenstate of A 1 and A 2 : Due to the lattice structure, exciting qubit 1 with the operator Z creates a single e particle at v 1 rather than creating a pair, whereas exciting qubit 5 with the operator X still creates a pair of m particles at the plaquettes associated with B 3 and B 5 . Refer to figure 3(b) for the particle locations. Starting with this excited state, there are (a) Applying X 1 and X 2 flips the eigenvalues of their shared B p . Hence, creating two defects on the same plquette annihilates the defects. (b) Net effect of annihilating the m particles while creating another is moving the particle. (c) Full operation which braids a m particle around e is X X X X 0 3 2 1 . After this braiding operation, the wave function gains a π phase. Since such braiding corresponds to exchanging the particles twice, this π phase demonstrates that the anyonic statistic of e and m particles is p 2. Note that it is not possible to exchange the two particles' positions once, since they live in different places (e on vertex and m on plquette). three possible loops to braid m particles as shown in figures 3(b) and (c): the trivial braiding operation = l X X X X 0 4 5 6 7 where a m particle braids around v 2 without an e particle, and the two non-trivial braiding operations = l X X X X X X 1 1 2 3 5 6 7 and = l X X X X 2 1 2 3 4 where a m particle braids around v 1 with an e particle. The wave function remains the same when the braiding operation is trivial; however, if the operation is non-trivial, the wave function picks up a π phase from the fractional statistics.
To experimentally demonstrate the path independence of anyonic braiding, we simulated anyonic physics manifested in the seven-qubit toric code using a liquid-state NMR quantum simulator. Since it is experimentally challenging to engineer the toric code Hamiltonian which involves four-body interactions, we took the state preparation approach: dynamically preparing the ground and excited states of the toric code Hamiltonian in a NMR system, instead of generating the toric code Hamiltonian and cooling the system.
The quantum circuit which simulates the anyonic physics is shown on the right in figure 3. The main idea is to prepare |y ñ g 7 and then create a superposition (| | ) y y ñ + ñ 2 mm mme , where |y ñ mm is a state with the pair of m particles, and |y ñ mme is a state with both the e particle and pair of m particles. If braided along the non-trivial paths such as l 1 or l 2 in which e circulates around |y ñ m, mme gains a π phase due to the fractional statistics; otherwise, |y ñ mme remains unchanged. By measuring the variation of the relative phase on |y ñ mme before and after the braiding, one can deduce whether the braiding path is trivial or not and,furthermore, demonstrate the path independence. The details are described as follows.
First, two Hadamard and six controlled-NOT (CNOT) gates are applied to prepare the ground state |y ñ g 7 of the seven-qubit toric code Hamiltonian from |  ñ 00 0 , as depicted in figure 3. Then, applying X 5 and Z X 1 5 on |y ñ g 7 generates , 4 mm g mme g 5 1 5 7 7 respectively. To create a superposition of the two, we apply Z X 1 5 since . When the anyons are braided along a non-trivial loop, the superposition picks up a relative phase on the |y ñ mme component. Finally, anyons are annihilated by reversing the creation operator Z X 1 5 in order to measure this relative phase. The system ultimately evolves to either the ground state |y ñ g 7 or the excited state |y ñ e 7 depending on different braiding paths. Therefore, we can experimentally demonstrate the path independence nature of anyonic braiding if the two phases obtained under the two different non-trivial loops l 1 and l 2 are the same.
The states corresponding to each step of the circuit shown on the right in figure 3 are where q a is the phase gained from the anyonic statistics for different loops l 0,1,2 . When the m moves around the trivial loop l 0 , the final state |y ñ d ends up at the ground state |y ñ g 7 (q = 0 a ), whereas when the m is moved around the non-trivial loops l 1 or |y ñ l , d 2 ends up at the excited state |y ñ e 7 (q p = a ). In order to demonstrate path independence of anyonic braiding experimentally, we need to implement the entire circuit and observe q a for different loops.

Experimental implementation in NMR
Our seven-qubit NMR processor is the per-13 C-labelled dichlorocyclobutanone derivative [32,33] dissolved in d6-acetone. The molecule consists of seven 13 C spins and the five 1 H spins. We denoted the seven nuclear spins of 13 C as qubits, while 1 H nuclei were decoupled throughout all experiments except for the initialization step to boost polarization on 13 C. The molecular structure is depicted in figure 4(a), where two nearest-neighbouring 13 Cs have stronger coupling strengths, implying the ability to implement a faster two-qubit gate. Therefore, by comparing the geometry of toric code qubits and the structure of nuclear spins, we mapped each toric code qubit to the nuclear spin in as shown figure 4(b). The natural Hamiltonian of this system is described as where n i is the chemical shift frequency of the ith spin, and J ij is the coupling strength between the ith and jth spins (refer to appendix A for values of the parameters). All experiments were conducted on a Bruker DRX 700 MHz spectrometer at room temperature. The experiment was divided into five main steps as shown in figure 4 (b), as follows: PPS initialization. We first utilized the cat-state method proposed in [31] to initialize the system to a labelled pseudo-pure state (PPS) state. It can be represented by a trace-zero deviation matrix of the form | | r = ñá Z 000000 000000 PPS C 7 , where C 7 is the labelling qubit, and Z C 7 describes the state of C 7 is s z . Two techniques were adopted before this initialization step to improve the signal-to-noise ratio (SNR). One is turning on 13 C and 1 H couplings temporarily at the very beginning and applying a SWAP gate between C 7 and H 5 , to achieve a ∼4 times higher polarization on C 7 . The other one is performing RF-selection [34] sequence to pick out a slice of the NMR sample which experiences much better radio-frequency (RF) homogeneity by randomizing the other part with worse RF homogeneity. Subsequently, the labelled PPS state was prepared using non-unitary transformations via gradient fields and phase cycles [31]. The total length of the initialization sequence is about 100 ms. Refer to appendix A for more details about the PPS initialization step.
Ground state preparation. Unlike the theoretical circuit on the right in figure 3, the implemented circuit in NMR prepared the ground state from˜| | r = ñá Z 000000 000000 PPS C 7 , rather than the required pure state | ñ 0000000 . Nevertheless, sincer PPS contains half of | ñ 0000000 C 7 and half of | ñ 0000001 C 7 , we can simply write the deviation matrix after the ground state preparation as: where |ỹ ñ g 7 results from | ñ 0000001 C 7 . Under perfect unitary transformation, | |ỹ y ñá g g 7 7 stays orthogonal to | | y y ñá g g 7 7 throughout the circuit, thus not interfering with the final result if it can be separated in the NMR spectra. In fact, the additional two CNOT gates in the beginning of the ground state preparation shown in figure 4(b) were specifically added to achieve this separation. However, in the presence of errors, | |ỹ y ñá g g 7 7 did slightly modify with the final result, as analysed in section 4.
The entire ground state preparation step was optimized by a 60 ms gradient ascent pulse engineering (GRAPE) pulse [35] based on a subspace approach [36]. The simulated fidelity of this pulse is over 0.99. Additionally, a special rectification method was used in the experiment to ensure that all of the GRAPE pulses acting on the spins were very close to theoretical expectations [37,38]. We performed modified stabilizer measurements after the ground state preparation step to verify the state. This step is explained in detail in appendix C.
Anyon creation, braiding and annihilation. These three parts shown in the emulation circuit (on the right of figure 3) are compressed together to simplify the circuit as they only involve single-qubit rotations. The trivial loop l 0 and non-trivial loops l l , 1 2 are all depicted in figure 4(b), and in each experiment only one loop was implemented. The three braiding operators were realized by 1 ms GRAPE pulses, respectively. In principle, after this stage we can determine q a by measuring coefficients of |y ñ g 7 and |y ñ e 7 in equation (8), but it does require many measurements in a seven-qubit system.
Measurement. This additional 'measurement' step is added to estimate q a with a few measurements, which allows us to measure diagonal elements of the final state and then extract the value of q a . It separates diagonal elements of |y ñ g 7 and |y ñ e 7 via basis transformation by evolving the state |y ñ d to . In this case  figure 3, respectively. The circuit contains five steps with the sequence length at the bottom of each step: the labelled PPS state preparation based on the cat-state method [31], and the detailed network can be found in the appendix A; preparation of the ground state |y ñ g 7 of the toric code Hamiltonian; anyonic manipulation including creation, braiding, and annihilation (all three braiding paths are shown here, but in each experiment we just implement one); measurement circuit which converts the state tomography to simpler diagonal elements measurement; readout pulse on C 7 to measure the required diagonal elements.
To evaluate q a , we estimated | | a 2 by measuring the diagonal elements of | | y y ñá p0 p0 and similarly | | b 2 by measuring the diagonal elements of | | y y ñá p1 p1 . Diagonal elements readout. Since the diagonal elements cannot be directly observed in NMR, we indirectly measured them by applying the readout pulse which rotates C 7 by π/2 around the y-axis. This readout pulse generated single coherences from the diagonal elements, and thus a detectable signal with distinct frequencies depending on the state of the other qubits (see appendix for detailed descriptions). In particular, the transitions relevant to | | a 2 and | | b 2 estimations are at four distinct frequencies centred around n 7 (resonant frequency of C 7 ): 61.25 Hz, 24.09 Hz, 32.24 Hz, and −4.93 Hz. Therefore, the real coefficients of the peaks at these specified frequencies can indirectly estimate the diagonal elements of interest. There is one assumption in the measurement of diagonal elements in the above method. The peaks are actually generated by the subtraction of two relevant diagonal elements after rotating C 7 by π/2 around the yaxis (see appendix). For example, the intensity of the peak at 61.25 Hz corresponds to | | | | ñá -ñá 0000000 0000000 0001000 0001000 , but we only need the value of the first term. So we assume that the latter term is 0 in order to get the value of the first term. We simulated the contributions from such small elements and found that this assumption should be good enough for the accurate estimation of q a .

Result and discussion
We measured the anyonic phases of the four different cases: where noBD and BD0 ideally have q = 0 a , and BD1 and BD2 have q p = a . GD and MM refer to the ground state preparation and measurement steps, respectively. Figure 5 shows the C 7 spectra of the labelled PPS and the above four cases. The experimental spectra agree qualitatively with our theoretical predictions. First, in theory, we expect to observe the same spectra for the noBD and BD0 cases and the same spectra for the two non-trivial braiding cases (BD1 and BD2) due to the path independent nature. From figure 5, it is clear that the spectra of noBD and BD0 match well, and also that BD1 and BD2 match well. Second, our experimental spectra matched well with the simulated spectra. In theory, the spectra resulting from the four cases are expected to show four peaks with equal height (two generated from | | y y ñá g g 7 7 and the other two from | |ỹ y ñá g g 7 7 ), which is a quarter of the labelled-PPS peak. The spectra shown in figure 5 qualitatively illustrate the expected behaviours. Third, recalling equation (15), we expect to observe no peaks at c=32.24 Hz and = - and q a are displayed in table 1 for all noBD, BD0, BD1 and BD2 cases. The experimental results show that the anyonic phases under the two different non-trivial braiding paths l 1 and l 2 agree within the errors: (153.9 ± 3.8)°and (151.4 ± 3.8)°. These experimental values clearly demonstrate path independence and the phase gained under the non-trivial paths compared to the cases of the trivial and no braiding paths ((17.4 ± 6.0)°and (12.1 ± 9.5)°, respectively). However, the experimental q a have discrepancies with the theoretical values, which are  0 for the trivial and no braiding paths, and  180 for the two non-trivial paths. For the non-trivial cases, this deviation is mostly attributed to the tiny peaks at a and b (figure 5), which result in | | a » 0.05 2 , because q a is highly sensitive to | | a 2 as it is small and in the denominator (equation (15)).
For instance, consider a theoretical case when | | a 2 is 0. In this case, regardless of a value of | | b q , a 2 is always π. Similarly, for the trivial and no braiding cases, the deviation of q a is mostly caused by the tiny peaks at c and d ( figure 5), resulting in | | b » 0.02 2 rather than the theoretical value of 0.
To investigate how the unwanted small peaks arise, we numerically simulated the NMR circuit starting from the ideal labelled PPS state using 99% fidelity unitaries calculated from the GRAPE pulses in the presence of the decoherence effect. The assumptions that we used to simulate decoherence are shown in appendix. The results of the simulation indicate that the errors increase the trivial loop phases to ∼20°whereas the non-trivial loop phases decrease to ∼160°, blurring the difference between the two. It should be noted that most of the phase deviation comes from the decoherence effect; simulating only the gate imperfection from 99% fidelity unitaries results in the non-trivial phases of ∼177°. Now we discuss the different sources of error in detail.
First, the error primarily comes from the decoherence effect, and the ground state and measurement pulses contribute the most in causing the biases in the q a determination. In particular, the ground state we prepared was the ground state of the seven-qubit toric code, not a ground state of our physical system. Therefore, the ground state preparation step is susceptible to decoherence, as there is no protection of the ground state by the energy gap in our NMR system. Second, to a much lesser extent, equation (15) no longer accurately determines the anyonic phase in the presence of gate imperfections. Therefore, to estimate the anyonic phase independent of imperfections of ground state and the measurement pulses, a different equation is required. However, it is difficult to find such an Figure 5. NMR spectra of C 7 after the labelled-PPS, BD0, BD1, BD2, and noBD cases. The experimental data are shown in blue, and the red spectra are the fit of the experimental spectra produced by the least-square method. The labelled-PPS state shows a single peak at the expected frequency. In theory, the PPS peak splits into the four distinct peaks with equal heights for BD0, BD1, BD2 and noBD cases, and the experimental spectra show that indeed the expected four peaks appear for all cases. However, the peak height is less than a quarter of the PPS peak due to decoherence effect. As expected, for the BD0 and noBD cases, the peaks at a and b are more dominant than peaks at c and d. Whereas for the BD1 and BD2 cases, the peaks at c and d are more dominant than peaks at a and b.   equation that is accurate and whose variables can be easily measured. Since the braiding operation is 1 ms, whereas the ground state and measurement pulses are 60 ms, the ground state and measurement pulse imperfections contribute more significantly to the q a determination. Moreover, we expect that gate imperfections are worse in experiments than in simulation, which could explain the <10  discrepancy between the simulation and experimental values after accounting for the other sources of error. Third, we also examined the effect of | |ỹ y ñá g g 7 7 on the q a determination through numerical simulations. We simulated two scenarios with one started from | |   ñá 00 0 00 0 and the other from the labelled PPS. As mentioned above, the one starting with the labelled PPS results in the non-trivial q a of ∼160°, whereas the one started from the pure state results in ∼150°. This signifies that the contribution from | |ỹ y ñá g g 7 7 cannot be neglected completely when both gate imperfections and decoherence effects are present.

Conclusions
We have successfully demonstrated path independence of anyonic braiding statistics by braiding two anyons under two different non-trivial paths in a seven-qubit NMR quantum simulator. The anyonic phases of the two non-trivial paths l 1 and l 2 agree within the errors: ( ) •  153.9 3.8 and ( ) •  151.4 3.8 for l 1 and l 2 , respectively. As references, the cases of no braiding and braiding along a trivial path are also implemented. We measured significantly smaller phases for these trivial cases compared to the non-trivial cases, confirming the extra phase acquired by the anyons in the non-trivial cases. The deviation of the anyonic phases from the theoretical value are well accounted for by the inherent errors of decoherence and imperfect gates. These contributions can be mostly attributed to the ground state preparation and measurement steps, as these steps are significantly longer than the braiding step. Other experimental schemes or setups where such a long preparation step can be prevented may be less prone to such errors. Moreover, the measurement step which is used to remarkably reduce the number of experiments in our NMR system may be eliminated in other settings.
As a step towards the realization of TQC, we do not simulate the many-body interactions in the toric code Hamiltonian but alternatively use a state preparation approach to simulating the toric code. This method is sufficient to simulate some particular anyonic properties such as the path independent nature shown in this paper; however, realizing fault-tolerant topological quantum computation would ultimately require engineering such Hamiltonians with many-body interactions. Fortunately, quantum simulation provides exponential speedup, outperforming classical computers as well as highly controllable systems instead of the natural intractable solid state systems. Hence, quantum simulation is a promising solution for creating and engineering the full toric code Hamiltonians [39][40][41][42] in the near future, and it may shed light on the goal of building a topological quantum computer in a fault-tolerant manner.
å å where γ is the gyromagnetic ratio of the nuclear spins,  is the2 2 12 12 identity matrix, and  » -10 5 represents the polarization of the system. Typically, g = 1 C and g = 4 H with a constant factor ignored. As the large identity matrix part does not evolve under unital operators (which is roughly the case in our experiment as the experimental time is far less than T 1 ) and it cannot be measured in NMR experiments, we can simply neglect the identity part and rewrite the input state as In the following calculations we only focus on this deviation density matrix assuming that the identity has no influence on the entire experiment. a. Rotate 13 C to å = X i 1 7 C i by a 1 ms p 2 GRAPE pulse around y-axis on 13 C channel, and then crush it by a 2 ms gradient pulse. The total length is 3 ms and the state at step a is r = å = Z 4 a i 1 5 H i . b. SWAP the signal of C 7 and H 5 by applying a 8 ms GRAPE pulse. This GRAPE pulse was designed via stateto-state approach and hence not a universal SWAP gate. The reason of implementing this SWAP operation is to improve the C 7 signal by four times in principle, which enables a much better SNR in experiment. The state at . c. Turn on the Waltz-16 decoupling sequence on 1 H channel. It averages out the signals of all 1 H spins and their interactions with the 13 C spins. In quantum information, this step is equivalent to reducing the 12-qubit system to seven qubits which only involve 13 C spins. Hence, the state at step c is r = Z 4 c C 7 . Compared to the input thermal equilibrium state of r 0 , the signal of C 7 has been boosted by four times.
d. RF-selection technique is used to pick out a sub-sample which has much better RF homogeneity. As the sample in NMR has some volume in centimeters, the RF pulse applied to the sample may have inhomogeneity. Some molecules located in the centre of the RF coil experience the ideal RF amplitude, while majority of molecules experience over-rotation or less-rotation for the sake of RF amplitude inhomogeneity along the sample size. Since NMR readout is an ensemble average, the large portion with bad homogeneity contributes a lot to the final signal and causes accumulated errors when multiple pulses are implemented. RF-selection sequence [34] is such a technique to randomize this inhomogenous portion to x−y plane while keeping the homogenous portion in the thermal equilibrium state, followed by a gradient pulse in z-direction to destroy all x−y plane signals. It is usually applied before the primary circuit, and the inhomogenous portion will stay at no-signal case during the following pulse sequence. A typical RF-selection sequence with 64 loops is . When the molecules feel perfect RF amplitude, their states remain as thermal equilibrium after this sequence. By contrast, when the molecules feel for example 4.5% error in RF amplitude, their states mostly evolve to x−y plane and thus be killed by the following gradient field. Note that although RF-selection Figure A1. Molecular structure of dichloro-cyclobutanone, where C 1 -C 7 form a seven-qubit system. The diagonal elements are chemical shifts (Hz), and the off-diagonal elements are scalar coupling strengths (Hz). T 1 and T 2 are the relaxation times (second) of the individual spins, respectively. All parameters are obtained on a Bruker DRX 700 MHz spectrometer at room temperature. enables a better SNR as the RF pulses are much more precise, the cost of this technique is the absolute loss of signal as many molecules have no contributions to the signal any longer.
In our experiment, we used a GRAPE pulse instead of the long sequence to realize this RF-selection technique. This GRAPE pulse was designed on a single-qubit system via the state-to-state approach, by setting two constraints: evolve Z to x−y plane when the RF inhomogeneity is more than 1%, or else do nothing to Z. After applying this GRAPE pulse on our seven-qubit system, we found the signal reduced to about 30% but the RF pulses were indeed much more homogeneous by running the Rabi oscillation experiment. The two gradients and p 2 rotations in step d are used to kill the minor signal of multi-coherence generated by the J-coupling evolution during the RF-selection sequence. The state at step d is the same as step c, but with some loss that r =´Z 30% 4 d C 7 . For convenience, we simply mark this state as Z C 7 . Compared to the original thermal equilibrium state, this new state gains signal boost from H 5 and owns much better RF homogeneity.
e. The main body of cat-state method [31] is implemented which creates the labelled PPS | | r = ñá Z 000000 000000 PPS C 7 from Z C 7 . It consists of three steps: encoding, phase cycling, and decoding. The detailed NMR sequence is shown in figure B1(b). Starting from Z C 7 , the system evolves to Z Z Z ...   figure 5 for its NMR spectrum. Hi by rotating all the carbon spins by p 2 around the y-axis (Y90) followed by a gradient field (GR); b: boosting the polarization of C 7 by exchanging the state of C 7 and H 5 (SWAP); c: decoupling the 1 H channel for the rest of experiments; d: the RF-selection targeted on C 7 ; and e: labelled PPS preparation. The above steps are repeated for seven times with different phases of y j to select the appropriate coherence. For simplicity, the rest of hydrogen spins are not shown in the figure. (b) Detailed sequence of step e in the above circuit. 7 . When these stabilizer operators evolve under the ground state preparation circuit shown in figure 3, one can reconstruct the stabilizer operators of the ground state of the seven-qubit toric code. However, since our circuit starts from | | ñá Z 000000 000000 C 7 , the S pps are modified to Z , Z , Z  Figure B2. Simulated non-trivial phase q a (loop 1) as a function of T 2 scaling factor η. To observe the decoherence effect on the anyonic phase, T 2 values are coherently varied for all the 13 C's by scaling the T 2 values shown in figure A1 by η. As expected, as T 2 gets longer, the anyonic phase approaches values closer to the expected value π. The phase of the other non-trivial loop 2 has almost the same behaviour compared to this one. For this particular simulation, T 2 effect is considered by assuming the initial state and readout stages are perfect.

19
Therefore, the experimentally prepared ground state have +1 expectation values of the above transformed operators S ground . We measured the expectation values of S pps of the labelled PPS state and the expectation values of S ground of the ground state. For the S pps measurements, a single readout pulse which rotates C 7 by p 2 around y-axis is sufficient to measure all six operators; whereas five different readout pulses are required (thus, five different measurements) to measure the S ground operators. The readout pulses are composed of the single qubit rotations that transform the product operator components of a density matrix corresponding to the S ground operators to the measurable product operators in C 7 spectra, which are a combination of X C 7 or Y C 7 and different Z C i , where i indicates the ith 13 C. For instance, the readout pulse required to measure the expectation value of the second (#2 in the above list) S ground operator is which rotates Y X Y X   figure C2 were fitted using the least-square method to estimate the coefficients of the peaks. The expectation values of the desired operators were evaluated by taking the appropriate linear combinations of the estimated coefficients [43]. The error bars were calculated by using the method of the error propagation with the initial standard deviations from the fitting procedure. Figure C2 shows the experimental and simulated spectra of C 7 after the ground state preparation which were measured by the five different readout pulses. These spectra were used to estimate the S ground operators.