Photonic simulation of Majorana-based Jones polynomials

Link or knot invariants [1, 2, 3] serve as a powerful tool to determine whether or not two knots are topologically equivalent. In 1984, Vaughan Jones identified a novel type of knot invariant, known as the Jones polynomial [4]. While being conceptually simple, these polynomials are sufficiently powerful in distinguishing between most topologically inequivalent knots. Currently, there is much interest in determining Jones polynomials as they are important in the field of low-dimensional topology and have applications in various disciplines such as DNA biology [5, 6] and condensed matter physics [7, 8]. Unfortunately, even approximating the value of Jones polynomials at certain non-lattice roots of unity falls within the #P-hard complexity class, with the most efficient classical algorithms requiring an exponential amount of resources [9].

Surprisingly, it has been demonstrated that quantum states of non-Abelian anyons can be used to efficiently evaluate Jones polynomials [10, 11]. Non-Abelian anyons have exotic statistics manifested in the evolution of their quantum state when two of them are exchanged or “braided” [12, 13]. These topological evolutions can be related to Jones polynomials in the following way. Consider a set of identical non-Abelian anyons corresponding to the SU(2)_k Chern-Simons theory with $k$ a positive integer [14]. We pair-create these anyons from the vacuum and we move them on the plane so they can span arbitrary $(2+1)$-dimensional paths with their worldlines. Closed worldlines that form a desired link, $L$, are constructed by demanding that the final anyons are also pairwise combined to the vacuum state. The quantum amplitude of such an anyonic evolution is proportional to the Jones polynomial of $L$ evaluated at the polynomial parameter $t=e^{\frac{2\pi i}{(k+2)}}$. The Jones polynomials can be considered as computational primitives for topological quantum computation. Universality is obtained for anyons that correspond to SU(2)_k Chern-Simons theory with $k\neq 2,4$. Anyons with $k=2,4$ generate the Pauli group, which can become universal with the addition of a simple dynamical gate [15]. The fault-tolerance of topological quantum computation is inherited by the topological invariance of the Jones polynomials against continuous deformations of the worldlines that do not change the topology of the corresponding link.

The SU(2)_k non-Abelian anyons are expected to emerge in strongly interacting systems such as fractional quantum Hall liquids [16, 17, 18]. Unfortunately, we are not yet in position to experimentally generate and manipulate SU(2)_k non-Abelian anyons. As an alternative, significant effort is dedicated in the experimental realisation of Majorana zero modes (MZMs) [19, 20, 21]. Even though the computational power of MZMs is equivalent to SU(2)₂ anyons and thus are not universal they are the most plausible candidate for experimentally realising non-Abelian statistics. Still, the practical realisation of specific quantum algorithms on this platform, such as the calculation of Jones polynomials, is yet to be achieved as the braiding operations using MZMs in an experiment remains a challenging and unresolved issue. As an alternative, quantum simulations of MZMs offer a platform to test the viability of these non-Abelian anyons, investigate their properties and determine the resources required for their full realisation.

In this study, we calculate Jones polynomial based on the braiding operation of Majorana zero modes. We experimentally simulate the realisation and manipulation of MZMs through our all-optical multi-state interferometers. By employing this photonic quantum simulator that harnesses two-photon correlations and non-dissipative imaginary-time evolution, we implement a sequence of Hamiltonians that induce specific braiding evolutions of MZMs. This approach enables us to perform two distinct braiding operations and apply them to compute the Jones polynomials for five typical links: the Hopf link, Trefoil knot, the Solomon link, Figure Eight knot and Borromean rings [22], as shown in Fig. 1a. The obtained quantum amplitudes demonstrate the topological non-equivalence of these links. Our research demonstrates the feasibility of simulating non-Abelian braiding operations and the subsequent evaluation non-trivial Jones polynomials, setting the stage for further exploration of non-Abelian anyons. Importantly, the topological encoding of quantum information in non-Abelian anyons and its manipulation by braiding them with each other guaranties the highly desired fault-tolerance of the algorithm, throughout its execution. Moreover, it showcases that our photonic quantum simulator, even with its limited scalability, facilitates the experimental realisation of complex topological quantum algorithms.

