Tutorial: Calibration refinement in quantum annealing

Quantum annealing has emerged as a powerful platform for simulating and optimizing classical and quantum Ising models. Quantum annealers, like other quantum and/or analog computing devices, are susceptible to nonidealities including crosstalk, device variation, and environmental noise. Compensating for these effects through calibration refinement or"shimming"can significantly improve performance, but often relies on ad-hoc methods that exploit symmetries in both the problem being solved and the quantum annealer itself. In this tutorial we attempt to demystify these methods. We introduce methods for finding exploitable symmetries in Ising models, and discuss how to use these symmetries to suppress unwanted bias. We work through several examples of increasing complexity, and provide complete Python code. We include automated methods for two important tasks: finding copies of small subgraphs in the qubit connectivity graph, and automatically finding symmetries of an Ising model via generalized graph automorphism. Code is available at https://github.com/dwavesystems/shimming-tutorial.

Quantum annealing (QA) [1,2] is a computing approach that physically realizes a system of Ising spins in a transverse magnetic field.A common application of QA is to find low-energy spin states of the Ising problem Hamiltonian Here, {σ z i } N i=1 ∈ {−1, 1} N is a set of Pauli z-operators, which can be thought of as a vector of classical ±1 Ising spins; h i denotes a longitudinal field (bias) on spin i, and J ij (used interchangeably with J i,j depending on context) denotes a coupling (quadratic interaction) between spins i and j.Minimizing H P is intractable, i.e., NP-hard [3].
QA adds to H P an initial driving Hamiltonian The ground state of H D , which is a uniform quantum superposition of all classical states, is easy to prepare.QA guides a time-dependent Hamiltonian H(s) from H D to H P , by linearly combining H D and H P as where s is a unitless annealing parameter ranging from 0 to 1. Unless stated, s is simply t/t a : time normalized by annealing time.The functions Γ(s) and J (s) define the annealing schedule: Γ(s) decreases toward 0 as a function of s, and J (s) increases as a function of s; Γ(0) J (0).An example is shown in Fig. 1.

B. CALIBRATION IMPERFECTIONS AND REFINEMENT
Quantum processing units (QPUs, in this case quantum annealers) are typically made available with a single onesize-fits-all calibration.Nonidealities in the calibration can arise from a number of sources.For example, small fluctuations in the magnetic environment can bias qubits in one direction or the other.So can crosstalk, in which a Hamiltonian term, e.g., a programmed coupler J ij , can cause an undesired perturbation in another Hamiltonian term corresponding to a physically nearby device, e.g. a bias field h i .
In short, no calibration is perfect.Oftentimes, in-depth studies of a single system (Ising model) or ensemble of systems (e.g., a set of realizations of a spin-glass model) can be improved by suppressing crosstalk and other nonidealities.This is achieved by "shimming": inferring statistical features of an ideal annealer, and tuning the Hamiltonian to produce these features.An ideal annealer, in this work, is defined simply as one that respects  The FM loop has a twofold-degenerate ground state (all spins up or all spins down) with no frustration; the frustrated loop has 2L ground states, each with one frustrated bond.When h = 0, all qubits trivially have zero average magnetization in an ideal annealer.
symmetries in the Hamiltonian-each qubit behaves identically, and each coupler behaves identically.Variations on the methods described herein have by now been used in many works [4][5][6][7][8][9][10][11].Often, when behavior of the system relies on precise maintenance of energy degeneracy between states, or energy splitting from the transverse field, results are highly sensitive to these tunings.Particularly for the simulation of exotic magnetic phases, calibration refinement is an essential ingredient of successful experiments.Yet, so far the discussion of these methods has mostly been relegated to supplementary materials.Our aim here is to provide an accessible guide that will encourage the use of these powerful but simple methods.
Specific visual demonstrations of the benefit of these methods "in the wild" include: • Frustrated 2D lattice, [4] Extended Data Fig. 7.

