Consistency tests of classical and quantum models for a quantum annealer

Recently the question of whether the D-Wave processors exhibit large-scale quantum behavior or can be described by a classical model has attracted significant interest. In this work we address this question by studying a 503-qubit D-Wave Two device in the “black box” model i.e., by studying its input-output behavior. Our work generalizes an approach introduced in Boixo et al. [Nat. Commun. 4, 2067 (2013)] and uses groups of up to 20 qubits to realize a transverse Ising model evolution with a ground-state degeneracy whose distribution acts as a sensitive probe that distinguishes classical and quantum models for the D-Wave device. Our findings rule out all classical models proposed to date for the device and provide evidence that an open-system quantum dynamical description of the device that starts from a quantized energy level structure is well justified, even in the presence of relevant thermal excitations and a small value of the ratio of the single-qubit decoherence time to the annealing time.

Figure 1

DW2 annealing schedules $A\left(t\right)$ and $B\left(t\right)$ along with the operating temperature of $T=17\mathrm{mK}$ (black dashed horizontal line). The large $A\left(0\right)/\left(k_{B}T\right)$ value ensures that the initial state is the ground state of the transverse-field Hamiltonian. The large $B\left(t_f\right)/\left(k_{B}T\right)$ value ensures that thermal excitations are suppressed and the final state reached is stable. Also shown are the attenuated $\alpha B\left(t\right)$ curves for (a) the value of $\alpha$ at which the intersection between $A\left(t\right)$ and $\alpha B\left(t\right)$ coincides with the operating temperature (blue dot-dashed curve) and (b) the largest $\alpha$ such that $\alpha B\left(t\right)$ remains below the temperature line for the entire evolution (blue dotted curve).

Figure 2

The eight-spin Ising Hamiltonian. The inner “core” spins (green circles) have local fields $h_i=+1$ [using the convention in Eq. (2)], while the outer spins (red circles) have $h_i=-1$. All couplings are ferromagnetic: $J_{ij}=1$ (black lines).

Figure 3

Schematic representation of the 12, 16, and 20 spin Hamiltonians used in our tests. Extensions to larger $N$ follow the same pattern, with $N/2$ qubits in the inner ring and $N/2$ qubits in the outer ring. Notation conventions are as in Fig. 2.

Figure 4

Numerically calculated evolution of the gap between the instantaneous ground state and the 17th excited state (which becomes the first excited state at $t=t_f)$, for the eight-spin Hamiltonian in Eq. (1), following the annealing schedule of the DW2 device (Fig. 1). The gap value is shown for some interesting values of $\alpha$ (see Fig. 1). The kinks are due to energy level crossings, as explained in Appendix pp4. A reduction in $\alpha$ results in a reduction of the size of the minimal gap and delays its appearance.

Figure 5

ME simulation for the time dependence of the probability of being in the lowest 17 energy eigenstates, for different values of $\alpha$. Simulation parameters are $t_f=20\mu \mathrm{s}$ (the minimal annealing time of the DW2) and $\kappa =1.27\times {10}^{-4}$, where $\kappa$ is an effective, dimensionless system-bath coupling strength defined in Appendix pp5-s3. The chosen value of $\kappa$ allows us to reliably probe the small $\alpha$ regime.

Figure 6

Distribution of the ground states for $N=8$ for (a) ME with no noise on $\{h_i,J_{ij}\}$, (b) SSSV with no noise on $\{h_i,J_{ij}\}$, and (c) SSSV with $\{h_i,J_{ij}\}$ noise using $\sigma =0.085$. The cluster states are labeled by their Hamming distance $H$ from the isolated state and by their multiplicity $M$ for a given value of $H$. The vertical axis is the final probability $p$ of a given $(H,M)$ set, divided by its multiplicity and the total ground-state probability. The data symbols $(\circ$, etc.) are the mean values of the bootstrapped [31] distributions, and the error bars are two standard deviations below and above the mean representing the 95% confidence interval. Note that the SSSV model prefers the $|11110000\rangle$ cluster state, whereas the ME gives a uniform distribution over all cluster states. SSSV parameters are $T=10.56\mathrm{mK}$ and $1\times {10}^5$ Monte Carlo step updates per spin (“sweeps”). The same parameters are used in all subsequent SSSV figures. These results do not include the cross-talk correction.

Figure 7

