Optimizing for periodicity: a model-independent approach to flux crosstalk calibration for superconducting circuits

Flux tunability is an important engineering resource for superconducting circuits. Large-scale quantum computers based on flux-tunable superconducting circuits face the problem of flux crosstalk, which needs to be accurately calibrated to realize high-fidelity quantum operations. Typical calibration methods either assume that circuit elements can be effectively decoupled and simple models can be applied, or require a large amount of data. Such methods become ineffective as the system size increases and circuit interactions become stronger. Here we propose a new method for calibrating flux crosstalk, which is independent of the underlying circuit model. Using the fundamental property that superconducting circuits respond periodically to external fluxes, crosstalk calibration of N flux channels can be treated as N independent optimization problems, with the objective functions being the periodicity of a measured signal depending on the compensation parameters. We demonstrate this method on a small-scale quantum annealing circuit based on superconducting flux qubits, achieving comparable accuracy with previous methods. We also show that the objective function usually has a nearly convex landscape, allowing efficient optimization.

Operating large-scale superconducting-circuits based quantum computers requires accurate calibration of the system parameters and control crosstalk [27,28].For flux control specifically, crosstalk arises due to the physical proximity between circuit elements and control lines, as well as reasons associated with the electromagnetic environment hosting the circuits, such as ground loops.For most large-scale superconducting circuits today, which are based on transmons, solving the calibration problem is often helped by the fact that transmons interact via the charge degree of freedom and the interaction strength is weak [29,30,31,32].For annealers based on flux qubits, calibration is more challenging because of the strong flux interaction between circuit elements, which makes it hard to directly measure the coupling between bias lines and flux loops [33,34,24].
In this work we introduce a new approach to calibrate crosstalk between flux biases for superconducting circuits.Relying on the fundamental property that superconducting circuits respond periodically to external flux biases [35,36], orthogonal control of flux biases can be obtained by optimizing for the periodic response of the circuit relative to each target bias coordinate, over the compensation parameters from all other bias lines.The periodicity analysis can be automated, allowing a closed loop optimization to be performed.Unlike conventional calibration methods which require either a simple model to describe the measurement signal, or high resolution scans, the periodicity optimization approach is circuit-model independent, and does not require high resolution scans.We demonstrate this approach on a small quantum annealing circuit, achieving a comparable accuracy with previous calibration methods on a subcircuit with three loops [33].
The paper is organized as follows.In Section 2 we introduce the notations used to describe the flux crosstalk problem and state the periodicity condition.In Section 3 we briefly review the principles behind earlier crosstalk calibration approach.In Section 4 we discuss the framework for crosstalk calibration based on maximizing periodicity, as well as the metric used to quantify periodicity.In Section 5 we discuss the experimental results of implementing this method on a quantum annealing circuit, followed by conclusions in Section 6.

Flux crosstalk and periodicity
The properties of superconducting circuits depend on the external flux biases of the superconducting loops.For a superconducting circuits with N flux bias loops, the external fluxes are usually controlled by N bias lines, which are mutually coupled to the flux loops.We denote the external flux bias in loop i, reduced by the flux quantum Φ 0 , as f i , and the corresponding bias line current as I i .The fluxes {f i } in all loops and currents {I i } on all bias lines can be written as vectors f and I respectively, and they are in general related by a linear transformation where M is the N × N mutual matrix describing the coupling between bias lines and flux loops, and f 0 is the vector of flux offsets arising from spurious sources.Often, and in particular in the context of our experiment, bias currents are controlled by voltage sources.For a more direct representation of the experiment, we will refer to the relation between fluxes and voltages, written as where V is the vector of voltages with each element controlling the corresponding element in I, and C = MR −1 with R a diagonal matrix consisting of the resistances between the voltage sources and the bias lines.From here onward we will work with voltage controls and the crosstalk matrix C.
Measurements on the superconducting qubits can be considered as a function mapping the flux biases to the signal R l , where l denotes a particular readout channel.Note that the number of readout channel is not restricted to the number of physical signal processing units; rather each channel corresponds to reading out the signal of an experiment, with a particular set of experimental parameters.The experiment could consist of one quadrature of a transmission measurement at a particular frequency or more complex experiments involving microwave excitations of the system.
Superconducting circuits respond periodically to external bias fluxes, with the period of one flux quantum.We denote the readout data as R, which is a vector with dimension of the number of readout channels.The periodicity condition can be formally stated as where m k is an integer.

