Adiabatic Spectroscopy and a Variational Quantum Adiabatic Algorithm

Preparing the ground state of a Hamiltonian is a problem of great significance in physics with deep implications in the field of combinatorial optimization. The adiabatic algorithm is known to return the ground state for sufficiently long preparation times which depend on the a priori unknown spectral gap. Our work relates in a twofold way. First, we propose a method to obtain information about the spectral profile of the adiabatic evolution. Second, we present the concept of a variational quantum adiabatic algorithm (VQAA) for optimized adiabatic paths. We aim at combining the strengths of the adiabatic and the variational approaches for fast and high-fidelity ground state preparation while keeping the number of measurements as low as possible. Our algorithms build upon ancilla protocols which we present that allow to directly evaluate the ground state overlap. We benchmark for a non-integrable spin-1/2 transverse and longitudinal Ising chain with $N=53$ sites using tensor network techniques. Using a black box, gradient-based approach, we report a reduction in the total evolution time for a given desired ground state fidelity by a factor of ten, which makes our method suitable for the limited decoherence time of noisy-intermediate scale quantum devices.


I. INTRODUCTION
Remarkable scientific progress in recent years has led to the first noisy intermediate-scale quantum (NISQ) [1] devices. Current NISQ devices are still limited by the number of qubits available, their gate fidelity and the maximum circuit depth. Much research effort is put into the investigation of algorithms for digital NISQ devices that could hold the promise of a (practical) quantum advantage [2,3]. Nevertheless, simulating large quantum many-body systems on these devices still remains challenging within the next years. Analog quantum simulators, however, are able to implement some quantum dynamics very efficiently. Furthermore, analog quantum simulators are especially well suited for probing universal features of quantum many-body systems. Several powerful concepts and experimental realizations of digital NISQ devices and analog quantum simulators have been put forward [4][5][6]. One particularly important problem is preparing the ground state of quantum many-body systems, relevant for a wide range of physics applications and closely related to combinatorial optimization tasks. Quantum computers already save an exponential memory cost compared to a classical computer in the representation of the quantum state, making a quantum device the natural choice to compute desired quantum states. Moreover, ground states are also of great importance for optimization problems as the solutions to combinatorial problems can be naturally encoded into a classical Hamiltonian. The ground state finding problem is known to be QMAcomplete [7], which translates roughly to the analogue * Corresponding author: Benjamin.Schiffer@mpq.mpg.de of NP-complete for a quantum computer. However, two different promising heuristic approaches have been established. First and most straightforward, there is the adiabatic approach. In order to prepare the ground state of a target Hamiltonian H T with a quantum adiabatic algorithm (QAA) [8], one takes an adiabatic path H(s) = (1 − s)H 0 + sH T . Starting from the ground state of a trivial Hamiltonian H 0 at s(0) = 0, the parameter s(t) is changed with time until the final value s(T ) = 1. Given a non-degenerate ground state along the path, the adiabatic algorithm is known to return the ground state for a sufficiently long preparation time T . However, T is a function of the spectral energy gap ∆(s) between the ground state and the first excited state along the adiabatic path and not known a priori [9][10][11][12]. Due to the limited decoherence times of current NISQ devices and analog quantum simulators, the preparation time T is of great relevance for the feasibility of the adiabatic approach. Different adiabatic paths can be constructed and a linear time schedule is generally not optimal in the sense of a minimal T . This is especially relevant when the gap ∆(s) becomes very small as it occurs in the presence of a quantum phase transition. Optimal adiabatic paths have been the subject of intensive research efforts both in the framework of shortcuts to adiabaticity as well as optimal control theory [13,14]. For the problem of an unstructured search, Grover-type speed-ups have been shown where the full spectral information about the problem is available [15,16]. Another approach to prepare the ground state is the quantum approximate optimization algorithm (QAOA) [17] which seeks to overcome the limitations of the QAA. It generalizes the QAA by splitting the total time into chunks, {T i }, where one alternates between H 0 and H T , and takes {T i } as the variational parameters. The QAOA includes the QAA in the sense that there arXiv:2103.01226v2 [quant-ph] 7 Apr 2022 exist T i = iT /L, where i = {1, . . . , L} and L sufficiently large, corresponds to the Trotter evolution of the QAA. However, it is expected that the QAOA can provide a speedup by choosing the parameters {T i } larger than in trotterized QAA. In fact, the QAOA belongs to a wider class of variational quantum algorithms (VQA) in the spirit of the variational principle which is widely used in physics. While the quantum computer is used to prepare the states and perform the measurements, the optimization of the parameters is carried out classically. In practice, the performance of VQA can be curbed since the number of measurements necessary to estimate an objective function may scale unfavorably, and due to the presence of plateaus in the energy landscape, including noise-induced barren plateaus [18][19][20][21]. In this work, we propose the concept of a variational quantum adiabatic algorithm (VQAA) to find optimal adiabatic paths for high-fidelity ground state preparation, thereby combining the strengths of adiabatic state preparation and VQA. Akin to the QAOA, the times {T i } are treated as the variational parameters. However, the evolution in the different chunks is performed adiabatically, similarly to the QAA. Here, the VQAA differs from the fixed Hamiltonians found in the QAOA. The VQAA allows for a significant acceleration compared to the QAA with a linear adiabatic path, yet requires fewer parameters and measurements than the QAOA. We discuss different approaches to find such a parametrized optimal adiabatic path. The approaches are suited to different resource requirements, some making use of ancilla protocols to estimate the ground state overlap.
These protocols rely on controlled unitary evolution which leads to many benefits in the quantum computing setting. Besides being central to the canonical quantum phase estimation algorithm [22], it also allows to access overlaps between initial and final states for a given unitary transformation by looking at the ancilla only and enables spectral projections on the state of the system if one enables postselection on the ancilla [23]. The ancilla techniques presented here may prove to be useful tools by themselves for other variational algorithms. This is especially true for our protocol with two ancilla qubits, as this constitutes a case of nontrivial distributed quantum computing. Interconnecting multiple quantum devices using a coherent link is a promising path forward for the field of quantum computing [24]. We benchmark the different algorithms with a quantum Hamiltonian and up to N = 100 qubits, which is what is expected for the new generation of NISQ devices or analog quantum simulators. We take a Hamiltonian that is non-trivial (non-integrable) but for which we can simulate the action of a quantum computer classically using tensor network techniques. In the case of a small gap in the adiabatic path, in the presence of a phase transition, we report significant reductions in the required evolution time to reach a given desired ground state fidelity. For a chain of N = 53 qubits and a target fidelity of 90%, the VQAA is able to reduce the total evolution time by a factor of 10 compared to a linear adiabatic path. Having mentioned before that the spectral gap ∆(s) is a priori unknown, obtaining knowledge about the spectral gap ∆(s) can be as hard a problem as finding the ground state itself. However, due to the interconnection between the spectral gap and the performance of adiabatic state preparation, we are able to propose a form of adiabatic spectroscopy to find spectral gap properties. Performing adiabatic sweeps on the system and being equipped with either the ancilla protocols for ground state estimation or backward-time evolution, a measure for both the position and the smallness of the spectral gap can be obtained. A related scheme using backwardtime evolution, albeit initializing the evolution in a superposition, has been suggested in [25]. Knowledge about the spectral gap may be applied to gain important insights about quantum many-body physics, obtain the phase diagram of a physical system or formulate the optimal path for adiabatic state preparation. The structure of this paper is as follows. In Section II we outline our main results. Then, we discuss our approach to adiabatic spectroscopy in Sec. III, and in Sec. IV, we present two protocols for estimating the ground state at a given point along the adiabatic path using one or two ancillas. In Sec. V, we introduce the general concept of a variational quantum adiabatic algorithm for ground state preparation and discuss different specific algorithms with different resource requirements (Sec. VI). After that, in Sec. VII, the model for benchmarking our algorithms is described and we present results. Finally, we comment on the number of measurements necessary in our approach (Sec. VII D) and the impact of noise on our algorithms (Sec. VIII). In the Appendices, we give the necessary background to make the paper self-contained: a theoretical description of the adiabatic algorithm and QAOA as well as Bayesian inference for Beta-Bernoulli models which are relevant for performing hypothesis testing.

