Counterdiabatic driving in the quantum annealing of the $p$-spin model: a variational approach

Finding the exact counterdiabatic potential is, in principle, particularly demanding. Following recent progresses about variational strategies to approximate the counterdiabatic operator, in this paper we apply this technique to the quantum annealing of the $p$-spin model. In particular, for $ p = 3 $ we find a new form of the counterdiabatic potential originating from a cyclic ansatz, that allows us to have optimal fidelity even for extremely short dynamics, independently of the size of the system. We compare our results with a nested commutator ansatz, recently proposed in P. W. Claeys, M. Pandey, D. Sels, and A. Polkovnikov, Phys. Rev. Lett. 123, 090602 (2019), for $ p = 1 $ and $ p = 3 $. We also analyze generalized $ p $-spin models to get a further insight into our ansatz.


I. INTRODUCTION
It is possible to solve a complex optimization problem by finding the ground state of a Hamiltonian H p that encodes the problem of interest. Adiabatic quantum computation [1] and quantum annealing [2,3] achieve this task by slowly dragging the system to the desired state, having first initialized it in the ground state of a known (simple) Hamiltonian H q . The adiabatic theorem assures that a time-dependent Hamiltonian H 0 (t), such that H 0 (0) = H q and H 0 (T) = H p , will lead towards the wanted solution if, at t = 0, the system is prepared in the ground state of H 0 (0) and if T is large enough to avoid transitions to excited states.
The ultimate bound to the efficiency of this approach is related to the presence, and to the scaling with system size, of vanishing small gaps in the instantaneous spectrum of H 0 (t). The evolution time must be tuned so that at t = T there is a large probability of finding the system in its ground state (longer annealing times obviously leave also the system more prone to errors and thermal noise [4]). According to the adiabatic criterion, this means that T must be longer than /∆ 2 , where ∆ is the minimal gap between the ground state and the first excited state. When the spectral gap closes in the thermodynamic limit, as in the case of a quantum phase transition [5], the evolution time T must grow with the system size as well, to compensate for the closure of the minimal gap. The hardness of a given optimization problem is directly reflected in the presence of exponentially small (with the system size) gaps during the adiabatic evolution.
Finding ways to reduce the computation time T, while keeping the fidelity of the time-evolved state (at t = T) close to the ground state of H p , is a central question in this context that still needs an adequate answer. In the last years, there have been several attempts to circumvent this problem, resorting to quantum optimal control approaches applied to the evolution of a complex system. A particularly appealing approach that received increasing attention recently is the implementation of shortcuts to adiabaticity [6][7][8] to quantum annealing.
A Shortcut To Adiabaticity (STA) [or Counterdiabatic Driving (CD)] amounts in finding a time-dependent Hamiltonian which evolve the state of a system as if the evolution was adiabatic. Given the bare Hamiltonian H 0 (t), it is possible to find a term H cd (t) such that that the dynamics governed by H(t) = H 0 (t) + H cd (t) has vanishing diabatic transitions between pairs of energy eigenstates at all times [7]. Since its initial formulation, STAs have found numerous applications. An overview of the field can be found in Ref. [9].
The exact CD Hamiltonian H cd (t), derived by Berry [7], can be expressed in terms of the exact instantaneous eigenstates and eigenvalues of H 0 (t). In a many-body system, H cd (t) becomes progressively nonlocal and hence not easy for a practical implementation, when H 0 (t) is close to a quantum phase transition [10,11]. More importantly from a conceptual point of view, the severe limitation in the use of the exact form of H cd (t) arises because i) exact eigenstates cannot be determined (this would correspond to a prior knowledge of the solution of the optimization), and ii) H cd (t) may be ill-defined due to small denominators near the quantum critical point.
A big step forward in overcoming these limitations has been achieved by Sels and Polkovnikov [12], who derived a variational principle to determine approximate forms of H cd . This new variational approach to STA solicits detailed scrutiny of its potentialities in speeding up quantum annealing protocols. Some very recent papers addressed this question [13][14][15]. Hartmann and Lechner [13] analyzed the approximate optimal counterdiabatic driving in the so called LHZ lattice gauge model [16], Claeys et al. [14] studied a non-integrable onedimensional Ising model, Hatomura [15] adopted a mean field approximation for the CD driving of the infinite-range Ising model.
The paper is organized as follows. In Section II, we introduce the basic concepts of STA and describe the variational procedure developed in Ref. [12]. In the same section, we show the approximate counterdiabatic operator adopting two different ansätze for the variational potential: The first one is based on the Nested Commutator (NC) approach [14], the sec-ond one is based on a simple Cyclic Ansatz (CA), the origin of this name will be clear in the following. In Section III, we discuss different properties of the p-spin model for p = 1, 2 and 3 in detail. This section contains a detailed discussion of the results we obtained. We compare the two approaches and show that, in the case p = 3, the CA outperforms the NC ansatz considerably, leading to an almost perfect recovery of the ground state in the final time T. This effect occurs even for very large system sizes, the scaling of the fidelity being almost independent of the number of spins. Given the excellent performance of the CA for the p-spin model, in Section IV, we move away from this highly-symmetric case and consider different generalizations that do not have constant all-to-all couplings. The question is to see if the CA has sufficient power and flexibility to be used in a generic situation. To this aim, we consider different finite-range and random variants of the p-spin model and compare the two variational ansätze also in these cases. Finally, in the Section V, we summarize our results and briefly discuss future directions.

