Onset of a quantum phase transition with a trapped ion quantum simulator

Abstract

A quantum simulator is a well-controlled quantum system that can follow the evolution of a prescribed model whose behaviour may be difficult to determine. A good example is the simulation of a set of interacting spins, where phase transitions between various spin orders can underlie poorly understood concepts such as spin liquids. Here we simulate the emergence of magnetism by implementing a fully connected non-uniform ferromagnetic quantum Ising model using up to 9 trapped ¹⁷¹ Yb⁺ ions. By increasing the Ising coupling strengths compared with the transverse field, the crossover from paramagnetism to ferromagnetic order sharpens as the system is scaled up, prefacing the expected quantum phase transition in the thermodynamic limit. We measure scalable order parameters appropriate for large systems, such as various moments of the magnetization. As the results are theoretically tractable, this work provides a critical benchmark for the simulation of intractable arbitrary fully connected Ising models in larger systems.

ARTICLE

Received 22 Mar 2011 | Accepted 2 Jun 2011 | Published 5 Jul 2011 DOI: 10.1038/ncomms1374

R. Islam1, E.E. Edwards1, K. Kim1, S. Korenblit1, C. Noh2, H. Carmichael2, G.-D. Lin3, L.-M. Duan3, C.-C. Joseph Wang4, J.K. Freericks4 & C. Monroe1

A quantum simulator is a well-controlled quantum system that can follow the evolution of a prescribed model whose behaviour may be difcult to determine. A good example is the simulation of a set of interacting spins, where phase transitions between various spin orders can underlie poorly understood concepts such as spin liquids. Here we simulate the emergence of magnetism by implementing a fully connected non-uniform ferromagnetic quantum Ising model using up to 9 trapped 171Yb + ions. By increasing the Ising coupling strengths compared with the transverse eld, the crossover from paramagnetism to ferromagnetic order sharpens as the system is scaled up, prefacing the expected quantum phase transition in the thermodynamic limit. We measure scalable order parameters appropriate for large systems, such as various moments of the magnetization. As the results are theoretically tractable, this work provides a critical benchmark for the simulation of intractable arbitrary fully connected Ising models in larger systems.

Onset of a quantum phase transition with a trapped ion quantum simulator

1 Joint Quantum Institute, University of Maryland Department of Physics and National Institute of Standards and Technology, College Park, Maryland 20742, USA. 2 Department of Physics, University of Auckland, Private Bag 92019, Auckland, New Zealand. 3 MCTP and Department of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA. 4 Department of Physics, Georgetown University, Washington DC 20057, USA. Correspondence and requests for materials should be addressed to R.I. (email: [email protected]).

NATURE COMMUNICATIONS | 2:377 | DOI: 10.1038/ncomms1374 | www.nature.com/naturecommunications

2011 Macmillan Publishers Limited. All rights reserved.

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1374

Cold atomic systems provide an ideal standard for quantum simulation1,2, by virtue of their ability to support many classes of interactions as well as their excellent quantum

coherence and readout properties. Spin chains with nearest-neigh-bour interactions have been simulated in neutral atoms stored in an optical lattice3,4, whereas long-range Ising models have been implemented with small numbers of trapped atomic ions510. The

quantum coherence in such atomic systems should allow the observation of quantum phase transitions (QPTs)11 that are driven by non-thermal parameters, like the transverse magnetic eld in the long-range quantum Ising model.

In the experiment reported here, we nd the ground state of the transverse eld Ising Hamiltonian for N interacting spin-1/2 systems:

H N J B

i j

= <

1

, s s s

where i is the Pauli matrix for the ith spin ( = x,y,z and i = 1,2, , N), Ji,j>0 is the ferromagnetic (FM) Ising coupling matrix, with

J = Ji,j and B is an external eective magnetic eld. Here we set Plancks constant h to unity, and use x,y,z for the Bloch sphere coordinates and X,Y,Z for the spatial coordinates throughout the paper. Our experiment is performed according to adiabatic quantum simulation protocol12, where the dimensionless coupling B/|J| is tuned slowly enough so that the system follows instantaneous eigenstates of the changing Hamiltonian810. As B/|J|, the ground state has all spins polarized along the magnetic eld, or is paramagnetic, along the Ising direction-x. In the other limit B/|J| = 0, the spins order according to the Ising couplings and the ground state is a superposition of FM states | and | where | and

| are eigenstates of x.

