# Non-perturbative Floquet engineering of the toric-code Hamiltonian and its
ground state

Francesco Petiziol<EMAIL_ADDRESS>Technische Universität Berlin,
Institut für Theoretische Physik, Hardenbergstraße 36, Berlin 10623, Germany
Sandro Wimberger  Dipartimento di Scienze Matematiche, Fisiche e Informatiche,
Università di Parma, Parco Area delle Scienze 7/A, 43124 Parma, Italy  INFN,
Sezione di Milano Bicocca, Gruppo Collegato di Parma, Parco Area delle Scienze
7/A, 43124 Parma, Italy  André Eckardt  Technische Universität Berlin,
Institut für Theoretische Physik, Hardenbergstraße 36, Berlin 10623, Germany
Florian Mintert  Blackett Laboratory, Imperial College London, London SW7 2AZ,
United Kingdom  Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400,
01328 Dresden, Germany

###### Abstract

We theoretically propose a quantum simulation scheme for the toric-code
Hamiltonian, the paradigmatic model of a quantum spin liquid, based on time-
periodic driving. We develop a hybrid continuous-digital strategy that
exploits the commutativity of different terms in the target Hamiltonian. It
allows one to realize the required four-body interactions in a nonperturbative
way, attaining strong coupling and the suppression of undesired processes. In
addition, we design an optimal protocol for preparing the topologically
ordered ground states with high fidelity. A minimal implementation of a
topological device and the tools to manipulate it are also discussed. The
proposed scheme finds natural implementation in architectures of
superconducting qubits with tuneable couplings.

##  I Introduction

The discovery of the fractional quantum Hall effect [1], together with later
theoretical formulations of quantum spin liquids [2], unraveled the existence
of phases of matter, eluding traditional classifications based on symmetry
breaking and local order parameters, which can be understood instead through
the notion of topological order [3, 4]. Topologically ordered states are
characterized by properties such as long-range entanglement, ground state
degeneracy on topologically nontrivial manifolds, intrinsic robustness to
perturbations, and they support quasiparticles with anyonic quantum statistics
[5, 6]. Such features, alongside their importance from a fundamental
condensed-matter physics perspective, make topologically ordered states
promising candidates for scalable quantum information processing [7, 8, 9].

The paradigmatic model in the understanding of topological order and of its
potential for quantum computation is Kitaev’s toric-code Hamiltonian [10].
This model describes spin-1/2 systems in a two-dimensional (2D) square
lattice, experiencing purely four-spin interactions. It underpins surface
codes for fault-tolerant quantum computing [11, 12, 7] and is tightly bound to
$\mathbb{Z}_{2}$ lattice gauge theory [13, 14, 15, 16] and string-net
condensation [17, 6]. The development of quantum simulators and quantum
processors has enabled to probe signatures of $\mathbb{Z}_{2}$ topological
order in experiments [18, 19, 20, 21, 22, 23]. For example, the toric-code
ground state has been recently prepared in a superconducting quantum processor
via a quantum-circuit-based approach [23, 24], without implementing the
background Hamiltonian. However, due to the difficulty to attain clean four-
spin interactions in synthetic systems [25, 26, 27], the realization of the
full Kitaev’s model is still a challenge.

Figure 1: Sketch of the spin lattice and of the quantum simulation approach,
which will be detailed in the following. Square symbols indicate spins, solid
lines represent tuneable pairwise couplings, and wavy purple lines represent
couplings that are driven in a given time substep $\tau$, with associated
angular frequency $\omega=2\pi/\tau$. In each substep $\tau$, different
subsystems are decoupled from each other and driven to produce effective four-
spin interactions (light-blue squares). At the end of the sequence (which is
periodic with period $T=4\tau$), a complete toric-code Hamiltonian is achieved
with high accuracy. The proposed scheme further allows one to prepare its
ground state(s) and to study anyonic quasiparticles (‘$\mathrm{e}$’ and
‘$\mathrm{m}$’).

In this work, we propose a scheme for the accurate quantum simulation of the
toric-code Hamiltonian via periodic driving (sketched in Fig. 1). Adopting the
perspective of recently proposed hybrid approaches [28, 29], we exploit the
interplay of techniques from continuous-drive Floquet engineering [30, 31, 32,
33] and digitization [34] to minimize undesired terms in the quantum simulated
Hamiltonian, achieving clean four-spin interactions. By taking advantage of
the commutativity of different four-spin terms characterizing the Hamiltonian,
we can reduce the problem of finding suitable control functions to individual
four-spin subsystems, which can be optimized based on numerically exact
methods, thus circumventing perturbative treatments. Moreover, this approach
provides, by construction, a method that can be applied straightforwardly to
arbitrary lattice size. Thanks to the fact that the proposed scheme exploits
the integrability of the target model, the level of accuracy attained for the
single plaquette is maintained when addressing the whole lattice, such that no
amplification of the quantum simulation error occurs. It is further shown how
the scheme can be adapted to implement four-spin entangling gates, which, in
turn, allow us to design systematic ground-state-preparation protocols
applicable in arbitrarily large systems. The achievement of topologically
ordered states in these protocols is verified by probing long-range
entanglement, detected via topological entanglement entropy, and through
explicit creation and braiding of anyons. A minimal, proof-of-concept
implementation is then proposed, which can be used either to study the
transition to the ordered phase, through Floquet-adiabatic passage, or as a
prototypical topological qubit.

Our scheme applies to a lattice of driven two-level systems and employs
single-spin control and time-periodic modulations of the nearest-neighbour
hopping. While such spin lattices are routinely realized in various quantum
simulation platforms, the necessity to modulate spin-spin hopping makes
architectures of superconducting qubits with tuneable-coupling a natural
framework for implementing our proposal [35, 36, 37, 38, 39, 40]. Indeed, this
type of control has been demonstrated and used for the realization of
efficient two-qubit gates [39, 40] and artificial gauge fields for microwave
photons [37].

The presentation is organized as follows. In Section II, the driven model and
the target effective Hamiltonian are introduced. In Section III.1, the driving
sequence which Floquet engineers the characteristic four-spin interactions of
the toric code is presented. The protocol to generate the target Hamiltonian
on the full lattice is then described in Sec. III.2. In Section IV, it is
shown how the toolbox developed can be used to prepare the toric code ground
state with high fidelity and, in Sec. V, signatures of topological order in
the state prepared are analysed. In Section VI, we propose a minimal
implementation of a proof-of-principle topological device.

##  II Driven and target system Hamiltonians

The Hamiltonian of the system describes two-level systems (spins) in a 2D
lattice, and has the form $\hat{H}(t)=\hat{H}_{1}(t)+H_{2}(t)$ with

$\hat{H}_{1}(t)=\sum_{\langle\alpha,\beta\rangle}g_{\alpha\beta}(t)(\hat{X}_{\alpha}\hat{X}_{\beta}+\hat{Y}_{\alpha}\hat{Y}_{\beta})\
,$ (1)

where the summation runs over nearest-neighbor pairs and
$\\{\hat{X}_{\alpha},\hat{Y}_{\alpha},\hat{Z}_{\alpha}\\}$ denote Pauli
matrices related to the $\alpha$th spin. The indices $\alpha$ and $\beta$
indicate pairs of coordinates $(i,j)$ in the 2D lattice. The geometry and
connectivity of the lattice is sketched in Fig. 1, and will be specified more
in detail in the following (see Fig. 4). The inter-qubit hopping will be
periodically modulated in time, $g_{\alpha\beta}(t)=g_{\alpha\beta}(t+T)$. The
Hamiltonian $\hat{H}_{2}(t)$ describes additional terms corresponding to
resonant single-qubit pulses, of the form
$\Omega_{\alpha}(t)\hat{\sigma}_{\alpha}$ for
$\hat{\sigma}_{\alpha}=\hat{X}_{\alpha},\hat{Y}_{\alpha}$ or
$\hat{Z}_{\alpha}$. The Hamiltonian (1) is written in the interaction picture
with respect to the qubit energies (see Appendix A) and can describe an array
of superconducting qubits with controllable coupling, as has been realized
with different strategies and architectures [35, 36, 37, 38, 39, 40].
Modulation of $g_{\alpha\beta}(t)$ is achieved, for example, by means of an
intermediate coupler and has been used to generated artificial magnetic fields
[37].

The target of the quantum simulation is to realize dynamics given by Wen’s
plaquette model [41] on a square lattice,

$\displaystyle\hat{H}_{w}=-\mathcal{J}\sum_{i,j}\hat{P}_{i,j}\ ,$
$\displaystyle\hat{P}_{i,j}=\hat{X}_{i,j}\hat{Z}_{i,j+1}\hat{Z}_{i+1,j}\hat{X}_{i+1,j+1}\
.$ (2)

Wen’s model directly maps to the more commonly studied Kitaev’s toric-code
Hamiltonian [10] via a local basis rotation for a subset of qubits. Indeed,
Kitaev’s model is recovered by transforming Eq. (2) through Hadamard gates
$(\hat{X}_{i,j}+\hat{Z}_{i,j})/\sqrt{2}$ on all sites with $i+j$ either even
or odd [42, 43]. Besides, Wen’s formulation has recently found renewed
attention as a promising surface code design against biased noise [44, 45].
Our goal is then to find control functions $g_{\alpha\beta}(t)$ and
$\Omega_{\alpha}$ such that the dynamics produced by the Hamiltonian (1), at
times matching multiples of a fundamental Floquet period $T$, reproduces as
accurately as possible the dynamics generated by $\hat{H}_{w}$ of Eq. (2).

##  III Floquet engineering scheme

###  III.1 Single-plaquette sequence and four-spin interactions

Figure 2: (a) Sketch of the driving scheme for realizing a four-spin term in a
single four-spin plaquette, involving oscillating fields at $\omega$ and
$2\omega$; (b) Components (larger than $10^{-7}\omega$ in magnitude) of the
effective plaquette Hamiltonian, confirming the achievement of a clean four-
spin interaction, in units of $\omega$.

The central difficulty in achieving the Hamiltonian $\hat{H}_{w}$ is the
realization of the four-spin interactions $\hat{P}_{i,j}$ of Eq. (2), and the
effective cancellation of natural single- and two-spin terms. This problem
will now be addressed at the level of a fundamental four-spin plaquette.
Considering a time step $\tau$ (which will be a submultiple of the overall
Floquet period $T$), we analyze the effective Hamiltonian $\hat{H}_{i,j}$
generated at time $\tau$ by the driven Hamiltonian $\hat{H}(t)$ restricted to
a single plaquette. The effective Hamiltonian is here defined by ($\hbar=1$)