II. QUANTUM ANNEALING AND COUNTERDIABATIC DRIVING
As already mentioned in the introduction, a typical annealing schedule can be formulated as follows. A (many-body) quantum system is governed by a time-dependent Hamiltonian of the form with the constraint that H 0 (0) = H(0) and H 0 (T) = H(T), which implies λ(0) = 0, λ(T) = 1. In the rest of the paper we will further impose that λ(0) = λ(T) = 0. A convenient parametrization that automatically includes all these requirements has the form The goal is to reach, as close as possible, the ground state of H(T) if the initial state is the ground state of H(0). Having introduced the instantaneous eigenstates of H 0 (λ(t)) as H 0 (λ(t)) | m (t) = m | m (t) , ideally we require that

II.1. Shortcut to adiabaticity
Achieving the condition in Eq. (3) is the goal of adiabatic quantum computation via quantum annealing. If one wants to speed up the evolution, an additional control on the dynamics is required. The proposals collectively named shortcuts to adiabaticity or counterdiabatic drivings [6][7][8]33], reviewed in Ref. [9], are based on the idea that the evolution of a quantum state can coincide with the instantaneous ground state of H 0 (λ(t)) if an additional time-dependent contribution H cd (t) is added to the Hamiltonian. The exact quantum state governed by the Hamiltonian H 0 (λ(t))+H cd (t) is given by |ψ(t) = | 0 (t) at any time (here, for simplicity, we considered the case in which one follows the ground state).
Berry [7] derived the exact CD potential to add to the bare Hamiltonian H 0 (t) to completely suppress diabatic transitions between pairs of energy eigenstates at all times. It reads Albeit exact, this potential is ill-defined at the quantum critical point due to exponentially small denominators. Furthermore, even in the cases where Eq. (4) can be analytically computed, the resultant operator can be highly nonlocal and impossible to implement on the available hardware. Finally, and most important for its application to quantum annealing, the Hamiltonian in Eq. (4) requires the knowledge of the exact eigenvalues and eigenvectors of the Hamiltonian at any time. A big step forward in offering a solution to these problems comes with the variational approach to counterdiabatic driving formulated in Refs. [12,34].

