Quantum Approximate Optimization of the Long-Range Ising Model with a Trapped-Ion Quantum Simulator

TL;DR

This paper demonstrates the implementation of a low-depth Quantum Approximate Optimization Algorithm on a trapped-ion quantum simulator, effectively estimating ground states of long-range Ising models with up to 40 qubits, showing promising scalability and error analysis.

Contribution

The study presents the first scalable implementation of QAOA on a trapped-ion system for long-range Ising models, including optimization, benchmarking, and error analysis.

Findings

QAOA performance remains stable as system size increases.

Runtime is approximately independent of the number of qubits.

High-fidelity measurements enable effective sampling and optimization.

Abstract

Quantum computers and simulators may offer significant advantages over their classical counterparts, providing insights into quantum many-body systems and possibly improving performance for solving exponentially hard problems, such as optimization and satisfiability. Here we report the implementation of a low-depth Quantum Approximate Optimization Algorithm (QAOA) using an analog quantum simulator. We estimate the ground state energy of the Transverse Field Ising Model with long-range interactions with tunable range and we optimize the corresponding combinatorial classical problem by sampling the QAOA output with high-fidelity, single-shot individual qubit measurements. We execute the algorithm with both an exhaustive search and closed-loop optimization of the variational parameters, approximating the ground state energy with up to 40 trapped-ion qubits. We benchmark the experiment with bootstrapping heuristic methods scaling polynomially with the system size. We observe, in agreement with numerics, that the QAOA performance does not degrade significantly as we scale up the system size, and that the runtime is approximately independent from the number of qubits. We finally give a comprehensive analysis of the errors occurring in our system, a crucial step in the path forward towards the application of the QAOA to more general problem instances.

Equations

$$H=\underset{H_A}{\underbrace{\sum _{i<j}{J}_{ij}{\sigma }_i^x{\sigma }_j^x}}+\underset{H_B}{\underbrace{B\sum _i{\sigma }_i^y}}.$$

$$|\vec{\beta },\vec{\gamma }⟩=\prod _{k=1}^p{e}^{-i{\beta }_k(H_B/J_0)}{e}^{-i{\gamma }_k(H_A/J_0)}|{\psi }_0⟩,$$

$$\eta \equiv \frac{E(\vec{\beta },\vec{\gamma })-E_{max}}{E_{gs}-E_{max}},$$

$$\mathrm{TVD}$$

$${\mathrm{D}}_{\mathrm{K}-\mathrm{L}}$$

$$T=O\left(\sum _{q=1}^{p}q{ϵ}_q^2\right)\le O(p^2)$$

$$\eta \equiv \frac{E(\vec{\beta },\vec{\gamma })-E_{max}}{E_{gs}-E_{max}},$$

$$E(\beta ,\gamma )=E_I+E_{II}+E_{III}$$

$$E_I=B\sum _{i=1}^N\prod _{k\ne i}\cos \left(2\gamma J_{ik}\right)$$

$$E_{II}=-\frac{\sin \left(4\beta \right)}{2}\sum _{i,j}{J}_{ij}\sin \left(2\gamma J_{ij}\right)\prod _{k\ne i,j}\cos \left(2\gamma J_{ik}\right)$$

$$E_{III}=-\frac{{\sin }^2\left(2\beta \right)}{4}\sum _{s=\pm 1,i,j}{J}_{ij}\prod _{k\ne i,j}\cos \left(2\gamma \left(J_{ik}+(-1{)}^s{J}_{jk}\right)\right)$$

$$\frac{{\sigma }_E}{\Delta }=\frac{\sqrt{⟨E^2{⟩}_{\beta ,\gamma }-⟨E{⟩}_{\beta ,\gamma }^2}}{\Delta }$$

$$⟨f{⟩}_{\beta ,\gamma }:=\underset{T_{\beta },T_{\gamma }\to \infty }{\lim }\frac{1}{4T_{\beta }{T}_{\gamma }}\underset{-T_{\beta }}{\overset{T_{\beta }}{\int }}\underset{-T_{\gamma }}{\overset{T_{\gamma }}{\int }}f(\beta ,\gamma )d\beta d\gamma$$

$$⟨E_I{⟩}_{\beta ,\gamma }^2=B^2\sum _{i,j}\prod _{k=1}^N\cos (2\gamma J_{ik})\cos (2\gamma J_{jk})$$

$$⟨E_I{⟩}_{\beta ,\gamma }^2\lesssim 4N{B}^2{\left(\frac{3}{8}\right)}^{(N-1)/2}$$

$$⟨E_{II}{⟩}_{\beta ,\gamma }^2$$

$$⟨E_{III}{⟩}_{\beta ,\gamma }^2$$

$$⟨E{⟩}_{\beta ,\gamma }^2\lesssim N{\left(\frac{3}{8}\right)}^{N/2}\left[8B^2+J_0^2\right]\sim O(N\cdot {\left(3/8\right)}^{N/2})$$

$$⟨0^N|{\tilde{H}}^{\otimes N}{e}^{-i\gamma H_z}{H}^{\otimes N}|0^N⟩=\frac{1}{2^N}\sum _{x,y}⟨y|i^{{\sum }_i{y}_i+\tilde{f}(x)}|x⟩,$$

$$U_1(\pi /2)=e^{-i\frac{\pi }{2}{\sigma }_i^{\theta }}{e}^{-i\pi {\sigma }_i^{3\theta }}{e}^{-i\frac{\pi }{2}{\sigma }_i^{\theta }}{e}^{-i\frac{\pi }{4}{\sigma }_i^x},$$

$$J_{ij}={\Omega }^2{\nu }_R\sum _m\frac{b_{im}{b}_{jm}}{{\mu }^2-{\nu }_m^2}.$$

$$J_{ij}$$

$$J_{ij}^{\ell }$$

$$J_{ij}\approx \frac{J_0}{r^{\alpha }},$$

$$J_{ij}\approx \frac{J_0}{r^{{\alpha }^{\prime }}}{e}^{-{\beta }^{\prime }(r-1)}$$

Full text

A promising near-term application of quantum devices is the production of highly entangled states with metrological advantage or with properties of interest for many-body physics and quantum information processing. One possible approach to produce useful quantum states is to use quantum devices to perform adiabatic quantum computing Kadowaki1998 ; Farhi2000 , which in some cases may provide an advantage over classical approaches Hastings2013 . However, adiabatic quantum computing has stringent adiabaticity requirements that hinder its applicability on existing quantum platforms that have finite coherence times Richerme2013 .

Alternatively, hybrid quantum-classical variational algorithms may approximately solve hard problems in realms such as quantum magnetism, quantum chemistry mcclean2016 , and high-energy physics Kokail2018 . This is because the key resource of quantum computers and simulators is quantum entanglement, which is exactly what makes these many-body quantum problems hard. In a hybrid variational algorithm, entangled states are functions of variational parameters that are iteratively optimized by a classical algorithm. One example is the Quantum Approximate Optimization Algorithm Farhi2014 , which consists of a “bang-bang” protocol that can provide approximate answers in a time-efficient way, using devices with finite coherence times and without the use of error-correction harrow2016qaoa ; Lloyd2018 ; Preskill2018 ; Zhou2018 ; Hastings2019 .

Similarly to adiabatic quantum computing, the QAOA protocol encodes the objective function of the optimization problem in a target spin Hamiltonian. The optimization steps of the QAOA are based on unitary evolution under the target Hamiltonian and a non commuting “mixing” operator. In general, the QAOA relies on a classical outer loop to optimize the quantum circuit, aided by physical intuition bapat2018 ; verdon2018 ; Jiang2017 ; Wang2018 or observed structure of the variational parameters Zhou2018 ; brady2019 ; crooks2018 ; Mbeng2019 , producing fast, low-depth circuits for approximate solutions. The QAOA has also been proposed as an efficient way to produce entangled quantum states, such as the ground states of critical Hamiltonians, which gives access to their corresponding energies Ho2018a ; Ho2018b .

