Dynamic Hamiltonian engineering of 2D rectangular lattices in a one-dimensional ion chain

Controlling the interaction graph between spins or qubits in a quantum simulator allows user-controlled tailoring of native interactions to achieve a target Hamiltonian. The flexibility of engineering long-ranged phonon-mediated spin-spin interactions in a trapped ion quantum simulator offers such a possibility. Trapped ions, a leading candidate for simulating computationally hard quantum many-body dynamics, are most readily trapped in a linear 1D chain, limiting their utility for readily simulating higher dimensional spin models. In this work, we introduce a hybrid method of analog-digital simulation for simulating 2D spin models and dynamically changing interactions to achieve a new graph using a linear 1D chain. The method relies on time domain Hamiltonian engineering through a successive application of Stark shift gradient pulses, and wherein the pulse sequence can simply be obtained from a Fourier series decomposition of the target Hamiltonian over the space of lattice couplings. We focus on engineering 2D rectangular nearest-neighbor spin lattices, demonstrating that the required control parameters scale linearly with ion number. This hybrid approach offers compelling possibilities for the use of 1D chains in the study of Hamiltonian quenches, dynamical phase transitions, and quantum transport in 2D and 3D. We discuss a possible experimental implementation of this approach using real experimental parameters.


I. INTRODUCTION
Dynamical evolution of interacting quantum manybody systems are often intractable with classical computation. Controlled studies are best done in quantum simulators [1][2][3][4] wherein the essential many-body dynamics is manifest but resides in an experimentally manageable configuration. Trapped ions [3] are among the most versatile platforms for quantum simulation, especially for simulating quantum spin systems owing to their inherent long-range interactions even when the ions are situated in a 1D topology. Even amongst the gallery of other physical implementations of quantum simulation of spin models, including ultracold bosonic and fermionic atoms in optical lattices [2,5,6], and superconducting circuits [7][8][9], the trapped ion platform offers the versatile and advantageous ability to manipulate individual spin-spin interactions, in principle, arbitrarily [10].
Long range spin-spin interactions are exceedingly simple to generate in ion trap quantum simulators, and additionally, can be controlled in their range, magnitude, and sign [11][12][13][14][15][16][17][18][19]. Leveraging phonon modes to build the inter-spin interactions makes a trapped ion system fully-connected and so inherently higher dimensional, allowing potentially the ability to probe a rich variety of physical phenomena, such as quantum transport and localization, topological insulators [20], the Haldane model [21], as well as in topological quantum computation following the Kitaev honeycomb model [22,23] and can be advantageous for quantum computing [24].
However, despite a few notable experiments and proposals [15,19,[25][26][27][28], most quantum simulations have been limited to one-dimensional (1D) chain of ions due to the constraints of radio-frequency ion traps [29]. While experimental efforts to broaden the number of ion traps with higher dimensional ion arrays are underway, significant experimental simplification is offered by leveraging existing 1D ion chains, especially considering remarkable progress, where N >100 ions have been trapped in a linear geometry [30,31], and prospects for still larger system sizes looking optimistic in the future. Additionally, existing simulators can be experimentally resource-intensive, operating on either analog [14-19, 30, 32-37] or digital [24,[38][39][40] quantum simulation protocols. Whereas an analog approach requires fine tuning of several experimental parameters for careful frequency, amplitude and phase modulation of optical beams, leading to control parameters that often scale quadratically with ion number; a digital approach requires individual beams for each ion and consequently stringent control over optical beam paths of the same number as there are ions.
In this work, we propose a hybrid method of analogdigital quantum simulation that can allow the dynamic engineering of a fully-connected 1D ion chain to, in prin-ciple, an arbitary 2D lattices (see Figure 1). When the target lattice contains certain symmetries, for instance in the case of engineering a 2D rectangular lattice, the quantum control required (O(N )) scales exceedingly favorably compared to other methods. The method relies on a repeated and stroboscopic application [41][42][43] of the full interaction HamiltoniansĤ int and -Ĥ int , and laser driven Stark shift gradients, allowing the time-domain engineering of the interaction graph. In an analogy to holography, the exact Hamiltonian engineering is efficent in the Fourier domain of couplings in the interaction graph, allowing a powerful means to engineer the target Hamiltonian while exploiting its inherent symmetries. Indeed, as we shall demonstrate, the time-sequence of Hamiltonians to be applied, can be simply read-off from a Fourier series expansion of the target Hamiltonian graph in an appropriate encoding space. Most importantly, since the engineered Hamiltonians can be dynamically modified, this opens several possibilities for studying quantum transport, dynamical phase transitions under a Hamiltonian quench [30,37,44,45], and thermalization [46] and manybody localization [47][48][49][50][51] in high dimensions. Figure 1 shows a schematic of the Hamiltonian engineering scheme. It works by removing (decoupling) interactions (forthwith "class B" interactions) that are absent in the target Hamiltonian graph, while appropriately weighting (engineering) the other interactions ("class A"), all by the global manipulation of all spins in the linear ion chain. Thus, the experimental implementation is considerably simpler than a fully digital simulation model, which requires individual two qubit gates on the ion chain. Practically, the global spin-spin interactions (±Ĥ int ) are realized by laser driven Mølmer-Sørensen couplings [11] and the single qubit phase gates (byĤ ext ) are realized by imprinting light shift (AC Stark shift) in the qubit frequency by an additional laser beam with an intensity gradient. The sign of the internal Hamiltonian (±Ĥ int ) can be flipped by changing the frequencies of global laser beams [19]. The scheme can be extended to other 2D lattice geometries, 3D lattices, and can potentially be adapted to other systems with long-range interactions and control over individual spins. Our approach therefore offers both a simplification of control parameters and favorable scaling with ion number, and offers compelling possibilities for exploiting the remarkable versatility of long-range coupled linear chain of ions for the generation of exotic engineered Hamiltonians.

II. LONG RANGE SPIN-SPIN INTERACTIONS IN TRAPPED IONS
Internal electronic or hyperfine states of trapped ions act as pristine spin-1/2 objects, with coherences extending to many minutes [52]. These long coherence times are useful for implementation of control pulse se-quences. In a typical radio-frequency Paul trap, multiple ions can be laser-cooled to a linear chain configuration in an anisotropic confinement potential. The Coulomb repulsion between ions results in collective phonon or vibrational normal modes. Off-resonant optical dipole forces exerted by laser beams or global microwave radiation with a strong spatial field gradient [53,54] can induce spin-phonon couplings, resulting in phonon-mediated spin-spin interactions [11,13,[55][56][57][58]. For example, a long range flip-flop or XY Hamiltonian, has been engineered in recent experiments [17,18], wherê S ± =Ŝ x ±iŜ y are the raising and lowering spin operators. The interactions in Eq. 1 can be tuned according to a power law, Here, 0 < α < 3 sets the range of interactions [12,58] and J 0 is the nearest neighbor coupling strength. Without loss of generality, we assume J 0 > 0. More generally, with full control over individual spin-phonon couplings, it has been shown that J ij can be arbitrarily programmed [10]. However, controlling all the spin-phonon couplings is a challenging experimental task and may not be scalable with current technologies. As a special case of Eq. (2), the range of interactions can be quasi-infinite range, α ≈ 0 when the field driving coherent spin-phonon coupling is tuned close to the center of mass (COM) phonon mode. Further, by changing frequency of laser beams, the sign of interactions [13,19] can be flipped (see Appendix A). For a fully connected system of N spins, there are couplings. To engineer the target interaction graph, we categorize these couplings into two classes, which we name class A and class B (see Fig. 1a). The couplings in class A are present in the target graph, while those in class B need to be removed. Further, the couplings in class A may need additional scaling depending on the target interaction graph. For example, in Fig. 1c, the vertical nearest neighbor bonds J ′ V in the target square lattice are obtained from the third neighbors (J i,i+3 ) in the original 1D chain, while the horizontal nearest neighbor bonds J ′ H are obtained from the nearest neighbors (J i,i+1 except J 34 ) in the 1D chain. In a square lattice, we require J ′ H = J ′ V and hence our engineering protocol must reduce the strength of the nearest neighbor couplings (except J 34 ) in the 1D chain to match the third neighbor couplings if α > 0 in Eq. 2. This is because we cannot enhance the strength of a coupling or engineer one that does not already exist in the original network. The protocol for Hamiltonian engineering is described in detail in section III. The interaction network can be modified by subjecting the system to a periodic sequence of free evolution under the native Hamiltonian ±Ĥint and single qubit (phase) gatesĤext acting simultaneously on all qubits. c) The average Hamiltonian of the system,Ĥ eff resembles the target HamiltonianĤTarget, here a 2 × 3 rectangular lattice, at discrete time steps, t = nTcyc (n = 1, 2, · · · ). The horizontal and the vertical couplings of the target Hamiltonian, J ′ H and J ′ V respectively, can be dynamically changed in a simulation. A square lattice is obtained when

III. METHOD
A. Cancellation of couplings by reversing unitary time evolution One way to effectively cancel spin-spin interactions is subjecting the system to periodic time evolutions alternating underĤ int and −Ĥ int for equal durations. The unitary time evolution underĤ int is reversed under −Ĥ int , returning to the initial state. In Fig. 2a, we illustrate this for N = 2 spins. The spins, when initialized in the product state |↑↓ will undergo coherent oscillation between |↑↓ and |↓↑ under the flip-flop Hamiltonian, Eq. (1). Here, |↑ and |↓ are the eigenstates ofŜ z . If the Hamiltonian is instead switched alternately between H int and −Ĥ int the spins go back to their initial state after each cycle of duration T cyc . Thus, when observed discretely at time t = nT cyc (n = 0, 1, 2, · · · ), the spins appear to be non-interacting. However, we must need additional ingredients in this protocol as we want to protect interactions in class A. Our protocol must distinguish between couplings that we want to cancel by reversing unitary time evolution and couplings that we want to rescale. This is done by imprinting separate 'phase tags' between these classes of couplings by a spatial field gradient applied simultaneously on all the spins [43]. The Hamiltonian from this external field is, where {ω i } depends on the nature of the field gradient. H ext can be engineered in experiments, e.g., by a laser beam with spatially inhomogeneous intensity pattern imprinting AC Stark shifts on the spins. If the external Hamiltonian, Eq. (3), is applied for a duration τ , the Hamiltonian in Eq. (1) is transformed tỗ with n being an integer. Thus, the interactions in class A evolve underĤ int = e −i(2n−1)πĤ int = −Ĥ int , whereas interactions in class B evolve underĤ int = e −i2nπĤ int = H int . The external field gradient and switching between H int and −Ĥ int can be combined in the time evolution cycle to preserve interactions in class A while canceling interactions in class B. Fig. 2b demonstrates the basic principle for N = 2 spins. One cycle in the time evolution consists of two subsequent blocks of pulses. The first block consists ofĤ ext for time τ followed byĤ int for time t, and the second block consists ofĤ ext for time τ followed by −Ĥ int for time t. Setting φ ij = (2n − 1)π preserves the interaction between the (i, j) pair of spins (class A), while φ ij = 2nπ cancels the interaction (class B).

B. Scaling of couplings by a multi-pulse sequence
To rescale the couplings {J ij }, the simple two pulse sequence can be replaced by a multi-pulse sequence acting as a filtering function, as schematically shown in Fig. 2c. Each block now consists of alternate pulses ofĤ ext for time τ i and ±Ĥ int for time t i . The total phase accumulated by the (i, j) spin-pair fromĤ ext must satisfy φ ij tot = ω ij τ tot = (2n − 1)π for it to not cancel under the reversal of unitary evolution by −Ĥ int . Here, τ tot = l k=1 τ k , l being the number of pulses ofĤ ext in each block. The unitary evolution of the system at t = T , the end of a multi-pulse block, is given bŷ (4) whereĤ s = ±Ĥ int . This equation can be rewritten in the following form whereĤ s k = e iĤext k j=1 τjĤ s e −iĤext k j=1 τj is the internal Hamiltonian in the rotating frame defined by successive frame transformation induced by the external gradient field. Eq. 5 is derived usingÎ = e −iĤextτj e iĤextτj . Upon applying the average Hamiltonian theory [41,59,60], an effective Hamiltonian for the full block can be defined aŝ so thatÛ Here we have assumed that l j=1 t j ≈ T , which is a valid approximation if τ tot ≪ T and can be realized with the external gradient field much stronger than the spinspin interactions. Eq. (7) implies that the the multipulse block can be thought of as a single pulse ofĤ ext applied for a duration of τ tot followed by a pulse ofĤ eff for a duration of T . Note that, in order for the average Hamiltonian formalism (Eq. (6)) to be valid, J 0 T ≪ 1.
Eq. (6) can be written in a more explicit form for the coupling J ′ ij ofĤ eff in terms of the J ij ofĤ int .
the total phase accumulated in the multi-pulse block. The realvalued scaling parameter β ij is explicitly given by, and can be engineered by choosing {t k , τ k } (for k = 1, 2, · · · , l) for a chosen {ω i }. If the multi-pulse block is repeated while all ±Ĥ int are replaced with ∓Ĥ int , the effective interactions at the end of the cycle, t = T cyc = 2T , becomes [43] and the couplings are rescaled. While, if φ ij tot = 2nπ, the couplings vanish automatically regardless of β ij , Eqs. (9)- (11) indicate that we only need to consider those pairs with φ ij tot = (2n−1)π in designing the scaling parameter β ij . As we will show in Section III D, the pulse sequence to engineer the target scaling factors (Eq. (8)) can be constructed from a Fourier series expansion of the target Hamiltonian in the space of the interactions.

C. Labeling and field gradient scheme
To map the interactions J ij in the 1D spin-chain onto the target 2D rectangular lattice, we categorize them into classes N d according to the distance d = |j − i| between the spins, with N 1 denoting the nearest neighbor couplings, N 2 denoting the next nearest-neighbor couplings and so on. Here, we ignore inhomogeneities due to the finite size effect, which we discuss in Section IV C. As seen in Fig. 3a, the N 1 and N m couplings form the horizontal and vertical bonds of an m ′ ×m rectangular lattice. That is class A = {N 1 , N m }, with the exception of couplings J km,km+1 , k = 1, ..., m ′ − 1, that form a toroidal linkage between the edges of the rectangular lattice and must be excluded from class A. Thus, we need an external gradient field {ω i } that assigns, is reversed and the spins return to the initial state at t = nTcyc (n = 1, 2, · · · ). Thus, the interactions are effectively canceled when the system is viewed at discrete time steps, t = nTcyc. b) Evolution 1 and Evolution 2 are controllably reproduced by the same pulse sequence, if we introduce an external gradient field (Eq. 3). In this scheme, one cycle in the time evolution consists of two subsequent blocks. The first block consists ofĤext for time τ followed byĤint for time t and the second block consists ofĤext for time τ followed by −Ĥint for time t. TheĤext pulse assigns a phase tag φ = ω0τ to the coupling between the ions, where ω0 = ω2 − ω1. If the assigned phase φ = (2n − 1)π, n being an integer number, Evolution 1 is reproduced and the interaction is retained. If φ = 2nπ, Evolution 2 is reproduced and the interaction is effectively canceled. c) To rescale the interactions, the two-pulse sequence is replaced by a multi-pulse sequence. Each block now consists of alternate blocks ofĤext and ±Ĥint for different durations. The duration of the multi-pulse sequence cycle Tcyc is 2T , where T is the duration of each block. In alternate blocks ±Ĥint are switched with ∓Ĥint [42].
The integer n does not have to be the same for all of the couplings. However, the problem of designing the filtering function simplifies with a small number of different n's, 2) distinct phase tags to J km,km+1 from other couplings in class A, 3) φ ij tot = 2nπ to as many couplings as possible in class B, 4) φ ij tot = (2n−1)π to couplings in class B for which it is not possible to assign φ ij tot = 2nπ after satisfying the constraints for couplings in class A. These couplings have to be rescaled to zero (β ij = 0) by the Fourier filtering function.
A semi-linear field gradient as shown in Fig. 3b satisfies the conditions above. We choose the external field profile to increase linearly with a constant slope ω 0 with added jumps of 2ω 0 (for even m) or 3ω 0 (for odd m) between the km th and (km + 1) th spins to address the toroidal linkages J km,km+1 . Hence, ω i,i+1 = ω 0 for N 1 except J km,km+1 , ω i,i+m = (m + 2)ω 0 for N m when m is odd, ω i,i+m = (m + 1)ω 0 for N m when m is even.  3)) provides the necessary conditions, described in the text, for efficiently canceling and scaling interactions to achieve the target graph. In this profile, ωi increases linearly with a constant slope ω0 with some added jumps of +2ω0 (for even m) and +3ω0 (for odd m). These jumps are designed between the km th and (km + 1) th ions, such that the torroidal couplings J km,km+1 can be canceled (k = 1, 2, . . . , m ′ − 1). See Fig. 9 in Appendix B for an example.
So, ω 0 τ tot = π satisfies the conditions above for all couplings within class A and some of the couplings in class B.  (12) is fit to the Fourier series target scaling parameters, βij in Eq. (10). The field gradient in Eq. (3) creates a phase of φ = φij = ωij τtot for all the couplings Jij that need to be rescaled by the same factor βij . For this target lattice and the field gradient of Fig. 3(b), there are four relevant phases, π, 3π, 5π, and 7π given the constraints set by Eq. (10) and Eq. (11). The points with negative phases are included to satisfy the symmetry condition of Eq. (12). All the nearest neighbors, N1, except J34, acquire a phase of φ = π. They are scaled by a factor of 1/3 α = 0.8 with respect to N3 to attain a square lattice geometry, if the interactions in the 1D chain is decaying with an exponent α = 0.2 in Eq. (2). The torroidal coupling φ = φ34 = 3π has to be rescaled to zero to be removed from the target graph. The Fourier coefficients directly give the duration of the pulses {tj } within a time cycle of the simulation.

