Intrinsic phonon effects on analog quantum simulators with ultracold trapped ions

Intrinsic phonon effects on analog quantum simulators with ultracold trapped ions

C.-C. Joseph Wang Department of Physics, Georgetown University, 37th and O Sts. NW, Washington, DC 20057, USA    J. K. Freericks Department of Physics, Georgetown University, 37th and O Sts. NW, Washington, DC 20057, USA
July 12, 2019
Abstract

Linear Paul traps have been used recently to simulate the transverse field Ising model with long-range spin-spin couplings. We study the intrinsic effects of phonon creation (from the initial phonon ground state) on the spin-state probability and spin entanglement for such quantum spin simulators. While it has often been assumed that phonon effects are benign because they play no role in the pure Ising model, they can play a significant role when a transverse field is added to the model. We use a many-body factorization of the quantum time-evolution operator of the system, adiabatic perturbation theory and exact numerical integration of the Schrödinger equation in a truncated spin-phonon Hilbert space followed by a tracing out of the phonon degrees of freedom to study this problem. We find that moderate phonon creation often makes the probabilities of different spin states behave differently from the static spin Hamiltonian. In circumstances in which phonon creation is minor, the spin dynamics state probabilities converge to the static spin Hamiltonian prediction at the cost of reducing the spin entanglement. We show how phonon creation can severely impede the observation of kink transitions in frustrated spin systems when the number of ions increases. Many of our results also have implications for quantum simulation in a Penning trap.

pacs:
05.30.Rt, 03.67.Lx, 37.10.Ty, 75.10.Jm

I Introduction

Complex states of matter like spin liquids are suspected to exist in quantum spin models with frustration due to geometry or due to the nature of the spin-spin interaction spinliquid_Balents (); spinliquid_Fisher1 (); spinliquid_Varney (). Spin liquids are complicated quantum many-body states that exhibit significant entanglement of their wave functions without symmetry breaking, and could also exhibit emergent quantum phenomena within their low-energy excitation spectra. Classical computation, such as exact diagonization and quantum Monte Carlo simulation, or conventional theories based on local order parameters fail to describe these systems without bias. For example, exact diagonalization studies are limited to small size lattices and hence usually have strong finite-size effects, while quantum Monte Carlo simulations can suffer from the sign problem or have a large computational expense to describe long-range interactions and hence cannot reach the low temperatures needed to see the predicted exotic phases.

Feynman proposed that one could use controlled quantum-mechanical systems with few quantum gates to simulate many-body problems feynman (); Seth () as an useful quantum computation before achieving universal quantum computation. In recent years, there has been significant success in trying to achieve this goal by quantum simulation of desired spin models through analogous cold atom systems emulators (); emulators1 (); emulators2 (). We focus here on one platform for performing analog quantum computation, the simulation of interacting quantum spins via manipulation of hyperfine states of ions in a linear Paul trap gate (); three-ion (); Kim (); Islam (); A.Friedenauer () although many ideas presented here can be generalized to adiabatic quantum state computation in the two dimensional Penning trap as well emulators2 (). In the Paul trap systems, clock states of the ions (states with no net -component of angular momentum) are the pseudospin states, which can be manipulated independently by a pseudospin-dependent force driven by laser beams. The lasers couple the pseudospin states to the lattice vibrations of the trapped ions, which leads to effective spin-spin interactions when the phonon degrees of freedom are adiabatically eliminated spin-spin-interactions (); Duan (); Kim () based on the idea of geometric phase gate Leibfried () or Mølmer-Sørensen gate Sorensen (). Theoretically, the analog ion-trap simulators can be described as nonequilibrium driven quantum systems with both spin and phonon degrees of freedom. Sufficiently small systems can be treated numerically in an exact fashion by truncating the phonon basis and taking into account all possible quantum states in the solution of the time-dependent Schrödinger equation. Experimentally, ion traps have been used to simulate the transverse-field Ising model with a small number of ions  A.Friedenauer (); three-ion (); Islam () based on simulated quantum annealing quantum annealing () (see Ref. kim_review, for a review). It has been known experimentally that moderate phonon creation is commonplace (on the order of one phonon per mode) Islam (), even when the system is cooled to essentially the phonon ground state prior to the start of the simulation. In addition, the role phonons play are intrinsic and essential for the mediated spin-spin interaction in trapped ion systems especially in the presence of noncommuting magnetic field Hamiltonian in addition to the spin Hamiltonian of interest. Therefore, an understanding of the role phonons play in the spin simulator is crucial to understanding its accuracy.

