Topological contextuality and anyonic statistics of photonic-encoded parafermions

Quasiparticle poisoning, expected to arise during the measurement of Majorana zero mode state, poses a fundamental problem towards the realization of Majorana-based quantum computation. Parafermions, a natural generalization of Majorana fermions, can encode topological qudits immune to quasiparticle poisoning. While parafermions are expected to emerge in superconducting fractional quantum Hall systems, they are not yet attainable with current technology. To bypass this problem, we employ a photonic quantum simulator to experimentally demonstrate the key components of parafermion-based universal quantum computation. Our contributions in this article are twofold. First, by manipulating the photonic states, we realize Clifford operator Berry phases that correspond to braiding statistics of parafermions. Second, we investigate the quantum contextuality in a topological system for the first time by demonstrating the contextuality of parafermion encoded qudit states. Importantly, we find that the topologically-encoded contextuality opens the way to magic state distillation, while both the contextuality and the braiding-induced Clifford gates are resilient against local noise. By introducing contextuality, our photonic quantum simulation provides the first step towards a physically robust methodology for realizing topological quantum computation.

Quasiparticle poisoning, expected to arise during the measurement of Majorana zero mode state, poses a fundamental problem towards the realization of Majorana-based quantum computation. Parafermions, a natural generalization of Majorana fermions, can encode topological qudits immune to quasiparticle poisoning. While parafermions are expected to emerge in superconducting fractional quantum Hall systems, they are not yet attainable with current technology. To bypass this problem, we employ a photonic quantum simulator to experimentally demonstrate the key components of parafermion-based universal quantum computation. Our contributions in this article are twofold. First, by manipulating the photonic states, we realize Clifford operator Berry phases that correspond to braiding statistics of parafermions. Second, we investigate the quantum contextuality in a topological system for the first time by demonstrating the contextuality of parafermion encoded qudit states. Importantly, we find that the topologically-encoded contextuality opens the way to magic state distillation, while both the contextuality and the braiding-induced Clifford gates are resilient against local noise. By introducing contextuality, our photonic quantum simulation provides the first step towards a physically robust methodology for realizing topological quantum computation.

I. INTRODUCTION
Non-Abelian anyonic systems, in which information is non-locally encoded, offer an attractive approach to faulttolerant quantum computation [1,2]. Anyonic quantum states are protected by topological order and are intrinsically immune to erroneous local perturbations [3]. Moreover, their non-Abelian statistics [4] naturally induce topologically quantum gates in the encoded space that are immune to a wide range of control errors. Presently, Majorana zero modes remain the most conspicuous realization of the non-Abelian anyons. Several experimental signatures of Majorana fermions have been confirmed indicating their existence [5,6]. Nevertheless, the realization of their non-Abelian statistics remains elusive. An important obstacle towards the experimental realization of Majorana braiding is that these systems suffer contamination by leaked quasiparticles from the environment, a mechanism known as quasiparticle poisoning [7][8][9][10][11]. This notorious mechanism dramatically reduces the coherence times of quantum states encoded in Majorana zero modes, thus creating significant challenges in the verification of their non-Abelian statistics by adiabatic braiding evolutions. * These two authors contributed equally to this work. † j.k.pachos@leeds.ac.uk ‡ jsxu@ustc.edu.cn § smhan@ustc.edu.cn ¶ cfli@ustc.edu.cn To overcome the above demerit of Majorana fermions, the use of more exotic non-Abelian anyons is proposed.
Parafermions [12] with Z n symmetry (n positive integer), a natural extension of Majorana fermions that possess Z 2 symmetry, provide symmetryprotected n-fold degenerate ground states which can be utilized to encode an n-level qudit. Remarkably, due to the existence of a composite topological charge for a pair of parafermions in superconducting fractional quantum Hall systems [13,14], the quasiparticle poisoning in the parafermions-based quantum computation can be substantially suppressed [15]. Similar to Majorana zero modes, braiding of parafermions only generates the Clifford operators and is not solely competent for universal quantum computation.
However, various approaches have been proposed to implement universal quantum computation with parafermions, for example, by obtaining Fibonacci anyon [16] from a nucleated phase of Z 3 parafermions [17], or by asymptotically simulating the non-Clifford gates, such as the π/8phase gate, with magic state distillation [18].
To determine if magic state distillation can be employed, one can use the onset of quantum contextuality [19] that benchmarks resource suitability for universal quantum computation [20]. Moreover, the contextuality encoded in anyons is topologically protected, making the process robust to noisy input states.
Realizations of parafermions have been proposed [21][22][23][24][25] by manipulating fractional quantum Hall states. The realization of such a proposal is technically challenging, and only very recently the experimental signatures of arXiv:2011.05008v2 [quant-ph] 11 Jul 2021 generic anyons have been observed [26,27]. Quantum simulations [28,29], on the other hand, have been very successful in demonstrating key ingredients of anyonic systems. Examples include the photonic implementation of anyonic braiding [30,31] and quantum gate [32] by virtue of Majorana zero modes. In this article, we apply the idea of photonic quantum simulation, together with a novel form of dense encoding, to demonstrate the two key ingredients of the universal quantum computation with Z 3 parafermions: the braiding to generate Clifford group and the contextuality in the topologically encoded space. The photonic simulator has the advantage of being virtually decoupled from the environment. Meanwhile, it supports the realization of accurate quantum gates, so it is ideal for implementing Berry phase evolutions that conclusively identify the anyonic statistics of parafermions.
When the parafermionic states are faithfully encoded with the beamsplitter network, their subsequent manipulations are robust against control errors due to their topological character. Moreover, our experiment comprises a quantum simulation of an intrinsically interacting system [33], the Z 3 parafermions chain, compared to the Majorana chain, which can be described by non-interacting fermions [34]. Due to the peculiarities of the photonic simulation, the encoding of the quantum states of the parafermion chain in an optical beamsplitter network lacks a tensor product structure, thus the size of the Hilbert space increases only polynomially with the employed photonic resources. Nevertheless, the recently empowered modulation [35] and detection [36] of hundreds of spatial modes paved the way to large-scale photonic quantum simulation. We envision that such photonic simulations will make it possible to implement prototype topological quantum algorithms in the near future [32].