$e^{-i\hat{H}_{i,j}\tau}=\mathcal{T}\exp\left(-i\int_{0}^{\tau}dt\hat{H}(t)\right)\
.$ (3)

We then search for $\tau$-periodic control functions $g_{\alpha\beta}(t)$
yielding $\hat{H}_{i,j}\simeq-\mathcal{J}\hat{P}_{i,j}$ for some coupling
parameter $\mathcal{J}$. In Sec. III.2, it will then be discussed how the
resulting scheme can be used to engineer $\hat{H}_{w}$ on the full lattice,
for an arbitrary system size. The starting point is a plaquette with coupling
connectivity and spin labelling $k=1,2,3,4$ chosen as in Fig. 2(a). The driven
Hamiltonian of Eq. (1), in the absence of single-spin drives
($\hat{H}_{2}=0$), conserves the total magnetization $\sum_{k}\hat{Z}_{k}$ in
the whole plaquette, while the target four-body term conserves only the total
number of excitations residing on sites 2 and 3 instead. This excitation
conservation is broken for $\hat{H}(t)$ by introducing a strong resonant
$\hat{X}$-drive on spins $1$ and $4$,
$\hat{H}_{2}=(\Omega_{1}\hat{X}_{1}+\Omega_{4}\hat{X}_{4})$ with amplitudes
much larger than that of the two-spin coupling $g_{kk^{\prime}}$(t). The
intuition is that this drive energetically penalizes processes induced by the
$\hat{Y}\hat{Y}$, as compared to $\hat{X}\hat{X}$, interaction. These extra
drive terms would not be needed if one is able to implement, separately,
tuneable couplings of either $\hat{X}\hat{X}$ or $\hat{Y}\hat{Y}$ type. The
Hamiltonian becomes

$\hat{H}(t)=\Omega_{1}\hat{X}_{1}+\Omega_{4}\hat{X}_{4}+\sum_{\langle
k,k^{\prime}\rangle}g_{kk^{\prime}}(t)(\hat{X}_{k}\hat{X}_{k^{\prime}}+\hat{Y}_{k}\hat{Y}_{k^{\prime}})\
.$ (4)

In order to inspect the range of operators that can possibly compose the
effective Hamiltonian generated by $\hat{H}(t)$, it is useful to study the
dynamical Lie algebra $\mathfrak{L}$ of the control system [46]. This is
defined as the space spanned by all possible nested commutators of the control
operators $\hat{X}_{1}$, $\hat{X}_{4}$,
$\hat{X}_{k}\hat{X}_{k^{\prime}}+\hat{Y}_{k}\hat{Y}_{k^{\prime}}$ entering
$\hat{H}(t)$. It characterizes entirely the set of evolution operators that
can be produced by the driven dynamics and the set of reachable states [46].
Hence, the effective Hamiltonian $\hat{H}_{i,j}$ for a plaquette can be
generically expressed in the form

$\hat{H}_{i,j}=\sum_{\hat{O}_{\ell}\in\mathfrak{L}}c_{\ell}\hat{O}_{\ell}\ ,$
(5)

for some (unknown) coefficients $c_{\ell}$. A linearly independent set of
operators $\hat{O}_{\ell}$ spanning $\mathfrak{L}$ is given in Table 1. For
each term, the number of commutators of the control operators needed to
produce it is also reported.

Table 1: Basis of operators spanning the dynamical Lie Algebra of the single-plaquette system. The left column indicates the minimal number of nested commutators of initial control operators needed to produce the corresponding terms in the right column.  # commutators | Operators
---|---
0 |  $\hat{X}_{1}$, $\hat{X}_{4}$, $\hat{X}_{1}\hat{X}_{2}+\hat{Y}_{1}\hat{Y}_{2}$,
|
$\hat{X}_{2}\hat{X}_{3}+\hat{Y}_{2}\hat{Y}_{3},\,\hat{X}_{3}\hat{X}_{4}+\hat{Y}_{3}\hat{Y}_{4}$
1 |  $\hat{Z}_{1}\hat{Y}_{2}$, $\hat{Y}_{3}\hat{Z}_{4}$, $\hat{X}_{1}\hat{Z}_{2}\hat{Y}_{3}-\hat{Y}_{1}\hat{Z}_{2}\hat{X}_{3}$,
| $\hat{X}_{2}\hat{Z}_{3}\hat{Y}_{4}-\hat{Y}_{2}\hat{Z}_{3}\hat{X}_{4}$
2 |  $\hat{Y}_{1}\hat{Y}_{2}$, $\hat{Y}_{3}\hat{Y}_{4}$, $\hat{Z}_{1}\hat{Z}_{2}\hat{X}_{3}$, $\hat{X}_{2}\hat{Z}_{3}\hat{Z}_{4}$,
|  $\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}$,
$\hat{Y}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Y}_{4}$,
3 |  $\hat{Z}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Y}_{4}$, $\hat{Y}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Z}_{4}$
4 | $\hat{Z}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Z}_{4}$

Building on intuition derived from high-frequency expansions for Floquet
systems and explained in Appendix B, the driving functions are chosen as
follows [and sketched in Fig. 2(a)],

$\displaystyle g_{13}(t)=g_{13}^{(0)}+g_{13}^{(1)}\cos(\omega t)\ ,$
$\displaystyle g_{23}(t)=g_{23}^{(0)}+g_{23}^{(1)}\cos(2\omega t)\ ,$
$\displaystyle g_{24}(t)=g_{24}^{(0)}+g_{24}^{(1)}\cos(2\omega t)\ ,$ (6)

with angular frequency $\omega=2\pi/\tau$. The choice of functions which are
time symmetric within one evolution period $\tau$,
$g_{kk^{\prime}}(t)=g_{kk^{\prime}}(\tau-t)$, together with the particular
structure of the dynamical Lie algebra, allows one to exclude from the
effective Hamiltonian some of the terms belonging to $\mathfrak{L}$, as
detailed in Appendix B.3. The remaining operators are linear combinations of
those corresponding to an even number of commutators in Table 1.

Although the high-frequency regime and related high-frequency expansions are
helpful for finding a promising ansatz for the functions $g_{kk^{\prime}}(t)$,
their perturbative character would lead to a very weak four-spin interaction,
and still sizeable higher-order unwanted terms. To circumvent this limitation,
one can take advantage of the fact that the single-plaquette dynamics can be
accessed in a numerically exact way [28], thanks to the small system size. In
particular, the driving amplitudes $g_{kk^{\prime}}^{(\ell)}$ are optimized
together with $\Omega_{j}$ to achieve a four-spin interaction
$\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}$ with a desired strength
$\mathcal{J}$. Details on the optimization procedure are given in Appendix C.
Different optimizations can be performed to have a set of different values of
$\mathcal{J}$ available. Since the Hamiltonian does not depend on $t$ and
$\omega$ separately, but only on their product, in the following, $\omega$ is
chosen as the unit of energy.

The magnitude of different terms in the effective Hamiltonian of a single
plaquette resulting from this procedure is shown in Fig. 2(b), for a set of
optimal parameters yielding a four-body term of chosen strength
$\mathcal{J}\tau=\pi/8$. Approximated values of these parameters are given in
Table 2. With this choice of $\mathcal{J}$, evolution for two Floquet periods
produces an entangling four-spin gate
$e^{i2\mathcal{J}\hat{P}_{i,j}\tau}=(\mathbb{1}+i\hat{P}_{i,j})/\sqrt{2}$,
which will be employed in Sec. IV for ground state preparation.

Table 2: Example sets of (rounded) optimal driving parameters in units of $\omega$. Row A corresponds to a four-spin interaction strength $\mathcal{J}\tau=\pi/8$, wheres row B corresponds to $\mathcal{J}\tau=\pi/50$. In row C, the parameters correspond to the three-spin interaction protocol described in Sec. VI.1 with $\mathcal{J}_{xzx}\tau=-\pi/50$. For parameters labeled by the index $j$, the values in parenthesis correspond to increasing values of $j$. | $\Omega_{j}$ | $g_{12}^{(j)}$ | $g_{23}^{(j)}$ | $g_{24}^{(j)}$
---|---|---|---|---
A | [12, 10] | [3.6$\cdot 10^{-4}$, 1.3] | [0.0034, 2.6] | [$-$9.4$\cdot 10^{-5}$, 0.4]
B | [14, 11] | [1.2, 2.6] | [5.4 , $0.0021$] | [0.0003, 0.0012]
| $\overline{\Omega}_{j}$ | [$g_{12}$, $g_{23}$] | $z_{j}$
C | [8.5, 9.5] | [0.3, 0.3] | [$-9.3\cdot 10^{-3}$, 0.8, 2.4]

As can be appreciated from Fig. 2(b), the four-body interaction is
successfully achieved and it is by far the dominant term in the effective
Hamiltonian, with other terms almost three orders of magnitude smaller. The
effective Hamiltonian for a plaquette can thus be written in the form

$\hat{H}_{i,j}=-\mathcal{J}\hat{P}_{i,j}+\varepsilon\hat{V}_{i,j}\ ,$ (7)

with a quantum simulation error $\varepsilon/\mathcal{J}\simeq 1.5\cdot
10^{-3}$ and with $\hat{V}_{i,j}$ containing the extra terms depicted in Fig.
2 with dimensionless prefactors $c_{\ell}/\varepsilon<1$.

Figure 3: Sensitivity of the single plaquette protocol to errors in the
control amplitudes. Different colors correspond to different maximal error
amplitude $\eta_{\mathrm{max}}$. The bars represent the operator with maximal
magnitude at a given weight (number of non-identity operators in the tensor
product), which is indicated on the top of the bar.

In order to benchmark the robustness of the effective four-spin Hamiltonian
against imperfections in the control parameters, the leading effective terms
in the presence of errors are explored in Fig. 3. In particular, errors are
introduced in all parameters $p$ in the first row of Table 2 according to

$p\longrightarrow p+\eta,$ (8)

where $\eta$ represents uniformly distributed numbers in the interval
$[-\eta_{\mathrm{max}},\eta_{\mathrm{max}}]$. The effective Hamiltonian is
then averaged over 1000 realizations of errors affecting all parameters
simultaneously. The resulting averaged effective Hamiltonian can be decomposed
in the form (5) with strings of Pauli operators $\hat{O}_{\ell}$ and
coefficients $c_{\ell}$. Figure 3 depicts the maximal value of the magnitudes
$|c_{\ell}|$ of different terms for a given weight (number of non-identity
terms in the Pauli string). One can appreciate that the effective four-body
interaction persists even in the presence of relatively strong imperfections.
For example, from Fig. 3 one can see that errors larger than $0.1\omega$ must
occur for the four-body term to be masked by other effective terms.

###  III.2 Full lattice sequence

