Quantum annealing simulation of out-of-equilibrium magnetization in a spin-chain compound

Geometrically frustrated spin-chain compounds such as Ca3Co2O6 exhibit extremely slow relaxation under a changing magnetic field. Consequently, both low-temperature laboratory experiments and Monte Carlo simulations have shown peculiar out-of-equilibrium magnetization curves, which arise from trapping in metastable configurations. In this work we simulate this phenomenon in a superconducting quantum annealing processor, allowing us to probe the impact of quantum fluctuations on both equilibrium and dynamics of the system. Increasing the quantum fluctuations with a transverse field reduces the impact of metastable traps in out-of-equilibrium samples, and aids the development of three-sublattice ferrimagnetic (up-up-down) long-range order. At equilibrium we identify a finite-temperature shoulder in the 1/3-to-saturated phase transition, promoted by quantum fluctuations but with entropic origin. This work demonstrates the viability of dynamical as well as equilibrium studies of frustrated magnetism using large-scale programmable quantum systems, and is therefore an important step toward programmable simulation of dynamics in materials using quantum hardware.

Geometrically frustrated spin-chain compounds such as Ca3Co2O6 exhibit extremely slow relaxation under a changing magnetic field. Consequently, both low-temperature laboratory experiments and Monte Carlo simulations have shown peculiar out-of-equilibrium magnetization curves, which arise from trapping in metastable configurations. In this work we simulate this phenomenon in a superconducting quantum annealing processor, allowing us to probe the impact of quantum fluctuations on both equilibrium and dynamics of the system. Increasing the quantum fluctuations with a transverse field reduces the impact of metastable traps in out-of-equilibrium samples, and aids the development of three-sublattice ferrimagnetic (up-up-down) long-range order. At equilibrium we identify a finite-temperature shoulder in the 1/3-to-saturated phase transition, promoted by quantum fluctuations but with entropic origin. This work demonstrates the viability of dynamical as well as equilibrium studies of frustrated magnetism using large-scale programmable quantum systems, and is therefore an important step toward programmable simulation of dynamics in materials using quantum hardware.

I. INTRODUCTION
Geometrically frustrated magnetic systems exhibit a variety of interesting dynamical and equilibrium properties. The calcium oxalate Ca 3 Co 2 O 6 , a canonical example of a geometrically frustrated spin-chain compound, can be extremely slow to relax from metastable configurations when subjected to a changing magnetic field. The magnetic Co 3+ ions align in ferromagnetic (FM) spin chains along the c axis that form a triangular antiferromagnetic lattice in the a-b plane. This compound shows unusual out-of-equilibrium magnetization curves, a phenomenon which has been studied both in situ [1][2][3][4][5] and in simulation with a variety of simulated 2D and 3D models [6][7][8][9]. Hysteresis in this system shows counterintuitive response to changes in a longitudinal field: increasing the field sweep rate can increase the response to the changing field, as measured by bulk magnetization. Simulations implicate the FM spin chains in the slow dynamics, and indicate the importance of intrachain physics in correctly understanding the magnetic behavior of Ca 3 Co 2 O 6 as a rare example of effective lowdimensional Ising-like triangular frustration. Some models have suggested quantum tunneling within the FM spin chains [10]; the influence of quantum fluctuations on relaxation in spin-chain compounds is therefore an important question.
Two decades ago, Brooke et al. explored the role of quantum fluctuations for the disordered spin glass * aking@dwavesys.com FIG. 1. Geometrically frustrated square-octagonal lattice viewed as a spin-chain antiferromagnet. Red and blue lines indicate antiferromagnetic and ferromagnetic Ising exchange with couplings J1 and J0 = −1.8J1 respectively, where ferromagnetic bonds form four-spin chains. As in Ca3Co2O6, the strong FM chains interact via weaker AFM couplings in a triangular geometry. LiHo x Y 1−x F 4 , annealing it to a low-energy state in a laboratory by attenuating thermal and quantum fluctuations [11]. Their finding that quantum annealing relaxed the system faster than thermal annealing was a major motivating factor in the development of programmable superconducting quantum annealing (QA) processors [12]. In turn, recent experiments have demonstrated that a class of quantum condensed matter sys-tems can be simulated using QA processors [13][14][15][16][17]. One such system is a geometrically-frustrated 2D Ising magnet ( Fig. 1), in which quantum fluctuations induce orderby-disorder (OBD) at zero longitudinal field. Upon increasing temperature, this ordering is suppressed via a sequence of two Berezinskii-Kosterlitz-Thouless (BKT) transitions, a phenomenon previously demonstrated using path-integral quantum Monte Carlo (QMC) simulations [18] and later using QA simulation to probe equilibrium and dynamical properties [14,15].

