# Observing and braiding topological Majorana modes
on programmable quantum simulators

Nikhil Harle  Department of Physics, Yale University, New Haven CT, 06520, USA
IBM Quantum, MIT-IBM Watson AI lab, Cambridge MA, 02142, USA  Oles Shtanko
IBM Quantum, IBM Research – Almaden, San Jose CA, 95120, USA  Ramis Movassagh
IBM Quantum, MIT-IBM Watson AI lab, Cambridge MA, 02142, USA  Google Quantum
AI, Venice Beach, CA, 90291, USA

###### Abstract

Electrons are indivisible elementary particles, yet paradoxically a collection
of them can act as a fraction of a single electron, exhibiting exotic and
useful properties. One such collective excitation, known as a topological
Majorana mode, is naturally stable against perturbations, such as unwanted
local noise, and can thereby robustly store quantum information. As such,
Majorana modes serve as the basic primitive of topological quantum computing,
providing resilience to errors. However, their demonstration on quantum
hardware has remained elusive. Here, we demonstrate a verifiable
identification and braiding of topological Majorana modes using a
superconducting quantum processor as a quantum simulator. By simulating
fermions on a one-dimensional lattice subject to a periodic drive, we confirm
the existence of Majorana modes localized at the edges, and distinguish them
from other trivial modes. To simulate a basic logical operation of topological
quantum computing known as braiding, we propose a non-adiabatic technique,
whose implementation reveals correct braiding statistics in our experiments.
This work could further be used to study topological models of matter using
circuit-based simulations, and shows that long-sought quantum phenomena can be
realized by anyone in cloud-run quantum simulations, whereby accelerating
fundamental discoveries in quantum science and technology.

It is a unique time in the history of science and engineering when we are
witnessing significant advances in the development of fully controllable,
coherent many-body quantum systems that contain dozens to hundreds of qubits
Chow  _et al._ (2021). Quantum simulators hold the promise of exponentially
outperforming classical computers, which would bring about a host of
applications beyond the reach of classical computers. Perhaps the most
promising application of these systems is the simulation of quantum many-body
systems Feynman (1982), which includes topological phases of matter Wen
(2017); Qi and Zhang (2011). In addition to their exotic nature, topological
quantum states are a promising route to fault-tolerant quantum computation
that is based on non-Abelian excitations such as Majorana fermions Aasen  _et
al._ (2016). Majorana fermions are exotic particles: each is its own
antiparticle, unlike an electron being distinct from its antiparticle
(positron). Despite the remarkable progress, the original proposal for the
realization of Majorana-based quantum memories on solid state devices Lutchyn
_et al._ (2010); Oreg  _et al._ (2010); Beenakker (2013) ultimately
encountered difficulties due to disorder and lack of control, as well as the
inability to separate Majorana modes from other trivial zero-energy states Lee
_et al._ (2014); Kayyalha  _et al._ (2020); Valentini  _et al._ (2021); Yu
_et al._ (2021); Saldaña  _et al._ (2021); Wang  _et al._ (2021). At the same
time, quantum simulators may help in this search with their unprecedented
levels of parameter control for a range of topological models Khemani  _et
al._ (2016); Liu  _et al._ (2013); Potter  _et al._ (2016).

Realization of topological phases hosting Majorana modes in bosonic multi-
qubit devices was first envisioned few decades ago Levitov  _et al._ (2001),
with subsequent theoretical developments You  _et al._ (2014); Backens  _et
al._ (2017). Since then signatures of topological modes were detected in
photonic experiments Kitagawa  _et al._ (2012); Cheng  _et al._ (2019); Xiao
_et al._ (2017) and programmable digital quantum information processors Smith
_et al._ (2019a); Tan  _et al._ (2019); Fauseweh and Zhu (2021); Bassman  _et
al._ (2021); Koh  _et al._ (2022); Smith  _et al._ (2019b); Neill  _et al._
(2021). While these devices are limited to non-equilibrium settings they still
are able to exhibit long-lived signatures of topological modes Wiekowski  _et
al._ (2018); Shtanko and Movassagh (2020). Some of these signatures were
analyzed in programmable processors with methods usually tailored to free-
fermionic models Azses  _et al._ (2020); Choo  _et al._ (2018); Zhang  _et
al._ (2022). However, the qualitative study of the properties of these
topological excitations remained a challenge. Braiding of the Majorana
fermions is yet another motivation as it provides the exchange statistics of
the topological excitations and is a necessary step for topological quantum
computation. While there has been progress in manipulation of toy Majorana
modes in photonics Xu  _et al._ (2016, 2018); Liu  _et al._ (2021) and
superconducting architectures Wootton (2017); Zhong  _et al._ (2016); Song
_et al._ (2018); Huang  _et al._ (2021), they were limited to a few qubit
systems, not a real topological phase. Thus, direct probing of the topological
modes and their manipulation remained an open problem.

Using existing noisy quantum hardware, we aim to perform quantitative
simulations of topological quantum matter. We recreate the state of one-
dimensional topological superconductor widely known for hosting a pair of
exotic “half-electron” Majorana modes at its boundaries. We show how to use
Fourier transformation of multi-qubit observables to reliably determine the
structure of Majorana modes. We also demonstrate how the detection of two-
point correlation functions make it possible to distinguish between trivial
and topological modes. Finally, we introduce and implement the fast
approximate swap (FAS): a general non-adiabatic method to approximately braid
Majorana fermions in one dimension. Unlike conventional adiabatic methods, it
allows implementation on the current generation of noisy quantum hardware.

Floquet engineering. Time-periodic (Floquet) systems had proven to be
particularly suitable for simulations on digital quantum processors. In
particular, when system Hamiltonian alternates between two or more local
Hamiltonians being sums of mutually commuting terms, this choice of quantum
dynamics provides a remarkable resource utilization. In this way, unlike
trotterized continuous dynamics, a constant-time Floquet dynamics can be
simulated using constant-depth circuit. While Floquet systems may be compared
in their form to rough trotterization of continuous dynamics, they exhibit a
wide variety of topological phases Rudner and Lindner (2020). The Floquet
topological phase may be quite robust despite the presence of disorder Shtanko
and Movassagh (2018).

Our focus is on the time-periodic Hamiltonian

$H(t)=\sum_{j=1}^{N-1}\Bigl{(}J(t)\,X_{j}X_{j+1}+\lambda(t)\,Z_{j}Z_{j+1}\Bigl{)}+h(t)\sum_{j=1}^{N}Z_{j},$
(1)

where $X_{j}$ and $Z_{j}$ are single-qubit Pauli operators,
$\\{J,\lambda,h\\}(t+T)=\\{J,\lambda,h\\}(t)$ is a set of time-periodic
parameters, $T$ is the time period, $N$ is the number of qubits.

We propose a protocol that divides a single driving period into three parts.
For simplicity, we consider the driving period acting from $t=0$ to $t=T$.
During the first part, from the start of the period to time $\tau_{1}$, we set
$h(t)=h$ and the other coefficients to zero, $J(t)=\lambda(t)=0$. Next, for
times in between $\tau_{1}$ and $\tau_{2}$, we set $J(t)=J$ and the rest of
the coefficients to zero. Lastly, between $\tau_{2}$ and the end of the period
$T$, we set the last term to be on, $\lambda(t)=\lambda$, and all other terms
to zero. Therefore, only one term in the Hamiltonian in Eq. (1) is active at
any given moment.

A quantum circuit can reproduce such a quantum dynamics protocol at discrete
times $t_{n}=nT$. At such times, the system’s state is described by the
wavefunction $|\psi_{n}\rangle=U_{F}^{n}|\psi_{0}\rangle$, here
$|\psi_{0}\rangle$ is the initial state and
$U_{F}=\exp(-i\int_{0}^{t}H(t^{\prime})dt^{\prime})$ is the Floquet unitary,

$\begin{split}\quad U_{F}=\prod_{j=1}^{N-1}e^{-i\varphi
Z_{j}Z_{j+1}}&\prod_{j=1}^{N-1}e^{-i\theta
X_{j}X_{j+1}}\prod_{j=1}^{N}e^{-i\phi Z_{j}},\end{split}$ (2)

where the gate angles are $\phi=h\tau_{1}$, $\theta=J(\tau_{2}-\tau_{1})$, and
$\varphi=\lambda(T-\tau_{2})$. The corresponding experimental protocol that
involves local single- and two-qubit gates is depicted in Fig. 1(a), where
each cycle corresponds to a single Floquet unitary.

