# Scalable and robust quantum computing on qubit arrays with fixed coupling

N.H. Le  Advanced Technology Institute and Department of Physics, University
of Surrey, Guildford GU2 7XH, United Kingdom  M. Cykiert  Advanced Technology
Institute and Department of Physics, University of Surrey, Guildford GU2 7XH,
United Kingdom  E. Ginossar  Advanced Technology Institute and Department of
Physics, University of Surrey, Guildford GU2 7XH, United Kingdom

###### Abstract

We propose a scheme for scalable and robust quantum computing on two-
dimensional arrays of qubits with fixed longitudinal coupling. This opens the
possibility for bypassing the device complexity associated with tunable
couplers required in conventional quantum computing hardware. Our approach is
based on driving a subarray of qubits such that the total multi-qubit
Hamiltonian can be decomposed into a sum of commuting few-qubit blocks, and
then efficient optimization of the unitary evolution within each block. Each
driving pulse can implement a target gate on the driven qubits, and at the
same time implement identity gates on the neighbouring undriven qubits,
cancelling any unwanted evolution due to the constant qubit-qubit interaction.
We show that it is possible to realise a universal set of quantum gates with
high fidelity on the basis blocks, and by shifting the driving pattern one can
realise an arbitrary quantum circuit on the array. Allowing for imperfect
Hamiltonian characterisation, we use robust optimal control to obtain
fidelities around 99.99% despite $1\%$ uncertainty in the qubit-qubit and
drive-qubit couplings, and a detuning uncertainty at 0.1% of the qubit-qubit
coupling strength. This robust feature is crucial for scaling up as parameter
uncertainty is significant in large devices.

## Introduction

Great progress has been achieved recently in various physical platforms for
quantum computing, most notably is the 54-qubit programmable superconducting
processor Arute  _et al._ (2019). High fidelity two-qubit gates were also
demonstrated for trapped ions Ballance  _et al._ (2016), neutral atoms in
optical tweezers Levine  _et al._ (2019), and spin qubits in silicon He  _et
al._ (2019); Huang  _et al._ (2019) and GaAs Nichol  _et al._ (2017). These
experimental implementations are based on tunable coupling between qubits
where the interaction is switched on only when the two qubit gates are
implemented. In solid state quantum computers, tunable couplers typically
involve more circuit elements and require their own external control for
tuning the magnitude of the interaction Arute  _et al._ (2019); Chen  _et al._
(2014); He  _et al._ (2019) leading to overheads in fabrication and wiring.
For solving real-world problems, a quantum computer needs a large number of
qubits Fowler  _et al._ (2012), and the complexity of the tunable couplers
adds to the technological difficulties in scaling up the device. In contrast,
fixed couplers do not require the extra components for controlling the
magnitude of the interaction, resulting in a substantial simplification of the
hardware architecture and hence a significant advantage for scaling up.

An important requirement for quantum computing with fixed coupling is the
ability to cancel the unwanted evolution due to the fixed interaction on
qubits where no gate is needed. In NMR quantum computing, where the qubits
have fixed longitudinal couplings, this is achieved by applying a series of
cleverly designed of refocusing pulses Vandersypen  _et al._ (2001); Tsunoda
_et al._ (2020). For large arrays of qubits these series become increasingly
complex, which is a bottleneck for scaling up Tsunoda  _et al._ (2020). In
this paper we describe a simple method for quantum computing on qubit arrays
with fixed coupling without refocusing pulses. Instead, we rely on a key
observation that, by driving a specific subarray, one can implement any gate
on the driven qubits, and _at the same time_ implement an identity operator on
all undriven qubits, effectively cancelling the unwanted evolution on these
undriven qubits. Any arbitrary quantum circuit can then be implemented by
changing the driven subarray between the time steps. An illustration of the
driving pattern and the implementation of gates is given in Fig. 1. Our method
can be scaled up to an arbitrarily large array in a straightforward manner,
opening an alternative pathway for a simplified quantum computer hardware
architecture based entirely on fixed coupling.

In principle, designing the subarray could be difficult. This is because
simulating a constantly interacting system of qubits is in general not
possible due to the exponential wall: the cost in memory and time increases
exponentially with the number of qubits. Thus, one cannot predict the unitary
gate implemented by a driving pulse. In our method this problem is avoided,
because the driven subarray can be chosen such that the total Hamiltonian of
the system can be decomposed into a sum of commuting blocks of only a few
qubits. Each block has a low dimensional Hilbert space, and thus its unitary
evolution can be simulated and optimized efficiently. This decomposition
exists when the qubit-qubit coupling term is longitudinal, i.e., diagonal in
the computational basis, for example, the ZZ interaction.

Figure 1:  a) An example of implementing quantum gates on an array of qubits
with fixed couplings. In any given step only a sub-array can be driven. This
sub-array can be chosen to satisfy specific requirements so that the array’s
Hamiltonian can be decomposed into commuting few-qubit blocks. Each drive can
implement a gate on the driven qubit, and through its combination with the
fixed qubit-qubit couplings also implements an identity operator on the
neighboring undriven qubits. In the next step a different subarray is driven
for implementing gates on a different set of qubits. Here C-X on two adjacent
qubits denotes the CNOT gate and $I$ the identity gate. b) Illustration of a
quantum circuit in our scheme. The key feature is that any idle interval
between gates can be used to realize with an active identity gate, for
preventing unwanted evolution due to the fixed couplings.

An appealing feature of our method is the robustness of the gates against
uncertainty in all the physical parameters of the array. By using robust
optimal control we find pulses for realizing gates with fidelities around
$99.99\%$ despite a $1\%$ uncertainty in all the qubit-qubit and drive-qubit
couplings, and a detuning uncertainty at 0.1% of the qubit-qubit coupling
strength. This robustness of the fidelity against uncertainties is crucial for
an architecture with entirely fixed couplers because it is not possible to
isolate a qubit or qubit pair for a precise measurement of the parameter
values, and hence there is always a significant residual uncertainty even
after the device characterisation process.

This paper is organised as follows: We first describe the key details of our
method, including the driving pattern to make the Hamiltonian decomposable,
the application of a universal set of gates, and the implementation of an
arbitrary circuit. Next, we show how to use optimal control to make the gates
robust against parameter uncertainty in the Hamiltonian. Finally, we discuss
potential physical realisations of our method.

## Results

We first describe our method for implementing an arbitrary quantum circuit on
qubit arrays with fixed longitudinal coupling. We consider a system of qubits
coupled by fixed nearest-neighbor longitudinal interaction, i.e., an
interaction that commutes with the bare qubit’s Hamiltonian. For simplicity we
choose the ZZ interaction, which has been realised experimentally for
superconducting qubits Johnson  _et al._ (2011); Collodo  _et al._ (2020).
When a subset of qubits is driven by external fields, the system’s Hamiltonian
is

