# Scalable quantum control and non-abelian anyon creation in the Kitaev
honeycomb model

Omar Raii  Blackett Laboratory, Imperial College London, SW7 2AZ, United
Kingdom  Center for Engineered Quantum Systems, Macquarie University, NSW
2109, Australia  Florian Mintert  Blackett Laboratory, Imperial College
London, SW7 2AZ, United Kingdom  Helmholtz-Zentrum Dresden-Rossendorf,
Bautzner Landstraße 400, 01328 Dresden, Germany  Daniel Burgarth  Center for
Engineered Quantum Systems, Macquarie University, NSW 2109, Australia

(2022)

###### Abstract

The Kitaev honeycomb model is a system allowing for experimentally realisable
quantum computation with topological protection of quantum information.
Practical implementation of quantum information processing typically relies on
adiabatic, i.e. slow dynamics. Here we show that the restriction to adiabatic
dynamics can be overcome with optimal control theory, enabled by an extension
of the fermionization of the Kitaev honeycomb model to the time-dependent
case. Moreover we present a quantum control method that is applicable to large
lattice models due to sub-exponential scaling.

##  I Introduction

Non-abelian anyons are the foundational theoretical tools for topological
quantum computation. These generalizations of bosons and fermions, which can
only exist in two-dimensional systems Kitaev (2006); Wilczek (1982), allow for
a topological form of gate implementation due to their non-trivial braiding
statistics Sarma  _et al._ (2006); Pachos (2012); Kitaev (2003). Logical
qubits are encoded non-locally and anyons may be braided around one another
and subsequently fused to carry out computations and measurements respectively
as part of the overall implementation of a quantum algorithm Pachos (2012);
Brennen and Pachos (2007).

The Kitaev honeycomb model is a notable example of a relatively simply defined
system with non-trivial topological order Pachos (2012); Kitaev (2003). The
simplicity of the model’s definition has led to multiple proposals in recent
years for experimental realization Knolle and Moessner (2019); Takagi  _et
al._ (2019). In this system, anyons manifest as vortices introduced into the
model which may be fused to create fermionic excitations corresponding to
anyonic fusion rules known as Ising anyons Lahtinen (2011); Génetay Johansen
and Simula (2021). Anyonic braiding itself and indeed any form of particles
being interchanged is generally assumed to be an adiabatic process Cheng  _et
al._ (2011) so that unwanted excitations may be safely suppressed.

In practice, however, the restriction to adiabatic dynamics is typically
conflicting with the requirement to realize all operations on a time-scale
that is short compared to the system’s coherence time Lai and Zhang (2020).
Quantum control has proven successful in speeding up adiabatic evolution in a
wide range of scenarios He  _et al._ (2016); Zhou  _et al._ (2017); Xu  _et
al._ (2019); Wang  _et al._ (2012); Simsek and Mintert (2021) suggesting its
suitability for anyon creation in topological systems. Common optimal control
techniques, however, are limited in their applicability to the Kitaev
honeycomb model. Due to the exponential scaling of composite quantum systems,
numerical simulations of the time-dependent Kitaev honeycomb model are only
possible for very small system sizes. fermionization of the Kitaev honeycomb
model allows for improved scaling and solving larger systems but thus far this
has been restricted to systems with time-independent Hamiltonians. The scope
of this article is to use quantum control in applying fermionization within
the context of a time-dependent version of the model. This will demonstrate
that optimal control does indeed provide access to faster-than-adiabatic anyon
creation.

Sec. II of this paper provides a brief overview of anyons in the Kitaev
honeycomb model to setup the operators and terminology required for the
control problem. Sec. III describes quantum control and pulse optimization for
anyon creation and sets up the key result on fermionization, which is proven
in the appendix. Sec. IV presents explicit numerical results of the optimal
control problem defined. An overall summary of results and conclusions are
presented in the final section.

##  II Anyons in the Kitaev honeycomb model

Numerous methods have been demonstrated which allow for solving the Kitaev
honeycomb model including that of Feng et al. Feng  _et al._ (2007) using the
Jordan-Wigner transformation as well as that of Kells et al. Kells  _et al._
(2009); Brennan2018TheGgeq2 which employs a more specific Jordan-Wigner type
mapping that describes the system in terms of hard-core bosons which are then
fermionized. Such methods have many advantages and have applications to other
lattice structures Kells  _et al._ (2010). The method most suitable for our
work in extending to a time-dependent model is Kitaev’s original Majorana
fermionization procedure and in this section we review this method Kitaev
(2006); Knolle (2016). This is followed by a demonstration of how vortex
creation is implemented within the model Pachos (2012); Lahtinen and Pachos
(2017) and how this corresponds to creation of non-abelian anyons.

###  II.1 Diagonalising the honeycomb

Although diagonalization of the Kitaev honeycomb model is not required for
solving the time-dependent control problem we define later, we still outline
its strategy, as the operators and terminology introduced will also play a
role in the time-dependent version. The model takes its name from its
hexagonal lattice geometry consisting of spin-1/2 particles located at the
vertices of hexagonal plaquettes, as shown in Fig. 1. It is defined by the
Hamiltonian

$\displaystyle H=-\sum_{\\{j,k\\}\in
N_{2}}J_{s}\sigma^{s}_{j}\sigma^{s}_{k}-K\sum_{\\{j,k,l\\}\in
N_{3}}\sigma^{x}_{j}\sigma^{y}_{k}\sigma^{z}_{l},$ (1)

where $N_{2}$ correspond to honeycomb edges and $N_{3}$ are certain triplets
described further below. There are three types $s=x,y,z$ of two-body nearest
neighbor Pauli interactions determined by the position of the edge in the
lattice, highlighted in three colors in Fig. 1. The three-body terms act
within each hexagonal plaquette in the following way: three adjacent spins
contribute to a three-body interaction term with the middle spin interacting
through the Pauli operator corresponding to the link pointing outwards from
the plaquette, while each of the two remaining spins interact through the
Pauli operator corresponding to their link to the middle spin. For example in
the plaquette highlighted in Fig. 1, one of the three-body interaction terms
would be $\sigma^{y}_{1}\sigma^{x}_{2}\sigma^{z}_{3}$, with $5$ similar terms
following clockwise along the hexagonal plaquette. While the two-body part of
the Hamiltonian allows for the model to be solved by a process of Majorana
fermionization, the three-body part preserves the solvability of the model
while also breaking time-reversal symmetry and it consequently gives the
system non-trivial topological order Kitaev (2006); Pachos (2012).

Figure 1: A honeycomb lattice showing the three kinds of interactions between
neighbouring spins which are on the vertices of hexagonal plaquettes. Red,
blue and green links indicate $\sigma^{x}\sigma^{x}$, $\sigma^{y}\sigma^{y}$
and $\sigma^{z}\sigma^{z}$ interactions respectively. The hexagonal plaquette
operators correspond to many conserved quantities in the system. A generic
plaquette whose spin sites have been numbered 1 to 6 is highlighted in grey.

For every hexagonal plaquette, a corresponding plaquette operator $W_{p}$ may
be defined which acts on every spin with the Pauli operator of the outward
pointing interaction, so that for example on the numbered gray plaquette in
Fig. 1 we have
$W_{p}=\sigma^{z}_{1}\sigma^{x}_{2}\sigma^{y}_{3}\sigma^{z}_{4}\sigma^{x}_{5}\sigma^{y}_{6}$.
Each plaquette operator squares to the identity so that its eigenvalues adopt
the values $\pm 1$ only. Since the plaquette operators all commute with the
Hamiltonian and with one another, the system Hilbert space is naturally
partitioned into simultaneous eigenspaces of all plaquette operators.

Negative plaquette eigenstates are known as vortices and by a well known
theorem Lieb (1994) it is known that the ground state eigenspace is in the no-
vortex sector Pachos (2012). Different vortex sectors relate to the presence
of anyons localized at the respective vortex plaquettes.

While solving the system remains an intractable exponential problem even after
restriction to one vortex sector, a mapping of the problem into a Majorana
fermionic Hamiltonian provides a pathway towards diagonalization. For a more
detailed breakdown of diagonalization see Kitaev (2006); Pachos (2012). The
mapped Majorana fermionic Hamiltonian is defined by replacing spin qubit sites
$j$ with two fermionic sites and their corresponding creation operators
$a^{\dagger}_{j,1}$ and $a^{\dagger}_{j,2}$. For each site $j$, the real and
imaginary parts of the two fermionic modes constitute a total of four Majorana
modes per site. The Majorana creation/annihilation operators are defined

$\displaystyle b^{x}_{j}$
$\displaystyle=a_{j,1}+a_{j,1}^{\dagger},~{}~{}~{}b^{y}_{j}=\frac{1}{i}(a_{j,1}-a_{j,1}^{\dagger}),$
(2) $\displaystyle b^{z}_{j}$
$\displaystyle=a_{j,2}+a_{j,2}^{\dagger},~{}~{}~{}c_{j}=\frac{1}{i}(a_{j,2}-a_{j,2}^{\dagger}).$
(3)