Figure 1: Circuit and phase diagram. (a) Schematics for an 8-qubit circuit
including the initialization, evolution, and measurement parts. The
initialization process involves single-qubit Hadamard gates (yellow, #1).
Evolution is composed of cycles consisting of $Z$-gates (green, #2),
$XX$-gates (blue, #3), and $ZZ$-gates (red, #4). Measurement provides the
expectation of operators $\gamma_{\mu}^{L,R}$ or $\gamma_{\mu}\gamma_{\nu}$
(see SI Section 1). (b) Phase diagram for $\lambda=0$, depicting four possible
phases, see text. (c) Experimentally measured Fourier component
$|F^{L}_{1}(\omega)|$ as a function of $\phi$ for fixed $\theta=\pi/8$ using a
21-qubit system, implemented on ibm_hanoi. The system exhibits transitions
from MZM to trivial phase and from trivial to MPM phase. Detected peaks
indicate the presence of Majorana modes at frequencies $\omega=0$ and
$\omega=\pi$.

The model has received considerable attention in the study of condensed matter
systems due to its alternative description in terms of spinless fermions. By
Jordan-Wigner transformation, the qubit Pauli operators can be transformed
into nonlocal Majorana fermion operators $\gamma_{\mu}$ satisfying
$\\{\gamma_{\mu},\gamma_{\nu}\\}=2\delta_{\mu\nu}$, where $\mu,\nu=1,\dots,2N$
Suzuki  _et al._ (2012). It is not a unique mapping; here we use two
equivalent Jordan-Wigner representations, denoted as $\gamma^{L,R}_{\mu}$ and
associated with the right and left boundaries. In these representations a
Majorana operator becomes a string of Pauli operators connected to one of the
boundaries. As we show in Section 1 of Supplementary Information (SI), the
expectation values of these operators can be obtained from single-qubit
measurements preceded by a series of two-qubit gates. We will not include the
superscripts for the Majorana operators when the choice of representation is
not important.

Figure 2: Detection of Majorana modes. Panels (a)-(c) show the absolute value
of the experimentally observed Majorana mode wavefunctions $|\psi^{s}_{\mu}|$
(dots) in comparison with its theoretical prediction (lines) for $N=10$
qubits. Wavefunctions are further normalized because under noise effects Eq.
(3) is inexact. Bottom panel illustrates the density function
$g(x,\omega)=|F^{L}_{2x-1}(\omega)|^{2}+|F^{R}_{2x-1}(\omega)|^{2}$, the
bright peaks show the frequency of the modes. (a) MZM extracted using
ibm_montreal device in the topological phase $\theta=\pi/4$, $\phi=\pi/8$, and
$\varphi=0$, using $D=11$ cycles (b) MPM extracted using the same device in
topological phase $\theta=\pi/4$, $\phi=3\pi/8$, and $\varphi=0$, using $D=11$
cycles (c) MZM wavefunction extracted using ibm_mumbai device for interacting
topological phase $\theta=\pi/4$, $\phi=\pi/16$, and $\varphi=\pi/16$, using
$D=21$ cycles. (d) Difference between trivial phase $\theta=\pi/16$ and
$\phi=\pi/4$ with two trivial boundary modes (blue, circles) and topological
phase $\theta=\pi/4$ and $\phi=\pi/16$ (red, squares) quantified by
$|T_{1,2x}|$ in Eq. (5), measured using ibm_toronto device. The result is
calculated as the average of 10 random initial states and $D=11$ cycles. The
error bars are one standard deviation. The expectation values used to generate
all figures are calculated by averaging over 8192 circuit runs.

In the case that $\lambda=0$, the Hamiltonian in Eq. (1) is non-interacting
and takes the simple quadratic form
$H(t)=\sum_{\mu,\nu=1}^{2N}h_{\mu\nu}(t)\gamma_{\mu}\gamma_{\nu}$, where
$h_{\mu\nu}$ is an antisymmetric Hermitian matrix. Due to its free fermionic
nature, dynamics generated by such a Hamiltonian are classically efficient to
simulate (see SI Section 2). In this regime, depending on the ratio between
$J$ and $h$, the system exhibits various phases including the symmetry-
protected topological phases Khemani  _et al._ (2016); Potter  _et al._
(2016), as summarized by the phase diagram shown in Fig. 1(b). Among these
four phases, there is one (shown in white) that is trivial and topologically
equivalent to a product state. There are three more topological phases. The
first phase is topologically equivalent to the static Kitaev chain (blue).
Under open boundary conditions, this phase exhibits two symmetry-protected
modes at zero quasi-energy called Majorana zero modes (MZM). The remaining
topological phases only occur in time-driven systems. For example, the second
phase (red) exhibits a pair of Majorana $\pi$ modes (MPM) occurring at quasi-
energy $\pi$ Liu  _et al._ (2013). The third phase (green) is distinct from
the rest and hosts both MZM and MPM. Majorana modes in non-interacting systems
manifest themselves by the presence of a pair of conserved boundary-localized
operators $\Gamma^{\omega}_{s}$ that satisfy
$U^{\dagger}_{F}\Gamma^{\omega}_{s}U_{F}=e^{-i\omega}\Gamma^{\omega}_{s}$
Jermyn  _et al._ (2014); Alicea and Fendley (2016) where index $s\in\\{L,R\\}$
defines right and left eigenmodes respectively, and $\omega\in\\{0,\pi\\}$. We
will skip the frequency index $\omega$ when the context is clear.

In the interacting case $\varphi\neq 0$, Majorana mode operators are not
conserved across the spectrum, i.e.
$U^{\dagger}_{F}\Gamma^{\omega}_{s}U_{F}-e^{-i\omega}\Gamma^{\omega}_{s}=O(\tau^{-1})$.
As a result, the observables associated with topological modes must decay with
characteristic lifetime $\tau$. As was shown in Ref. Shtanko and Movassagh
(2020), if the bulk has vanishing dispersion, for small interaction angles
$\varphi$ the lifetime diverges as $\tau\propto\mathcal{O}(\exp(c/\varphi))$,
where the constant $c$ depends on the details of interaction. In practice, the
lifetime may exceed dozens of Floquet cycles even if the bulk has finite
dispersion and interactions are not too strong. This approximate conservation
of Majorana modes leads to the persistent signal for some local observables
when the rest reach infinite-temperature values. The primary goal of this work
is to use this long-lived signal to restore the structure of the modes from
the experiment. In this case we look for Majorana modes of the form
$\Gamma_{s}=\sum_{\mu=1}^{2N}\psi^{s}_{\mu}\gamma_{\mu},$ where
$\psi^{s}_{\mu}$ are real-valued wavefunctions. We also develop a method to
distinguish trivial and topological modes.

Finally, we illustrate the exchange of Majorana modes and verify that the
exchange results in the desired change of phase of the wavefunction.
Conventionally, such an exchange is modeled by a slow adiabatic implementation
of the unitary map $\mathcal{E}_{\rm ex}(\cdot)=U^{\dagger}_{\rm
ex}(\cdot)U_{\rm ex}$, where $U_{\rm
ex}=\exp(-\frac{\pi}{4}\Gamma_{L}\Gamma_{R})$. Such a map provides
$\mathcal{E}_{\rm ex}(\Gamma_{R})=\Gamma_{L}$ and $\mathcal{E}_{\rm
ex}(\Gamma_{L})=-\Gamma_{R}$. While it is possible to carry out this procedure
for one-dimensional Floquet systems Bomantara and Gong (2018); Bauer  _et al._
(2019), it might require quantum circuits with depths beyond what is available
on noisy devices. Below we show an alternative way to perform such an exchange
on a noisy quantum hardware.

Majorana wavefunctions. Our first objective is to detect the presence of
Majorana modes and measure the details of their structure using Fourier
transformation Neill  _et al._ (2021); Wiekowski  _et al._ (2018). We assume
that there are no other eigenmodes with zero or $\pi$ frequencies. In this
case, we can use the asymptotic formula (see SI Section 3)

$\begin{split}&\psi^{L}_{\mu}(\omega)=F^{L}_{\mu}(\omega)/\sqrt{F^{L}_{1}(\omega)},\quad\psi^{R}_{\mu}(\omega)=F^{R}_{\mu}(\omega)/\sqrt{F^{R}_{2N}(\omega)},\end{split}$
(3)

where $\omega\in\\{0,\pi\\}$ is the mode frequency, the positivity of
$F^{L}_{1}(\omega)$ and $F^{R}_{2N}(\omega)$ is proven in SI Section 3, and

$F^{s}_{\mu}(\omega)=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{i\omega
n}\langle\psi_{0}|U_{F}^{{\dagger}n}\gamma^{s}_{\mu}U^{n}_{F}|\psi_{0}\rangle,$
(4)

with $|\psi_{0}\rangle=\bigotimes_{i=1}^{N}|+\rangle_{i}$ being the product
state of eigenstates of Pauli operator $X$ with eigenvalue one, and
superscript $s\in\\{L,R\\}$ conforming with the representation of the Majorana
operator. The order in the limit is important: one first takes the limit over
the number of qubits $N$, and then the limit over the number of cycles $D$.

In spite of the fact that the true limit cannot be reached experimentally, we
measure the quantities $F^{s}_{\mu}(\omega)$ approximately using the largest
available $N$ and $D$. The values of $D$ must not exceed the Majorana mode
lifetime $\tau$ such that $D/\tau\ll 1$. First, we initialize the qubits in
the product state $|\psi_{0}\rangle$ and apply an $n$-cycle circuit as shown
in Fig. 1(a) for $n=0,\dots,D-1$. For each circuit, we determine the
expectation of $\gamma^{R,L}_{\mu}$. In the last step, we estimate the
approximate value of $F_{\mu}(\omega)$ by summing up the results for each
$n$-cycle circuit with corresponding Fourier coefficients $e^{i\omega n}/D$.

As an example, Fig. 1(c) shows the function $|F^{L}_{1}(\omega)|$ and its use
in detecting Majorana modes and topological phases. The plots illustrate the
dependence of this function on angle $\phi$ for the fixed $\theta=\pi/8$ and
are similar to differential conductance spectra found in solid-state
experiments Beenakker (2013). The function is equal to the topological mode
density at the boundary, $F^{L}_{1}(0)=(\psi^{L}_{1})^{2}$. In particular, we
observe a strong signal for $\omega=0$ in the topological phase for value
$\phi=0$, as it indicates the presence of the left MZM. Strong peaks also
appear at frequencies $\omega=\pm\pi$ indicating the presence of MPM for
$\phi=\pi/2$. The peaks’ intensities decrease in the bulk for intermediate
angles. For $\theta=\pi/8$ the boundary signal disappears at $\phi=\pi/8$ and
$3\pi/8$ as the system transitions into the trivial phase.

Next, the values of $F^{s}_{\mu}(\omega)$ for $\mu>1$ help us recover Majorana
wavefunctions $\psi^{s}_{\mu}$. Plots in Fig. 2(a)-(c) illustrate the
normalized absolute values of wavefunctions corresponding to MZMs and MPMs in
both non-interacting ($\varphi=0)$ and interacting ($\varphi=\pi/16)$ regimes.
The results for the non-interacting regime are in a good agreement with the
theoretical prediction. In the interacting regime, where we add an extra set
of noisy two-qubit ZZ gates in each Floquet cycle, we expect to see a visibly
higher level of noise in the resulting wavefunction as can be seen in Fig. 2c.
More data assessing the device’s performance is presented in the Supplementary
Information.

Detecting trivial modes. Majorana modes may not be the only modes responsible
for zero-frequency signals Lee  _et al._ (2014); Kayyalha  _et al._ (2020);
Valentini  _et al._ (2021); Yu  _et al._ (2021); Saldaña  _et al._ (2021);
Wang  _et al._ (2021). In this work, we demonstrate that quantum simulators
can be used to distinguish unpaired Majorana zero modes from the other
topologically trivial localized excitations. Topological Majorana $\pi$ modes
can be treated similarly. We use a generalized notation
$\Delta_{k}=\sum_{\nu}\psi_{k\nu}\gamma_{\nu}$ for both zero-frequency trivial
and topological Majorana modes, $[\Delta_{k},U_{F}]=0$, and $\psi_{k\nu}$ are
real wavefunctions that are localized at the boundaries. In contrast to
Majorana modes residing at opposite boundaries, any pair of trivial modes must
always be localized near the same position. Below we assume that the effect of
disorder on the localization of the wavefunction is negligible.

We examine the two-point correlation function ($\omega=0$)

$T_{\mu,\nu}=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}\langle\tilde{\psi}_{0}|U_{F}^{{\dagger}n}\gamma_{\mu}\gamma_{\nu}U^{n}_{F}|\tilde{\psi}_{0}\rangle,$
(5)

where
$|\tilde{\psi}_{0}\rangle=|\psi_{a}\rangle|s_{2}\rangle|s_{3}\rangle\dots|s_{N-2}\rangle|\psi_{a}\rangle$,
where $|s_{i}\rangle$ are random states in $Z$-basis with eigenvalues
$s_{i}=\pm 1$, and $|\psi_{a}\rangle=\cos a|0\rangle+i\sin a|1\rangle$ for
$a\in[0,\pi]$ being a phase. For simplicity, we consider the non-interacting
case $\lambda=0$. Then the value of the correlation function for $\mu=1$ and
$\nu=2$ is (see SI Section 4)

$T_{1,2}=i\cos
2a\lim_{N\to\infty}\sum_{kk^{\prime}}(\psi^{2}_{k1}\psi^{2}_{k^{\prime}2}-\psi^{2}_{k2}\psi^{2}_{k^{\prime}1}).$
(6)

If there is only one pair of topological modes separated by the system size,
then $T_{1,2}=0$. Indeed, in this case
$\sum_{kk^{\prime}}\psi^{2}_{k1}\psi^{2}_{k^{\prime}2}-\psi^{2}_{k2}\psi^{2}_{k^{\prime}1}=(\psi^{R}_{1})^{2}(\psi^{L}_{2})^{2}-(\psi^{R}_{2})^{2}(\psi^{L}_{1})^{2}\propto
O(2^{-\Theta(N)})$. A pair of trivial localized states at the left boundary,
however, would result in $T_{1,2}>0$. At the same time, $T_{1,2N}$ is nonzero
for both cases, while in the middle of the system, i.e. $T_{1,x}=0$ for
$x=2cN$ and $1>c>0$. As a consequence, correlation function indicates the
presence of zero-frequency modes but has a different structure for trivial and
Majorana modes.

Figure 3: Braiding. Comparison of normalized original wavefunction in Eq. (3)
for the left (a) and right (c) modes and braided wavefunction in Eq. (10) for
the left (b) and right (d) modes with the theoretically estimated angle
$\alpha_{0}=0.263127\pi$. We use the 5-qubit system on $ibm\\_hanoi$ device
with the parameters $\phi=\pi/16$, $\theta=\pi/4$, and $\varphi=0$ and maximum
number of cycles $D=11$, averaged over 30 experiments each with 8192 shots.
Error bars are one standard of deviation. Experimental data are represented by
points, whereas theoretical predictions are represented by lines. Plots
illustrate that modes acquire a relative minus sign after braiding
$\tilde{\psi}^{L}_{\mu}=\psi^{R}_{\mu}$ and
$\tilde{\psi}^{R}_{\mu}=-\psi^{L}_{\mu}$.

In order to illustrate this method, we consider two examples of non-
interacting systems, $\lambda=0$. In the first example, we use the Hamiltonian
in the topological phase $(\theta=\pi/4,\phi=\pi/16)$. We compare this case to
a trivial system with a slightly modified Hamiltonian. In particular, we set
to zero the $XX$-term and $Z$-term for the first and the last qubits, thus
decoupling them from the rest of the system (see SI Section 5). The state of
the rest of the qubits is governed by the Floquet evolution in Eq. (2) with
parameters $(\theta=\pi/16,\phi=\pi/4)$. This modification mimics a possible
error when some of the links between the qubits are dysfunctional. The
modification produces two trivial full-electron modes at opposite boundaries,
which is equivalent to four non-topological Majorana modes
$\Delta_{1}=\gamma_{1}$, $\Delta_{2}=\gamma_{2}$, $\Delta_{3}=\gamma_{2N-1}$,
and $\Delta_{4}=\gamma_{2N}$. Using only the observables in Eq. (4), it is
difficult to distinguish between these modes and topological modes. However,
if we measure the sequence $|T_{1,2x}|$ for $x=1,\dots,N$ for a random
configuration of the initial state, it shows an important difference. As shown
in Fig. 2(d), the curve for trivial case is characterized by two peaks at
$x=1$ and $x=N$, while topological system has only one peak around $x=N$.
Thus, observation of a single peak provides reliable evidence distinguishing
topological Majorana modes from the other possible trivial modes.

Braiding Majorana modes. Finally, we introduce a method for braiding the
Majorana modes, which we call Fast Approximate Swap (FAS). Here we examine the
parametrized map

$\mathcal{E}_{\alpha}(\cdot):=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}U_{n\alpha}^{\dagger}(\cdot)\,U_{n\alpha}\;,$
(7)

where $U_{n\alpha}=U_{F}^{n}\exp(-\alpha\gamma_{1}\gamma_{2N})U_{F}^{n}$, and
$\alpha\in[0,\pi]$ is a real parameter. This quantum channel is equivalent to
selecting the unitary $U_{n\alpha}$ for $n=0,\dots D-1$ with uniform
probability $1/D$.

Let us assume that the system is reflection-symmetric such that the localized
modes satisfy $\psi^{L}_{1}=\psi^{R}_{2N}=\xi$ and $\xi^{2}\geq 1/2$. Then, by
setting the angle $\alpha_{0}=\arcsin(1/\sqrt{2}\xi)$, the action of the map
on topological Majorana operators is

$\mathcal{E}_{\alpha_{0}}(\Gamma_{R})=p\Gamma_{L},\quad\mathcal{E}_{\alpha_{0}}(\Gamma_{L})=-p\Gamma_{R},$
(8)

where $p=\sqrt{2\xi^{2}-1}\leq 1$ (see SI Section 5). This procedure
constitutes approximate FAS method of braiding that aims to replace the
conventional adiabatic process. This method applies in both the interacting
and non-interacting regimes.

We also establish the effect of proposed braiding map on Majorana operators in
absence of localization in non-interacting limit $\lambda=0$,

$\begin{split}&\mathcal{E}_{\alpha_{0}}(\gamma_{\mu})=p(\psi^{R}_{\mu}\Gamma_{L}-\psi^{L}_{\mu}\Gamma_{R}).\end{split}$
(9)

This allows us to detect the relative phase of Majorana fermions after
braiding. The braided mode wavefunction can be defined similarly to Eq. (3) as

$\tilde{\psi}^{s}_{\mu}=\frac{1}{\mathcal{N}}\langle\psi_{0}|\mathcal{E}_{\alpha_{0}}(\gamma^{s}_{\mu})|\psi_{0}\rangle,$
(10)

where $\mathcal{N}$ is normalization coefficient. After braiding, assuming
that the system is reflection-symmetric, we expect the braided wavefunctions
to satisfy $\tilde{\psi}^{L}_{\mu}=\psi^{R}_{\mu}$ and
$\tilde{\psi}^{R}_{\mu}=-\psi^{L}_{\mu}$. This behavior is illustrated in Fig.
3, where we compare wavefunctions in Eq. (3) and Eq. (10).

Our braiding procedure depends on the parameter $\alpha_{0}$ that is generally
unknown without prior access to the system. Although here we calculated it
analytically, it may be difficult to find this angle theoretically for generic
Hamiltonians, in which case it would be necessary to rely on experimental
data. For instance, one can evaluate the angle using the measured Majorana
wavefunction. In SI Section 5, we discuss an alternative method of finding the
proper value of $\alpha_{0}$.

Discussions and outlook.In this work, we propose a framework for detecting,
verifying, and braiding Majorana modes on near-term programmable quantum
simulators by employing the Floquet dynamics. This scheme can be generalized
to the continuous evolution of static Hamiltonians by replacing the discrete
Fourier transformation in our work by its continuous version. It would have
been possible to run our experiments on larger qubit devices. However,
Majorana modes exist at the boundaries rather than the bulk, and our current
experiments are sufficient to make conclusive statements about the detection
and braiding of the topological Majorana modes.

The finite lifetime of the Majorana modes is attributed to natural tendency of
Floquet systems to “heat up”. Adding disorder such as randomization of phases
in $Z$ gates, i.e. $\phi\to\phi+\delta_{i}$, where $\delta_{i}\in[-W,W]$, can
reduce the heating because of the many-body localization (MBL) phenomenon
Abanin  _et al._ (2019). However, such a simplistic scheme may require
disorder values $W$ that can cause transition into the trivial phase. Avoiding
phase transition would require finding good model parameters Decker  _et al._
(2020) or using more sophisticated techniques Shtanko and Movassagh (2020).

This work illustrates the power of synthetic near-term qubit-based quantum
computers for demonstrating and studying topological phases of electronic
systems. Indeed, if we neglect noise, the observed dynamics of bosonic system
can perfectly simulate fermionic topological phases if the measurements are
made in the non-local qubit basis. Unlike solid state devices, however, the
Majorana modes in this work are subject to decoherence because noise breaks
the parity symmetry protecting the topological phase in fermionic systems.
This is a serious drawback for using then in topological quantum computation.
Nonetheless, with improvements in coherence time, the role of noise can be
sufficiently reduced as to make this type of quantum simulation useful in
studying topological quantum matter. Solid state systems such as nanowire
devices Lutchyn  _et al._ (2010); Oreg  _et al._ (2010) can be studied through
continuous-time local Hamiltonian simulations. Floquet systems similar to
those studied in this work, in the limit of large frequency, are equivalent to
such simulations. A potential model of a nanowire could incorporate a qubit
“ladder” representing the two spin values (up and down) and local gates that
account for hopping, spin-orbital coupling, and density-density interactions.

This work can be extended to regimes beyond the current classical simulation
capabilities. By using devices with higher connectivity, it is possible to
study generic two-dimensional materials with a broader variety of topological
phases and to explore new possibilities for topological quantum computation.
Further, the method studied for extracting the Majorana modes may be extended
to the study of local integrals of motion in many-body localized systems,
similarly to the proposal in Ref. Chandran  _et al._ (2015).

Acknowledgements. We thank Frank Pollmann and Bela Bauer for helpful
discussions. We would also like to thank Sergey Bravyi, Zlatko Minev, and
Sarah Sheldon for their help in preparing this publication. The research was
partly supported by the IBM Research Frontiers Institute. This research was
also supported in part by the National Science Foundation under Grant No. NSF
PHY-1748958.

Code availability. The code to run the experiment and access the data
presented in this study is publicly available on GitHub using the link:
https://github.com/IBM/observation-majorana.git

## References

  * Chow  _et al._ (2021) J. Chow, O. Dial, and J. Gambetta, “Ibm quantum breaks the 100-qubit processor barrier,”  (2021).
  * Feynman (1982) R. P. Feynman, Int. J. Theor. Phys. 21, 467 (1982).
  * Wen (2017) X.-G. Wen, Rev. Mod. Phys. 89, 041004 (2017).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Aasen  _et al._ (2016) D. Aasen, M. Hell, R. V. Mishmash, A. Higginbotham, J. Danon, M. Leijnse, T. S. Jespersen, J. A. Folk, C. M. Marcus, K. Flensberg, and J. Alicea, Phys. Rev. X 6, 031016 (2016).
  * Lutchyn  _et al._ (2010) R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ (2010) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Beenakker (2013) C. Beenakker, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).
  * Lee  _et al._ (2014) E. J. H. Lee, X. Jiang, M. Houzet, R. Aguado, C. M. Lieber, and S. De Franceschi, Nat. Nanotech. 9, 79 (2014).
  * Kayyalha  _et al._ (2020) M. Kayyalha, D. Xiao, R. Zhang, J. Shin, J. Jiang, F. Wang, Y.-F. Zhao, R. Xiao, L. Zhang, K. M. Fijalkowski, _et al._ , Science 367, 64 (2020).
  * Valentini  _et al._ (2021) M. Valentini, F. Peñaranda, A. Hofmann, M. Brauns, R. Hauschild, P. Krogstrup, P. San-Jose, E. Prada, R. Aguado, and G. Katsaros, Science 373, 82 (2021).
  * Yu  _et al._ (2021) P. Yu, J. Chen, M. Gomanko, G. Badawy, E. Bakkers, K. Zuo, V. Mourik, and S. Frolov, Nat. Phys. 17, 482 (2021).
  * Saldaña  _et al._ (2021) J. C. E. Saldaña, A. Vekris, L. Pavešič, P. Krogstrup, R. Žitko, K. Grove-Rasmussen, and J. Nygård, arXiv:2101.10794 (2021).
  * Wang  _et al._ (2021) D. Wang, J. Wiebe, R. Zhong, G. Gu, and R. Wiesendanger, Phys. Rev. Lett. 126, 076802 (2021).
  * Khemani  _et al._ (2016) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016).
  * Liu  _et al._ (2013) D. E. Liu, A. Levchenko, and H. U. Baranger, Phys. Rev. Lett. 111, 047002 (2013).
  * Potter  _et al._ (2016) A. C. Potter, T. Morimoto, and A. Vishwanath, Phys. Rev. X 6, 041001 (2016).
  * Levitov  _et al._ (2001) L. Levitov, T. Orlando, J. Majer, and J. Mooij, arXiv:cond-mat/0108266 (2001).
  * You  _et al._ (2014) J. You, Z. Wang, W. Zhang, and F. Nori, Sci. Rep. 4, 1 (2014).
  * Backens  _et al._ (2017) S. Backens, A. Shnirman, Y. Makhlin, Y. Gefen, J. E. Mooij, and G. Schön, Phys. Rev. B 96, 195402 (2017).
  * Kitagawa  _et al._ (2012) T. Kitagawa, M. A. Broome, A. Fedrizzi, M. S. Rudner, E. Berg, I. Kassal, A. Aspuru-Guzik, E. Demler, and A. G. White, Nat. Commun. 3, 1 (2012).
  * Cheng  _et al._ (2019) Q. Cheng, Y. Pan, H. Wang, C. Zhang, D. Yu, A. Gover, H. Zhang, T. Li, L. Zhou, and S. Zhu, Phys. Rev. Lett. 122, 173901 (2019).
  * Xiao  _et al._ (2017) L. Xiao, X. Zhan, Z. Bian, K. Wang, X. Zhang, X. Wang, J. Li, K. Mochizuki, D. Kim, N. Kawakami, _et al._ , Nat. Phys. 13, 1117 (2017).
  * Smith  _et al._ (2019a) A. Smith, M. Kim, F. Pollmann, and J. Knolle, npj Quantum Inf. 5, 1 (2019a).
  * Tan  _et al._ (2019) X. Tan, Y. X. Zhao, Q. Liu, G. Xue, H.-F. Yu, Z. D. Wang, and Y. Yu, Phys. Rev. Lett. 122, 010501 (2019).
  * Fauseweh and Zhu (2021) B. Fauseweh and J.-X. Zhu, Quantum Inf. Process. 20, 1 (2021).
  * Bassman  _et al._ (2021) L. Bassman, M. Urbanek, M. Metcalf, J. Carter, A. F. Kemper, and W. A. de Jong, Quantum Sci. Technol. 6, 043002 (2021).
  * Koh  _et al._ (2022) J. M. Koh, T. Tai, Y. H. Phee, W. E. Ng, and C. H. Lee, npj Quantum. Inf. 8, 16 (2022).
  * Smith  _et al._ (2019b) A. Smith, B. Jobst, A. G. Green, and F. Pollmann, arXiv:1910.05351 (2019b).
  * Neill  _et al._ (2021) C. Neill, T. McCourt, X. Mi, Z. Jiang, M. Y. Niu, W. Mruczkiewicz, I. Aleiner, _et al._ , Nature 594, 508 (2021).
  * Wiekowski  _et al._ (2018) A. Wiekowski, M. M. Maśka, and M. Mierzejewski, Phys. Rev. Lett. 120, 040504 (2018).
  * Shtanko and Movassagh (2020) O. Shtanko and R. Movassagh, Phys. Rev. Lett. 125, 086804 (2020).
  * Azses  _et al._ (2020) D. Azses, R. Haenel, Y. Naveh, R. Raussendorf, E. Sela, and E. G. Dalla Torre, Phys. Rev. Lett. 125, 120502 (2020).
  * Choo  _et al._ (2018) K. Choo, C. W. von Keyserlingk, N. Regnault, and T. Neupert, Phys. Rev. Lett. 121, 086808 (2018).
  * Zhang  _et al._ (2022) X. Zhang, W. Jiang, J. Deng, K. Wang, J. Chen, P. Zhang, W. Ren, H. Dong, S. Xu, Y. Gao, _et al._ , Nature 607, 468 (2022).
  * Xu  _et al._ (2016) J.-S. Xu, K. Sun, Y.-J. Han, C.-F. Li, J. K. Pachos, and G.-C. Guo, Nat. Commun. 7, 1 (2016).
  * Xu  _et al._ (2018) J.-S. Xu, K. Sun, J. K. Pachos, Y.-J. Han, C.-F. Li, and G.-C. Guo, Sci. Adv. 4, eaat6533 (2018).
  * Liu  _et al._ (2021) Z.-H. Liu, K. Sun, J. K. Pachos, M. Yang, Y. Meng, Y.-W. Liao, Q. Li, J.-F. Wang, Z.-Y. Luo, Y.-F. He, D.-Y. Huang, G.-R. Ding, J.-S. Xu, Y.-J. Han, C.-F. Li, and G.-C. Guo, PRX Quantum 2, 030323 (2021).
  * Wootton (2017) J. R. Wootton, Quantum Sci. Technol. 2, 015006 (2017).
  * Zhong  _et al._ (2016) Y. P. Zhong, D. Xu, P. Wang, C. Song, Q. J. Guo, W. X. Liu, K. Xu, B. X. Xia, C.-Y. Lu, S. Han, J.-W. Pan, and H. Wang, Phys. Rev. Lett. 117, 110501 (2016).
  * Song  _et al._ (2018) C. Song, D. Xu, P. Zhang, J. Wang, Q. Guo, W. Liu, K. Xu, H. Deng, K. Huang, D. Zheng, S.-B. Zheng, H. Wang, X. Zhu, C.-Y. Lu, and J.-W. Pan, Phys. Rev. Lett. 121, 030502 (2018).
  * Huang  _et al._ (2021) H.-L. Huang, M. Narożniak, F. Liang, Y. Zhao, A. D. Castellano, M. Gong, Y. Wu, S. Wang, J. Lin, Y. Xu, H. Deng, H. Rong, J. P. Dowling, C.-Z. Peng, T. Byrnes, X. Zhu, and J.-W. Pan, Phys. Rev. Lett. 126, 090502 (2021).
  * Rudner and Lindner (2020) M. S. Rudner and N. H. Lindner, Nat. Rev. Phys. 2, 229 (2020).
  * Shtanko and Movassagh (2018) O. Shtanko and R. Movassagh, Phys. Rev. Lett. 121, 126803 (2018).
  * Suzuki  _et al._ (2012) S. Suzuki, J. Inoue, and B. Chakrabarti,  _Quantum Ising Phases and Transitions in Transverse Ising Models_, Lecture Notes in Physics (Springer Berlin Heidelberg, 2012).
  * Jermyn  _et al._ (2014) A. S. Jermyn, R. S. K. Mong, J. Alicea, and P. Fendley, Phys. Rev. B 90, 165106 (2014).
  * Alicea and Fendley (2016) J. Alicea and P. Fendley, Annu. Rev. Condens. Matter Phys. 7, 119 (2016).
  * Bomantara and Gong (2018) R. W. Bomantara and J. Gong, Phys. Rev. Lett. 120, 230405 (2018).
  * Bauer  _et al._ (2019) B. Bauer, T. Pereg-Barnea, T. Karzig, M.-T. Rieder, G. Refael, E. Berg, and Y. Oreg, Phys. Rev. B 100, 041102 (2019).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Decker  _et al._ (2020) K. S. C. Decker, C. Karrasch, J. Eisert, and D. M. Kennes, Phys. Rev. Lett. 124, 190601 (2020).
  * Chandran  _et al._ (2015) A. Chandran, I. H. Kim, G. Vidal, and D. A. Abanin, Phys. Rev. B 91, 085425 (2015).
  * (53) “IBM Quantum Compute Resources,” https://quantum-computing.ibm.com/services.