$\displaystyle\mathcal{H}(t)=$
$\displaystyle-\sum_{j}\frac{\omega_{j}}{2}\sigma_{j}^{z}+\sum_{j\in\mathcal{L}}d_{j}E_{j}(t)\sigma^{x}_{j}+\sum_{jk}J_{jk}\sigma_{j}^{z}\sigma_{k}^{z},$
$\displaystyle E_{j}(t)=$ $\displaystyle\
E^{x}_{j}(t)\cos\left(\nu_{j}t\right)+E^{y}_{j}(t)\sin\left(\nu_{j}t\right),$

where $\nu_{j}$ is the frequency of the drive on the $j$-th qubit,
$E^{x}_{j}(t)$ and $E^{y}_{j}(t)$ the two quadratures of the field, $d_{j}$
the $j$-th qubit’s dipole matrix element, and $\mathcal{L}$ the driven subset.
Typically, the qubit’s transition energy, $\omega_{j}$, is much larger than
the interaction, $J_{jk}$, and hence $\ket{0,0,...,0}$ is the ground state of
the undriven Hamiltonian and can be initialised by cooling.

In the frame rotating with the qubits’ frequencies, described by the unitary
transformation $U_{0}(t)=e^{i\sum_{j}\frac{\omega_{j}}{2}\sigma^{z}_{j}t}$,
the Hamiltonian in the rotating wave approximation is

$\displaystyle H(t)$
$\displaystyle\approx\sum_{j\in\mathcal{L}}\frac{1}{2}\left[\Omega_{j}(t)\sigma^{x}_{j}+\Omega^{\prime}_{j}(t)\sigma^{y}_{j}\right]+\sum_{jk}J_{jk}\sigma_{j}^{z}\sigma_{k}^{z},$
(1)

where

$\displaystyle\Omega_{j}(t)=\Omega^{x}_{j}(t)\cos(\delta_{j}t)+\Omega^{y}_{j}(t)\sin(\delta_{j}t),$
$\displaystyle\Omega^{\prime}_{j}(t)=\Omega^{y}_{j}(t)\cos(\delta_{j}t)-\Omega^{x}_{j}(t)\sin(\delta_{j}t).$
(2)

Here, $\Omega^{x,y}_{j}(t)\equiv d_{j}E^{x,y}_{j}(t)$ is the Rabi frequency
and $\delta_{j}\equiv\nu_{j}-\omega_{j}$ is the detuning.

### Hamiltonian decomposition into commuting blocks

Computing the unitary evolution of a many-body Hamiltonian like $H(t)$ is in
general intractable due to the exponential complexity of the wave function,
unless one can decompose the Hamiltonian into a sum commuting few-qubit
blocks, i.e., $H(t)=\sum_{l}H_{l}(t)$, where all the $H_{l}(t)$ are mutually
commuting. The unitary evolution after a time duration $T$ is then
$U(T)=\prod_{l}U_{l}(T)$ where
$U_{l}(T)=\mathcal{T}e^{-i\int_{0}^{T}dt^{\prime}H_{l}(t^{\prime})}$ and
$\mathcal{T}$ is the time-ordering operator. $U_{l}(T)$ can be efficiently
computed since it involves only a few qubits. The $U_{l}(T)$ factors are also
mutually commuting and can be seen as parallel quantum gates applied on
separate qubit blocks.

Figure 2:  a) An arbitrarily large 2D array of coupled qubits in a honeycomb
array. Each link represents a ZZ coupling term. A yellow (cyan) vertex
represents a driven (undriven) qubit. The total Hamiltonian can be decomposed
into a sum of identical four-qubit blocks that commute with each other
(enclosed in the dashed circles). An undriven qubit is shared by three
neighboring blocks. (b) The pattern of driving for implementing a two-qubit
gate, resulting in a central row of identical six-qubit blocks (enclosed by
the dashed rectangles). The rest of the array can be decomposed into the four-
qubit blocks as in (a).

We find that the simplest geometry that allows the decomposition of $H(t)$
into few-qubit commuting blocks is a honeycomb array of qubits with nearest
neighbor ZZ coupling, as shown in Fig. 2. We consider an alternating driving
pattern where only the subarray colored in yellow in Fig. 2a is driven, then
$H(t)=\sum_{j\in\mathcal{L}}H_{j}(t)$ where $\mathcal{L}$ is the driven
subarray and

$\displaystyle
H_{j}(t)=\frac{1}{2}\left[\Omega_{j}(t)\sigma^{x}_{j}+\Omega^{\prime}_{j}(t)\sigma^{y}_{j}\right]+\sum_{k\in\text{NB}_{j}}J_{jk}\sigma_{j}^{z}\sigma_{k}^{z},$
(3)

where $\text{NB}_{j}$ is the set of the three nearest neighbors of qubit $j$.
Each of the block Hamiltonians, $H_{j}$, has only four qubits and they commute
with each other. This can be seen in a more transparent way by the graphical
representation in Fig. 2a. Each link in the graph represents a $ZZ$ coupling
term; a yellow vertex represents a single qubit driving term, and a cyan
vertex represents an undriven qubit and hence nothing in $H(t)$. All the ZZ
terms commute with each other. A yellow vertex does not commute only with the
three vertices connected to it, because $\sigma^{x}$ and $\sigma^{y}$ do not
commute with $\sigma^{z}$. Thus, the total Hamiltonian can be expressed as a
sum of commuting four-qubit blocks enclosed by the dash circles in Fig. 2a.
Note that each block has only one driven qubit in the center. The qubits at
the intersection of two neighboring blocks must not be driven for the
commutativity to hold. We will show below that it is possible to implement a
single qubit gate on the driven qubit without changing the state of the
undriven qubits, at the end of the gate, despite the permanent ZZ interaction
in the block.

The driving pattern need to be modified slightly for implementing two-qubit
gates. In a conventional device with tunable couplers the qubit-qubit
interaction is turned on only when a two-qubit entangling gate is applied. In
our case the ZZ coupling is always on, and in general it entangles all the
qubits at all times. However, we find that it is still possible to implement a
specific two-qubit entanging gate, for example the CNOT gate, between two
targeted qubits by driving both. Turning on the drives on two neighboring
qubits results in the pattern of Fig. 2b where the central row is built from
identical six-qubit blocks. The rest of the array can be driven in the
alternating pattern as before. The reader may wonder why the six-qubit blocks
are required for the entire central row when only one two-qubit gate is
needed. This is necessary for applying the identity operators on all undriven
qubits for cancelling the actions of the fixed ZZ coupling, which requires
that any undriven qubit must have at least one neighboring driven qubit (more
details below). If a link, $J_{jk}\sigma^{z}_{j}\sigma^{z}_{k}$, is not
connected to any driven qubit, then it commutes with every other terms in the
Hamiltonian, and its contribution to the total unitary evolution is simply the
factor $e^{-iJ_{jk}\sigma^{z}_{j}\sigma^{z}_{k}T}$, which cannot be cancelled
due to the absence of control.

### Applying gates using optimal control

