Entanglement in a quantum annealing processor

Entanglement lies at the core of quantum algorithms designed to solve problems that are intractable by classical approaches. One such algorithm, quantum annealing (QA), provides a promising path to a practical quantum processor. We have built a series of scalable QA processors consisting of networks of manufactured interacting spins (qubits). Here, we use qubit tunneling spectroscopy to measure the energy eigenspectrum of two- and eight-qubit systems within one such processor, demonstrating quantum coherence in these systems. We present experimental evidence that, during a critical portion of QA, the qubits become entangled and that entanglement persists even as these systems reach equilibrium with a thermal environment. Our results provide an encouraging sign that QA is a viable technology for large-scale quantum computing.


I. INTRODUCTION
The past decade has been exciting for the field of quantum computation. A wide range of physical implementations of architectures that promise to harness quantum mechanics to perform computation have been studied [1][2][3]. Scaling these architectures to build practical processors with many millions to billions of qubits will be challenging [4,5]. A simpler architecture, designed to implement a single quantum algorithm such as quantum annealing (QA), provides a more practical approach in the near term [6,7]. However, one of the main features that makes such an architecture scalable, namely, a limited number of low-bandwidth external control lines [8], prohibits many typical characterization measurements used in studying prototype universal quantum computers [9][10][11][12][13][14]. These constraints make it challenging to experimentally determine whether a scalable QA architecture, one that is inevitably coupled to a thermal environment, is capable of generating entangled states [15][16][17][18]. A demonstration of entanglement is considered to be a critical milestone for any approach to building a quantum computing technology. Herein, we demonstrate an experimental method to detect entanglement in subsections of a quantum annealing processor to address this fundamental question.

II. QUANTUM ANNEALING
QA is designed to find the low-energy configurations of systems of interacting spins. A wide variety of optimization problems naturally map onto this physical system [19][20][21][22]. A QA algorithm is described by a time-dependent Hamiltonian for a set of N spins, i ¼ 1; …; N, where the dimensionless H P is and σ x;z i are Pauli matrices for the ith spin. The energy scales Δ and E are the transverse and longitudinal energies of the spins, respectively, and the biases h i and couplings J ij encode a particular optimization problem. The timedependent variation of Δ and E is parametrized by s ≡ t=t f with time t ∈ ½0; t f and total run (anneal) time t f . QA is performed by first setting Δ ≫ E, which results in a ground state into which the spins can be easily initialized [6]. Then Δ is reduced and E is increased until E ≫ Δ. At this point, the system Hamiltonian is dominated by H P , which represents the encoded optimization problem. At the end of the evolution, a ground state of H P represents the lowest energy configuration for the problem Hamiltonian and thus a solution to the optimization problem.

III. QUANTUM ANNEALING PROCESSOR
We have built a processor that implements H S using superconducting flux qubits as effective spins [6,7,23,24]. Figure 1(a) shows a photograph of the processor. Figure 1(c) shows the circuit schematic of a pair of flux qubits with the magnetic flux controls Φ x qi and Φ x ccjj . The annealing parameter s is controlled with the global bias Φ x ccjj ðtÞ (see Appendix A for the mapping between s and Φ x ccjj and a description of how Φ x qi is provided for each qubit). The strength and sign of the inductive coupling between pairs of qubits is controlled with magnetic flux Φ x co;ij that is provided by an individual on-chip digital-to-analog converter for each coupler [8]. The parameters h i and J ij are thus in situ tunable, thereby allowing the encoding of a vast number of problems. The time-dependent energy scales ΔðsÞ and EðsÞ are calculated from measured qubit parameters and plotted in Fig. 1(d). We calibrate and correct the individual flux qubit parameters in our processor to ensure that every qubit has a close-to-identical Δ and E (the energy gap Δ is balanced to better than 8% between qubits and E to better than 5%). See Appendix A for measurements of these energy scales. The interqubit couplers are calibrated as described in Ref. [25]. The processor we study here is mounted on the mixing chamber of a dilution refrigerator held at temperature T ¼ 12.5 mK.

