Towards analog quantum simulations of lattice gauge theories with trapped ions

Zohreh Davoudi,1,2 Mohammad Hafezi,3,4 Christopher Monroe,3,5 Guido Pagano,3,5,6 Alireza Seif,3 and Andrew Shaw1 1Maryland Center for Fundamental Physics and Department of Physics, University of Maryland, College Park, Maryland 20742, USA 2RIKEN Center for Accelerator-based Sciences, Wako 351-0198, Japan 3Joint Quantum Institute and Department of Physics, University of Maryland, College Park, Maryland 20742, USA 4Department of Electrical and Computer Engineering and Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland 20742, USA 5Joint Center for Quantum Information and Computer Science, University of Maryland, College Park, Maryland 20742, USA 6Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA


I. INTRODUCTION
The invariance of physical systems under local transformations of fields leads to fundamental constraints on how matter fields interact, and introduces new bosonic degrees of freedom, the gauge fields. Gauge field theories coupled to matter are responsible for a wide range of phenomena in nature, and permeate condensed matter, nuclear, and particle physics. In the case of gauge theories comprising the Standard Model (SM) of particle physics, progress in perturbative tools has enabled predictions for high-energy experiments at the Large Hadron Collider [1]. Furthermore, progress in nonperturbative tools has led to theoretical input for precision experiments in search of violations of fundamental symmetries in nature, and to predicting hadronic excitations and their internal structure [2]. Nonetheless, the computational complexity of such studies grows significantly with the system size. In the strong-coupling regime, in which nonperturbative Monte Carlo sampling of quantum vacuum configurations is a common practice, questions such as the nature of the phase diagram of finite-density systems and the real-time dynamics of matter remain unanswered. It is therefore essential to explore a broader set of computational approaches, including those based on quantum simulation and quantum computation [3][4][5], to tackle these problems.
While the idea of simulating a quantum system using another quantum system with a higher level of control dates back to Feynman [6], only the experimental advancements in recent years have enabled powerful and sizable quantum simulations to become a reality. As in the case of classical computations, digital computations on quantum platforms may be the ultimate solution to all computational problems, including quantum simulations of physical systems. However, in the era of noisy intermediate-scale quantum (NISQ) computing [7], the number of high-fidelity operations that can be performed on a device can be highly constrained by the short coherence time of the quantum state. As a result, the digitalization of complex dynamics, such as those associated with gauge field theories [8][9][10][11], can be limited to small system sizes and short evolution times. It is therefore important to seek alternative approaches in the NISQ era. An interesting possibility is offered by analog simulations, in which the native Hamiltonian of the controlled quantum system is engineered to be mapped to that of the target system. The quantum operations are then naturally implemented once the system is prepared to evolve according to the desired Hamiltonian.
Among the most compelling platforms for analog simulations of quantum systems, including those governed by gauge theories, are cold neutral atoms in optical lattices [12][13][14][15][16][17][18][19], optical tweezers [20,21], and trapped ions [22,23]. Simple lowdimensional field theories such as relativistic Dirac fermions, 1+1D 1 and 2+1D scalar and fermionic quantum electrodynamics (QED), and non-Abelian SU (2) and SO(3) gauge theories have been studied in this context, and proposals exist to map the desired lattice Hamiltonians (or their approximated forms) to that of the engineered Hamiltonian of neutral atoms in optical lattices . Recent implementations of simple static and dynamical gauge theories with neutral atoms in optical lattices [46][47][48][49][50][51][52], however, demonstrate the challenge of simulating more phenomenologically relevant gauge theories. Given the current size of controlled quantum systems, only a small number of degrees of freedom can be studied, leading to unavoidable truncations in the Hilbert space of a gauge theory that lives in a continuous infinite-volume spacetime. Such a limitation is present in other digital and analog quantum platforms as well. It is nonetheless important that theoretical developments in formulating and mapping gauge theories for a quantum simulation proceed alongside the continual experimental progress that aims to significantly improve capabilities and capacities of simulating platforms.
Trapped ions provide a pristine platform for quantum simulations [22]. Given the extremely high level of control enabled by laser-cooled and localized ions confined by electromagnetic fields, exceedingly high fidelities in state preparation and measurement, all-to-all entangling capability enabled through control over the excitations of the motional normal modes, and scalability potential of such systems, this architecture has become a primary candidate for digital quantum computations in recent years [9,[53][54][55][56][57][58][59][60][61][62][63][64][65][66]. A unique feature of the trapped ion architecture is that global addressing of the ions using a few laser beams allows the realization of tunable longrange spin-spin interactions in the chain. With no need for individual addressability, systems of a few tens of ions have been successfully realized, and analog simulations of sizable quantum spin systems are made possible [67][68][69][70][71][72][73][74][75]. More complex quantum many-body systems, such as those described by gauge field theories, require either some degree of individual addressing or higher-order spin interactions among different species, as put forward in several proposals for simulating the relativistic Dirac equation [5,[76][77][78] a quantum field theory of scalar fields [79], and 1+1D QED [9,80]. A milestone in quantum simulations of lattice gauge theories (LGTs) using trapped ions was achieved in Ref. [9], where the real-time dynamics of 1+1D QED in a system of four trapped ions in a linear trap was made possible through a digital protocol, but the number of operations required for a Trotterized procedure prevented a long evolution time to be achieved in the presence of noise. While fully analog proposals exist for simulating simple low-dimensional LGTs [80], none have been implemented so far due to technical limitations.
It is important to classify gauge field theories of interest in terms of whether analog simulation of their dynamics is feasible given current technology. It is also essential to investigate whether fully analog implementations can circumvent the accumulated noise due to digitalization [81], and whether the noise in an analog setup can be effectively mitigated. Finally, it would be beneficial to assess the practicality of existing ideas, and to develop new proposals for extending the quantum toolkit of trapped ions, to enable a one-to-one mapping between the engineered Hamiltonian of the ionlaser quantum system to the dynamics of a fermionic system coupled to gauge degrees of freedom (bosons). This paper is a first step in addressing these questions. Here we focus on identifying goals that can be achieved in the near term, by specifying, in detail, practical proposals for a range of accessible gauge theories.
The gauge field theories studied in this paper are examples of the theories whose discretized formulations can be mapped entirely to systems with spin-1 2 degrees of freedom. One example is 1+1D quantum electrodynamics (Schwinger model): this model has similarities to quantum chromodynamics in 3+1D, including exhibiting a nontrivial vacuum. The second example is 2+1D Abelian Chern-Simons theory coupled to matter fields: this model is an example of a topological gauge theory with applications in many areas of physics. The last example we consider is 2+1D Z 2 gauge theory with a nontrivial phase diagram on a lattice, including exhibiting confinement. We discuss the mapping of these theories to spin systems, and present experimental protocols for realizing these interactions in current and near-term ion-trap systems. In order to provide a reference for upcoming implementations in the case of the Schwinger model, detailed examples for four and eight fermion-site theories will be presented.
A linear chain of trapped ions is often viewed as a platform for simulating spin-1 2 systems in 1+1D. However, once such a system is augmented with individual addressing, it offers far more possibilities for quantum simulations of arbitrary spin systems, including those in higher dimensions. Such proposals have been put forth in Ref. [82] and are explicitly taken advantage of in the current work to simulate the dynamics of the LGTs mentioned above. We also demonstrate the accessibility of nearly perfect nearest-neighbor interactions by simply controlling the lasers' phases and intensities on each ion, and demonstrate the sensitivity of the evolution to the imperfections of the engineered Hamiltonian in the case of the Schwinger model. By controlling intensities, phases, and frequencies of laser beams addressing each ion, a highly accurate mapping to spin-spin Hamiltonians with arbitrary interaction profiles is enabled. An important feature of the protocols devised in this work is a thorough optimization procedure that maximizes the closeness to the desired Hamiltonian, while simultaneously minimizes errors stemming from residual couplings to motional excitations. The proposed experimental scheme will have applications beyond the examples discussed and is a general protocol for realizing interesting spin systems described by a Heisenberg Hamiltonian in arbitrary dimensions.
We must emphasize that the sole goal of this work is to engineer effective Hamiltonians in an ion-laser system to realize, with high accuracy, the Hamiltonians of several lattice gauge theories. What is achieved is a protocol for applying the fundamental unitary operation e −iHt in a quantum simulation, which is required for studying a range of dynamical correlation functions in the theory, as well as state preparation in a class of protocols. The ion-laser system, once tuned to reproduce the Hamiltonian of the desired system, can be used to study all quantities of interest in the theory when combined with a full quantum-simulation scheme.
The paper is organized as follows. Section II includes details of ion-laser Hamiltonian considered in the scheme of this work and presents the effective Hamiltonian obtained, its range of validity, and the associated undesired contributions that must be minimized subsequently. The two associated Appendices A and B offer details on a particular experimental platform, and a scheme that eliminates an unwanted bias term in engineering the effective Hamiltonian. The full evolution operator is further detailed in Appendix C. Section III presents the example of the lattice Schwinger model, its purely spin representation, and explicit experimental proposals for simulating four and eight fermion-site theories. The former case is implemented with a single detuning for each set of the lasers used, while the latter takes advantage of a multifrequency, multiamplitude scheme, requiring a thorough optimization of interaction couplings. Additional results on the eight fermionsite theory are presented in Appendix E. The results of the numerical evaluation of the full evolution operator up to the order considered are presented in another associated Appendix (Appendix F) as well as in the Supplemental Material [83]. Section IV presents examples of LGTs in higher dimensions and their dual spin representation, along with discussions on their amenability to the quantum simulation scheme of this work. We conclude in Sec. V by highlighting the differing features of the scheme presented here compared with the previous work, the significance of the results obtained, and future extensions that may enable addressing a wider class of gauge theories.

II. 1D CHAIN OF TRAPPED IONS AND ENGINEERED EFFECTIVE INTERACTIONS
Consider N ions confined in a radio-frequency Paul trap [84]. The "qubit" in this system can be encoded in two stable internal levels of the ion, denoted in the following as |↑ and |↓ . These states are separated in energy by an angular frequency ω 0 (with Planck's constanth = 1 here and in the rest of the paper). Coherent operations on spin degrees of freedom are realized through stimulated Raman transitions using two laser beams with a momentum-vector difference k. The physics of ion-laser interactions and the single and two-qubit manipulations in an ion trap is well known [23,67,69,[85][86][87]. However, the involved evolution of the system under multiple pairs of Raman beams, which are needed for engineering the Hamiltonians of models considered here, requires a few technical novelties, and warrants a dedicated discussion which will follow in this section. For clarity in the presentation, further details of the proposed scheme and a number of involved analytical forms will be provided in the appendices.

A. Devised scheme and intrinsic Hamiltonian
The ion-laser interaction Hamiltonian for a system of N trapped ions can be written as [69] Index L in Eq. (1) runs over n L pairs of Raman beams.
(i) L is the Rabi frequency associated with the laser L. ϕ (i) L denotes the phase difference between the two lasers in each pair of Raman beams, ω (i) L is the difference in their angular frequency, namely, the beatnote frequency, and k (i) L is the difference in their momentum k-vector. In general, each ion is addressed with multiple pairs of Raman beams individually [hence the superscript (i) on quantities], requiring both amplitude and frequency control of the beams. Such individual addressing of the ions is widely used in digital ion-trap platforms and can be ported to analog platforms in upcoming experiments. r (i) denotes the displacement vector of ion i from its equilibrium position. The Pauli matrices σ (i) act on the quasispin of ion i, and α 0 , α 1 , α 2 , and α 3 are constants related to the spin-dependent forces on the two states of the qubit [69] and are controlled by the intensity, geometry, and polarization of the laser beams; see Appendix B for further details.
We assume that the confining potential is sufficiently stronger along the transverse axes of the trap so that the ions form a 1D crystal in space. With appropriate anharmonic axial confinement forces, the ions can be nearly equally spaced [88,89], with a typical spacing between adjacent ions of a few micrometers. Due to the long-range Coulomb force among the ions and the common trapping potential applied, the motion of the ions can be described in terms of a set collective normal modes. Then r (i) in Eq. (1) can be expressed in terms of phononic degrees of freedom, whose excitation energies are quantized in units of the normal-mode frequencies of the system. For the Hamiltonians of gauge theories considered in this work, it is necessary to introduce multiple pairs of bichromatic Raman beams directed at each ion, such that each pair couples to only one set of the three independent sets of normal modes. Such a scheme can be achieved with N individual beams and three global beams. Each of the individual beams will have three frequencies 2 that are tuned sufficiently apart such that each frequency will drive the qubit only by pairing with one of the global beams. This setup will allow to tune independently Hamiltonians acting along orthogonal directions of the Bloch sphere with negligible undesired cross couplings as shown below. The chosen directionality of the beams can ensure that each global-individual pair will result in a net k-vector along one of the three orthogonal principal axes of the trap, X, Y , and Z; see Fig. 1. 3 Here X and Y denote the most-confined directions in the trap, which will have the same normal-mode spectra for symmetric traps commonly used. These will be denoted as transverse directions. The least-confined direction is denoted as Z and is named the axial direction. Consider now the ion-laser system in the interaction picture, in which all excitations arising from the free Hamiltonian are rotated away by frequencies of the order of ω 0 , ω T m , and ω A m . 4 a m (a † m ) annihilates (creates) a phonon excitation of the transverse normal mode m with angular frequency ω T m along the X direction of the trap, k I = k IX . Similarly, b m and c m (b † m and c † m ) are, respectively, the phonon annihilation (creation) operators for the axial normal modes along the Z direction, k II = k IIẐ , and the transverse normal modes along the Y direction, k III = k IIIŶ . 5 The corresponding normal mode frequencies are denoted as ω A m and ω T m . Different superscripts are introduced to distinguish the transverse and axial normal modes which have different spectra. Finally, in the Lamb-Dicke regime where k (i) r (i) 1, and when the laser frequencies are chosen such that all transitions except for those near the first sideband transitions 6 are far 4 Although the axial modes are generally low in frequency, such a rotating-frame approximation is still valid as long as lasers' detunings from these modes remain small compared to the sideband Rabi frequencies of the axial motion. 5 At this point, such assignments of a given set of normal modes to one of the Hamiltonians in Eqs. (3)-(5) appear arbitrary. The rationale behind the choices made will become clear in applications of the scheme to nearest-neighbor Hamiltonians considered in this work; see Sec. III. 6 The n th blue (red) sideband transition for mode m adds (removes) n quanta of motion each with frequency ω m . off-resonant, the three sets of Raman-beam pairs at each ion induce the laser-ion Hamiltonians of the form where σ (i) , and the tilde over the Hamiltonians implies the use of the rotated frame described above. Here it is assumed that On the other hand, for the Hamiltonian H III , it is assumed that |μ III | ω 0 where μ III ≡ − ω III = ω III . Further, two distinct Raman-beam phase differences are assigned to each of the red (unprimed) and blue (primed) detuned frequencies of the beam. η (i) m is the (normalized) normal-mode eigenvector components between ion i and mode m, and M denotes the mass of the ion. Similarly, η (i) for the axial and transverse modes, respectively. For each pair of Raman beams L, the same k L vector is applied at the location of each ion. α 1 = 1 2 and α 2 = 0 correspond to the well-known Molmer-Sorenson scheme, already applied in a number of experiments. In order to eliminate a bias σ z interaction arising from H III , it is essential that α 0 is set to zero. With the scheme presented in Appendix B, it is shown that one can achieve this requirement by tuning the Raman-beam frequencies and polarization vectors. We further set α 3 = 1 4 for consistency between the effective spin-spin couplings arising from H I , H II , and H III . 7 Now by setting the phases of the blue-and redsideband detuned beams to ϕ (i)  7 There will be no ambiguity in the overall constants in the Hamiltonian. Rescaling these coefficients by a constant means the Rabi frequencies must be rescaled accordingly so that the expected strength of the state-dependent force is produced on a given ion, and with given choices of the internal levels for the qubit.
H III can be seen to be proportional to σ (i) x , σ (i) y , and σ (i) z , respectively.
Finally, an effective longitudinal magnetic field can be introduced at the location of each ion by another N sets of beams inducing a Stark shift to be tuned to the desired value of the magnetic field. Alternatively, a B z field can be generated with the existing sets of Raman beams, i.e., by shifting the frequency of red-and blue-detuned beams by B (i) z . This can be seen by noting that if the rotating frame that led to Eqs. (3)-(5) is assumed to rotate with the Hamiltonian H 0 + is generated, but at the cost of the following change: z to the laser detuning in the first and second occurrences of μ I in Eq. (3), respectively. Similarly, μ II must be replaced by z in the first and second occurrences in Eq. (4), respectively. The laser detuning μ III , on the other hand, remains unchanged. Note that this scheme requires a frequency control, as the detunings are now generally different at the location of each ion.

B. Time evolution and effective Hamiltonian
With the Hamiltonians in Eqs. (3)-(6), an evolution operator can be formed by applying a Magnus expansion, taking into account all contributions up to and including O(η 2 , ηB z ) in the exponent: where The definitions of the rest of the functions in Eqs. (7)-(10) are provided in Appendix C. When B (i) z = 0, all contributions proportional to phonon creation and annihilation operators in the exponent in Eq. (7) are bounded in time, provided that μ I = μ II = μ III . As a result, an effective Heisenberg model can be achieved when so that the terms linear in time in Eq. (7) (those proportional to χ (α) i, j ) dominate the evolution. In such a limit, (15)], and other contributions will be subdominant. For practical (noisy) implementations, one needs to minimize the spin-phonon entanglement arising from the first term in the exponent in Eq. (7) at early times. This is achieved with |η (i) i,m (t ) in Eqs. (8) and (9) develop an oscillatory time dependence but with a linear growth in the magnitude of its amplitude. These terms are proportional to Assuming that the magnetic field is comparable in size to the effective spin-spin couplings, such contaminating terms do not severely impact the desired evolution as long as Unfortunately, this condition limits the size of (effective) magnetic fields that can be studied in models considered below. Nonetheless, a range of interesting possibilities can still be explored.
Under the conditions described above, the time-evolution operator in Eq. (7) can be approximated as where As a result, the individual-addressing scheme proposed here enables analog quantum simulations of a rather generic Heisenberg spin model. The spin-spin coupling matrices in Eq. (12) are derived from discussions above (see also Appendix C) and read Here is the recoil frequency of the ion given the lasers L = I, II, III.
It is worth noting that despite the case of a usual Molmer-Sorenson transition where the starting Hamiltonian is proportional to σ x , the Magnus expansion in the scheme described above is not cut off at any order in the Lamb-Dicke parameter, due to the nonzero commutation of Pauli operators in Eqs. (3)- (6). It is therefore important to ensure that not only | 1 as stated before, but also This guarantees that contributions from the pth-sideband transitions are suppressed compared to the first-sideband transitions. These conditions are easier to satisfy for transverse modes than the axial modes. This is because the axial modes have lower frequencies, and their corresponding Lamb-Dicke parameters are larger. Finally, one notes that coherent operations on a single spin correspond to the zeroth-order terms in Eq. (1) in the Lamb-Dicke limit, and with ω (i) L = ω 0 . Hence, the laser frequencies applied must be far detuned from such "carrier transitions" of the ions.

III. OPTIMIZED SPIN-SPIN HAMILTONIANS IN AN ION TRAP: 1+1D SCHWINGER MODEL
A unique test bed for exploring theoretical and experimental proposals for quantum simulations of gauge theories is the 1 + 1D QED, i.e., the Schwinger model. It is an Abelian gauge theory, hence avoiding complexities of its non-Abelian counterparts. It is also a low-dimensional theory, allowing numerical and experimental studies of its approximate dynamics with finite resources. Despite these simplifications in the formulation, the theory exhibits rich properties, similar to those seen in more complex theories such as QCD. In particular, phenomena such as confinement and spontaneous symmetry breaking arise in the model. The spontaneous creation of electron-positron pairs in the time evolution of the "vacuum" exhibits a clear signature of such nontrivial dynamics. Since the time evolution of quantum states is, in general, a computationally intractable problem with classical Monte Carlo methods, addressing such a problem using a quantum simulation platform is of significant value; see Refs. [8,9] for digital implementations.
The strong-coupling dynamics of the Schwinger model can be studied through nonperturbative LGT methods. In the staggered formulation of Kogut and Susskind [90,91], the (scaled) lattice Hamiltonian takes the form where n ( † n ) is a one-component fermion field that creates (annihilates) an electron on the odd site while annihilates (creates) a positron on an even site. Due to this distinction, there is a staggered mass term in the Hamiltonian, with the fermion (scaled) mass μ. θ n is the U (1) gauge potential with the corresponding gauge link e iθ n originating at site n. The latter is introduced in the Hamiltonian to render the fermion hopping (kinetic) term gauge invariant. The pair creation and annihilation in the theory originates from this term. The corresponding electric field at site n is denoted as L n (with the operator relation [θ n , L m ] = iδ n,m ), which adds a contribution to the Hamiltonian due to the energy stored in the electric field. The Hamiltonian in Eq. (16) is written in units of ag 2 /2, where a denote the lattice spacing and g is the original fermion-gauge field coupling. The dimensionless parameters x and μ are related to dimensionful parameter g (with mass dimension one) and the original mass m via: x = 1/(ag) 2 and μ = 2m/(ag 2 ). 8 The familiar Jordan-Wigner transformations n = l<n (iσ (l ) z )σ (n) − and † n = l<n (−iσ (l ) z )σ (n) + can be applied to Eq. (16) in order to map the fermionic degrees of freedom to those of a qubit. A unique feature of the lattice Schwinger model with open boundary condition is that the remaining degrees of freedom that are bosonic, namely, gauge links and electric field, can be entirely eliminated in favor of new spin-spin interactions. Explicitly, by performing gauge transformations σ (n) ± → l<n e ±iθ l σ (n) ± , and further imposing the Gauss's law L n − L n−1 = 1 2 [σ (n) z + (−1) n ], the Hamiltonian becomes [9,92,93] Here 0 is the electric field flux into the first lattice site which can be set to zero without loss of generality. To make explicit the mapping of this Hamiltonians to that of the Hamiltonian of the ion-laser system in our proposed scheme, Eq. (12), one can note that Eq. (17) can be rewritten as where Hamiltonian in Eq. (12) is equal to the Schwinger Hamiltonian in Eq. (18). This is an optimization problem that can be solved straightforwardly provided that multiple laser frequencies are used with each set of beams each with a corresponding Rabi frequency [82]. With n μ L number of beatnote frequencies on each pair of lasers L, the total number of free parameters is N n μ L , while the the number of independent nonzero elements in each J i, j coupling matrix is N (N − 1)/2. Empirically, it is seen that a solution to the optimization problem can be achieved with n μ L N. It is, however, conceivable that in the first generation of experiments planned, only the amplitude control of Raman beams will be a reality. As a result, we first focus on experimental proposals that do not require a frequency control.

A. A single-detuning and multiamplitude scheme
With a single beatnote frequency on each pair of Raman beams, the Schwinger Hamiltonian on small lattices can be realized with good accuracy. For this example, an ion trap consisting of 171 Yb + ions will be considered. The specifications of this system are presented in Appendix A. Consider the case of N = 4, and further set the values of the parameters of the Schwinger Hamiltonian to x = 6 and μ = 1. The Hamiltonian H (xx) can be achieved by first noting that a certain detuning from the CM transverse mode with the same amplitude on each ion produces the coupling matrix shown in Fig. 2(a). This matrix can be systematically turned into a nearest-neighbor form: the slope of the decline in the strength of nearest-neighbor couplings from the center of the chain can be determined and be used to systematically adjust the Rabi frequencies in such a way that an equal strength is achieved on all J i, j with |i − j| = 1, as demonstrated in Fig. 2(b). The most accurate nearest-neighbor Hamiltonian achieved with this procedure presents a ∼3% contamination on the non-nearest-neighbor elements and no contamination on the nearest-neighbor elements.
As mentioned in Sec. II, the H (yy) effective Hamiltonian is chosen to arise from the Raman beams that address the axial modes of the ions. If the transverse modes were to be addressed, the Raman beams would have to be detuned from the modes by the same amount as those for the H (xx) Hamiltonian, as these appear with the same coupling in the Schwinger Hamiltonian. This, however, would cause the dynamics to deviate from the effective Heisenberg model in Eq. (12), given the nonzero commutations between H I and H II in Eqs. (3) and (4), generating phonon-dependent terms that grow (or decline) linearly with time. Such contaminations are circumvented by producing the H (yy) Hamiltonian with the Raman beams that couple to the axial modes. Note that the axial modes have a very different frequency spectrum compared with the transverse modes. The same procedure as for the H (xx) mapping can be used to find the values of the laser beatnote and Rabi frequencies that generate a nearest-neighbor interaction with these modes; see Fig. 3. As discussed at the end of Sec. II, a critical check is to ensure the higher-sideband contributions to the applied Molmer-Sorenson scheme are not significant given the low normal-mode frequencies in the axial direction and given the laser frequencies applied. It can be shown that the largest contribution from these higherorder sidebands is only a few percent of the contribution from the first sideband and will be ignored in the current proposal.
An effective H (zz) Hamiltonian that matches that of the Schwinger model can be achieved with a single beatnote frequency and by addressing the other set of transverse normal modes of the ions. Here the values shown in Fig. 4 allow the J i, j coupling to be tuned to the desired values with below-percent accuracy. However, in contrast with the case of nearest-neighbor Hamiltonians, the procedure that finds the adjusted Rabi frequencies for H (zz) is not systematic, making it challenging to generalize such an ad hoc tuning procedure to a higher number of ions. Finally, an effective H (z) Hamiltonian can be induced using N sets of Raman beams with their Stark shift tuned to reproduce H (z) of the Schwinger Hamiltonian in Eq. (22). The values of the effective magnetic fields that are required given the chosen parameters of the model are depicted in Fig. 5.
It is crucial to verify that the laser parameters found in such a mapping do not violate the conditions enumerated in the previous section, and the true dynamics is that dictated by the effective Heisenberg Hamiltonian in Eq. (12). This check can be done by a numerical evaluation of all contributions to the exponent of the evolution operator in Eq. (7), up to and including O(η 2 , ηB). Here we assume that the experiment can be initiated in a state with zero phonon occupation in all modes. The results of this investigation are shown in Fig. 14 of Appendix F for the first ion, and in the Supplemental Material [83] for the rest of the ions. As shown, the dominant source of error is related to the nonzero commutations of H B in Eq. (6) with H I and H II in Eqs. (3) and (4), introducing effective magnetic fields along the x and y spin axes. These are a small fraction of the desired field along the z direction, but are, however, dependent upon the phonon occupation in the system.
Hamiltonians of the lattice Schwinger model for a larger number of fermion sites can be shown to be accessible through the single-frequency and multiamplitude scheme described, but deviations from the exact Hamiltonian can be significant. For N = 10 and the nearest-neighbor Hamiltonian with transverse modes, the best parameters found give rise to errors as high as ∼20% in the non-nearest-neighbor elements. To investigate the effect of inexact Hamiltonians on the dynamics of the Schwinger model, we have studied the four fermion-site Schwinger model is then considered. The quantity of interest here is the vacuum persistence amplitude (VPA), defined as the (square) of the overlap of the state of the system at time t, |ψ (t ) with the "vacuum" (a state in the physical sector of the theory with no net electron-positron pair), |ψ (vac) . This quantity is plotted for select times in Figs. 6(b)-6(e) for all the 20 inexact Hamiltonians used in the evolution. A procedure is described to estimate a mean and uncertainty band from the most accurate Hamiltonians employed. Nonetheless, as is seen in Fig. 6(a), during certain times, the estimate of VPA deviates significantly from the expected result, and this feature is amplified at longer times.
This observation promotes adopting a multifrequency and multiamplitude scheme, 9 as proposed previously in Ref. [82] in the context of quantum simulation of the Ising model on two-dimensional lattices. With this scheme, mapping of the effective Hamiltonian of the ion-laser system to that of the Schwinger model can be achieved with unprecedented accuracy, as is shown in the following.  [83]. These values are obtained from a numerical study and do not represent experimental findings.

B. A multifrequency and multiamplitude scheme
The extension of the formalism presented in Sec. II to a multifrequency scheme is straightforward. For example, the effective spin-spin coupling engineered by Raman pairs I generalizes to where n μ I is the number of beatnote frequencies, and where each detuning μ I,m is associated with the Rabi frequency (i) I,m . 10 Similarly, the J (yy) i, j and J (zz) i, j coupling matrices can be obtained by replacements μ II → μ II,m , (i) II → (i) II,m , μ III → μ III,m , and (i) III → (i) III,m , where a summation over m is assumed. For J (yy) i, j , one must replace ω T m with ω A m . More generally, the full time evolution operator in Eq. (7) can be constructed by performing the changes described in the ionlaser Hamiltonians in Eqs. (3)-(5). This introduces additional 10 We remind that the effective spin-spin Hamiltonian arises from a bichromatic pair of Raman beams, one detuned by −μ I,m (reddetuned) and one by μ I,m (blue-detuned) from the carrier transition, see discussions after Eq. (5).
off-resonant terms that would scale as the number of beatnote frequencies introduced. One therefore needs to ensure that the cumulative effect of such terms remains negligible compared with the desired effective Heisenberg Hamiltonian. Figure 7 demonstrates the success of this scheme in an accurate generation of the long-range part of the Schwinger Hamiltonian, H (zz) , for the case of N = 8 ions. Here the corresponding optimization problem is solved (see Appendix D for details), and the desired effective spin-spin Hamiltonian is achieved with errors that are comparable with the machine precision. The laser frequencies are fixed such that μ I, ), with f s = −0.5, and where m runs from 1 to n III = 7, see the lower-right plot of Fig. 7. 11 The corresponding Rabi frequencies at the location of each ion are plotted in the upper-right plot of Fig. 7. As is seen, a perfect agreement between J (zz) i, j and that in the Schwinger model with x = 6 and μ = 1 is achieved. The reason for choosing a large value of the coupling x in the original theory is to minimize 11 In the convention of this work, the normal-mode frequencies are ordered in a set from the highest value to the lowest value. Therefore for the axial mode, ω A N denotes the c.m. mode, while for the transverse mode, the c.m. mode is ω T 1 . Because of this convention, the normal-mode eigenvectors b (i) m must be ordered accordingly for the transverse and axial modes.  3) and (5). Note that the desired effective B z field in the Schwinger Hamiltonian grows with N even in the limit μ = 0. Hence, in order to keep the undesired contribution small compared with the effective Hamiltonian, the strength of the nearest-neighbor terms is taken to be stronger by setting x = 6. As is shown in Appendix F for the first ion, and in the Supplemental Material [83] for the rest of the ions, all the contributions to the exponent in the full time-evolution operator (up to the order considered) are small (and mostly bounded) compared with those that constitute the Hamiltonian of the Schwinger model. The laser parameters for a nearly exact engineering of H (xx) , H (yy) , and H (z) are shown in Figs. 11-13 in Appendix E. It must be noted that the optimization problem in all cases is solved under two constraints: (1) the sum of Rabi frequencies at the location of each ion is less than or equal to 2 MHz and (2) the contribution to the full evolution from the first-order terms, those proportional to coefficients α (x) i,m , α (y) i,m , and α (z) i,m in Eqs. (8)-(10), remains below 0.5 at several random times up to 1 ms, see Appendix D.
To summarize, we have provided detailed experimental protocols for a fully analog simulation of the Schwinger model for given parameters with (1) a scheme that requires only individual amplitude and phase control of the laser beams and engineers an approximate Schwinger Hamiltonian and (2) a scheme that takes advantage of individual amplitude, phase, and frequency control and engineers the desired Hamiltonian with great accuracy (up to errors associated with the difference between the full ion-laser evolution and the effective Heisenberg model, which are nonetheless assured to remain negligible in the schemes proposed). It is clear that the second scheme can be easily applied to any number of ions at the cost of introducing a multitude of laser frequencies, the number of which grows with the number of ions. This can be already achieved with current technologies for up to ∼30 ions, and most importantly is scalable, as it involves a linear growth in the complexity of the classical control hardware of the experiment.
In the following, other examples of LGTs whose dynamics can be mapped onto a spin- 1 2 system will be discussed. The goal is to only point out the potential of an ion-trap quantum simulator in addressing more complex spin systems by providing examples of relevant gauge theories. Explicit scenarios for given ion-trap architectures are straightforward to obtain, following optimization strategies presented for the case of the Schwinger model.

