Reduction of the energy-gap scaling by coherent catalysis in models of quantum annealing

Non-stoquastic drivers are known to improve the performance of quantum annealing by reducing first-order phase transitions into second-order ones in several mean-field-type model systems. Nevertheless, statistical-mechanical analysis shows that some target Hamiltonians still exhibit unavoidable first-order transitions even with non-stoquastic drivers, making them difficult for quantum annealing to solve. Recently, a mechanism called coherent catalysis was proposed by Durkin [Phys. Rev. A \textbf{99}, 032315 (2019)], in which he showed the existence of a particular point on the line of first-order phase transitions where the energy gap scales polynomially as expected for a second-order transition. We show by extensive numerical computations that this phenomenon is observed in a few additional mean-field-type optimization problems where non-stoquastic drivers fail to change the order of phase transition in the thermodynamic limit. This opens up the possibility of using coherent catalysis to search for exponential speedups in systems previously thought to be exponentially slow for quantum annealing to solve.


I. INTRODUCTION
Quantum annealing (QA) [1][2][3][4][5][6][7][8][9] is a method for solving combinatorial optimization problems by making use of the tunneling effects induced by quantum fluctuations to search for the optimal solution. An important challenge facing QA is that for certain classes of problems one encounters a first-order phase transition during the annealing process. At a first-order transition, the energy gap ∆ = E 1 − E 0 between the energies of the ground and first-excited states of the system, E 0 and E 1 , generally decreases exponentially with system size, resulting in an exponentially long annealing time, according to the adiabatic theorem of quantum mechanics [10][11][12]. One such example is the ferromagnetic p-spin model (p-spin ferromagnet) whose target Hamiltonian is given by (1) where the power p ≥ 2 is an integer, N is the total number of spins, J α = N −1 N i=1 σ α i , and σ α i (α = x, y, z) is the α-component of the Pauli matrix for the ith spin. Jörg et al. [13] showed that, if QA is performed with just a transverse field J x , the system undergoes a firstorder phase transition for p > 2 and is therefore hard to solve although the solution is trivially ferromagnetic. It was later found that, by introducing non-stoquastic drivers of the form +J x q (q ≥ 2) into the annealing Hamiltonian, many of the first-order phase transitions that plagued a simple transverse field annealing in the p-spin ferromagnet can be avoided and be replaced by gentler second-order ones whose annealing times are expected to scale only polynomially with the system size * patrickkyw@gmail.com [14][15][16][17][18]. An intractable exception remains, however, for the case of p = 3 where mean-field analysis, which is valid in the thermodynamic (large-size) limit, suggests that first-order transitions cannot be avoided even with a non-stoquastic driver [14,15,17]. It is also known that the p-spin ferromagnet with p ≥ 4 becomes difficult to solve even with a non-stoquastic driver in the presence of random longitudinal fields [19].
In a recent development, Durkin [20] examined the p = 3 case in detail and found that it is in fact actually possible to reduce the computation time from exponential to polynomial as a function of the system size with the help of a non-stoquastic driver +J x 2 if one carefully controls the system parameters as functions of the system size. By exact numerical diagonalization and mapping of the problem to an effective single-particle problem, he showed that, along the curve of minimum energy gap that separates the paramagnetic and ferromagnetic phases in the phase diagram, there exists a particular point where the gap scales polynomially with the system size. At this point on the phase transition curve, kinetic effects due to the transverse antiferromagnetic interaction (non-stoquastic driver) delocalizes the ground-state wavefunction into a form that simultaneously overlaps the two minima of the classical potential, thereby softening the first-order phase transition into an effectively second-order one. This process is termed 'coherent catalysis' by the author, and opens up the possibility of searching for quantum speedups in systems which are predicted by the mean-field theory to exhibit only first-order transitions.
In Ref. [20], the use of coherent catalysis to change the gap scaling from exponential to polynomial in QA was demonstrated exclusively for the p = 3 model. Since this result is a consequence of intricate balance among various terms in the Hamiltonian, a natural questions arises whether or not the mechanism of coherent catalysis applies to other problems. The goal of the present paper is to answer this question. The models we have chosen to study here are difficult for QA in the sense that the conventional mean-field theory predicts that a first-order phase transition is unavoidable even in the presence of non-stoquastic drivers. Our approach is computational and we shall focus mainly on computing the energy gap by diagonalizing the Hamiltonian numerically and obtaining its scaling. Our results indicate that the mechanism of coherent catalysis works also in our problems.
The rest of the paper is organized as follows. Sections II to IV are devoted to three different models. In Sec. II we extend the Hamiltonian studied in Ref. [20] to include many-body transverse interactions. This is to verify that the result reported in Ref. [20] is not just peculiar to the case of +J x q with q = 2 but holds for higher powers of q as well. In Sec. III we go beyond the simple 3-spin ferromagnet and study a target Hamiltonian whose energy landscape exhibits a barrier and a local minimum. In Sec. IV we study the p-spin ferromagnet in a random longitudinal field. A previous mean-field study of this model has shown that the random longitudinal field induces a first-order transition that is not lifted by a nonstoquastic driver even for p ≥ 4 where the difficulty can be removed by a non-stoquastic driver in the absence of random field [19]. Section V discusses and concludes the paper.