(1)(1)

i

i j x

i xj

y

We characterize the magnetic order in the system by measuring various correlation functions between all N spins, including the probability of FM occupation and the second and fourth moments of the total magnetization. We compare the results with theory, which itself may become intractable for non-uniform Ising couplings as the number of spins grows beyond 2030 (ref. 13), and even NP complete for a fully connected frustrated Ising model14. This

experiment is thus an important benchmark for large-scale quantum simulation.

ResultsEngineering the quantum Ising Hamiltonian. We represent each spin-1/2 system by the hyperne clock states 2S1/2|F = 0,mF = 0

and |F = 1,mF = 0 of 171Yb + separated by vHF = 12.642819 GHz (in a real magnetic eld of ~4 Gauss dening the quantization axis), which are denoted by the eigenstates z and z of z, respectively.

These states are detected by standard spin-dependent resonant uorescence on the cycling 2S1/2 to 2P1/2 transition at 369.5 nm using a photomultiplier tube15. The ions are trapped along the Z axis of a three-layer linear Paul trap (Fig. 1a) with centre of mass (CM) vibrational frequencies of vX = 4.748, vY = 4.300, and vZ = 1.002 MHz along the X,Y (transverse) and Z (axial) directions, respectively16.

The modes of motion along X are cooled to near their vibrational ground states and within the LambDicke regime.

O-resonant laser beams address the ions globally, driving stimulated Raman transitions between the spin states and also imparting spin-dependent forces exclusively in the X-direction, as depicted in Figure 1b9,10,17. The Raman beams contain a carrier beat-note

at frequency vHF , which provides an eective uniform transverse magnetic eld B. Raman beatnotes at frequencies vHF are near

motional sidebands and generate a spin-dependent force. The radio frequency phase dierence between the carrier beatnote and the

i

a

b

Z

X

400 m

Z

X

250 m

200 m

DC

RF

DC

c

/2.1 PMT

Spins initialized along y

Spins measured along x

Raman 1 Raman 2

Figure 1 | Experimental schematics. (a) Schematic of the three-layer linear radio frequency (Paul) trap, with the top and bottom layers carrying static potentials and the middle one carrying radio frequency. (b) Two Raman beams globally address the 171Yb + ion chain, with their wave-vector difference (

k) along the transverse (X) direction of motion, generating the Ising couplings through a spin-dependent force. The same beams generate an effective transverse magnetic eld by driving resonant hyperne transitions. A CCD image showing a string of nine ions (not in present experimental condition)

is superimposed. A photomultiplier tube (PMT) is used to detect spin-dependent uorescence from the ion crystal. (c) Outline of quantum simulation protocol. The spins are initially prepared in the ground state of Biiy, then the Hamiltonian 1 is turned on with starting eld B0|J| followed by an adiabatic exponential ramping to the nal value B, keeping the Ising couplings xed. Finally the x-component of the spins are detected.

NATURE COMMUNICATIONS | 2:377 | DOI: 10.1038/ncomms1374 | www.nature.com/naturecommunications

2011 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1374

ARTICLE

Carrier 9 8 7 6 5 4 3 107 KHz

30 KHz

(N=9)

21

Ji,j / N

(KHz)

Fluoresscence from

ions (arb. unit)

1.0

0.5

0.0

63 KHz

(N =2)

1

9

2

3

8 7 6 6

5 5

4

4

i

1 2 3 4 5

Raman beatnote from carrier (MHz)

0

3

7

8

j

2

1

9

Figure 2 | Motional modes and Ising couplings. Transverse26 vibrational modes are used in the experiment to generate Ising couplings according to equation (2). (a) Raman sideband spectrum of vibrational normal modes along transverse X-direction for nine ions, labelled by their index m. The two highest frequency modes at v1 (CM mode) and n n n

2 12 2

= z (tilt mode) occur at the same position independent of the number of ions. The dotted and the dashed lines show beatnote detunings of v1 + 30 KHz and v1 + 63 KHz used in the experiment for N = 9 and N = 2 ions respectively. Carrier transition, weak excitation of transverse -Y and axial -Z normal modes and higher order modes are faded (light grey) for clarity. (b) Theoretical

Ising-coupling pattern (equation (2)) for N = 9 ions and uniform Raman beams. The main contribution follows from the uniform CM mode, with inhomogeneities given by excitation through the other nearby modes (particularly the tilt mode). Here J1,1 + r~1/r0.35 (r1), as found out empirically. For

