Rapid counter-diabatic sweeps in lattice gauge adiabatic quantum computing

We present a coherent counter-diabatic quantum protocol to prepare ground states in the lattice gauge mapping of all-to-all Ising models (LHZ) with considerably enhanced final ground state fidelity compared to a quantum annealing protocol. We make use of a variational method to find approximate counter-diabatic Hamiltonians that has recently been introduced by Sels and Polkovnikov [Proc. Natl. Acad. Sci. 114, 3909 (2017)]. The resulting additional terms in our protocol are time-dependent local on-site y-magnetic fields. These additional Hamiltonian terms do not increase the minimal energy gap, but instead rely on coherent phase maximization. A single free parameter is introduced which is optimized via classical updates. The protocol consists only of local and nearest- neighbor terms which makes it attractive for implementations in near term experiments.


INTRODUCTION
A fundamental limitation in adiabatic quantum computation (AQC) is subject to the adiabatic theorem [1][2][3] which states that a physical system follows its instantaneous eigenstate if the rate of change of a timedependent Hamiltonian is much smaller than the energy gap between its lowest eigenstates. This inevitable poses a speed limit for any algorithm based on AQC, such as solving combinatorial optimization problems by quantum annealing [4][5][6][7][8][9][10][11][12][13][14][15]. As the minimal energy scales with the system size, the question whether a possible speedup in quantum annealing exists is thus still open.
However, for general spin glass models that solve optimization problems, counter-diabatic Hamiltonians contain complicated non-local k-body interactions that are hard to achieve in experiment. Recently, a variational principle to find approximate counter-diabatic protocols for arbitrary Hamiltonians has been introduced by Sels and Polkovnikov [61]. An open challenge is to make use * w.lechner@uibk.ac.at of these local counter-diabatic protocols in many-body systems with non-local all-to-all connectivity.
In this work, we present a counter-diabatic Hamiltonian consisting of single on-site local magnetic fields (σ x , σ y and σ z ) and 4-body interactions between neighboring qubits (σ z σ z σ z σ z ) to solve all-to-all connected combinatorial optimization problems. The scheme is based on the recently introduced encoding of optimization problems in a lattice gauge model (LHZ) [62] where the optimization problem is fully determined by local magnetic fields and problem-independent interactions among nearest-neighbor qubits. We present an approximate counter-diabatic protocol in LHZ employing additional local single-body magnetic fields (σ y ) and derive the analytic expression for the time-dependent protocol as a function of local properties of each spin. For small systems, we demonstrate numerically the implementation of rapid sweeps with large ground state fidelity and small excess energies compared to an adiabatic protocol. A free control parameter for further improvement of the efficiency in a hybrid quantum-classical iterative update is introduced. The efficiency varies as a smooth function of this free parameter even for short sweep times which allows for a simple variational optimization of the parameter.
We note that the improved fidelity does not stem from an increase of the minimal energy gap between ground and excited state and cannot be understood in incoherent quantum annealing or with path integral Monte Carlo methods. It is rather the result of full quantum phase coherence. Thus, the results further encourage the current efforts to build next generation quantum annealing experiments in the fully coherent regime.

