## 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.

## Abstract

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 ^{171}Yb^{+} 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 (*1*–*3*). 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 (*10*–*14*) 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.

## Model

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

(1)where the Planck constant *h* is set to 1, are the spin-1/2 Pauli operators for the *i*th spin (*i *= 1,2,…*N*), *B* is the effective transverse magnetic field, and *J _{ij}* > 0 is the Ising coupling between spins

*i*and

*j*, falling off with the lattice spacing |

*i*−

*j*| approximately as

where 0 < α < 3 (*8*).

For *B* >> *J _{ij}* on all pairs, the spins are polarized along the effective transverse magnetic field in the ground state |↑

*↑*

_{y}*↑*

_{y}*…〉 of the Hamiltonian in Eq. 1, where |↑*

_{y}*〉 denotes a spin along the +*

_{y}*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* >> *J _{ij}*, 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: ξ ≡ 5

^{1/α}.) 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*/

*J*

_{0}and has the same spin order as the ground state near

*B*/

*J*

_{0 }= 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 ( 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.

## Experiment

The spins are realized in a collection of ^{171}Yb^{+ }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 |↓

*〉, separated by the hyperfine frequency ν*

_{z}_{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 |↑

*〉 along the effective transverse magnetic field. The Hamiltonian (Eq. 1) is then switched on with an initial field*

_{y}*B*

_{0}≈ 5

*J*

_{0}, where

*J*

_{0}(∼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 or 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 〈... 〉).

## 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 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*) *g _{s }*= 〈(

*m*−〈

_{s}*m*〉)

_{s}^{4}〉/〈(

*m*−〈

_{s}*m*〉)

_{s}^{2}〉

^{2}, scaled to to remove finite size effects, varies from to as the paramagnetic state gives way to AFM order (see Fig. 2). Here 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 , allowing a direct probe of spin order for each experimental realization. The Fourier transform of this correlation function is the structure function , where 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 *C*_{1,}* _{r}* between a chosen edge spin and the other spins are presented in Fig. 3B at various stages in the ramp

*B*/

*J*

_{0 }= 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*/

*J*

_{0}is lowered, however, a zig-zag pattern emerges, with negative (positive) correlations between spins separated by odd (even) lattice spacings. For

*B*/

*J*

_{0}≈ 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*/

*J*

_{0}≈ 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 2

^{10}= 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*/

*J*

_{0}≈ 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*/

*J*

_{0}= 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.

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*/*J*_{0} ≈ 0.01 for four different ranges of interactions. In Fig. 4A, we show the measured structure function *S*(*k*) at wave vector valuesfrom 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/*J*_{0}. 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*(*E _{i}*) 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 particle 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.

## 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 *m _{y}* of the total transverse spin operator at 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 〈

*S*〉 = 10, but finite detection efficiency broadens the distribution. At the lowest value of the field (

_{y}*B*/

*J*

_{0}≈ 0.01), the transverse magnetization is distributed near 〈

*S*〉 = 0, as the spins are presumably ordered along the Ising

_{y}*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 〈

*S*〉 = 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 〈

_{y}*S*〉 ≈ 70% of its initial value, obtained by numerically integrating the Schrödinger’s equation without any decoherence.

_{y}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 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 5*J*_{0} 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 |↓* _{y}*↓

*↓*

_{y}*…〉 and ramp the field down as before (*

_{y}*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*).

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

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵See supplementary materials on
*Science*Online. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵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*). **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.