## Simulating correlated electron systems

Correlated electron systems are generally difficult to simulate because of limited capabilities of computational resources. Harris *et al.* used a D-Wave chip based on a large array of superconducting elements to simulate the phases of a complex magnetic system. They tuned the amount of frustration within the lattice and varied the effective transverse magnetic field, which revealed phase transitions between a paramagnetic, an ordered antiferromagnetic, and a spin-glass phase. The results compare well to theory for this spin-glass problem, validating the approach for simulating problems in materials physics.

*Science*, this issue p. 162

## Abstract

Understanding magnetic phases in quantum mechanical systems is one of the essential goals in condensed matter physics, and the advent of prototype quantum simulation hardware has provided new tools for experimentally probing such systems. We report on the experimental realization of a quantum simulation of interacting Ising spins on three-dimensional cubic lattices up to dimensions 8 × 8 × 8 on a D-Wave processor (D-Wave Systems, Burnaby, Canada). The ability to control and read out the state of individual spins provides direct access to several order parameters, which we used to determine the lattice’s magnetic phases as well as critical disorder and one of its universal exponents. By tuning the degree of disorder and effective transverse magnetic field, we observed phase transitions between a paramagnetic, an antiferromagnetic, and a spin-glass phase.

Can quantum physics be efficiently simulated on a computer? What kind of computer would one need to do so? Feynman posed these questions to motivate the concept of a probabilistic quantum mechanical computer that could be used to simulate nature (*1*). Whereas the implementation of such computers seemed daunting in Feynman’s time, subsequent progress in the laboratory has made small-scale quantum simulation possible today (*2*–*9*). Although the field of quantum simulation is still in its infancy, the scale and the control of engineered quantum systems has reached the point where they can begin to address fundamental problems in physics.

In this Report, we describe the use of an engineered quantum system to investigate the properties of a particular quantum magnetism problem known as the transverse field Ising model (TFIM). The TFIM is often viewed as the archetypal model for order-disorder phase transitions in quantum systems (*10*). A diverse collection of condensed matter systems, such as ferroelectrics (*11*) and magnetic spin glasses (*12*), can be described by the TFIM. Although this model can be analytically solved in one dimension, it is generally only amenable to numerical study via computationally intensive means, such as quantum Monte Carlo (*13*), in higher dimensions. We have simulated the TFIM on three-dimensional (3D) simple cubic (SC) lattices up to dimensions 8 × 8 × 8 by embedding the model on a D-Wave quantum processing unit (QPU) (D-Wave Systems, Burnaby, Canada) (*14*). By tuning the amount of disorder within the lattice and varying the effective transverse magnetic field, we demonstrate phase transitions between a paramagnetic (PM), an ordered antiferromagnetic (AFM), and a spin-glass (SG) phase. The experimental results compare well with theory for this particular SG problem, thus validating the use of a probabilistic quantum computer to simulate materials physics. This represents an important step forward in the realization of integrated quantum circuits at a scale that is relevant for condensed matter research.

The system studied here is a SC lattice of Ising spin-1/2 particles subject to a global transverse field splitting (*10*). The single-spin energy splitting induced by the transverse field is denoted by Γ, and the magnitude of the nearest-neighbor interaction energy is denoted by J. To induce glassy behavior, let the sign of the interaction vary randomly, with *p* and 1 – *p* denoting the probability that a particular interaction *J _{ij}* be –1 (ferromagnetic) and +1 (AFM), respectively. The Hamiltonian for this system can be written as(1)where and are the Pauli

*z*and

*x*operators, respectively, acting on spin

*i*, and the notation indicates summing over all nearest neighbor pairs. A schematic phase diagram for this system is depicted in Fig. 1A. For Γ = 0, this is a well-studied classical system, with phase transitions induced by thermal fluctuations. The critical disorder

*p*

_{c}= 0.222 (

*15*) marks the value of

*p*below which one observes an ordered AFM phase, whereas above this value one observes a SG phase. Critical temperatures at

*p*= 0 (

*16*) and at

*p*= 0.5 (

*17*), where

*k*

_{B}is the Boltzmann constant, indicate the locations of PM-AFM and PM-SG thermal phase transitions, respectively, at Γ = 0. As Γ is increased, phase boundaries evolve because of the increased role that quantum fluctuations play over thermal fluctuations. The critical transverse field energies at

*p*= 0 (

*18*) and at

*p*= 0.5 (

*19*) indicate locations of quantum phase transitions between the PM-AFM and PM-SG phases, respectively, at

*T*= 0.

