How to experimentally evaluate the adiabatic condition for quantum annealing

We propose an experimental method for evaluating the adiabatic condition during quantum annealing (QA), which will be essential for solving practical problems. The adiabatic condition consists of the transition matrix element and the energy gap, and our method simultaneously provides information about these components without diagonalizing the Hamiltonian. The key idea is to measure the power spectrum of a time domain signal by adding an oscillating field during QA, and we can estimate the values of the transition matrix element and energy gap from the measurement output. Our results provides a powerful experimental basis for analyzing the performance of QA.


I. INTRODUCTION
The adiabatic theorem is a crucial result in quantum mechanics, first introduced by Ehrenfest in 1916 [1].Later, Born and Fock proved a more modern version of the theorem in 1928.The theorem states that if an initial state is prepared in the ground state of the Hamiltonian, it will remain in the ground state as long as the change in the Hamiltonian is slow enough.Since Born and Fock's proof in 1928, there have been numerous studies that have improved and expanded the theorem, including more rigorous formulations [2] and extensions to open systems [3][4][5].
An essential application of this theorem is quantum annealing (QA).This was originally proposed by Apolloni et al. in 1989 [6].The original proposal aimed to improve the simulated annealing utilizing the quantum effects of tunneling.However, an alternative approach was subsequently presented [7,8], where the Hamiltonian changes over time.In this approach, a ground state of the transverse-field Hamiltonian is prepared, and the Hamiltonian is gradually changed to the target problem Hamiltonian.The adiabatic theorem guarantees that if the alteration of the Hamiltonian is gradual enough, the final state will be the ground state of the problem Hamiltonian.
One of the problems in QA is that there is no known efficient method for checking whether the adiabaticity is satisfied or not.In principle, if we can diagonalize the Hamiltonian, we can use an approximate version of the adiabatic conditions are given as follows [22][23][24][25]: for all s and m, where T ann denotes the annealing time, s = t/T ann denotes the time normalized by T ann , t denotes the time, |m(s)⟩ (|0(s)⟩) denotes the m-th excited (ground) state, Ḣ(s) denotes the s derivative of the instantaneous Hamiltonian at a time s and E m (s) (E 0 (s)) denote the eigenenergies of the m-th excited (ground) state (see Appendix A).Throughout this paper, we consider a dimensionless time s normalized by T ann .These conditions are obtained by an argument that considers only the first order perturbation expansion and neglects higher order terms [26], and so are not mathematically rigorous.In particular, conditions (1) are not known to be sufficient for adiabaticity.However, when the interest is in the qualitative properties of the computation time, these conditions are widely used, and so we adopt them as the adiabatic conditions in our paper.
In the case of applying QA to practical problems, it is unworkable to diagonalize the Hamiltonian with using a classical computer.Consequently, we cannot directly apply the adidabatic conditions (1) to check whether the dynamics is adiabatic or not.Experimental methods have been proposed to measure the energy gap [16,27,28], which corresponds to the denominator in Eq. (1).However, to our knowledge, no studies have been conducted to measure the numerator of the adiabatic condition (1), i.e., the size of the transition matrix element of the time derivative of the Hamiltonian.
In this paper, we propose a method for simultaneously measuring the numerator and denominator of Eq. ( 1).This method involves utilizing an oscillating field during quantum annealing to induce a Rabi oscillation between the ground and excited states.By performing Fourier transformation on a time domain signal, we obtain a power spectrum and extract relevant information from the data.These steps enable us to evaluate the values of the numerator and denominator of the adiabatic condition (1).
The remainder of this paper is organized as follows.In Sec.II, we review QA.In Sec.III, we introduce our method for simultaneously measuring the values of the transition matrix element and the energy gap, based on an analytical calculation using some approximations.In Sec.IV, we describe numerical simulations (with noise) performed to quantify the performance of our method in realistic cases.Finally, in Sec.V, we summarize our results and discuss possible directions for future work.

II. REVIEW OF QA
To review the conventional QA, we consider the following Hamiltonian: where H D is a driver Hamiltonian, H P is a problem Hamiltonian, and f (s) is a schedule function satisfying the condition Here and in the following we make the choice Due to the condition (3), the Hamiltonian at s = 0 is the driver Hamiltonian H D and the Hamiltonian at s = 1 is the problem Hamiltonian.After obtaining a ground state of the driver Hamiltonian, we let the state evolve by the annealing Hamiltonian from s = 0 to s = 1.According to the adiabatic theorem, if the annealing time T ann is sufficiently large, the state after QA becomes a ground state of the problem Hamiltonian.