C. INFERRING STATISTICAL FEATURES: QUBIT AND COUPLER ORBITS
The approach described in this tutorial can be stated simply and generically: In theory, two observables of a QPU output are expected to be identical due to symmetries in the Ising model being studied.In experiment they can differ systematically.We tune Hamiltonian terms to reduce these differences.
In this work we only consider one-and two-spin observables-spin magnetizations and frustration probabilities-in part because they can be fine-tuned easily using the available programmable terms in the QPU.A call to the QPU typically results in a number of classical samples, which we set to 100 for all examples.From these samples we can compute a magnetization for each spin s i , and a frustration probability for each coupler J i,j ; f i,j is the observed probability of the coupler having a positive contribution to the energy in H P .This raises the first question: how do we identify observables that should be identical in expectation?The answer is: through symmetries of the Ising model under spin relabeling and gauge transformation1 .We understand and formalize these symmetries-and automate their detection-through graph isomorphisms (especially automorphisms) and generalizations thereof [12].The symmetries we find and exploit here are a subset of all possible symmetries.
These symmetries admit two types of equivalence relations on an Ising model H P : one on the qubits, and one on the couplers.We call the equivalence classes qubit orbits and coupler orbits, respectively.We use notation O(s i ) for a qubit orbit containing spin s i , and O(s i , s j ) for a coupler orbit containing coupler (s i , s j ).We define them as having the following properties guaranteed by symmetry in an ideal annealer: • All qubits in the same orbit have the same expected magnetization.
• All couplers in the same orbit have the same expected frustration probabilities.
Due to spin-flip symmetries, each qubit and coupler orbit can additionally have up to one nonempty orbit that is opposite.
• If qubit orbits O(s i ) and O(s j ) are opposite, then -If O(s i ) = −O(s j ), then h i = −h j and, in an ideal annealer, m i = −m j .
• If coupler orbits O(s i , s j ) and O(s k , s ) are opposite, then -J i,j = J k, and, in an ideal annealer, f i,j = f k, .
We will sometimes overload notation, conflating O(s i ) with O(i), and O(s i , s j ) with O(i, j).
Qubit and coupler orbits are related to, but not identical to, automorphism orbits of an auxiliary graph.In particular, qubit and coupler orbits are not unique: Putting each qubit and each coupler in a separate orbit is sufficient to meet the definition, but does not provide any useful information.We seek large orbits that satisfy the requirements.
Note that in the commonly arising situation where h i = 0 on all qubits, each qubit orbit is its own opposite, so all qubits have m i = 0.The analogous situation does not exist for couplers, because we do not consider symmetries between pairs of qubits with zero coupling between them.Two simple examples are shown in Fig. 2: a frustrated loop and an unfrustrated loop.In each case, all qubits are expected to have magnetization and all couplers are expected to have the same probability of frustration, but this is less obvious in the frustrated case than in the ferromagnetic case.
Having defined qubit and coupler orbits, we now consider how to find them.

Automorphisms of the signed Ising model
Let (h, J) denote an Ising model with fields h = {h i |v i ∈ V } and J = {J i,j |e i,j ∈ E}, with an underlying graph G = (V, E) with vertex and edge sets V and E. We construct a signed Ising model S(h, J) as follows: • For each spin v i ∈ V , S(h, J) has two spins v i and vi , with fields h i and −h i respectively.
• For each coupler e i,j = (v i , v j ) ∈ E, S(h, J) has four couplers: two couplers (v i , v j ) and (v i , vj ) with coupling J i,j , and two couplers (v i , v j ) and (v i , vj ) with coupling −J i,j .
Informally, we simply replace each spin with two: itself and its negation, and replace each coupler with four couplers with appropriate parity-based sign flipping.Fig. 3 shows an example of this construction applied to a four-spin Ising model.
Our aim is to find large qubit and coupler orbits for (h, J), and we will begin by finding the automorphism group of S(h, J), which can be considered as a vertex-and edge-labeled graph.The automorphism group of S(h, J) naturally generates one equivalence relation defining qubit orbits, and another equivalence relation defining coupler orbits (see Fig. 4)2 .Our orbits of S(h, J) immediately give us orbits of (h, J), constructed by simply discarding the qubits and couplers that do not exist in (h, J).
There is more usable information held in the orbits of S(h, J).First, we can combine coupler orbits of S(h, J) such that for each coupler e i,j ∈ E, (v i , v j ) and (v i , vj ) are in the same orbit, and (v i , v j ) and (v i , vj ) are in the same orbit.Second, we can then easily derive opposite orbits: These orbits are already very useful, but we can combine some to make even larger orbits.As demonstrated in the example in Fig. 3, in S(h, J) the couplers between pairs (v i , v j ) and (v i , vj ) are not necessarily automorphic.Ising model (h, J) Auxiliary Ising model S(h, J)  Orbits of (h, J) Orbits of S(h, J) Orbits of signed and original Ising model.By computing automorphism groups of the edge-and vertex-labeled graph of the signed Ising model (Fig. 3, right) we can construct orbits of qubits and couplers that should behave identically by symmetry in S(h, J) (left).Here, vertex and edge labels indicate orbits.By identifying equivalent orbits (e.g., coupler orbits 0 and 5) and reducing back to the original Ising model (h, J), we determine effective qubit and coupler orbits of (h, J), and their opposite relations (right).
Yet they are clearly equivalent under a flip of all spins.Thus we combine the coupler orbits containing these two couplers.Likewise, the same applies to couplers between pairs (v i , vj ) and (v i , v j ).This is all demonstrated in the accompanying code example0 1 orbits.py,and shown in Fig. 4. We now consider how to exploit orbits to improve performance in quantum annealers, building up a set of tools in the following worked examples.
For our first example of calibration refinement, we study the ferromagnetic loop in which each coupling J i = J i,i+1 is equal and each field h i is zero.In this case, by rotation, it is obvious that all qubits are in the same orbit and all couplers are in the same orbit.Furthermore, the orbit containing all qubits is its own opposite.Thus we will perform two refinements: First, we will balance each qubit at zero magnetization m i ≈ 0. Second, we will balance the couplings so that each coupler is frustrated with approximately equal probability.
Since the FM loop has no frustrated bonds in the ground state, the latter condition is only interesting if we sample excited states.To ensure abundant excitations, we study a reasonably long loop with weak couplings: L = 64 and J ij = −0.2.