II. BRAIDING OF PARAFERMION EDGE ZERO MODES
Our starting model is a parafermion chain with two Z 3 parafermions (α ja and α jb ) at each site j. The parafermion operators satisfy the commutation relation α ja α ka = ωα ka α ja , α jb α kb = ωα kb α jb and α ja α kb = ωα kb α ja for j < k, where ω = e 2πi/3 . For simplicity and clarity, throughout this article we focus on a chain with three sites, j = 1, 2, 3, and manipulate it in order to obtain the parafermion statistics as a Berry phase evolution. We start from the configuration subjecting to the following Hamiltonian This chain has two Z 3 parafermion edge zero modes (α 1a and α 3b ) emerging at the edges of the chain, as illustrated in the middle panel of Fig. 1. The term "zero" mode is coined because these modes are free operators vanishing from the Hamiltonian, thus contributing zero energy shift when populated. The ground state of H Pf 0 is 3-fold degenerate and can encode a non-local qutrit. For convenience, we choose the basis |ψ Pf l (l = 0, 1, 2) as the eigenstates of the parity operator Q = k α † ka α kb with eigenvalues ω l .
The non-trivial statistics of parafermions is manifested by the emergence of different phase factors between the above eigenstates when the parafermion edge zero modes, α 1a and α 3b , are exchanged [37]. To braid these two parafermion edge zero modes, we can adiabatically cycle between the Hamiltonians H 0 , H 1 and H 2 , where Here, the ferromagnetic term α 1a α † 1b in (2) resembles the effect of a local external field, while the other terms are analogous to a paramagnetic interaction between neighbor sites. Therefore, all employed Hamiltonians have a clear physical interpretation with components that are amenable in quantum Hall systems. We also note that, in our configuration the external fields do not overlap with the interacting parts of the chain that can freely support parafermion edge zero modes. In the case they overlap the maximum tolerable external field is determined by the chirality of the interactions, represented by the phase factor of the α jb α † (j+1)a terms [12]. Our choice of π/6 makes the parafermion edge zero modes most robust.
We expound the effect of the adiabatic evolution, H Pf 0 → H Pf 1 → H Pf 2 → H Pf 0 , on the parafermion edge zero modes, in the left panel of Fig. 1. First, by turning off the interaction between sites 1 and 2, we swap the parafermion edge zero mode initially on site 1a and the auxiliary parafermion at site 2a. Second, we interact the solitary site 1 with site 3, so the parafermion edge zero mode initially situated on site 3b is driven to site 1a which is now at the end of the interaction chain; finally, we transfer the parafermion edge zero mode from site 2a to site 3b by resetting the system's couplings to their configuration. Modulo the braiding with the auxiliary parafermion which has no effect, the evolution results in the exchange of the two parafermion edge zero modes α 1a and α 3b . As the braiding operator and the parity operator can be simultaneously diagonalized [37], the braiding induces a Berry phase acting on the degenerate subspace spanned by |ψ Pf 0 , |ψ Pf 1 and |ψ Pf 2 .
A. Simulation of parafermionic dynamics with a bosonic spin chain Realization of parafermions based on the fractional quantum Hall system is highly non-trivial and still intractable with current technologies.
Fortunately, through the Fradkin-Kadanoff transformation [38]  Braiding of parafermion edge zero modes. Left: braiding two parafermion edge zero modes located on a parafermionic chain's two remote sites using an auxiliary parafermionic site. The taichi symbols with green and blue edges denote two parafermion edge zero modes, the silver disk with pink edge denotes the auxiliary parafermion and the solid disks denotes idle parafermionic sites that do not participate in the braiding procedure. The braiding operation yields a Clifford gate on the encoded parafermions qutrit. Middle: physical realization of the braiding operation on an interacting parafermion chain. By adiabatically tuning the system Hamiltonian to H Pf 0 , H Pf 1 , H Pf 2 and again H Pf 0 (see Eqs. (1) and (2)), the system evolves to different configurations that result in two parafermion edge zero modes being exchanged. Note that between H1 and H2 the green parafermion undergoes a type-2 Reidemeister move which does not affect the braiding process. The curled lines denote the interaction between different sites, and the rounded rectangles pair the parafermions that belong to a single site. Right: Applying the Fradkin-Kadanoff transformation on the chain yields a spin-1 quantum Potts model with Z3 symmetry. Its evolution can be effectively simulated by a series of imaginary-time evolutions which repeatedly project the system onto the ground state of new Hamiltonians. In each subpanel, the spin states corresponding to the empty sphere have an eigenvalue of 1 on the basis denoted as a subindex of the ket symbol, the spheres with upright arrows have the same eigenvalues, and the spheres with tilted arrows have spin eigenvalues that sum up to the sphere with inverted arrow (modulo 3; see the definition of the eigenstates after Eq. (3)). The dashed horizontal lines link the configurations in the three different pictures at each stage of the evolution. model that describes a chiral spin-1 chain H S [12]. As the Fradkin-Kadanoff transformation is non-local, it scrambles the topological order of the parafermions, so the corresponding spin-1 system is only symmetry protected, i.e., it is only naturally resilient against noise with Z 3 symmetry. To reproduce adiabaticity during braiding we artificially constrain the dynamics of our photonic system into the ground state subspace of the spin chain by suppressing the wavefunction leaked into excited states. However, as the spectra of the parafermions and the spin-1 chains are identical, the statistical behavior, such as the Berry phase resulting from the braiding of parafermion edge zero modes, can be identically produced by the equivalent manipulations of the chiral spin-1 chain.
The Fradkin-Kadanoff transformation maps the parafermionic Hamiltonians (1) and (2) into spin Hamiltonians with three sites of the form: where the subindex corresponds to the three sites j = 1, 2, 3 of the chain. As the Hamiltonians H S and H Pf have identical spectra, we know that the adiabatic evolution H S 0 → H S 1 → H S 2 → H S 0 of the spin system shall reproduce the Berry phase resulting from the adiabatic evolution of H Pf that corresponds to the parafermion edge mode braiding. It is convenient to define three sets of bases for each site j = 1, 2, 3 of the chain, namely |k σ , |k τ and |k χ , which are eigenstates of the operators σ j , τ j and σ j τ j , respectively, with eigenvalues ω k . In terms of these states, the correspondence between the basis vectors of the parafermionic and spin picture in the degenerate spaces are as follows: The representation of the parafermions states in terms of spin states makes it possible to realize them with a photonic simulator, as we shall see in the following subsection.
B. Berry phase and dense encoding The first task of our quantum simulation is to obtain the Berry phase evolution that corresponds to the statistics of parafermions. The Berry phases induced by braiding two parafermion edge zero modes are given in terms of projections by where l = 0, 1, 2 parametrizes the degenerate eigenstates of H S 0 , while Π i projects the wavefunction on the ground states of the spin-1 Hamiltonian H S i (i = 1, 2) given in (3). Note that in general φ B,l is not the same for different l's indicating the non-Abelian character of the parafermions. To realize this Berry phase, we implement the cyclic imaginary time evolution operator given by [30] The parameter t is chosen to be large, so each term e −H S t acts as a projector Π onto the space of ground states of H S ; thus (6) implements the Berry phase (5). A vital characteristic of the imaginary-time evolution method is that we do not need to know the ground state of every involved Hamiltonian in order to carry out the simulation; this is because whenever two terms in Hamiltonian H S commute, the imaginary time evolution operator (6) can be further disassembled into two consecutive projections corresponding to each of the terms. Furthermore, as the system is finally projected back onto H S 0 's ground states, the first term of H S 1 and H S 2 in (3) can be completely omitted. As such, we only need to project the system three times onto the eigenstates of the operators σ j , τ j and σ j τ j , which are |k σ |k σ |k σ , |0 τ |k σ |k σ and |k σ |l τ |l + k mod 3 χ , respectively, with k, l ∈ {0, 1, 2}. This simulation procedure is shown in the right panel of Fig. 1.
The encoding of the parafermionic Berry phase in the imaginary-time evolution of a spin system can be further simplified before proceeding to the actual experiment. To elucidate this point, we note that although the set of vectors |k σ |j τ |j + k χ (j, k = 0, 1, 2) lives in a nine-dimensional Hilbert space and it needs a qudit with d = 9 to be encoded, projecting from this subspace onto the degenerate space of H S 1 or H S 0 only yields three independent coefficients. More concretely, the destination amplitude of a basis rotation from the eigenstates of H S 1 or H S 0 to |k σ |j τ |j + k χ is solely determined by k, regardless of j, so the projection onto different j's can be assimilated into the same physical mode, and only three orthogonal modes are required throughout the simulation. In other words, by densely encoding different eigenstates of H S 2 , it is sufficient to simulate the braiding of a parafermionic chain through on-demand imaginary-time evolution on a physical qutrit (d = 3). Our dense encoding method can be directly generalized to the quantum simulation of Z n parafermions.