In this work, we employ a collection of interacting trapped-ion qubits to experimentally implement a specific instance of the QAOA, which is native to our quantum hardware. We focus on both the energy minimization of the quantum Hamiltonian and the combinatorial optimization of the corresponding classical problem. Both problems are encoded in the transverse field anti-ferromagnetic Ising Hamiltonian with long-range interactions:

[TABLE]

Here we set the reduced Planck’s constant $\hbar=1$, $\sigma_{i}^{\gamma}$ ($\gamma=x,y,z$) is the Pauli matrix acting on the $i^{\text{th}}$ spin along the $\gamma$ direction of the Bloch sphere, $J_{ij}>0$ is the Ising coupling between spins $i$ and $j$, which, in our case, falls off as a power law in the distance between the spins, and $B$ denotes the transverse magnetic field. It is well-known Koffel2012 that the Hamiltonian (1) exhibits a quantum phase transition for anti-ferromagnetic interactions with power law decay. One of the goals of this work is to find an approximation of the ground state energy both at the critical point $(B/J_{0})_{c}$, where $J_{0}$ is the average nearest-neighbour coupling, and in the case of $B=0$, optimizing the QAOA output for the classical Hamiltonian $H_{A}$. The realization of the QAOA entails a series of unitary quantum evolutions (see Fig. 1) under the non-commuting Hamiltonians $H_{A}$ and $H_{B}$ (defined under Eq. (1)) that are applied to a known initial state $|\psi_{0}\rangle$. The state obtained after $p$ layers of the QAOA is:

[TABLE]

where the evolution times (or, henceforth, “angles”) $\beta_{k}$ and $\gamma_{k}$ are variational parameters used in the $k$-th QAOA layer to minimize the final energy $E(\vec{\beta},\vec{\gamma})=\langle\vec{\beta},\vec{\gamma}|H|\vec{\beta},\vec{\gamma}\rangle$.

In order to implement the quantum optimization algorithm, each spin in the chain is encoded in the 2S1/2 $|F=0,m_{F}=0\rangle\equiv|\!\!\downarrow\rangle_{z}$ and $|F=1,m_{F}=0\rangle\equiv|\!\!\uparrow\rangle_{z}$ hyperfine “clock” states of a 171Yb+ ion (see Supplementary). In this work, depending on the number of qubits and measurements required, we employ two different quantum simulation apparatus to run the QAOA, which will herein be referred to as system 1  Kim2009 and system 2  Pagano2019 (see Supplementary). Both systems are based on a linear rf Paul trap where we store chains of up to $N=40$ ions and initialize the qubits in the ground state of $H_{B}$, namely the product state $|\!\!\uparrow\uparrow\cdots\uparrow\rangle_{y}\equiv|+\rangle^{\otimes N}=|\psi_{0}\rangle$, where $|\!\!\uparrow\rangle_{y}\equiv(|\!\!\uparrow\rangle_{z}+i|\!\!\downarrow\rangle_{z})/\sqrt{2}$ and $B$ is assumed to be negative. The unitary evolution under $H_{A}$ is realized by generating spin-spin interactions through spin-dependent optical dipole forces implemented by an applied laser field. This gives rise to effective long-range Ising couplings that fall off approximately as $J_{ij}\approx J_{\textrm{0}}/|i-j|^{\alpha}$ Porras2004 . The power-law exponent $\alpha\sim 1$ and the interaction strengths vary in the range $J_{0}/2\pi=$(0.3-0.57) kHz, depending on the system size and the experimental realization (see Supplementary for details). The unitary evolution under $H_{B}$ is generated by applying a global rotation around the $y$-axis of the Bloch sphere.

After each run of the algorithm, we perform a projective measurement of each spin in the $x\,(y)$ basis to measure $\langle H_{A}\rangle$ ($\langle H_{B}\rangle$) (see Fig. 1). Measurements in the $x$ and $y$ bases are carried out by performing a $\pi/2$ rotation about the $y$($x$)-axis of the Bloch sphere, illuminating the ions with resonant laser light, and collecting the $\sigma_{i}^{z}$-dependent fluorescence on a camera with site-resolved imaging. The energy is calculated by combining the measurements of the two-body correlators $\langle\sigma_{i}^{x}\sigma_{j}^{x}\rangle$ and the total magnetization along the $y$ axis $\sum_{i}\langle\sigma^{y}_{i}\rangle$, where the indices $i,j$ range from 1 to $N$. We benchmark the experimental outcome $E(\vec{\beta},\vec{\gamma})$ with the ground state $E_{gs}$ of the target Hamiltonian (see Eq. 1) calculated numerically with exact diagonalization or Density Matrix Renormalization Group (DMRG) Jaschke2018 . In order to quantify the performance of the QAOA, we use the dimensionless quantity

[TABLE]

where $E_{max}$ is the energy of the highest excited state. This choice maps the entire many-body spectrum to the $[0,1]$ interval. In the following we show that the best experimental performance $\eta^{*}$ is close to the theoretical performance $\eta_{th}$, which itself is less than unity for a finite number $p$ of QAOA layers.

Quantum Hamiltonian Optimization - We first focus on the $p=1$ optimization of the full quantum problem, where two variational parameters ($\gamma$ and $\beta$) are used to minimize the energy of the Hamiltonian (1). In this case, the time-evolved one- and two-point correlation functions can be efficiently computed dylewsky2016 ; Hadfield2018 . This leads to a general formula for the energy expectation under a state produced by the $p=1$ QAOA that is used to compute the theoretical performance of the algorithm (see Supplementary). In Fig. 2a we show an experimental exhaustive search over the parameter space $\{\gamma,\beta\}$ and compare it to the theoretical performance of the algorithm, showing good agreement for $N=20$ qubits. We also compare the performance of our algorithm as a function of $B/J_{0}$ with the expected QAOA performance $\eta_{th}$ (see Fig. 2b).

As shown in Ref. Koffel2012 , for transverse field greater than the critical value, the ground state is a low entanglement paramagnet, whereas below the critical point the ground state is an entangled superposition of anti-ferromagnetic states. We locate this critical point at $|B/J_{0}|=0.31$ for 20 qubits by computing the half-chain entanglement entropy $S_{L/2}=-\mathrm{Tr}(\rho_{L/2}\log{\rho_{L/2}})$ of the ground state numerically, where $\rho_{L/2}$ is the half-chain reduced density matrix. As shown in Fig. 2b, while the experimental performance is $\eta>94\%$ when $|B/J_{0}|$ is above the critical point, the gain relative to the initial state $|\psi_{0}\rangle$ is modest. On the other hand, below the critical point, the target state is more entangled, which allows for a larger experimental performance gain, at the expense of a reduced absolute performance. In order to quantitatively assess the gain over the finite initial state performance, we introduce a performance natural scale based on the quantity $\sigma_{\eta}(J_{0},B,N)$, namely the standard deviation around the mean performance achieved implementing a QAOA algorithm with random angles (see Supplementary for details). For $N=20$ and $B/J_{0}\sim-0.3$, $\sigma_{\eta}\sim 2\times 10^{-3}$. Our experimental performance at the critical point $\eta^{*}$ is more than 20$\sigma_{\eta}$ away from the initial state. On the other hand, the discrepancy between the ideal and experimental performance can be explained by taking into account our noise sources in the numerics (see Fig. 2c and the Combinatorial Optimization section below).

