Quantifying the efficiency of state preparation via quantum variational eigensolvers

Recently, there has been much interest in the efficient preparation of complex quantum states using low-depth quantum circuits, such as Quantum Approximate Optimization Algorithm (QAOA). While it has been numerically shown that such algorithms prepare certain correlated states of quantum spins with surprising accuracy, a systematic way of quantifying the efficiency of QAOA in general classes of models has been lacking. Here, we propose that the success of QAOA in preparing ordered states is related to the interaction distance of the target state, which measures how close that state is to the manifold of all Gaussian states in an arbitrary basis of single-particle modes. We numerically verify this for the ground state of the quantum Ising model with arbitrary transverse and longitudinal field strengths, a canonical example of a non-integrable model. Our results suggest that the structure of the entanglement spectrum, as witnessed by the interaction distance, correlates with the success of QAOA state preparation. We conclude that QAOA typically finds a solution that perturbs around the closest free-fermion state.


I. Introduction
In recent years, algorithms involving a hybrid quantum-classical procedure for cost function minimization have attracted much attention. 1,2 Among these is the Quantum Approximate Optimization Algorithm (QAOA), which employs an alternating operator ansatz for solving optimization problems that are mappable to the problem of finding the ground state of a classical Ising-type Hamiltonian. 2 These ansatz circuits are of great interest since they have been shown to successfully approximate or even exactly prepare states at remarkably low circuit depths. This makes them amenable to implementation using the currently available "noisy intermediate scale" quantum computers, 3 potentially enabling useful applications in the near future before full-fledged quantum computers, featuring robust error correction, may become operational. It has also been proven that QAOA circuits can perform universal quantum computation for certain classes of Hamiltonians. 4,5 QAOA was originally proposed to tackle classical optimization problems, such as the MaxCut problem, 6 and several others. 7-10 More recently, it has been pointed out that QAOA could also serve as a tool for exactly preparing quantum many-body states, such as the GHZ state, the ground state of the Ising model at the critical point, for both short-range 11 and long-range interactions, 12,13 the ground state of the toric code, 11 the ground state of the two-dimensional Hubbard model, 14 and the thermofield double states. 15,16 In this paper we focus on the latter type of applications of QAOA in the context of preparing ground states of non-integrable quantum Hamiltonians.
Some analytical results on state preparation in the classical Ising model have been established. It was shown that in one dimension, the uniform nearest-neighbor Ising model in the absence of a magnetic field can be reduced to a system of pseudospins. 17 This reduction was used the prove a conjecture 2 that the ground state of this model can be prepared using a circuit with depth linear in the system size, 17 and the associated bounds on the best attainable fidelity for QAOA circuits below this depth were derived. 18 At the same time, a deeper understanding of QAOA that would, for example, allow for the systematic construction of a circuit to prepare a given quantum state, and to predict the circuit depth needed to reach a certain accuracy, is still lacking.
In its original formulation, QAOA involved the alternating application of a "mixer" Hamiltonian and a "problem" Hamiltonian, 2,8 where different choices for the mixer Hamiltonian can be made. 2,19 A different approach is to split the problem Hamiltonian into a number of terms and alternate between the application of those; in this way, QAOA can be seen as the digitization plus Trotter splitting of a quantum annealing protocol. 18 This justifies that, in this latter case, as the circuit depth increases, states can be prepared with increasing accuracy. It does not explain, however, why some states can be prepared with very high accuracy using low circuit depths. Thus, there remain open questions about the inner workings of QAOA. Some of the difficulties in developing a deeper understanding of QAOA, as well as its numerical implementations, stem from the fact that the optimization landscape is generally riddled with local minima and other issues. [20][21][22] Different techniques, such as different optimization methods [23][24][25] or machine learning, [26][27][28][29] have attempted to address these problems. Heuristics for producing a starting set of optimization angles have also been developed. 30,31 In this paper, we address the problem of using QAOA to prepare the ground state of some quantum Hamiltonian H that depends on one or more tunable parameters H(h 1 , h 2 , . . .), such that H is non-integrable for generic values of h 1 , h 2 , etc. We are interested in predicting the relative success of QAOA state preparation across the phase diagram defined by the parameters h 1 , h 2 , . . .; moreover, our aim is to relate the success of preparation to some physical property of the target state. While standard quantities such as the von Neumann (entanglement) entropy of the target state poorly correlate with the success of QAOA, we find that the success of state preparation correlates with the property called interaction distance. 32 The latter can be evaluated from the eigenvalue spectrum of the system's (reduced) density matrix. Our findings are numerically supported by examples of non-integrable quantum Ising models.
The remainder of this paper is organized as follows. Secs. II and III contain a brief overview of the QAOA variational ansatz and interaction distance, respectively. In Sec. IV, we introduce an alternating operator protocol for the Ising model in transverse and longitudinal fields, and we compare the success of the QAOA ground-state preparation with interaction distance across the phase diagram of the model. In Sec. V we provide analytical arguments for the numerically-observed correlation between QAOA and interaction distance, while in Sec. VI we analyze the optimization landscape for these models. Our conclusions are presented in Sec. VII, while Appendices contain generalizations of our results; in particular Appendix A contains the results for a variant of the Ising model which features interactions between nearestneighbour triplets of spins. This model has a critical line in the universality class of the Potts model, 33 and its ground state is much harder to prepare than that of the quantum Ising model.