Supplementary Information for
“Observing and braiding topological Majorana modes
on programmable quantum simulators”

Nikhil Harle,1,2 Oles Shtanko3, and Ramis Movassagh2,4

1Department of Physics, Yale University, New Haven CT, 06520, USA
2IBM Quantum, MIT-IBM Watson AI lab, Cambridge MA, 02142, USA
3IBM Quantum, IBM Research – Almaden, San Jose CA, 95120, USA
4Google Quantum AI, Venice Beach, CA, 90291, USA

##  I Section 1: Hardware setting

We replicate the Floquet dynamics in Eq. (2) on IBM Qiskit using the circuit
in Fig. 1(a). We transpile the circuit on Qiskit using native gates and run it
on the IBM quantum hardware. In particular, for the special value
$\theta=\pi/4$ ($\varphi=\pi/4$), each $XX$-gate ($ZZ$-gate) requires a CNOT
gate in combination with single qubit gates. For other non-zero angle values,
two-qubit gates require two CNOTs in combination with other single-qubit
gates. Therefore, to reduce the depth, part of the experiment is designed to
investigate the case $\theta=\pi/4$.

Figure S1: Configuration of IBM hardware. (a) Layout of the ibm_hanoi device.
Sequence A of 21 qubits was used to generate the frequency-resolved boundary
oscillations shown in Fig. 1(c); sequence B of 10 qubits was used to perform
the braiding experiment shown in Fig. 3. (b) Layout for the ibm_montreal
device. Sequence C of 10 qubits was used to reproduce the Majorana mode
tomography in Fig. 3(a),(b) (C). Layout for the ibm_mumbai device. Sequence D
of 10 qubits used to generate Majorana modes tomography in fig. 2(c). Layout
for the ibm_toronto facility. Sequence D of 10 qubits used to generate the
topological/non-topological mode separation experiment in Fig. 2(d). Figure
S2: Experimental data. The figure shows the raw data for the expectations of
the operators $X_{0}$ (top) and $Z_{0}Z_{1}X_{2}$ (bottom) for noiseless
simulation (dashed black curve), the experimental data (circles, blue), and
rescaled experimental data (squares, orange), where the rescaling takes the
form $\exp(\Gamma n)$, where $n$ is the index of the Floquet cycle and
$\Gamma$ is the compensated decay rate. This rescaling accounts for the decay
of Majorana modes that does not contribute to the measured wavefunction (see
discussion below Eq. (S.47)). Panel (a) shows the results for MZM in Fig. 2(a)
(here, we use $\Gamma=0.0328$ for rescaling), Panel (b) shows the results for
MPM in Fig. 2(b) (using $\Gamma=0.0376$), and Panel (c) shows the results for
MZM in the presence of $ZZ$ gates shown in Fig. 2(c) (using $\Gamma=0.120$) .

