Phase transitions in a programmable quantum spin glass simulator

See allHide authors and affiliations

Science  13 Jul 2018:
Vol. 361, Issue 6398, pp. 162-165
DOI: 10.1126/science.aat2025

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


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 (29). 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 Jij be –1 (ferromagnetic) and +1 (AFM), respectively. The Hamiltonian for this system can be written asEmbedded Image(1)where Embedded Image and Embedded Image are the Pauli z and x operators, respectively, acting on spin i, and the notation Embedded Image 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 pc = 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 Embedded Image at p = 0 (16) and Embedded Image at p = 0.5 (17), where kB 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 Embedded Image at p = 0 (18) and Embedded Image at p = 0.5 (19) indicate locations of quantum phase transitions between the PM-AFM and PM-SG phases, respectively, at T = 0.

Fig. 1 The bond-disordered TFIM on a simple cubic lattice.

(A) Schematic phase diagram. Axes correspond to scaled transverse magnetic field splitting Embedded Image, disorder p, and scaled temperature Embedded Image, where J is the Ising interaction strength. The experiments reported on here probed this space on a curved manifold (translucent surface). Any single experiment for fixed p involves annealing along a trajectory on this manifold (black dotted curve) that starts in the PM phase, and a critical point is encountered where the trajectory penetrates a phase boundary (silver point). (B) Illustration of a ground state at Γ = 0 for an 8 × 8 × 8 lattice instance for p = 0.1 < pc. Ferromagnetic (Jij = –1) and AFM (Jij = ±1) bonds are indicated by gold and silver cylinders, respectively. Ising spin orientations ↑ and ↓ are represented by red and blue spheres, respectively.

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 Embedded Image and Embedded Image. For fixed temperature T ≈ 12 mK (kBT/h ≈ 250 MHz), the QPU can be used to explore the phase space on a curved manifold in Embedded Image–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 pc, 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 Jij, 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) (2022). 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 Embedded Image, where N = L3 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)Embedded Image(2)for each state Embedded Image, where the vector of integers (ai, bi, ci) corresponds to the position of site i in the SC lattice (relative to an arbitrary origin), and si is the state of the spin at site i. A summary plot of Embedded Image versus p, where Embedded Image 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, Embedded Image monotonically decreases with increasing p, and the slope near p = pc becomes increasingly steep with larger L. However, this quantity does not strictly go to zero for p > pc because of finite size effects (23). To determine the experimentally observed critical point Embedded Image, we calculated the Binder cumulant (24) Embedded Image(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 pc 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 Embedded Image. A collapse of the data are shown in Fig. 2C for a posited scaling function Embedded Image, from which we obtain a critical exponent 1/vexp = 0.85 ± 0.05. These experimental results are in good agreement with the numerical results of pc = 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).

Fig. 2 Determination of the critical disorder pc and AFM susceptibility χAFM versus quantum annealing parameter s.

(A) Disorder-averaged AFM magnetization Embedded Image versus p for three different lattice sizes L × L × L. (B) Binder cumulant g versus p. For each L, the results near pc have been independently fit to second-order polynomials (solid curves). The critical disorder is given by the mutual crossing point, which is estimated to be Embedded Image (dashed vertical line). (C) Finite-size scaling collapse of the experimental results. Assuming that g is a function of scaled disorder Embedded Image, the best collapse was obtained for 1/vexp = 0.85 ± 0.05. (D) χAFM versus s for p < pc. (E) χAFM versus s for p > pc. (F) Normalized χAFM versus s for p > pc.

We next turned to the identification of phase transitions between the PM, AFM, and SG phases as a function of Embedded Image and Embedded Image. 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 ≡ δmAFMh|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 < pc = 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 sc. For p > pc = 0.22, sc 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(sc) versus s in Fig. 2F.

The peaks observed for p < pc suggest an incipient divergence, which reliably indicates a second-order phase transition between the PM and AFM phases. For p > pc, 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 Embedded Image, as given in Eq. 2, is an order parameter that one can use to easily identify the AFM phase. However, Embedded Image 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 asEmbedded Image(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 < sc, which is consistent with the PM phase in a finite-sized system. For the AFM example in Fig. 3A, P(q) bifurcates for s > sc, 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 Z2 symmetry (28). For the SG example in Fig. 3B, P(q) broadens in the vicinity of sc 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).

Fig. 3 Example spin-spin overlap probability distributions P(q).

(A) Example PM-to-AFM phase transition for a p = 0.1 instance at four representative values of s. The results for s = 0.1 have been magnified ×10. (B) Example PM-to-SG phase transition for a p = 0.5 instance at four representative values of s.

The final analysis presented here is the phase diagram projected on to the Embedded Image–plane for the manifold depicted in Fig. 1A. To convert sc inferred from linear susceptibility peak positions to Embedded Image, 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 pc from the literature (15), as well as values of Embedded Image (18) and Embedded Image(19) for an infinitely large system at T = 0. The experimentally determined phase diagram projection yields a value of pc 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 Embedded Image–plane. Results for p = 0 are shown in Fig. 4B. Decreasing or increasing J moves the phase boundary toward the classical transition point Embedded Image (16) or quantum critical point Embedded Image (18), respectively. Fitting the rightmost edge of the phase boundary to a power law Embedded Image yielded an experimental quantum critical point Embedded Image and exponent aexp = 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.

Fig. 4 Embedded Image and Embedded Image phase diagrams.

(A) Embedded Image phase diagram computed from χAFM peak positions (black points with error bars). All measurements performed at T = 12.0 mK mK. Thick dashed black line indicates experimentally determined critical disorder Embedded Image, and the thin dashed black lines indicate error estimates for that quantity. Critical transverse fields Embedded Image and Embedded Image for an infinite system at T = 0 and critical disorder pc were taken from the literature. (B) χAFM peak positions for varying relative coupling strength J translated into the Embedded Image–plane for p = 0. χAFM(s) data are provided in fig. S13. Annealing trajectories for varying relative coupling strength are indicated by dashed curves (supplementary materials), and the corresponding transition points are indicated by color-coordinated points, with error estimate along each trajectory. Trajectory and critical point denoted as 1.00× correspond to the experimental manifold depicted in Fig. 1A and the p = 0 point in (A), respectively. Trajectories and points labeled 0.33×, 0.50×, and 0.67× indicate lower J, whereas those labeled 1.33×, 1.67×, 2.00×, 2.33×, and 2.67× indicate higher J. Best fit of phase boundary and extrapolation to T = 0 is indicated by solid black curve. Dashed black curves indicate error bounds on the best fit model. Open circles represent the theoretical results (30) for an infinite system.

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

Materials and Methods

Figs. S1 to S15

Table S1

References (3444)

Movie S1

Data Files

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.

Stay Connected to Science

Navigate This Article