III. OUR METHOD FOR EVALUATING THE ADIABATIC CONDITION
We will now present a technique to experimentally determine the numerator and denominator of the left-hand side of Eq. ( 1) for a given time s 1 , using the Hamiltonian defined in Eq. (2).In this scenario, the Hamiltonian in Eq. ( 1) is the Hamiltonian for quantum annealing H conv defined by Eq. (2).We introduce the Hamiltonian H(s), which comprises the driver Hamiltonian H D , the problem Hamiltonian H P , and an external driving Hamiltonian H ext (s) with strength λ(s) and frequency ω as follows.
H ext (s) = λ(s) Ḣconv (s 1 ) cos (ωT ann (s − s 1 )) (7) Here, A(s) is the schedule function that modulates the weight of H D and H P in H QA (s).We plot A(s) and λ(s) as functions of time in Fig. 1, where we note that A(s) satisfies A(s) = f (s) when 0 ≤ s ≤ s 1 .
Our experimental protocol proceeds as follows.Firstly, we prepare the ground state of the driver Hamiltonian |0(s = 0)⟩.Secondly, we slowly vary the Hamiltonian H QA (s) from s = 0 to s = s 1 by setting λ(s) = 0, allowing the system to evolve under this Hamiltonian adiabatically.Thirdly, at s = s 1 , we introduce a driving term by setting λ(s) = λ and fixing A(s) = f (s 1 ), and we let the system evolve for s 1 < s ≤ s 1 + τ /T ann .Fourthly, we terminate the driving at s = s 1 +τ /T ann by setting λ(s) = 0, and gradually vary the Hamiltonian from H QA (s 1 ) to H D for s 1 + τ /T ann < s ≤ 2s 1 + τ /T ann , allowing the system to evolve adiabatically.Finally, we measure the probability of the system occupying the m-th excited state |m(s = 0)⟩ of the driver Hamiltonian using projective measurements, which we denote as p 0,m (ω, s 1 , τ ).We repeat these steps multiple times, varying ω, s 1 and τ .We emphasize the importance of adiabaticity during the second and fourth steps, while it is not necessary for the third step.
Let us explain how to realize H ext (s) in the third step of the actual experiment.We have The driver Hamiltonian and problem Hamiltonian can be decomposed using the Pauli operators as follows: where O i (O ′ j ) denote the Pauli matrices and h i (h ′ j ) denotes a time-independent coefficient.Hence, we obtain Thus, if we can temporarily change the coefficient of the Pauli matrices to a cosine function, it is possible to realize the Hamiltonian H ext (s).As the problem Hamiltonian usually contains two-body interaction terms, we must change the interaction coupling strength.Such a technique has been developed for superconducting circuits [29].
Here, we describe the dynamics of the system in the third step of our scheme, which is crucial for measuring the adiabatic condition.We begin by describing a simplified scenario in which the dynamics is adiabatic in the second and fourth steps, and we will consider more general cases later.For simplicity, we omit the expression of "(s 1 )" to mention H QA (s 1 ) or Ḣconv (s 1 ) in the remainder of this section.In our proposal, the measurements are performed while sweeping the time period τ ; hence, we treat τ as a variable in the remainder of this section unless mentioned otherwise.
Let us diagonalize H QA as follows: where E i ≤ E j is satisfied for i < j.By moving to a rotating frame, we can express the state of the system as follows: and the Hamiltonian in the rotating frame is expressed as Note that we set ℏ = 1 throughout this paper.Here, we assume that the transition frequency between the ground state and the m-th excited state is close to the frequency of the external driving field.Then, we set r as the ratio between |E m − E 0 | and ω as follows: where E 0 denotes the energy of the ground state.The second term in Eq. ( 14) becomes ⟨i| Ḣconv |j⟩ e ir(Ei−Ej )τ cos ωτ |i⟩ ⟨j| .
Here, we adopt the rotating wave approximation (RWA) [30].The coefficient |i⟩ ⟨j| in Eq. ( 16) includes an oscillatory component: If r|E i − E j | = ω is satisfied, one of the terms in Eq. ( 17) becomes time-independent while the other term has a high-frequency oscillation.Owing to the condition of Eq. ( 15), we have at least two time-independent terms, (i, j) = (m, 0) and (0, m), which remain after RWA.We assume a condition Then, all terms except (i, j) = (0, m) and (i, j) = (m, 0) are dropped, and the Hamiltonian ( 16) can be simplified as . Therefore, the effective Hamiltonian Eq. ( 14) can be expressed as These calculations indicate that if the initial state is prepared in a subspace spanned by the ground state and m-th excited state, the system's dynamics will be confined to this subspace.Notably, projecting out the states except |m(s = s 1 )⟩ and |0(s = s 1 )⟩ results in an effective Hamiltonian with the same structure as the single-qubit Hamiltonian that induces Rabi oscillations.A known analytical formula that characterizes tha Rabi oscillation without decoherence involves two parameters: detuning and Rabi frequancy, and details of the behavior of Rabi oscillations in a single-qubit system are presented in Appendix B. By using this analytical formula, we can fit the data obtained from our method and acquire information about the transition matrix element | ⟨m| Ḣ|0⟩ | and the energy gap E m − E 0 .
To observe the oscillation experimentally, we need to construct a projective measurement of |m⟩ ⟨m| in the rotating frame.In our idea, the fourth and fifth steps enable us to construct a projective measurement |m⟩ ⟨m| in the laboratory frame effectively, provided the dynamics in the fourth step is adiabatic.If the state |ψ(τ )⟩ is an eigenstate of the Hamiltonian H QA , the change in the frame only results in a global phase.Therefore, as long as the second step and fourth step are adiabatically performed, p 0,m (ω, s 1 , τ ) is approximately described as follows: where α(ω) is a function of ω given by and Ω ana (ω) is a hyperbolic curve on the ω − Ω plane that is represented by We obtain the right-hand side of Eq. ( 19), which is independent of s 1 , under an assumption that the adiabatic condition is satisfied at the second and fourth steps.However, if there are non-adiabatic transitions, the probability p 0,m has a dependence on s 1 .
In the aforementioned discussion, a ground state of the driver Hamiltonian is assumed to be prepared in the first step, and we perform a projective measurement into the m-th excited state in the fifth step.Meanwhile, if we prepare the k-th excited state in the first step and perform a projective measurement into the l-th excited state in the fifth step, we can obtain the hyperbolic curve as Ω The details of these derivations are presented in Appendix B 1.
The adiabatic condition described in Eq. ( 1) is valid only when we can consider that the effect of the nonadiabatic transitions is weak.Therefore, throughout our paper (except in the Appendix), we assume that the effect of non-adiabatic transitions is negligible.We will discuss how the non-adiabatic transitions affect the spectroscopic measurements in our methods later.
Let us explain how to specify the values of |E m − E 0 | and | ⟨m| Ḣconv |0⟩ | by using our method.We repeat these by sweeping ω, and we can find an optimal value of ω = |E m − E 0 | to minimize the frequency of the Rabi oscillation; this corresponds to the energy gap ∆.Furthermore, the Rabi frequency with the optimal Ω observed in our method corresponds to the numerator in Eq. ( 1).Thus, our estimated transition matrix element | ⟨m| Ḣ|0⟩ | est and our estimated energy gap ∆ est are given by respectively.Here, Ω exp (ω) is the angular frequency of the Rabi oscillation obtained experimentally, which is analytically considered to be expressed by Eq. ( 22).
In actual experiments, owing to some imperfections, p 0,m (ω, t 1 , τ ) cannot be fully explained by the analytical formula Eq. ( 19), which was derived under ideal conditions (see Fig. 11 in the Appendix).To find the relevant frequency of Ω exp (ω) in the dynamics, we perform a Fourier transformation and obtain a power spectrum that is defined by If p 0,m (ω, s 1 , τ ) is expressed as Eq. ( 19), the power spectrum is given by Therefore, in the actual experiment, we define the peak with a positive frequency in the spectrum as Ω exp (ω), and we expect to satisfy Ω exp (ω) ≃ Ω ana (ω) in the power spectrum; this allows us to use the formulas of Eqs. ( 23) and ( 24).Thus, we can estimate the values of the transition matrix element | ⟨m| Ḣconv |0⟩ | and the energy gap ∆ using our method.

