Analytical solution for low energy state estimation by quantum annealing to arbitrary Ising spin Hamiltonian

We point to the existence of an analytical solution to a general quantum annealing (QA) problem of ﬁnding low energy states of an arbitrary Ising spin Hamiltonian H I by implementing time evolution with a Hamiltonian H ( t ) = H I + g ( t ) H t . We will assume that the nonadiabatic annealing protocol is deﬁned by a speciﬁc decaying coupling g ( t ) and a speciﬁc mixing Hamiltonian H t that make the model analytically solvable arbitrarily far from the adiabatic regime. In speciﬁc cases of H I , the solution shows the possibility of a considerable quantum speedup of ﬁnding the Ising ground state. We then compare predictions of our solution to results of numerical simulations, and argue that the solvable QA protocol produces the optimal performance in the limit of maximal complexity of the computational problem. Our solution demonstrates for the most complex spin glasses a power-law energy relaxation with the annealing time T and uncorrelated from H I annealing schedule. This proves the possibility for spin glasses of a faster than ∼ 1 / log β T energy relaxation.

We point to the existence of an analytical solution to a general quantum annealing (QA) problem of finding low energy states of an arbitrary Ising spin Hamiltonian HI by implementing time evolution with a Hamiltonian H(t) = HI + g(t)Ht. We will assume that the nonadiabatic annealing protocol is defined by a specific decaying coupling g(t) and a specific mixing Hamiltonian Ht that make the model analytically solvable arbitrarily far from the adiabatic regime. In specific cases of HI , the solution shows the possibility of a considerable quantum speedup of finding the Ising ground state. We then compare predictions of our solution to results of numerical simulations, and argue that the solvable QA protocol produces the optimal performance in the limit of maximal complexity of the computational problem. Our solution demonstrates for the most complex spin glasses a power-law energy relaxation with the annealing time T and uncorrelated from HI annealing schedule. This proves the possibility for spin glasses of a faster than ∼ 1/ log β T energy relaxation.

I. INTRODUCTION
Quantum Annealing (QA) computing that has been implemented in hardware [1][2][3][4] searches for the ground state of an arbitrary classical Ising spin Hamiltonian H I (σ 1 , . . . , σ N ) by mapping the Ising spins to the zprojections σ k z of quantum spins 1/2 (qubits) and implementing an evolution with the Hamiltonian where f (t) is monotonically increasing from zero to a finite value and r(t) is monotonically decreasing from a finite value to zero functions of time; H t is the initial Hamiltonian whose ground state is easy to prepare, and k̸ =s a ks σ k z σ s z + k̸ =s̸ =r a ksr σ k z σ s z σ r z + . . . .
(2) The number of different terms in (2) can be exponentially large as H I can have arbitrary k-local terms that couple k spins directly with different coefficient a {k} .
The configuration of the ground state of H I encodes the solution of the target computational problems. Allowing only binary couplings in (2), this already includes many NP-complete problems [5][6][7][8], which means that many important QA problems that are usually formulated with a different from (2) target Hamiltonian, can be mapped to the model (1) with only a polynomial overhead. The integer number factorization and the Grover algorithm can be also formulated as QA problems with some H I [9,10]. Finally, in appendix VI A, we discuss that the result that we will present applies also to QA with different from Ising target Hamiltonians.
Today, accessible hardware for a large number, over 100, qubits uses only heuristic approaches to QA [11], for which the operator H t and the annealing schedule, f (t) and r(t), in (1) are not specifically tuned to the choice of H I . The QA protocol is chosen then mainly for simplicity to implement it in practice. Still, H t must not commute with H I , and have a large gap between the lowest eigenvalue and the rest of its spectrum. According to the adiabatic theorem, if the time-dependent parameters change sufficiently slowly, the system remains in the instantaneous ground state and thus transfers to the ground state of H I as t → ∞. Measuring the qubit polarizations σ k z , k = 1, . . . , N , we then obtain the desired configuration of Ising spins that minimize H I .
In real heuristic QA experiments, time is restricted by the coherence time of qubits, so the adiabatic regime is practically never achievable. Given the widths ∆E I of the energy band of H I , it is possible to perform a pseudo-adiabatic evolution with T ≫ 1/∆E I , where T is the achievable QA time. However, given that the gap between nearest levels of H I is generally δ ∼ ∆E I /2 N , i.e., exponentially smaller than ∆E I , the practical situation usually corresponds to the nonadiabatic regime T ≪ 1/δ. Hence, the experimentally accessible QA computing is inspired by a phenomelogical assumption that there are computational problems whose partial solutions, i.e., the low Ising spin energy states can be obtained during the nonadiabatic QA process faster than during classical computations. If this assumption is correct, the quantum coherent evolution can be used in combination with incoherent classical annealing during a longer time.
Whether this is true or not is hard to verify either numerically or analytically because we deal with a driven and nonadiabatic many-body dynamics. We still do not have definite answers on how quickly the useful information is gained during nonadiabatic QA computations, and whether there can be quantum algorithms that outperform classical computations during time that is accessible in practice.
The needed intuition can be gained from physics using the similarity of the complex Ising Hamiltonians with spin glass systems that correspond to randomly chosen couplings between spins. The glass phase appears at low temperatures and corresponds to extremely slow energy relaxation. Indeed, classical annealing simulations of spin glasses generally show a logarithmic residual energy dependence on time T of the temperature decay from a finite value to zero [12,13]: The transition to the glass phase is also expected for QA but the scaling of the residual energy with QA time is not clear. On one hand, quantum tunneling is more efficient than thermal fluctuations when overcoming spikes of a potential barrier. On the other hand, such barrier spikes can be bypassed in the multidimensional phase space of many qubits, whereas stochastic fluctuations are more efficient for transiting over shallow but broad potential barriers. Moreover, disordered quantum systems show purely quantum effects, such as many-body localization, that resist propagation of information inside a system. An example of this behavior is found in gammamagnets [14] -the models of arbitrarily many interacting spins that resist flipping even a single spin in response to arbitrarily strong and fast magnetic fields. Thus, there are arguments both in favor and against QA in comparison with classical annealing performance.
Early numerical studies found that QA leads to an inverse power of the logarithmic decay (3) as well, where T is the time of the QA protocol, but with a larger power β, and hence outperforms classical annealing [15,16]. However, later studies [17] claimed that this behavior might be a numerical artifact caused by time discretization, and the improvement of QA reduces only to a small finite offset in the time-continuum limit. If the system passes into a glassy phase, there are analytical arguments showing that QA has no advantage over classical annealing at all [18]. In any case, if slow energy relaxation (3) describes QA of spin glasses in the pseudo-adiabatic regime generally, the heuristic QA method becomes impractical for computations, apart from niche applications that avoid the spin glass behavior.
Here, we propose an approach to the question of QA performance arbitrarily far from the adiabatic limit. Namely, we introduce complete analytical solution of the Schrödinger equation with the Hamiltonian (1) having arbitrary H I part. The control protocol is not arbitrary but it depends on an arbitrary parameter, g, that describes deviation from the adiabatic regime of computation. This exact solution provides an insight into the accuracy of nonadiabatic computations for arbitrarily complex H I , which may not be possible to study by any alternative method. In particular, we will present a QA protocol that in the pseudo-adiabatic regime leads to a monotonic ε res ∼ 1/T power-law relaxation scaling with time T of QA, without any signature of a transition to a glass phase for arbitrary H I . This proves the superiority of QA over classical simulated annealing explicitly, and suggests that the energy relaxation can differ in classical and quantum spin glasses strongly, when it is assisted by external time-dependent fields.

