Simulating the Shastry-Sutherland Ising Model using Quantum Annealing

Frustration represents an essential feature in the behavior of magnetic materials when constraints on the microscopic Hamiltonian cannot be satisfied simultaneously. This gives rise to exotic phases of matter including spin liquids, spin ices, and stripe phases. Here we demonstrate an approach to understanding the microscopic effects of frustration by computing the phases of a 468-spin Shastry-Sutherland Ising Hamiltonian using a quantum annealer. Our approach uses mean-field boundary conditions to mitigate effects of finite size and defects alongside an iterative quantum annealing protocol to simulate statistical physics. We recover all phases of the Shastry-Sutherland Ising model -- including the well-known fractional magnetization plateau -- and the static structure factor characterizing the critical behavior at these transitions. These results establish quantum annealing as an emerging method in understanding the effects of frustration on the emergence of novel phases of matter and pave the way for future comparisons with real experiments.

Frustration represents an essential feature in the behavior of magnetic materials when constraints on the microscopic Hamiltonian cannot be satisfied simultaneously. This gives rise to exotic phases of matter including spin liquids [1], spin ices [2], and stripe phases [3]. Here we demonstrate an approach to understanding the microscopic effects of frustration by computing the phases of a 468-spin Shastry-Sutherland Ising Hamiltonian using a quantum annealer.
Our approach uses mean-field boundary conditions to mitigate effects of finite size and defects alongside an iterative quantum annealing protocol to simulate statistical physics. We recover all phases of the Shastry-Sutherland Ising model -including the well-known fractional magnetization plateau -and the static structure factor characterizing the critical behavior at these transitions. These results establish quantum annealing as an emerging method in understanding the effects of frustration on the emergence of novel phases of matter and pave the way for future comparisons with real experiments.
Materials simulation using quantum annealing (QA) relies on fundamentally different principles as compared to conventional computing methods. Whereas traditional simulated annealing uses thermal excitations and the ergodic principle to find the energetic ground state [4], quantum annealing uses an additional transverse magnetic field as a quantum tuning parameter that drives transitions between quantum states [5,6]. The dynamics induced by the transverse field allows QA to explore the energy landscape, potentially faster than classical approaches [7].
We use a quantum annealer realized by a two-dimensional lattice of superconducting qubits that implements a transverse-field Ising model (TFIM) defined by [6] where the amplitudes A(s) and B(s) determine the relative strength of the transverse field and Ising terms, and are controlled through the dimensionless parameter s ∈ [0, 1]. The amplitudes A(s) and B(s) are smooth functions of s with A(0) >> B(0) and A(1) << B (1). The parameters J ij and h i define tunable Ising interactions and a per-qubit tunable longitudinal magnetic field, respectively, which enables the realization of various magnetic systems and lattice geometries within the processor. The available couplers are arranged in a "Chimera" graph [8], used by several recent examples to validate materials simulations 2 using QA. This includes the simulation of the three-dimensional (3D) spin-glass transition [9] and the demonstration of the Berezinskii-Kosterlitz-Thouless (BKT) transition in a twodimensional TFIM [10].
We investigate geometric frustration within an Ising model represented by the Hamilto- defined over the Shastry-Sutherland lattice [11], as shown in Fig. 1(a). In the Shastry-Sutherland Ising model, geometric frustration originates from the nearest-neighbor J 1 and the next-nearest neighbor J 2 interactions acting on the square frustrated grid. In particular, this model exhibits multiple phases including a trivial ferromagnetic (FM) phase when h z J 1 , J 2 , a Néel anti-ferromagnetic (AFM) phase when J 1 h z , J 2 , and a more interesting AFM spin-dimer phase when J 2 h z , J 1 [12]. However, for J 1 ≈ J 2 > 0, a highly non-trivial phase arises in which spins order to form a 1/3 magnetization plateau. This plateau arises from a 6-fold degenerate solution in the model above a critical longitudinal field strength [12,13].
Magnetization plateaus in Shastry-Sutherland models have generated significant interest historically. Notably, the Shastry-Sutherland Ising model and its variations have been proposed to account for the fractional magnetic phases observed in frustrated magnets such as the rare-earth tetraborides (RB 4 ) [13][14][15], although there is some debate on the role of additional terms beyond equation (2). The exact solubility of the Ising model and its simple square topology makes a compelling testbed to evaluate the efficacy of QA for addressing fundamental materials science and exploring a new pathway for the cross-examination of theory and experimental data.
We use QA to simulate the macroscopic behaviors of the Shastry-Sutherland Ising model by sampling the low-energy states across all phases of the model. Our approach goes beyond conventional QA approaches by using an iterative simulation protocol known as Quantum Evolution Monte Carlo (QEMC) chaining [7,10] to sample the ground-state manifold of the frustrated Ising Hamiltonian. We demonstrate this method across the full range of Hamiltonian parameters to recover the complete phase diagram. In order to simulate accurately the bulk behavior, we also apply a technique for mean-field boundary conditions [16,17] to mitigate finite-size effects as well as the presence of defects in the quantum annealer. We then recover both the phase diagram and the static structure factor for the model Hamiltonian from the simulated ensemble of spin states. These results validate the use of programmable quantum annealing as an avenue to understanding frustration in a magnetic Hamiltonian and the potential to simulate the macroscopic behaviors of bulk materials suitable for future experimental comparison.
The geometry of the Shastry Sutherland Ising model was embedded into 1,872 superconducting qubits in a Chimera graph on a D-Wave 2000Q quantum annealer [18]. The embedding of the Hamiltonian into the processor uses strongly-coupled ferromagnetic chains to encode each logical spin site [19]. The embedding is described in Fig. 1b, where each logical spin maps to a cyclic chain of 4 physical spins, each logical dimer bond to 8 physical internal couplers, and each logical square bond to one physical external coupler. This so-called "halfcell" embedding exploits the symmetry of the Shastry-Sutherland lattice to ensure dimer couplers map to internal couplers and square couplers map to external couplers. The largest model obtained with the half-cell embedding has 468 logical spins as shown in Fig. 1c.
The simulation of the embedded Hamiltonian was susceptible to the defects and finite-size effects shown in Fig. 1c, which can undermine estimates of bulk properties. We introduced mean-field boundary conditions in the embedded Hamiltonian in order to better approximate the equilibrium statistics of the thermodynamic limit. For an ideal, translationally invariant model, the magnetization at the edge of the embedded lattice reflects the bulk system far from any other interface or boundary. We leveraged independent per-qubit control over the longitudinal magnetic fields on each spin to mimic this expected behavior. The longitudinal field applied to each spin site on the boundary was iteratively adjusted to ensure homogeneous magnetization between interior and boundary qubits (see Methods).
Convergence of these boundary fields was obtained using gradient descent methods with a logical constraint that the boundary magnetic fields must match the sign of the uniform bulk magnetic field. ferromagnetic phase to the plateau phase was found at J 1 = 1, J 2 = 1, h z ≈ 1.8, and the transition from the plateau phase to the ferromagnetic phase was found at J 1 = 1, J 2 = 1, h z ≈ 5.
In addition to the correct order parameter, we observe the correct microscopic ordering expected in the different phases. In Fig. 3c we obtain the typical Néel AFM ordering with zero average magnetization. In Fig. 3e,f we display two of the six observed degenerate spin structures found within the plateau phase, each with m = 1/3. Shown in Fig. 3d  analysis of diffuse neutron scattering spectra can be utilized here to characterize the spatial correlations in the ensemble of magnetic orders.
In order to quantify the spatial correlations, we calculated the static structure factor S( q) defined in terms of the two-point correlation function as where R ij is the relative position of two spins and q is the wave vector in reciprocal space. The static structure factor provides a Fourier decomposition of the spatial correlations encoded by the spin system and quantifies the ordering of different magnetic phases. The calculated correlations for the Néel AFM, dimer AFM, and 1/3 plateau phase are shown in Fig. 4a-(c) and agree with theoretical expectations [13]. Because the static structure factor is directly measured by the neutron diffuse scattering spectrum, we anticipate similar QA results for appropriately modified Hamiltonians will enable comparisons with real Shastry-Sutherland magnets in the future.
Unlike the standard forward-annealing use of QA, QEMC chains allow us to probe ordering and statistical convergence near critical regimes [7] including, for example, the first-order transition from the Néel AFM to the 1/3 plateau. However, while theory predicts a sharp discontinuity in the magnetization, our annealing simulations produce a more gradual change as shown in Fig. 4e. This smoothing has many origins, including disorder in J 1 , J 2 , defects and finite size effects, as well as non-negligible transverse-field effects persisting near the 10 boundary conditions, we simulated the low-energy manifold of the Shastry-Sutherland Ising Hamiltonian over a range of model parameters. By accurately recovering the complex phase diagram for this model system as well as the static structure factors for each phase, our results indicate that QA can provide a powerful utility for understanding the emergence of macroscopic behavior from statistical physics. This work indicates that simulations using a quantum annealer of more complex Hamiltonians and those of real materials is no longer a matter of a conceptual, but rather a technical, challenge. Some of these challenges may be addressed by improved quantum annealing hardware as well as more sophisticated methods of materials simulation [21].