IV. NUMERICAL ANALYSIS
We perform numerical simulations to evaluate the effectiveness of our method.In the previous section, we derived an analytical formula under the following assumptions: I.The time evolution is adiabatic in both step 2 and step 4.
II.The rotating wave approximation holds.
III.The time evolution in step 3 only involves the ground state and the m-th excited state.
IV.There is no decoherence.
However, these assumptions may not be met in actual experiments and we perform numerical simulations to examine the validity of our method under different conditions, as summarized in Table I.Condition I is only satisfied when the process in steps 2 and 4 is completely adiabatic.In cases A and D from Table I, we use diagonalization to prepare the ground state of H QA (s 1 ).In the remaining cases, we solve the timedependent Schroedinger equation with specific annealing times to prepare the ground state of H QA (s 1 ).
Condition III is naturally satisfied for a single-qubit system, while it is violated for a system with two or more qubits.Thus, in cases A, B, and C, condition III is satisfied, whereas in cases D, E, and F, condition III is violated.Condition IV is satisfied if we solve a time-dependent Schroedinger equation of the system, as with cases A, B, D, and E. Meanwhile, we consider the effect of decoherence by solving the master equation in cases C and F.

A. Settings and methods for all cases
Here, we introduce some conditions that are common throughout our numerical analysis.

Schedule function
For the schedule function A(s) in Eq. ( 5), we use where T ann is the annealing time.In actual experiments, this value is typically around 10 to 100 µs, and the typical energy scale of the Hamiltonian is of the order of GHz [31].We take the schedule function ( 27) as A(s) = 1 − s up to s 1 , and we evaluate the adiabatic condition at time s 1 according to our method.Hence, for our simulation, ḢQA (s 1 ) is given by for any t 1 .

Strength λ
The Rabi frequency can be controlled by changing the strength λ.If the decoherence is negligible, we set λ to be as small as possible, because RWA is valid only when the Rabi frequency is much smaller than the energy gap ∆.Meanwhile, when there is decoherence, the choice of λ is not straightforward.As we decrease λ, the decoherence becomes more relevant and RWA becomes more valid.Therefore, the following condition should be satisfied: where T c is the coherence time.In our simulation, we set λ = 0.05.

Time evolution and measurement process
In real experiments, we need to consider decoherence.To address this, we use the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation [32], for cases C and F, where L n is the Lindblad operators.
In the fifth step, we assume that an ideal projective measurement into the state |m(t = 0)⟩ can be performed.It is worth mentioning that, in the actual experiment, this projective measurement corresponds to σ x on all qubits.By using a post processing with a classical computer, we can obtain the projection probability for not only m = 1 but also all m.However, the non-adiabatic transitions between the ground state and the first excited state is considered as the most relevant part.Actually, as long as | ⟨1| Ḣ|0⟩ | is similar to or larger than | ⟨m| Ḣ|0⟩ | for m ≥ 2, the non-adiabatic transitions between the ground state and the first excited state is more relevant than the others, Thus, for the numerical simulations, we consider a case of m = 1 in this paper.(See Fig. 2)