Since with this mapping the Hilbert space is enlarged, a projection is
required to obtain vectors that correspond to states in the original Hilbert
space of the model. This requires the stabilizer projector Kitaev (2006);
Pachos (2012)

$\displaystyle P_{D}=\prod_{j}^{N}\left(\frac{1+D_{j}}{2}\right)$ (4)

with $D_{j}=b^{x}_{j}b^{y}_{j}b^{z}_{j}c_{j}$.

The newly mapped Hamiltonian

$\displaystyle H_{f}=\frac{iJ}{2}\sum_{\\{j,k\\}\in
N_{2}}\hat{u}_{jk}c_{j}c_{k}+\frac{iK}{2}\sum_{\\{j,l\\},\\{k,l\\}\in
N_{2}}\hat{u}_{jl}\hat{u}_{kl}c_{j}c_{k},$ (5)

is defined in terms of Majorana operators where we may also define link
operators $\hat{u}_{jk}=ib^{s}_{j}b^{s}_{k}$ in a system of $L$ total links
and $N$ spin qubits.

The eigenvalues $\pm 1$ of link operators allow for further partitioning of
each vortex sector into link sectors. To this end, we can define a
corresponding link sector projector

$\displaystyle P_{u}=\prod_{\\{j,k\\}\in
N_{2}}^{L}\frac{1+u_{jk}\hat{u}_{jk}}{2},$ (6)

which amounts to a tuple of chosen eigenvalues $u_{jk}\in\\{\pm 1\\}$.

Picking a certain link sector corresponds to fixing a gauge for a specific
vortex sector and leads to a quadratic fermionic Hamiltonian $P_{u}H_{f}P_{u}$
that is easily diagonalized Blaizot and Ripka (1986). The trivial gauge would
consist of setting all link eigenvalues to $u_{jk}=1$, alongside the
constraints imposed by their anti-symmetry $\hat{u}_{jk}=-\hat{u}_{kj}$. This
amounts to defining an orientation which for concrete purposes we define as
follows: a positive orientation on an $x$-link is directed from the bottom-
left qubit to the top-right one ($j$ to $k)$, for a $y$-link it is directed
from the bottom-right to the top-left, and for a $z$-link it is directed from
top to bottom.

The Hamiltonian becomes diagonal in a certain quasiparticle basis
$H=\sum_{\omega>0}\omega_{j}b_{j}^{\dagger}b_{j}-E_{g}$ and the ground state
is the quasiparticle vacuum state. As previously stated, for physical states
of the original Hamiltonian, the states will need to be projected; thus
eigenstates of the spin Hamiltonian $\ket{\Psi}$ are related to eigenstates of
the quadratic fermionic Hamiltonian $\ket{\psi}$ by

$\displaystyle\ket{\Psi}=P_{D}P_{u}\ket{\psi}.$ (7)

The energy gap in the system corresponding to the vortex energy $\delta$ is
simply the difference in ground state energy $E_{g}$ between differing vortex
sectors. In general this is the energy gap between pairs of vortices but in
exceptional examples, such as with the single plaquette lattice, this gap
corresponds to a single vortex.

Additionally the fermionic excitation gap is equivalent to the smallest
fermionic energy $\omega_{j}$ within a chosen vortex sector.

##  III Anyon creation as an optimal control problem

Let us first describe the usual adiabatic approach. Vortex creation is
implemented in terms of the time-dependent Hamiltonian

$\displaystyle H(t)$ $\displaystyle=H+\frac{t}{T}H_{\text{control}}^{j,k},$
(8)

comprised of the original honeycomb Hamiltonian as a drift combined with a
control Hamiltonian

$\displaystyle H_{\text{control}}^{j,k}$
$\displaystyle=2J_{jk}\sigma^{s}_{j}\sigma^{s}_{k}+2K\sum_{j,k\in\\{a,b,c\\}}\sigma^{x}_{a}\sigma^{y}_{b}\sigma^{z}_{c}.$
(9)

where $\sigma^{x}_{a}\sigma^{y}_{b}\sigma^{z}_{c}$ are those three-body terms
such that two of $a,b$ and $c$ are $j$ and $k$. There are four such three-body
terms for a chosen $j,k$ pair due to the definition of the three-body
interactions described in the previous section. This amounts to gradually
reversing the sign of a specific $s$-link, as well as the sign of the nearby
three-way interactions, using a linear time-dependence, with steepness and
therefore adiabaticity determined by the duration $T$ of the adiabatic
protocol. As we will see in the numerical examples, $T$ needs to be very large
to obtain a good fidelity.

We now want to set up anyon creation as an optimal control problem in the hope
that we can obtain similar fidelities in much shorter times compared to the
adiabatic evolution. To this end, generalize the time-dependence of Eq. (8) as

$\displaystyle H(t)$ $\displaystyle=H+f(t)H_{\text{control}}^{j,k},$ (10)

where $f(t)$ is the _ramp function_ defined such that $f(0)=0$ at the initial
point in time $t=0$ and such that $f(T)=1$ at the final point in time $t=T$.
Typically, $f(t)$ is assumed to be piecewise smooth or piecewise constant.

A typical figure of merit to be maximized is the state fidelity $\mathcal{F}$
defined in terms of an initial state $\ket{\Phi_{0}}$, the propagator
$U[f(t),T]$ induced by the time-dependent Hamiltonian $H(t)$ and a target
state that is meant to be created. In the present case, the initial state
$\ket{\Phi_{0}}$ would usually be the ground state of the Kitaev honeycomb
model whereas the target state $\ket{\Phi_{\text{target}}}$ is a state with an
additional vortex-pair created.

While optimizing for such a state fidelity is a generally successful approach,
it has two flaws when it comes to the Kitaev honeycomb model. Since
topological stability only arises for large lattice size, any practical
application of the model requires a vast number of qubits. Evaluating the
time-evolution operator therefore requires numerics in exponentially large
spaces. Secondly, even if the ground state is solvable analytically in the
free fermion picture, we would have to translate it back into the spin picture
to evaluate $\mathcal{F}$, which is again exponentially hard. In the
following, we will resolve both problems to obtain a scalable optimization
method.

###  III.1 Time-dependent fermionic picture

Let us first describe the time-independent case. In the fermionic picture, the
quadratic Majorana Hamiltonian can be written in the most general form with a
matrix $J_{jk}$ that incorporates all interaction factors $J$ as
$H=\frac{i}{2}\sum_{jk}J_{jk}c_{j}c_{k}$. When written in terms of full
fermionic creation and annihilation operators this is:

$\displaystyle H=\frac{1}{2}\alpha^{\dagger}M\alpha$ (11)

where $\alpha=(a_{1},\ldots,a_{N},a_{1}^{\dagger},\ldots,a_{N}^{\dagger})^{T}$
and the Hermitian matrix $M$

$\displaystyle M=\begin{pmatrix}\mu&\nu\\\ -\nu^{*}&-\mu^{*}\end{pmatrix}$
(12)

may be defined in terms of a Hermitian matrix $\mu$ and an antisymmetric
matrix $\nu$.

A canonical transformation $T$ can then be found so that
$TMT^{-1}=\text{diag}\\{\omega,-\omega\\}$ where $\omega$ is a diagonal
$2N$-by-$2N$ matrix. This allows for the Hamiltonian to be diagonalised in
terms of quasiparticle modes Blaizot and Ripka (1986).

We now consider how, in the fermionic picture, we may calculate the fidelity
between a state evolved from an initialized state by a _time-dependent
Hamiltonian_ towards a target state. To this end, we write a ground state of
$H$ as $\ket{\Phi_{0}}\equiv A\ket{\text{vac}}$, with the vacuum state
$\ket{\text{vac}}$ satisfying the relation $a_{j}\ket{\text{vac}}=0\;\forall
j$. The operator $A$ is some appropriately chosen function of creation and
annihilation operators. In the Appendix, we show that the state fidelity in
the fermionic picture is given by

$\displaystyle\mathcal{F}(t)$
$\displaystyle=|\bra{\Phi_{\text{target}}}kP_{D}P_{u}~{}V[f(t),t]A\ket{\text{vac}}|^{2}$
(13)
$\displaystyle=|\bra{\Phi_{\text{target}}}kP_{D}P_{u}A(t)V[f(t),t]\ket{\text{vac}}|^{2}.$
(14)

Here the projector $P_{u}$ is given by Eq. (6), $V[f(t),t]$ is the evolution
operator corresponding to the quadratic Hamiltonian $P_{u}H(t)P_{u}$ and
$P_{D}$ is given by Eq. (4), while $k$ is a real number which depends on the
specific lattice (see Appendix for specific examples). $A(t)$ is a Heisenberg
picture operator $A(t)\equiv V[f(t),t]AV[f(t),t]^{\dagger}$.

In analogy to Eqs (11) and (12) it is useful to decompose $P_{u}H(t)P_{u}$ as

$\displaystyle P_{u}H(t)P_{u}=P_{u}\frac{1}{2}\alpha^{\dagger}M(t)\alpha$ (15)

Since $A$ depends on annihilation and creation operators, we may write it as
$A(\alpha)$. It can then be shown Blaizot and Ripka (1986) that