II.2. A variational approach to counterdiabatic driving
A variational principle for the counterdiabatic potential H cd can be derived as follows. If one parametrizes the CD potential as with It is possible to introduce the operator [12,34]: whose diagonal elements in the energy eigenbasis do not depend on A * λ , and whose off-diagonal elements in the energy eigenbasis are zero if A * λ = A λ (i. e., if G λ = −F ad ). Therefore, the true counterdiabatic potential A λ minimizes the Hilbert-Schmidt norm of G λ , i. e., the operator distance between G λ and −F ad : In the next Section we are going to apply the above mentioned variational approach to the p-spin model. As in any variational approach, the next important step is an educated guess of (in this case) the operator A λ . We will consider two variational forms of the operator A λ : 1. An operator A (l) λ as defined in Ref. [14], expressed in the form of NC: It turns out that this form can be manipulated efficiently. Important details of the method are discussed in Appendix A. In principle, the larger is l, the more accurate is the approximation of the CD potential.
2. For p = 1, 3, we also study a new form, that we name CA, where the corresponding variational operator A CA λ has a particularly symmetric form when applied to the model considered here. It is expressed by Details on this ansatz are discussed in Appendix B.
In the next sections, we scrutinize these two possibilities testing to which extent it is possible to optimize the annealing schedule in these cases. The new variational ansatz, especially tailored to this specific model, leads to exceptional good fidelities for the p-spin model.

III. FERROMAGNETIC p-SPIN MODEL
The ferromagnetic p-spin model is defined by the Hamiltonian where S x , S y , S z are total spin operators: For odd p, the ground state of H p is ferromagnetic and nondegenerate, whereas for even p there are two degenerate ferromagnetic ground states, related by a Z 2 symmetry. In the following, we choose Γ = J = 1 and = 1. We consider three prototypical cases of this model [17][18][19][20][21][22][23][24][25][26][27][28][29][30]35]: 1. p = 1: in this case, the qubit system acts as a single spin S = n/2. The operators O 2k−1 are all proportional to S y , as shown below. Thus, one variational parameter is sufficient to recover the limit l → ∞ of the NC ansatz. We will numerically show that this corresponds to the exact counterdiabatic potential, up to numerical errors.
2. p = 2: the system exhibits a second-order quantum phase transition, where the minimal gap ∆ scales as ∆ ∼ n −1/3 . In this case, the CD operator derived within the nested commutator ansatz improves the success probability of quantum annealing by increasing the fidelity as a function of the order l. We show that the number of NC required grows with the number of qubits n.
3. p = 3: the system exhibits non degenerate ground state.
It shows a first-order QPT. The exponent p = 3 is the smallest odd integer for which the p-spin model has this property. Even for this (relatively) simple case, a large number l of NC in Eq. (9) is needed in order to have a significant improvement in the success probability of quantum annealing. Moreover, this number increases with the system size. However, in the following we show that the cyclic ansatz yields an almost perfectlyefficient and size-independent counterdiabatic driving in the symmetry subspace with maximum spin.
Note that the p-spin model Hamiltonian is SU(2) invariant, hence states with different total spin S 2 belong to disconnected subspaces. If we apply a quantum annealing protocol using the Hamiltonian of Eq. (11) as target Hamiltonian and H q = − i σ x i as starting Hamiltonian, we can work out the whole procedure within the maximum spin subspace, whose dimension is N = n + 1. This will allow us to perform numerical simulations for large system sizes. Moreover, while in the definition of S l (ì α) of Eq. (8) the traces are evaluated over the 2 n basis states of the full Hilbert space, in this specific case we can define another family of functionalsS l (ì α), in which the traces are restricted to the N = n + 1 states of the maximal spin subspace. In this subspace,S l (ì α) obeys variational equations analogous to those for S l (ì α).
We will show that both the nested commutator ansatz and the CA can be studied in the whole Hilbert space or in the maximum spin subspace, with different performances. These two approaches are compared in Appendix C.