III. A PHOTONIC PARAFERMION SIMULATOR
The spin-1 chain system given in (3) can be directly encoded into an off-the-shelf photonic simulator. In particular, we encode the ground states of the homomorphic spin-1 chain into different spatial modes of a single photon, which constitutes a qutrit system. Since we only encode the information in the degenerated ground space, we can cast imaginary-time evolution (e −H S t ) corresponding to different Hamiltonians in the photonic setup to enforce the encoded spin system to their ground spaces. The consecutive implementation of photonic imaginary-time evolution is the core procedure of our dissipative quantum simulation.
The optical simulator that realizes generic preparation, coherent evolution and the imaginary time evolution is shown in Fig. 2. The beam displacers, which are birefringent crystals separating the horizontally-and vertically-polarized propagating beams, are employed to prepare the photonic path qutrit. Before and between two sets of beam displacers, half-wave plates and quarterwave plates are inserted for the on-demand control of the complex amplitude of the three modes. In principle, a 3 4 -  Parafermion photonic simulator. The qutrit state, |ψ Pf , of the parafermion edge zero modes, is encoded on photonic paths. The preparation of the qutrit state in the encoder is implemented using arrays of beam displacers. The halfwave plates and quarter-wave plates are exploited to control the complex amplitude of each mode. The photonic modes are then subjected to the required braiding evolution gates consisting of two phase gates P inducing phase unitaries between the three modes and a rotation gate R inducing a generic SU(3) rotation. Then, the architecture of beam displacers is again exploited in the decoder to recover the status of the parafermion edge zero modes by tomographic measurements. In the subplots with dashed boxes, we give the detailed photonic implementation of the P and R gates, where the individual modulation of the nine optical modes inside R is realized by additional amplitude filtering (not shown). Throughout the setup, we insert the temporal compensators to align the optical lengths of the modes so that the interference visibility is maximized and tilt them to introduce the necessary phase compensations. The brightness of the photon source in this experiment was about 100K counts per second on the maximum detection basis, and the integration time was 10 seconds per data point.
dimensional Hilbert space is required to fully support the complete evolution of the Z 3 parafermions braiding. By employing the method presented in [32], which restricts to the ground state's dynamics, 3 3 independent optical modes are required. We additionally employ the novel dense encoding presented in Sec. II B to further suppress the required resources in the simulator to only 3 2 modes. This encoding extends the range of operations that can be simulated with our architecture and facilitates the realization of parafermions braiding evolutions.
The basic building blocks of the Berry phase evolution are the qutrit projectors Π i , which implement the imaginary-time evolution. When the projector and the current Hamiltonian have the same set of eigenstates, then the projector is represented by a phase unitary realized by a non-dissipative phase gate P, which induces relative phases between the three optical modes; a P gate is constructed using two quarter-wave plates with their optical axis fixed at 45 • , sandwiching three adjustable half-wave plates acting on each path independently. When the projector does not share the same eigenstates with the current Hamiltonian, it must be implemented dissipatively. An arbitrary imaginary-time evolution is realized by a rotation gate R comprised of first a basis rotation and a subsequent dissipation of complex amplitudes of the eigenstates with higher energies. The R gate operates by first horizontally separating three vertically displaced optical modes with beam displacers twice in a row; the resulting nine complex amplitudes are then individually controlled; subsequently, these modes are vertically merged to obtain a horizontally arranged qutrit, with some photons dissipated as the result of imaginary-time evolution. This procedure is shown in the subpanels of Fig. 2.
Finally, in the decoding stage, another pair of beam displacers combines the three modes back. The additional half-and quarter-wave plates, together with a polarizing beam splitter, are employed for qutrit tomographic measurement, as shown in Fig. 2. As the measurement relies on interference between different optical modes and the photons used in our experiment have a coherence length of about 200µm, the total optical length of each evolution channel must be aligned, so the differences between them are much smaller than the coherence length in order to preserve the temporal coherence and maximize interference visibility.
In our experimental configuration, we introduce additional temporal compensators to compensate for the path differences introduced by the wave plates and optical misalignments and slightly tilt them to ensure the correct relative phase between different paths.

