Quantum error mitigation in quantum annealing

Open quantum model

In this appendix, we use open quantum modeling to explain ZNE. We write the total Hamiltonian as

where H_S and H_B are system and bath Hamiltonians, respectively. We write the interaction Hamiltonian in a general form as

where \({{\mathcal{O}}}^{\alpha }\) is an operator acting on the system and Q^α is the noise operator acting on the corresponding environment. The operators \({{\mathcal{O}}}^{\alpha }\) are typically Pauli matrices, e.g., \({\sigma }_{\alpha }^{z}\) for flux noise and \({\sigma }_{\alpha }^{y}\) for charge noise. For static noise, such as control error, Q^α can be considered as a constant random energy.

In open quantum systems, the state of the system is described by the reduced density matrix ρ(t). We choose the preferred basis to be the eigenstates \(\left\vert n\right\rangle\) of the system Hamiltonian H_S with eigenvalues E_n. The Bloch-Redfield equation for the reduced density matrix ρ(t) is written as³⁷

where ω_nm = E_n − E_m and the Redfield tensor is defined as

with

Here, \({{\mathcal{O}}}_{nm}^{\alpha }\equiv \left\langle n\right\vert {{\mathcal{O}}}^{\alpha }\left\vert m\right\rangle\) and the noise spectral density is defined as

For the case of uncorrelated heat baths we have

For a time-dependent system Hamiltonian, the generalized Bloch-Redfield equation becomes⁴⁵

where

In the energy basis, we write the density matrix as a linear vector with 2^2N elements:

The master equation becomes a matrix equation:

where

is a diagonal matrix made of energy differences and \(\hat{M}\) and \(\hat{R}\) are the basis-rotation and relaxation matrices corresponding to M_nmkl and R_nmkl, respectively. A formal solution to this equation is

The time-ordering operator \({\mathcal{T}}\) is introduced to take care of non-commuting \(\hat{\gamma }(t)\), \(\hat{M}(t)\), and \(\hat{R}(t)\) at different times. Changing the integration variable to s = τ/t_a, we get

The term \(\gamma \hat{R}(s)\) in the exponent takes care of all relaxation and dephasing processes during the evolution, with γ appearing as a common factor in all time-scales. In the limit of small γ when those decay processes are not very strong, it is possible to assume this term is small and expand the exponent. Using Taylor expansion, we can write

where

is the noise-free contribution to the density matrix and

captures the effect of noise to order μ in perturbation. Turning back the density matrices from vector form to matrix form we obtain (5).

In some cases, only a small number of energy states get occupied during the annealing process. Therefore, we only need to incorporate relaxation and dephasing processes between those states. The experiment reported in ref. ³³ is one such case. The evolution is limited to the lowest two energy states, during which the relaxation process is extremely slow. This allows expansion in powers of those relaxation rates only, although other decay processes could be much faster. As a result zero-T extrapolation worked for time scales millions of times longer than the expected single-qubit decoherence time of the qubits.

Practical energy-time rescaling

The goal of this section is to introduce a relation between energy and time rescaling so that the critical dynamics remain unchanged. We first write the time-dependent Hamiltonian in a dimensionless form with one dimensionless parameter τ_Q, which measures the speed of passing the critical point. The same quantity is used in some theory papers (see e.g., refs. ^42,43) and rederived in “Methods—Critical and post-critical dynamics of the transverse field Ising chain”, hence our formalism allows direct comparison with these analytical results.

We start by writing the Hamiltonian (1) as

where

is the dimensionless transverse field. The corresponding Schrödinger equation reads

Introducing a dimensionless time

we can write (39) in dimensionless form

where

is the dimensionless Hamiltonian. Note that \(\tilde{t}\in [0,{\tilde{t}}_{a}]\) with \({\tilde{t}}_{a}={t}_{a}\mathop{\int}\nolimits_{0}^{1}{\mathcal{J}}(s)ds\) being the dimensionless annealing time.

For fast annealing, only the speed of passing the critical point determines the critical dynamics. We define the critical point s = s_c by

where g_c is the critical value of g. Using linear expansion near the critical point we write

