Non-Stoquastic Interactions in Quantum Annealing via the Aharonov-Anandan Phase

We argue that a complete description of quantum annealing (QA) implemented with continuous variables must take into account the non-adiabatic Aharonov-Anandan geometric phase that arises when the system Hamiltonian changes during the anneal. We show that this geometric effect leads to the appearance of non-stoquastic terms in the effective quantum Ising Hamiltonians that are typically used to describe QA with flux-qubits. We explicitly demonstrate the effect of these geometric interactions when QA is performed with a system of one and two coupled flux-qubits. The realization of non-stoquastic Hamiltonians has important implications from a computational complexity perspective, since it is believed that in many cases QA with stoquastic Hamiltonians can be efficiently simulated via classical algorithms such as Quantum Monte Carlo. It is well-known that the direct implementation of non-stoquastic interactions with flux-qubits is particularly challenging. Our results suggest an alternative path for the implementation of non-stoquastic interactions via geometric phases that can be exploited for computational purposes.

Introduction.-It is well known that the solution of computational problems can be encoded into the ground state of a time-dependent quantum Hamiltonian. This approach is known as adiabatic quantum computation (AQC) [1][2][3], and is universal for quantum computing [4] (for a review of AQC see Ref. [5]). Quantum annealing (QA) is a framework that incorporates algorithms [6][7][8] and hardware [9][10][11][12][13] designed to solve computational problems via quantum evolution towards the ground states of final Hamiltonians that encode classical optimization problems, without necessarily insisting on universality or adiabaticity.
QA thus inhabits a regime that is intermediate between the idealized assumptions of universal AQC and unavoidable experimental compromises. Perhaps the most significant of these compromises has been the design of stoquastic quantum annealers. A Hamiltonian H is stoquastic with respect to a given basis if H has only real nonpositive offdiagonal matrix elements in that basis, which means that its ground state can be expressed as a classical probability distribution [14,15]. Typically, one chooses the computational basis, i.e., the basis in which the final Hamiltonian is diagonal. The computational power of stoquastic Hamiltonians has been carefully scrutinized, and is suspected to be limited in the ground-state AQC setting [5]. E.g., it is unlikely that ground-state stoquastic AQC is universal [16]. Moreover, under various assumptions ground-state stoquastic AQC can be efficiently simulated by classical algorithms such as quantum Monte Carlo [14,15,17,18], though certain exceptions are known [19,20].
One is thus naturally motivated to consider a departure from the stoquastic setting. Indeed, this is the setting of proofs of the universality of AQC and of various specific results that use non-stoquasticity to improve upon the performance of a stoquastic Hamiltonian [21][22][23][24][25][26]. For example, it is known that non-stoquastic interactions that are turned on temporarily during QA can modify the annealing path so that it encounters a polynomially small gap rather than an exponentially small one, by replacing a first order quantum phase transition by a second order one [27].
Introducing non-stoquastic interactions is especially important in the physical implementation of QA devices, in order to allow them to escape the trap of efficient classical simulability. As a case in point, and setting aside heavily debated concerns about whether such devices are sufficiently quantum [28][29][30][31][32][33], despite intense efforts [32,[34][35][36][37][38][39][40] there is currently no evidence of an example where stoquastic QA hardware such as the D-Wave devices [10][11][12] delivers a quantum speedup over all possible classical algorithms. This is true even though in this setting QA is not limited to ground state evolution. However, the implementation of non-stoquastic interactions is technologically challenging, at least with superconducting flux-qubits. The D-Wave devices, for example, have a scalable design that can only implement the Hamiltonian of the quantum transverse field Ising model, which is stoquastic. To remedy this would require additional couplings between the flux qubits, to realize, e.g., σ x ⊗ σ x interactions, in addition to the existing σ z ⊗ σ z interactions. This greatly complicates the design of the current generation of superconducting circuits.
Rather than attempt to introduce non-stoquasticity via new components, here we revisit the assumptions that lead to the derivation of the Hamiltonian generating the effective time evolution of a continuous-variable system, such as inductively coupled flux qubits. We show that non-stoquastic terms arise naturally as a non-adiabatic geometric phase, due to the Aharonov-Anandan effect [41,42]. In other words, nonstoquastic terms are in fact present all along when QA is performed in systems of inductively coupled flux qubits. We study these geometric terms in detail, and show that the geometric effect is amplified in proportion to the inverse gap of the flux Hamiltonian, appearing and then disappearing during the anneal. The geometric effect can thus be considered as a type of non-stoquastic catalyst [5], with the potential to lead to quantum speedups [21][22][23][24][25][26][27]. The presence of geometric terms also has clear experimental consequences. Geometric effects should be taken into account in the validation of current and future QA devices. Conversely, the same devices could be used to perform experimental measurements of non-trivial geometric phases.
General formalism and the geometric term.-For concreteness, we assume that QA is performed by implementing a time-dependent Hamiltonian that can be generically expressed as follows: Such Hamiltonians can be realized with superconducting fluxqubits [43], in which case the continuous variables φ = {φ i } n i=1 are the magnetic fluxes trapped by the n flux-qubits, and E C represents a charging energy. The term P (φ, t) is a time-dependent potential that controls both the anneal and the interactions between qubits. We assume that the lowest 2 n energy levels of the Hamiltonian (1) are separated from the rest of the spectrum by an energy gap. At sufficiently low temperatures and for a slowly variable potential P (φ, t), the system of Eq. (1) is then effectively confined to an N = 2 n dimensional Hilbert space H N (t).
To proceed, we follow the approach introduced by Anandan in the discussion of non-adiabatic non-Abelian geometric phases [42]. We first consider the time-evolved Ndimensional subspace H N (t), which is defined by the map where dt is the unitary timeevolution operator (T denotes time ordering), and H(t) is the Hamiltonian (1). Let {|ψ a (0) } N a=1 denote an orthonormal basis for H N (0). Thus, H N (t) has a basis whose orthonormal elements obey the Schrödinger equation: i∂ t |ψ a (t) = H(t)|ψ a (t) . ( We now define another (arbitrary) orthonormal basis |ψ a (t) for H N (t), such that |ψ a (0) = |ψ a (0) . Then there exists an N × N unitary transformation W (t) between the two bases such that |ψ b (t) = N a=1 W ab (t)|ψ a (t) . An arbitrary state |ω(0) = a ω a (0)|ψ a (0) ∈ H N (0) thus evolves according to: can thus be considered as describing the effective evolution of the state |ω(t) inside the subspace H N (t) in the frame rotating with the basis |ψ a (t) (we derive the evolution equation of this basis in the SM, section A). By substituting the basis transformation into Eq. (3) one easily finds: whereH(t) and G(t) are the N × N matrices defined by the following matrix elements, respectively: Let us now assume that H(t) depends on t only via the invertible and differentiable "annealing schedule" s ≡ κ −1 (τ ), where τ ≡ t/t f , and t f denotes the final time. Then it follows directly from Eq. (5b) that G(t)dt = G(s)ds (see the SM, section B), which shows that G(s) is a geometric term, i.e., it depends only on the schedule s (and not on its parametrization) [44,45]. Consequently, W (t f ) = T exp −i 1 0 H eff (s)ds with the dimensionless effective Hamiltonian where from hereon a dot denotes d/ds, and we set s ≡ τ for simplicity. Note that the geometric term is negligible only in the adiabatic limit t f → ∞. As we shall see, is responsible for the appearance of non-stoquastic terms when t f is finite. Application to QA with superconducting flux qubits.-Let {|a(s) } denote the N lowest energy instantaneous eigenstates of the original Hamiltonian (1), i.e., H(s)|a(s) = E a (s)|a(s) . We identify the earlier {|ψ a (s) } with these eigenstates, so that W (s) describes the unitary evolution in the subspace rotating with the instantaneous eigenbasis of H(s). In this basis, the QA process is described by the dimensionless effective Hamiltonian (6) with: In QA applications, the Hamiltonian (1) is designed such that its N = 2 n lowest energy levels can be put into 1-to-1 correspondence with the energy levels of a transverse field Ising model with n qubits. Thus, there exists a unitary V (s) such that, up to a term proportional to the identity matrix [32]: where the profile functions A(s) and B(s) are particular to the qubit Hamiltonian (see SM, section C),H X = − n i=1 σ x i is the usual transverse-field driver, and is the problem Hamiltonian whose ground state encodes the answer to the optimization problem of interest. The unitary V (s) is the transformation between the energy eigenbasis and the computational basis |ψ C a (s) . Note that the geometric term transforms as a geometric connection [42] (see SM, section D): Quantum annealing of a system of n flux-qubits controlled by the Hamiltonian (1) is then described by the following effective Hamiltonian in the computational basis: The first term, proportional to t f , is the usual Hamiltonian discussed in literature in the context of QA. The second term has a geometric origin and is non-vanishing for finite annealing times t f . The geometric term is non-stoquastic. To see this, note that the original Hamiltonian (1) is real, and thus there is a basis choice in which the energy eigenbasis states |a(s) have only real amplitudes. Consequently, the geometric term G(s) in Eq. (7) is then a purely imaginary Hermitian matrix. The transformed geometric term G C (s) [Eq. (10)] is also purely imaginary in the computational basis, since the basis change of Eq. (8) can be performed with a real unitary (i.e., orthogonal) matrix V (s). Therefore, including only interactions up to two-body terms, the geometric term can be written in the most general form as follows: (12) Quantum annealing with geometric terms: one qubit.-We now apply the results obtained so far to QA with one qubit. For concreteness we focus on the C-shunt flux-qubit with three Josephson junctions. This qubit can be described by the following Hamiltonian [46][47][48] (see SM, section E): where φ is the flux (in units of Φ 0 /2π) trapped in the superconducting ring, E S is the charging energy of the shunting capacitor, and φ x CJJ (s) and φ x (s) are external fluxes used to control the anneal.
We next construct the effective Hamiltonian of Eq. (11). We start by numerically computing the ground state |0(s) and first excited state |1(s) of the flux-qubit Hamiltonian Eq. (13) [examples are given in the SM, Figs. 5(b)-5(d)], from which we numerically obtain the effective Hamiltonian [Eq. (6)] in the instantaneous energy eigenbasis: (we have ignored a term proportional to the identity ma- Fig. 1(a)] and g(s) ≡ 1(s)|∂ s |0(s) . We now transform to the computational basis using the unitary V (s) = exp i 2 arctan A(s) B(s) σ y , after which the Hamiltonian of Eq. (14) becomes This has the form of the most general single-qubit Hamiltonian: H eff,C 1 (s) = α∈x,y,z c α (s)σ α . It can be checked easily that (see SM, section F): By defining the computational basis as the basis of states with well-defined persistent current, the annealing schedule functions A(s) and B(s) [shown in Fig. 1(a)] can be determined in terms of the external fluxes φ x CJJ and φ x of Eq. (13b) (see the SM, section C). Thus, Eq. (16b) shows that these flux parameters in turn control the properties of the geometric term g y . The profile function g y (s) for the geometric term, also shown in Fig. 1(a), is non-vanishing towards the middle of the adiabatic evolution, when the gap closes. The effects of the geometric term are shown in Fig. 1(b), obtained by numerically solving for the corresponding unitary evolution. The effects become significant when the gap is small. There is also a significant effect on the final ground state population.
Quantum annealing with geometric terms: two qubits.-To study the interacting case, we consider two inductivelycoupled compound Josephson junction flux qubits [46] (see SM, section E): where the flux-qubit Hamiltonian is given by and the interaction potential is explicitly written as: Proceeding as in the single qubit case, we start from the Hamiltonian (17) to numerically computeH and G appearing in Eq. (7), and construct the effective Hamiltonian in the instantaneous energy eigenbasis [Eq. (6)]. Once again, the effective Hamiltonian can be expressed in the computational basis by (numerically) finding a unitary V such that the final Hamiltonian reads: This is the usual transverse field Ising model, plus an additional geometric term, whose general form is given in Eq. (12). Figure 2(a) shows the gaps ∆ k,0 = E k − E 0 , where {E k } 3 k=0 are the ordered eigenvalues of H eff,C

2
, and the components of the geometric term G C (s, h 1 , h 2 , J 12 ) in the computational basis as defined in Eq. (12). Fig. 2(a) shows that, in general, the geometric functions g y i (s), g xy ij (s), g zy ij (s) are all non-vanishing, and their magnitude grows as the ground state gap shrinks. In particular, as in the single-qubit case, geometric effects introduce a non-stoquastic contribution to the driver term which is non-vanishing towards the middle of the anneal. Fig. 2(b) shows the effect of the geometric interactions in the annealing of the system of two coupled qubits under consideration. As in the single qubit case, we see that ignoring the geometric terms results in consistently different final populations in the computational basis.
Dependence on the annealing time.-The contribution of the geometric terms is inherently a non-adiabatic effect, arising when diabatic transitions populate the excited states. In Fig. 3 we show the fidelity | ψ C 1,G (s)|ψ C 1,no-G (s) | 2 between the time-evolved states with and without geometric terms, as a function of the total annealing time. As expected, the effect of the geometric terms increases with decreasing annealing time and it is reflected in the decreasing fidelity. Figure 3 shows that the total annealing time over which this contribution is significant is on the order of a few nanoseconds, for parameters relevant for current QA devices. While this is significantly shorter than the typical microsecond timescale of current QA experiments using the D-Wave devices, this is the case for the one and two-qubit cases which we have analyzed here. Since, as is clear from Eq. (16b) and from Fig. 2, the magnitude of the geometric term grows in inverse proportion to the gap, we expect it to become more significant for multi-qubit problems whose gap dependence can be inverse polynomial or even exponential. This effect is already visible in Fig. 3, which shows that the fidelity for the two-qubit system tends to be smaller than the one-qubit system for larger annealing times.
Conclusions.-We have shown that even the simplest implementation of QA with flux-qubits induces effective Hamiltonians with non-stoquastic interactions arising from a geometric phase. The appearance of such interactions is ubiquitous when the Hamiltonian of a continuous-variable system is changed over time, and the evolution is non-adiabatic, due to the appearance of the Aharanov-Anandan effect. Since arbitrarily small gaps are inevitable in QA for hard optimization problems, non-adiabatic evolutions are ultimately inescapable. We thus argue that, similarly, the geometric effects studied here are unavoidable and relevant in practical applications of QA. Moreover, since these geometric effects give rise to non-stoquastic terms in the effective Hamiltonian, they provide a natural and desirable mechanism for avoiding classically efficient simulation of QA. This may point to the possibility of "quantum supremacy" experiments with QA devices featuring fewer than 100 physical qubits [49][50][51][52]. (A1) Differentiation yields: Thus, the evolution equation satisfied by the basis element |ψ a (t) is where we usedẆ = −iH eff (t)W (t) and unitarity, and de-finedH eff (t) ≡ W † (t)H eff (t)W (t).

Appendix B: Proof that G is a geometric term
Let us show explicitly that G is a geometric term, i.e., that G(t)dt = G(s)ds where s = κ −1 (t/t f ). Note that since we assumed that κ is invertible and differentiable, we can write t = t f κ(s), so that dt = t f dκ(s) where in the second line we used the assumption that H(t) depends on t only via s, which means, by Eq. (5a), that the same must be true of ψ a (t). In this section we closely follow Ref. [32] to briefly describe how the annealing profile functions A(s) and B(s) can be calculated. The main observation is that the control flux φ x , i.e., the bias between the two potential wells, is a small perturbation for the potential of H 1 in Eq. (18). The eigenstates |0(s) and |1(s) of the unperturbed Hamiltonian with φ x = 0 are used to define the states of the computational basis: These symmetric and antisymmetric combinations have opposite and well-defined circulating persistent current along the whole anneal, with the persistent current operator defined as I p = E L φ. This justifies the use of | ↑(s) and | ↓(s) as states of the computational basis.
We can now expand the flux qubit Hamiltonian H 1 in Eq. (18) to first order in φ x to get: Evaluating the Hamiltonian above in the basis (| ↑(s) , | ↓(s) ) then gives (up to a term proportional to identity) the Hamiltonian: where The profile functions above are completely controlled via the external fluxes φ CJJ (s) and φ x (s).
The schedule of the flux φ x (s) and thus of the profile function B(s) is further constrained in the case of multi-qubit interactions. The interaction potential is given by Eq. (19), whose expectation value in the computational basis is given by: from which it follows that the interaction term in the effective Hamiltonian is given by To ensure the same annealing schedule for the local and interaction terms, the control field φ x is then chosen to be proportional to the persistent current φ x (s) = E M E −2 L I p (s). This implies that B(s) = E M E −2 L I p (s) 2 . We have considered an annealing schedule linear in the control field φ CJJ (s) [see Fig. 4(a)]. The corresponding values for the control flux φ x (s) and profile functions A(s) and B(s) are shown in Fig. 4(b) and are computed using the prescription described in this section, using the numerical methods of the next section. Figure 6 shows the ratios between the exact gaps computed by numerical diagonalization of Eqs. (13) and (17), and the gaps computed by diagonalizing the Hamiltonians Eqs. (15) and (19) (without geometric terms), when the profile functions A(s) and B(s) are computed as described in this section. Due to the perturbative expansion, the method described in this section only approximately recovers the exact gaps. The relative error is largest (up to 2%) in the middle of the anneal, when the gaps close.

b. C-shunt Flux-Qubit
As in the CJJ case, we will treat φ x as a small perturbation. Expanding the Hamiltonian (13) to first order in φ x gives: The Hamiltonian above has the following representation: with For consistency with the CJJ example, we have taken φ x (s) ≡ E M E −2 L I p (s) 2 /I p (s) with E −1 M E 2 L = 10 4 GHz, such that B(s) is proportional to the square of the persistent current.
Appendix D: Proof that G transforms as a geometric connection Let us show that G transforms as in Eq. (10), i.e., that Recall that we transform from the basis |a(s) (e.g., the energy eigenbasis) to the new basis |ψ C a (s) (e.g., the computational basis) using the unitary V (s). Thus, . From now on we drop the explicit s-dependence for simplicity. By definition, Using these two expressions, we have: i.e., G C = V GV † −iV V † , which gives the desired result after using unitarity to writeV V † = −VV † .

Appendix E: Flux-Qubit Hamiltonians
The Langrangian of a superconducting circuit with Josephson junctions is generically written as follows [43]: where the first term is the kinetic term, the second is the junction energy, and the third is the induction energy. The φ c are the phase differences at each capacitance. The charging energy is E C,c ≡ (2e) 2 /C c , where C c is the c-th capacitance and e the electron charge. The junction energies are given by E J,j ≡ I c,j (Φ 0 /2π), where I c,j is the critical current of the j-th junction. The induction energies are given by E L,l ≡ (Φ 0 /2π) 2 /L l , where L l is the induction of the l-th loop. The conjugate momenta are p c = ∂L/∂φ c = E −1 C,cφ c . The circuit is quantized by promoting the momenta to opera-tors: p c → −i∂ φc . The corresponding Hamiltonian is:

c. Compound Josephson Junction (CJJ) Flux-Qubit
The basic design of a CJJ flux-qubit [46] is shown in Fig. 7(a). We assume the same charging and junction energies for the two Josephson junctions and a negligible inductance of the small loop. The Hamiltonian for this device can then be written as:  i.e., the sum of the charging, junction and induction energies of the circuit. By defining φ ≡ (φ L − φ R )/2 and φ CJJ = φ R + φ L we have ∂ φ L,R = ∂ φCJJ ± 1/2∂ φ , from which we obtain: where we have neglected the term − 1 2 (2E C )∂ 2 φCJJ since the small loop inductance gives φ x CJJ = φ CJJ , i.e., the flux φ CJJ is locked to the external flux φ x CJJ . The equation above reduces to Eq. (18) with the redefinition φ x CJJ → 2π − φ x CJJ .

d. Capacitively Shunted (C-shunt) Flux-Qubit
The basic design of a C-shunt flux-qubit [46,47] involves four Josephson junctions and a large shunting capacitance and is shown in Fig. 7(b). We start by writing the kinetic term: where the first term comes from the junctions while the last term is the shunting capacitor energy with E S = (2e) 2 /C S .
Combining the last two equations yields Eq. (16b) reported in the main text.

Appendix G: Numerical Methodology
All the "static" quantities, e.g., the quantities shown in Figs. 1(a), 2(a), 4, 5 and 6, are determined at a given value of the schedule parameter s by first numerically computing the wave functions |0(s) and |1(s) [see, e.g., Fig. 4 and Fig. 5].
We first discretized the flux-qubit Hamiltonian Eq. (1). For example, the one-qubit Hamiltonian Eq. (13a) is reduced to the following L-dimensional system: where the first term is the discretized Laplacian and the continuous flux φ is discretized as φ i = φ −L +i(φ L −φ −L )/(L− 1), i = 0, . . . , L − 1, with L being the size of the mesh. Similarly we can discretize the two-qubit Hamiltonian Eq. (17) to obtain an L 2 -dimensional system. We used L = 600, which was sufficient for numerical convergence. We numerically computed all the static functions on a mesh of 100 points for the annealing parameter s. Once all the static quantities where computed, the "dynamic" quantities of Figs. 1(b), 2(b) and 3 were computed by solving the Schrödinger equations resulting from the effective one-and two-qubit Hamiltonians Eqs. (14), (15) and (19). We used a standard ode45 solver provided with Matlab.