A. Observation of parafermionic Berry phase in the simulator
To investigate the effect of the parafermions braiding operations with the photonic simulator, we utilize nine generic qutrit states, |ψ j , j = 1, ..., 9 encoded in the three degenerate ground states of H S 0 and then implement the simulated braiding evolution on them. Quantum state tomography is then employed to identify the resulting states. The tomographic results, together with the theoretical predictions, are plotted in Fig. 3(a). By translating the braiding evolution of the parafermions chain into three consecutive basis rotations and projections onto corresponding ground states, we find that two P gates and one R gate are required. From this procedure, we experimentally obtain that an extra phase of δφ B,2 = 2.061 ± 0.128 emerges between |ψ S 2 and |ψ S 0 , while the phase between |ψ S 1 and |ψ S 0 is insignificant, taking the value δφ B,1 = 0.066 ± 0.104. These Berry phases correspond to 2π/3 and 0, respectively, in agreement with the expected Z 3 parafermion braiding matrix diag(1, 1, ω) [37], with an average fidelity of 94.9%. Note that φ B,l is not the same for different l's emphasizing the non-Abelian character of the parafermions. The data in our experiment were recorded with single photon avalanche detectors; in order to estimate the counting errors, we assume a Poisson distribution for the counting statistics and resample the coincidence event numbers for 100 times to calculate the 1σ standard deviations of desired quantities.
Quantum process tomography is a reliable tool for the complete characterization of the braiding dynamics. From the tomography, we deduce the process matrix χ jk of the evolution, spanned by the Gell-Mann basis {λ 1 , ..., λ 8 } and the identity I, with maximum likelihood estimation. The resulting process matrix corresponds to a relative phase of 2π/3 between |ψ S 2 and the other two states, as shown in Fig. 3(b). By comparing the result from the braiding quantum simulation with the expected theoretical value, we obtain the fidelity of the operation to be 93.4%. Therefore, the braiding of the parafermions is faithfully simulated by a series of imaginary time evolution operators realized by mesoscale linear optical gates [30].
We now comment on the sources of error present in the experimental simulation of parafermion braiding. The imaginary-time evolution of the system dissipates all optical modes that correspond to excited parafermionic states. As a result all the recorded errors occur inside the ground degenerate space, and are solely caused by experimental imperfections. From the quantum process tomography we see that the decoherence due to experimental imperfections plays a crucial role in the reduction of our measured fidelity. This arises in the photonic experiment from the different times the photons take to travel along different paths, causing the final interference fringes to have less than 100% visibility. To a lesser degree, the decrease of fidelity is also  (a) We employ nine states |ψj , j = 1, ..., 9, to experimentally implement quantum state tomography before and after the braiding process. The top, green Bloch sphere shows the The black and white points denote the theoretical initial and final states, respectively. In the top subplot, no significant relative phases can be seen in agreement with the theoretically predicted zero Berry phase between states |0 and |1 resulting from parafermion braiding. In the bottom subplot, each pair of points is connected by a circular arc of 2π/3 central angle, which is the Berry phase induced between|0 and |2 by braiding the parafermions. The location of the points is calculated using the expectation values of the Gell-Mann matrices {λ1, ..., λ8} and the identity matrix, I. (b) The real (top) and imaginary (bottom) parts of the process matrix elements, χ jk , expressed in the Gell-Mann basis. The edge and filling of the cuboids show the theoretical values and experimental results from quantum process tomography, respectively. caused by misalignment in the optical elements, element imperfections, and counting statistics.

IV. PARAFERMION CONTEXTUALITY FOR QUANTUM COMPUTATION
It is plausible that our experimental simulation of Z 3 parafermions braiding can be extended to more complicated systems with several chains that encode topologically encoded qutrits in their parafermion edge zero modes just like the Majorana zero modes in [32]. In this case, braiding between arbitrary pairs of parafermion edge zero modes can generate all Clifford operators in the encoded space.
As circuits with Clifford gates can be efficiently simulated by classical computers [39], exclusively parafermions braiding evolutions are not universal for quantum computation. A highly promising method for rectifying this problem is to employ magic state distillation on non-classical states generated e.g. by dynamical evolutions [5] or by employing half-fluxons [15]. Magic state distillation not only complements the Clifford circuit for producing universal quantum computation, but also has a fault-tolerant character. State-dependent quantum contextuality is a necessary criterion that identifies if magic state distillation is applicable in a circuit or not [20]. In other words, quantum contextuality and Clifford operations bootstrap the full power of universal quantum computation.
In this section, we experimentally demonstrate that a parafermionic universal quantum computation architecture is possible based on the fault-tolerant procedures of braiding and magic state distillation, as summarized in Fig. 4(a). Assuming that we can generate quantum states that are not classically simulable, e.g., by activating interactions between parafermion edge zero modes, then magic state distillation can be applied at an arbitrary stage of the parafermions-based quantum computation to provide the non-Clifford gates required for universality.
In the first subsection, we experimentally determine the contextuality properties of various topologically encoded states that enable magic state distillation. Specifically, we demonstrate that contextuality remains unaffected by the application of braiding-induced Clifford gates. In the second subsection, we elucidate that the parafermion encoded contextuality inherits the faulttolerant feature from topological order. In the last subsection, we use another contextuality criterion to selftest the encoded state, which best describes its merit of local noise resilience. Therefore, we demonstrate that fault-tolerance extends to all the necessary elements of our parafermionic quantum computation architecture.