QUANTUM COUNTER-DIABATIC ANNEALING
Quantum annealing aims at solving optimization problems which can be translated into Ising spin glasses with logical spins that can take the values ±1 [4,5]. Finding the minimum energy of the spin glass is thus equivalent to determining the solution of the optimiza-tion problem [63]. The problem is cast into the form i whereσ z i is the z-Pauli matrix for the i-th logical spin and N the total number of logical spins. The magnetic fieldsb i and interaction matrixJ ij between the two sites i and j fully parametrize the system. Starting in the ground state of a trivial initial state, for example H i = N i=1 h iσ x i , the ground state of the problem Hamiltonian H p and thus the solution of our optimization problem is obtained by sufficiently slowly transferring the trivial initial state into the ground state of the problem Hamiltonian via the protocol where λ(t) is a smoothly varying parameter in time with λ(t = 0) = 1 and λ(t = τ ) = 0 at the beginning and end of the sweep, respectively. If the running time τ is infinitely large, the quantum system remains in its initial instantaneous eigenstate for all times during the sweep.
To overcome this limitation from the adiabatic condition, the counter-diabatic expression as a means of fast transitionless driving has firstly been mentioned by Demirplak and Rice [27][28][29] and later by Berry [30]. The basic idea of counter-diabatic driving (CD) is to evolve the system as where H 0 is the Hamiltonian of interest, A λ the adiabatic gauge potential andλ a free control parameter which in general are dependent on time.
Let us highlight the idea of the counter-diabatic Hamiltonian H CD and the role of the adiabatic gauge potential A λ and control parameter λ, respectively. For further considerations, we will consider λ as a singlecomponent parameter, although in general it can be a multi-component vector λ. A state |ψ evolves under a time-dependent Hamiltonian H(t) ≡ H(λ(t)) as i ∂ t |ψ = H(λ(t)) |ψ and |ψ = U † (t) |ψ in a rotating frame with respect to the unitary transformation U † (t). The Hamiltonian in the rotating frame has the form Solving this equation for A λ results in high order k-body interactions which are not realistic in current experimental implementations. Recently, a variational principle method to schematically apply counter-diabatic terms to arbitrary Hamiltonians has been introduced by Sels and Polknovnikov [61]. Here, the exact adiabatic gauge potential is approximated with an appropriate ansatz A * λ . Solving Eq.(4) for the adiabatic gauge potential A * λ is equivalent to minimizing the Hilbert-Schmidt norm of the Hermitian operator with respect to A * λ where we seek for the minimum of the In turn, minimizing the operator distance is equivalent to minimizing the action associated with the approximate adiabatic gauge potential A * λ (see Ref. [61] and [64] for more details), that is, where δ denotes the partial derivative.

COUNTER-DIABATIC DRIVING IN THE LHZ ARCHITECTURE
In the recently introduced lattice gauge model (LHZ) [62], the physical qubits describe the relative configuration of two logical spins taking values 1 for parallel (i.e. ↑↑, ↓↓) and 0 for antiparallel (↑↓, ↓↑) alignment. The time-dependent Hamiltonian in LHZ can be written in the form 1 where σ z k is the z-Pauli matrix for the k-th physical qubit and all local fields h k , J k and constraints C l depend on some tuning parameter λ which in turn depends on time. The first two sums run over all K = N (N − 1)/2 physical qubits where h k andJ ij → J k are controllable local fields that act on physical qubits. In the third sum, C l are 4body constraints constructed by closed loops of logical spins emerging due to the increased number of degrees of freedom from N logical to K = N (N − 1)/2 physical qubits. To account for this, K − N + 1 four-body constraints among nearest neighbors on a square lattice are introduced. The indices (l, n), (l, w), (l, s) and (l, e) denote the northern, western, southern and eastern physical qubit of the constraint l, respectively. The time-dependent protocols for h k , J k and C l are chosen such that h k vanishes at time t = τ with protocol function λ(t) and J k and C l vanish at time t = 0 with a protocol 1 − λ(t) where Here, τ is the sweep time and λ 0 and λ f the corresponding values for initial and final time, respectively. The time derivative of this protocol iṡ As the ansatz for the adiabatic gauge potential A λ of the LHZ Hamiltonian (8), we choose The additional local magnetic field (σ y ) is introduced for each physical qubit including auxiliary qubits in the lower row fixed. This ansatz is imaginary; thus it breaks instantaneous time-reversal symmetry and adds a new degree of freedom to the system. The operator (5) in LHZ reads where the prime stands for the derivative with respect to λ. Pauli matrices are traceless; therefore, we simply compute the Hilbert-Schmidt norm by adding up squares of coefficients in front of independent single spin terms of operator (12), that is, where 2 K is the dimension of the Hilbert space. The goal is to find the values of the local magnetic field strengths with minimal action in Eq.(6) corresponding to a minimum in operator distance between exact and approximate gauge potential. The optimal approximate solution for the adiabatic gauge potential A * λ is found by computing the derivative of the action with respect to α k and applying Eq. (7), to obtain where the sum in the denominator runs over all nearest neighbor constraints C k,n of the k-th physical qubit. 2 Note, that this solution for the adiabatic gauge potential is exact for any constraint strength C l equal to zero, as it is just a generator of local spin-rotations in the x-z plane. The potential A * λ also vanishes, if either h k = 0 or J k = 0 for all physical qubits implying that the leading contribution to A λ actually comes from the 4body interaction terms. For completeness, we can include 4-body interaction terms in our ansatz (see Appendix). The experimental implementation of the resulting terms is challenging and we will focus on the local solutions in this work.
The resulting local CD Hamiltonian in the lattice gauge model has the form where Equations (15) and (16), together with the variational optimization of the parameter λ f in Eq.(10) are the main result of this work. The complete protocol reads as follows: 1. Initial State: Prepare the ground state of the trivial driver Hamiltonian H i and set an appropriate sweep time τ .
2. CD sweep: The local fields h k (t), J k (t) and constraint strengths C l (t) are driven according to protocol (9) with initial and final values h k (0) = 1, h k (τ ) = 0, J k (0) = 0, J k (τ ) = 1, C l (0) = 0 and C l (τ ) = 1 for k, l in the set of all physical qubits and constraints, respectively. Implement the protocol of the magnetic field strength in y-direction as in Eq. (16).
3. λ f values: The parameter λ f corresponds to a global factor for the magnetic field strength and can be optimized as a variational parameter from iterating the sweep in step 2 and minimization of the final energy.