Construction of Ωexp(ω)
In our method, we calculate the probability of a projection into the first excited state p 0,1 (τ ), and we use the power spectrum P (Ω) to determine the function Ω exp (ω) as explained in the previous section.In this case, we expect to observe a peak at Ω = Ω ana (ω) in the power spectrum.To determine the function Ω exp (ω), we fix ω and maximize the height of the power spectrum by sweeping Ω so that we can determine the position of the resonance peak as follows: Finally, by sweeping ω, we can obtain the function Ω exp (ω).When we sweep Ω, it is crucial to choose an appropriate range.First, we explore the frequency range Ω > 0. As indicated by Eq. ( 26), three peaks emerge.However, to evaluate the adiabatic condition, our focus lies solely on the positive-frequency peak, because the negative frequency peak contains the same information as its positive counterpart, while the zero frequency peak lacks relevant information.Also, we should consider only the frequency range of Ω ≪ ω because we use RWA to derive the analytical formula of Eq. ( 19), which is valid only for Ω ≪ ω.
Even if we restrict the frequency range, we may not find a correct peak for several reasons.We discuss the case in which such a problem occurs, and we present a possible solution to overcome such a problem at least for some cases.

B. Single-qubit cases (A, B, and C)
We examine the single-qubit cases (A, B, and C).For these cases, the driver Hamiltonian H D and the problem Hamiltonian H P are given by respectively.In our simulation, we fixed ω 1 = 1 GHz and g = 0.4 GHz.

Case A
We set the parameters T ann and s 1 as follows: T ann = 10, 30, 100, 300, 1000 ns, As shown in Fig. 3, our estimated values (dots in the figure) are in good agreement with the theoretically expected values (lines in the figure).Indeed, the relative error in the estimation of the transition matrix element | ⟨1| Ḣ|0⟩ | (the energy gap E 1 − E 0 ) is at most 0.99 % (0.071 %).These errors are small compared to the resolution owing to the discretization performed while processing the data.The estimation error of the transition matrix element (energy gap) is 0.9 (0.1) times smaller than the resolution.As shown in Fig. 3, we confirm that the adiabatic condition (1) is reasonably satisfied.

Case B
Next, the effect of non-adiabatic transitions in steps 2 and 4 is studied for case B. Similar to case A, we can accurately measure both the transition matrix element | ⟨1| Ḣ|0⟩ | and the energy gap (E 1 − E 0 ) for case B, and (see Fig. 4) the relative error of the transition matrix element (energy gap) is at most 2.2 % (0.7 %) and 0.77 (0.93) relative to the resolution.
For the single-qubit case, our scheme is robust against the non-adiabatic transitions.Actually, we consider cases with T ann = 1, 2, 4, and 8 ns (see Fig. 5), and these results show that a shorter annealing time does not impair the performance of our methods.
We show that, as long as RWA is valid, the power spectrum contains a peak corresponding to a frequency of Ω(ω) (see Appendix B 2).Thus, we can accurately estimate the transition matrix element and energy gap using Eqs.( 23) and ( 24) for the single-qubit case without decoherence.
FIG. 4. Top (bottom): estimated value of the transition matrix element (energy gap) in case B (single qubit, incomplete adiabaticity, and no decoherence).For the solid lines and dots, we use the same notation as that in Fig. 3.

Case C
In case C, to consider decoherence, we employ the GKSL master equation, and we select the Lindblad operator as where κ denotes the decay rate.We fix κ = 2.5 × 10 −3 ns −1 , which is a typical value for a superconducting flux qubit [33].
The results are shown in Fig. 6.The relative error of the transition matrix element | ⟨1| Ḣ|0⟩ | (the energy gap ∆) is at most 2.1 % (0.05 %), which is 0.77 (0.023) times smaller than the resolution.
These errors are as small as those in cases A and B, indicating the robustness of our method against decoherence.This resilience stems from the fact that decoherence primarily impacts the width rather than the position of the peaks in the power spectrum.Consequently, accurate estimation of the transition matrix element and energy gap remains achievable even in the presence of weak decoherence.

C. Two-qubit cases (D, E, and F)
In the two-qubit cases, the problem and driver Hamiltonians are given by respectively.Here, we set ω 1 = 1.0 GHz, ω 2 = 1.1 GHz, g 1 = 0.5 GHz, g 2 = 0.3 GHz, and g 3 = 0.For these cases, we select the parameters T ann and s 1 as follows.

Case D
In this case, we can accurately measure the transition matrix element | ⟨1| Ḣ|0⟩ | and the energy gap (E 1 − E 0 ) as shown in Fig. 7.The relative error of the transition matrix element (energy gap) is at most 3.5 % (0.04 %), which is 0.55 (0.99) times smaller than the resolution.FIG. 6. Top (bottom): estimated value of the transition matrix element (energy gap) in case C (single qubit, incomplete adiabaticity, and decoherence).For the solid lines and dots, we use the same notation as that in Fig. 3.
Despite not satisfying condition III for considering two qubits in this case, the dynamics can be effectively confined within a two-level system, ensuring the accuracy of our method, especially when the Rabi frequency is low.

Case E
In case E, the relative error of the transition matrix element (energy gap) is at most 3.5 % (1.2 %), which is 0.55 (0.99) times smaller than the resolution, as shown in Fig. 8.
In the case of weak non-adiabatic transitions, it is possible to estimate both the transition matrix element and the energy gap with high accuracy even for the two-qubit case.Meanwhile, as described in detail in Appendix C, in the case of strong non-adiabatic transitions, the power spectrum contains peaks other than the one that we want to use in our estimation.We discuss a possible solution for this problem in Appendix C.