larger detunings, the range of the interaction falls off even faster with distance, approaching the limit Ji,j~1/|i j|3 for v15.

average beatnote of the sidebands is set to /2 so that the magnetic eld is transverse to the Ising couplings (equation (1))9,10. We sup

press direct sideband (phonon) excitation by ensuring that the beatnote detuning is sufficiently far from any normal mode frequency17. This requires that | vm|i,mi, where i,m is the Lamb

Dicke parameter of the ith ion and mth normal mode at frequency vm (with v1 = vX denoting the CM mode), and i = g2i/ is the carrier

Rabi frequency on the ith ion. Here gi is the single photon Rabi frequency of ith ion and is the detuning of the Raman beams from the 2S1/2 2P1/2 transition.

This results in an Ising interaction between the spins with control parameter that dictates the form of the coupling matrix17.

J N

i j i j

m

N i m j m m

m

80 s to a nal value B, keeping the Ising couplings xed. We nally measure the spins along the Ising (x) direction by coherently rotating the spins through /2 about the Bloch y-axis before uorescence detection. We repeat the experiment ~1,000 N times for a system of N spins and generate a histogram of uorescence counts and t to a weighted sum of basis functions to obtain the probability distribution P(s) of the number of spins in state (|), where s = 0,1, , N, as described in the Methods section.

Extraction of order parameters from measured probabilities. We can generate several magnetic order parameters of interest from the distribution P(s), showing transitions between dierent spin orders. One order parameter is the average absolute magnetization (per site)

along the Ising direction, m N N s P s

h h n

m n

In the experiment ~2.7 THz, i~370 kHz and we expect Ji,j/1 kHz for the beatnote detuning such that v14i,1i, as shown in Figure 2a. This beatnote corresponds to 63 KHz blue of the CM mode frequency for 2 ions and 30 KHz for 9 ions, as the Lamb Dicke parameter i,m~1/N. This maintains roughly the same level of virtual phonon excitation for any number of ions. The expected Ising coupling pattern for a uniformly illuminated ion chain is shown in Figure 2b for N = 9 ions and the couplings are dominated by uniform contribution of the CM mode. The non-uniformity in the Ising couplings arises from other vibrational modes and variation in i across the ion chain for gaussian Raman beams with a waist of ~70 m along the ion chain and ~6 m perpendicular to the ion chain used in the experiment. For N = 9 ions the chain is ~14 m long, and the variation in i is ~2%.

Experimental protocol. In the experiment, we follow the highest excited state of the Hamiltonian-H (refs 8,17), which is formally equivalent to the ground state of Hamiltonian H (equation (1)). It proceeds as follows (Fig. 1c). We cool all the X-transverse modes of vibration to near their ground states, and deep within the Lamb Dicke regime by standard Doppler and Raman sideband cooling procedures. We initialize the spins to be aligned to the y-direction of the Bloch sphere by optically pumping to |zzz and then

coherently rotating the spins through /2 about the Bloch x-axis with a carrier Raman transition. Next, we switch on the Hamiltonian H with an eective magnetic eld B0~5|J| so that the spins are prepared predominantly in the ground state. Then, we exponentially ramp down the eective magnetic eld with a time constant of

(2)(2)

W W

1

, , .

, =

=

N

=

=

2 2

x s

1 2

0

| | ( ) . Hamiltonian 1

has a global time-reversal symmetry of {ix ix, iz iz,

iyiy}, and this does not spontaneously break for a nite system, necessitating the use of average absolute value of the magnetization per site along the Ising direction as the relevant order parameter. For a large system, this parameter shows a second-order phase transition, or a discontinuity in its derivative with respect to B/|J|. On the other hand, the fourth-order moment of the magnetization or

Binder cumulant g N s P s N s P s

= sN sN

0 4 0 2 2

2 2

( ) ( )/( ( ) ( )) (refs 18,19)

becomes a step function at the QPT and should therefore be more sensitive to the phase transition. We illustrate this point by plotting the exact ground-state order in the simple case of uniform Ising couplings for a moderately large system (N = 100) in Figure 3a. Here we scale the two-order parameters properly to account for trivial nite size eects, as described in the Methods section. The scaled magnetization and Binder cumulant are denoted by m