The continuous-driving protocol for the single plaquette presented in Sec.
III.1 represents the fundamental building block that is used, in the
following, to construct the full Hamiltonian on the whole lattice, by means of
a Suzuki-Trotter sequence [47, 48]. Although the present scheme can be applied
to generic lattice topologies, for experimental practicality we will focus on
planar lattices. The connectivity, scaling up the single-plaquette one, is
shown in Fig. 4(a), and is equivalent to a (distorted) square-lattice
geometry. Some additional links [dashed in Fig. 4(a)], making up triangular
cells, are needed for implementing the mixed boundary terms discussed in the
following.

Figure 4: (a) Lattice connectivity, equivalent to a distorted square lattice.
Some triangular cells appear when including links (dashed lines) needed to
implement mixed boundaries. (b) Trotter sequence to realize the Hamiltonian on
the whole lattice. In each of the four steps of duration $\tau$, a group
$G_{k}$ of disconnected plaquettes is driven (blue colour), while couplings
connecting such plaquettes are turned off. The full lattice Hamiltonian
progressively builds up (light blue).

A hybrid “Floquet-Trotter” quantum simulation strategy is particularly suited
for the problem treated in this work and can attain a dramatic improvement as
compared to straightforward digitization or high-frequency perturbative
Floquet engineering. This is related to the interplay of two ingredients,
namely, the fact that (i) at the single-plaquette level, continuous drives can
be numerically optimized yielding very clean four-spin interactions within one
time substep $\tau$, and that (ii) such terms approximately commute on
different sites, due to the properties of the toric code, thus yielding very
small digitization errors when considering multiple plaquettes [28]. Hence,
once the driving parameters for one plaquette are optimized, achieving the
fundamental four-spin interaction with a small error $\varepsilon$ [see Eq.
(7)], scaling to arbitrary system size is immediate without dramatic error
amplification and without needing further parameter optimization.

A possible first-order Suzuki-Trotter sequence features four steps and is
depicted in Fig. 4(b). In each step, corresponding to one driving substep
$\tau$, the coupling along lattice links connecting different blue-colored
plaquettes are turned off by setting $g_{\alpha\beta}(t)=0$, while the single-
plaquette protocol depicted in Fig. 2(a) is applied to each dark-blue coloured
plaquettes, engineering the desired four-body terms. These terms add up
progressively in the sequence, building up the Hamiltonian of Eq. (2) at time
$T=4\tau$. The propagator resulting from the Trotter sequence can be written
in the form

$\hat{U}(T)=\prod_{k=1}^{4}\prod_{(i,j)\in{G}_{k}}\hat{U}_{ij}(\tau)\ ,$ (9)

where $\hat{U}_{ij}(\tau)=e^{-i\tau\hat{H}_{i,j}}$ represent single-plaquette
propagators, while ${G}_{k}$ are the four groups of plaquettes involved in the
sequence and depicted in Fig. 4(b). Since each group $G_{k}$ contains
disconnected plaquettes, the corresponding effective Hamiltonian is simply the
sum of the effective Hamiltonians of the individual plaquettes involved,

$\hat{H}_{k}=\sum_{(i,j)\in G_{k}}\hat{H}_{i,j}\ .$ (10)

The product of propagators $\hat{U}_{k}=e^{-i\hat{H}_{k}\tau}$ belonging to
two different groups do not commute instead and can be estimated by means of
the exponential product formula,

$\prod_{k}\hat{U}_{k}=\exp\left({-i\tau\sum_{k}\hat{H}_{k}+\frac{\tau^{2}}{2}\sum_{j<k}[\hat{H}_{j},\hat{H}_{k}]+\dots}\right).$
(11)

Using the fact that the desired terms $\hat{P}_{i,j}$ contained in
$\hat{H}_{k}$ via the $\hat{H}_{i,j}$ commute with each other, the unwanted
sub-leading-order terms in the exponential in Eq. (11) are of the order of the
error terms $\propto\epsilon$ in Eq. (7). The effective Hamiltonian at the end
of the sequence can thus be estimated as

$\hat{H}_{\mathrm{eff}}=\hat{H}_{w}+\widehat{\mathrm{Err}},$ (12)

where $\widehat{\mathrm{Err}}$ is at most of order $\varepsilon$,

$\displaystyle\widehat{\mathrm{Err}}=\varepsilon\sum_{i,j}\hat{V}_{i,j}-$ (13)
$\displaystyle\frac{i\varepsilon\mathcal{J}\tau}{2}\sum_{k<k^{\prime}}\sum_{\begin{subarray}{c}(i,j)\in
G_{k}\\\ (i^{\prime},j^{\prime})\in
G_{k^{\prime}}\end{subarray}}\left([\hat{P}_{i,j},\hat{V}_{i^{\prime},j^{\prime}}]+[\hat{V}_{i,j},\hat{P}_{i^{\prime},j^{\prime}}]\right)+\dots\
.$

Thanks to the precision in the single-plaquette effective Hamiltonian achieved
through numerical optimization, the Trotter error is very small, independently
of the system size. Thus, here we are using the integrability of the target
Hamiltonian in order to reduce errors in our hybrid Floquet-Trotter approach.

Figure 5: Weight vs magnitude $|c_{\ell}|$ of the operators composing the
effective Hamiltonian for a 5-by-5 system. The desired four-spin and two-spin
boundary terms (exemplified in the inset) are at least two orders of magnitude
larger than undesired terms. The red colour of markers has non-unit opacity,
such that more colour-intense markers appear as a result of overlapping
markers.

To quantitatively inspect the magnitude of undesired terms arising from the
combination of Floquet engineering and Trotterization in a larger system, we
consider a 5-by-5 lattice with mixed boundary conditions as depicted in the
inset of Fig. 5. Such a geometry is interesting, since it features two
degenerate, topologically ordered ground states that can be used to define a
topologically protected qubit [23]. The boundary terms, involving
$\hat{X}_{\alpha}\hat{Z}_{\beta}$ terms, are Floquet engineered with the
driving scheme

$g_{\alpha\beta}\cos(\omega
t)(\hat{X}_{\alpha}\hat{X}_{\beta}+\hat{Y}_{\alpha}\hat{Y}_{\beta})+\lambda\sin(\omega
t)\hat{Y}_{\beta}\ ,$ (14)

where $g_{\alpha\beta}$ and $\lambda$ are numerically optimized to give the
same coupling $\mathcal{J}$ as that engineered for the plaquette interactions
in a time substep $\tau$. The Floquet-Trotter realization of the Hamiltonian
also allows us to develop a practical route to numerically compute the
effective Hamiltonian for such a large system. Once the effective Hamiltonian
$\hat{H}_{i,j}$ for a single plaquette (or boundary term) is determined
numerically by diagonalizing the corresponding end-of-period propagator
$\hat{U}_{ij}(\tau)=e^{-i\hat{H}_{i,j}\tau}$, we compute the full lattice
Hamiltonian via Eq. (9) and iteratively using Baker-Campbell-Hausdorf formula
for estimating the Hamiltonian, where terms in the series are included until
convergence is reached. To avoid evaluating the commutators via explicit
matrix multiplications, we develop a semi-analytic representation of the
Hamiltonian (Appendix D), which allows one to compute the commutators
analytically instead.

Considering a decomposition of the Hamiltonian into Pauli strings with
associated coefficients $c_{\ell}$, the magnitude $|c_{\ell}|$ of the
different terms in the effective Hamiltonian is shown in Fig. 5, where they
are organized according to their weight (namely, the number of non-identity
operators in the tensor product). The desired terms (four-body and boundary
two-body) are almost three orders of magnitude larger than unwanted terms. The
leading error terms appear at order $|c_{\ell}|/\omega\sim 10^{-4}$ and
correspond to two-spin $\hat{X}\hat{X}$ and $\hat{Y}\hat{Y}$ terms contained
already in $\hat{V}_{i,j}$ [first line of Eq. (13)]. The renowned topological
properties of the toric code are in general jeopardized when strings of Pauli
operators spanning the whole system and connecting different boundaries are
present [10]. From Fig. 5(b), one can further see that terms of such a weight
are very weak, as expected since they can be created only via deeply nested
commutators of single-plaquette effective Hamiltonians $\hat{H}_{i,j}$. For
this reason, their magnitude will also decay with the system size, since more
and more commutators will be needed to create a similar string. The
preservation of the toric-code topological properties for the quantum
simulated model will be inspected more in detail Sec. V.

Before moving on, we would like to stress that our hybrid continuous-digital
approach to Floquet engineering allows us to realize four-body terms, which in
the regime where standard high-frequencies expansions work, would correspond
to very small third-order contributions.

##  IV Ground state preparation

By exploiting the possibility to realize the plaquette operators
$\hat{P}_{i,j}$ with high precision and the related optimal driving parameters
found, an efficient ground-state preparation protocol can be formulated, which
can be applied to systems of arbitrary size. This is optimal in the sense that
the time needed to prepare the target state scales linearly with $N$ for an
$N\times N$ lattice, thus saturating Lieb-Robinson bounds on the preparation
of topological ordered states via local Hamiltonians [49]. The protocol
features a sequence in which different plaquette operators $\hat{P}_{i,j}$ are
Floquet engineered in a way such that evolution for two time steps $\tau$
implements a four-spin entangling gate, with the help of additional single-
qubit $\hat{Z}_{i,j}$ rotations. A ground state $\ket{G}$ of $\hat{H}_{w}$
satisfies the “stabilizer” constraint $\langle\hat{P}_{i,j}\rangle=1$ on all
plaquettes operators, and can be written as an equal superposition of closed
strings [10],

$\ket{G}=\prod_{i+j\,\mathrm{even}}\frac{(1+\hat{P}_{i,j})}{\sqrt{2}}\ket{\psi_{0}},$
(15)

where $\ket{\psi_{0}}=\mathcal{H}\ket{0}^{\otimes N}$, with $\mathcal{H}$ the
simultaneous Hadamard gate $(\hat{X}_{i,j}+\hat{Z}_{i,j})/\sqrt{2}$ on odd
sites $(i,j)$, namely such that $i+j$ is odd [the bottom-left spin has
coordinate (1,1)]. The sum in Eq. (15) involves even plaquettes. Even (odd)
plaquettes are defined as those whose bottom-left site is an even (odd) site.
Note that the ordering of factors in the product does not matter since all
factors commute. The ground state obtained has the same parity of
noncontractible loops along the $x$ or $y$ direction as $\ket{\psi_{0}}$. In
this case, for example,
$\bra{G}\hat{X}_{L}\ket{G}=\bra{\psi_{0}}\hat{X}_{L}\ket{\psi_{0}}=1$, where
$\hat{X}_{L}$ is a loop of alternating $\hat{X}_{i,j}$ and $\hat{Z}_{i,j}$
operators spanning the whole system along the $y$ direction, which commutes
with the Hamiltonian $\hat{H}_{w}$.

