Universal digital quantum simulation with trapped ions

Supplementary Online Material

Experimental details

For each simulation a string of 40Ca+ ions is loaded into a linear Paul trap. Qubits are encoded in two internal states of each ion. We use the meta-stable (lifetime $\sim$ 1 s) $|\!\uparrow\rangle=|D_{5/2},m_{j}=3/2\rangle$ state and the $|\!\downarrow\rangle=|S_{1/2},m_{j}=1/2\rangle$ ground state, where $m_{j}$ is the magnetic quantum number, for the experiments with 2 ions and $|\!\uparrow\rangle=|D_{5/2},m_{j}=-1/2\rangle$, $|\!\downarrow\rangle=|S_{1/2},m_{j}=-1/2\rangle$ 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 $S_{1/2}\leftrightarrow P_{1/2}$ transition at 397 nm. Optical pumping and resolved-sideband cooling on the $|\!\downarrow\rangle\leftrightarrow|D_{5/2},m_{j}=5/2\rangle$ transition prepare the ions in the state $|\!\downarrow\rangle$ 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 $S_{1/2}\leftrightarrow P_{1/2}$ transition on a photomultiplier tube and/or CCD camera. Instances where fluorescence light is detected correspond to the ion being in the state $|\!\downarrow\rangle$, instances where it is not correspond to the ion being in state $|\!\uparrow\rangle$. 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. $P_{0}$ (probability for 0 ions being bright), $P_{1}$ (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 (?).

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.

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

where $\sigma_{\phi}{=}\cos{\phi}~{}\sigma_{x}+\sin{\phi}~{}\sigma_{y}$ and $\sigma^{j}_{k}$ denotes the $k$-th Pauli matrix acting on the $j$-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 $\approx$ 30 $\mu s$.

In our experiments the coupling between the $j^{th}$ ionic qubit and a laser beam at a frequency that is resonant with the qubit transition is well described by the Hamiltonian $H_{3}{=}\hbar\Omega\sigma_{\phi}^{j}$. Here, $\Omega$ is the Rabi frequency which represents the coupling strength between the ion and the laser field. $O_{3}$ is realised via such a resonant laser pulse using the global beam, where $\theta_{3}{=}\Omega t$, $\phi_{3}$ is given by the laser phase and $t$ is the duration of the pulse i.e. $O_{3}{=}e^{-i\sum_{i}H_{3}t}$.

The interaction between an ionic qubit and a beam at a frequency detuned from resonance by $\Delta\gg\Omega$ is given by the Hamiltonian $H_{2}{=}\hbar\Omega^{2}/(4\Delta)\sigma_{z}^{i}$, describing an AC-Stark effect caused by off-resonant coupling to the $S_{1/2}$-$D_{5/2}$ transition and to very off-resonantly driven dipole transitions. $O_{2}$ is realised by such a detuned laser pulse using the global beam, with $\theta_{3}{=}\Omega t$ and $t$ is the duration of the pulse i.e. $O_{2}{=}e^{-i\sum_{i}H_{2}t}$. This same interaction is used to create $O_{1}$ by using an addressed beam.

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

Two-spin simulations

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 $O_{2}$ and $O_{4}$ 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 $O_{2}(\pi)$ pulses took $\approx 30\mu$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 $\pi/(4\sqrt{2}),\pi/(8\sqrt{2}),\pi/(12\sqrt{2}),\pi/(16\sqrt{2})$, respectively. We note that these fractions relate to a point in the evolution of the simulated system where a maximally entangled state is created ($\pi/2\sqrt{2}$). The pulse length $t_{MS}$ of the $O_{4}$ operations is set by the laser detuning $\delta$ (see section 1.1). For figures 1A i. - iv. detunings were chosen that yield operation times of 120, 60, 40 and 30$~{}\mu$s 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 $O_{3}$ operation) at $\approx 600~{}\mu$s.

1.2 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 $O_{2}(\pi/16)$ and $O_{4}(\pi/16,0)$ operations, of 10 $\mu$s and 30 $\mu$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 (?, ?). 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 (?).

1.3 Higher-order Trotter approximation

Consider a Hamiltonian with two terms $H{=}A{+}B$. A first-order Trotter approximation is:

which has errors on the order of $t^{2}/n$. A second-order Trotter-Suzuki approximation (?, ?) is:

which has errors on the order of $t^{3}/n$. 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 $B{=}J{=}1$, using both first- and second-order approximations. Figure 3 shows results and details. For the first-order simulation we use building blocks $O_{4}(\pi/8,0)$ and $O_{2}(\pi/8)$, which is seen to poorly reproduce the ideal evolution of the initial state $|\!\uparrow\uparrow\rangle$. For the second-order simulation we split the $O_{2}(\pi/8)$ operator into two pieces (each $O_{2}(\pi/16)$) 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 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 $\theta=\pi/4$. As shown in Figure 2 in the main text, simulations were built from $O_{2}(\pi/16)$, $O_{4}(\pi/16,0)$, $O_{4}(\pi/16,\pi)$, $O_{3}(\pi/4,0)$ operations, of duration 10, 30, 30, 5 $\mu$s respectively.

Simulations with more than two spins

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

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 $O_{4}$ also couples all pairs of spins (ionic qubits) with equal strength. Each digital time step of a simulation requires the sequence $O_{4}O_{2}$. 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 $O_{2}$ operation in the digital sequence.