Case F
In case F, we select the Lindblad operator as follows.
FIG. 7. Top (bottom): estimated value of the transition matrix element (energy gap) in case D (two qubits, complete adiabaticity, and no decoherence).Even when two qubits are used, we can estimate both the transition matrix element and the energy gap with high accuracy.For the solid lines and dots, we use the same notation as that in Fig. 3.
Here, κ denotes the decay rate.For the numerical simulations, we chose κ = 2.5 × 10 −3 ns −1 .We plot the estimated transition matrix element and energy gap against s 1 , and we demonstrate that our method is accurate except for two points, s 1 = 0.3 and s 1 = 0.9 for T ann = 100 ns, as shown in Fig. 9.In the former case, as shown in the power spectrum (see Fig. 10 (a)), where ω is smaller than 0.575 or larger than 0.65, a low-frequency (Ω < 0.02) peak exists, and the height of this peak is greater than that of the target peak at the same ω.
As shown in Eq. ( 26), strictly speaking, a peak around Ω ≃ 0 should exist in the spectrum, and this peak has a finite width owing to decoherence so that we can observe this in case F. Therefore, if we naively adopt our method described in Eq. ( 31), we generate an inappropriate Ω exp (ω) and obtain incorrect estimated values of the transition matrix element and energy gap.
To identify the target peak in the presence of decoherence and non-adiabatic conditions, we employ a modified approach outlined as follows.Initially, we assess the value of Ω not only for the highest peak but also for the second-and third-highest peaks at each ω.These pairs of values (ω, Ω) then constitute candidates for the data in the estimated function Ω exp (ω).Subsequently, FIG. 8. Top (bottom): estimated value of the transition matrix element (energy gap) in case E (two qubits, incomplete adiabaticity, and no decoherence).For the solid lines and dots, we use the same notation as that in Fig. 3.
we attempt to fit the data using the analytical formula presented in Eq. ( 22).In the third step, we eliminate data that cannot be adequately fitted by the analytical formula.Finally, we designate data successfully fitting the analytical formula as the target peaks.
In the former case (s 1 = 0.3), after using this modified method, the relative error of the transition matrix element (energy gap) is 1.0% (0.2%), and the ratio to the resolution is 0.07 (0.26).Thus, our modified method is effective for this case, as shown in Fig. 10 (a).
However, in the latter case (s 1 = 0.9), the decoherence is so strong that the target peak nearly disappears, and we cannot identify the target peak anymore, as shown in Fig. 10 (b).