A. Contextuality and its interplay with braiding
Quantum contextuality [19] is one of the most profound non-classical features of quantum theory. It states that an observable cannot have a predefined value regardless of the specific orthonormal basis used to measure it. The discrepancy between quantum and classical theories, where the observables have contextualindependent values, is revealed by the violation of non-contextual hidden-variable (NCHV) inequalities [40]. Given a set of projective measurements, a NCHV inequality can always be formulated in the light of a recent graph-theoretic approach [41]. In particular, when the set of projectors corresponds to the qudit stabilizing group, the resulting NCHV inequalities bound the set of classically simulable states. As a result, a state that violates the inequality is suitable for quantum computation with magic state distillation [20].
The construction of this NCHV inequality relies on the stabilizer formalism. Explicitly, the Weyl-Heisenberg displacement operators for a qutrit system are defined as Each of these operators has a spectrum of {1, ω, ω 2 }. For convenience, we denote the list of displacement operators D = {D 0,1 , D 1,0 , D 1,1 , D 1,2 }, whose eigenstates span a complete set of mutually unbiased bases. By exploiting the above notation, a witness of contextuality can be constructed from the operators where Π rj j is the projector of the eigenstate ω rj of the j-th element in D. The vector r is in turn defined as r = xa + zb with a = {1, 0, 1, 2}, b = −{0, 1, 1, 1} and {x, z} ∈ {0, 1, 2}. For simplicity, we also adopt the shorthand notation A xz = A xa+zb . There are in total 9 witness operators, and the one with the maximum expectation value determines the suitability of a quantum state for magic state distillation. In particular, the NCHV inequality that identifies the magic state distillation resource states reads where ρ is the quantum state to be considered. Because the projectors Π rj j of stabilizing operators form a closed group under Clifford operations, the witnesses M of contextuality are also unaffected by Clifford gates applied on the state ρ.
Experimentally, we apply the contextuality test (9) to the nine photonic states ρ i = |ψ i ψ i |, i = 1, ..., 9 shown in Fig. 3(a), and determine the value of M for these states before and after braiding with the quantum simulator in order to determine their contextuality properties. Fig. 4(b) illustrates the dynamics of the contextuality witness operators A ij , with i, j = {0, 1, 2} during braiding, which form a closed subgroup, acquiring each another's values. Direct comparison between these two groups of results shows that the values of the contextuality witnesses M are almost not affected by the process of braiding, as shown in Fig. 4(c). Consequently, the Clifford operations induced by braiding parafermions do not affect the state-dependent contextuality required for magic state distillation and can be safely performed at all stages of the computation.