RESULTS
In the first part of this section, we discuss the protocol for a single random instance of J k in Hamiltonian (15). In the second part, we present the statistics from an ensemble of 100 randomly chosen instances.

Single instance
Now, let us first describe the results for the local counter-diabatic Hamiltonian (15) using an example with N = 4 logical and thus K = 6 physical qubits plus 2 fixed auxillary qubits in LHZ. The drive fields are bounded to λ(t) ∈ [−10 J, 10 J]. A measure for the efficiency of the model is the probability to find the system in its ground state for different times during a counter-diabatic sweep. We consider the instantaneous ground state fidelity squared F 2 (t) = | ψ(t)|φ 0 | 2 , where φ 0 is the instantaneous ground state and ψ(t) the state of the system at time t. Figure 1 depicts the instantaneous ground state fidelity squared F 2 (t) of a single instance of Hamiltonian (15) with randomly uniformly distributed values between −1 and +1 during a whole counter-diabatic sweep with running time τ = 1 /J. For intermediate times during the sweep, the instantaneous ground state fidelity drops rapidly and then increases to a value of around t/τ = 0.28, whereas for the naive annealing protocol (9) rapidly decreases and then stays at a value of around 0.01. The inset in Figure 1 depicts the free control parameter λ f in protocol (10) which can dramatically enhance the performance of the counter-diabatic Hamiltonian (15). The distribution of the final ground state fidelity squared F 2 (τ ) for multiple values of λ f close to its optimal value of 1.04 J is shown in the inset of Figure 1 makes the local CD Hamiltonian (15) stable against perturbations and allows for an experimental implementation of the iterative variational update. Figure 2 depicts the values in front of σ y k for each physical qubit during this counter-diabatic sweep with the same parameters as described above. The terms σ y 6 and σ y 7 correspond to the auxiliary physical qubits in the LHZ architecture. The strengths of the two auxiliary local  (8). The minimal energy gap ∆Emin between ground and first excited state shifts from around t ≈ 0.46τ for the naive annealing to t ≈ 0.56τ for the counter-diabatic Hamiltonian. Note, that even though the fidelity increases considrably, the minimal gap of the counterdiabatic protocol is smaller compared to that of the annealing protocol.
magnetic fields are identical as their final value is fixed for both to 10 J.
In adiabatic protocols, the minimal energy gap is the fundamental limitation for the sweep time. To analyse the influence of the minimal energy gap on the CD Hamiltonian (15) during the sweep, we compare the energy spectra of the counter-diabatic Hamiltonian (15) and the annealing Hamiltonian (8) using the parameters described above (see Figure 3). The minimal energy gap ∆E min = E 1 − E 0 between ground and first excited state shifts from around t ≈ 0.46τ for the naive annealing Hamiltonian (8) to t ≈ 0.56τ for the counter-diabatic Hamiltonian (15) while the minimal gap even decreased.