Time evolution (according to the ME) of energy eigenstate populations for $\alpha =0.1$ and $\kappa =8.9\times {10}^{-4}$ (this relatively large value was chosen here since it results in increased thermal excitation/relaxation). $P_i$ denotes the population of the $i\mathrm{th}$ eigenstate, with $i=1$ being the instantaneous ground state. The energy eigenstate that eventually becomes the isolated ground state is $i=6$ (dashed red line). This state acquires more population at the end of the evolution than the other 16 eigenstates that eventually become the cluster (solid purple line). (Inset) The difference of the population ratio between the open-system and the closed-system evolution, $\Delta (P_{\mathrm{I}}/P_{\mathrm{C}})={(P_{\mathrm{I}}/P_{\mathrm{C}})}_{\mathrm{open}}-{(P_{\mathrm{I}}/P_{\mathrm{C}})}_{\mathrm{closed}}$. The deviation from closed-system dynamics starts at $t/t_f\approx 0.4$, when the $i=6$ eigenstate becomes thermally populated at the expense of the lowest five eigenstates.

Figure 8

Statistical box plot of the probability divided by the multiplicity of being in a given state with Hamming distance $H$ and multiplicity $M$. Shown are the ME isolated state and cluster states for $N=8$, with 512 noise realizations applied to the $h$'s and $J$'s [with distribution $\mathcal{N}(0,0.06)]$ at $\alpha =1$. The isolated state $\left(H=0,M=1\right)$ is suppressed, while the cluster states are, on average, equally populated. The red bar is the median, the blue box corresponds to the lower and upper quartiles, respectively, the segment contains most of the samples, and the $+$'s are outliers [32]. The horizontal axis label indicates the Hamming distance from the isolated state and the multiplicity of the cluster states at each value of $H$. States that are equivalent up to ${90}^{\circ }$ rotations are grouped together. For example, there are four rotationally equivalent cluster states that have two adjacent outer qubits pointing down, while the other two are pointing up. Only the $H=6$ case splits into two rotationally inequivalent sets.

Figure 9

(a) Ground-state populations for DW2, $N=8$ . Legend: $(H,M)$, corresponding to Hamming distance from the isolated state and multiplicity, respectively. Error bars represent the 95% confidence interval. (b) Cluster-state populations for the ME, $\chi =0.015,N=8$. Solid lines correspond to the results with no noise on the $\{h_i,J_{ij}\}$'s, while the data points include Gaussian noise with mean 0 and standard deviation $\sigma =0.025$ for 100 noise realizations. The error bars represent the 95% confidence interval. The DW2 data from (a) is also plotted as the shaded region representing a 95% confidence interval with the dashed lines corresponding to the mean. The inset shows the behavior for the noiseless ME for small $\alpha$. (c) Cluster-state population for SSSV for $\chi =0.035,\sigma =0.085,N=8$ with the DW2 data plotted as in (b). In contrast to Fig. 6, both the ME and the noisy SSSV model include the cross-talk correction [Eq. (6)], with $\chi$ chosen to optimize the fit for the cluster-state populations at $\alpha =1$. (d) Only the cluster states with Hamming distance 4 and 8 from the isolated state are shown for DW2, the ME $\left(\chi =0.015\right)$, SSSV $\left(\chi =0.035,\sigma =0.085\right)$, $N=8$, and SSSV from panels (b) and (c) in order to highlight their differences. Panel (e) displays the same for $N=20$ (excluding the ME, which is too costly to simulate at this scale). (f) The isolated-state populations for DW2 and SSSV $\left(\chi =0.035,\sigma =0.085\right)$, ME $\left(\chi =0.015\right)$, and noisy ME $\left(\chi =0.015,\sigma =0.025\right)$, $N=8$, which highlights the qualitative agreement between the models and DW2. Experimental data were collected using the in-cell embeddings strategy described in Appendix pp3. The embedding and gauge-averaging strategies are also discussed in Appendix pp3. The color coding of states is consistent across all panels.

Figure 10

Subset of the first excited-state populations for (a) DW2 for $N=8$; (b) SSSV for $\chi =0.035,\sigma =0.085,N=8$; and (c) ME for $\chi =0.015,\sigma =0.025,N=8$. In (b) and (c), the simulations include qubit cross-talk correction with $\chi$ chosen as in Fig. 9 to optimize the fit for the cluster-state populations at $\alpha =1$. Panels (d) and (e) are for $N=20$ and $\chi =0.035,\sigma =0.085,N=20$, respectively. The $\Pi$ symbol denotes all permutations. The SSSV model does not reproduce the correct ordering. The error bars represent the 95% confidence interval.