The organization of this paper is as follows. In Sec. II, we describe the microscopic Hamiltonian for the ion-trap-based simulators and then show how one can factorize the time-evolution operator into a pure phonon term, a coupled spin-phonon term, a pure spin-spin interaction term, and a complicated term that primarily determines the degree of entanglement of the spins. Next we use adiabatic perturbation theory to determine how adiabatic state evolution can be used to reach a complicated, potentially spin-liquid-like ground state, and detail under what circumstances the evolution is not adiabatic (diabatic). In Sec. III, we show numerical comparison studies in various relevant circumstances based on a direct integration of the time-dependent Schrödinger equation, including both spin and phonon degrees of freedom (the latter in a truncated basis). In Sec. IV, we conclude with discussions and possible experimental limitations and improvements.

Ii Theory

ii.1 Microscopic Hamiltonian

When ions are placed in a linear Paul trap three-ion (); Kim (); MichaelJohanning () with harmonic trapping potentials, they form a nonuniform (Wigner) lattice, with increasing interparticle spacing as one moves from the center to the edge of the chain. The ions vibrate in all three spatial dimensions about these equilibrium positions James () with normal modes. Two hyperfine clock states (relatively insensitive to external magnetic field fluctuations because the -component of total angular momentum is zero) in each ion will be the pseudospins (and are split by an energy difference ). Hence, the bare Hamiltonian including the pseudospin and motional degrees of freedom for the ion chain is given by

 H0=∑jℏω02σzj+∑ανℏωαν(a†ανaαν+12), (1)

where is the Pauli spin matrix at the th ion site and the second term is the phonon Hamiltonian with the phonon creation operator of the normal mode along the three spatial directions . The notation refers to the pseudospin orientation in the Bloch sphere. The th spatial component of the th ion displacement operator is related to the th phonon normal mode amplitude (unit norm eigenvector of the dynamical matrix) and the th phonon creation and annihilation operator via with the mass of the ion and the normal-mode frequency.

A laser-ion interaction is imposed to create a spin-dependent force on the ions by using bichromatic laser beams to couple these clock states to a third state via stimulated Raman transitions PJLee (). Effectively, this process is equivalent to an off-resonant laser coupling to the two clock states by a small frequency detuning determined by the frequency difference of the bichromatic lasers. The ions are crystallized along the easy axis (-axis) of the trap with hard axes in the - and -directions where the transverse phonons lie. Then coupling the Raman lasers in the transverse direction minimizes effects of ion heating and allows for an identical spin axis for each ion Duan (). By accurate control of the locked phases of the blue detuned and red detuned lasers with similar Rabi frequencies, an effective laser-ion Hamiltonian PJLee (); Wineland () along the spin direction can be engineered in the Lamb-Dicke limit :

 HLI(t)=−ℏN∑j=1ΩejΔk⋅δ^Rj(t)σxjsin(μt), (2)

in which the effective Rabi frequency is generated by one effective blue-detuned beam and one red-detuned beam simultaneously (refer to Appendix A for details).

In experiments, one uses adiabatic quantum state evolution to evolve the ground state from an easily prepared state to the desired complex quantum state that will be studied. For spin models generated in an ion trap, it is easy to create a fully polarized ferromagnetic state along the -direction via optical pumping, and then apply a spin rotation (for instance, with a pulsed laser) to reorient the ferromagnetic state in any direction (usually chosen to be the -direction). Then, if one introduces a Hamiltonian with a magnetic field in the direction of the polarized state, it is in the ground state of the system. By slowly reducing the magnitude of the field and turning on the spin Hamiltonian of interest, one can adiabatically reach the ground state of the Hamiltonian (at least in principle). Hence, one has an additional Zeeman term with a spatially uniform time-dependent effective magnetic field coupled to the different Pauli spin matrices as

 HB(t)=N∑j=1B(t)⋅^σj (3)

(the magnetic field has units of energy here, since we absorbed factors of the effective magnetic moment into the definition of B, and it can also be expressed in units of frequency if we use units with ). Note that we are using an unconventional sign for the coupling to the magnetic field, since this is the sign convention often used in the ion-trapping community. In this case, the ground state of the magnetic field Hamiltonian has the spins aligned opposite to that of the field, while the highest-energy state has them aligned with the field. The magnetic field is made in the -direction by directly driving a resonant radio-frequency field with frequency between the two hyperfine states to implement the spin flips Wineland () or by indirect Raman coupling through a third state to effectively couple the two hyperfine states PJLee (). The full Hamiltonian is then .

ii.2 Factorization of the time evolution operator

We solve for the quantum dynamics of this time-dependent Hamiltonian by calculating the evolution operator as a time-ordered product and operating it on the initial quantum state . For the adiabatic evolution of the ground state, we start our system in a state with the spins aligned along the magnetic field and the system cooled down so that there are no phonons at time : . This means we will be following the highest excited spin state of the system, as described in more detail below. While it is also possible to examine incoherent effects due to thermal phonons present at the start of the simulation, we do not do that here, and instead focus solely on intrinsic phonon creation due to the applied spin-dependent force.