II. Quantum Approximate Optimization Algorithm
Various names for quantum-classical variational algorithms have been proposed in the literature, depending on the context and the specific implementation. 1, 2,11,14,34,35 Among the first such algorithms is the Variational Quantum Eigensolver (VQE) -proposed in the context of quantum chemistry 1,34 -for preparing approximate eigenstates and calculating eigenvalues of a given Hamiltonian. The QAOA 2 introduced the alternating operator ansatz, which we review below. We refer to the general class of variational quantum-classical algorithms based on the alternating operator ansatz simply as QAOA, with the understanding that it can be seen as a specialization of VQE for this particular class of ansätze.
As previously mentioned, QAOA is a variational algorithm based on a classical optimization routine which performs a minimization over a parametrized family of quantum circuits. The goal of this minimization is to find the circuit which, starting from some initial state |ψ initial , best prepares the target quantum state, |ψ target . This family of quantum circuits is defined by a set of operators H 1 , H 2 , . . . , H M , and takes the alternating operator "bang-bang" form, 36 defined by the unitary The circuit is parameterized by the set of variational angles θ = (θ 1,1 , ..., θ p,M ).
A sketch of the QAOA protocol is given in Fig. 1(a). The algorithm starts with some chosen initial state |ψ initial and an initial set of values for the circuit parameters. The initial state |ψ initial is, in principle, arbitrary, but it should be sufficiently easy to prepare (e.g., a product state or some low-entangled state). The target state ψ target is often assumed to be the ground state of some Hamiltonian H, sometimes called the "problem Hamiltonian". As mentioned previously, there is freedom in the choice of the operators {H j } j∈{1,...,M } . Unlike the original formulation, 2 for the problems considered in this paper, we choose the operators by splitting the problem After preparing the initial state, the simulator performs the quantum evolution After the state |ψ(θ) is obtained, a cost function is measured. In what follows, we assume all states to be normalized, and the cost function may be defined as the expectation value of the energy It is often more convenient to use the rescaled relative energy 37 where E min , E max are the extremal eigenenergies in the spectrum of H. The relative energy is bounded between 0 and 1, such that = 0 corresponds to finding an exact ground state. Alternatively, if |ψ target is known, the cost function can be taken to be the quantum infidelity, which is similarly bounded between 0 and 1. Although evaluating fidelity in experiment is impractical or even impossible, it is often useful in numerical simulations. Once the value of the cost function is measured, it is passed back to the optimization algorithm running on the classical computer. This algorithm returns a new set of angles, which are passed again to the quantum simulator, and the process repeats itself until the optimization algorithm running on the classical computer halts.
In Ref. 11 it was observed that the ground state of the transverse-field Ising model with periodic boundary conditions could be prepared exactly (i.e., with f = 1) in precisely p = N/2 steps, where N is the total number of spins. This was done using the same M = 2operator QAOA protocol originally proposed for the MaxCut problem. 2 This was a surprising result, given that the ground state of the Ising model can be very complex depending on the magnitude of the transverse magnetic field. For instance, at the critical value of the field, the excitation gap goes to zero and the ground state displays logarithmically diverging von Neumann entropy (VNE) of entanglement as a function of subsystem size. 38 This example demonstrates that the success of QAOA protocol is not determined by the VNE of the target state. One of the main results of the present paper is to show that a different quantum-information measure called interaction distance 32 serves as an error-metric for the quality of QAOA state preparation. In the following section, we briefly introduce and review the properties of interaction distance (see also Ref. 39).

