We propose a protocol for time-domain engineering of the interaction graph between trapped ion spins in a quantum simulator. The method relies on a repeated and stroboscopic application38,39,40 of the full interaction Hamiltonians \(\hat{H}_{{\mathrm{int}}}\) and \(-\hat H_{{\mathrm{int}}}\), and laser driven Stark shift gradients. A quantum simulator can be tuned such that the inter-spin interactions in a one-dimensional chain of trapped ions has the form,

$$\hat H_{{\mathrm{int}}} = \mathop {\sum}\limits_{i < j} {J_{ij}} \left( {\hat S_i^ + \hat S_j^ - + \hat S_i^ - \hat S_j^ + } \right),$$ (1)

where \(\hat S_i^ \pm = \hat S_{x_i} \pm i\hat S_{y_i}\) are the raising and lowering spin operators acting on spin i. The long-range couplings are characterized by

$$J_{ij} \approx \frac{{J_0}}{{\mid i - j\mid ^\alpha }},$$ (2)

where 0 < α < 3 sets the range of interactions.7,41 Here, J 0 is the nearest-neighbor coupling strength. Without the loss of generality, we assume J 0 > 0. Our goal is to modify the interaction profile in Eq. (2) such that the 1D spin-chain can be mapped onto rectangular lattices. We achieve this with the help of an external gradient field, represented by the Hamiltonian,

$$\hat H_{{\mathrm{ext}}} = \mathop {\sum}\limits_{i = 1}^N {\omega _i} \hat S_{z_i},$$ (3)

where \(\hat S_z\) is the z component of the spin-1/2 operator. If the external Hamiltonian, \(\hat H_{{\mathrm{ext}}}\), is applied for a duration τ, the Hamiltonian in Eq. (1) is transformed to \(\hat{\tilde H}_{{\mathrm{int}}} = e^{i\hat H_{{\mathrm{ext}}}\tau }\hat H_{{\mathrm{int}}}e^{ - i\hat H_{{\mathrm{ext}}}\tau } = \mathop {\sum}

olimits_{i < j} {J_{ij}} \hat S_i^ + \hat S_j^ - e^{i\omega _{ij}\tau } + h.c.\), where ω ij = ω i − ω j (note \(\omega _{ij} \gg J_{ij}\)). Thus the phase tags ϕ ij = ω ij τ appear in \(\hat{\tilde H}_{{\mathrm{int}}}\). These phase tags give us a handle to modify the couplings J ij . 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 single and two qubit gates on the ion chain.

Our scheme consists of a periodic application (with a cycle time of T cyc ) of a multi-pulse sequence on all the qubits simultaneously, as shown in Fig. 1b. Each cycle consists of two multi-pulse blocks with alternating pulses of \(\pm \hat H_{{\mathrm{int}}}\) (with duration {t k }) and \(\hat H_{{\mathrm{ext}}}\) (with duration {τ k }). \(\pm \hat H_{{\mathrm{int}}}\) in the first block are substituted by \(\mp \hat H_{{\mathrm{int}}}\) in the second block. As explained in the Supplementary Materials, the average Hamiltonian at the end of each cycle retains the form of Eq. (1), with modified couplings \(J_{ij}^\prime\) that depend on \(\{ t_k,\tau _k,\phi _{ij}^{(k)}\}\), where \(\phi _{ij}^{(k)} = \omega _{ij}\tau _k\). \(J_{ij}^\prime\) is related to J ij by,

$$J_{ij}^\prime (T_{{\mathrm{cyc}}}) = \frac{{J_{ij}}}{2}[\beta _{ij}(1 - e^{i\phi _{ij}^{({\mathrm{tot}})}})],$$ (4)

where

$$\beta _{ij} = \frac{2}{{T_{{\mathrm{cyc}}}}}\left[ { \pm t_1e^{i\phi _{ij}^{(1)}} \pm t_2e^{i(\phi _{ij}^{(1)} + \phi _{ij}^{(2)})} \cdots \pm t_le^{i\phi _{ij}^{({\mathrm{tot}})}}} \right],$$ (5)