Figure 3:  a) A star graph with undriven qubits on the boundary and a driven
subet of one or more qubits in the center. Each undriven qubit must be
connected to at least one driven qubit. By utilizing optimal control on the
driven subset it is possible to apply a unitary of the form
$U_{\mathcal{C}}\otimes I_{\mathcal{B}}$ where $U_{\mathcal{C}}$ is a target
unitary on the central driven subset, and $I_{\mathcal{B}}\equiv I\otimes
I\otimes\dots\otimes I$ is the identity in the Hilbert space of the undriven
qubits on the boundary (see text). b) Examples of applying a Hamadard gate on
the four-qubit block and CNOT (CX) gate on the six-qubit block of Fig. 2. In
both cases the identity operators are applied on all undriven qubits.

We now describe how to apply targeted gates on the driven qubits while at the
same time apply the identity operators on the neighboring undriven qubits.
Note that the four and six qubit blocks in Fig. 2a and 2b have the form of a
star graph where only a central subset of qubits is driven, as depicted in
Fig. 3a. Our method lies in the key numerical finding that, for such a star
graph, it is possible to use optimal control algorithm to find pulse shapes,
$\Omega_{j}^{x,y}(t)$ where $j\in\textit{driven subset}$, to implement a
unitary operation of the type $U_{\mathcal{C}}\otimes I_{\mathcal{B}}$, where
$U_{\mathcal{C}}$ is a unitary acting on the driven subset, and
$I_{\mathcal{B}}$ the identity matrix acting on the undriven subset at the
boundary. The net effect is that the gate $U_{\mathcal{C}}$ is applied to the
driven subset while the rest remains unchanged. If the driven subset has one
(two) qubit, then $U_{\mathcal{C}}$ is a single qubit (two-qubit) gate.

Obviously the qubits on the boundary are acted on by the ZZ interactions, and
hence their states are changed, during the pulse, but by choosing the right
shape one can use the combined effect of the central driving term and the ZZ
connectors to ensure that the identity operators are applied at then end of
the pulse, removing the ZZ interactions in a stroboscopic fashion. This can be
partly understood by looking at the Baker–Campbell–Hausdorff expansion. For
the four-qubit block with the Hamiltonian of Eq. (3) , for example, the
unitary evolution in a small time step is

$\displaystyle e^{-iH_{j}(t)dt}$
$\displaystyle=e^{-i\left(H^{d}_{j}(t)+H^{int}_{j}\right)dt}$
$\displaystyle\approx
e^{-iH_{j}^{d}(t)dt}e^{-iH^{int}_{j}dt}e^{-[H^{d}_{j}(t),H^{int}_{j}]dt^{2}/2},$
(4)

where
$H^{d}_{j}(t)\equiv(1/2)\left[\Omega_{j}(t)\sigma^{x}_{j}+\Omega^{\prime}_{j}(t)\sigma^{y}_{j}\right]$
is the driving term,
$H^{int}_{j}\equiv\sum_{k\in\text{NB}_{j}}J_{jk}\sigma_{j}^{z}\sigma_{k}^{z}$
the ZZ interaction terms, and $[H^{d}_{j}(t),H^{int}_{j}]$ the commutator of
the two. While the first term is responsible for applying a gate on the driven
qubit, the second is the unwanted evolution due to the ZZ interaction. Since
$[H^{d}_{j}(t),H^{int}_{j}]=\sum_{k\in\text{NB}_{j}}iJ_{jk}\left[\Omega_{j}(t)\sigma_{j}^{x}-\Omega^{\prime}_{j}(t)\sigma_{j}^{y}\right]\sigma_{k}^{z}$,
it is obvious that the third term allows partial control of the undriven
qubits, labelled by $k$, through shaping $\Omega_{j}^{x,y}(t)$; and we find
that this is sufficient for undoing the evolution due to the ZZ interaction.

Using optimal control we are able to obtain pulses for realising the unitary
operator $U_{\mathcal{C}}\otimes I_{\mathcal{B}}$ with maximum fidelity,
$F=1$, up to numerical precision, where $U_{\mathcal{C}}$ is the Hadamard,
$\pi/8$ and the _direct_ identity gate on the the single driven qubit of the
four-qubit block (we use “ _direct_ identity gate" to refer to an identity
gate applied on an driven qubit to differentiate it from the identity
operators applied on the undriven qubits). The same was achieved where
$U_{\mathcal{C}}$ is the CNOT gate and the direct two-qubit identity gate,
$I\otimes I$, on the two driven qubits of the six-qubit block of Fig. 3b. Note
that in all of these examples the identity operators are applied on the
undriven qubits. These one and two qubit gates form a unversal set, i.e., a
set from which any muti-qubit unitary can be approximated from with arbitrary
precision, _allowing the implementation of an arbitrary quantum circuit_
Nielsen and Chuang (2000). More details of the optimal control algorithm and
pulse shapes are given later where we discuss the robustness of the these
gates against parameter uncertainty in the Hamiltonian.

### Implementing quantum circuits

An implementation of quantum computation on the array is illustrated in Fig. 1
where the driven subarray is varied from one step to the next to apply the
target gate on the right qubits. As can be seen in Fig. 1 the identity
operators are applied on the undriven qubits at every step. Note that the
commuting blocks are not fixed, but are changed constantly during the
execution of a quantum circuit, depending on where the gates are applied and
whether they are single or two-qubit gates.

We show in Fig. 4 a simple example of how the driven/undriven subarrays and
the blocks are varied during the implementation of a simple quantum circuit.
Consider the following sequence of gates on two qubits, denoted by A and B: a
two-qubit gate on A and B, followed by a single qubit gate on A, and then a
single qubit gate on B. Note how the driven/undriven subarrays and the blocks
are changed at each step. In step 1 both qubits, A and B, are driven in a six-
qubit block, in step 2 only qubit A is driven in a four-qubit block and qubit
B now becomes an undriven qubit, and in step 3 qubit B is driven and qubit A
undriven. At each step gates can also be implemented in parallel on the driven
qubits other than A and B. This parallel processing helps reduce the number of
steps in a computation. If there is no gate on a driven qubit at a given step
one simply applies the direct identity gate to keep its state unchanged. The
undriven qubits are always subjected to the identity operators at all steps.
Following this example it is straightforward to derive the driving pattern for
an arbitrary quantum circuit.

Figure 4:  An illustration of how the blocks are changed/shifted during an
implementation of the simple quantum circuit described in the text.

In our method the undriven subarray is crucial for the Hamiltonian
decomposition into commuting blocks, but this means that there are always
qubits that have no gate at a given step, leading to an overhead in the number
of steps compared with conventional quantum computation with tunable couplers.
The exact amount of this overhead depends on the gate configuration of a
circuit. However, as shown in Fig. 2, the number of undriven qubits is no more
than half of the total qubits, and it can be shown that the overhead is at
worst a factor of 2, which does not change the computational complexity of a
quantum algorithm.