Statistical Ensemble
Let us now examine the counter-diabatic sweeps for an ensemble of 100 randomly chosen instances J k . Figure 4 depicts the averaged final ground state fidelities squared and excess energies ∆E = E − E 0 , respectively, where averages are taken from fixed protocols τ and uniformly distributed J k instances in (15).
In the quench limit τ → 0, a final ground state fidelity squared for the counter-diabatic Hamiltonian (15) is 0.59. In the annealing protocol (8)  factor of 150. On the other hand, the corresponding excess energies of the counter-diabatic Hamiltonian is 2.9 J; whereas for the naive annealing Hamiltonian they are approximately 28 J. For long running times, the final ground state fidelities squared and excess energies for both, counter-diabatic (16) and naive annealing (9) protocols, converge towards the same value as the amplitude of the added y-magnetic field becomes negligible compared to the annealing problem Hamiltonian for long running times due the inversely proportional dependence ofλ(t) on τ .
We note, that there is an experimental limitation in implementing counter-diabatic protocols. This is a trade-  Figure 5. Energy scaling. The energy scaling of the counter-diabatic term (16) and protocol (9) for different running times is shown.
off between the obtained increase in speed and energetic cost of our implemented CD Hamiltonian (15) and its feasible applicability in the experiment. Figure 5 depicts the energy scaling of the counter-diabatic Y k (t) term (16) and derivative of the executed protocol (10) for different fast protocols τ . For different sweep times τ , the control parameterλ(t) scales with the factor 1/τ . Considering the counter-diabatic protocol (16), enhancing the performance of the counter-diabatic Hamiltonian (15) by one order of magnitude corresponds to increasing the energy of the eigenstates during intermediate times by a factor of 100.

CONCLUSION AND OUTLOOK
We have introduced an approximate optimal counterdiabatic driving protocol for the LHZ lattice gauge model architecture from a variational principle. Using an experimentally accessible local ansatz for the adiabatic gauge potential A λ , we derived a counter-diabatic Hamiltonian that consists of local fields only and enables counterdiabatic Hamiltonian driving even for all to-all connected spin glass problems. The counter-diabatic term Y k (λ, t) added in Eq.(15) is a local magnetic field in y-direction which depends on some free tuning parameter λ f . This enables a hybrid classical-quantum algorithm where λ f is updated from measurements after the quantum process.
Furthermore, we demonstrated a large increase in final ground state fidelity and decrease in excess energy with the counter-diabatic Hamiltonian (15). The CD driving keeps the system close to its ground state and dramatically enhances the performance of quantum annealing protocols to solve optimization problems using the lattice gauge model for Ising spins (LHZ). The increase in final ground state fidelity and decrease in excess energy, respectively, do not emerge due to an increase in the minimal energy gap, thus does not just follow adiabatic theorem and Landau-Zener's formula. Instead, the position of the minimal energy gap shifts such that we optimize the total phase θ for the eigenstates of our counter-diabatic Hamiltonian H CD,LHZ . Therefore, we propose finding an optimal approximate gauge potential and CD protocol in the mapping of lattice gauge model as a means of coherent phase maximization. We anticipate that due to coherent evolution we can extend the number of logical and physical spins as in this example with just K = 6 plus two auxiliary physical qubits.
Remarkably, the ratio between final ground state fidelity of the counter-diabatic protocols and the annealing protocols increases with N (see Appendix for the example with N = 3 physical qubits for comparison). While for N = 3 we see an improvement by one order of magnitude it increases to a relative improvement of two orders of magnitude for N = 4. This is an encouraging result which we will study in detail if this trend continues for larger systems.
The energies of the σ y terms in Eq.(16) scale with 1/τ 2 as shown in Figure 5. Therefore, we expect an experimental limit for these counter-diabatic terms. The relevant and accessible regime is the one where the energy scale of the σ y term is comparable to that of the σ x term. In this regime, i.e. τ ≈ 10 −1 /J, the resulting speedup for solving optimization problems with counter-diabatic terms compared to annealing is three orders of magnitude in the example above.
As a future direction, the counter-diabatic Hamiltonian (15) of the LHZ architecture may be applied to quantum approximate optimization algorithms (QAOA) where the system is sequentially quenched with unitaries, that is, we can combine the speedup of the counterdiabatic protocol in the quench limit τ → 0 with the unitary quenches in QAOA [65,66] which may result in improved efficiency of this method.
The CD Hamiltonian (15) is non-stoquastic [67] and thus cannot be effectively solved with classical algorithms such as path integral Monte Carlo methods or simulated in stoquastic quantum annealing devices. We thus hope that our work contributes as a possible application to the current efforts in building next generation quantum annealing experiments with full quantum coherence.
Here we describe in detail the derivation of the adiabatic gauge potential in Eq. (15). Evolving a state |ψ according to the Schrödinger equation i ∂ t |ψ = H(λ(t)) |ψ with a time-dependent Hamiltonian H(λ(t)) ≡ H(t) in a rotating frame |ψ = U † |ψ leads tõ where the adiabatic gauge potential reads Going back to the laboratory frame, that is removing the tildes, and the fact that the gauge potential eliminates the off-diagonal terms of the moving Hamiltonian, i.e.
[∂ λH ,H] = 0, we get where the first element in the commutator is precisely the operator G(A λ ) (5) in the maintext with = 1.
Using the ansatz A * λ = i α i σ y i and computing the operator G(A * λ ), leads to the commutator 2C l (α l,n σ x l,n σ z l,w σ z l,s σ z l,e + α l,w σ z l,n σ x l,w σ z l,s σ z l,e + α l,s σ z l,n σ z l,w σ x l,s σ z l,e + α l,e σ z l,n σ z l,w σ z l,s σ x l,e ).