IV. ANALOG SIMULATIONS OF SYSTEMS IN HIGHER DIMENSIONS WITH A 1D CHAIN OF IONS
With a generic Heisenberg model and an effective magnetic field engineered in Sec. II, it is clear that a wide range of couplings among spins can be tailored, as was demonstrated for the case of the Schwinger model. In particular, as seen in Sec. III, the H (αα) with α = x, y, z does not have to be necessarily nearest neighbor or of any particular form, as the multifrequency, multiamplitude scheme of this work allows an arbitrary J i, j to be produced. This observation implies that spin systems in higher spatial dimensions can be engineered as well, as was also noted in Ref. [82]. One only needs to map the points on a 2D or 3D lattice to a linear chain of ions along with their corresponding couplings. Of course, with a fixed number of ions in a given experiment, this means that the finite-size effects in the dynamics of the system under study will be larger, such as in the case of square and cubic lattices the spatial extent of the system will be N 1/2 and N 1/3 , respectively. Nevertheless, this possibility implies that a linear quantum system can be used as a platform for analog simulations of theories in any dimension, bringing the versatility of such an analog platform closer to its digital counterpart.

A. 2+1D Abelian Chern-Simons theory coupled to fermions
As an example of an interesting field theory in 2 + 1D, consider the Chern-Simons theory coupled to fermions. This theory is of broad impact on a range of problems in theoretical physics, from the theory of the integer and fractional quantum Hall effects to knot theory and parity anomalies in quantum field theory; see Ref. [94] for a review. Since the theory is topological in the continuum, the construction of a discretized counterpart of the theory turned out to be nontrivial as a lattice has explicit reference to a given coordinate system and metric. However, it has been shown [95,96] that one can still formulate a U (1) LGT that retains gauge invariance on arbitrary 2D planar graphs, has no local excitations (hence is topological), and in the long-wavelength limit approaches the Chern-Simons theory in the continuum. As is discussed in Ref. [96], a lattice formulation of the Chern-Simons theory is invaluable in investigations of fractional Chern insulators that occur in given lattice geometries. As a result, it is interesting to ask if a quantum-simulation protocol for this theory can be devised on the simulating platform of this work.
A known result [95] in the context of the generalized Jordan-Wigner transformation in higher dimensions is Fradkin's proof of equivalence between the spin-1 2 XY model on a 2D Bravais lattice and a Chern-Simons theory in 2+1D coupled to fermions, provided that the strength of the Chern-Simons' term in the Lagrangian density, is θ = 1 2π [95]. Here time is assumed to be continuous while spatial coordinates are defined on a square lattice, i.e., x = (t, n) where n is a vector whose components are integer multiples of the lattice spacing. 12 μ, ν = 0, 1, 2 with the zeroth direction being the time direction, a is a complex spinless fermion field, A μ is the gauge field, D μ = ∂ μ − iA μ is the covariant derivative, F μν is the field-strength tensor: and μνλ is the Levi-Civita symbol. Note that the A 0 field does not have any dynamics and can be set equal to zero with the choice of a temporal gauge. The physical sector of the theory, i.e., states that satisfy the Gauss's law, can be identified from the condition δS δA 0 = 0, where S is the action. These states then correspond to those for which It is also clear that the Hamiltonian of the theory vanishes in the absence of matter fields, which is a desired feature of the topological theory. In the presence of matter fields, the Hamiltonian corresponding to Eq. (24) is [a † (n)e iA j (n) a(n +n j ) + H.c.].
Note that the time dependence of the fields is now implicit considering the Hamiltonian equations of motion. As is shown in Ref. [95], the gauge links can be eliminated from the Hamiltonian with the use of Gauss's law, at the cost of changing the equal-time commutation relation of fermions. This is in fact a great advantage since when θ = 1 2π (or in general when 1 2θ is an odd multiple of π ), the new commutation relations are those of hardcore bosons, i.e., the spin-1 2 matrices. As a result, this procedure can be realized as a 2D generalization of the familiar Jordan-Wigner transformation. Explicitly, by performing the transformations a → e iφ a ≡ã and a † → a † e −iφ ≡ã † , where A j (n) ≡ φ j (n +n i ) − φ j (n), one arrives at where the following identifications are assumed: σ (n) + = a † (n), σ (n) − =ã(n), and σ (n) z = 1 − 2a † (n)a(n). Equation (26) clearly corresponds to an XY spin model. Note that a parameter h could be introduced to control the magnitude of the hopping term in the Hamiltonian.
To perform an analog simulation of such a 2D XY model within the scheme presented in Sec. II requires optimizing a matrix by performing a multifrequency, multiamplitude Molmer-Sorenson scheme using the transverse and axial normal modes of motion. For a 4 × 4 lattice in the target theory, a system of N = 16 ions can be used as is shown in Fig. 8, along with the required J i, j matrix. Obtaining the laser frequencies and amplitudes is a straightforward optimization process, as detailed in the previous section, and in fact machine precision accuracy can be achieved, as demonstrated in Ref. [82] for similar geometry and coupling profiles. Finally, we should remark that the full Hamiltonian in such a 2+1D Abelian LGT must include the energy stored in electric and magnetic fields, giving rise to the Maxwell-Chern-Simons theory [98,99]. 13 Aside from the question of what is the proper formulation of a discretized Maxwell-Chern-Simons theory, one needs to account for the full dynamics of the gauge fields by mapping them to those in an ion-trap quantum-simulation platform, which is beyond the scope of the present work.
B. 2+1D pure Z 2 lattice gauge theory Z N gauge theories are discrete Abelian gauge theories that given their simple underlying symmetry have long served as a testbed for gaining deeper perspectives on gauge theories. Despite their simple structure, they can have nontrivial phase diagrams exhibiting, e.g., a confining phase. In fact, since Z 3 is the center of the SU(3) group, the confinement in the Yang-Mills theory is attributed to the Z 3 symmetry. These gauge theories have been the focus of numerous theoretical and experimental proposals for quantum simulation of gauge theories, in particular using neutral atoms in optical lattices [36,52,101]. An interesting feature of Z N is its duality with spin models. This connection has been developed over decades [102], starting from Wegner's demonstration of such a duality for the case of a Z 2 LGT [103], and has inspired similar duality constructions for non-Abelian gauge theories such as SU(N) [104]. Further, recent work has suggested that the 4D Z 2 LGT provides a complete model for all classical spin models and all Abelian discrete LGTs [105,106].
The example that will be presented here is a 2+1D Z 2 LGT that is dual to a 2D Ising model and is therefore amenable to the quantum simulation protocol of this work. The Hamiltonian of the 2+1D Z 2 LGT can be expressed with a pair of conjugate spin operators {σ x (l ), σ z (l )}, where σ x (l ) = e iπE (l ) and σ z (l ) = e iA(l ) . Here l denotes a link on the 2D spatial lattice, A(l ) is the gauge field evaluated on link l with A(l ) = {0, π}. E (l ) is the corresponding "electric field" with E (l ) = {0, 1}. Note that in order to keep the presentation simple, we have not used bold-faced quantities for the two-dimensional vectors A(l ) and E (l ), as their directionality on the 2D plane is implicit from the directionality of the link arguments. The lattice Hamiltonian of such a pure gauge theory consists of "electric" and "magnetic" terms: Here the first (second) sum runs over all links (plaquettes) on the 2D lattice, and open boundary conditions are assumed. A plaquette is defined as the product of four gauge links staring from the lower-left corner and moving counterclockwise; see Fig. 9(a). The Hamiltonian in Eq. (27) remains invariant under a local gauge transformation which flips the sign of σ z on links sharing site n, but does not affect σ x on links sharing the same site. The Gauss's law corresponding to this symmetry defines the physical sector of the theory, namely, states for which the eigenvalue of the Gauss's law operator G(n) = n σ x (l n ) is unity, where l n denotes all the four links that meet at point n.
To establish a duality relation with the 2D Ising model, the gauge invariance can be taken into account to (1) fix the gauge conveniently such that σ z on all links along one of the spatial directions is set to unity and (2) use the operator identity G(n) = 1 in the physical Hilbert space of the theory to replace σ x along the same space direction as in (1) with those along the other direction. These two steps inspire the replacements σ z (l 1 )σ z (l 2 )σ z (l 3 )σ z (l 4 ) → σ x (p), and l l σ x (l ) → σ z (p) (which is allowed as the new {σ x , σ z } set has the same commutation relations as the original set). In the first replacement rule, p denotes the plaquette formed by links l 1 , l 2 , l 3 , l 4 , and in the second rule, it denotes the plaquette whose left bottom corner is the point at whichl starts. The product is over all links prior to and including linkl, and tilde is used to denote the space dimension for which the gauge remains unfixed. It is now easy to see that in terms of the new spin operators, the Hamiltonian in Eq. (27) can be written as where in the last line, n refers to points on the "dual" lattice defined by the center of spatial plaquettes in the original lattice; see Fig. 9. p, p in the first line denotes the nearestneighbor plaquettes. For further details on the expected phase diagram of the theories at different coupling regimes, see, e.g., Ref. [102]. The duality between Eqs. (27) and (28) allows to simulate the dynamics of a Z 2 LGT in 2+1D using a chain of ions in 1D whose interactions are tailored to correspond to the Ising Hamiltonian, as discussed in the previous example of this section. The correspondence between the original 2D lattice, the dual lattice, and the chain of ions is depicted in Fig. 9. Engineering the nearest-neighbor σ z ⊗ σ z interactions was detailed in Sec. II, and the additional global transverse magnetic field can be easily introduced by performing singlequbit rotations, with an angle determined by the coupling λ in the original theory.

