Phonon mediated quantum spin simulator employing a planar ionic crystal in a Penning trap

# Phonon mediated quantum spin simulator employing a planar ionic crystal in a Penning trap

C.-C. Joseph Wang Department of Physics, Georgetown University, Washington, DC 20057, USA Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA    Adam C. Keith Department of Physics, Georgetown University, Washington, DC 20057, USA Department of Physics, University of Colorado at Boulder, CO 80309    J. K. Freericks Department of Physics, Georgetown University, Washington, DC 20057, USA
July 12, 2019
###### Abstract

We derive the normal modes for a rotating Coulomb ion crystal in a Penning trap, quantize the motional degrees of freedom, and illustrate how they can by driven by a spin-dependent optical dipole force to create a quantum spin simulator on a triangular lattice with hundreds of spins. The analysis for the axial modes (oscillations perpendicular to the two-dimensional crystal plane) follow a standard normal-mode analysis, while the remaining planar modes are more complicated to analyze because they have velocity-dependent forces in the rotating frame. After quantizing the normal modes into phonons, we illustrate some of the different spin-spin interactions that can be generated by entangling the motional degrees of freedom with the spin degrees of freedom via a spin-dependent optical dipole force. In addition to the well-known power-law dependence of the spin-spin interactions when driving the axial modes blue of phonon band, we notice certain parameter regimes in which the level of frustration between the spins can be engineered by driving the axial or planar phonon modes at different energies. These systems may allow for the analog simulation of quantum spin glasses with large numbers of spins.

###### pacs:
37.10.Ty , 03.67.Lx, 75.10.Jm, 75.10.Nr

## I Introduction

Richard Feynman motivated the idea of using analog quantum simulation as a means to describe the behavior of complex quantum-mechanical systems feynman (). This idea has been experimentally realized for the transverse-field Ising model in cold ion traps in one-dimension two-ion (); Kim1 (); Kim2 (); Edwards (); Islam (), in neutral atoms driven to a nonequilibrium state Markus (), and with preliminary results also available for cold ion traps in two-dimensions John (). The tunability and control of cold atom and cold ion systems are unmatched in traditional (real-material) condensed matter systems and therefore they are an ideal platform for examining idealized condensed matter systems and searching for nontrivial quantum ground states. In addition, the isolation of the system from the environment with long decoherence times (in milliseconds) also offers tremendous opportunities for the study of non-equilibrium coherent many-body dynamics in real time such as thermalization Marcos (); dynamic QS () and quench dynamics Santos (); kastner () in cold atoms.

Cold ions in traps have been proposed theoretically as potential emulators for effective spin models spin-spin-interactions (). For the case of linear ion chains in Paul traps, the experimental protocol has been well established Kim2 (); Edwards (); Islam (). However, there are still difficulties in scaling up Paul-trap based systems to more than approximately sixteen ions while maintaining good quantum control, although next generation traps are planned for up to potentially 100 ions. Theoretically, phonon modes in linear Paul traps have been well understood for a decade James () and are essential for the realization of the spin models in the system (although effects of micromotion on two-dimensional crystals in Paul traps are more complicated paultrap_2d ()). On the other hand, ion crystals in Penning traps can be easily stabilized with several hundreds of ions in a two dimensional structure and can potentially simulate quantum phases beyond what modern computers can simulate. Recently, a feasibility study has shown that it is possible to generate spin-spin interactions in these systems that decay like a tunable power law for long distances John (). Similar to cold ions in linear Paul traps, a theoretical understanding of phonon modes is needed as a prerequisite for the adiabatic state evolution of these systems.

In this paper, we provide a complete theory for the oscillatory normal modes in a Penning trap, tuned to parameters similar to those used in current experiments. We then quantize these phonons and show how they can be used with a spin-dependent optical dipole force to generate effective spin-spin interactions. Recently, others have shown how to find equilibrium positions for similar trapping potentials Chinese (), and have evaluated the phonon spectra using different methods from the ones we employ Jake ().