A. FINDING MULTIPLE EMBEDDINGS OF A SMALL ISING MODEL
Code reference: embed loops.py.
The first task is to find a copy of the FM loop in the qubit connectivity graph A QPU of the QPU being used.This is an embedding-a mapping of spins of an Ising model to qubits in a QPU.In an Advantage processor, a 64-qubit loop can be embedded many times on disjoint sets of qubits, so we can run many copies in parallel for a richer and larger set of measurements.
To find these embeddings we use the Glasgow graph solver [14], which has been incorporated into the embedding finding module minorminer [15].To make the embedding search faster, we raster-scan across 2 × 2 blocks of unit cells in the QPU's Pegasus graph [16], then greedily construct a large set of non-intersecting embeddings.The file embed loops.pyprovides a code example that finds multiple disjoint copies of a 64-qubit loop in A QPU .

B. BALANCING QUBITS AT ZERO
Code reference: example1 1 fm loop balancing.py.
We will use simple parameters for the experiment, running 1 µs anneals and drawing 100 samples for each QPU call.We set auto scale=False to ensure that the QPU will not automatically magnify the energy scale.
In D-Wave's annealing QPUs, each qubit s i can be biased toward −1 or +1 in two ways: first, with a programmable longitudinal field h i as in Eq. 1; second, with a programmable flux-bias offset (FBO) Φ i [17,18].In the quantum annealing Hamiltonian (3), the bias conferred by the h i term is scaled by J (s), meaning that it changes as a function of s.The FBO Φ i , in contrast, confers a constant bias that is independent of s.We prefer to mitigate biases using FBOs, in part because they are programmed independently of h i .
We employ an iterative gradient descent method for minimizing |m i | with a step size α Φ .For a given iteration we consider the observed magnetization m i = s i .If m i < 0 we adjust the FBO to push s i toward +1; if m i > 0 we adjust the FBO to push s i toward −1.This is done by updating for each qubit after each iteration, where m is the average observed magnetization across all qubits.In this case, we can simply replace m with 0 since h i = 0 for all qubits.In Figure 6 we show the resulting FBOs for a single copy of the 64-qubit chain, as well as magnetization statistics.We show experiments for three choices of α Φ .One (flux 1e−4, in units of Φ 0 ) is too large, and creates oscillations in Φ i and m i .One (1e−6) is too small, and takes many iterations to converge.One (1e−5) is in between, and performs well.The choice of step size is a common concern in gradient descent applications, and we will consider automatic tuning of α Φ in a later section.For best results, we should at a minimum ensure: • The calibration refinement appears to have converged to the vicinity of a fixed point.
• The parameters do not oscillate wildly.Balancing qubits and couplers in a FM chain with flux-bias offsets and coupler adjustments.This experiment is similar to that shown in Fig. 6, but with αJ > 0 for the last 100 iterations.Couplers remain distributed about the average value of J = −0.2.
When seeking evidence that qubit bias is improved by the FBOs, we should not just look at qubit statistics over a single QPU call, since fluctuations can be large.Rather, we should look at the average magnetization of a qubit over multiple calls, which indicates systematic bias.The middle row of Fig. 6 shows the average magnetization of each qubit across the first and last ten iterations.For each step size, the shim results in a significant improvement in variation of m i from one qubit to another.However, the standard deviation among qubit magnetizations for individual iterations shows that the case α Φ = 1e−4 causes broad spreading of biases, so we need to be careful with our step sizes.
Having balanced qubits at zero with linear terms with an FBO shim, we now address homogenizing the spin-spin correlations on adjacent qubits, which by symmetry should be equal for all coupled pairs.The couplings J i,i+1 are all nominally −0.2; we will fine-tune the couplings in the vicinity of this value.This is similar to how we fine-tuned the FBOs, but with the added constraint that we do not change the average coupling.
For a given iteration we take the observed probability f i,i+1 of the coupler being frustrated: Let f denote the average frustration across all couplers in all disjoint embeddings of the chain-in general, we will compute f across all couplers in the union of a coupler's orbit and its opposite orbit.We then adjust couplings based on the residual frustration f i,i+1 − f : Fig. 7 shows data for the same experiment as Fig. 6, but with the "coupler shim" added, with α J = 0.001.To show the effect of the two shims, we run 100 iterations with α Φ = 0 and α J = 0, then 100 iterations with α Φ = 1e−5 and α J = 0, then 100 with α Φ = 1e−5 and α J = 0.001.In this particular case, the coupler shim is small but some systematic signals can be seen.We will show more impactful cases later in the tutorial.