III. Interaction distance
Given some density matrix ρ, the interaction distance 32 of ρ is defined as where F is the manifold of Gaussian density matrices σ, Here, H being quadratic means that it is a free-particle Hamiltonian, e.g., in second quantization, H = c † hc for some matrix h and some set of creation and annihilation operators {c † j }, {c j }, with j ∈ {1, ..., N }. D F , as defined in Eq. (5), measures distinguishability between a given density matrix ρ and the set of all free-particle density matrices, σ. Physically, the density matrix ρ can represent a thermal state of the system, in which case it is the standard Boltzmann-Gibbs density matrix when the system is in thermodynamic equilibrium at some temperature β = 1/T . On the other hand, ρ can also be a reduced density matrix which represents the subsystem A of a larger system in the pure state |ψ . In this case, ρ is obtained as the partial trace ρ A := TrĀ |ψ ψ| over the degrees of freedom of the sub-systemĀ complementary to the subsystem A.
The reduced density matrix is useful for characterizing properties of |ψ , such as the entanglement entropy of A, Since ρ A is readily available in numerical simulations, in what follows we focus on D F evaluated with respect to the reduced density matrix of the model's ground state.
There is a crucial simplification in evaluating D F as written in Eq. (5), which was shown in Ref. 32 using results from Ref. 40. The minimization over F is equivalent to where ρ k denote the eigenvalues of ρ in descending order (normalized such that k ρ k = 1), and where n (k) j ∈ {0, 1} is the occupancy number on the jth site of the kth element of a Fock basis with the energy j . The normalization Z ensures that k σ k = 1, and we assume that σ k are in the same (descending) order as ρ k , which is necessary to achieve a minimum in Eq. (8). 40 The utility of Eq. (8) is that the value of D F (ρ) can thus be determined solely from the information of the spectrum {ρ k }, also known as the "entanglement spectrum". 41 Comparing Eq. (5) with Eq. (8), we see that the minimization over all matrices σ ∈ F was traded for a minimization over scalars { j }. The latter is a much simpler optimization problem, as the number of parameters is expected to scale linearly with the system size. Thus, the problem becomes numerically tractable, as the computational complexity is only polynomial in system size N once the spectrum {ρ k } is known. 32 Note that D F is strictly bounded 0 ≤ D F ≤ 1, 42 and states that have D F = 0 can be expressed as Gaussian states in terms of some free-particle modes as in Eq. (9). This is, of course, true for unentangled (product states) in the computational basis, but it is also the case for certain entangled states. An example is the ground state of the Ising model in the transverse field, as discussed in the following section. Interestingly, unlike the lower bound, D F does not seem to saturate its upper boundit was conjectured that D F ≤ 3 − 2 √ 2. 42 Physical states that realize this upper bound of D F were identified as ground states of certain types of parafermion chains. 42 These states do not have a particularly high value of VNE, but the structure of their entanglement spectrum is as distinct as possible from that of free fermions, in the sense of Eq. (8).