The organization of the paper is as follows. In Sec. II, we formulate the theory for the normal modes of cold ions in Penning traps. We determine the ground state of the crystal configuration by finding the static equilibrium of the system in the rotating frame of a rotating crystal. The normal modes are then found by harmonic expansion around the equilibrium state. Because the Penning trap creates a rotating crystal, the normal modes must be solved in the rotating frame, and the longitudinal (planar) modes require the solution of a quadratic eigenvalue problem because they have velocity-dependent forces. The phonon Hamiltonian is then derived by quantizing the normal modes. In Sec. II, we also discuss quantum simulation based on spin-dependent forces arising from laser-ion dipole interactions. The effective Ising models mediated by the axial phonon modes and the planar phonon modes are also described. In Sec. III, we show the numerical results for phonon frequencies and for spin-spin interactions for parameters relevant for current experiments. In Sec. IV, we provide our conclusions.

## Ii Theoretical formulation

We consider hundreds of spins, realized via the two hyperfine states , of ions, localized in a two-dimensional plane. The setup of cold ions in Penning traps is briefly illustrated in Fig. 1. Further details can be found in the published literature John (); John2 (). In current experiments, there also is a small amount of impurities such as which forms via collision of the Beryllium ions with hydrogen molecules that are in the background gas. Hence a real system will have different mass particles (which will tend to separate from one another as the centrifugal force pushes the heavier particles to the outer regions of the crystal). In this paper, we discuss the clean limit only, and ignore such defects. The effect of defects will be studied elsewhere.

Given the above considerations, we form the ion Lagrangian in the laboratory frame of reference as

 L=N∑j=1[12m|˙rj|2−eϕj+eAj⋅˙rj], (1)

in which is the total number of ions, is the positive unit charge of and is the mass of a ion, is the ion spatial configuration in cylindrical coordinates, is the total scalar potential acting on the th ion and is the vector potential in the symmetric gauge for the axial magnetic field . The potential , includes the harmonic trapping from the electrodes, the rotating wall potential, and the Coulomb interaction between the ions. It satisfies

 ϕj(t)=V0[z2j−12r2j]+VWr2jcos[2(θi+Ωt)]+kee22∑k≠j1rkj, (2)

in which is the amplitude of the static potential from the Penning trap electrodes, is the amplitude of the rotating wall potential, is the rotating wall angular frequency , is the Coulomb force constant, and is the inter-particle distance between the th and th ions. With the application of the quadrupole rotating wall potential , the ion potential is time dependent in the laboratory frame of reference. However, in the co-rotating frame with the angular speed , the corresponding effective potential becomes time independent and the problem can be understood as an equilibrium problem dubin_oneil (). The geometric configuration of the ions under rotation in the steady state is simplified by finding the (static) equilibrium ion configuration in the rotating frame. In general, the coordinate transformation between the rotating frame and the laboratory frame is described by the rotational operator with the angle of rotation :

 (xj(t)yj(t))=Rz(θ(t))(xRj(t)yRj(t)) (3)

where the superscript denotes the rotating frame of reference. The operator satisfies

 Rz(−Ωt)=(cos(Ωt)sin(Ωt)−sin(Ωt)cos(Ωt)), (4)

with . Accordingly, in the rotating frame, the rotating wall potential does not depend on time and is given by by direct substitution. As a consequence, the effective potential energy in the rotating frame is then given by the expression

 eϕRj=eV0zRj2+12(eBzΩ−mΩ2−eV0)rRj2+eVW(xRj2−yRj2)+ke22∑k≠j1rRjk (5)

which now includes the potential terms from the centrifugal and Lorentz forces as well.

Ions trapped in a Penning trap do not always crystallize into a two-dimensional plane. The following approximate confinement criterion has to be satisfied to guarantee it is energetically favorable for ions to stay in the plane with :

 β1=2eV0(eBzΩ−mΩ2−eV0)≫1. (6)

In addition, when the rotating wall potential vanishes (), the system is confined in the plane only if the radial potential is attractive. In this case, the confinement criterion for zero rotating wall potential is given by

 β2=eBzΩ−mΩ2−eV0>0. (7)

For nonzero rotating wall potentials, the ion saddle potential caused by the rotating wall potential breaks rotational symmetry and leads to a confinement of ions only when the coefficients of the trapping potential along the weak trapping axis is positive. In other words, the following criterion has to be satisfied to have a confined equilibrium state for all ions:

 β3=12(eBzΩ−mΩ2−eV0)−e|VW|>0. (8)