11
FIG. S1. Half-cell embedding into the D-Wave quantum annealer with unit cell as depicted in Fig. 1 (b)

METHODS
The D-Wave 2000Q quantum annealer is constructed of superconducting flux-qubits coupled together in the Chimera topology [18]. The quantum annealer has two types of couplers: internal unit cell couplers, which control interactions within a unit cell, and external couplers, which mediate interactions between unit cells. Here we present an embedding, shown in Fig. 1 and Fig. S1, which places one logical dimer onto each Chimera unit cell and realizes square bonds between dimers through inter-cell Chimera couplers. This embedding leverages the natural structure of the quantum annealer architecture to reinforce the symmetries of the logical model.

Chi-compensation
The radio-frequency SQUID flux qubits in the D-Wave 2000Q system do not perfectly implement the TFIM Hamiltonian. Specifically, each qubit mediates an effective coupling between neighboring qubits regardless of the existence of a physical coupler. Additionally there is leakage of an applied h bias from a qubit to its neighbors. The strength of these additional effects are dependent on the normalized background susceptibility where M AF M is the maximum available AFM mutual inductance and χ q is the physical qubit susceptibility. These effects lead to a first-order modification to the embedded hs and Js as shown below.
For the low-noise D- here it was sufficient to derive expressions of the impact of these interactions on the logical problem, which are given below. We used these expressions to determine the correct input such that the embedded system corresponded to the correct logical problem instance.

