Quantum annealing (QA; Kadowaki and Nishimori, 1998; Johnson et al., 2011) 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 as follows:

Here, {σiz}i=1N∈{-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 HP is intractable, i.e., NP-hard (Barahona, 1982).

QA adds to HP to an initial driving Hamiltonian as follows:

The ground state of HD, which is a uniform quantum superposition of all classical states, is easy to prepare. QA guides a time-dependent Hamiltonian H(s) from HD to HP by linearly combining HD and HP as follows:

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). Units are GHz, convertible to Joules by multiplication by ℏ (reduced Planck constant). An example is shown in Figure 1.

Quantum processing units (QPUs, in this case quantum annealers) are typically made available with a single one-size-fits-all calibration. Non-idealities 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. Moreover, 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 non-idealities. 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 study, is defined simply as one that respects symmetries in the Hamiltonian—each qubit behaves identically and each coupler behaves identically.

Variations on the methods described herein have been used in many studies (King et al., 2018, 2021a,b,c, 2022, 2023; Kairys et al., 2020; Nishimura et al., 2020). Often, when behavior of the system relies on precise maintenance of energy degeneracy between states, or energy splitting from the transverse field, the results are highly sensitive to these tunings. Particularly for the simulation of exotic magnetic phases, calibration refinement is an essential ingredient of successful experiments. However, so far the discussion of these methods has mostly been relegated to Supplementary material. Here, our aim 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:

The tutorial is organized as follows. In the remainder of this section, we introduce concepts that form the bases of the QPU calibration procedure. In Section 2, we illustrate the essence of our approach through a toy example. In Section 3, we extend the method and improve calibration efficiency by exploiting symmetries in a given model. In Section 4, we introduce a non-trivial system to demonstrate additional concepts useful for realistic applications. Collectively, these sections provide a comprehensive walkthrough of the calibration procedure. In Section 5, we survey additional methods for narrower use cases of the QPU. Finally, we provide practical tips and considerations in Section 6 and conclude the tutorial in Section 7.

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 theory, symmetries in an Ising model admit identical expectation values of observables.¹ However, empirical averages over many realizations of these observables (from a QPU) may differ systematically.² We tune Hamiltonian terms to reduce the discrepancy between the expected value and observed averages. For example, given a Hamiltonian consisting of a single qubit with no bias, HP=0(s1), we expect to observe a mean spin of 0 for s₁. However, this observed quantity may deviate from 0 systematically; we, thus, attempt to correct this deviation by perturbing the Hamiltonian.

In this study, 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 as follows:

For each spin s_i, a frustration probability is as follows:

For each coupler J_i,j; f_i,j is the observed probability of the coupler having a positive contribution to the energy in HP.

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 relabelling and gauge transformation (discussed below).³ We understand and formalize these symmetries—and automate their detection—through graph isomorphisms (especially automorphisms) and generalizations (Godsil and Royle, 2001). We understand and formalize these symmetries—and automate their detection—through graph isomorphisms (especially automorphisms) and generalizations. We briefly introduce the concept of graph automorphisms and orbits below (see Godsil and Royle, 2001 for a more complete treatment). Notably, the symmetries we find and exploit here are a subset of all possible symmetries.

Given a graph G = (V, E) with vertex and edge sets V, E, a graph automorphism is a mapping π:V ↦ V such that (π(u), π(v)) ∈ E if and only if (u, v) ∈ E. Intuitively, a graph automorphism is an adjacency-preserving relabelling of vertices. Two vertices u, v ∈ V are said to belong in the same vertex orbit if there exists an automorphism mapping u to v (or v to u). If u, v ∈ V belong in the same vertex orbit, the edges incident to u or v also belong to the same edge orbit.

We now relate the definitions of graph automorphisms and orbits back to our goal of detecting and exploiting symmetries. These symmetries admit two types of equivalence relations on an Ising model HP: one on the qubits and the other on the couplers. We call the equivalence classes qubit orbits and coupler orbits, respectively. We use notation O(si) for a qubit orbit containing spin s_i, and O(si,sj) 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:

Formally, a set O is said to be a qubit orbit if si,sj∈O, then m_i = m_j. Similarly, a set O is said to be a coupler orbit if (si,sj),(sk,sl)∈O, then f_i,j = f_k,l.

For example, consider the Hamiltonian HP=hs1-hs2+hs3 for h ≠ 0. The two independent spins s₁, s₃ can be trivially relabeled (permuted) by each other, thus the two qubits belong to the same qubit orbit; s₂ belongs in its own qubit orbit as it does not have the same magnetization as s₁, s₃. Now, let us consider the Hamiltonian HP,2=HP+Js1s2+Js2s3 for J ≠ 0. In this case, (s₁, s₂) and (s₂, s₃) exist in the same coupler orbit because s₁, s₃ can be swapped while preserving the couplers in HP,2; the permutation preserves adjacency structures. As a non-example, let us consider the Hamiltonian HP,3=HP+Js1s2-2Js2s3. In this case, (s₁, s₂) and (s₂, s₃) no longer exist in the same coupler orbit because swapping s₁, s₃ no longer preserves the couplers in HP,3.

Due to spin-flip symmetries, or spin reversal transformations (SRTs; described below), each qubit and coupler orbit can additionally have up to one non-empty orbit that is opposite.

An SRT, as its name suggests, flips the sign of a spin. For example, consider the Hamiltonian with a single qubit HP=hs. An SRT transforms on s yields an identical Hamiltonian HP=-hs~ where s~=-s. Similarly, for a Hamiltonian consisting of both biases and coupling terms such as HP=h1s1+h2s2+J1,2s1s2, we can apply an SRT on one (or multiple) variable(s) to obtain HP=-h1s~1+h2s2-J1,2s~1s2=h1s1-h2s~2-J1,2s1s~2=-h1s~1-h2s~2+J1,2s~1s~2, where s~1=-s1,s~2=-s2.

In the earlier example HP=hs1-hs2+hs3, qubit s₂ belongs to the orbit opposite of s₁, s₃'s orbit. We will sometimes overload notation, conflating O(si) with O(i) and O(si,sj) 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.

Notably, 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 Figure 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. In each case, all qubits and couplers have, respectively, identical magnetization and frustration probabilities. The unfrustrated case is trivially true. The frustrated scenario is less obvious but can be verified by computing frustration probabilities for each edge.

Having defined qubit and coupler orbits, we now consider how to find them.

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 (McCreesh et al., 2020), which has been incorporated into the embedding finding module minorminer (D-Wave, 2023). To make the embedding search faster, we raster-scan across 2 × 2 blocks of unit cells in the QPU's Pegasus graph (Boothby et al., 2020) and then greedily construct a large set of non-intersecting embeddings. The file embed_loops.py provides a code example that finds multiple disjoint copies of a 64-qubit loop in A_QPU.

Code reference: example1_1_fm_loop_balancing.py.

We will use simple parameters for the experiment, running 1 μs anneals forward anneals (where s increases linearly in time as t/t_a) 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 Equation 1; second, with a programmable flux-bias offset (FBO) Φ_i (Harris et al., 2009; D-Wave, 2022). 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 as follows:

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.

The choice of a step size α_Φ has a strong influence on the convergence of the calibration procedure. 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 1 × 10⁻⁴, in units of Φ₀) is too large and creates oscillations in Φ_i and m_i. One (1 × 10⁻⁶) is too small and takes many iterations to converge. One (1 × 10⁻⁵) 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 (Section 4.4). For best results, we should ensure:

When seeking evid ence 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 Figure 6 shows the average magnetization of each qubit across the first and last 10 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 αΦ=1×10-4 causes broad spreading of biases, so we need to be careful with our step sizes.

Code reference:

example1_2_fm_loop_correlations.py.

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 as follows:

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 fi,i+1-f̄:

Figure 7 shows data for the same experiment as Figure 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 αΦ=1×10-5 and α_J = 0, then 100 with αΦ=1×10-5 and α_J = 0.001. In this particular case, the coupler shim is small but some systematic signals can be observed. We will show more impactful cases later in the tutorial.