Now we discuss _readout_ for the qubit array which obviously takes a finite
duration. At the end of the computation the drives are turned off and the
system evolves according to the background Hamiltonian
$H_{0}=\sum_{j,k}J_{j,k}\sigma_{j}^{z}\sigma_{k}^{z}$. The final wavefunction
can be expanded in the computational basis,
$\ket{\psi_{f}}~{}=~{}\sum_{j_{1},...,j_{N}=0,1}c_{j_{1},...,j_{N}}\ket{j_{1},...,j_{N}}$.
Since $H_{0}$ consists of exclusively $\sigma^{z}$ terms,
$\ket{j_{1},...,j_{N}}$ is an eigenstate of $H_{0}$ and the evolution under
$H_{0}$ introduces only addition phases in the coefficients
$c_{j_{1},...,j_{N}}$. Since the solution of a computation is usually encoded
in the bit string of the final wavefunction, as is the case for the Shor’s
algorithm for factoring Shor (1997) and HHS algorithm for solving linear
systems Harrow  _et al._ (2009), the additional phase will not change the
result.

### Robustness against parameter uncertainty

$\Delta J/\bar{J}$ | $\Delta\alpha$ | $\Delta\delta/\bar{J}$ | $H\otimes I_{\mathcal{B}}$ | $\pi/8\otimes I_{\mathcal{B}}$ | $I\otimes I_{\mathcal{B}}$ | CX$\otimes I_{\mathcal{B}}$ | $I_{2}\otimes I_{\mathcal{B}}$
---|---|---|---|---|---|---|---
0% | 0% | 0% | 10 | 10 | 10 | 10 | 10
0.1% | 0.1% | 0.1% | 5.6 (4.5) | 5.5 (4.6) | 5.5 (4.7) | 4.4 (4.1) | 4.6 (4.3)
1% | 1% | 0.1% | 5.6 (2.8) | 5.5 (2.7) | 5.6 (2.7) | 3.7 (2.4) | 4.2 (2.3)
5% | 5% | 0.1% | 5.3 (1.4) | 5.3 (1.2) | 5.4 (1.1) | 3.2 (1.1) | 4.0 (0.9)
Table 1: Maximum worst-case fidelities. The figures shown are the exponents of
the infidelity, i.e., $-\log_{10}\left(1-\mathcal{F}_{\text{max}}\right)$,
which is the number of nines in $\mathcal{F}_{\text{max}}$. The fidelities are
calculated for the single qubit gates on the four-qubit block, and the two-
qubit gates on the six-qubit block of Fig. 2. Four levels of uncertainty in
the qubit-qubit coupling strength and control amplitude are considered, while
the uncertainty in the detuning is kept at $0.1\%$ of the coupling strength.
The maximum amplitude of the Rabi frequencies is constrained to less than
$10\bar{J}$. $I$ and $I_{2}$ are the direct identity gates for one and two
qubits, respectively. The gate duration is $T=2\pi/\bar{J}$, divided into
$M=100$ time bins. The figures inside the bracket are results obtained with
non-robust optimization.

We now describe the optimal control algorithm for maximising the fidelity of
the unitary $\mathcal{U_{C}}\otimes\mathcal{I_{B}}$ and how to make this
unitary robust against parameter uncertainty in the Hamiltonian. We divide the
pulse duration, $T$, into $M$ time bins of interval $\Delta t$. In each time
bin the field amplitudes are kept constant. The set of the Rabi frequencies
form the control vector, $\bm{c}=\\{\Omega_{jn}^{\mu}:1\leq n\leq
M;\mu=x,y;j\in\mathcal{C}\\}$, where $\mathcal{C}$ is the driven subset. The
unitary evolution $U_{\mathcal{G}}(T)$ of the star graph $\mathcal{G}$ is then
a function of $\bm{c}$. Each qubit-qubit coupling, $J_{jk}$, and detuning,
$\delta_{j}$, is allowed to vary independently in the uncertainty intervals
$[\bar{J}-\Delta J/2,\bar{J}+\Delta J/2]$, and
$[\bar{\delta}-\Delta\delta/2,\bar{\delta}+\Delta\delta/2]$, respectively. The
uncertainty in the Rabi frequencies can be caused by that in the dipole-matrix
elements, or a slow drift in the drive leading to changes in the field
amplitudes from one experiment to the next. This can be modelled by replacing
$\Omega_{j}^{x,y}(t)$ in $H_{\mathcal{G}}$ by $\alpha_{j}\Omega_{j}^{x,y}(t)$
where $\alpha_{j}$ is a dimensionless parameter that varies in the interval
$[1-\Delta\alpha/2,1+\Delta\alpha/2]$. Now the unitary $U_{\mathcal{G}}(T)$
also depends on $J_{jk}$, $\alpha_{j}$, and $\delta_{j}$. The qubit-qubit
coupling and dipole-matrix element of a qubit have to be measured or estimated
at the point of fabrication, so their values can change substantially after
more qubits are added to the array due to additional interaction or
experimental drift. In contrast, we find that it is possible to determine the
frequency of every qubit in the completed array with typically very precise
spectroscopic measurement. In the array the resonant frequency of each qubit
is shifted due to the ZZ interactions, but there exists a procedure of one and
two-photon absorption measurements that can be combined to cancel these shifts
and obtain the bare qubit frequency, $\omega_{j}$ (see Method). Thus, we
assume that the driving fields are tuned to resonance, $\bar{\delta}=0$, with
residual detuning uncertainties much smaller than the average qubit-qubit
interaction, $\Delta\delta_{j}/\bar{J}=0.1\%$, which is typical for
superconducting qubits foo . This small detuning uncertainty should be
achievable in most physical realisations of qubits owing to the high accuracy
of spectroscopic measurements.

Denote the set of these uncertain parameters as $\bm{v}$, then the robust
optimal control problem can be defined as a max-min optimization problem: We
find an optimal control that maximizes the minimum fidelity over $\bm{v}$,

$\displaystyle\mathcal{F}_{\text{max}}=\max_{\bm{c}}\mathcal{F}(\bm{c}),\quad\mathcal{F}(\bm{c})=\min_{v\in\mathcal{V}}F(\bm{c},\bm{v}),$
(5)

where $\mathcal{V}$ is the hypercube containing the possible values of
$\bm{v}$, and $F(\bm{c},\bm{v})$ the gate fidelity based on the trace distance

$\displaystyle
F(\bm{c},\bm{v})=\left|\frac{1}{D}\,\text{tr}\\!\left[U_{\mathcal{G}}^{\dagger}(T)\left(U_{\mathcal{C}}\otimes
I_{\mathcal{B}}\right)\right]\right|^{2},$ (6)