III.1. Results for p = 1
The p = 1 case is trivial, it is however useful to set the stage for further calculations. The Hamiltonian of the p-spin model for p = 1 is It is easy to see that Thus, the ansatz of Eq. (9) contains only one variational parameter (l = 1) and is Details on the minimization of the corresponding quadratic action is given in Appendix A. We numerically simulate the dynamics of the p-spin system in the symmetric spin sector, for T = 1, for system sizes ranging from n = 10 to 100. The time evolution operator time interval [0, T] into N t steps of dt and using the following rule: The probability of being in the instantaneous ground state is and the fidelity is i. e., the probability of being in the ground state at the annealing time t = T. In the absence of the CD term, the fidelity F is very small in the analyzed cases for this choice of T (F ≈ 1 × 10 −3 for n = 10, F < 1 × 10 −15 for n = 50 and above). The scaling of the fidelity as a function of the system size is summarized in Fig. 1, up to n = 100. This clearly shows that the ansatz of Eq. (14) indeed, in this simple case, gives the exact counterdiabatic potential of Eq. (4).
In this case, the fidelity is given by F = P n/2 (T) + P −n/2 (T). For p = 2, we simulate the quantum annealing up to a final time T = 1 and study the scaling of F as a function of the system size in the maximum spin subspace and for different orders of approximation of the counterdiabatic operator A (l) λ of Eq. (9).
In particular, in Fig. 2 we report the scaling of F as a function of n, for l = 1, 3 and 8, compared to bare quantum annealing with no CD terms. These results are obtained using a time In the bare annealing case, the fidelity rapidly goes to zero as we increase the number of qubits n. This is easily understood, as increasing the number of qubits n the energy levels become more and more dense and, for fixed T, the dynamics quickly leaves the adiabatic regime. The starting paramagnetic state is metastable for all t and the p-spin system occupies this state for the whole dynamics. At the annealing time t = T, the system state would be where |n/2 − w are eigenstates of S z with eigenvalues σ w = n/2 − w. Only two terms of this sum contribute to the fidelity: A single variational parameter (l = 1) yields good improvement for small systems, but eventually the fidelity goes to zero for large n. By increasing the number of variational parameters, the improvement in F can even be of several orders of magnitude for small systems. However, the general trend is that, for large n, this improvement still gives a fidelity very close to zero and the proposed variational ansatz seems to be inefficient in the thermodynamic limit.

III.3. Results for p = 3
For p = 3, the ground state is nondegenerate. In this specific case, we discovered an alternative ansatz, which we named CA, yielding strikingly large fidelities in the symmetric sector, almost independent of the system size. The CA is shown in Eq. (10). In Appendix B, we show that for p = 3 it takes the compact form having only three variational parameters.  We start this section showing the scaling of the fidelity as a function of the system size. Due to the first-order QPT, we expect that the fidelity scales as We ask whether the CD driving can change this scaling law or not. In Fig. 3, we show the fidelity F as a function of the system size for T = 1 and we compare the bare quantum annealing with the nested commutator ansatz (l = 1, 3 and 8) and the CA. We perform a fit of the results conjecturing the exponential behavior of Eq. (21) even in the presence of CD. The coefficients φ and γ are summarized in Table I, comparing the standard quantum annealing, the NC ansatz (for several orders l), and the CA.
Note that the exponent γ in the CA is three orders of magnitude smaller than both the unitary case and the nested commutator ansatz. Moreover, in the latter case, we observe that the fidelity grows with increasing order l only for small system sizes, whereas for larger systems the fidelity shows a maximum as a function of l and then decreases. In fact, the exponent γ for l = 8 is larger than that for l = 3.
To summarize, in the presence of a CD the scaling of F with n remains exponential, but the coefficient γ is reduced with respect to the bare case. Moreover, the CA yields an almost constant fidelity up to n = 100, and a large one (F > 1/2) up to n = 1000, providing a robust mechanism to counteract the exponentially vanishing spectral gap for macroscopic systems.
In Fig. 4, we show the time evolution of the ground state probability P gs (t) for T = 1, for n = 10 [panel (a)] and n = 20 [panel (b)] in the maximal spin subspace. The blue dashed line indicates the fidelity of bare quantum annealing with no CD. The lightest green line corresponds to a CD dynamics with l = 8. Darker lines are for all orders starting from l = 1 (the lighter is the line, the larger is l). The red dotdashed line represents the CA. Fig. 4 clearly shows that the nested commutator ansatz of Eq. (9) can improve the fidelity of quantum annealing for p = 3. However, the CA yields an even larger fidelity. In comparison, a similar fidelity could be reached only going beyond order l = 8. However, increasing the number of NC the improvement in the fidelity is gradually smaller, and we guess that in order to achieve results similar to the one obtained using the CA an unpractical large number of nested commutator would be required. Moreover, as n grows, more and more variational parameters are required to achieve a similar level of fidelity, whereas the CA requires only three variational parameters. Hence, in this particular case, the CA is extremely efficient and outperforms the other known approximation schemes. This is even more evident increasing the number of qubits. For instance, Fig. 4(b) shows the same results for n = 20. Here, the fidelity for l = 8 is significantly smaller than the previous case (F ≈ 0.20 versus F ≈ 0.68). By contrast, the fidelity of the CA is almost unchanged with respect to n = 10, while the ground state probability at intermediate times is still affected by the system size.