Strictly speaking the MZMs are not described by a Chern-Simons gauge theory. Nevertheless, the elementary braiding evolutions of Majorana fermions differ from those of SU(2)₂ anyons by an overall constant phase factor and they both have quantum dimension $d=\sqrt{2}$. Hence, we can mathematically relate the expectation values, $\bra{\phi_{0}}U_{M}\ket{{\phi_{0}}}$, of arbitrary Majorana evolutions $U_{M}$ with the Jones polynomial $V_{L}(i)$ corresponding to the Chern-Simons theory with $k=2$, in the following way: $V_{L}(i)=e^{{9i\pi\over 8}{(w(L)-\chi(L))}}{2}^{n-1\over 2}\langle\phi_{0}|U_{M}|\phi_{0}\rangle$, where $\chi(L)=\chi_{\rm ac}(L)-\chi_{\rm c}(L)$ with $\chi_{\rm ac}(L)$ and $\chi_{\rm c}(L)$ the number of anti-clockwise and clockwise exchanges of the link $L$, respectively, and $U_{M}$ describes the Majorana braiding evolutions that span the link with their worldlines. For the links $L$ considered here obtained from the Markov trace of braiding evolutions $U_{M}$ we always have $\chi(L)=w(L)$, resulting in the simple expression

Note that in the conventional algorithm given in [24], Hadamard test is required to measure the amplitude, i.e. the $Re\bra{\phi_{0}}U_{M}\ket{{\phi_{0}}}$ and $Im\bra{\phi_{0}}U_{M}\ket{{\phi_{0}}}$. However, for the SU(2)₂ anyons the Jones polynomial $V_{L}(i)$ is always a real number (see Supplemental Material for the details). Moreover, to simplify the complexity of our photonic experiment we restrict to the calculation of $|\langle\phi_{0}|U_{M}|\phi_{0}\rangle|$, which gives us the value $|V_{L}(i)|$ of the Jones polynomial up to an overall sign. While this reduces the predictive power of the Jones polynomials, our experiment can still distinguish between several of the targeted links, as we shall see in the following.

We want to realise a set of braiding operations that correspond to a wide family of knots and links, while they are still feasible to implement with our photonic simulator. For that we consider three pairs of Majorana zero modes positioned at the endpoints of three Kitaev chains [25]. Moreover, we choose the size of the chains to be large enough so MZMs are never brought on the same site during the braiding evolution, which would expose the topologically encoded information to local potentials. To achieve fault-tolerance against local potentials at all times during the performance of our algorithm we employ ten fermion sites, as shown in Fig. 1b. Each Tai Chi-like dual-colored sphere represents a Majorana fermion, and two adjacent spheres together constitute a normal fermion. Here, canonical operators $c_{j}$ and $c_{j}^{\dagger}$ can be used for describing these fermions, with position $j$=1,…,10. The first chain is composed of $j$=1, 2; the second chain is composed of $j$=4, 5, 6 and the third chain is composed of $j$=8, 9, 10. Sites $j$=3 and $j$=7 corresponds to the connections between chains that facilitate the braiding. The Hamiltonian of the initial state in the Kitaev chains can be written as

where $\gamma_{ja}=c_{j}+c_{j}^{\dagger}$ and $\gamma_{jb}=i(c_{j}^{\dagger}-c_{j})$. The Majorana operators $\gamma_{m}$ satisfies the relations $\gamma_{m}^{\dagger}=\gamma_{m}$ and $\gamma_{l}\gamma_{m}+\gamma_{m}\gamma_{l}=2\delta_{lm}$ for $l,m=1a,1b,...10a,10b$. Operators $\gamma_{1a}$, $\gamma_{2b}$, $\gamma_{4a}$, $\gamma_{6b}$, $\gamma_{8a}$ and $\gamma_{10b}$ are not in the initial Hamiltonian $H_{M_{0}}$, thus $[H_{M_{0}},\gamma_{j}]=0$ for $j=1a,2b,4a,6b,8a,10b$. Hence, these Majorana modes have zero energy, corresponding to six endpoint MZMs, denoted from A to F in Fig. 1b.

After completing the two sequences of Hamiltonians of braiding operation $\sigma_{1}$ and $\sigma_{2}^{-1}$ (see Methods for the details), the theoretically resulting braiding evolutions in the logical basis are given by $L_{\sigma_{1}}=\frac{1}{\sqrt{2}}(\mathbb{1}\otimes\mathbb{1}+i\sigma_{x}\otimes\sigma_{x})\otimes\mathbb{1}$ and $L_{\sigma_{2}^{-1}}=\frac{1}{\sqrt{2}}\mathbb{1}\otimes(\mathbb{1}\otimes\mathbb{1}-i\sigma_{y}\otimes\sigma_{x})$. In $L_{\sigma_{1}}$($L_{\sigma_{2}^{-1}}$), the last (first) $\mathbb{1}$ is introduced due to the chain that is not involved in the braiding. According to (2) we can now utilize the combination of the braiding operations $\sigma_{1}$ and $\sigma_{2}^{-1}$ to realise the braiding matrix $U_{M}$ that corresponds to the desired knots and links shown in Fig. 1a.