In time-dependent perturbation theory, one rewrites the evolution operator in the interaction picture with respect to the time-independent part of the Hamiltonian. This procedure produces an exact factorized evolution operator

 U(t,t0)=e−iℏHph(t−t0)UI(t,t0) (4)

which is the first step in our factorization procedure [the first factor is called the phonon evolution operator .] [Note that is time independent and it is multiplied by the factor in the exponent.] The second factor is the evolution operator in the interaction picture, which satisfies an equation of motion given by , with since . The only difference between and is that the phonon operators are replaced by their interaction picture values: and .

We now work to factorize the evolution operator further. Motivated by the classic problem on driven harmonic oscillators  gottfried (), we factorize the interaction picture evolution operator via with defined by (we call the factor on the left the phonon-spin evolution operator and the one on the right is the remaining evolution operator). The key step in this derivation is that the multiple commutator satisfies . This fact greatly simplifies the analysis below.

The equation of motion for the remaining evolution operator satisfies

 iℏ∂∂t¯U(t,t0)=¯H(t)¯U(t,t0), (5)

in which the operator is given by the expression

 ¯H(t)=eiℏWI(t)[−iℏ∂t+HB(t)+VI(t)]e−iℏWI(t). (6)

The operator can then be expanded order by order as

 eiℏWI(t)VI(t)e−iℏWI(t) = VI(t)+iℏ[WI(t),VI(t)] (7) eiℏWI(t)HB(t)e−iℏWI(t) = N∑j=1{B(t)⋅^σj+ \lx@stackreliℏ[WI(t),B(t)⋅^σj]+…}Residual terms, eiℏWI(t)iℏ∂te−iℏWI(t) = VI(t) (9) + 12iℏ[WI(t),VI(t)]

where we used the facts that and . Explicit calculations then yield

 VI(t) = −N∑j=1N∑ν=1∑αℏΩjηανbανj(aανe−iωανt+a†ανeiωανt) (10) × sin(μt)σxj,

in which is the Lamb-Dicke parameter for the phonon mode with the th component of the laser momentum . When the terms in Eq. (II.2) vanish, virtually excited phonons will be shown to play no role on the spin-state probabilities as a function of time, but in the presence of a transverse field, due to the noncommuting nature of quantum operators, phonon creation can significantly affect the spin-state probabilities. This fact has not been considered in detail before, and involves one of the most important results of our work, as detailed below.

ii.3 Ising spin models and Mølmer-Sørensen gate for vanishing transverse magnetic field

With a vanishing transverse magnetic field, the Hamiltonian can be greatly reduced to the spin-only Hamiltonian . Because the spin operators in and commute, one can exactly derive the following Ising spin Hamiltonian spin-spin-interactions (); Duan ()

 Hspin(t)=N∑j,j′=1Jjj′(t′)σxjσxj′. (11)

Then the expression for the spin exchange interaction is

 Jjj′(t)=ℏ2N∑ν=1ΩjΩj′η2XνbXνjbXνj′μ2−ω2Xν[ωXν−ωXνcos2μt−2μsinωXνtsinμt], (12)

which can be uniformly antiferromagnetic () or ferromagnetic () for the instantaneous ground state of the Hamiltonian when the laser detuning is detuned close to center of mass mode frequency . However, the interaction can also be inhomogeneous and frustrated when the laser is detuned in between phonon modes with details depending on the properties of the nearby phonon modes , , and .

The Mølmer-Sørensen gate Sorensen () was originally proposed to disentangle phonon effects from the spins in ion-trap quantum computing. It was discovered that because the phonons are harmonic, one could operate on the spins in such a way that the phonon state is unmodified after the gate operation (irrespective of the initial population of phonons). But one needs to keep in mind that this gate has no transverse field present, which can modify it because the transverse field operator does not commute with the Ising Hamiltonian. We begin our discussion by assuming that the laser is closely detuned to the transverse center of mass mode with angular frequency () and the addressing laser intensity for each ion is uniform and moderate (). In this situation, the time dependent term with the frequency in the interaction can be neglected. Therefore, the interaction and the operator , which are proportional to the collective spin operator , can be reduced to the following forms

 VI(t) = iℏηCMΩ2√NSx[e−i(ωCM−μ)taCM+h.c.], (13) WI(t) = ℏηCMΩ2√N(ΩCM−μ)Sx (14) × [(1−e−i(ωCM−μ)t)aCM+h.c.]

where is the Lamb-Dicke parameter for the center of mass mode and is chosen.