We investigate the performance of the $p=1$ QAOA algorithm as a function of the number of qubits. For each system size, we ensure that the spin-spin couplings $J_{ij}$ have the same dependence on the qubit distance $|i-j|$ by varying the trap parameters (see Supplementary). As shown in the inset of Fig. 2d, the half-chain entanglement entropy as a function of system size $N$ exhibits a peak located at $B/J_{0}\sim-0.33$, displaying the onset of the phase transition as $N$ tends to infinity. For all system sizes, we optimize the algorithm by performing a scan of the interaction angle $\gamma$ and applying discrete variations of the mixing angle $\beta$ around the optimal value predicted by the theory. In Fig. 2d we compare the optimal experimental and theoretical performances $\eta$ for different system sizes from 20 up to 40 qubits for fixed $B/J_{0}\sim-0.3$. We observe experimentally that the QAOA yields a similar performance as a function of number of qubits even if the algorithm runtime stays approximately constant as the number of qubits increases. Numerically, we found that the performance $\eta$ scales polynomially with $N$ and with the number of layers $p$ (see Supplementary). Assuming extrapolation to higher numbers of qubits holds, this scaling, combined with a polynomial-time search heuristic, suggests that for any desired energy threshold $\epsilon$, our approach allows us to approximate the energy to a degree $\eta>1-\epsilon$ in time and number of layers that scale as $\text{poly}(N,1/\epsilon)$.

We experimentally perform a search for the optimal $p=2$ QAOA performance using 20 qubits. Unlike the $p=1$ case, there is no known analytic formula to efficiently compute the energy. However, exploiting relationships between optimal angles as a function of increasing $p$, we use a bootstrapping heuristic (see Supplementary for details) that allows the experiment to identify a set of optimal angles faster than a global parameter search. The bootstrapping heuristic computes a guess for optimal angles at $p$ given optimal angles at lower $p$. A local optimizer, such as the greedy gradient descent described below, is then needed to take this guess to the true optimum. Our new heuristic method allows us to find variational parameters in time that scales polynomially with the number of layers and sublinearly in the number of qubits (when used in conjunction with the quantum device).

We start from the optimal guess and perform a fine scan of $\gamma_{2}$, while varying $\gamma_{1},\beta_{1}$ and $\beta_{2}$ in larger steps. The result is shown in Fig. 2d, where we plot the performances $\eta$ as a function of $\gamma_{2}$ for every set of parameters used in the experiment. Fig. 2d shows also a colour plot of all the optimal energies found as a function of the other three parameters $\gamma_{1},\beta_{1}$ and $\beta_{2}$. The $p=2$ QAOA performance with 20 qubits $\eta^{*}=(93.9\pm 0.3)\%$ is in agreement with the $p=1$ performance in system 2, taken with the same parameters (see Fig. 2c). This indicates that decoherence and bit-flip errors (see Supplementary) accumulated during longer evolution times are already balancing out the $2\%$ expected performance gain of one additional optimization layer.

As a brute force approach is inefficient, we implement a closed-loop QAOA by interfacing the analog trapped-ion quantum simulator with a greedy gradient-descent algorithm to optimize the measured energy. In the $p=1$ QAOA, we can visualize the optimization trajectory on the theoretical performance surface as shown in Fig. 3. Starting from a guess $(\beta^{(0)},\gamma^{(0)})$, we measure the approximate local gradient by performing the energy measurements in two orthogonal directions $\beta^{(0)}+\delta\beta$ and $\gamma^{(0)}+\delta\gamma$ to compute the new guess $(\beta^{(1)},\gamma^{(1)})$, where we measure the new energy on the quantum simulator. As shown in Fig. 3, the algorithm converges after about 10 iterations. Compared to an exhaustive search, the gradient descent uses fewer queries to the quantum simulator and is therefore more robust to slow drifts in the experimental system. For this reason, we are able to achieve a better performance compared to the exhaustive search method.

Combinatorial Optimization - We further explore the performance of the trapped-ion system by investigating the combinatorial optimization of the classical Hamiltonian $H_{A}$ (see Eq. (1) with $B=0$) approximately sampling the output of the $p=1$ QAOA, using high-fidelity, single-shot measurement of all the qubits. It has been proven, under reasonable complexity-theoretic assumptions, that no classical algorithm can efficiently sample exactly from a sufficiently general class of $p=1$ QAOA circuits harrow2016qaoa . Recent results Bremner2016 ; Bouland2019 suggest that this could also hold in the case of approximate sampling (see Supplementary). In this case, by measuring in the $x$ basis, it is possible to sample the probability distribution of all the $2^{N}$ eigenstates $|x_{i}\rangle$ of the Hamiltonian $H_{A}$. We performed the experiment with 12 qubits so that we can both compute the expected QAOA theoretical output and also experimentally over-sample the Hilbert space of all the possible $2^{12}=4096$ possible outcomes. In Fig. 4a we show on a log scale the QAOA eigenstates probability distribution using the optimal variational parameters ${\beta^{*},\gamma^{*}}$ and compare the experimental eigenstate histogram with the exact diagonalization prediction of the QAOA output state, sorting the eigenstates according to their energies.

However, sampling from the full QAOA output distribution is a daunting task, since the experimental outcome is extremely sensitive to fluctuations in the Hamiltonian parameters and to experimental errors caused by detection and phonon-assisted bit-flip events and unwanted effective magnetic fields along the $z$ direction of the Bloch sphere caused by uncompensated light shift (see also Supplementary). Given our measured experimental parameters, we can calculate the effect of these errors on the quantum evolution, resulting in a good agreement with the experimental outcome, as shown in Fig. 4a.

Another useful way to compare numerics and experimental data is to implement the coarse-graining procedure of the Hilbert space proposed in Ref. Wang2015 . After sorting in decreasing order the observed states according to their experimental probability, we iteratively group the states into “bubbles” of Hamming distance $L$ around the most probable state, producing a coarse-grained dataset. We then apply the same coarse-graining to the theoretical probability distribution and plot the comparison in Fig. 4b. In this procedure the Hamming distance radius is varied to ensure that each bubble contains a comparable number of experimental shots, leading to bubbles of average Hamming distance $\bar{L}=2.5$. In order to quantitatively compare the coarse-grained experiment and the theory, we use two different metrics, namely the total variation distance (TVD) and the Kullback-Leibler divergence ($\mathrm{D_{K\rm{-}L}}$), defined as:

[TABLE]

where $p_{i}(q_{i})$ is the experimental (theoretical) probability of observing the $i$-th outcome. As shown in Fig. 4c, when the system is in the initial state, it is closer to a uniform probability distribution since $|\psi_{0}\rangle$ is an equal superposition of all the eigenstates of $H_{A}$. On the other hand, as the $\gamma$ parameter is scanned, we observe a net decrease of both TVD and $\mathrm{D_{K-L}}$ between the experiment and the numerical minimum, in agreement with the decrease in energy, computed by measuring one and two-body correlators.

The variational quantum algorithm reported here, with up to 40 trapped-ion qubits, is the largest ever realized on a quantum device. We approximate the ground state energy of a non-trivial quantum Hamiltonian showing almost constant time scaling with the system size. Single-shot high-efficiency qubit measurements in different bases give access to the full distribution of bit-strings that is difficult or potentially impossible to model classically. With the addition of individual control over the interactions between qubits as well as improvements to fidelity and system size, the variational quantum-classical hybrid approach can be employed in this experimental platform to give insight into quantum chemistry Kandala2017 ; Hempel2018 ; Nam2019 and hard optimization problems Rigetti2017 , such as Max-SAT or exact cover Farhi2001 , or be used for the production of highly entangled states of metrological interest Giovannetti2011 .

.1 Acknowledgements

We acknowledge illuminating discussions with L. Duan, Y. Wu, B. Fefferman and S. Wang. The DMRG simulations were performed using Open Source Matrix Product States Jaschke2018 . This work is supported by the ARO and AFOSR QIS and Atomic and Molecular Physics Programs, the AFOSR MURIs on Quantum Measurement/Verification and Quantum Interactive Protocols, the IARPA LogiQ program, the ARO MURI on Modular Quantum Systems, the ARL Center for Distributed Quantum Information, the NSF QIS program and the NSF Physics Frontier Center at JQI. A.B., A. Des., F.L., and A.V.G. were mainly supported by the DoE BES QIS program (award No. DE-SC0019449), with additional support by the NSF PFCQC program, DoE ASCR Quantum Testbed Pathfinder program (award No. DE-SC0019040), ARO MURI, AFOSR, NSF PFC at JQI, and ARL CDQI. S. J. and A.B. were supported by the DOE HEP QuantISED Program, and Quantum Algorithms Teams (QOALAs) programs.