II. SPIN-CHAIN MODEL AND RELATED COMPOUNDS
In this work we consider the lattice studied in Refs. [14,15] as a spin-chain antiferromagnet (Fig. 1) and simulate it under a longitudinal magnetic field, in addition to the transverse field that induces quantum fluctuations. The Hamiltonian of the system is where ν ∈ {1, 2, 3, 4} is the spin index of each vertical chain, and ν(j, l) is the spin index of the chain j that interacts with a spin in the chain l. The intra-chain interaction is ferromagnetic (J 0 > 0) and the inter-chain interaction is antiferromagnetic (J 1 > 0). The Hamiltonian parameters can be expressed as B(t) = J (t)H(t), J 0 (t) = J (t)I 0 , J 1 (t) = J (t)I 1 and Γ(t) to indicate the time-dependent control over strengths of the Ising, longitudinal field, and transverse field terms respectively, where I 0 = 1.8, I 1 = 1 and J (t) ≥ 0. Note that the ratio between the FM and the AFM Ising interactions is timeindependent J 0 /J 1 = I 0 /I 1 = 1.8. H(t) ranges from 0 to 2. Thus we can alternatively express the time-dependent Hamiltonian as where J (t) and Γ(t) are the time-dependent Ising and transverse energy functions comprising the annealing schedule, familiar in the field of quantum annealing [13]. The choices of I 0 = 1.8 and I 1 = 1 ensure that-as in related materials of interest-spin chains have few breaks at equilibrium, and the lattice can be fully saturated with H = 2. While the spin Hamiltonian H cannot be directly applied to any real material, it can still be used to model the behavior of known quantum magnets. In addition to  [14,27]. [19][20][21][22][23], the compound TmMgGaO 4 has recently become the subject of intense study [24][25][26][27][28]. Here the magnetic Tm 3+ ions form a perfect triangular lattice. These moments can be described in terms of an effective spin-1/2 variable because the two lowest crystal field (CF) levels are separated from the rest by an energy gap that is much bigger than the Ising-like exchange interaction between different magnetic moments. Given that the Tm 3+ ion has total angular momentum J = 6, there is no Kramers degeneracy and the two lowest energy CF levels are singlets [26,27]. The energy splitting between the two singlets corresponds to an intrinsic transverse field acting on the spin-1/2 variables and the low-energy physics of Tm 3+ is described by a transverse field Ising model (TFIM) on a triangular lattice [26,27], which coincides with the effective low-energy model for H in the limit Γ/J 0 1. In this limit, the low-energy degrees of freedom of each chain are also spin-1/2 variables that interact via the effective Hamiltonian:

III. EQUILIBRIUM PROPERTIES
As in [14], we program the QA processor to realize a cylindrical lattice on 1800 spins; in this experiment defects reduce the number of working spins to 1764. To suppress boundary effects and better simulate the thermodynamic limit, we tune longitudinal field and coupling terms on a per-device basis to empirically homogenize qubit magnetizations and coupled spin-spin correlations, respectively [14,29,30].
We estimate statistics of the system at thermal equilibrium under a fixed longitudinal magnetic field B(t) = B(0) that ranges from 0 to 2, using the "quantum evolution Monte Carlo" (QEMC) method introduced in Ref. [14]. In this mode of operation the system is initialized in a random classical spin state, and repeatedly exposed to both quantum and thermal fluctuations, in this case for 500 µs per exposure. Between exposures, the fluctuations are quenched as the system is destructively projected to the computational σ z basis, and a classical state is read (see appendix). Because the Ising system is repeatedly quenched and relaxed, the observed output does not precisely reflect the system at equilibrium. Rather, prior to each readout the quantum and thermal fluctuations are turned off, and a small amount of local relaxation occurs during this process-most obviously in the erasure of single-spin excitations. Bulk observables reflecting long-range order or magnetization are relatively protected from this distortion [15]. One can therefore think of the equilibrium estimates in two ways: first, as a perturbed observation of the system at equilibrium, or second, as an observation of a periodically-driven Floquet system. Figure 2 shows the thermodynamic phase diagram obtained with classical (Γ = 0) and quantum (Γ = 0) Monte Carlo, and the transitions that are extracted from the QA processor for the QEMC runs. In the classical Ising limit (Γ = 0), the phase diagram only includes two phases. The ordered phase at low enough temperature and longitudinal field B = 0 corresponds to a threefold degenerate three-sublattice up-up-down ferrimagnetic (FIM) state that induces a 1/3 plateau in the magnetization curve M (B)/M sat (M = jν σ z jν ) and can be characterized by the complex order parameter: where m 1 , m 2 , and m 3 are individual sublattice magnetizations [18]. The three ground states of the 1/3 plateau have ψ equal to −(2/ √ 3)e k2πi/3 for k = 0, 1, 2. Thus to measure the onset of FIM long-range order (LRO) in our simulations of finite size lattices, we use the normalized cubic invariant The second phase is simply the paramagnetic (PM) state that becomes a fully polarized state along the longitudinal field direction at T = 0. As shown in Figure 2, the FIM phase survives for Γ = 0. In addition, there are two low-temperature phases identified in the B = 0 case: a critical phase and an ordered phase, in which ψ concentrates around six values rather than three [31]. This is also a three-sublattice ordering in which one of the two spin-up sublattices of the FIM state is polarized along the x-direction. The critical phase disappears for finite B with the two BKT transitions that mark its upper and lower boundary. For B Γ, the system transitions directly from the PM state to a sixfold degenerate state upon cooling, with the spins of the third sublattice canted along the longitudinal field direction. Unfortunately, the base temperature of the QA processor is not low enough to access this ordered phase [14].
The phase transition points between the FIM and the PM phases can be determined with the QA processor (blue and red squares in Figure 2) by measuring the crossing point of m FIM = 1/2. The transition is accompanied by a rather abrupt change of the magnetization that produces a peak in the longitudinal susceptibility dM/dB. As shown in Fig. 3, this peak splits into two upon increasing the temperature and the transverse field. The lower field peak still signals the metamagnetic phase transition between the FIM and the PM phases. The higher field peak results from the combined effect of quantum and thermal fluctuations. As such, this second peak in dM/dB, or the corresponding magnetization shoulder in M (B) (see Fig. 3), can be used as an experimental fingerprint of the presence of a transverse field. In other words, the double-peak structure can be used to identify materials that are realizations of quantum Ising models.
To understand the origin of the second peak, it is convenient to start from theΓ = 0 limit ofH. In this classical limit,H has an extensive ground state degeneracy at the effective saturation fieldB sat = 4B sat = 6J, which is lifted by the transverse field: the ground space S is generated by spin states |φ j that do not contain any pair of nearest-neighbor spins anti-aligned with the longitudinal field B. To first order inΓ, the new ground state is obtained by diagonalizing theΓ-term restricted to S. Given that theΓ-term is irreducible on S and all the off-diagonal terms are semi-negative defined, The extensive ground state entropy reappears at a finite temperature T Γ and induces a magnetization shoulder near the middle of the magnetization ramp that is clearly visible in Fig. 3 (b). A numerical estimate of the magnetization value that maximizes the number of states |φ j gives M m /M sat 0.583, which is in good agreement with the value of the magnetization shoulder shown in Fig. 3 (b). Further evidence for the entropic origin of the shoulder is provided in the appendix. Based on the entropy argument, the width of the shoulder is expected to be proportional to T and to Γ (note that δB Γ ). In other words, the shoulder is the combined effect of quantum fluctuations, that stabilize the low energy manifold S over a finite longitudinal field interval of order Γ, and thermal fluctuations that select the submanifold of states with magnetization M M m . We note that an intermediate regime between the FIM and the fully polarized phases has been detected with neutron scattering measurements of TmMgGaO 4 [28]. However, it is difficult to separate the relative roles of the intrinsic quantum and thermal fluctuations of the clean limit versus the spatial fluctuations induced by the significant amount of disorder present in TmMgGaO 4 .
Ca 3 Co 2 O 6 comprises a triangular lattice of ferromagnetic (FM) Ising chains coupled by weak antiferromagnetic (AFM) exchange interactions. This compound exhibits longitudinal field-induced magnetization steps whose heights depend on the field sweep history and rate [19][20][21][22][23]. These out-of-equilibrium magnetization steps, which appear at regular magnetic field intervals, originate from the extensive ground state degeneracy of the triangular Ising model and by the need to overcome an energy barrier to connect different ground states by a sequence of individual spin flips (local dynamics). As it was pointed out in Ref. [32], the regular spacing between the magnetization steps can be explained by the regular spacing between the molecular fields produced the six nearest-neighbor spins that surround the spin to be flipped. Experimental studies of the quantum effect of an external transverse field in Ca 3 Co 2 O 6 are challenging because of the very small effective gyromagnetic factor in the direction perpendicular to the chains.