Figure 11

Time-dependence of the negativity [Eq. (9)] for (a) a closed-system evolution and (b) an open-system evolution of $N=8$ qubits (modeled via the ME with $\kappa =1.27\times {10}^{-4})$, as a function of $\alpha$. The rapid decay of negativity for small $\alpha$ in the open-system case signals a transition to classicality. However, for large $\alpha$ the closed- and open-system negativities are similar, suggesting that the system is quantum in this regime. The apparent jaggedness of the closed-system plot near $\alpha =0$ is due to our discretization of $\alpha$ in steps of 0.01.

Figure 12

Qubits and couplers in the DW2 device. The DW2 “Vesuvius” chip consists of an $8\times 8$ two-dimensional square lattice of eight-qubit unit cells with open boundary conditions. The qubits are each denoted by circles, connected by programmable inductive couplers as shown by the lines between the qubits. Of the 512 qubits of the device located at the University of Southern California used in this work, the 503 qubits marked in green and the couplers connecting them are functional.

Figure 13

Embedding according to the random parallel embeddings strategy for eight spins. An example of 15 randomly generated different parallel copies of the eight-spin Hamiltonian. Our data collection used a similar embedding with 50 different copies.

Figure 14

Embedding according to the random parallel embeddings strategy for 16 spins. An example of 10 randomly generated different parallel copies of the 16-spin Hamiltonian. Our data collection used a similar embedding with 93 different copies.

Figure 15

Embedding according to the “in-cell embeddings” strategy. An example of a randomly generated in-cell embedding.

Figure 16

Some representative autocorrelation tests at $t_f=20\mu \mathrm{s}$ showing the standard error of the mean $\Delta x$ as a function of binning size, for different values of $\alpha$: (a) $\alpha =0.1$, (b) $\alpha =0.35$, (c) $\alpha =1$. Each curve is the result of the binning test for a different state. The relatively flat lines for all states suggest that there are no significant autocorrelations in the data.

Figure 17

(a) Time dependence of the gap between the 18th excited state and the instantaneous ground state for different values of the energy scale factor $\alpha$. (b)–(d) Time dependence of the lowest 56 energy eigenvalues for different values of $\alpha$ [(b) $\alpha =1/7$, (c) $\alpha =2/7$, (d) $\alpha =3/7]$. Note that the identity of the lowest 17 energy eigenvalues changes over the course of the evolution.

Figure 18

Simulated annealing shows quantitatively similar behavior for various annealing schedule. The schedules are exponential $T\left(k\right)=T\left(0\right){\left[T\left(K\right)/T\left(0\right)\right]}^{k/K}$, linear $T\left(k\right)=T\left(0\right)+\frac{k}{K}[T\left(K\right)-T\left(0\right)]$, and constant $T\left(k\right)=T\left(K\right)$, with $K=1000$, with $k_{B}T\left(0\right)/\hslash =8$ GHz and $k_{B}T\left(K\right)/\hslash =0.5$ GHz.

Figure 19

Numerical results distinguishing the quantum ME and classical SA, SD, and SSSV models. (a) Results for the ratio of the isolated-state population to the average population in the cluster states $\left(P_{\mathrm{I}}/P_{\mathrm{C}}\right)$, and (b) the ground-state probability $\left(P_{\mathrm{GS}}\right)$, as a function of the energy scale factor $\alpha$, at a fixed annealing time of $t_f=20\mu \mathrm{s}$. The error bars represent the 95% confidence interval. Two striking features are the “ground-state population inversion” between the isolated state and the cluster (the ratio of their populations crosses unity), and the manifestly nonmonotonic behavior of the population ratio, which displays a maximum. At the specific value of the system-bath coupling used in our simulations $\left(\kappa =1.27\times {10}^{-4}\right)$, it is interesting that the ME underestimates the magnitude and position of the peak in $P_{\mathrm{I}}/P_{\mathrm{C}}$ but qualitatively matches the experimental results shown in Fig. 9, capturing both the population inversion and the presence of a maximum even in the absence of noise. In contrast to the ME results shown, the SA, SSSV, and SD results for the population ratio are not in qualitative agreement with the ME. Specifically, all three classical models miss the population inversion and maximum seen for the ME. Simulation parameters can be found in Appendix pp5.

Figure 20