Part III
Worked Example: Frustrated loop Code reference: example2*.py.Take the ferromagnetic loop considered in the previous example, and flip the sign of a single coupler J 1,2 .It is again obvious that all spins should have zero average magnetization, since there is no symmetry-breaking field (i.e., h i = 0 everywhere).Less obvious is the fact that we can have two coupler orbits: one containing all FM couplers, and one containing the AFM coupler, and they are opposite.Consequently, every coupler should be frustrated with equal probability in an ideal annealer.
We can derive this fact as follows.Flipping the sign of s 2 , and the sign of both couplers incident to it, is a gauge transformation and, as such, will not change the probability of any coupler being frustrated in an ideal annealer.The result of this gauge transformation is again a frustrated loop with a single AFM coupler J 2,3 ; note that this is equivalent to the original loop by a cyclic shift of qubit labels.From this we can infer that J 1,2 and J 2,3 should have the same frustration probability; repeating this argument tells us that all couplers should have the same frustration probability.
Coupler orbits of frustrated loops.The code example2 1 frustrated loop orbits constructs three disjoint frustrated loops and programmatically generates qubit and coupler orbits.All qubits are in the same orbit.There are two signed coupler orbits, O1 and O2, and in this example they form an opposite pair, meaning that a coupler in O1 and a coupler in O2 have opposite sign (J = −1 and J = 1 in this case) but equal probability of frustration in an ideal annealer.
For more complicated examples, we would prefer to find such statistical identities programmatically as described in Section C. We do this in the file example2 1 frustrated loop orbits.pyby computing automorphisms of an auxiliary graph.The result is a mapping O of qubits and couplers to orbits.If spins s i and s j satisfy O(s i ) = O(s j ), then in an ideal annealing experiment m i = m j .Likewise, if couplers s i s j and s k s satisfy O(s i s j ) = O(s k s ), then they have identical frustration probabilities f i,j = f k, .The code also gives us a mapping of orbits to "opposite" orbits, such that if spins s i and s j are in opposing orbits, m i = −m j , and if couplers s i s j and s k s are in opposing orbits then f i,j = f k, and J i,j = −J k, .Running the code on three disjoint copies of a frustrated six-qubit loop tells us that all AFM couplers are in one coupler orbit O 1 = {s i , s j | O(s i s j ) = 1}, and all FM couplers are in its opposite, O 2 = −O 1 (see Fig. 9).
We point out an obvious but useful fact: If we are using multiple embeddings of an Ising model, then all copies of a given qubit are in the same orbit, and all copies of a given coupler are in the same orbit.Here we use disjoint embeddings, but they need not be disjoint: the embeddings could overlap, and be annealed in separate calls to the QPU.