B. Observation of fault-tolerant contextuality
The omnipresent local noise in a non-topological quantum system may diminish the observed contextuality, necessitating more copies of resource states to achieve a target fidelity in magic state distilation [18], or even completely hinder its implementation when contextuality is not present. In contrast, quantum states encoded in parafermionic systems are topologically protected. It is thus expected that the contextuality properties of these states are topologically immune enabling the application of magic state distillation even in the presence of local noise.
To demonstrate this aptitude of parafermion systems, we apply local noise to both a topologically-encoded and a non-topological qutrit, and compare their response to noise. In the topological system, the flip of a single parafermion is prevented by parity conservation, so only simultaneous flips of two sites are allowed. As a result, the parafermionic system is susceptible to two possible forms of local noise: local hopping errors, which occur simultaneously on two adjacent sites, and local phase errors, which act on single sites. In order to characterize their behavior, it is convenient to define the Fock parafermion operators, C j , associated with the canonical parafermions operators α ja and α jb by Then, the local hopping errors and the local phase errors can be expressed as C † k C k+1 and, C † k C k , respectively [42], occurring with a certain probability.
We first consider analytically the effect of hopping noise, C † 1 C 2 , between sites 1 and 2. Through the Fradkin-Kadanoff transformation, the corresponding operator in the spin picture reads where the superscript S indicates that the operator is expressed in the spin picture. If we apply the hopping error with probability p on the three computational basis states we find: where δ ij is the Kronecker-δ function. Hence, applying the neighborhood hopping operator only causes the parafermions to leak out of the degenerate subspace and has no effect on the topological qutrit up to a state renormalization. Similarly, the local phase noise on the parafermionic system with probability q takes the form with Direct calculation shows that the result in (13) still holds; namely, the dephasing noise always gives rise to an increase in system's energy and is thus corrected during imaginary time evolution.
To understand the protection exhibited by the topologically encoded qutrit, we compare it to a normal qutrit state subject to noise. We act on this state with two error models that resemble the errors in the parafermionic system, namely, the probabilistic flip error caused by the shift operator, τ , occurring with probability p, and dephasing errors caused by the clock operator, σ, occurring with probability q, where τ and σ are given in (3). Their corresponding effects on a state ρ can be described by the superoperators and While the qutrit error operators τ and σ are in principle different to the operators X and Z acting on the parafermion qutrits, they correspond to the type of errors expected to emerge in typical physical realizations.
To quantify the effect of errors, we exploit the behavior of the contextuality witness M given in (9), which also quantifies the usefulness of the state in magic state distillation. For concreteness, we start with the state which is an almost optimal state for magic state distillation, encoded in the topological or non-topological qutrit. Then, we numerically analyze the response of the contextuality witness M with respect to joint hoppingphase noise in the parafermionic system and flip-phase noise in the non-topological system. These errors are described by the composed noise superoperators X S • Z S and T • Σ, respectively. Fig. 5(a) shows that the presence of flip-phase noise quickly destroys the resource for magic state distillation in non-topological system. In contrast, the parafermionic qutrit is still free from any containment penetrating the degenerate subspace, as shown in Fig. 5(b). We now experimentally investigate the response of the parafermion edge zero mode qutrit state to the local noise in our photonic simulator and compare it with the effect of noise acting on a normal qutrit. We initially encode the state (18) in our simulator. For simplicity, we restrict to hopping or flip errors as the full set of errors studied theoretically requires a larger Hilbert space than what can be supported by our simulator. This noise is simulated in our setup shown in Fig. 2 by removing the two P gates and keeping only the R gate. In particular, for the parafermionic qutrit we separate each of the three optical modes into two parts according to p, the probability of error to occur, and then project back into the computational basis with an imaginary-time evolution, which dissipates the modes outside the ground state subspace. For the non-topologically protected qutrit, we use the R gate to emulate the behavior of the superoperator T . The crucial difference in this case is that the modes in the excited subspace are no longer dissipated, resulting in the lift of the symmetry protection mechanism.
The experimental results for the contextuality witness M against noise are shown in Fig. 5(c). From these data we see that for the parafermionic system, the witness does not decrease with respect to any noise strength, which is attributed to the topological protection. In particular, the parafermion encoded state remains an almost optimal resource state with M 0.580 ± 0.013. Indeed, when hopping noise acts on a parafermionic system, it only decreases the population of uncontaminated parafermions quasiparticles leaving the encoded states unaffected. Consequently, the parafermion system is suitable for magic state distillation even in the presence of noise.
The behavior of the normal qutrit under the presence of noise is rather different. From the experimental data, shown in Fig. 5(c), we see that, for a normal qutrit, the value of M initially decreases for increasing error probability p denoting that its quality as a resource for magic state distillation is diminished. At p = 2 3 , the inequality (9) can no longer be violated and the noise has completely deprived the capability of the non-topological state for magic state distillation-based quantum computation. However, for sufficiently large error probabilities, p > 2 3 , the value for M increases again with p, because the maximum function in (9) is obtained for another operator in A r . This behavior can be best understood from the dynamics of the states in the qubit Hilbert space. As is visualized in Fig. 5(d), an exceeding amount of noise may eventually drive the state to the vicinity of another resource state of magic state distillation. We conclude that, unlike its non-topological counterpart, the qutrit encoded in the parafermion zero modes is able to protect contextuality and thus facilitates quantum computation with magic state distillation even for arbitrary levels of noise. C. Self-testing a topologically-protected qutrit in the face of local noise In the previous section, we used contextuality as an indicator for the applicability of magic state distillation. We now exploit contextually as a tool for self-testing a quantum system.
The term 'self-testing' means that the observed correlations can only be explained, up to local isometries, by a specific state and set of measurements. In other words, a successful self-testing reveals a unique relationship between the state and the measurements. It has been shown that some NCHV inequalities from quantum contextuality are robust self-  (9). The purple points denote two of the magic states, the point marked I denotes the maximally mixed state, and the orange (blue) point represent the final non-topological (parafermionic) state at p → 1. The pink (cyan) area corresponds to the states which are (not) applicable for magic state distillation. The trajectory of the non-topological qutrit under flip noise T is illustrated by the orange arrow, whereas the hopping noise X S has almost no effect on the parafermionic qutrit.
testings: if the choices of projective measurements are fixed, then the maximum violation of the NCHV inequality automatically singles out the form of the tested quantum state [43]. The Klyachko-Can-Binicioğlu-Shumovsky (KCBS) inequality [44] is such an example. Here we use the quantum violation of the KCBS inequality to self-test the encoded parafermionic qutrit in the face of local noise in order to further demonstrate its trait of noise-resilience.
The KCBS inequality is obtained by applying the graph-theoretic approach of contextuality [41] on a pentagon and it has an elementary form. Consider a set of measurements B i , i ∈ {1, ..., 5}, where B i and B i mod 5 +1 are orthogonal. We define P |ψ (B i = 1) as the probability of observing measurement outcome B i = 1 when the system is prepared in |ψ . Then, any NCHV theory assigning dichotomic values 1 for "true" and 0 for "false" regardless of context, satisfies However, a quantum qutrit system admits a maximum witness of contextuality K of √ 5 ≈ 2.236. Hence, the KCBS inequality can be used as a robust selftesting: the maximum violation of the inequality can be achieved by only a single quantum state given the known, optimal set of projective measurements. Moreover, the difference between the experimentally obtained value and the maximum value of K lower bounds the trace distance between the optimal and actual input state.
We first derive theoretically the response of contextuality witness K to noise. In this subsection we again focus on the hopping and flip noises whose realizations are experimentally friendly.
Without loss of generality, let the initial qutrit state be |ψ q 2 = (0, 0, 1) T , and the five optimal measurement settings be B k = |k k|, with Then, under local disturbance X S and T on topological, |ψ S 2 , and non-topological, |ψ q 2 , qutrit, respectively, the value of K in (19) becomes So for topologically encoded states, the violation of KCBS inequality is not affected by local noise, whereas for the normal qutrit states, the violation of (19) is sensitive to the probability p of local hopping errors.
Using the KCBS inequality, we investigate the resilience of contextuality granted by the topological qutrits. The experimental results are shown in Fig. 6(a). The value of K obtained from our measurements on the topological qutrit stays above 2.199 even under the effect of strong local noise. This value significantly violates the prediction of NCHV and confidently shows the existence of quantum contextuality. The proximity of the obtained K with its theoretical maximum value guarantees that the trace distance of our parafermions states before and after perturbation does not exceed the order of O( √ 0.037) [43]. In comparison, when the flip noise with probability p is applied on a normal qutrit, K quickly falls below 2 as the error probability p increases and the violation of the KCBS inequality fails. As a complement to our self-testing method, the process matrix χ jk of noise operator acting on topological and trivial qutrit for p = 2 3 is given in Fig. 6(b) and (c). For the topologically encoded qutrit the process matrix is the identity matrix with fidelity 96.0% indicating the immunity of the parafermions states to local noise. On the other hand, the process matrix of the normal qutrit is significantly disturbed.
The resilience of the parafermion qutrit is due to the increase of the system's energy caused by the application of noise and thus it is suppressed by the excitation gap. Equivalently, in the photonic parafermion quantum simulation, the effect of erroneous channels is dissipated by the imaginary time evolution. The observed noise-immunity clearly showcases the protection offered by topological order in the task of universal quantum computation. The observed reduction in the fidelity of the topologically encode qutrit is due to experimental imperfections. During our demonstration, such imperfections originate from phase mismatch of different photon paths or misalignment and imperfections of the photonic elements cause, on average, an overall reduction in fidelity of 4%.