V. CONCLUSION AND OUTLOOK
In this paper, we took on the question of how to best leverage the current technologies in ion-trap analog quantum simulators to engineer the Hamiltonian of gauge field theories. Towards this goal, gauge theories that can be experimentally realized in such platforms in the near future are enumerated and are shown to be amenable to a particular quantumsimulation scheme devised in this work. Experiments that will take advantage of the proposals of this work are being planned. The highlights of the scheme presented, and its promising applications, can be summarized as the following: (1) N sets of laser beams are used to address individual ions in a 1D chain. With the addition of three global laser beams, the Hamiltonian of a Heisenberg model can be engineered. Certain orientations and frequencies of the beams compared with each other (see Fig. 1 spin-spin interactions to be generated with negligible couplings among different Raman processes. Each set of lasers couples to one set of normal modes of motion (two transverse and one axial), allowing arbitrary spinspin couplings to be engineered. Our scheme is inspired by that presented in Ref. [67] but does not require an asymmetric trap in the transverse directions, as long as one is interested in a Heisenberg XYZ and XXZ models (see the example of the Schwinger model in Sec. III).
(2) The experimental scheme of this work offers the capability of engineering a range of interesting dynamics with a single beatnote frequency for each set of the lasers, denoted as μ L with L = I, II, III, but with tunable phases and with Rabi frequencies (i) L at the location of each ion. Moreover, introducing a frequency control to the system, as is common in the digital ion-trap platforms, allows arbitrary spin-spin Hamiltonians to be engineered with unprecedented accuracy.
(3) The frequency control allows an effective local magnetic field to be engineered via asymmetrically shifting the frequency of the red-and blue-detuned Raman beams, eliminating the need for introducing another N laser beams to induce local Stark shifts on the ions.
(4) Engineering an arbitrary Heisenberg Hamiltonian is enabled in this work by a thorough optimization procedure that minimizes the contributions arising from unwanted couplings to phonon excitations, contributions that drive the dynamics away from the effective spin-spin Hamiltonians. This is a crucial requirement for a reliable quantum simulation that is addressed for the first time in this work. The purely spin formulation of the lattice Schwinger model exists and corresponds to a Heisenberg XXZ model with both short-and long-range interactions, and with an effective local magnetic field. The optimization procedure described above was applied to this example with N = 8, and can be scaled straightforwardly to any number of ions. (5) In this work, equal-size nearest-neighbor couplings along the spin axesx andŷ are achieved through coupling to transverse and axial modes of the motion, respectively, eliminating any significant undesired coupling between the two resulting interacting Hamiltonians in the evolution given the Raman-beam detunings required. This feature does not demand the use of a strong effective magnetic field to induce such nearest-neighbor interactions [70,71,107], with its known limitations [108]. Although it may be challenging to implement such a scheme in larger chains of ions with low axial normal-mode frequencies, ideas such as that proposed in Ref. [109] may allow a scalable scheme in future investigations.
(6) Another feature of the proposed scheme is a high degree of flexibility in tuning the spin-spin interaction couplings of arbitrary forms along each axis of the qubit independently. This feature, which for example is not offered in single Molmer-Sorenson schemes [110], is shown to be particularly useful for engineering the Hamiltonians of gauge theories considered in this work. (7) The high level of control allows quantum simulation of models in higher dimensions. Two interesting examples of lattice gauge theories presented in this work (see Sec. IV) are Abelian Chern-Simons theory coupled to matter, and a Z 2 pure gauge theory, both in 2+1D, whose dynamic can be mapped to a planar Ising model with nearest-neighbor interactions. Such capability opens up the possibility of analog quantum simulations of systems beyond what has been possible to date.
A few directions can be recognized as natural extensions of the ideas presented in this paper. These include the following: (1) There are a range of methods that lead to a truncated angular-momentum representation of the gauge degrees of freedom in LGTs, such as the quantum link models [36,38,44], or the use of a tensor-network construction in Abelian gauge theories coupled to matter [44,111]. With the manipulation of a larger number of internal levels of the ions, the approach advocated in this paper can be applied to engineer interactions of spin systems with s 1 2 . An experimental realization of a spin Hamiltonian with s = 1 is presented in Ref. [112], and can be extended to allow quantum simulation of select gauge theories in spin-1 representations.
(2) For a wide range of phenomenologically interesting lattice gauge theories for which a purely spin representation does not exist, it is essential to extend the toolkit of iontrap analog simulation to leverage the control over phononic degrees of freedom. This will require further technological advancement on the experimental front, as well as new proposals for engineering gauge and gauge-matter interactions in a highly controlled spin-phonon system.

ACKNOWLEDGMENTS
We are grateful to Jiehang Zhang for his encouragement during the early stages of this interdisciplinary collaboration.

APPENDIX A: EXPERIMENTAL SPECIFICATIONS OF THE TRAPPED ION SYSTEM CONSIDERED FOR EXAMPLES OF THIS WORK
In order to provide explicit protocols in the examples provided in Sec. III and Appendix B, the ion-trap system that is considered is assumed to share similar features as those realized in Refs. [72][73][74]. Nonetheless, the general procedure for obtaining these protocols can be identically applied to systems containing other species of ions, and exhibiting different laser characteristics.
Consider N 171 Yb + ions confined in a radio-frequency Paul trap [84]. The "qubit" in this system has been commonly encoded in a magnetically insensitive clock state of  171 Yb + . However, for the quantum simulations of the gauge theories considered in this study, magnetically sensitive hyperfine levels |F = 0, m F = 0 and |F = 1, m F = −1 will be needed; see Fig. 10. The former (latter) level corresponds to s z = − 1 2 ( 1 2 ) component of a quasispin operator. These are split in energy by a corresponding frequency ν 0 ≡ ω 0 /2π = 12.642819 GHz + 310.8B 2 0 Hz/G 2 , where B 0 denotes an external magnetic field [113]. Highly efficient state initialization and readout are performed using a laser tuned to 369.5 nm, which strongly couples the ground 2 S1/2 and excited 2 P1/2 states.
For the Paul trap considered in the proposals of this work, ν A = 0.713 MHz and ν T = 4.1351 MHz, where ν A and ν T are the axial and transverse frequencies of the confining potential, respectively. The axial and transverse normal-mode frequencies in such a trap are tabulated in Table I Hamiltonian. Such a bias magnetic field introduces a significant error to the desired evolution. Any attempt to null out such a local magnetic field with additional sets of lasers will cause further nonzero commutations with the H I and H II Hamiltonians, which are generally non-negligible given the strength of the bias magnetic field. 14 It is therefore important to investigate solutions that eliminate the term proportional to α 0 in the native Hamiltonian in Eq. (5). One such solution relies on tuning the polarizations and detuning of the Raman beams used to produce the H III Hamiltonian such that the 14 Such a bias magnetic field term is discussed in Ref. [69].
spin-dependent force acting on the state |↑ is negative to that on the state |↓ : F ↑ = −F ↓ . This then sets α 0 = 0, which is the choice used in our proposal in Sec. II. To demonstrate this solution, we consider the example of 171 Yb + ; however, the same approach can be taken to find schemes that work for other ion traps as well.
As mentioned in Appendix A, the qubit is encoded in the magnetically sensitive |↑ ≡ |F = 0, m F = 0 and |↓ ≡ |F = 1, m F = −1 hyperfine 2 S1/2 states of 171 Yb + . Consider a set of Raman beams with frequencies ω r and ω b , detuned from 2 P1/2 manifold by . In order to produce a spin-dependent force as discussed in Sec. II, the beams have to be detuned from each other by the motional mode's frequency ω m , that is ω = ω b − ω r = ω m ; see Fig. 10. In order to find appropriate polarizations and detuning that allow a pure σ z Hamiltonian, three quantities must be calculated in this scheme: (1) the Stark shift induced by red and blue lasers in the Raman pair, (2) the spontaneous emission rate from excited states, and (3) the spin-dependent force on the qubit. Quantity (3) must be studied to deduce the conditions under which F ↑ = −F ↓ , while at the same time quantity (1) must be ensured to vanish, and quantity (2) must be minimized.
The reduced matrix element α J ||d||αJ is related to the spontaneous emission rate γ between states with J and J quantum numbers for an atom coupled to free space: where c 0 is a number that depends on the transitions. For simplicity, in the following we assume that the 2 P1/2 and the 2 P3/2 states have the same c 0 and γ .

Stark shift
In the limit where γ , the Stark shift for |m S = |↑ , |↓ is given by [115] δ Stark (m S ) = 1 4 where i is the detuning from the states that are virtually occupied, and E j is the electric-field amplitude. Using Eq. (B1), the net Stark shift is found to be As is evident, by choosing |b − | 2 + |r − | 2 = |b + | 2 + |r + | 2 , the net shift can be set to zero.

Spontaneous emission
The spontaneous emission rate can be evaluated using [115] where P m S is the probability of being in the m S ground state. Under the constraint that sets Eq. (B4) to zero, one finds that As is seen, with the choice = ( √ 2 − 1)ω F one is close to a local minimum of the spontaneous emission rate.

Spin-dependent force
Finally, the spin-dependent force can be found by considering the resonant two-photon Raman Rabi rate [115] where ϕ r and ϕ b are the phases of the red-and bluedetuned beams, respectively. With ϕ ≡ ϕ b − ϕ r = 0 and = ( √ 2 − 1)ω F , one finds that In order to satisfy the condition (↓) = − (↑) or in turn F ↑ = −F ↓ , 15 a choice for the polarization vectors iŝ Of course, these analytical solutions rely on the approximations that were made throughout these calculations, such as equal spontaneous emission rate from all the excited states considered. When precise values of the physical parameters in the system are input, the optimal values for the parameters can still be evaluated numerically using the formalism outlined. See also Ref. [116] for a similar approach in achieving the condition F ↑ = −F ↓ .

APPENDIX C: DETAILS OF THE LASER-ION EVOLUTION OPERATOR
In this Appendix, the explicit forms of the functions appeared in Eqs. (7)-(10) of the main text will be provided. The following frequency parameters are used: while the rest of the parameters/functions are already defined in Sec. II: 15 Note that the spin-dependent force is related to the Rabi frequency via F m S = k (m S ). In order to find the Rabi frequencies in the multifrequency scheme of Sec. III B, the following optimization procedure was implemented. A cost function is defined as where { (i) m } are the set of Rabi frequencies at ion i corresponding to a detuning from mode m . J targ.
i, j is the target spinspin coupling matrix, e.g., corresponding to either the nearestneighbor or the long-range interactions in the Schwinger Hamiltonian in Eqs. (19)(20)(21) for N number of ions. J i, j ({ (i) m }) in the multifrequency scheme is given in Eq. (23) for a given set of lasers. The laser indices are suppressed in the following discussion. The number of beatnote frequencies and their values are fixed (they can be chosen by running the optimization routine for select values and find the optimal values). These can also be treated as variables to be simultaneously optimized along with Rabi frequencies, but it was found that the optimization is much more robust when fixed values of beatnote frequencies were used.
The cost function is then minimized with respect to variables { (i) m } using a numerical routine, such as Mathematica's NMinimize, under the following conditions. First, one assures that the Rabi frequencies obtained do not allow the sum of the maximum magnitude of the first-order contaminating terms in Eqs. (C4)-(C6) exceed a chosen value , so the evolution remains close to the desired one. Hereñ is the number of beatnote frequencies used for a single laser. A sample of random times were picked in the α i,m (t ) function to approximate the maximum amplitude. was set to at most 30% of the maximum matrix element in J targ. in the examples shown (although this sum rarely exceeded a few percent in all cases). Second, the sum of Rabi frequencies at each ion i is set to be less than 2 MHz to conform to the current experimental limits,˜n The optimization for the case of the single-frequency scheme is a special instance of the problem described above wheñ n = 1.

APPENDIX E: ENGINEERED HAMILTONIAN OF THE SCHWINGER MODEL WITH N = 8 IONS
The multifrequency, multiamplitude scheme presented at the end of Sec. III describes the engineering of the long-range Hamiltonian of the Schwinger model in the eight fermionsite theory; see Fig. 7. The same optimization procedure can be adopted to engineer the nearest-neighbor Hamiltonians in the same theory using sets of laser beams that address transverse (for H (xx) ) and axial (for H (yy) ) normal modes of motion. The associated results, as well as the required effective magnetic field that produces H (z) , are depicted in Figs. 11-13. Associated numerical values are presented in the Supplemental Material [83]. These values are obtained from a numerical study and do not represent experimental findings.

APPENDIX F: NUMERICAL EVALUATION OF LASER-ION EVOLUTION
In order to confirm that the evolution of laser-ion systems in the scheme proposed in this work follows that of a Heisenberg spin model with a magnetic field, the exponent  7)] can be numerically evaluated for each set of laser beatnote and Rabi frequencies found. Here we assume that the ions are in their motional ground state, which can be achieved in current ion-trap experiments. The results of this evaluation are plotted, respectively, in Figs. 14 and 15 for the case of the Schwinger-model parameters with N = 4 and N = 8 that were studied in Sec. III. These figures correspond to the evolution of the first ion in the chain and the results for the rest of the ions are included in the Supplemental Material [83]. To interpret these plots, note that the quantities that are plotted are contributions to the exponent of the evolution operator as a function of time t in milliseconds (ms) and the following: