Universal digital quantum simulation with trapped ions

# Universal digital quantum simulation with trapped ions

A digital quantum simulator is an envisioned quantum device that can be programmed to efficiently simulate any other local system. We demonstrate and investigate the digital approach to quantum simulation in a system of trapped ions. Using sequences of up to 100 gates and 6 qubits, the full time dynamics of a range of spin systems are digitally simulated. Interactions beyond those naturally present in our simulator are accurately reproduced and quantitative bounds are provided for the overall simulation quality. Our results demonstrate the key principles of digital quantum simulation and provide evidence that the level of control required for a full-scale device is within reach.

While many natural phenomena are accurately described by the laws of quantum mechanics, solving the associated equations to calculate properties of physical systems, i.e. simulating quantum physics, is in general thought to be very difficult [1]. Both the number of parameters and differential equations that describe a quantum state and its dynamics, grow exponentially with the number of particles involved. One proposed solution is to build a highly controllable quantum system that can efficiently perform the simulations [2]. Recently, quantum simulations have been performed in several different systems [3, 4, 5, 6, 7, 8, 9, 10, 11, 12], largely following the analog approach [2] whereby an analogous model is built, with a direct mapping between the state and dynamics of the simulated system and those of the simulator. An analog simulator is dedicated to a particular problem, or class of problems.

A digital quantum simulator [13, 14, 2, 15] is a precisely controllable many-body quantum system on which a universal set of quantum operations (gates) can be performed, i.e. a quantum computer [16]. The simulated state is encoded in a register of quantum information carriers, and the dynamics are approximated with a stroboscopic sequence of quantum gates. What makes this device so special is that it can, in principle, be reprogrammed to efficiently simulate any local quantum system [13] and is therefore referred to as a universal quantum simulator. Furthermore, there are known methods to efficiently correct for and quantitatively bound experimental error in large-scale digital simulations [17].

We report on digital simulations using a system of trapped ions. We focus on simulating the full time evolution of networks of interacting spin-1/2 particles, which are models of magnetism [18] and exhibit rich dynamics. We do not use error correction, which has been demonstrated separately in our system [19], and must be included in a full-scale fault-tolerant digital quantum simulator.

The central goal of a quantum simulation is to calculate the time-evolved state of a quantum system . In the case of a time-independent Hamiltonian the form of the solution is . A digital quantum simulator can solve this equation efficiently for any local quantum system [13], i.e. where contains a sum of terms that operate on a finite number of particles, due to interaction strengths that fall of with distance for example. In this case the local evolution operators can be approximated using a fixed number of operations from a universal set. However, these terms do not generally commute . This can be overcome using the Trotter approximation [20, 13], , for integer , which is at the heart of the digital quantum simulation algorithm. For finite the Trotter error is bounded and can be made arbitrarily small. The global evolution of a quantum system can therefore be approximated by a stroboscopic sequence of many small time-steps of evolution due to the local interactions present in the system. The digital algorithm can also be applied to time-dependent Hamiltonians and open quantum systems [13, 16, 21, 15].

Our simulator is based on a string of electrically trapped and laser-cooled calcium ions (see [22]). The and Zeeman states encode a qubit in each ion. Simulated states are encoded in these qubits and manipulated by laser pulses that implement the operation set: ; ; ; . Here and denotes the -th Pauli matrix acting on the -th qubit. is an effective qubit-qubit interaction mediated by a common vibrational mode of the ion string [23]. Recent advances have seen the quality of these operations increase appreciably [24]. For our simulations, we define dimensionless Hamiltonians , i.e. such that and the system evolution is quantified by a unitless phase .