D. Fourier filtering of interactions
The interactions that were not automatically be canceled by our chosen field gradient can be rescaled to their target value by designing a Fourier filter. The N 1 couplings in class A have to be rescaled to match the N m couplings for α = 0 in Eq. (2). In addition, the couplings in class B for which φ ij tot = (2n− 1)π have to be rescaled to zero. The scaling is performed through a Fourier filter function F (φ) that produces the desired scaling parameter β ij in Eq. (10). The filter function F (φ) is defined as, The filter function F (φ) should satisfy F (φ ij tot ) = β ij for all couplings with a phase φ ij tot = (2n − 1)π.
By comparing Eq. (8) with Eq. (12), we construct a multi-pulse sequence for implementing the Fourier filtering function F (φ), with the following features: 1. The number of single qubit phase gates (byĤ ext ) l in each block is l = 2i + 1, where the number of Fourier terms in Eq. (12) is i + 1.
2. The pulse sequence in each block is anti-symmetric about the centralĤ ext pulse. That is t j = t l−j and ±Ĥ int j = ∓Ĥ int l−j for j = 1, ..., l − 1. This leads to the cancellation of all even order correction terms to the average Hamiltonian in Eq. (6).
3. The time intervals {t j } are proportional to the coefficients in the Fourier filter function : t j = T |a j |/2 for j = 1, . . . , l − 1 and t l = T |a 0 |, with the constraint that l j=0 |a j | = 1. This constraint can be relaxed for an efficient search for the Fourier coefficients at the expense of reducing all the couplings in the target lattice by a global rescaling factor. Numerically, we find that β i,i+m = 0.7 allows us to find efficient solutions for up to N = 100. 4. A negative coefficient (a j < 0) in Eq. (12) is implemented by anĤ ext pulse followed by −Ĥ int .