V. CONCLUSIONS AND DISCUSSION
We have proposed an experimental method to assess adiabaticity in QA by evaluating adiabatic conditions.Our approach uses an oscillating field to induce Rabi oscillations, providing insights into the energy gap and the transition matrix element of the time derivative of the Hamiltonian.To validate our method, we performed numerical simulations, considering non-adiabatic transitions and decoherence effects.The results confirm the robustness of our method against these experimentally FIG. 9. Top (bottom): estimated value of the transition matrix element (energy gap) in case F (two qubits, incomplete adiabaticity, and decoherence).In this case, we have significant estimation errors for a few points.For the solid lines and dots, we use the same notation as that in Fig. 3. inevitable problems.
Our method is valuable for determining an optimal annealing schedule and Hamiltonian form to improve the performance of QA.When a phase transition takes place, the performance of QA is degraded.Some methods have been proposed to address this issue in specific cases [34][35][36].To apply these methods, we need to change the annealing schedule and form of the Hamiltonian.However, a potential problem is that we cannot easily optimize the annealing scheduling and form of the Hamiltonian for general problems if we do not know whether the adiabatic conditions are satisfied.Meanwhile, by using our methods to evaluate the adiabaticity of the dynamics, we can select an appropriate annealing schedule and form of the Hamiltonian when we try to solve practical optimization problems using QA.Also, we discuss possible experimental implementations of our proposal.In the currently available Dwave quantum annealing machine, it is not possible to perform our proposal because we cannot perform microwave pulses to induce the Rabi oscillation in the annealer.However, the recently proposed method, called the spin-lock quantum annealing [12,37], is compatible with the requirement of our method.In the spin-lock quantum annealing, we use superconducting qubits for a gate type quantum computer, and there is an exper- imental demonstration of the spin-lock with superconducting qubits by using microwave pulses [38].Let us also examine the validity of the parameters used in our numerical simulations.During the spin lock, the system is in the rotating frame and so we can set ω 1 = 1.0 GHz and ω 2 = 1.1 GHz as the detuning between the microwave frequency and qubit resonance.The reported coherence time of the superconducting qubit [39] is much longer than 1/κ = 400 ns which is used in our simulations.Also, it is possible to realize a coupling strength of g 1 = 0.5 GHz and g 2 = 0.3 GHz by using inductive coupling between the superconducting qubits [40].We set the Rabi frequency λ to be around tens of MHz, which is also available in the current experiment [41].
It is important to note that our method remains applicable even if the number of qubits increases.Even in such scenarios, there must be finite energy gaps between the eigenenergies.Our approach is designed to effectively operate in such situations (as long as the energy gap is finite) where long-lived qubits are available for quantum annealing.
In QA, the energy gap could approach to zero as we increase the size of the system.In such a case, we will find a highly degenerate spectrum.In our method, we sweep s 1 and investigate the adiabatic condition for several values of s 1 .Since we sweep the value of s 1 from 0 to another value in our method, we can study the adiabatic condition during QA before the energy gap closes.In this case, we will recognize that the energy gap becomes smaller as we increase s 1 , and this lets us know the existence of the energy gap closing.We could adopt several strategies to enlarge the energy gap in this case.For example, twisted field [42,43], counterdiabatic term [44,45], nonstoquastic Hamiltonian [46][47][48], inhomogeneous driving magnetic field [49] are helpful for such a purpose.It is considered as difficult to find such additional fields to resolve the problem of the energy gap closing.Importantly, our method is useful to find such additional fields for the following reasons.When we introduce the additional fields, the adiabatic conditions will be changed, and we can experimentally measure the adiabatic conditions from the spectrum.So, if introducing the additional fields improves the adiabatic conditions, we can detect it.In principle, by measuring the adiabatic conditions, we can variationally change the strength and the direction of the additional fields to optimize the adiabatic conditions.The detailed research about this is left for future work.In the revised manuscript, we explained this point.
Finally, we comment on the adiabatic condition itself.In general, it has not been proved that the condition ( 1) is sufficient to achieve the adiabaticity.In addition, more sophisticated criteria have been proposed [26], and it was shown that higher order derivative of the annealing Hamiltonian could affect the adiabaticity.In our numerical examples, we show that the condition (1) actually provides an upper bound of the population of the excited state due to the non-adibatic transitions in Appendix E. However, if we need to know the information of the higher order derivative of the Hamiltonian, we can use a modified version of our method.For example, if we are interested in the value of | ⟨m| Ḧ|0⟩ |, we can replace Ḣ(s 1 ) in Eq. ( 7) with Ḧ(s 1 ).We leave a detailed study of this for future work.Secondly, although the conventional adiabatic condition in Eq. ( 1) is derived from the unitary evolution, it is possible to generalize the adiabatic theorem to open quantum systems [4].The aim of our method is to know the value of Eq. ( 1), which is different from the adiabatic condition in the open quantum systems.We also leave the extension of our method to the adiabatic condition in open quantum systems for future work.
Based on the discussion in the main text, we consider a Rabi oscillation between the states |k⟩ and |l⟩ (for example, the Hamiltonian in Eq. ( 18) corresponds to a case with k = 0 and l = m).We have where σ z = |l⟩ ⟨l| − |k⟩ ⟨k|, ∆ = E l − E k , and λ = λ ⟨l| ḢQA |k⟩.This coincides with the conventional Hamiltonian to induce the Rabi oscillation, where λ denotes the Rabi frequency and (1 − r)∆ denotes the detuning.For this Hamiltonian, the operation e iθσz is adopted for rotation about the z-axis by an appropriate angle θ.Then, the Hamiltonian (B1) becomes where σ x = |k⟩ ⟨l| + |l⟩ ⟨k|.We can rewrite Eq. (B2) with a unitary operator U diag as where λ′ satisfies By setting r = ω/∆, we obtain λ′ = Ω (k,l) ana (ω) in Eq. ( 22).FIG.11.Plot of the Rabi oscillation.We use a single-qubit Hamiltonian (32), with Tann = 30 and t1 = 0.3 Tann.The oscillation period is the longest when the angular frequency ω of the external field coincides with the energy gap ∆ = 0.74, indicating that the oscillation period is shorter for small angular frequencies (ω = 0.715) and large angular frequencies (ω = 0.785).
To observe the Rabi oscillations, we calculate the amplitude as follows: using the relations Finally, we obtain Owing to the form of Eq. (B4), the minimum value of λ′ is given by s = 1; then, the angular frequency of the Rabi oscillation λ′ is | λ| = λ| ⟨l| Ḣ|k⟩ |.In this case, U diag = I and the amplitude is maximized.Thus, the Rabi oscillation is a type of resonance phenomenon whose resonant angular frequency is the energy gap ∆.When the dynamics in the second and fourth steps is adiabatic, both the initial state and the final state are energy eigenstates.Meanwhile, Fig. 11 shows that the Rabi oscillations were obtained numerically using the Hamiltonian induced by Eq. ( 32) directly, without the approximation described in this section.It is generally a simple sinusoidal curve; however, there are slight fine oscillations when the details are observed.

Dynamics with non-adiabatic conditions
We explain the dynamics of our system when the non-adiabatic transitions occur in the second and fourth steps.Here, for simplicity, we assume that RWA is valid.Owing to the non-adiabatic transitions, in general, our prepared state at time t 1 is not an energy eigenstate; it is a superposition state of the energy eigenstates.Similarly, the non-adiabatic transition also occurs from t 1 + τ to 2t 1 + τ .We will explain that these non-adiabatic transitions cause high-frequency oscillations in our scheme.
In step 3, we prepare a state |ψ(0)⟩ and let this state evolve by the Hamiltonian for time t.The state we obtain at time t denotes |ψ(t)⟩.Combining the unitary evolution in step 4 with the projective measurement performed in step 5, this process can be considered as a projective measurement |ϕ(t)⟩⟨ϕ(t)| on the state |ψ(t)⟩.
An overlap between the states is given by ⟨ϕ(t)|ψ(t)⟩ = ⟨ φ(t)| ψ(t)⟩, where we define (B9) Using RWA, the final transition amplitude is calculated as where we use the effective Hamiltonian described in Eq. ( 18).As we prepare a superposition of different energy eigenstates by a non-adiabatic transition and then perform a projective measurement onto another superposition of the energy eigenstates, the difference between the energy eigenvalues affects the oscillation.