We begin by simulating an Ising system of two interacting spin-1/2 particles, which is an elementary building block of larger and more complex spin models, and was one of the first systems to be simulated with trapped ions following an analog approach [6, 25]. The Hamiltonian is given by . The first bracketed term represents the interaction of each spin with a uniform magnetic field in the -direction and the second an interaction between the spins in an orthogonal direction. The interactions do not commute, giving rise to non-trivial dynamics and entangled eigenstates. Each spin is mapped directly to an ionic qubit (, ). The dynamics are implemented with a stroboscopic sequence of and gates, representing the magnetic field and spin-spin evolution operators, respectively. We first simulate a time-independent case which couples the initial state to a maximally entangled superposition of and (Fig 1A). The simulated dynamics converge closer to the exact dynamics as the digital resolution is increased. The overall simulation quality is quantified using quantum process tomography (QPT) [26], yielding a process fidelity of % at the finest digital resolution used. In [22] we show how higher-order Trotter decompositions can be used to achieve more accurate digital approximations with fewer operations.

We now consider a time-dependent case where increases linearly from 0 to 4 during a total evolution . In the limit , spins initially prepared in the paramagnetic ground state of the magnetic field () will evolve adiabatically into the anti-ferromagnetic ground state of the final Hamiltonian: an entangled superposition of the eigenstates and . As a demonstration, we approximate the continuous dynamics, for , using a stroboscopic sequence of 24 and gates, and measure the populations in the basis (Fig 1B). The evolution closely follows the exact case and an entangled state is created ( tangle [27]). Full quantum state reconstructions are performed after each digital step, yielding fidelities between the ideal digitised and measured state of at least , and overlaps with the instantaneous ground state of no less than . Note that the observed oscillation in expectation values is a diabatic effect, as excited states become populated.

More complex systems with additional spin-spin interactions in the (‘XY’ model) and (‘XYZ’ model) directions can be simulated by reprogramming the operation sequence. The dynamics due to an additional spin-spin interaction in the -direction is simulated by adding another operation to each step of the Ising stroboscopic sequence (with ). A third spin-spin interaction in the -direction is realised by adding an gate sandwiched between a pair of operations set to rotate the reference frame of the qubits. In the simulated dynamics of the initial state under each model, for a fixed digital resolution of and up to 12 trotter steps (Fig 2), up to 24, 48 and 84 gates are used for the Ising, XY and XYZ simulations, respectively. This particular initial state is chosen because the ideal evolution is different for each model. The results show close agreement with the exact dynamics and results from QPT after four digital steps yield process fidelities, with the exact unitary evolution, of 88(1)%, 85(1)% and 79(1)%, for the Ising, XY and XYZ respectively. With perfect operations the Trotter error would be less than 1% in each case. Note that while analog simulations of Ising models have previously been demonstrated in ion traps [6, 8], XY and XYZ models have not.

The digital approach allows arbitrary interaction distributions between spins to be programmed. For three-spin systems, we realise various interactions that give rise to the dynamical evolutions of the initial state (Fig 3). Fig 3A shows a system supporting interactions between all spin pairs with equal strength, and between each spin and a transverse field. The initial state couples equally to , and , while the strength of the field determines the amplitude and frequency of the dynamics. For the case shown () an equal superposition of the coupled states is periodically created (an entangled W state [28]). Fig 3B shows how non-symmetric interaction distributions can be programmed, using sequences of and to add spin-selective interactions. The interaction between one spin pair is dominant. Due to this broken symmetry, one coupled state () is populated faster than the others, yielding more complex dynamics than in the symmetric case. Fig 3C demonstrates the ability to simulate -body interactions; specifically . A clear signature is observed: direct coupling between eigenstates and , periodically producing an entangled GHZ state [28]. Many-body spin interaction of this kind are an important ingredient in the simulation of systems with strict symmetry requirements [29] or spin models exhibiting topological order [30]”. Measurements in other bases and simulations of nearest-neighbour and many-body interactions with a transverse field using over 100 gates are presented in [22].