where $D$ is the dimension of the Hilbert space of $\mathcal{G}$,
$I_{\mathcal{B}}$ the identity matrix of the subset of undriven qubits on the
boundary, and $U_{\mathcal{C}}$ the target unitary that we want to apply on
the central driven subset. In numerical computation one chooses a set of
sampling points, $\bm{v}_{i}$ , in $\mathcal{V}$, and find the minimum
fidelity in this set. We found that when the uncertainties are all smaller
than $5\%$ and $F(\bm{c},\bm{v})$ is larger than $99\%$ its minimum over
$\bm{v}$ always lies at one of the extreme points of $\mathcal{V}$, i.e., one
of the corners of the hypercube. Thus, we can redefine
$\mathcal{F}(\bm{c})\equiv\min_{\bm{v}_{i}\in\mathcal{X}}F(\bm{c},\bm{v}_{i})$
where $\mathcal{X}$ is the discreet set of the extreme points of
$\mathcal{V}$. This drastically reduces the number of sampling points where
$F(\bm{c},\bm{v})$ has to be computed. There are $2^{n_{u}}$ extreme points in
a hypercube of $n_{u}$ uncertain parameters. For example, the 6-qubit block in
Fig. 2b has 9 uncertain parameters, five $J$s, two $\alpha$s, and two
$\delta$s, so we need to compute $2^{9}\equiv 512$ values of
$F(\bm{c},\bm{v}_{i})$ for any given $\bm{c}$.

We develop a numerical computation for gradient-based robust optimal control
that can handle systems with up to 12 qubits and multiple uncertain
parameters. Starting with a random initial guess for $\bm{c}$, we use
gradients to identify a step $\delta\bm{c}$ to maximize
$\min_{i}\nabla_{\bm{c}}{F(\bm{c},\bm{v}_{i})}.\delta\bm{c}$, the first order
increment, so that $F(\bm{c},\bm{v}_{i})$ is increased for all $\bm{v}_{i}$.
In this way $\mathcal{F}(\bm{c})$, the minimum fidelity over $\bm{v}_{i}$, can
be increased to a value very close to one. In practice we find that one can
also raise $\mathcal{F}(\bm{c})$ by maximizing the average fidelity,
$\sum_{i=1}^{n_{\mathcal{X}}}F(\bm{c},\bm{v}_{i})/n_{\mathcal{X}}$, with
gradient ascent, which in some cases works faster than trying to increase
$F(\bm{c},\bm{v}_{i})$ for all $\bm{v}_{i}$. For the $n$-th time bin the time
evolution is calculated by the mid-point rule
$e^{-iH_{\mathcal{G}}(t_{n}-\Delta t/2)\Delta t}+O(\Delta t^{3})$ which is
accurate when the time step is small. The matrix exponentiation is sped up by
using the Krylov subspace method on sparse matrices; and the gradient computed
from a simple and efficient second-order formula de Fouquieres  _et al._
(2011).

For the four-qubit block of Fig. 2a, we derive optimal pulses to realize
$U_{\mathcal{C}}\otimes I_{\mathcal{B}}$ where $U_{\mathcal{C}}$ is the
Hadamard and $\pi/8$ gates Nielsen and Chuang (2000). And for the six-qubit
block of Fig. 2b we want $U_{\mathcal{C}}$ to be the CNOT gate. These three
gates form a universal set where any multi-qubit unitary can be approximated
from with arbitrary precision Nielsen and Chuang (2000), allowing universal
quantum computing on the array. For a system with exclusively fixed coupling,
it is also necessary to implement the identity gate that keeps all the qubits
in a block unchanged despite the permanent interaction.

In Table 1 we show the robust fidelities obtained for the universal set and
the identity gates at various levels of uncertainty. The gate duration is
divided into $M=100$ time bins. For $1\%$ uncertainty the fidelity is higher
than 99.999% for the single qubit gates and exceeds 99.98% for the two-qubit
gates. Even if the uncertainty is as high as 5% five-nines fidelities are
still achieved for the single qubit gates, and above 99.94% for the two-qubit
gates. This can be improved by increasing the number of initial guesses,
relaxing the constraints, or raising the number of control variables. The
optimal pulse shape for the Hadamard gate is shown in Fig. 5, and the pulse
shapes for the other gates in Table 1 is given in the Supplemental Materials.
We choose $T=2\pi/\bar{J}$ but this specific value of the duration is not
essential; the same order of magnitude is achieved for the fidelities when $T$
is changed by 10%. In order to see the effectiveness of robust optimal control
we also calculate the fidelities with non-robust optimal control: We first
neglect all the uncertainties and optimize the fidelity for the ideal case
where $J_{j}=\bar{J}$, $\alpha_{j}=1$ and $\delta_{j}=0$ for all $j$, then we
use the obtained optimal control, $\bm{c}_{\text{ideal}}$, to calculate the
minimum fidelity in the hypercube,
$\mathcal{F}=\min_{v\in\mathcal{V}}F(\bm{c}_{\text{ideal}},\bm{v})$. The
results are shown in the parentheses of Table 1. robust optimization improves
the fidelities by two to three nines when the uncertanties are significantly
large ($1\%$ and $5\%$).

Figure 5:  Optimal pulse for the Hadamard gate at $1\%$ uncertainty in Table
1.

## Discussion

A promising physical realisation of our model is the superconducting qubits
including the flux and transmon qubits. A direct ZZ interaction between flux
qubits can be realised by coupling the qubits inductively, as demonstrated in
quantum annealers Johnson  _et al._ (2011); Tennant  _et al._ (2021). There is
another interesting scheme based on the inductive longitudinal coupling of the
flux qubits with a common bus cavity Billangeon  _et al._ (2015); Didier  _et
al._ (2015), which can be scaled up to 2D arrays Richer and DiVincenzo (2016).
In addition, a cross-Kerr ZZ interaction has been recently demonstrated for
transmon qubits using a flux-tunable coupler Collodo  _et al._ (2020). While
flux qubits are very good two-level systems and hence our results are
applicable, for transmons the leakage to higher excited states must be
accounted for in the optimal control algorithm Motzoi  _et al._ (2009). The ZZ
coupling is also the natural interaction in a nuclear magnetic resonance
quantum computer Vandersypen  _et al._ (2001), for which sophisticated pulse-
shaping is available Yang  _et al._ (2020), making it a good test bed for our
model.

In experimental realisations, the pulse duration, $T$, is limited by the
coherence time, $T_{2}$, of the qubits. The decoherence rate of a block of
$N^{\prime}$ qubits is enhanced by a factor of approximately $N^{\prime}$ in
the worst case, giving rise to a lifetime of
$T_{2}^{\prime}=T_{2}/N^{\prime}$. The fidelity of a multi-qubit unitary on
the block is then bounded by $F\leq 1-T/T_{2}^{\prime}\equiv
1-N^{\prime}T/T_{2}$. Therefore, to achieve a fidelity $F$ for the four-qubit
and six-qubit blocks in Table 1 the pulse duration needs to be shorter than
$T_{2}(1-F)/N^{\prime}$ where $N^{\prime}=4$ and 6, respectively.