with \(\phi _{ij}^{({\mathrm{tot}})} = \mathop {\sum}

olimits_{k = 1}^l {\phi _{ij}^{(k)}}\). Our choice of the pulse sequence in each block results in real-valued β ij . The total phase accumulated in one block, \(\phi _{ij}^{({\mathrm{tot}})}\) has to be an integer multiple of π for \(J_{ij}^\prime\) to be real-valued (Eq. (4)). For the most efficient use of this protocol, we use a labeling scheme for mapping the ions in the 1D chain to the target lattice and the external field gradient {ω ij } such that \(\phi _{ij}^{({\mathrm{tot}})} = 2n\pi\) (n = 0, 1, 2, ⋯), and hence \(J_{ij}^\prime = 0\) for as many class B couplings as possible. This is achieved by choosing a semi-linear gradient field, which can be engineered in experiments (see Methods section). The remaining couplings in class B are set to zero, versus the class A couplings which are scaled to their target values by appropriately choosing \(\{ t_k,\tau _k = \phi _{ij}^{(k)}/\omega _{ij}\}\) in Eq. (5). In a target square lattice, we require equal horizontal and vertical nearest-neighbor bonds, \(J_H^\prime = J_V^\prime\) as seen in Fig. 1a.

In an analogy to holography, the target Hamiltonian can be engineered by a Fourier expansion of the target couplings in the domain of phases imparted by the external field over the duration T = T cyc /2 of a multi-pulse block. The Fourier ‘filtering’ function

$$F(\phi ) = a_0 + \mathop {\sum}\limits_{i^{\prime} = 1}^i {a_{i^{\prime} }} \cos (i^{\prime} W\phi )$$ (6)

is chosen such that \(F(\phi _{ij}^{({\mathrm{tot}})}) = \beta _{ij}\), where \(\phi _{ij}^{({\mathrm{tot}})} = (2n - 1)\pi\), and W is a fitting parameter. The coefficients {t k } in Eq. (5) are simply obtained from the Fourier coefficients {a i }, as explained in Methods. The Fourier engineering allows a powerful means to achieve the target Hamiltonian efficiently while exploiting its inherent symmetries. Most importantly, by dynamically modifying the Fourier filtering function, the engineered Hamiltonians can be dynamically modified. This opens several possibilities for studying quantum transport, dynamical phase transitions under a Hamiltonian quench,24,30,42,43 and thermalization44 and many-body localization45,46,47,48,49 in high dimensions.