II. SOLVABLE MODEL
First, let us show that the original model (1) can be rewritten in the form of a scattering problem that depends on a single time-dependent parameter g(t). In the Schrödinger equation, we switch to a new time variable Here, f (τ ) is positive, so s(t) is a single-valued function, which is growing monotonically with t. Moreover, since both f (t) and r(t) are changing monotonically with t, they are single valued functions of s: f (s) ≡ f (t(s)) and r(s) ≡ r(t(s)). Using that in (4), we find that (4) is mathematically equivalent to equation Note that since f (s) → 0 as s → 0, the initial conditions now become and since r(s) decays to zero as s → ∞, so does the redefined coupling g(s). Thus, the QA problem in (1) is equivalent to the QA model with the Hamiltonian where g(t) is decaying from an infinite value to zero. Next, if the goal is to study the accuracy of computations, one needs the probabilities of nonadiabatic excitations that are produced during QA starting from the ground state as t → 0. Here, we point to the fact that there is a fully solvable model that provides all excitation probabilities for evolution (6) with an arbitrary H I . This model has g(t) and H t , which satisfy the basic requirements for a QA protocol. Namely, and H t is the projection operator onto the state with all spins pointing along x-axis: This transverse Hamiltonian has been considered for QA problems previously in relation to the adiabatic Grover algorithm [10].
As t → 0 + , state |ψ 0 ⟩ is the ground state of H with an energy Since all the other eigenstates of H t are zero, |E 0 | is also the leading order energy gap to the rest of the spectrum of H as t → 0. Let |n⟩ be the state of an arbitrary configuration of all the spins with definite projections along the z-axis. For this state, where is the dimension of Hilbert space of N spins-1/2's. Thus, the matrix form of H t in the computational basis has identical exponentially small but nonzero entries. Let us also introduce the Ising energies where we reserve n = 0 for the ground state of H I , and assume that the state indices are chosen so that We will call n in ε n the number of excitations, because this index tells how many basis states have smaller Ising energy than the given state. Let a 0 (t), . . . , a N −1 (t) be the amplitudes of the basis states in the Schrödinger equation solution: For our QA protocol, the Schrödinger equation is given by a k , n = 0, . . . , N − 1.
(12) The solvability of Eqs. (12) follows from the fact that, after the Laplace transform, the N coupled equations reduce to a single first order ordinary differential equation in the Laplace transform of v, which can always be solved analytically (Appendix VI B). This model is a special case of a model that was solved by one of us [19]. Algebraic properties of this model were also mentioned in Refs. [20,21], but the relation of its solution to the QA problem has not been discussed before.
The analytical solution gives a simple formula for the probabilities of excitation numbers at the end of the evolution. If as t → 0 + the system is in the ground state |ψ 0 ⟩, the probability to produce n excitations as t → ∞ is given by Note that the final state probabilities do not depend on the particular expressions for the eigenstates |n⟩, and in this sense tell nothing about the ground state of H I . However, Eq. (13) gives complete information about the performance of the given QA protocol. For example, the probability to obtain the ground state is given by and the average number of excitations is These expressions simplify for a large number of interacting qubits N ≫ 1, for which N is exponentially large, and we can disregard p N in comparison to p. For g ≫ 1 we find p N ≪ 1, and P n follows the geometric distribution, with To provide an intuition about the properties of the distribution (13), we also note that if the energy dispersion of H I were linear, i.e., if ε n = nδ, then the distribution (13) would be a Gibbs distribution where 1/Z is a normalization factor and As the dimensionless parameter g is growing, the effective temperature (17) of the final excitation distribution is decreasing.

