Quantum Approximate Optimization
with a TrappedIon Quantum Simulator
Abstract
Quantum computers and simulators may offer significant advantages over their classical counterparts, providing insights into quantum manybody systems Feynman1982 () and possibly solving exponentially hard problems Nielsen2011 (), such as optimization Farhi2000 () and satisfiability Farhi2001 (). Here we report the first implementation of a shallowdepth Quantum Approximate Optimization Algorithm (QAOA) Farhi2014 () using an analog quantum simulator to estimate the ground state energy of the transverse field Ising model with tunable longrange interactions. First, we exhaustively search the variational control parameters to approximate the ground state energy with up to 40 trappedion qubits. We then interface the quantum simulator with a classical algorithm to more efficiently find the optimal set of parameters that minimizes the resulting energy of the system. We finally sample from the full probability distribution of the QAOA output with singleshot and efficient measurements of every qubit.
It remains an open question whether quantum computers or simulators can efficiently solve classically hard combinatorial problems. One possible approach is to use quantum devices to perform quantum annealing Kadowaki1998 (); Farhi2000 (), which in some cases may provide an advantage over classical approaches Santoro2002 (); Hastings2013 (). However, quantum annealing has stringent adiabaticity requirements that hinder its applicability on existing quantum platforms that have finite coherence times Richerme2013 ().
Alternatively, hybrid quantumclassical variational algorithms may approximately solve hard problems in realms such as quantum magnetism, quantum chemistry mcclean2016 (), and highenergy physics Kokail2018 (). This is because the key resource of quantum computers and simulators is quantum entanglement, which is exactly what makes these manybody 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 “bangbang” protocol that can provide an approximate answer in a timeefficient way, using devices with finite coherence times and without the use of errorcorrection Preskill2018 (); Zhou2018 ().
Similarly to quantum annealing, 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. The QAOA is theoretically motivated from complexity theory harrow2016qaoa (); Lloyd2018 () and optimal control theory mcclean2016 (); Yang2017 (); bapat2018 (), and may be seen more intuitively as replacing infinitesimal (Trotterized) Hamiltonian time evolution by an appropriate sequence of finite time steps. The QAOA relies on a classical outer loop to optimize the quantum circuit, aided by physical intuition bapat2018 (); verdon2018 () or observed structure of the variational parameters Zhou2018 (); brady2019 (); crooks2018 (), producing fast, lowdepth circuits for approximate solutions.
In this work, we employ a collection of interacting trappedion qubits to experimentally implement a specific instance of the QAOA. We focus on the energy minimization of a transverse field antiferromagnetic Ising Hamiltonian with longrange interactions Zhang2017a ():
(1) 
Here we set Planck’s constant , () is the Pauli matrix acting on the spin along the direction of the Bloch sphere, is the Ising coupling between spins and , and denotes the transverse magnetic field. It is wellknown Koffel2012 () that the Hamiltonian (1) exhibits a quantum phase transition for antiferromagnetic interactions with power law decay. One of the goals of this work is to find an approximation of the ground state energy at the critical point , where is the average nearestneighbour coupling. In this case, the realization of the QAOA entails a series of unitary quantum evolutions (see Fig. 1) under the noncommuting Hamiltonians and that are applied to a known initial state . The state obtained after layers of the QAOA is:
(2) 
where the angles and are variational parameters used in the th QAOA layer to minimize the final energy .
In order to implement the quantum optimization algorithm, each spin in the chain is encoded in the S and hyperfine “clock” states of a Yb ion (see Methods). 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 Methods). Both systems are based on a linear rf Paul trap where we store chains of up to ions and initialize the qubits in the ground state of , namely the product state , where and is assumed to be negative. The unitary evolution under is realized by generating spinspin interactions through spindependent optical dipole forces implemented by an applied laser field. This gives rise to effective tunable longrange Ising couplings that fall off approximately as Porras2004 (). The powerlaw exponent and the interaction strengths vary in the range (0.30.57) kHz, depending on the system size and the experimental realization (see Methods for details). The unitary evolution under is generated by applying a global rotation around the axis of the Bloch sphere.
After each run of the algorithm, we perform a projective measurement of each spin in the basis to measure () (see Fig. 1). Measurements in the and bases are carried out by performing a rotation about the ()axis of the Bloch sphere, illuminating the ions with resonant laser light, and collecting the dependent fluorescence on a camera with siteresolved imaging. The energy is calculated by combining the measurements of the twobody correlators and the total magnetization along the axis , where the indices range from 1 to . We benchmark the experimental outcome with the ground state 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
(3) 
where is the energy of the highest excited state. This choice maps the entire manybody spectrum to the interval. In the following we show that the best experimental performance is close to the theoretical performance , which itself is less than unity for a finite number of QAOA layers.
We first focus on the case, where only two variational parameters ( and ) are used to minimize the energy. In this case, the timeevolved one and twopoint correlation functions can be efficiently computed dylewsky2016 (). This leads to a general formula for the energy expectation under a state produced by the QAOA that is used to compute the theoretical performance of the algorithm (see Methods). In Fig. 2a we show an experimental exhaustive search over the parameter space and compare it to the theoretical performance of the algorithm, showing good agreement for qubits. We also compare the performance of our algorithm as a function of with the expected QAOA performance (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 antiferromagnetic states. We locate this critical point at for 20 qubits by computing the halfchain entanglement entropy of the ground state numerically. As shown in Fig. 2b, while the experimental performance is when is above the critical point, the gain relative to the initial state 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.
We investigate the performance of the QAOA algorithm as a function of the number of qubits. For each system size, we ensure that the spinspin couplings have the same dependence on the qubit distance by varying the trap parameters (see Methods). As shown in Fig. 2c, the halfchain entanglement entropy as a function of system size exhibits a peak located at , displaying the onset of the phase transition as tends to infinity. For all system sizes, we optimize the algorithm by performing a scan of the interaction angle and applying discrete variations of the mixing angle around the optimal value predicted by the theory. In Fig. 2c we compare the optimal experimental and theoretical performances for different system sizes from 20 up to 40 qubits for fixed . We conclude from both the numerics and experimental data that the QAOA yields a similar performance as a function of number of qubits, even if the algorithm runtime changes negligibly as the number of qubits increases.
We experimentally perform a search for the optimal QAOA performance using 20 qubits. Unlike the case, there is no known analytic formula to efficiently compute the energy. However, exploiting relationships between optimal angles as a function of increasing , we use a bootstrapping heuristic (see Methods for details) that allows the experiment to identify a set of optimal angles faster than a global parameter search. We start from the optimal guess and perform a fine scan of , while varying and in larger steps. The result is shown in Fig. 2d, where we plot the performances as a function of 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 and . The QAOA performance with 20 qubits is in agreement with the performance in system 2, taken with the same parameters (see Fig. 2c). This indicates that decoherence and bitflip errors (see Methods) accumulated during longer evolution times are already balancing out the expected performance gain of one additional optimization layer.
As a brute force approach is inefficient, we implement a closedloop QAOA by interfacing the analog trappedion quantum simulator with a greedy gradientdescent algorithm to optimize the measured energy. In the QAOA, we can visualize the optimization trajectory on the theoretical performance surface as shown in Fig. 3. The procedure goes as follows: (i) we start from a guess and measure its energy. (ii) We measure the approximate local gradient by performing the energy measurements in two orthogonal directions and . (iii) The algorithm computes the new guess . The new energy is then measured on the quantum simulator. Steps  are iterated until convergence of the algorithm. As shown in Fig. 3, the algorithm converges after about 10 iterations. Compared to an exhaustive search, the gradient descent uses a fewer number of 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.
We further explore the performance of our system by investigating higherorder observables in the output of the QAOA, using highfidelity singleshot measurement of all the qubits. It has been proven, under reasonable complexitytheoretic assumptions, that no classical algorithm can efficiently sample exactly from a sufficiently general class of QAOA circuits harrow2016qaoa (). Recent results Bremner2016 (); Bouland2019 () suggest that this could also hold in the case of approximate sampling (see Methods). In order to approximately sample from the QAOA, we optimize the algorithm to solve for the ground state of alone (see Eq. (1)), namely in the absence of a transverse field. In this case, by measuring in the basis, it is possible to sample the probability distribution of all the eigenstates of the Hamiltonian . We performed the experiment with 12 qubits so that we can both compute the expected QAOA theoretical output and also experimentally oversample the Hilbert space of all the possible possible outcomes. In Fig. 4a we show the raw QAOA outcome using the optimal variational parameters and compare the experimental eigenstate histogram with the exact diagonalization prediction of the QAOA output state, showing qualitative agreement.
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 phononassisted bitflip events (see Methods). We can take into account these errors by implementing the coarsegraining 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 around the most probable state, producing a coarsegrained dataset. We then apply the same coarsegraining 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 .
In order to quantitatively compare the coarsegrained experiment and the theory, we estimate two different metrics, namely the total variation distance (TVD) and the KullbackLeibler divergence (), defined as:
(4)  
(5) 
where is the experimental (theoretical) probability of observing the th outcome. As shown in Fig. 4c, when the system is in the initial state, it is closer to a uniform probability distribution since is an equal superposition of all the eigenstates of . On the other hand, as the parameter is scanned, we observe a net decrease of both TVD and , in agreement with the decrease in energy, computed by measuring one and twobody correlators.
The variational quantum algorithm reported here, with up to 40 trappedion qubits, is the largest ever realized on a quantum device. We approximate the ground state energy of a nontrivial quantum Hamiltonian showing almost constant time scaling with the system size. Singleshot highefficiency qubit measurements in different bases give access to the full distribution of bitstrings that is difficult or potentially impossible to model classically. With the addition of individual control over the interactions between qubits, the variational quantumclassical hybrid approach can be employed in this experimental platform to give insight into quantum chemistry Kandala2017 (); Hempel2018 (); Nam2019 () and hard optimization problems, such as MaxSAT 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, the NSF PFCQC Program, the DOE BES QIS Program, the DOE HEP QuantISED Program, the DOE ASCR Quantum Testbed Pathfinder and Quantum Algorithms Teams (QOALAs) programs, and the NSF Physics Frontier Center at JQI.
.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)
References
 (1) Feynman, R. P. Simulating physics with computers. International Journal of Theoretical Physics 21, 467–488 (1982).
 (2) Nielsen, M. A. & Chuang, I. L. Quantum Computation and Quantum Information: 10th Anniversary Edition (Cambridge University Press, New York, NY, USA, 2011), 10th edn.
 (3) Farhi, E., Goldstone, J., Gutmann, S. & Sipser, M. Quantum Computation by Adiabatic Evolution. arXiv eprints quant–ph/0001106 (2000). eprint quantph/0001106.
 (4) Farhi, E. et al. A quantum adiabatic evolution algorithm applied to random instances of an npcomplete problem. Science 292, 472–475 (2001).
 (5) Farhi, E., Goldstone, J. & Gutmann, S. A Quantum Approximate Optimization Algorithm. arXiv eprints arXiv:1411.4028 (2014). eprint 1411.4028.
 (6) Kadowaki, T. & Nishimori, H. Quantum annealing in the transverse Ising model. Phys. Rev. E 58, 5355–5363 (1998). eprint condmat/9804280.
 (7) Santoro, G. E., Martoňák, R., Tosatti, E. & Car, R. Theory of quantum annealing of an ising spin glass. Science 295, 2427–2430 (2002).
 (8) Hastings, M. B. & Freedman, M. H. Obstructions To Classically Simulating The Quantum Adiabatic Algorithm. arXiv eprints arXiv:1302.5733 (2013). eprint 1302.5733.
 (9) Richerme, P. et al. Experimental performance of a quantum simulator: Optimizing adiabatic evolution and identifying manybody ground states. Phys. Rev. A 88, 012334 (2013).
 (10) McClean, J. R., Romero, J., Babbush, R. & AspuruGuzik, A. The theory of variational hybrid quantumclassical algorithms. New Journal of Physics 18, 023023 (2016).
 (11) Kokail, C. et al. SelfVerifying Variational Quantum Simulation of the Lattice Schwinger Model. arXiv eprints arXiv:1810.03421 (2018). eprint 1810.03421.
 (12) Preskill, J. Quantum Computing in the NISQ era and beyond. Quantum 2, 79 (2018).
 (13) Zhou, L., Wang, S.T., Choi, S., Pichler, H. & Lukin, M. D. Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on NearTerm Devices. arXiv eprints arXiv:1812.01041 (2018). eprint 1812.01041.
 (14) Farhi, E. & Harrow, A. W. Quantum supremacy through the quantum approximate optimization algorithm. arXiv preprint arXiv:1602.07674 (2016).
 (15) Lloyd, S. Quantum approximate optimization is computationally universal. arXiv eprints arXiv:1812.11075 (2018). eprint 1812.11075.
 (16) Yang, Z. C., Rahmani, A., Shabani, A., Neven, H. & Chamon, C. Optimizing variational quantum algorithms using pontryagin’s minimum principle. Physical Review X 7, 1–8 (2017). eprint 1607.06473.
 (17) Bapat, A. & Jordan, S. Bangbang control as a design principle for classical and quantum optimization algorithms. arXiv preprint arXiv:1812.02746 (2018).
 (18) Verdon, G., Pye, J. & Broughton, M. A universal training algorithm for quantum deep learning. arXiv preprint arXiv:1806.09729 (2018).
 (19) Brady, L., Bapat, A. & Gorshkov, A. in preparation (2019).
 (20) Crooks, G. E. Performance of the quantum approximate optimization algorithm on the maximum cut problem. arXiv preprint arXiv:1811.08419 (2018).
 (21) Zhang, J. et al. Observation of a manybody dynamical phase transition with a 53qubit quantum simulator. Nature 551, 601–604 (2017).
 (22) Koffel, T., Lewenstein, M. & Tagliacozzo, L. Entanglement entropy for the longrange ising chain in a transverse field. Physical Review Letters 109 (2012).
 (23) Kim, K. et al. Entanglement and tunable spinspin couplings between trapped ions using multiple transverse modes. Physical Review Letters 103 (2009).
 (24) Pagano, G. et al. Cryogenic trappedion system for large scale quantum simulation. Quantum Science and Technology 4, 014004 (2019).
 (25) Porras, D. & Cirac, J. I. Effective Quantum Spin Systems with Trapped Ions. Phys. Rev. Lett. 92, 207901 (2004).
 (26) Jaschke, D., Wall, M. L. & Carr, L. D. Open source matrix product states: Opening ways to simulate entangled manybody quantum systems in one dimension. Computer Physics Communications 225, 59 – 91 (2018).
 (27) Dylewsky, D., Freericks, J., Wall, M., Rey, A. & FossFeig, M. Nonperturbative calculation of phonon effects on spin squeezing. Physical Review A 93, 013415 (2016).
 (28) Bremner, M. J., Montanaro, A. & Shepherd, D. J. Averagecase complexity versus approximate simulation of commuting quantum computations. Phys. Rev. Lett. 117, 080501 (2016).
 (29) Bouland, A., Fefferman, B., Nirkhe, C. & Vazirani, U. On the complexity and verification of quantum random circuit sampling. Nat. Phys. 15, 159–163 (2019).
 (30) Wang, S.T. & Duan, L.M. Certification of Boson Sampling Devices with CoarseGrained Measurements. arXiv eprints arXiv:1601.02627 (2016). eprint 1601.02627.
 (31) Kandala, A. et al. Hardwareefficient variational quantum eigensolver for small molecules and quantum magnets. Nature 549, 242 EP – (2017).
 (32) Hempel, C. et al. Quantum chemistry calculations on a trappedion quantum simulator. Phys. Rev. X 8, 031022 (2018).
 (33) Nam, Y. et al. Groundstate energy estimation of the water molecule on a trapped ion quantum computer. arXiv eprints arXiv:1902.10171 (2019). eprint 1902.10171.
 (34) Giovannetti, V., Lloyd, S. & Maccone, L. Advances in quantum metrology. Nature Photonics 5, 222 EP – (2011).
 (35) Farhi, E., Goldstone, J. & Gutmann, S. A Quantum Approximate Optimization Algorithm Applied to a Bounded Occurrence Constraint Problem. arXiv eprints arXiv:1412.6062 (2014). eprint 1412.6062.
 (36) Wineland, D. et al. Experimental issues in coherent quantumstate manipulation of trapped atomic ions. J. Res. Natl. Inst. Stand. Technol. 103, 259–328 (1998).
 (37) Olmschenk, S. et al. Manipulation and detection of a trapped hyperfine qubit. Phys. Rev. A 76, 052314 (2007).
 (38) Brown, K. R., Harrow, A. W. & Chuang, I. L. Arbitrarily accurate composite pulse sequences. Phys. Rev. A 70, 052318 (2004).
 (39) Sørensen, A. & Mølmer, K. Quantum computation with ions in thermal motion. Phys. Rev. Lett. 82, 1971–1974 (1999).
 (40) James, D. F. V. Quantum dynamics of cold trapped ions with application to quantum computation. Applied Physics B 66, 181 (1998).