x and g respectively. In Figure 3a, we also plot the exact ground-state order parameters for N = 2 and N = 9 spins. In Figure 3b-d, we present data for these two-order parameters, as B/|J| is varied in the adiabatic quantum simulation. Figure 3b shows the scaled magnetization, m

x for N = 2 to N = 9 spins that depict the sharpening of the crossover curves from para-magnetic to ferromagnetic spin order with increasing system size. The linear time scale indicates the exponential ramping prole of the (logarithmic) B/|J| scale. Figure 3c,d compares the two extreme system sizes in the experiment, N = 2 and N = 9, and clearly shows the increased steepness for larger system size. The scaled magnetization m

x is suppressed by ~25% (Fig. 3b,c) and the scaled Binder cumulant g is suppressed by ~10% (Fig. 3d) from unity at B/|J| = 0,

= =

NATURE COMMUNICATIONS | 2:377 | DOI: 10.1038/ncomms1374 | www.nature.com/naturecommunications

2011 Macmillan Publishers Limited. All rights reserved.

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1374

1.0

1.0

Order parameters

(Perfect adiabatic theory)

0.8

Binder cumulant, g

Magnetization, mx

N =2

N=9

N =100 (Uniform Ising) N =2

N =9

N =100 (Uniform Ising)

N =2 DataN =2 Theory

N=9 Theory

N=9 Data

0.5

0.0

0.01

0.6

m x

0.4

0.2

1 6

B /|J|

0.0

0.01

0.1

1 6

B /|J|

B /|J|

9

8

7

N

6

5

1.0

x 0.5

0.0

6

4

1.0

3

m

2

0.8

0.6

g

0

1

0.4

0.2

N =2 DataN =2 Theory

N =9 Theory

N =9 Data

0.1

0.25

0.0

0.01 0.01

0.1

1 6

t (ms)

0.1

0.5

B /|J|

Figure 3 | Experimental results of adiabatic quantum simulation. (a) Theoretical values of order parameters m

x and g are plotted vs B/|J| for N = 2 and N = 9 spins with non-uniform Ising couplings as used in the experiment in the case of a perfectly adiabatic time evolution. The order parameters are calculated by directly diagonalizing Hamiltonian 1. Order parameters are also calculated for a moderately large system (N = 100) with uniform Ising couplings, to show the difference between the behaviours of m

x and g. In case of uniform Ising couplings, the effective ground-state manifold reduces to N + 1 dimensions in the total spin basis. The scaled Binder cumulant g approaches a step function near the transition point B/|J| = 1 unlike the scaled magnetization m

x, making it experimentally suitable to probe the transition point for relatively small systems. (b) Scaled magnetization, mx vs B/|J| (and simulation time) is plotted for N = 2 to N = 9 spins. As B/|J| is lowered, the spins undergo a crossover from a paramagnetic to ferromagnetic phase. The crossover curves sharpen as the system size is increased from N = 2 to N = 9, prefacing a QPT in the limit of innite system size. The oscillations in the data arise because of imperfect initial state preparation and non-adiabaticity due to nite ramping time. The (unphysical) three-dimensional background is shown to guide eyes. (c) Magnetization data for N = 2 spins (circles) is contrasted with N = 9 spins (diamonds) with representative detection error bars. The data deviate from unity at B/|J| = 0 by ~20%, predominantly due to decoherence from spontaneous emission in Raman transitions and additional dephasing from Raman beam intensity uctuation, as discussed in the text. The theoretical time evolution curves (solid line for N = 2 and dashed line for N = 9 spins) are calculated by averaging over 10,000 quantum trajectories (Methods). (d) Scaled Binder cumulant (g) data and time evolution theory curves are plotted for N = 2 and N = 9 spins. At B/|J| = 0 the data deviate by ~10% from unity, due to decoherence as mentioned before.

predominantly due to decoherence from o-resonant spontaneous emission and additional dephasing due to intensity uctuations in Raman beams, during the simulation.

We compare the data shown in Figure 3c,d with the theoretical evolution, taking into account experimental imperfections and errors discussed below, including spontaneous emission to the spin states and states outside the Hilbert space, and additional decoherence. The evolution is calculated by averaging 10,000 quantum trajectories. This takes only 1 minute on a single computing node for N = 2 spins and approximately 7 h, on a single node, for N = 9 spins. Extrapolating from this calculation suggests that averaging 10,000 trajectories for N = 15 spins would require 24 hours on a 40-node cluster, indicating the inefficiency of classical computers to simulate even a small quantum system.