The above heuristic criteria narrow down the relevant parameter regime in which a stable planar ion crystal formation is feasible for quantum simulation. Of course, the accurate quantitative prediction of the stability of the ion crystals must also include the Coulomb interaction between the ions. The precise criterion can only be found by solving for the equilibrium positions and showing that all normal modes have real frequencies, so that the equilibrium is globally stable dubin1 ().

Based on the invariance of the Lagrangian, we derive the Lagrangian in the rotating coordinates from Eq. (1) as

 LR=N∑j=1[12m|˙rRj|2−eBeff[Ω]2(˙xRjyRj−˙yRjxRj)−eϕRj] (9)

in which the second term will produce the Lorentz force term with the effective magnetic field for the th ion in the rotating coordinates with the magnitude which depends on the rotational angular frequency . The modification of the magnetic field is due to the Coriolis force for the moving ions in the rotating frame. There is no Coriolis force present for ions when the ions are in equilibrium in the rotating frame but the Coriolis force can have effects on the oscillating normal mode motion away from equilibrium. The centrifugal force originates from the term in the effective potential energy in Eq. (5).

### ii.1 Stable configurations and normal modes in the rotating frame

It has been well established experimentally that ions in Penning traps reach different static equilibrium states in the rotating frame of reference (steady state in the laboratory frame of reference) for different values of the adjustable rotating angular frequency given by a rotating wall potential  John3 (). Therefore, we can find the stable spatial configuration of the ions by determining the local minima of the effective potential energy in the rotating frame of reference. In general, such a minimization problem is difficult to solve for the absolute minimum in two or higher dimensions, because many different configurations can be competitive for the lowest energy state and there can be large potential barriers between them. To find an experimentally viable solution, we are guided by the experiments themselves, which typically observe the ions arranged in an approximate triangular lattice in a single plane.

Our strategy is to construct the trial configurations based on a “closed shell” construction analogous to finding stable electronic configurations in an atom as detailed in the numerical discussions. Even though we cannot guarantee that the stable configuration we find is the global minimum, the NIST experimental systems seem to find the same local minimum which validates our approach. In further support of the strategy, we have successfully predicted the phonon mode spectra and spin-spin interaction observed in the NIST experiments John (); Brian ().

In the following, we discuss the collective phonon excitations about the previously determined equilibrium configuration. Based on the equilibrium configuration, normal mode dynamics of ions near the equilibrium structure can be fully captured by the full system Lagrangian expanded around this minimal ion spatial configuration to quadratic orders by a Taylor series. We use the following notation for the ion coordinates and for the ion velocities in the rotating frame. Since the equilibrium configuration is the local minimum of the classical action of the system , the first order terms in and do not contribute to the expansion by construction. Therefore, the system Lagrangian can be expanded to quadratic order as

 L=L0+12N∑j,k=1[(δRj⋅∂∂Rj)(δRk⋅∂∂Rk)+(˙δRj⋅∂∂˙Rj)(˙δRk⋅∂∂˙Rk)+2(δRj⋅∂∂Rj)(˙δRk⋅∂∂˙Rk)]L∣∣0 (10)

where the first term is the Lagrangian from the equilibrium state and the second term is the Lagrangian due to fluctuations away from the equilibrium state. One can isolate the planar motions of the ions in the plane from the axial motions along the axis based on the consideration that the Lorentz force due to the external magnetic fields acts only in the crystal plane and the potentials are separable in cylindrical coordinates, hence there is no harmonic coupling between planar and axial degrees of freedom (anharmonic effects can couple normal modes together but are ignored here). Therefore, we can study the axial and planar modes independently.

The Lagrangian which governs the axial motions of the ions is described as

 LAph=12N∑j=1mδ˙Rzjδ˙Rzj−12N∑j,k=1KzzjkδRzjδRzk, (11)

in which the symmetric stiffness matrix is a real Hermitian matrix and the matrix elements are given by

 Kzzjk=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩2eV0−N∑l=1kee2R0lj3    j=k,l≠jkee2R0jk3    j≠k,⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (12)

where is the distance between two ions in equilibrium (in the rotating frame). Similarly, we derive the Lagrangian for the planar normal modes

 LPph=12N∑j=1m|˙δRj|2−12N∑j,k=1∑αβ∈(x,y)KαβjkδRαjδRβk+N∑j=1eBeff[Ω]2[δRxjδ˙Ryj−δRyjδ˙Rxj], (13)