5.
We choose the phase gates to be of equal duration τ , except the central phase gate in each block which has to have a duration of τ ′ = τ tot − (l − 1)τ . τ can be read-off from W = τ /τ tot .
In Fig. 4, we show the Fourier filter function fit for engineering a 2 × 3 square lattice when α = 0.2 in Eq.
(2). The nearest neighbor couplings N 1 , except the toroidal linkage (J 34 ) have a phase of φ ij tot = π. The N 3 couplings accumulate a phase of 5π. Hence we require F (π)/F (5π) = 1 3 α = 0.80 such that N 1 couplings (except J 34 ) are equal to N 3 . The toroidal linkage J 34 (φ ij tot = 3π) and N 5 couplings (φ ij tot = 7π) are scaled to zero. We have introduced a global rescaling factor to all the target couplings, by setting F (5π) = 0.7. The global rescaling of all the target couplings ensures an efficient Fourier fit with minimum number of parameters for at least up to N = 100 spins.
In Table II, the Fourier fit parameters for engineering a 2 × 3 (listed in Row 1) and a 3 × 3 square lattices (listed in Row 2) are given for α = 0.2.

A. Numerical simulations
Applying the tools described above, we successfully reproduce the spin-spin interaction graph of the target square lattice at evolution times t = nT cyc , n = 1, 2, · · · , starting with the initial long range couplings with an exponent α = 0.2. Figs. 5a and 5b illustrate the results for a 2 × 3 and a 3 × 3 square lattices, respectively. The engineered interaction matrix formed by the couplings {J ij } matches the target interaction matrix of the 2D square lattice with an RMS error of < 0.1%. Here we define the RMS error as ij (Target) denote the couplings in the ideal lattice. In Fig. 5, we also compare spin dynamics under the engineered lattice (red dots) with that of the ideal target (green curve), and find excellent agreement. Here the systems are initially prepared in |ψ 0 =|↑ 1 ↓ 2 ↓ 3 ↓ 4 ↓ 5 ↓ 6 for N = 6 and |ψ 0 =|↑ 1 ↓ 2 ↓ 3 ↓ 4 ↓ 5 ↓ 6 ↓ 7 ↓ 8 ↓ 9 for N = 9 at t = 0, when the pulse sequence is turned on. The probability of the system being in |ψ 0 (Figs. 5a(iii) and 5b(iii)) follows the expected dynamics of the ideal 2D square lattice. These numerical simulations were performed using the time dependent master equation solver based on QUTIP python package [61,62].
The near-perfect matching of the target and engineered spin dynamics in Fig. 5 indicates small intrinsic errors. These errors arise from the Fourier fit in Eq. (12) and numerical rounding error in the time interval values {t i , τ i }. However, additional error may creep into an experimental realization due to imperfect single qubit gates. In our numerical simulation, we find that RMS error between the target and engineered interaction matrices scales linearly with single qubit phase error, with 1% error in single qubit phase (in each pulse ofĤ ext ) contributing to approximately 1.2% error in J ′ ij .
A crucial feature of the protocol presented here is the ability to dynamically change the Hamiltonian within the same symmetry class, by changing the Fourier coefficients in Eq. (12). This enables simulation of many-body dynamics, such as quantum quench experiments that are hard to simulate numerically. As an example, we show a quench from two decoupled chains of 3 spins each into the 2 × 3 square lattice plaquette in Fig. 6. To simulate the decoupled chains, N 3 couplings are set to zero (see Fig. 4) in estimating the Fourier filtering function and hence the multi-pulse sequence. The spin-spin correlations between the previously uncoupled chains start to build up after the quench. We show the engineered dynamics (red dots) of the two point correlation functions C 12 (t) between spins 1 and 2 and C 14 (t) between spins 1  vs (i, j). The couplings are normalized to J0 (ii) The engineered interaction graph closely resembles the target interaction graph of the square lattices with an RMS error (defined in the text) of < 0.1%. The engineered couplings J ′′ ij = J ′ ij /max{J ′ ij } are shown vs (i, j) in the color plot. (iii) The evolution of the engineered lattice (red dots) is compared with the evolution of the target lattice (green curve). The system is initially prepared in (a) |ψ0 =|↑1↓2↓3↓4↓5↓6 and (b) |ψ0 =|↑1↓2↓3↓4↓5↓6↓7↓8↓9 state. The probability of measuring the state of the system in its initial state is shown over time. The evolution of the system under multi-pulse sequence replicates the evolution of the target lattice. and 4. They follow closely to the ideal target dynamics (green line).