The individual operators $(1+\hat{P}_{i,j})$ are not unitary, and thus cannot
be achieved by a Hamiltonian evolution. Nonetheless, an equivalent result is
obtained, starting from state $\ket{\psi_{0}}$, by employing a suitably
ordered product of operators

$\displaystyle\hat{\mathcal{U}}_{i,j}^{A}=(1-i\hat{A}_{i,j})/\sqrt{2}=\exp(-i\pi\hat{A}_{i,j}/4)\
,$
$\displaystyle\hat{\mathcal{U}}_{i,j}^{B}=(1-i\hat{B}_{i,j})/\sqrt{2}=\exp(-i\pi\hat{B}_{i,j}/4)\
,$ (16)

where

$\displaystyle\hat{A}_{i,j}=\hat{X}_{i,j}\hat{Z}_{i,j+1}\hat{Z}_{i+1,j}\hat{Y}_{i+1,j+1}\
,$
$\displaystyle\hat{B}_{i,j}=\hat{Y}_{i,j}\hat{Z}_{i,j+1}\hat{Z}_{i+1,j}\hat{X}_{i+1,j+1}\
.$ (17)

The construction of this protocol is discussed in Appendix E. The ground state
$\ket{G}$ is obtained, for example, by applying ordered products of
$\hat{A}_{i,j}$ (or $\hat{B}_{i,j}$) operators along bottom-left to top-right
diagonals of the lattice involving even plaquettes, which are ordered from
bottom-left to top-right for $\hat{A}_{i,j}$, and vice versa for
$\hat{B}_{i,j}$. The choice between using operators $\hat{A}_{i,j}$ or
$\hat{B}_{i,j}$ for a given diagonal depends on the boundary: a boundary term
must always contain one $\hat{Y}$ operator, and hence the right boundary (and
thus the whole corresponding diagonal) necessitates operators $\hat{B}_{i,j}$,
whereas the left boundary needs $\hat{A}_{i,j}$.

Figure 6: Ground state preparation sequence for a 5-by-4 lattice. Dark and
light purple tiles indicate an operation $\hat{\mathcal{U}}_{i,j}^{A}$ and
$\hat{\mathcal{U}}_{i,j}^{B}$, respectively, each involving two time steps
$\tau$ and two additional single-qubit rotations according to Eq. (18).

Exploiting the Floquet-engineered plaquette operators, the gates
$\hat{\mathcal{U}}_{i,j}^{A}$ and $\hat{\mathcal{U}}_{i,j}^{B}$ can be
realized by combining the single-plaquette Floquet protocol with two single-
qubit rotations. Indeed, introducing the notation
$\hat{\mathcal{U}}_{i,j}^{P}=e^{i\frac{\pi}{4}\hat{P}_{i,j}}$ and
$\hat{\mathcal{Z}}_{i,j}^{\pi/4}=e^{-i\frac{\pi}{4}\hat{Z}_{i,j}}$, we find

$\displaystyle\hat{\mathcal{U}}_{i,j}^{A}=\hat{\mathcal{Z}}_{i+1,j+1}^{\pi/4,\dagger}\hat{\mathcal{U}}_{i,j}^{P}\hat{\mathcal{Z}}_{i+1,j+1}^{\pi/4}\
,$
$\displaystyle\hat{\mathcal{U}}_{i,j}^{B}=\hat{\mathcal{Z}}_{i,j}^{\pi/4,\dagger}\hat{\mathcal{U}}_{i,j}^{P}\hat{\mathcal{Z}}_{i,j}^{\pi/4}\
.$ (18)

To obtain the $\pi/4$ rotation needed, the driving parameters have been
optimized to obtain $|\mathcal{J}|\tau=\pi/8$, such that evolution for two
steps $2\tau$ implements the desired gate with $\pi/4$ angle. The single-qubit
gates will typically be much faster than the time step $\tau$, such that one
can approximately consider them to be instantaneous on the timescale of
$\tau$.

The sequence for a $5\times 4$ lattice is shown in Fig. 6. The state
$\ket{G_{\mathrm{eff}}}$ approximating $\ket{G}$ prepared using the Floquet
gates for such a system has an energy matching the exact toric-code ground-
state energy up to an error $\sim 10^{-5}\omega$. Since different diagonals
can be addressed in parallel, the overall number of four-spin gates needed
scales with the number of plaquettes tiling the longest diagonal of the
lattice. This, in turn, scales linearly in the system size $N$, for a $N\times
N$ system. Therefore, the ground state preparation protocol attains “quantum
speed limits” on the preparation of topologically ordered states [49], being
in this sense optimal.

##  V Probing entanglement and ayons

The robustness of topological order in the toric code ground state against
perturbations is the paradigmatic example of topological protection [10, 50,
15, 51, 52, 53]. Weak local perturbations can in general only slightly deform
the ideal eigenstates, without spoiling their topological properties.
Moreover, the ground state and anyonic excitations can still be efficiently
manipulated by means of unperturbed string operators [51]. The goal of this
Sec. is to show that, despite the fact that the quantum simulation procedure
also produces non-local higher-order contributions, approximate ground states
prepared with the procedure described in Sec. IV indeed feature and enable the
detection of topological order. To this end, typical signatures of topological
order are probed in a 5-by-4 lattice, namely the presence of long-range
entanglement and the existence of quasiparticles with non-trivial mutual and
exchange statistics.

Figure 7: (a) Topological entanglement entropy for two different partitions of
the system, where the subsets A, B, C are depicted in orange, purple and blue,
respectively. (b) von Neumann entropy of the reduced state including the spins
in red, which satisfies $(n-1)\log(2)$ where $n$ is the number of even
plaquettes crossed (indicated in dark blue). (c) Specific anyon creation and
braiding procedures and (d) dyon exchange, studied starting from the effective
ground state $\ket{G_{\mathrm{eff}}}$ prepared through the Floquet scheme of
Sec. IV.

To inspect the presence of entanglement as predicted in the toric code, the
topological entanglement entropy $S_{\mathrm{topo}}$ [54, 55, 56] is reported
in Fig. 7(a) for two different partitions of the system. It is defined as [54,
23]

$S_{\mathrm{topo}}=S_{A}+S_{B}+S_{C}-S_{AB}-S_{BC}-S_{AC}+S_{ABC},$ (19)

where $S_{X}=\mathrm{tr}\rho_{X}\log\rho_{X}$ is the von Neumann entropy of
the system state reduced to subsystem $X$. The subsystems $A$, $B$ and $C$ are
depicted in orange, purple and blue color, respectively, in Fig. 7(a). As
indicated in the figure, $S_{\mathrm{topo}}$ is very close to the value of the
ideal toric code, $-\log 2$, which characterizes $\mathbb{Z}_{2}$ topological
order, and differs only by $\sim 10^{-4}$ in all cases. Additionally, it is
verified that, for a single subsystem depicted in Fig. 7(b), the von Neumann
entropy scales like $(n-1)\log(2)$, where $n$ is the number of even plaquettes
crossed.