IV. Preparing the ground state of the non-integrable quantum Ising model
In this section we present our main findings on the correlation between the success of QAOA state preparation and the interaction distance of the target state. As a toy model, we consider the one-dimensional quantum Ising model in the presence of both transverse and longitudinal fields, (10) where X i and Z i are the standard Pauli matrices on site i, and we assume periodic boundary conditions (PBCs) by identifying sites j + N ≡ j. The model is either ferromagnetic (FM) or antiferromagnetic (AFM) depending on whether the coupling of the first term is chosen to be +1 or −1, respectively. The Ising models in Eqs. (10) serve as a useful laboratory for studying a number of phenomena in condensed matter physics. [43][44][45]  The properties of the ground state of the model in Eq. (10) are insensitive to the sign of the FM/AFM coupling in the absence of the longitudinal field h z . However, once h z > 0, the phase diagram is substantially different for the two models. The FM model has a critical point at (h x , h z ) = (1, 0), while the AFM model has a critical line connecting the point (h x , h z ) = (1, 0) with the point (h x , h z ) = (0, 2). The critical line is not known analytically, but it has been determined numerically using density-matrix renormalization group simulations in Ref. 46. In both cases, the limit of purely transverse field (h z = 0) is particularly important. Along this line, the Hamiltonian is diagonal when written in terms of free fermions after performing a combination of the Jordan-Wigner and Bogoliubov transforms. 47 The phase diagrams of the FM and AFM models diagnosed by the value of D F in their ground state are shown in Fig. 2(a)-(b), respectively. The ground state of the Hamiltonian in Eq. (10) is obtained numerically using exact diagonalization, and its entanglement spectrum is computed by partitioning the system into two equal halves. From the entanglement spectrum, D F is evaluated by numerical optimization following Eq. (8) 48 . For both models, D F is found to be zero (to machine precision) when h z = 0, regardless of the value h x . Away from this line, D F is a sensitive indicator of interaction effects and changes by many orders of magnitude depending on the location in the phase diagram. For example, in the FM model, D F exhibits a sharp peak just off the free Ising critical point, (h x = 1, h z = 0). While the Ising critical point is described by a free Ising conformal field theory 33 and thus it has D F = 0, the properties of this CFT change dramatically once h z field is introduced. 49 This is consistent with the fermionic picture, where the h z field introduces long-range interaction between fermions after the Jordan-Wigner transformation, which makes the system's ground state highly interacting. Somewhat surprisingly, away from the critical point, the value of D F sharply decays to values as low as ∼ 10 −7 , even though the interaction is comparable in magnitude to other terms in the Hamiltonian. This implies that there are large regions of the phase diagram where the ground state of the system is effectively freefermion-like, even though the Hamiltonian itself is "interacting". On the other hand, the AFM model features a critical line that extends from the free Ising critical point (h x = 1, h z = 0). While D F = 0 at (h x = 1, h z = 0), the value of interaction distance progressively increases along the critical line towards the interior of the phase diagram -see Fig. 2 Next, we explore how to prepare the ground state of Eq. (10) using QAOA for arbitrary values of fields h x and h z . To this end, we have found it necessary to employ a M = 3-step QAOA protocol from Eq. (2) with the operators which satisfy H 1 + H 2 + H 3 = H, see the illustration in Fig. 1(b). The initial state of the protocol is taken to be the ground state of H 2 , i.e., all spins polarized along x-direction, |ψ init = | → . . . → . Since Pauli matrices are involutions, these operators H j satisfy e −i(θ+π)Hj = ±e −iθHj , so in what follows we restrict the representation of the angles θ to the [0, π) interval. We further restrict the angle θ p,1 associated with H 1 to the [0, π 2 ) interval, since The initial guesses for the angles were determined sequentially as p is increased, following the method in Appendix B1 of Ref. 31. For minimizations involved in both QAOA and D F we use a basinhopping algorithm with a Metropolis acceptance criterion, 50 as implemented in the Python package scipy.optimize.basinhopping. This is a global strategy that performs multiple minimizations, taking as the initial condition for the next minimization the stochastically perturbed result of the previous one. This allows us to avoid the local minima associated with the rugged landscapes of both QAOA and D F , as discussed further in Sec. VI. This, however, was not enough to completely eliminate local minima, and all the data presented here required two additional rounds of minimization. Each of these consisted in running the basinhopping algorithm across the phase diagram again, this time using as initial value for each point the optimal values of each of the adjacent points from the previously obtained data, and keeping the minimum value found.
Note that our protocol in Eqs. (11)-(13) is a generalization of the one considered in Ref. 11, which was restricted to the purely transverse field (h z = 0) and made use of a M = 2-step ansatz with only H 1 , and H 2 . In that case, both the Hamiltonian and the protocol conserve the total fermion parity, generated by P = i X i . This symmetry is broken once the z-field is introduced and the ground state acquires a non-zero magnetisation ψ| i Z i |ψ = 0. While it is easy to come up with a twostep protocol that does not conserve parity, we have not been able to find one that accurately prepares the ground state for general values of (h x , h z ), thus we introduced a third operator into the QAOA protocol.
In Figs. 2(c)-(d) we present results of the QAOA protocol across the phase diagram h x − h z . The color scale in Fig. 2(c)-(d) shows the infidelity 1 − f obtained after fixed p = N 2 steps of QAOA. We observe that this metric of success of ground state preparation looks remarkably similar to the behaviour of D F in Figs. 2(a)-(b). In particular, we recover f = 1 when h z = 0, 11 while the QAOA no longer finds an exact ground state when h z > 0. Nevertheless, it approximates the ground state very closely when D F is small. Once again, it is easy to see that that in this case there is no clear relation between QAOA's 1 − f and the VNE of the ground state. For example, in the FM model, the VNE should be largest at the critical point; and, since adding h z opens a gap in the spectrum, increasing this parameter should reduce the VNE, as its scaling changes from logarithmic divergence with system size to an area law. However, from the point of view of QAOA, we find precisely the opposite: it is harder to prepare the state with some small amount of h z compared to h z = 0.
Examining the optimal angles found at each point of the phase diagram of both the FM and the AFM Ising models when running the protocol in Eqs.(11)-(13), we found no continuous variation of the angles across the phase diagram of the kind in, e.g., Ref. 51. However, we found that the optimal angles θ 3,j had a striking tendency to be very close to multiples of π 2 (see Fig. 7 in the Appendix). This suggests that the Hamiltonian H 3 has a restricted role in the evolution, and that the symmetry could perhaps be broken in a simpler way. This property could be exploited by having the initial guess be close to multiples of π 2 through an ansatz, or by giving higher weight to regions close to these two points (0 and π 2 ) in the minimization algorithm. In Fig. 3, we show a scatter plot of D F vs. 1 − f from the data extracted from phase diagrams such as in Fig. 2, but using different numbers of QAOA steps p, indicated in the legend. In both FM and AFM models, we expect correlation between D F and 1 − f around p = N 2 . As p → ∞, we expect states to be exactly prepared and this correlation to break down, as in this limit QAOA protocols satisfying H = H 1 + ... + H M should have the same power as quantum annealing with an arbitrary schedule. 18 As p → 1, we expect that the variational method, in general, is not powerful enough for a correlation to emerge. However, in the FM model, we see that D F and 1 − f are correlated even at lower p. The minimization landscape at higher p becomes considerably more complex, which results in a larger spread of the data (see Sec. VI for further discussion of the minimization landscape). We compute the Pearson correlation coefficients for the data in Fig. 3 and plot them in Fig. 10 of the appendix; as expected, these seem to peak around p = N 2 .