Fig 4A shows the observed dynamics of the 4-spin state under a long-range Ising-type interaction. The rich structure of the dynamics reflects the increased complexity of the underlying Hamiltonian: oscillation frequencies correspond to the energy gaps in the spectrum. This information can be extracted via a Fourier transform of the data (see [22]). Specific energy gaps could be targeted by preparing superpositions of eigenstates via an initial quasi-adiabatic digital evolution [10], or study the complex non-local correlations generated by this model. Fig 4B shows the observed dynamics for our largest simulation: a six spin many-body interaction, which directly couples the states and , periodically producing a maximally entangled GHZ state.

Direct quantification of simulation quality for more than two qubits is impractical via QPT: for three qubits, expectation values must be measured for 1728 experimental configurations, and this increases exponentially with qubit number ( for six qubits). However, the average process fidelity () can be bounded more efficiently [31]. We demonstrate this for the 3- and 6-spin simulations of Fig 3C and 4B, respectively. Bounds of (3 spins) and (six spins) are obtained at , using 40 and 512 experimental configurations respectively [22]. The largest system for which a process fidelity has previously been measured is 3 qubits [32]. Note that a different measure of process quality is given by the worst-case fidelity, over all input states, and may be better used to assess errors in future full-scale fault-tolerant simulations. Regardless of the measure used, the error in large-scale digital simulations built from finite-sized operations can be efficiently estimated. Each operation can be characterized with a finite number of measurements, then the error in any combination can be chained [33]. In order to exploit this the number of ionic qubits on which our many-qubit operators can act must be restricted.

The dominant effect of experimental error can be seen in Fig 3B and 4B: the dynamics damps due to decoherence processes. Laser frequency and ambient magnetic field fluctuations are far from the leading error source: in the absence of coherent operations, qubit lifetimes are over an order of magnitude longer (coherence times  ms) than the duration of experiments ( ms). The current leading sources of error, which limit both the simulation complexity and size, are thought to be laser intensity fluctuations [22]. This is not currently a fundamental limitation and, once properly addressed, should enable an increase in simulation capabilities.

The digital approach can be combined with existing tools and techniques for analog simulations to expand the range of systems that can be simulated. In light of the present work, and current ion trap development [34], digital quantum simulations involving many tens of qubits and hundreds of high-fidelity gates seems feasible in coming years.

Acknowledgements
We thank Wolfgang Dür, Alán Aspuru-Guzik and Michael Brownnutt for discussions. We acknowledge financial support by the Austrian Science Fund (FOQUS), the European Commission (AQUTE), the Institut für Quanteninformation GmbH, IARPA, and two Marie Curie International Incoming Fellowships within the 7th European Community Framework Programme.

## 1 Experimental details

For each simulation a string of Ca ions is loaded into a linear Paul trap. Qubits are encoded in two internal states of each ion. We use the meta-stable (lifetime  1 s) state and the ground state, where is the magnetic quantum number, for the experiments with 2 ions and , for the experiments with more ions. These states are connected via an electric quadrupole transition at 729 nm and an ultra-stable laser is used to perform qubit operations. At the start of each experiment the ions are Doppler cooled on the transition at 397 nm. Optical pumping and resolved-sideband cooling on the transition prepare the ions in the state and in the ground state of the axial center-of-mass vibrational mode.

The internal states of the ions are measured by collecting fluorescence light on the transition on a photomultiplier tube and/or CCD camera. Instances where fluorescence light is detected correspond to the ion being in the state , instances where it is not correspond to the ion being in state . Each simulation experiment (consisting of cooling, state preparation, simulated time evolution and detection) is repeated many times to obtain enough statistics (at least 200 times per data point). The photomultiplier measures the collective fluorescence state of the ion string, from which the probability for any number of ions in the string to be bright can be extracted i.e. (probability for 0 ions being bright), (probability for any 1 ion to be bright) etc. More information can be obtained by using the CCD camera. By defining regions of interest on the camera sensor, the fluorescence of each ion can be measured individually, from which the logical state of any individual ion can be extracted. Measurements in other bases are achieved by applying single-qubit operations to the ions before measurement, to map eigenstates of the desired observable to the logical eigenstates.