III. CHARACTERISTIC ANNEALING TIMES
Our choice of H t is hard to implement experimentally. The currently studied QA systems use a slowly changing transverse magnetic field with where σ k x are Pauli x-operators acting in space of individual spins. In later sections, we will argue that the model with schedule g(t) in (7) and H t from (8) is, for a certain large subclass of H I , optimal. Therefore, its  Table I. The colored vertical lines mark the corresponding times at which the excitations reach the halfway into their saturation, which verifies the same effective annealing rate for the protocols with equal τa.
solution can be used to learn about the entire strategy of using nonadiabatic QA for finding low-energy states. In order to explain the consequences of our solution, we must first introduce a method to compare performance of different QA protocols with g(t) ∼ 1/t α and different H t , but the same H I and the computation time T .
There is an additional time scale that characterizes the speed of QA. The operator g(t)H t has a bounded spectrum. Due to the exponentially large Hilbert space, this spectrum must have a high density region at some distance ∆E t from the ground state of g(t)H t . The Ising part H I also has a characteristic energy scale ∆E I , that is, the bandwidth of its spectrum ( Fig. 1, left panel). Since H I and g(t)H t do not commute, the resonant nonadiabatic transitions between the ground level of g(t)H t and the dense region of its spectrum become most probable near the time τ a , when the operators H I and g(τ a )H t become comparable (Fig. 1, left and middle panels), i.e., For example, for our solvable model (see Appendix VI C) where τ I = 1/∆E I is the characteristic time of dephasing that can be induced by the Ising part H I . We will call τ a the annealing time, in contrast to the total evolution time T that we will call computation time.
Any QA protocol must pass through the moment (19). Hence, τ a can always be defined consistently. We will say that two different protocols with power-law decays of g(t) and the same H I and T , have the same speed of QA if they also have the same τ a . The practically interesting values of τ a are restricted to the range The first inequality in (20) follows from the fact that the case of τ a < τ I corresponds to a strongly nonadiabatic regime, for which the gap in the spectrum of g(t)H t closes faster than the characteristic interaction rates of H I . We will say that one of the compared protocols is better if it produces fewer excitations, ⟨n⟩, when T /τ a = const ≫ 1 and the same characteristic times, τ I and τ a , are set for the different protocols. If a protocol is optimal, i.e., outperforms all other protocols at some imposed conditions on the QA schedule and for a certain class of H I , it must remain optimal after time-rescaling, t → λt, in the Schrödinger equation, because the latter merely means the change of timecounting procedure. It has been recently proved [? ] that if such a protocol exists, it must correspond to a power-law decay of the coupling: g(t) ∼ t a . We will use this result because it strongly restricts the class of the schedules that should be tested in order to prove the optimality. Here we also note that the solvable protocol has g(t) ∼ 1/t, which means that it may be optimal for some class of H I , which we will identify later.