IV. FINITE RANGE AND RANDOM INSTANCES
The efficiency of the CA could depend on the peculiarities of the p-spin model, i. e., spin symmetry and infinite-range interactions. However, it is difficult to prove this statement theoretically, therefore in the following we will try to limit the range of the interactions and to break spin symmetry to gain some insights into this problem. In the absence of spin symmetry, we need to extend the analysis presented in Sec. III.3 to the whole Hilbert space. Of course, this process is exponentially more demanding, since now we have to consider all the 2 n computational basis states. Therefore, in this section, we will limit our analysis to small systems, i. e., n = 3 to 8.
The total Hilbert space can be decomposed as the direct sum of the eigenspaces of the maximum total spin operator, corresponding to S = n/2, and of its orthogonal: In principle, all traces in Eq. (A6) have to be evaluated over H . However, we have already discussed that the variational procedure can be easily restricted to the interesting subspace. Using a parameter 0 ≤ η ≤ 1, we can choose whether to evaluate the traces in the whole Hilbert space or rather in the symmetric subspace. For any operator O, we replace The case η = 0 is analogous to Section III.3. Minimizing in the whole Hilbert space corresponds to choosing η = 1/2. The comparison between the two approaches is discussed in Appendix C. Finally, η = 1 is the minimization in the orthogonal subspace H ⊥ n/2 . Here focus our attention to the cases in which η = 0 and η = 1/2.

IV.1. Finite-range p-spin model
To highlight the infinite-range nature of the p-spin model, it is more convenient to rewrite its Hamiltonian [Eq. (11)] in the following form: where i j = 1, . . . , n for all j. As {σ z i , σ z j } = 2δ i, j , the Hamiltonian in Eq. (24) is a polynomial function of order p of Pauli operators, containing terms of orders P = p, p − 2, . . . , 1 (odd p). Each term represents an infinite-range P-body interaction between qubits, with uniform coupling constant J.
A possible way for turning this infinite-range p-spin model into a finite-range model is by weighting J with the "distance" between the qubits involved in the p-body term. This can be easily understood by considering the simple case p = 2, where the Hamiltonian would be In this case, we can replace J → J i, j = J/|i − j | ν and build a finite-range version of the p-spin model, where the exponent ν determines how punctual the interactions between the qubits are: ν = 0 is the infinite-range model and ν → ∞ represents nearest-neighbor interacting qubits. The same reasoning holds for any value of p. In particular, we can always replace J by J i 1 ,...,i p = J/dist(i 1 , . . . , i p ) ν . Here, we propose to consider the following form for the distance function: otherwise. (26) The parameter Z = (p 3 − p)/6 is a normalization factor chosen so that dist(1, 2, . . . , p) = 1. For p = 3, that normalization factor is Z = 4. We also note that this choice of the distance function does not allow finite-range models for n = 3 with p = 3, as in that case J i 1 ,i 2 ,i 3 = 1 for all combinations of indices. Of course, for all other values of n, this procedure breaks the spin symmetry of the p-spin model, therefore we will work in the whole Hilbert space and consider η = 1/2.
In Fig. 5, we show the ground state probability P gs (t) as a function of time, for n = 8 and p = 3. The green solid line corresponds to NC while the red dotdashed line corresponds to the CA. The left-hand (right-hand) panel refers to ν = 1 (ν = 10). This figure can be compared with the right-hand panel of Fig. 8, corresponding to ν = 0. Moving from the infinite-range model to the finite-range one, we note that the efficiency of both ansätze, NC and CA, is improved, as both curves are pushed upwards. In particular, for ν = 10, the fidelity F in the CA case is F ≈ 0.995, comparable with the nested commutator ansatz with l = 4. Can we conclude that the reason why the CA works so well cannot be the fact that the model is infinite-range as it works even better without this feature? It is difficult to compare to the p-spin model as 1) we are not working in the symmetric subspace and 2) the two models have different spectra (different gap, different time of minimal gap, see the unitary dynamics), however the results of this section motivated us to go even deeper and to analyze different (random) instances and check our ansatz also in that case.