IV. FERROMAGNETICALLY COUPLED INSTANCES
The experiments reported herein focus on one of the eight-qubit unit cells of the larger QA processor as indicated in Fig. 1(a). The unit cell is isolated by setting all couplings outside of that subsection to J ij ¼ 0 for all experiments. We then pose specific H P instances with strong ferromagnetic (FM) coupling J ij ¼ −2.5 and h i ¼ 0 to that unit cell, as illustrated in Figs. 1(e) and 1(f). These configurations produce coupled two-and eight-qubit systems, respectively. The Hamiltonian (1) describes the behavior of these systems during QA.
Typical observations of entanglement in the quantum computing literature involve applying interactions between FIG. 1. (a) Photograph of the QA processor used in this study. We report measurements performed on the eight-qubit unit cell indicated. The bodies of the qubits are extended loops of Nb wiring (highlighted with red rectangles). Interqubit couplers are located at the intersections of the qubit bodies. (b) Electron micrograph showing the cross section of a typical portion of the processor circuitry (described in more detail in Appendix A). (c) Schematic diagram of a pair of coupled superconducting flux qubits with external control biases Φ x qi and Φ x ccjj and with flux through the body of the ith qubit denoted as Φ qi . An inductive coupling between the qubits is tuned with the bias Φ x co;ij . (d) Energy scales ΔðsÞ and EðsÞ in Hamiltonian (1) calculated from a rf SQUID model based on the median of independently measured device parameters for these eight qubits. See Appendix A for more details. (e),(f) The two-and eight-qubit systems studied were programmed to have the topologies shown. Qubits are represented as gold spheres, and interqubit couplers, set to J ¼ −2.5, are represented as silver lines. qubits, removing these interactions, and then performing measurements. Such an approach is well suited to gatemodel architectures (see, e.g., Ref. [11]). During QA, however, the interaction between qubits is determined by the particular instance of H P , in this case, a strongly ferromagnetic instance, and cannot be removed. In this way, systems of qubits undergoing QA have much more in common with condensed-matter systems, such as quantum magnets, for which interactions cannot be turned off. Indeed, a growing body of recent theoretical and experimental work suggests that entanglement plays a central role in many of the macroscopic properties of condensedmatter systems [26][27][28][29][30][31][32]. Here, we introduce other approaches to quantifying entanglement that are suited to QA processors. We establish experimentally that the two-and eight-qubit systems, comprising macroscopic superconducting flux qubits coupled to a thermal bath at 12.5 mK, become entangled during the QA algorithm.
To illustrate the evolution of the ground state of these instances during QA, a sequence of wave functions for the ground state of the two-qubit system is shown in Fig. 2. A similar sequence could be envisioned for the eight-qubit system. We consider these systems subject to zero biases, h i ¼ 0. For small s, Δ ≫ 2jJ ij jE, and the ground state of the system can be expressed as a product of the ground states of the individual qubits: Fig. 2 Fig. 2(b)]. The state jþi is the maximally entangled Bell [or Greenberger-Horne-Zeilinger (GHZ), for eight qubits] state [17]. As s → 1, the energy gap g between the ground and first excited states approaches g ≡ ðE 2 − E 1 Þ ∝ ΔðsÞ N = ½2jJ ij jEðsÞ N−1 and vanishes as ΔðsÞ → 0 [ Fig. 2(c)]. At some point late in the evolution, g becomes less than k B T, where T characterizes the temperature of the thermal environment to which the qubits are coupled. At this point, we expect the system to evolve into a mixed state of jþi and j−i, and the entanglement will vanish with g for sufficiently long thermalization times. At the end of QA, s ¼ 1, Δ ∼ 0, and Hamiltonian (1) predicts two degenerate and localized ground states, namely, the FM-ordered states j↑…↑i and j↓…↓i.