where τ_Q is the dimensionless quench time scale and \({\tilde{t}}_{c}=\tilde{t}({s}_{c})\). It is possible to express Kibble-Zurek dynamics only in terms of τ_Q. In other words, two systems with the same τ_Q are expected to yield the same statistics as long as only critical dynamics affect their evolution. For a slower evolution where the dynamics outside the critical region affect the results, details of the schedule become important.

We expand the schedule near the critical point as

where Γ_c ≡ Γ(s_c), \({\Gamma }_{c}^{{\prime} }\equiv {[{\partial }_{s}\Gamma (s)]}_{s = {s}_{c}}\), etc. Expanding g(s) near s_c, we obtain

On the other hand

Therefore,

Comparing with (44), we obtain

where

defines a time scale that only depends on the annealing schedule. For a 1D spin chain with a problem Hamiltonian

the critical point is the point where \({\Gamma }_{c}={\mathcal{J}}({s}_{c})\), thus g_c = 1. On the other hand, for 2D and 3D problems we have g_c > 1 and, therefore, we expect the critical point to occur earlier in the schedule.

In order for the two systems to have the same critical dynamics, they need to have the same τ_Q. Equation (50) requires that t_a should be scaled the same way as t_Q in (51). If both Γ(s) and \({\mathcal{J}}(s)\) are scaled the same way such that H_S → H_S/κ, we obtain t_Q → κt_Q, therefore we need λ = κ, as expected for simple rescaling.

For the more practical case of only rescaling H_P, it is equivalent to \({\mathcal{J}}(s)\to {\mathcal{J}}(s)/\kappa\) without changing Γ(s). Obviously, the critical point will be shifted to a new point \({s}_{c}^{\kappa }\) which is the solution to \(\Gamma ({s}_{c}^{\kappa })={g}_{c}{\mathcal{J}}({s}_{c}^{\kappa })/\kappa\). Thus, one needs to calculate the new \({t}_{Q}^{\kappa }\), obtained from (51) at \({s}_{c}^{\kappa }\). The correct time rescaling coefficient will therefore be:

For the case of the 1D transverse field Ising problem, we can write

where \({\Gamma }_{c}^{\kappa {\prime} }={\Gamma }^{{\prime} }({s}_{c}^{\kappa })\), etc. We should emphasize that the above analysis is based on the assumption that Kibble-Zurek dynamics can be completely characterized by a single time-scale. In reality, there are other time scales that are relevant to post-critical dynamics, which may explain the small t_a-dependence in Fig. 8a.

Critical and post-critical dynamics of the transverse field Ising chain

Quantum Ising chain

The transverse field quantum Ising chain is⁴⁶

where we assume periodic boundary conditions, \({\overrightarrow{\sigma }}_{N+1}\,=\,{\overrightarrow{\sigma }}_{1}\), and even N for definiteness. For N → ∞ the quantum critical point at \({\Gamma }_{c}/{{\mathcal{J}}}_{c}=1\) separates the paramagnetic and ferromagnetic phases. Note that the programmed ferromagnetic coupling (J_ij = −1.8/κ) is absorbed into the definition of \({\mathcal{J}}\) for purposes of this section. After the Jordan-Wigner transformation:

introducing fermionic annihilation operators c_n, the Hamiltonian (55) becomes

Here the projectors

define subspaces with even/odd numbers of fermions as the parity, related to \({\prod }_{n}{\sigma }_{n}^{x}\), is a good quantum number.

are the corresponding reduced Hamiltonians. Fermionic operators c_n in H^± satisfy (anti-)periodic boundary conditions: c_N+1 = ∓ c₁.

The ground state has even parity for any non-zero Γ. For a time evolution that begins in the ground state, the state remains in the even parity subspace. Terms relevant to H⁺ can be expressed by an anti-periodic Fourier transform

where pseudomomenta take half-integer values

The Hamiltonian can thereby be written in additive form,

where

Each Hamiltonian H_k lives in a 4-dimensional subspace but its ground state, and also the state excited during the ramp, belongs to a 2-dimensional subspace spanned by