Next, we encode the initial logical state $|\phi_{0}\rangle=|0_{1}0_{2}0_{3}\rangle_{l}$ of the three chains in terms of spin states. For convenience, we express the spin states of the system in terms of the eigenvectors $\{|x\rangle,|\bar{x}\rangle\}$, $\{|y\rangle,|\bar{y}\rangle\}$ and $\{|z\rangle,|\bar{z}\rangle\}$ of the corresponding Pauli operators $\sigma^{x}$, $\sigma^{y}$ and $\sigma^{z}$, with eigenvalues $\{1,-1\}$, respectively. Then the transformation is given by: $|0_{1}\rangle_{l}=(|x_{1}x_{2}\rangle+|\bar{x}_{1}\bar{x}_{2}\rangle)/\sqrt{2}$, $|1_{1}\rangle_{l}=(|x_{1}x_{2}\rangle-|\bar{x}_{1}\bar{x}_{2}\rangle)/\sqrt{2}$, $|0_{2}\rangle_{l}=(|x_{4}x_{5}x_{6}\rangle-|\bar{x}_{4}\bar{x}_{5}\bar{x}_{6}\rangle)/\sqrt{2}$, $|1_{2}\rangle_{l}=(|x_{4}x_{5}x_{6}\rangle+|\bar{x}_{4}\bar{x}_{5}\bar{x}_{6}\rangle)/\sqrt{2}$, $|0_{3}\rangle_{l}=(|x_{8}x_{9}x_{10}\rangle-|\bar{x}_{8}\bar{x}_{9}\bar{x}_{10}\rangle)/\sqrt{2}$ and $|1_{3}\rangle_{l}=(|x_{8}x_{9}x_{10}\rangle+|\bar{x}_{8}\bar{x}_{9}\bar{x}_{10}\rangle)/\sqrt{2}$. Thus, the initial logical state $|\phi_{0}\rangle$ is represented in terms of spins $\ket{\phi_{0}^{\prime}}=(|x_{1}x_{2}\rangle+|\bar{x}_{1}\bar{x}_{2}\rangle)/\sqrt{2}\otimes(|x_{4}x_{5}x_{6}\rangle-|\bar{x}_{4}\bar{x}_{5}\bar{x}_{6}\rangle)/\sqrt{2}\otimes(|x_{8}x_{9}x_{10}\rangle-|\bar{x}_{8}\bar{x}_{9}\bar{x}_{10}\rangle)/\sqrt{2}$.

Rather than using single photons to encode the desired spin states we improve the efficiency of our photonic system by employing two-photon correlations [27, 28, 29]. Methods provide the details of two-photon correlations encoding. Fig. 2 shows our experimental setup. A 404 nm continuous wave diode laser is used to pump a type-II periodically poled KTiOPO₄ (PPKTP) to generate the two correlated photons through spontaneous parametric down conversion [30]. One of the two photons is sent to side A, and the other photon is sent to side B. On each side, spatial modes of the photon are manipulated via beam displacers (BDs), which are birefringent crystals capable of spatially separating light beams based on their horizontal or vertical polarization. Here, two types of beam displacers are introduced: BD3 with beams separated by 3 mm and BD6 with beams separated by 6 mm. Additionally, half-wave plates (HWPs) and quarter-wave plates (QWPs) are employed to regulate the polarization of the photon, i.e., exchanging the basis between Pauli operators $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$.

After completing the initial state preparation, we need to perform the braiding evolution according to their corresponding sequences of Hamiltonians. The braiding can be implemented through a series of imaginary-time evolution (ITE) operators (see Supplemental Material for the details), and the whole process for $\sigma_{1}$ and $\sigma_{2}^{-1}$ can be simplified as $e^{-H_{0}t}e^{-H_{3}t}e^{-H_{2}t}e^{-H_{1}t}\ket{\phi_{0}}=e^{-\sigma_{3}^{z}}e^{\sigma_{4}^{x}\sigma_{5}^{x}}e^{-\sigma_{4}^{z}}e^{-\sigma_{3}^{x}\sigma_{4}^{y}}e^{\sigma_{2}^{x}\sigma_{3}^{x}}\ket{\phi_{0}}$ and $e^{-H_{0}t}e^{-H^{\prime}_{1}t}\\ e^{-H^{\prime}_{2}t}e^{-H^{\prime}_{3}t}e^{-H^{\prime}_{4}t}e^{-H^{\prime}_{5}t}\ket{\phi_{0}}=e^{\sigma_{4}^{x}\sigma_{5}^{x}}e^{-\sigma_{7}^{z}}e^{\sigma_{8}^{x}\sigma_{9}^{x}}\\ e^{-\sigma_{8}^{z}}e^{\sigma_{9}^{x}\sigma_{10}^{x}}e^{-\sigma_{7}^{x}\sigma_{8}^{y}}e^{\sigma_{5}^{y}\sigma_{6}^{z}\sigma_{7}^{x}}e^{-\sigma_{4}^{z}}\ket{\phi_{0}}$, respectively. Previously [31, 15], we selected an extended evolution time $t$ in order to diminish the influence of the high-energy excited states (labeled as $\ket{e}$), thereby allowing only the ground states (labeled as $\ket{g}$) with the lowest energy to persist following each evolution step.