For the simulation of quantum systems, Majorana operators must be encoded
using qubits. We use the Jordan-Wigner transformation to implement this
encoding. In particular, we define left representation as
$\gamma^{L}_{2k-1}=\mathcal{Z}^{L}_{k}X_{k}$ and
$\gamma^{L}_{2k}=\mathcal{Z}^{L}_{k}Y_{k}$, where
$\mathcal{Z}^{L}_{k}=\prod_{i=1}^{k-1}(-Z_{i})$ are $Z$-string operators and
$k=1,\dots,N$. This representation is equivalent to the most common
convention. Alternatively, the right representation is
$\gamma^{R}_{2k-1}=\mathcal{Z}^{R}_{k}Y_{k}$ and
$\gamma^{R}_{2k}=-\mathcal{Z}^{R}_{k}X_{k}$, where
$\mathcal{Z}^{R}_{k}=\prod_{i=k+1}^{N}(-Z_{i})$. A more traditional approach
for accessing these operators experimentally is to measure each qubit inside
the string in the basis ($X$, $Y$, or $Z$), returning $\pm 1$ values for each
qubit. In this case, we can use the product of the obtained results as the
measured value. The disadvantage of the conventional method is that it has
lower precision due to accumulated measurement errors. Therefore, we adopt a
scheme that involves only one measurement, as described below.

In order to measure $\gamma^{L}_{\mu}$, we use the expressions

$\begin{split}U^{\dagger}_{k}X_{1}U_{k}=(-1)^{k-1}\mathcal{Z}^{L}_{k}X_{k},\qquad
V^{\dagger}_{k}Y_{1}V_{k}=(-1)^{k-1}\mathcal{Z}^{L}_{k}Y_{k},\end{split}$
(S.1)

where

$U_{k}:=U^{YX}_{1}\dots U^{YX}_{k-1},\qquad V_{k}:=U^{XY}_{1}\dots
U^{XY}_{k-1},$ (S.2)

and can be expressed as a product of two-qubit unitaries
$U^{YX}_{j}:=\exp(-i\frac{\pi}{4}Y_{j}X_{j+1})$ and
$U^{XY}_{j}:=\exp(i\frac{\pi}{4}X_{j}Y_{j+1})$. According to these
expressions, to measure the string operator
$\gamma^{L}_{2k-1}=\mathcal{Z}^{L}_{k}X_{k}$, we apply the gates $U^{YX}_{j}$
consecutively and in the reverse order for $j=k-1,\dots,1$ as shown in Fig.
1(a). Next, we measure the first qubit ($j=1$) in X-basis. Similarly, to
measure the operator $\gamma^{L}_{2k}=\mathcal{Z}^{L}_{k}Y_{k}$, we perform
similar gate sequence but with unitaries $U^{XY}_{j}$ and measure the qubit
$j=1$ in $Y$-basis. Finally, the measurement of $\gamma^{R}_{\mu}$ can be done
by mirroring the entire circuit in Fig. 1(a) upside down.

Next, to evaluate $T_{1,2k}$ defined in Eq. (5) we need to probe the operator
$i\gamma_{1}\gamma_{2k}=(-1)^{k}Y_{1}Z_{2}\dots Z_{k-1}Y_{k}$. In order to
measure the Pauli string $Y_{1}Z_{2}\dots Z_{k-1}Y_{k}$, we note that

$Y_{1}Z_{2}\dots
Z_{k-1}Y_{k}=(-1)^{k-1}G_{1}^{\dagger}\mathcal{Z}^{L}_{k}Y_{k}G_{1}$ (S.3)

where $G=HS^{\dagger}$, where $S={\rm diag}\\{0,i\\}$ is $S$-gate and $H$ the
Hadamard gate. Thus, the procedure of measuring this operator is the same as
$\mathcal{Z}^{L}_{k}Y_{k}$ with the difference that in the latter we apply
$G_{1}$ before the series of $U^{XY}$ gates.

Finally, applying the unitary $\exp\left(-\alpha\gamma_{1}\gamma_{2N}\right)$,
which is necessary for generating $U_{n\alpha}$ in Eq. (7), can be implemented
in a similar manner. The implementing circuit consists of the series of gates
$U^{XY}_{j}$ for $j=2,\dots,N-1$, followed by $\exp\left(-i\alpha
Y_{1}Y_{2}\right)$, then by the series of $U^{XY{\dagger}}_{j}$ in reverse
order.

Results are obtained using IBM quantum hardware ibm . The experiments are
performed on four different 27-qubit devices: ibm_hanoi, ibm_montreal,
ibm_mumbai, and ibm_toronto. The number of qubits utilized for each experiment
vary; Fig. S1 shows the chosen subsets. For example, we perform experiments
represented by Fig. 2 using 10 qubits as the smallest system size that
exhibits an overlap between unpaired Majorana modes which is smaller than the
effect of noise. In contrast, braiding experiments are performed on 5 qubits
as the effect of noise is stronger due for deeper circuits. The depth of the
circuits are chosen to be between 11-21 cycles for most experiments.

##  II Section 2: eigenmodes and their properties

Analyzing the eigenmodes of a system can be an effective way to describe its
dynamics. In this section, we introduce and analyze the properties of
eigenmodes for time-periodic (Floquet) dynamics generated by the unitary
operator $U_{F}$ in Eq. (2). We prove three propositions concerning the
eigenmodes of the system, which we then use in the following sections.

######  Proposition 1.