II. THE 3-SPIN FERROMAGNET WITH MANY-BODY TRANSVERSE INTERACTIONS
The first model we study is the 3-spin ferromagnet with many-body transverse interactions given by the Hamiltonian where −J z 3 is the target Hamiltonian, Γ and κ are parameters controlling the strengths of the transverse field J x and the non-stoquastic driver term +J x q , respectively, and the integer power q ≥ 2 generates the q-body transverse interactions. We have used the above parameterization by κ and Γ following the convention adopted in Ref. [20]. We have also chosen the Hamiltonian to be of order N 0 in the present section for consistency with Ref. [20].
When q = 2, H I reduces to the model studied in Ref. [20]. For QA, one starts from Γ = 1 and κ arbitrary and tries to find a path in the Γ-κ parameter space to reach (Γ, κ) = (0, 1) where the target classical Hamiltonian is recovered. A mean-field analysis of this system shows that the Γ-κ plane is separated into two regions by a firstorder transition line [15]. Although the precise location of this line depends on q, qualitatively, it lies within the curved orange-yellow segment of the heat map shown in Fig. 1(a).   (a) Density plot of the logarithm of the inverse of the energy gap, log 10 (1/∆), as a function of the parameters κ and Γ for q = 3 and N = 50. The location of the saddle point on this 'gap landscape' is shown as a green circle. Equipotential contours of log 10 (1/∆) are overlaid to aid in the visualization of the saddle point. Annealing starts from Γ = 1, κ arbitrary and ends at Γ = 0, κ = 1. (b) The energy gap evaluated at the saddle point, denoted ∆D, decreases polynomially as the system size N increases. Results for q = 2 (studied in Ref. [20]) and higher powers q = 3 and 4 are shown. The solid black line corresponding to each q is obtained by fitting a straight line through all the data points.
According to the mean-field theory, during the annealing process the system must cross this line and undergo a first-order phase transition, and therefore the gap is naively expected to decrease exponentially as a function of the system size. This problem with general q was touched upon in Ref. [20] but was not analyzed in detail except for the case of q = 2. Another reason that we study this model is to confirm the accuracy of our numerical method to discern exponential and polynomial gap closings, which will turn out to be serious also in other problems to be discussed later.
The Hamiltonian H I commutes with the total angular momentum operator J x 2 + J y 2 + J z 2 . Since quantum annealing starts from the state with only the transverse field term Γ = 1 in the Hamiltonian Eq. (2), one can simply diagonalize H I in the sector with total angular momentum N/2. The result is summarized in Fig. 1(a) as the density plot of the logarithm of the energy gap as a function of Γ and κ for the case of q = 3 with N = 50. Qualitatively, the 'gap landscape' here is similar to that of q = 2 reported in Ref. [20]. The narrow orange-yellow curve slanting down the middle is where the energy gap is smallest, and separates the ferromagnetic phase on the left from the paramagnetic phase on the right. In the limit N → ∞, the minimum gap curve asymptotically approaches the first-order phase transition curve given by the mean-field theory [15]. It therefore seems impossible to avoid a first-order transition to reach the final state at Γ = 0, κ = 1 from the initial state at Γ = 1, κ arbitrary, suggesting an exponentially long computation time.
An outstanding feature about this gap landscape is the existence of a saddle point along the minimum gap curve where the gap is a maximum as pointed out for the case q = 2 in Ref. [20]. In Fig. 1(a), this point is indicated by a green circle. It was reported for q = 2 in Ref. [20] that the gap evaluated at this point obeys a different, polynomial, scaling compared to the other points with exponential scaling along the minimum gap curve. To see if the same happens in the present case of q = 3, we tracked the saddle point as N increases, which is known to scale as κ c ∝ 1/ √ N for κ 1 for q = 2. Numerical identification of a saddle point in the very steep gap landscape is a highly non-trivial task, and we used the Newton-Raphson algorithm to find the point (Γ D , κ D ) where the first derivatives (with respect to Γ and κ) of the function log 10 (1/∆) vanishes. The gradient and Hessian required are also computed numerically, by taking finite differences of log 10 (1/∆) along the Γ and κ directions. As the initial condition for the Newton-Raphson algorithm, it is preferable to choose a point which is near to the solution (such as the saddle point from the previous N ).
The results for the scaling of the gap evaluated at the saddle point, denoted as ∆ D , are shown in Fig. 1(b) for q = 2, 3, and 4 (from N = 50 to 500). The three solid lines shown are obtained by fitting the data points (from all N ) to a polynomial function of the form c 1 N −c2 (c 1 , c 2 : fitting parameters). It is seen that ∆ D decreases polynomially with increase in system size for all three q values. The scaling result we obtained for q = 2 also agrees with the analytic result 3.4 × N −2 derived in Ref. [20] (which holds in the large N limit) [21]. Elsewhere along the minimum gap curve, the energy gap closes exponentially fast with N , in accordance with the prediction of the mean field theory. This has been discussed for the case of q = 2 in Ref. [20], and which we also verified for the cases of q = 3 and 4 based on our numerical calculations. We have thus confirmed that the method of coherent catalysis works also in this slightly different problem and that our numerical technique succeeds to adequately follow the saddle point in the gap landscape where the valley of minimum gap has a very sharp, edge-like, shape.