The SC lattice described above was embedded into the native graph of a prototype QPU (*14*). The QPU and the embedding are described in the supplementary materials. Lattices up to size 8 × 8 × 8 fit into the QPU, an example of which is shown in Fig. 1B. Within this experimental platform, the coupling energy J and transverse field splitting Γ are monotonic functions of a dimensionless control parameter 0 ≤ *s* ≤ 1 so that and . For fixed temperature *T* ≈ 12 mK (*k*_{B}*T*/*h* ≈ 250 MHz), the QPU can be used to explore the phase space on a curved manifold in –space, as indicated in Fig. 1A. Details concerning the functional form of this manifold can be found in the supplementary materials.

To experimentally determine the critical disorder *p*_{c}, we first worked in a classical regime by noting that the final Hamiltonian at *s* = 1 is given by the classical form of Hamiltonian (Eq. 1) with Γ(*s* = 1) → 0. Ensembles of up to 1000 random instances, with each instance corresponding to a specific set of *J _{ij}*, were generated for a series of values of

*p*and lattice sizes

*L*×

*L*×

*L*. Each instance was programmed into the QPU, and low-energy states of the final classical Hamiltonian were obtained by slowly evolving from

*s*= 0 to

*s*= 1 in a process known as quantum annealing (QA) (

*20*–

*22*). Experimental details are available in the supplementary materials. QA was performed 1000 times per instance, and the final states of every spin were recorded. In order to approximately estimate ground-state properties, the results were then filtered in order to identify the lowest-energy state returned by the QPU for each instance. In cases in which the lowest observed energy was degenerate, all such states were kept. The result was a set of low-energy states , where

*N*=

*L*

^{3}is the number of lattice spins, for each value of

*p*that could then be used to determine disorder-averaged quantities, as described below.

To detect the AFM-SG phase transition as a function of *p*, we calculated the absolute value of the AFM magnetization (also known as staggered magnetization)(2)for each state , where the vector of integers (*a _{i}*,

*b*,

_{i}*c*) corresponds to the position of site

_{i}*i*in the SC lattice (relative to an arbitrary origin), and

*s*is the state of the spin at site

_{i}*i*. A summary plot of versus

*p*, where denotes averaging over disorder, is shown in Fig. 2A. Data are shown for three system sizes

*L*×

*L*×

*L*, for

*L*∈ [4, 6, 8]. For all system sizes, monotonically decreases with increasing

*p*, and the slope near

*p*=

*p*

_{c}becomes increasingly steep with larger

*L*. However, this quantity does not strictly go to zero for

*p*>

*p*

_{c}because of finite size effects (

*23*). To determine the experimentally observed critical point , we calculated the Binder cumulant (

*24*) (3)versus

*p*for

*L*∈ [4, 6, 8], because this quantity should be scale invariant at a critical point. The results of this analysis are shown in Fig. 2B, with independent fits of the data near

*p*

_{c}to second-order polynomials for each value of

*L*. From these fits, we estimate the mutual crossing point for all three system sizes to be . A collapse of the data are shown in Fig. 2C for a posited scaling function , from which we obtain a critical exponent 1/

*v*

_{exp}= 0.85 ± 0.05. These experimental results are in good agreement with the numerical results of

*p*

_{c}= 0.222 and 1/

*v*= 0.9 reported in literature (

*15*). This observation confirms that the embedded classical problem is in the same universality class as the 3D cubic lattice problem in the limit of vanishing transverse field (

*25*).

We next turned to the identification of phase transitions between the PM, AFM, and SG phases as a function of and . Unlike the above section, we worked in the quantum regime, where Γ is nonzero. The PM-to-AFM phase transition is known to be of second order (*18*), thus motivating the use of a linear susceptibility measurement that would diverge at a critical point in the thermodynamic limit *L* → ∞. Because the ordered phase at *p* = 0 is AFM, we applied a global longitudinal magnetic field that was uniform in magnitude but alternating in sign between nearest-neighbor sites (±*h*). We refer to the quantity χ_{AFM} ≡ δ*m*_{AFM}/δ*h*|_{h}_{→0} as the AFM or staggered linear susceptibility. χ_{AFM} was measured as a function of QA control parameter 0 < *s* < 1 for small sets of randomly chosen instances for a series of values of *p*. Disorder-averaged traces of χ_{AFM} versus *s* for 8 × 8 × 8 lattices are presented in Fig. 2, D to F. For *p* < *p*_{c} = 0.22, this quantity yields obvious peaks that broaden and shift to increasing *s* with increasing *p* in Fig. 2D. In what follows, we refer to the peak positions as *s*_{c}. For *p* > *p*_{c} = 0.22, *s*_{c} remains fixed to within experimental resolution, but the peak amplitude continues to decrease with increasing p, as shown in Fig. 2E. Furthermore, the peak broadens until *p* = 0.4, beyond which the breadth remains fixed. This latter observation is made more clear in the scaled plot of χ_{AFM}(*s*)/χ_{AFM}(*s*_{c}) versus *s* in Fig. 2F.

The peaks observed for *p* < *p*_{c} suggest an incipient divergence, which reliably indicates a second-order phase transition between the PM and AFM phases. For *p* > *p*_{c}, the absence of such a divergence is in agreement with the results of numerical simulations of the PM-to-SG phase transition in a related system (*26*). We consider this to be only part of the evidence needed to confirm that a PM-to-SG phase transition has been observed.

The AFM magnetization , as given in Eq. 2, is an order parameter that one can use to easily identify the AFM phase. However, in both the PM and SG phases. The latter two phases may be distinguished by using an order parameter proposed by Edwards and Anderson (*27*), hereafter referred to as the spin-spin overlap *q*, which is defined as(4)where *A* and *B* are independently determined states for a given instance of disorder. Given a statistically large number of such states, one can construct the probability distribution *P*(*q*), whose characteristics dramatically differ between the PM, AFM, and SG phases. Examples of experimentally determined *P*(*q*) for an AFM (*p* = 0.1) and a SG (*p* = 0.5) instance of size 8 × 8 × 8 at four values of *s* are shown in Fig. 3, A and B, respectively. In both cases, the results reveal a Gaussian-shaped *P*(*q*) centered at *q* = 0 for *s* < *s*_{c}, which is consistent with the PM phase in a finite-sized system. For the AFM example in Fig. 3A, *P*(*q*) bifurcates for *s* > *s*_{c}, and all of the weight becomes consolidated in the two sharp peaks at *q* = ±1 above the transition. This is the signature of an ordered phase, in which there is a single thermodynamic state up to *Z*_{2} symmetry (*28*). For the SG example in Fig. 3B, *P*(*q*) broadens in the vicinity of *s*_{c} and then develops multiple peaks that are approximately symmetric around *q* = 0 and continuous throughout. This is consistent with the existence of many competing thermodynamic states, as typically characterizes a disordered SG phase (*28*).

The final analysis presented here is the phase diagram projected on to the –plane for the manifold depicted in Fig. 1A. To convert *s*_{c} inferred from linear susceptibility peak positions to , we used the lattice spin energy scale versus *s* mapping given in the supplementary materials. The results of this mapping are shown in Fig. 4A along with the phases identified based on the order parameter analysis presented above. For reference, we have also included *p*_{c} from the literature (*15*), as well as values of (*18*) and (*19*) for an infinitely large system at *T* = 0. The experimentally determined phase diagram projection yields a value of *p*_{c} at Γ = 0 that agrees with numerical results. As of the time of writing, we had not established an experimental method to resolve any possible reentrant behavior of the AFM-SG phase boundary in the (*p*, Γ)–plane, as is suspected to exist in the (*p*, *T*)–plane (*29*). Critical transverse fields are systematically lower than the *L* → ∞ results at *T* = 0 because we are constrained to measure on a particular manifold, as depicted in Fig. 1A. The intersection between this experimentally accessible manifold and the phase boundaries can be varied by altering the annealing trajectory, which can be achieved by adjusting the magnitude of the interaction energy J by an overall scale factor (with further details given in the supplementary materials). For fixed *p*, one can measure the χ_{AFM} peak position, similar to that for *p* = 0 depicted in Fig. 2D, as a function of relative scaling of J. Doing so samples the PM-AFM phase boundary in the –plane. Results for *p* = 0 are shown in Fig. 4B. Decreasing or increasing J moves the phase boundary toward the classical transition point (*16*) or quantum critical point (*18*), respectively. Fitting the rightmost edge of the phase boundary to a power law yielded an experimental quantum critical point and exponent *a*_{exp} = 0.7 ± 0.3. We have also included the theoretical results of Elliott and Wood (*30*) for an infinite-sized cubic lattice in Fig. 4B.

The experiments reported on here demonstrate how an in situ programmable TFIM in an integrated circuit can be used as a materials physics simulator. The ability to manipulate individual spins and bonds allows one to explore order parameters in ways that are not possible in bulk condensed matter systems. This work opens many potential avenues of research within the field of quantum magnetism, with the possibility of studying lattices with built-in geometrical frustration (*31*, *32*), lattice gauge theories (*33*), and defect/domain formation.

## Supplementary Materials

www.sciencemag.org/content/361/6398/162/suppl/DC1

Materials and Methods

Figs. S1 to S15

Table S1

Movie S1

Data Files

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We thank the technical staff of D-Wave Systems for making these experiments possible; W. Bernoudy for illustrations and animations of the embedding; and A. King, J. Raymond, A. P. Young, and H. G. Katzgraber for useful discussions and for reviewing this manuscript.

**Funding:**This work was wholly funded by D-Wave Systems.

**Authors contributions:**All authors contributed to the conceptualization of the quantum processing unit used in this study. K.B., P.B., C.E., R.M., N.T., and M.V. provided resources by generating the detailed design of the circuit. S.H., E.L., N.L., T.O., and J.Y. provided resources by directly participating in and directing the fabrication of the circuit. R.L., T.M., and R.N. provided resources by preparing and qualifying the circuit. R.H., Y.S., A.J.B., M.R., F.A., C.D., E.H., M.W.J., T.L., I.Pa., I.Pe., G.P.-L., C.R., L.S., and J.W. provided resources by calibrating and characterizing the circuit. R.H. and Y.S. performed the investigation, validation, and data curation. R.H., Y.S., M.H.A., and A.S. contributed to the formal analysis. R.H. and Y.S. wrote the original draft of the manuscript. R.H., Y.S., A.J.B., F.A., M.H.A., E.H., E.L., T.L., and A.S. contributed to the revised manuscript.

**Competing interests:**All authors are employees and holders of stock options in D-Wave Systems. R.H. and K.B. are inventors of U.S. patent application no. 15/881260.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper or the supplementary materials.