Practically, the global spin–spin interactions (\(\pm \hat H_{{\mathrm{int}}}\)) are realized by laser driven Mølmer–Sørensen couplings6 and the single qubit phase gates (by \(\hat H_{{\mathrm{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 (\(\pm \hat H_{{\mathrm{int}}}\)) can be flipped by changing the frequencies of global laser beams.50 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 a favorable linear 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.

Numerical simulations

We successfully reproduce the spin–spin interaction graph of target 2D square lattices at discrete evolution times t = nT cyc (n = 1, 2, ⋯). We demonstrate control over dynamically changing the vertical and horizontal nearest neighbour bonds (\(J_{\mathrm{V}}^\prime\) and \(J_{\mathrm{H}}^\prime\) in Fig. 1a).

Figure 2a, b illustrates the results for a 2 × 3 and a 3 × 3 square lattice, respectively. These results are obtained using the labeling and field gradient schemes presented in Fig. 5 in Methods section. The semi-linear field gradient in Fig. 5 assigns \(\phi _{ij}^{({\mathrm{tot}})} = 2n\pi\) to most coupling in Class B and \(\phi _{ij}^{({\mathrm{tot}})} = (2n - 1)\pi\) to all Class A couplings. Decoupling interactions in Class B with an \(\phi _{ij}^{({\mathrm{tot}})} = (2n - 1)\pi\) and rescaling Class A interactions (e.g., \(J_{\mathrm{H}}^\prime = J_{\mathrm{V}}^\prime\) in Fig. 1a) are accomplished through the Fourier filtering functions F(ϕ) specified by coefficients {a i } (see Eq. (6) and Table 1). The first and second row of the table correspond, respectively, to {a i } for engineering a 2 × 3 and a 3 × 3 square lattice. The target lattices are engineered through application of multi-pulse sequences constructed from the given Fourier coefficients taking the steps discussed in Methods section. The Fourier expansion coefficients for engineering an m × m square lattice (\(m = \sqrt N\)) up to m = 7 are given in Supplementary Materials. In Fig. 3, we show that the number of pulses in each cycle scales linearly with the number of ions.

Fig. 2 Numerical results for engineering a a 2 × 3 square lattice and b 3 × 3 square lattice. (i) The interaction graphs in the 1D ion chains with α = 0.2 are shown, along with a 2D color plot of the couplings J ij vs (i, j) in Eq. (2). The couplings are normalized to J 0 . (ii) The engineered interaction graph closely resembles the target interaction graph of the square lattices with an RMS error of <0.1%. The engineered couplings \(J_{ij}^{\prime\prime} = J_{ij}^\prime /{\mathrm{max}}\{ J_{ij}^\prime \}\) 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 Full size image

Table 1 Fourier expansion coefficients for a 2 × 3 and a 3 × 3 square lattice when α = 0.2 (see Eq. (2)) Full size table

Fig. 3 Scaling of number of pulses with system size. The number of pulses (includes \(\hat H_{{\mathrm{ext}}}\) and \(\pm \hat H_{{\mathrm{int}}}\)) in a time cycle, T cyc required for engineering m × m square lattices scales linearly with the number of ions, N = m2 in a 1D chain Full size image

On the target square lattices, spins interact with strength \(J_{ij}^\prime\), which is non-zero only when the ion-pairs (i, j) are nearest neighbors in the target lattice. Since the interactions are decaying with distance in the original 1D chain, for α > 0 in Eq. (2), \(J_{{\mathrm{ij}}}^\prime\) can be at most J 0 /mα for an m′ × m square lattice.

The engineered interaction matrix formed by the couplings \(\{ J_{ij}^\prime \}\) matches the target interaction matrix of the 2D square lattice with an RMS error of <0.1%. Here we define the RMS error as \(\sqrt {\mathop {\sum}

olimits_{ij} {(J_{ij}^\prime - J_{ij}^\prime ({\mathrm{Target}}))^2} } /\mathop {\sum}

olimits_{ij} {|J_{ij}^\prime ({\mathrm{Target}})|}\). In Fig. 2, 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. Here, |↑ i 〉 and |↓ i 〉 are the eigenstates of \(\hat S_z\) for the ith spin. The probability of the system being in |ψ 0 〉 (Fig. 2a(iii), b(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 the QuTip python package.51,52

The near-perfect matching of the target and engineered spin dynamics in Fig. 2 indicates small intrinsic errors in the dynamical Hamiltonian engineering protocol. However, additional errors 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 \(\hat H_{{\mathrm{ext}}}\)) contributing to ~1.2% error in \(J_{ij}^\prime\).

A crucial feature of the protocol presented here is the ability to dynamically change the Hamiltonian within the same symmetry class, by changing the pulse sequence obtained from the Fourier decomposition of the target interaction profile. 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 three spins each into a 2 × 3 square lattice in Fig. 4. To simulate the decoupled chains, N 3 couplings are set to zero (see Fig. 6 in Methods section) 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 and 4. They follow closely to the ideal target dynamics (green line).

Fig. 4 Dynamical changing of the target Hamiltonian. The Fourier filtering function, F(ϕ) in Eq. (6) 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 three spins each into a 2 × 3 square lattice. a The correlation measure between spins 1 and 4 defined by \(C_{14}(t) = \langle S_z^1S_z^4\rangle - \langle S_z^1\rangle \langle S_z^4\rangle\). b The correlation measure between spins 1 and 2 defined by \(C_{12}(t) = \langle S_z^1S_z^2\rangle - \langle S_z^1\rangle \langle S_z^2\rangle\). 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 Full size image

Proposal for experimental implementation

The experimental implementation of the multi-pulse scheme can be achieved in a trapped ion system such as 171Yb+ trapped in a radio-frequency (Paul) trap. When the confining potential is sufficiently anisotropic, laser-cooled ions will form a linear chain.53 The hyperfine states |↓〉 ≡ |2S 1/2 , F = 0, m F = 0〉 and |↑〉 ≡ |2S 1/2 , F = 1, m F = 0〉 form the spin (qubit) states for this ion.54

The inter-spin interactions in Eq. (1) can be simulated12,13 by shining the ions with a laser beam that imparts optical dipole forces, using the Mølmer–Sørensen scheme.6 When the frequencies of the laser beams (referred to as the Mølmer–Sørensen detuning in this manuscript) are appropriately chosen to off-resonantly excite the center-of-mass collective vibrational (phonon) modes of the ion chain, a long-range interaction in the form of Eq. (2) can be obtained. The global sign of \(\hat H_{{\mathrm{int}}}\) can be flipped by changing the Mølmer–Sørensen detuning, with additional laser beams improving the accuracy as discussed in Supplementary Materials.