(Existence of eigenmodes) For any unitary $U_{F}$, there exists a complete set
of eigenoperators (eigenmodes) $\Delta_{b}$ and real eigenfrequencies
$\omega_{b}$, $1\leq b\leq 4^{N}$, such that

$\begin{split}&U_{F}^{\dagger}\Delta_{b}U_{F}=e^{-i\omega_{b}}\Delta_{b},\\\
&{\rm
Tr\,}(\Delta_{b}^{\dagger}\Delta_{b^{\prime}})=2^{N}\delta_{bb^{\prime}},\end{split}$
(S.4)

where $N$ is the number of qubits and $\delta_{bb^{\prime}}$ is the Kronecker
delta.

###### Proof.

Consider $\mathcal{P}=\\{P_{\alpha}:\alpha=1,\dots,d\\}$ as a complete set of
$2^{N}\times 2^{N}$ Pauli basis operators, where $d=4^{N}$ is the dimension of
the space they span. The unitary transformation of each basis operator is
given by

$U^{\dagger}_{F}P_{\alpha}U_{F}=\sum_{\beta=1}^{d}E_{\alpha\beta}P_{\beta},$
(S.5)

where $E_{\alpha\beta}:=2^{-N}{\rm
Tr\,}(U_{F}^{\dagger}P_{\alpha}U_{F}P_{\beta})$ are matrix elements of a real
orthogonal matrix,

$E_{\alpha\beta}=E^{*}_{\alpha\beta},\quad EE^{T}=E^{T}E=I.$ (S.6)

To prove that $E$ is orthogonal, we use Eq. (S.5) to express

$\begin{split}U^{\dagger}_{F}P_{\alpha}P_{\alpha^{\prime}}U_{F}=\sum_{\beta,\beta^{\prime}=1}^{d}E_{\alpha\beta}E_{\alpha^{\prime}\beta^{\prime}}P_{\beta}P_{\beta^{\prime}}.\end{split}$
(S.7)

Next, we take the trace for the two sides of Eq. (S.7) and use the
orthogonality of Pauli matrices, ${\rm
Tr\,}(P_{\alpha}P_{\beta})=2^{N}\delta_{\alpha\beta}$. Using the fact that
$E_{\alpha\beta}=E^{*}_{\alpha\beta}$, we get that

$\sum_{\beta=1}^{d}E_{\alpha\beta}E^{T}_{\beta\alpha^{\prime}}=\sum_{\beta=1}^{d}E_{\alpha^{\prime}\beta}E^{T}_{\beta\alpha}=\delta_{\alpha\alpha^{\prime}}.$
(S.8)

Since the matrix $E$ is orthogonal, it has an orthonormal set of eigenstates
and eigenvalues that lie on the unit circle in the complex plane:

$v^{T}_{b}E=e^{-i\omega_{b}}v^{T}_{b}.$ (S.9)

Using this set, we construct the operators

$\Delta_{b}=\sum_{\alpha=1}^{d}v_{b\alpha}P_{\alpha},$ (S.10)

where $v_{b\alpha}$ stands for the $\alpha$-th entry of $v_{b}$. These
operators satisfy

$\begin{split}U_{F}^{\dagger}\Delta_{b}U_{F}&=\sum_{\alpha=1}^{d}v_{b\alpha}U_{F}^{\dagger}P_{\alpha}U_{F}=\sum_{\alpha,\beta=1}^{d}v_{b\alpha}E_{\alpha\beta}P_{\beta}\\\
&=e^{-i\omega_{b}}\sum_{\beta=1}^{d}v_{b\beta}P_{\beta}=e^{-i\omega_{b}}\Delta_{b},\end{split}$
(S.11)

which proves the first equation in Eq. (S.4). Furthermore, the orthogonality
of the eigenstates of the operator $E$ leads to the expression

$\begin{split}{\rm
Tr\,}(\Delta^{\dagger}_{b}\Delta_{b^{\prime}})&=\sum_{\alpha,\beta=1}^{d}v^{*}_{b\alpha}v_{b\beta}{\rm
Tr\,}(P_{\alpha}P_{\beta})\\\
&=2^{N}\sum_{\alpha=1}^{d}v^{*}_{b\alpha}v_{b^{\prime}\alpha}=2^{N}\delta_{bb^{\prime}},\end{split}$
(S.12)

proving the second equation in Eq. (S.4) and thus completing our proof.∎

The spectrum of the system can contain degeneracies. If we define $B$ as the
full set of eigenmodes, then $B_{\omega}=\\{b_{1},\dots,b_{m}\\}$ is defined
as the subset of eigenmodes with the same frequency $\omega_{b_{k}}=\omega$
for $1\leq k\leq m$ in the limit $N\to\infty$. We define
$\overline{B_{\omega}}=B\setminus B_{\omega}$ as the set of eigenmodes with
frequencies different from $\omega$.

######  Proposition 2.

For the dynamics in Eq. (2), non-trivial eigenmodes $\Delta_{b}\neq I$ and
Majorana operators satisfy

$\Delta_{b}=\sum_{\mu=1}^{2N}w_{\mu
b}\gamma_{\mu}+C_{\perp},\quad\gamma_{\mu}=\sum_{b\in B}w^{*}_{\mu
b}\Delta_{b},$ (S.13)

where $w_{\mu b}$ are complex-valued coefficients and $C_{\perp}$ is an
operator that satisfies ${\rm Tr\,}(\gamma_{\mu}C_{\perp})=0$ for all $\mu$.

###### Proof.

First, we prove the second equation. Since $\Delta_{b}$ form a complete basis,
we have

$O=\frac{1}{2^{N}}\sum_{b\in B}{\rm Tr\,}(O\Delta^{\dagger}_{b})\Delta_{b}.$
(S.14)

Then, the second part of Eq. (S.13) follows from this expression by setting
$w^{*}_{\mu b}:=2^{-N}{\rm Tr\,}(\gamma_{\mu}\Delta^{\dagger}_{b})$. Any many-
body operator can be written as a decomposition of Majorana operators,

$\Delta_{b}=\sum_{\mu=1}^{2N}C^{(1)}_{b,\mu}\gamma_{\mu}+\sum_{\mu\nu=1}^{2N}C^{(2)}_{b,\mu\nu}\gamma_{\mu}\gamma_{\nu}+\dots\;.$
(S.15)

Multiplying both sides by the operator $\gamma_{\mu}$ and taking the trace, we
get

${\rm Tr\,}(\Delta_{b}\gamma_{\mu})=2^{N}C^{(1)}_{b,\mu}.$ (S.16)

That is, we derive the expression

$C^{(1)}_{b,\mu}=w_{\mu b}.$ (S.17)

which leads to Eq. (S.13) and completes our proof. ∎

As a useful tool for analytical derivations, we introduce the Fourier channel

$\mathcal{F}_{\omega}(\cdot):=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{i\omega
n}U_{F}^{{\dagger}n}(\cdot)U^{n}_{F},$ (S.18)

where we first take $N$ and then $D$ to infinity. The action of this channel
is expressed in terms of the eigenmodes of the system by the following
Proposition.

######  Proposition 3.

The action of the Fourier map in Eq. (S.18) can be expressed as

$\mathcal{F}_{\omega}(O)=\lim_{N\to\infty}\frac{1}{2^{N}}\sum_{b\in
B_{\omega}}{\rm Tr\,}(O\Delta^{\dagger}_{b})\Delta_{b}.$ (S.19)

###### Proof.

Using the decomposition in Eq. (S.14), we have

$\mathcal{F}_{\omega}(O)=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{i\omega
n}\sum_{b\in B}\frac{1}{2^{N}}{\rm
Tr\,}(O\Delta^{\dagger}_{b})U_{F}^{{\dagger}n}\Delta_{b}U^{n}_{F}.$ (S.20)

Now, using Proposition 1, we have

$\begin{split}\mathcal{F}_{\omega}(O)=\lim_{N,D\to\infty}\frac{1}{2^{N}D}\sum_{n=0}^{D-1}\sum_{b\in
B}e^{i(\omega-\omega_{b})n}{\rm
Tr\,}(O\Delta^{\dagger}_{b})\Delta_{b}.\end{split}$ (S.21)

Using the property