B. Proposal for experimental implementation
The experimental implementation of the multi-pulse scheme can be achieved in a trapped ion system such as 171 Yb + trapped in a radio-frequency (Paul) trap. When the confining potential is sufficiently anisotropic, lasercooled ions will form a linear chain [63]. The hyperfine states |0 = | 2 S 1/2 , F = 0, m F = 0 and |1 = | 2 S 1/2 , F = 1, m F = 0 form the qubit states for this ion [64].
The flip-flop HamiltonianĤ int in Eq. (1) can be simulated [17,18] by global Raman laser beat-notes using the Mølmer-Sørensen scheme [11] for generating phononmediated spin-spin interactions. When the Mølmer-Sørensen detuning is tuned close to center of mass phonon mode, a long range interaction in the form of Eq. (2) can be obtained. The global sign ofĤ int can be flipped by changing the Mølmer-Sørensen detuning, with additional detunings improving the accuracy (see Appendix A for details).
The field gradient inĤ ext can be implemented with laser beams imprinting AC Stark shifts, and by spatially modulating the laser intensity using a Spatial Light Modulator (SLM) or an Acousto Optic Deflector (AOD). Another potentially easier experimental implementation will be to combine a global tightly focused laser beam with additional relatively low power beams created by an AOD. The global beam propagating along the axis of the ion chain can be focused before hitting the ions, such that its intensity varies linearly on the ion chain. The jumps in the gradient field can be added by beams created by an SLM or AOD and shining on the ions from the transverse direction. A 100 mW laser beam propagating along the ion chain, detuned from the 171 Yb + 2 S 1/2 − 2 P 1/2 resonance by 10 5 natural line-widths, and  (12) can be updated at each time cycle to realize a dynamically changing target Hamiltonian. For example, here we perform a quench from two decoupled chains of 3 spins each into a 2 × 3 square lattice. a) the correlation measure between Spin 1 and 4 defined by the correlation measure between Spin 1 and 2 defined by C12(t) = S 1 z S 2 z − S 1 z S 2 z . The spin-spin correlations between the previously uncoupled chains start to build up after the quench, indicated by the dashed line. The engineered dynamics follow closely to the ideal target dynamics.
focused to approx. 2 microns will create a two-photon differential AC Stark shift gradient (ω 0 ) of approximately 1 MHz. Thus, the total time (τ tot ) that the Stark shifting beam is shining on the ions in a time cycle can be limited to a few microseconds, minimizing spontaneous emission errors.