Flux-bias offset calibration
The D-Wave 2000Q system is calibrated as a quantum annealer to perform well under a variety of input types, and the unique problem here permits the re-calibration of the device for improved performance. This calibration should be done to reinforce the expected symmetries in the system, and in the absence of an applied longitudinal field the system is invariant under a spin-flip operation. Using this symmetry we determine flux-bias offsets for each qubit with a gradient descent method such that the average qubit magnetization is 13 zero. The measurement of this magnetization is performed via forward anneals in this work but is generally a function of anneal schedule.

Mean-field boundary conditions
Due to the nature of the half-cell embedding, any missing qubit or coupler on the D-Wave 2000Q chip leads to a missing qubit or coupler in the logical SS graph, respectively. Due to the finite size of the chip and the presence of these defects appropriate boundary conditions much be chosen in order to recover the expected solutions in the thermodynamic limit.
Here we implement a form of mean-field boundary conditions [16,17] which are found via an optimization over the longitudinal magnetic fields of qubits on external and internal boundaries. Specifically, we seek to determine the appropriate magnetic fields such that the magnetization of boundary qubits is equivalent to the magnetization of bulk qubits and is defined as where h bound. is the vector of longitudinal magnetic fields on the boundary and m bulk , m bound. are the magnetization of the bulk and boundary qubits, respectively. The sign constraint is used to restrict the boundary fields to physically realistic local environments.
The optimization is done via a simple gradient descent method that updates each boundary longitudinal field according to where σ z (i) is the average magnetization of qubit i over all samples and 0 ≤ δh ≤ 0.05 is the step size. We typically see convergence of this parameterization up to statistical noise in at most 500 iterations.

Quantum Evolution Monte Carlo Chains
In quantum annealing the goal is to evolve a system adiabatically from an eigenstate of the transverse field at s = 0 to an eigenstate of the Ising model at s = 1. However for systems composed of chains of ferromagnetically coupled qubits there may be a value of s at which the spins of the chain fail to evolve meaningfully and freeze [22]. This freeze out point

15
Here we utilize the method of King et al. and chain together a sequence of QEMC steps [10]. In a reverse anneal the system is initiated in a classical state at s = 1, follows an anneal schedule, and then is read out to obtain a new classical state. The anneal schedule used here is shown in Fig. S2(b) and is quite similar to that used by King et al. [10]: rapidly reverse-anneal to a point in the anneal schedule with non-negligible transverse field, pause to allow the system to thermalize and populate low energy intermediate states, and forward anneal to finish the protocol. When repeated, this protocol resembles a Markov chain Monte Carlo and enables the determination of ground states through an iterative refinement of solutions shown in Fig. S2(c). Just as in a conventional Markov chain simulation, we discard the first half of the iterations as burn-in steps while the system equilibrates and then sample from the last half of the chain. The size of chain length depends greatly on the pause point s p as can be seen in Fig. 2(b).
All experiments performed in this work used the same reverse anneal time, t r = 1 µs and used a chain length of 100, so that 50 steps were used for chain equilibration and 50 states are used for statistical sampling. For the system under study the time scale of dynamics is assumed to be much larger than the time scale of the anneal and readout, therefore a rapid forward anneal permits more time to pause and thermalize while not populating high energy states during the forward anneal. Finally, the initial state to begin each QEMC chain in this work was chosen to be the, FM, all spin-up state. We found little dependence of the final state of the chain on the initial state if s p was small enough such that the transverse field was significantly strong (s p ≤ 0.6), as can be seen in Fig.2(b). 16