As a second signature, the possibility to create and manipulate anyons in the
same 5-by-4 system is probed, verifying mutual and exchange statistics through
explicit braiding and exchange. The toric code possesses three types of gapped
quasiparticle excitations, that appear as violations of the stabilizer
constraint ($\langle\hat{P}_{i,j}\rangle=-1$) occurring at the end of open
strings of Pauli operators. These are so-called electric charges
($`\mathrm{e}$’) and magnetic vortices (‘$\mathrm{m}$’), associated to a
violation of even or odd plaquettes, respectively (the distinction is
arbitrary [57]), and dyons (‘$\epsilon$’, a compound
$\epsilon=\mathrm{e}\times\mathrm{m}$ quasiparticle). The former are bosons,
whilst the latter are fermions. However, if an ‘$\mathrm{e}$’ charge is
braided around a ‘$\mathrm{m}$’ charge or vice versa, the wavefunction
acquires a $\pi$ phase, and thus these particles exhibit anyonic (semionic)
mutual statistics. In the presence of weak perturbations, these quasiparticles
are not exact eigenstates of the system, and will propagate with an effective
finite mass [10]. Nonetheless, their characteristic properties are predicted
to persist, until the perturbation becomes strong enough to disrupt the
topological phase [10, 51, 50, 52, 16].

We verify these features through the processes depicted in Fig. 7, where the
creation and braiding of an electric charge around a magnetic vortex [7(c)],
and the exchange of two dyons $\epsilon$ [7(d)], are depicted. These processes
are realized in the system by applying strings of single-spin operators (shown
in Fig. 7) on top of the ground state $\ket{G_{\mathrm{eff}}}$ prepared with
the protocol of Sec. IV. In the exact toric code, both operations result in
the wavefunction acquiring a phase $\pi$ due to the semionic nature of
‘$\mathrm{e}$’ and ‘$\mathrm{m}$’, in the first case, and to the fermionic
nature of ‘$\epsilon$’, in the second. This is obtained almost exactly in the
system studied here, where the phases obtained are $\pi$ up to numerical
precision (resulting in a sign change $s$ deviating from $-1$ by $\sim
10^{-4}$). The bosonic exchange statistics of ‘$\mathrm{e}$’ and
‘$\mathrm{m}$’ was also verified (not shown).

These results give strong evidence that the state prepared is topologically
ordered and features the defining properties of the $\mathbb{Z}_{2}$
topological phase very accurately. This further confirms that higher order
contributions to the effective Hamiltonian have very little impact on the
system properties, providing an additional signature of successful Floquet
engineering.

##  VI Proof-of-principle device

In the earlier Sections, tools were developed for quantum simulating a toric-
code Hamiltonian, preparing its ground state, and it has been verified that
the prepared state indeed exhibits clear signatures of topological order.
Based on this toolbox, we next propose a full-fledged proof-of-principle
realization, that can be used both to study transitions into the spin-liquid
phase and to implement a prototypical topological qubit. For this purpose, the
nine-spin geometry shown in Fig. 8(a) is considered. The system features two
degenerate ground states, which define the logical qubit. Topological order
implies that these ground states cannot be distinguished via local
observables, nor be turned one into the other by weak local perturbations. A
string of single-spin errors needs to cross the entire system, involving at
least three spins, to produce a logical error. The connectivity considered is
shown in Fig. 8(a), and the Floquet-Trotter sequence realizing the
Hamiltonian, developed following the prescription of Sec. III.2, is shown in
Fig. 8(b). The two ground states are distinguished by logical operators
$\hat{X}_{L}$ and $\hat{Z}_{L}$ which span the system vertically and
horizontally, commute with the Hamiltonian and anticommute with each other. A
possible choice, which will be adopted in the following, is
$\hat{X}_{L}=\hat{X}_{1,2}\hat{Z}_{2,2}\hat{X}_{3,2}$ and
$\hat{Z}_{L}=\hat{Z}_{2,1}\hat{X}_{2,2}\hat{Z}_{2,3}$. Logical states
$\ket{0_{L}}$ and $\ket{1_{L}}$ can then be defined according to

$\hat{Z}_{L}\ket{0_{L}}=+\ket{0_{L}},\quad\hat{Z}_{L}\ket{1_{L}}=-\ket{1_{L}}\
.$ (20)

For example, $\hat{X}_{L}$ is depicted in Fig. 8(a). For a chosen four-spin
coupling $\mathcal{J}=0.01\omega$, the quasienergy spectrum of the nine-qubit
device is shown in Fig. 8(c) [black colour], where it is also compared with
the ideal spectrum (red colour). The quasienergy spectrum reproduces the ideal
spectrum faithfully, featuring almost flat bands that acquire a slight
dispersion due the weak higher-order terms in the effective Hamiltonian.
Ground states can be prepared with the protocol of Sec. IV. By taking
advantage of the small size of the system, we also propose in the following
additional operations that allow one to manipulate the system, although their
scaling to larger system is non-obvious.

Figure 8: (a) Device connectivity, minimal-qubit geometry, and $\hat{X}_{L}$
logical operator. (b) Trotter sequence realizing $\hat{H}_{w}$. (c)
quasienergy spectrum of the effective Hamiltonian (black), compared to the
exact spectrum of the ideal Hamiltonian (red).

###  VI.1 Realizing a logical operator

For a large lattice, the implementation of logical operators $\hat{X}_{L}$ and
$\hat{Z}_{L}$ into the Hamiltonian is generally not feasible, since they are
long Pauli strings spanning all across the system. However, for the small
prototypical device analyzed in this Sec., the operator $\hat{X}_{L}$ spans
three sites “only” and can also be realized via Floquet engineering, by
adapting earlier protocols for achieving three-spin interactions [28]. In
particular, considering the three spins in the central column, labeled from 1
to 3, the protocol Hamiltonian involves periodic driving of the central site
and reads

$\displaystyle\hat{H}_{X_{L}}(t)=$
$\displaystyle\overline{\Omega}_{1}\hat{X}_{1}+\overline{\Omega}_{3}\hat{X}_{3}+g_{12}(\hat{X}_{1}\hat{X}_{2}+\hat{Y}_{1}\hat{Y}_{2})$
$\displaystyle+g_{23}(\hat{X}_{2}\hat{X}_{3}+\hat{Y}_{2}\hat{Y}_{3})+z(t)\hat{Z}_{2}\
.$ (21)

The $\tau$-periodic drive $z(t)=z(t+\tau)$ is chosen to feature multiple
harmonics,

$z(t)=z_{1}+z_{2}\cos(\omega t)+z_{3}\cos(2\omega t)\ .$ (22)

The parameters $\overline{\Omega}_{j}$, $g_{ij}$, $z_{j}$ are numerically
optimized simultaneously to obtain an effective Hamiltonian
$\mathcal{J}_{xzx}\hat{X}_{L}=\mathcal{J}_{xzx}\hat{X}_{1}\hat{Z}_{2}\hat{X}_{3}$,
with a chosen magnitude $\mathcal{J}_{xzx}$. A set of optimal parameters is
given in Table 2. The capability to implement a logical operator at the
Hamiltonian level is particularly useful since it allows one to distinguish,
and actively manipulate, protected qubit states. For example, it can be used
to induce Rabi oscillations between the states $\ket{0_{L}}$ and
$\ket{1_{L}}$. In the following, it will be exploited to adiabatically prepare
logical states $\ket{\pm_{L}}=(\ket{0_{L}}\pm\ket{1_{L}})/\sqrt{2}$,
satisfying $\hat{X}_{L}\ket{\pm_{L}}=\pm\ket{\pm_{L}}$.

###  VI.2 Adiabatic manipulation

While a system’s ground state can be efficiently prepared with the protocol of
Sec. IV, an alternative approach can also be used based on adiabatic parameter
variation at the level of the effective Hamiltonian [58, 32]. Although
adiabatic— and even optimized diabatic [59]— preparation in interacting
Floquet systems is challenging, due to the interplay of both nonadiabaticities
and Floquet heating [32], this approach is particularly interesting in view of
studying transitions into or from the spin liquid phase. We move first steps
in this direction by discussing the passage into a spin-liquid “puddle” in the
prototypical nine-spins device. In this system, this method is also
interesting since it provides a versatile qubit manipulation toolbox. The
protocol performs an adiabatic interpolation from an initial effective
Hamiltonian

$\hat{H}_{\mathrm{in}}=-R\sum_{\alpha}\hat{X}_{\alpha}+\hat{H}_{w}\pm\mathcal{J}\hat{X}_{L}\
,$ (23)

with $R\gg\mathcal{J}$, to a final effective Hamiltonian

$\hat{H}_{\mathrm{fin}}=\hat{H}_{w}\pm\mathcal{J}\hat{X}_{L}\ .$ (24)

At the initial time, the system is prepared with all spins aligned in the
trivial magnetized state $\ket{-}$ ($-1$ eigenstate of $\hat{X}_{\alpha}$),
which is approximately the ground state of $\hat{H}_{\mathrm{in}}$. If the
preparation is successful, at the final time the system will be in the ground
state of $\hat{H}_{\mathrm{fin}}$. This will be a logical state
$\ket{\pm_{L}}$, depending on the sign of the $\hat{X}_{L}$ term (which can be
controlled with the drive). The role of $\hat{X}_{L}$ is indeed to distinguish
between the otherwise degenerate ground states of $\hat{H}_{w}$. For both
cases, the final state will be a topologically ordered state, such that the
adiabatic preparation realizes a finite-size precursor of a topological phase
transition [50, 15, 60, 52].

Figure 9: (a) Quasienergy spectrum as a function of $R$. (b) Adiabatic state
preparation. The dark red curve refers to the left $y$ axis and represents the
squared fidelity with respect to the exact toric-code ground state
$\ket{-_{L}}$. The darkblue dashed curve refers to the right $y$ axis and
depicts the corresponding growth of the entanglement entropy for the subsystem
depicted in the inset.

The Floquet-Trotter sequence of Sec. III.2 implementing $\hat{H}_{w}$ is
“always-on” during the protocol, which has duration $t_{f}$. The control
parameters are optimized in order to achieve a weaker four-body coupling
strength $\mathcal{J}=0.01\omega$: in this way, the quasienergy spectrum of
$\hat{H}_{w}$ fits entirely within a single Floquet-Brillouin zone,
simplifying the application of the adiabatic principle to the driven system
[58, 32]. The fundamental Trotter step features one additional substep (thus,
$T=5\tau$), as compared to the sequence of Sec. III.2, which is used to
implement $\hat{X}_{L}$ according to the driving pattern of Sec. VI.1. The
$\hat{X}_{\alpha}$ terms in $\hat{H}_{\mathrm{in}}$ are obtained by modifying
the magnitude of single-spin terms in Eq. 4 when driving each plaquette during
the four-step sequence, $\Omega_{k}\to\Omega_{k}-R$. The interpolation takes
place in a total time $t_{f}$ by ramping down $R\to R(t)=Rr(t)$, with ramping
function $r(t)$ such that $r(0)=1$ and $r(t_{f})=1$. Note that, at
intermediate times, the effective Hamiltonian will not be simply a linear
combination of $\hat{H}_{\mathrm{in}}$ and $\hat{H}_{\mathrm{fin}}$ only, but
additional terms will appear due to the interplay of the ramped $\hat{X}_{j}$
terms with the Floquet-Trotter drive. The quasienergy spectrum as a function
of $R$ for the driven system is depicted in Fig. 9(a).

A successful Floquet-adiabatic state preparation is obtained in $200$ Trotter
steps using, for instance, $r(t)=\arctan(15t/t_{f})/\arctan(15)$, as shown in
Fig. 9(b). The latter Figure reports the instantaneous fidelity
$\lvert\braket{\psi(t)}{-_{L}}\\!\lvert^{2}$ between the evolving state
$\ket{\psi(t)}$ and the exact ground state $\ket{-_{L}}$ (red curve, left
ordinates). The corresponding growth in time of the topological entanglement
entropy is also reported (dashed blue, right ordinates), for the system
partition shown in the inset, and further signals the achievement of
topological order.

##  VII Outlook

The possibility to Floquet engineer a clean toric-code Hamiltonian as proposed
in this work, with the added flexibility in the control of coupling and
operators defining the effective Hamiltonian, provides a full-fledged quantum
simulation strategy of $\mathbb{Z}_{2}$ lattice gauge theory [16, 15] and
string-net condensation [41, 17, 6]. Moreover, the proposed Floquet scheme for
implementing Wen’s Hamiltonian (2) provides an interesting perspective for
investigating related quantum computational models. For example, it could be
used to study quantum computation based on topological defects. Topological
defects such as holes, which consist in the removal of a subset of stabilizer
operators $\hat{P}_{i,j}$ from the Hamiltonian, admit the presence of anyons
residing in the hole, which do not cost any energy (since the stabilizers in
the hole have been removed from the Hamiltonian). The possibility to create,
annihilate, and braid anyons, together with the holes potentially hosting
them, enables universal quantum computation [61]. While anyon manipulation
needs only single-spin operations, once the vacuum state is prepared, creating
and moving holes is typically a more challenging task [61, 7]. This is
particularly simple in the Floquet-Trotter scheme proposed here, where
stabilizers $\hat{P}_{i,j}$ can be removed by simply avoiding to drive the
corresponding plaquette in the sequence of Sec. III.2. Similarly, a hole can
be restored by simply reintroducing the corresponding drive in the subsequent
Trotter step. To displace an occupied hole, one can combine such operations
with single-spin operations to move the hosted anyon along with the hole.
These fundamental operations then constitute a promising toolbox for
implementing a prototypical Floquet quantum computation based on holes.

###### Acknowledgements.

F.P. and A.E. acknowledge support from the Deutsche Forschungsgemeinschaft
(DFG) via the Research Unit FOR 2414 under Project No. 277974659.

## References

  * Tsui  _et al._ [1982] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Two-Dimensional Magnetotransport in the Extreme Quantum Limit, Phys. Rev. Lett. 48, 1559 (1982).
  * Kalmeyer and Laughlin [1987] V. Kalmeyer and R. B. Laughlin, Equivalence of the resonating-valence-bond and fractional quantum Hall states, Phys. Rev. Lett. 59, 2095 (1987).
  * Wen [1989] X.-G. Wen, Vacuum degeneracy of chiral spin states in compactified space, Phys. Rev. B 40, 7387 (1989).
  * Wen [1990] X.-G. Wen, Topological Orders in Rigid States, Int. J. Mod. Phys. B 4, 239 (1990).
  * Wen [2017] X.-G. Wen, Colloquium: Zoo of quantum-topological phases of matter, Rev. Mod. Phys. 89, 041004 (2017).
  * Wen [2004] X.-G. Wen,  _Quantum Field Theory of Many-Body Systems_ (Oxford University Press (Oxford), 2004).
  * Fowler  _et al._ [2012] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Surface codes: Towards practical large-scale quantum computation, Phys. Rev. A 86, 032324 (2012).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-Abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Terhal [2015] B. M. Terhal, Quantum error correction for quantum memories, Rev. Mod. Phys. 87, 307 (2015).
  * Kitaev [2003] A. Kitaev, Fault-tolerant quantum computation by anyons, Ann. Phys. 303, 2 (2003).
  * [11] S. B. Bravyi and A. Y. Kitaev, Quantum codes on a lattice with boundary, arXiv:quant-ph/9811052 .
  * Dennis  _et al._ [2002] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, Topological quantum memory, J. Math. Phys. 43, 4452 (2002).
  * Wegner [1971] F. J. Wegner, Duality in Generalized Ising Models and Phase Transitions without Local Order Parameters, J. Math. Phys. 12, 2259 (1971).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Rev. Mod. Phys. 51, 659 (1979).
  * Hamma and Lidar [2008] A. Hamma and D. A. Lidar, Adiabatic Preparation of Topological Order, Phys. Rev. Lett. 100, 030502 (2008).
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Progr. Phys. 80, 016502 (2016).
  * Wen [2003a] X.-G. Wen, Quantum order from string-net condensations and the origin of light and massless fermions, Phys. Rev. D 68, 065003 (2003a).
  * Lu  _et al._ [2009] C.-Y. Lu, W.-B. Gao, O. Gühne, X.-Q. Zhou, Z.-B. Chen, and J.-W. Pan, Demonstrating Anyonic Fractional Statistics with a Six-Qubit Quantum Simulator, Phys. Rev. Lett. 102, 030502 (2009).
  * Luo  _et al._ [2018] Z. Luo, J. Li, Z. Li, L.-Y. Hung, Y. Wan, X. Peng, and J. Du, Experimentally probing topological order and its breakdown through modular matrices, Nature Phys. 14, 160 (2018).
  * 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, Demonstration of Topological Robustness of Anyonic Braiding Statistics with a Superconducting Quantum Circuit, Phys. Rev. Lett. 121, 030502 (2018).
  * Andersen  _et al._ [2020] C. K. Andersen, A. Remm, S. Lazar, S. Krinner, N. Lacroix, G. J. Norris, M. Gabureac, C. Eichler, and A. Wallraff, Repeated quantum error detection in a surface code, Nature Phys. 16, 875 (2020).
  * 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 (2021).
  * Satzinger  _et al._ [2021] K. J. Satzinger, Y.-J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, C. Gidney, I. Aleiner, F. Arute, K. Arya, J. Atalaya, and et al., Realizing topologically ordered states on a quantum processor, Science 374, 1237 (2021).
  * [24] Y.-J. Liu, K. Shtengel, A. Smith, and F. Pollmann, Methods for simulating string-net states and anyons on a digital quantum computer, arXiv:2110.02020 .
  * Dai  _et al._ [2017] H.-N. Dai, B. Yang, A. Reingruber, H. Sun, X.-F. Xu, Y.-A. Chen, Z.-S. Yuan, and J.-W. Pan, Four-body ring-exchange interactions and anyonic statistics within a minimal toric-code Hamiltonian, Nature Phys. 13, 1195 (2017).
  * Paredes and Bloch [2008] B. Paredes and I. Bloch, Minimum instances of topological matter in an optical plaquette, Phys. Rev. A 77, 023603 (2008).
  * Sameti  _et al._ [2017] M. Sameti, A. Potočnik, D. E. Browne, A. Wallraff, and M. J. Hartmann, Superconducting quantum simulator for topological order and the toric code, Phys. Rev. A 95, 042330 (2017).
  * Petiziol  _et al._ [2021] F. Petiziol, M. Sameti, S. Carretta, S. Wimberger, and F. Mintert, Quantum Simulation of Three-Body Interactions in Weakly Driven Quantum Systems, Phys. Rev. Lett. 126, 250504 (2021).
  * [29] N. H. Le, M. Cykiert, and E. Ginossar, Scalable and robust quantum computing on qubit arrays with fixed coupling, arxiv.2110.07737 .
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Periodically Driven Quantum Systems: Effective Hamiltonians and Engineered Gauge Fields, Phys. Rev. X 4, 031027 (2014).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Adv. Phys. 64, 139 (2015).
  * Eckardt [2017] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017).
  * Oka and Kitamura [2019] T. Oka and S. Kitamura, Floquet Engineering of Quantum Materials, Annu. Rev. Condens. Matter Phys. 10, 387 (2019).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).
  * Grajcar  _et al._ [2006] M. Grajcar, Y.-x. Liu, F. Nori, and A. M. Zagoskin, Switchable resonant coupling of flux qubits, Phys. Rev. B 74, 172505 (2006).
  * Chen  _et al._ [2014] Y. Chen, C. Neill, P. Roushan, N. Leung, M. Fang, R. Barends, J. Kelly, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, A. Megrant, J. Y. Mutus, P. J. J. O’Malley, C. M. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, M. R. Geller, A. N. Cleland, and J. M. Martinis, Qubit Architecture with High Coherence and Fast Tunable Coupling, Phys. Rev. Lett. 113, 220502 (2014).
  * Roushan  _et al._ [2017] P. Roushan, C. Neill, A. Megrant, Y. Chen, R. Babbush, R. Barends, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, P. J. O’Malley, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, E. Kapit, H. Neven, and J. Martinis, Chiral ground-state currents of interacting photons in a synthetic magnetic field, Nature Phys. 13, 146 (2017).
  * Yan  _et al._ [2018] F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Tunable Coupling Scheme for Implementing High-Fidelity Two-Qubit Gates, Phys. Rev. Applied 10, 054062 (2018).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, and Y. e. a. Chen, Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * [40] D. Weiss, H. Zhang, C. Ding, Y. Ma, D. I. Schuster, and J. Koch, Fast high-fidelity gates for galvanically-coupled fluxonium qubits using strong flux modulation, arXiv:2207.03971 .
  * Wen [2003b] X.-G. Wen, Quantum Orders in an Exact Soluble Model, Phys. Rev. Lett. 90, 016803 (2003b).
  * Nussinov and Ortiz [2009] Z. Nussinov and G. Ortiz, A symmetry principle for topological quantum order, Ann. Phys. 324, 977 (2009).
  * Brown  _et al._ [2011] B. J. Brown, W. Son, C. V. Kraus, R. Fazio, and V. Vedral, Generating topological order from a two-dimensional cluster state using a duality mapping, New J. Phys. 13, 065010 (2011).
  * Bonilla Ataides  _et al._ [2021] J. P. Bonilla Ataides, D. K. Tuckett, S. D. Bartlett, S. T. Flammia, and B. J. Brown, The XZZX surface code, Nat. Commun. 12, 2172 (2021).
  * [45] Google Quantum AI, Suppressing quantum errors by scaling a surface code logical qubit, arXiv:2207.06431 .
  * D’Alessandro [2007] D. D’Alessandro,  _Introduction to Quantum Control and Dynamics_ , Chapman & Hall/CRC Applied Mathematics & Nonlinear Science (CRC Press, 2007).
  * Lloyd [1996] S. Lloyd, Universal Quantum Simulators, Science 273, 1073 (1996).
  * Nielsen and Chuang [2004] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and quantum information_ (Cambridge University Press (Cambridge, UK), 2004).
  * Bravyi  _et al._ [2006] S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-Robinson Bounds and the Generation of Correlations and Topological Quantum Order, Phys. Rev. Lett. 97, 050401 (2006).
  * Trebst  _et al._ [2007] S. Trebst, P. Werner, M. Troyer, K. Shtengel, and C. Nayak, Breakdown of a Topological Phase: Quantum Phase Transition in a Loop Gas Model with Tension, Phys. Rev. Lett. 98, 070602 (2007).
  * Bravyi  _et al._ [2010] S. Bravyi, M. B. Hastings, and S. Michalakis, Topological quantum order: Stability under local perturbations, J. Math. Phys. 51, 093512 (2010).
  * Dusuel  _et al._ [2011] S. Dusuel, M. Kamfor, R. Orús, K. P. Schmidt, and J. Vidal, Robustness of a Perturbed Topological Phase, Phys. Rev. Lett. 106, 107203 (2011).
  * Jamadagni  _et al._ [2018] A. Jamadagni, H. Weimer, and A. Bhattacharyya, Robustness of topological order in the toric code with open boundaries, Phys. Rev. B 98, 235147 (2018).
  * Kitaev and Preskill [2006] A. Kitaev and J. Preskill, Topological Entanglement Entropy, Phys. Rev. Lett. 96, 110404 (2006).
  * Levin and Wen [2006] M. Levin and X.-G. Wen, Detecting Topological Order in a Ground State Wave Function, Phys. Rev. Lett. 96, 110405 (2006).
  * Hamma  _et al._ [2005] A. Hamma, R. Ionicioiu, and P. Zanardi, Bipartite entanglement and entropic boundary law in lattice spin systems, Phys. Rev. A 71, 022315 (2005).
  * Bombin [2010] H. Bombin, Topological Order with a Twist: Ising Anyons from an Abelian Model, Phys. Rev. Lett. 105, 030403 (2010).
  * Eckardt  _et al._ [2005] A. Eckardt, C. Weiss, and M. Holthaus, Superfluid-Insulator Transition in a Periodically Driven Optical Lattice, Phys. Rev. Lett. 95, 260404 (2005).
  * [59] O. Raii, F. Mintert, and D. Burgarth, Scalable quantum control and non-abelian anyon creation in a honeycomb lattice, arXiv:2205.10114 .
  * Yu  _et al._ [2008] J. Yu, S.-P. Kou, and X.-G. Wen, Topological quantum phase transition in the transverse Wen-plaquette model, Europhys. Lett. 84, 17004 (2008).
  * Wootton [2012] J. R. Wootton, Quantum memories and error correction, J. Mod. Opt. 59, 1717 (2012).
  * Eckardt and Anisimovas [2015] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a Floquet-space perspective, New J. Phys. 17, 093039 (2015).
  * Burum [1981] D. P. Burum, Magnus expansion generator, Phys. Rev. B 24, 3684 (1981).
  * Brinkmann [2016] A. Brinkmann, Introduction to average Hamiltonian theory. I. Basics, Concepts Magn. Reson. Part A 45A, e21414 (2016).
  * Blanes  _et al._ [2009] S. Blanes, F. Casas, J. Oteo, and J. Ros, The Magnus expansion and some of its applications, Phys. Rep. 470, 151 (2009).

##  Appendix A Interaction picture Hamiltonian

The Hamiltonian of the system in the Schrödinger picture describes a lattice
of two-level systems with tuneable nearest-neighbour hopping rate and subject
to resonant qubit control,

$\displaystyle\hat{H}(t)=$
$\displaystyle\sum_{\alpha}\Big{[}\frac{\omega_{\alpha}}{2}\hat{Z}_{\alpha}+\Omega_{\alpha}(t)\sin(\omega_{\alpha}t+\phi_{\alpha})\hat{Y}_{\alpha}\Big{]}$
$\displaystyle+\sum_{\langle\alpha,\beta\rangle}g_{\alpha\beta}(t)\cos(\omega_{\alpha\beta}t)(\hat{\sigma}_{\alpha}^{+}\hat{\sigma}_{\beta}^{-}+\hat{\sigma}_{\alpha}^{-}\hat{\sigma}_{\beta}^{+})\
,$ (25)

where $\hat{\sigma}_{\alpha}^{\pm}=(\hat{X}_{\alpha}\pm i\hat{Y}_{\alpha})/2$,
$\omega_{\alpha\beta}=\omega_{\alpha}-\omega_{\beta}$, $\Omega_{\alpha}(t)$
are the envelopes of the resonant qubit pulses, $g_{\alpha\beta}(t)$ a hopping
modulation of frequency $\omega$, and $\langle\alpha,\beta\rangle$ indicates
nearest neighbours. Indices $\alpha$ and $\beta$ denote pairs of lattice
coordinates. The characteristic frequencies of $\Omega_{\alpha}(t)$ and
$g_{\alpha\beta}(t)$ will be much smaller than the qubit frequencies
$\omega_{\alpha}$ and their differences $\omega_{\alpha\beta}$. The
Hamiltonian (25) can describe an architecture of superconducting qubits with
tuneable coupling and additional single-qubit microwave control, such as the
one implemented in Ref. [37]. In the interaction picture with respect to the
perturbed problem, the Hamiltonian 25 reads

$\displaystyle\hat{H}(t)=$
$\displaystyle\sum_{\alpha}\frac{\Omega_{\alpha}(t)}{2}(e^{i\phi_{\alpha}}\hat{\sigma}_{\alpha}^{+}+e^{-i\phi_{\alpha}}\hat{\sigma}_{\alpha}^{-}$
$\displaystyle-e^{i2\omega_{\alpha}t+i\phi_{\alpha}}\hat{\sigma}_{\alpha}^{+}-e^{-2i\omega_{\alpha}t-i\phi_{\alpha}}\hat{\sigma}_{\alpha}^{-})$
$\displaystyle+\sum_{\langle\alpha,\beta\rangle}g_{\alpha\beta}(t)\big{[}\hat{\sigma}_{\alpha}^{+}\hat{\sigma}_{\beta}^{-}+\hat{\sigma}_{\alpha}^{-}\hat{\sigma}_{\beta}^{+}$
$\displaystyle+e^{2i\omega_{\alpha\beta}t}\hat{\sigma}_{\alpha}^{+}\hat{\sigma}_{\beta}^{-}+e^{-2i\omega_{\alpha\beta}t}\hat{\sigma}_{\alpha}^{-}\hat{\sigma}_{\beta}^{+}\big{]}\
.$ (26)

The second and fourth line involve terms oscillating at frequencies
$\omega_{\alpha}$ and $\omega_{\alpha\beta}$, which are assumed to be much
larger than the frequencies in $\Omega_{\alpha}(t)$ and $g_{\alpha\beta}(t)$.
Hence, they can neglected in rotating-wave approximation, and a Hamiltonian of
the form of Eq. (1) is obtained. Note that, if single-qubit control was
possible for resonant qubits, without crosstalk issues, the cosinusoidal
modulation in Eq. (25) could directly implement the modulations of Eq.
(III.1), without needing the additional separation of timescales
$\omega_{\alpha\beta}\gg\omega$.

##  Appendix B Single plaquette driving scheme

We motivate in this Section the choice of the driving functions of Eq. III.1.

###  B.1 Oscillating part

First, we explore the problem from the point of view of a high-frequency limit
and related high-frequency expansions. This helps one to build some intuition
about the choice of the driving fields, which provides the base for proceeding
with numerically exact methods for designing the final protocol. The desired
four-spin term $\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}$ can be written
as the commutator of at least three two-spin operators,

$\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}\propto[\hat{X}_{1}\hat{X}_{2},[\hat{X}_{2}\hat{X}_{3}+\hat{Y}_{2}\hat{Y}_{3},\hat{X}_{3}\hat{X}_{4}]]\
.$ (27)

Double-commutator terms of this form appear in a Magnus expansion [31, 62, 32]
generated by the driven Hamiltonian $\hat{H}(t)$ at third order. In
particular, such terms have, e.g., the form

$\displaystyle[\hat{H}_{-m},[\hat{H}_{m-n},\hat{H}_{n}]]/\omega^{2},$
$\displaystyle[\hat{H}_{-m},[\hat{H}_{0},\hat{H}_{m}]]/\omega^{2}\ ,\dots$
(28)

where $\hat{H}_{m}=1/\tau\int_{0}^{\tau}e^{im\omega t}\hat{H}(t)$. Choosing
$\hat{H}_{0}=0$, since we want lower-order two-spin terms to be zero, leaves,
for example, the first term. A choice of values of $m$ and $n$ that makes such
a term non-zero is, for instance, $m=-1$, $n=2$. Then, comparing Eqs. (27) and
(28), one can start guessing a drive Hamiltonian of the form

$g(\omega t)\hat{X}_{1}\hat{X}_{2}+g(-3\omega
t)(\hat{X}_{2}\hat{X}_{3}+\hat{Y}_{2}\hat{Y}_{3})+g(2\omega
t)\hat{X}_{3}\hat{X}_{4},$ (29)

with periodic $g(t+\tau)=g(t)$. Now, since the goal is to make third-order
terms in the expansion to become of leading order, the second-order must
vanish. This can be achieved by exploiting a property of the Magnus expansion,
namely that all even order terms vanish if the driven Hamiltonian is time-
symmetric within the integration step $\tau$, $\hat{H}(t)=\hat{H}(\tau-t)$
[63, 64]. We will elaborate more on this point in subsection B.3. It is thus
convenient to choose $g(t)$ to be time-symmetric. We thus obtain driving
functions (using the labelling of Fig. 2)

$\displaystyle g_{13}(t)=g_{13}\cos(\omega t),\quad
g_{23}(t)=g_{23}\cos(3\omega t),$ $\displaystyle\quad
g_{24}(t)=g_{24}\cos(2\omega t),$ (30)

which yield the desired four-body term with strength $\mathcal{J}\propto
g_{13}g_{23}g_{24}/\omega^{2}$. The ansatz (30) is the starting point for the
numerical optimization, where testing slightly different harmonics yields in
the end optimal parameters with $\omega$ and $2\omega$ harmonics as in Eq.
(III.1). In the next subsection, we discuss the intuition behind the
introduction of the static parameters $\Omega_{1}$, $\Omega_{4}$,
$g_{kk^{\prime}}^{(0)}$.

###  B.2 Constant part

The Hamiltonian of Eq. (29) has a form similar to the Hamiltonian of Eq. (1),
except that it also involves two-qubit operators not conserving the total spin
magnetization, $\sum_{j}\hat{Z}_{j}$. If one is experimentally able to
modulate magnetization-nonconserving couplings like $\hat{X}\hat{X}$ and
$\hat{Y}\hat{Y}$ separately, then the ansatz (29) would be sufficient.
However, the Hamiltonian (1), as compared to Eq. (29), also contains terms
$\hat{Y}_{1}\hat{Y}_{2}$ and $\hat{Y}_{3}\hat{Y}_{4}$, which would lead to the
formation of unwanted effective terms
$\hat{Y}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Y}_{4}$, in addition to
$\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}$. These terms are undesired
since they do not commute with $\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}$
on different plaquettes, and they are thus detrimental. Their suppression in
the effective Hamiltonian can be achieved by introducing single-spin $\hat{X}$
drives on spin 1 and 4. The Hamiltonian (29), including the additional
$\hat{Y}\hat{Y}$ terms, then becomes

$\displaystyle\Omega_{1}\hat{X}_{1}+\Omega_{4}\hat{X}_{4}+g_{13}(t)(\hat{X}_{1}\hat{X}_{3}+\hat{Y}_{1}\hat{Y}_{3})$
$\displaystyle+g_{23}(t)(\hat{X}_{2}\hat{X}_{3}+\hat{Y}_{2}\hat{Y}_{3})+g_{24}(t)(\hat{X}_{2}\hat{X}_{4}+\hat{Y}_{2}\hat{Y}_{4}),$
(31)

By choosing $\Omega_{1}=\Omega_{4}=M\omega$ with integer $M$ and moving to a
rotating frame with respect to the $\hat{X}_{1}$ and $\hat{X}_{4}$ terms, the
Hamiltonian transforms to

$\displaystyle g_{13}(t)(\hat{X}_{1}\hat{X}_{2}+e^{iM\omega
t\hat{X}_{1}}\hat{Y}_{1}e^{-iM\omega t\hat{X}_{1}}\hat{Y}_{2})$
$\displaystyle+g_{23}(t)(\hat{X}_{2}\hat{X}_{3}+\hat{Y}_{2}\hat{Y}_{3})$
$\displaystyle+g_{24}(t)(\hat{X}_{2}\hat{X}_{4}+\hat{Y}_{2}e^{iM\omega
t\hat{X}_{4}}\hat{Y}_{4}e^{-iM\omega t\hat{X}_{4}})\ .$ (32)

Note that, thanks to the choice of $\Omega_{1}$ and $\Omega_{4}$ as a multiple
of $\omega$, the rotating frame coincides with the non-rotating frame at
$t=\tau$. One can now observe that the terms multiplied by exponential
functions oscillate at frequencies $\sim\pm(M\pm k)\omega$ with $k=1,2$, and
will thus produce terms in a high-frequency expansion that are damped at least
by a factor $1/(M\pm k)\omega$. By choosing $M$ sufficiently large, these
terms can be made smaller than the effective four-body
$\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4}$ term. A very large $M$ can
become impractical, and this is an additional reason for including the
parameters $\Omega_{j}$ in the optimization.

Finally, the Hamiltonian of Eq. (31) will likely generate, at higher orders in
the effective Hamiltonian, terms composed of the same operators which enter
the Hamiltonian itself. In the numerical optimization, small static
corrections $g_{kk^{\prime}}^{(0)}$ are also included as free parameters for
the purpose of cancelling these terms, thus obtaining in conclusion the ansatz
(III.1) for the control functions.

###  B.3 Effective reduction of the dynamical Lie algebra via time-symmetric
drive

In this subsection, we describe how the choice of time-symmetric drives,
together with the specific algebraic structure of the control problem, leads
to a situation where not all operators belonging to the dynamical Lie algebra
can actually contribute to the effective Hamiltonian. The latter property is
particularly advantageous, since it drastically reduces the number of
potential unwanted operators, while easing the numerical optimization of
parameters.

As the first step, we observe that the dynamical Lie algebra $\mathfrak{L}$,
spanned by the operators in Table 2, admits a Cartan decomposition [46]
$\mathfrak{L}=\mathfrak{h}\oplus\mathfrak{p}$, such that

$[\mathfrak{h},\mathfrak{h}]\subseteq\mathfrak{h},\quad[\mathfrak{h},\mathfrak{p}]\subseteq\mathfrak{p},\quad[\mathfrak{p},\mathfrak{p}]\subseteq\mathfrak{h}\
,$ (33)

with subspaces $\mathfrak{h}$ and $\mathfrak{p}$ spanned by

$\displaystyle\mathfrak{h}=\mathrm{span}\\{$
$\displaystyle\hat{Z}_{1}\hat{Y}_{2},\,\hat{Y}_{3}\hat{Z}_{4},\,\hat{X}_{1}\hat{Z}_{2}\hat{Y}_{3},\,\hat{Y}_{1}\hat{Z}_{2}\hat{X}_{3},\,$
$\displaystyle\hat{X}_{2}\hat{Z}_{3}\hat{Y}_{4},\,\hat{Y}_{2}\hat{Z}_{3}\hat{X}_{4},\,\hat{Z}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Y}_{4},$
$\displaystyle\hat{Y}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Z}_{4}\\}\ ,$ (34)
$\displaystyle\mathfrak{p}=\mathrm{span}\\{$
$\displaystyle\hat{X}_{1},\,\hat{X}_{4},\,\hat{X}_{1}\hat{X}_{2},$
$\displaystyle\hat{X}_{2}\hat{X}_{3},\,\hat{X}_{3}\hat{X}_{4},\
\hat{Z}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Z}_{4},$
$\displaystyle\hat{Y}_{1}\hat{Y}_{2},\,\hat{Y}_{3}\hat{Y}_{4},\,\hat{Z}_{1}\hat{Z}_{2}\hat{X}_{3},\,\hat{X}_{2}\hat{Z}_{3}\hat{Z}_{4},$
$\displaystyle\hat{X}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{X}_{4},\,\hat{Y}_{1}\hat{Z}_{2}\hat{Z}_{3}\hat{Y}_{4}\\}\
.$ (35)

The dynamics $\hat{U}(t)$ of the system governed by the Hamiltonian
$\hat{H}(t)$ can be formally expressed via a Magnus expansion [65],
$\hat{U}(t)=e^{\hat{M}(t)}$, characterized by the Magnus exponent
$\hat{M}(t)=\sum_{n=1}^{\infty}\hat{M}_{n}(t)$. The first terms read

$\displaystyle\hat{M}_{1}(t)$ $\displaystyle=-i\int_{0}^{t}\\!dt^{\prime}\
\hat{H}(t^{\prime})\ ,$ (36) $\displaystyle\hat{M}_{2}(t)$
$\displaystyle=\frac{(-i)^{2}}{2}\
\int_{0}^{t}\\!dt_{1}\int_{0}^{t_{1}}\\!\\!dt_{2}\
[\hat{H}(t_{1}),\hat{H}(t_{2})]\ .$ (37)

In particular, the $n$th term $\hat{M}_{n}(t)$ in the Magnus exponent involves
$n-1$ nested commutators of $\hat{H}(t)$,

$\Big{[}\hat{H}(t_{1}),\big{[}\hat{H}(t_{2}),\dots[\hat{H}(t_{n-1}),\hat{H}(t_{n})]\dots\big{]}\Big{]}\
,$ (38)

evaluated at different times. Since $\hat{H}(t)$ is a linear combination of
elements of $\mathfrak{p}$, it holds that $\hat{H}(t)\in\mathfrak{p}$. Hence,
from the properties (33) of the Cartan decomposition, one can conclude by
induction that

$\hat{M}_{2n}(t)\in\mathfrak{h},\quad\hat{M}_{2n+1}(t)\in\mathfrak{p}\ .$ (39)

Finally, choosing driving functions that are time-symmetric within the time
interval $\tau$ makes even terms $\hat{M}_{2n}(\tau)$ in the Magnus expansion
to vanish [63, 64], as mentioned in Sec. B.1. As a result, the effective
Hamiltonian is restricted to be a linear combination of operators belonging to
$\mathfrak{p}$ only, rather than to the whole Lie algebra $\mathfrak{L}$.

##  Appendix C Numerical optimization

The numerical optimization is performed by starting with different random
choices of drive parameters, computing the end-of-period propagator
$\hat{U}_{ij}(\tau)$ generated by the driven dynamics and minimizing

$\|\hat{U}_{ij}(\tau)-e^{i\tau\mathcal{J}\hat{P}_{i,j}}\|$ (40)

with a gradient-descent algorithm. Well converged results are obtained by
first setting $g_{\alpha\beta}^{(0)}=\Omega_{\alpha}=0$ and optimizing
$g_{\alpha\beta}^{(1)}$. Then, all parameters are optimized together starting
from initial random guesses that are stochastically perturbed variants of the
optimal $g_{\alpha\beta}^{(1)}$ found at the first step.

##  Appendix D Semi-analytical representation of operators

In order to efficiently represent and manipulate the operators for large
systems in numerical simulations, in particular for computing commutators, we
employ the following semi-analytical representation. A generic multi-qubit
operator

$\hat{O}=\sum_{\begin{subarray}{c}\hat{\sigma}_{\alpha_{j}}\in\\\
\\{\hat{X}_{\alpha_{j}},\hat{Y}_{\alpha_{j}},\hat{Z}_{\alpha_{j}},\hat{\mathbb{1}}_{\alpha_{j}}\\}\end{subarray}}O_{\sigma_{\alpha_{1}},\dots,\sigma_{\alpha_{N}}}\hat{\sigma}_{\alpha_{1}}\otimes\dots\otimes\hat{\sigma}_{\alpha_{N}},$
(41)

is represented by storing (i) the position of non-identity terms, (ii) the
type of Pauli matrix at each position (as a string), (iii) the coefficient
$O_{\sigma_{\alpha_{1}},\dots,\sigma_{\alpha_{N}}}$. For example, an operator
$\hat{O}=\omega_{1}\hat{Z}_{1}+\omega_{2}\hat{Z}_{2}+J\hat{X}_{1}\hat{X}_{2}$
is represented as

$\displaystyle O=\\{$ $\displaystyle[1\to`Z\text{'},\omega_{1}],$
$\displaystyle[2\to`Z\text{'},\omega_{2}],$
$\displaystyle[(1,2)\to(`X\text{'},`X\text{'}),J]\\}\ .$ (42)