Appendix A Methods
a.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 operators, typically a problem (or cost) Hamiltonian that encodes the cost function on the diagonal in (say) the basis, and a transverse term that generates transitions between bit strings, such that the initial state evolves into an approximate ground state of .
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 (). We show in this work that 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 state preparation in a scalable fashion and with a low overhead on quantum resources.
First, we discuss how to discover optimal QAOA1 schedules, i.e., QAOA schedules for .
a.1.1 Qaoa,
Despite its apparent simplicity, the QAOA (or QAOA1) can be a powerful state preparation ansatz. For example, hardnessofsampling results are known for QAOA1 circuits harrow2016qaoa (), closely mirroring the hardness of sampling from instantaneous quantum polynomial (IQP) circuits. 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 (). For a twolocal transverse field spin Hamiltonian as in Eq. 1, this leads to a formula for the energy expectation under a state produced by the QAOA1, starting from the product state (see Supplement for the result and derivation). These formulas are applicable to many cases of interest in quantum state preparation and optimization. Importantly, the time complexity to compute the formula is in the worst case, making it tractable to optimize the QAOA1 protocols for large spin chains.
a.1.2 Qaoa,
The general analytical formula for 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 for each layer . For layers, this is an optimization on a 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 , form smooth curves for any . While this observation points to a deeper theoretical mechanism at play, it doesn’t 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 , denote the special optimal angles by , which we can also think of as a pair of angle curves (as a function of step index ). As is varied, we may think of these minima as a family. We numerically find that this family exhibits the following desirable features (for sufficiently large):