$\lim_{D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{i\omega n}=\delta_{\omega,0},$
(S.22)

we get the expression

$\mathcal{F}_{\omega}(O)=\lim_{N\to\infty}\frac{1}{2^{N}}\sum_{b\in
B}\delta_{\omega-\omega_{b},0}{\rm Tr\,}(O\Delta^{\dagger}_{b})\Delta_{b}.$
(S.23)

Next, using the definition of $B_{\omega}$, we obtain the statement of this
Proposition. ∎

If the system is non-interacting, i.e. $\lambda=0$, the eigenmodes can be
efficiently expressed in the free fermionic representation. The action of the
Floquet unitary on the Majorana operators is then the linear map

$U^{\dagger}_{F}\gamma_{\mu}U_{F}=\sum_{\nu=1}^{2N}u_{\mu\nu}\gamma_{\nu},$
(S.24)

where $u$ is a $2N\times 2N$ unitary matrix

$u=\exp(4\theta h_{xx})\exp(4\phi h_{z}),$ (S.25)

and $h_{z}$ and $h_{xx}$ are

$\begin{split}&h_{z}=-\frac{1}{2}\sum_{k=0}^{N-1}\Bigl{(}|2k+1\rangle\langle
2k+2|-{\rm h.c.}\Bigl{)},\\\
&h_{xx}=\frac{1}{2}\sum_{k=0}^{N-2}\Bigl{(}|2k+2\rangle\langle 2k+3|-{\rm
h.c.}\Bigl{)}.\end{split}$ (S.26)

Using the solution of the eigenproblem

$\psi^{T}_{k}u=e^{-i\omega^{0}_{k}}\psi^{T}_{k},$ (S.27)

we define the set of single-fermion modes

$\Delta^{0}_{k}=\sum_{\mu=1}^{2N}\psi_{k\mu}\gamma_{\mu},\quad\\{\Delta^{0{\dagger}}_{k},\Delta^{0}_{k^{\prime}}\\}=2\delta_{kk^{\prime}}.$
(S.28)

Since the unitary matrix in Eq. (S.25) is real-valued, the complex conjugate
for both sides in Eq. (S.27) gives us another, opposite-frequency solution

$(\psi^{*}_{k})^{T}u=e^{i\omega^{0}_{k}}(\psi^{*}_{k})^{T}.$ (S.29)

This conclusion proves that the spectrum of the problem is zero-symmetric,
i.e. for each mode $k$ there exists an orthogonal mode $\sigma(k)$ such that
$\psi_{\sigma(k)}=\psi_{k}^{*}$ and $\omega_{\sigma(k)}^{0}=-\omega_{k}^{0}$.
This mathematical property is a reflection of the physical particle-hole
symmetry.

For the non-interacting case, each many-body eigenmode in Eq. (S.4) and its
frequency can be expressed as

$\Delta_{b}=\Delta^{0}_{m_{1}}\dots\Delta^{0}_{m_{K}},\qquad\omega_{b}=\sum_{i=1}^{K}\omega^{0}_{m_{i}}.$
(S.30)

Each mode appears only once in the product, i.e. $m_{i}\neq m_{j}$, and the
number of single-particle modes involved is $K\in\\{1,\dots,2N\\}$. A crucial
advantage of the free fermion representation is its efficiency in evaluating
the eigenmodes on a classical computer by diagonalizing the unitary operator
in Eq. (S.25).

Finally, we define the sets of approximate modes.

######  Definition 1.

An orthonormal operator set $\\{\Delta^{\prime}_{b}\\}$ is called an
$\epsilon$-approximate set of eigenmodes if it satisfies

$\begin{split}&U_{F}^{\dagger}\Delta^{\prime}_{b}U_{F}=e^{-i\omega_{b}}\Delta^{\prime}_{b}+\delta
O,\quad\|\delta O\|\leq\epsilon,\\\ \end{split}$ (S.31)

where $\|\cdot\|$ is the operator norm.

Approximate modes behave similarly to regular eigenmodes up to times
$\tau\propto 1/\epsilon$.

Let us consider the realistic case $D<\infty$ and in the presence of noise. In
the experiment, we implement the map

$\tilde{\mathcal{F}}_{D,\omega}(\cdot):=\frac{1}{D}\sum_{m=0}^{D-1}e^{i\omega
n}\mathcal{E}^{n}_{F}(\cdot),$ (S.32)

where $\mathcal{E}_{F}$ is a general map that includes unitary evolution and
the effect of noise. The effect of this map can be approximated by

$\mathcal{E}_{F}(\Delta_{b})\approx e^{i\omega_{b}-\Gamma_{b}}\Delta_{b}$
(S.33)

where $\Gamma_{b}$ is the decay rate that includes both the effect of the
noise and the finite lifetime of the approximate integral of motion
$\Delta_{b}$. This effect is visible in Fig. S2, especially in the presence of
$ZZ$-gates. Then, the effect of noise and finite depth can be understood as

$\tilde{\mathcal{F}}_{D,\omega}(O)\approx\sum_{b\in
B}f_{D}(\omega-\omega_{b},\Gamma_{b}){\rm
Tr\,}(O\Delta^{\dagger}_{b})\Delta_{b}.$ (S.34)

where the coefficients are

$f_{D}(\omega,\Gamma)=\frac{1}{D}\frac{1-e^{i\omega D-\Gamma
D}}{1-e^{i\omega-\Gamma}}.$ (S.35)

This function generates the peak broadening and effective attenuation factor
as $|f_{D}(\omega,\Gamma)|\leq 1$.

##  III Section 3: Majorana eigenmodes

The topological Majorana modes $\Gamma_{s}$, $s\in\\{L,R\\}$, satisfy

$U^{\dagger}_{F}\Gamma_{s}U_{F}=e^{-i\omega_{M}}\Gamma_{s}+O(1/\tau),$ (S.36)

where $\omega_{M}$ are the Majorana mode frequencies,
$\omega_{M}\in\\{0,\pi\\}$ in the limit $N\to\infty$, and $\tau$ is the
lifetime. The Majorana modes are commonly described as “strong” if
$\tau\to\infty$, while the case $\tau<\infty$ corresponds to “weak” modes.

For the topological Majorana modes, both strong and weak, we use the following
form

$\Gamma_{s}=\sum_{\mu=1}^{2N}\psi^{s}_{\mu}\gamma_{\mu},$ (S.37)

where $\gamma_{\mu}$ are physical Majorana operators and $\psi^{s}_{\mu}$ are
real Majorana wave functions. In the non-interacting topological phase,
$\lambda=0$, there always exists a pair of solutions for the unitary in Eq.
(S.25) which satisfies

$\sum_{\mu=1}^{2N}u_{\mu\nu}\psi^{s}_{\mu}=\pm\psi^{s}_{\nu}$ (S.38)

corresponding to MZMs (plus sign) and MPMs (minus sign). Therefore, non-
interacting regimes are always characterized by strong Majorana modes. On the
other hand, in the interacting regime, i.e. $|\lambda|>0$, the existence of
such operators is not guaranteed. However, for certain values of the angle
$(\theta,\phi)$ such modes are weak eigenmodes characterized by large
$\tau\propto\exp(c/\lambda)\gg 1$ Shtanko and Movassagh (2020). Even in the
ideal case where we ignore the noise in the actual experiment, as will become
clear in Eq. (S.42), the lifetime $\tau$ limits the maximum depth $D$ such
that $D/\tau\to 0$. In practice, this connection means that the Fourier
components evaluated in Eq. (4) must be compared with the experimental results
obtained for $\tau\gg D\gg 1$.

For both strong and weak modes, the wavefunctions $\psi^{s}_{\mu}$ must be
properly normalized, as follows from the condition $\Gamma_{s}^{2}=1$. Indeed,

$\Gamma_{s}^{2}=\sum_{\nu,\nu^{\prime}=1}^{2N}\psi^{s}_{\nu}\psi^{s}_{\nu^{\prime}}\gamma_{\nu}\gamma_{\nu^{\prime}}=\sum_{\nu=1}^{2N}(\psi^{s}_{\nu})^{2}=1,$
(S.39)

where we used the statistics of the Majorana operators,
$\\{\gamma_{\mu},\gamma_{\nu}\\}=2\delta_{\mu\nu}$.

Next, we analyze how to detect Majorana modes using the Fourier components in
Eq. (4) generated from the experiment. To do this, we express these components
using the Fourier map in Eq. (S.18) as

$\begin{split}F^{s}_{\mu}(0):=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}\langle\psi_{0}|&U_{F}^{{\dagger}n}\gamma_{\mu}U^{n}_{F}|\psi_{0}\rangle\\\
&=\langle\psi_{0}|\mathcal{F}_{0}(\gamma_{\mu})|\psi_{0}\rangle.\end{split}$
(S.40)

In the following section we will show how to distinguish true MZM modes from
other trivial edge oscillations. For now, let us focus on zero-frequency MZMs
and assume that they are the only pair of zero-frequency modes that have a
non-zero overlap with the Majorana operators $\gamma_{\mu}$. Then, from
Proposition 3, this action is given by

$\begin{split}\mathcal{F}_{0}(\gamma_{\mu})&=\sum_{s\in\\{L,R\\}}\psi^{s}_{\mu}\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}U_{F}^{{\dagger}n}\Gamma_{s}U^{n}_{F}\\\
&\quad+\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}\frac{1}{2^{N}}\sum_{b\in\overline{B_{0}}}w^{*}_{b\mu}e^{-i\omega_{b}n}\Delta_{b}.\end{split}$
(S.41)

where $\overline{B_{0}}$ is the complementary to zero frequency subspace
$B_{0}$. Under the limit over $D$ the last term of this expression vanishes.
At the same time, the first term can be simplified with Eq. (S.36) as

$\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}U_{F}^{{\dagger}n}\Gamma_{s}U^{n}_{F}=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}(\Gamma_{s}+O(n/\tau))$
(S.42)

Assuming that $\lim_{D\to\infty}D/\tau=0$, we get

$\begin{split}\mathcal{F}_{0}(\gamma_{\mu})=\lim_{N\to\infty}\sum_{s\in\\{L,R\\}}\psi^{s}_{\mu}\Gamma_{s}\equiv\lim_{N\to\infty}(\psi^{L}_{\mu}\Gamma_{L}+\psi^{R}_{\mu}\Gamma_{R}),\end{split}$
(S.43)

Using Eqs. (S.43) and (S.37), we transform Eq. (S.40) into

$F^{\alpha}_{\mu}(0)=\lim_{N\to\infty}\sum_{s\in\\{L,R\\}}\sum_{\mu^{\prime}=1}^{2N}\psi^{s}_{\mu}\psi^{s}_{\mu^{\prime}}\langle\psi_{0}|\gamma^{\alpha}_{\mu^{\prime}}|\psi_{0}\rangle.$
(S.44)

where we restored the Pauli representation index for Majorana operators. Using
$\langle\psi_{0}|\gamma^{L}_{\mu}|\psi_{0}\rangle=\delta_{\mu,1}$ and
$\langle\psi_{0}|\gamma^{R}_{\mu}|\psi_{0}\rangle=\delta_{\mu,2N}$, we rewrite

$\begin{split}&F_{\mu}^{L}(0)=\lim_{N\to\infty}(\psi^{L}_{\mu}\psi^{L}_{1}+\psi^{R}_{\mu}\psi^{R}_{1}),\\\
&F^{R}_{\mu}(0)=\lim_{N\to\infty}(\psi^{L}_{\mu}\psi^{L}_{2N}+\psi^{R}_{\mu}\psi^{R}_{2N}).\end{split}$
(S.45)

Taking into account that the eigenmodes are exponentially suppressed away from
the respective boundaries of the system, i.e.
$\psi^{R}_{1}\sim\psi^{L}_{2N}\sim 2^{-\Theta(N)}$, we obtain the expressions

$\begin{split}&F_{\mu}^{L}(0)=\lim_{N\to\infty}\psi^{L}_{\mu}\psi^{L}_{1},\\\
&F^{R}_{\mu}(0)=\lim_{N\to\infty}\psi^{R}_{\mu}\psi^{R}_{2N},\end{split}$
(S.46)

from which we derive

$\begin{split}&\psi^{L}_{\mu}=F^{L}_{\mu}(0)/\sqrt{F^{L}_{1}(0)},\\\
&\psi^{R}_{\mu}=F^{R}_{\mu}(0)/\sqrt{F^{R}_{2N}(0)}.\end{split}$ (S.47)

The evaluation can be done in a similar way for MPM. The final expression is
the same as Eq. (3) in the main text.

In the case of finite circuit depth and in the presence of noise, we use Eq.
(S.34) to modify Eq. (S.47) as

$\psi^{L,R}_{\mu}\to\sqrt{|f(\omega,\Gamma)|}\,\psi^{L,R}_{\mu},$ (S.48)

where $\omega\sim 2^{-O(N)}$ is the frequency of the Majorana mode for a
system of finite size and $\Gamma$ is the effective decay rate. Since
$|f(\omega,\Gamma)|\leq 1$, the approximation given by Eq. (S.33) leads to a
uniform damping of the wavefunction as seen in Eq. (S.48). Normalizing the
wavefunction eliminates this effect. However, more general noise models would
result in more complex effects.

##  IV Section 4: Two-point correlation function

In this section we will assume a non-interacting scenario where $\lambda=0$.
For simplicity, we will focus on zero-frequency modes, although
$\pi$-frequency modes can be treated similarly. First, we will establish a
connection between the two-point function in Eq. (5) and an expectation of the
Fourier map given by

$T_{\mu\nu}=\langle\tilde{\psi}_{0}|\mathcal{F}_{0}(\gamma_{\mu}\gamma_{\nu})|\tilde{\psi}_{0}\rangle,$
(S.49)

where we have chosen the initial product state in the form of the product
state
$|\tilde{\psi}_{0}\rangle=|\psi_{a}\rangle|s_{2}\rangle\dots|s_{N-1}\rangle|\psi_{a}\rangle$,
where $|\psi_{a}\rangle=\cos a|0\rangle+i\sin a|1\rangle$, and $|s_{i}\rangle$
are arbitrary states in the $Z$ basis, $s_{i}\in\\{0,1\\}$. Next, we express
the Majorana operators using free fermion eigenmodes as

$\gamma_{\mu}=\sum_{k=1}^{2N}\psi^{*}_{k\mu}\Delta^{0}_{k},$ (S.50)

where the wavefunctions $\psi_{k\mu}$ are defined in Eq. (S.27). This
decomposition allows us to write

$\begin{split}\mathcal{F}_{0}&(\gamma_{\mu}\gamma_{\nu})=\\\
&=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}\sum_{k,k^{\prime}=1}^{2N}\psi^{*}_{k\mu}\psi^{*}_{k^{\prime}\nu}e^{-i(\omega^{0}_{k}+\omega^{0}_{k^{\prime}})n}\Delta^{0}_{k}\Delta^{0}_{k^{\prime}}\\\
&=\lim_{N\to\infty}(A_{\mu\nu}+B_{\mu\nu}),\end{split}$ (S.51)

where the variable $A_{\mu\nu}$ denotes the contribution from the set of zero-
frequency modes ($B_{0}$), while $B_{\mu\nu}$ represents the contribution from
pairs of modes with opposite frequencies from the complementary set
$\overline{B_{0}}$,

$\begin{split}&A_{\mu\nu}=\sum_{k,k^{\prime}\in
B_{0}}\psi_{k\mu}\psi_{k^{\prime}\nu}\Delta^{0}_{k}\Delta^{0}_{k^{\prime}},\\\
&B_{\mu\nu}=\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi^{*}_{\sigma(k)\nu}\Delta^{0}_{k}\Delta^{0}_{\sigma(k)}\\\
&\qquad\qquad=\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi_{k\nu}\Delta^{0}_{k}\Delta^{0{\dagger}}_{k}\end{split}$
(S.52)

The remaining terms vanish in the limit $D\to\infty$. Here $\sigma(k)$
represents the opposite frequency mode with respect to mode $k$. We also take
advantage of the fact that we can always choose the zero-frequency modes $k\in
B_{0}$ to be real-valued, $\psi_{k\mu}=\psi^{*}_{k\mu}$, while the remaining
modes satisfy the relations $\psi_{\sigma(k)\mu}=\psi^{*}_{k\mu}$,
$\omega_{\sigma(k)}^{0}=-\omega_{k}^{0}$, and
$\Delta^{0}_{\sigma(k)}=\Delta^{0{\dagger}}_{k}$.

Using these notations, the target two-point correlation function can be
expressed as