A. Computational convergence rates
The analytical solution says that the probability to find the ground state configuration is growing linearly with τ a , however, starting from an exponentially small value. Thus, if we assume that g = τ a /τ I ≫ 1, then Hence, in order to make P 0 ∼ 1, we need the QA time The theory of simulated QA has previously produced various bounds on the rate of change of the coupling [13,22,23]. The simulated QA is a Monte-Carlo algorithm, which performance dependence on N and T can be different from the performance of the physical QA but both algorithms are interesting to compare. According to Ref. [22], to guarantee the convergence of the simulated QA for binary couplings in the Ising Hamiltonian, as t → ∞, to O(1) ground state probability, the field should change as where ξ is exponentially small for large N . Our solution agrees with this estimate. It shows the convergence of QA computing to the ground state in the adiabatic limit, during a finite non-polynomial in N annealing time (21). However, for a fair comparison, the result in Ref. [22] must be extended to the limit of maximal complexity of (2). At least the fact that the number of terms in H I can be exponentially large adds an extra exponential overhead on the Monte-Carlo algorithms, such as the simulated QA. The result (21) also shows that the generally exponentially hard computational problem requires exponentially large calculation time for a precise solution. Hence, computational difficulties reemerge in some form in different computational approaches. For specific problems, this annealing time can be generally obtained by the gap analysis and fine-tuning of the protocol for a specific H I . For example, if the minimal gap over the ground state scales as ∼ 1/ √ N , this imposes the same constraint for the annealing time τ a ∼ N . However, we stress that the gap analysis for complex H I can be very challenging, and a proper choice of the annealing protocol, g(t) and H t , requires individual tuning [24,25]. In contrast, our analytic solution applies for all H I with a fixed simple form of the annealing protocol.
The time estimate (21) can be compared to the one for a classical search algorithm that would identify the ground state of the diagonal matrix H I . If the entries of H I are random, there is no other way but to compare all eigenvalues, which requires N computational steps. Using this analogy, Eq. (21) suggests that τ I can be considered as an analog of the single computation time step and τ a is the analog of the full computation time in the classical search algorithms.
The modern attempts to develop QA hardware are largely based on a heuristic assumption that at moderate QA rates we can obtain considerable reduction of computational error rate even when the true ground state cannot be found. This effect would correspond to fast suppression of the average number of excitations, which for N ≫ 1 is given by As expected, ⟨n⟩ decreases with the growing annealing time τ a but nonexponentially and starting from an exponentially large initial value. The average excitation number ⟨n⟩ is related to the average energy For example, ⟨n⟩ is proportional to this energy for constant density of states. For spin glasses with random H I , the density of states is smooth and broad, as we show in (Fig. 1, right panel). In this case, for a broad range of annealing times, ⟨n⟩ and the average energy after QA are linearly related: ε res ∼ ⟨n⟩δ, where δ = ∆E I /N is the characteristic distance between nearest energy levels. Then, Eq. (23) means a surprising fact that the energy relaxation as a function of the annealing time follows a power law: rather than a logarithmic relaxation with growing τ a , which is found in classical annealing of spin glasses. The deviations are expected for truly slow QA because the density of states near the ground level generally decreases in comparison to the bulk of the spectrum. However, any power-law energy dispersion near the ground level, ε n ∼ n α , leads to a power law rather than logarithmic residual energy dependence on 1/τ a after averaging over the distribution (13).
Let us now discuss the fact that, formally, the computation time T in the solvable model is infinite but in practice it has to be finite. Let us set T to be proportional to τ a . The same scaling then would be found for the dependence of ⟨n⟩ on T if the deviation of the QA result at finite T from the exact solution is suppressed by a small parameter τ a /T . Numerically, we always found that ⟨n⟩ saturates for T > τ a close to the T → ∞ value, up to corrections of some order of τ a /T (Fig. 1, right  panel).
The following analytical arguments show that, indeed, a sudden termination of the protocol at finite T ≫ τ a produces a negligible difference from our analytical prediction. Using the Landau-Zener formula, the nonadiabatic transitions may not be suppressed during t > T for the states within the energy difference δε 2 ∼ |d∆E t /dt| ≤ g/T 2 . For spin glasses with smooth density of states the introduced deviations from ⟨n⟩ are suppressed, at least, by a factor O(τ I /T ), which has the same dependence on T as the ⟨n⟩ dependence on τ a but the factor 1/T is much   Table I. The Hamiltonian HI takes the form (2) with the coupling coefficients independently drawn from the standard normal distribution. The main figure and inset show the adiabatic (large g) and nonadiabatic (small g) regimes in log-log and semi-log scales, respectively. The solvable protocol (red points) always outperforms the other protocols for the same g.
smaller. For example, if we set τ a /T ∼ 0.01, then the deviations from the analytical prediction for ⟨n⟩ should not exceed ∼ 1%. Thus, we find the scaling assuming that τ a /T = const ≪ 1. Equation (26) is the main result of our article. We showed analytically that QA with the solvable protocol does not lead to a logarithmically slow relaxation for arbitrarily complex H I . In fact, the exact solution does not show any sharp changes in the relaxation curve, which are expected for the transition to a glass phase.
Below, we discuss other properties of the solvable protocol, which should be of interest for the heuristic QA hardware developments.

B. Degenerate ground state
The exponentially large QA time is needed for the solvable protocol to obtain the ground state only if this state is nondegenerate. We consider now the case with the ground state degeneracy. Let ε 1 = . . . = ε M −1 = ε 0 . Summing the first M equations in (12), we then find that the superposition |m⟩ is coupled to any state |n⟩, where n ≥ M , with a larger coupling g √ M . All other orthogonal superpositions of the Ising ground states then decouple and have zero probability to be at the end of the evolution. This means, in particular, that all Ising ground state configurations have equal probabilities to be found correctly.
The solvable model in Appendix B of Ref. [19] (see also appendix VI A) is applied even when all N states are coupled to each other with different independent N parameters. Thus, the modification of the effective coupling to state |+⟩ is still described by the exact solution in [19], which leads to the probability of the final state |+⟩: whereas the probabilities of the energy excitations do not change. This gives us an estimate for the time to prepare the state |+⟩ with probability P + ∼ 1: If M is large, e.g., this leads to an exponential speedup for extracting nonlocal information that can be obtained from measurements on the prepared superposition |+⟩.
For example, suppose that all excitation energies of H I are random positive and ε 0 = . . . = ε M −1 = 0 appear periodically, so that, when sorted in the known standard computational basis, they correspond to the eigenstates |x 0 + rT ⟩, where x 0 and T are integers, such that x 0 < T ∼ log a e N ; r = 0, 1, 2, . . ., and N /T is also an integer. This corresponds to M ∼ N /log a e N , so during the QA time of an order τ a ∼ τ I log a e N the solvable protocol prepares a state of the qubits as a symmetric superposition: The Quantum Fourier Transform then can be used to change this state into a superposition of the states |k⟩, where k is the integer multiple of N /T . Finding only two different k, one can then find their greatest common divisor by classical means, and thus determine the period T faster than by classical means. The possibility to solve the period finding problem on a qunatum computer is an essential ingredient in many quantum algorithms, such as the Shor's factorization algorithm. An important step in such algorithms is to find a symmetric superposition of equal energy eigenstates of a quantum function that has a high degeneracy of eigenstates in the entire phase space. Such a function can be usually encoded in the target Hamiltonian H I and thus one of its eigenstates can be found using QA. However, it is clear from our solution why such algorithms are hard to implement with the heuristic protocols, such as with the transverse field (18). This field couples different Ising ground states with the higher Ising energy states differently. Hence, even if we assume that the ground state can be prepared quickly, it will appear generally in a nonsymmetric superposition where the coefficients C r have not only different absolute values but also different phases which depend on all parameters of H I . Hence, further manipulations, such as making the Quantum Fourier Transform, may not provide a desired effect on this state, which is needed to complete the algorithm.