Shimming
We can now approach the frustrated loop similarly to the unfrustrated loop: all qubits should have average magnetization zero, and all couplers should be frustrated with the same probability.Again, tuning FBOs and individual couplings helps to reduce bias in the system.This example shows how to exploit orbits for our shim.
There is one detail worth pointing out.In Eqs. ( 6) and ( 8), the terms m and f can be computed as averages over an orbit.If we are dealing with opposing qubit orbits O q and −O q , we can simply use m = 0, as we do in the first example.For opposing coupler orbits O c and −O c , we can compute f across the union of the two orbits.In this case, that means that f is the average frustration probability across all couplers.Fig. 10 shows the results of shimming FBOs and couplings for 165 parallel embeddings of a 16-qubit frustrated loop, using nominal coupling strength |J i | = 0.9.Here, both components of the shim show a marked improvement of statistical homogeneity.Taking moving means for 10 iterations at a time, we see that both σ m (standard deviation of qubit magnetization) and σ f (standard deviation of coupler frustration probability) decrease as a result of turning on the FBO shim and the coupler shim, respectively.

Finding orbits of an arbitrary Ising model
Code reference: example2 3 buckyball orbits.py.
Here we present an example of an Ising model that is read from a text file and run through our orbit-finding code.The user may want to edit this code to analyze other Ising models of interest.
Consider another antiferromagnetic Ising model (J ij = 1) with a Buckyball graph as its underlying structure and no linear fields (h i = 0).We apply the same methodology described in Section C to find its orbits.Fig. 11 visualizes the Buckyball model with its orbits labelled by text, as well as its signed Ising counterpart with coupling values encoded by colour.FIG.11.An antiferromagnetic Ising model with a Buckyball graph structure.The node and edge colours encode the resulting qubit and edge orbits respectively.All qubits are in the same orbit since the graph is vertex-transitive.There are only two coupler orbits: those couplers sitting between two hexagons, and those sitting between a hexagon and a pentagon.Code reference: example3*.py.
In the previous examples we demonstrated several key methods: • Finding qubit and coupler orbits.
We can now apply these tools to a nontrivial system: the triangular antiferromagnet (TAFM).This is a classic example of a frustrated 2D spin system.Moreover, the addition of a transverse field to a TAFM leads to order-by-disorder at low temperature [19,20].For this and other reasons, including qualitative similarity to real materials, the TAFM has been simulated extensively using quantum annealers [4,10].We will use it as an example to showcase several concepts in calibration refinement for quantum simulation: • Truncating and renormalizing Hamiltonian terms.
• Simulating logical versus embedded systems.
• Simulating an infinite system versus faithfully simulating boundary conditions.

A. EMBEDDING AS A SQUARE LATTICE
Code reference: example3 1 tafm get orbits.py.
In D-Wave's Advantage systems, we can minor-embed the TAFM using two-qubit FM chains.First, we will embed a 12 × 12 square lattice with cylindrical boundary conditions, then we ferromagnetically couple pairs of qubits with a strong coupling J F M .The cylindrical boundaries are very helpful in providing rotational symmetries that we can exploit in our calibration refinement methods (as in the 1D chains already studied).
The provided code uses the Glasgow subgraph solver to find embeddings of the 12 × 12 square lattice, but note that this can take several hours.For larger square lattices, up to 32 × 32 or even larger depending on the location of inoperable qubits, one can inspect embeddings of smaller lattices and generalize the structure, since subgraph solvers are unlikely to be efficient at that size.We proceed with 10 disjoint 12 × 12 embeddings generated by the code.
In this example we will set AFM couplers to J AFM = 0.9, and all FM couplers to J FM = −2 * J AFM .Since FM couplers are very rarely frustrated in this system, we will only shim the AFM couplers.

B. ANNEALING WITH AND WITHOUT SHIMMING
Code reference: example3 2 tafm forward anneal.py.As in the previous example, we will compare performance of three methods: no shim, FBO shim only, and FBO and coupler shims together.We perform 800 iterations, turning on the FBO shim after 100 iterations and the coupler shim after 300 iterations.Fig. 13 shows data for this experiment, and we can see that as with the frustrated loop example, shimming improves statistical homogeneity of magnetizations and frustration.Note, however, that there is no appreciable impact on the average magnitude of the order parameter |ψ| .This will change when we vary boundary conditions (see Fig. 16).

C. MANIPULATING ORBITS TO SIMULATE AN INFINITE SYSTEM
Code reference: example3 2 tafm forward anneal.py.
The shim shown in Fig. 13 used coupler orbits for the square lattice with cylindrical boundaries, which are naturally different for couplers that are different distances from the boundary, or different orientations with respect to the boundary (and to FM chains).But what if we want to simulate, to the extent possible, an infinite TAFM?In that system, a coupler's probability of frustration is independent of its orientation and position, unlike in the square-lattice embedded system.We can simulate this case by putting all AFM couplers in one orbit, and all FM couplers in a second orbit, and proceeding as before.The coupler orbits no longer reflect the structure of the programmed Ising model, but rather the structure of the Ising model we wish to simulate.
Results for the "infinite triangular" shim are shown in Fig. 14.This experiment is performed just like the previous The coupler shim deviates significantly from nominal values (note axis scale), and has not converged even after 500 iterations.