$\displaystyle V[f(t),t]A(\alpha)V^{\dagger}(t)=A(W[f(t),t]\alpha)$ (16)

where $W(t)$ is the $2N$-by-$2N$ the time-ordered product solving the
differential equation

$\dot{W}[f(t),t]=iM[f(t)]W[f(t),t].$ (17)

This solves the problem of an exponentially sized evolution operator, as
$W(t)$ scales linearly in the system size. We will refer to calculations in
this picture as the Heisenberg picture, since it is directly obtained from the
Heisenberg equations of motions of $\alpha$. However, the problem of
expressing the target and initial state in the spin picture remains. This will
be tackled next.

###  III.2 Heisenberg fidelity as optimization target

In the previous paragraph we showed that the evolution is fully determined by
the Heisenberg picture of the quadratic Hamiltonian. If we knew a good target
evolution, rather than target state, we could therefore free ourselves from
the state picture and obtain all quantities directly in the Heisenberg
picture. The key idea here is to get back to the adiabatic evolution to find
such good evolution. We phrase such evolution directly in the Heisenberg
picture. To this end, let $W_{\text{ad}}$ be the solution of Eq. (17) for the
adiabatic ramp Eq. (8) with some suitably large duration $T_{\text{ad}}$. This
can be computed efficiently without having to refer to states. We define a
corresponding Heisenberg fidelity

$\displaystyle\mathcal{F}_{H}=\frac{1}{2N}\left|\Tr(W_{\text{ad}}^{\dagger}W[f(t),T])\right|$
(18)

This quantity obtains its maximum of $1$ if and only if the evolutions match
up to a phase and it can be used for efficient numerical optimization. In
Appendix B, we show that

$\displaystyle
1-\mathcal{F}_{H}\geq\frac{1}{32N^{3}}\left(1-\sqrt{\mathcal{F}}\right).$ (19)

This shows that $\mathcal{F}_{H}$ is a good surrogate for $\mathcal{F}$ and
may be optimized instead. In order to appropriately compare very high
fidelities and plot them logarithmically we give our results in terms of
infidelity $\mathcal{I}_{H}=1-\mathcal{F}_{H}$.

##  IV Numerical results

Here, we use the QuTiP Johansson  _et al._ (2013) implementation of the
gradient-ascent pulse engineering (GRAPE) algorithm Khaneja  _et al._ (2005);
Rowland and Jones (2012) using the limited memory version of the
Broyden–Fletcher–Goldfarb–Shanno method (BFGS) as an optimizer with exact
gradients Broyden (1970); Fletcher (1970); Goldfarb (1970); Shanno (1970). The
optimization takes place over piecewise-constant functions, which means that
the number of time-steps becomes an additional parameter of our numerics. This
algorithm begins with a random choice of initial pulse. Using BFGS, exact
gradients of the infidelity are calculated and parameters within the pulse are
varied as the infidelity goes down in the optimization space. This is done by
calculating second derivatives of the infidelity with respect to the
parameters, using this to approximate a parabola, and descending to the minima
of this approximation until reaching a tolerable minimum infidelity. The
limited memory version of this algorithm (L-BFGS) does not require calculation
of the entire Hessian matrix of second derivatives and so is not
computationally intensive Liu and Nocedal (1989).

###  IV.1 Optimized non-adiabatic pulse in a simple lattice

While the timescale of anyon creation through adiabatic evolution can be very
long, if instead of using linear ramps we use non-linear time-dependence in
the Hamiltonian which have been specially designed, then we can achieve high
fidelities at shorter timescales. The well-tested aforementioned GRAPE
algorithm is used to develop such time-dependent control functions also known
as pulses. The explicit Hamiltonians used in our optimizations incorporate
interaction factors $J$ of magnitude 1 with the effective magnetic field
interaction term equaling $J/100$, and thus all times are considered to be in
units of $J^{-1}$ to be applicable to the general case.

The results of using this procedure for a single plaquette of 6 spins are
depicted in Fig. 2 showing infidelities as function of the ramp time $T$. The
infidelities obtained with a linear ramp are depicted in blue. There is a
slight improvement with increasing ramp time, but the fidelity of about $90\%$
achieved with a ramp time of $T=1$ is only a very small improvement compared
to the initial fidelity at time $T=0$. This is consistent with an estimate
based on the spectral gap condition Van Dam  _et al._ (2001) that implies ramp
times $T\gg{\delta}^{-2}=3.480\ldots$ are required for high-fidelity
operations. With a value of $J=1$ this spectral gap is calculated from the
vortex gap of $\delta=|-4J-(-2\sqrt{3}J)|\approx 0.536$. This is a narrower
gap than the first fermionic excitation energy penalty which is $2J$. While in
general the fermionic excitation gap is a function of $J$, and $K$ and so can
be narrower than the vortex gap, in this special case of the single plaquette
lattice this is not the case.

The behaviour with optimized ramps, depicted in orange, is fundamentally
different. In the range $T<0.2$ there is a much faster decrease of infidelity
with increasing ramp time than in the case of linear ramps. This decrease is a
bit less pronounced in the range $0.2\lesssim T\lesssim 0.4$, but for ramp
time $T>0.4$ this decrease becomes increasingly pronounced with increasing
ramp time. For ramp times $T\gtrsim 0.8$, there is a rapid drop in the
infidelity, and for ramp times exceeding the threshold value of $T=T_{d}\simeq
0.85$, the deviations between the numerically obtained infidelities and the
ideal value of $0$ are consistent with noise due to finite numerical accuracy.

It is by no means surprising that even with an optimized ramp a finite ramp
time is required to reach perfect fidelities. This is due to the constant part
in the system Hamiltonian Eq. (10) that defines a natural time-scale of the
system. This effect is also referred to as the quantum speed limit Burgarth
_et al._ (2022); Deffner and Campbell (2017) and we will refer to the
threshold value $T_{d}$ of ramp durations at which fidelities drop to values
close to their ideal value as the drop time. Apart from limitations imposed by
finite-dimensional parametrization of the ramp function, the numerical
optimization routine and numerical accuracy, this drop time coincides with the
minimal duration required to reach perfect fidelity.

Figure 2: Logarithmic comparison of infidelity between initial and adiabatic
target state for linear control pulses and optimized control pulses, at
various timescales, for the simplest 6-qubit lattice. 100 time steps were used
in both cases. We can see more clearly the dramatic improvement in fidelity at
approximately $T=0.8$. The minimum infidelity reached by optimized pulses, on
the order of $10^{-9}$ is reached at time $T\approx 0.9$, many orders of
magnitude less than the time to reach this infidelity with the linear ramp,
which is at $T\approx 1350$. Figure 3: 6-qubit lattice: Heisenberg
infidelities based on differing targets (circles) are compared with their
corresponding state fidelities (squares). These are compared overall with
state fidelities achieved from linear ramp pulses (triangles). Different
colors represent different values of $T_{\text{ad}}$, which are the timescales
of 100 (red, e.g. upper squares), 200 (orange, middle squares) and 300 (green,
lower squares) that define the target unitary $W_{\text{ad}}$. All optimized
pulses are comprised of $200$ time steps. There is a significant improvement
at times between $T=2$ and $T=4$. Thereafter there is effectively perfect
fidelity with the presence of numerical noise. It confirms that both
Heisenberg and state infidelities improve markedly at the same time $T$ and
this is an improvement on the fidelities achieved with a linear ramp pulse.

The example of a single plaquette with $6$ qubits is also a good test-case to
compare optimization of state-fidelity and Heisenberg fidelity. Fig. 3 depicts
the Heisenberg fidelity obtained with various ramp functions as function of
the ramp time $T$.

The black triangles represent state fidelity data obtained with linear ramps.
Consistently with Fig. 2, there is only a moderate decrease of the infidelity
with increasing ramp time. The circles correspond to ramp functions optimized
for Heisenberg fidelity, and the different colors correspond to different
chosen adiabatic target times with $T=100$, $200$ and $300$ for red, orange
and green respectively. Similarly to the observations in Fig. 2 there is a
clear drop of the infidelities at a drop time $T_{d}\simeq 3$. The fact that
the numerically observed drop time is essentially the same in all three cases
indicates that the drop time is not dependent on the length of the adiabatic
target time that is chosen.

The squares depict the state infidelity obtained from the ramp functions that
had been optimized for Heisenberg infidelity, in order to see numerical
evidence that when implementing our procedure, good fidelity is achieved in
the one case which ensures a good infidelity in the other.

Here also is a clear drop of the infidelities and it occurs at the same drop
time as for Heisenberg infidelities. The fidelities obtained for ramps with a
longer ramp time than the drop time $T_{d}$, however, are not merely limited
by numerical accuracy, but they are indeed finite. Their exact value depends
on the parametrization of the ramp function, with finer parametrizations
resulting in lower infidelities. Since state fidelity and Heisenberg fidelity
are not strictly equivalent, it is not surprising, that a ramp that is
optimized for one of these fidelities does not yield the optimally achievable
value of the other fidelity. The results in Fig. 2, however, clearly show that
ramps optimized for Heisenberg fidelity result in high state fidelities and,
in particular, in infidelities that are between $3$ and $4$ orders of
magnitude lower than infidelities obtained with linear ramps.

