Scalable and robust quantum computing on qubit arrays with fixed coupling

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 [1]. High fidelity two-qubit gates were also demonstrated for trapped ions [2], neutral atoms in optical tweezers [3], and spin qubits in silicon [4,5] and GaAs [6]. 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 [1,4,7] leading to overheads in fabrication and wiring. For solving real-world problems, a quantum computer needs a large number of qubits [8], 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 [9,10]. For large arrays of qubits these series become increasingly complex, which is a bottleneck for scaling up [10]. 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 evolu-tion 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.
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 FIG. 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.
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 [11,12]. When a subset of qubits is driven by external fields, the system's Hamiltonian is where ν 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 L the driven subset. Typically, the qubit's transition energy, ω j , is much larger than the interaction, J jk , and hence |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) = where Ω j (t) = Ω x j (t) cos(δ j t) + Ω y j (t) sin(δ j t), Here, Ω x,y j (t) ≡ d j E x,y j (t) is the Rabi frequency and δ j ≡ ν j − ω 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) = l H l (t), where all the H l (t) are mutually commuting. The unitary evolution after a time duration T is then and 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.
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) = j∈L H j (t) where L is the driven subarray and where 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 σ x and σ y do not commute with σ 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 σ z j σ 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 σ z j σ z k T , which cannot be cancelled due to the absence of control.