Translation-based approach to flux crosstalk calibration
Most previously developed approaches to flux crosstalk calibration assume that one can identify a particular readout channel l that depends on the external flux in a single loop i, R l ({f k }) ≈ R l (f i ).This allows estimating the coupling coefficient from bias line j to loop i, C ij , by measuring the translation of R l as a function of V j .For this reason we denote such calibration method as the translation-based approach.When a simple model for R l (f i ) exists, the method becomes particularly effective as one only requires measurements at a few voltage bias values to extract the coupling parameters, and the model can be fitted to the data to obtain the coupling C ij .This is the case for many calibration methods used in tunable transmons, with the readout channel being the frequency of the transmon or its readout resonator [29,31], or the Ramsey phase shift [30,31,32].However, this method would only work if the circuit elements interact weakly, and each superconducting loop can be sufficiently decoupled from the other loops.We also note that the work presented in Ref. [37] uses an optimizationbased crosstalk calibration approach, however, this too relies on simplifying the full superconducting circuit to an effective description in terms of weakly coupled bosonic modes.
Recently, an iterative version of the translation-based approach was developed to tackle the issue of strong circuit interactions [33].The idea is that in the first iteration of the flux crosstalk calibration, one can assume the readout signal depends on only one flux bias, which changes due to the coupling from the bias lines alone and not due to inductive coupling to other loops.In each new iteration, the control coordinates to be swept become the estimated flux coordinates from the previous iteration, which allows one to gradually decouple the different superconducting loops and converge towards the true crosstalk.However, due to the strong circuit interactions, which are hard to model, this method often requires both high resolution scans and a large number of iterations to calibrate the crosstalk accurately.

Periodicity maximization approach
In this section we introduce the periodicity optimization approach to crosstalk calibration.We will first discuss the framework used to treat the calibration task as an optimization problem.Then we will discuss the measurement and analysis required to quantify periodicity.

Crosstalk calibration as an optimization problem
The task of crosstalk calibration is to obtain estimates of the coupling matrix C and independent control of the external flux biases.This is equivalent to finding N independent control coordinates, such that the circuit responds periodically to changes in each of them.To do this, we break the calibration task into N independent optimization problems, as described below.
We start by introducing initial estimates of the crosstalk and flux offsets, given by C init and f init 0 .Introducing them makes it convenient to discuss the optimization with or without prior knowledge on the same footing.When no prior knowledge is available, the initial estimates are identity and zeros for the crosstalk matrix and flux offsets respectively.The initial estimates allow us to define the initial control coordinates f init , The initial control coordinates f init are related to the actual fluxes f via the residual crosstalk and flux offsets, where where the compensation matrix O ′ has elements O ′ jk with There are N − 1 non-zero elements in the matrix O ′ , denoted as {Ω ji }.These are the compensation parameters to be optimized when calibrating the ith control coordinate.The objective of the optimization problem is to maximize the periodicity of the measured signal when varying the ith coordinate of the trial flux, f ′ i .This can be done by performing measurements varying f ′ i , and quantifying the periodicity using the metric discussed in Sec.4.2.A schematic for one iteration of the optimization step is shown in figure 1(a).
The maximum periodicity of the signal is achieved when compensation parameters satisfy specific relations relative to the residual crosstalk C res .To see this, consider the relation between the trial control fluxes and the actual fluxes, which follows from Eq. (4,5,6), where It can be seen that when the following condition is satisfied one has where f i , f ′ i , f ′ 0,i and C res ij are elements of f , f ′ , f ′ 0 and C res respectively.The relations between the actual fluxes f and trial control fluxes f ′ given by Eqs.(10,11) indicate that when the ith control flux f ′ i is being varied, only the ith actual flux f i changes.In other words, the residual crosstalk from the ith control coordinate to other coordinates l = i is completely cancelled out by setting the compensation parameters {Ω ji } satisfying Eq. 9. Since the circuit response is periodic to each flux f i with period 1, the circuit also responds periodically with respect to f ′ i , with period [(C res ) −1 ] ii .Hence, optimizing the periodicity for the ith coordinate gives the optimized compensation parameters approximately satisfying Eq. 9, and they are denoted as Ω ′ ji .After completing the optimization for all N control coordinates, we obtain N (N − 1) optimized compensation parameters, and another N parameters corresponding to the periods of the N control coordinates.Using Eq. 9, estimates for residual crosstalk matrix can be obtained, which we denote as C res′ .