II. MAIN RESULTS
Our main results include (1) protocols for eigenstate closeness estimation, (2) a proposal for adiabatic spectroscopy, as well as a (3) concept for variational quantum adiabatic algorithms including a black box gradient-based method.
(1) Provided the possibility to implement controlled unitary evolution on a quantum state, where the dynamics are controlled by a single ancilla qubit, information about the closeness to the next eigenstate can be extracted from the ancilla. We show that for a quantum state |ψ = j ψ j |φ j and Hamiltonian H = j E j |φ j φ j |, we need to obtain α := j |ψ j | 2 exp(−iE j τ ) by measuring the ancilla. Then, for suitable τ , we can deduce that the state |ψ is an eigenstate of H only if |α| = 1.
By a selfconsistent argument and suitable construction, the closest eigenstate may be identified with the ground state.
(2) The profile of the spectral gap ∆(s) between the ground state and first excited state energy of an interpolating Hamiltonian H(s) is closely connected to the evolution time T required for adiabatic state preparation.
We analyze the evolution time T (s) , the adiabatic path is a linear interpolation from H(0, 1, 0) to H(J, 1, 1). The minima of the different curves correspond to position and smallness of the respective spectral gap. In the lower plot, −∆(s) −2 is shown as a comparison, obtained with DMRG. The match with the Landau-Zener scaling is not exact as higher order corrections from adiabatic perturbation theory are playing a role.
necessary to prepare the ground state of H(s) with a given target fidelity. The ground state fidelity is computed either by time-evolving both forward and backward and measuring the fidelity with the initial product state, or by making use of the ancilla protocol which is generally expected to give superior results. After obtaining the data points for T (s), we interpolate the curve and compute the derivative. Corresponding to physical intuition, the curve of the evolution time features a strong increase when a small gap is crossed. The derivative ∂T /∂s can then provide the position and also a measure for the smallness of the spectral gap ( Fig. 1).
(3) In our work, we present the concept of the VQAA and give specific algorithms which attempt to improve over QAA. The main idea is to optimize the adiabatic path s(t) by performing a moderate number of measurements. In our setup, s(t) depends on a set of parameters, which 10 0 10 1 10 2 10 3 T are chosen through optimization procedures. The key ingredient is to use the overlap with the ground state as the figure of merit. The algorithms presented in this work make use of different protocols to estimate this overlap, and are also distinguished in the way the optimization is performed: either by intending to remain in the ground state at intermediate steps in the adiabatic path, or by optimizing for a maximum ground state overlap at the end of the state preparation (like in QAOA). The main feature of the method presented in this paper is that it requires fewer measurements. Our approach is well suited to be used in NISQ devices or analog quantum simulators by reducing the required preparation time and thus avoiding decoherence. Without knowing the adiabatic spectrum beforehand, the optimal adiabatic paths yield the desired ground state with high fidelity at only a fraction of the total evolution time of a nonoptimized adiabatic algorithm (Fig. 2).

III. ADIABATIC SPECTROSCOPY
We seek to gain insights about the adiabatic spectrum of the Hamiltonian The Hamiltonian H(s) describes the adiabatic evolution from the ground state of a trivial Hamiltonian H 0 at s = 0 to the ground state of H T at s = 1 where s = t/T is parametrized time. The probability that the system transitions out of the ground state in the course of the evolution is intimately connected with the size of the spectral gap ∆(s) = E 1 (s) − E 0 (s) between the ground state and the first excited state energy. For each s i in a set of data points {s i }, i ∈ {1, . . . , r} and r being the resolution of the spectroscopy, an efficient root finding algorithm is used (e.g. using a bisection algorithm on the time variable T i ) to find the evolution time T i which reaches a given target ground state overlap O T (see Suppl. for a simple algorithmic description). If a small gap is crossed in the adiabatic path, e.g. at s * , this will correspond to a large increase of the {T i } around and after this value s * . Hence, if we observe such a rise in the required evolution times {T i }, we conclude that a local minimum in the spectral gap ∆(s * ) must be present. In this manner, the position of a small spectral gap can be obtained. Moreover, the smallness of ∆(s * ) is related to the steepness of the increase in the {T i } around s * . For the Landau-Zener model, we establish the relation ∂T (s)/∂s ∼ 1/∆(s) 2 around the minimal gap (App. A). Since the Landau-Zener approach is a toy model to qualitatively study the property of the spectral gap in adiabatic algorithms, we presume that a similar scaling could persist in a more general sense.
In order to obtain a measure for the ground state overlap for given s i , we propose two different approaches. The first approach is not making use of the ancilla protocol and is thus very simple to implement. The initial ground state |ψ 0 at s = 0 is adiabatically time-evolved forward with evolution time T j , implemented by the unitary operator U 0→si ( T j ). Here, the T j denote the probing values in the search. We seek to determine the forward time so that we obtain the evolution time T i which succeeds in reaching the given overlap O T . Next, the interactions are reversed in a backward-time evolution from s = s i to s = 0 implemented by U 0←si ( T B j ). The backward time T B j T j is chosen to be larger than the forward time evolution. This allows for a trivial measurement of the ground state overlap at s = 0 as the state |ψ 0 is a product state. Once a T j is found for which O j ≈ O T , we set T i := T j and proceed with the next data point s i+1 . Note that it is also sensible to use a large constant evolution time for the backward path. We extend this approach by noting that if a small spectral gap ∆(s * ) has been found and we continue to probe for those s i for which s i > s * , it is economical to make use of the spectral information already obtained.
In the spirit of the VQAA, the adiabatic schedule used in the spectroscopy should be altered so that the evolution around s * is performed very slowly. In this manner, multiple gaps in the adiabatic spectrum could be investigated. This ancilla-free approach provides the least stringent requirements on the experimental setup in a NISQ framework. However, with some additional effort, there is the possibility with to directly estimate the ground state overlap without going back to the initial state at s = 0. This second approach leverages on the ability to obtain a measure for the eigenstate closeness through an ancilla protocol. For every data point in {s i }, we obtain the ground state overlap for different T j directly: with |GS si the ground state at s i . Now just like in the first approach, we search for the T j such that O j ≈ O T and set T i := T j . Hence, this form of adiabatic spectroscopy can provide a tool in order to experimentally gain knowledge about the adiabatic spectrum. We simulate this technique classically with N = 53 qubits for different spectral profiles and present the results in (Fig. 1) and Sec. VII B.

IV. PROTOCOL FOR GROUND STATE CLOSENESS ESTIMATION
A very important ingredient of our algorithms is the ability to estimate the overlap with a non-degenerate ground state.
In this section we introduce two suitable protocols using ancillas and comment briefly on their benefits. The ancilla protocols only require the ability to implement unitary evolution controlled on a single ancilla and to perform measurements on the ancilla. Controlled unitary evolution is a valuable ingredient for quantum computing, e.g. in the well-known quantum phase estimation algorithm [22]. Using ancilla measurements for ground state preparation has been investigated previously and can be used to construct spectral projection operators similar to the quantum Zeno effect if one allows for postselection on the ancilla [23]. In our work, we utilize the ancilla in order to extract information about the eigenstate closeness which can then be used in the different applications of the VQAA in a highly versatile manner.

A. Single-ancilla protocol
We first present a protocol for eigenstate closeness estimation using a single ancilla in the |+ state [26]. The quantum circuit shown in (Fig. 3) can be used to compare the overlap of an initial quantum state |ψ with the state after the unitary evolution |ψ evo For a fixed Hamiltonian H = j E j |φ j φ j | with the unitary we write the normalized state |ψ = j ψ j |φ j in the eigenbasis of H. Then, we analyze the quantum circuit for this choice of U . We obtain for the density matrix of the ancilla qubit after the controlled unitary evolution (App. B 1) For a quantum state |ψ that is an eigenstate, the rank of ρ a will be 1. Due to the specific structure of ρ a , only one off-diagonal matrix element is needed in order to determine the rank of ρ a . We note with (Eqn. B8) that For suitable τ , |ψ evo is an eigenstate only if |α| = 1. The time τ needs to be chosen so that the complex summands of α with non-vanishing amplitude do not have approximately equal phases. In this unlikely case of matching phases, we would see constructive interference so that |α| = 1 could be true even if |ψ evo is not an eigenstate. Visualizing the summands of α on a complex plane (Fig. 4), this becomes rather intuitive. The choice for τ is related to the spectrum of H. For the sake of this argument we assume that the overlap with the ground state and first excited state are the only other non-vanishing overlaps. Then, an arbitrary τ would correspond to choosing an l ∈ Z in τ = πl/∆ at random (App. B 3). For l 1, the probability of choosing an odd value of l is approximately 1/2. Therefore, by testing several random values of τ ∈ O(∆ −1 ) it is possible to deduce information about the system whether it is in a mixed state or an eigenstate with high confidence (cf. the Chernoff-Hoeffding bound in Sec. VII D). Through a self-consistent argument we can conclude that the main contribution to α comes from the ground state, provided we remain nearly adiabatic throughout the path. We can bound the maximal error in |α(τ )| 2 by computing the average for up to large values of an uniformly distributed τ as Then, we obtain FIG. 4. Exemplary eigenstate clock featuring the complex summands of α for τ =1. Here, the term corresponding to the ground state |ψ0| 2 exp(−iE0τ ) and the first excited state term |ψ1| 2 exp(−iE1τ ) lie quite closely together in phase which leads to problematic constructive interference. However, as the pointers rotate with their respective eigenenergy, they will be well-separated for suitable larger τ . Note that in this instance, the summands corresponding to higher excited eigenstates are taken very small and not visible.
for the ground state population |ψ 0 | 2 . In practice, terms in (Eqn. 7) corresponding to higher eigenenergies will be rather small and destructively interfere with each other. Therefore, this bound is not tight and much smaller errors are expected in an experiment (cf. App. B 4). Our protocol takes inspiration from semi-classical approaches to the quantum phase estimation algorithm [27,28]. However, these algorithms also seek to determine to energy. As we argued above, the ground state overlap is better suited as the cost function for ground state finding algorithms. Therefore we devised this protocol that is oblivious to the energy value at any given point while being very simple to implement.