More detailed information about our experimental setup and techniques can be found in the Ph.D. thesis of Gerhard Kirchmair [35].

Quantum simulations are carried out in two different ion traps. Both are based on linear Paul traps of similar design and working with only slightly different experimental parameters. Trap 1 has trapping frequencies  1.2 MHz and 3 MHz in the axial and radial directions and operates at a magnetic field of 4 G. Trap 2 has trapping frequencies  1.1 MHz and  3 MHz in the axial and radial directions and operates at a magnetic field of 3.4 G.

Trap 2 offers the possibility of trapping larger strings of ions due to an improved vacuum quality and all experiments involving more than two ions were carried out in this trap.

### 1.1 Universal operation set

In this section we explain how the set of operations used in the experiments are implemented. To recap, the operations are

 O1(θ,i) = exp(−iθσiz) (1) O2(θ) = exp(−iθ∑iσiz) (2) O3(θ,ϕ) = exp(−iθ∑iσiϕ) (3) O4(θ,ϕ) = exp(−iθ∑i

where and denotes the -th Pauli matrix acting on the -th qubit.

Each operation is implemented using one of two laser beam paths, which impinge on the ion string from different directions. One beam illuminates all ions equally and can be used to perform global spin rotations on all ions. This is referred to as the ‘global beam’. The direction of this beam has a large overlap with the axis of the ion string, such that it can couple to the axial motional state. The other beam is tightly focussed, impinges at a 90 degree (68 for trap 2) angle to the ion string and can be used to address ions individually. This is referred to as the ‘addressed beam’. The particular ion illuminated by the addressed beam can be changed using an electro-optic deflector in 30 .

In our experiments the coupling between the ionic qubit and a laser beam at a frequency that is resonant with the qubit transition is well described by the Hamiltonian . Here, is the Rabi frequency which represents the coupling strength between the ion and the laser field. is realised via such a resonant laser pulse using the global beam, where , is given by the laser phase and is the duration of the pulse i.e. .

The interaction between an ionic qubit and a beam at a frequency detuned from resonance by is given by the Hamiltonian , describing an AC-Stark effect caused by off-resonant coupling to the - transition and to very off-resonantly driven dipole transitions. is realised by such a detuned laser pulse using the global beam, with and is the duration of the pulse i.e. . This same interaction is used to create by using an addressed beam.

is an effective qubit-qubit (spin-spin) interaction generated via a Mølmer-Sørensen type interaction [23]. For this, the global beam is used with a bichromatic laser pulse, whose frequencies are detuned from the carrier by , where is the frequency of the axial centre of mass vibrational mode (1.2 MHz) and is between 10 and 100 kHz depending on the simulation. After a time , this generates the unitary , with , and given by the sum of the phases of the two light fields [36]. The unitaries in equations 1-4 form a universal set of operations [37]—any arbitrary unitary qubit evolution can be implemented using only these operations.

We note that by choosing to be a far detuned beam, there are no phase stability requirements between the global and addressed laser paths.

## 2 Two-spin simulations

### 2.1 Ising demonstration

#### Time-independent dynamics

Figure 1A, in the main text, shows results from digital simulations of a time-independent case of a two-spin Ising system, for increasing levels of digital approximation. Each simulation corresponded to a stroboscopic sequence of and operations, where the former simulates an interaction with an external field, and the latter an orthogonal spin-spin interaction. The laser power was set such that pulses took s. The shorter phase evolutions required for each simulation were achieved by varying the pulse length to the correct fraction of this time. For figures 1A i. - iv. these fractions were , respectively. We note that these fractions relate to a point in the evolution of the simulated system where a maximally entangled state is created (). The pulse length of the operations is set by the laser detuning (see section 1.1). For figures 1A i. - iv. detunings were chosen that yield operation times of 120, 60, 40 and 30s respectively. The power of the laser was adjusted to realise the required phase evolution in each case. Varying the operation lengths in this way enabled the total simulation time to be kept constant (up to small changes due to the operation) at s.

The caption of figure 1A quotes quantum process fidelities for the lowest and highest resolution digital simulations. These are calculated after performing full quantum process tomography, which allows the quantum process matrix to be reconstructed for a given point in the simulation. For details on process tomography we refer to references  [26, 38]. In summary we input a complete set of states into the simulation, and for each measure the output in a complete basis. A maximum-likelihood reconstruction algorithm is used to determine the most likely quantum process to have produced our data. Monte-Carlo error analysis is then employed to estimate the uncertainty in the process. Figure 1 shows the experimentally reconstructed quantum process matrices after a simulated phase evolution of . The matrices clearly demonstrate the improved simulation quality at higher digital resolution, and that our simulations are of high quality across the full Hilbert space.

#### Time-dependent dynamics

Figure 1B in the main text shows results from a digital simulation of a time-dependent two-spin Ising model. The spins are first prepared in the ground state of the external magnetic field, the simulated dynamics corresponds to slowly increasing the strength of the spin-spin interaction such that the state evolves to an approximation of the joint ground state, which is highly entangled. The continuous dynamics are approximated by an 8 step digital simulation built from and operations, of 10 s and 30 s duration respectively.

Figure 2 reproduces and extends the data and details in the main text, showing experimentally reconstructed density matrices at all 9 stages of the digital simulation (including the initial state). These matrices are constructed via a full quantum state tomography [26, 38]. Maximum-likelihood tomography is used to assure a physical state and Monte-Carlo analysis is used to estimate errors in derived quantities. The fidelity and entanglement properties quoted in Figure 1B are calculated from these states. The fidelity is between the measured and ideal digital case, assuming perfect operations. The entanglement is quantified by the tangle which can be readily calculated from the density matrices [27].

#### Higher-order Trotter approximation

Consider a Hamiltonian with two terms . A first-order Trotter approximation is:

 U(t)=(e−iAt/ne−iBt/n)n+O(t2/n) (5)

which has errors on the order of . A second-order Trotter-Suzuki approximation [39, 16] is:

 U(t)≈(e−iBt/2ne−iAt/ne−iBt/2n)n+O(t3/n) (6)

which has errors on the order of . By splitting the second evolution operator into two pieces and rearranging the sequence, a closer approximation to the correct dynamics is achieved. Practically this means that a more accurate digital approximation can be achieved at the expense of more operations, but for the same total phase evolution for each step. To illustrate this concept we performed a digital simulation of the two-spin Ising model for , using both first- and second-order approximations. Figure 3 shows results and details. For the first-order simulation we use building blocks and , which is seen to poorly reproduce the ideal evolution of the initial state . For the second-order simulation we split the operator into two pieces (each ) and rearranged the sequence according to Eq. 6, thereby achieving a much more accurate simulation. Since each evolution operator is simulated directly with our fundamental operations, the higher-order approximation can be employed with little overhead. However, in the more general case where operators must be constructed, such as in our simulation of the XYZ model, there is some finite overhead with simulating any evolution operator regardless of how short it is. In this case there will be a trade-off between increased digital resolution offered by higher-order approximations and additional experimental error introduced through using more operations.

### 2.2 Digital simulations of the XY and XYZ models

Figure 2 in the main text shows results of digital simulations of time-independent instances of the Ising, XY and XYZ models. Figure 4 in this document shows experimentally reconstructed process matrices of these simulations after 4 digital time steps. This corresponds to a simulated phase evolution of . As shown in Figure 2 in the main text, simulations were built from , , , operations, of duration 10, 30, 30, 5 s respectively.

## 3 Simulations with more than two spins

### 3.1 Ising-type models

#### Long-range

Our basic set of operations is well suited to simulating the Ising model with long-range interactions:

 H=J∑i

which corresponds to a system with interactions between each pair of spins with equal strength J, and a transverse field of strength B. This is because the effective interaction underlying also couples all pairs of spins (ionic qubits) with equal strength. Each digital time step of a simulation requires the sequence . In Figure 5 we give a more complete set of results for the simulations of three and four spin cases shown in Figures 3A and 4A, of the main text. Specifically, time dynamics measured in a complementary basis and results for different strength transverse fields are shown. The simulated transverse field strength is adjusted by varying the phase evolution of each operation in the digital sequence.

#### Aysmmetric and nearest-neighbour

While our operation is best suited for simulating symmetric interactions between all pairs of qubits (spins), it is possible to engineer interactions that break this symmetry. For this we make use of refocussing techniques in the spirit of nuclear magnetic resonance quantum computing as described in [37]. In particular, we can use the pulse sequence to exclude ion from the long-range spin-spin interaction. In principle, sequences of this form can be repeated to simulate any arbitrary spin-spin coupling network.

Two examples with increasing difficulty, in terms of the number of operations required, are given in Figure 6. In Figure 6A the system considered has an asymmetric interaction between the spins: specifically, the interaction strength between one pair is three times larger than any other. A subset of these results is shown in Figure 3B in the main text. Each digital step is constructed from four operations: the first () simulates the evolution due to an interaction between all spin pairs with equal strength for a phase , the next three operations simulate the evolution due to an interaction between one pair of spins for a phase . The overall effect is equivalent to evolving the system for a phase due to the desired asymmetric Hamiltonian.

Figure 6B considers a spin system with nearest-neighbour interactions. The large number of operations required for each digital step causes the simulated dynamics to damp due to decoherence processes in the operations themselves (largely the results of laser intensity fluctuations). Clearly these simulations require significantly more gate operations than the long-range Ising model. From an experimental point of view it might therefore be advantageous to use a different set of universal operations for simulating such systems. For this, spin-spin interactions between neighbouring ions could be realised using lasers focussed on pairs of ions. This will be the subject of future work.

### 3.2 3-body interaction with additional transverse field

In the main text we presented simulations of a three-body interaction (Figure 3C). The circuit decomposition for three-body interactions is a special case of a general scheme to simulate -body spin interactions, which is derived and discussed in detail in [40]. We now give simulation results with an additional transverse field. This is particularly challenging as a large number of operations, many of which have large fixed phase evolutions, are required for each digital step. Note that the three-body interaction alone can be simulated for any phase evolution using only three operations, as shown in Figure 3C in the main text. However, the Trotter approximation must be employed in the case of an additional transverse field, costing 4 operations (3 for the three-body interaction and one for the magnetic field interaction) for each digital step of the total phase evolution. Figure 7 shows results, first for the case where the field is zero but following a stroboscopic approach with fixed operation settings, and second with a non-zero field. A coarse digital resolution of is chosen so as to observe some dynamics before decoherence mechanisms equally distribute population among each possible spin state.

Note that here, for the first time, we simulate a transverse field using instead of . This is because the three-body interaction that we simulate is , and the correct transverse field axis is therefore . An alternative approach would be to use two extra pulses on spin 1 to rotate the axis of its spin-spin interaction to the basis, at the expense of two more operations for each digital step.

### 3.3 Process bounding method

Quantum process tomography enables a complete reconstruction of the experimental quantum process matrix [26, 16] from which any desired property, such as the process fidelity, can be calculated. However, the number of measurements required grows exponentially with the qubit number. In an ion trap system expectation values must be estimated to reconstruct the process matrix of an qubit process. This number is already impractical for processes involving more than two qubits: it simply takes too long to carry out the measurements with sufficient precision, while maintaining accurate control over experimental parameters.

In [31] it is shown that the overall process fidelity can be bounded without reconstructing the process matrix, and with a greatly reduced number of measurements. In summary, the technique requires classical truth tables to be measured for two complementary sets of input basis states. The two sets, and , are complementary if for all . Conceptually this means that a measurement in one basis provides no information about the outcome of a subsequent measurement in the other basis. In this way there is no redundancy in these measurements and maximal information is returned.

For a unitary quantum process a truth table shows the probability for measuring the ideal output state ( and ) for each input state in a basis set. If we define the fidelity (overlap) of truth-table with its ideal case as , then the process fidelity is bound above and below in the following way:

 F1+F2−1≤Fp≤Min{F1,F2} (8)

It is useful to note that the truth table fidelity is equivalent to the average output state fidelity. Therefore, for an qubit process, the requirements are to prepare two complementary sets of input states, and to measure the probability of obtaining the correct output state (the state fidelity) in each case. The technique is highly dependent on the particular process to be characterised: the challenge is to choose complementary sets of states that can be accurately prepared and for which the output state fidelities can be accurately measured.

#### Process bounding the 3-body interaction

We bounded the process fidelity of the 3-body operation for and , considered in Figure 3C of the main text. As the first basis set we chose the 8 separable eigenstates, i.e. (where we now use the conventional qubit state notation for simplicity, and ). These states can be created experimentally using a sequence of coherent laser pulses that include both global and addressed beams. The inverse of this pulse sequence, followed by fluorescence detection, is used to effectively perform a projective measurement in this eigenstate basis. From this measurement the average output state fidelity can be calculated directly. The results of these measurements, for , are presented in Table 1.

For the second basis set we chose the 8 states (where ). These input states are mapped to entangled output states by . We then map each output state to a GHZ-like state with the ideal form , using an additional set of local operations. The fidelity of an experimentally produced state (density matrix) with can be derived from 5 expectation values. The first two are the probability for finding the qubits all in () and the probability for finding them all in (). These can be estimated directly from fluorescence measurements. The last three are the parity of the state (spin correlations) measured in different bases. For this purpose we apply an operation to the output state, for three different values of and then measure the parity via fluorescence measurement and post-processing. It is straightforward to show that the fidelity of an experimentally produced state with the GHZ-like state is given by

 F(ρ,ψideal)=cos2(θ)P000+sin2(θ)P111+cosθsinθ33∑i=1αiq(ϕi) (9)

where is the measured value of the parity at and depending on whether the ideally expected parity is at a minimum or at a maximum.

The measurement results for this second set of input states are presented in Table 2. Together with the results of Table 1 and Equation 8 the process fidelity can be bound to . This procedure was repeated for , yielding very similar results: .

#### Process bounding the 6-body interaction

We bounded the process fidelity of the 6-body operation for , shown in figure 4B of the main text. Our method is conceptually equivalent to that for the 3-body case described above. As the first basis set of input states we chose the 64 separable eigenstates and directly measured in this basis to extract the output state fidelities. These results are split between Table 3 and Table 4. For the second set we chose a complementary basis which evolve into entangled states that are locally equivalent to GHZ states. In the 6-qubit case 2 populations and 6 parities are required for the state fidelity. These results are split between Tables 5 and 6. The fidelity of an experimentally produced state with a GHZ-like state of qubits () is given by

 F(ρ,Ψ)=cos2(θ)P00...0+sin2(θ)P111...1+cosθsinθnn∑i=1αiq(ϕi) (10)

where the values of are equally spaced by and alternately correspond to parity maxima () and minima (). The requirement to measure the parity at different angles for an -qubit state reflects the increasing number of possible entanglement partitions with .

In total therefore expectation values have to be measured to bound these many-body processes: for the basis that becomes entangled and; for the separable eigenbasis. This compares well with the required for full quantum process tomography. For three qubits this means 40 instead of 1728 and for six qubits 512 instead of .

Interestingly, further analysis of the data in Tables 5 and 6 shows that decoherence of the GHZ states is an error source. The 64 states can be separated into 4 groups, determined by the magnitude of the difference in the number of 0’s and 1’s in the input state. The possible values are 6, 4, 2, 0. Input states with a difference of 0 (e.g. ) are converted, by , to states that are eigenstates of rotations on any or all qubits. This makes them free of decoherence effects due to fluctuating magnetic fields, for example. Input states with the maximum difference ( and ) are converted to states that are maximally sensitive to these kinds of rotations and errors - by a factor of 6 times more than a single qubit [41]. The effect is to reduce the coherence between the populations, which would reduce the parity amplitude while keeping the population values the same. We should expect the average parity amplitude (over the 6 measurements) for the groups 6, 4, 2 and 0 to be progressively better in that order. The results are consistent: the average absolute parity amplitudes for groups 6, 4, 2, 0 are 0.58(2), 0.67(1), 0.71(1) and 0.76(1), respectively, while the average total populations are 0.78(5), 0.83(2), 0.81(1) and 0.83(2), respectively.

### 3.4 Fourier transform to extract energy gaps

Oscillation frequencies in the time evolution of observables are energy gaps in underlying Hamiltonian. A Fourier transform can extract this information. We now give a brief example for one of the observables measured in the 4-ion long-range Ising simulation (Figure 5C i. and Figure 4A), which shows the richest dynamics of all our simulations.

Figure 5D shows the spectrum of the ideal Hamiltonian and the probability distribution of the initial state used for the simulation, amongst the energy levels. Three of the nine energy levels are populated, therefore at most 3 energy gaps (oscillation frequencies) can be observed in the dynamics. However, the observed spectral amplitudes in the Fourier transform depend not only on the population distribution of the initial state, but also the observable coupling strengths and trace length (total simulated phase evolution). Figure 5D shows the Fourier transform of the black data trace in Figure 4A of the main text (and Figure  5C), which represents the total probably of finding all combinations of two spins up and two down. One of the three fundamental frequencies is clearly resolved.

### 3.5 Error sources in gate operations

#### Laser-ion coupling strength fluctuations

For previous work on error sources in our quantum operations we refer the reader to [42, 24, 36]. A conclusion from this work is that fluctuations in the laser-ion coupling strength are a significant source of experimental error. These fluctuations can be introduced by laser-intensity noise or thermally occupied vibrational modes, for example. Since the phase angles of the , and operations are all proportional to this error source is particularly important in our simulations.

We measured the laser intensity fluctuation, at the entrance into the ion-trap vacuum vessel, using a fast photo-diode. A slow oscillation in the intensity by between 1 and 2% was observed over periods of several minutes. The precise value in the range varies on a daily basis. This corresponds to a coupling-strength fluctuation (intensity is proportional to ) of approximately 1%. This measurement is an underestimate of intensity fluctuations at the point of the ions, due to possible beam-pointing fluctuations and wave-front aberrations.

The effect of coupling strength fluctuations on our digitised simulations was modelled, and compared to one of our key results: the time dynamics of a two-spin Ising system at the highest digital resolution used (shown in Figure 1A, panel iv, in the main text). The model makes the assumption that is constant for each simulation sequence (), as supported by our photo-diode measurements, but varies from sequence to sequence (in experiments expectation values are calculated by averaging a large number of repeated experimental sequences taken minutes apart). Fluctuations are incorporated by post-mixing a large number of simulated sequences, with each subjected to noise randomly sampled from a gaussian distribution with a standard deviation .

Figure 8 shows the results: the measured coupling strength fluctuation of 1% qualitatively reproduces the observed damping of the spin dynamics, while a much closer fit is obtained for a larger fluctuation of 2%. There are a large number of other errors sources that could contribute to deviations between the observed simulations and the ideal, as discussed in [24, 36].

#### Frequency shifts in simulated dynamics

The 3-spin transverse Ising model results, presented in Figure 3A in the main text, exhibit a slight frequency shift compared to the ideal digitised case. These effects could easily be the result of errors made in the setting-up/optimising of the gate operations required for the sequence. Figure 9 shows how the observed frequency mismatch would be expected if the phase angle of the operation used in each trotter step is set incorrectly by only . The sensitivity to these effects suggests the need for future work on developing even more accurate methods to optimise gate operations in the lab than are currently employed [24].