V. Relation between QAOA and interaction distance
In Sec. IV, we have established numerically a correlation between D F and the success of QAOA protocols. This suggests that the protocol's success depends on how close to being Gaussian (in the sense of Eq. (6)) the target ground state is. In this section, we support these numerical observations by analytic arguments.
As mentioned in Sec. IV, the angle θ 3,j associated with the operator Eq. (13) is found to cluster around either 0 or π 2 . Now, note that a shift in π/2 of the θ 3,j part of the protocol results in an overall parity flip, as easily seen from the following sequence of identities: where | → , | ← denote eigenstates of X. This implies that, if we have the freedom of choosing either |→ ... → or |← ... ← as the initial state, we can restrict, without loss of generality, all angles θ 3,j to an interval of length π/2. Since these angles are clustered around 0 and π/2 as shown in Fig. 7 of the Appendix, they can all be mapped to be close to 0.
Next, note that since both initial states are product states, they are also Gaussian states. Moreover, the evolution under the unitaries generated by H 1 and H 2 maps Gaussian states into Gaussian states, while the evolution under H 3 spoils this property. However, for θ 3,j close to 0, the evolution under H 3 introduces only a small, perturbative deviation from a Gaussian state. This heuristically accounts for the high correlation of QAOA success with interaction distance of the target state, as the states prepared by QAOA are close to being free. As p gets larger, more perturbations are possible and the success of the QAOA increases. At a fixed p, the success of QAOA is related to the distance of the target state from the Gaussian state manifold.
We conclude that there is a practical limitation to the "natural" QAOA protocol proposed in Sec. IV, which was obtained as a Trotter splitting of the Hamiltonian (10) into its translation invariant components: the protocol is unable to prepare ground states that are far from being Gaussian (as measured by interaction distance). This limitation is fundamentally related to the probability spectrum of the target state, i.e., the eigenvalue spectrum of its reduced density matrix. Indeed, when performing QAOA using as a cost function the relative entropy 52 between the probability spectra of the trial state and of the target state, one finds heat maps and scatter plots similar to those in Fig. 2 (data not shown). Thus, there is a correlation between QAOA and D F , even though the former minimizes the overlap of two vectors, while the latter one employs a minimization using the probability spectrum of the subsystem's reduced density matrices. Note that if one performs QAOA using the relative entropy as the cost function and then computes the overlap between the optimal state and the target state, one finds that it is, in general, close to zero. On the other hand, if one takes the QAOA optimal states obtained using overlap as a cost function and computes its relative entropy with respect to the target state, we obtained heat maps similar to those in Fig. 2. This indicates that there is some "degeneracy" in states having the same relative entropy as the optimal state (with respect to overlap) in the space of states accessible through QAOA; however, this degeneracy is easily lifted when minimizing using overlap.