C. Discussions
The Hamiltonian engineering protocol presented here makes efficient use of Fourier filtering to achieve a desired spin-spin interaction topology. For an arbitrary target lattice, O(N 2 ) Fourier coefficients, hence number of pulses, will be needed. However, in presence of common symmetries between the target lattice and the external field gradient, the number of pulses are reduced drastically. For example, the rectangular lattices presented here need O(N ) pulses (Fig. 7), as our chosen field gradient (Fig. 3b) creates the same phase tag for the class N d , except a small number of (O( √ N )) torroidal linkages. The number of Fourier coefficients are further reduced, approximately by a factor of 2, by using ±Ĥ int instead ofĤ int only. This is because all of the interactions for which φ ij = 2nπ in a time cycle are canceled automatically (Eq. (11)) and hence are not required to be included in estimating the Fourier filtering function.
The engineered interactions in the target 2D lattice will become weaker for a given J 0 in the original 1D chain, as the system size increases. This is because of the following reasons. The average Hamiltonian theory employed here works when J 0 T cyc ≪ 1. Due to the linear scaling of the number of pulses in a cycle with N (Fig. 7), T cyc is expected to scale linearly, and hence the initial coupling J 0 has to be reduced linearly with increasing N . We may have to further scale some couplings down to match the target interaction graph, as the interactions are decaying with distance in the original 1D lattice. For example, the N 1 couplings (except the torroidal linkages) have to be scaled down by a factor of 1/m α compared to the N m couplings for an m ′ × m square lattice, as demonstrated in Fig. 3 and Fig. 4. For the results presented here with N = 6 and N = 9 ions, we have chosen α = 0.2, which is experimentally realizable in current systems, and provides sufficiently strong target interactions. However, for α > 0, the couplings will scale down with increasing N , and hence longer simulation times are necessary as the system size increases. We estimate a target coupling of 2π ×300 Hz in a 2×3 square lattice. Since the separation between the vibrational normal modes in the ion chain will decrease with increasing N , the coupling strength J 0 will have to scale down accordingly in order to avoid direct excitation of phonons that limit the validity of a spin-only Hamiltonian. While trapped ion qubits have long single qubit coherence time [14], scaling the simulation to a large number of spins where classical computation of dynamics may be intractable will require isolating experimental noise sources, such as intensity fluctuations of the global Mølmer-Sørensen and Stark shifting laser beams, and drifts in the collective phonon mode frequencies.
In a small chain of ions, the couplings are inhomogeneous. The errors due to the inhomogeneity can be mitigated by using an anharmonic trapping potential for the ions [65] and spatially modulating the global Mølmer-Sørensen laser beams to increase the homogeneity of the couplings. The errors can also be reduced at the expense of increasing the number of pulses within a cycle and using a field gradient that breaks the symmetry between interactions belonging to a class N d .
As the number of ions N increases, the spacing between the vibrational modes decreases. This will make it harder to engineer −Ĥ int with a single global Mølmer-Sørensen beam that is detuned close to the center of mass mode. The accuracy of engineering −Ĥ int is enhanced by introducing additional global beams with Mølmer-Sørensen detuning near the neighboring phonon modes (see Appendix A). Engineering −Ĥ int for very large N will require either a large number of global Mølmer-Sørensen beams to cancel the effect of multiple modes, or reducing the overall intensity of laser beams resulting in a reduc- Engineering all the interactions via Fourier filtering alone (by usingĤ int only instead of ±Ĥ int ) may become experimentally preferable at the expense of a longer pulse sequence (and longer T cyc ).

