Research Article

Emergence and Frustration of Magnetism with Variable-Range Interactions in a Quantum Simulator

See allHide authors and affiliations

Science  03 May 2013:
Vol. 340, Issue 6132, pp. 583-587
DOI: 10.1126/science.1232296

Magnetic Frustration

The study of magnetic frustration has a long history in solid-state physics, but cold-atom systems now offer the possibility of simulating the problem with exquisite control. Islam et al. (p. 583) study a system of 16 trapped ions, using the Coulomb interactions between the ions to simulate exchange interactions present in magnetic systems. The measured spin correlations provide a window into the behavior of the system, as the effective magnetic field and the range of the interactions are tuned.


Frustration, or the competition between interacting components of a network, is often responsible for the emergent complexity of many-body systems. For instance, frustrated magnetism is a hallmark of poorly understood systems such as quantum spin liquids, spin glasses, and spin ices, whose ground states can be massively degenerate and carry high degrees of quantum entanglement. Here, we engineer frustrated antiferromagnetic interactions between spins stored in a crystal of up to 16 trapped 171Yb+ atoms. We control the amount of frustration by continuously tuning the range of interaction and directly measure spin correlation functions and their coherent dynamics. This prototypical quantum simulation points the way toward a new probe of frustrated quantum magnetism and perhaps the design of new quantum materials.

Predicting the behavior of many-body quantum materials such as frustrated magnets is difficult because the number of relevant configurations scales exponentially with the number of particles (13). Feynman proposed the use of a quantum simulator for the task. Here, interactions are engineered in a “standard” quantum system to illuminate the physics behind the real material (4). Cold atoms are excellent standards for quantum simulations of magnetism, with the ability to tailor frustrated magnetic interactions and measure the individual atomic spins (5, 6). Neutral atomic systems are typically limited to nearest-neighbor interactions (7), although geometrically frustrated interactions can be realized in certain optical lattice geometries (8). The strong Coulomb interaction between cold atomic ions (9) has led to the realization of long-range Ising couplings between individual trapped ion spins (1014) and the observation of spin frustration and quantum entanglement in the smallest system of three spins (15).

Here, we report the implementation of variable-range antiferromagnetic (AFM) Ising interactions with up to 16 atomic ion spins, using optical dipole forces. We directly measure the emergence and frustration of magnetic ordering through spatially resolved imaging of the ions. The spins are initially polarized along a strong effective magnetic field transverse to the Ising couplings, and when the field is lowered, the spins order themselves according to the characterstics of the Ising interactions. We can increase the level of frustration by increasing the range of the interactions, which results in a more equitable balance of competing interactions and a suppression of magnetic order. The quantum coherence in the system is characterized by reversing the transverse field back to its initial value and comparing the resulting state with the initial state. These experiments present simulations in a nontrivial quantum system that approaches a complexity level at which it becomes difficult or impossible to calculate the spin dynamics.


We simulate the quantum transverse Ising model with long-range AFM interactions in a one-dimensional spin chain, given by the Hamiltonian

Embedded Image (1)

where the Planck constant h is set to 1, Embedded Image are the spin-1/2 Pauli operators for the ith spin (i = 1,2,…N), B is the effective transverse magnetic field, and Jij > 0 is the Ising coupling between spins i and j, falling off with the lattice spacing |i j| approximately as

Embedded Image (2)

where 0 < α < 3 (8).

For B >> Jij on all pairs, the spins are polarized along the effective transverse magnetic field in the ground state |↑yyy …〉 of the Hamiltonian in Eq. 1, where |↑y〉 denotes a spin along the +y direction of the Bloch sphere. As the ratio of B to the Ising couplings is reduced, the system crosses over to an ordered state dictated by the form of the Ising couplings, and the spectrum of energy levels depends on the range of the interactions. For any finite α > 0, we find by direct calculation that the staggered AFM states |↑↓↑↓ …〉 and |↓↑↓↑ …〉 constitute the doubly degenerate ground state manifold at B = 0, with the degeneracy arising from the time reversal or the global spin flip symmetry of the Hamiltonian. Here |↑〉 and |↓〉 are spin states oriented along the Ising or x direction of the Bloch sphere. Thus, the system exhibits nearest-neighbor AFM or Néel ordering at sufficiently low temperatures. When the interactions are uniform over all pairs of spins (α → 0), the system becomes maximally frustrated and the excitation gap (Fig. 1A) closes, leading to a finite entropy density in the ground state. In this case, any spin configuration with a net magnetization of zero (1/2) for even (odd) numbers of spins belongs to the ground state.

Between the limits B = 0 and B >> Jij, the energy spectrum features a minimum gap, whose position and size depends on the degree of frustration in the system, or the interaction range. (The interaction range is defined as the number of lattice spacings ξ where the interaction falls off to 20% of the nearest-neighbor Ising coupling: ξ ≡ 51/α.) Figure 1A shows a few low-lying energy states of the Hamiltonian in Eq. 1 for an interaction range of ξ = 5 (corresponding to α = 1). The first excited eigenstate merges with the ground state for small B/J0 and has the same spin order as the ground state near B/J0 = 0. The critical gap Δc between the ground and the first coupled excited state determines the adiabaticity criterion (16). Figure 1B compares the position (dotted line) and size (solid line) of the critical gap of the Hamiltonian for various ranges. As the range and hence the amount of frustration increases, the critical field is pushed toward zero, and the gap closes. To observe the effects of frustration, reflected in the density of states near the ground state, we quench the system by ramping the effective transverse magnetic field faster than the critical gap (Embedded Image to populate the lowest coupled excited states. The observed spin order depends on the resulting degree of excitation and hence on the level of frustration.

Fig. 1 Theoretical energy spectrum and critical gap in the long-range antiferromagnetic Ising model (Eq. 1) for N = 10 spins.

(A) Low-lying energy states for α = 1 (characteristic range of ξ = 5 sites) as a function of the dimensionless parameter B/J0. The spacing between the ground state at E = E0 and the first coupled excited state (black lines) reaches a bottleneck at a critical value Bc/J0 with critical gap Δc. (B) Theoretical dependence of Bc/J0 (dotted line) and Δc/J0 (solid line) on the range of the interaction. As the interaction range increases, the competing long-range couplings make it easier to create excitations and the critical gap is reduced, so a relatively small effective transverse field can break the spin ordering. Both parameters approach zero asα → 0 or ξ → ∞. Present experiments are performed with parameters in the shaded region.


The spins are realized in a collection of 171Yb+ ions confined in a linear radiofrequency (Paul) trap, with the effective spin-1/2 system represented by two hyperfine “clock” states within each ion |↑z〉 and |↓z〉, separated by the hyperfine frequency νHF = 12.642819 GHz (17). The variable-range AFM Ising interactions are generated by applying off-resonant spin-dependent optical dipole forces (11) that drive stimulated Raman transitions between the spin states while modulating the Coulomb interaction between the ions in a controlled way (18). The effective magnetic field is generated by simultaneously driving coherent transitions between the spin states with a π/2-phase shift with respect to the dipole force beams. At any time, we measure the state of the spins by illuminating the ions with resonant radiation and collecting state-dependent fluorescence on an imager with site-resolving optics (17). From this information, we can extract all spin correlation functions (18).

The quantum simulation begins by optically pumping all spins to the |↓z〉 state and then coherently rotating them all about the x axis of the Bloch sphere to initialize each spin in state |↑y〉 along the effective transverse magnetic field. The Hamiltonian (Eq. 1) is then switched on with an initial field B0 ≈ 5 J0, where J0 (∼1 kHz) is the average nearest-neighbor Ising coupling, thus preparing the spins in the ground state of the initial Hamiltonian with a fidelity better than 97%. The effective magnetic field is ramped down exponentially in time with a time constant of 400 μs to a final value B of the transverse field, but no longer than 2.4 ms overall, to avoid decoherence effects present in the system (see Fig. 2). At this point, the Hamiltonian is switched off, freezing the spins for measurement. We then measure the x or y component of each spin Embedded Image or Embedded Image by first rotating our measurement axes with an appropriate global π/2 pulse similar to the initialization procedure, before capturing the spin-dependent fluorescence on the imager. The experiments are repeated 2000 to 4000 times to measure expectation values of certain spin operators and correlation functions (denoted by 〈... 〉).

Fig. 2 Spin order versus speed of ramp, for N = 10 spins.

The spins are initialized with B/J0 = 5 and the transverse field is ramped exponentially down for six time constants, and the experiment is repeated for various values of time constant τ. Symbols (solid line) indicate the scaled staggered Binder cumulant Embedded Image , versus the total duration 6τ of the ramp measured (expected from theory). As the Hamiltonian evolves more slowly, the observed spin order shows more ground state order, and less excitation for ramp times under ∼ 2.5 ms. For longer times, the spins become disordered, implying external decoherence in the system. Here (and in the following figures) the error bars include statistical fluctuations and estimated detection uncertainties.

Order Parameters and Correlation Functions

From these measurements, we can construct order parameters appropriate for observing low-energy AFM states. Various moments of the staggered magnetization operator Embedded Image distinguish between paramagnetic and AFM order, and also quantify spin flip excitations. In particular, the normalized fourth moment of this magnetization, known as the Binder cumulant (19) gs = 〈(ms−〈ms〉)4〉/〈(ms−〈ms〉)22, scaled to Embedded Image to remove finite size effects, varies from Embedded Image to Embedded Image as the paramagnetic state gives way to AFM order (see Fig. 2). Here Embedded Image is the Binder cumulant in a perfect paramagnetic state of the N spins. We can also form any correlation function of the spins such as the two-point correlation Embedded Image, allowing a direct probe of spin order for each experimental realization. The Fourier transform of this correlation function is the structure function Embedded Image , where Embedded Image is the average correlation over r sites in the chain. The structure function shows spin order versus wave number k, with S(k = π) singling out the presence of the nearest-neighbor Néel AFM order.

Figure 3 shows the onset of AFM ordering in the quantum simulation of the frustrated transverse field Ising model in a system of 10 spins. Two-point spin correlations C1,r between a chosen edge spin and the other spins are presented in Fig. 3B at various stages in the ramp B/J0 = 5, 0.25, and 0.01. For larger transverse magnetic fields, there are no appreciable correlations between the spin components along the Ising direction. As the ratio of B/J0 is lowered, however, a zig-zag pattern emerges, with negative (positive) correlations between spins separated by odd (even) lattice spacings. For B/J0 ≈ 0.01, the nearest-neighbor spin correlation reaches about 60%, and the correlation length (defined to be the distance at which the absolute correlation drops to 1/e of the nearest neighbor value) reaches about six lattice sites. The effective field was ramped exponentially down from B/J0 ≈ 5 with a time constant of 400 μs in this experiment. This ramping is not slow enough to be adiabatic, and the diabatic effects prevent the spin ordering from reaching a perfect AFM phase. Figure 3C shows the measured probabilities of all 210 = 1024 possible spin states measured along the Ising x direction, sorted in binary order with spin |↓〉 ≡ 0 and |↑〉 ≡ 1. The net detection fidelity of each spin is ~97%, after postfiltering the measurements based on calibrating the known detection errors for each spin (18, 20). The initial paramagnetic phase shown in red (B/J0 ≈ 5) exhibits a roughly uniform probability of 1/1024 ≈ 0.1% for each state (the residual peaks in the red trace are consistent with detection errors). The spin-ordered phase shown in black (B/J0 = 0.01) displays the emergence of the two AFM states, each with an occupation probability of about 8.5%. Other prominent peaks correspond to single spin-flips and other low-lying excitations from the two ground states.

Fig. 3 Quantum phase transition from a paramagnet to an antiferromagnet (Eq. 1), with Jij ≈ |i j|−1.05 in a system of 10 spins.

The exponent (α = 1.05) is estimated from the average couplings between spins in this inhomogeneous system. (A) Image of 10 trapped ions, with a distance of 22 μm between the first and last ion. (B) Measured two-point correlation function between a chosen spin (on the edge) and the other spins separated by r lattice sites, Embedded Image, averaged over 4000 experiments for each value of the parameter B/J0. For B/J0 = 5, the spins are initially polarized along the transverse y field with little correlation along the Ising x direction. As the field is reduced, the spins cross over to predominantly AFM states |↑↓↑↓ …〉 and |↓↑↓↑ …〉, resulting in alternating signs in the two-point correlations with separation. The solid lines are shown to guide the eye. (C) Measured occurrence probability of all 210 = 1024 states at B/J0 ≈ 5 (paramagnetic state, red trace) and B/J0 ≈ 0.01 (AFM phase, black trace). The states are listed in binary order, with spin |↓〉 ≡ 0 and |↑〉 ≡ 1. The residual peaks in the red trace are consistent with detection errors biased toward states with many |↑〉 spins such as 127, 511, and 1023. The two tall peaks in the black trace at 341 (0101010101) and 682 (1010101010) correspond to two Néel-ordered staggered AFM states, shown with camera images of these cases and contributing ∼17% to the population.

In Fig. 4, we probe the frustration in the system for various ranges of interactions. Here, we look at the spin order achieved in the quantum simulation when the external magnetic field is ramped down to B/J0 ≈ 0.01 for four different ranges of interactions. In Fig. 4A, we show the measured structure function S(k) at wave vector valuesEmbedded Imagefrom the measured two-point correlation functions. To directly compare the different interaction ranges, we choose the same external magnetic field ramp time constant, τ = 0.4/J0. As the range of interaction increases, the ground state AFM order (given by the structure function peak at k = π) disappears, reflecting increased occupation of the excited states as the frustration grows. Figure 4B displays the observed distribution of energy P(Ei) for α = 1.05 (shorter range) and α = 0.76 (longer range) power-law exponents, along with the cumulative energy distribution function. The eigenenergies of each configuration are calculated using Eq. 1 with B = 0. For the longer-range interactions, excitations are more prevalent, and the energy gap between the ground and the first excited state is reduced, both of which are signatures of increasing frustration in the system. The observed final entropy per particleEmbedded Image is seen to increase from 0.832 to 0.903 as the interaction range grows from α = 1.05 (ξ = 4.6 sites) to α = 0.67 (ξ = 11 sites), which is also a signature of the increased frustration in the system. As a reference, the paramagnetic state distribution shows an entropy per particle of 0.959, which is slightly less than unity because of detection errors.

Fig. 4 (A) Structure function S(k) versus wave vector k for various ranges of AFM interactions, for B/J0 = 0.01 in a system of N = 10 spins.

The increased level of frustration for the longer-range interactions reduces the observed antiferromagnetic spin order. The detection errors may be larger than shown here for the longest range of interactions, owing to spatial crosstalk from their closer spacing. (B) Distribution of observed states in the spin system, sorted according to their energy Ei that was previously calculated by diagonalizing Eq. 1 with B = 0. Data are presented for two ranges (red for α = 1.05 and blue for α = 0.76). The dashed lines indicate the cumulative energy distribution functions for these two ranges.

Quantum Coherence

The above measurements of the state distribution concern only the diagonal components or populations of the density matrix. To characterize the quantum coherence in the simulation, we retrace the external magnetic field back to its initial value after ramping it down to almost zero and measure each spin along the transverse (y) magnetic field. Figure 5 shows the distribution of the measured value my of the total transverse spin operator Embedded Imageat three different times: first at the beginning of the simulation, second when the transverse field has been ramped down to nearly zero, and third after the transverse field returns to its initial value (α = 0.92 for these data). The initial state is ideally a delta function at 〈Sy〉 = 10, but finite detection efficiency broadens the distribution. At the lowest value of the field (B/J0 ≈ 0.01), the transverse magnetization is distributed near 〈Sy〉 = 0, as the spins are presumably ordered along the Ising x direction. When the external field is ramped back to its initial state, the distribution of the total spin returns toward the initial distribution, with an average magnetization that is approximately 〈Sy〉 = 68(4)% of its initial value, indicating that some level of quantum coherence is maintained throughout the simulation. This number is in agreement with the theoretically estimated average magnetization of 〈Sy〉 ≈ 70% of its initial value, obtained by numerically integrating the Schrödinger’s equation without any decoherence.

Fig. 5 Quantum coherence in the simulation, in a system of N = 10 spins.

Probabilities of different values of the total spin component along y in the initial polarized state (red), when the transverse field is ramped to ≈0.01 J0 (black), and when the transverse field is reversed back to its initial value (green). After reversal, the y magnetization returns to ∼68(4)% of its initial value, indicating quantum coherence in the evolution. The trajectory of the transverse field (B, in green) and all the Ising couplings (J, in blue) are shown in the inset. α = 0.92 for these data.

To probe decoherence in the simulation, we repeat the experiment with various ramping speeds of the effective magnetic field. In Fig. 2, we plot the AFM order parameter Embedded Image versus the total duration for the experiment for a long-range coupling (α = 1.05) for N = 10 spins. Each data point represents the spin order achieved after ramping the magnetic field B down exponentially from 5J0 for a total duration of six time constants. The AFM order grows with slower ramping, as expected, for up to τ = 400 μs. But we also observe a saturation and then decay in the spin order, which might indicate the presence of decoherence in the system at long times. During the simulation, spontaneous Raman scattering from the optical beams is expected to occur at a rate of less than 6 s−1 per spin (21), which is consistent with separate measurements of the spin relaxation from a single spin and is therefore not expected to contribute to decoherence given the time scales in the experiment. The phonon population is expected to be well under 10% for all the data presented here (22). A principal source of decoherence appears to be the intensity fluctuations in the Raman beams, from beam pointing instabilities and fluctuations in the optical power.

Maintaining adiabaticity while ramping the magnetic field down is more difficult when the Ising interactions in the experiment are frustrated, because the relevant energy gaps are smaller. To directly compare frustrated versus nonfrustrated systems, we execute quantum simulations of both long-range AFM and ferromagnetic (FM) Ising models in a system of N = 16 spins (Fig. 6). For the FM experiment (Fig. 6C), we initialize the spins in the highest-energy state with respect to the transverse field |↓yyy …〉 and ramp the field down as before (10, 15). For the AFM experiment with the same ramp rate, we find that the best AFM nearest-neighbor correlation (Fig. 6A) is only 30(3)%, corresponding to a staggered magnetization of about 30(2)%, whereas the simulation of the FM model shows a clear FM spin order across the chain (Fig. 6B), reaching 73(3)% magnetization. We have also observed a level of 70(10)% FM magnetization emerging in N = 18 spins (18).

Fig. 6 Magnetic ordering in N=16 spins.

(A) Image of 16 trapped ions, with a distance of 30 μm between the first and last ion. (B) Pair correlation function measured at various stages of the quantum simulation, for B/J0 = 5 (red) and B/J0 = 0.01 (blue) in an AFM coupling falling off with distance as Jij ∼ |i j|−1 among N = 16 spins. Small amounts of staggered order are seen, tempered by the small gaps and frustration in the low-energy states. (C) In contrast, for all FM couplings (again with Jij ∼ |i j|−1), the gaps are large and clear FM order is seen. Here the measured distribution of magnetization is plotted. The paramagnetic phase of 16 spins is indicated in red, and after the field is ramped to nearly zero, the distribution clearly bifurcates, indicating population weighted heavily toward the FM states |↓↓↓…〉 and |↑↑↑…〉, indicated in blue. The resulting magnitude magnetization is ~73%.

The interacting spin system that we study is approaching a complexity level at which it becomes difficult or impossible to predict its behavior. Static properties such as the ground state order or the excited state energies can be calculated for hundreds of spins using Monte Carlo methods (23); however, the calculation of dynamics of fully connected spin models is currently limited to approximately N = 30 spins (24, 25). Small ion trap quantum simulators such as that reported here may soon reach this milestone with technical upgrades in the hardware, including lower vacuum chamber pressures to prevent collisions with the background gas, better stability of the optical intensities, and higher optical power so that fluctuations in the beam inhomogeneities can be suppressed.

Supplementary Materials

Supplementary Text

Fig. S1

References (2628)

References and Notes

  1. See supplementary materials on Science Online.
  2. Quantum Monte Carlo algorithms can be used to calculate static equilibrium properties of the transverse field Ising model for large numbers of interacting spins, so ground states and static correlation functions can indeed be calculated for large systems (23). However, the calculation of dynamics and nonequilibrium behavior of quantum spin models is not currently feasible for these Monte Carlo approaches, and in the presence of frustrated long-range interactions, the general behavior of such systems requires exact diagonalization. Other techniques such as the density matrix renormalization group become too difficult with long-range interactions. Because the size of the Hilbert space grows exponentially with the number of spins, sparse-matrix techniques such as the Lanczos method must therefore be used. Current state-of-the-art work on such systems is limited to sizes on the order of 30 to 35 spins (24).
  3. Acknowledgments: We thank E. Demler, L. Duan, D. Huse, K. Kim, P. Richerme, R. Sensarma, and P. Zoller for critical discussions. This work is supported by the U.S. Army Research Office (ARO) award W911NF0710576 with funds from the Defense Advanced Research Projects Agency Optical Lattice Emulator Program, ARO award W911NF0410234 with funds from the Intelligence Advanced Research Projects Activity, and the NSF Physics Frontier Center at the Joint Quantum Institute. J.K.F. was also supported by the McDevitt bequest at Georgetown.
View Abstract

Stay Connected to Science

Navigate This Article