Here, to avoid the dissipation of the excited states, we design a new scheme inspired by the concept of quantum cooling [32]. The quantum circuit is shown in Fig. 2a. By inserting a control-not (CNOT) gate in the traditional quantum cooling circuit, the non-dissipative ITE is achieved. This can be realised using a designed polarization-dependent Sagnac interferometer (SI). Methods provide the details of the realisation of the CNOT gate. This process emulates quantum cooling by effectively reducing the excited states to their ground states, resulting in zero photon loss during the ITE evolution. The ratio of reflected to transmitted vertical polarization photons exceeds 500:1, ensuring the accuracy of the non-dissipative ITE evolution. This represents a significant advancement compared to our previous experimental approaches, where approximately half of the photons were discarded in each step of the basis transformation (after $n$ steps, $n=1,2,3...$ the total number of photons will decay to $1/2^{n}$ of the initial photon numbers) and thus making it challenging to extend the evolution to multiple steps. This enhancement significantly improves our ability to conduct experiments involving multiple evolutionary steps, thereby facilitating the experimental simulation of complex physical processes.

For braiding operation $\sigma_{1}$, to realise the non-dissipative evolution according to its spin Hamiltonians, we built five cascaded SIs on side A while keeping the device on side B stationary. For $\sigma^{-1}_{2}$, due to the more complex Hamiltonians, multiple BD3 and BD6 were used alongside the cascaded SIs on side A, which is more experimentally challenging compared with $\sigma_{1}$. (The experimental details for these two evolutions can be found in Supplemental Material).

Upon completion of the specific evolution operation, we need to reconstruct the quantum states during the evolution thus ensuring the effectiveness. The states have $2^{3}$ basis from $\ket{000}_{meas}$ to $\ket{111}_{meas}$, which corresponds to a three-qubit quantum state. We thus project the final state to 64-measurement basis composed of $\ket{0}_{meas},\ket{1}_{meas},\ket{+}=(\ket{0}_{meas}+\ket{1}_{meas})/\sqrt{2},\ket{R}=(\ket{0}_{meas}-i\ket{1}_{meas})/\sqrt{2}$ to reconstruct the density matrix of the quantum state such as $|\phi_{0}^{\prime}\rangle$ according to the standard three-qubit quantum state tomography [33]. Specifically, it requires combining the beams from sides A/B into a single beam. This process involves the use of BDs and HWPs to recombine all the equal-optical path spatial modes on both sides, thus transforming the spatial modes encoded information to polarization information. The relative phase between different states in the interferometer is adjusted by inserting multiple fused silica plates (4×20×0.15 mm, with reflectivity less than 0.2%, not depicted in Fig. 2). The analysis of the polarization information is accomplished through the combination of a QWP, an HWP, and a polarization beam splitter (PBS). (More experimental details for implementing quantum state tomography can be found in Supplemental Material).

In addition to the reconstruction of the quantum states, we also perform quantum process tomography [34] to characterize the braiding operation $\sigma_{1}$. We employ the experimentally intuitive basis $I,X,Y$ and $Z$ representing the identity, $\sigma^{x}$, $\sigma^{y}$ and $\sigma^{z}$, respectively. During the $L_{\sigma_{1}}$ braiding evolution the third chain does not participate so we only need to present the reconstruction of the operator $(II+iXX)/\sqrt{2}$. The experimental results are presented in Fig. 3c (real part) and d (imaginary part). The process fidelity of braiding operation $\sigma_{1}$ is 0.983 $\pm$ 0.014.