Quantifying periodicity
The objective function for the optimization is the periodicity of the readout data with respect to f ′ i .To measure the periodicity, the readout data is collected sweeping a large enough range of f ′ i to cover a few periods, while keeping {f ′ j =i } fixed.The readout data is denoted as R l (f ′ i,s ), where s = 1, 2, . . ., m goes over the values of f ′ i taken during the sweep and m is the total number of f ′ i steps.Readout data from different channels is first normalized, by applying the operation For each control coordinate i For each loop i, the optimization parameters are elements of a trial compensation matrix O ′ , which defines the trial flux coordinates f ′ .Then the measurement is done by sweeping f ′ i and the periodicity of the measurement signal is determined.If the periodicity is high, the compensation parameters give good estimates of the ratio between the crosstalk matrix elements, otherwise, the compensation parameters are updated and the optimization is repeated until periodicity is high.(b) Schematic of the subcircuit of the device measured, with the tunable flux qubit on the left (purple), the quantum flux parametron (QFP) in the middle (yellow), and the tunable resonator on the right (grey).In addition, the qubit and the QFP are each coupled to a fixed-frequency resonator (grey).All resonators are coupled to a joint feedline (red).External flux biases in the loops are controlled via the on-chip bias lines (blue).

Measurement
where R l is the average of the readout data from channel l over all values of f ′ i taken.The periodicity can be quantified by first computing the correlation of the signal and its own with a translation of t steps along the f ′ i coordinates.Defining the translated signal as the correlation is where δ is the step size of the f ′ i sweep and t is an integer for the translation considered.The R l , R l,t refer to averages of the readout data over S for a particular readout channel l.From the definition of correlation, we have the range of ρ ∈ [−1, 1], with 1 for perfect correlation, −1 for perfect anti-correlation, and 0 for no correlation.
The correlation for a periodic signal is largest when the translation is an integer multiple of the period.However, since the period of the signal is in general not commensurate with the step size δ, we fit the following function around the maximum of ρ i  where τ can take non-integer values and τ max i corresponds to the period of the signal.The fit parameter ρ max i could be identified with the periodicity of the measurement signal.However, to be more precise, we choose to do another measurement where the sweep range is shifted from the original measurement by . The correlation between this signal and the original one is then computed and denoted as P , which is the objective function used in the optimization.

Experimental implementation
We implement the optimization procedure outlined above on subcircuits of a small prototype coherent quantum annealer.The device consists of two coupled tunable capacitively-shunted flux qubits [38], fabricated using a three-stack process in Lincoln Laboratory [39].Each qubit is coupled to a quantum flux parametron (QFP), which is in turn coupled to a flux-tunable resonator for readout.The QFP acts as an amplifier for the flux signal of the qubit, hence ensuring high-fidelity readout in the qubit flux basis, which is critical for quantum annealing applications [40].In addition, each qubit and QFP is inductively coupled to a fixed-frequency resonator to assist crosstalk calibration.A schematic of one qubit unit cell consisting of the qubit, the QFP and the tunable resonator is shown in Fig. 1(b).The full two-qubit system, including its readout circuits have been calibrated using the iterative translation-based method [33].The procedure is briefly reviewed in Appendix A and the crosstalk matrix obtained via this method is denoted as the reference crosstalk matrix, C ref .