Truncating and renormalizing couplers
In this code example (and others) we use an important method in the coupler shim: truncation.Programmed couplings must be in the range [−2, 1], so AFM couplers must remain less than 1, which is 1.11 * J AFM .Therefore, when couplers go out of range, we truncate them to within the range.To avoid persistent shrinking of the couplings due to truncation, we renormalize to the correct average coupling value (0.9) before truncation-this prevents cumulative shrinkage over many iterations.

Better initial conditions
Looking at the data, we can see that the most reduced couplers are those on the boundary.This suggests that if we want to simulate the infinite TAFM, we should start with a thoughtful setting of couplers.In this case, setting the AFM couplers on the boundary to J AFM /2 reduces the need to shim enormously.This makes sense, since doing so maximizes the ground-state degeneracy of the classical system, as previously noted [4].This shim is shown in Fig. 15.The experiment is performed just like the previous one, but with the parameter We can see that now, the coupler shim only deviates a few percent from nominal, at most.-15, we plot the evolution of the average magnitude ψ , as well as complex histograms of ψ (showing only data for one of the ten embeddings) before and after shimming.

Complex order parameter
Order in the TAFM can be characterized by a complex order parameter ψ, which we define now.Let c : S → {0, 1, 2} be a 3-coloring of the spins of the TAFM, mapping them onto three sublattices so that no two coupled spins are in the same sublattice (this coloring is unique, up to symmetries).Then for a spin state S we can define where c i = c(s i ) and i = √ −1.Due to symmetries among the sublattices arising from the cylindrical boundary condition, as well as up-down symmetry of spins since h = 0, we expect sixfold rotational symmetry (among other symmetries) in the distribution of ψ in an ideal annealer.Thus ψ can serve as a good indicator of any biases in the system, as well as global ordering.
We can use ψ to compare the "embedded finite" shim and "triangular infinite" shim, as seen in Fig. 16.Although we are simply forward-annealing the system, and therefore not sampling from the mid-anneal Hamiltonian, we expect the same characteristic ring histogram-without a peak near ψ = 0-that is seen in the quantum system (cf.[4] Fig. 3c).This is seen only after the "triangular infinite" shim.We mainly attribute this to the halving of the boundary couplings.In all cases, the shim improves the theoretically expected sixfold rotational symmetry of ψ.
It is often difficult or impractical to determine appropriate step sizes a priori.Here we demonstrate a simple method for adapting step sizes based on statistics of the shim.Note that due to noise in the QPU's surrounding environment, there is no well-defined asymptote or steady state for a shim.However, we can loosely assume that such a state exists: we expect high-frequency fluctuations in the environment to be small compared to low-frequency fluctuations and static cross-talk.
If the step size is sufficiently small and we are sufficiently close to the steady state, we can expect fluctuations of the Hamiltonian terms (FBOs, couplers, or fields) to behave like unbiased random walks.In an unbiased random walk with position x(t) at time t = 0, 1, . .., the probability distribution of x(t) approaches the normal distribution with mean 0 and variance t.
If the shim is far from the steady state and has a relatively small step size, the random walks will be biased in one direction, and thus the variance of fluctuations will grow superlinearly in t.Finally, if the step size is very large, then it will tend to overshoot the steady state, and oscillate.This leads to variance of fluctuations growing sublinearly in t.Thus we can periodically adjust the step size of a shim as follows, using a 20-iteration lookback and a tuning term ε = 0.1:

Part V A survey of additional methods
We have provided detailed demonstrations and free-standing Python implementations for several worked examples.These cover the basics of calibration refinement.Here we discuss some additional methods that have been used successfully in recent works.