Similarly, for braiding operation $\sigma_{2}^{-1}$, we also characterize its performance according to the same procedure in $\sigma_{1}$. After the braiding operation $L_{\sigma_{2}^{-1}}=\mathbb{1}\otimes(\mathbb{1}\otimes\mathbb{1}-i\sigma_{y}\otimes\sigma_{x})/\sqrt{2}$, the initial state $|\phi_{0}\rangle=|000\rangle_{l}$ evolves to the final state $|\phi^{\sigma^{-1}_{2}}_{f}\rangle=(|000\rangle_{l}+|011\rangle_{l})/\sqrt{2}$. The experimental results are presented in Fig. 3e (real part) and f (imaginary part) with a fidelity of 0.964 $\pm$ 0.028. The average fidelity of the final states according to the Hamiltonians during the braiding evolution is 0.978 $\pm$ 0.015. Moreover, the reconstructed process matrix $(II-iYX)/\sqrt{2}$ (the first chain does not participate) is shown in Fig. 3g (real part) and h (imaginary part) with a fidelity of 0.980 $\pm$ 0.015. The experimental results of the state tomography and process tomography demonstrate that we have realised both braiding operations, $\sigma_{1}$ and $\sigma_{2}^{-1}$, successfully with high fidelity for each step of the evolution.

In our all-optical platform the imperfections of the experimental results stem from flawed interference, imprecise wave plates and fluctuation in the count of photons. Among these three factors, imperfect interference is the primary reason for the incomplete alignment between the experimental results and theoretical values. During the experimental realisation of $\sigma_{1}$ the main sources of error are imperfections in the five cascaded Sagnac interferometers on Side A, which are experimentally challenging. Meanwhile, the fidelity of both state and process tomography of ${\sigma_{2}^{-1}}$ have decreased slightly compared to $\sigma_{1}$. This is because of the additional multiple beam splitters on side A of the setup required for implementing the ${\sigma_{2}^{-1}}$ evolution. For instance, when we implement the ITE process for evolving the final state corresponding to $H^{\prime}_{5}$ to the final state corresponding to $H^{\prime}_{4}$, operation $e^{\sigma_{5}^{y}\sigma_{6}^{z}\sigma_{7}^{x}}$ needs to be applied to the final state of $H^{\prime}_{5}$. This basis transformation corresponds to experimentally creating eight beams on side A, and then recombining them to four beams which needs interference processes. This process poses significant experimental challenges, especially when these beams have passed through a Sagnac interferometer before recombining. Despite this technical complexity the fidelity of braiding operation ${\sigma_{2}^{-1}}$ still remains exceptionally high, demonstrating the excellent characteristics of our device and its ability to demonstrate Majorana-based Jones polynomial calculations.

The resulting values of the Jones polynomial $V_{L}(i)$ for the realised links $L$ are shown in Fig. 4. In particular, we obtain that the values of the Jones polynomials for the Hopf link, Trefoil knot, Solomon link, Figure Eight knot and Borromean rings are 0.074 $\pm$ 0.008, 1.000 $\pm$ 0.011, 1.369 $\pm$ 0.013, 0.976 $\pm$ 0.042 and 1.888 $\pm$ 0.058, respectively. All of these experimentally obtained values match well the theoretical predictions, as shown in Fig. 4. Mathematically, even if the Jones polynomial values of two knots are identical, it does not necessarily mean that they are topologically equivalent [35]. Notably, the Trefoil knot and the Figure Eight knot are topologically distinct despite sharing the same value, $|V_{L}(i)|=1$. This limitation underscores the difficulty in discerning specific distinct links from their corresponding Jones polynomials when they are evaluated at a specific value $t$ [23]. Nevertheless, we can clearly distinguish the links whose values of the Jones polynomials, $V_{L}(i)$, are different. Such an advance can greatly contribute in the fields of statistical physics [36], molecular synthesis technology [37] and integrated DNA replication [38], where intricate topological links and knots emerge frequently.

The topological encoding of information with Majorana anyons is equivalent to encoding information in a quantum error correcting code (QECC) [23]. The added value of manipulating this logical information by braiding anyons is that the information remains encoded in the QECC at all times during the performance of the logical gates. This required a minimum size of ten fermionic sites to fault-tolerantly encode the manipulations between three pairs of anyons, as shown in Fig. 1. Hence, the prototype algorithm we performed for evaluating Jones polynomials is encoded in a QECC that protects it against local errors acting on the fermionic sites at all times. We leave the demonstration of fault-tolerance of our encoded quantum algorithm to a future investigation [15]. The versatility of our platform enables the simulation of a broader range of links, paving the way for the implementation of scalable quantum computation based on fault-tolerant Majorana fermions in various platforms, including superconducting circuits [39], integrated silicon chips [40], and Rydberg atoms [41].