Control parameter protocol
The function λ(t) as in Eq. (9) in the maintext is vanishing in first and second order derivatives at the beginning and end of the sweep, respectively, i.e.λ(t = 0) = λ(t = 0) =λ(t = τ ) =λ(t = τ ) = 0. In the philosophy of counter-diabatic driving, this protocol (9) behaves like an adiabatic protocol in the beginning and end of the sweep and accelerates during intermediate times. Figure  6 depicts the protocol λ(t) and its derivativeλ(t) for the case of λ 0 = 0 and λ f = 1.04 J, respectively.

Fidelity Distribution
In the maintext, we have seen that for the quench limit τ → 0, the mean final ground state fidelity squared for the case of K = 6 physical qubits is around 0.59. Figure 7 shows the distribution of the reached maximal final ground state fidelities squared for all 100 J k instances. The distribution is roughly Gaussian with most of the instances having a final ground state fidelity squared between 0.50 and 0.60.

N=3
With the aim to compare different system sizes, we consider the minimal example for LHZ with N = 3 logical qubits. Figure 8 depicts the final ground state fidelities squared and excess energies for 100 randomly uniformly chosen J k instances for an Ising spin glass in a LHZ lattice gauge model with K = 3 physical qubits and one constraint C 1 where we have added one auxiliary qubit in the lower row to obtain a 4-body constraint. In the quench limit τ → 0, we achieve a final ground state fidelity squared of around 0.93 for the counter-diabatic Hamiltonian (15) in the maintext and 1/2 4 = 0.0625 for the annealing Hamiltonian (8), respectively, which gives an enhancement of a factor of around 15. The excess energy of the counter-diabatic Hamiltonian on the other hand is around 0.4 J; whereas for the annealing case it stays at around 13.2 J which gives an enhancement (a) (b) Figure 8. Ground state fidelities. The CD Hamiltonian (15) in the LHZ architecture with K = 3 physical qubits and parameters C = −2 J, h k = 1 J and randomly chosen J k interaction matrices over 100 instances undergoes a counter-diabatic sweep. (a) shows the final ground state fidelity squared during different fast protocols τ ; (b) shows the excess energy during different fast protocols. The blue circles are associated with the counter-diabatic Hamiltonian (15) and the red circles with the naive annealing Hamiltonian (8). For both plots we have optimized the bounded parameters λ f ∈ [−10 J, 10 J]. of a factor of around 33. Even in the quench limit, the probability to prepare the ground state of the CD Hamiltonian (15) is finite.
The enhancement of the counter-diabatic protocol compared to annealing for different sweep times is measured via the ground state fidelity squared as in the main text. Figure 9 depicts the ratio between counter-diabatic and naive annealing final ground state fidelities squared.
For the quench limit τ → 0, the ratio F 2 C (τ )/F 2 n (τ ) between counter-diabatic and naive annealing ground state fidelities squared is around 15 and decreases to a value of 1, whereas the excess energies for the counter-diabatic Hamiltonian (15) are just around 0.03 of the naive annealing Hamiltonian (8) and increase to a value of 1.  Figure 9.
Relative improvement of the counterdiabatic protocol. An ensemble of 100 J k instances undergo counter-diabatic sweeps for different sweep times τ . The ratio between the final ground state fidelity squared F 2 (τ ) of counter-diabatic (15) and naive annealing case (8) is shown. The blue line depicts the ratio in fidelity and the orange in excess energy.
Again like in first order, minimizing the action leads to the optimal solution. The derivative of T r[G 2 (A λ )] with respect to all parameters α k , β l , γ l , δ l and l and solving the linear equation system leads to the optimal, yet very unhandy, optimal solution and counter-diabatic Hamiltonian in 2nd order.