A. SHIMMING A SYSTEM IN A UNIFORM MAGNETIC FIELD
Certain Ising models in a uniform magnetic field are of interest to physicists, and these have been simulated in quantum annealers both at equilibrium [5] and out of equilibrium [22].If we want to simulate an infinite system, we would ideally study a large system with no missing spins, and with fully periodic boundaries.However, this is often not possible, so we wish to make the magnetization m i independent of the spin's position relative to the boundary (although it may depend on the spin's position in a unit cell of the lattice being simulated).
To deal with this, we can shim individual longitudinal field terms, h i , such that all spins of a given type (i.e. in the same position of the unit cell) are in the same orbit.We can then shim all h i terms for each simulated field magnitude h that we want to study.Then the average value of h i is forced to remain at h throughout the shim, perhaps with an adjustment arising from boundary spins.
To shim the case h = 0, we use FBOs (as in the worked examples) instead of tuning h i .We can additionally ensure that each h i is a locally smooth function of h by adding a smoothing term.For example, if h i has values h − i and h + i for the next lower and higher values of h being simulated, we can make the adjustment for some small constant > 0.

B. SHIMMING AN ISING MODEL WITH NO SYMMETRIES
In Ref. [9], a qubit spin ice was implemented using a checkerboard Ising model.The system had open boundary conditions and missing spins due to inoperable qubits, so no geometric symmetries were available.However, due to the rich automorphism group of the qubit connectivity graph (ignoring unused qubits), it was possible to generate many distinct embeddings of the same system, using different mappings of qubits to spins.Therefore we could simulate a collection of distinct embeddings (in this case, 20) and shim in the same way we did in the worked examples.The only difference is that in the qubit spin ice example, the embeddings are not disjoint and therefore must be sampled from using separate calls to the QPU.However, once we have a set of samples from each embedding, we can analyze the data as though the embeddings are disjoint, whether or not this is actually the case.The benefit remains the same: by simulating with 20 distinct embeddings, we get qubit and coupler orbits of size at least 20.

C. SHIMMING A COLLECTION OF RANDOM INPUTS
In Ref. [11], shimming was used to study spin-glass ensembles-collections of random problems with certain parameters.As we have seen, we can spend hundreds of iterations shimming a single problem, and this becomes impractical when studying ensembles of thousands of instances.
The approach used was to exploit a common symmetry: all problems in the ensembles had h = 0. Shimming the couplers was abandoned as being impractical for such a large set of inputs.Shimming FBOs, however, is straightforward.By cycling through 300 spin-glass realizations using the same set of qubits and couplers, simulating each realization several times, it is possible to combine the work and arrive at a good set of FBOs that mitigates the majority of systematic offsets.

D. SHIMMING ANNEAL OFFSETS FOR FAST ANNEALS
As described in the Supplementary Materials to Ref. [10], D-Wave quantum annealing processors have recently demonstrated the capacity to anneal much faster than currently generally available, at an anneal time of 10 nanosec-onds or less [10,11].This speed exceeds the ability of the control electronics to synchronize the annealing lines (eight in the Advantage processor, four in D-Wave 2000Q™) satisfactorily.Therefore, frustration statistics can be used to infer which lines are out of sync with the others, and in which direction.Anneal offsets, which allow individual qubits to be annealed slightly ahead of or behind other qubits, were used to synchronize the qubits on each annealing line.These fast anneals are not currently generally available, but they may be in the future.

A. MAKING CALIBRATION REFINEMENT MORE EFFICIENT
As we have seen, shimming can take many iterations to converge.Naively repeating the process across many combinations of parameters (e.g., annealing time, energy scale, etc.) can be extremely time consuming.However, there are ways to improve the efficiency of the process.Here we outline some important things to bear in mind.

Adjustments are often continuous functions of other parameters
If we determine a set of adjustments for a given experiment, then slightly vary some parameters of the experiment, we can generally expect that the adjustments will not change much.For example, FBOs and coupling adjustments are expected to vary smoothly as functions of annealing time, energy scale, and various perturbations to the system (for example the ratio between FM and AFM couplers in an embedded triangular antiferromagnet).An important example is the annealing parameter s, in cases where we simulate a system at 0 s 1 ([4, 8, 10] etc).As an example of how this can help speed up a shim, if we double the annealing time, FBOs and coupling adjustments will remain relatively stable.Thus, rather than starting our shim anew from the nominal Hamiltonian, we can start from an adjusted Hamiltonian that was determined using similar parameters.One could go further than this, and extrapolate or interpolate based on multiple values.

Predictable adjustments should be programmed into the initial Hamiltonian
As shown in Figs. 14 and 15, starting with halved boundary couplings can immediately bring the couplings close to their converged values.If we are aware of such adjustments, using them as initial conditions can make shims converge far faster.