$T_{\mu\nu}=\lim_{N\to\infty}\Bigl{(}\langle\tilde{\psi}_{0}|A_{\mu\nu}|\tilde{\psi}_{0}\rangle+\langle\tilde{\psi}_{0}|B_{\mu\nu}|\tilde{\psi}_{0}\rangle\Bigl{)}.$
(S.53)

Next, we use Eq. (S.13) to express the expectation value

$\begin{split}&\langle\tilde{\psi}_{0}|\Delta^{0}_{k}\Delta^{0}_{k^{\prime}}|\tilde{\psi}_{0}\rangle=\sum_{\mu,\mu^{\prime}=1}^{2N}\psi_{k\mu}\psi_{k^{\prime}\mu^{\prime}}\langle\tilde{\psi}_{0}|\gamma_{\mu}\gamma_{\mu^{\prime}}|\tilde{\psi}_{0}\rangle,\\\
&\langle\tilde{\psi}_{0}|\Delta^{0}_{k}\Delta^{0{\dagger}}_{k^{\prime}}|\tilde{\psi}_{0}\rangle=\sum_{\mu,\mu^{\prime}=1}^{2N}\psi_{k\mu}\psi^{*}_{k^{\prime}\mu^{\prime}}\langle\tilde{\psi}_{0}|\gamma_{\mu}\gamma_{\mu^{\prime}}|\tilde{\psi}_{0}\rangle.\end{split}$
(S.54)

The expected values on the right-hand side of these equations can be expressed
in terms of the chosen product state $|\tilde{\psi}_{0}\rangle$ as

$\begin{split}\langle\tilde{\psi}_{0}|\gamma_{\mu}\gamma_{\mu^{\prime}}|\tilde{\psi}_{0}\rangle=\delta_{\mu\mu^{\prime}}&+iC_{1}(\delta_{\mu
1}\delta_{\mu^{\prime}2N}-\delta_{\mu 2N}\delta_{\mu^{\prime}1})\\\
&+iC_{2}(\delta_{\mu 1}\delta_{\mu^{\prime}2}-\delta_{\mu
2}\delta_{\mu^{\prime}1})\\\ &+iC_{2}(\delta_{\mu
2N-1}\delta_{\mu^{\prime}2N}-\delta_{\mu
2N}\delta_{\mu^{\prime}2N-1}),\end{split}$ (S.55)

where we used the notations

$\begin{split}&C_{1}=\langle\psi_{a}|Y|\psi_{a}\rangle^{2}\prod_{i=2}^{N-1}\langle
s_{i}|(-Z_{i})|s_{i}\rangle=\,(-1)^{N+S}\sin^{2}(2a),\\\
&C_{2}=\langle\psi_{a}|Z|\psi_{a}\rangle=\cos 2a.\end{split}$ (S.56)

where $S=\sum_{i=2}^{N-1}s_{i}$. We use these expressions to obtain

$\begin{split}\langle\tilde{\psi}_{0}|\Delta^{0}_{k}\Delta^{0}_{k^{\prime}}|\tilde{\psi}_{0}\rangle&=\sum_{\mu=1}^{2N}\psi_{k\mu}\psi_{k^{\prime}\mu}+iC_{1}(\psi_{k1}\psi_{k^{\prime}2N}-\psi_{k2N}\psi_{k^{\prime}1})\\\
&\qquad+iC_{2}(\psi_{k1}\psi_{k^{\prime}2}-\psi_{k1}\psi_{k^{\prime}2}\\\
&\qquad\qquad+\psi_{k2N-1}\psi_{k^{\prime}2N}-\psi_{k2N}\psi_{k^{\prime}2N-1}).\end{split}$
(S.57)

and

$\begin{split}\langle\tilde{\psi}_{0}|\Delta^{0}_{k}\Delta^{0{\dagger}}_{k^{\prime}}|\tilde{\psi}_{0}\rangle&=\sum_{\mu=1}^{2N}\psi_{k\mu}\psi^{*}_{k^{\prime}\mu}+iC_{1}(\psi_{k1}\psi^{*}_{k^{\prime}2N}-\psi_{k2N}\psi^{*}_{k^{\prime}1})\\\
&\qquad+iC_{2}(\psi_{k1}\psi^{*}_{k^{\prime}2}-\psi_{k1}\psi^{*}_{k^{\prime}2}\\\
&\qquad\qquad+\psi_{k2N-1}\psi^{*}_{k^{\prime}2N}-\psi_{k2N}\psi^{*}_{k^{\prime}2N-1}).\end{split}$
(S.58)

For real-valued zero frequency modes, we use the orthogonality condition

$\forall k\in
B_{0}:\qquad\sum_{\mu=1}^{2N}\psi_{k\mu}\psi_{k^{\prime}\mu}=\delta_{kk^{\prime}}.$
(S.59)

As a result, the expression for the contribution of the zero frequency modes
to the two-point correlation function can be written as

$\begin{split}\langle\tilde{\psi}_{0}|&A_{\mu\nu}|\tilde{\psi}_{0}\rangle=\sum_{k,k^{\prime}\in
B_{0}}\psi_{k\mu}\psi_{k^{\prime}\nu}\langle\tilde{\psi}_{0}|\Delta^{0}_{k}\Delta^{0}_{k^{\prime}}|\tilde{\psi}_{0}\rangle\\\
=&\sum_{k\in B_{0}}\psi_{k\mu}\psi_{k\nu}\\\ &+iC_{1}\Bigl{(}\sum_{k\in
B_{0}}\psi_{k\mu}\psi_{k1}\sum_{k^{\prime}\in
B_{0}}\psi_{k^{\prime}\nu}\psi_{k^{\prime}2N}-(\mu\leftrightarrow\nu)\Bigl{)}\\\
&+iC_{2}\Bigl{(}\sum_{k\in B_{0}}\psi_{k\mu}\psi_{k1}\sum_{k^{\prime}\in
B_{0}}\psi_{k^{\prime}\nu}\psi_{k^{\prime}2}-(\mu\leftrightarrow\nu)\Bigl{)}\\\
&+iC_{2}\Bigl{(}\sum_{k\in B_{0}}\psi_{k\mu}\psi_{k2N-1}\sum_{k^{\prime}\in
B_{0}}\psi_{k^{\prime}\nu}\psi_{k^{\prime}2N}-(\mu\leftrightarrow\nu)\Bigl{)},\end{split}$
(S.60)

where $(\mu\leftrightarrow\nu)$ is the same term as before with the indices
$\mu$ and $\nu$ swapped.

Similarly, the second part has the form

$\begin{split}\langle\tilde{\psi}_{0}|B_{\mu\nu}&|\tilde{\psi}_{0}\rangle=\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi_{k\nu}\\\
&+\Bigl{(}iC_{1}\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi_{k1}\psi_{k\nu}\psi^{*}_{k2N}\\\
&+iC_{2}\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi_{k1}\psi_{k\nu}\psi^{*}_{k2}\\\
&+iC_{2}\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi_{k2N-1}\psi_{k\nu}\psi^{*}_{k2N}-(\mu\leftrightarrow\nu)\Bigl{)}.\end{split}$
(S.61)

Since we assume that the bulk of the system is delocalized, the wavefunctions
of single-fermion modes with non-zero frequency must satisfy
$\psi_{k\mu}\propto O(N^{-1/2})$. This means that all terms in Eq. (S.61)
except the first have $O(N^{-1})$ scaling and therefore vanish in the limit
$N\to\infty$. Therefore, combining the contributions in Eqs. (S.60) and (S.61)
we get

$\begin{split}T_{\mu,\nu}=\delta_{\mu\nu}&+\lim_{N\to\infty}\Bigl{[}iC_{1}\sum_{k\in
B_{0}}\psi_{k\mu}\psi_{k1}\sum_{k^{\prime}\in
B_{0}}\psi_{k^{\prime}\nu}\psi_{k^{\prime}2N}\\\ &+iC_{2}\sum_{k\in
B_{0}}\psi_{k\mu}\psi_{k1}\sum_{k^{\prime}\in
B_{0}}\psi_{k^{\prime}\nu}\psi_{k^{\prime}2}\\\ &+iC_{2}\sum_{k\in
B_{0}}\psi_{k\mu}\psi_{k2N-1}\sum_{k^{\prime}\in
B_{0}}\psi_{k^{\prime}\nu}\psi_{k^{\prime}2N}\\\
&-(\mu\leftrightarrow\nu)\Bigl{]}.\end{split}$ (S.62)

Next, we analyze the behavior of this function for a few characteristic values
of $\mu$ and $\nu$. First we consider the combination $\mu=1$ and $\nu=2$ and
assume that the wavefunctions $\psi_{k\mu}$ for zero modes $k\in B_{0}$ are
strongly localized, i.e. $\psi_{k\mu}\psi_{k\nu}\leq 2^{-c|\mu-\nu|}$, $c>0$.
Then the correlation function takes the form

$T_{1,2}=iC_{2}\sum_{k\neq k^{\prime}\in
B_{0}}(\psi^{2}_{k1}\psi^{2}_{k^{\prime}2}-\psi^{2}_{k2}\psi^{2}_{k^{\prime}1}).$
(S.63)

In particular, if the only modes are Majorana modes, then
$\psi^{L}_{1}\psi^{R}_{2}\sim\psi^{L}_{1}\psi^{R}_{2}\sim O(2^{-\Theta(N)})$,
and this quantity vanishes in the thermodynamic limit. In contrast, if there
is a localized state at the boundary of the system, there exists a pair $b\neq
b^{\prime}$ for which $T_{1,2}=O(1)$.

At the same time, the correlation function for $\mu=1$ and $\nu=2N$ has the
form

$T_{1,2N}=iC_{1}\sum_{b,b^{\prime}\in
B_{0}}(\psi^{2}_{k1}\psi^{2}_{k^{\prime}2N}-\psi^{2}_{k2N}\psi^{2}_{k^{\prime}1}).$
(S.64)

This expression does not vanish for both Majorana and trivial modes. For other
points $\mu$ and $\nu$ far from the boundaries, the two-point correlation
function vanishes.

To model the trivial system with approximate zero-energy localized boundary
eigenmodes, we consider the Hamiltonian

$H(t)=\sum_{j=1}^{N-1}J_{j}(t)\,X_{j}X_{j+1}+\sum_{j=1}^{N}h_{j}(t)Z_{j},$
(S.65)

with the same protocol. The couplings and fields for the bulk qubits are the
same, $J_{j}T=\pi/16$ for $j\neq 1,N-1$ and $h_{j}T=\pi/4$ for $j\neq 1,N$. At
the same time, we set $J_{1}=J_{N-1}=h_{1}=h_{N}=0$. Due to decoupling of the
boundary qubits, Majorana operators corresponding to these qubits (i.e.
$\gamma_{1}$, $\gamma_{2}$, $\gamma_{2N-1}$, and $\gamma_{2N}$) are integrals
of motion.

##  V Section 5: Braiding map

In this section, we provide a rigorous proof of the properties of the braiding
map in Eq. (7) in the main text. First, we formulate Lemma 1, which
establishes the action of the map on MZM operators.

######  Lemma 1.

Suppose Eq. (S.43) holds and $\psi^{L}_{1}=\psi^{R}_{2N}=\xi$, $\xi^{2}\geq
1/2$. Then for the angle $\alpha_{0}=\arcsin(1/\sqrt{2}\xi)$ the action of the
map in Eq. (7) on MZM operators is

$\mathcal{E}_{\alpha_{0}}(\Gamma_{R})=p\Gamma_{L},\quad\mathcal{E}_{\alpha_{0}}(\Gamma_{L})=-p\Gamma_{R},$
(S.66)

where $p=\sqrt{2\xi^{2}-1}$.

Lemma 1 leads to Eq. (8) in the main text. This lemma applies only to systems
where MZMs are the only zero-frequency modes overlapping with single-fermion
operators, as manifested by Eq. (S.43). At the same time, it applies to a
generic setting including the interacting case $|\lambda|>0$.

Next, we formulate Lemma 2, which gives the expression for the map action on
Majorana operators.

######  Lemma 2.

Under conditions of Lemma 1, the map in Eq. (7) satisfies

$\mathcal{E}_{\alpha_{0}}(\gamma_{\mu})=p(\psi^{R}_{\mu}\Gamma_{L}-\psi^{L}_{\mu}\Gamma_{R})+\delta
C,$ (S.67)