B. Entangled-ancillas protocol
Building upon the single-ancilla protocol, we introduce a protocol using two identical quantum systems with one ancilla each. The protocol is motivated by the intuition that the single-ancilla protocol gathers information about the complex phase of ψ|ψ evo which is without practical use to us. Instead we would like to estimate |α| directly, which this protocol achieves. We require the possibility to conduct an entangling measurement (i.e. a Bell measurement) between the two ancillas of the two systems. Such an entangling measurement could be performed with a microwave quantum link between two superconducting circuits [24]. This protocol constitutes an instance of distributed quantum computing for such a quantum network. Moreover, it is well-suited to existing NISQ devices lacking an all-to-all connectivity where the adiabatic evolutions could be implemented on separate but connected sub-graphs [29]. Our goal is to determine the purity of the ancilla. In general, there is the relation ρ pure ⇔ Tr[ρ 2 ] = λ 2 1 + λ 2 2 for density matrices, with λ 1 and λ 2 the eigenvalues of ρ. We write the density matrix and its square as where matrix elements irrelevant for our protocol are denoted with a dot. With the Bell state |Φ − = (|00 − |11 ) / √ 2 we construct a Bell measurement. The diagonal matrix elements of ρ 2 a may then be attained by considering a composite system where the second system has controlled negative time evolution (implementable by changing the sign of H). Then, the density matrix of the ancilla of the second system effectively corresponds to the transpose of the density matrix of the first ancilla ρ a2 = ρ T a1 . The composite system gives If |ψ is in an eigenstate, the density matrix of the ancilla ρ a is pure and the expectation of the |Φ − measurement is allowing for very low variance measurements when |ψ is in the vicinity of the ground state. This protocol is especially well-suited for hypothesis testing (cf. Suppl.).

V. VARIATIONAL QUANTUM ADIABATIC ALGORITHMS
Preparing the ground state of a Hamiltonian is a problem of great significance in physics with deep implications in the field of combinatorial optimization. While adiabatic state preparation is known to return the ground state for sufficiently long preparation times only, variational quantum algorithms require a very large number of measurements in the training phase. We present a toolbox for variational quantum adiabatic algorithms (VQAA). Our objective is to prepare the ground state of a problem Hamiltonian H T with high fidelity while keeping the number of measurements in the process as low as possible. We aim at finding an optimized profile for the adiabatic evolution H(s) from the ground state of H 0 to the ground state of H T . For reasons of completeness, we refer to the Appendix (see Suppl.) for a treatment of the adiabatic algorithm and its implementation. In order to find an optimized adiabatic evolution, we choose a positive resolution L ∈ N for this velocity profile and split up the adiabatic path into L chunks. Then, for every chunk i ∈ {1, . . . , L} an optimal adiabatic evolution time T i needs to be determined (Fig. 5).
In every optimization task, keeping the number of parameters to be optimized as low as possible is paramount. This is because every new parameter yields an additional cost in the number of repetitions necessary to make all the estimations which are required in the optimization of that parameter. In a VQAA, two different options are possible. One could distribute the total time budget T for the evolution from s = 0 to s = 1 onto the L chunks (e.g. evenly spaced) and optimize the chunk lengths. Or alternatively, one could distribute the chunks in a given way and optimize the T i . For both options, we present suitable specific algorithms with different resource requirements. An essential ingredient for the algorithms is the ability to obtain information about the closeness of the quantum state at a given point in the adiabatic path with the ground state.

VI. PRESENTATION OF THE ALGORITHMS
The concept for the VQAA allows flexibility in the question of whether the chunk lengths or the chunk evolution times T i are the parameters to be optimized. Also, different specific classical algorithms can be used for the optimization process. Here, we present three different algorithms for finding the chunk lengths for fixed total evolution time and one algorithm for finding the chunk evolution times for flexible total T . These algorithms have different resource requirements and can make use of different ground state closeness protocols. The results stated in the abstract have been obtained using the gradient-based black box optimization for fixed total time (Sec. VI A 3).

A. Fixed total time optimization
Our proposal for the VQAA allows to set a maximum total time for the adiabatic evolution, which is well-suited for the limitations of current quantum devices. Here, the total adiabatic evolution time T is allocated evenly between the chunks of the adiabatic path, so that The chunk lengths become the variational parameters to be optimized, effectively controlling the density of adiabatic steps.

Ancilla-free optimization for fixed total time
In order to optimize | ψ(s = 1)|GS |, we try to keep the loss in fidelity in the ground state overlap along the adiabatic path as small as possible. Chunk lengths are initialized with equal lengths so that the end positions of each chunk are at reproducing what we call naive QAA. The chunk lengths ares i = s i − s i−1 , with s 0 = 0. Then, an adiabatic evolution is performed from the initial trivial product state |ψ 0 = |ψ(s = 0) = |− ⊗N up to the end points of each chunk.
This adiabatic evolution is then time-reversed (at the same speed) by changing the sign in the unitaries.
The total forward and backward-time evolution between s = 0 and s = s i are described by W 0→si and W 0←si , respectively, with the backward evolution being slower than the forward evolution. By going back to the initial product state |− ⊗N , the implementation of this protocol is rather simple and allows for low variance ground state overlap measurements without ancillas. We denote the individual adiabatic evolution operators for chunk j from s j−1 to s j in time T j with the unitaries V sj−1,sj (T j ), so that we write Now, we compute for all i ∈ {1, . . . , L} with O 0 = 1. The consecutive ratios of the overlap are These R i correspond to the drop in ground state overlap with each next chunk along the adiabatic path. In order to find chunk lengths which correspond to a smooth decrease of the ground state overlaps, chunks i where the drop in ground state overlap is larger than the average are made smaller and vice-versa (where R j is below average, the jth chunk is made larger). The sum of the chunk lengths is kept normalized to one ( L i=1s i = 1). Then, new values O i are computed and the procedure repeats until convergence. Clearly, what is optimized here is not exactly the instantaneous ground state overlap because the reversed adiabatic evolution will accumulate extra phases distorting the results slightly. Nevertheless, this method can be a useful compromise between a protocol which is very simple to execute and can still yield improved results of an optimized adiabatic routine (Sec. VII C 2).

Forward only evolution for fixed total time
Improving over the ancilla-free algorithm, this algorithm makes use of the ancilla-based ground state closeness protocol (Sec. IV A) in order to estimate the ground state overlap at the point s i in the adiabatic path without backward-time evolution. Here, |GS si is the instantaneous ground state at point s i of the adiabatic path. The rest of the procedure is analogue to the previous algorithm. This algorithm enables us to try to keep quantum state along the adiabatic path close to the instantaneous ground state for fixed total adiabatic runtime.

Black box optimization for fixed total time
We present a black box optimizer routine which makes use of the gradient to find optimized chunk lengths {s i }.
Here, the chunk lengths are optimized in a quantumclassical feedback loop similar to typical variational quantum algorithms.
In our approach, the vector containing the chunk lengths is used as the input for a quantum black box (Fig. 6). The chunk lengths remain normalized to 1. Within the quantum black box, an optimized adiabatic evolution is implemented according to the current chunk length vector, and the output of the black box is the final ground state overlap at s = 1. The ground state overlap may be obtained, e.g. by making use of our proposed one-ancilla protocol (Sec. IV A). This value of the ground state overlap is fed into a classical optimizer which updates the input vector. Our cost function is the ground state overlap which we seek to maximize. We use the (quasi-Newtonian) bounded limited memory BFGS (L-BFGS-B) algorithm [30] or the gradient-free Nelder-Mead (or downhill simplex) algorithm [31] for the classical optimization. To make the setting of the optimizer more realistic for an experimental set-up, we fix the relative step size in L-BFGS-B to be larger than 1% of the chunk lengths {s i }. The feedback loop is repeated until convergence or until another suitable termination criterion is reached, e.g. a desired ground state fidelity.

Quantum Black Box
Ground state overlap Update FIG. 6. Quantum classical feedback loop for optimizing the chunk lengths of an adiabatic evolution while keeping the total time fixed. A parameter vector containing the chunk lengths is used as input for a quantum black box. The black box implements an optimized adiabatic evolution accordingly and outputs the ground state overlap. Making use of a classical optimizer, the chunk lengths vector is updated in order to maximize the ground state overlap.

B. Target fidelity profile with flexible total time
The concept of the VQAA also allows for flexible total time optimization. Here, the goal is to find an optimized adiabatic evolution which gives a final fidelity as close as possible to a given threshold the adiabatic path is split up into L chunks of length s/L with respective evolution times T i . The state |GS is the ground state at the end of the adiabatic path at s = 1. By smoothly interpolating from θ 0 at s = 0 to θ L at s = 1 the intermediate target thresholds θ i are set. Starting with the first chunk, we now search for the evolution time T 1 which suffices so that | ψ(s 1 = 1/L)|GS | θ 1 . Here, we can make use of hypothesis testing which is highly efficient in the number of measurements (cf. Suppl.). Search algorithms such as bisection methods feature exponentially fast convergence. The {T i } define the optimized adiabatic evolution. This algorithm is then able to follow a target fidelity profile with resolution L as closely as possible. We note that due to the limited decoherence time of NISQ devices, a maximal value for the T i can be fixed.