Although a honeycomb array is the focus of this paper, the qubits can be
arranged in any physical shape that has the same connectivity, for example a
square array with each qubit connected to only three nearest neighbors.
Moreover, there exist driving patterns that satisfy the conditions of
commutitativity and robust control for other geometries such as square arrays
and one-dimensional chains (see Supplemental Materials). One can also envisage
a hybrid architecture where large clusters of fixed couplers are connected
with tunable couplers, keeping the number of required tunable couplers low.
Such a modular structure can help ease the technological difficulties in
scaling quantum computers.

To conclude, we find that it is feasible to implement quantum computing with
accurate operations on a 2D qubit array with exclusively fixed couplers. The
quantum gates are robust against significant uncertainty in the qubit’s
frequency, qubit-qubit and drive-qubit coupling caused by fabrication
imperfection and/or slowly fluctuating fields. Both 1D and 2D geometries are
possible with a robust optimization process capable of handling the minimal
cluster size required for each scheme. Our proposal shows that scalability can
be accelerated with simplified hardware architecture based on fixed
longitudinal coupling schemes, thus motivating further development of this
coupling in different physical platforms.

## Method

### Calculating fidelity and gradient

We first describe how the fidelity and its gradients are calculated with the
midpoint rule. The time duration is divided into $M$ equal time bins with
$t_{0}=0$ and $t_{M}=T$. The field amplitudes are kept constant during each
time bin. The Hamiltonian of a star graph, $\mathcal{G}$, at the midpoint of
the $n-$th interval from $t_{n-1}$ to $t_{n}$ is

$H_{\mathcal{G},n}=\sum_{j\in\mathcal{C}}\alpha_{j}\left[\Omega_{jn}\sigma^{x}_{j}+\Omega^{\prime}_{jn}\sigma^{y}_{j}\right]+\sum_{j,k\in\mathcal{G}}J_{jk}\sigma_{j}^{z}\sigma_{k}^{z},$

where $\alpha_{j}$ is the dimensionless factor introduced to model the
uncertainty in the Rabi frequencies (see the Main Text), $\mathcal{C}$ is the
driven subset in the center of the graph and $J_{jk}\neq 0$ only for nearest
neighbors, and

$\displaystyle\Omega_{jn}$
$\displaystyle=\Omega^{x}_{jn}\cos(\delta_{j}(t_{n}-\Delta
t/2))+\Omega^{y}_{jn}\sin(\delta_{j}(t_{n}-\Delta t/2)),$
$\displaystyle\Omega^{\prime}_{jn}$
$\displaystyle=\Omega^{y}_{jn}\cos(\delta_{j}(t_{n}-\Delta
t/2))-\Omega^{x}_{jn}\sin(\delta_{j}(t_{n}-\Delta t/2)),$

where $\Omega^{x,y}_{jn}$ are the Rabi frequencies of the driving field on the
$j$-th driven qubit during the interval from $t_{n-1}$ to $t_{n}$. They are
the elements of the $2MN_{\mathcal{C}}\times 1$ control vector, $\bm{c}$,
where $N_{\mathcal{C}}$ is the number of qubits in the driving subset. The
unitary evolution from $t_{n-1}$ to $t_{n}$,
$U_{n}=e^{-iH_{\mathcal{G},n}\Delta t}+O(\Delta t^{3})$, is computed with the
expm function in Matlab. For an efficient calculation of the fidelity and the
gradients we compute and store all the $U_{n}$, and then obtain the forward
and backward unitary propagation operators de Fouquieres  _et al._ (2011);
Motzoi  _et al._ (2011), defined by

$\displaystyle U^{f}_{n}$ $\displaystyle=U_{n}U_{n-1}\dots U_{1},$
$\displaystyle U^{b}_{n+1}$ $\displaystyle=U_{M}U_{M-1}\dots U_{n+1},$

using the recursive relations $U^{f}_{n}=U_{n}U^{f}_{n-1}$ and
$U^{b}_{n+1}=U^{b}_{n+2}U_{n+1}$. Then the fidelity is

$\displaystyle
F(\bm{c},\bm{v})=\left|\frac{1}{2^{N_{\mathcal{G}}}}\,\text{tr}\\!\left[{U_{M}^{f}}^{\dagger}\left(U_{\mathcal{C}}\otimes
I_{\mathcal{B}}\right)\right]\right|^{2},$

where $N_{\mathcal{G}}$ is the number of qubits in the star graph.

For computing the gradients, we note that

$\displaystyle
H_{\mathcal{G},n}=\sum_{\mu=x,y}\sum_{j\in\mathcal{C}}\Omega^{\mu}_{jn}K^{\mu}_{jn}+\sum_{j,k\in\mathcal{G}}J_{jk}\sigma_{j}^{z}\sigma_{k}^{z},$

where

$\displaystyle K^{x}_{jn}$
$\displaystyle=g_{j}\left[\sigma_{j}^{x}\cos(\delta_{j}(t_{n}-\Delta
t/2))-\sigma_{j}^{y}\sin(\delta_{j}(t_{n}-\Delta t/2))\right],$ $\displaystyle
K^{y}_{jn}$
$\displaystyle=\alpha_{j}\left[\sigma_{j}^{y}\cos(\delta_{j}(t_{n}-\Delta
t/2))+\sigma_{j}^{x}\sin(\delta_{j}(t_{n}-\Delta t/2))\right].$

One can show that the derivative of $U_{n}\equiv e^{-iH_{\mathcal{G},n}\Delta
t}$ with respect to $\Omega^{\mu}_{jn}$ is de Fouquieres  _et al._ (2011)

$\displaystyle\frac{\partial
U_{n}}{\partial\Omega^{\mu}_{jn}}=\left\\{-i\Delta tK^{\mu}_{jn}+\frac{\Delta
t^{2}}{2}\left[H_{\mathcal{G},n},K^{\mu}_{jn}\right]\right\\}U_{n}+O(\Delta
t^{3}),$

where $\left[H_{\mathcal{G},n},K^{\mu}_{j,n}\right]$ is a commutator. It
follows that the derivative of $U_{M}^{f}\equiv U^{b}_{n+1}U_{n}U^{f}_{n-1}$
is

$\displaystyle\frac{\partial
U_{M}^{f}}{\partial\Omega^{\mu}_{jn}}=U^{b}_{n+1}\left\\{-i\Delta
tK^{\mu}_{jn}+\frac{\Delta
t^{2}}{2}\left[H_{\mathcal{G},n},K^{\mu}_{jn}\right]\right\\}U^{f}_{n},$

and from this it is straight forward to compute the gradient of the fidelity.
The most computationally expensive part of the calculation is the matrix
exponentiation for obtaining $U_{n}$, which is done $M$ times.

### Robust optimisation