IV. OUT-OF-EQUILIBRIUM BEHAVIOR
We simulate relaxation of the system under field B(t) that either increases from B(0) = 0 or decreases from a saturating field B(0) = 2J 1 . Fig. 4a shows the timedependent Hamiltonian terms for a 100 µs field sweep, which follows a 100 µs anneal to the desired values of Γ and J 1 (cf. [13] Fig. S13). Since QA readout is achieved after projection to the σ z basis following a rapid quench of both Γ and J 1 , each value of B/J 1 must be simulated individually, as opposed to taking multiple measurements from a single sweep of B. We probe sweep rates of 1 µs, 10 µs, and 100 µs for low and high transverse field (Γ/J 1 = 0.51 and 1.06); J 1 = 1.89 GHz is roughly 4.5 times the effective qubit temperature T = 10 mK. As the sweep of B becomes slower, observations approach equilibrium values, which we estimate using a quantum evolution Monte Carlo (QEMC) protocol in which the system is relaxed iteratively [14], in contrast to the singleshot measurements. Fig. 4b shows magnetization hysteresis curves for different sweep rates. As in investigations of Ca 3 Co 2 O 6 , we observe distinct out-of-equilibrium behavior, including metastable overshooting of the equilibrium 1/3 plateau. Increasing the transverse field Γ/J 1 from 0.51 to 1.06 significantly reduces the signatures of metastability in the magnetization curves.
The 1 µs sweep of B results in almost no LRO when Γ/J 1 = 0.51 despite the order seen at thermal equilibrium (Fig. 4c). Throughout the curve, increasing Γ/J 1 to 1.06 appears to hasten the development of LRO by over an order of magnitude. In this sense it is clear that quantum fluctuations suppress metastable trapping in the out-of-equilibrium experiment. Like in the case of Ca 3 Co 2 O 6 [33], the out-of-equilibrium simulation of 1 µs sweep also shows three steps, in addition to the above-mentioned shoulder. It is important to note that H is far from being a realistic model for Ca 3 Co 2 O 6 . Besides the fact that Ca 3 Co 2 O 6 is a 3D material, a more realistic Hamiltonian should include multiple competing inter-chain exchange interactions that induce a commensurate-incommensurate transition at finite temperature [6]. Nevertheless, as it was shown by Kudasov [32], a simplified 2D triangular Ising model can explain the observation of three out-of-equilibrium steps, so it not surprising that a time-evolution controlled by H can reproduce the three steps. Near-term QA processors will be able to simulate more flexible geometries [34], including longer spin chains.
The structure of output states can shed light on the effect of quantum fluctuations at equilibrium and out-ofequilibrium. Fig. 5 shows representative output states from the QA simulation. Fig. 5a and b show equilibrium states at B/J 1 = 1.48 for Γ/J 1 = 0.51 and Γ/J 1 = 1.06, respectively. We see clear evidence of the transverse field's role in suppressing FIM LRO while remaining in the S manifold. In Fig. 5c and d we see outof-equilibrium states for the fastest sweep of B/J 1 from 2 to 1. In this case we see far more FIM LRO in the high-Γ case (in line with Fig. 4c), and many competing FIM domains in the low-Γ case. At the domain boundaries we see a variety of trapped defects, including up-down-down plaquettes and broken spin chains with nearly zero average magnetization, that contribute to the overshooting of the 1/3 plateau seen in Fig. 4b. The presence of broken chain defects indicates that the out of equilibirum dynamics is not fully captured by the low-energy effective model [32].
It is not clear to what extent this trapping is escaped versus avoided as the system passes through the phase transition either as B increases from 0 or decreases through B c ≈ 3 2 J 1 . Although path-integral Monte Carlo (PIMC) simulation is not expected to correctly reproduce open-system relaxation dynamics of complex frustrated systems [15,35], such simulations suggest that the acceleration arising from the transverse field is dynamical, and not reflective of a new path through the phase transition that avoids metastable configurations (see appendix).