.2 Author Contributions

G.P., A.K., P.B., W.L.T., H.B.K., K.S.C., A.D., P.W.H. and C.M. all contributed to experimental design, construction, data collection and analysis. A.B., L.T.B., F.L., A.D., C.B., S.J. and A.V.G. contributed to the theory for the experiment. All authors contributed to this manuscript.

.3 Author Information

The authors declare competing financial interests: details are available in the online version of the paper. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to G.P. (pagano@umd.edu)

I Supplementary material

I.1 Quantum Approximate Optimization Algorithm (QAOA)

The QAOA is an approximate optimization algorithm first introduced in 2014 by Farhi et al. Farhi2014 , and has since enjoyed growing interest. The QAOA uses alternating evolutions under two non-commuting operators, typically a problem (or cost) Hamiltonian $H_{A}$ that encodes the cost function on the diagonal in (say) the $\sigma^{x}$ basis, and a transverse term $H_{B}=-\sum\limits_{i=0}^{N}\sigma^{y}_{i}$ that generates transitions between bit strings, such that the initial state $|+\rangle_{y}^{\otimes N}$ evolves into an approximate ground state of $H_{A}$.

Practically, the most valuable feature of the QAOA seems to be its “learnability” via a classical outer loop optimizer, where the discovery of the evolution angles in the optimal QAOA schedule is achieved via the discovery of structure in the angle sequences Zhou2018 ; brady2019 ; crooks2018 . These patterns are seen quite generally across local Hamiltonian problems, and while steps towards a theory describing optimal QAOA sequences have been taken brady2019 , several questions surrounding it remain open. Regardless, the structure in optimal QAOA schedules may be harnessed to implement approximate state preparation in a scalable manner and with a low overhead on quantum resources. We present a new heuristic method that helps achieves this goal.

First, we discuss how to discover optimal QAOA1 schedules, i.e., QAOA schedules for $p=1$.

I.1.1 QAOA, $p=1$

Despite its apparent simplicity, the $p=1$ QAOA (or QAOA1) can be a powerful state preparation ansatz. For example, hardness-of-sampling results are known for QAOA1 circuits harrow2016qaoa , closely mirroring the hardness of sampling from instantaneous quantum polynomial (IQP) circuits (see next section for details). Furthermore, it is known that the performance of the QAOA1 for certain combinatorial optimization problems can be competitive with the best classical algorithms for the same problems Farhi2014maxE3lin2 . Another desirable feature of the QAOA1 for local spin Hamiltonians is the tractability of computing energy expectation values, as observed in Farhi2014 . A very similar result has also been known in the setting of quantum dynamics dylewsky2016 ; Hadfield2018 . For a two-local transverse field spin Hamiltonian as in Eq. (1) in the main text, this leads to a formula for the energy expectation under a state produced by the QAOA1, starting from the product state $|+\rangle^{\otimes N}$. These formulas are applicable to many cases of interest in quantum state preparation and optimization. Importantly, the time complexity to compute the formula is $O(N^{3})$ in the worst case, making it tractable to optimize the QAOA1 protocols for large spin chains.

I.1.2 QAOA, $p>1$

The general analytical formula for $p=1$ does not extend to the case where we apply the QAOA for more than one layer. Here, we must turn to classical numerical methods to find the optimal QAOA angles $\beta_{i},\gamma_{i}$ for each layer $i$. For $p$ layers, this is an optimization on a $2p$-dimensional space that grows exponentially with the depth of the circuit. However, numerics done here and in crooks2018 ; Zhou2018 have identified the existence of minima that exhibit patterns in the optimal QAOA angles, namely that the angles, when plotted as a function of their index $i$, form smooth curves for any $p$. While this observation points to a deeper theoretical mechanism at play, it does not directly simplify the optimization problem, since we must still search over all approximately smooth sequences of the angles. Zhou et al. Zhou2018 have exploited the smoothness of the functions by carrying out searches in the Fourier domain. Here, we follow a different route that arises from some novel observations of these family of minima.

For each $p$, denote the special optimal angles by $\left\{\left(\boldsymbol{\beta}^{*(p)},\boldsymbol{\gamma}^{*(p)}\right)\right\}_{p}$, which we can also think of as a pair of angle curves (as a function of step index $i$). As $p$ is varied, we may think of these minima as a family. We numerically find that this family exhibits the following desirable features (for $p$ sufficiently large):

The angles are non-negative, small and bounded. 2. 2.

For $p$ sufficiently large, the two angle sequences ${\beta}^{*(p)}$ and ${\gamma}^{*(p)}$ are approximately smooth. 3. 3.

The angle sequence ${\beta}^{*(p)}$ (and correspondingly, ${\gamma}^{*(p)}$) when viewed as a function on the normalized time parameter $s_{i}=\frac{i-1}{p-1}$, is convergent in the parameter $p$. In other words, as $p$ is increased, the angle sequences ${\beta}^{*(p)}$ and ${\gamma}^{*(p)}$ approach a smooth, asymptotic curve (See Fig. S1.) 4. 4.

The energy expectation $E(\beta^{*(p)},\gamma^{*(p)})$ approaches the global minimum as $p\rightarrow\infty$, and hence this family is asymptotically optimal.

The significance of the first point is that in experimental settings, large evolution times are infeasible to implement due to decoherence, so these minima correspond to practicable QAOA protocols. The third and fourth points suggest an inductive algorithm where a locally optimal schedule for a given $p$ may be discovered using the optimal schedule for $p-1$ as a prior.

Point 3 in the above list is a novel observation that allows us to construct a heuristic that is efficiently scalable for large $p$. The main idea behind this construction is that the minimal angle curves for a larger $p$ may be guessed from the optimal curve of a smaller $p^{\prime}<p$ by interpolation.

Using the above points, we use a bootstrapping algorithm to find the optimal angle sequences, ${\beta}^{*(p)}$ and ${\gamma}^{*(p)}$, for a given $p$, as described below. Let $q=1,\ldots,p$ denote an intermediate angle index. Then:

For $q=1$, use an analytic formula to find ${\beta}^{*(1)}$ and ${\gamma}^{*(1)}$. 2. 2.

For $q=2$, choose an initial guess of ${\beta}^{(2)}=\left({\beta}^{*(1)},{\beta}^{*(1)}-0.2\right)$ and ${\gamma}^{(2)}=\left({\gamma}^{*(1)},{\gamma}^{*(1)}+0.2\right)$. 3. 3.

Perform a local optimization of $\beta^{(2)}$ and $\gamma^{(2)}$ in order to find $\beta^{*(2)}$ and $\gamma^{*(2)}$. 4. 4.

Repeat the next steps (5-7) for $q=3,\ldots,p$. 5. 5.

Create interpolating functions through the angle sequences, $\beta^{*(q-1)}$ and $\gamma^{*(q-1)}$, using the normalized time $s_{i}=\frac{i-1}{q-2}$ as the independent parameter (we use a linear interpolation for $q=3$ and cubic for $q>3$). 6. 6.

Choose the initial guesses for $\beta^{(q)}$ and $\gamma^{(q)}$ by sampling the interpolating function from (5) at evenly spaced points separated by a normalized time distance of $\Delta s=1/(q-1)$. 7. 7.

Perform a local optimization of $\beta^{(q)}$ and $\gamma^{(q)}$ in order to find $\beta^{*(q)}$ and $\gamma^{*(q)}$.