Figure 4: 10-qubit lattice: Heisenberg infidelities between optimized
unitaries and an optimized target are shown (orange), as compared with
infidelities between the target and a unitary defined by a linear ramp control
(blue). Each optimized pulse consists of 200 time steps. The drop time is
$T\approx 10$, later than for the 6-qubit lattice. Figure 5: 30-qubit lattice
with open boundary conditions and pulses with 400 time steps. Heisenberg
infidelities between optimized unitaries and an optimized target are shown.
The spectral gap condition requires that $T_{\text{ad}}\gg 30$ calculated
analogously to the 6 and 10-qubit lattices. The results show that $T_{d}$ is
within the order of magnitude of 30 rather than many orders of magnitude
larger as required by adiabaticity, highlighting the improvement in timescale
achieved.

###  IV.2 Optimized pulses for larger lattices

Since numerical optimizations of state fidelity rapidly become infeasible with
increasing system size, the subsequent examples for larger systems feature
only Heisenberg fidelities with ramps that are optimized for this Heisenberg
fidelity. Fig. 4 shows the Heisenberg infidelity as function of ramp time for
linear ramps (blue) and for optimized ramps (orange) for a lattice consisting
of two adjacent plaquettes made up of 10-spins. Similarly to the cases
discussed above, there is a clearly identifiable drop time $T_{d}$, but its
value $T_{d}\simeq 10$ is larger than in the examples of smaller systems. The
abscissa depicts that ramp time on the log-scale, highlighting that linear
ramps with durations exceeding the drop time by many orders of magnitude are
required to achieve any sizable decrease in infidelity. The specific adiabatic
timescale is a consequence of the energy gap between vortex pairs. This is
calculated in the same way as for the single-plaquette system and with the
parameters being set at $J=1$ and $K=J/100$ this gap is $\delta=0.375\ldots$,
requiring a timescale of $T\gg\delta^{-2}\approx 7.070$.

The inset depicts a zoomed-in look into the domain around the drop time. It
highlights that, on top of the rapid drop of infidelity there is also a finite
noise level. When we compare the optimization results of a system made up of
10 spin-qubits and one with 30-qubits, whose optimized infidelities are shown
in Fig. 5, we see again the marked increase in drop time that is achieved. The
geometry of the 30-qubit spin lattice is identical to the lattice systems
presented in Fig. 7 and is subject to open boundary conditions. As detailed in
Appendix A, lattices with periodic boundary conditions require up to four
times as much optimization to be carried out compared with the same geometry
without periodicity therefore showing results for open boundary conditions is
the most efficient choice given limited computational resources. Here the
spectral gap condition, calculated analogously as with the 6 and 10-qubit
systems, requires the adiabatic timescale for this system to be
$T\gg\delta^{-2}\approx 30$. In the $30$-qubit lattice the vortex gap remains
smaller than the lowest fermionic energy gap but for larger systems this may
not necessarily be the case. Additionally the smallest gap may not be the
vortex gap as compared with the fermionic gap in scenarios where vortices are
created from a non-zero vortex sector.

Figure 6: Drop time (Above) and minimum required time steps (Below) compared
with lattice size by number of qubits. These quantities are defined as the
earliest time $T$ and minimum number of time steps in an optimized piecewise-
constant pulse respectively where Heisenberg infidelity drops below $10^{-6}$.
They are compared with quadratic growth curves (blue) to show sub-exponential
growth.

The observation that the drop time increases with the system size is depicted
more systematically in Fig. 6. Generally, the practically achievable drop time
depends on the parametrization of the ramp function and the number of free
parameters that can be optimized. Only for a ramp function with sufficiently
many piecewise constant elements is the drop time independent of the number of
free parameters. Fig. 6 depicts that shortest achievable drop time together
with the minimal number of piecewise constant elements in the ramp function
that is required to achieve this drop time. In both sub-plots the scaling is
consistent with a quadratic increase (blue line) with system size.

Within the validity of extrapolation from numerically accessible system sizes,
there is thus a very moderate scaling with system size, highlighting that
numerical optimizations based on the framework of fermionization can be
performed efficiently.

##  V Conclusions

It is well known that the Kitaev honeycomb model is a useful system for the
encoding of quantum information using non-abelian anyons. Our use of quantum
control techniques has allowed us to create anyons not only in adiabatic
timescales, but to overcome the necessity for such a restriction through the
use of gradient ascent pulse-engineering to construct non-linear ramps whose
timescales are many orders of magnitude faster than adiabatic. These ramps
achieve very high fidelities at these much faster timescales, and thus are
more useful for the encoding of quantum information and quantum algorithm
implementation given that decoherence times entail a preference for shorter
timescales.

Along with the drawback of requiring long timescales with adiabatic control,
the other chief drawback we would encounter with non-adiabatic quantum control
is the difficulty in completely determining the dynamics of honeycomb lattices
with large numbers of spin qubits. The method we have developed allows us to
overcome this obstacle by solving a matrix control problem where the matrices
grow linearly in the number of lattice spin-qubits, allowing for the
implementation of non-linear pulse-engineering and optimization of a related
matrix or Heisenberg fidelity. This, alongside the analytic expression
directly comparing state fidelity and Heisenberg fidelity allows us to be
assured of the success of these optimized ramps and allows us to have
confidence in the successful implementation of our procedure for use in
experimental realization of the model. Given the various avenues suggested for
experimental implementation of the Kitaev honeycomb model such as in solid
state materials Takagi  _et al._ (2019) and with polar molecules in optical
lattices Micheli  _et al._ (2006), as well as the success in observing
topological order in quantum spin liquids using the Rydberg blockade mechanism
Verresen  _et al._ (2021); Semeghini  _et al._ (2021) we believe our control
methods would be a useful for realization within a setup that allows for
varying spin-qubit interactions in a time-dependent manner.

Furthermore we have seen that there is indeed sub-exponential growth in the
computational difficulty of our control problem with growing system size and
so carrying out our methods on lattices on the order of hundreds of qubits is
possible, as necessary for scalable topological robustness.

## Acknowledgements

OR acknowledges funding from EPSRC Quantum Systems Engineering Hub and support
from MQCQE. DB acknowledges funding by the Australian Research Council
(project numbers FT190100106, DP210101367, CE170100009). Financial support
from the QuantERA ERANET Cofund in Quantum Technologies implemented within the
European Union’s Horizon 2020 Programme under the project Theory-Blind Quantum
Control TheBlinQC and from EPSRC under the grant EP/R044082/1 is gratefully
acknowledged

##  Appendix A Dynamics of the time-dependent Hamiltonian in the fermionic
picture

Previous work on the Kitaev honeycomb model has primarily been focused on
time-independent Hamiltonians, yet for quantum control we require time-
dependent Hamiltonians. Examples of such investigations of a time-dependent
Kitaev honeycomb model have included those focused on periodic driving Fulga
_et al._ (2019); Molignini  _et al._ (2021) which have used Floquet theory,
allowing for periodic driving to be mapped to an effectively time-independent
system Verdeny  _et al._ (2016). Other studies have used Jordan-Wigner
transformations without projections to look at specific time-dependent
behavior such as the Kitaev honeycomb model with a quenched magnetic field
Nasu2019NonequilibriumLiquids.

Here we show that in regular lattices with open, periodic or half-periodic
lattices, the fermionization procedure with projections is still possible in
the time-dependent case.

First let us recall the projections onto a link sector $u$

$\displaystyle P_{u}=\prod_{\\{jk\\}}^{L}\frac{1+u_{jk}\hat{u}_{jk}}{2},$ (20)

defined by an $L$-tuple of link eigenvalues $u_{jk}\in\\{\pm 1\\}$; the
projector onto a vortex sector $w$

$\displaystyle P_{w}=\prod_{j}^{P}\frac{1+w_{j}\hat{w}_{j}}{2}$ (21)

defined by a $P$-tuple of plaquette eigenvalues $w_{j}\in\\{\pm 1\\}$; and the
projector onto the physical subspace of the fermionic space given by

$\displaystyle P_{D}=\prod_{j}^{N}\frac{1+D_{j}}{2}\ ,\mbox{ with }\
D_{j}=b^{x}_{j}b^{y}_{j}b^{z}_{j}c_{j}.$ (22)

Since

$\displaystyle\hat{w}_{j}=\prod_{\\{k,l\\}\in w_{\beta}}\hat{u}_{kl},$ (23)

$P_{u}$ and $P_{w}$ commute, and moreover $w$ is fully determined by $u$. We
denote this relationship as $w=\omega(u)$ and thus have
$P_{w}P_{u}=\delta_{w,\omega(u)}P_{u}$ and

$\displaystyle P_{w}=P_{w}\sum_{u}P_{u}=\sum_{u:\omega(u)=w}P_{u}.$ (24)