V. CONCLUSIONS
We have demonstrated a large-scale out-of-equilibrium simulation of a frustrated transverse field Ising model using a programmable superconducting QA processor. When quantum fluctuations are small, we observe the signature metastable magnetization curves also seen in experimental measurements of spin-chain compounds. Introducing larger quantum fluctuations produces both dynamical and equilibrium effects on the system. Outof-equilibrium simulations are a promising application for quantum simulation: standard quantum Monte Carlo methods-be they single-spin Glauber dynamics or imaginary-time cluster update dynamics (e.g. Swendsen-Wang)-cannot faithfully reproduce the the real timeevolution of a quantum system [15,35,36]. Although evolution of a quantum system in the presence of a thermal bath is fundamental to QA, the relationship between QA dynamics and time-evolution of a quantum Ising model has not been established experimentally in detail for large systems. Moreover, the spin-chain system we have simulated is only qualitatively similar to real compounds-future hardware generations with greater control and more flexible Hamiltonian terms and geometry [34,37] will allow deeper and more realistic simulations of this nature.
Although independent control of transverse and longitudinal magnetic fields is difficult to achieve for in situ crystal studies, a laboratory reproduction of our findings in Ca 3 Co 2 O 6 would be of great interest. This would mark the first example of a quantum simulation preceding confirmation in a physical sample, an important milestone in quantum simulation for materials science. ; the transverse field modulates long-range order. a, at Γ/J1 = 0.51 the system adopts long-range FIM order and the system is dominated by plaquettes of a single phase, with mFIM large. b, with Γ/J1 = 1.06, the system is characterized by competing FIM domains, with saturated plaquettes at the interfaces. c-d, out-of-equilibrium samples, with H decreasing from 2 to 1 over 0.5 µs. c, with Γ/J1 = 0.51 the system is chaotic, with many small FIM domains. The interfaces of these domains host trapped excitations, many in the form of broken four-qubit chains with magnetization zero. d, with Γ/J1 = 1.06 the system forms large FIM domains during the sweep of H. The interfaces between these domains host plaquettes of magnetization +1 and −1/3. For estimating equilibrium statistics, we perform "quantum evolution Monte Carlo", forming a chain of reverse anneals. Each reverse anneal starts from a classical state with Γ = 0, then increases Γ, dwells, and decreases Γ to 0, arriving at a new classical state. In the easy case (left) the protocol involves a simple dwell. In the slow-relaxing case (right) Γ is attenuated from a higher value before the dwell. Each chain repeats the reverse anneal 100 times, producing 100 classical output samples that ideally represent projections of the system to the σ z basis. In both cases, B(t)/J1 is held constant (shown here for example B(t)/J1 = 1), in contrast to the out-of-equilibrium protocols shown in Fig. 4. of [15] for a different superconducting processor of similar design. Similarly, the leading qubit-to-spin nonideality of background susceptibility is compensated in the coupling energies as in Methods of [15].
In this work we study a single lattice programmed into the qubit architecture; this has the same 1800-qubit size as studied in [14] but due to inoperable qubits only uses 1764 spins. Longitudinal fields H = B/J 1 can be programmed either negative or positive; we run experiments for H = 0, 0.04, . . ., 2 in both directions and combine symmetric results.
With each programming of the quantum processing unit (QPU) we draw 100 samples. Out-of-equilibrium experiments are performed using the protocol described in the main text and shown in Fig. 4; these protocols require the "H-gain" control. Equilibrium estimates are generated using quantum evolution Monte Carlo (QEMC) [14,15,30], a method that, analogous to Markov-chain Monte Carlo, forms a Monte Carlo chain of reverse annealing protocols that dwell at the desired Hamiltonian rather than annealing to it. The QEMC (equilibrium estimating) calls also draw 100 samples, of which we discard the first 50 as QEMC burn-in.We dwell for 500 µs per step, except in the case where Γ/J 1 < 0.77, where relaxation is slow. In this case we reverse anneal to an intermediate value of Γ/J 1 , then anneal "forward" (i.e., attenuating Γ) to the target value of Γ/J 1 over 250 µs, then dwell at the target Γ/J 1 for the remaining 250 µs. Fig. 7 shows the waveforms for a single reverse annealing step in each of these cases.
In practice the relative strengths of Γ and J 1 are not controlled independently, but rather through single annealing parameter s; the functions Γ(s) and J 1 (s) are fixed functions of the QA processor, which allows timedependent control of s(t) (cf. [15] Methods). Varying Γ for fixed J 1 therefore requires tuning Γ(s) and J 1 (s) via s, then compensating by adjusting the individual programmed coupling terms.
All reported data are taken over 120 individual programmings, representing 12, 000 samples for out-ofequilibrium and 6, 000 samples for equilibrium. Reported bulk averages (magnetization, susceptibility, FIM order) are taken over non-boundary spins to reduce the effect of open boundaries and defects. eters calibration to compensate for device variation and boundary conditions [14,15,30,38] (specifically cf. [14] Extended Data Fig. 7). The idea of per-spin tuning of longitudinal fields in the presence of open boundaries goes back to the 1970s [29], and although it now seems obsolete for computer simulations, here it allows us to better simulate the thermodynamic limit of a large system in a uniform magnetic field [30], in spite of the fact that limited qubit connectivity prevents us from implementing fully periodic boundaries.
In contrast to previous studies on cylindrical lattices [14,15] that attempted to faithfully simulate a finite system with partially open boundaries, here we are interested in the thermodynamic limit. We therefore use the same approach of iteratively tuning J ij terms to homogenize spin-spin correlations. However, instead of doing this over small rotational symmetry groups of the cylinder, we homogenize correlations and magnetizations across the entire lattice.
When homogenizing the coupling terms we tune terms using H = B/J 1 = 0 observations and propagate the couplings to all simulated H values. When homogenizing the per-spin field strengths we use all observations and employ smoothing of the Hamiltonian terms so that relative field strengths vary smoothly as a function of H.

Appendix B: Entropic origin of magnetization shoulder
Further evidence for the entropic origin of the shoulder is shown in Fig. 6, where we compute the average local degeneracy of observed states, defined as the number of states reachable from a given state by moving a single down-spin chain to a neighboring site. This local entropy is found to peak at the phase transition, with peak height that increases with both quantum and thermal fluctuations, in the form of changing Γ and T respectively.