The resulting angles $\beta^{*(p)}$ and $\gamma^{*(p)}$ should be at least a good local minimum of the energy expectation value and approaches the global minimum as $p\to\infty$.

The $q=2$ interpolation in step 2 is based on our observation that the $\beta$ angles tend to curve down at the end and the $\gamma$ angles tend to curve up.

An important feature of our algorithm is that its asymptotic runtime is expected to be efficient in $p$. This feature is predicated on the previous result that the angle curves are generally convergent as $p$ tends to infinity. The argument proceeds as follows: if we assume a maximal deviation of the initial guess for layer $q$ to be $\epsilon_{q}\geq 0$, then the total $l_{2}$-norm distance between the initial guess and the optimized curve is no greater than $\epsilon_{q}\sqrt{q}$, by the Cauchy-Schwarz inequality. Therefore, the local search algorithm is confined to a ball of radius at most $\epsilon_{q}\sqrt{q}$, and for a fixed error tolerance, the convergence time for a standard local optimizer is $O(\epsilon_{q}^{2}q)$. Summing over convergence times for all from $q=1,\ldots,p$, we have

[TABLE]

The last inequality above comes about as follows: while the summand depends on the convergence rate of the sequence $\left\{\epsilon_{q}\right\}_{q=1}^{p}$, it is upper bounded by $O(q)$ for a converging set of paths and an initial error $\epsilon_{1}$ of order 1. The latter is true since our angle search domain is bounded and independent of $N$. Therefore, the sum is no greater than $O(p^{2})$. In practice, even faster runtimes are possible. Therefore, the bootstrap algorithm exploits the structure of the special minima and provides a scalable route to multi-step QAOA for the long-range TFIM. In fact, as discussed in the supplement and in brady2019 , there is mounting numerical evidence that the path approach applies across a very general variety of models on discrete as well as continuous systems.

I.1.3 Convergence in $N$

In the previous sections, we introduced a bootstrap algorithm that is asymptotically efficient in the number of layers $p$. However, in order to be fully scalable the algorithm must also be scalable in the system size $N$. This may not be possible in general (say for random spin models), as the optimized angles for a particular small system may have no bearing on the angles for a larger system. However, for the long-range TFIM, and indeed any translationally-invariant model with a well-defined notion of metric and dimension arising from the functional form of the coupling coefficients $J_{ij}$, it is reasonable to expect that the optimized angles depend on system size in a predictable way. This is indeed the case for the long-range TFIM. There, it can be seen that the angle curves for varying $N$ appear similar in shape. Usefully, the curves also appear to be convergent to an idealized curve for a hypothetical continuous, long-range spin chain. Once again, this feature suggests that the optimized QAOA angle curves for small systems may be used as initial guesses for larger systems within the same Hamiltonian family.

While it is not clear (due to numerical limitations) how fast the curves converge, we argue that the rate should be weakly dependent (or independent) of the system size $N$. For a given coupling function (such as inverse power-law) that decays as a function of distance, we define a characteristic length scale, which may be called the skin depth $\delta$, that is the number of sites from the boundary that the coupling is a factor of $e$ smaller than the nearest-neighbour value. In other words, we define $\delta$ such that $J_{i,i+\delta}\sim J_{i,i+1}/e$. Clearly, $\delta$ is independent of the system size $N$ and depends only on the parameters of the coupling function. For instance, for the long-range TFIM, $\delta\sim e^{1/\alpha}$. As $N$ tends to infinity, the fractional skin depth $\delta/N$ then “falls away” and becomes vanishing with respect to the bulk region of the chain. Now, we make the assumption that any deviations in the optimal QAOA schedules from $N$ to $N+1$ arise from change in the fractional skin depth, which is reasonable for a translationally invariant model. The incremental change in the fractional skin depth from $N$ to $N+1$ is $\delta/N-\delta/(N+1)\sim O(1/N^{2})$. Therefore, if the change in the optimal QAOA curves $\epsilon_{N}$ (in, say, $l_{1}$-norm distance) is a smooth function of the the fractional skin depth, then we expect it to vary as $\epsilon_{N}\sim 1/\text{poly}(N)$. Therefore, the total running time of a bootstrap from small system sizes to a given size $N$ should be $O\left(\sum\limits_{k=1}^{N}1/\text{poly}(N)\right)$ which is sub-linear in $N$. Combining this observation with the convergence in $p$, we see that for a given Hamiltonian family, optimized QAOA angle curves for small $p$ may be used as a rubric for the optimization for longer circuit depths. Furthermore, if the Hamiltonian is translationally-invariant with decaying interactions, the optimized QAOA schedules are expected to scale with $N$ as well. Therefore, the state preparation procedure under the QAOA for such a Hamiltonian family is scalable in circuit “volume”, for a wide range of Hamiltonian parameters (Fig. S2). This is our main theoretical contribution in this work.

I.1.4 Scaling of $\eta$ in $p,N$

Our performance parameter $\eta$, defined as

[TABLE]

measures how close (in energy) the prepared state is to the ground state of the system. As described in previous sections, the optimal angle curves for QAOA appear to converge to a smooth, hypothetical curve, as a function of $p$ as well as $N$. We show that under the assumption that such a curve exists, there is a fast heuristic for finding optimal angles for any finite $p$ that is time-efficient in $p$ and the number of spins $N$ (when used in conjunction with the quantum device). In this section, we show that not only is the search efficient, but the quality of the optimum is numerically seen to improve with $p,N$ as well.

In Fig. S3, we show the result of the numerical study. We chose as the target Hamiltonian an idealized transverse field Ising model with inverse power-law couplings, with the power $\alpha=1.1$ chosen to closely mimic the experimental Hamiltonian. The number of spins was varied from $N=8$ to $20$. Via DMRG, the critical value of the transverse field for a finite chain can be located by maximizing the von Neumann entropy at half-cut. This was done independently for each value of $N$. Then, using our heuristic, we located the optimal angle curve, and computed $\eta$ for the final state prepared using this angle sequence, for each $N$. The plot shows the trend of $1/(1-\eta)$ with $N$, for a range of $p=0,1,2,3,4,5$, with [math] corresponding to a trivial protocol where the initial state is returned. While the number of spins could not be extended beyond $20$ due to computational limitations, the trend is clear. We see that $1/(1-\eta)$ grows linearly with $N$ and $p$ (inset). While the linear trend in $N$ is encouraging, we similarly expect the inverse spectral gap (and indeed, the density of low-lying states) to increase with $N$. Empirically for the target Hamiltonian, we observe a gap scaling of $\sim 1/N^{2}$. Assuming the density of low-lying states scales similarly, this suggests that the squared overlap with the ground state should fall off with $N$. Numerics confirm this expectation and indicate a scaling of the squared overlap of $|\langle\psi|\psi_{0}\rangle|^{2}\sim p/N$.

The linear scaling with $p$ for both the energy and fidelity metric, combined with a polynomial-time search heuristic, suggests that for any desired energy (or probability) threshold $\epsilon$, our approach allows us to approximate the state to within $1-\epsilon$ (in energy or fidelity) in time and number of layers that scale as $\text{poly}(N,1/\epsilon)$.

I.1.5 Characteristic scale for $\eta$

The figure of merit $\eta$ characterizes how close the final state is to the ground state of the system. At $\eta=0$, the system is in the highest excited configuration, while $\eta=1$ corresponds to a perfectly prepared ground state. QAOA, starting from the initial state $|+\rangle^{\otimes n}$, gives a state with figure of merit $\eta\in[0,1]$, from the initial value of $\eta_{0}$. The difference between the final $\eta$ and $\eta_{0}$ indicate the success of our QAOA protocol.

While $\eta$ is normalized to the range $[0,1]$, differences in $\eta$ are still somewhat arbitrary. In long-range Ising models with a transverse field $\eta_{0}$ is not [math] but typically greater than $0.5$, making the difference in $\eta$ an unsatisfactory metric of success. Therefore, in addition to the initial and final $\eta$, we must provide a characteristic scale for $\eta$ that quantifies the typical deviation from $\eta_{0}$. A natural choice is the standard deviation of $\eta$ for QAOA with random angles.