III. A TARGET HAMILTONIAN EXHIBITING AN ENERGY BARRIER
One may ask if the very simple structure of the target Hamiltonian −J z 3 with a single trivial minimum may be the main culprit of the success of coherent catalysis. We therefore study a model with a local minimum separated from the global minimum by an energy barrier in the target Hamiltonian. The total Hamiltonian is given by where g(J z ) is the target Hamiltonian, and s and λ are parameters controlling the strengths of the transverse field J x and the antiferromagnetic interaction (nonstoquastic driver) +J x 2 , respectively. The parameterization by s and λ follows the convention of Refs. [14,15,17]. The function g(m) is given by and its graph is shown in Fig. 2 [22]. The parameters a and b are the coordinates of the local minimum and local maximum, respectively. The height of the energy barrier measured relative to the potential at m = a is (b − a) 3 . For comparison, the graph of the target Hamiltonian of the 3-spin ferromagnet is shown in the inset. The latter is a monotonic function with only one minimum solution at m = 1. From an algorithmic point of view, finding the optimal minimum of g(m) should be more non-trivial than that of −m 3 due to the presence of the local minimum. For the QA of H II , one starts from s = 0, λ arbitrary and finishes at (λ, s) = (1, 1) where the target is recovered. A mean-field analysis of this model once again predicts that the system must encounter an unavoidable first-order phase transition during the annealing process [23]. In the following, we examine this model from the perspective of coherent catalysis.
A. Saddle point with polynomial scaling of the energy gap The Hamiltonian H II again commutes with the total angular momentum operator, so the numerical methods used in Sec. II are applicable to this model. Figure 3(a) shows the plot of log 10 (1/∆) as a function of s and λ for the parameter values a = 0.6, b = 0.8 with N = 100. An interesting observation is that, unlike the first model, we were able to identify two saddle points on the gap landscape: One inside Box (i) and another inside Box (ii). We first discuss the left one in Box (i) which has a larger gap. Figure 3 We next discuss the right saddle point lying inside Box (ii) of Fig. 3(a). Figure 3(c) shows a close-up view of the gap landscape in its vicinity. This saddle has a smaller gap compared to the left one and lies between two minima on the landscape. The right (blue) curve in the inset of Fig. 4 shows how the location of this saddle evolves on the λ-s plane as N increases. The blue squares in fitting parameters) through the data points. It is seen that ∆ s decreases exponentially with increase in system size N . Hence, the computation time remains exponential if the annealing path were to pass through the right saddle point.