C. Effectiveness of the solvable protocol in the limit of maximal complexity of HI
The annealing protocol in our solvable model is unbiased in the sense that the amplitudes a n (t) (12) do not depend on the specific structure of the basis states. This is not the case for the protocol with a transverse field [26], which couples directly only to the basis states whose net spin polarization differ by ±1. Our protocol is also unbiased in the sense that degenerate ground state configurations as a symmetric superposition couple to the other states equally, which results into equal probabilities to find such ground states of H I at the end (section IV B).
Moreover, the statistical learning theory [27] says that direct approaches, which avoid the work with irrelevant information, should be favorable for learning algorithms. This is partly addressed by our finding that the final state probabilities obtained by solving Eq. (12) are independent of the precise values of ε k , i.e., the transition probability to any state |n⟩ depends only on how many other states have smaller Ising energies. For example, the probability to find the ground state does not depend on the choice of H I at all. This independence of the scattering probabilities of certain basic parameters is shared by all integrable models with time-dependent Hamiltonians [20] but is not expected otherwise. Hence, it must be unique for g(t) ∼ 1/t annealing protocol because other g(t) are not among the known solvable models with arbitrary H I . This property means that our solvable proto- col does not produce irrelevant information about specific values of ε k , as needed because only the ordering of these eigenvalues matters for finding good approximations to the ground energy.
Such properties altogether are unique among the possible QA protocols, which suggests that the solvable protocol, for some types of problems, could be favorable. Due to the universality of the analytical solution, if true, this should be true for the most complex form of H I . Thus, let H I be the sum of all possible terms in (2) with independent random coefficients a {k} . Such a stiff limit reduces to the problem of identifying the minimal value from an unsorted array of independent random energies ε n that are sampled from some distribution. For instance, for Gaussian random coupling coefficients, ε n form a Gaussian distribution as well (Fig. 1, left panel). Such a construction of complex H I does not favor any particular ground state spin configuration and even any systematic correlations between excited states. Hence, it is expected that the low energy states are estimated faster with a maximally unbiased QA protocol, which is our solvable protocol.
To test this hypothesis, we employ the result in [? ] that allows us only to compare the performance of the solvable protocol with a family of the protocols with a power law decay of the coupling, g(t) ∼ 1/t α , and identical for each protocol fully random H I , as well as the annealing time τ a and T /τ a . First, we note that the protocols with α < 1 produce definitely worse than ⟨n⟩ ∼ 1/τ a scaling for the excitations if we set T /τ a = const. This follows from the fact that even in the adiabatic approximation the term H t /t α mixes any Ising eigenstate with other states within the window of energy ε ∼ 1/t α . Hence, a sudden termination of such protocols at finite time T cannot resolve the states within the energy window that scales as 1/T α , which decays slower than 1/T .
For α ≥ 1, we resort to numerical investigation. Figure 2 compares numerically calculated final ⟨n⟩ for different protocols at N = 12 and the Hamiltonian (2) with randomly chosen all possible couplings. For large g, which we define for all protocols as g ≡ τ a /τ I , the excitation number decays as a power law. For any g and N , our analytically solvable model (Protocol 1) always outperforms the other protocols, although all of them show scaling similar to 1/g for large g. In numerous other tests (not shown), we found that all non-powerlaw schedules, e.g., with g(t) decaying exponentially, had much worse performance for the same values of τ I , τ a , and T , in agreement with [? ]. Figure 3 also shows the data that we used to extrapolate the results to larger N . For such interpolations, we always found that the solvable protocol produced smaller residual energy for the fully random Hamiltonian H I . Hence, as far as we could test numerically and extrapolate our results, the solvable protocol was, indeed, optimal for our comparison criteria and the most complex form of H I .
An alternative argument for the optimality of the solvable protocol for fully random H I follows from the estimate (23), which says that the performance of this protocol is actually the same as in the classical Monte-Carlo search. Indeed, a random search for the lowest eigenvalue has probability n max /N per step to pick up an eigenvalue from the first n max excitations. Hence it takes time τ ∼ N τ step /n max to find an eigenvalue with 0 ≤ n ≤ n max , where τ step is the time of one eigenvalue of H I computation and its comparison to a previously found lowest value. This is precisely the estimate of Eq. (23), in which we identify τ a with τ , τ I with τ step and ⟨n⟩ with n max . Since our QA protocol has the same convergence rate as the classical Monte-Carlo search of the completely unsorted array, any improvement over its performance on H I with all random entries, either for the full or the partial search, would mean the quantum supremacy that does not rely on hints such as the oracle in the Grover algorithm, which is believed to be impossible.
Thus, our protocol gives an explicit example of heuristic QA computations leading to the same performance as for one of the known classical algorithms. This includes all possible H I with nondegenerate spectrum, and all possible time restrictions. As our QA protocol, the unbiased random search Monte-Carlo is the preferable choice for searching through a completely random array but then by classical means. This raises a question whether many other heuristic approaches, such as using the practically most accessible QA protocols without correlating them with a desired task, or post-processing the final state as in the case of the ground state degeneracy, have also the same performance for all possible tasks as certain classi-cal algorithms.