Discussion

A faithful quantum simulation requires an excellent understanding of errors and their scaling, especially when the underlying problem is otherwise intractable. We characterize errors in the current simulation by plotting the observed parameter P(FM) = P(0) + P(N) for N = 2 to N = 9 spins in Figure 4. Theoretically P(FM) = 2/2N when there is no ferromagnetic order that is, P(FM) = 0.5 for N = 2 spins and exponentially goes down to 0.004 for N = 9 spins, and unity, when there is perfect ferromagnetic order. Because P(FM) involves only two of the 2N basis states, it is more sensitive to errors compared with the order parameters m

x and g. For instance, at B/|J| = 0

in Figure 3b-d and Figure 4a-d, we nd that m

x and g do not change appreciably with system size, but P(FM) degrades to ~0.55 for N = 9 spins from ~0.9 for N = 2. In Figure 4 we compare the data with theory results that include experimental sources of diabatic errors.

There are several primary sources of experimental error. Diabaticity due to nite ramping speed and error in initialization is estimated to suppress P(FM) by ~3% for N = 2 to ~8% for N = 9. This also gives rise to oscillations seen in the data (Figs 3bd, 4ad). A major source of error is the spontaneous emission from Raman beams that amounts to a ~10% spontaneous emission probability per spin in 1 ms. Spontaneous emission dephases and randomizes the spin state and loosely behaves like a spin temperature in this system, though the spins do not fully equilibrate with the bath and the total probability of spontaneous emission increases linearly during the quantum simulation. In addition, each spontaneous emission event populates other states outside of the Hilbert space of each spin with a probability of 1/3. Spontaneous emission errors grow with increasing system size, which also suppresses P(FM) order with increasing N, as seen in Figure 4a-d. We theoretically estimate the suppression of P(FM) due to diabaticity and spontaneous emission, together by averaging over quantum trajectories to be ~7% for N = 2 spins and ~26% for N = 9 spins. Intensity uctuations on the Raman beams during the simulation modulate the AC Stark shi on the spins, and dephase the spin states, which causes additional diabaticity and degrades the nal ferromagnetic order. When we introduce a theoretical dephasing rate of 0.3 per ms per ion (Methods) in

NATURE COMMUNICATIONS | 2:377 | DOI: 10.1038/ncomms1374 | www.nature.com/naturecommunications

2011 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1374

ARTICLE

1.00.80.60.40.20.0

1.00.80.60.40.20.0

1.0

P(FM)

0.80.60.40.20.0

P(FM)

N=2 Data N=3 Data

N=2 Theory N=3 Theory

N=4 Data N=5 Data

N=4 Theory N=5 Theory

N=8 Data N=9 Data

N=8 Theory N=9 Theory

1 6 0.1 0.01

B/|J|

1

Diabaticity

Spontaneous

emission

Dephasing

P(FM)

1.00.80.60.40.20.0

6 0.1 0.01

B/|J|

P(FM)

N=6 Data N=7 Data

N=6 Theory N=7 Theory

6 0.1 0.01

B/|J|

1

6 0.1 0.01

B/|J|

1

Figure 4 | Suppression of P(FM) with increasing number of spins. (ad) Ferromagnetic order P(FM) = P(0) + P(N) is plotted vs B/|J| for N = 2 to N = 9 spins. The circles are experimental data, and the lines are theoretical results including decoherence and imperfect initialization. As this quantity includes only two of 2N basis states, random spin-ips and other errors degrade it much faster than the magnetization and Binder cumulant. The representative detection error bars are shown on a few points for each N. The P(FM) reduces from ~0.9 to ~0.55 as the system size is increased from two to nine.

The principle contribution to this degradation is decoherence, predominantly due to spontaneous emission from intermediate 2P1/2 states in the Raman transition and additional dephasing, primarily due to intensity uctuations in Raman beams. Shown in d is an estimated breakdown of the suppression of P(FM) from various effects for N = 9 spins. Non-adiabaticity due to nite ramping speed, and spontaneous emission and additional dephasing due to uctuating Raman beams suppress P(FM) by ~8%, ~18% and ~24% respectively from unity (B/|J|0).

the quantum trajectory computation the predicted suppression of P(FM) increases to ~9% for N = 2 and ~50% for N = 9.