where the norm of the correction operator is bounded as

$\begin{split}&\|\delta
C\|\leq\frac{1}{\xi^{2}}\max_{\mu}\sqrt{\sum_{\nu}\kappa_{\mu\nu}^{2}},\\\
&\kappa_{\mu\nu}=\lim_{N\to\infty}\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\psi^{*}_{k\nu}(\psi^{2}_{k1}+\psi^{2}_{k2N}),\end{split}$
(S.68)

where $\|\cdot\|$ is the spectral norm.

This conclusion leads us to Eq. (9) when the bulk modes are delocalized.
Indeed, in this case $\kappa_{\mu\nu}=O(N^{-1})$, therefore the absolute value
in Eq. (S.68) scales as $O(N^{-1/2})$ and the correction vanishes in the
thermodynamic limit $N\to\infty$.

As an alternative to the theoretical prediction in Lemma 1, we explore the
possibility of finding the correct angle $\alpha_{0}$ by optimization. To
illustrate this method, we run the circuit for multiple angles and find the
optimal value of $\alpha_{0}$. To estimate the braiding quality, we propose a
cost function that favors a braided wavefunction if it is located at the
opposite boundary of the chain. In particular, for braiding the left
eigenmode, our cost function is given by

$\begin{split}\mathcal{L}(\alpha_{0})=\sum_{x=1}^{N}\Bigl{(}|\tilde{\psi}^{L}_{2x-1}|^{2}&+|\tilde{\psi}^{L}_{2x}|^{2}\Bigl{)}(N-x)^{2},\end{split}$
(S.69)

where $\tilde{\psi}^{L,R}_{\mu}$ are the braided wavefunctions in Eq. (10) for
the angle $\alpha_{0}$. The cost functions for different runs on the same
device are shown in Fig. S3 with the same parameters as in Fig. (3). We use a
simple polynomial approximation to find the optimal angle $\alpha_{0}$
corresponding to the minimum of the approximation of the curve. In our
experiment, the optimal value ($0.256934\pi$) differs slightly from the
theoretical value ($0.263127\pi$). This difference is due to the presence of
disorder and noise, which modify the original Hamiltonian dynamics.

Figure S3: Finding the angle by optimization. Cost function in Eq. (S.69)
derived from experimental data (dots) approximated by a polynomial
$I(\alpha_{0})=a_{1}\alpha_{0}^{2}+a_{2}\alpha_{0}+a_{3}$ (solid curve), where
$a_{i}$ are coefficients. The dotted vertical line shows the theoretically
predicted value of the angle $\alpha_{0}$. We use a 5-qubit system on the
$ibm\\_hanoi$ device with parameters $\phi=\pi/16$, $\theta=\pi/4$, and
$\varphi=0$ and maximum number of cycles $D=11$.

Below we provide the proofs for both Lemma 1 and Lemma 2.

Proof of Lemma 1. As a first step, we use the commutativity
$[U_{F},\Gamma_{s}]=0$ for MZM operators and the decomposition in Eq. (S.37)
to rewrite

$\begin{split}\mathcal{E}_{\alpha}(\Gamma_{s})&=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}U^{{\dagger}n}_{F}V^{\dagger}(\alpha)\Gamma_{s}V(\alpha)U^{n}_{F}\\\
&=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}\sum_{\mu=1}^{2N}\psi^{s}_{\mu}U^{{\dagger}n}_{F}V^{\dagger}(\alpha)\gamma_{\mu}V(\alpha)U^{n}_{F}.\end{split}$
(S.70)

where we define $V(\alpha):=\exp(-\alpha\gamma_{1}\gamma_{2N})$. Now we can
use the decomposition
$V(\alpha)=\cos\alpha\,I-\sin\alpha\,\gamma_{1}\gamma_{2N}$ and commutation
relation between Majorana operators to express

$\begin{split}V^{\dagger}(\alpha)\gamma_{\mu}V(\alpha)=\gamma_{\mu}&-\delta_{\mu
1}(2\sin^{2}\alpha\gamma_{1}+\sin 2\alpha\gamma_{2N})\\\ &-\delta_{\mu
2N}(2\sin^{2}\alpha\gamma_{2N}-\sin 2\alpha\gamma_{1}).\end{split}$ (S.71)

Combining this result with Eq. (S.70), we get

$\begin{split}\mathcal{E}_{\alpha}(\Gamma_{s})=\Gamma_{s}&-(2\psi^{s}_{1}\sin^{2}\alpha-\psi^{s}_{2N}\sin
2\alpha)\mathcal{F}_{0}(\gamma_{1})\\\
&-(2\psi^{s}_{2N}\sin^{2}\alpha+\psi^{s}_{1}\sin
2\alpha)\mathcal{F}_{0}(\gamma_{2N}),\end{split}$ (S.72)

where we use the notation for the Fourier map from Eq. (S.18). Then, using the
property in Eq. (S.43), we get the expression

$\begin{split}\mathcal{E}_{\alpha}(\Gamma_{s})=\Gamma_{s}&-\lim_{N\to\infty}\psi^{L}_{1}(2\psi^{s}_{1}\sin^{2}\alpha-\psi^{s}_{2N}\sin
2\alpha)\Gamma_{L}\\\
&-\lim_{N\to\infty}\psi^{R}_{2N}(2\psi^{s}_{2N}\sin^{2}\alpha+\psi^{s}_{1}\sin
2\alpha)\Gamma_{R}.\end{split}$ (S.73)

Inserting the particular values $s\in\\{L,R\\}$ and taking into account that
$\psi^{L}_{2N}\sim\psi^{R}_{1}\sim O(2^{-\Theta(N)})$ vanish as $N\to 0$, we
get

$\begin{split}&\mathcal{E}_{\alpha}(\Gamma_{R})=(1-2(\psi^{R}_{2N})^{2}\sin^{2}\alpha)\Gamma_{R}+\sin
2\alpha\,\psi^{L}_{1}\psi^{R}_{2N}\Gamma_{L},\\\
&\mathcal{E}_{\alpha}(\Gamma_{L})=(1-2(\psi^{L}_{1})^{2}\sin^{2}\alpha)\Gamma_{L}-\sin
2\alpha\,\psi^{L}_{1}\psi^{R}_{2N}\Gamma_{R}.\end{split}$ (S.74)

By putting $\psi^{L}_{1}=\psi^{R}_{2N}=\xi$, according to the Lemma’s
conditions, and choosing the angle

$\alpha\to\alpha_{0}=\arcsin\frac{1}{\sqrt{2}\xi},$ (S.75)

expression in Eq. (S.74) converts into the statement of the Lemma. This step
concludes our proof.

Proof of Lemma 2. Assuming that Majorana modes are the only single-fermion
modes with zero frequency, we can write

$\gamma_{\mu}=\psi^{L}_{\mu}\Gamma_{L}+\psi^{R}_{\mu}\Gamma_{R}+\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\Delta^{0}_{k},$
(S.76)

where $\overline{B_{0}}$ is the set of single-fermion modes in Eq. (S.28)
whose frequencies are distinct from zero. Then, the action of the target map
on Majorana operator is

$\begin{split}\mathcal{E}_{\alpha_{0}}(\gamma_{\mu})&=\psi^{L}_{\mu}\mathcal{E}_{\alpha_{0}}(\Gamma_{L})+\psi^{R}_{\mu}\mathcal{E}_{\alpha_{0}}(\Gamma_{R})+\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\mathcal{E}_{\alpha_{0}}(\Delta^{0}_{k})\\\
&=p(\psi^{R}_{\mu}\Gamma_{L}-\psi^{L}_{\mu}\Gamma_{R})+\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}\mathcal{E}_{\alpha_{0}}(\Delta^{0}_{k}),\end{split}$
(S.77)

where we used Lemma 1 to express the action of the target map on topological
Majorana operators $\Gamma_{L,R}$. The last term, in turn, can be evaluated
using the explicit form of the map and Eq. (S.28),

$\begin{split}\mathcal{E}_{\alpha}&(\Delta_{k}^{0})=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{-i\omega^{0}_{k}n}U_{F}^{{\dagger}n}V^{\dagger}(\alpha)\Delta^{0}_{k}V(\alpha)U_{F}^{n}\\\
&=\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}\sum_{\nu=1}^{2N}\psi_{k\nu}e^{-i\omega^{0}_{k}n}U^{{\dagger}n}_{F}V^{\dagger}(\alpha)\gamma_{\nu}V(\alpha)U^{n}_{F}.\end{split}$
(S.78)

Using the result in Eq. (S.71), we get

$\begin{split}\mathcal{E}_{\alpha}(\Delta_{k}^{0})=&\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{-2i\omega^{0}_{k}n}\Delta^{0}_{k}\\\
&-(2\psi_{k1}\sin^{2}\alpha-\psi_{k2N}\sin
2\alpha)\mathcal{F}_{-\omega^{0}_{k}}(\gamma_{1})\\\
&-(2\psi_{k2N}\sin^{2}\alpha+\psi_{k1}\sin
2\alpha)\mathcal{F}_{-\omega^{0}_{k}}(\gamma_{2N}),\end{split}$ (S.79)

where we use the notation for the Fourier map from Eq. (4). The first term of
this expression vanishes for the set $k\in\overline{B_{0}}$. In turn, the
action of the Fourier map is

$\begin{split}\mathcal{F}_{-\omega_{k}}(\gamma_{\mu})&=\sum_{k^{\prime}\in\overline{B_{0}}}\psi^{*}_{k^{\prime}\mu}\Delta^{0}_{k^{\prime}}\lim_{N,D\to\infty}\frac{1}{D}\sum_{n=0}^{D-1}e^{-i(\omega^{0}_{k}+\omega^{0}_{k^{\prime}})n}\\\
&=\psi_{k\mu}\Delta_{k}^{0{\dagger}},\end{split}$ (S.80)

where we used the fact that the wave-functions corresponding to opposite
frequencies, i.e. $\omega^{0}_{k^{\prime}}=-\omega^{0}_{k}$, satisfy
$\psi_{k^{\prime}\mu}=\psi_{k\mu}^{*}$ and
$\Delta^{0}_{k^{\prime}}=\Delta_{k}^{0{\dagger}}$. Then, Eq. (S.79) takes the
form

$\begin{split}\mathcal{E}_{\alpha_{0}}(\Delta^{0}_{k})&=-2\sin^{2}\alpha_{0}(\psi_{k1}^{2}+\psi_{k2N}^{2})\Delta_{k}^{0{\dagger}}\\\
&=-\frac{1}{\xi^{2}}(\psi_{k1}^{2}+\psi_{k2N}^{2})\Delta_{k}^{0{\dagger}}.\end{split}$
(S.81)

Inserting this expression into Eq. (S.77), we finally get

$\begin{split}\mathcal{E}_{\alpha_{0}}(\gamma_{\mu})&=p(\psi^{R}_{\mu}\Gamma_{L}-\psi^{L}_{\mu}\Gamma_{R})\\\
&\quad-\frac{1}{\xi^{2}}\sum_{k\in\overline{B_{0}}}\psi^{*}_{k\mu}(\psi_{k1}^{2}+\psi_{k2N}^{2})\Delta_{k}^{0{\dagger}}\\\
&=p(\psi^{R}_{\mu}\Gamma_{L}-\psi^{L}_{\mu}\Gamma_{R})\\\
&\quad-\frac{1}{\xi^{2}}\sum_{k\in\overline{B_{0}}}\sum_{\nu=1}^{2N}\psi^{*}_{k\mu}\psi^{*}_{k\nu}(\psi_{k1}^{2}+\psi_{k2N}^{2})\gamma_{\nu}\\\
&=-\frac{1}{\xi^{2}}\sum_{\nu}\kappa_{\mu\nu}\gamma_{\nu},\end{split}$ (S.82)

where $\kappa_{\mu\nu}$ are real numbers due to the fact that
$\psi_{k\mu}=\psi^{*}_{k^{\prime}\mu}$ for symmetric pairs
$\omega_{k}^{0}=-\omega_{k^{\prime}}^{0}$ in $\overline{B_{0}}$. The maximum
singular value of the operator in Eq. (S.82) is $\Lambda^{\rm
max}_{\mu}=\xi^{-2}\sqrt{\sum_{\nu}\kappa_{\mu\nu}^{2}}$. This result leads us
to the statement of the Lemma.