$P_{w}$ commutes with $P_{D}$, $P_{u}$ and the time-dependent fermionic
Hamiltonian $H_{f}(t)$, so it will suffice to restrict ourselves to a single
plaquette sector $w$. While $P_{D}$ is the projection that determines
physicality, $P_{u}$ will turn the fermionized Hamiltonian into a quadratic
(and thereby easy to solve) one. A difficulty arises from the fact that
$P_{D}$ _does not_ commute with $P_{u}$. We can however find another useful
relationship between these projectors. Let $N$ be the number of qubits of the
original spin lattice and $\\{\Gamma_{k}|k=1,\ldots,2^{N}\\}$ be the set of
all possible products of stabilizer operators $D_{i}$ on the qubits, without
repetition, where we take an arbitrary but fixed order. For our fixed $w$
consider the corresponding pre-image $\omega^{-1}(w)$ of link sectors. We
define an equivalence relationship $\sim$ on this set by $u\sim
v:\Leftrightarrow\exists k:P_{u}=\Gamma_{k}P_{v}\Gamma_{k}.$ Let $\kappa$ be
the number of equivalence classes. Let us choose an arbitrary but fixed set of
representatives $u_{1},\ldots,u_{\kappa}$ and define
$P_{\bar{u}}=\sum_{i=1}^{\kappa}P_{u_{i}}.$

Upon expanding $P_{D}$ in terms of stabilizers, we obtain

$\displaystyle
P_{D}P_{\bar{u}}P_{D}=\frac{1}{2^{N}}P_{D}\sum_{i=1}^{\kappa}\sum_{k=1}^{2^{N}}\Gamma_{k}P_{u_{i}}\Gamma_{k}.$
(25)

To understand the right hand side better, we make a counting argument.
Firstly, it follows from the anticommutation relationships between link
operators and the $D_{j}$ that the $\Gamma_{k}P_{u_{i}}\Gamma_{k}$ are again
link projectors. Since we sum over all $\Gamma_{k}$ and by definition of the
equivalence classes, we know that every $u\in\omega^{-1}(w)$ appears at least
once on the right hand side, and that there are no overlaps between the
classes of fixed $i$. Furthermore the equality $\Gamma P_{u}\Gamma=P_{u}$
holds if and only if $\Gamma=1$ or $\Gamma=\prod_{i}D_{i}$. Therefore, there
are $2^{N-1}$ distinct projectors for each $i$. It follows that

$\displaystyle
P_{D}P_{\bar{u}}P_{D}=\frac{1}{2^{N}}2P_{D}\sum_{u:\omega(u)=w}P_{u}=\frac{1}{2^{N-1}}P_{D}P_{w}.$
(26)

From the commutativity relation $[H_{f}(t),P_{D}]=0$ it also follows that

$\displaystyle
V(t)P_{D}P_{w}=2^{N-1}VP_{D}P_{\bar{u}}P_{D}=2^{N-1}P_{D}P_{\bar{u}}V(t)P_{\bar{u}}P_{D},$
(27)

where $V(t)$ is the propagator corresponding to $H_{f}(t)$. Hence, the
evolution can be computed in the subspace $P_{\bar{u}}$. To conclude the
argument, we need to know the value of $\kappa$, as this determines how many
link sectors we need to consider. As long as $\kappa$ is not exponential, we
can efficiently simulate the dynamics.

To this end, we need another counting argument. To simplify the analysis, we
only consider three different regular lattice types dubbed open (o), periodic
(p) and half-periodic (h), and find their corresponding values of $\kappa$.

To do this, we first find relationships for the number of link operators $L$,
the number of plaquettes $P$, and the number of qubits $N$ for the various
lattices. Simple but rather tedious counting of such regular lattices shows
that $L-P=N-1$ in the open case and $L-P=N$ in the other two cases. Next,
compute the size of $\omega^{-1}(w).$ We show in the lemma below that
$|\omega^{-1}(w)|_{o,h}=2^{L-P}$ and $|\omega^{-1}(w)|_{p}=2^{L-P+1}.$ Since
each equivalence class has exactly $2^{N-1}$ elements, we have to have

$|\omega^{-1}(w)|=2^{N-1}\kappa$

such that $\kappa_{o}=1$, $\kappa_{h}=2$,$\kappa_{p}=4$.

######  Lemma A.1.

For all $w$, $|\omega^{-1}(w)|_{o,h}=2^{L-P}$ and
$|\omega^{-1}(w)|_{p}=2^{L-P+1}.$

###### Proof.

Consider first the case of a lattice with open boundary conditions. Since
there are no boundary constraints, all possible configuration of plaquette
eigenvalues $\\{\Lambda_{j}\\}$ are possible, and

$\displaystyle|\\{\Lambda_{j}\\}|_{o}=2^{P}.$ (28)

From the above, for each $u$, we can find $2^{N-1}$ other $v$ with
$\omega(v)=\omega(u)$ by conjugation with $\Gamma_{k}$. Since in the open
lattice the number of qubits $N$ follows the relation $N-1=L-P$, once can
conclude that for each $w$, the inequality $|\omega^{-1}(w)|_{o}\geq 2^{L-P}$
holds. Since there are by definition $2^{L}$ different link sectors, we have

$\displaystyle 2^{L}=\sum_{w}|\omega^{-1}(w)|_{o}\geq 2^{L-P}2^{P}=2^{L}$ (29)

so equality holds and the statement follows.

Next, consider a lattice with full periodic boundary conditions. Any link
eigenvalue change leads necessarily to exactly two plaquette eigenvalues being
flipped. Therefore only even numbers of vortices may ever be present, and the
number of plaquette eigenvalue configurations is

$\displaystyle|\\{\Lambda_{j}\\}|_{p}=2^{P-1}.$ (30)

Now, for each $u$, we can find $2^{N-1}$ other $v$ with $\omega(v)=\omega(u)$
by conjugation with $\Gamma_{k}$, but for each of these we can find 4
inequivalent link sectors. Since $L-P=N$ for periodic boundaries, we obtain

$\displaystyle 2^{L}=\sum_{w}|\omega^{-1}(w)|_{p}\geq 4\times
2^{L-P-1}2^{P-1}=2^{L}$ (31)

and the statement follows again.

Finally, in the half-periodic case we only have 2 inequivalent link sectors,
but $|\\{\Lambda_{j}\\}|_{p}=2^{P}$ is twice as big as in the periodic case,
and the same counting argument holds.

∎

#### Inequivalent sectors: Periodic boundary conditions

(a)

(b)

Figure 7: Simple sets of $y$ (above) and $z$ (below) link flips that preserve
the vortex sector while being gauge inequivalent to the trivial link sector.

(a)

(b)

Figure 8: Topological interpretation of the set of $y$ and $z$ flips that
preserve vortex sector but are gauge inequivalent to the trivial link sector
which topologically amounts to no loop at all.

We now consider the case of a lattice with periodic boundary conditions in
both the $a$ and $b$ directions (that is, both vertically and horizontally as
on a torus). We consider the no-vortex sector and and any arbitrary link
sector associated with it which we call $P_{u_{0}}$. We have already shown
that $P_{u_{0}}$ is gauge equivalent, that is, equivalent up to conjugation by
$\Gamma_{k}$ operators, with $2^{N-1}$ link sector projectors.

Acting on a non-zero eigenstate of $P_{u_{0}}$ with any stabilizer operator
$D_{j}$ will change three link sector eigenvalues at a time, due to $D_{j}$
overlapping with the three types of link that contain spin site $j$, the same
is therefore also true of $\Gamma_{k}$ operators. On the other hand acting on
spin sites with a Pauli operator such as $Z_{j}$, which in the fermionic
picture is

$\displaystyle\tilde{Z}_{j}=ib^{z}_{j}c_{j}$ (32)

will flip the sign of only the $z$-link corresponding to that site and no
others. In general a Pauli operator $\tilde{\sigma}^{\alpha}$ will only flip
the sign of the $\alpha$-link containing spin site $j$. This is due to the
fact that

$\displaystyle\\{\tilde{\sigma}_{j},\hat{u}_{jk}\\}$
$\displaystyle=\\{ib^{\sigma}_{j}c_{j},ib_{j}^{\sigma}b_{k}^{\sigma}\\}=0$
(33) $\displaystyle\left[\tilde{\sigma}_{j},\hat{u}_{kl}\right]$
$\displaystyle=\left[ib^{\sigma}_{j}c_{j},ib_{k}^{\sigma}b_{l}^{\sigma}\right]=0,~{}~{}~{}j\notin\\{k,l\\}.$
(34)

Thus we have that

$\displaystyle\tilde{Z}_{\alpha}P_{u_{j}}=P_{u_{k}}\tilde{Z}_{\alpha}$ (35)

and thus

$\displaystyle P_{u_{j}}=\tilde{Z}_{\alpha}P_{u_{k}}\tilde{Z}_{\alpha}$ (36)

where in general $\omega(u_{j})\neq\omega(u_{k})$.

Flipping a single link eigenvalue will change the vortex sector as a pair of
vortices are introduced on adjacent plaquettes. Taking $\zeta$ to be an
arbitrary product of $\tilde{Z}$ Pauli operators, then in order to have

$\displaystyle P_{u_{j}}=\zeta P_{u_{k}}\zeta$ (37)

such that $\omega(u_{j})$ does equal $\omega(u_{k})$, we require $\zeta$ to
consist of a pair of link-flipping $\tilde{Z}$ operators for each plaquette.