Imperfect spin detection efficiency contributes ~510% error in P(FM). Fluorescence histograms for P(0) and P(1) have a ~1% overlap (in detection time of 0.8 ms) owing to o-resonant coupling of the spin states to the 2P1/2 level. This prevents us from increasing detection beam power or photon collection time to separate the histograms. Detection error in the data include uncertainty in tting the observed uorescence histograms to determine P(s), intensity uctuations and nite width of the detection beam.

The role of phonons in the results of the quantum simulation is investigated both experimentally and numerically. The average number of centre of mass phonons excited during the simulation is numerically found to be always under 1.5 for N = 9 spins, and even lower for N < 9. We perform another set of experiments with v1 = 63 KHz for all N, which amounts to reducing the phonon excitation as the number of spins N is increased, as the LambDicke parameter i,m~1/N. We do not note any appreciable dierence (beyond the margin of experimental errors) in spin population with the results reported here. In the presence of an eective magnetic eld, phonon modes are coherently populated and generally exhibit spin-motional entanglement. However, these phonons do not alter the spin ordering and hence preserve spin-spin correlation, even if the entanglement between spin states is partly destroyed when tracing over the phonon states. This will hence degrade entanglement characterization beyond a few spins using the standard GHZ type witness operators.

In this experiment, we have qualitatively observed sharpening of crossover curves that indicates the onset of a QPT as the system size increases. This scheme can be scaled up to a larger number of spins where it is possible to quantitatively estimate nite size eects, for example, scaling in critical exponents near the phase transition point20. The primary challenges in experimenting with larger system sizes include the spontaneous emission as described above,

the requirement of larger optical power to maintain the same level of Ising couplings, and non-adiabatic eects due to shrinking gap between ground and rst excited state21,22. One solution is to implement a high power laser with a detuning far from the 2P energy levels, which would minimize spontaneous emission while maintaining the same level of Ising couplings. This would also allow versatility in varying the Ising interaction (together with the eective external eld) during the simulation, as the dierential AC stark shi between spin states is negligible for a sufficiently large detuning. The coherence time increases in the absence of spontaneous emission, allowing for a longer simulation time necessary to preserve adiabaticity as the system grows in size. Recently Raman transitions have been driven using a mode-locked high power pulsed laser at a wavelength of 355 nm, which is optimum for 171Yb + wherein the ratio of dierential AC Stark shi to Rabi frequency is minimized and spontaneous emission probabilities per Rabi cycle are < 10 5 per spin23.

With this system, it is possible to engineer dierent Ising coupling patterns by controlling the Raman beatnote detuning and observe interesting spin ordering such as with antiferromagnetic long-range couplings leading to frustration9,17, and phase transi

tions some of which can be very sharp, or of rst order24. Long-

range interactions in this spatially one-dimensional system allow for simulating multidimensional spin models by selectively exciting vibrational modes using multiple Raman beatnote detunings. With additional laser beams, this scheme can potentially simulate more complicated and higher spin-dimensional Hamiltonians like xy and xyz models, which map onto nontrivial quantum Hamiltonians such as the Bose-Hubbard Hamiltonian.

Methods

Detection of spin states. The spin states are detected by spin-dependent uorescence signals collected through f/2.1 optics by a photomultiplier tube. Spin state |z is resonantly excited by the 369.5 nm detection beam and uoresces from 2P1/2

states, emitting Poisson-distributed photons with mean ~12 in 0.8 ms. This state appears as bright to PMT. The detection light is far o-resonant to spin state |z

NATURE COMMUNICATIONS | 2:377 | DOI: 10.1038/ncomms1374 | www.nature.com/naturecommunications

2011 Macmillan Publishers Limited. All rights reserved.

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1374

and this state appears dark to the PMT. However, due to weak o-resonant excitation bright state leaks onto dark state, altering the photon distribution25. Unwanted scattered light from optics and trap electrodes also alter the photon distribution. We construct the basis function for s bright ions by convolution techniques, and include a 5% uctuation in the intensity of detection beam, which is representative of our typical experimental conditions. We then t the experimental data to these basis functions, and obtain probabilities P(s) at each time step ti in the experiment.

Mean photon counts for dark (mD) and bright (mB) states are used as tting parameters so as to minimize the error residues.