Applying gates using optimal control
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, Ω x,y j (t) where j ∈ driven subset, to implement a unitary operation of the type U C ⊗ I B , where U C is a unitary acting on the driven subset, and I B the identity matrix acting on the undriven subset at the boundary. The net effect is that the gate U C is applied to the driven subset while the rest remains unchanged. If the driven subset has one (two) qubit, then U 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 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 UC ⊗ IB where UC is a target unitary on the central driven subset, and IB ≡ I ⊗ I ⊗ · · · ⊗ 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. small time step is 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.
it is obvious that the third term allows partial control of the undriven qubits, labelled by k, through shaping Ω x,y j (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 C ⊗ I B with maximum fidelity, F = 1, up to numerical precision, where U C is the Hadamard, π/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 C is the CNOT gate and the direct twoqubit identity gate, I ⊗ 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 [13]. 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.
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 Hamil- The final wavefunction can be expanded in the computational basis, Since H 0 consists of exclusively σ z terms, |j 1 , ..., j N is an eigenstate of H 0 and the evolution under H 0 introduces only addition phases in the coefficients c j1,...,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 [14] and HHS algorithm for solving linear systems [15], the additional phase will not change the result.

Robustness against parameter uncertainty
We now describe the optimal control algorithm for maximising the fidelity of the unitary U C ⊗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 ∆t. In each time bin the field amplitudes are kept constant. The set of the Rabi frequencies form the control vector, c = {Ω µ jn : 1 ≤ n ≤ M ; µ = x, y; j ∈ C}, where C is the driven subset. The unitary evolution U G (T ) of the star graph G is then a function of c. Each qubit-qubit coupling, J jk , and detuning, δ j , is allowed to vary independently in the uncertainty intervals [J − ∆J/2,J + ∆J/2], and [δ − ∆δ/2,δ + ∆δ/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 Ω x,y j (t) in H G by α j Ω x,y j (t) where α j is a dimensionless parameter that varies in the interval [1 − ∆α/2, 1 + ∆α/2]. Now the unitary U G (T ) also depends on J jk , α j , and δ 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, ω j (see Method). Thus, we assume that the driving fields are tuned to resonance,δ = 0, with residual detuning uncertainties much smaller than the average qubit-qubit interaction, ∆δ j /J = 0.1%, which is typical for superconducting qubits [16]. 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 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 v, where V is the hypercube containing the possible values of v, and F (c, v) the gate fidelity based on the trace distance where D is the dimension of the Hilbert space of G, I B the identity matrix of the subset of undriven qubits on the boundary, and U C the target unitary that we want to apply on the central driven subset. In numerical computation one chooses a set of sampling points, v i , in V, and find the minimum fidelity in this set. We found that when the uncertainties are all smaller than 5% and F (c, v) is larger than 99% its minimum over v always lies at one of the extreme points of V, i.e., one of the corners of the hypercube. Thus, we can redefine F(c) ≡ min vi∈X F (c, v i ) where X is the discreet set of the extreme points of V. This drastically reduces the number of sampling points where F (c, v) has to be computed. There are 2 nu extreme points in a hypercube of n u uncertain parameters. For example, the 6-qubit block in Fig. 2b has 9 uncertain parameters, five Js, two αs, and two δs, so we need to compute 2 9 ≡ 512 values of F (c, v i ) for any given c. We develop a numerical computation for gradientbased robust optimal control that can handle systems with up to 12 qubits and multiple uncertain parameters. Starting with a random initial guess for c, we use gradients to identify a step δc to maximize min i ∇ c F (c, v i ).δc, the first order increment, so that F (c, v i ) is increased for all v i . In this way F(c), the minimum fidelity over v i , ∆J/J ∆α ∆δ/J H ⊗ IB π/8 ⊗ IB I ⊗ IB CX⊗IB I2 ⊗ IB 0% 0% 0% 10 10 10 10 10 0.1% 0.1% 0.1% 5.6 (4.5) 5.5 (4.6) 5. 5 Fmax), which is the number of nines in Fmax. 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 10J. I and I2 are the direct identity gates for one and two qubits, respectively. The gate duration is T = 2π/J, divided into M = 100 time bins. The figures inside the bracket are results obtained with non-robust optimization.
can be increased to a value very close to one. In practice we find that one can also raise F(c) by maximizing the average fidelity, with gradient ascent, which in some cases works faster than trying to increase F (c, v i ) for all v i . For the n-th time bin the time evolution is calculated by the mid-point rule e −iH G (tn−∆t/2)∆t + O(∆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 [17].
For the four-qubit block of Fig. 2a, we derive optimal pulses to realize U C ⊗ I B where U C is the Hadamard and π/8 gates [13]. And for the six-qubit block of Fig. 2b we want U 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 [13], 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 I 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 I is given in the Supplemental Materials. We choose T = 2π/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 =J, α j = 1 and δ j = 0 for all j, then we use the obtained optimal control, c ideal , to calculate the minimum fidelity in the hypercube, F = min v∈V F (c ideal , v). The results are shown in the parentheses of Table I. robust optimization improves the fidelities by two to three nines when the uncertanties are significantly large (1% and 5%).

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 [11,18]. There is another interesting scheme based on the inductive longitudinal coupling of the flux qubits with a common bus cavity [19,20], which can be scaled up to 2D arrays [21]. In addition, a cross-Kerr ZZ interaction has been recently demonstrated for transmon qubits using a flux-tunable coupler [12]. 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 [22]. The ZZ coupling is also the natural interaction in a nuclear  Table I. magnetic resonance quantum computer [9], for which sophisticated pulse-shaping is available [23], 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 qubits is enhanced by a factor of approximately N in the worst case, giving rise to a lifetime of T 2 = T 2 /N . The fidelity of a multiqubit unitary on the block is then bounded by F ≤ 1 − T /T 2 ≡ 1−N T /T 2 . Therefore, to achieve a fidelity F for the four-qubit and six-qubit blocks in Table I the pulse duration needs to be shorter than T 2 (1 − F )/N where N = 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 onedimensional 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.

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, G, at the midpoint of the n−th interval from t n−1 to t n is where α j is the dimensionless factor introduced to model the uncertainty in the Rabi frequencies (see the Main Text), C is the driven subset in the center of the graph and J jk = 0 only for nearest neighbors, and Ω jn = Ω x jn cos(δ j (t n − ∆t/2)) + Ω y jn sin(δ j (t n − ∆t/2)), Ω jn = Ω y jn cos(δ j (t n − ∆t/2)) − Ω x jn sin(δ j (t n − ∆t/2)), where Ω 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 2M N C × 1 control vector, c, where N C is the number of qubits in the driving subset. The unitary evolution from t n−1 to t n , U n = e −iH G,n ∆t + O(∆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 [17,24], defined by Then the fidelity is where N G is the number of qubits in the star graph.
One can show that the derivative of U n ≡ e −iH G,n ∆t with respect to Ω µ jn is [17] ∂U n ∂Ω µ jn = −i∆tK µ jn + where H G,n , K µ j,n is a commutator. It follows that the 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, F(c), defined in Eq. 5. The first is based on sequential convex programming [25]. We start with initial guesses c 0 and u 0 for the control vector and the upper limit (trust region) of the step, respectively. Then a step |δc| < u 0 is found to maximize min i ∇ c F (c, v i ).δc, i.e., to maximize the minimum first-order increment over v i . This ensures all the fidelities at the extreme points of the hypercube, 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 ∇ c F (c, v i ).δc is positive then we increase the upper bound 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, using a quasi-Newton method. Obviously this does not guarantee that the worst-case fidelity over 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 op-timizeF(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. Summing all the euqations in Eqs. (11) and (12) Finally, subtracting Eq. (10) from Eq. (13), the dispersive shifts are cancelled and one arrives at 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 [26].

Code Availability
The Matlab code for this work are available without restriction [26].

SUPPLEMENTAL MATERIALS
Pulse shapes for the gates in Table I 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.
FIG. 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) π/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.

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.
FIG. 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.