1.2 Aysmmetric and nearest-neighbour

While our $O_{4}$ 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 (?). In particular, we can use the pulse sequence $O_{4}(\theta/2,\phi)O_{1}(\pi/2,n)O_{4}(\theta/2,\phi)$ to exclude ion $n$ 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 ($O_{4}$) simulates the evolution due to an interaction between all spin pairs with equal strength for a phase $\theta$, the next three operations simulate the evolution due to an interaction between one pair of spins for a phase $2\theta$. The overall effect is equivalent to evolving the system for a phase $\theta$ 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.

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 $n$-body spin interactions, which is derived and discussed in detail in (?). 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 $\pi/4$ 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 $O_{3}$ instead of $O_{2}$. This is because the three-body interaction that we simulate is $\sigma_{z}^{1}\sigma_{x}^{2}\sigma_{x}^{3}$, and the correct transverse field axis is therefore $\sum_{i}\sigma_{y}^{i}$. An alternative approach would be to use two extra pulses on spin 1 to rotate the axis of its spin-spin interaction to the $x$ basis, at the expense of two more operations for each digital step.

3 Process bounding method

Quantum process tomography enables a complete reconstruction of the experimental quantum process matrix (?, ?) 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 $12^{n}$ expectation values must be estimated to reconstruct the process matrix of an $n$ 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 (?) 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, $\{\psi_{i}\}$ and $\{\phi_{i}\}$, are complementary if $|\langle\psi_{i}|\phi_{i}\rangle|^{2}=1/N$ for all $i$. 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 $U$ a truth table shows the probability for measuring the ideal output state ($\{U\psi_{i}\}$ and $\{U\phi_{i}\}$) for each input state in a basis set. If we define the fidelity (overlap) of truth-table $i$ with its ideal case as $F_{i}$, then the process fidelity $F_{p}$ is bound above and below in the following way:

It is useful to note that the truth table fidelity is equivalent to the average output state fidelity. Therefore, for an $n$ qubit process, the requirements are to prepare two complementary sets of $2^{n}$ 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.

We bounded the process fidelity of the 3-body operation $U_{3}(\theta)=e^{-i\theta\sigma_{z}^{1}\sigma_{x}^{2}\sigma_{x}^{3}}$ for $\theta{=}\pi/4$ and $\pi/8$, considered in Figure 3C of the main text. As the first basis set we chose the 8 separable eigenstates, i.e. $|0\rangle|++\rangle_{x},|0\rangle|+-\rangle_{x},...,|1\rangle|--\rangle_{x}$ (where we now use the conventional qubit state notation for simplicity, and $|\pm\rangle_{x}=(|0\rangle\pm|1\rangle)/\sqrt{2}$). 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 $\theta=\pi/4$, are presented in Table 1.

where $q(\phi_{i})$ is the measured value of the parity at $\phi_{i}$ and $\alpha_{i}=\pm 1$ 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 $0.850(8)\leq F_{process}\leq 0.908(6)$. This procedure was repeated for $\theta{=}\pi/8$, yielding very similar results: $0.839(9)\leq F_{p}\leq 0.909(7)$.

3.2 Process bounding the 6-body interaction

We bounded the process fidelity of the 6-body operation $U_{6}(\theta)=e^{-i\theta\sigma_{y}^{1}\sigma_{x}^{2}\sigma_{x}^{3}\sigma_{x}^{4}\sigma_{x}^{5}\sigma_{6}^{3}}$ for $\theta{=}\pi/4$, 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 $n$ qubits ($\Psi{=}\cos{\theta}|0\rangle^{\otimes n}+\sin{\theta}|1\rangle^{\otimes n}$) is given by

where the $n$ values of $\phi$ are equally spaced by $\pi/n$ and alternately correspond to parity maxima ($\alpha=+1$) and minima ($\alpha=-1$). The requirement to measure the parity at $n$ different angles for an $n$-qubit state reflects the increasing number of possible entanglement partitions with $n$.

In total therefore $2^{n}(n+1)+2^{n}$ expectation values have to be measured to bound these many-body processes: $2^{n}(n+1)$ for the basis that becomes entangled and; $2^{n}$ for the separable eigenbasis. This compares well with the $12^{n}$ required for full quantum process tomography. For three qubits this means 40 instead of 1728 and for six qubits 512 instead of $2,985,984$.

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. $|000111\rangle$) are converted, by $U_{6}$, to states that are eigenstates of $\sigma_{z}$ 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 ($|000000\rangle$ and $|111111\rangle$) 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 (?). 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.

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.

5 Error sources in gate operations

For previous work on error sources in our quantum operations we refer the reader to (?, ?, ?). A conclusion from this work is that fluctuations in the laser-ion coupling strength $\Omega$ 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 $\theta$ of the $O_{4}$, $O_{1}$ and $O_{2}$ operations are all proportional to $\Omega^{2}$ 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 $\Omega^{2}$) 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 $\Omega$ is constant for each simulation sequence ($\approx 1~{}ms$), 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 $\delta\Omega/\Omega$.

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 (?, ?).

5.2 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 $\theta$ of the $O_{4}$ operation used in each trotter step is set incorrectly by only $1\%$. 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 (?).

References and Notes