The best tting at time step ti is obtained for the parameters {mD,i,mB,i}. These parameters uctuate at dierent time steps of the quantum simulation, primarily due to uctuations in the intensity of detection beam and background scatter, and also due to uncertainties in a multivariate tting. The tting errors are propagated to the spin-state probabilities P(s) using Monte Carlo method of error analysis, as follows. We extract P(s) and compute the order parameters at time step ti with mean dark and bright state counts chosen randomly from a Gaussian distribution with mean {m

D,mB} and standard deviations {mD,mB} respectively. Here mD and ,mB are averages of mD,i and mB,i respectively over dierent time steps ti. Similarly mD and mB are standard deviations of mD,i and mB,i respectively. By repeating this process ~400 times we generate a histogram of each order parameter and t the histogramsto a Gaussian distribution. The standard deviation of the distribution is chosen to represent the random error due to tting in that order parameter. The uncertaintyin amount of uctuation of the detection beam power during the experiment is conservatively included in the error analysis by repeating the tting process for a range of uctuations. The nite width of the detection beam is taken care of by modelling the Gaussian beam having a three-step intensity prole with appropriate intensity ratios.

Scaled order parameters. To characterize the spin orders we use dierent order parameters in the experiment, namely the average absolute magnetization per site (mx) and the Binder cumulant (g). When the spins are polarized along the y-direction of the Bloch sphere, the distribution of total spin along x-direction is Binomial and approaches a Gaussian (with zero mean) in the limit of N. For system

size of N, mx takes on theoretical value of m N C N n

x N N n

in

the perfect paramagnetic phase (B/|J|) and unity in the other limit of B/|J| = 0. In Figure 3b,c, we rescale mx to m

x = (m0x,N mx)/(m0x,N 1) which should ideally be zero in perfect paramagnetic phase and unity in perfect ferromagnetic phase for any N. This accounts for the trivial nite size eect due to the dierence between Binomial and Gaussian distribution. Similarly the Binder Cumulant g is scaled to g = (g0N g)/

(g0N 1) in Figure 3d, where g0N = 3 2/N is the theoretical value of g for B/|J|.

Quantum Monte-Carlo simulations. Quantum trajectories are generated by numerically integrating the Schrdinger equation, with Hamiltonian (1), while simultaneously executing quantum jumps to account for spontaneous emission and decoherence. Spontaneous emission from ion i either localizes the spin of the ion, projecting it into 2S1/2|F = 0,mF = 0 (spin state |z) or 2S1/2|F = 1,mF = 0 (spin state |z), or projects the ion into 2S1/2|F = 1,mF = 1, in which case ion i is factored out of the Schrdinger evolution, though it is counted as spin up at the time of measurement. Decoherence (dephasing) is modelled by the quantum jump operator x; thus a jump for ion i, |

ix|, introduces a phase shi between the spin states | and | (in x-basis). Jump rates are taken to be xed and equal for all ions. Note that a decoherence jump rate of decoh leads to decay of the spin coherence at rate 2decoh. To determine the entangled state of the spin ensemble aer a spontaneous emission, for example, from ion i, we assume that the ground-state conguration before emission,

| | | | | |

, , , ,

1 N N n

0 , | |

0

= 2 2

=

+

( )+ +

z i

j i

i j z j i j z

a b g d jj

i j z j i j z j z i

j i

( ),

is mapped, by the far-detuned Raman beams, into a very small excited-state contribution to the overall system entangled state,

l a g b d

| [( )| ( )| ],

/ , , , ,

2 1 2

P i

j i

+ + +

i j i j z j i j i j z j

with 1 proportional to the amplitude of the Raman beams and inversely proportional to their detuning. The (unnormalized) state aer the emission is

|? [( )| ( )| ],

, , , ,

+ + +

i

i j i j z j i j i j z j

a g b d

where |?i is |zi, |zi, or the factored state 2S1/2|F = 1,mF = 1.

References

1. Feynman, R. Simulating physics with computers. Int. J. Theor. Phys. 21, 467488 (1982).

2. Lloyd, S. Universal quantum simulators. Science 273, 1073 (1996).3. Simon, J. et al. Quantum simulation of an antiferromagnetic spin chains in an optical lattice. Nature 472, 307312 (2011).

4. Struck, J. et al. Quantum simulation of frustrated magnetism in triangular optical lattices arXiv:1103.5944.

5. Porras, D. & Cirac, J. I. Eective quantum spin systems with trapped ions. Phys. Rev. Lett. 92: 207901 (2004).