There are two parameter regimes where phonon effects disappear. In the weak-coupling regime , the operator almost vanishes and the time evolution operator is solely determined by the spin-only Hamiltonian because the phonon dynamics are adiabatically eliminated. Any extra phonon state redistribution takes a long time to be experimentally observable and therefore phonon effects are under control. Outside the weak-coupling regime, one can also prevent phonon effects by preparing spin states determined by the spin Hamiltonian at a particular waiting time interval . The idea is to choose the time interval such that the operator is periodic with integer cycles so that . Therefore, the initial phonon state at the start of the simulation will be revived at the time intervals as can be clearly seen from Eq. (14) when the we start with the phonon ground state for example, but it is generally true for any occupancy of the phonon states.

ii.4 Effective spin Hamiltonian with a transverse field

In the presence of a transverse magnetic field, the ideas of the Mølmer-Sørensen gate are modified. While it is tempting to claim that the residual spin-phonon terms in the magnetic field are irrelevant in the Lamb-Dicke limit , it is difficult to quantify this if the residual terms are relevant in the presence of the time dependent magnetic field which can be large in magnitude. In fact, phonon effects often modify the time evolution of the spin states when a transverse field is present. However, one can say that in cases where the integral of the field over time is small (which occurs when the field is small, or when it is rapidly ramped to zero) or when the field lies along the -direction only (or vanishes), the residual terms in Eq. (II.2) are irrelevant. For large detuning (weak-coupling regime), where or , the residual terms are always higher-order perturbations with respect to the leading transverse magnetic field term in the course of the quantum simulation. On those occasions, one can also consider the system as described by only the quantum Ising spin model in a transverse magnetic field. In general though, we need to determine how large the residual terms are, which often can only be done with numerical calculations. We illustrate this for a number of different cases below. The residual terms are

 Hres(t)=exp[iWI(t)/ℏ]HB(t)exp[−iWI(t)/ℏ]−HB(t). (15)

The equation of motion for can be written as

 iℏ∂∂t¯U(t,t0) = {i2ℏ[WI(t),VI(t)]+HB(t)+Hres(t)} (16) × ¯U(t,t0).

We perform the final factorization by writing where the spin-evolution operator satisfies and the entangled evolution operator satisfies

 iℏ∂∂tUent(t,t0)=U†spin(t,t0)Hres(t)Uspin(t,t0)Uent(t,t0). (17)

The spin evolution operator becomes

 Uspin(t,t0)=Ttexp⎡⎣−iℏ∫tt0dt′⎛⎝N∑j,j′=1Jjj′(t′)σxjσxj′+B(t′)N∑j=1σyj⎞⎠⎤⎦=Ttexp[−iℏ∫tt0dt′Hspin(t′)], (18)

which is the third factor for the evolution operator of the Ising model in a transverse field and we define in the exponent. The spin exchange terms as given in Eq. (12) include a time-independent exchange interaction between two ions and a time-dependent exchange interaction . The time-independent term can be thought of as the effective static spin-spin Hamiltonian that is being simulated, while the time-dependent terms can be thought of as diabatic corrections, which are often small in current experimental set-ups, but need not be neglected. For simplicity, we set the initial time .

The entanglement evolution operator is a complicated object in general, but it simplifies when one can approximate the operator as for the special situations discussed at the end of the last subsection. In general, this evolution operator involves a coupling of spins to phonons in all directions and has a very complicated time dependence. If one evaluates the first few terms of the series for the time-ordered product, one finds it involves multispin interactions, spin-phonon coupling, and spin-exchange interactions in all spatial directions. But the net weight of all of the terms is governed by the integral of the magnetic field over time, so if that integral is small, then this factor will also be small. Therefore, the adiabatic elimination of phonons based on Mølmer-Sørensen gate Sorensen () can be justified only in the case of a vanishing transverse magnetic field. With a constant magnetic field, the entanglement between spins and phonons can be periodic so that phonon effects can continue to be nulled at integer multiples of the appropriate period Leibfried (); Sorensen (). But such a procedure would be more complicated than the standard gate, and is not relevant for adiabatic state creation simulations, so we won’t discuss it further here. From a mathematical standpoint, because the entanglement evolution operator is on the far right of the factorization, it’s main effect is to modify the state from an initial spin state in direct product with the phonon vacuum to a state that will typically involve some degree of entanglement between phonon and spin degrees of freedom.

We can use this factorization to show that in cases where the spin-entanglement evolution operator can be approximated by the unit operator, then phonons have no observable effects on the probability of product states (regardless of the number of coherent phonons created during the simulation), so this result is similar in spirit to the original Mølmer-Sørensen gate, but is different because it holds in the presence of a transverse field and requires no special times for periodic variations to recur. To do this, we need one final identity. We further factorize the entangled phonon-spin evolution operator into the product with the spin operator defined to be and its complex conjugate is , while the function satisfies .