D. Avoiding the Bound
The limit of fully random H I represents the largest class of all possible computational problems (6). Classical optimization algorithms usually trade between good and bad performance in different applications, which is known as the "no-free-lunch" property. Although similar results are not known for QA, it is expected that the effectiveness of the solvable protocol for the big class of the most complex problems generally means that there are protocols that outperform it on simpler problems with more structured H I . Below, let us show several examples in support of this hypothesis.

Nonadiabatic Grover's algorithm
A well known example of a problem with simpler H I is the one that is solvable by the Grover algorithm. It prepares the ground state of an operator H I that has all but one zero eigenvalues, whereas the ground state energy is −1. Let η k = ±1, where the sign depends on whether this ground state has the k-th spin, respectively, up or down. Then, H I for the Grover's problem can be written as In comparison to the most complex version of (2), this Hamiltonian is much simpler. It depends only on N sign parameters, and it has considerable symmetry: changes of these parameters do not affect the spectrum of H G I . It is, indeed, known that the ground state of H G I can be found by adiabatic QA during time that scales only as N 1/2 [10]. Achieving this adiabatically requires a very fine-tuned choice of the schedule g(t). However, if our solvable protocol is not optimal for the structured problems there must be protocols that achieve better estimates for the ground state for the Grover's problem also beyond the adiabatic regime, and such protocols may not need to be very complex.
Let us show that this expectation is true. Consider the QA Hamiltonian where H t is given by (8). Due to the degeneracy of eigenvalues of H G I , the evolution equation (12) reduces to two coupled differential equations for the amplitude a 0 of the ground state and the normalized sum of the other amplitudes: Namely, The initial conditions, as t → 0 + , correspond to a 0 = 1/ √ N ≈ 0 and, hence, a + ≈ 1. The protocol that makes P 0 ≡ |a 0 | 2 ∼ 1 is obtained by immediately setting the schedule to a constant value and then letting the system evolve at such conditions during time One can verify that this makes P 0 ≈ 1 by noting that Eqs. (32) with condition (33) are equivalent to the evolution equations for a spin 1/2 in a transverse magnetic field, which rotates this spin. Condition (33) is needed to remove the component of this field that points along the spin axis. Time T corresponds to a rotation angle that switches between orthogonal states of this spin. Unlike the time of the solvable protocol with g(t) = −g/t, which scales as T ∼ N , the time in (34) scales as ∼ √ N , which is expected for the Grover's computational problem.
This efficient protocol to solve the Grover's problem is fine-tuned for H G I and cannot show good performance on other tasks. Identifying such algorithms for heuristic computations requires additional optimization steps, e.g., using the methods of the machine learning [28], which would correlate the annealing protocol to a given structured H I . Such methods, however, become inefficient in the limit of maximal complexity with fully random H I because of the emergence of the barren plateau [29]. We leave the question whether they can produce a more efficient protocol than the solvable one in this limit open.

Models with limited connectivity
Another example corresponds to the systems with small connectivity between qubits in H I . It is expected then that a QA protocol that emphasizes interactions without many direct spin flips can achieve a better performance, such as the protocol induced by the decaying transverse field.
To test this hypothesis, we performed simulations for H I with limited connectivity ranges, i.e., a rangek Hamiltonian is of form (2) but only contains terms with at most k simultaneously coupled spins. This allows the control of the problem complexity by tuning the connectivity range. Our numerical simulations (Fig. 4) show that, for finite size systems of up to 12 spins and the transverse field protocol (18), the final excitation numbers always scale as a power law of g, i.e., ⟨n⟩ ∼ 1/g α ∼ 1/τ α a , and α increases with the decrease of the connectivity range of H I . At k = 2, α reaches the value 2. Figure 4 demonstrates the convergence of the performance to the universality domain of the solvable protocol with increasing complexity. In agreement with our expectations, as far as we could see numerically, the protocol with the decaying transverse field produced better performance on the structured problems than the solvable protocol, in agreement with the "no-free-lunch" property.
All such numerical tests had to be restricted to a small number of spins and a finite set of different protocols. Hence, our claims about optimality of the solvable protocol for the stiffest QA problems and the "no-free-lunch" consequences remain conjectures. Our numerical results only mean that such conjectures will be hard to disprove having access to the modern computers. If such conjectures are accepted as true, the performance of our solvable QA protocol can then be used for estimates of the lower bound of the time that is needed to achieve a desired accuracy of QA computation for the most complex problems at given time restrictions. Moreover,"no-freelunch" arguments suggest that for any less complex H I there should be QA protocols that outperform the solvable one. This gives a reference for characterizing efficiency of arbitrary problem-specific QA protocol.