6. Deng, X.- L., Porras, D. & Cirac, J. I. Eective spin quantum phases in systems of trapped ions. Phys. Rev. A 72 (2005).

7. Taylor, J. M. & Calarco, T. Wigner crystals of ions as quantum hard drives. Phys. Rev. A 78: 062331 (2008).

8. Friedenauer, A., Schmitz, H., Glueckert, J. T., Porras, D. & Schaetz, T. Simulating a quantum magnet with trapped ions. Nature Phys. 4, 757761 (2008).

9. Kim, K. et al. Quantum simulation of frustrated ising spins with trapped ions. Nature 465, 590593 (2010).

10. Edwards, E. E. et al. Quantum simulation and phase diagram of the transverseeld ising model with three atomic spins. Phys. Rev. B 82: 060412 (2010).

11. Sachdev, S. Quantum magnetism and criticality. Nature Phys. 4, 173185 (2008).

12. Farhi, E. et al. A quantum adiabatic evolution algorithm applied to random instances of an np-complete problem. Science 292, 472 (2001).

13. Lanczos, C. An iteration method for the solution of the eigenvalue problem of linear dierential and integral operators. J. Res. Nat. Bur. Standards 45, 255282 (1950).

14. Barahona, F. On the computational complexity of Ising spin glass models. J. Phys. A: Math. Gen. 15, 32413253 (1982).

15. Olmschenk, S. et al. Manipulation and detection of a trapped Yb+ hyperne qubit. Phys. Rev. A 76: 052314 (2007).

16. Deslauriers, L. et al. Zero-point cooling and low heating of trapped 111Cd+ ions.

Phys. Rev. A 70: 043408 (2004).17. Kim, K. et al. Entanglement and tunable spin-spin couplings between trapped ions using multiple transverse modes. Phys. Rev. Lett. 103: 120502 (2009).

18. Binder, K. Finite size scaling analysis of Ising model block distribution functions. Physik B 43, 119 (1981).

19. Binder, K. Critical properties from monte carlo coarse graining and renormalization. Phys. Rev. Lett. 47, 693696 (1981).

20. Fisher, M. E. & Barber, M. N. Scaling theory for nite-size eects in the critical region. Phys. Rev. Lett. 28, 15161519 (1972).

21. Dusuel, S. & Vidal, J. Continuous unitary transformations and nite-size scaling exponents in the lipkin-meshkov-glick model. Phys. Rev. B 71: 224420 (2005).

22. Caneva, T., Fazio, R. & Santoro, G. E. Adiabatic quantum dynamics of the lipkin-meshkov-glick model. Phys. Rev. B 78: 104426 (2008).

23. Campbell, W. C. et al. Ultrafast gates for single atomic qubits. Phys. Rev. Lett. 105: 090502 (2010).

24. Lin, G.- D., Monroe, C. & Duan, L.- M. Sharp phase transitions in a small frustrated network of trapped ion spins ArXiv 1011. 5885.

25. Acton, M. Detection and Control of Individual Trapped Ions and Neutral Atoms Ph.D. thesis (University of Michigan, 2008).

26. Zhu, S.- L., Monroe, C. & Duan, L.- M. Trapped ion quantum computation with transverse phonon modes. Phys. Rev. Lett. 97: 050505 (2006).

Acknowledgements

We thank D. Huse for help with theoretical understanding of the quantum Ising model and nite-size scaling. This work is supported under Army Research Office under Award No. 911NF0710576 with funds from the DARPA Optical Lattice Emulator Program, IARPA under ARO contract, the NSF Physics at the Information Frontier Program, the European Program on Atomic Quantum Technologies, and the NSF Physics Frontier Center at JQI.

Author contributions

R.I., E.E.E., K.K., S.K. and C.M. designed and performed the experiments, and analysed the data. C.N., H.C., G.-D.L., L.-M.D., C.-C.J.W and J.K.F provided theoretical support. C.N. and H.C. performed the quantum trajectory calculations. All the authors contributed to the manuscript.

Additional information

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/

How to cite this article: Islam, R. et al. Onset of a quantum phase transition with a trapped ion quantum simulator. Nat. Commun. 2:377 doi: 10.1038/ncomms1374 (2011).

j i

NATURE COMMUNICATIONS | 2:377 | DOI: 10.1038/ncomms1374 | www.nature.com/naturecommunications

2011 Macmillan Publishers Limited. All rights reserved.