IV.2. Random p-spin model
Previously, we showed that the CA can be efficient also for models featuring finite-range interactions, as in that case we can even get larger fidelities than those of the infinite-range ferromagnetic p-spin model. In this section, we will address another question. Starting from the Hamiltonian of Eq. (24), here we randomly suppress some of the coupling constants J with a certain probability. In this way, we can build a family of infinite-range models, where the full-connectivity of the p-spin model is progressively lost.
The resulting Hamiltonian is identical to that in Eq. (24), but the coupling constant J is replaced by a random variable K satisfying K = J with probability P J ; 0 with probability 1 − P J .
This model is the usual infinite-range ferromagnetic p-spin model when P J = 1. For any P J 0, 1, this model breaks the spin symmetry and we have to work in the whole Hilbert space, with η = 1/2.
For several choices of P J , we performed dynamical simulations for M randomly generated instances of this infinite-range random p-spin model and measured the fidelity, both with and without counterdiabatic terms. We focus here on the case n = 5 and p = 3, however we obtained qualitatively similar results also for larger system sizes.
In Fig. 6, we show the fidelity distributions for P J = 0.1, 0.3, 0.5, 0.7 and 0.9. We divided the fidelity interval [0.30, 0.45] into N b = 100 bins and counted the occurrences over M = 900 repetitions of the dynamics.
The peaks of the distributions are equally spaced, which implies that the mean fidelity F linearly depends on P J . According to the Landau-Zener formula, in a two-level approximation around the avoided crossing the mean fidelity would be As a consequence F ≈ ∆ 2 ∝ P J which implies that ∆ ∼ P 1/2 J . As P J < 1, all randomly generated instances have smaller gaps than the infinite-range ferromagnetic p-spin model with P J = 1, therefore the corresponding fidelity is always smaller than that of the original model in the absence of counterdiabatic terms. This is shown in Fig. 6, using a black dashed line to highlight the fidelity for P J = 1.
In Fig. 7, we show the same fidelity distributions in the presence of counterdiabatic driving. The left-hand panel is for the NC ansatz of order l = 3. The right-hand panel is for the CA. Here, we consider N b = 100 bins for the fidelity interval [0, 1]. Except for a few minor differences, such as the presence of outliers around F = 0 for P J = 0.1, the two plots look similar. However, we note that the mean fidelity in the cyclic case is smaller than that in the NC case for P J ≤ 1/2, while the opposite is true for P J > 1/2. In both cases, the presence of counterdiabatic driving allows for a significantly As opposed to the case with no CD driving, here there are some instances showing larger fidelity than that for P J = 1. This is highlighted in Fig. 7, where the black dashed line indicate the fidelity of the case P J = 1. This evidence confirms that the efficiency of counterdiabatic driving does not entirely depend on the spectral properties of the analyzed model. In fact, even if the average minimal gap is smaller than that for P J = 1, there are instances of the P J < 1 case where quantum annealing with CD is more efficient than the case P J = 1 with CD.