For QAOA1 with evolution angles $\beta,\gamma$, it is possible to estimate the standard deviation analytically as a function of the underlying model parameters $B$ and $J_{0}$ and on the number of qubits $N$. This derives from the analytical formula for the energy expectation $E(\beta,\gamma)$ which can be stated as follows:

[TABLE]

where

[TABLE]

where the Hamiltonian has long-range power law couplings $J_{ij}\sim\frac{1}{|i-j|^{\alpha}}$ (with $J_{ii}=0$), and a transverse field of strength $B$. Then, our goal is to compute the standard deviation (normalized by the spectral bandwidth $\Delta:=E_{max}-E_{gs}$),

[TABLE]

which gives us the characteristic scale for $\eta$. We define the average $\langle\cdot\rangle_{\beta,\gamma}$ as

[TABLE]

In the limit, the average is precisely the constant term of the Fourier transform of $f$. Since the function is a sum of trigonometric monomials, its moments over the angle variables $\beta,\gamma$ can be computed analytically term by term. We will need the following properties of the coupling function:

(Symmetry) Since the inverse power law only depends on distance between nodes, we have $J_{ij}=J_{(2j-i)j}$ In other words, the inverse power-law is symmetric under a lateral flip (or “mirroring”). We assume a finite, open chain, and therefore couplings $J_{ij}$ with $|j-i|>N-j$ do not have an image under mirroring. 2. 2.

(Incommensurateness) The coupling strengths $J_{ij}$ are, in general, mutually indivisible irrational numbers whose sums and differences are also irrational and mutually distinct, e.g. for $i\neq j,k\neq l$, $J_{ik}\pm J_{jk}\neq J_{il}\pm J_{jl}$ (with a very small set of exceptions due to, say, symmetry).

The mean $\langle E\rangle_{\beta,\gamma}$ consists of three parts corresponding to the terms $E_{I},E_{II},E_{III}$. Performing the $\beta$ integral first, we see that $\langle E_{II}\rangle_{\beta,\gamma}=0$. Next, we may argue that in products of the form $\prod\limits_{k}\cos(2\gamma J_{ik})$, the cosine factors are of degree one if they have no mirror images, and degree two otherwise. The only way to have a non-zero expectation is if all terms are systematically paired up by mirroring, so that the overall product is quadratic in a product of cosines. For the summand in $E_{I}$, this can only happen if $N$ is odd and $i$ is exactly at the center of the chain, in which case the average is $B/2^{(N-1)/2}$. When $N$ is even, the mean is 0. Finally, for general $i,j$ the last term is zero by property 2, since the cosines are generically incommensurate and therefore barring very few exceptions, most phases do not cancel out. However, in the special case that $i,j$ are mirror images, i.e. $i=N-j$, we have perfectly paired terms when $N$ is even (and one unpaired term at $k=\lfloor{N/2}\rfloor$ when $N$ is odd). Counting all occurrences of this case, the mean is approximately $\frac{1}{2^{N/2+1}}\sum\limits_{i=1}^{N}J_{i(N-i)}\lesssim NJ_{0}/2^{N/2}$ where $J_{0}$ is the nearest-neighbor coupling in the chain. Note that asymptotically in $N$, $\langle E\rangle_{\beta,\gamma}\sim O(N/2^{N/2})$ which approaches [math] in the infinite $N$ limit.

Next, we estimate the term $\langle E\rangle_{\beta,\gamma}^{2}$. By the orthogonality of trigonometric polynomials in $\beta$, we first have that $\langle E\rangle_{\beta,\gamma}^{2}=\langle E_{I}\rangle_{\beta,\gamma}^{2}+\langle E_{II}\rangle_{\beta,\gamma}^{2}+\langle E_{III}\rangle_{\beta,\gamma}^{2}$. Therefore, we estimate each term separately. As before, we require that the cosines pair up so that their phases can cancel. First, we have

[TABLE]

Each summand is a product of $2N$ cosines, and only survives averaging if every cosine is paired. This happens exactly when either $i=j$ or $i=N-j$ (There is also the “disconnected” contribution that cancels with the mean). In each case, the squared cosines give a factor of $1/2$ from averaging. Moreover, using mirror symmetry we can have fourth powers of some of the cosines, which give a factor $3/8$ from averaging. In all, the mean (minus the disconnected part) is no greater than

[TABLE]

A similar reasoning for $E_{II},E_{III}$ give us the following estimates:

[TABLE]

Finally, this gives

[TABLE]

Therefore, we see that the standard deviation $\sigma_{\eta}=\sigma_{E}/\Delta\sim\frac{\sqrt{8B^{2}+J_{0}^{2}}}{\Delta}\cdot N^{1/4}\left(3/8\right)^{N/4}$, which is exponentially suppressed for large $N$. For $N=20$ ions, we have $N^{1/4}\cdot\left(3/8\right)^{N/4}\sim 0.02$. While this is already small, the normalization $\frac{\sqrt{8B^{2}+J_{0}^{2}}}{\Delta}$ will have an additional linear $N$ factor in the denominator, making the scale for $\eta$ about $0.002$. Therefore, a typical final QAOA performance of $\eta\gtrsim 0.95$ is several standard deviations above a typical $\eta_{0}\sim 0.85$.

I.2 Evidence for hardness of sampling from general QAOA circuits

In this section we expand upon previous work Farhi2014 that gives evidence for exact sampling hardness of QAOA circuits, using the techniques of Refs. Bremner2016 ; Bouland2019 to give evidence for hardness of approximate sampling. First we relabel the bases $Y\rightarrow X\rightarrow Z$ so that the $p=1$ experiment is equivalent to preparing a state $|\psi_{0}\rangle=|\!\!\uparrow\rangle^{\otimes N}_{x}$, evolving under a Hamiltonian $H_{z}$ diagonal in the computational basis, followed by a uniform rotation $\tilde{H}=e^{-i\beta\sum_{i}\sigma^{x}_{i}}$ and measurement in the computational basis. Following Ref. Farhi2014 , it suffices to consider QAOA circuits with $\beta=\pi/4$. The output state is $\tilde{H}^{\otimes N}e^{-i\gamma H_{z}}H^{\otimes N}|0^{N}\rangle$ for some cost function $C$ diagonal in the computational basis.

I.2.1 Generalized gap of a function