QA superiority beyond adiabatic limit
Let us now return to the question whether QA computations in the nonadibatic regime can provide a better performance, in terms of scaling with the number of qubits, than the adiabatic quantum computations for the same problem. Our solvable protocol, as well as the nonadiabatic Grover protocol in section IV D 1, do not show this feature, as their performances scale equally with the adiabatic QA. Generally, this may not be true.
Here, we note that there is one more solvable model of QA that can be used to explore the scaling of τ a (N ) for a specific simple H I : Consider that is subject to a nonlocal constraint N k=1 σ k z = 0. Let us assume that |ε k | are of the order ε. The ground state of H ε I has N/2 spins pointing up. They correspond to the smaller half of ε k values. The other N/2 spins point down. Here, H I is parametrized by only N numbers ε k . Naturally, a wise algorithm should not look through all 2 N eigenvalues of H I but rather learn those parameters.
Due to the constraint, the ground state of (35) has zero total qubit polarization. To find this state, one can use FIG. 4. The scaling rate, α, obtained by fitting numerically obtained excitation numbers with ⟨n⟩ ∼ 1/τ α a , for various interaction ranges of the Hamiltonian HI and the transverse field protocol (Protocol 2 in Table I). The inset shows a comparison for the typical power-law decay for such fits of the data for protocols 1 and 2 with N = 12 spins under a range-2 (only binary spin-spin couplings) Hamiltonian.
the protocol with H t that also has the ground state with zero initial total spin [30]: As t → 0, the ground state energy of H(t) is separated from the dense region of g(t)H t near zero energy by ∆E a ∼ gN (N −1)

2t
, and the H ε I bandwidth scales linearly with N : ∆E I ∼ εN . The exact solution of this model was found in Ref. [30]. It says that the ground state is determined if g ≈ 1.
Using our definition of the annealing time in section III, we can now compare the performance of such QA computations with the performance of classical algorithms for the same problem. We find for the model (36) that g ≈ 1 corresponds to The same solution in Ref. [30] also shows that if we need only a partial search by allowing a fraction α ≪ 1 of mistakes, i.e., allowing αN spins pointing in a wrong direction, then it is sufficient to choose g ∼ 1/(N α), i.e., the computation time reduces by a factor ∼ 1/(N α), so in our notation τ a ∼ O(1/α).
If we can call H ε I as an oracle that returns an eigenvalue of any given Ising spin configuration in one time step, then learning the ground state classically here is equivalent to sorting the array of numbers ε k , which can be done, either fully or partially but for the same memory restrictions, in ∼ N logN steps. The partial QA solution thus has a better N -scaling than both the best available classical algorithm and the complete solution in the adiabatic limit. This example supports the speculations that a hybrid approach that involves a moderately fast QA step combined with a subsequent classical relaxation may improve the search for the true ground state.

V. ESTIMATES FOR PHYSICAL TIME OF COMPUTATION
The hardware to implement heuristic QA with simple annealing protocols exists [31][32][33][34][35][36], but its tests on specific problems gave contradictory results. There are claims for superior performance of QA in some instances [37], but achieving scalable quantum supremacy [11] using QA is still far from conclusive.
Let us estimate the performance of our solvable protocol at the current level of technology. The coupling energy of a single qubit to the rest of the quantum processor is physically restricted to some value ϵ max . For example, for a superconducting qubit, a coupling larger than the superconducting gap may produce unwanted excitations outside the qubit phase space. The bandwidth for H I is then restricted by ∆E I < ϵ max N . Hence, τ I for N qubits is restricted by If we assume ϵ max = 10GHz as the upper bound for a superconducting qubit, then to find the ground state of only 20 qubits, from (20), we need at least time τ a ∼ 0.1µs, which is the typical upper bound on coherence time of such qubits. The required computation time τ a is growing exponentially with extra qubits, so chances to solve an optimization problem for more than 20 qubits with the modern level of quantum technology are quickly vanishing.
Let us also estimate how many spins, n c , are expected to find their correct directions with respect to the true ground state for the solvable protocol. Assuming that the N − n c spins end up in the random directions, using (23) we then find ⟨n⟩ ≈ 2 N −nc . In practice, τ a is restricted by the decoherence time τ dec , which leads to the bound Assuming, optimistically, ε max ∼ 10 GHz, N = 1000, and the decoherence time τ dec ∼ 10 −3 s, from (37) we find that only about n c ≈ 35 qubits can correctly find their ground state directions due to the purely quantum annealing effects. The true number is likely much smaller because we used optimistic estimates of the physical parameters, while quantum supremacy requires at least n c ∼ 50. One practical advantage of the solvable protocol that may justify the efforts to implement it in hardware may follow from the complexity to retrieve the H I eigenvalues. Namely, when the sorting problem is encoded in the Hamiltonian of spin projection operators the direct classical algorithm requires additional computation of eigenvalues of H I at each step, which can be exponentially long on its own for the most complex H I , but is not required during QA. To exploit this resource, one should create a small processor, with only ∼ 25 high quality qubits, but with H I that depends on ∼ 2 25 different coupling parameters.