This representation is convenient for manipulating Hamiltonians for large
systems since it allows one to perform matrix multiplications without needing
to store and multiply the matrices explicitly. This is particularly
advantageous in the case of Floquet-engineered effective Hamiltonians, since
they are typically not sparse, because of all the higher-order effective terms
produced by the Floquet drives. For example, the commutator of two operators
$\hat{O}_{1}$ and $\hat{O}_{2}$ is computed by simply selecting shared non-
identity sites between $O_{1}$ and $O_{2}$ and providing the analytical result
of the commutator of such terms.

##  Appendix E Construction of the ground-state preparation protocol

The aim of this Section is to explain how the ground state preparation
protocol proposed in Sec. IV can indeed produce the same results of Eq. (15).
This is based on three observations:

1) Since $\hat{X}$ and $-i\hat{Y}$ give the same result if applied to
$\ket{0}$, $\hat{X}\ket{0}=-i\hat{Y}\ket{0}=\ket{1}$, it is also true that

$(1+\hat{P}_{i,j})\ket{\psi_{0}}=(1-i\hat{A}_{i,j})\ket{\psi_{0}}=(1-i\hat{B}_{i,j})\ket{\psi_{0}}\
.$ (43)

2) $\hat{A}_{i,j}$ commutes with all other operators
$\hat{A}_{i^{\prime},j^{\prime}}$ and $\hat{B}_{i^{\prime},j^{\prime}}$ except
for $\hat{A}_{i+1,j+1}$ and $\hat{B}_{i-1,j-1}$; this implies that the
ordering of the product along diagonals spanning the system in direction
bottom-right (BR) to top-left (TL) do not matter, whereas the ordering of
products along diagonals in direction bottom-left (BL) to top-right (TR) do
matter. The correct order can be evicted from the next point.