Statistical box plot of the average $z$ component for all core qubits (main plot) and all outer qubits (inset) at $t=t_f=20\mu \mathrm{s}$. (a) The SD model. The data are taken for 1000 runs with Langevin parameters $k_{B}T/\hslash =2.226\mathrm{GHz}$ (i.e., 17 mK, to match the operating temperature of the DW2) and $\zeta ={10}^{-3}$. In Appendix pp10 we show that the results do not depend strongly on the choice of $\zeta$. (b) The SSSV model. The data are taken for 1000 runs with parameters $k_{B}T/\hslash =1.382\mathrm{GHz}$ (i.e., $T=10.56$ mK, as in Ref. [11]) and $5\times {10}^5$ sweeps.

Figure 21

“Forced” and strongly or weakly “decohering” SSSV models. Shown are the results for the ratio of the isolated-state population to the average population in the cluster states $\left(P_{\mathrm{I}}/P_{\mathrm{C}}\right)$ as a function of the energy scale factor $\alpha$, for $N=8$ and at a fixed annealing time of $t_f=20\mu \mathrm{s}$. The error bars represent the 95% confidence interval. For reference the plot also includes the curves for the ME from Fig. 19. Additional parameters for the modified SSSV models: $g^2\eta ={10}^{-6}$ for the strongly decohered model and $g^2\eta =2.5\times {10}^{-7}$ for the weakly decohered and forced models.

Figure 22

Trace-norm distance of the ME, SA, SD, SSSV, weakly and strongly decohered SSSV, and forced SSSV states from the $T=17$ mK Gibbs state at $t_f=20\mu \mathrm{s}$ and $N=8$. The error bars represent the 95% confidence interval. Three regions are clearly distinguishable for the ME: (1) $1\ge \alpha \gtrsim 0.3$, where $\mathcal{D}$ is decreasing as $\alpha$ decreases; (2) $0.3\gtrsim \alpha \gtrsim 0.1$, where $\mathcal{D}$ is increasing as $\alpha$ decreases; (3) $0.1\gtrsim \alpha \ge 0$, where $\mathcal{D}$ is again decreasing as $\alpha$ decreases. Both SA and SD lack the minimum at $\alpha \approx 0.3$.

Figure 23

Master-equation results for the populations of the 17 Ising ground states, with $\alpha =1,\left|h\right|=0.981\left|J\right|,t_f=20\mu \mathrm{s}$, and $\kappa =1.27\times {10}^{-4}$. The cluster states split by HD from the isolated state (bottom curve), in agreement with the experimental results shown in Fig. 24.

Figure 24

Statistical box plot of the gauge-averaged ground-states population (a),(c) before and (b),(d) after optimization of $J_{ij}$ as per Table 1, for (a),(b) $N=8$ and (c),(d) $N=16,\alpha =1$ and $t_f=20\mu \mathrm{s}$. Only the $N=8,H=6$ case splits into two rotationally inequivalent sets. Note the clear step structure in the cluster states $\left(H>0\right)$ in (a),(c), while in (b),(d) the population of the cluster states is fairly equalized (less so in the $N=16$ case since Table 1 is optimized for $N=8)$. (a) Data taken with the random parallel embeddings strategy. (b) Data taken using the in-cell embeddings strategy, with the optimized values of the couplings given in Table 1. The same optimization removes the step structure from data taken with the random parallel embeddings strategy (not shown).

Figure 25

Ratio of isolated-state population to average cluster-state population as a function of the energy scale factor $\alpha$, for $t_f=20\mu \mathrm{s}$ and $N=8$. Shown are the ratios calculated with both uncorrected and corrected values of $J$ (as per Table 1), the latter tuned to flatten the steps seen in the population of the cluster states. Error bars represent the standard error of the mean value of the ratio estimated using bootstrapping.

Figure 26

(a) Noisy SSSV, offset $h=0.97J$; (b) ME, offset $h=0.985J$; (c) perturbation theory [Eq. (H2)] offset $h=0.975J$; with a population ordering correction (offset of $h$ vs $J)$ at $\alpha =1$. The error bars represent the 95% confidence interval. The SSSV model now has the right ordering of the cluster states but clearly disagrees with the DW2 result [Fig. 9] for $\alpha \lesssim 0.3$, near where the isolated state has its maximum. The ME result is in qualitative agreement with the DW2 result except that the cluster-state populations do not equalize for small $\alpha$, which is a consequence of not including the $\alpha$ dependence of the offset.

Figure 27