Optimization of a Subcircuit with Three Flux Biases
For a proof-of-principle demonstration of the periodicity optimization approach, we start with a subcircuit consisting of just the QFP and the tunable resonator.The subcircuit has strong coupling due to the large persistent current of the QFP and the resonator, which makes it very time-consuming to calibrate using the translation-based method.The three flux biases in the subcircuit are denoted as QFPZ, QFPX and TR.The optimization starts with the initial crosstalk C init = C ref .This allows setting the qubits and couplers outside the subcircuit in a flux bias such that they are decoupled from the measured subcircuit.Using the reference crosstalk also allows systematic investigation of the performance of the optimization relative to particular initial conditions and bounds on the trial compensation parameters.We have also demonstrated the optimization starting with C init given by a single iteration of the translation-based method, which is discussed in Appendix B.
For the readout channel, we choose to measure transmission through the feedline that is coupled to the resonators.Both the fixed-frequency and tunable resonators are measured, each at six different readout frequencies.The readout frequencies are chosen to be around the bare resonator frequencies and the step size is around their resonance linewidth.The flux bias sweep range is chosen to cover about two periods and the step size is about 20 mΦ 0 .The other flux biases not being swept are set to values which avoid the flux-insensitive bias points of the QFP and tunable resonator.This is needed to avoid the tunable resonator and the QFP coincidentally being in fluxinsensitive spots, which would make the measurement signal insensitive to crosstalk.Such settings can be achieved without accurate initial estimates of the crosstalk or flux offsets.
As examples for the measurement and analysis at a single step in the optimization, we show the transmission measurement results at the start and the end of the optimization for the QFPZ control periodicity in figure 2(a,b,d,e).It is clear that the measurement signal is more periodic after the optimization.This is also reflected in the maximum correlations with respect to translations of the signal, which are shown in figure 2(c,f).
We use primarily an optimization algorithm based on Bayesian optimization [41], which is a global optimizer suited for black-box optimization with objective functions which are expensive to evaluate.The algorithm uses a Gaussian process to approximate the objective function, which is called the prior.At each step, the optimizer samples the distribution at a new point in the parameter space, which is probabilistically chosen according to the prior to improve upon the existing samples while minimizing the uncertainties of the prior [42].The Gaussian process is then updated according to the Bayesian inference rule, and is used as the prior for the next iteration.The compensation parameters Ω ji 's are bounded to within [−0.2, 0.2], and the optimization is initialized with evaluations at 20 random points in the parameter space.The bounds correspond to typical levels of crosstalk in large-scale devices [33].We defer to Sec. 5.2 for a discussion of how the bounds and initial conditions could affect the optimization.In figure 3(a), the trial compensation parameters and the periodicity is plotted versus the optimization step.It can be seen that the optimum parameters have been found after about 40 iterations.There is a drop in the periodicity near step 50, even when the compensation parameters remain mostly unchanged.This is due to the presence of hysteresis in the QFP (see Appendix C for further discussion).In figure 3  method [33].This shows that the crosstalk matrix obtained by the optimization method is comparable to the crosstalk matrix obtained by the iterative calibration method.We also note that the differences are much smaller than the flux sweep step size, which shows that the method does not require high resolution scans to be accurate.As a result of this, the optimization-based measurements required fewer data as compared to one iteration of the translation-based method.
Using the same optimizer setting but starting the optimization with estimated crosstalk from one iteration of the translation-based method, the estimated crosstalk obtained converged with similar level of accuracy, as compared to starting the optimization with the reference crosstalk matrix, obtained from multiple iterations of translation-based method.The result is presented in Appendix B.