B. DAMPING SHIM TERMS
It is sometimes useful to gently encourage a shim to remain close to the nominal values, for example to prevent drifting Hamiltonian terms.This issue can be particularly important near a phase transition, where statistical fluctuations can be very large.Drift can be suppressed by adding a damping term to the shim.For example, we can set a damping constant 0 ≤ ρ ≤ 1, and after every iteration we can move each coupler J ij towards its nominal value Ĵij : Doing this can discourage random fluctuations, but can also lead to undercompensation of biases.It is only recommend to use damping when the shim is otherwise badly behaved.
Part VII

Conclusions
In this document we have presented several basic examples that introduce the value of calibration refinement or "shimming" in quantum annealing processors.These methods should be applied to any detailed study of quantum systems in a quantum annealer, and will generally provide a significant improvement to the results.Depending on the sensitivity of the system under study, these methods can mean the difference between an unsuccessful experiment and an extremely accurate simulation.
We have provided fully coded examples in Python, which should be easy to generalize and adapt.As part of these examples, we include methods for embedding many copies of a small Ising model in a large quantum annealing processor.This is a valuable and straightforward practice that can enormously improve both the quantity and the quality of results drawn from a single QPU programming.
Another important perspective, which has been introduced here for the first time, is the notion of constructing an auxiliary Ising model and using automorphisms of it to infer qubit and coupler orbits automatically.We encourage users to experiment with this method and report on any challenges or benefits found.
The examples in this document are written for use in an Advantage processor, but are not specific to that model, or even to D-Wave quantum annealers in general.These results may prove useful in diverse analog Ising machines, both quantum and classical.

FIG. 1 .
FIG.1.Annealing schedule for Hamiltonian 3 in a D-Wave™ Advantage™ processor.Γ(s) and J (s) control the magnitude of quantum fluctuations and the Ising energy scale, respectively.These values vary slightly from one processor to another.

FIG. 2 .
FIG. 2. (a) Ferromagnetic loop (periodic 1D chain) on L spins.(b) Frustrated loop, with one antiferromagnetic coupler.The FM loop has a twofold-degenerate ground state (all spins up or all spins down) with no frustration; the frustrated loop has 2L ground states, each with one frustrated bond.When h = 0, all qubits trivially have zero average magnetization in an ideal annealer.

FIG. 6 .
FIG. 6. Balancing qubits in a FM chain with flux-bias offsets.Iterative correction of qubit biases is demonstrated using three step sizes αΦ for 100 iterations: 10 −4 (left), 10 −5 (middle), 10 −6 (right).Step size is set to zero for the first 10 iterations.(Top) evolution of flux-bias offsets for 64 qubits in a FM chain.(Middle) Qubit magnetization averaged over first 10 iterations and last 10 iterations.(Bottom) Standard deviation of qubit magnetizations per iteration.
FIG.7.Balancing qubits and couplers in a FM chain with flux-bias offsets and coupler adjustments.This experiment is similar to that shown in Fig.6, but with αJ > 0 for the last 100 iterations.Couplers remain distributed about the average value of J = −0.2.

FIG. 10 .
FIG.10.Shimming a frustrated loop.300 iterations are performed.A flux-bias offset shim is used after iteration 100, and a coupler shim is used after iteration 200.Nominal couplings are ±0.9.The third panel shows the standard deviation of qubit magnetizations taken as a moving mean over 10 iterations, σm.The fourth shows the corresponding quantity σ f for frustration probability.
FIG.12.A 12 × 12 square lattice with cylindrical boundary conditions (periodic top/bottom).Contracting two-qubit FM chains into single spins results in a triangular antiferromagnet.

FIG. 13 .FIG. 14 .
FIG.13.Shimming an embedded cylindrical triangular antiferromagnet.800 iterations are performed.A flux-bias offset shim is used after iteration 100, and a coupler shim is used after iteration 300.For clarity, we only show FBOs for 12 qubits, and couplings for 12 couplers in the same orbit.Standard deviation of frustration probabilities, σ f , is computed for the couplers in each orbit, and the average over all orbits is taken.

FIG. 16 .
FIG.16.Complex order parameter ψ.For the three shims shown in Figs.13-15, we plot the evolution of the average magnitude ψ , as well as complex histograms of ψ (showing only data for one of the ten embeddings) before and after shimming.
Construction of signed Ising model.To detect exploitable symmetries, we search for automorphisms of an auxiliary Ising model in which each spin is duplicated into itself and its negation; each coupler is then expanded to four copies of itself, two of them negated.Automorphisms of the auxiliary Ising model can be detected by conversion into an equivalent automorphism-finding problem on an edge-labeled graph.Here, vertex labels indicate the identities of spins, and show how each spin is duplicated for the signed Ising model.