(a) DW2, offset $h=0.981J$; and (b) noisy SSSV, offset $h=0.988J$, model with a population equalizing correction at $\alpha =1$. (c) Noisy SSSV model with offset $h=0.94J$ chosen to equalize the cluster-state populations at $\alpha =0.2$. The error bars represent the 95% confidence interval.

Figure 28

Subset of the first excited-state populations for (a) perturbation theory [as in Eq. (H2)] offset $h=0.975J$, (b) offset $h=0.97J$, (c) offset $h=0.988J$, and (d) SSSV offset $h=0.94J$. Panel (f) shows the case with qubit cross-talk discussed in Sec. 3. The $\Pi$ symbol denotes all permutations. Whereas the perturbation theory result for a QA Hamiltonian [Eq. (H2)] reproduces the correct ordering, none of the three SSSV cases shown does. These three SSSV cases were chosen to optimize the fit for the cluster-state populations. The error bars represent the 95% confidence interval.

Figure 29

Ratio of the isolated-state population to the average population in the cluster $\left(P_{\mathrm{I}}/P_{\mathrm{C}}\right)$ as a function of the energy scale factor $\alpha$, for two different values of $t_f$, and $N=8$. The inset shows the ME results. Data were collected using the “in-cell embeddings” strategy (see Appendix pp3 for details). Error bars are one standard deviation above and below the mean.

Figure 30

Numerically calculated instantaneous energy gap between the ground and first excited states for the 8, 12, and 16 spin Hamiltonians. The gap vanishes since the first excited state becomes part of the $2^{N/2}+1$-fold degenerate ground-state manifold at $t=t_f$.

Figure 31

Ratio of the isolated-state population to the average population in the cluster states $\left(P_{\mathrm{I}}/P_{\mathrm{C}}\right)$ as a function of the energy scale factor $\alpha$, for different values of $N$, at a fixed annealing time of $t_f=20\mu \mathrm{s}$. The nonmonotonic dependence of the population ratio on $\alpha$ is observed for all values of $N$. The growth of the $P_{\mathrm{I}}/P_{\mathrm{C}}$ peak with increasing $N$ is consistent with the discussion presented in Appendix pp9. The increasingly large error bars are due to the smaller amount of data collected as $N$ grows. For $N=12,16,20$ data were collected using the “random parallel embeddings” strategy and for $N=40$ using the “designed parallel embedding” strategy (see Appendix pp3 for details). Error bars are one standard deviation above and below the mean.

Figure 32

Evolution of a core (blue) and outer (green) spin with $t_f=20\mu s$, subject to the O(3) SD model with $\alpha =1$. All spins start with $M^x=1,M^y=M^z=0$, i.e., point in the $x$ direction. (a) Closed-system case given by Eq. (J17). (b) Open-system case given by Eq. (J19). Rapid oscillations at the beginning of the evolution in (a) are because the initial conditions used are not exactly the ground state of the system [because of the finite $B\left(0\right)]$. In (b), Langevin parameters are $k_{B}T/\hslash =2.226\mathrm{GHz}$ and $\zeta ={10}^{-6}$.

Figure 33

Distribution of $M^z$ at the end of the evolution for all (a) core spins and (b) all outer spins for $\alpha =1$. Langevin parameters are $k_{B}T/\hslash =2.226\mathrm{GHz}$ and $\zeta ={10}^{-3}$. Data collected using 1000 runs of Eq. (J19).

Figure 34

Statistical box plot of the $z$ component for all core qubits (main plot) and all outer qubits (inset) at $t=t_f=20\mu \mathrm{s}$. The data is taken for 1000 runs of Eq. (J19) with Langevin parameters $k_{B}T/\hslash =2.226\mathrm{GHz}$ (to match the operating temperature of the DW2) and (a) $\zeta ={10}^{-1}$ and (b) $\zeta ={10}^{-5}$. The $\zeta ={10}^{-3}$ is shown in Fig. 20. This illustrates that the results do not depend strongly on the choice of $\zeta$.

Figure 35

Distribution of $M^z$ at the end of the evolution for (a) all core spins and (b) all outer spins (i.e., the values of all core spins and all outer spins are included in each respective box plot) for $\alpha =2/7$. Langevin parameters are $k_{B}T/\hslash =2.226\mathrm{GHz}$ and $\zeta ={10}^{-3}$. Data collected using 1000 runs of Eq. (J19). Note that the $t_f$-axis scale is not linear.