Optimization Landscape
After demonstrating that the optimizations converge with high accuracy to the expected compensation parameters, we examine the structure of the optimization problem.First we looked at how periodicity changes as the compensation parameters deviate from the optimized values.We define the distance from the optimized compensation parameters as and plot the periodicities measured during the optimization versus Ω i in figure 4.
It can be seen that for all of the loops measured, when Ω i 0.001, the periodicity function plateaus at about 0.99, This suggests given the current set of readout channels, the optimized compensation parameters would allow us to control each bias coordinate to 1 mΦ 0 accuracy over one flux quantum range.The sharp peak for the QFPX loop is likely due to hysteresis of the QFP, which can be avoided by choosing a different set of independent flux control coordinates (see Appendix C).
When Ω i 0.1, the periodicity P ≈ 0. This means that the sampled trial compensation parameters are only informative when they satisfy Ω i 0.1.Hence, the optimization method is likely only efficient when initial crosstalk is known to within 10% accuracy, relative to the diagonal coupling elements.Various sources of estimation could provide such accuracy, such as one single iteration of the translationbased calibration method [33], measurement on different copies of the same device, or potentially careful electromagnetic simulation of the device.
We further characterize the landscape of the periodicity function by directly measuring it.This is done by first updating the initial crosstalk matrix with the optimized parameters, according to and then doing measurement in the updated f init coordinates.For each loop, the periodicity is measured sweeping one trial compensation parameter, while keeping the other at zero.This measured periodicities are plotted in figure 5(a).It can be seen that the periodicity is mostly a smooth function of the compensation parameters with a single maximum.There are two features outstanding.First, the periodicity relative to the QFPX loop has rugged landscape.This is likely due to the QFP becoming hysteric and not responding to flux bias variations fast enough compared to the experiment time.The hysteresis is caused by the discontinuous change in the ground state wavefunction of the QFP at the flux bias symmetry points.This can be systematically avoided by choosing a different set of linearly independent flux bias coordinates, along which the ground state wavefunction changes smoothly (see Appendix C for more detailed discussion).Second, the periodicity maxima for the compensations to TR loop are slightly deviated from zero.The reason for this still requires further investigation.One possibility could be that the periodicity function, under the measurement setting used, is sensitive to drifts in flux offsets, which could occur between the optimization measurement and the landscape measurement.The offsets therefore need to be kept track of in future implementations of the optimization, otherwise the accuracy of the crosstalk calibration based on periodicity optimization could be limited.The periodicity along the TR flux bias is also measured sweeping a two-dimensional grid of values for the trial compensations Ω QFPZ,TR , Ω QFPX,TR , over the range of [−0.1, 0.1] .The result is plotted in figure 5(b).It confirms that the periodicity is a smooth function over the entire range, and has a single maximum at around (0, 0).Such characteristics of the objective function means the optimization problem is likely convex in general.This opens the possibilities of using optimization algorithms that approximates and make use of the local gradients [43,44,45].We successfully implement one such optimization method, called simultaneous perturbation stochastic approximation (SPSA) [43] and the result is discussed in Appendix D.

Robustness of the periodicity metric to different measurement ranges
It is desirable that the periodicity optimization method is robust for a broad range of measurement settings, in particular when limited measurement data is available.In this section, we evaluate the effectiveness of the periodicity metric if the transmission measurement is done with a smaller range of flux bias or fewer frequency values than the measurements presented above.
To obtain the data with smaller flux bias ranges, we take the measurement done at every step of the optimization in section 5.1, which has a flux bias range of ∼ 2 Φ 0 , truncate the data to a specified smaller range, and then perform the periodicity analysis on the truncated data.Similar to figure 4, we plot in figure 6(a) the periodicity versus the distance from the trial compensation parameters to the optimum compensation parameters, for data truncated to approximately [1.2, 1.4, 1.6, 1.8] Φ 0 flux bias range respectively.It can be seen that with decreasing flux bias ranges, the maximum periodicity P reached is lower.This is because the overall contrast between the signal at a particular trial flux value R l (f ′ i,s ) and the average signal R l for a given readout channel l is decreased, reducing the value of the correlation function as given by Equation 14.Specifically, for bias range smaller or equal to ∼ 1.4 Φ 0 , the periodicity metric could not distinguish changes in the compensation parameters on the same order of accuracy as with the full data.However, as long as the bias range is larger than or equal to ∼ 1.6 Φ 0 , the cost landscape remains qualitatively the same as the full data, with the maxima at |Ω QFPX | 0.003.The analysis shows that the periodicity maximization does not require multiple periods of bias range, which allows good flexibility in designing the circuit.
To evaluate the effectiveness of the periodicity metric with fewer frequency values, we again take the data presented in Sec.5.1 and truncate it to a specified smaller range.It is found that reducing the number of frequency values measured to one per resonator would not significantly alter the cost landscape in the majority of cases.In figure 6(b) we show the periodicity versus the distance between the trial compensation parameters to the optimum compensation parameters, evaluated using only one particular frequency for each resonator measurement.It can be seen that when the first frequency value is used, the periodicity achieved is lower, and the landscape becomes noisy near Ω QFPX 0.01.In contrast, when the third or the fifth frequency value is used, a higher periodicity could be achieved, at a distance |Ω QFPX | 0.003.‡ The periodicity analyzed using the full data corresponds to an "averaged" value of the periodicity obtained using different frequencies, in line with the definition of the correlation function given in equation 14.This suggests that the periodicity maximization is not particularly sensitive to the readout frequencies chosen for the measurement, and likely does not require prior fine-tuning.