V. MEASUREMENTS
In order to experimentally verify the change in spectral gap in the two-and eight-qubit systems during QA, we use qubit tunneling spectroscopy (QTS), as described in more detail in Ref. [33] and Appendix B. QTS allows us to measure the eigenspectrum and level occupation of a system during QA by coupling an additional probe qubit to the system. We perform QTS on the two-and eight-qubit systems shown in Figs. 1(e) and 1(f). Figures 3(a) and 3(b) show the measured energy eigenspectrum for the two-and eight-qubit systems, respectively, as a function of s. The measurements are initial tunneling rates of the probe qubit, normalized by the maximum observed tunneling rate. Peaks in the measured tunneling rate map the energy eigenstates of the system under study [33]. As the system evolves (increasing s), ΔðsÞ in Hamiltonian (1) decreases and the gap between ground and first excited states closes. The spectroscopy data in Fig. 3(a) reveal two higher-energy eigenstates. We observe a similar group of higher-energy excited states for the eight-qubit system in Fig. 3(b). Note that g closes earlier in the QA algorithm for the eight-qubit system as compared to the two-qubit system. In all of the panels of Fig. 3, solid curves indicate the theoretical energy levels predicted by Hamiltonian (1) using the measured 2. An illustration of entanglement between two qubits during QA with h i ¼ 0 and J < 0. We plot calculations of the two-qubit ground-state wave-function modulus squared in the basis of Φ q1 and Φ q2 , the flux through the bodies of q 1 and q 2 , respectively. The color scale encodes the probability density with red corresponding to high probability density and blue corresponding to low probability density. We use Hamiltonian (1) and the energies in Fig. 1(d) for the calculation. The four quadrants represent the four possible states of the two-qubit system in the computation basis. We also plot the single-qubit potential energy (U versus Φ q1 ) calculated from measured device parameters. (a) At s ¼ 0 (Δ ≫ 2jJjE ∼ 0), the qubits weakly interact and are each in their ground state ð1= ffiffi ffi 2 p Þðj↑i þ j↓iÞ, which is delocalized in the computation basis. The wave function shows no correlation between q 1 and q 2 and, therefore, their wave functions are separable. (b) At intermediate s (Δ ∼ 2jJjE), the qubits are entangled. The state of one qubit is not separable from the state of the other, as the ground state of the system is approximately jþi ≡ ðj↑↑i þ j↓↓iÞ= ffiffi ffi 2 p . A clear correlation is seen between q 1 and q 2 . (c) As s → 1, Δ ≪ 2jJjE, and the ground state of the system approaches jþi. However, the energy gap g between the ground state (jþi) and the first excited state (j−i) is closing. When the qubits are coupled to a bath with temperature T and g < k B T, the system is in a mixed state of jþi and j−i and entanglement is extinguished.
ΔðsÞ and EðsÞ. The agreement between the experimentally obtained spectrum and the theoretical spectrum is good.
The data presented in Figs. 3(a) and 3(b) indicate that the spectral gap between ground and first excited state decreases monotonically with s when all h i ¼ 0. Under these bias conditions, these systems possess Z 2 symmetry between the states j↑…↑i and j↓…↓i. The degeneracy between these states is lifted by finite ΔðsÞ. To explicitly demonstrate that the spectral gap at h i ¼ 0 is due to the avoided crossing of j↑…↑i and j↓…↓i, we perform QTS at fixed s as a function of a "diagnostic" bias h i ≠ 0 that was uniformly applied to all qubits, thus sweeping the systems through degeneracy at h i ¼ 0. As a result, either the state j↑…↑i or j↓…↓i becomes energetically favored, depending upon the sign of h i . Hamiltonian (1) predicts an avoided crossing, as a function of h i , between the ground and first excited states at degeneracy, where h i ¼ 0, with a minimum energy gap g. The presence of such an avoided crossing is a signature of ground-state entanglement [14,35]. For large gaps, g > k B T, there is persistent entanglement at equilibrium (see Refs. [18,26,28,29,31]).
We experimentally verify the existence of avoided crossings at multiple values of s in both the two-and eight-qubit systems by using QTS across a range of biases h i ∈ f−4; 4g. In Fig. 3(c), we show the measured spectrum of the two-qubit system at s ¼ 0.339 up to an energy of 6 GHz for a range of bias h i . The ground states at the far left and far right of the spectrum are the localized states j↓↓i and j↑↑i, respectively. At h i ¼ 0, we observe an avoided crossing between these two states. We measure an energy gap g at zero bias, h i ¼ 0, between the ground state and the first excited state, g=h ¼ 1.75 AE 0.08 GHz, by fitting a Gaussian profile to the tunneling rate data at these two lowest-energy levels and subtracting the centroids. Here, h (without any subscript) is the Planck constant. Figure 3(d) shows the two-qubit spectrum later in the QA algorithm, at s ¼ 0.351. The energy gap has decreased to g=h ¼ 1.21 AE 0.06 GHz. Note that the error estimates for the energy gaps are derived from the uncertainty in extracting the centroids from the rate data. We discuss the actual source of the underlying Gaussian widths (the observed level broadening) below. For both the two-and eight-qubit Spectroscopic data for two-and eight-qubit systems plotted in false color (color indicates normalized qubit tunnel spectroscopy rates). A nonzero measurement (false color) indicates the presence of an eigenstate of the probed system at a given energy (ordinate) and s (abscissa). Panel (a) shows the measured eigenspectrum for the two-qubit system as a function of s. Panel (b) shows a similar set of measurements for the eight-qubit system. The ground-state energy E 1 has been subtracted from the data to aid in visualization. The solid curves indicate the theoretical expectations for the energy eigenvalues using independently calibrated qubit parameters and Hamiltonian (1). We emphasize that the solid curves are not a fit, but rather a prediction based on Hamiltonian (1) and measurements of Δ and E. The slight differences between the high-energy spectrum prediction and measurements are due to the additional states in the rf SQUID flux qubits. A full rf SQUID model that is in agreement with the measured high-energy spectrum is explored in the Supplemental Material [34]. Panels (c) and (d) show measured eigenspectra of the two-qubit system versus h 1 ¼ h 2 ≡ h i for two values of annealing parameter s, s ¼ 0.339 and s ¼ 0.351 from left to right, respectively. Note the avoided crossing at h i ¼ 0. Panels (e) and (f) show analogous measured eigenspectra for the eight-qubit system (with  EðsÞ. Again, the agreement between the experimentally obtained spectra and the theoretical spectra is good. For the early and intermediate parts of QA, the energy gap g is larger than temperature, g ≫ k B T, for both the two-and eight-qubit systems. We expect that if we hold the systems at these s, then the only eigenstate with significant occupation will be the ground state. We confirm this by using QTS in the limit of long tunneling times to probe the occupation fractions. Details are provided in Appendix C. Figures 4(a) and 4(b) show the measured occupation fractions of the ground and first excited states as a function of s for both the two-and eight-qubit systems. The solid curves show the equilibrium Boltzmann predictions for T ¼ 12.5 mK and are in good agreement with the data.
The width of the measured spectral lines is dominated by the noise of the probe device used to perform QTS [33]. The probe device is operated in a regime in which it is strongly coupled to its environment, whereas the system qubits we study are in the weak coupling regime. The measured line widths, which do not increase with system size, therefore, do not represent the intrinsic width of the two-and eight-qubit energy eigenstates. During the intermediate part of QA, the ground and first excited states are clearly resolved. The ground state is protected by the multiqubit energy gap g ≫ k B T, and these systems are coherent. At the end of the annealing trajectory, the gap between the ground state and first excited state shrinks below the probe qubit line width of 0.4 GHz. An analysis of spectroscopic data, presented in the Supplemental Material [34], shows that the intrinsic energy levels remain distinct until later in QA. The interactions between the two-and eight-qubit systems and their respective environments represent small perturbations to Hamiltonian (1), even in the regime in which entanglement is beginning to fall due to thermal mixing.

VI. ENTANGLEMENT DETECTION
The tunneling spectroscopy data show that, midway through QA, both the two-and eight-qubit systems have avoided crossings with the expected gap g ≫ k B T and have ground-state occupation P 1 ≃ 1. While observation of an avoided crossing is evidence for the presence of an entangled ground state [35], we can make this observation more quantitative with entanglement measures and witnesses.
For a large part of the QA algorithm, the two-and eightqubit systems are in their ground states with high occupation fraction. We, therefore, begin the analysis with a susceptibility-based witness W χ , which detects groundstate entanglement. This witness does not require explicit knowledge of Hamiltonian (1), but requires a nondegenerate ground state, confirmed with the avoided crossings shown in Fig. 3, and high occupation fraction of the ground state, confirmed early in QA by the measurements of P 1 ≃ 1 shown in Fig. 4. We then perform measurements of all available linear cross susceptibilities where hσ z i i is the expectation value of σ z i for the ith qubit andh j ¼ Eh j is a bias applied to the jth qubit. The measurements are performed at the degeneracy point (in the middle of the avoided crossings), where the classical contribution to the cross susceptiblity is zero.
From these measurements, we calculate W χ as defined in Ref. [35] (see Appendix D for more details). A nonzero value of this witness detects ground-state entanglement, and global entanglement in the case of the eight-qubit system (meaning every possible bipartition of the eightqubit system is entangled). Figures 4(c) and 4(d) show W χ for the two-and eight-qubit systems. Note that for two qubits at degeneracy, W χ coincides with ground-state concurrence. These results indicate that the two-and eight-qubit systems are entangled midway through QA. Note also that a susceptibility-based witness has a close analogy to susceptibility-based measurements of nanomagnetic systems that also report strong nonclassical correlations [29,31].
The occupation fraction measurements shown in Fig. 4 indicate that midway through QA, the first excited state of these systems is occupied as the energy gap g begins to approach k B T. The systems are no longer in the ground state but, rather, in a mixed state. To detect the presence of mixed-state entanglement, we need knowledge about the density matrix of these systems. Occupation fraction measurements provide measurements of the diagonal elements of the density matrix in the energy basis. We assume that the density matrix has no off-diagonal elements in the energy basis (they decay on time scales of several nanoseconds). We relax this assumption below. Populations P 1 and P 2 plotted in Figs. 4(a) and 4(b) indicate that the system occupies these states with almost 100% probability. This means that the density matrix can be written in the form ρ ¼ P 2 i¼1 P i jψ i ihψ i j, where jψ i i represents the ith eigenstate of Hamiltonian (1).
We use the density matrix to calculate standard entanglement measures, Wootters' concurrence C [18] for the two-qubit system, and negativity N [16,36] for the two-and eight-qubit system. For the maximally entangled two-qubit Bell state, we note that C ¼ 1 and N ¼ 0.5. Figure 4(c) shows C as a function of s. Midway through QA, we measure a peak concurrence C ¼ 0.53 AE 0.05, indicating significant entanglement in the two-qubit system. This value of C corresponds to an entanglement of formation E f ¼ 0.388 (see Refs. [16,18] for definitions). This is comparable to the level of entanglement, E f ¼ 0.378, obtained in Ref. [11]. Because concurrence C is not applicable to more than two qubits, we use negativity N to detect entanglement in the eight-qubit system. For N > 2, N A;B is defined on a particular bipartition of the system into subsystems A and B. We define N to be the geometric mean of this quantity across all possible bipartitions. A nonzero N indicates the presence of global entanglement. Figures 4(c) and 4(d) show the negativity calculated with measured P 1 and P 2 (and with the measured Hamiltonian parameters Δ and EJ ij ) as a function of s for the two-and eight-qubit systems. The eight-qubit system has nonzero N for s < 0.315, thus indicating the presence of mixed-state global entanglement. Both concurrence C and negativity N decrease later in QA, where the first excited state approaches the ground state and becomes thermally occupied. The experimental values of these entanglement measures are in agreement with the theoretical predictions (solid lines in Fig. 4). The error bars in Figs. 4(c) and 4(d) represent uncertainties in the measurements of P 1 ðsÞ, P 2 ðsÞ, ΔðsÞ, and EðsÞ.
As stated above, the calculation of C and N relies on the assumption that the off-diagonal terms in the density matrix decay on time scales of several nanoseconds. We remove this assumption and demonstrate entanglement through the use of another witness W AB , defined on some bipartition A-B of the eight-qubit system. The witness, described in Appendix D, is designed in such a way that Tr½W AB σ ≥ 0 for all separable states σ. When Tr½W AB ρðsÞ < 0, the state ρðsÞ is entangled. Measurements of populations P 1 and P 2 provide a set of linear constraints on the density matrix of the system ρðsÞ. We then obtain an upper bound on Tr½W AB ρðsÞ by searching over all ρðsÞ that satisfy these linear constraints. If this upper bound is < 0, then we have shown entanglement for the bipartition A-B [37]. Figure 5 shows the upper limit of the witness Tr½W AB ρðsÞ for the eight-qubit system. We plot data for the bipartition that gives the median upper limit. The error bars are derived from a Monte Carlo analysis wherein we used the experimental uncertainties in Δ and J to estimate the uncertainty in Tr½W AB ρ. We also plot data for the two partitions that give the largest and smallest upper limits. For all values of the annealing parameter s, except for the last two points, upper limits from all possible bipartitions of the eight-qubit system are below zero. In this annealing range, the eightqubit system is globally entangled.

VII. CONCLUSIONS
In summary, we provide experimental evidence for the presence of quantum coherence and entanglement within subsets of qubits inside a quantum annealing processor during its operation. Our conclusion is based on four levels of evidence: (a) the observation of two-and eight-qubit avoided crossings with a multiqubit energy gap g ≫ k B T; (b) the witness Wχ, calculated with measured cross susceptibilities and coupling energies, which reports ground-state entanglement of the two-and eight-qubit systems [note that these two levels of evidence do not require explicit knowledge of Hamiltonian (1)]; (c) the measurements of energy eigenspectra and equibrium occupation fractions during QA, which allow us to use Hamiltonian (1) to reconstruct the density matrix, with some weak assumptions, and calculate concurrence and negativity (these standard measures of entanglement report nonclassical correlations in the two-and eight-qubit systems); (d) the entanglement witness W AB , which is calculated with the measured Hamiltonian and with constraints provided by the measured populations of the ground and the first excited states (this witness reports global entanglement of the eight-qubit system midway through the QA algorithm).
The observed entanglement is persistent at thermal equilibrium, an encouraging result as any practical hardware designed to run a quantum algorithm will be inevitably coupled to a thermal environment. The experimental techniques we discuss provide measurements of energy levels, and their populations, for arbitrary configurations of Hamiltonian parameters Δ, h i , J ij during the QA algorithm. The main limitation of the technique is the spectral width of the probe device. Improved designs of this device will allow much larger systems to be studied. Our measurements represent an effective approach for exploring the role of quantum mechanics in QA processors and ultimately to understanding the fundamental power and capability of quantum annealing.

ACKNOWLEDGMENTS
We thank C. Williams, P. Love, and J. Whittaker for useful discussions. We acknowledge F. Cioata and P. Spear for the design and maintenance of electronics control systems, J. Yao for fabrication support, and D. Bruce, P.

Chip description
The experiments discussed in herein were performed on a sample fabricated with a process consisting of a standard Nb/AlOx/Nb trilayer, a TiPt resistor layer, planarized SiO 2 dielectric layers, and six Nb wiring layers. The circuit design rules include a minimum linewidth of 0.25 μm and 0.6-μm-diameter Josephson junctions. The processor chip is a network of densely connected eight-qubit unit cells that are more sparsely connected to each other (see Fig. 1 for photographs of the processor). We report measurements made on qubits from one of these unit cells. The chip is mounted on the mixing chamber of a dilution refrigerator inside an Al superconducting shield and temperature controlled at 12.5 mK.

Qubit parameters
The processor facilitates quantum annealing of compound-compound Josephson junction rf SQUID flux qubits [38]. The qubits are controlled via the external flux biases Φ x qi and Φ x ccjj , which allow us to treat them as effective spins (see Fig. 1). Pairs of qubits interact through tunable inductive couplings [25]. The system can be described with the time-dependent QA Hamiltonian, where σ x;z i are Pauli matrices for the ith qubit, i ¼ 1; …; N. The energy scales Δ and E are the transverse and longitudinal energies of the spins, respectively, and the unitless biases h i and couplings J ij encode a particular optimization problem. We defineh i ≡ Eh i andJ ij ≡ EJ ij . We map the annealing parameter s for this particular chip to a range of Φ x ccjj with the relation where t f is the total anneal time. We implement QA for this processor by ramping the external control Φ The energy scale E ≡ M eff jI p q ðsÞj 2 is set by the s-dependent persistent current of the qubit jI p q ðsÞj and the maximum mutual inductance between qubits M eff ¼ 1.37 pH [8]. The transverse term in Hamiltonian (A1), ΔðsÞ, is the energy gap between the ground and first excited state of an isolated rf SQUID at zero bias. Δ also changes with annealing parameter s. Φ x qi ðtÞ is provided by a global external magnetic flux bias along with local in situ tunable DAC that tunes the coupling strength of this global bias into individual qubits and thus allows us to specify individual biases h i . The coupling energy between the ith and jth qubit is set with a local in situ tunable DAC that controls Φ x co;ij . The main quantities associated with a flux qubit, Δ and jI p q j, primarily depend on macroscopic rf SQUID parameters: junction critical current I c , qubit inductance L q , and qubit capacitance C q . We calibrate all of these parameters on this chip as described in Refs. [6,8]. We calibrate all interqubit coupling elements across their available tuning range from 1.37 pH to −3.7 pH, as described in Ref. [25]. We correct for variations in qubit parameters with on-chip control as described in Ref. [8]. This allows us to match jI p q j and Δ across all qubits throughout the annealing trajectory. Table I shows the median qubit parameters for the devices studied here. Figure 6 shows measurements of Δ and jI p q j versus s for all eight qubits. Δ is measured with single-qubit Landau-Zener measurements from s ¼ 0.515 to s ¼ 0.658 [39] and with qubit tunneling spectroscopy from s ¼ 0.121 to s ¼ 0.4 [33]. The resolution limit of qubit tunneling spectroscopy and the bandwidth of our external control lines during the Landau-Zener measurements prevent us from characterizing Δ between s ¼ 0.4 and s ¼ 0.5, respectively. jI p q j is   6. (a) ΔðsÞ versus s. We show measurements for all eight qubits studied in this work. We use a single-qubit Landau-Zener experiment to measure Δ=h < 100 MHz [39]. We use QTS to measure Δ=h > 1 GHz [33]. The red line shows the theoretical prediction for a rf SQUID model employing the median qubit parameters of the eight devices. The vertical black line separates coherent (left) and incoherent (right) evolution as estimated by analysis of single-qubit spectral line shapes. (b) jI p q jðsÞ versus s. We show measurements for all eight qubits studied in this work. We use a two-qubit coupled flux measurement with the interqubit coupling element set to 1.37 pH [8]. The red line shows the theoretical prediction for a rf SQUID model employing the median qubit parameters of the eight devices. measured by coupling a second probe qubit to the qubit q i with a coupling of M eff ¼ 1.37 pH and measuring the flux M eff jI p qi ðsÞj as a function of s. jI p q j is matched between qubits to within 3% and ΔðsÞ is matched between qubits to within 8% across the annealing region explored in this study.

APPENDIX B: QUBIT TUNNELING SPECTROSCOPY
QTS allows one to measure the eigenspectrum of an Nqubit system governed by Hamiltonian H S . Details on the measurement technique are presented elsewhere [33]. For convenience in comparing with this reference, we define a qubit energy bias ϵ i ≡ 2h i . Measurements are performed by coupling an additional probe qubit q P , with qubit tunneling amplitude Δ P ≪ Δ, jJj, to one of the N qubits of the system under study, for example, q 1 . When we use a coupling strengthJ P between q P and q 1 and apply a compensating bias ϵ 1 ¼ 2J P to q 1 , the resulting system plus probe Hamiltonian becomes For one of the localized states of the probe qubit, j↑i P , for which an eigenvalue of σ z P is equal to þ1 (i.e., the probe qubit in the right well), the contribution of the probe qubit is exactly canceled, leading to H SþP ¼ H S , with composite eigenstates jn; ↑i ¼ jni ⊗ j↑i P and eigenvalues E R n ¼ E n , which are identical to those of the original system without the presence of the probe qubit. Here, jni is an eigenstate of the Hamiltonian H S (n ¼ 1; 2; …; 2 N ).
For the other localized state of the probe qubit, j↓i P , when this qubit is in the left well, the ground state of H SþP is jψ L 0 ; ↓i ¼ jψ L 0 i ⊗ j↓i P , with eigenvalueẼ L 0 ¼ E L 0 þ ϵ P , where jψ L 0 i is the ground state of H S − 2J P σ z 1 and E L 0 is its eigenvalue. We choose jJ P j ≫ k B T such that the state jψ L 0 ; ↓i is well separated from the next excited state for ferromagnetically coupled systems, and, thus, system plus probe can be initialized in this state to high fidelity.
Introducing a small transverse term, − 1 2 Δ P σ x P , to Hamiltonian (B1) results in incoherent tunneling from the initial state jψ L 0 ; ↓i to any of the available jn; ↑i states [40]. A bias on the probe qubit ϵ P changes the energy difference between the probe j↓i P and j↑i P manifolds. We can thus bring jψ L 0 ; ↓i into resonance with any of jn; ↑i states (whenẼ L 0 ¼ E R n ), allowing resonant tunneling between the two states. The rate of tunneling out of the initially prepared state jψ L 0 ; ↓i is thus peaked at the locations of jn; ↑i.
The measurement of the eigenspectrum of an N-qubit system thus proceeds as follows. We couple an additional probe qubit to one of the N qubits (say, to q 1 ) with coupling constantJ P . We prepare the (N þ 1)-qubit system in the state jψ L 0 ; ↓i by annealing from s ¼ 0 to s ¼ 1 in the presence of large bias ϵ pol < 0 on all the system and probe qubits. We then adjust s for the N-qubit system to an intermediate point s Ã ∈ ½0; 1, such that Δ ≫ k B T=h, and adjust s for the probe qubit to s P ¼ 0.612, such that Δ P =h ∼ 1 MHz (here, h is the Planck constant). We assert a compensating bias ϵ 1 ¼ 2J P to this qubit. We dwell at this point for a time τ, complete the anneal s → 1 for the system plus probe, and then readout the state of the probe qubit. Figure 7 summarizes these waveforms during a typical QTS measurement. We perform this measurement for a range of τ, which allows us to measure an initial rate of tunneling Γ from jψ L 0 ; ↓i to jψ; ↑i. We repeat this measurement of Γ for a range of the probe qubit bias ϵ P . Peaks in Γ correspond to resonances between the initially prepared state and the state jn; ↑i, thus allowing us to map the eigenspectrum of the N-qubit system.
For the plots in the main paper, measurements of Γ are normalized to [0, 1] by dividing Γ by the maximum value across a vertical slice to give a visually interpretable result. Figure 8(b) shows a typical raw result in units of μs −1 .
We pose ferromagnetically coupled instances of the form with J ij < 0 for two-and eight-qubit subsections of the QA processor. Figure 8(a) shows typical measurements of Γ for FIG. 7. Typical waveforms during QTS. We prepare the initial state by annealing probe and system qubits from s ¼ 0 to s ¼ 1 in the presence of a large polarization bias ϵ pol . We then bias the system qubit q 1 (to which the probe is attached) to a bias ϵ 1 and the probe qubit to a bias ϵ P . With these biases asserted, we then adjust the system qubits' annealing parameter to an intermediate point s Ã and the probe qubit to a point s P and dwell for a time τ. Finally, we complete the anneal s → 1 and readout the state of the qubits.
a two-qubit subsection at several biases h i and at s ¼ 0.339 (J P < 0). We assemble multiple measurements to produce the spectrum shown in Fig. 8(b).

APPENDIX C: EQUILIBRIUM DISTRIBUTION OF SYSTEM
In addition to the energy eigenspectrum, QTS also provides a means of measuring the equilibrium distribution of an N-qubit system with a probe qubit. Suppose we are in the limit jJ P j ≫ k B T, such that there is only one accessible state in the j↓i P manifold: jψ L 0 i ⊗ j↓i P . As described above, the other available states in the system are the composite eigenstates jni ⊗ j↑i P in the j↑i P manifold, where jni is an eigenstate of the N-qubit system without the probe qubit attached. Energy levels E R n of the j↑i P manifold coincide with the energy levels E n of the system, E R n ¼ E n , even in the presence of coupling between the probe qubit and the system. We make the assumption that the population of an eigenstate depends only on its energy. Degenerate states have the same population.
Let P L represent the probability of finding the probe plus system in the state jψ L 0 i ⊗ j↓i P and let P R n represent the probability of finding the probe plus system in the state jni ⊗ j↑i P . At any point in the probe plus system evolution, we expect As described in the previous section, we can alter the energy of jψ L 0 i ⊗ j↓i P with the probe bias ϵ P . Based on the spectroscopic measurements of the N-qubit eigenspectrum, we can choose an ϵ P such that jψ L 0 i ⊗ j↓i P and jni ⊗ j↑i P are degenerate. Since the occupation of the state depends on its energy, we expect that, after long evolution times, these two degenerate states are occupied with equal probability, P L ðϵ P ¼ E n Þ ¼ P R n . Aligning the state jψ L 0 i ⊗ j↓i P with all possible 2 N states jni ⊗ j↑i P , we obtain a set of relative probabilities P R n . These relative probabilities characterize the population distribution in the system since they are uniquely determined by the energy spectrum E n . However, For better interpretability, we have subtracted off a baseline energy with respect to (a) such that the ground and first excited levels are symmetric about zero. Note the avoided crossing at h i ¼ 0. The peak tunneling rate Γ ∼ jΔ P hψ L 0 jnij 2 [33]. The solid black and white curves plot the theoretical expectations for the energy eigenvalues using independent measurements shown in Fig. 6 and Hamiltonian (1). as follows from Eq. (C1), the set P R n is not properly normalized. The probability distribution of the system itself is given by where P 2 N n¼1 P n ¼ 1. At every eigenenergy, ϵ P ¼ E n , the denominator of Eq. (C2) can be found from Eq. (C1), so that the population distribution of the system P n has the form Thus, the probability P n to find the system of N qubits in the state with energy E n can be estimated by measuring P L at ϵ P ¼ E n and using Eq. (C3).
Measurements of P L proceed as they do for the spectroscopy measurements. The system plus probe is prepared in jψ L 0 ; ↓i. We then adjust ϵ P ¼ E n and an annealing parameter s for the N-qubit system to some intermediate point, and also s P ¼ 0.612 for the probe qubit, such that Δ P =h ∼ 1 MHz. We dwell at this point for a time τ ≫ 1=Γ, complete the anneal s → 1, and then readout the state of the probe qubit. We typically investigate a range of τ to ensure that we are in the long evolution time limit in which P L is independent of τ. We use P L measured with τ ¼ 7041 μs to estimate P 1 and P 2 . The Supplemental Material [34] contains typical data used for these estimates.

APPENDIX D: ENTANGLEMENT MEASURES
AND WITNESSES

Definition of entanglement
A pure state jΨi of a system S consisting of two parts A and B, S ¼ A∪B, is entangled [16] with relation to this bipartition if the state jΨi cannot be represented as a product of states jΨ A i and jΨ B i describing the subsystems A and B: jΨi ≠ jΨ A i ⊗ jΨ B i. An open quantum system is characterized by a density matrix ρ. An open system S is entangled [16] relative to the bipartition S ¼ A∪B if its density matrix ρ cannot be written as a convex sum of product states ρ k Here, fρ k A g and fρ k B g are sets of density matrices for the components A and B, respectively; w k ≥ 0, P k w k ¼ 1. If there is no bipartition for which the system is entangled, the state is completely separable or unentangled. If the system is entangled for all possible bipartitions, the state is globally entangled.
2. W χ : A susceptibility-based, ground-state entanglement witness For a bipartion of the system into two parts, A and B, we define a witness R AB as where χ ij is a cross susceptibility,J ij ¼ EJ ij , and N AB is a number of nonzero couplings, J ij ≠ 0, between qubits from the subset A and the subset B (see Ref. [35]). We note that at low temperature, T ¼ 12.5 mK, the measured susceptibility χ ij ðTÞ almost coincides with the ground-state susceptibility χ ij ðT ¼ 0Þ since contributions of excited states to χ ij ðTÞ are proportional to their populations, P n ≪ 1, for n > 1. We analyze a deviation of the measured susceptibility from its ground-state value in the Supplemental Material [34]. To characterize global entanglement in the system of N qubits, we introduce a witness W χ , which is given by a bounded geometrical mean of witnesses R AB calculated for all possible partitions of the whole system into two subsystems. Here, N p is a number of such bipartitions, in particular, N p ¼ 127 for the eightqubit ring.

Thermal density matrix ρ
Systems coupled to a thermal bath are in a mixed state when the energy gap between the ground and first excited states approaches k B T. Most mixed-state entanglement measures and witnesses require knowledge of the density matrix ρ. The density matrix can be measured using quantum state tomography. However, this approach is limited to a small number of qubits [41]. As an alternate approach, we consider an N-qubit system in a thermal (stationary) state [26]. The stationary system described by Hamiltonian (1) can be characterized by a density matrix that is diagonal in the energy basis, ρ E ¼ diagfP 1 ; P 2 ; …; P 2 N g. The off-diagonal elements of this matrix disappear on a very short decoherence time scale (within a few nanoseconds) in the quantum annealing processor analyzed in Ref. [6]. In thermal equilibrium, the occupation probability P μ ¼ e −E μ =k B T =Z is the Boltzmann distribution over the eigenstates jμi with energies E μ , such that Hjμi ¼ E μ jμi, T is temperature, and Z ¼ P μ e −E μ =k B T is the partition function. The density matrix ρ ¼ P ρMNjMihNj in the computation basis jMi, where ρMN ¼ P μ P μ hMjμihμjNi, is obtained from ρ E by a unitary transformation that depends on the parameters of the Hamiltonian H. The parameters of Hamiltonian (1), namely, the energy scale E, the qubit biases h i , tunneling amplitudes Δ, coupling constants J ij , as well as probabilities P μ , can be independently measured. This allows us to restore the density matrix ρ. We emphasize that the stationary-state entanglement, or thermal entanglement, is robust and does not decay with time [26].

Concurrence C
For an open system described by a thermal density matrix ρ, concurrence C is defined as [16,18] where λ 4 > λ 3 > λ 2 > λ 1 are the roots of the matrix R ¼ ρðσ y 1 ⊗ σ y 2 Þρ Ã ðσ y 1 ⊗ σ y 2 Þ; (D4) fσ x i ; σ y i ; σ z i g are the Pauli matrices for the ith qubit, and ρ is the density matrix of the system in the computation basis.
An entanglement of formation E f is another measure of two-qubit entanglement. This measure is a monotonic function of concurrence C [16,18], where F ðxÞ ¼ −xlog 2 ðxÞ − ð1 − xÞlog 2 ð1 − xÞ is the entropy function.

Negativity N
Negativity is a measure that provides a sufficient, but not necessary, condition for entanglement of an arbitrary number of qubits [36]. A nonzero value of negativity detects entanglement. To calculate negativity, we find all bipartitions of the system. For two qubits, there is only one bipartition. An eight-qubit system can be bipartitioned into 127 (¼ 8 þ 28 þ 56 þ 35) possible combinations of two subsystems, A and B. In the case when the state of the system is separable, its density matrix ρ should retain all properties of the true density matrix after the partial transposition of ρ with respect to the subsystem A or to the subsystem B [16,36]. In particular, the partially transposed density matrix ρ T A should not have negative eigenvalues. The negativity, N ¼ ð1=2Þ P i ðjλ i j − λ i Þ, is proportional to the sum of all negative eigenvalues λ i of the matrix ρ T A , thus quantifying a degree of entanglement of the subsystem A and the subsystem B. For the eight-qubit system, we analyze negativities for all 127 partitions, N 1=7 , N 2=6 , N 3=5 , N 4=4 , and calculate their geometrical average [42], or the global negativity, N ðρÞ ¼ ðN 1=7 N 2=6 N 3=5 N 4=4 Þ 1=127 . Here, the negativity N 1=7 is equal to the product of 8 negativities for all possible partitions of the eight-qubit system into subsystems of one and seven qubits, and so on. Nonzero values of the global negativity mean that all possible subsystems of the whole eight-qubit system are globally entangled. Note that the negativity of the maximally entangled GHZ state, N GHZ , is equal to 1=2.

Entanglement witness W AB
Consider Hamiltonian (1) with measured parameters. This Hamiltonian describes a transverse Ising model having N qubits. The ground state jψ 1 i of this model is entangled with respect to some bipartition A − B of the N-qubit system. We can form an operator jψ 1 ihψ 1 j T A , where T A is a partial transposition operator with respect to the A subsystem [16]. Let jϕi be the eigenstate of jψ 1 ihψ 1 j T A with the most negative eigenvalue. We can form a new operator W AB ¼ jϕihϕj T A . This operator can serve as an entanglement witness (it is trivially positive on all separable states).