At this stage, we have factorized the evolution operator into four main terms, each term being an evolution operator evolving the system from time to time . We have explicit values for the first three factors, but the last term (the entanglement evolution operator) can be quite complicated; we have also described situations where the exponent of that term is small and can be neglected. In this case, the probabilities to observe any of the product states with a quantization axis along the Ising axis ( for the ionic spins) is unaffected by the presence of an arbitrary number of real excited phonons (which are excited by the phonon-spin evolution operator). Using the fundamental axiom of quantum mechanics, the probability to observe a product spin state starting initially from the phonon ground state and not measuring any of the final phonon states involves the trace over all possible final phonon configurations

 Pβ(t)=∞∑nX1=0...∞∑nXν=0...∞∑nXN=0|⟨β|⊗⟨nX1..nXν..nXN|U(t,t0)|0⟩⊗|Φ(0)⟩|2 (19)

where is any initial spin state (it need not be a product state) and denotes the number of phonons excited in the th mode in the -direction. The operator in the matrix element entangles the phonons and the spins, so we evaluate the matrix element in two steps: (1) first we evaluate the phonon part of the operator expectation value, and then (2) we evaluate the spin part. Note that since the pure phonon factor of the phonon evolution operator is a phase factor, it has no effect on the probabilities when evaluated in the phonon number operator basis, so we can drop that factor. Next, the term gives 1 when operating on the phonon vacuum to the right, so it can be dropped. We are thus left with three factors in the evolution operator. One involves exponentials of the phonon creation operator multiplied by spin operators (and is essentially a coherent-state excitation for the phonons with the average phonon excitation number determined by the spin state being measured), one involves products of spin operators that resulted from the factorization of the coupled phonon-spin evolution operator factor, and one is the pure spin evolution factor . The two remaining factors that appear on the left involve only spin operators, and hence the product state basis is an eigenbasis for those operators. This fact allows us to directly evaluate the expression in Eq. (19). We expand the evolution of the initial state at time in terms of the product-state basis with denoting each of the product state basis vectors and is a (complex) number. Using the fact that the product states satisfy the eigenvalue equation with eigenvalues (for ) or (for ), we arrive at the expression for the probability

 Pβ(t)=|cβ(t)|2exp⎡⎣−∑νjj′γνXj(t)γ∗νXj′(t)mxjβmxj′β⎤⎦
 ×∞∑nX1...nXN=0N∏ν=1⎡⎢⎣1nXν!⎧⎨⎩∑jj′γνXj(t)γ∗νXj′(t)mxjβmxj′β⎫⎬⎭nXν⎤⎥⎦. (20)

We used the matrix element

 ⟨nX1..nXν..nXN|e−i∑νΓ∗Xν(t)a†Xν|0⟩
 =(−i)nX1+...nXν+..nXN√nX1!..nXν!..nXN!Γ∗X1(t)nX1...Γ∗Xν(t)nXν...Γ∗XN(t)nXN (21)

in the derivation. The summations become exponentials, which exactly cancel the remaining exponential term and finally yield , which is what we would have found if we evaluated the evolution of the spins using just the spin evolution operator and ignoring the phonons altogether. Hence, the coherently excited phonons have no observable effects on the probability of product states for the transverse-field quantum Ising model when we can neglect the entanglement evolution operator. If we do not measure the probability of product states, then the terms from the coupled spin-phonon evolution operator remain spin operators, and one can show that the probabilities are changed by the phonons. In other words, it is because the spin-phonon evolution operator is diagonal in the product space basis for phonons and spins that allows us to disentangle the phonon and spin dynamics. In cases where this cannot be done, we expect the phonon and spin dynamics to remain entangled. In other words, phonons have observable effects on any spin measurements which introduces spin operators away from the Ising quantization axis such as most entanglement witness operator measurements. Finally, we may ask what does the entanglement evolution operator do to the system? It is difficult to find any simple analytic estimates of the effect of this term, but it acts on the initial state which has the spins aligned along or opposite to the magnetic field and has no phonons. During the evolution of that operator, new terms will be created which involve entanglement of spin states with states that have created phonons. If the amplitude of those extra terms is small, they will not have a large effect, but if it is not, then one has no other recourse but to examine the full problem numerically, which is what we do next. First we examine a perturbation-theory treatment, and then we consider the full numerical evolution of the system.

ii.5 Diabatic effects from time-dependent spin-spin interaction

One may have noticed that the spin evolution operator was not the evolution of a static Ising spin model. There were additional time-dependent factors in the evolution operator which arose from the additional time dependence of the exchange operators that was inherited by the phonons when they were “adiabatically” removed from the problem. In this section, we use adiabatic perturbation theory (reviewed in Appendix B) to analyze the effect of those extra time-dependent terms on the spin evolution of the system. In an adiabatic quantum simulation, one initially prepares the system in a certain pure state of the initial Hamiltonian with the occupation and the probability amplitudes in all other states vanishing []. Thereafter, the probability amplitudes to be excited into the other states can be approximated by

 αm(t)≈∫tt0dt′⟨m(t′)|∂t′H(t′)|n(t′)⟩Em(t′)−En(t′)ei[θm(t′)−θn(t′)], (22)