Appendix C: Results for strong non-adiabatic transitions
In the main text, we considered a case in which the non-adiabatic transition is not relevant.In this section, we investigate the performance of our scheme when we increase the effect of the non-adiabatic transitions in case E.
We plot the spectrum by setting T ann = 3 and s 1 = 0.9, as shown in Fig. 12 (a).Here, as a visual guide, we plot a blue line corresponding to the analytical curve of Eq. ( 22), where we use the actual values of | ⟨1| Ḣ|0⟩ | and (E 1 − E 0 ), and this is the target peak in the spectrum.The magnified view is shown in Fig. 12 (b).
We observe unexpected peaks at frequencies of around Ω = 0.38 in the spectrum, and their height is more significant than that of the target peak.Therefore, if we naively adopt our method described in Eq. ( 31), we obtain incorrect estimated values of the transition matrix element and energy gap.
To understand the origin of these unexpected peaks, we perform analytical calculations in order to obtain res-onant frequencies in the spectrum with non-adiabatic transitions in Appendix B 2. Although the peaks at Ω = ω − Ω (0,1) ana (ω), ω, ω + Ω (0,1) ana (ω) should exist, we could not observe them owing to the restricted range of Ω in the spectrum, as mentioned in Appendix B 2. Meanwhile, if there is a non-negligible population of the second excited state induced by the non-adiabatic transitions, we can observe peaks at frequencies of Ω = ω − Ω ana (ω).To see these points, we plot the spectrum in Fig. 13, and we actually observe three such peaks.Moreover, when ω is far from the energy difference (E 2 − E 1 ), the peak at a frequency of ω − Ω  In case E, we plot the power spectrum P (ω, Ω), where we set Tann = 3 and t1/Tann = 0.9 in (a).The magnified view is shown in (b).The blue and yellow lines represent the target peaks obtained by diagonalization.The peak we want is seen at the correct position when Ω is enlarged in the small region; however, it is lower than the high-frequency peak around Ω = 0.38.Even when we observe peaks other than the target peaks, there is a way to estimate the transition matrix element and energy gap.As mentioned in the main text (in Section IV C 2), we adopt a modified method to identify the target peak, which is useful for this case as well.The positions of the target peak are expected to be fitted by the analytical formula in Eq. ( 22).Thus, if we fail to fit the peaks, we can guess that such peaks do not correspond to the peak from the Rabi oscillation.Indeed, the peak Ω ≃ 0.38 cannot be well fitted by Eq. ( 22).Meanwhile, if we focus on the peaks with frequencies of around Ω = 0.01, as shown in Fig. 12 (b), we can fit these peaks by the analytical formula; hence, we can accurately estimate the transition matrix element and energy gap.
Furthermore, we plot the spectrum by setting T ann = 3 and s 1 = 0.7 in Fig. 14 (a).Here, as a visual guide, we plot a yellow line corresponding to the analytical curve of Eq. ( 22), which are the target peaks.Here, the target peaks as well as other peaks are observed.These come from higher-order perturbations, which can be observed for a larger Rabi frequency (see Appendix D).We can distinguish these secondary peaks from the target peaks as follows.
First, the secondary peaks are usually smaller than the target peaks.As shown in Fig. 14, except for a few points, the target peaks are the highest in this frequency region.Second, from the fitting results by Eq. ( 22), we can distinguish the target peaks from the secondary peaks (for example, the slope of the target peaks is given by dΩana(ω) dω ≃ 1 for large ω, while the slope of the secondary peaks is twice as large).Third, the λ-dependence of the peak height is different.The height of the target peaks scales as λ 2 , while that of the secondary peaks scales as λ 4 , and this lets us identify the target peaks by sweeping λ.As shown in Fig. 14, we plot the spectrum by select-ing a smaller λ, and we show that the secondary peak becomes nearly invisible compared to the target peak.FIG.14. Plot of the power spectrum in case E with Tann = 3 and t1/Tann = 0.7.Here, as a visual guide, we plot a yellow line corresponding to the analytical curve of Eq. ( 22), which are the target peaks.(a) We select λ/Tann = 0.05.We observe the target peaks as well as the peaks due to higher-order perturbations.(b) We select λ/Tann = 0.01.Compared to the case with λ/Tann = 0.05, the peaks caused by higher-order perturbations are smaller.

Appendix D: Perturbative approach of this method
In the main text, we used RWA, and the Hamiltonian was effectively transformed into a simple two-dimensional one.When we consider higher-order perturbations, there can be other peaks without the |E 1 − E 0 | = ω condition.
We calculate the transition amplitude ⟨f |U (t, 0)|i⟩, where U (t, 0) is a unitary operator expressing the time evolution from time 0 to time t, and |i⟩ and |f ⟩ are the initial state and the final state, respectively.
First, we describe the Hamiltonian of Eq. ( 7) in the so-called interaction picture: Then, the transition amplitude is given by where the tilde symbol implies that the states are in the interaction picture.As the Schroedinger picture coincides with the interaction picture at t = 0, we have | ĩ⟩ = |i⟩.We assume that λ is small, and we calculate the transition amplitude up to the second order of the perturbation theory.
In our method, we assume that the dynamics is adiabatic in the second and fourth steps.In this case, we can set |i⟩ = |0⟩ and |f ⟩ = |m⟩, and we evaluate the quantity of p 0,m (t) = | ⟨m|U (t, 0)|0⟩ | 2 .