C. Smaller energy barrier in the target Hamiltonian
The gap landscape of this model is dependent on the parameters values a and b chosen. To illustrate it, consider a = 0.7, b = 0.78, which gives a smaller energy barrier in the target Hamiltonian compared to the case of a = 0.6, b = 0.8. We found that although the left saddle point still exists, this time the right saddle point does not. The scaling of the energy gap at this left saddle is calculated and the result is plotted in Fig. 4 using green circles. The solid black line passing through the circles is obtained by fitting a polynomial function through the data points. It is seen that the gap again decreases polynomially with increase in N . The rate of decrease is also slower compared to the previous case of a = 0.6, b = 0.8, which is consistent since a target Hamiltonian with a smaller barrier should be easier to solve.

IV. THE 3-SPIN FERROMAGNET IN RANDOM LONGITUDINAL FIELD
The third model we study is the 3-spin ferromagnet in a random longitudinal field. This is a more difficult problem than the case without a random longitudinal field discussed in Ref. [20] and Sec. II because the non-stoquastic driver fails to reduce first-order transitions into second order even in the p-spin ferromagnet with p ≥ 4 [19]. Remember that the reduction of the order of transition is achieved by non-stoquastic driver for the p(≥ 4)-spin model in the absence of a random longitudinal field [14,15,17].

The target Hamiltonian is
where h 0 is the strength of the longitudinal field given by The random variables ξ i are each +1 or −1 with probability 1/2. When h 0 = 0, we recover the 3-spin ferromagnet, which QA can solve in polynomial time using coherent catalysis. In this section, we shall discuss the case where exactly N/2 of the ξ i are +1 (we consider even N in this section). Our annealing Hamiltonian is where the parameters s, λ as well as the starting and ending points of the annealing schedule are the same as for the previous model in Sec. III. According to an earlier mean-field study, the system must undergo an unavoidable first-order phase transition during the annealing process [19].
Unlike the two previous models, the Hamiltonian H III does not commute with the total angular momentum operator due to the longitudinal field L z . To avoid diagonalizing the entire Hamiltonian matrix of dimension 2 N , we note that H III exhibits permutation symmetry among the spin indices with the same value of ξ i . We renumber the indices i such that i = 1 to N/2 now have ξ i = +1, and i = N/2 + 1 to N now have ξ i = −1, which is possible because of the symmetry of the first term −J z 3 in the target Hamiltonian Eq. (5) under the permutation of spin index. The Hamiltonian can be rewritten as In the tensor product A ⊗ B, the operator A(B) lies in the Hilbert space of the spins whose ξ i is +1(−1).Ĩ denotes the identity operator, The Hamiltonian H III commutes with the total angular momentum operators of each subspace, αJ α 2 ⊗Ĩ andĨ ⊗ αJ α 2 . To obtain just the ground and first excited states, one can diagonalize H III using the product basis |N/4, m A ⊗ |N/4, m B where |j, m is an angular momentum eigenstate with total quantum number j and projection quantum number m. The dimension of the Hamiltonian matrix that we need to diagonalize is therefore (N/2 + 1) 2 instead of 2 N . In our calculations, we diagonalized H III using the Lanczos algorithm and computed only the ground and first excited states.
A. Absence of coherent catalysis when h0 = 1 The valley in the gap landscape in the random-field model is very sharp and numerical computations require extremely careful steps. Figure 5(a) shows the density plot of the logarithm of the energy gap log 10 ∆ as a function of s and λ for h 0 = 1 with N = 50. For this model, we cannot identify the saddle point on the gap landscape just by visual inspection because the minimum gap curve is sandwiched within an extremely narrow region of the λ-s parameter space [24]. To locate the saddle point, we first compute the minimum gap curve [25], which is shown in Fig. 5(a) by the broken red curve. The gap evaluated along the curve, denoted as log 10 [∆ min (λ)], is shown in Fig. 5(b). It is seen to be rich in structures such as local maxima and cusp-shaped minima. We shall focus on the global maximum indicated by a blue circle, which is also a saddle point. Inset (i) shows a close-up view of the gap in the vicinity of the maximum. To verify that the maximum is a saddle point, we took the gradient ∇ = (∂/∂λ, ∂/∂s) of the function log 10 ∆(λ, s) along the minimum gap curve [26]. Inset (ii) shows that the Euclidean norm of the gradient, ||∇(log 10 ∆)||, vanishes at the maximum. The scaling of the energy gap evaluated at this saddle point, denoted as ∆ s , is plotted in  Fig. 6(a) using solid red circles. The solid black curve passing through the circles is obtained by fitting an exponential function through the data points. It is seen that ∆ s decreases exponentially with increase in N . Hence, it is not possible to reduce the gap scaling from exponential to polynomial using coherent catalysis when h 0 = 1.
B. Crossover from exponential to polynomial scaling between h0 = 1 and 0 As mentioned earlier, when h 0 = 0 the target Hamiltonian Eq. (5) reduces to the 3-spin ferromagnet whose energy gap at the saddle point exhibits polynomial scaling. This is shown in Fig. 6(a) using solid blue triangles; the solid red line passing through the triangles is obtained by fitting a polynomial function through the data points. The scaling of the saddle point gap for h 0 = 0.25, 0.50, and 0.75 is also shown in Fig. 6(a). Figure 6(b) shows a detailed view of the scaling curves, with values of h 0 close to 0. The gap scaling seems to change from polynomial to exponential as h 0 is increased from 0 to 1.
To quantify this observation, let us fit the scaling data of the gap ∆ s of each h 0 to a polynomial curve of the form log 10 ∆ s = c log 10 N + d The fitting parameters c and d obtained are shown in the inset of Fig. 7(a). How well the data points are fitted by the polynomial function Eq. (9) is given by the sum of the square of residuals [27], which is shown in Fig. 7. It is seen that for h 0 0.25 the polynomial function fits the data very well; as h 0 increases beyond 0.25 the residual increases rapidly, indicating a deterioration of the fit. For comparison, we repeated the curve fitting, this time to an exponential function of the form The fitting parametersc andd obtained are shown in the inset of Fig. 7(b), and the sum of the square of residuals are shown in Fig. 7. It is seen that although for h 0 0.5 the exponential form fits the data well, for h 0 < 0.5 the residual increases rapidly, indicating a deterioration of the fit for small values of h 0 . To conclude, we found evidence that there exist a range of values of h 0 where the scaling of the saddle point gap is likely to be polynomial. Coherent catalysis works to reduce computational complexity in this difficult problem.