VI. Minimization landscape
In this section, we explore the minimization landscape of the optimization problem studied in Sec. IV for the AFM Ising model (we reached qualitatively similar conclusions in the FM model). The target state in the cost function is taken to be the ground state of the Hamiltonian at a set of representative points in the (h x ,h z ) phase diagram, S = { (0.1, 0.1), (0.1, 2), (1, 1), (2, 0.1), (2, 2)}. These points are drawn from regions of both "hard" and "easy" state preparation according to Fig. 2. Here, we use the rescaled relative energy, defined in Eq. (3), instead of the quantum infidelity. In Fig. 4 we first look at the probability distribution function for log in selected points S. We generate a sample of 10 4 initial θ angles, drawn from a uniform distribution in the [0, π[ interval. The distribution of log gives us insight about the structure of the landscape. A sharply-defined distribution of log is only obtained in the case where h x = h z = 0.1, with the peak at close to 0. The mean of the distribution shifts to large values of upon approaching the critical line, e.g., at h x = h z = 1. In addition to the shift of the mean, the distribution also develops multiple peaks corresponding to local minima. At other points in the phase diagram, such as h x = h z = 2, the separate minima form a smooth curve with larger variance. Finally, in some cases like h x = 2, h z = 0.1, we observe a clear bimodal distribution of the minima. Thus, the distribution of minima varies considerably across the phase diagram and, generally, has multiple peaks.
A systematic investigation of the nature of the landscape of a related minimization problem was performed in Refs. 21 and 22 using a discretized adiabatic state preparation protocol. In these works, the behavior of the minimization landscape was examined as a function of the total allowed time for the protocol. It was found that there are distinct "phases" associated with different intervals for the total allowed time. Particularly, at intermediate times, there is a glassy phase presenting with multiple clusters of minima where the minimization becomes difficult. Following Refs. 21 and 22, we have probed the nature of the minimization landscape in our models and using our QAOA protocol when the total time T (θ) = p i=1 M j=1 θ i,j is restricted. We impose this restriction in two different ways. First, we allow T (θ) to be less than or equal to some maximum total time T ≤max , which can be easily achieved by constraining the allowed interval for each θ i,j angle in our protocol. The second method is to demand T (θ) to be exactly equal to a given total time T =max . The results of these two approaches are contrasted in Fig. 5(a) and (b).
In Fig. 5(a) we see that, as expected, as T ≤max in-creases, decreases. Perhaps surprisingly, this occurs in a very clear step-wise fashion, suggesting that there are discrete values of T ≤max that show significant improvement in state preparation. By contrast, in Fig. 5(b) we see that as T =max increases, the behaviour of is more erratic, indicating that there are discrete, optimum values of T for which states can be prepared under this restriction. This shows that the protocol can not accommodate non-optimal values T =max , that is, there is no way for the protocol to continuously "stall" and wait, "wasting time" so as to emulate the last optimal value of T =max . The protocol can, however, "stall" in discrete values of π, due to the symmetry in the angles. A consequence of this seems to be the peaks and troughs pattern in the graphs in Fig. 5(b), which show some irregular pattern. This contrasts with the results in Ref. 21 and 22, which display an almost monotonically increasing success in state preparation as T =max increases.
Next, we took 500 random angle samples restricted to either T =max or T ≤max and ran QAOA with target state coming from the ground state at each of the representative points in S. In order to plot the highdimensional minimization landscape, we have used t-SNE, 53 a dimensionality-reduction algorithm for data visualization that embeds high dimensional data in a space with lower dimension while preserving the relative position of the data points. Performing t-SNE on these samples, we find that, for T ≤max , T =max < 1, there exists some clustering of minima, although some of the clusters are significantly less compact than others -see Fig. 5(c). For T ≤max , T =max > 1, the clustering rapidly disappears, first for the T ≤max restriction and then for the T =max restriction -an example of the latter is shown in Fig. 5(d) for T =max = 8. This indicates that QAOA, which usually does not place restrictions on the values of θ angles and therefore implicitly operates in the large-T regime, does not display a glassy phase in its minimization landscape as found for a different protocol in Refs. 21 and 22.

VII. Conclusions
In this paper we have investigated the preparation of ground states of non-integrable quantum models using QAOA. Our motivation was to identify physical properties of a state that can be associated with its preparation, thereby allowing us to bound the relative success of QAOA. While this task appears challenging for rigorous analytical treatment, we have numerically demonstrated a correlation between interaction distance and the success of the QAOA protocols in two variants of the quantum Ising model. This suggests that, in these models, states which are far from free, as measured by interaction distance, are harder to prepare, i.e., in order to prepare states with larger interaction distance, QAOA needs higher values of p to achieve the same degree of success as for states with lower interaction distance and lower p. We have also performed an analysis of the landscape associated with the QAOA optimization problem. We have found that there are several local minima associated with this landscape, though they are spread out and show no distinctive clustering. Limiting the total allowed QAOA time did not alter this landscape significantly for total time T 1.
One of the applications of our results is that theoretical insight into the closest free states representing the target state can be gained by using the experimentally obtained QAOA ansatz and setting the small θ 3,j angles to be zero. The absence of the glassy phase in the minimization landscape implies that the natural QAOA protocols constructed here do not lead to a NP-hard optimization problem and the time to find optimal angles should scale polynomially with the system size.

VIII. Acknowledgements
We acknowledge useful discussions with Chris Self, Christopher Turner, and Marin Bukov. This work was supported by EPSRC grant EP/R020612/1. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data.
Here we demonstrate that our findings from the main text also apply in a different model featuring three-spin interactions. The model is defined by the Hamiltonian where, again, X i and Z i are the standard Pauli matrices on site i, and we assume PBCs. The critical behaviour of this model is in the same universality class as the twodimensional classical three-state Potts model. The phase diagram of the model has been mapped out in Ref. 54 (see also Ref. 55 for further generalisations of the model). In the h x > 0, −3 ≤ h z < 0 region, it contains a critical line connecting (h x , h z ) = (0, −3) to (h x , h z ) = (1, 0). Below that critical line, there exists a threefold ground state degeneracy, while above it the ground state is unique.
The motivation for studying the model in Eq. (A1) is that its ground state is expected to be more strongly interacting and have a higher value of D F . For example, unlike the FM/AFM models in Eq. (10), the ground state of the model in Eq. (A1) cannot be obtained in closed form using the Jordan-Wigner and/or Bogoliubov transformation (apart from classical line, h x = 0). Moreover, the 3-fold ground state degeneracy in the ordered phase gives rise to an approximate 3-fold degeneracy of the entanglement spectrum, as generally found in "symmetryprotected topological phases". 56 This can be understood by picking a point (h x = 0, h z = −1), where the exact ground state of the system (with zero momentum under

translation) is given by
The corresponding entanglement spectrum is given by This is the type of entanglement spectrum that gives D F = 1 6 , a value close to the upper bound 3 − 2 √ 2. 42 The same spectrum is obtained in the Z 3 parafermion model at its fixed point. 57 An approximate 3-fold degeneracy in the entanglement spectrum persists throughout the ordered phase of the model, thus we expect the ground state throughout this phase to be more difficult to prepare using QAOA compared to the disordered phase.
The comparison between D F and QAOA for the model in Eq. (A1) is shown in Fig. 6. The QAOA protocol was chosen such that H 2 is defined as in Eq. (12) and H 3 as in Eq. (13), but here we use Note that this protocol also satisfies H 1 + H 2 + H 3 = H.
We have found that, like the two-spin Ising model, the success of the protocol also correlates well with interaction distance, as we see in the top row of Fig. 6. Here, as in Section IV, minimizations are done using a basinhopping algorithm. Moreover, we find correlation between D F and 1 − f for several values of p, as shown in the bottom row of Fig. 6. As before, the data in the bottom row of Fig. 6 was obtained by sampling across the entire phase diagram in the top row of Fig. 6. It worth noting that we can prepare the ground state in Eq. (A2) exactly by choosing the protocol while the initial state is the ground state of H 2 in the sector with magnetisation −N/3, as this is the sector where the states {|↑↓↓ ... , |↓↑↓ ... |↓↓↑ ... } live. It can be verified that this protocol prepares the exact ground state in Eq. (A2) in N/2 steps. Moreover, supplementing the protocol with a third operator, H 3 = − j X j , leads to good results across the entire phase with the 3-fold ground-state degeneracy. However, the infidelity 1 − f of the latter protocol does not capture the phase transition in a way that our protocol [Eqs. (12), (13), (A3)] does. Moreover, the initial state is more difficult to prepare in this case, unlike the product state of spins in our protocol.  Similar to the models studied in Sec. IV, we found no continuous variation of angles in the three-spin Ising model, and the angles θ 3,j tended to be close to multiples of π/2 (see Fig. 7). However, in this case the heuristic arguments of of Sec. V do not directly apply as the Gaussianity of the protocol is broken by the triple spin interaction term (A3). It is an interesting open problem to analytically explain the approximate Gaussianity in this case.

B. Additional data
This section contains some additional results that support the conclusions in the main text. Fig. 7 shows the distribution of the angle θ 3,j associated to the Hamiltonian (13) in the QAOA protocol. As claimed in Sections IV-V, these angles tend to be close to 0 or π 2 , which is clearly seen in the figure when j (and consequently p) is sufficiently large.
In Figures 8 and 9, we compare the heat maps depicting the QAOA infidelity across the phase diagrams of the FM and the AFM Ising models, respectively, for different protocol durations p. Comparing these heat maps to that of interaction distance in Fig. 2, we can clearly see that the correlation between QAOA and interaction distance is most pronounced at p = N 2 . In light of our analysis in Sec. V, this makes sense because we know that p = N 2 is necessary to prepare the free Ising ground state (h z = 0). Moreover, we see that 1 − f decreases monotonically as p increases; this is also expected, as in the limit p → ∞ we expect to be able to exactly prepare the states, as explained in Section IV.
Finally, we compute the Pearson correlation coefficients for the data in Fig. 6 and plot them in Fig. 10. As in Section IV, these seem to peak around p = N 2 . We see that the Pearson coefficient peaks around p = N 2 . For the FM and AFM Ising models, the coefficient remains high across the plotted range of p.