In order to return to the no-vortex sector, vortices need to be annihilated
and so the vortex string must form a closed loop. Algebraically speaking,
while each $\tilde{\sigma}_{j}$ will commute with all overlapping and non-
overlapping $D_{j}$ operators, as

$\displaystyle\left[\tilde{\sigma}_{j},D_{j}\right]$
$\displaystyle=\left[ib^{\alpha}_{j}c_{j},b_{j}^{x}b_{j}^{y}b_{j}^{z}c_{j}\right]$
(38) $\displaystyle=0$ (39)

and

$\displaystyle\left[\tilde{\sigma}_{j},D_{k}\right]$
$\displaystyle=\left[ib^{\alpha}_{j}c_{j},b_{k}^{x}b_{k}^{y}b_{k}^{z}c_{k}\right]$
(40) $\displaystyle=0,$ (41)

such Paulis will not commute with overlapping plaquette operators $W_{p}$ as
for a spin site $j$, plaquette operators will act with Paulis $\tilde{X}$ or
$\tilde{Y}$ and thus there is anticommutation.

Figure 9: Clearly
$\left[Z_{6},W_{p}\right]=\left[Z_{6},Z_{1}X_{2}Y_{3}Z_{4}X_{5}Y_{6}\right]\neq
0$ but
$\left[Z_{6}Z_{2},W_{p}\right]=\left[Z_{6}Z_{2},Z_{1}X_{2}Y_{3}Z_{4}X_{5}Y_{6}\right]=0$.

In order for there to be commutation with all plaquette operators, there needs
to be $\tilde{Z_{j}}$ operators acting on two spin sites per plaquette, as
shown in Fig. 9. The simplest example of a product of $\tilde{Z}_{j}$
operators that commutes with all link operators is shown explicitly and
schematically in Fig. 7(b) and Fig. 8(b) respectively.

Figure 10: A lattice where the grey spin sites have been acted on with $D_{j}$
operators and link eigenvalues on yellow links have been flipped. It is not
possible to flip only a single line of $z$-links using $\Gamma_{k}$ operators
alone.

Starting with link sector $P_{u_{0}}$ and using only $\Gamma_{k}$ operators,
it is not possible to flip the sign of only a single row of $z$-links as in
general the action of a $D_{j}$ operator flips the sign of all three types of
link $x,y$ and $z$. In order to keep all $x$ and $y$ links unflipped, products
of $D_{j}$ operators must act on the sites of all spins that correspond to
those links, leaving at least a pair of rows with flipped eigenvalues on
$z$-links, as demonstrated in Fig. 10. We therefore have found another link
sector $P_{u_{z}}$ whose corresponding vortex sector is the no-vortex sector
which cannot be reached from $P_{u_{0}}$ by conjugation with $\Gamma_{k}$
operators. Algebraically this can be represented by

$\displaystyle\nexists~{}\Gamma_{k}~{}\text{such
that}~{}\Gamma_{k}P_{u_{0}}\Gamma_{k}=P_{u_{z}}$ (42)

where

$\displaystyle P_{u_{z}}:=\zeta P_{u_{0}}\zeta$ (43)

and $\zeta$ is a horizontal, topologically closed loop of Pauli $\tilde{Z}$
operators. This is because when commuting $\zeta$ through $P_{u_{0}}$, only
the signs of a row of $z$-link eigenvalues are changed. If we call this row of
links $\rho$ then

$\displaystyle\zeta P_{u_{0}}$
$\displaystyle=\zeta\prod_{jk}\frac{1+u^{\alpha}_{jk}\hat{u}^{\alpha}_{jk}}{2}$
(44)
$\displaystyle=\prod_{\alpha=z,{jk}\in\rho}\frac{1-u^{\alpha}_{jk}\hat{u}^{\alpha}_{jk}}{2}\prod_{\text{etc.}}\frac{1+u^{\alpha}_{jk}\hat{u}^{\alpha}_{jk}}{2}\zeta$
(45) $\displaystyle=P_{u_{k}}\zeta.$ (46)

However with $\Gamma_{k}$ operators, in order to commute them through
$P_{u_{0}}$ and have only $z$-link eigenvalues flipped, this requires a
$\Gamma_{k}$ consisting of $D_{j}$ operators acting on a row of spin sites in
the manner of Fig. 10 and thus a minimum of two rows of $z$-link eigenvalues
are flipped. Thus for all $\Gamma_{k}$

$\displaystyle\Gamma_{k}P_{u_{0}}\Gamma_{k}\neq\zeta P_{u_{0}}\zeta.$ (47)

We can now repeat the above process, with the same reasoning, starting with
the link sector $P_{u_{z}}$ rather than $P_{u_{0}}$, but this time with Pauli
$\tilde{Y}$ operators to flip links as shown in Fig. 7(a) and Fig. 8(a). This
would have given another link sector unreachable via actions of $\Gamma_{k}$
operators, just as $P_{u_{z}}$ was unreachable from $P_{u_{0}}$. This would
have corresponded to flipping link signs vertically rather than horizontally.
Pauli $\tilde{X}$ operators can be written, up to a phase, as simply products
of $\tilde{Z}$ and $\tilde{Y}$ Paulis, thus no new link sectors can be found
by use of the above operations with $\tilde{X}$ Paulis. Therefore for a
lattice defined by $a$ plaquettes per row and made up of $b$ rows, if we
define

$\displaystyle\zeta$ $\displaystyle=\tilde{Z}_{j_{1}}...\tilde{Z}_{j_{a}}$
(48) $\displaystyle\chi$ $\displaystyle=\tilde{Y}_{j_{1}}...\tilde{Y}_{j_{b}}$
(49)

then the four link sector operators

$\displaystyle P_{u_{0}}$ , (50) $\displaystyle P_{u_{z}}$
$\displaystyle=\zeta P_{u_{0}}\zeta,$ (51) $\displaystyle P_{u_{y}}$
$\displaystyle=\chi P_{u_{0}}\chi,$ (52) $\displaystyle P_{u_{x}}$
$\displaystyle=\zeta\chi P_{u_{0}}\chi\zeta$ (53)

each define gauge-inequivalent link sectors in a particular vortex sector.

#### Inequivalent sectors: Half-periodic boundary conditions

For a system with half-periodic boundary conditions, depending on whether
periodicity is in the horizontal or vertical direction, only a single row or
‘column’ of flipped signs would be required to find a link sector projector
not reachable from $P_{u_{0}}$ by actions of $\Gamma_{k}$ operators. There
would therefore be two sets of $2^{N-1}$ link sector projectors corresponding
to two equivalence classes in each vortex sector. Accounting for all $2^{P}$
vortex sectors there would be a total of

$\displaystyle 2^{P}(2\times 2^{N-1})=2^{N+P}=2^{L}$ (54)

and so we have accounted for all link sectors.

##  Appendix B Proof of Heisenberg-infidelity state-infidelity bound

This section contains the proof of Eq (19) which relates state infidelity and
Heisenberg infidelity. State infidelity is defined between two states with
differing dynamics $U_{j}~{}(j=1,2)$ applied to an initial spin state
$\ket{\Psi_{I}}$. These dynamics are induced by spin Hamiltonians $H_{s,j}$
whose counterparts in the quadratic Majorana fermionic picture
$\tilde{H}_{\text{ferm},j}=\alpha^{\dagger}M_{j}\alpha$ (see Section. II.1)
define orthogonal matrices $\mathcal{O}_{j}$ generated by the matrices $M_{j}$
in a relation which satisfies the differential equation
$\dot{\mathcal{O}}_{j}=iM_{j}\mathcal{O}_{j}$.

The relation to be proven reads

$\displaystyle\mathcal{I}_{H}(\mathcal{O}_{1},\mathcal{O}_{2})\geq\frac{1}{4d^{3}}\left(1-\sqrt{\mathcal{F}_{s}(\Psi_{1},\Psi_{2}})\right)$
(55)

with the state fidelity

$\displaystyle\mathcal{F}_{s}(\Psi_{1},\Psi_{2})\
=\left|\langle\Psi_{1}|\Psi_{2}\rangle\right|^{2}\
=\left|\bra{\Psi_{I}}U_{2}^{\dagger}U_{1}\ket{\Psi_{I}}\right|^{2},$ (56)

and the Heisenberg infidelity $\mathcal{I}_{H}=1-\mathcal{F}_{H}$ with

$\displaystyle\mathcal{F}_{H}(\mathcal{O}_{1},\mathcal{O}_{2})=\frac{1}{d}\left|\Tr(\mathcal{O}_{1}^{\dagger}\mathcal{O}_{2})\right|\
,$ (57)

and the matrix dimension $d$.

The full proof consists of Eq. (58) and the series of inequalities Eq. (59) to
(61)