in which is the effective magnetic field in the rotating frame with angular frequency . The real Hermitian stiffness matrix satisfies the relation . We derive the stiffness matrix for the planar modes as:

 Kxxjk=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩mΩ2+eΩBeff[Ω]−eV0+2eVW−kee2N∑l=1R0jl2−3(x0j−x0l)2R0jl5    j=k,l≠jkee2R0jk2−3(x0j−x0k)2R0jk5    j≠k,
 Kyyjk=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩mΩ2+eΩBeff[Ω]−eV0−2eVW−kee2N∑l=1R0jl2−3(y0j−y0l)2R0jl5    j=k,l≠jkee2R0jk2−3(y0j−y0k)2R0jk5    j≠k,
 Kxyjk=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩3kee2N∑l=1(x0j−x0l)(y0j−y0l)R0jl5    j=k,l≠j−3kee2(x0j−x0k)(y0j−y0k)R0jk5    j≠k (14)

where the equilibrium configuration is represented by the ion coordinates for . Notice that the off-diagonal matrix elements have nonzero values, indicating the collective nature of the planar motional degrees of freedom, which couple the motion in the two coordinate directions together.

### ii.2 Axial phonon modes

Phonon modes are the quantized modes corresponding to the classical normal modes discussed earlier. In this section, we solve the classical normal mode problem and then find the quantized Hamiltonian of the axial phonons, which is most relevant for the description of conventional quantum simulation in a Penning trap.

By minimizing the action , the classical equation of motion for the ion displacement along the axial () axis can be written as

 mδ¨Rzj+N∑k=1KzzjkδRzk=0, j=1,2,…N. (15)

The axial normal modes are found by direct substitution into Eq. (15) of the eigenvector solution after initial time where the eigenvectors are -tuples of real numbers with unit norm. The following eigenvalue equation must therefore be satisfied

 N∑k=1[mω2zνδjk−Kzzjk]bzνk=0, j,ν=1,2,…N (16)

where the axial eigenvectors are the eigenvectors of the stiffness matrix (whose eigenvalues are ) but with different eigenvalues (since all eigenvalues are real, the frequencies are either real or pure imaginary). Note that this eigenvalue problem is the simplest quadratic eigenvalue problem, but because one can immediately solve it by this mapping onto a linear Hermitian eigenvalue problem, we need not discuss this point further here. For the planar phonons, the analysis is more complicated, as shown below.

The planar crystal system is stable against a one-to-two plane transition when all the eigenfrequencies are real, which is analogous to case of stability against the well-known zigzag transition for ions in linear Paul traps. If the eigenfrequencies are imaginary, then the two-dimensional planar equilibrium that was found is an unstable equilibrium, and will not form in experiment. A detailed numerical analysis of this instability can tell us the parameter regimes where a single sheet of ion crystals is energetically stable dubin1 ().

The procedure for the quantization of axial normal modes is the standard procedure. One identifies the generalized coordinates and momentum for each phonon mode as canonically conjugate variables, which satisfy the relation given by the Poisson bracket . The quantization of phonon modes can be implemented by promoting the relation to the commutation relations for operators. To identify the canonically conjugate variables for the axial motional degrees of freedom, the Lagrangian for the axial phonon modes can be recast into the form in Eq. (17) (as long as the system is in the parameter regime in which the planar crystal configuration is stable with Im for all modes ). We let be the displacements, where the quantum dynamics due to the phonon mode is given by the real generalized coordinate . The Lagrangian becomes

 LAph=12N∑ν=1m[˙ξ2ν−ω2zνξ2ν]. (17)

The corresponding generalized momentum for the phonon mode is then given by the expression

 PAν=∂LLph∂˙ξν=m˙ξν. (18)

As a consequence, the Hamiltonian for the axial phonon modes is represented by the summation over independent harmonic modes

 HAph=N∑ν=1⎡⎣PA2ν2m+12mω2zνξ2ν⎤⎦. (19)

Therefore, the second quantized form of the Hamiltonian for the axial phonon modes simply becomes

 ^HAph=N∑ν=1ℏωzν(^nzν+12), (20)

in which is the number operator for the phonon mode and the phonon creation and annihilation operators are related to the generalized coordinates via

 azν = √mωzν2ℏ(ξν+imωzνPAν), a†zν = √mωzν2ℏ(ξν−imωzνPAν) (21)