Code reference:

example2_1_frustrated_loop_orbits.py.

We can derive this fact as follows. Flipping the sign of s₂, 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; 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.

For more complicated examples, we would prefer to find such statistical identities programmatically as described in Section 1.3. We do this in the file

example2_1_frustrated_loop_orbits.py

by 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(si)=O(sj), then in an ideal annealing experiment m_i = m_j. Similarly, if couplers s_is_j and s_ks_ℓ satisfy O(sisj)=O(sksℓ), 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_is_j and s_ks_ℓ are in opposing orbits, 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 O1={si,sj∣O(sisj)=1}, and all FM couplers are in its opposite, O2=-O1 (see Figure 9).

We point out an obvious but useful fact: If we are using multiple embeddings of an Ising model, 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.

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

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. Figure 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 Figure 16).

Code reference: example3_2_tafm_forward_anneal.py.

The shim shown in Figure 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 Figure 14. This experiment is performed just like the previous one, but with the parameter

shim['type']='triangular_infinite'

instead of

shim['type']='embedded_finite'.

The coupler shim deviates significantly from nominal values (note axis scale), and has not converged even after 500 iterations.

Code reference: example3_2_tafm_forward_anneal.py.

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 act as though 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:

In the example code example3_2_tafm_forward_anneal.py, this method is applied by setting

This check is done every iteration, but this is not necessary.

Adaptive step sizes are so far a largely unexplored research area, and various approaches could be taken. Using different step sizes for each orbit is certainly worth exploring; note in Figure 14 that different coupler orbits have hugely varying deviations from the mean. More general frameworks like “Adam” (Kingma and Ba, 2014) could also be useful in this context.

Certain Ising models in a uniform magnetic field are of interest to physicists, and these have been simulated in quantum annealers both at equilibrium (Kairys et al., 2020) and out of equilibrium (King et al., 2021a). 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). In a typical experiment, we want to measure a system under an average field h̄ for each value in an increasing set of equally spaced values {h̄(1),h̄(2),…,h̄(m)}.

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̄(k) that we want to study, and denote each individual term hi(k). After each iteration we renormalize the fields so the average value 1N∑i=1Nhi(k) remains equal to h̄(k) 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. To shim the case h̄(k)=0, we use FBOs (as in the worked examples) to set a zero point instead of tuning hi(k). In doing so, we can compensate for time-dependent flux drift, particularly when determining the location (in h̄) of a phase transition.

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 hi- and hi+ for the next lower and higher values of h̄ being simulated, we can make the adjustment.

We can additionally ensure that each sequence {hi(1),…,hi(m)} is a smooth function of h̄ by adding a smoothing term.⁵ For example, we can add make an adjustment in two steps. First, set all

for intermediate values of k and some smoothing constant 0 < ϵ < 1. Second, set all hi(k)←h~i(k).

In King et al. (2021b), 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.

In King et al. (2023), 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.

As described in the Supplementary material to King et al. (2022), D-Wave quantum annealing processors have recently demonstrated the capacity to anneal much faster than currently generally available, at an anneal time of 10 nanoseconds or less (King et al., 2022, 2023). 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.

The shimming techniques we have discussed in this tutorial are efficient to perform and are simple to justify theoretically. However, we note a couple of its limitations here. First, there exists a potential computational bottleneck in the full shimming paradigm. The method partly relies on identifying qubit and coupler orbits. Identifying these orbits can be reduced to the problem of identifying the automorphism group of a graph, which can usually be done quickly in practice even though no general polynomial-time algorithm is known. Second, shims inherently exhibit time-dependent fluctuations due to noise in the environment. These fluctuations tend to be smaller than the shim terms themselves, so do not significantly diminish the potential benefit of calibration refinement.

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.

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 toward its nominal value Ĵ_ij:

Doing this can discourage random fluctuations, but can also lead to under compensation of biases. It is only recommend to use damping when the shim is otherwise badly behaved. In practice, a suitable value for ρ is determined through trial and error, i.e., by assessing whether shims have converged (see Section 2 for discussion on shim convergence).