VII. BENCHMARKING AND RESULTS
For benchmarking purposes, we choose a non-trivial problem where we are able to simulate the evolution of a quantum computer with and without noise on a classical computer. Simulating the dynamics of a large quantum system on a classical computer is, in general, very hard as the number of coefficients necessary for the classical description increases exponentially with the system size N . For this reason, we choose a gapped one-dimensional Hamiltonian as our system. In the algorithms we propose, the state is always close to the ground state, i.e. has few excitations only. The ground states of gapped one-dimensional systems can be described efficiently using a matrix product state (MPS) ansatz. We give a brief outline of the tensor network techniques [32] we use to simulate our algorithms classically in the Appendix (cf. Suppl.). Therefore, even though MPS cannot approximate time evolution in general, it is ideally suited for our problem since our states have few excitations.

A. Model
For benchmarking a simple, yet non-integrable quantum model, we use the translationally invariant Ising model with transverse and longitudinal fields, hereafter referred to as the ZZXZ model adiabatic evolution from the trivial Hamiltonian H 0 = i hσ x i stays entirely within the paramagnetic phase, therefore no phase transition is crossed [34,35]. This is different for choices, e.g. J = {2, 3, 5, 7}, where the adiabatic paths crosses a second-order phase transition from the paramagnetic phase into the antiferromagnetic phase.
The adiabatic spectrum following a linear interpolation from H 0 to H T is discrete for the lowest energy eigenstates (for a discussion on smooth reparametrizations of the path, see Suppl.).

B. Results for adiabatic spectroscopy
We benchmark the adiabatic spectroscopy for the ZZXZ model on a qubit chain with N = 53 sites. For our model, the spectral gap ∆(s) has been obtained using DMRG methods (Fig. 8). In our simulations, we compute the ground state fidelity directly using tensor contractions, however, in an experiment, this information would be gathered using our ancilla protocols (Sec. IV). In an experimental setting, every time we would like to make a measurement to estimate |α(τ )|, we can choose to evolve uniformly at random for different values of τ between 0 and some very large K. Under these assumptions a simple analytical formula for the expectation value E 2 of |α(τ )| 2 can be given (App. B 4) and a target ground state fidelity of |ψ 0 | = 0.8 translates to In the adiabatic spectroscopy, we obtain the evolution times T (s) required to reach this target for a given value of s (Fig. 9). Clearly, there is a strong increase in T around the position of the minimal spectral gap for

C. Results for VQAA
Here we compare the benchmarking results obtained for the different algorithms and interpret their respective behavior. In particular we would like to highlight that there are two regimes where an optimized adiabatic path improves over non-optimized (or naive) QAA. When a phase transition is crossed in an adiabatic evolution, as intuition would suggest, optimal adiabatic paths concentrate the majority of the evolution time around the position of the spectral gap. We also examine the case without a phase transition (J = 1) and benchmark for systems with up to N = 100 qubits. Here, we observe that rotations in the low-energy eigenspace can significantly improve the performance of the quasiadiabatic evolution. Moreover, we find that leaving the instantaneous ground state can pay off in finding a better final ground state overlap. This seems to be especially significant for short total evolution times.

Black box optimization for fixed total time
The best results for fixed time are achieved with the black box algorithm. Even for large system sizes, this algorithm can be able to improve significantly over naive QAA for fixed total time. Good results are usually achieved for a handful of chunks already and the performance of the algorithm does not improve much further for more chunks (which would correspond to a higher resolution in the adiabatic velocity profile). This is the case both when a small spectral gap needs to be crossed, but also in the absence a phase transition. In the latter case, the black box algorithm very often returns an adiabatic path including one or even several very small chunks. These rotations can be understood in the following way. The set-up of the VQAA with fixed time allows some chunk lengths i to become very small. In these small chunks, however, there is still an amount of time T i spent. The evolution in such a chunk in the limit ofs i → 0 is implemented by a unitary which effectively changes the local phases of |ψ(s i ) thereby physically changing the quantum state. As |ψ(s i ) in our algorithm is expected to have a considerable ground state population and also small populations in the lowest excited states, the change in the local phases of |ψ(s i ) corresponds to a rotation in the low energy eigensector. Some rotations eventually become beneficial for the performance of the adiabatic routine. Rotations in the low energy eigensector have also been considered in the context of the Eigenstate Thermalization Hypothesis [36].
In the case of a phase transition, the improvement in the evolution time T between naive QAA and an optimized adiabatic evolution is much larger than without a small gap. For very few chunks only, target fidelities of over 90% in a system of 53 qubits were reached at around T ≈ 100. With a non-optimized adiabatic path, this would have required very long (≈ factor 10 longer) preparation times as depicted in (Fig. 2). The classical optimizer used for obtaining the data in this figure is the Nelder-Mead algorithm. In our simulations, Nelder-Mead proved to be more robust at large T at the expense of requiring more measurements. Being a gradient-free method, we expect it to behave better than the gradient-based L-BFGS-B for noisy systems. In (Fig. 10), the evolution of the chunk positions for an instance of the black box algorithm is shown for 53 qubits and J = 3. Here, the total time is T = 5 and L-BFGS-B was chosen as the classical optimizer due to its faster convergence in the classical simulations. The adiabatic spectrum for J = 3 features a small gap around s = 0.34. This can be approximately captured with three chunks only. Over a few iterations, the center chunk becomes smaller around the position of the gap, so that more time in the adiabatic evolution is spend there.
In the case without a phase transition, for J = 1, we find reductions of the total evolution time by a factor of over 3 when comparing with naive QAA. Here, we benchmark the black box optimization routine for large values of T and 100 qubits (Fig. 11). We note that for T > 40, the classical L-BFGS-B algorithm was in most cases not able to find an optimized adiabatic path. It was made sure that the reason for this was not due to memory-limitations in the optimizer, the relative step size or tolerance values for the termination of the algorithm. Increasing the number of chunks does not generally help in finding better optimized paths. However, sequential initializations of the optimizers can improve the performance. Instead of starting with naive QAA, the optimizer then begins with the chunk lengths of the previous, shorter-time optimizer instance. A typical black box algorithm instance (for J = 1, i.e. no small gap) is shown in (Fig. 12) where the small chunks can be easily observed around s ≈ 0.4. Indeed, the ground state fidelity is maximized and the energy minimized in five iterations only. In (Fig. 13), the fidelity between |ψ(s) and the ZZXZ ground state or the instantaneous ground state of H(s), respectively, is shown. The exact ground states have been obtained using DMRG methods. A discontinuity in the path of the optimized QAA is clearly visible for s ≈ 0.4 due to the implemented rotation. As the black box optimizer only has access to the ground state overlap at s = 1, it is agnostic to the actual curve of the optimized QAA for values below s < 1. It is interesting to observe that in fact leaving the instantaneous ground state can lead to better final results at s = 1, also discussed for adiabatic evolutions in [37] and in a related way in the context of diabatic transitions in QAOA in [38]. This feature can also be observed in a black box optimizer instance for a smaller system in (Fig. 14) where the optimized QAA curve gives a final fidelity of approximately 0.999 and recovers most of its ground state fidelity in the last 20% of the adiabatic evolution path only.

Converging fidelity ratios for fixed total time
Besides the black box approach to VQAA, we have presented other algorithms which aim to stay close to the instantaneous ground state along the adiabatic path. Both the ancilla-free and the one-ancilla method seek for convergence in the fidelity ratios between consecutive chunks. The one-ancilla method is cleaner in the sense that it directly uses the overlap and does not accumulate extra transitions and phases on the backward path (Sec. VI A 2). The one-ancilla method can achieve a smooth adiabatic path which remains as close as possible In this instance, albeit for N = 10, we see near perfect ground state preparation for T = 10 without increasing the total time budget. Even though the instantaneous ground state fidelity with H(s) along the adiabatic path is smaller, the final ground state fidelity is 0.999. to the ground state at all times for fixed total evolution time. However, the ancilla-free method finds adiabatic paths which effectively implement a rotation in the low energy eigensector. This property is also observed in the gradient-based black box method (Sec. VII C 1) and seems to be very relevant for preparing the ground state on NISQ devices. Therefore, which method will perform better is likely quite model-dependent for these two algorithms.

Target fidelity profile for flexible total time
The flexible total time algorithm with a given target fidelity with the ground state at the end of the adiabatic evolution is well-suited in the kind of instances in which staying close to the ground state is desired. As we observe in (Fig. 15) the naive QAA will occasionally leave the ground state leading to oscillations in the instantaneous ground state fidelity. Considering the Instantaneous ground state fidelity for QAA optimized with the flexible total time algorithm compared to naive QAA with the same time budget. The target fidelity at s = 1 was set to θ = 0.99 and linearly interpolated along the adiabatic path. The search interval for the time spent in each chunk was upper bounded such that Ti ≤ 20 ∀i.
adiabatic path split-up into several chunks, due to the changing Hamiltonian H(s) different evolution times {T i } for each chunk are necessary in order to achieve a given fidelity with the ground state at the end of the adiabatic evolution. Values for the {T i } strongly depend on the spectral gap ∆(s) between the ground state and the excited state as well as on the Berry connections (cf. Suppl.). We observe that our algorithm is able to reduce the total adiabatic evolution time compared to naive QAA in this toy example. This is because it uses the time available in a more economic way by spending much of the evolution time only when required to stay close to the ground state. One useful property of this algorithm is the fact that it is self-verifying in a sense that the hypothesis testing at the end of every chunk guarantees with high confidence that the ground-state fidelity is larger than a given value. Further developments of this algorithm can be envisioned where the {T i } are optimized to follow a more complicated profile. Moreover, the results of this algorithm serve as a good initial point for the gradientbased black box algorithm for fixed total evolution time. Here, we have set the chunk lengths {s i } to equal values. Adaptions of this algorithm with unevenly spaced chunks are also possible. Including additional chunks with very small chunk length could possibly provide a better performance of the algorithm by making it possible to implement rotations in the low energy eigensector at the expense of an increased number of measurements (Sec. VII C 2).