The phonon creation and annihilation operators are related to the quantum fluctuations of ions along the axial direction via

 δ^Rzj=∑νbzνj√ℏ2mωzν[^azν+^a†zν]. (22)

It is obvious that the Hamiltonian is invariant in the rotating frame and the laboratory frame since each phonon mode oscillates along the rotation axis and hence they are not influenced by the overall rotation of the coordinate system.

### ii.3 Planar phonon modes

So far, the planar phonon modes have not been utilized for quantum simulation due to the complexity of the normal modes as well as the complication introduced by the rotation of the ion crystals as observed in the laboratory frame. We plan to give a detailed explanation of the nature of the planar phonon modes and hope this will stimulate further insights toward quantum simulation. Our procedure is based on previous work on Coulomb crystals in a magnetic field, but not in a rotating frame japanese (); russian_old (); chen_thesis (); baiko (). Let us expand the Lagrangian in terms of the eigenmodes of the hermitian stiffness matrix in Eq. (14), which are labeled by the index with (note that these eigenvectors are not the normal mode eigenvectors for the planar motion, they are a convenient basis to use for the normal mode equations of the planar motion, except for the one case noted below). These eigenvectors satisfy

 ∑kβKαβjkbβνk=mων02bανj, (23)

where the spatial indices are the site indices for each ion and or , denote the directions of the spatial displacements of the phonons in the plane, are the real eigenvalues of the real-symmetric matrix for the stable ion configuration. These eigenvectors do become the eigenvectors for the planar modes only in the case when (which occurs when the rotating frequency is equal to half the cyclotron frequency), where the eigenvalue problem simplifies to the traditional form for a phonon. This fact has been used for an analysis of the phonons in a Penning trap Jake (). The basis vectors can be chosen to be real and they are also orthonormal . Therefore, the component of the displacement of the th ion in the th direction, , can be expanded in this basis as

 δRαj(t)=2N∑ν=1ζν(t)bανj, j=1,2,⋯N. (24)

where is the dynamic variable. The Lagrangian for the planar normal modes in this basis can be rewritten following Eq. (13) as

 LPph=12m[˙ζ2ν−ων02ζν2]−12ζν˙ζν′Tνν′, (25)

where the Einstein summation convention is used for repeated indices and the matrix elements are related to the antisymmetric matrix by the unitary transformation in which (where and all others vanish) due to the velocity dependent term in Eq. (13). Based on the Lagrangian , we can identify that the canonical momentum for mode is related to the mechanical momentum by

 Pν=∂LPph∂˙ζν=Πν−12ζν′bαν′jTαβjkbβνk=Πν−12ζν′Tν′ν (26)

Finally, we arrive at the Hamiltonian for the planar modes as

 HPph=Pν˙ζν−LPph=12mΠν2+12mων02ζν2. (27)

Notice that the Hamiltonian is not diagonal in the canonical conjugate variables and but is diagonal in the variables and . We need to perform a canonical transformation in order to diagonalize the Hamiltonian in the new canonical conjugate variables and we then can proceed with the standard quantization procedure (as in the axial phonon case).

Our strategy is to first quantize the Hamiltonian with respect to the canonical variables and then choose a nonstandard linear combination of coordinate and momentum to create the appropriate raising and lowering operators (this strategy was described by Baiko baiko ()). Hence we use the following fundamental commutation relations

 [ζν,ζν′]=0, [Pν,Pν′]=0, [ζν,Pν′]=iℏδνν′, (28)

which lead to the following commutation relations for the mechanical momenta

 [Πν,Πν′] = [Pν,12ζ¯νbα¯νjTαβjkbβν′k]+[12ζ¯νbα¯νjTαβjkbβνk,Pν′] (29) = −iℏTνν′≠0,

in which the antisymmetric property is used. To diagonize the Hamiltonian in the conventional second quantized form with collective mode frequencies for the stable ion configurations, we look for the phonon creation operator in the presence of the magnetic field as the unconventional superposition of the operators and

 a†λ=ανλΠν+βνλζν,aλ=αν∗λΠν+βν∗λζν, (30)

with and (possibly complex) numbers. The operators and must satisfy the commutation relation As a consequence, the Hamiltonian is required to satisfy the commutation relation

 [HPph,a†λ]=ℏωλa†λ (31)