for later times, as long as the transition amplitudes are much smaller than one during the time evolution. This is the main expression we will use to evaluate the diabatic effects due to the time-dependent exchange interactions . Here are the instantaneous eigenstates of the spin Hamiltonian with a static exchange interaction and are the corresponding dynamic phases given by the integrals , and we assume there are no degeneracies in the instantaneous spectrum as a function of time.

Let us briefly describe the experimental protocol for a typical trapped ion quantum simulator restricted to the spin-only Hamiltonian defined in Eq. (18). The system is initially prepared in a spin-polarized state along (or opposite to) the direction of the transverse magnetic field by optical pumping followed by a spin rotation. The spin-only Hamiltonian is then turned on with a maximum effective transverse magnetic field followed by an exponential ramping down of the magnetic field to a final value at time ( is the exponential ramping time constant for the decay of the magnetic field). After evolving to time , the projection of the spin states along the -axis of the Ising Hamiltonian is taken to find the probability to be in a particular spin state at time (in actual experiments another pulse is applied to rotate the -axis to the -axis where the measurement is made).

If the system is perfectly adiabatic during the evolution, the outcome of the quantum state would be the highest excited state of the Ising Hamiltonian if the simulation starts out in the highest excited state of the magnetic field Hamiltonian at time , which corresponds to the spins aligned along the -axis. This procedure is theoretically identical to the ground state passage of the spin polarized state to the ground state of the negative of the Ising Hamiltonian with the system Hamiltonian being modified as  Islam (). In a typical trapped ion quantum simulator, the frequency is sufficiently far from any phonon frequencies such that the condition holds to avoid the heating of the system away from the initial phonon vacuum state during the simulation. In addition, the maximum magnetic field strength is much larger than the time independent exchange interactions to ensure the system initially starts in an eigenstate of the initial Hamiltonian. To optimize the adiabaticity of the simulation, the ramping time constant for the magnetic field has to be chosen to be much greater than the largest characteristic time scale of the system, which is shown below to be the minimum of the inverse of the frequencies .

We now discuss the effects of the time-dependent exchange interaction . For concreteness, we will follow the highest energy state, starting from the spin state aligned along the direction of the magnetic field. Starting with the expression for the transition probability amplitude in Eq. (22), we find the dominant diabatic transition is to the state with the minimum energy difference with the initial spin polarized state , assuming the matrix element in the numerator does not depend too strongly on , which is true when . At the initial time , where the Ising couplings can be shown to always vanish, all of the spin states with one spin flipped along the -axis are degenerate. This degeneracy will be broken by the Ising Hamiltonian at finite time . Due to spin-spin interaction in , the states along the y axis of the Bloch sphere called , have nonzero matrix components with the lowest energy gap with respect to the initial spin state .

To approximately evaluate the transition amplitude from the initial state to the two spin-flipped states , we do not actually need to know the state . The only relevant information we need is that it is one of the two spin-flipped states which tells us what the denominator is. Hence, we can approximate

 ⟨m(t′)|∂t′HIs(t′)|n(t′)⟩Em(t′)−En(t′)≈−⟨m|∂t′HIs(t′)|n⟩2B(t). (23)

Using similar reasoning, we approximate as

 Δθ(t′)≈−∫t′0dt2Bmℏe−tτ=−2ωBτ(1−e−t′τ), (24)

in which is the magnetic angular frequency. The operator consists of modes with frequencies , , and with the time derivative given by

 ∂t′Jj,j′(t′)=ℏN∑ν=1ΩjΩj′η2XνbXνjbXνj′μ2−ω2XνμωXν
 ×[sin2μt′−cosωXνt′sinμt′−μωXνsinωXνt′cosμt′].
 ≈ℏ2N∑ν=1ΩjΩj′η2XνbXνjbXνj′μ2−ω2XνμωXν(1−μωXν)sin(ωXν−μ)t′. (25)

The last approximate expression is derived by keeping the contribution from the slow mode and dropping the high frequency modes , and because the detuning is closely detuned to certain phonon modes in the quantum simulation.

As a consequence, the probability amplitude induced by a single phonon mode is given by the expression

 αIsm,ν(t)=−i8∑jj′[μΩjΩj′η2XνbXνjbXνj′ωB(μ+ωXν)]⟨m|σxjσxj′|n⟩f(t). (26)

The function can be approximated in experiments [when is much larger than at slow ramping ] as

 f(t)≈2i(ωXν−μ)[eiΔθ(t)+tτcos(ωXν−μ)t′∣∣t0]. (27)