D. Cost of implementation
In every variational algorithm, the number of measurements needed in order to obtain satisfactory results is of upmost importance as it directly determines the feasibility of the approach. Here, we discuss the number of measurements necessary for each respective algorithms presented in this paper. We focus on the number of ground state overlap evaluations necessary in our classical simulations. To obtain the ground state overlap from |α|, a small overhead is required (Sec. IV A). Both the ancilla-free and one-ancilla fixed time algorithms converge with very few iterations even for large system sizes (N = 100). The resolution of the velocity profile in the adiabatic evolution is the number of chunks chosen. In our case, the number of ground state overlap evaluations was fairly low with #GS overlaps = #iterations · #chunks 100 (24) which sufficed for our simulations. The number of chunks are the number of parameters to be estimated and for the flexible total time algorithm this corresponds to the number of bisection searches required. As these search algorithms converge exponentially fast, roughly 10-20 ground state overlap estimations per search are usually sufficient for highest accuracies of the optimized {T i } values. In this algorithm, the two-ancilla algorithm is used which is suitable for hypothesis testing. Hypothesis testing converges exponentially fast as well, asking for about another 10 measurements for each ground state overlap estimation. Exemplary numbers for hypothesis testing are given in the Appendix (cf. Suppl.).
The black box algorithm depends on ground state overlap evaluations in order to estimate the gradient. The number of ground state overlap evaluations necessary in the optimization process is directly related to the number of chunks. Already for very few iterations, good results can be obtained with this method, so that for five chunks, we achieved good results for N = 100 and T = 1 with significantly less than 50 ground state overlap evaluations in total (App. Fig. 20). For a random variable X with values in [a, b],∆ = b − a and independent and identically distributed samples, the Chernoff-Hoeffding inequality [39] gives an upper bound on the probability to find the measurement deviating more than from its expectation value µ For Pauli measurements, we can have either +1 or −1 as results, so∆ = 2 and the number of measurements depends on the η and needed. Setting the precision to (1−(GS fidelity))/20 and the failure probability so that 1 in 2 experiments is successful with all estimations within the deviation , we estimate the number of measurements necessary to be of the order of 10 3 . The number of ground state fidelities required to reach very high fidelities > 0.9 at larger T is naturally larger than what stated above for T = 1. In (Fig. 2), the maximum number of ground state fidelity evaluations was set to 1000, while typically few hundred evaluations sufficed to find an adiabatic path with maximal ground state fidelity up to 10 −8 relative accuracy (i.e. the relative accuracy of the optimization process). Without doubt, this accuracy will be unattainable in current experiments and much less ground state evaluations give already very good results. In fact, we observed that for N = 53 in the case of a phase transition and T ≈ 100, when obtaining less than 200 ground state evaluations, the Nelder-Mead method will yield results of practically the same quality. These estimations suggest considerably low number of measurements even for optimizing the adiabatic evolution of large quantum systems. While the number of measurements seems to be the most relevant figure of merit to assess the cost of the methods presented here, we also include a short discussion of the number of gates required on different architectures of quantum devices. On analogue quantum simulators, for instance, the ancilla-free optimization method can be implemented natively with the only overhead being the additional parametrized adiabatic sweeps to find optimized adiabatic paths. In a gate-model architecture, for a qubit-chain with N = 53 sites and only allowing for next-neighbour interaction, we upper bound the total number of CNOT gates for one unit of time to be around 120 CNOTs/N . Here, we consider a decomposition of the unitary gates in the trotterized evolution of the numerical simulation. Actual gate counts will be significantly lower because the circuit will be optimized for the respective experimental hardware platform. In the single-ancilla protocol, an upper bound on the number of CNOTs per unit of time (τ = 1) including the required SWAPs for the chain topology is at around 2100 CNOTs/N . We note that with τ scaling as the inverse of the spectral gap ∆(s), the cost of the protocol is not excessively demanding, especially at the end of the adiabatic sweep (s = 1) where the gap is large. A more detailed discussion of the gate count is provided in the Supplement.

A. Noise in adiabatic quantum computation
The noise in current quantum devices severely limits the performance of many quantum algorithms [40]. Therefore, we discuss some important properties of noisy quantum adiabatic algorithms. A general inherent robustness of adiabatic evolution has already been established for some time [41], here we focus on a few points that are especially important to our method. In a gate-model quantum algorithm without error correction, a flipped qubit will in the worst case render the whole quantum computation nonsensical. This is quite different in an adiabatic algorithm as for physical instances low energy spectral lines are rare (App. Fig. 21a). While a flipped qubit in the preparation of the ground state in an adiabatic algorithm can also lead to a quantum state orthogonal to the ground state, the energy, however, of this orthogonal excited state will still be a very good approximation to the ground state energy. Intuitively, one bit flip corresponds to a single excitation of the system. We can therefore assume that errors increase the energy only by O(1) for fixed time, when flipped qubits are rare. Also, the position s 0 in the adiabatic path, where a flipped qubit occurs, is not critical. (App. Fig. 21b). This may seem somehow surprising, but it can be explained because a perfect adiabatic evolution suppresses all transitions between the eigenstates. It does not only apply to the ground state but also to excited states. Therefore, a noise-induced excitation will in the regime of an adiabatic evolution not lead to further deviations from the ground state energy. However, for time evolutions that are faster than an adiabatic evolution, which is in general the case on quantum devices limited in coherence time, we would generally expect light-cone spreading of noise through the spin-chain. Yet, in our simulations this was not observed to be problematic for the performance of the VQAA. In general, the noise behavior of adiabatic algorithms is encouraging as it suggests very benign noise features in these kind of algorithms making them a suitable candidate for NISQ devices.

B. Impact of noise in the presented algorithms
Here, we include a qualitative discussion about the expected performance of the algorithms presented in this paper in the presence of noise. We expect the adiabatic spectroscopy to be quite robust to noise as the information obtained using this method relies on multiple data points and a rather distinctive feature in the T (s) curve resulting from a small spectral gap. Noise effects will become stronger towards the end of the adiabatic evolution as noise accumulates in the circuit. However, as the results of this spectroscopy are quite pronounced, a qualitative description of the gap is likely only slightly impaired by moderate noise in the circuit. Regarding the VQAA algorithms, when considering noise, there are two main points to consider. First, noise can substantially impede the training phase of the algorithm, when the parameters of an optimal adiabatic path are being searched for. Second, in order to prepare the ground state with a desired high fidelity, an adiabatic evolution time T that is only a fraction of the T required for naive QAA suffices with an optimal adiabatic path. This may help strongly in suppressing errors. Following a target profile is especially tricky when noise comes into play. This is because noise strongly alters the required target profile. Therefore, concerning the target profile method, we do not expect this method to be very robust to noise, especially when finite size effects are playing a role. In the presence of noise, optimizing with regard to the final ground state overlap instead, is thus advisable. For the black box method, in classical simulations with noise, we observed the gradient-based training to be not too well-behaving. Convergence to optimized paths that improve over naive QAA was often impossible even when only few bit flips occurred in the quantum circuit. Classical optimization routines which are more robust towards noise seem to be asked for here, i.e. a classical optimizer that is combined with an error mitigation technique so that noisy outputs of the quantum black box can be corrected. We note, that gradient-free optimization methods are expected to be more robust to noisy environments and could provide better performance in an experimental setup than a gradient-based method. For this reason, we made use of the COBYLA method [42]. In (Fig. 16), it can be observed that small amounts of noise significantly impair the training process of an optimized adiabatic path. For small noise strengths p, however, the results can improve even over naive QAA. For simulating the noisy quantum circuit, 100 MPS samples were taken. We refer to the Appendix for more details on the noise model (cf. Suppl.). Besides noise in the adiabatic evolution, noise also is present in the measurement process. For both the oneancilla method and the two-ancilla method, we have benchmarked the black box routine numerically. The shot noise simulation is performed by finite sampling of m independent measurements and taking the average of these measurements. From the numerical data, we can conclude that for a large system around 10.000 measurements can reduce sufficiently reduce the shot noise. This is one order of magnitude larger than the estimations made earlier in (Sec. VII D). The benchmark for the one-ancilla method is shown in the main text, the plot for the two-ancilla method can be found in the Appendix. A further comment shall be specifically addressing the noise in the black box algorithm which makes use of a gradient to optimize the adiabatic path. In a recent work, it has been shown that the gradient in a variational quantum algorithm vanishes exponentially in the number of qubits N when the number of layers scales as poly(N ) [18]. These noise-induced barren plateaus severely hinder the scalability of variational quantum algorithms on NISQ devices. In our work, however, the algorithms either optimize the adiabatic path for fixed total evolution time (which includes the gradient-based black box algorithm), or have a maximum time budget in the case of the target fidelity profile algorithm. Thus, the (maximum) circuit depth is fixed in our approach which makes the results on noise-induced barren plateaus not directly applicable.