First-order perturbation and the Rabi oscillation
Let us consider the first order of the perturbative term in Eq. (D2) as follows: The absolute square of Eq. (D3) includes terms with frequencies of 2ω, E m − E 0 ± ω.This result is consistent with the analytical result in Eq. ( 22) in the limit of small λ, which is used to predict the resonance at Ω = |E m − E 0 − ω|.The absolute square of this amplitude includes various modes; however, one of them is E m − E 0 − 2ω.Thus, the probability function p 0,m (τ ) includes an oscillation with a frequency of Ω = 2(ω − (E m − E 0 )/2).Importantly, for this frequency, we have dΩ dω = 2, while we have dΩana dω = 1 for our analytical formula in Eq. (22).The adiabatic condition (1) is believed to be related to success of quantum annealing but it is not proven to be sufficient.The adiabatic condition Eq. (1) comes from the equation (A11), and the final transition amplitude to the excited states is approximately given by (A15).It is worth mentioning that, since a perturbative theory is used, the Eq.(A15) is not exact.Although we cannot show the validity of the Eq.(A15) for any quantum annealing, we can show that the Eq.(A15) actually provides us with an upperbound of the population of the excited states for our numerical examples in section IV.According to Eq. (A15), the population of the n-th excited state is evaluated by  We perform a numerical simulation of the quantum annealing defined by (2) with the two qubit model (35), and calculate the populations of the first, the second, and the third excited states at the time s = 1.Then, we plot the dependence of these populations on the annealing time T ann .(see Fig. 15) Fig. 15 shows that sum of p n defined by (E1) is actually larger than the populations obtained by the quantum annealing for the large T ann .Also, we confirm that the populations are inversely proportional to the square of the annealing time T ann in the limit of the large T ann .If T ann is small, the approximation used in (A10) is violated and p n does not provide the upper bound of the population.Thus, the adiabatic condition we used ( 1) is expected to be valid for the large T ann .

Appendix F: Numerical results for larger systems
To establish the applicability of our method for larger systems, we adopt a nine-qubit system with the following Hamiltonian.
where σ i x (σ i z ) is the Pauli X (Z) matrix acting on ith qubit, J i,i+1 is a coupling strength between nearest neighbor qubits, h i is a strength of the external field.Actual values of these parameters are shown in Table II.
Let us explain how we can apply our method to this system.During the implementation of QA with the driver Hamiltonian and problem Hamiltonian described above, a notable observation emerges where the success probability suddenly decreases at a specific point in the process.(see Fig. 16).We apply our method to measure the adiabatic condition at the point, and we obtain a power spectrum depicted as Fig. 17.As depicted in Fig. 17, while the spectrum displays numerous curves, the curve corresponding to the desired signal is discernible, and so we can fit this by using the analytical expression of Eq. (22).Therefore, we can estimate the value of the right-hand side of the Eq. ( 1) from the spectrum, and this information is valuable when we try to optimize the schedule of QA to maximize the success probability.22) with diagonalization of the Hamiltonian Eq. (F1).We select λ/Tann = 0.01 in this simulation.

FIG. 1 .
FIG. 1. Plot of the scheduling function A(t) and strength of the external driving field λ for our protocol.The dotted line shows f (s), which is the scheduling function of the conventional QA.

FIG. 2 .
FIG. 2. The actual adiabatic conditions (1) of our simulated two-qubit cases.m = 1 is the largest at almost all s.

FIG. 3 .
FIG. 3. (Top) Estimation of the transition matrix element in case A (single qubit, complete adiabaticity, and no decoherence).(Bottom) Estimation of the energy gap in case A. The solid lines represent the solution obtained by diagonalization of the Hamiltonian, and the dots represent the estimated values obtained from our method by numerical simulation.

FIG. 5 .
FIG. 5. Top (bottom): estimated value of the transition matrix element (energy gap) in case B (single qubit, incomplete adiabaticity, and no decoherence) with a shorter annealing time such as Tann = 1, 2, 4, 8 ns.For the solid lines and dots, we use the same notation as that in Fig. 3.

FIG. 10 .
FIG. 10.Plot of the power spectrum P (ω, Ω) for Tann = 100 for case F (two qubits, incomplete adiabaticity, and decoherence).The horizontal axis represents the angular frequency of the driving field and the vertical axis represents the Fourier frequency.(a) We use s1 = 0.3.The yellow line is obtained by fitting Eq. (22) to the plot.(b) We use s1 = 0.9.The dotted line represents the exact value obtained by diagonalization.

FIG. 12 .
FIG. 12.In case E, we plot the power spectrum P (ω, Ω), where we set Tann = 3 and t1/Tann = 0.9 in (a).The magnified view is shown in (b).The blue and yellow lines represent the target peaks obtained by diagonalization.The peak we want is seen at the correct position when Ω is enlarged in the small region; however, it is lower than the high-frequency peak around Ω = 0.38.

2 .
Second-order perturbation Let us consider three energy eigenstates of H QA as |0⟩ , |I⟩ , |m⟩, and we assume that |i⟩ = |0⟩ and |f ⟩ = |m⟩.The third term of the right-hand side of Eq. (D2) is given

FIG. 15 .
FIG. 15.Plot of the population of being one of the excited states (solid line) and sum of pn(1) (Dotted line) against the annealing time Tann for n = 1, 2, 3.For the large Tann, we confirm that pn(1) actually gives an upper bound of the population in these cases.

FIG. 16 .
FIG.16.Success probability of our nine-qubit annealing process with respect to s. Around s = 0.4, success probability is suddenly worsen.
We will demonstrate this point below.Let us assume that we are interested in only two states, |k⟩ and |l⟩.In this case, we can approximate the Hamiltonian as H QA ≃ ∆ 2 σ z = |l⟩⟨l| − |k⟩⟨k|.Furthermore, we can use Eq.(B2) for the effective Hamiltonian.The transition amplitude (B10) is calculated as