where \(\left\vert {0}_{k}\right\rangle\) is a state without fermions: \({c}_{\pm k}\left\vert {0}_{k}\right\rangle =0\). The eigenstates satisfy stationary Bogoliubov-de Gennes (BdG) equations:

where eigenvalues ε are minus eigenenergies of Hamiltonians in Eq. (63). There are two eigenenergies corresponding to ε = ± ε_k, where

is the quasiparticle dispersion.

A fermionic kink operator⁴³ that annihilates a kink on a bond connecting sites n and n + 1 is defined

With the Jordan-Wigner transformation (56), (67) becomes

and its quasimomentum representation is

In the ferromagnetic Ising chain at g = 0, Eq. (65) has two eigenstates with eigenvalues \(\varepsilon =2{\mathcal{J}}\) and \(\varepsilon =-2{\mathcal{J}}\), respectively,

The ground state is a no-kink vacuum, \({\gamma }_{\pm k}\left\vert {\rm{GS}}\right\rangle =0\), and the excited state is a pair of kinks: \(\left\vert {\rm{ES}}\right\rangle ={\gamma }_{-k}^{\dagger }{\gamma }_{k}^{\dagger }\left\vert {\rm{GS}}\right\rangle\).

D-Wave QPU ramp

The transverse field Γ(s) and the ferromagnetic coupling \({\mathcal{J}}(s)\) are ramped with

where the time t runs from 0 to the annealing time t_a. The states in (64) evolve according to the time-dependent Bogoliubov-de Gennes equations:

Near the critical point, where \(\Gamma ({s}_{c})={\mathcal{J}}({s}_{c})\), the ramps can be linearized as

Here the prime means d/ds and \({\Gamma }_{c}={{\mathcal{J}}}_{c}\). The approximation is accurate for long enough t_a when the KZ excitation is completed near \({s}_{c}+\hat{s}\) soon after s_c, where \(\hat{s}\propto {t}_{a}^{-1/2}\). In the same regime, only small quasimomenta k ≪ 1 are excited in the range proportional to \(\pm {t}_{a}^{-1/2}\). To the leading nontrivial order in s − s_c and k, Eq. (73) becomes the Landau-Zener (LZ) problem,

where

The LZ excitation probability is

where the dimensionless KZ quench time is defined as

The final kink density is

The last equality defines the KZ length including not only the universal KZ power law, \({\tau }_{Q}^{-1/2}\), but also the non-universal numerical factor. Up to a phase, the final state at the end of the ramp is

and the overall state reads

Here \(\left\vert {\rm{GS}}\right\rangle\) is a vacuum for the kink annihilation operators, γ_k, i.e. the fully polarized ferromagnetic ground state.

The Gaussian state (81) is fully characterized by its quadratic correlators: normal

and anomalous

We approximate⁴²

and expand the final phase as⁴³

The anomalous correlator (83) becomes

where α = (19/20)²(3/2)³π = 9.57 and the length

The connected kink correlator can be assembled as^42,43

It has two scales of length: the shorter KZ length ξ and the longer l. The latter depends on the parameter φ⁽²⁾ in (85).

Derivation of φ
⁽²⁾

On the one hand, from the asymptote of the exact solution of (75), expressed by the Weber functions, one obtains the phase⁴⁶

accurate when Δ_kτ² ≫ 1. Here Γ is the gamma function. With the support of (77) we can further approximate⁴³

where γ_E is the Euler gamma constant. On the other hand, after the Landau-Zener excitation is completed around \({\Delta }_{k}{\hat{\tau }}^{2}=A\), where \(A={\mathcal{O}}(1)\), φ_k continues to grow as a dynamical phase:

where \(2{\epsilon }_{k}(\tau )=\sqrt{1+{\Delta }_{k}^{2}{\tau }^{2}}\) is the energy gap in (75). With (76), the crossover condition \({\Delta }_{k}{\hat{\tau }}^{2}=A\) becomes

This defines the parameter \(s={s}_{c}+\hat{s}\) after which the phase is approximated by the dynamical one in (91). The parameter A should be neither too small, to avoid the LZ excitation regime, nor too large, for (75) to remain accurate. In the next step we go beyond the approximate (73) that is accurate near the critical point only.