V. CONCLUSIONS
To summarize, we have applied the variational approach (developed in Ref. [12]) to find an approximate CD potential for the quantum annealing of the p-spin model. This approach is very promising [34] and has been recently applied to other kinds of model Hamiltonians with p-body interactions [13]. Our results confirm these expectations. The p-spin model possesses a spin symmetry, hence we first restricted our approach to the maximum spin subspace, where the starting and final states lay. We used two different ansätze for the variational potential: one is based on the NC hypothesis [14], and another based on a cyclic ansatz [see Eq. (10)]. We focused our attention on the case p = 3, which is the first odd integer showing a first-order quantum phase transition with exponentially small gap closure in the thermodynamic limit. In this case, both ansätze give substantial improvement in the fidelity with respect to the bare dynamics (i. e., standard quantum annealing). However, the CA seems to be much more efficient, in particular when the number of qubits grows. For both ansätze, we can also perform the variational optimization in the whole Hilbert space. However, in this case, the efficiency of the CA is reduced, as discussed in Appendix C. This is likely due to the fact that the corresponding optimal counterdiabatic potential also addresses diabatic transitions between pairs of energy eigenstates that are already uncoupled by the spin symmetry. This redundancy harms the performances of the CA in the relevant symmetry subspace, where the dynamics occurs. By contrast, restricting the traces to the subspace with maximal spin, we can readily implement the spin symmetry of the p-spin model and obtain better performances. We have also shown new results concerning modified model Hamiltonians accounting for short-range interactions and random instances. These results may be relevant for the future implementation of experimentally viable counterdiabatic operators on the available quantum processors. original ansatz. The exact counterdiabatic driving breaks time reversal, thus an odd number of σ y operators is required at all orders. However, in this approach, it is a priori unknown which operators have to be included to have a sufficiently accurate description of the counterdiabatic potential.
In Ref. [14], the authors proposed the following approximate gauge potential: For ease of notation, we have defined Now the operator G λ of Eq. (7) reads thus the corresponding action is we can define ì α = (α 1 , . . . , α l ) and recast S l as a quadratic polynomial is diagonal We can introduce the matrix U that diagonalizes C, i. e., D = U T CU. Then, S l is rewritten as with ì α = U T ì α and ì B = U T ì B. The stationary point is The exact CD driving is recovered in the limit l → ∞. In fact, the matrix elements of A (l) λ in the energy basis are to be compared with the exact gauge potential: The proposed variational ansatz is a power-series approximation of the factor ( ω mr ) −1 ≡ ( m − r ) −1 and will generally fail near ω mr = 0. Compared with the "heuristic" approach presented above, this NC ansatz has the obvious advantage that it can be improved by including higher-order terms and more variational parameters. The downside is that operators O 2k−1 are nonlocal and as difficult to implement as the original proposal by Berry [7]. As opposed to the exact CD potential, however, these operators are well-defined around the quantum critical point.
In this manuscript, we apply the counterdiabatic scheme to the ferromagnetic p-spin model. In the case p = 1, we only need l = 1 (l is the order of nested commutators) to recover the exact CD potential. For l = 1, the quadratic action S 1 (α) In the same way, it is possible to prove that Tr O 2 2 = 64n(n + 1)(n + 2) 12 1 − 2λ + 2λ 2 , so that independently of the system size.

Appendix B: Cyclic Ansatz
In Eq. (10), we propose a general form for the CA. In this section, we will work out explicitly all the relevant terms in the case p = 3. Starting from Eq. (10), we first explicitly evaluate all the terms: Then, using the angular momentum commutation rules, we can simplify the previous expression to A CA λ = α 1 S y + α 3 S 3 y + α (S x S y S z + h. c.), where we have defined α = α xyz = −α zyx = −α yxz = α zxy = −α xzy = α yzx and we have used the fact that (S x S y S z + h. c.) = (S z S x S y + h. c.) = (S x S z S y + h. c.). With this assumption, the CA requires just three variational parameters. In this appendix, we compare the efficiency of our algorithms both for η = 0 and η = 1/2, i. e., calculating the traces in the maximum spin subspace or in the whole Hilbert space, respectively. In Fig. 8, we show the instantaneous ground state probability as a function of time, for T = 1 and n = 8. The green line is for the NC ansatz of orders l = 8 and the red dotdashed line is the CA. The left (right) panel refers to the case η = 0 (η = 1/2).
Comparing the two panels, we notice that the ground state probability for the nested commutator ansatz is almost independent of η. By contrast, the CA outcome is strongly affected by η. In fact, for η = 0, the fidelity is very close to F = 1, while, for η = 1/2, the fidelity drops to F ≈ 0.7. Also the instantaneous value of P gs (t) has a different behavior depending on η. In fact, for η = 0, the fidelity drops before the avoided crossing and then immediately recovers following a nonmonotonic behavior, with a maximum at the avoided crossing (see also Fig. 4). By contrast, for η = 1/2, the ground state probability shows a minimum after the avoided crossing and then grows, with no maximum. Similar results also hold for other values of the system size n.