V. CONCLUSIONS AND OUTLOOK
In this work, we proposed a hybrid analog-digital simulation protocol that leverages the long range spin-spin interactions in a linear chain of trapped ions to engineer a 2D lattice topology. This protocol efficiently engineers certain target interaction topology, such as 2D square lattices, by Fourier engineering analogous to holography. Our work opens up opportunities to identify classes of interaction graphs that can be simulated efficiently.
The scheme needs only global beams creating spin-spin interactions and a single laser beam with intensity gradient. Unlike a full digital quantum simulator, this scheme does not require individual precise two qubit gates. An attractive feature of this protocol is the dynamical engineering of Hamiltonians, as we have demonstrated by a quench experiment. The dynamical engineering enables investigation of a range of many-body phenomena of active research interest, such as quantum quench and transport and dynamical phase transitions. The protocol is feasible with a moderately large system of tens of spins in existing trapped ion experiments.
is canceled. This can be accomplished by applying a second pair of Raman beams with µ 2 = ω tilt − δ 2 , reddetuned from the tilt mode and chosen so that δ 2 < δ ′ 1 . This brings in another contribution from the tilt mode that scales as − (η tilt Ω 2 ) 2 /δ 2 . One can optimize for the lasers parameters, i.e., Rabi frequencies and detunings, so that the coupling profiles forĤ int and −Ĥ int match as closely as possible. That is the difference between the corresponding coupling profiles defined by In Fig. 8 we show a set of experimental parameters for a chain of 6 ions when α = 0.2. To implement +Ĥ int , a Raman beat-note with a global sideband Rabi frequency η COM Ω 1 = 2π × 18 kHz and δ 1 = 2π × 55 kHz is bluedetuned from the COM mode (Fig. 8a). This results in the coupling profile J ij [δ 1 ] with the nearest neighbor coupling constant J 0 ≈ 2π × 0.520 kHz. To implement −Ĥ int , two Raman beat-notes with δ ′ 1 = 2π × 45 kHz and δ 2 = 2π × 12.3 kHz are red-detuned from the COM and tilt modes (Fig. 8a). The sideband Rabi frequencies are taken to be η COM Ω ′ 1 = 2π × 15 kHz and η tilt Ω 2 = 2π × 4.1 kHz, respectively, where η COM Ω ′ 1 is associated with µ ′ 1 . These parameters result in the coupling profile J ij [δ ′ 1 , δ 2 ]. The bar chart corresponding to Fig. 8b. This bar chart indicates up to 1.9% error when switching be-tweenĤ int to −Ĥ int during the experiment.
Appendix B: LABELING/FIELD GRADIENT SCHEME FOR 9 IONS Fig. 9a illustrates the labeling scheme for engineering a 3 × 3 square lattice from a 1D chain of 9 ions. The class A interactions in the ion chain network consists of N 1 and N 3 interactions. The {J 34 , J 56 ∈ N 1 } should be excluded as they form toroidal linkages at the horizontal edges of the square lattice and should be set to zero to engineer the target geometry. Fig. 9b illustrates the external field gradient scheme proposed for engineering a 3 × 3 square lattice. The external field profile increases linearly across the ion chain with two jumps of +3ω 0 to assign a distinct tag to J 34 and J 56 interactions. Similar to the 6-ion chain, one can adjust ω 0 τ tot = π so that φ A ij tot = (2n − 1)π and φ B ij tot = 2nπ except for N 5 , N 7 ∈ class B.  9. The labeling/field gradient scheme for engineering a 3×3 square lattice. a) The ions in the 1D chain can be labeled in a way so that only N1 and N3 couplings are required for engineering the square lattice. The N1 couplings form the horizontal bonds of the square lattice while N3 couplings form the vertical bonds. b) The proposed field gradient scheme for engineering a 3 × 3 square lattice from a 1D chain of 9 ions. The external field profile corresponding to ωi increases linearly along the chain of ions with two jumps of +3ω0 to assign distinct phase tags to J34 and J56 forming the toroidal linkage between the edges of the lattice.