V. DISCUSSION AND CONCLUSIONS
Topological systems are powerful candidates for future realizations of fault-tolerant quantum computation. A handicap present in both Majorana and parafermionbased systems for quantum computation is the operations created by braiding always belong to the Clifford group, which is not a universal set of quantum gates. To resolve this issue, we resort to magic state distillation as a pathway to acquiring universality.
We have simulated two main elements of parafermions-based quantum computation under this architecture: the parafermions braiding that results in topological Clifford gates and the topologically-protected contextuality that enables the use of magic state distillation.
While contextuality is a vital resource per se for various quantum information tasks like one-way quantum computation [45] and communication [46], its interplay with topological system and noise resilience is not yet well understood. Our work may serve as a trailblazer for encoding contextuality in a real topological system and protecting it from environmental noise, thus benefiting the community of quantum information processing.
From a pragmatic perspective, simulation of a topological system with a realistic quantum simulator provides not only some novel insights on topological quantum computing, but also a testbed for investigating fault-tolerance in the encoded system. For example, quantum teleportation of Majorana zero modes with active error correction has been demonstrated on a superconducting quantum simulator [47]. Even though our photonic simulator is not exponentially scalable due to the use of a single photon for the encoding of the whole system, our method can be directly extended to several parafermions chains and further detect the contextuality of two qudits as proposed in [20]. Thus, it is plausible to expect the simulation of a few parafermions chains in the near future, where small topologically encoded algorithms can be performed in a fault-tolerant way [32]. and No. AHY060300) and the EPSRC (Grant No. EP/R020612/1).
The data that support the findings of this study are available from the authors upon request.
Appendix A: Theoretical Details

Parafermions from the Fradkin-Kadanoff Transformation
A pair of isolated Majorana zero modes can appear at the end of a chain of suitably coupled fermions [3]. A wellknown mathematical equivalence, namely, the Jordan-Wigner transformation, can subsequently connect the fermionic chain to a chain of spin-1/2 particles. The Jordan-Wigner transformation reads: Where γ km denotes the creation and annihilation operator of the Majorana fermion modes m = a, b at the k-th site. Perhaps the most notable feature of the Jordan-Wigner transformation is its non-locality. The constructed fermionic operator is defined with respect to an entire chain of spin operators located on its left hand sites. Due to the non-locality of the transformation, the local properties of the fermionic system are drastically different from those of the spin-1/2 system. Specifically, local noise on the spin system is no longer restricted to one site in the fermionic system. As a result the topological protection is removed.
Similarly to Majorana chains their generalization to Z n parafermion chains (where n = 2 corresponds to Majorana fermions) can be transformed to higher spin chain systems via the Fradkin-Kadanoff transformation. First, the shift and clock operator, τ and σ, are defined as generalizations to the Pauli matrices σ x and σ z to higher dimensions: with ω = e 2πi/n .
The two operators have the commutation relation of στ = ωτ σ.
In analogue to the Jordan-Wigner transformation, the pair of parafermion annihilation operators on the k-th site α ka and α kb is defined with the spin operators as: To demonstrate the effect of the transformation, substituting (S3) into the initial definitions of braiding Hamiltonians (main text formula (1)) yield the braiding Hamiltonians in the spin picture (main text formula (6)). The inverse transformation from parafermion picture to spin picture is: To demonstrate the effect of the inverse transformation consider the parity operator Q = k α † ka α kb . From the first subequation of (S4) one finds that the form of parity operator in the spin picture is The Fradkin-Kadanoff transformation gives rise to exotic commutation relations for the parafermion operators; they neither commute nor anti-commute. Instead, α ka α la = ω sgn(l−k) α la α ka and α kb α lb = ω sgn(l−k) α lb α kb holds. Furthermore, the parafermion Hermitian conjugate relation α † ka = α n−1 ka implies that unlike Majorana fermions, parafermions are not their own antiparticles. Table. SI gives a comparison between Majorana fermions emerging as Majorana edge zero modes and parafermions as parafermion edge zero modes.

Basis Rotation
We use three sets of basis in the spin picture. Namely, Majorana zero modes Parafermion edge zero modes Symmetry

Ground states
Twofold degenerate when J > f > 0 n-fold degenerate when zero modes exist Ground states crosstalk Exponentially suppressed w.r.t. L − 1 Exponentially suppressed w.r.t. L − 1.

Involutory Yes No
Parity Allowed states Eigenstates of (−1) F Eigenstates of ω Q Homologous spin system Ising model Quantum Potts model The eigenstates with index k have eigenvalue ω k with respect to their corresponding operator. Let us span a qutrit |φ on the three sets of bases |φ = A |0 σ + B |1 σ + C |2 σ = a |0 τ + b |1 τ + c |2 τ = α |0 χ + β |1 χ + γ |2 χ , then, the coefficients in each case are linked by: To check the result of the basis rotation, we start from an arbitrary superposition of spin state composed of three ground states of H S 0 , |φ 0 = k 00 |0 σ |0 σ |0 σ + k 01 |1 σ |1 σ |1 σ + k 02 |2 σ |2 σ |2 σ . In the first step of adiabatic evolution, the system is projected on H S 1 's ground state, the new ground state is |φ 1 = |0 τ (k 10 |0 σ |0 σ + k 11 |1 σ |1 σ + k 12 |2 σ |2 σ ). Use the relations of basis rotation, we conclude that The components that do not belong to the ground state of H S 1 have been discarded. In the second step of adiabatic evolution, the system is projected on H S 2 's ground state, which is realized in two steps. Firstly, projecting on the first term of H S 2 , which has no effect because the system is already in its ground state. Secondly, projecting on the second term of H S 2 yields the new ground state, which has the form of 3 p,q=1 k 2pq |p σ |q τ |p + q mod 3 χ . By using the relations of basis rotation, we conclude that As the results for q = 0, 1 and 2 are identical we can omit the subscript q and rewrite the states as In the third step of the adiabatic evolution, the system is projected on H S 0 's ground state, the new ground state is |φ 3 = k 30 |0 σ |0 σ |0 σ + k 31 |1 σ |1 σ |1 σ + k 32 |2 σ |2 σ |2 σ . Projecting the |. τ |. χ terms in the parentheses back to the σ-basis yields: (S11) Experimentally, the three basis rotation are realized separately with an identity phase gate P 1 =    and an inverse Fourier transform gives the result of braiding in the eigenbasis of Q: (S13) Appendix B: Experimental Details