A quasiparticle dispersion that follows from the exact (73) can be expanded in powers of k²:

By the logic of (91) the final phase at s_f = 1 becomes

and we can identify the coefficients in (85) as

where

After isolating a logarithmic divergence of the integral near s_c, that appears with increasing annealing time t_a and decreasing \(\hat{s}\), we obtain:

Neglecting terms linear in \(\hat{s}\propto {t}_{a}^{-1/2}\) and higher, we arrive at

Here \(\ln \xi\) is a universal contribution to the phase scrambling from near the critical point and f is a non-universal one from the post-critical dynamics that does not depend on ξ but only on the schedule:

In general f(κ) depends on the scaling parameter κ in the transformation \({\mathcal{J}}\to {\mathcal{J}}/\kappa\). This is verified by the following examples.

Γ-linear ramp—As in refs. ^42,43 we assume a linear Γ(s) = Γ₀(1 − s) and a constant \({\mathcal{J}}(s)={\Gamma }_{0}/\kappa\). With s_c = 1 − 1/κ and s_f = 1 we obtain f(κ) equal to

It does not depend on κ.

Linear ramp—Here we assume that both parameters are linear functions: \(\Gamma (s)={\Gamma }_{0}+(s-{s}_{0}){\Gamma }_{0}^{{\prime} }\) and \({\mathcal{J}}(s)={{\mathcal{J}}}_{0}+(s-{s}_{0}){{\mathcal{J}}}_{0}^{{\prime} }\), where \({{\mathcal{J}}}_{0}={\Gamma }_{0}\) and \({s}_{f}={\Gamma }_{0}/(-{\Gamma }_{0}^{{\prime} })\). After the rescaling, \({\mathcal{J}}(s)\to {\mathcal{J}}(s)/\kappa\), the critical point moves to \({s}_{c}={s}_{0}+{\Gamma }_{0}(1-1/\kappa )/({{\mathcal{J}}}_{0}^{{\prime} }/\kappa -{\Gamma }_{0}^{{\prime} })\) and we obtain

For κ → 0 it diverges to + ∞ sending l → ∞ and making the anomalous part of the correlator negligible in this limit.

D-Wave QPU ramp—In order to avoid numerical instabilities in the integral (100), where the integrand is a finite difference of two functions that diverge for \(s\to {s}_{c}^{+}\), the annealing schedules \({\mathcal{J}}(s)\) and Γ(s) are approximated by smooth degree-10 polynomials in the range s = 0.025. . . 0.45. Furthermore, the integral is regularized by replacing its lower limit s_c with \({s}_{c}+\hat{s}\), where the parameter A in (97) is treated as a regulator: A → 0⁺. The coherent phase-scrambling length (87) as a function of the scaling parameter κ is plotted in Fig. 4, where it is converged for A = 10⁻⁸. The connected kink-kink correlator is shown in Fig. 5. Despite the theory being developed with corrections of order 1/ξ, and assuming the large system size limit, it quite accurately describes the results of exact simulation of the Bugoliubov-de Gennes equations for the experimental regime. Dependence on t_a is not insignificant through the experimental regime. Furthermore, the length scale l increases slightly as a function of κ in the range [1, 2.57] matching the experiment. Extrapolating backward in an experiment or simulation suppresses the length scale. We have assumed this is a small effect for purposes of experiments. Indeed, it seems to be small compared to the scale of the correction achieved, which is consistent with the statistical variation, with κ being dominated by noise rather than the post-critical (closed system) dynamical evolution.

The phase-scrambling length l(87) (divided by ξ) as a function of the scaling parameter κ. A broad maximum is located near κ = 4.

Top: the correlation for different annealing times t_a at the scaling parameter κ = 1. Bottom: the correlation for different scaling parameters κ at t_a = 7 ns. Theory has small deviations from the simulated (sim.) dynamics, consistent with corrections of O(1/ξ).

Special care is needed in the limit of small λ (κ), where the problem-energy scale diverges—particularly given that we consider dynamics initiated in a ground state at non-zero H_P. For practical purposes, we might avoid these technicalities by considering extrapolation only to ‘small enough’ λ such that the statistical error is no longer dominated by ‘correctable’ noise (but by sampling error, systematic control errors, etc).