with the normal mode frequency for the planar phonon. By direct substitution of Eq. (30) into Eq. (31), the coefficients and must satisfy the following coupled eigenvalue equations

 ωλανλ = −imβνλ−imTνν′αν′λ (32) ωλβνλ = imων02ανλ. (33)

Eq. (33) can be solved for in terms of via . Substituting into Eq. (32), we find the coefficients satisfy a quadratic eigenvalue problem (QEP)

 (mω2λδνν′+iωλTνν′−mων02δνν′)αν′λ=0 (34)

for the eigenvalue which is chosen to be a nonnegative frequency. Since QEP are not so common in physics, we discuss how to solve them in the appendix (see also Ref. QEP, ). In the simple form in Eq. (34), we can map the QEP onto a conventional linear hermitian eigenvalue problem in twice as many dimensions, which then allows us to use orthonogonality of eigenvectors and completeness to derive a number of nontrivial identities of the eigenfunctions. Details can be found in the appendix. For example, if we consider two different positive eigenvalues and , the corresponding eigenvectors satisfy the following orthogonality relation

 2N∑ν=1(ωλωλ′+ων02)αν∗λανλ′=0, λ≠λ′, (35)

which is derived in the appendix. The normalization of is fixed by the following commutation relation

 [aλ,a†λ′] = −iℏαν∗λTνν′αν′λ′+ℏmων02(1ωλ+1ωλ′)αν∗λανλ′ (36) =

which has been simplified by using Eqs. (34) and (35).

We are interested in representing the operator in terms of the raising and lowering operators and so we can relate the displacement in the plane to the collective mode operators and . Using Eq. (30), is related to and via

 αν∗λa†λ−ανλaλ = (αν∗λβνλ−ανλβν∗λ)ζν, (37) = 2imων02ωλ|ανλ|2ζν, = iℏζν, (38)

in which the relation , Eq. (65) and Eq. (66) are used. Hence,

 ζν=−iℏ∑λ(αν∗λa†λ−ανλaλ). (39)

This expression is crucial for the discussion of the planar modes in a quantum simulation. The component of the ion displacement associated with phonon coherent state in Heisenberg picture can be evaluated as

 ⟨δRλβj⟩ = 2ℏ∑λ:ωλ>0Re{iανλbανjϕλeiωλt} (40) = −2ℏ∑λ:ωλ>0|ϕλ||αβνj|sin(ωλt+δλ) (41)

where is the complex eigenvalue for the coherent state () , is the phase angle for mode , and is the average phonon occupation for mode , is the projection of state at site along the orientation.

### ii.4 Effective spin Hamiltonian with axial phonon modes

We discuss the generation of effective spin-spin couplings from an optical dipole force within the context of the Penning trap experiment at NIST, which uses ions of Beryllium (modifications for other systems would be straightforward, but are omitted here). Consider the hyperfine qubit states = and = in a ion. The qubit level splitting in the presence of a magnetic field of magnitude Tesla (due to Zeeman effects) is approximately GHz. Even though the nuclear spin is not a spin singlet, the nuclear spin dynamics can be completely frozen out by optically pumping the nuclear spin to the single nuclear spin state throughout the duration of an experiment John (); John2 (); John3 (), which is what we assume is done here.

A spin-dependent force along the axis of the Bloch sphere is generated through a gate PJ (), which is implemented by a coupling of the qubit states with the excited states ( and ) via two off-resonant laser beams with different angular frequencies and wavevectors respectively. The AC Stark shift due to the interference of two linear-polarized beams for an ion labeled by generates the following spin-dependent force (after adjusting the polarization of the laser beams) John ():

 ^HOD=−N∑j=1FOδkzcos(δk^zj−μt)σZj, (42)

where the transverse wavenumber is determined by the magnitude of the wavevector difference between the two beams, is the angle between the beam with respect to the tangential orientation of the two-dimensional crystal plane, is the beatnote frequency given by the frequency difference of the two beams, and is the displacement of the transverse phonons along the axial direction . The Pauli spin operator shows that the force is equal in magnitude, but opposite in direction for each of the two spin states of the Beryllium ion (hence it is known as a spin-dependent force). In the Lamb-Dicke limit, we have , so the interaction can be further reduced to the following form

 HOD=−N∑j=1FOδ^Rzjsin(μt)σZj (43)