$\displaystyle\sqrt{2d}\sqrt{\mathcal{I}_{H}}$
$\displaystyle=\|\mathcal{O}_{1}-\mathcal{O}_{2}\|_{F}$ (58)
$\displaystyle\geq\frac{1}{2d}\|\Phi_{U_{1}}-\Phi_{U_{2}}\|_{\lozenge}$ (59)
$\displaystyle\geq\frac{1}{2d}\min_{\varphi}\|U_{1}-e^{i\varphi}U_{2}\|_{\text{op}}$
(60) $\displaystyle\geq\frac{1}{\sqrt{2}d}\sqrt{1-\sqrt{\mathcal{F}_{s}}}$
(61)

that will be discussed separately in the following subsections.

Eq. (58) is expressed in terms of the Frobenius norm

$\displaystyle\|A\|_{F}:=\sqrt{\Tr(A^{\dagger}A)}$ (62)

for any operator $A$. Eq. (59) is expressed in terms of the operator norm

$\displaystyle\|A\|_{\text{op}}$
$\displaystyle:=\sup\left\\{\frac{\|Ax\|}{\|x\|}:x\in V^{d}\text{ with }x\neq
0\right\\}.$ (63)

The diamond norm for a quantum channel $\Phi$ in Eq. (60) is given by

$\displaystyle\|\Phi\|_{\lozenge}:=\max_{\rho}\|(\Phi\otimes\mathds{1})\rho\|_{1}$
(64)

where $\|\circ\|_{1}=\Tr\sqrt{A^{\dagger}A}$ is the trace norm Benenti and
Strini (2010), and the maximization is taken over all density matrices in a
space of dimension corresponding to the size of the quantum channel.

###  B.1  (58)

The Frobenius norm of the difference between two orthogonal operators
$\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ reduces to

$\displaystyle\|\mathcal{O}_{1}-\mathcal{O}_{2}\|_{F}$
$\displaystyle=\sqrt{2d-\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})-\Tr(\mathcal{O}_{2}^{\mathrm{T}}\mathcal{O}_{1})}\
$ (65)
$\displaystyle=\sqrt{2d-2\Re{\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})}},$
(66)

where $d$ is the dimension of $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$. It thus
depends on the real part of $\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})$
and not on its absolute value as it is the case for $\mathcal{F}_{H}$. As we
will show in the following, however, in the present case, the object
$\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})$ is real and positive, so
that Eq. (58) is indeed satisfied.

####  B.1.1 Proof that $\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})$ is
real

The orthogonal matrices $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ satisfy the
differential equation $\dot{\mathcal{O}}_{j}=iM_{j}\mathcal{O}_{j}$ with
generally time-dependent generators $iM_{j}$. Since the $M_{j}$ are purely
imaginary, the generators $iM_{j}$ are purely real. Together with the initial
condition $\mathcal{O}_{j}(0)={\mathds{1}}$, such that $\mathcal{O}_{j}(0)$ is
real, this implies that $\mathcal{O}_{j}(t)$ for $j=1,2$ is real for all
times. Consequently the overlap
$\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})$ is also real.

####  B.1.2 Proof that $\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})$ is
non-negative

Since Majorana fermions move between only one of two pairs of fermionic sites
per spin site, the full space that $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ act
on, can be divided into two subspaces $H_{x}$ and $H_{y}$ of equal dimension
$d/2$. Both $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ are given as a direct sum
of the identity ${\mathds{1}}$ in $H_{x}$ and orthogonal operators
$\tilde{\mathcal{O}}$ in $H_{y}$. The complete trace is thus given as the sum
of the two traces $\Tr_{x}{\mathds{1}}$ and
$\Tr_{y}\tilde{\mathcal{O}}_{1}^{\mathrm{T}}\tilde{\mathcal{O}_{2}}$. The
latter trace can also be expressed as the sum over the eigenvalues
$\lambda_{j}$ of
$\tilde{\mathcal{O}}_{1}^{\mathrm{T}}\tilde{\mathcal{O}_{2}}$. This results in
the relation

$\displaystyle\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})=\Tr_{x}{\mathds{1}}+\Tr_{y}\tilde{\mathcal{O}}_{\text{1}}^{\mathrm{T}}\tilde{\mathcal{O}_{2}}=\frac{d}{2}+\sum_{j=1}^{\frac{d}{2}}\lambda_{j}\
.$ (67)

Since the trace $\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})$ is purely
real, the imaginary parts $\imaginary\lambda_{j}$ add up to zero. Because all
the eigenvalues $\lambda_{j}$ are phase factors, i.e.
$\lambda_{j}=\exp(i\varphi_{j})$ with $\varphi_{j}$ real, the real parts of
the $\lambda_{j}$ satisfy the inequality $\real\lambda_{j}\geq-1$, so that the
relation

$\displaystyle\Tr(\mathcal{O}_{1}^{\mathrm{T}}\mathcal{O}_{2})=\frac{d}{2}+\sum_{j=1}^{\frac{d}{2}}\lambda_{j}\geq
0$ (68)

is indeed given.

###  B.2 from (58) to (59)

The Frobenius and operator matrix norms satisfy the inequality
$\|A\|_{F}\geq\|A\|_{\text{op}}$ Golub and Van Loan (1996). For the present
case, this implies the relation

$\|\mathcal{O}_{1}-\mathcal{O}_{2}\|_{F}\geq\|\mathcal{O}_{1}-\mathcal{O}_{2}\|_{\text{op}}\
.$ (69)

Together with the relation

$\displaystyle\|\mathcal{O}_{1}-\mathcal{O}_{2}\|_{\text{op}}\geq\frac{1}{2d}\|\Phi_{U_{1}}-\Phi_{U_{2}}\|_{\lozenge}\
,$ (70)

that is proven in Eq. (171) in Oszmaniec  _et al._ (2022a) this yields the
desired inequality.

###  B.3 from (59) to (60)

The required inequality

$\displaystyle\|\Phi_{U_{1}}-\Phi_{U_{2}}\|_{\lozenge}$
$\displaystyle\geq\min_{\varphi}\|U_{1}-e^{i\varphi}U_{2}\|_{\text{op}}$ (71)

is proven in Eq. (2.1) in Oszmaniec  _et al._ (2022b)

###  B.4 from (60) to (61)

Following the definition of the operator norm (Eq. (63)), the operator norm in
Eq. (60) satisfies the inequality

$\displaystyle\|U_{1}-e^{i\varphi}U_{2}\|_{\text{op}}\geq\|(U_{1}-e^{i\varphi}U_{2})\ket{\Psi}\|_{2}$
(72)

for any state vector $\ket{\Psi}$ in the spin picture, and, as such, in
particular for the initial state $\ket{\Psi_{I}}$ of the dynamics. That is,
the inequality

$\displaystyle\|U_{1}-e^{i\varphi}U_{2}\|_{\text{op}}$
$\displaystyle\geq\|\ket{\Psi_{1}}-e^{i\varphi}\ket{\Psi_{2}}\|_{2}$ (73)
$\displaystyle=\sqrt{(\bra{\Psi_{1}}-e^{-i\varphi}\bra{\Psi_{2}})(\ket{\Psi_{1}}-e^{i\varphi}\ket{\Psi_{2}}})$
(74)
$\displaystyle=\sqrt{2-e^{i\varphi}\innerproduct{\Psi_{1}}{\Psi_{2}}-e^{-i\varphi}\innerproduct{\Psi_{2}}{\Psi_{1}}}$
(75)
$\displaystyle=\sqrt{2-2|\langle\Psi_{1}|\Psi_{2}\rangle|\cos(\varphi+\theta)}$
(76)

holds, where the overlap of $\ket{\Psi_{1}}$ and $\ket{\Psi_{2}}$ can be
defined in terms of state infidelity and a phase

$\displaystyle\langle\Psi_{1}|\Psi_{2}\rangle=|\langle\Psi_{1}|\Psi_{2}\rangle|e^{i\theta}=\sqrt{\mathcal{F}_{s}}e^{i\theta}.$
(77)

When Eq. (76) is minimized over all phases $\varphi$ this gives the required
result

$\displaystyle\min_{\varphi}\|U_{1}-e^{i\varphi}U_{2}\|_{\text{op}}$
$\displaystyle\geq\sqrt{2}\sqrt{1-\sqrt{\mathcal{F}_{s}(\Psi_{1},\Psi_{2}})}.$
(78)