We use two algorithms to raise the worst-case fidelity, $\mathcal{F}(\bm{c})$,
defined in Eq. 5. The first is based on sequential convex programming Allen
_et al._ (2017). We start with initial guesses $\bm{c}_{0}$ and $\bm{u}_{0}$
for the control vector and the upper limit (trust region) of the step,
respectively. Then a step $|\delta\bm{c}|<\bm{u}_{0}$ is found to maximize
$\min_{i}\nabla_{\bm{c}}F(\bm{c},\bm{v}_{i}).\delta\bm{c}$, i.e., to maximize
the minimum first-order increment over $\bm{v}_{i}$. This ensures all the
fidelities at the extreme points of the hypercube, $\mathcal{V}$, are
increased. The above optimization problem can be solved by sequential convex
programming (SCP). We used the YALMIP toolbox and SPDT3 package in Matlab for
this purpose. If a step can be found such that
$\min_{i}\nabla_{\bm{c}}F(\bm{c},\bm{v}_{i}).\delta\bm{c}$ is positive then we
increase the upper bound $\bm{u}_{0}$ by 1.15, otherwise we decrease it by 2.
We chose these factors as they give the fastest convergence in our numerical
tests. The procedure is repeated until either the maximum iteration is reached
or the step’s upper bound drops below a small tolerance. The second approach
is to simply maximize the average fidelity,

$\displaystyle\bar{\mathcal{F}}(\bm{c})=\sum_{i=1}^{n_{\mathcal{X}}}F(\bm{c},\bm{v}_{i})/n_{\mathcal{X}},$
(7)

using a quasi-Newton method. Obviously this does not guarantee that the worst-
case fidelity over $\bm{v}$ is increased, as the mean can be increased without
increasing the minimum value in the set. However, we found that in our
calculations the worst-case fidelity is always improved substantially when we
maximize the average fidelity. We optimize $\bar{\mathcal{F}}(\bm{c})$ using
the interior-point method implemented in Matlab’s fmincon function, where the
Hessian is computed from the exact gradients with the BFGS approximation.

In our numerical tests the first algorithm is more sensitive on the initial
guesses of the control parameters. For the two-qubit gates in Table I the
computation is very expensive and hence it is not practical to run the
optimization with too many initial guesses. We find that for the same running
time the second algorithm gives higher fidelities, and the results in Table I
are obtained with it.

### Spectroscopic measurement of qubit frequency

Figure 6:  The frequency of qubit 1 can be determined by one and two-photon
absorption measurements on the four qubits of the numbered unit cell (see
text).

In this section we propose a procedure using one and two-photon absorption
measurements for determining the frequency of every qubit in the 2D honeycomb
array. Initially all qubits in the array are in the ground state, $\ket{0}$.
We want to measure the frequency of qubit 1 in Fig. 6. Due to the ZZ
interactions with neighboring qubits causing a dispersive shift the one-photon
absorption peak of qubit 1 is at

$\displaystyle\omega^{p}_{1}=\omega_{1}-2(J_{12}+J_{13}+J_{14}).$ (8)

The shift is the energy cost for flipping three bonds from $\ket{00}$ to
$\ket{01}$. Similarly, the one-photon absorption peaks for the other numbered
qubits are

$\displaystyle\omega^{p}_{2}$
$\displaystyle=\omega_{2}-2J_{21}-2\sum_{\begin{subarray}{c}k\in\text{NB}_{2}\\\
k\neq 1\end{subarray}}{J_{2k}},$ $\displaystyle\omega^{p}_{3}$
$\displaystyle=\omega_{3}-2J_{31}-2\sum_{\begin{subarray}{c}k\in\text{NB}_{3}\\\
k\neq 1\end{subarray}}J_{3k},$ $\displaystyle\omega^{p}_{4}$
$\displaystyle=\omega_{4}-2J_{41}-2\sum_{\begin{subarray}{c}k\in\text{NB}_{4}\\\
k\neq 1\end{subarray}}J_{4k},$ (9)

where $\text{NB}_{j}$ is the set of nearest neighbors of qubit $j$. The last
sums are the dispersive shifts caused by the ZZ intreaction with neighbors
that are not qubit 1. Summing all the equations in Eqs. (8) and (Spectroscopic
measurement of qubit frequency), and using $J_{kl}=J_{lk}$, one obtains

$\displaystyle\sum_{j=1}^{4}\omega^{p}_{4}=\sum_{j=1}^{4}\omega_{j}-4\sum_{j=2}^{4}J_{1k}-2\sum_{j=2}^{4}\sum_{\begin{subarray}{c}k\in\text{NB}_{j}\\\
k\neq 1\end{subarray}}J_{jk}.$ (10)

If both qubits 1 and 2 are subjected to the the same field, for example by
splitting the field in a Mach–Zehnder interferometer set-up with one qubit in
each arm, one can observe the two-photon transition
$\ket{0}_{1}\ket{0}_{2}\rightarrow\ket{1}_{1}\ket{1}_{2}$. Now the coupling
energy of the bond between qubit 1 and 2 are not changed, and the two-photon
resonance happens at the frequency

$\displaystyle 2\omega^{p}_{12}$
$\displaystyle=\omega_{1}+\omega_{2}-2(J_{13}+J_{14})-2\sum_{\begin{subarray}{c}k\in\text{NB}_{2}\\\
k\neq 1\end{subarray}}{J_{2k}},$ (11)

Similarly, the two-photon resonance frequency for the other qubit pairs are

$\displaystyle 2\omega^{p}_{13}$
$\displaystyle=\omega_{1}+\omega_{3}-2(J_{12}+J_{14})-2\sum_{\begin{subarray}{c}k\in\text{NB}_{3}\\\
k\neq 1\end{subarray}}{J_{3k}},$ $\displaystyle 2\omega^{p}_{14}$
$\displaystyle=\omega_{1}+\omega_{4}-2(J_{12}+J_{13})-2\sum_{\begin{subarray}{c}k\in\text{NB}_{4}\\\
k\neq 1\end{subarray}}{J_{4k}}.$ (12)

Summing all the euqations in Eqs. (11) and (Spectroscopic measurement of qubit
frequency), one obtains

$\displaystyle
2\sum_{j=2}^{4}\omega^{p}_{1j}=2\omega_{1}+\sum_{j=1}^{4}\omega_{j}-4\sum_{j=2}^{4}J_{1k}-2\sum_{j=2}^{4}\sum_{\begin{subarray}{c}k\in\text{NB}_{j}\\\
k\neq 1\end{subarray}}J_{jk}.$ (13)

Finally, subtracting Eq. (10) from Eq. (13), the dispersive shifts are
cancelled and one arrives at

$\displaystyle\omega_{1}=\sum_{j=2}^{4}\omega^{p}_{1j}-\frac{1}{2}\sum_{j=1}^{4}\omega^{p}_{j},$

which gives the bare frequency of qubit 1 in terms of the one and two-photon
resonant frequencies. In most qubit platforms spectroscopic measurements are
precise and thus the uncertainty in the qubit’s frequency can be small.

### Data Availability

The data for this work are available without restriction zen .

### Code Availability

The Matlab code for this work are available without restriction zen .