after dropping a term with no dependence. In this limit, the spin-dependent force is spatially uniform.

The effective spin Hamiltonian and the time dependent evolution of the wave function from the Hamiltonian have been rigorously studied when the system is cooled to low temperature to start from the phonon vacuum spin-spin-interactions (); Phonons () immediately before the onset of the quantum simulation. The evolution of the entangled spin-phonon states are captured by the time evolution operator

 U(t,t0)=exp[−iHTph(t−t0)/ℏ]exp[−i^WI(t)/ℏ]Uspin(t,t0) (44)

with the operator determined by where the interaction is the optical dipole interaction term expressed in the Heisenberg picture of the bare phonon states (also called the interaction picture) via

 ^VI(t)=exp[i^HTph(t−t0)/ℏ]HODF(t)exp[−i^HTph(t−t0)/ℏ], (45)

in which is the (time-independent) bare planar phonon Hamiltonian given in Eq. (20). As a consequence, we find the effective spin Hamiltonian is dictated by an Ising spin Hamiltonian

 Uspin(t,t0)=Ttexp⎡⎣−iℏ∫tt0dt′⎛⎝N∑j,j′=1Jjj′(t′)σZjσZj′⎞⎠⎤⎦ (46)

with time-dependent exchange integrals.

These spin-spin exchange integrals are given by spin-spin-interactions ()

 Jjj′(t) = F2O4mN∑ν=1bzνjbzνj′μ2−ω2zν (47) × [1+cos2μt−2μωzνsinωzνtsinμt].

The detailed phonon mode properties , , and determine the values of the different spin-spin interactions. The sign of the effective spin-spin interaction depends on the redness or blueness of the laser detuning with respect to each phonon mode and the range of the interaction can be tuned depending upon the magnitude of the laser detuning away from all the axial phonon modes. In our numerical discussion, we will show these effects with direct numerical calculations.

Note that the analysis becomes more complicated if there is a transverse magnetic field present in addition to the spin-dependent optical dipole force. This is needed for quantum simulations that employ adiabatic state creation. Nevertheless, we do not discuss the evolution of the system further here when there is a transverse magnetic field present. Details of the case of the Paul trap can be found in Ref. Phonons, .

### ii.5 Spin-dependent interaction by planar modes

Let us now discuss the situation where the momentum transfer from the two laser beams lies in the crystal plane. One complication of this set-up is due to the rotation of the crystal with respect to the laser beams in the laboratory frame of reference. Examining the planar-mode coupling is also important in order to estimate the errors for a quantum simulation mediated through the axial phonon modes when the momentum difference of the two Raman beams is not oriented precisely along the axis due to laser alignment errors. Nevertheless, one can also involve the planar modes on purpose for quantum simulation, which may be useful when the planar modes are far away from the axial modes in energy. To simplify our discussion, let us choose the orientation of the planar spin-dependent forces due to the off-resonant laser beams to be along the -axis in the laboratory frame.

Analogous to the earlier discussion for axial phonon modes, the AC Stark shift along the -axis for the planar phonon modes is associated with the effective momentum transfer due to the photons via

 ^HOD=−N∑j=1FOδkxcos(δkx^xj(t)−μt)σZj. (48)

Notice that the ion coordinates in the laboratory frame are determined by the rigid-body rotation of the ion crystals with the axis chosen along the orientation of the effective wavevector in the laboratory frame and fluctuations in the laboratory frame. Hence, in the Lamb-Dicke limit , the optical dipole interaction is described by

 ^HOD=N∑j=1FOδ^rxj(t)sin(δkxx0j(t)−μt)σZj (49)

where the phase for the planar modes is modulated by the rotation of the crystal for the ion coordinates in the laboratory frame given by with the static phase determined by the equilibrium configuration in the rotating frame with respect to the orientation of the spin-dependent force acting along the direction at the initial time . Hence, the optical dipole interaction in Eq. (49) can be expanded in harmonics of . We expect the effective spin-spin interaction due to the planar modes to be sensitive to the rotational angular frequencies and the laser detuning . The value of the spin-spin interaction is sensitive to the structure of the different planar modes. Using the standard partial wave expansion Partial (), we can expand the function in harmonics of as

 sin(δk