Experimental details

Several methods are presented in this section, alongside code examples sufficient for the reproduction of the main-text experiments in problem-Hamiltonian rescaling, and additional experimental data.

For simulation and the evaluation of λ(κ), we use the schedules shown in Fig. 6. These are developed based upon single and coupled-pair qubit experimental measurements, a parametric rfSQUID model, and a parametric (persistent-current based) mapping from an rfSQUID to Ising model (online documentation is available at https://docs.dwavesys.com/docs/latest/index.html. This includes QPU-specific properties (schedules and parameters)). This methodology leaves room for uncertainty and systematic errors^18,19 (online documentation is available at https://docs.dwavesys.com/docs/latest/index.html. This includes QPU-specific properties (schedules and parameters)) that might explain some deviation in critical and post-critical phenomena of experimental results.

The schedules of the experimental and GA QPUs used in experiments parameterized by Γ(s) and \({\mathcal{J}}(s)\). The system Hamiltonian is \(H(s)/h=-\frac{1.8}{\kappa }{\mathcal{J}}(s){\sum }_{i=0}^{L-1}{\sigma }_{i}^{z}{\sigma }_{i+1 \,(mod \,L)}^{z}-\Gamma (s){\sum }_{i=0}^{L-1}{\sigma }_{i}^{x}\), and h is the Planck constant and σ^x,y are the Pauli matrices and L is the system size L = 232 (experimental), 1024 (general access).

Zero-noise extrapolation was performed using linear extrapolation in T (Fig. 1), quadratic extrapolation in λ (Figs. 2 and 7), and linear extrapolation in λ (Figs. 3, 8, 9) to zero; in each case, we used least-square methods.

Error mitigation of kink density on the experimental Advantage2 system. (left) Kink density is measured for various annealing times t_a and coupling energies J = −1.8/κ. An empirical rescaling of time is performed by collapsing fits of \(\rho \propto {t}_{a}^{-1/2}\) for t_a < 10. (center) The collapse achieved by rescaling, circles represent the extrapolated points at fixed t_a/λ. Center(Inset) ZNE is performed in λ, using a quadratic fit, Eq. (17). (right) λ(κ) obtained from collapse (circles) compared to those obtained from Eq. (15) (triangles).

Mitigation of kink correlation on the experimental Advantage2 system. (Left) Simulation results are obtained for \({\mathcal{C}}(r)\); see “Methods—Simulation of dynamics” for details. Dashed line indicates a degree-10 polynomial that is replicated on the other panels to the right to facilitate comparison between theory and experiment. (Center-left) Unmitigated QA data (κ = 1) show correct qualitative behavior, but with a significantly flattened peak. (Center-right) Results obtained via linear ZNE to λ = 0 show good agreement with simulation, extrapolated curves are shifted to smaller r as is also observed for the GA Advantage2 prototype. (Right) We show example extrapolations for varying lattice distances r, with t_a = 7.34 ns.

See description for Fig. 3. The kink density (top) is weakly dependent on the calibration refinement; there is small suppression in the kink density. The kink-correlation with calibration refinement at t_a/λ(κ) = 8 ns (center) has enhanced peak-height, consistent with reduced disorder. We can combine calibration refinement with the collapse method for t_a/λ(κ) = 8 ns, whereby λ kink-density is matched for all κ shown (described “Methods—Experimental details—Experimental Advantage2 prototype methodology at a scale L = 232”). This further enhances the peak height, but agreement with the 1-parameter theory (best fit l/ξ shown) is less convincing (at this t_a); one can also observe that the data is clustered rather than linearly spaced in κ, indicating the perturbative (linear) approach is inappropriate.

Data and simulation curves shown in Fig. 3 are interpolated with a fitted four-parameter functional form, allowing reasonable agreement,

qualitatively matching the theory of “Methods—Critical and post-critical dynamics of the transverse field Ising chain”. Theory predicts, up to corrections O(1/ξ), parameters f = 9.31ξ/l³, 3π/l², 1, 2π/ξ².

Experimental Advantage2 prototype methodology at a scale L = 232

Quantum annealing experiments (Figs. 1, 7, and 8) were performed using a prototype D-Wave Advantage2 processor, with 232 flux qubits coupled in a one-dimensional periodic chain. Refinements of calibration are performed as previously detailed in the supplementary materials of ref. ¹⁸, and qualitatively matched to the description of the section “Experimental details—Calibration refinement method”⁴⁷.

Average kink densities are measured over 20 programmings; 100 samples are taken in each programming. Kink correlations require more robust statistics, and are measured over 1000 programmings; 1000 samples are taken in each programming. All error bars represent 95% confidence intervals generated from bootstrap resampling of the 20 or 1000 programmings.

Qubit temperatures reported in Fig. 1 are measured by the population slope for uncoupled qubits subjected to a varying longitudinal field.

In the remainder of this section, we apply a problem-Hamiltonian rescaling on the same experimental Advantage2 system for which temperature extrapolation was performed (section “Zero-temperature extrapolation”), with a method assuming no knowledge of the schedule. This is an alternative to the schedule-derived method (15), and most appropriate when the schedule is not known with high confidence.

Figure 7. (left) shows measurements of kink density for different values of J = −1.8/κ as a function of annealing time. As expected, reducing the energy scale increases the number of kinks. Moreover, it increases the relative influence of the thermal environment as evidenced by the increasingly obvious deviation from \(n\propto {t}_{a}^{-1/2}\) at high t_a. Nonetheless, it is possible to rescale the horizontal axis t_a → t_a/λ, using λ as a free parameter, to collapse all of the data within some part of the diagram (the coherent annealing regime at small t_a, specifically we choose t_a < 10). We assume that this kink-density statistic is error-free, thereby eliminating the possibility of correction by extrapolation, but we obtain λ(κ) that can be applied at larger t_a for kink density (or at any t_a for other statistics such as C(r)). The resulting collapse is shown in Fig. 7(center) and the values of λ(κ) obtained by this collapse are plotted in Fig. 7(right) together with theoretical predictions using (15). Note that the collapse method will fail if static errors, such as Hamiltonian misspecification, significantly impact the coherent regime. The agreement between the empirical collapse and the annealing-schedule-inferred theoretical values of λ(κ) in the inset suggests that kink-density is insensitive to such errors, and conversely demonstrates that the schedule is relatively accurate with respect to the evaluation of (15) at experimental κ.

Given λ(κ) we can consider behavior outside of the collapsed regime, at large t_a. Since the dominant noise is thermal, one should extrapolate using quadratic fits per Eq. (17). The inset in Fig. 7b illustrates this extrapolation for three different values of t_a/λ. The black circles in the main panel show the result of this extrapolation. As in Fig. 1, the extrapolated points agree with the exact \(n\propto {t}_{a}^{-1/2}\) behavior of larger t_a. Thus one can successfully infer noise-free expectation values over a much broader range of annealing time than the ostensible coherent annealing regime.

We can also consider the kink correlation (18). For purposes of this section, we assume the kink-correlation with r normalized to ξ is weakly dependent on t_a (as shown in Figs. 3 and 5) through the narrow range of t_a studied, and we extrapolate to a curve averaged on t_a. Figure 8 shows \({\mathcal{C}}(r)\) as a function of r normalized to the correlation length ξ = 1/n. Figure 8(left) illustrates the result obtained from simulation (BdG and DMRG are in agreement, as described in section “Simulation of dynamics”). Data is fitted with a 10th-order polynomial over the plot-range. Figure 8(center-right) illustrates the ZNE result, showing a significant improvement in alignment with the exact theoretical predictions (depicted by the solid lines), compared to unmitigated results (κ = 1, Fig. 8(center-left)). \({\mathcal{C}}(r)\) calculated by extrapolation, depicted in Fig. 8(center-right) exhibits a similar systematic deviation to that shown in Fig. 3, perhaps consistent with a reduced phase-scrambling-length scale (reduced post-critical dynamics). In Fig. 8(right), a few examples of extrapolation at different values of r are displayed, and unlike the case of thermal noise (Fig. 7), linear extrapolation works well, consistent with the correction of non-thermal noise sources.

GA Advantage2 prototype methodology at a scale L = 1024

In this section, we present details of the experimental method applied to the GA Advantage2_prototype2.3 QPU with supporting data. Code implementations are presented based on open-source repositories, principally the python package dwave-ocean-sdk (code examples rely upon dwave-ocean-sdk, an open-source Python repository hosted at https://github.com/dwavesystems/dwave-ocean-sdk)^48,49. A connection to a QPU in the Leap quantum cloud service is established by DWaveSampler() and can be configured to call any GA QPU: the QPU used in these experiments is an Advantage2 prototype. The 1024 spin-variables i must be mapped to QPU-specific qubits for programming purposes via an embedding. We require a set of qubits programmable as a ring, which amounts to solving a subgraph isomorphism problem⁴⁹. This can be done with open-source tools typically in less than a second as follows Models with L > 1024 can be embedded, as above using additional time, or exploiting graph-specific insights. Embedding of smaller graphs is quicker, and can allow for parallelization. The ability to sample embeddings is useful to mitigate any locality-dependent errors. This can be achieved by creating the source and/or target graphs with shuffled variable and/or edge ordering before calling the deterministic subgraph search. Standard errors presented are determined with respect to variance between estimators on different programmings.

Given an anneal duration target t_a (in microseconds, for κ = 1), and a rescaling κ, we can program and sample the QPU as lam(κ) is the function (15), calculated from the published schedule.

For the kink-density experiments, we drew samples for five embeddings, a total of 5 programming, and 5000 samples. For the kink-correlation experiments we drew 16,000 samples for the same 5 embeddings, a total of 80 programmings and 80,000 samples. The kink density (inverse correlation length) fluctuates slightly as a function of κ, and the zero-noise extrapolated value differs from the simulation result. This can be accounted for by small errors in the schedule; the empirical method of the the section “Experimental details—Experimental Advantage2 prototype methodology at a scale L = 232” by which curves are collapsed can mitigate for such misspecification in principle. In Fig. 3, we calculate C(r) for fixed (programmed) t_a at each κ, using the empirical kink correlations, and kink density (which is non-constant in κ).

Calibration refinement method

A refinement of the calibration, dependent on programmed Hamiltonian, is possible through exploiting symmetries. Coded examples of the method have been published⁴⁷. Specifically, the transverse field Ising model is translationally invariant, and parity (defined by the operator \({\prod }_{i}{\sigma }_{i}^{x}\)) is a conserved quantity. As a consequence \({C}_{i}=\langle {\sigma }_{i}^{z}{\sigma }_{i+1}^{z}\rangle =\bar{C}=1-2n\), and \({m}_{i}=\langle {\sigma }_{i}^{z}\rangle =0\), respectively for all i in the noiseless limit. Deviations are a calibration failure that can be corrected by adjustment of the programmed J and flux bias (ϕ) to minimize sum square residuals: \({\mathcal{O}}(\phi )={\sum }_{i}{m}_{i}^{2}\) and \({\mathcal{O}}(J)={\sum }_{i}{({C}_{i}-\bar{C})}^{2}\). Stochastic local search is sufficient, given that the residuals are smooth functions of the calibration parameters about the correct calibration:

Here, {m, C, n} are sampling-based estimates on the t^th programming, and α(t) is a learning rate.

For the presented data, we specify the learning rate as follows. In the first T/2 iterations, the learning rate is heuristically adapted so that the objective approaches a static value controlled by sampling error. On iterations for which \({{\mathcal{O}}}_{\phi }\) (\({{\mathcal{O}}}_{J}\)) grows, the learning rate is decreased by a factor 1.1; otherwise it grows by the same factor. A fixed value of α_J ≲ 1, or α_ϕ ~ 10μΦ₀, can also work well, but adaptation is useful for efficiency and stability, particularly at higher susceptibility parameterizations (lower kink rate). In the second T/2 iterations the learning rate decays polynomially as ~ 1/t, so that we mitigate for sampling noise, and obtain theoretical guarantees on convergence (up to non-static, but low-frequency variation beyond the experimental time scale, which is practically small).

The calibration is found to be converged for our purposes at T = 64. Given the refined calibration, we collect data per section “Experimental details—GA Advantage2 prototype methodology at a scale L = 1024”. In Fig. 9 we show the impacts on kink density and kink correlation. The impacts on the kink density are small: a modest decrease in the kink density is observed. The kink-correlation peak value is enhanced: this is consistent with reduced dephasing or control errors^18,43. Whilst the result is improved, the calibration refinement does not substantially change the conclusions relative to experiments performed with the baseline calibration.

Similarly, we can use the method of the section “Experimental details—Experimental Advantage2 prototype methodology at a scale L = 232” to collapse data with matched ξ = 1/n. We can take the collected data, and interpolate C(r) so that kink density is matched at t_a = 7ns, and extrapolate this collapsed data, as shown in Fig. 9. We combine this method with the calibration refinement, as was used on the experimental Advantage2 system. The result is a further enhancement of the peak C(r) height, which is compatible with reduced noise. However, the agreement with a coarsening theory based only on l/ξ is not improved (see departure from the fit curve). The data clearly shows signs of clustering in κ, indicating the linear (perturbative) fitting method has failed. Failure of the linear extrapolation ZNE method to match the coarsening theory at small t_a seems to occur in conjunction with C(r) data clustering (as a function of κ), the origin of which requires further study.

Simulation of dynamics

Simulation results are obtained both by integration of the Bugolioubov-de Gennes equations (73), and matrix product state (MPS) dynamics. There is agreement between these methods. Results of simulation through the experimental range are shown in Fig. 10, and in other experiments for the GA Advantage2 system used in experiments. Results are well described by the analysis of IV C. In the experimental range the kink density is described by Kibble-Zurek theory; the kink correlation can be understood after accounting for post-critical dynamics with an additional length scale.

Simulation results for the kink correlation are shown for the GA Advantage2 system through relevant experimental parameters as a function of annealing time (top) and problem Hamiltonian rescaling (bottom). The peak value and position of the kink correlation curve grow larger with t_a, and to a lesser extent with κ. The curves are well characterized by a kink density (determined by the Kibble-Zurek theory), and a phase-scrambling-length scale l, as described in section “Critical and post-critical dynamics of the transverse field Ising chain”; best fit values for l are shown. A two-length scale model predicts slightly smaller peak values than the simulation, but this discrepancy is consistent with corrections of order O(1/ξ, 1/L), and corrections vanish at larger scales as expected.

MPS dynamics are here simulated by the time-dependent variational principle (TDVP) with a two-site update^50,51 using ITensor library⁵². We also compared the results against the W_II method introduced in ref. ⁵³ using the TenPy library⁵⁴ and the local Krylov method⁵⁵ (for details about this method in MPS language see ref. ⁵⁶) implemented by the DMRG++ library⁵⁷. Overall, we found that TDVP emerges as the most efficient method, providing converged results using the largest time step (dt=0.01 ns). For the lattice geometry, we adopted the same mapping from periodic to open chain used in ref. ¹⁸ as shown in Fig. 11a. We found a much better performance of the numerical simulations when using a local Hilbert space of two qubits (see Fig. 11b) for the MPS wave-functions.

a Site ordering for MPS simulations of a periodic chain with one-qubit local Hilbert space. b Site ordering with a two-qubit local Hilbert space.

For the results shown in the main part of the manuscript, TDVP simulations were performed with a bond dimension D = 32 and a time step of 0.01 ns, with a maximum truncation error of 10⁻¹⁰.

Results can also be obtained, without the complications of a bond-dimension parameter, by a simple first- or second-order (Euler or Runge-Kutta) integration of the Bugoliobov de-Gennes equations at step size dt ~ 0.01 ns, as described in the section “Critical and post-critical dynamics of the transverse field Ising chain”. A DMRG method is also implemented validating the implementation. A strength of DMRG is the ability to efficiently generalize for inclusion of fields (h_i), beyond first-neighbor couplings, or even explicit modeling of a bath, as analyzed in early works¹⁸.