V. SUMMARY AND DISCUSSIONS
We have studied the possibility of avoiding an exponentially closing energy gap implied by first-order phase transitions in quantum annealing using the mechanism of coherent catalysis recently proposed by Durkin [20]. The models studied here are difficult for QA in the sense that even with non-stoquastic drivers, the mean-field theory predicts that a first-order phase transition is unavoidable during the annealing process. Coherent catalysis, on the other hand, predicts the existence of a saddle point on the energy gap landscape where the exponential gap scaling is softened into polynomial. Our purpose is to investigate, via numerical calculations, whether this phenomenon is not specific to the 3-spin ferromagnet but also occurs in other systems as well. In all the three models studied, we have found the effectiveness of coherent catalysis, i.e. the existence of a saddle point in the phase diagram where the gap scales polynomially though the thermodynamic phase transition is of first order. It was argued in Ref. [20] that coherent catalysis represents intrinsic quantum effects in the mean-field-type models, and thus the present results may be understood as addi-tional examples of quantum enhancement in optimization algorithms.
In this work, we have approached coherent catalysis purely from an empirical (i.e. numerical) perspective. Some theoretical aspects of the mechanism, however, remain unclear. For instance, we have seen from the models in Secs. III and IV that a saddle point on the gap landscape can exhibit polynomial or exponential gap scaling, not just the former. Hence, the existence of such a point does not automatically imply the reduction of the scaling, and it would be interesting to understand the conditions under which polynomial scaling will occur. Another issue concerns the type of models that has been investigated. Both here and in Ref. [20], coherent catalysis was demonstrated mainly using mean-field type models. An interesting question would be if it can occur in finitedimensional or spin-glass systems as well. We hope to address some of these questions in our future efforts.