VI. DISCUSSION
Finding the ground state of an arbitrary Ising spin Hamiltonian is generally an exponentially hard computatinal problem. Even harder, it seems, is to study a dynamics with a time-dependent quantum Hamiltonian that implements quantum annealing computation in the nonadiabatic regime. Nevertheless, we showed that a fully solvable model for the most general case of Ising spin interactions exists.
In other branches of physics, integrable many-body models have been very influential -often not for a particular experimental application but for the opportunity to understand the behavior of complex matter in the regimes unreachable to numerical simulations. Similarly, our exact solution produces an insight into both spin glass physics and quantum computing from an original perspective. Thus, we used it to set new limits on the computation precision and proved the better relaxation scaling of the residual energy for quantum over classical annealing computations.
Numerically, we found considerable evidence to our conjecture that in the limit of the maximal complexity of the computational problem our solvable QA protocol outperforms other protocols for arbitrary QA rate at identical conditions for the time of computation. Given also the "no-free-lunch" property of algorithms, this leads to a new conjecture that more structured computational problems can be solved by certain QA protocols faster than in our solvable model. We provided the arguments in support of this conjecture too. Hence, our analytical solution can serve as a reference for the performance that can be achievable in the nonadiabatic regime for arbitrary H I .
A currently discussed technical question, besides improving quantum coherence, is how to redesign the interqubit connections and the annealing protocol in order to improve heuristic QA [36]. It is often stated that the performance can improve if one-to-many qubit couplings are implemented in the Ising Hamiltonian, and if the annealing protocol has a simpler spectrum in order to make it less biased and thus reduce the effects of resonances that are specific to H t . Our results show that such approaches may not lead to the boost of the performance. In fact, the solvability of our model follows from a high symmetry that makes the solvable protocol maximally unbiased. We showed that this provides the advantage, over other protocols, only for the tasks with the maxi-mal complexity but not for more structured Ising spin Hamiltonians. Hence, by implementing stiffer problems in a hardware, e.g., by adding one-to-many qubit connections and preparing less biased QA protocols, we may only bring the complexity of the QA computations closer to the domain of our model's superiority.
Our findings suggest that the quantum annealing superiority, for a specific problem, over all classical algorithms should be searched either in small size processors but with combinatorially complex interactions in H I or among relatively simple-structured H I , with a polynomial number of parameters but a transverse part g(t)H t that is tailor-made for this specific computational task. It is thus important to understand how the QA performance depends on the correlations between H I and H t , and on the prepared correlations in the initial state for quantum annealing. In addition, due to the lack of the logarithmically slow relaxation, the solvable protocol may find applications in problems that are traditionally studied by simulated classical annealing, such as finding the ground states of interacting classical spins in micromagnetic calculations [38].

ACKNOWLEDGMENTS
This work was carried out under the support of the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, Condensed Matter Theory Program. B.Y. also acknowledges partial support from the Center for Nonlinear Studies.

A. Solution for QA model with arbitrary target Hamiltonian
The annealing problem is sometimes formulated so that the target Hamiltonian, H 0 , is different from the Ising Hamiltonian: Let us show that our protocol with g(t) = g/t and H t given by (8) is still solvable in the sense that we can write the probabilities of the final eigenstates of H 0 in terms of the parameters of H 0 . Suppose that U is the unitary operator that diagonalizes H 0 , i.e., H I = U H 0 U † is a diagonal matrix. The latter means that it can be written in the Ising form (2), and we can define the basis states |n⟩ where n is the index of the excitation, as in the main text. Let us define the state |ψ ′ ⟩ = U |ψ 0 ⟩, |ψ 0 ⟩ ≡ | →, →, . . . , →⟩.
In the basis |n⟩, the entire Hamiltonian has the form It is now almost the same as in the problem considered in the main text but the state |ψ ′ ⟩ is dependent on the matrix U . Hence, the matrix elements of the mixing part are given by (gH t ) nm = g n g * m , where g n = √ g⟨n|U |ψ 0 ⟩.
Thus, unlike the model in the main text, the mixing Hamiltonian gH t depends on N generally different parameters that depend on the eigenstates of H 0 via the matrix elements of U . Nevertheless, the most general form of the model that was solved in Appendix B of Ref. [19] includes this particular case. Thus, if we define the probabilities p n ≡ e −2π|gn| 2 then Eq. (13) for the excitation probabilities (see also Eq. (B13) in Ref. [19]) is extended to Returning to the original problem (6) in the main text, it follows from (40) that knowledge of a unitary transformation U H I U † , such that its action increases the overlap of the ground state with the state |ψ 0 ⟩, can be used to increase the probability to find the ground state.

B. Solution of the model
Following steps from Appendix B in Ref. [19], we perform Laplace transformation a n (t) = A e −st b n (s) ds, where A is a contour in the complex plane such that the integrand vanishes when A originates and escapes to infinity (Fig. VI1). Substituting (41) into (12), we find a first order differential equation with a simple solution for b n (s), which we substitute to (41) to find a n (t) = c A e −st −s + ε n N −1 n=0 (−s + ε n ) ig/N ds, where c is a normalization constant that is fixed by the initial conditions. Following [19], as t → ∞ this integral is evaluated using the saddle point method and suitable deformation of A into the paths that go around the branch cuts in Fig. VI1. This results in the analytical expression for a n (t → ∞) in terms of Gamma-function of the parameters. The excitation probability is then obtained from P n = |a n (t → ∞)| 2 , and using the properties of the Gamma-function.
C. Setting parameters of protocols to compare their performance First, we note that H with H t in (8) and H 0 t in (18) have the same ground states both as t → 0 + and t → ∞. For both of them, the maximum of the density of states is at zero energy. Hence, for H t , ∆E t = g(t), and for H 0 t , ∆E 0 t = N g(t), where N is the number of spins. If, for the analytically solvable protocol with H t , we choose the time dependent form g(t) = g/t and fix the quench parameter g, then the annealing time is given by g/τ a = ∆E I , or where τ I = 1/∆E I . This also gives the meaning to the parameter g, that is, the ratio of the annealing time and the characteristic time of the dephasing by H I . For the transverse field protocol (8) with H 0 t , the same annealing time τ a in [43] is achieved if we set g 0 (t) = g N t .
Similar arguments for g 0 (t) ∼ 1/t 2 lead to g(t) = −g/(at 2 ), where a = ∆E I /g, as listed in Table I.