3) In a product of two noncommuting terms $(1-i\hat{A}_{i+1,j+1})$ and
$(1-i\hat{A}_{i,j})$, the noncommuting parts satisfy

$\displaystyle-\hat{A}_{i+1,j+1}\hat{A}_{i,j}=$
$\displaystyle-i\hat{X}_{i,j}\hat{Z}_{i+1,j}Z_{i,j+1}\hat{Z}_{i+1,j+1}\hat{Z}_{i+1,j+1}$
$\displaystyle\otimes\hat{Z}_{i+2,j+1}\hat{Z}_{i+1,j+2}\hat{Y}_{i+2,j+2}\ .$
(44)

According to observation 1 and the fact that $\hat{Z}\ket{0}=\ket{0}$, the
operator in Eq. (44) acts on $\ket{\psi_{0}}$ in exactly the same way as
$\hat{P}_{i+1,j+1}\hat{P}_{i,j}$, such that

$(1-i\hat{A}_{i+1,j+1})(1-i\hat{A}_{i,j})\ket{\psi_{0}}=\\\
(1+\hat{P}_{i+1,j+1})(1+\hat{P}_{i,j})\ket{\psi_{0}}.$ (45)

Generalizing, the product of $n$ terms $-i\hat{A}_{i,j}$ along a diagonal
ordered from BL to TR has the same action on $\ket{\psi_{0}}$ as the
equivalent product of $\hat{P}_{i,j}$.