The main idea behind proving exact sampling hardness is to examine a particular output amplitude, say the amplitude of the $|0^{N}\rangle$ basis state. In Ref. Bremner2016 , the output state after a so-called IQP circuit (which only differs from the one here in that the final rotation is a global Hadamard $H^{\otimes N}$ instead of $\tilde{H}^{\otimes N}$) has an amplitude proportional to a quantity known as the gap of a Boolean function, $\mathrm{gap}(f)=\sum_{x:f(x)=0}1-\sum_{x:f(x)=1}1$, the difference in the number of inputs that map to 1 and the number of inputs that map to 0 under $f$. Finding the gap of a general function is a $\mathsf{GapP}$-complete problem. This is a very hard problem since the class $\mathsf{GapP}$ includes $\mathsf{\#P}$, which in turn includes the whole of $\mathsf{NP}$. The authors of Ref. Bremner2016 prove that the gap of a degree-3 polynomial over $\mathbb{Z}_{2}$, $f$, may be expressed as an output amplitude of an IQP circuit. They also show that the finding the gap of such functions $f$ is still $\mathsf{GapP}$-complete. Following Ref. Bremner2016 , we examine the $|0^{N}\rangle$ output amplitude of a QAOA state:

[TABLE]

where now we define the function $\tilde{f}$ to have the range $\mathbb{Z}_{4}$ and the Hamiltonian $H_{z}$ satisfies $e^{-i\gamma H_{z}}|x\rangle=i^{\tilde{f}(x)}|x\rangle$ for a computational basis state $|x\rangle$. The output amplitude is thus proportional to a ‘generalized gap’ $\mathrm{ggap}(f):=\sum_{x:f(x)=0}1+i\sum_{x:f(x)=1}1+i^{2}\sum_{x:f(x)=2}1+i^{3}\sum_{x:f(x)=3}1$ of a function $f(x)=\tilde{f}(x)+\mathrm{wt}(x)$, where $\mathrm{wt}(x)$ is the Hamming weight of $x$. This modified function $f(x)$ is also a degree-3 polynomial over $\mathbb{Z}_{4}$. Note that this restriction to degree-3 comes from the fact that the gates $Z$, $CZ$ and $CCZ$ are universal for classical computation (indeed, the Toffoli alone is universal for classical computation) and there is a natural degree-3 polynomial coming from this construction. The quantity we have defined, $\mathrm{ggap}(f)$, can be easily shown to be $\mathsf{GapP}$-hard to compute, by reducing $\mathrm{gap}$ to $\mathrm{ggap}$. This suffices for exact sampling hardness assuming the polynomial hierarchy ($\mathsf{PH}$) does not collapse.

I.2.2 Approximate sampling hardness

For approximate sampling hardness, we need two other properties, namely anti-concentration and a worst-to-average case reduction. Anti-concentration of a circuit roughly says that the output probability is sufficiently spread out among all possible outcomes so that not many output probabilities are too small. We choose a random family of QAOA circuits by choosing $H_{z}$ such that the function $f(x)$ is a degree-3 polynomial $\sum_{i,j,k}a_{i,j,k}x_{i}x_{j}x_{k}+\sum_{i,j}b_{i,j}x_{i}x_{j}+\sum_{i}c_{i}x_{i}$ with uniformly random weights $b_{i,j}$ and $c_{i}\in\mathbb{Z}_{4}$. Anti-concentration then follows from the Paley-Zygmund inequality and Lemma 4 of the Supplemental Material of Ref. Bremner2016 (with $r=s=4$).

Finally, we need to show that the problem of approximating the generalized gap is average-case hard. Currently, no scheme for quantum computational supremacy has achieved this, and the best known result in this direction is in Ref. Bouland2019 , where the authors show a worst-to-average case reduction for the problem of exactly computing an output probability of a random quantum circuit. The authors remark that their techniques may be extended to any distribution parametrized by a continuous variable. In principle, we have such a parameter $\gamma$ available here, which continuously changes the parameters $b_{i,j}$ and $c_{i}$. However, we have only shown anti-concentration when the weights $b_{i,j}$ and $c_{i}$ are chosen from a finite set. It remains to be seen whether one can have the property of anti-concentration and average-case hardness holding at the same time for some specific QAOA output distribution.

I.3 Trapped-ion experimental systems

In this work two quantum simulators have been used, referred to as system 1 and 2. System 1 Kim2009 is a room-temperature ion-trap apparatus, consisting of a 3-layer linear Paul trap with transverse center-of-mass (COM) motional frequency $\nu_{\text{\rm COM}}=4.7$ MHz and axial center-of-mass frequencies $\nu_{x}$ ranging from $0.39$ to $0.6$ MHz depending on the number of trapped ions. In this system Langevin collisions with the residual background gas in the ultra high vacuum (UHV) apparatus are the main limitation to ion chain lifetime Wineland1998 . These events can melt the crystal and eject the ions from the trap because of rf-heating or other mechanisms.

System 2 Pagano2019 is a cryogenic ion-trap apparatus based on a linear blade trap with four segmented gold coated electrodes. The trap is held at 6.5 K in a closed cycle cryostat, where differential cryo-pumping reduces the background pressure at low $10^{-12}$ Torr level, which allows for long storage times of large ion chains. For this reason system 2 has been used to perform the QAOA with a large number of qubits (Fig. 2b) or when a large number of measurements was required (Fig. 4). The two transverse trap frequencies are $\nu^{y}_{\rm COM}=4.4$ MHz and $\nu^{z}_{\rm COM}=4.26$ MHz, and the axial frequency ranges from $0.27$ to $0.46$ MHz.

I.3.1 State preparation

The qubit is initialized by applying resonant 369.5 nm light for about 20 $\mu$s to optically pump into the $|\!\!\downarrow\rangle_{z}$ state. To perform global rotations in the Bloch sphere, we apply two far-detuned, non-copropagating Raman beams whose beatnote is tuned to the hyperfine splitting $\nu_{0}=12.642821$ GHz of the clock states ${}^{2}S_{1/2}|F=0,m_{F}=0\rangle$ and ${}^{2}S_{1/2}|F=1,m_{F}=0\rangle$ encoding the qubit Olmschenk2007 . State preparation in our implementation of the QAOA requires qubit initialization in the $|\!\!\downarrow\rangle_{z}$ state by optically pumping the ions and then a global rotation into the $|\!\!\uparrow\rangle_{y}$ state using stimulated Raman transitions. We detect the state of each ion at the end of each experimental sequence using state-dependent fluorescence, with single site resolution. In order to improve the accuracy of global qubit rotations, we employ a composite pulse sequence based on the dynamical decoupling BB1 scheme Brown2004 . This allows us to compensate for inhomogeneity due to the Raman beam’s Gaussian profile and achieve nearly $99\%$ state preparation fidelity. The BB1 four pulse sequence is:

[TABLE]

where after the first $\pi/2$ rotation $e^{-i\frac{\pi}{4}\sigma_{i}^{x}}$, three additional rotations are applied: a $\pi$-pulse along an angle $\theta=\textrm{cos}^{-1}(-1/16)=93.6\degree$, a $2\pi$-pulse along $3\theta$, and another $\pi$-pulse along $\theta$. The axes of these additional rotations are in the $x$-$y$ plane of the Bloch sphere with the specified angle referenced to the $x$-axis.

I.3.2 Generating the Ising Hamiltonian

We generate spin-spin interactions by employing a spin dependent force with a pair of non-copropagating 355 nm Raman beams, with a wavevector difference $\Delta k$ aligned along the transverse motional modes of the ion chain. The two off-resonant Raman beams are controlled using acousto-optic modulators which generate two interference beatnotes at frequencies $\nu_{0}\pm\mu$ in the Mølmer-Sørensen configuration Molmer1999 . In the Lamb-Dicke regime, the laser-ion interaction gives rise to the effective spin-spin Hamiltonian in Eq. (1) in the main text, where the coupling between the $i$-th and $j$-th ion is:

[TABLE]

Here $\Omega$ is the Rabi frequency, $\nu_{R}=h\Delta k^{2}/(8\pi^{2}M)$ is the recoil frequency, $\nu_{m}$ is the frequency of the $m$-th normal mode, $b_{im}$ is the eigenvector matrix element for the $i$-th ion’s participation to the $m$-th normal mode $(\sum_{i}|b_{im}|^{2}=\sum_{m}|b_{im}|^{2}=1)$ James1998 , and $M$ is the mass of a single ion.

Differently from system 1, where the wavevector difference $\Delta k$ of the Raman beams is aligned along one of the principal axes of the trap, in system 2 the spin-spin interaction stems from the off-resonant coupling to both families of transverse normal modes. Eq. (20) is then generalized to:

[TABLE]

where $\nu^{\ell}_{R}$ is the recoil frequency given by the projection of the Raman wavevector $\Delta k$ along the two transverse principal axes of the trap $\ell=y,z$. We infer an angle $\vartheta\sim 40^{o}$ between $\Delta k$ and the $z$ principal axis (see inset in Fig. S4a) from the ratio between the resonant spin-phonon couplings to the two transverse COM modes. Before every experiment, we perform Raman sideband cooling on both the COM and the two nearby tilt modes for both transverse mode families.

As we scale up the number of qubits (see Fig. 2c in the main text), we vary the axial confinement in order to maintain a self-similar functional form of the spin-spin interaction (see Fig. S4b). For the data in Fig. 2c in the main text, we set the detuning to $\delta=\mu-\omega^{y}_{\rm COM}=2\pi\times 45$ kHz and the axial frequency to $\nu_{x}=0.46,0.37,0.36,0.31,0.27$ MHz, for $N=20,25,30,35,40$ respectively. For the data in Fig. 4 in the main text, the detuning is $\delta/2\pi=45$ kHz and the $\nu_{x}=0.54$ MHz.

I.3.3 Fitting Ising Couplings to Analytic Form

By directly measuring trap parameters and spin-phonon couplings, we can calculate the spin-spin interaction matrix $J_{ij}$ with Eqs. (20) and (21). However, in order to efficiently compute the ground state energy of the Hamiltonian in Eq. (1) (see main text) for $N\gtrsim 25$ using DMRG, we approximate the Ising couplings using a translational invariant analytic function of the ion separation $r=|i-j|$. For $N<20$ the spin-spin coupling $J_{ij}$ between the two qubits at distance $r$ is well approximated by a power law decay:

[TABLE]

where, as stated in the main text, $J_{0}$ is the average nearest-neighbor coupling and $\alpha$ is the power law exponent Porras2004 . However for larger system sizes, this approximation fails to capture the actual decay of the interaction matrix.

In order to use the DMRG algorithm to accurately compute the ground state energies, we developed a compound function to better fit our couplings. This function is a product of a power law decay and an exponential decay parametrized by $J_{0}$, $\alpha^{\prime}$ and $\beta^{\prime}$:

[TABLE]

As seen in Fig. S5, this functional form fits well the exact Ising couplings even for a chain of 40 ions, while both a power law and a pure exponential fit diverge significantly.

I.3.4 State Detection

We detect the ion spin state by globally rotating all the spins into the measurement basis with a composite BB1 $\pi/2$ pulse as described above, to rotate the $x$ or $y$ basis into the $z$ basis), followed by the scattering of resonant laser radiation on the 2S${}_{1/2}|F=1\rangle\leftrightarrow^{2}$P${}_{1/2}|F=0\rangle$ cycling transition (wavelength near 369.5 nm and radiative linewidth $\gamma/2\pi\approx 20$ MHz). If the atom is projected in the $|\!\!\uparrow\rangle_{z}$ “bright” state, it fluoresces strongly, while if projected in the $|\!\downarrow\rangle_{z}$ “dark” state it fluoresces almost no photons because the laser is far from resonance Olmschenk2007 .