We therefore reach the conclusion that the probability amplitude is given by

 αIsm,ν(t)=−14∑j,j′⟨m|σxjσxj′|n⟩μΩjΩj′η2XνbXνjbXνj′ωB(ω2Xν−μ2)eiΔθ(t)+tτ
 ×[1−cos(ωXν−μ)t]. (28)

We note that diabatic effects manifested in due to time-dependent Ising couplings grow exponentially in time as signifying that the theory is only accurate for short times. To suppress the diabatic effects, the criterion that has to hold for all phonon modes is

 ∣∣ ∣∣14∑ν∑j,j′⟨m|σxjσxj′|n⟩μΩjΩj′η2XνbXνjbXνj′ωB(ω2Xν−μ2)etτ∣∣ ∣∣≪1. (29)

Based on this expression, when the laser is closely detuned to one of the phonon resonance frequencies , the transition probability between states caused by becomes large (diabatic). In addition, a stronger magnetic field is required to suppress the diabatic transitions with smaller detuning . This is supported by the following numerical discussion in section III C. Notice that the above expression should be a reasonable estimation as long as the condition applies at time after the beginning of the quantum simulation. We can estimate the maximal time for which holds. The cutoff is set by where is the absolute value for the maximum exchange interaction between the spins. As a result, the cut-off time is proportional to the ramping time constant with a logarithmic factor given approximately by

 tc=τln[BmJmax]. (30)

In the parameter regime where , in which our theory holds, can be extended somewhat beyond the ramping time constant . Based on our numerical discussion, the diabatic effects are largest when the magnetic field is ramped through the transition from paramagnetic state to other targeted spin states, which is also accompanied by larger phonon creations due to the shrinkage of the spin gap near the transition.

Iii Numerical results

In this section, we focus on showing the circumstances where quantum emulators can or cannot be described by the transverse-field Ising model with high fidelity. Since our goal is to understand under what circumstances the effect of the phonons is small, we consider different cases for the time-evolution of the system including various detunings and initial transverse magnetic field strengths.

To isolate different effects, we compare two spin-only models in the presence of the ramping magnetic fields with the theoretically exact spin-phonon model based on numerical diagonization. The first is the ideal spin model which considers the evolution of the system with a static Ising model (spin-exchange coefficients are the time-averaged exchange coefficients) and a time-dependent magnetic field. While one might think this is a purely adiabatic model, it has some diabatic effects, since the fully polarized state is not generically the ground state of the Ising plus magnetic field Hamiltonian because the (static) Ising exchange interactions are nonzero at the initial time. Hence, one can invoke a sudden approximation to the system initially, and find that the initial state is a superposition of different energy eigenstates. In addition, the magnetic field varies in time and hence can cause additional diabatic effects due to its derivative with respect to time.

The second is the effective spin model which involves, essentially, evolution of the spin system according to the spin evolution operator only in Eq (18). Hence it has the static Ising Hamiltonian, the time-dependent Ising interactions and the time-varying transverse magnetic field. This model can have its Schrödinger equation solved in a spin-basis only, since all phonon effects are neglected except virtual phonon excitations.

The third model is the exact spin-phonon model, where we evolve the system according to the original spin-phonon Hamiltonian expanded in the Lamb-Dicke limit [Eq. (2)]. The only approximation used in this last model is the cutoff for the phonons. The strategy we use is to numerically integrate the Schrödinger equation using a direct product basis which involves a spin state in direct product with a phonon state. We do this because the Hamiltonian only connects states that differ by plus or minus one phonon number, and hence is block sparse in this basis. The spin states are chosen to include all possible Ising spin states for the number of ions in the trap. The phonon basis is chosen to have a cutoff of a maximal phonon excitation. The maximal cutoff is always chosen to be larger than the average occupancy of the phonons in each normal mode of the ion chain. Of course, we expect more phonons to be excited into the phonon modes closest to the beatnote frequency of the lasers, so the cutoffs that are chosen will vary from one normal mode to another. For example, we often find we can set the phonon cutoff to be one for some of the phonon modes far from the driving frequency of the spin-dependent force.

To facilitate our discussion, we define the root-mean square average of the fully connected Ising interaction for ions as

 Jrms= ⎷∑j,j′|Jj≠j′|2N(N−1), (31)

in which the static Ising interaction is given by the static term in Eq. (12) and the integer indices both range from to .

iii.1 Symmetry and experimental protocols

Let us discuss the symmetry of the spin-only system first, which is relevant for the exact diagonization of the spin-only Hamiltonian. There is one spatial inversion symmetry () in the ion chain, since the equilibrium ion positions are distributed symmetrically about the origin in the trap and all phonon modes involve symmetric or antisymmetric displacements of corresponding ion positions. There is also a spin reflection symmetry (, , and ) in the spin-only models with a transverse magnetic field . This spin-reflection symmetry preserves all commutation relations of the spin operators and leaves the Hamiltonian invariant. Therefore, there are four symmetry sectors for the eigenstates of the spin model (even space, even spin; even space, odd spin; odd space, even spin; and odd space, odd spin).