IX. DISCUSSION
In this paper, we present a toolkit for quantum adiabatic computation.
This toolkit includes a proposal for adiabatic spectroscopy, ancilla protocols for estimating the ground state overlap as well as the VQAA as a flexible yet powerful framework for variationally optimized adiabatic paths for high-fidelity ground state preparation. The adiabatic spectroscopy offers a straightforward approach to obtain information about an adiabatic spectrum. Our method relies on protocols to evaluate the closeness of an N -qubit quantum state |ψ to an eigenstate using controlled unitary evolution. By remaining sufficiently adiabatic throughout the evolution, a self-consistent argument applies, enabling us to identify this eigenstate with the ground state. We note, that the error in the ancilla protocols presented is smallest when being close to an eigenstate. The protocol which we propose requires the ability to implement controlled time evolution on a single ancilla qubit. Current technology already meets the requisites of our protocol: Conditional dynamics have been explored in the context of trapped ion simulators and Rydberg atom arrays [43,44] and they can be implemented efficiently in a gate model. A natural requirement for the spectroscopy is that the decoherence time of the quantum device is not a limiting factor to determine the evolution time required in order to reach the target overlap. The applications of the adiabatic spectroscopy go beyond obtaining the spectral gap for a given Hamiltonian H(s). By aggregating spectral information from several adiabatic paths which cut through the phase diagram of a target Hamiltonian H T , rich properties of quantum many-body systems might be acquired. We note that the numerical results presented in (Fig. 1) support our argument that this technique is suitable to derive information about the spectral gap. However, the relation ∂T (s)/∂s ∼ 1/∆(s) 2 , which is obtained from the Landau-Zener model, is only a first order approximation.
Improvements of the quantitative validity of the adiabatic spectroscopy are left for future work and might build upon the rich literature on adiabatic perturbation theory [45]. We shall now discuss the VQAA from two perspectives. First, the perspective of VQAA as a quantum algorithm for optimal adiabatic paths, and second, VQAA as a variational quantum algorithm which requires only few measurements. Adiabatic quantum computation is known to prepare the ground state of target Hamiltonian H T for sufficiently long preparation time T . However, T scales as a function of the minimal spectral energy gap ∆(s). For general H T , the best rigorous bound on T has a inverse cubic gap dependence T = O((min s ∆(s)) −3 ) [11]. Due to the limited decoherence time of NISQ devices, large evolution times necessary for high-fidelity ground state preparation can be a difficulty. If it was possible to reduce T to fit into the coherent time frame of a quantum device, it would become possible to adiabatically prepare ground states, e.g. the solution of an optimization problem, that have remained unattainable before. The question of finding an optimal path for the adiabatic evolution has been in the focus of research efforts for several years already and for the unstructured search problem, the Grover-type speed-up has been recovered using an optimized adiabatic sweep [15]. However, the position and size of the spectral gap are, in general, a priori unknown, and obtaining the spectral properties of the adiabatic path can be as hard a problem as preparing the ground state. Therefore, it remains a challenge how an optimized adiabatic path can be obtained when no or only little spectral information is available. One recent work employed techniques from reinforcement learning to find an optimal adiabatic path [46]. In our work, we phrase the problem of finding an optimized adiabatic path as a problem to be solved through a variational quantum algorithm with step-wise adiabatic velocity profile. Turning now towards a discussion of the VQAA as a variational algorithm, we begin by noting that the ground state overlap can be a suitable cost function for quantumclassical feedback loops. Variational approaches such as the quantum approximate optimization algorithm (QAOA) [17] sparked intensive research interest in recent years.
The number of measurements necessary to estimate an objective function scales with O( −2 ), where is the maximum error that can be tolerated in the optimization process. The proposal of the VQAA aims to reduce the number of measurements by requiring relatively fewer parameters that need to be optimized and by considering a cost function with a low variance. Variational approaches for preparing a ground state generally make use of the energy as a cost function and that energy is estimated via the local observables that compose the Hamiltonian [47]. As the actual ground state will generally not be an eigenstate of the Hamiltonian local terms (e.g. if there is frustration), a low variance in the estimates cannot be guaranteed. Moreover, the orthogonal eigenstates in the low-energy sector typically yield similar energy values, hindering convergence to the ground state. If we were indeed able to directly measure in the eigenbasis of the Hamiltonian at a given point in the adiabatic path, we would seek to exploit the property of proximity to an eigenstate, which is inherent to adiabatic algorithms. The textbook approach for this problem would be a quantum phase estimation (QPE) algorithm, and direct implementation requires an ancilla overhead [22]. Recent semi-classical approaches are able to use a single ancilla only by utilising post-processing schemes [27,28].
For the VQAA, we suggest the overlap with the ground state as a figure of merit. In the case of the adiabatic algorithm, the optimal value of some other figures of merit, such as the energy, are not directly accessible. Therefore, we present two protocols to evaluate the closeness to an eigenstate using controlled unitary evolution. The entangled ancilla protocol offers the possibility to perform low-variance measurements by harnessing the power of hypothesis testing when being close to an eigenstate. We note that in the case of a small spectral gap, e.g. when a phase transition is crossed, the ground state overlaps are in general very small and special care is needed to extract useful information with the ancilla protocol. In the limit of very large depth, QAOA has the possibility to recover a trotterized adiabatic evolution. Therefore, a black box VQAA algorithm bears some similarities to QAOA. Several key differences are remarked though. First, analog to the evolution times in a (trotterized) adiabatic evolution, the unitaries in QAOA feature angles as parameters. However, for a quantum cost function H T , the optimized angles could be too large for the decoherence limit of the NISQ device which the algorithms is supposed to be implemented on. On the contrary, limiting the maximum angles could be rather problematic for the performance of QAOA. This issue does not arise in the black box VQAA as the total evolution time has been fixed. Second, even for large system sizes with 53 or 100 qubits, VQAA significantly improves over naive QAA for a very small number of parameters L only. In QAOA, it is generally expected that deep circuits are necessary to obtain good results for large systems. Finally, even for very small L, the performance of VQAA is lower bounded by the QAA with a linear time profile. This is true for QAOA only when the angles are initialized akin to a trotterized adiabatic evolution which requires large depth QAOA.

X. SUMMARY AND OUTLOOK
We seek to combine the best from two worlds by combining the strengths of the adiabatic and the variational approaches. We present a toolbox for VQAA building upon ancilla-based methods to evaluate the ground state fidelity at any point in the adiabatic path. Our approach only obtains information about the proximity to the ground state and is deliberately oblivious to the actual value of the energy throughout the adiabatic path. Due to the small parameter space and the ground state overlap as our cost function, the number of measurements necessary in the optimization of the adiabatic evolution is dramatically lower than for typical variational quantum algorithms such as QAOA. On the whole, our work suggests that a further exploration of NISQ algorithms based on variational adiabatic concepts is indicated. For instance, there seems to be room for sequential VQAA algorithms. In some instances, when going towards larger T , the black box VQAA does not improve in a strictly monotonous manner. By reusing information from shorter T , the optimized paths for larger T might be improved and obtained more rapidly. In a similar direction, it would be interesting to see how the information gathered through adiabatic spectroscopy will be used for the optimal adiabatic path in an experimental setting. If it was possible to use adiabatic spectroscopy to eliminate or drastically reduce the quantum-classical training phase of VQAA, this would surely be a large advancement in NISQ algorithms for ground state preparation. Also, it remains to be seen, what bounds can be established on how quickly and closely VQAA can find an adiabatic path which is optimal. Furthermore, there is increasing research interest in error mitigation techniques for NISQ devices [48,49], and a further exploration of how VQAA might benefit from these techniques appears promising. More generally, the protocols for estimating the ground state overlap presented in this work could make a wide range of new, exciting quantum algorithms for ground state preparation possible. This might include the opportunity for an algorithm to find optimized adiabatic paths using techniques from reinforcement learning or a combination of the protocols with techniques such as projected measurements and the quantum Zeno effect [50]. Besides, in the regime where the time evolution is not strictly adiabatic, high final ground states fidelities might be achieved by not starting in the ground state of the initial Hamiltonian but in an appropriate superposition in the low energy sector of the initial Hamiltonian. Note added. -After the completion and submission of this manuscript, some of the variational adiabatic methods presented here have been tested experimentally using large arrays of neutral Rydberg atoms to prepare the ground state of the maximum independent set problem on unit disk graphs with disorder [51]. We seek to better understand adiabatic spectroscopy as presented in the main text. Therefore, we would like to obtain a qualitative relation between the spectral gap ∆(s) along the adiabatic path and the evolution time T (s * ) required to prepare the ground state of H(s * ) with a given target fidelity. In order to do so, we turn towards the well-known Landau-Zener (LZ) model which describes a simple two-level system [52,53]. The model Hamiltonian is given as . (A1) With tan θ = g/λ(t), we write the eigenvectors as |a = sin(θ/2) − cos(θ/2) and |b = cos(θ/2) sin(θ/2) . (A2) The eigenenergies are given as E ± = ± λ(t) 2 + g 2 and the coupling is assumed to be a linear function in time λ(t) = δt. This implies that the minimum of the spectral gap (the avoided level-crossing) is found at t = 0, i.e. an adiabatic evolution parametrized by s from 0 to the position of a small gap s * is understood to be mapped onto the LZ evolution from very small initial t to t = 0. A perturbative approach yields the probability to find the system in the excited state at t f after initializing in the ground state at t i [54]. Assuming the beginning of the adiabatic evolution at t i = −∞, we are given 6 (A4) where we identified t with t f . Now, we set a target fidelity A 2 + := |α + (t)| 2 and assume that the total evolution time T scales as T ∼ 1/δ: As we are interested in the change of T , we compute the time-derivative of Ṫ where we used (Eqn. A5) in (Eqn. A6). Because of ∆(t) = 2E + (t) in the LZ model, we obtainṪ ∼ 1/∆(t) 2 as an approximation to the scaling ofṪ up to coefficients and possible corrections.
Appendix B: Calculations for the one-ancilla and entangled-ancillas protocol 1. Single-ancilla protocol The combined system after the unitary evolution can be written as Denoting the quantum state after the unitary evolution as |ψ evo = U |ψ , we note the following relationships regarding Pauli measurements of the ancilla For a fixed Hamiltonian H = j E j |φ j φ j | with the unitary U |φ j = e −iHτ |φ j = e −iEj τ |φ j , we write the state |ψ = j ψ j |φ j in the eigenbasis of H with |ψ normalized ( j |ψ j | 2 = 1). The total quantum system can be expressed as which we now use to calculate the density matrix of the ancilla qubit ρ ancilla = Tr non-ancilla (|ζ ζ|) = 1 2 j,m

