{'fraction_non_alphanumeric': 0.029938247479090126, 'fraction_numerical': 0.014330753797649757, 'mean_word_length': 4.537440484778531, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Observation of topological phenomena in a programmable lattice of 1,800 qubits

The celebrated work of Berezinskii, Kosterlitz and Thouless in the 1970s revealed exotic phases of matter governed by topological properties of low-dimensional materials such as thin films of superfluids and superconductors. Key to this phenomenon is the appearance and interaction of vortices and antivortices in an angular degree of freedom---typified by the classical XY model---due to thermal fluctuations. In the 2D Ising model this angular degree of freedom is absent in the classical case, but with the addition of a transverse field it can emerge from the interplay between frustration and quantum fluctuations. Consequently a Kosterlitz-Thouless (KT) phase transition has been predicted in the quantum system by theory and simulation. Here we demonstrate a large-scale quantum simulation of this phenomenon in a network of 1,800 in situ programmable superconducting flux qubits arranged in a fully-frustrated square-octagonal lattice. Essential to the critical behavior, we observe the emergence of a complex order parameter with continuous rotational symmetry, and the onset of quasi-long-range order as the system approaches a critical temperature. We use a simple but previously undemonstrated approach to statistical estimation with an annealing-based quantum processor, performing Monte Carlo sampling in a chain of reverse quantum annealing protocols. Observations are consistent with classical simulations across a range of Hamiltonian parameters. We anticipate that our approach of using a quantum processor as a programmable magnetic lattice will find widespread use in the simulation and development of exotic materials.

The celebrated work of Berezinskii, Kosterlitz and Thouless in the 1970s 1,2 revealed exotic phases of matter governed by topological properties of low-dimensional materials such as thin films of superfluids and superconductors. Key to this phenomenon is the appearance and interaction of vortices and antivortices in an angular degree of freedom-typified by the classical XY model-due to thermal fluctuations. In the 2D Ising model this angular degree of freedom is absent in the classical case, but with the addition of a transverse field it can emerge from the interplay between frustration and quantum fluctuations. Consequently a Kosterlitz-Thouless (KT) phase transition has been predicted in the quantum system by theory and simulation [3][4][5] . Here we demonstrate a large-scale quantum simulation of this phenomenon in a network of 1,800 in situ programmable superconducting flux qubits arranged in a fully-frustrated square-octagonal lattice. Essential to the critical behavior, we observe the emergence of a complex order parameter with continuous rotational symmetry, and the onset of quasi-long-range order as the system approaches a critical temperature. We use a simple but previously undemonstrated approach to statistical estimation with an annealing-based quantum processor, performing Monte Carlo sampling in a chain of reverse quantum annealing protocols. Observations are consistent with classical simulations across a range of Hamiltonian parameters. We anticipate that our approach of using a quantum processor as a programmable magnetic lattice will find widespread use in the simulation and development of exotic materials.
Richard Feynman's vision of simulating quantum systems with a quantum computer 6,7 has motivated the field of quantum information since its inception. In the absence of large-scale programmable universal quantum computers, advances in quantum simulation are bound by available technology. Still, remarkable progress has been made using near-term approaches such as ultracold Rydberg atoms 8 , superconducting qubits 9 , trapped ions 10,11 , and quantum dots 12 , among others 13 . Due to the possibility of noise-tolerant applications, quantum simulation has been identified as a poten-tial area of commercial value for near-term quantum computing technologies 14,15 . Quantum annealing (QA) processors [16][17][18] can be used to simulate systems in the transverse field Ising model (TFIM) described by the Hamiltonian where h i are longitudinal fields, J ij are coupling terms, σ x i and σ z i are Pauli matrices acting on the ith spin, and Γ is the transverse field. Evolution of this Hamiltonian in a low-temperature environment allows sampling of lowenergy solutions, with applications to optimization and machine learning 14,19,20 . This in turn allows estimation of equilibrium statistics of various phases dictated by the parameters h i , J ij and Γ at a given temperature T .
Geometrically frustrated magnets 21 are systems described by spin Hamiltonians with competing terms that cannot be minimized simultaneously. These systems give rise to a broad spectrum of exotic phases of matter. The TFIM exhibits a particular type of topological phenomenon on certain frustrated lattices: a KT phase transition separates a disordered paramagnetic phase from a phase in which complementary topological defectsvortices and antivortices-form bound pairs, resulting in polynomial decay of correlations 1,2,4 . This phase transition has been predicted in theory and simulation, but to our knowledge never observed experimentally in the TFIM. Other KT phase transitions have, however, been observed in many systems including superconducting films, superfluid helium films, trapped atomic gases, and hybrid tin-graphene Josephson junction arrays 10,22,23 .
Here we exhibit this phenomenon in a large-scale programmable quantum simulation. We perform this simulation using a superconducting QA processor consisting of 2,048 rf-SQUID flux qubits fabricated as an integrated circuit (Fig. 1a), in which Hamiltonian terms h i and J ij are specified with programmable on-chip control circuitry [16][17][18] . The topology of available nonzero coupling terms J ij consists of a regular bipartite grid. On this programmable substrate we define a sequence of fully-frustrated square-octagonal lattices with cylindrical boundary condition (Fig. 1b), with the largest system using 1,800 qubits (Supplemental S-II B). This lattice features chains of four ferromagnetically (FM) coupled qubits. In the large FM coupling limit, the low-energy description of this lattice follows the same Landau-Ginzburg-Wilson (LGW) theory as the widely studied 3-5,24-26 triangular antiferromagnetic (AFM) lattice (Fig. 1c). We therefore use the same analytic machinery even well outside the large FM coupling, lowtemperature limit.
Central to this analysis is a mapping from each plaquette in the triangular lattice to a 2D pseudospin that can be represented by a complex number ψ j . This gives the local angular degree of freedom in which topological features, including vortices and antivortices, can emerge. To first order perturbation in the small-Γ limit, an AFM triangle has six degenerate ground states, each consisting of an up-spin (|↑ ), a down-spin (|↓ ), and a spin aligned with the transverse field (|→ = (|↑ + |↓ )/ √ 2). These six ground states map to complex six-state pseudospins via the mapping As with a single AFM triangle, the triangular lattice has, to first order perturbation, a six-fold degenerate clock ground state: One third of the spins align with the transverse field, one third take spin up and one third take spin down. This order by disorder 4 , in which the addition of quantum fluctuations brings order to a highly degenerate classical ground state manifold, naturally divides the spins into three sublattices. Their magnetizations m 1 , m 2 , m 3 , defined as the average of σ z i over spins in the sublattice, give a complex order parameter 26 In the absence of open boundary effects, ψ is simply the average of pseudospin ψ j over all plaquettes in the system. Further details are given in Supplemental S-I. In the square-octagonal lattice, we use the same order parameter with the three sublattices defined according to the natural mapping: two FM-coupled spins must be in the same sublattice, and two AFM-coupled spins must be in different sublattices. Fig. 2a shows ψ j evaluated on each plaquette for an output state from the QA processor. Fig. 2b shows the phase diagram of the squareoctagonal lattice in the Γ-T plane, which resembles that of the triangular lattice 5 (Supplemental S-IV). At high temperature the system is in a disordered paramagnetic phase. As T drops, for Γ < 1.76, the system enters a critical phase at a KT phase transition, in which vortices and antivortices form bound pairs. At even lower temperature the system transitions into a clock ordered phase, where ψ concentrates around the six low-temperature clock states.
Having described the lattice and its phase diagram, we now turn to the task of simulating the physics with the QA processor. The processor Hamiltonian is parameterized by a time-dependent annealing parameter s that varies between 0 and 1, and for this particular lattice can be written as where J ij is 1.0 for AFM couplers and −1.8 for FM couplers. Details of implementing this Hamiltonian are given in Supplemental S-II. As s increases, J(s) increases and Γ(s) decreases: J(0) ≈ 0, J(1) T , Γ(0) T , Γ(1) ≈ 0. Fig. 2b shows the set of points in the Γ-T plane realized by varying s at the operating temperature of 8.4±0.2 mK; points above or below this line can be reached by either rescaling J ij or changing the QA operating temperature.
To estimate the distribution of ψ we use the QA processor as a quantum evolution update operator in a Monte Carlo method, where each evolution follows a reverse annealing protocol (Fig. 3a, details in Supplemental S-II A). This protocol involves initializing the system in a classical state at s = 1, quickly reducing s to a target value between 0.20 and 0.30 over approximately 3.5 µs, and pausing for 65 ms before quenching back to s = 1 over 1 µs and reading the classical output. Repetition of this protocol resembles a Markov chain by which we generate Monte Carlo estimates of ψ. This is a departure from the traditional use case of QA, in which s is steadily increased from 0 to 1. Typical of Markov chain Monte Carlo methods, we discard the first half of samples in each experiment as burn-in to minimize correlation with the initial state; a further motivation is that we observe around 0.1 mK of cooling at low temperatures, as heat from the QA programming cycle dissipates. We estimate statistics using both an ordered initial state and a random intial state (m = |ψ| large and small, respectively). Deviations between the two estimators are small and included in the error bars in subsequent plots.
We experimentally probe values of s between 0.20 and 0.30 at physical temperature T between 8.4 mK and 21.4 mK. Fig. 3b shows the expectation m as a function of s for the 1,800-spin lattice (lattice width L = 15) at 8.4 mK, compared with estimates from continuoustime path-integral quantum Monte Carlo (QMC) simulation as detailed in Supplemental S-III. In agreement with QMC, QA shows a peak in m close to where the QA schedule cuts through the KT phase in the large system limit, with a small shift caused by the open boundary condition and finite size effects. For s ≥ 0.24, QA shows good agreement with QMC. For smaller s, dynamics of the system are fast and the system orders during the 1 µs quench, leading to an overestimate of m ; results are almost unchanged between s = 0 and s = 0.2.
The complex order parameter ψ is plotted in Fig. 3c for two points along the QA schedule. As s increases from 0.20 to 0.26, both QA and QMC show the emergence of order: ψ concentrates around a ring with the characteristic U (1) symmetry of the 2D XY model despite the discrete Z 2 symmetry of the Ising model 5,24,25,27 .
In the thermodynamic limit, the KT phase transition is marked by a change of correlation decay from exponential to power-law at the critical temperature. In the critical phase, correlations C ij decrease with distance x ij as C ij ∝ x −b ij with an exponent b bounded above by the universal 2D XY critical exponent η = 1/4 28 . At finite sizes, the transition is expected to be broadened, and the apparent transition point can be affected by boundary conditions. To probe the onset of critical behavior we study correlations in the complex field ψ j on the plaquettes of the largest lattice on 1,800 spins; correlation is measured as C ij = Re(ψ i ψ * j ) on plaquettes along the periodic dimension, halfway between the open boundaries. Fig. 4a shows C ij as a function of distance for two temperatures at s = 0.26, with good agreement between QMC and QA. The absolute shift between QMC and QA correlations is explained by evolution during the quench (Supplemental S-V). As temperature decreases, correlations increase, approaching power-law decay: quality of power-law fit on distances 1 to 5-where boundary effects are minimal-improves as temperature decreases. This is shown in Fig. 4b for s = 0.25 and 0.26, near the peak of m in s (Fig. 3b). The exponent b of this power-law regression decreases with temperature ( Fig. 4c), reaching a minimum b ≈ 0.35 for QA and b ≈ 0.32 for QMC.
In the vicinity of the critical region, the same onset of power-law behavior is expected for m as a function of size with a halved exponent; i.e., m ∝ L −b/2 . This is shown in Fig. 4d, where QA shows agreement with QMC over a range of system sizes. Scaling exponents of b ≈ 0.34 and b ≈ 0.29 for QA and QMC respectively are close to the values extracted from phase correlations in Fig. 4c.
We  Power-law correlation decay near the KT crossover. a, Correlations in the complex field ψj, calculated as Cij = Re(ψiψ * j ) and measured between plaquettes along the periodic dimension far from the open boundary. Data for QMC and QA are shown near the critical temperature (8.4 mK) and above the critical temperature (21.4 mK). In the thermodynamic limit these correlations should decay with a power-law exponent close to η = 1/4 near the critical temperature. We perform power-law regression fit on distances 1 to 5, where boundary effects are minimal. b, As temperature is reduced, root-meansquared regression error drops, indicating the onset of power-law decay. c, In both QA and QMC, the power-law exponent b moves towards the expected critical value of η = 1/4 as temperature decreases. d, Decay of m versus L is expected to go from exponential in the paramagnetic phase to power-law with exponent b/2 near the KT crossover.
der parameter with rotational symmetry that responds as expected to changes in Γ and T , and the onset of powerlaw scaling of correlations. Taken together, these constitute the experimental observation of topological order in a frustrated 2D transverse field Ising model, as theoret-ically predicted 3,4 and simulated in QMC 5 . Agreement with QMC over a range of system sizes and Hamiltonian parameters-independently extracted with no fitting parameters-validates the flux qubit implementation of the transverse field Ising model at large scales.
The reverse annealing technique used in this work promises to greatly expand the utility of quantum annealing processors. The ability for a QA processor to reverse anneal from an input state allowed us to validate convergence of our statistical estimators to a steady state, but its importance is far more general: quantum evolution of an input state is crucial to the implementation of hybrid quantum-classical algorithms 29 .
This programmable magnetic material has allowed us to study the many-body dynamics of matter that is time consuming to simulate classically and difficult to implement physically. The methods used here can be applied to explore other exotic phases of matter in a variety of lattices. Several future developments will improve this approach. First, greater qubit connectivity will allow more flexibility in lattice geometry and realization of fully periodic boundary conditions. Second, faster projective readout will enable more accurate sampling from systems with fast dynamics, and any advances in noise and control error will improve the accuracy of these simulations. And finally, with the addition of non-stoquastic couplings 30 , quantum annealing processors could foreseeably simulate systems of which classical simulations are intractable even at modest scale. a hybrid tingraphene Josephson junction array," Nature Physics 10, 380-386 (2014

DATA AVAILABILITY
The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Quantum processor and experimental methods
The quantum annealing processor used in this work is a D-Wave 2000Q system operated in a non-annealing mode. Operating temperatures ranging between 8.4±0.2 mK and 21.4±0.2 mK are measured with a cryostat thermometer and independently verified using single-qubit susceptibility measurements. In the reverse anealing protocol we initialize the system at s = 1 with a classical state loaded into the qubits. We then reduce s quickly, bringing the system to the quantum model we want to simulate. We let this model relax for a fixed evolution time of 65 ms before quickly quenching s to 1.
Each call to the processor consists of 50 reverse anneals. A single initial state (clock or random) is loaded for the first anneal, and each subsequent anneal is initialized with the final state of the previous anneal. This allows us to perform quantum evolution Monte Carlo to estimate equilibrium statistics of the target system. From each processor call, we discard the first 25 of 50 samples as described in the main text. Data shown in the main text represent 480 processor calls (with the exception of Fig. 3c, which represents 720 processor calls) divided equally between clock initial state and random intial state, each of which is further divided equally between two physical embeddings of the lattice (Supplemental S-II B).
We make two adjustments to the QA schedule to correct deviation between the rf-SQUID flux qubits and the two-level Ising spins 17 . First we shift the schedule according to the deviation between the spectra of 2-level and 8level models in a 4-SQUID FM-coupled chain; the result is a nonlinear shift in s of approximately 0.02. Second we adjust the programmed coupling strengths to compensate for the background susceptibility χ b 31 . This effect results in a next-nearest-neighbor coupling between any two qubits that are coupled to a third qubit. The magnitude of χ b is a function of s and is determined by independent measurements. We compensate for this effect by adjusting couplings to homogenize effective AFM couplings and effective FM couplngs (see Supplemental S-II).
We maintain adjustments to the general-purpose calibration to minimize the effect input-specific crosstalks. This involves maintaining fine flux-bias offsets and coupling energy adjustments that respect the spin-flip and rotational symmetries in the cylindrical lattice (see Supplemental S-II D). These symmetries are independent of temperature and transverse field, so these adjustments do not result in overfitting.

Classical Monte Carlo methods
For all classical QMC experiments in this work we use continuous-time path-integral Monte Carlo 32 with Swendsen-Wang 33 updates. For a given spin system, a model is parameterized by T and Γ, where J = 1.
Here we describe methods for simulations and results described in the main text; further experiments and methods are discussed in the Supplemental Material. We employ parallel tempering using a series of Hamiltonians that vary in T , with either Γ or Γ/T fixed. We also employ a specialized chain update, in which every second Monte Carlo sweep attempts flipping of all spins in a four-qubit FM-coupled chain, rather than one spin at a time. We find that this speeds up convergence by up to three orders of magnitude.
Convergence for a set of Hamiltonians is determined via standard error on m 2 and convergence of absolute sublattice magnetizations to at most 0.04 from a fullymagnetized classical clock initial state, which is then checked for consistency with the same experiment given a random initial state. Convergence of the Binder cumulant and moments of the order parameter are tested with respect to initialization in several different types of clas-sical ground state. Error is determined self-consistently between five independent runs. Pseudospin correlations are computed from 5120 classical states projected from QMC in ten independent runs. Experimental lengths are up to 2 22 Monte Carlo sweeps.

Statistical methods
All error bars show 95% bootstrap confidence intervals over 1000 bootstrap samples. For QA results in Fig. 3b and Fig. 4, independent confidence intervals are generated for the Monte Carlo estimators with ordered initial state and random initial state, and the union of these intervals is shown. For every data point analyzed in Fig. 3c we also show the equivalent data points under reversal of all spins and/or rotation (but not flipping) of the cylinder. Thus we show 27,000 effective samples for QA and 30,720 effective samples for QMC.

Phase diagram of square-octagonal lattice
We determined the phase diagram of the squareoctagonal lattice based on the statistics of L×L instances with toridal boundary conditions with L between 3 and 21. Following Isakov and Moessner 5 we use two methods to determine the KT critical temperatures: the critical exponent of m and universal scaling collapse. To investigate the quantum critical point at Γ c ≈ 1.76, we show a crossing and universal collapse for the Binder cumulant for L between 3 and 15 with inverse temperature β(L, Γ) = 3.5L/Γ. Our results are further discussed and presented in Supplemental S-IV. Supplemental Materials: Observation of topological phenomena in a programmable lattice of 1,800 qubits

S-I. PSEUDOSPIN PHASE AND TOPOLOGICAL FEATURES IN THE TFIM
In the TFIM, the triangular AFM lattice and fully-frustrated square-octagonal lattice admit a mapping to the 2D XY model that is essential to the topological phenomena we observe. Here we provide a description of this mapping, which has its origins in the study of the stacked magnet 26 and has since been described in the quantum case 3-5,24,25 . We begin with a description of the XY model and the upper KT phase transition intersected by the QA schedule.

A. Two-dimensional XY model and the vortex unbinding transition
The 2D XY model describes an interaction of classical spins. Each spin is a unit vector described by an angle θ i . The interaction between spins is given by the XY Hamiltonian where the sum is taken over all coupled pairs of spins. The ground state of this system is one in which all spins are aligned (see Fig. S1); note that the system described by (S1) has continuous O(2) or U (1) symmetry. In a ground state, the angle θ does not change as we traverse a closed path in the plane in a clockwise direction. Given a perturbation of a ground state, θ may fluctuate along this path but will not wind in a full rotation. A clockwise or anticlockwise winding is known as a vortex or an antivortex respectively; these topological defects imply an excitation that diverges with the radius R of the encircling path: where r is the lattice spacing. The entropy S vortex also grows proportionally to log R r since there are (R/r) 2 possible locations of the vortex in an area of R 2 . A competition of these two terms leads to a sign-variable expression for the free energy contribution of a vortex: The critical vortex unbinding temperature T 2 at which the Kosterlitz-Thouless phase transition occurs is given by the sign change of ∆F : Above this temperature, the free energy favors isolated vortices; below this temperature vortices and antivortices are attracted to one another and appear in bound pairs.

B. Pseudospins in geometrically frustrated quantum magnets
The mapping from the TFIM to the XY model arises from the interaction of frustration and quantum fluctuations. We first consider the effect of a perturbative transverse field on a single AFM triangle. In the classical case, the AFM triangle has six ground states. In each, one of two floppy spins can be flipped without changing the energy. Fig. S2 shows two such states differing by a flip of spin 1. The addition of a perturbative transverse field Γ lifts this degeneracy, allowing a superposition of these two states in which spin 1 aligns with the transverse field. We denote this superposition by |→ = (|↑ + |↓ ) √ 2 and note that this lowers the energy from −J to −J − Γ. To first order perturbation, this quantum ground state is six-fold degenerate, with the three spins taking values |↑ , |↓ , |→ in any permutation. We map these six states to the complex plaquette pseudospins ψ j as described by (2): ψ j = σ z 1 + σ z 2 e 2πi/3 + σ z 3 e 4πi/3 . This is shown in Fig. S3 for the state |→↑↓ . Pseudospins for six quantum and six classical states are shown in Fig. S3b. Also shown are the classical excited states |↑↑↑ and |↓↓↓ , which have ψ j = 0. As vanishing points of the pseudospin phase, these magnetized triangles correspond to vortices and antivortices in the lattice.
As shown in Fig. S4a, the trianglular AFM lattice can be tiled with a single clock state, resulting in a pseudospin that is constant across the lattice. This tiling naturally divides the sites of the lattice into three sublattices, which we use to determine the complex order paramater ψ as in (3). The introduction of a twist in the pseudospin phase, as shown in Fig. S4b, leads to a state where two of the triangles contain no spin aligned with the transverse field. This excitation can be described in terms of an effective XY model on the dual lattice, with effective XY coupling proportional to Γ. This is consistent with the linear scaling of T 2 in Γ seen in the small-Γ limit 34 . Fig. S5 shows a vortex in the triangular lattice, and in the corresponding state in the square-octagonal lattice. To compute the pseudospin of a plaquette in the square-octagonal phase we use all 12 spins in the three FM-coupled chains intersecting the plaquette, whether it is a square or an octagon.

S-II. QUANTUM PROCESSOR AND EXPERIMENTAL METHODS
The quantum annealing processor used in this work is a D-Wave 2000Q system operated in a non-annealing mode. Operating temperatures ranging between 8.4 ± 0.2 mK and 21.4 ± 0.2 mK are measured with a cryostat thermometer and independently verified using single-qubit susceptibility measurements.

A. Annealing schedule and reverse annealing protocol
We focus on the region of the annealing schedule between s = 0.20 and s = 0.30 (Fig. S6a). In quantum annealing, the annealing parameter s is typically increased smoothly from 0 to 1, bringing the system from a single-well superposition at s = 0 to a low-temperature classical system at s = 1 (Fig. S6b). In our experiments we take a different approach, using a reverse annealing protocol. In this protocol we initialize the system at s = 1 with a classical state loaded into the qubits. We then reduce s quickly, bringing the system to the quantum model we want to simulate. We let this model relax for a fixed evolution time of 65 ms before quickly quenching s to 1 (Fig. S6c).
Each call to the processor consists of 50 reverse anneals. A single initial state (clock or random) is loaded for the first anneal, and each subsequent anneal is initialized with the final state of the previous anneal (Fig. S6c). This allows us to perform quantum evolution Monte Carlo to estimate equilibrium statistics of the target system. From each processor call, we discard the first 25 of 50 samples as described in the main text.
When Γ is small, the rf-SQUID flux qubits faithfully model two-level systems. Since our experiments are performed in a region of the schedule where Γ is large, behavior deviates significantly from a two-level system but is well described by an eight-level model 17 . To account for this we make two adjustments. First we shift the schedule according to the deviation between the spectra of 2-level and 8-level models in a 4-SQUID FM-coupled chain. Second we adjust the programmed coupling strengths to compensate for the background susceptibility effect described in Supplemental S-II C.

B. Embedding the square-octagonal lattice
The cylindrical topology of the square-octagonal lattice is realized as two square sheets coupled together at the top and bottom ( Figure S8). In the processor used, the graph of available couplings is a 16 × 16 grid of Chimera unit cells 18 ; only a 7 × 7 block of cells is shown in Fig. S8b. Each FM-coupled chain of four qubits consists of two qubits in one unit cell, and one qubit in each of two other unit cells. If the size of the unit cell were doubled, it would be possible to embed a toroidal square-octagonal lattice with fully periodic boundary conditions. The halving of the antiferromagnetic couplers on the open boundaries of the cylinder ensures that the classical ground state space reflects that of the toroidal lattice.

C. Background susceptibility and compensation
The rf-SQUID flux qubits in the QA processor provide an imperfect implementation of two-level Ising spins. The primary deviation from ideal is that early in the anneal, every Ising spin mediates a coupling between its neighboring spins. This introduces next-nearest-neighbor couplings, meaning that for any spins i, j, k, there is an effective coupling term χ b J ij J jk added to J ik ; the strength is denoted χ b , representing normalized background susceptibility: χ b = M AFM χ q where M AFM is the maximum available AFM mutual inductance and χ q is the physical qubit susceptibility. In the QA experiments, we compensate for χ b based on the fact that FM couplers are seldom frustrated in the systems studied. We consider the application of next-nearest-neighbor terms as a function f applied to the coupling Hamiltonian J QA supplied to the processor, and seek a choice of J QA such that f (J QA ) approximates αJ for some constant α > 0. This gives us the following constraints on the classical Hamiltonians: 1. Any two chains have the same total coupling between them in f (J QA ) and αJ. In the absence of boundary conditions, this coupling is α.
2. If a chain has a single break (domain wall), the total strength of frustrated couplers is the same in f (J QA ) and αJ. In the absence of boundary conditions, this coupling is −1.8α.
We use an iterative method to find J QA given an average AFM coupling strength ρ; it follows that α is a function of ρ and χ b (in turn, χ b is a function of s). For ρ = 0.95 the value of α ranges monotonically from 1.32 at s = 0.25 to 1.13 at s = 0.4.

D. Calibration refinement
The QA processor was calibrated as a general-purpose quantum annealer, with the goal of good performance across a variety of inputs. Given our focus on one particular input type, we can improve performance by exploiting structure and symmetries in the input. After compensating for background susceptibility as described above, we exploit symmetries in the lattice.
The first symmetry is spin-flip invariance: since there is no longitudinal field h, the expected magnetization of each spin in the absence of systematic biases should be zero. To maintain this degeneracy we apply flux-offset biases to qubits in a gradient descent method, significantly improving the distribution of qubit magnetizations (Fig. S9a). The required offset is a function of annealing schedule, so we maintain an independent compensation for each schedule based on forward annealing results.
The second symmetry is rotational symmetry about the periodic dimension. For lattice width L (with 4·L·2L spins), each coupler is in an equivalence class of 2L rotationally equivalent couplers that should be frustrated equiprobably.
More generally, two couplers should have the same statistics if there is a graph automorphism over the lattice that maps one to the other. We make fine adjustments to individual coupling terms to tighten the frustration distributions of these equivalence sets (Fig. S9b).
There are no further trivial symmetries in the square-octagonal lattice; both the triangular lattice with semi-open boundary and the square-octagonal lattice with fully periodic boundary have richer automorphism groups that might be exploited in future research. The symmetries that we use for calibration refinement apply at every point in the phase diagram, so the calibration refinement used here does not result in overfitting.

E. Spin-bath polarization
The persistent current flowing in the qubit bodies during the QA protocol produces a magnetic field that can partially align or polarize the spin system. For long evolution times, the polarized environment can produce sampleto-sample correlations, biasing the QA towards previously achieved spin configurations. To reduce sample-to-sample correlations, we introduce a pause of 10 ms between anneals to give the spin-bath time to depolarize. However in the reverse annealing protocol we start in a polarized state, and the QA system will be biased towards the initial state by polarization in the spin-bath. This leads to the small systematic difference in estimators of m using a clock initial state versus a random initial state.

S-III. CLASSICAL MONTE CARLO METHODS
For all classical QMC experiments in this work we use continuous-time path-integral Monte Carlo 32 with Swendsen-Wang 33 updates. For a given spin system, a model is parameterized by T and Γ, where J = 1.
When seeking equilibrium estimates of a model, we employ parallel tempering using a series of models that vary in T , with either Γ or Γ/T fixed. We also employ a specialized chain update, in which every other Monte Carlo sweep attempts flipping of all spins in a four-qubit FM-coupled chain, rather than one spin at a time. We find that this speeds up convergence by up to three orders of magnitude.
Convergence for a set of models is determined via standard error on m 2 and convergence of absolute sublattice magnetizations to at most 0.04 from a fully-magnetized clock initial state, which is then checked for consistency with the same experiment given a random initial state. Convergence of the Binder cumulant and moments of the order parameter are tested with respect to initialization in several different types of classical ground state. Error is determined self-consistently between independent runs; error bars on QMC results are 95% bootstrap confidence intervals.

S-IV. PHASE DIAGRAM OF SQUARE-OCTAGONAL LATTICE
Here we provide evidence for the claimed phase diagram of the square-octagonal lattice. We show critical behavior and give estimates of the upper and lower critical temperatures bounding the critical KT region in the phase diagram of the square-octagonal lattice, and give an estimate of the quantum critical point at T = 0. We follow the methods used by Isakov and Moessner 5 on the triangular lattice. We study L × L lattices (on 4L 2 spins) with fully periodic (toroidal) boundary conditions with L ranging from 3 to 21. Our results are presented in Figure S10.

A. KT critical temperatures
For a range of transverse fields we use two methods to determine the lower (T 1 ) and upper (T 2 ) critical temperatures using properties of the 2D XY universality class.
First, we use the expected power-law decay of correlations-and consequently m-as system size increases. Within the critical region and for fixed Γ, m L is expected to scale as L −η/2 , with critical exponents η 1 = 1/9 and η 2 = 1/4. From this we determine T 1 and T 2 using a power-law fit on values of L between 6 and 21 (Fig. S10).
Second, we additionally determine T 1 and T 2 by fitting our data to universal scaling curves. For the lower transition, we have where m 0 is a universal function, a is a nonuniversal constant, and t = (T 1 − T )/T 1 is the residual temperature approaching T 1 from below. For the lower transition, we have where χ 0 is a universal function, a is a nonuniversal constant, and t = (T − T 2 )/T 2 is the residual temperature approaching T 2 from above. We expect from universality that in these fits, b ≈ (1/9)/2 and c ≈ 7/4. For the upper transition, we obtain good collapses for Γ between 0.7 and 1. 3. These values are in fairly good agreement with values of T 1 determined from η, as in the triangular lattice. 5 For the lower transition, fitting is a generally poor and we only obtain a reasonably convincing collapse for Γ ∈ [1.1, 1.2] by discarding data for L < 15. Values of T 1 deviate significantly from those given by η. A study of larger instances would clarify the picture.

B. Quantum critical point
Again extending from the triangular lattice, we anticipate a quantum phase transition consistent with the 3D XY universality class 5,25,26 with dynamical exponent z = 1 and universal exponent ν = 2/3. To estimate the quantum critical point, we measure the normalized Binder cumulant We do so for system sizes up to L = 15 using inverse temperature β(L, Γ) = 3.5L z /Γ. The Binder cumulant U crosses near a critical point Γ c ≈ 1.76. Scaling of U in the vicinity of the quantum critical point collapses as U L = L 1/ν (Γ − Γ c )/Γ c using ν = 2/3.

S-V. EFFECT OF QUENCH
The QA simulation generally provides shows good agreement with QMC in estimates of m , consistent with the formation and annihilation of vortex-antivortex pairs. However, these excitations appear far less often in QA output states than in projected QMC states. Here we show that this can be explained by evolution during the 1 µs QA quench. As quantum and thermal fluctuations are reduced, tightly-bound vortex-antivortex pairs are annihilated. Figure S11 compares the mean residual classical energy per spin of QA and QMC. QMC projected states have many more excitations than QA states. We model a local classical quench in QMC output by repairing classical excitations at the single-qubit level and four-qubit chain level. After this classical quench is applied to QMC output, residual energies of QMC states resemble those of QA states. Effect on ψ is minimal. This is consistent with the hypothesis that similar annihilation of defects occurs during the QA quench. We therefore expect that faster QA quench should lead to greater population of classical excitations-vortices and antivortices-in QA output.