If the static Ising couplings are all negative (positive), the spin ground state is ferromagnetic (antiferromagnetic) and the highest spin eigenstate is the opposite, namely antiferromagnetic (ferromagnetic). We will focus on a detuning to the blue of the center-of-mass phonon mode. In this case, all spin exchange couplings are positive and the ground state is antiferromagnetic, while the highest excited state is ferromagnetic. We will examine the adiabatic state evolution of the highest eigenstate. With all the respected discrete symmetries, we can construct the symmetric and antisymmetric ferromagnetic states of the spin-only Hamiltonian as or which is in the (even, even) or (even, odd) sectors, respectively.

The experimental protocol is to prepare the system initially in a spin polarized state , [which is the highest eigenstate of the transverse magnetic field ], by optical pumping and a coherent spin rotation and then to gradually turn off the magnetic field with an exponential ramp while keeping the spin-dependent laser force in the -direction on during simulation time through stimulated Raman transitions between the spin states. According to adiabatic evolution, if the quantum state is initially prepared in the highest eigenstate of the field-only Hamiltonian , the outcome of the quantum simulation will adiabatically follow the corresponding highest eigenstate of the Hamiltonian (Ising spin Hamiltonian), if there are no level crossings (which does not occur in this system). In the case with positive static Ising coupling, the ferromagnetic highest energy eigenstate is the symmetrical ferromagnetic entangled state (the so-called GHZ state) , when .

There are two intrinsic errors which can impede the quantum simulation in trapped ions. The first is diabatic effects which occur primarily when either parameters in the Hamiltonian are changed too rapidly in time, or when energy gaps in the instantaneous eigenvalue spectrum become to small. The second is the error induced by phonons in the presence of time-dependent transverse magnetic fields. For example, the phonon-spin Hamiltonian does not have spin-reflection symmetry because it is linear in the operators, and hence the spin-phonon interaction breaks this symmetry. One consequence of this is to couple the symmetric and antisymmetric ferromagnetic states which is likely to reduce the spin entanglement of the GHZ state. (Other errors such as phonon decoherence effects due to spontaneous emission are not considered here.)

iii.2 Probability and GHZ state entanglement measurements

Current experiments use atomic cycling transitions to measure the spin state of the ion (which clock state the ion is in), and do not measure the phonons excited in the system. Hence, the experimental observables are the probability of a spin-polarized state after tracing out phonons in the tensor product of the spin-phonon Hilbert space as mentioned above when discussing Eq. (19). If one performs rotations about the Bloch sphere prior to making the measurement of the probabilities, then one can also measure a number of different spin-entanglement witness operators.

A spin-entanglement witness operator (for a target entanglement state) Witness () is a mathematically constructed observable that has a negative expectation value when the system is entangled. No witness operator can measure general entanglement, but instead a witness operator is constructed to measure a specific type of spin entanglement. For example, the witness operator for an -ion chain can be constructed as Witness ()

 WGHZ=3I−2[SGHZN1+I2+N∏k=2SGHZNk+I2] (32)

with the stabilizing spin operators expressed in terms of the Pauli spin operators by

 SGHZN1=N∏k=1σyk,SGHZNk=σxk−1σxk. (33)

Notice that the target spin polarized state in this paper is along the Ising axis in the Bloch sphere instead of the axis. Therefore, we modified the original expression Witness () to our problem by the transformation and . Based on the above construction, GHZ state entanglement measurements can be detected by the observable with the density matrix constructed by pure states or mixed states during the quantum simulation. For a perfect GHZ state entanglement, one can show that the entanglement witness operator satisfies (refer to Appendix C). Any deviation from perfect GHZ entanglement would lead to a value greater than . Note that this is one of the few cases of a witness operator where the degree of entanglement is correlated with the magnitude of the expectation value of the witness operator.

iii.3 Transition to the ferromagnetic state when detuned blue of the center-of-mass phonon mode

The systems we consider range from to which is far from the thermodynamic limit. The quantum phase transition (QPT) due to the discontinuity of the ground-state wave function in the thermodynamic limit only manifests itself as a state avoiding crossing in the energy spectrum, which is adiabatically connected to the QPT at large . The system parameters and the higher set of transverse phonon modes, which belongs to the higher branch of two transverse motional degrees of freedom, for different numbers of ions are summarized in Table I. The trapping parameters are given by the aspect ratio and the CM mode frequency along the transverse (tight) axis. The axial (easy) trapping frequency is given by the product of the aspect ratio and the CM mode frequency . The choice of these parameters comes from trap parameters and typical operating regimes of the ion-trap experiment at the University of Maryland. Most results are robust with moderate changes of parameters and our choices do not intimate that fine tuning of parameters is needed to achieve the results we show.