Entangled-ancillas protocol
Our goal is to determine the purity of the ancilla. In general, there is the relation ρ pure ⇔ Tr[ρ 2 ] = λ 2 1 + λ 2 2 for density matrices, with λ 1 and λ 2 the eigenvalues of ρ.
We write the density matrix and its square as where matrix elements irrelevant for our protocol are denoted with a dot. With the Bell state we construct the Bell measurement operator The diagonal matrix elements of ρ 2 ancilla may then be attained by considering a composite system where the second system with controlled backward-time evolution (implementable by changing the sign of H). Then, the density matrix of the ancilla of the second system effectively corresponds to the transpose of the density matrix of the first ancilla ρ ancilla2 = ρ T ancilla1 . With (Eqn. B13), the composite system then gives If |ψ is in an eigenstate, the density matrix of the ancilla ρ ancilla is pure and the expectation of the |Φ − measurement is

Choosing suitable time values in the ancilla protocol
We have argued that a measure for the eigenstate closeness is given by (B20) using the one-ancilla protocol. For suitable τ , |ψ evo is an eigenstate of the fixed Hamiltonian H only if |α| = 1. The time τ needs to be chosen so that the complex summands of α with non-vanishing amplitude do not have approximately equal phases. In this unlikely case of matching phases, we would see constructive interference so that α = 1 could be true even if |ψ evo is not an eigenstate. Visualizing the summands of α on a complex plane (Fig. 18), this becomes rather intuitive. The choice for τ is related to the spectrum of H. It is reasonable to assume and confirmed in our simulations that a quantum state in our algorithm has the largest overlap with the ground state. As transitions to high energy excited states are extremely rare, for this argument we assume that the overlap with the first excited state was the only other non-vanishing overlap. Then, we simply have In the case of destructive interference which corresponds to with k ∈ Z and ∆ = E 1 − E 0 . An arbitrary τ would correspond to choosing an l ∈ Z in τ = πl/∆ at random. For l 1, the probability of choosing an odd value of l is approximately 1/2. Therefore, by testing several random values of τ ∈ O(∆ −1 ) it is possible to deduce information about the system whether it is in a mixed state or an eigenstate with high confidence.

Bound for the single-ancilla protocol
We extend the discussion about suitable values of τ . As argued above, an answer is generally dependent on For τ =1 the summand corresponding to the ground state and the summand corresponding to the first excited state lie quite closely together in phase which leads to problematic constructive interference (Fig. 4). However, as the pointers rotate with their respective eigenenergy for changing τ , they are well-separated for τ = 10. Then, we find destructive interference in the summation which is necessary for the explanatory power of the protocol. Note that in this instance, the summands corresponding to higher excited eigenstates are very small and not visible.
the spectral properties of H, which are not available a priori. We shall provide a bound for the single ancilla method without any knowledge about the spectrum of H or the populations of the different eigenstates. A byeffect of this generality is that the bound is not tight. This is because in a realistic setting of a (quasi-)adiabatic evolution, the eigenstate populations of higher excites are expected to be very small. We consider the expectation value [55] of α(τ ) In the limit of large K, we obtain where the spectral dependence has entirely averaged out.
Such an E 2 corresponds to what one would observe in the laboratory. We note that |ψ 0 | 2 was the ground state overlap of |ψ . As |ψ is normalized, we can write i.e. for all i which do not correspond to the ground state. Then, multiplying with |ψ i | 2 i>0 yields As the latter inequality holds for all i > 0, we obtain For the expectation value E 2 we can now write down the inequality Solving for |ψ 0 | 2 , we obtain So far we have not assumed anything about the populations. Letting the ground state population |ψ 0 | 2 be the largest of the eigenstate populations, we have as a lower bound for |ψ 0 | 2 for E 2 ≥ 1/2. Also, (Eqn. B34) implies that in this case the smaller ground state populations are therefore upper bounded as For E 2 < 1/2, no non-trivial bound for |ψ 0 | 2 (apart from |ψ 0 | 2 ∈ [0, 1]) can be given in the limit of very large K (with this approach).

Appendix C: Plots on the number of measurements and noise
We include plots showing the number of ground state overlaps for an instance of the black box algorithm for N = 100 and five chunks. Also, we present an estimate of the number of measurements needed using the Chernoff-Hoeffding inequality (Fig. 20). For a discussion of the inherent robustness of adiabatic algorithm, we show a the energy density of states for N = 12 for the ZZXZ model (Fig. 21a) and an analysis of the relative energy error due to a noisy gate at different positions in the adiabatic evolution (Fig. 21b).    Here, the adiabatic ramps have been smoothed using the reparametrization λ(s). The results for smoothened naive QAA are worse than in Fig. 2 leading to a larger separation between VQAA and QAA.