The angles are nonnegative, small and bounded.

For sufficiently large, the two angle sequences and are approximately smooth.

The angle sequence (and correspondingly, ) when viewed as a function on the normalized time parameter , is convergent in the parameter . In other words, as is increased, the angle sequences and approach a smooth, asymptotic curve (See Figure S1.)

The energy expectation approaches the global minimum as , 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 may be discovered using the optimal schedule for 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 . The main idea behind this construction is that the minimal angle curves for a larger may be guessed from the optimal curve of a smaller by interpolation.
Utilizing the above points, we use a bootstrapping algorithm to find the optimal angle sequences, and , for a given , as described below. Let denote an intermediate angle index. Then:

For , use an analytic formula to find and .

For , choose an initial guess of and .

Perform a local optimization of and in order to find and .

Repeat the next steps (57) for .

Create interpolating functions through the angle sequences, and , using the normalized time as the independent parameter (we use a linear interpolation for and cubic for ).

Choose the initial guesses for and by sampling the interpolating function from (5) at evenly spaced points separated by a normalized time distance of .

Perform a local optimization of and in order to find and .
The resulting angles and should be at least a good local minimum of the energy expectation value and approaches the global minimum as .
The interpolation in step 2 is based on our observation that the angles tend to curve down at the end and the angles tend to curve up.
An important feature of our algorithm is that its asymptotic runtime is expected to be efficient in . This feature is predicated on the previous result that the angle curves are generally convergent as tends to infinity. The argument proceeds as follows: if we assume a maximal deviation of the initial guess for layer to be , then the total distance between the initial guess and the optimized curve is no greater than . Therefore, the local search algorithm is confined to a ball of radius at most , and for a fixed error tolerance, the convergence time for a standard local optimizer is . Summing over convergence times for all from , we have
(6) 
The last inequality above comes about as follows: while the summand depends on the convergence rate of the sequence , it is upper bounded by for a converging set of paths and an initial error of order 1. The latter is true since our angle search domain is bounded and independent of . Therefore, the sum is no greater than . In practice, even faster runtimes are possible. Therefore, the bootstrap algorithm exploits the structure of the special minima and provides a scalable route to multistep QAOA for the longrange 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.
a.1.3 Convergence in
In the previous sections, we introduced a bootstrap algorithm that is asymptotically efficient in the number of layers . However, in order to be fully scalable the algorithm must also be scalable in the system size . 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 longrange TFIM, and indeed any translationallyinvariant model with a welldefined notion of metric and dimension arising from the functional form of the coupling coefficients , it is reasonable to expect that the optimized angles depend on system size in a predictable way. This is indeed the case for the longrange TFIM. There, it can be seen that the angle curves for varying appear similar in shape. Usefully, the curves also appear to be convergent to an idealized curve for a hypothetical continuous, longrange 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 . For a given coupling function (such as inverse powerlaw) that decays as a function of distance, one may define a skin depth that is the number of sites from the boundaries that the coupling is, say, a factor of smaller than the nearestneighbour value. Clearly, is independent of the system size and depends only on the parameters of the coupling function. For instance, for the longrange TFIM, . As tends to infinity, the fractional skin depth 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 to 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 to is . Therefore, if the change in the optimal QAOA curves (in, say, distance) is a smooth function of the the fractional skin depth, then we expect it to vary as . Therefore, the total running time of a bootstrap from small system sizes to a given size should be which is sublinear in . Combining this observation with the convergence in , we see that for a given Hamiltonian family, optimized QAOA angle curves for small may be used as a rubric for the optimization for longer circuit depths. Furthermore, if the Hamiltonian is translationallyinvariant with decaying interactions, the optimized QAOA schedules are expected to scale with 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 (Figure S2). This is our main theoretical contribution in this work.
a.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 so that the experiment is equivalent to preparing a state , evolving under a Hamiltonian diagonal in the computational basis, followed by a uniform rotation and measurement in the computational basis. Following Ref. Farhi2014 (), it suffices to consider QAOA circuits with . The output state is for some cost function diagonal in the computational basis.
a.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 basis state. In Ref. Bremner2016 (), the output state after a socalled IQP circuit (which only differs from the one here in that the final rotation is a global Hadamard instead of ) has an amplitude proportional to a quantity known as the gap of a Boolean function, , the difference in the number of inputs that map to 1 and the number of inputs that map to 0 under . Finding the gap of a general function is a complete problem. This is a very hard problem since the class includes , which in turn includes the whole of . The authors of Ref. Bremner2016 () prove that the gap of a degree3 polynomial over , , may be expressed as an output amplitude of an IQP circuit. They also show that the finding the gap of such functions is still complete. Following Ref. Bremner2016 (), we examine the output amplitude of a QAOA state:
(7) 
where now we define the function to have the range and the Hamiltonian satisfies for a computational basis state . The output amplitude is thus proportional to a ‘generalized gap’ of a function , where is the Hamming weight of . This modified function is also a degree3 polynomial over . Note that this restriction to degree3 comes from the fact that the gates , and are universal for classical computation (indeed, the Toffoli alone is universal for classical computation) and there is a natural degree3 polynomial coming from this construction. The quantity we have defined, , can be easily shown to be hard to compute, by reducing to . This suffices for exact sampling hardness assuming the polynomial hierarchy () does not collapse.
a.2.2 Approximate sampling hardness
For approximate sampling hardness, we need two other properties, namely anticoncentration and a worsttoaverage case reduction. Anticoncentration 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 such that the function is a degree3 polynomial with uniformly random weights and . Anticoncentration then follows from the PaleyZygmund inequality and Lemma 4 of the Supplemental Material of Ref. Bremner2016 () (with ).
Finally, we need to show that the problem of approximating the generalized gap is averagecase 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 worsttoaverage 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 available here, which continuously changes the parameters and . However, we have only shown anticoncentration when the weights and are chosen from a finite set. It remains to be seen whether one can have the property of anticoncentration and averagecase hardness holding at the same time for some specific QAOA output distribution.
a.3 Trappedion experimental systems
In this work two quantum simulators have been used, referred to as system 1 and 2. System 1 Kim2009 () is a roomtemperature iontrap apparatus, consisting of a 3layer linear Paul trap with transverse centerofmass motional frequency MHz and axial centerofmass frequencies ranging from to 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 rfheating or other mechanisms.
System 2 Pagano2019 () is a cryogenic iontrap 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 cryopumping reduces the background pressure at low 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 MHz and MHz, and the axial frequency ranges from to MHz.
a.4 State preparation
The qubit is initialized by applying resonant 369.5 nm light for about 20 s to optically pump into the state. To perform global rotations in the Bloch sphere, we apply two fardetuned, noncopropagating Raman beams whose beatnote is tuned to the hyperfine splitting GHz of the clock states and encoding the qubit Olmschenk2007 (). State preparation in our implementation of the QAOA requires qubit initialization in the state by optically pumping the ions and then a global rotation into the state using stimulated Raman transitions. We detect the state of each ion at the end of each experimental sequence using statedependent 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 state preparation fidelity. The BB1 four pulse sequence is:
where after the first rotation , three additional rotations are applied: a pulse along an angle , a pulse along , and another pulse along . The axes of these additional rotations are in the  plane of the Bloch sphere with the specified angle referenced to the axis.
Generating the Ising Hamiltonian
We generate spinspin interactions by employing a spin dependent force with a pair of noncopropagating 355 nm Raman beams, with a wavevector difference aligned along the transverse motional modes of the ion chain. The two offresonant Raman beams are controlled using acoustooptic modulators which generate two interference beatnotes at frequencies in the MølmerSørensen configuration Molmer1999 (). In the LambDicke regime, the laserion interaction gives rise to the effective spinspin Hamiltonian in Eq. 1, where the coupling between the th and th ion is:
(8) 
Here is the Rabi frequency, is the recoil frequency, is the frequency of the th normal mode, is the eigenvector matrix element for the th ion’s participation to the th normal mode James1998 (), and is the mass of a single ion.
Differently from system 1, where the wavevector difference of the Raman beams is aligned along a principal axis of the trap, in system 2 the spinspin interaction stems from the off resonant coupling to both families of transverse normal modes. Eq. (8) is then generalized to:
(9) 
where is the recoil frequency given by the projection of the Raman wavevector along the two radial principal axes of the trap . We infer an angle between and the principal axis (see inset in Fig. S3a) from the ratio between the spinphonon couplings to the two radial center of mass modes. Before every experiment, we perform Raman sideband cooling on both the centerofmass 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 selfsimilar functional form of the spinspin interaction (see Fig. S3b). For the data in Fig. 2c, we set the detuning to kHz and the axial frequency to MHz, for respectively.
a.4.1 Fitting Ising Couplings to Analytic Form
By directly measuring trap parameters and spinphonon couplings, we can calculate the spinspin interaction matrix with Eqs. (8) and (9). However, in order to efficiently compute the ground state energy of Hamiltonian (1) for using DMRG, we approximate the Ising couplings using a translational invariant analytic function of the ion separation . For the spinspin coupling between the two qubits at distance is well approximated by a power law decay:
(10) 
where, as stated in the main text, is the average nearestneighbor coupling and 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 , and :
(11) 
As seen in Fig. S4, this functional form fits the exact Ising couplings well even for a chain of 40 ions, while both a power law and a pure exponential fit diverge significantly.
a.5 State Detection
We detect the ion spin state by globally rotating all the spins into the measurement basis with a composite BB1 pulse as described above, to rotate the or basis into the basis), followed by the scattering of resonant laser radiation on the SP cycling transition (wavelength near 369.5 nm and radiative linewidth MHz). If the atom is projected in the “bright” state, it fluoresces strongly, while if projected in the “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 errors include offresonant optical pumping of the ion between spin states during detection as well as detector crosstalk between adjacent ions, readout noise, and background counts.
In system 2 the individual ion rangeofinterests (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 m pkpk drift due to thermal expansion/contraction of the cryostat. Ion separations range from 1.5 m to 3.5 m depending on the trap settings and the distance from the chain center, and are always much larger than the resolution limit of the diffractionlimited imaging system. We utilize the predetermined 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 crosstalk. We estimate crosstalk to be dominated by fluorescence from nearestneighbor, which can cause a dark ion to be erroneously read as bright.
a.6 Error sources
The fidelity of the quantum simulation is limited by experimental noise. A significant type of error that causes the system to depart from the ideal evolution is “bitflip errors”. These errors are caused by offresonant excitation of motional modes of the ion chain, which causes residual spin motionentanglement. When the motion is traced out at the end of the measurement this results in an unwanted bitflip. The probability of this error to occur on the th ion Kim2009 () is proportional to , where (see Eq. 8) and is the beatnote detuning from the th normal mode. We trade off a lower error for a weaker spinspin coupling by choosing a such that . By considering all the offresonant contributions of all the modes (see Fig. S5), we estimate the phonon error to cause about bitflip 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 phononassisted bitflip 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 440 experimental repetitions assuming noise spectral density falling as . Assuming a reasonable relative standard deviation and over the timescale required to average over quantum projection noise, we end up estimating an average bitflip probability (see Fig. S5a). Another source of bitflip errors is imperfect detection. Offresonant pumping limits our average detection fidelity to 97% for system 1 (2). A detection error is equivalent to a random bitflip event so the two errors will sum up.
Slow experimental drifts over the time scale of a few hours required for datataking also cause averaging over different Ising Hamiltonian parameters. These fluctuations are related to laser light intensity noise at the ions’ location and drifts of the trap frequency. Other less important noise sources are related to offresonant Raman scattering errors during the Ising evolution (estimated in per ion) and RF heating of the transverse COM motional mode of the ion chain in system 1. 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. This is equivalent to phasenoise on the Raman beams, which leads to dephasing of the qubits.
In Fig. S5b, we plot the experimentally measured energy as a function of , and the corresponding theoretical curves with and without incorporating errors. Using the bitflip probability estimated from our model considering phonons and detection errors, we get a better agreement with the experimental data, showing that we have a good understanding of the noise sources in our system.