Construction of Quantum Gates
To implement the basis rotation, two sets of quantum gates, the phase gate P, and the rotate gate R, are employed. P adds different phases to different paths which together represent a photonic qutrit; R realizes a general transformation on a qutrit.
Control of the relative phases in the P gate is realized by rotating the three independent half-wave plates between the two bulk quarter-wave plates acting on all of the three paths. Because the polarization states of photons in three paths are the same before entering the P gate, without loss of generality, we assume them to be horizontal. The operations of the wave plates are best described by the Jones matrices. More precisely, the polarization states of a photon is denoted by |H = 1 0 for horizontal polarization and |V = 0 1 for vertical polarization. A half-wave plates is a birefriengent crystal, which adds an extra π phase on photons whose polarization is perperdicular to its fast axis. The effect of half-wave plates with its fast axis at θ with respect to the horizontal axis can be represented by hwp(θ) = cos 2θ sin 2θ sin 2θ − cos 2θ . A quarter-wave plates is a birefriengent crystal, which adds an extra π/2 phase on photons whose polarization is perperdicular to its fast axis. The effect of quarter-wave plates with its fast axis at θ with respect to the horizontal axis can be represented by qwp(θ) = cos 2 θ + i sin 2 θ (1 − i) sin θ cos θ (1 − i) sin θ cos θ sin 2 θ + i cos 2 θ . In the P gate, the quarter-wave plates are both aligned with their fast axes at 45 • with respect to the horizontal axis.
The effect of the P gate on one of the photon paths can be expressed as: So a relative phase of 2θ is added to this path when the half-wave plates is set at θ angle. Using the wave plates, individual phase can be precisely adjusted.
Control of the parameters in the R gate is realized by reshuffling the three optical modes and adjust the configuration of the new modes from the old ones. The three modes are first modulated by three independent half-wave plates set at different angles, so there polarization state is shifted from |H to hwp(θ) · |H = cos 2θ sin 2θ . Next, a beam displacer  Calculated phase change during the braiding evolution for each data point. The blue bars stand for the relative phase between |ψ S 0 and |ψ S 2 , with a theoretical value of 2π/3, and the red bars stand for the relative phase between |ψ S 0 and |ψ S 1 , with a theoretical value of 0. (c) Calculated state fidelity between experimental result and theoretical prediction for each data point.

Quantum Process Tomography
Braiding two parafermions ideally induces a unitary transformation that is diagonal on the parity operator's eigenbasis.
Taking the experimental error and decoherence into consideration, the final form of the evolution is a completely positive trace-preserving map between initial and final states.
To completely characterize this map, we perform quantum process tomography and compare the acquired process matrix with its theoretical value to show the precision of applied operations. The effect of a general completely positive trace-preserving map E acting on a density matrix ρ can be expressed as ρ → E(ρ) = jk χ jkÊj ρÊ † k . Here,Ê is an orthonormal basis of the operators and χ is the process matrix. For a unitary evolution U , the theoretical value of χ can be linked by χ jk = c j c * k , where j c jÊj = U .
In our experiment, we reconstruct a completely positive trace-preserving map from a qutrit to a qutrit. The set of preparation and measurement basis are chosen to be |0 , |1 , |2 , where the set of eigenkets |0 , |1 and |2 represents different path modes in the photonic simulator. The basis of Kraus operators are chosen to be the Gell-Mann matrices plus the identity matrix. Explicitly, they read:  Top: real and imaginary parts of sample density matrices for the data points corresponding to p = 2/3 hopping probability from quantum state tomography. The edge (filling) give the theoretical (experimental) values. Bottom: calculated state fidelity between experimental result and theoretical prediction for each data point. In the graph on the right side, the gray and colored bars stand for fidelity with respect to the theoretical state when the effect of flip errors are taking and not taking into account in advance, respectively.

More Results for Braiding
In Fig. S1, we present more experimental results on the photonic simulation of the PF braiding. Specifically, in subplot (a), three density matrices (corresponding to the quantum states |ψ 3 , |ψ 6 , and |ψ 9 in the main text) obtained after braiding are plotted together with their theoretical values. Subplot (b) shows the calculated phase change during the braiding evolution for each data point. The average increase of relative phase between |ψ S 0 and |ψ S 2 is 118.06 • , and the average increase of relative phase between |ψ S 0 and |ψ S 1 is 3.80 • . The calculated fidelities for each pair of theoretical and experimental states are plotted in (c). The values of fidelity range between 0.916 and 0.968.

More Results for Contextuality and Noise Resilience
In Fig. S2, we present more experimental results on the observation of contextuality in a parafermionic encoded system and the testing of the parafermionic system's robustness against local noise. In the top row, the density matrices obtained after applying a hopping error which happens with a probability of 2/3 is plotted in the colored box, together with the theoretical values, plotted in the edge box. The graph with green bars stands for result for parafermionic qutrits, and the graph with red bars stands for result for non-topological qutrits. In the bottom row, we present the calculated state fidelity between experimental result and theoretical prediction at each flip probability. The average fidelity with respect to the undisturbed state for the topological qutrit is 0.98, while for the non-topologically encoded state the fidelity substantially decreases as the flip probability increase, and drops below 0.12 when the probability is 1. These results are plotted in the colored boxes. On the other hand, the fidelity for the trivial state with respect to the theoretical value when the hopping error has been taken into consideration is as high as 0.99.
In a contextuality experiment, the orthogonality of measurement basis has to be checked in order to satisfy the requirement of compatibility and not open the signaling loophole. Table. SII gives the measured detection probability when the preparation state is compatible with the target detection basis, so the detection probabilities on these bases are ideally 0. These measured detection probabilities are of the order of 10 −3 in accordance with the requirement of compatibility.