Optimization of a Subcircuit with Five Flux Biases
To understand the feasibility of the periodicity optimization on larger devices, we implement the procedure incorporating the qubit that is directly coupled to the QFP.The qubit loops are denoted as QZ and QX.In figure 7(a) we show the four compensation parameters and the periodicities with respect to QFPZ loop versus the optimization steps.It is noted that an increased number of initial evaluations, 50, is required for the optimization algorithm to approach the reference compensation parameters.The optimizations for other loops in the system did not approach the reference compensation parameters with the same optimizer hyperparameters.One possible explanation for the relative success of the QFPZ loop periodicity ‡ Using the fourth frequency value yields similar results to using the first, using the zeroth or the second yield similar results to using the third or fifth.They are not shown in the figure for visual clarity.(b) 1 minus the periodicity, analyzed with data at a chosen frequency for each resonator, versus the distance from the trial compensation parameters to the optimized compensation parameters.In both panels, the control flux coordinate to be optimized is QFPX, and the grey line corresponds to the analysis using the full data (exactly the same as figure 5 maximization, compared to the other loops, is that the QFPZ loop is special, both because of its large persistent current and it being directly coupled to most other loops in the subcircuit (except QX).The effectiveness of the optimization in larger circuits can potentially be improved by exploring different readout channels and optimization algorithms, which we didn't pursue in this proof-of-principle work.

Discussion and conclusion
There are a few points worth noting when comparing the iterative translationbased approach and the optimization-based approach.In terms of the number of measurement done at different bias settings, the optimization-based method requires about 3 times less measurements than the translation-based approach to achieve comparable accuracy.The majority of the saving comes from the fact that in the translation-based method, the QFP calibration requires a few high-density 2D scans of its X and Z loop control, while the optimization-based method only requires 1D scans.The translation-based approach cannot use 1D scans to calibrate the crosstalk between the QFP X and Z loop, as the response of the circuit to one of them is highly dependent on the other.This could be considered as an extreme case of two strongly coupled loops (though technically they belong to the same circuit element).Therefore, we expect the optimization-based method could reduce the number of measurements required for calibrating strongly-interacting circuits or multi-loop components, such as the Josephson phase-slip qubit [46].
On the other hand, we note that the translation-based approach is more robust than the optimization-based method.The latter requires prior knowledge of the crosstalk matrix to bound the parameter space of the optimization.When the number of optimization parameters increases, the size of the parameter space over which the periodicity metric is insensitive to changes in trial parameters.This leads to an increase in the number of initial evaluations required for the Bayesian algorithm to learn the cost landscape, as we see in section 5.4.Therefore, we expect that more initial knowledge of the crosstalk matrix is required in order to implement the current Bayesian optimization-based periodicity maximization to larger devices.In summary, we introduced a flux crosstalk calibration method based on the fundamental periodicity of superconducting circuits, without relying on any underlying circuit model.The method is successfully demonstrated on a coupled QFP-tunableresonator system, with an accuracy that is comparable to the previously developed iterative translation-based calibration method.Although the current implementation of the optimization approach is limited when used to calibrate devices with larger number of loops, it can already be utilized as a subroutine for calibrating parts of a larger system.Such hybrid calibration strategy is particularly useful for strongly interacting systems such as the quantum annealing circuits investigated here, where translation-based approach alone requires larger number of iterations and highresolution measurements.
The landscape measurement shows that the problem is nearly convex within some bounds on the optimization parameters.This points to exploring other optimization algorithms, such as momentum based optimizations [44,45] to speed up the convergence, which is crucial for extending the optimization-based calibration to larger devices.Another attractive future direction could be adaptive measurements, where different experiment parameters can be used to give different optimization landscapes.For example, an optimization landscape with a broad maximum could afford large tolerance to the initial guess of the crosstalk matrix, while optimization landscape with a narrow maximum could lead to higher accuracy for the optimized results.
The iterative approach is applied to the two-qubit circuit described in the main text.In figure A1 we show how the measured crosstalk and flux offsets converge towards identity and zeros.In figure A2 we show the final crosstalk matrix from the iterative procedure for the qubit, QFP and tunable resonator, on which the periodicity optimization approach is implemented.

Appendix B. Optimization initialized with single iteration of translation-based calibration
In this section we describe the results obtained by performing the periodicity optimization, starting from the estimated crosstalk of one iteration of the translationbased approach.In figure B1(a) we show the estimated crosstalk matrix obtained by a single iteration.This can be compared with the reference matrix elements plotted in figure A2.It can be seen that after a single iteration, the estimated crosstalk still deviates from the reference values, by as large as ∼ 10%.
In figure B1(b), we plot the deviation between the reference crosstalk matrix and the estimated crosstalk matrix after the optimization.Most of the deviation is about or less than 3 × 10 −3 .This is comparable to the accuracy of the results discussed in the main text, which starts the optimization directly from C ref .The only exception is the QFPZ diagonal element, which corresponds to its period.This is probably due to the hysteresis of the QFP, which can be resolved by repeating the QFPZ periodicity measurement at a different QFPX biasing point.versus the iteration number.The orange bar is the median, the black box corresponds to the lower and upper quartiles, the segments contain the 5th to 95th percentiles of the data and the dots are outliers.not tunnel to the persistent current state with lower energy.To solve this problem, we can replace the Z and X biases with By keeping fZ approximately fixed when sweeping fX , one can avoid switching the sign of the bias between the two persistent current states, and hence avoiding the need for tunneling to occur for the QFP to respond to changes in flux biases.Using the new flux bias coordinates, the periodicity along fX increases to 2 (keeping fZ fixed).

Appendix D. Optimization with SPSA
In this section we discuss the optimization results using an alternative optimizer called the Simultaneous Perturbation Stochastic Approximation (SPSA) [43].This algorithm approximates the gradient of the objective function by measuring the finite difference The difference is about twice as large as compared to the results obtained using Bayesian optimization.We expect the results to improve by using better hyper-parameters for the optimization, such as the magnitude of the perturbation, which would likely remove the parameter oscillations near the end of the optimization.

Figure 1 .
Figure 1.(Color Online) (a)Schematic representation of the optimization step.For each loop i, the optimization parameters are elements of a trial compensation matrix O ′ , which defines the trial flux coordinates f ′ .Then the measurement is done by sweeping f ′ i and the periodicity of the measurement signal is determined.If the periodicity is high, the compensation parameters give good estimates of the ratio between the crosstalk matrix elements, otherwise, the compensation parameters are updated and the optimization is repeated until periodicity is high.(b) Schematic of the subcircuit of the device measured, with the tunable flux qubit on the left (purple), the quantum flux parametron (QFP) in the middle (yellow), and the tunable resonator on the right (grey).In addition, the qubit and the QFP are each coupled to a fixed-frequency resonator (grey).All resonators are coupled to a joint feedline (red).External flux biases in the loops are controlled via the on-chip bias lines (blue).

Figure 2 .
Figure 2. (Color online) Transmission versus probing frequency ωp and trial flux coordinate f ′ QFPZ , through the fixed-frequency (a,d) and tunable (b,e) resonator, at the first (top) and last (bottom) step in the optimization.The plots on panels (c, f) show the corresponding correlation versus translation, with inset showing the absolute value linear fit around the maxima.
(b), the landscape of the objective function, predicted by the final Gaussian process model is shown, together with markers for the parameter values sampled during the optimization.The minimum is at around (Ω QFPZ QFPX , Ω TR QFPX ) = (0, 0), which is the expected optimum compensation parameter given C init = C ref ≈ C and hence C res ≈ I.In figure3(c), we show the difference between elements of C res′ with the identity matrix.The magnitudes of the elements are all below 3 × 10 −3 , which is about the error of the iterative

Figure 3 .
Figure 3. (Color online) (a) Trial compensation parameters (left axis), Ω QFPZ QFPX (dashed line), Ω QFPZ QFPX (solid line) and periodicity P (red dots, right axis) versus optimization step.(b) Gaussian process model of periodicity versus the compensation parameters.The cross markers correspond to parameters sampled by the optimizer and the gray scale of the markers indicates the sequence at which they are sampled, with darker color markers being sampled later.(c)Difference between the estimated residual crosstalk matrix C res′ and the identity matrix.

Figure 4 .
Figure 4. 1 minus the periodicity P versus the distance from the trial compensation parameters to the optimized compensation parameters for the three loops, QFPZ (left), QFPX (center), TR (right)

Figure 5 .
Figure 5. (Color online) (a) Measured (blue dots) periodicity versus deviation of the compensation parameters from the optimized value for each of the six offdiagonal parameters in the 3 × 3 matrix, with a quadratic fit (red curve) on top.(b) Measured periodicity along the f ′ TR bias versus deviation of the compensation parameters Ω QFPZ,TR , Ω QFPX,TR from their optimized values.

Figure 6 .
Figure 6.(Color online) (a) 1 minus the periodicity, analyzed with data truncated to a smaller bias range, versus the distance from the trial compensation parameters to the optimized compensation parameters.(b) 1 minus the periodicity, analyzed with data at a chosen frequency for each resonator, versus the distance from the trial compensation parameters to the optimized compensation parameters.In both panels, the control flux coordinate to be optimized is QFPX, and the grey line corresponds to the analysis using the full data (exactly the same as figure5(b).The spikes near |Ω QFPX | < 0.001 are due to the hysteretic behavior of the QFP.Repeating the analysis for the other two flux coordinates QFPZ and TR yields similar plots.
Figure 6.(Color online) (a) 1 minus the periodicity, analyzed with data truncated to a smaller bias range, versus the distance from the trial compensation parameters to the optimized compensation parameters.(b) 1 minus the periodicity, analyzed with data at a chosen frequency for each resonator, versus the distance from the trial compensation parameters to the optimized compensation parameters.In both panels, the control flux coordinate to be optimized is QFPX, and the grey line corresponds to the analysis using the full data (exactly the same as figure5(b).The spikes near |Ω QFPX | < 0.001 are due to the hysteretic behavior of the QFP.Repeating the analysis for the other two flux coordinates QFPZ and TR yields similar plots.

Figure 7 .
Figure 7. (Color online) Optimization including the qubit loops.(a) Trial compensation parameters (left axis) from QZ, QX, QFX and TR to the QFPZ loop and periodicity P (right axis) versus optimization steps.(b) 1 minus the periodicity P versus the distance from the trial compensation parameters to the optimized compensation parameters for the QFPZ loop.

Figure A1 .
Figure A1.For the full two-qubit device, statistical box plots of the diagonal (top), off-diagonal (middle) coupling coefficients in C (n)′ , and the flux offsets (bottom) in f (n)′ 0

Figure A2 .
Figure A2.Final crosstalk matrix C ref obtained from the iterative procedure for the qubit, QFP and tunable resonator subcircuit, on which the periodicity optimization is applied.

3 )
Figure B1.(Color online) Top: C init as given by the first iteration of the translation-based calibration method.Bottom: The deviation between the optimized crosstalk matrix C ′ and the reference crosstalk matrix C ref , defined as C ref (C ′ ) −1 − I.The optimized crosstalk matrix C ′ = C res′ C init is obtained by starting the optimization with the matrix given in the top panel.