Theoretical aspects of an adiabatic evolution
Here we outline the main concepts of a theoretical treatment of adiabaticity. When considering transitions from the ground state to the excited states, the spectral gap between the ground state and the first excited state is of the greatest importance. Notably, the accumulation of relative phases plays a crucial role as well. We derive the differential equations governing the time evolution of a quantum state. Here, we follow the calculations given in [61], yet it is instructive to use parametrized time s = t/T instead of real time t. Then, the time-dependent and time-independent Schödinger equations with = 1 are given by respectively. The |k(s) , k ∈ [0, n], are the eigenvectors of H(s) and the eigenenergies E n (s) are for simplicity of exposition assumed to be non-degenerate and distinct. E 0 (s) is the ground state energy at a given s. We use the representation for a quantum state writing φ k = s 0 E k (τ )dτ , and insert it into the time-dependent Schrödinger equation (Eqn. S8) where the dot denotes the derivative with respect to parametrized time: The second term on the left hand side and the right hand side of (Eqn. S12) vanish as |ψ ( (S16) Then, the system of differential equations for the coefficients to describe the quantum state is given bẏ The i m(s)|k(s) are also referred to as the Berry connections [62].
withċ m (0) = δ nm , i.e. at the beginning of the evolution the whole population is in the nth eigenstate and will strictly remain in this eigenstate subspace if the evolution is perfectly adiabatic. For n = 0 we have the special case of the QAA: if a (perfect) adiabatic evolution starts in the ground state, the initial state will remain in the instantaneous ground state subspace until the target Hamiltonian is reached at s = 1.

Variational algorithms and parametrized circuits
The concept of a parametrized circuit is to write a trial wave function which approximates the desired quantum state. This trial state is prepared as a product of many unitary operators, of which each depends on a classical variational parameter acting upon an initial trivial state. The operator preparing this state can be written as details of the target Hamiltonian and appropriate parameter initialization can be necessary, albeit is ultimately not sufficient due to the dependence of the algorithm performance on the representative power of the chosen ansatz U (ϕ).

The quantum approximate optimization algorithm
The quantum approximate optimization algorithm (QAOA) [17] is a highly popular example of a parametrized ansatz and much research has been done at investigating the properties of this algorithm [38]. We outline QAOA here: For an initial Hamiltonian H 0 = N j=1 σ x j , the state |+ ⊗N is usually chosen to be the initial state. However, the ground state of H 0 is |− ⊗N and it would correspond to the initial state of an adiabatic quantum algorithm. Starting with an initial spin state we write down the QAOA ansatz where β = (β 1 , ..., β p ), γ = (γ 1 , ..., γ p ) with integer p ≥ 1. H 0 and H T are non-commuting Hamiltonians. H 0 = N j=1 σ x j is also called the mixing Hamiltonian and H T is the problem Hamiltonian. In the usual case, H T is diagonal in the computational basis and encodes combinatorial optimization problems, e.g. MAX-CUT, 3SAT, etc. The energy cost function is determined with a quantum device by repeated measurements in the computational basis. Using a classical computer, the energy is then variationally minimized towards a local minimum E p (γ * , β * ).
Appendix F: Simulating the adiabatic evolution with matrix product states

Tensor network techniques
In order to describe a general N -qubit quantum state with d spin degrees of freedom, there are d N coefficients c ii,i2,...,i N in |ψ = 0≤i1,i2,...,i N <d c ii,i2,...,i N |i 1 , i 2 , . . . , i N (S1) that we need to account for. The |i k d i=1 are forming an orthonormal basis of local Hilbert space (for 1 ≤ k < N ). However, and to our great benefit, ground states of gaped local one-dimensional Hamiltonians obey an area law regarding their entanglement entropy [64]. In fact, area law states occupy an exponentially small corner of the manybody Hilbert space only, which makes a more economical description of the quantum state possible. If we consider the coefficients c ii,i2,...,i N being a very large tensor which can be decomposed into a chain of tensors of order 3, we obtain an alternative approach. We describe the quantum state by a matrix product state (MPS) ansatz for open boundary conditions by The complex matrices B i k k are (D × D)-dimensional with the exception of the boundary tensors which are vectors. D is the dimension of the virtual bonds in the MPS representation, limiting the amount of entanglement that can be represented by the MPS ansatz. Describing the state as a matrix product state only requires N D 2 d parameters which is clearly preferable if a moderate virtual bond dimension D suffices. Therefore, MPS methods are an extremely successful ansatz to describe such quantum states obeying the area law such as aforementioned ground states [65]. In general, MPS are not well-suited to simulate the dynamics of a quantum system, due to the increase of entanglement in the system. However, in our algorithm which describes a (quasi-)adiabatic dynamic, the quantum states are very close to the ground state of the system. This making MPS a perfect candidate for classically simulating quantum adiabatic algorithms. For the simulation of time evolution we employ the time-evolving block decimation algorithm (TEBD) [66] which is the computationally cheapest method of several time-evolution methods for tensor networks and it is also well-suited for Hamiltonians with local interactions. This being said, other powerful time-evolution methods have been developed in recent years, such as MPO W I,II , TDVP and the global Krylov method [67]. An appealing property of the ZZXZ model is the suppressed the light cone spreading due to real time confinement [68]. This leads to a very slow increase in entanglement when time evolving the system, making this model a good candidate for simulations with MPS. The maximum bond dimension of the MPS could be kept at modest values while maintaining very good accuracy in the MPS representation.

Simulating noisy quantum circuits
In order to simulate noisy quantum circuits, two main approaches are feasible. A noisy quantum channel can be represented by working with density matrices instead of quantum states. This is because an MPS cannot represent a mixed state. By considering Matrix Product Density Operators instead of MPS, such noisy quantum channels can be efficiently dealt with. We have chosen the second approach where noisy circuits are simulated through statistically averaging over an ensemble of MPS with the noise tensors N being unitary matrices. Each of the MPS | ψ k will be different due to the probabilistic nature of noise. In the limit of a large MPS sample size n, both approaches will yield the same results. A discrete noise model is considered, where p is the noise strength and probability that a noise event occurs. Such a noise event can be either a bit flip (σ x ), a phase flip (σ z ) or a combination of both. The noise channel for discrete noise is The distinct noise tensors {N i } are contracted with respective tensors {B i } for i ∈ {1, . . . , N } before the first and after every N -site unitary (Fig. S3). As an example, in the discrete noise model with noise strength of p = 10 −4 , N = 100 qubits, a time evolution for total time T = 5 and 16 discrete steps per time unit, we expect on average 10 −4 · 100 · (1 + 5) · 16 = 0.96 ≈ 1 noise event in the entire quantum circuit. of others. One straightforward way to implement a variational quantum adiabatic routine is by making use of the ancilla-free method on an analogue quantum simulator. As quantum simulators do not require any Trotter overhead for the time evolution, the cost of the optimization procedure is directly given by the overhead in evolution time. Our numerical estimations in (Fig. 19) show that around 20-30 ground state fidelity estimations can suffice for a significant improvement of the performance of the adiabatic algorithm. In the case with the highest cost, the ancilla-free method is used only at the end of the adiabatic evolution. The backward-time evolution leads to an additional factor of 2 in the cost of this method.
There are different approaches to implementing the time evolution in the VQAA on a gate-model architecture. One straightforward estimate of the cost shall be given here by counting the gates required to implement the trotterized evolution as in the TEBD for the numerical simulation for a large one-dimensional quantum system. The time evolution on the qubit chain with N = 53 sites is decomposed into two-qubit unitaries. We chose M = 16, K = 2 and the Suzuki-Trotter order to be 2 (cf. Suppl.) and every unitary can be decomposed into 3 CNOTs and 7 arbitrary single-qubit gates [69]. Thus, per unit of time, around 120 CNOTs/N and 270 single-qubit gates / N are required. Next, we would like to count the gates required for the single-ancilla method and focus on the number of CNOTs. For the ancilla method, every two-qubit unitary from our consideration above is now controlled by the auxiliary qubit. Hence, we require 3 · 6 = 18 CNOTs for the three Toffoli gates [70] and at most 7 · 3 = 21 CNOTs to control the single qubit rotations, resulting in around 1530 CNOTs/N per unit of time as an upper bound. As the gate count for the ancilla method depends crucially on the value of τ , we make a remark on how the ancilla method would be used in practice. As explained in (App. B 3), the time τ for the controlled time evolution depends on the spectral gap of the respective Hamiltonian H(s). While one might be interested in the instantaneous ground state overlap along the adiabatic path, the ground state overlap at s = 1 is of most interest. Note that while the spectral gap ∆ along the path might be very small for hard problem instances, this is not necessarily the case at the end of the path. The ZZXZ model is a good example (Fig. 8) where ∆ is of the order of one at s = 1. A chain of trapped ions are an experimental system corresponding to our numerical analyses where no SWAP gates are required to realize the circuit. Superconducting-qubit devices, however, usually do not have an all-to-all connectivity and rather feature a two-dimensional architecture. We provide a worst-case-connectivity estimate for the number of SWAPs by doing the count for a chain of N qubits where only neighbouring qubits are connected. The ancilla qubit can be moved through the chain next to the two-site unitary to be applied with N − 2 SWAPs per brick-wall layer. Also, three extra SWAPs are necessary per three-qubit controlled unitary so that the single qubit rotation to be controlled is next to the ancilla qubit. All-in-all, the overhead due to SWAP gates is at 580 additional CNOTs/N for swapping operations on this chain topology per unit of time.
We stress again that these estimates are concerned with counting the CNOTs that would have been required to implement the gates used in the numerical simulation on a quantum device. The actual gate count will be significantly lower if the circuit is optimized with regard to the experimental realization instead.

Appendix H: Bayesian inference in hypothesis testing
We give a brief summary on how to use Bayesian inference in order to do highly efficient hypothesis testing using our entangled-ancillas protocol (Sec. IV B) to decide whether the ground state overlap is larger than a given threshold value. In an experimental setup, we would like to only obtain as few samples as needed. Therefore, after each measurement, we would like make use of a stopping criterion to decide whether we need to continue sampling or are able to terminate. Unfortunately, stopping rules are quite cumbersome to deal with in a frequentist approach. And fortunately, Bayesian inference turns out to be a most convenient tool.
The measurement outcomes of a quantum expectation value such as the ground state overlap with our protocol are assumed to behave as an independent and identically distributed (i.i.d.) random Bernoulli variable. The total measurement data is X = {x 1 , ..., x N } with a Beta prior θ = Beta(a i , b i ). The prior prevents overfitting for few data points. The Beta distribution is normalized and given as The mean and the variance of the Beta distribution are var[X] = ab (a + b) 2 (a + b + 1) .
The Beta distribution is the conjugated prior for the Bernoulli distribution, implying that the posterior after one measurement is again a Beta distribution with new parameters [71]: When sampling in the experiment, our posterior becomes the prior for the next measurement until our stopping criterion is fulfilled. We merely need to keep track of the 0's and 1's or spin-ups and spin-downs that are measured.
p(θ|X) = Beta(a N , b N ), (S8) If we are expecting mostly zeros when close to the ground state, we might initialize our prior with a = 10 and b = 2. If the algorithm asks for a ground state overlap greater than H 0 (Zero hypothesis: | ψ|ψ evo | > H 0 ), samples are gathered until a stopping criterion is met and H 0 is either accepted or rejected. A decision upon the hypothesis is hard when H 0 is very close to the mean as many samples will be necessary to distinguish them. Setting up an interval [H 0 − , H 0 + ] where a decision cannot be made overcomes this problem. Then, H 0 can be updated to a different value and the sampling is restarted. The following outcomes of hypothesis testing are possible: Identifying a 0-measurement (perfect ground state overlap) with increasing a by 1 and a 1-measurement with increasing b by 1, the a N and b N contain the sampling data after N samples. When close to a ground state, we expect mainly a's to be counted and a distribution strongly biased towards p = 1. A suitable stopping criterion might be an alpha error probability of under 5%. A sample case is shown in (Fig. S4).