## References

  * Arute  _et al._ (2019) F. Arute _et al._ , Nature 574, 505 (2019).
  * Ballance  _et al._ (2016) C. Ballance, T. Harty, N. Linke, M. Sepiol, and D. Lucas, Physical Review Letters 117, 060504 (2016).
  * Levine  _et al._ (2019) H. Levine _et al._ , Physical Review Letters 123, 170503 (2019).
  * He  _et al._ (2019) Y. He, S. K. Gorman, D. Keith, L. Kranz, J. G. Keizer, and M. Y. Simmons, Nature 571, 371 (2019).
  * Huang  _et al._ (2019) W. Huang _et al._ , Nature 569 (2019), 10.1038/s41586-019-1197-0.
  * Nichol  _et al._ (2017) J. M. Nichol, L. A. Orona, S. P. Harvey, S. Fallahi, G. C. Gardner, M. J. Manfra, and A. Yacoby, npj Quantum Information 3 (2017), 10.1038/s41534-016-0003-1.
  * Chen  _et al._ (2014) Y. Chen _et al._ , Physical Review Letters 113, 220502 (2014).
  * Fowler  _et al._ (2012) A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Physical Review A 86, 032324 (2012).
  * Vandersypen  _et al._ (2001) L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, M. H. Sherwood, and I. L. Chuang, Nature 414, 883 (2001).
  * Tsunoda  _et al._ (2020) T. Tsunoda, G. Bhole, S. A. Jones, J. A. Jones, and P. J. Leek, Physical Review A 102, 032405 (2020).
  * Johnson  _et al._ (2011) M. W. Johnson _et al._ , Nature 473, 194 (2011).
  * Collodo  _et al._ (2020) M. C. Collodo, J. Herrmann, N. Lacroix, C. K. Andersen, A. Remm, S. Lazar, J.-C. Besse, T. Walter, A. Wallraff, and C. Eichler, Phys. Rev. Lett. 125, 240502 (2020).
  * Nielsen and Chuang (2000) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ , Vol. 1 (Cambridge University Press, 2000).
  * Shor (1997) P. W. Shor, SIAM Journal on Computing 26, 1484 (1997).
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Physical Review Letters 103, 150502 (2009).
  * (16) For superconducting qubits with $\sim 10$ kHz frequency drift and $\sim 10$ MHz coupling, $\Delta\delta/J\sim$ 0.1%.
  * de Fouquieres  _et al._ (2011) P. de Fouquieres, S. Schirmer, S. Glaser, and I. Kuprov, Journal of Magnetic Resonance 212, 412 (2011).
  * Tennant  _et al._ (2021) D. M. Tennant _et al._ , “Demonstration of long-range correlations via susceptibility measurements in a one-dimensional superconducting Josephson spin chain,”  (2021), arXiv:2111.04284 [quant-ph].
  * Billangeon  _et al._ (2015) P.-M. Billangeon, J. S. Tsai, and Y. Nakamura, Physical Review B 91, 094517 (2015).
  * Didier  _et al._ (2015) N. Didier, J. Bourassa, and A. Blais, Physical Review Letters 115, 203601 (2015).
  * Richer and DiVincenzo (2016) S. Richer and D. DiVincenzo, Physical Review B 93, 134501 (2016).
  * Motzoi  _et al._ (2009) F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm, Physical Review Letters 103, 110501 (2009).
  * Yang  _et al._ (2020) X.-D. Yang, C. Arenz, I. Pelczer, Q.-M. Chen, R.-B. Wu, X. Peng, and H. Rabitz, Physical Review A 102, 062605 (2020).
  * Motzoi  _et al._ (2011) F. Motzoi, J. M. Gambetta, S. T. Merkel, and F. K. Wilhelm, Physical Review A 84, 022307 (2011).
  * Allen  _et al._ (2017) J. L. Allen, R. Kosut, J. Joo, P. Leek, and E. Ginossar, Physical Review A 95, 042325 (2017).
  * (26) Data, to be uploaded to zenodo.

### Acknowledgements

This work is supported by the UK Hub in Quantum Computing and Simulation, part
of the UK National Quantum Technologies Programme with funding from UKRI EPSRC
grant EP/T001062/1, and the EPSRC strategic equipment grant no. EP/L02263X/1.

### Author contributions

N.H.L and E.G developed the method for quantum computation. N.H.L and M.C
worked on the optimisation algorithms. E.G and N.H.L designed the project. All
authors wrote the paper and discussed its implications.

### Additional information

Competing interests: The authors declare no competing interests.

## Supplemental Materials

### Pulse shapes for the gates in Table 1

In Figs. S1 and S2 below we provide the pulse shapes for the single qubit and
two-qubit gates in Table I. All the axis units and plot legends are the same
as in Fig. 5.

Figure S1:  Pulse shapes for the single qubit gates in Table I. a) Hadamard
gate at 1% uncertainty (left) and 5% uncertainty (right) in both the control
amplitude and qubit-qubit coupling. b) $\pi/8$ gate at 1% uncertainty (left)
and 5% uncertainty (right). c) active identity gate, I , at 1% uncertainty
(left) and 5% uncertainty (right). The detuning uncertainty is kept at 0.1% of
the qubit-qubit coupling strength in all cases.

Figure S2:  Pulse shapes for the two-qubit gates in Table I. a) CNOT gate at
1% uncertainty (top row) and 5% uncertainty (bottom row) in both the control
amplitude and qubit-qubit coupling. b) active two-qubit identity gate
$\text{{I }}\otimes\text{{I }}$, at 1% uncertainty (top row) and 5%
uncertainty (bottom row). The detuning uncertainty is kept at 0.1% of the
qubit-qubit coupling strength in all cases. Note that for qubit 1 the two
quadratures of the pulses have the same shape, hence they overlap in the plot.

### Commuting blocks for one dimensional chains and square arrays

In the main text we focus on the honeycomb array because it is the 2D graph
with the smallest number of nearest neighbors (3), resulting in smallest
possible building blocks in 2D. Here we describe how to drive 1D chains and 2D
square arrays so that the Hamiltonian can be decomposed into small commuting
building blocks. The result is illustrated in the graphs of Fig. S3 where the
yellow vertices are driven qubits and cyan vertices undriven ones. For the
square array the building block for implementing a target two-qubit gate is
thirteen-qubit block where the two-qubit gate can be applied to any two
nearest neighbors in the five-qubit driven subset. We could not find a simpler
block due to the requirements that both qubits subjected to the two-qubit gate
have to be driven, and that for robust control every ZZ link has to be
connected to at least one driven qubit. We found that any driving pattern with
a smaller block results in at least one link which is not connected to any
driven qubit.

Figure S3:  a) Decomposable driving pattern for a chain: One can implement any
single qubit gates on the driving qubit (yellow) of the three-qubit block and
any two-qubit gate between the two driving qubits of the four-qubit block. b)
Decomposable driving pattern for a square array: single qubit gates can be
implemented on the driving qubit of the five-qubit blocks. Two-qubit gates are
more complicated in this case due to the condition for robust control (see
text), requiring a thirteen-qubit block. A two-qubit gate can be implemented
on any two driven qubits in the block.