In both systems the fluorescence of the ion chain is imaged onto an Electron Multiplying Charge Coupled Device (EMCCD) camera (Model Andor iXon Ultra 897) using an imaging objective with 0.4 numerical aperture and a magnification of 90x for both systems. The fluorescence of each ion covers roughly a 7x7 array of pixels on the EMCCD. After collecting the fluorescence for an integration time of 0.65 (1) ms for system 1 (2), we use a binary threshold to determine the state of each ion, discriminating the quantum state of each ion with near 98$\%$ (97$\%$) accuracy in system 1 (2). The residual $2\,(3)\%$ errors include off-resonant optical pumping of the ion between spin states during detection as well as detector cross-talk between adjacent ions, readout noise, and background counts.

In system 2 the individual ion range-of-interests (ROIs) on the camera are updated with periodic diagnostic images, acquired by applying a nearly resonant cooling laser for 50 ms so that each ion fluoresces strongly regardless of its state. The signal to background noise ratio in the diagnostic shots is larger than 100, yielding precise knowledge of the ions’ center locations and taking into account the slow $\sim 2\,\mu$m pk-pk drift due to thermal expansion/contraction of the cryostat. Ion separations range from 1.5 $\mu$m to 3.5 $\mu$m depending on the trap settings and the distance from the chain center, and are always much larger than the resolution limit of the imaging system. We utilize the pre-determined ion centers to process the individual detection shots and optimize the integration area on the EMCCD camera to collect each ion’s fluorescence while minimizing cross-talk. We estimate cross-talk to be dominated by fluorescence from nearest-neighbor, which can cause a dark ion to be erroneously read as bright.

I.3.5 Error sources

The fidelity of the quantum simulation is limited by experimental noise that causes the system to depart from the ideal evolution and that can have several sources that are reviewed below. One important error source is off-resonant excitation of motional modes of the ion chain, which causes residual spin motion-entanglement. When the motion is traced out at the end of the measurement this results in a finite probability of an unwanted bit-flip. The probability of this error to occur on the $i$th ion Kim2009 is proportional to $p_{i}\sim\sum_{m=1}^{N}\left(\eta_{im}\Omega/\delta_{m}\right)^{2}$, where $\eta_{im}=b_{im}\sqrt{\nu_{R}/\nu_{\rm COM}}$ (see Eq. (20)) and $\delta_{m}=\mu-\omega_{m}$ is the beatnote detuning from the $m$-th normal mode. We trade off a lower error for a weaker spin-spin coupling by choosing a $\delta_{\rm COM}$ such that $\left(\eta_{\rm COM}\Omega/\delta_{\rm COM}\right)^{2}\lesssim 1/10$. By considering the off-resonant contributions of all the modes (see Fig. S6), we estimate the phonon error to cause about $1\%$ bit-flip per ion. Additionally, this effect is amplified by fluctuations in the trap frequency and laser light intensity at the ions’ location, increasing the probability of a phonon-assisted bit-flip event. To take this into account, we included slow drifts and fluctuations of the trap frequency and of the laser power on the timescale of 500 experimental repetitions assuming noise spectral density falling as $1/f$. Given our typical trap frequency and laser power fluctuations, we assume a relative standard deviation $\Delta\Omega/\Omega\sim 2\%$ and $\Delta\delta_{\rm COM}/\delta_{\rm COM}\sim 9\%$ over the timescale required to average over quantum projection noise and we end up estimating an average bit-flip probability $p_{i}\sim 9\%$ (see Fig. S6a). Moreover laser intensity, beam steering and trap frequency slow drifts over the time scale of a few hours required for data-taking cause averaging over different Ising parameters $J_{0}$. In particular, beam steering fluctuations create an imbalance between the red and blue $\nu_{0}\pm\mu$ beatnotes at the ions, producing an effective $B_{z}$ noisy field, that has been measured to be as high as $0.3J_{0}$. To take into account these drifts, we calculated several evolutions sampling from a gaussian distribution of values of $B_{z}$ and $J_{0}$, using as a variance the standard deviations ($\sigma_{J_{0}}=0.02J_{0}$ and $\sigma_{B_{z}}=0.3J_{0}$) observed in the experiment. Another source of bit-flip errors is imperfect detection. Off-resonant pumping limits our average detection fidelity to $98\%($97%$)$ for system 1 (2). A detection error is equivalent to a random bit-flip event so the two errors will sum up. A specific source of noise in system 2 is mechanical vibrations at 41 Hz and 39 Hz due to residual mechanical coupling to the cryostat Pagano2019 . This is equivalent to phase-noise on the Raman beams, which leads to dephasing of the qubits. Other less important noise sources are related to off-resonant Raman scattering errors during the Ising evolution (estimated in $7\cdot 10^{-5}$ per ion) and RF heating of the transverse COM motional mode of the ion chain in system 1.

In Fig. S6b, we plot the experimentally measured energy as a function of $\gamma$, and the corresponding theoretical curves with and without incorporating errors. Using the time dependent average bit-flip probability evolution that we estimated from our error model considering phonons and detection errors and averaging over slow drifts in experimental parameters $J_{0}$ and $B_{z}$, we get a good agreement with the experimental data (see also Fig. 2c in the main text, where the same parameters have been used), showing that we have a good understanding of the noise sources in our system.