## References

  * Kitaev (2006) A. Yu. Kitaev, “Anyons in an exactly solved model and beyond,” Annals of Physics 321, 2–111 (2006).
  * Wilczek (1982) Frank Wilczek, “Quantum mechanics of fractional-spin particles,” Physical Review Letters 49, 957–959 (1982).
  * Sarma  _et al._ (2006) S. Das Sarma, Michael Freedman, and Chetan Nayak, “Topological quantum computation,” Physics Today 59, 32–38 (2006).
  * Pachos (2012) Jiannis K. Pachos,  _Introduction to topological quantum computation_ (Cambridge University Press, 2012).
  * Kitaev (2003) A. Yu. Kitaev, “Fault-tolerant quantum computation by anyons,” Annals of Physics 303, 2–30 (2003).
  * Brennen and Pachos (2007) Gavin K. Brennen and Jiannis K. Pachos, “Why should anyone care about computing with anyons?” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 464, 1–24 (2007).
  * Knolle and Moessner (2019) J. Knolle and R. Moessner, “A Field Guide to Spin Liquids,” Annual Review of Condensed Matter Physics 10, 451–472 (2019).
  * Takagi  _et al._ (2019) Hidenori Takagi, Tomohiro Takayama, George Jackeli, Giniyat Khaliullin, and Stephen E. Nagler, “Concept and realization of Kitaev quantum spin liquids,” Nature Reviews Physics 2019 1:4 1, 264–280 (2019).
  * Lahtinen (2011) Ville Lahtinen, “Interacting non-Abelian anyons as Majorana fermions in the honeycomb lattice model,” New Journal of Physics 13, 075009 (2011).
  * Génetay Johansen and Simula (2021) Emil Génetay Johansen and Tapio Simula, “Fibonacci Anyons Versus Majorana Fermions: A Monte Carlo Approach to the Compilation of Braid Circuits in SU(2) k Anyon Models,” PRX Quantum 2, 10334 (2021).
  * Cheng  _et al._ (2011) Meng Cheng, Victor Galitski, and S. Das Sarma, “Nonadiabatic effects in the braiding of non-Abelian anyons in topological superconductors,” Physical Review B 84, 104529 (2011).
  * Lai and Zhang (2020) Hon Lam Lai and Wei Min Zhang, “Decoherence dynamics of Majorana qubits under braiding operations,” Physical Review B 101, 195428 (2020).
  * He  _et al._ (2016) Shuang He, Shi Lei Su, Dong Yang Wang, Wen Mei Sun, Cheng Hua Bai, Ai Dong Zhu, Hong Fu Wang, and Shou Zhang, “Efficient shortcuts to adiabatic passage for three-dimensional entanglement generation via transitionless quantum driving,” Scientific Reports 6, 1–11 (2016).
  * Zhou  _et al._ (2017) Brian B. Zhou, Alexandre Baksic, Hugo Ribeiro, Christopher G. Yale, F. Joseph Heremans, Paul C. Jerger, Adrian Auer, Guido Burkard, Aashish A. Clerk, and David D. Awschalom, “Accelerated quantum control using superadiabatic dynamics in a solid-state lambda system,” Nature Physics 13, 330–334 (2017).
  * Xu  _et al._ (2019) Kebiao Xu, Tianyu Xie, Fazhan Shi, Zhen Yu Wang, Xiangkun Xu, Pengfei Wang, Ya Wang, Martin B. Plenio, and Jiangfeng Du, “Breaking the quantum adiabatic speed limit by jumping along geodesics,” Science Advances 5 (2019).
  * Wang  _et al._ (2012) W. Wang, S. C. Hou, and X. X. Yi, “Adiabatic evolution under quantum control,” Annals of Physics 327, 1293–1303 (2012).
  * Simsek and Mintert (2021) Selwyn Simsek and Florian Mintert, “Quantum control with a multi-dimensional Gaussian quantum invariant,” Quantum 5, 409 (2021).
  * Feng  _et al._ (2007) Xiao Yong Feng, Guang Ming Zhang, and Tao Xiang, “Topological characterization of quantum phase transitions in a spin-1/2 model,” Physical Review Letters 98, 087204 (2007).
  * Kells  _et al._ (2009) G. Kells, J. K. Slingerland, and J. Vala, “Description of Kitaev’s honeycomb model with toric-code stabilizers,” Physical Review B - Condensed Matter and Materials Physics 80, 125415 (2009).
  * Brennan and Vala (2018) John Brennan and Jirí Vala, “The Kitaev honeycomb model on surfaces of genus g $\geq$ 2,” New Journal of Physics 20, 053023 (2018).
  * Kells  _et al._ (2010) G. Kells, D. Mehta, J. K. Slingerland, and J. Vala, “Exact results for the star lattice chiral spin liquid,” Physical Review B - Condensed Matter and Materials Physics 81, 104429 (2010).
  * Knolle (2016) Johannes Knolle,  _Dynamics of a Quantum Spin Liquid_, Springer Theses (Springer International Publishing, Cham, 2016).
  * Lahtinen and Pachos (2017) Ville T. Lahtinen and Jiannis K. Pachos, “A short introduction to topological quantum computation,” SciPost Physics 3, 021 (2017).
  * Lieb (1994) Elliott H. Lieb, “Flux Phase of the Half-Filled Band,” Physical Review Letters 73, 2158–2161 (1994).
  * Blaizot and Ripka (1986) Jean-Paul Blaizot and Georges Ripka,  _Quantum theory of finite systems_ (MIT Press, Cambridge Mass., 1986).
  * Johansson  _et al._ (2013) J. R. Johansson, P. D. Nation, and Franco Nori, “QuTiP 2: A Python framework for the dynamics of open quantum systems,” Computer Physics Communications 184, 1234–1240 (2013).
  * Khaneja  _et al._ (2005) Navin Khaneja, Timo Reiss, Cindie Kehlet, Thomas Schulte-Herbrüggen, and Steffen J. Glaser, “Optimal control of coupled spin dynamics: Design of NMR pulse sequences by gradient ascent algorithms,” Journal of Magnetic Resonance 172, 296–305 (2005).
  * Rowland and Jones (2012) Benjamin Rowland and Jonathan A. Jones, “Implementing quantum logic gates with gradient ascent pulse engineering: principles and practicalities,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 4636–4650 (2012).
  * Broyden (1970) C. G. Broyden, “The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations,” IMA Journal of Applied Mathematics 6, 76–90 (1970).
  * Fletcher (1970) Roger Fletcher, “A new approach to variable metric algorithms,” The Computer Journal 13, 317–322 (1970).
  * Goldfarb (1970) Donald Goldfarb, “A family of variable-metric methods derived by variational means,” Mathematics of Computation 24, 23–26 (1970).
  * Shanno (1970) D. F. Shanno, “Conditioning of quasi-Newton methods for function minimization,” Mathematics of Computation 24, 647–656 (1970).
  * Liu and Nocedal (1989) Dong C. Liu and Jorge Nocedal, “On the limited memory BFGS method for large scale optimization,” Mathematical Programming 1989 45:1 45, 503–528 (1989).
  * Van Dam  _et al._ (2001) Wim Van Dam, Michele Mosca, and Umesh Vazirani, “How powerful is adiabatic quantum computation?” Annual Symposium on Foundations of Computer Science - Proceedings , 279–287 (2001).
  * Burgarth  _et al._ (2022) Daniel Burgarth, Jeff Borggaard, and Zoltán Zimborás, “Quantum distance to uncontrollability and quantum speed limits,” Physical Review A 105, 042402 (2022).
  * Deffner and Campbell (2017) Sebastian Deffner and Steve Campbell, “Quantum speed limits: from Heisenberg’s uncertainty principle to optimal quantum control,” Journal of Physics A: Mathematical and Theoretical 50, 453001 (2017).
  * Micheli  _et al._ (2006) A. Micheli, G. K. Brennen, and P. Zoller, “A toolbox for lattice-spin models with polar molecules,” Nature Physics 2006 2:5 2, 341–347 (2006).
  * Verresen  _et al._ (2021) Ruben Verresen, Mikhail D. Lukin, and Ashvin Vishwanath, “Prediction of Toric Code Topological Order from Rydberg Blockade,” Physical Review X 11, 031005 (2021).
  * Semeghini  _et al._ (2021) G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, “Probing topological spin liquids on a programmable quantum simulator,” Science 374, 1242–1247 (2021).
  * Fulga  _et al._ (2019) I. C. Fulga, M. Maksymenko, M. T. Rieder, N. H. Lindner, and E. Berg, “Topology and localization of a periodically driven Kitaev model,” Physical Review B 99, 235408 (2019).
  * Molignini  _et al._ (2021) Paolo Molignini, Albert Gasull Celades, R. Chitra, and Wei Chen, “Crossdimensional universality classes in static and periodically driven Kitaev models,” Physical Review B 103, 184507 (2021).
  * Verdeny  _et al._ (2016) Albert Verdeny, Joaquim Puig, and Florian Mintert, “Quasi-periodically driven quantum systems,” Zeitschrift fur Naturforschung - Section A Journal of Physical Sciences 71, 897–907 (2016).
  * Benenti and Strini (2010) Giuliano Benenti and Giuliano Strini, “Computing the distance between quantum channels: Usefulness of the Fano representation,” Journal of Physics B: Atomic, Molecular and Optical Physics 43, 215508 (2010).
  * Golub and Van Loan (1996) Gene H. Golub and Charles F. Van Loan,  _Matrix computations_ , 3rd ed. (The Johns Hopkins University Press, 1996) pp. 56–57.
  * Oszmaniec  _et al._ (2022a) Michał Oszmaniec, Ninnat Dangniam, Mauro E.S. Morales, and Zoltán Zimborás, “Fermion Sampling: A Robust Quantum Computational Advantage Scheme Using Fermionic Linear Optics and Magic Input States,” PRX Quantum 3, 020328 (2022a).
  * Oszmaniec  _et al._ (2022b) Michał Oszmaniec, Adam Sawicki, and Michał Horodecki, “Epsilon-Nets, Unitary Designs, and Random Quantum Circuits,” IEEE Transactions on Information Theory 68, 989–1015 (2022b).