On the power of coherently controlled quantum adiabatic evolutions

It is not possible to provide a closed form solution to the Schrödinger equation for the case of time-dependent Hamiltonians in general. It is customary in such cases to express the time evolution operator, which is the formal solution to

$$\begin{eqnarray}&&\frac{\partial U(t,0)}{\partial t}=-{\rm i}H(t)U(t,0),\end{eqnarray} \tag{ 1 }$$

as

$$\begin{eqnarray}&&U(t,0)=\mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{t}H(\tau ){\rm d}\tau }}:=\mathop{{\rm lim} }\limits_{r\to \infty }{{\mathop{\prod }\limits_{j=0}^{r-1}}^{{}}}{{{\rm e}}^{-{\rm i}H(jt/r)t/r}}.\end{eqnarray} \tag{ 2 }$$

A wide array of approximation methods exist for U(t), including the Magnus expansion [17], Dyson series [18], Floquet theory [19], the Landau–Zener formula [20] and the adiabatic approximation.

The adiabatic approximation is widely used to approximate quantum dynamics in cases where rate of change of the time-dependent Hamiltonian is slow relative to an appropriate power of the minimum eigenvalue gap. In essence, the approximation states that if you prepare a system in an eigenstate of the Hamiltonian and evolve sufficiently slowly then the quantum system will evolve to the corresponding instantaneous eigenstate throughout the evolution. This lack of excitation throughout the process makes it analogous to reversible adiabatic processes in statistical mechanics. This analogy is not exact since the change in von Neumann entropy is also zero for any unitary process, and so the 'adiabatic' moniker persists for largely historical reasons.

Since the adiabatic approximation requires slow evolution, it is useful to consider how the approximation error scales as the speed of the transition from the initial to the final Hamiltonian decreases. This makes it natural to parameterize time via the variable s where

$$\begin{eqnarray}&&s=t/T,\end{eqnarray} \tag{ 3 }$$

and T is the total time for the adiabatic passage. While an adiabatic evolution occurs on $t\in [0,T]$, $s\in [0,1]$ regardless of the actual duration of the evolution. This means that if the Hamiltonian is re-parameterized as H(s), then we can increase T to make the evolution slower without fundamentally changing the form of the evolution.

We need to introduce some further notation before we can discuss the adiabatic approximation in greater detail. First we define $\left| n(s) \right\rangle$ to be the instantaneous eigenvectors of the time-dependent Hamiltonian,

$$\begin{eqnarray}&&H(s)\left| n(s) \right\rangle ={{E}_{n}}(s)\left| n(s) \right\rangle ,\end{eqnarray} \tag{ 4 }$$

and we make no assumptions about the ordering of E_n (i.e. we do not assume that ${{E}_{0}}\leqslant {{E}_{1}}$). Also, for notational simplicity, we define $\left| g(s) \right\rangle :=\left| {{E}_{0}}(s) \right\rangle$. We refer to this state as $\left| g(s) \right\rangle$ because it will represent the ground state in many practical examples of adiabatic QIP. The eigenvalue gaps will also be key to our analysis and so we use the following notation for them:

$$\begin{eqnarray}&&{{\gamma }_{\mu ,\nu }}(s)\;:={{E}_{\mu }}(s)-{{E}_{\nu }}(s).\end{eqnarray} \tag{ 5 }$$

The adiabatic approximation is often expressed in many different ways. The simplest of these states that

$$\begin{eqnarray}&&U(1,0)\left| g(0) \right\rangle \approx {{{\rm e}}^{-{\rm i}\int _{0}^{1}{{E}_{0}}\left( \xi \right)T{\rm d}\xi }}\left| g(1) \right\rangle +O\left( \frac{1}{T} \right).\end{eqnarray} \tag{ 6 }$$

In general the adiabatic approximation holds if

$$\begin{eqnarray}&&T\gg \frac{{{{\rm max} }_{s}}(\parallel \;{{\partial }_{s}}H(s){{\parallel }^{a}}+\parallel \;\partial _{s}^{2}\;H(s){{\parallel }^{b}}+\parallel \;\partial _{s}^{3}\;H(s){{\parallel }^{c}})}{{{{\rm min} }_{s}}{{\gamma }_{g,1}}{{(s)}^{d}}},\end{eqnarray} \tag{ 7 }$$

for integers $a,b,c$ and d that depend on the properties of the Hamiltonian [21–24]. A common misconception is that the adiabatic approximation holds if

$$\begin{eqnarray}&&T\gg \frac{{{{\rm max} }_{s}}\parallel \;{{\partial }_{s}}H(s)\parallel }{{{{\rm min} }_{s}}{{\gamma }_{0,1}}{{(s)}^{2}}},\end{eqnarray} \tag{ 8 }$$

although this criteria is appropriate for sufficiently slow evolutions under smoothly varying Hamiltonians [21, 22]. On the other hand, it is also known to fail in cases where there are resonant transitions which can occur in problems where $\parallel \partial _{s}^{2}H(s)\parallel$ is an increasing function of T [25, 26]. This can cause higher-order terms in the adiabatic perturbation series, which are typically negligible, to dominate the error in the approximation.

We refer to such results as zeroth-order adiabatic theorems, because they provide an estimate of the error that is correct to the zeroth order in ${{T}^{-1}}$, meaning that they simply tell you that the error is zero if the adiabatic process is infinitely slow. In order to show that we can combine different adiabatic evolutions to cancel the error, we need to have a sharper adiabatic condition that approximates the error to at least $O(1/T)$. It is necessary for us to use a first-order adiabatic approximation, which provides us with the error in the adiabatic approximation correct to $O(1/{{T}^{2}})$ [21]:

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} (\mathbb{I}-\left| g(1) \right\rangle \left\langle g(1) \right|)U\left( 1,0 \right)\left| g\left( s \right) \right\rangle \\ \quad =\left. \mathop{\sum }\limits_{n\ne g}{{{\rm e}}^{-{\rm i}\int _{0}^{1}{{E}_{n}}(\xi )T{\rm d}\xi }}\frac{\left\langle \dot{n}(s)|g(s) \right\rangle {{{\rm e}}^{{\rm i}\int _{0}^{s}{{\gamma }_{g,n}}(\xi ){\rm d}\xi T}}}{-{\rm i}{{\gamma }_{g,n}}\left( s \right)T} \right|_{s=0}^{1}\left| n\left( 1 \right) \right\rangle +O\left( \frac{1}{{{T}^{2}}} \right) \\ \end{array}.\end{eqnarray} \tag{ 9 }$$

This result can easily be found through the use of path integral methods [27–29]. Upper bounds on the magnitude of the sum of all $O(1/{{T}^{2}})$ terms are given in [21].

Equation (9) tells us something surprising: the leading order contribution to the error in the adiabatic approximation does not depend on the minimum gap but rather on the eigenvalue gap at the beginning and the end of the evolution, which motivates taking $\dot{H}(s)=0$ or equivalently $\left\langle \dot{n}(s)|g(s) \right\rangle =0$ on the boundary as per boundary cancellation methods [10]. The apparent contradiction posed by (9) is easily resolved. Adiabatic conditions like (7) give criteria for the convergence of the adiabatic perturbation series of $U(T,0)$ in powers of $1/T$, and equations such as (9) give a truncated expression for the power series. This means that after a critical evolution speed, the error in the adiabatic approximation no longer depends on the minimum gap, whereas the error depends crucially on the minimum gap before this point. We refer to the regime where the minimum gap dictates the error as the Landau–Zener regime and the regime where it does not as the adiabatic regime.

Similarly, the first-order adiabatic theorem relies on several conditions outlined in [21]. First, the Hamiltonian must be twice differentiable and three times piecewise differentiable with all such derivatives upper-bounded by a constant. Second, the system must already be in the adiabatic regime (i.e. the $\Theta (1/T)$ contribution to the error is much greater than the sum of all $O(1/{{T}^{2}})$ contributions). Third, we require that the norm of the Hamiltonian be upper-bounded by a constant for all times during the evolution. These criteria guarantee the validity of (9).

A common way to reduce errors in both the Landau–Zener regime and the adiabatic regime is to change the path used in the adiabatic evolution. The most frequently used adiabatic path, known as linear interpolation, is

$$\begin{eqnarray}&&H(s)=(1-s){{H}_{0}}+s{{H}_{1}},\end{eqnarray} \tag{ 10 }$$

where H₀ is the initial Hamiltonian and H₁ is the final Hamiltonian. There are, of course, many ways that one could imagine transitioning from the initial Hamiltonian to the final Hamiltonian. Each of these ways represents a particular 'adiabatic path' and (10) is known as the linear adiabatic path. More generally we could consider a path of the form

$$\begin{eqnarray}&&H(s)=(1-f(s)){{H}_{0}}+g(s){{H}_{1}},\end{eqnarray} \tag{ 11 }$$

where $f(0)=g(0)=0$ and $f(1)=g(1)=1$. Such paths can be extremely important for adiabatic quantum computing because they allow the evolution to slow down through, or even avoid, parts of the evolution that contribute substantially to the error; however, here we assume the simple case of $g(s)=f(s)$. We do not require that the range of f be $[0,1]$ here. In fact, some of the adiabatic paths that we consider will attain negative values and values greater than 1.

Other examples of non-linear paths include local adiabatic evolution, which seeks to minimize the error in the Landau–Zener regime by choosing the evolution speed to be smallest near the minimum gap. Boundary cancellation methods, on the other hand, choose paths that minimize the error in the adiabatic regime by choosing the evolution speed to be zero at s = 0 and s = 1. These two strategies are seemingly orthogonal. At present there is no known method that combines the best features of local adiabatic paths and the paths yielded by boundary cancellation methods. Our work provides a way to achieve this, thereby illustrating that controlled adiabatic evolution affords greater power than conventional adiabatic evolution.

The central idea behind our approach is to use a gadget that was recently proposed in [14] to non-deterministically implement the weighted average of two or more adiabatic evolutions. This idea of using controlled adiabatic evolutions and measurement has been recently explored by Itay Hen [15] and is also used in holonomic quantum computing [16]; however, these results do not consider using coherent control and measurement to suppress diabatic errors. The gadget that we use for this averaging process is given in figure 1. The circuit in figure 1 probabilistically implements linear combinations of unitary operations, as seen through the following argument:

$$\begin{eqnarray}\begin{array}{llr} \left| \psi \right\rangle \left| 0 \right\rangle & \to & \left| \psi \right\rangle \left( {\rm cos} \theta \left| 0 \right\rangle +{\rm sin} \theta \left| 1 \right\rangle \right) \\ {} & \to & {\rm cos} \theta {{U}_{A}}\left| \psi \right\rangle \left| 0 \right\rangle +{\rm sin} \theta {{U}_{B}}\left| \psi \right\rangle \left| 1 \right\rangle \\ {} & \to & \left( {{{\rm cos} }^{2}}\theta {{U}_{A}}+{{{\rm sin} }^{2}}\theta {{U}_{B}} \right)\left| \psi \right\rangle \left| 0 \right\rangle +{\rm sin} \theta {\rm cos} \theta \left( {{U}_{B}}-{{U}_{A}} \right)\left| \psi \right\rangle \left| 1 \right\rangle . \\ \end{array}\end{eqnarray} \tag{ 12 }$$

We then see that if the ancilla register is measured to be 0 then the circuit performs a weighted combination of U_A and U_B on the state $\left| \psi \right\rangle$; otherwise, the circuit implements the difference between the two operators:

$$\begin{eqnarray}&&p(0)\geqslant 1-\left\|\left( {{U}_{A}}-{{U}_{B}} \right)\left| \psi \right\rangle \right \|{{}^{2}}.\end{eqnarray} \tag{ 13 }$$

The generalization to cases where multiple U_A and U_B are used is trivial: it simply involves increasing the number of qubits used to control the overall rotation [29]. Such circuits, or variants thereof, are also used in [30, 31].

For the case of adiabatic evolution, we know that to the zeroth order

$$\begin{eqnarray}\begin{array}{rcl} {{U}_{A}}({{T}_{A}},0)\left| g(0) \right\rangle & = & {{{\rm e}}^{-{\rm i}\int _{0}^{1}{{E}_{0}}({{f}_{A}}(s)){\rm d}s{{T}_{A}}}}\left| g(1) \right\rangle +O(1/T), \\ {{U}_{B}}({{T}_{B}},0)\left| g(0) \right\rangle & = & {{{\rm e}}^{-{\rm i}\int _{0}^{1}{{E}_{0}}({{f}_{B}}(s)){\rm d}s{{T}_{B}}}}\left| g(1) \right\rangle +O(1/T), \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where $T={\rm max} \{{{T}_{A}},{{T}_{B}}\}$. This means that, to leading order, both U_A and U_B generate the same evolution up to a global phase and hence we expect the success probability to be high if the phases picked up by $\left| g \right\rangle$ under both evolutions are comparable.

Rather than choosing different paths that apply the same phase to $\left| g(s) \right\rangle$, we counter-rotate the evolution of each eigenstate by including an additional phase to each unitary. This affords us much greater freedom to choose adiabatic paths for U_B and U_A. In particular, we choose these phases such that

$$\begin{eqnarray}\begin{array}{rcl} {{U}_{A}}({{T}_{A}},0)\left| g(0) \right\rangle & = & {{{\rm e}}^{{\rm i}\int _{0}^{1}{{E}_{0}}({{f}_{A}}(s)){\rm d}s{{T}_{A}}}}\left( \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}H({{f}_{A}}(s)){\rm d}s{{T}_{A}}}}\left| g(1) \right\rangle \right), \\ {{U}_{B}}({{T}_{B}},0)\left| g(0) \right\rangle & = & {{{\rm e}}^{{\rm i}\int _{0}^{1}{{E}_{0}}({{f}_{B}}(s)){\rm d}s{{T}_{B}}}}\left( \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}H({{f}_{B}}(s)){\rm d}s{{T}_{B}}}}\left| g(1) \right\rangle \right). \\ \end{array}\end{eqnarray} \tag{ 15 }$$

Note that in certain highly symmetric cases, such as those considered in section 5, these phase factors are identical and can therefore be omitted.

We see from the choices of phases in (15) that (13) gives the failure probability of the linear combination $O(1/T)$ in the limit of large T. This means that the failure probability will typically be extremely small for adiabatic processes. Even if a failure is observed, the gadget in figure 1 informs the user that a failure has occurred and the state preparation process can be repeated until success is obtained. We see from numerical experiments that the failure probability of these circuits has a near-negligible impact on the cost of coherently controlled adiabatic state preparation in the adiabatic regime.

Generalization of these ideas to cases where more than two unitary evolutions are averaged is straightforward and is discussed in detail in [14]. We present the two-unitary case explicitly here since the majority of our results focus on averaging two different adiabatic evolutions.

Our first approach is a generalization of the strategy employed by Wiebe and Babcock in [32], which suppresses the dominant transition in the adiabatic passage for adiabatic paths satisfying

$$\begin{eqnarray}&&{{\left. \frac{\left\langle \dot{n}(s)|g(s) \right\rangle }{{{\gamma }_{g,n}}(s)} \right|}_{s=0}}={{\left. \frac{\left\langle \dot{n}(s)|g(s) \right\rangle }{{{\gamma }_{g,n}}(s)} \right|}_{s=1}},\end{eqnarray} \tag{ 16 }$$

by choosing the evolution time T appropriately. Our strategy is to suppress a single transition, not by choosing a single time and requiring a symmetry condition as per [32], but by interfering the adiabatic evolution with a dual evolution, as suggested in figure 2. This allows such errors to be suppressed for any evolution time and any primary path. We also provide a method for suppressing the two most significant diabatic transitions in appendix A.

Zoom In Zoom Out Reset image size

We wish to choose, for fixed H_A, an adiabatic path that quadratically suppresses the transition $\left| g(0) \right\rangle \to \left| e(1) \right\rangle$ where $\left| e(s) \right\rangle$ is any given instantaneous eigenstate of H(s) that is orthogonal to the $\left| g(s) \right\rangle$. From (9) and (15), we see that if we combine ${{U}_{A}}({{T}_{A}},0)$ with ${{U}_{B}}({{T}_{B}},0)$ and achieve a successful measurement outcome, then we obtain a result proportional to

$$\begin{eqnarray}&&({\rm cos} {{(\theta )}^{2}}{{U}_{A}}({{T}_{A}},0)+{\rm sin} {{(\theta )}^{2}}{{U}_{B}}({{T}_{B}},0))\left| g(0) \right\rangle =\left| g(1) \right\rangle +O(1/T)\;:=\left| \phi \right\rangle .\end{eqnarray} \tag{ 17 }$$

So to leading order, the linear combination will give the correct result. Then, using (9), it is clear that

$$\begin{eqnarray}&&\begin{array}{l} \left| e(1) \right\rangle \left\langle e(1) \right.\left| \phi \right\rangle \propto {{{\rm cos} }^{2}}(\theta ) \\ \qquad \times \left[ \frac{\left\langle \dot{e}_{1}^{A}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,e}^{A}(1){{T}_{A}}}-{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,{{e}^{A}}}^{A}\left( \xi \right){{T}_{A}}{\rm d}\xi }}\frac{\left\langle {{{\dot{e}}}^{A}}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,e}^{A}\left( 0 \right){{T}_{A}}} \right] \\ \qquad +{{{\rm sin} }^{2}}(\theta )\left[ \frac{\left\langle {{{\dot{e}}}^{B}}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,e(1){{T}_{B}}}^{B}}-{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,e}^{B}\left( \xi \right){{T}_{B}}{\rm d}\xi }}\qquad \frac{\left\langle {{{\dot{e}}}^{B}}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,e}^{B}\left( 0 \right){{T}_{B}}} \right] \\ \qquad +O(1/{{T}^{2}}). \\ \end{array}\end{eqnarray} \tag{ 18 }$$

This transition can therefore be canceled, to $O(1/{{T}^{2}})$, by choosing $\theta ,{{T}_{B}}$ and f_B such that the weighted average of the diabatic transitions to the state $\left| e \right\rangle$ is zero:

$$\begin{eqnarray}\begin{array}{rcl} 0 & = & {{{\rm cos} }^{2}}(\theta )\left[ \frac{\left\langle \dot{e}_{1}^{A}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,e}^{A}(1){{T}_{A}}}-{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,{{e}^{A}}}^{A}(\xi ){{T}_{A}}d\xi }}\frac{\left\langle {{{\dot{e}}}^{A}}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,e}^{A}(0){{T}_{A}}} \right] \\ {} & {} & +{{{\rm sin} }^{2}}(\theta )\left[ \frac{\left\langle {{{\dot{e}}}^{B}}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,e}^{B}(1){{T}_{B}}}-{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,e}^{B}\left( \xi \right){{T}_{B}}d\xi }}\frac{\left\langle {{{\dot{e}}}^{B}}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,e}^{B}(0){{T}_{B}}} \right] \\ \end{array}\end{eqnarray} \tag{ 19 }$$

where $\left\langle \dot{e}(s)|g(s) \right\rangle =\frac{\left\langle e(s)\left| \dot{H}(s) \right|g(s) \right\rangle }{{{\gamma }_{e,g}}(s)}$. Thus it is reasonable to expect that this condition can be met by choosing $\theta$ and f_B properly. The remaining question is, how can this be done in practice? We provide two strategies for finding f_B for any fixed f_A such that these errors cancel to leading order.

3.1. Partially anti-symmetric combination

Our first method chooses the paths f_A and f_B to satisfy an anti-symmetric condition on the derivatives at the beginning and end of the evolution. This approach is most useful in cases where it is desirable for f_B to be as similar to f_A as possible. When optimizing these paths, it is important to note that although f_A and f_B are arbitrary interpolation functions that describe the adiabatic paths, they are constrained to obey

$$\begin{eqnarray}&&{{f}_{A}}(0)={{f}_{B}}(0)=0,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{f}_{A}}(1)={{f}_{B}}(1)=1.\end{eqnarray} \tag{ 21 }$$

Furthermore, let us choose f_B such that its derivatives are symmetric with f_A at s = 0 and anti-symmetric at s = 1

$$\begin{eqnarray}\begin{array}{rcl} {{\left. {{{\dot{f}}}_{B}}\left( s \right) \right|}_{s=0}} & = & {{\left. {{{\dot{f}}}_{A}}\left( s \right) \right|}_{s=0}}, \\ {{\left. {{{\dot{f}}}_{B}}\left( s \right) \right|}_{s=1}} & = & -{{\left. {{{\dot{f}}}_{A}}\left( s \right) \right|}_{s=1}}. \\ \end{array}\end{eqnarray} \tag{ 22 }$$

Then using (19), (22) simplifies to

$$\begin{eqnarray}\begin{array}{rcl} \left( \frac{{{{\rm cos} }^{2}}(\theta )}{{{T}_{A}}}-\frac{{{{\rm sin} }^{2}}(\theta )}{{{T}_{B}}} \right)\frac{\left\langle \dot{e}(1)|g(1) \right\rangle }{{{\gamma }_{g,e}}\left( 1 \right)} & = & \left( \frac{{{{\rm cos} }^{2}}(\theta )}{{{T}_{A}}}{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,e}^{A}\left( \xi \right){{T}_{A}}{\rm d}\xi }} \right. \\ {} & {} & +\left. \frac{{{{\rm sin} }^{2}}(\theta )}{{{T}_{B}}}{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,e}^{B}\left( \xi \right){{T}_{B}}{\rm d}\xi }} \right)\frac{\left\langle \dot{e}(0)|g(0) \right\rangle }{{{\gamma }_{g,e}}\left( 0 \right)}. \\ \end{array}\end{eqnarray} \tag{ 23 }$$

Equation (23) can be satisfied for any f_A and T_A by setting

$$\begin{eqnarray}&&{{T}_{B}}=\frac{\int _{0}^{1}\gamma _{g,e}^{A}\left( \xi \right){\rm d}\xi {{T}_{A}}+\left( 2n+1 \right)\pi }{\int _{0}^{1}\gamma _{g,e}^{B}\left( \xi \right){\rm d}\xi }\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\theta ={\rm arctan} \left( \sqrt{\frac{{{T}_{B}}}{{{T}_{A}}}} \right).\end{eqnarray} \tag{ 25 }$$

This solution reduces to that of [32] in the limit as ${{T}_{A}}\to 0$; however, a non-trivial secondary path will always be needed if the symmetry condition demanded by [32] is not held.

An important consequence of taking the derivatives to be negative at s = 1 is that there exists $s^{\prime}$ such that $f(x)\gt 1$ for all $x\in (s^{\prime} ,1)$. This is a consequence of the fact that H(s) is twice differentiable and hence $f_{B}^{\prime }$ is continuous, from which the result directly follows from the mean value theorem. Thus f_b does not monotonically approach 1 as $s\to 1$, but rather it overshoots the value and then reverses direction to end the evolution at s = 1. Such reversals of direction are analogous to the backwards time steps used in Trotter–Suzuki methods and, although non-traditional, are not necessarily problematic for adiabatic evolution.

We see from this discussion that controlled adiabatic paths can be used to suppress diabatic errors in ways that are impossible using traditional adiabatic optimization strategies. In particular, for any optimization strategy, such as local adiabatic evolution, we can always find a second path to add to the primary path to suppress a chosen transition to one order higher. These ideas can also be generalized to suppress more than one transition; however, finding a closed form solution is difficult in such cases. We discuss generalizing this method to simultaneously suppress two diabatic transitions in appendix A. A drawback of this approach is that it cannot be used for arbitrary small times because 24 forces a difference at least $\frac{\pi }{\int _{0}^{1}{{\gamma }_{g,e}}\left( \xi \right)}$ between evolution times. We address this issue below by providing a method that does not require a shift in time, but requires a more substantial deformation to the primary adiabatic path.

3.2. Completely anti-symmetric combination

An alternative approach is to set the derivatives for the second path to be completely antisymmetric:

$$\begin{eqnarray}\begin{array}{rcl} {{\left. {{{\dot{f}}}_{B}}\left( s \right) \right|}_{s=0}} & = & {{\left. -{{{\dot{f}}}_{A}}\left( s \right) \right|}_{s=0}} \\ {{\left. {{{\dot{f}}}_{B}}\left( s \right) \right|}_{s=1}} & = & {{\left. -{{{\dot{f}}}_{A}}\left( s \right) \right|}_{s=1}}. \\ \end{array}\end{eqnarray} \tag{ 26 }$$

Pluging (26) into (19), we obtain

$$\begin{eqnarray}\begin{array}{lcl} {} & {} & \left( \frac{{{{\rm cos} }^{2}}(\theta )}{{{T}_{A}}}-\frac{{{{\rm sin} }^{2}}(\theta )}{{{T}_{B}}} \right)\frac{\left\langle \dot{e}(1)|g(1) \right\rangle }{\gamma _{g,e}^{B}\left( 1 \right)} \\ {} & = & \left( \frac{{{{\rm cos} }^{2}}(\theta )}{{{T}_{A}}}{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,e}^{A}\left( \xi \right){{T}_{A}}{\rm d}\xi }}-\left. \frac{{{{\rm sin} }^{2}}(\theta )}{{{T}_{B}}}{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,e}^{B}\left( \xi \right){{T}_{B}}{\rm d}\xi }} \right) \right. \\ {} & {} & \times \frac{\left\langle \dot{e}(0)|g(0) \right\rangle }{\gamma _{g,e}^{B}\left( 0 \right)}. \\ \end{array}\end{eqnarray} \tag{ 27 }$$

The error is suppressed when (25) holds and

$$\begin{eqnarray}&&{{T}_{B}}=\frac{\int _{0}^{1}\gamma _{g,e}^{A}\left( \xi \right){{T}_{A}}{\rm d}\xi +2n\pi }{\int _{0}^{1}\gamma _{g,e}^{B}\left( \xi \right){\rm d}\xi }.\end{eqnarray} \tag{ 28 }$$

In other words, the gap integrals for both paths must be equivalent modulo $2\pi$. This removes the difficulty with offsetting one of the times. However, in this case, the path f_B both begins and ends the evolution by moving backwards. Alternatively, we can modify one boundary from each path. This backwards motion at s = 0 means that there exists $\delta \gt 0$ such that the range of f_B(s) is within $[-\delta ,1+\delta ]$. Again, this use of backward evolution is atypical of conventional approaches to adiabatic evolution where the additional evolution time/speed required by backwards evolution would tend to be detrimental. In contrast, such backwards evolutions can lead to substantial reductions in the cost for coherently controlled adiabatic evolution.

3.3. Interpolation

There are many ways that these requirements can be satisfied by a dual path to f_A. The way that we satisfy these requirements is by adding a smooth polynomial continuation of f_A about s = 1 that allows the derivative to loop around and attain the opposite value. This interpolation must have piecewise continuous third derivatives in order to guarantee that the $O(1/{{T}^{2}})$ terms will remain sub-dominant in the limit of large T. This naturally leads to a quartic interpolation that takes the following form for a partially anti-symmetric combination:

$$\begin{eqnarray}{{f}_{B}}(s)=\left\{ \begin{array}{ccccccccccccccc} {{f}_{A}}(s) & s\lt 1-\Delta \\ e{{s}^{4}}+d{{s}^{3}}+c{{s}^{2}}+bs+a & s\geqslant 1-\Delta \\ \end{array} \right.,\end{eqnarray} \tag{ 29 }$$

where Δ is a free parameter that controls when f_B switches from the original adiabatic path f_A to the polynomial interpolation. The parameters are then set by requiring

$$\begin{eqnarray}\begin{array}{rcl} {{f}_{B}}(1) & = & {{f}_{A}}(1)=1 \\ {{f}_{B}}(1-\Delta ) & = & {{f}_{A}}(1-\Delta ) \\ {{{\dot{f}}}_{B}}(1) & = & -{{{\dot{f}}}_{A}}(1) \\ {{{\dot{f}}}_{B}}(1-\Delta ) & = & {{{\dot{f}}}_{A}}(1-\Delta ) \\ {{\ddot{f}}_{B}}(1-\Delta ) & = & {{\ddot{f}}_{A}}(1-\Delta ). \\ \end{array}\end{eqnarray} \tag{ 30 }$$

We could also equally well have chosen the backwards evolution to start at s = 0 rather than s = 1. Although seemingly arbitrary, this choice can have a substantial impact on the error depending on whether the gap is larger at s = 0 and s = 1. We also make use of this fact later in section 5 where we exploit this fact to suppress every transition simultaneously for Hamiltonians that satisfy a certain symmetry property.

The case of fully anti-symmetric boundaries is similar, except now two polynomial interpolations are needed:

$$\begin{eqnarray}{{f}_{B}}(s)=\left\{ \begin{array}{ccccccccccccccc} {{f}_{A}}(s) & \Delta /2\lt s\lt 1-\Delta /2 \\ e{{s}^{4}}+d{{s}^{3}}+c{{s}^{2}}+bs+a & s\leqslant \Delta /2 \\ e^{\prime} {{s}^{4}}+d^{\prime} {{s}^{3}}+c^{\prime} {{s}^{2}}+b^{\prime} s+a^{\prime} & s\geqslant 1-\Delta /2 \\ \end{array} \right.,\end{eqnarray} \tag{ 31 }$$

where Δ is a free parameter that controls how rapidly f_B switches from the original adiabatic path, f_A, to the polynomial interpolation. The parameters are then set by requiring

$$\begin{eqnarray}\begin{array}{rcl} {{f}_{B}}(0) & = & {{f}_{A}}(0)=0 \\ {{f}_{B}}(1) & = & {{f}_{A}}(1)=1 \\ {{f}_{B}}(1-\Delta /2) & = & {{f}_{A}}(1-\Delta /2) \\ {{f}_{B}}(\Delta /2) & = & {{f}_{A}}(\Delta /2) \\ {{{\dot{f}}}_{B}}(1) & = & -{{{\dot{f}}}_{A}}(1) \\ {{{\dot{f}}}_{B}}(0) & = & -{{{\dot{f}}}_{A}}(0) \\ {{{\dot{f}}}_{B}}(1-\Delta /2) & = & {{{\dot{f}}}_{A}}(1-\Delta /2) \\ {{{\dot{f}}}_{B}}(\Delta /2) & = & {{{\dot{f}}}_{A}}(\Delta /2) \\ {{\ddot{f}}_{B}}(\Delta /2) & = & {{\ddot{f}}_{A}}(\Delta /2) \\ {{\ddot{f}}_{B}}(1-\Delta /2) & = & {{\ddot{f}}_{A}}(1-\Delta /2). \\ \end{array}\end{eqnarray} \tag{ 32 }$$

In particular, the coefficients in (31) can then be found by substituting (31) into (32).

It then follows that for any fixed path f_A, we can choose f_B such that the dominant transition is suppressed to $O(1/{{T}^{2}})$. This opens the possibility that our error suppression methods may allow adiabatic state preparation to be performed using less evolution time (or equivalently, fewer gates on a quantum computer) than existing methods. However, local adiabatic evolution is known to be optimal for performing adiabatic Groverʼs search [8, 9], so we cannot expect that the algorithm will outperform all existing adiabatic algorithms in every time regime. We will see below that although our method does not outperform local adiabatic evolution for short times, it can come very close to matching its performance while giving substantially reduced error for slow evolutions.

We focus our numerical results on the case of adiabatic Groverʼs search. The Hamiltonian for adiabatic Groverʼs search is

$$\begin{eqnarray}&&H(f(s))=(1-f(s))({\mathbb{1}}-\left| {{+}^{\,\otimes n}} \right\rangle \left\langle {{+}^{\,\otimes n}} \right|)+f(s)({\mathbb{1}}-\left| 0 \right\rangle \left\langle 0 \right|),\end{eqnarray} \tag{ 33 }$$

and sufficiently slow evolution of this Hamiltonian causes the initial eigenstate $\left| {{+}^{\,\otimes n}} \right\rangle$ to transition to the marked state $\left| 0 \right\rangle$ as per Groverʼs search. Local adiabatic evolution is known to be optimal for adiabatic Groverʼs search [8, 9], meaning that the quadratic speedup over classical algorithms is attained for the adiabatic path.

The path for local adiabatic evolution, for cases where the search space is N-dimensional, is

$$\begin{eqnarray}&&f(s)=\frac{\sqrt{N-1}-{\rm tan} \left[ {\rm arctan} (\sqrt{N-1})(1-2s) \right]}{2\sqrt{N-1}}.\end{eqnarray} \tag{ 34 }$$

Our goal is to compare the cost of performing this adiabatic quantum algorithm using local adiabatic evolution to the cost incurred by using our methods. We choose f_A to be the path given by local adiabatic evolution, whereas f_B is taken to be a continuation of the local adiabatic path that satisfies (22) or (26) in all of the following numerical examples.

There are several ways that the cost of an adiabatic algorithm can be measured. The most straightforward method is to compare the time required for the evolution. Although this cost metric is appropriate in cases where the norm of the Hamiltonian is fixed, it is not appropriate for comparing different adiabatic evolutions because the energy required for both paths may differ substantially. Since there is a duality between energy and time in quantum mechanics, a fast evolution that requires a lot of energy may be dynamically equivalent to a slow evolution that requires little energy. Thus we need to consider not just the time but also the energy. For this reason, we use the following cost metric for the case where we combine j evolutions (where $j\geqslant 1$ ):

$$\begin{eqnarray}&&{\rm Cost}=\mathop{{\rm max} }\limits_{j}\left\{ \int _{0}^{1}\parallel {{H}_{j}}(s)\parallel {\rm d}s{{T}_{j}} \right\}/P(0),\end{eqnarray} \tag{ 35 }$$

where $P(0)$ is the success probability of the gadget that is given by (13). Here we implicitly assume that the cost of the rotations and the control logic is negligible and that each of the evolutions can be implemented in parallel. These assumptions may not hold in general, but they are appropriate for quantum computer simulations of such adiabatic evolutions because the query complexity of such evolutions depends on the maximum evolution time chosen rather than the total evolution time. This point will be made clear in section 7.

We see in figure 3 that including the second adiabatic path with partially anti-symmetric boundary conditions (as per section 3.1) to LAE yields comparable performance to LAE for short evolutions and also provides the improved scaling of boundary cancellation methods in the adiabatic regime. In particular, the second path follows the interpolation strategy of (29): it follows LAE (i.e. (34)) until s = 0.8 and then smoothly transitions to a fourth-order polynomial. Unlike the method of [32], this gives superior scaling over a discrete set of points, although it does enforce a minimum evolution time, as discussed in section 3.1. An important drawback of this method is that there is a manifest lack of symmetry in the derivatives in the adiabatic regime for this method. This means that the adiabatic interference effects that appear in the LAE and boundary cancellation paths will not appear here [32]. Note that if the Search Hamiltonian did not have symmetric derivatives or spectrum, the adiabatic interference effects would not appear and so they are an artifact of having a highly structured test Hamiltonian.

Zoom In Zoom Out Reset image size

Figure 4 tells a similar story. In that case we use fully anti-symmetric boundary conditions and add a second path that interpolates between LAE and polynomial evolution as per (31) with $\Delta =0.2$. This also corresponds to evolution under LAE for 80% of the time. The value $\Delta =0.2$ was chosen arbitrarily and in practice this value could be optimized to further reduce error. Unlike the case in figure 3, adiabatic interference patterns are again visible in the adiabatic regime because the two polynomials used to create the fully anti-symmetric boundary conditions between the two paths at s = 0 and s = 1 ensure that the derivatives are the same at the boundary, thereby allowing such interference effects to emerge again; this causes the error to be substantially reduced on a discrete subset of points, as per [32]. As a consequence, we can clearly see that our method substantially outperforms both methods for a range of evolutions with cost ranging from $[50,100]$, due in part to the presence of adiabatic interference effects that are absent from the boundary cancellation method.

Zoom In Zoom Out Reset image size

In both of the cases considered, our methods are less effective at suppressing errors in the adiabatic regime than boundary cancellation methods. This is because the $O(1/{{T}^{2}})$ terms in the error in the adiabatic approximation also depend on $\dot{H}(s)$. Such terms are zero for boundary cancellation methods and so we generically expect from the triangle inequality that boundary cancellation will lead to less error in this regime. An important point to note is that although these test cases do not outperform LAE for fast evolutions or boundary cancellation methods for slow evolutions, they can outperform both methods for evolutions that operate at an intermediate speed. This implies that these methods are not just a compromise between the two approaches: they also provide superior scaling in a region that is badly addressed by existing adiabatic optimization methods.

We now consider suppressing errors for Hamiltonians where $H(0)$ and $H(1)$ have the same spectra. Although restrictive, this condition is satisfied in many natural problems [3, 8, 32]. This symmetry is very useful because it guarantees that two adiabatic interpolations exist between H₀ and H₁ such that the amplitudes for every state orthogonal to $\left| g(1) \right\rangle$ that arise due to f_A are equal and opposite to those that arise under f_B. This means that the linear combination will simultaneously suppress diabatic leakage into every state. In contrast, the methods discussed in sections 3.1 and 3.2 guarantee this only for a single (but arbitrarily chosen) transition. First let us assume that the following conditions are met for f_A(s) and f_B(s)

$$\begin{eqnarray}&&\int _{0}^{s}\gamma _{g,n}^{A}\left( {{f}_{A}}\left( \xi \right) \right){\rm d}\xi =\int _{0}^{s}\gamma _{g,n}^{B}\left( {{f}_{B}}\left( \xi \right) \right){\rm d}\xi \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{{\left. \frac{{{\left\langle \dot{n}({{f}_{A}}(s))|m({{f}_{A}}(s)) \right\rangle }_{A}}}{\gamma _{g,n}^{A}\left( {{f}_{A}}(s) \right)} \right|}_{s=0}}=-{{\left. \frac{{{\left\langle \dot{n}({{f}_{B}}(s))|g({{f}_{B}}(s)) \right\rangle }_{B}}}{\gamma _{g,n}^{B}\left( {{f}_{B}}(s) \right)} \right|}_{s=0}}\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{{\left. \frac{{{\left\langle \dot{n}({{f}_{A}}(s))|m({{f}_{A}}(s)) \right\rangle }_{A}}}{\gamma _{g,n}^{A}\left( {{f}_{A}}(s) \right)} \right|}_{s=1}}=-{{\left. \frac{{{\left\langle \dot{n}({{f}_{B}}(s))|g({{f}_{B}}(s)) \right\rangle }_{B}}}{\gamma _{g,n}^{B}\left( {{f}_{B}}(s) \right)} \right|}_{s=1}}\end{eqnarray} \tag{ 38 }$$

for all states $\left| n \right\rangle \ne \left| g \right\rangle$. We will see that these conditions can always be met if the spectrum of H(s) is symmetric about $s=1/2$.

Such conditions do not naturally arise for all adiabatic passages but there are many examples where such Hamiltonians are natural. A natural example is the search Hamiltonian; however, such an application is trivial because the quantum dynamics occurs within a two-dimensional subspace. Other examples occur in adiabatic gates [15, 33, 34] and holonomic quantum computing [16, 35].

After substituting (36)–(38) into (18), we find

$$\begin{eqnarray}\begin{array}{rcl} |\left| n(1) \right\rangle \left\langle n(1) \right|\left| \phi \right\rangle | & = & \left| {{{\rm cos} }^{2}}(\theta )\left[ \frac{\left\langle \dot{n}_{1}^{A}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,n}^{A}\left( 1 \right){{T}_{A}}}-{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,{{n}^{A}}}^{A}\left( \xi \right){{T}_{A}}{\rm d}\xi }}\frac{\left\langle {{{\dot{n}}}^{A}}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,n}^{A}\left( 0 \right){{T}_{A}}} \right] \right. \\ {} & {} & -\left. {{{\rm sin} }^{2}}(\theta )\left[ \frac{\left\langle {{{\dot{n}}}^{A}}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,n}^{A}\left( 1 \right){{T}_{B}}}-{{{\rm e}}^{+{\rm i}\int _{0}^{1}\gamma _{g,n}^{A}\left( \xi \right){{T}_{B}}{\rm d}\xi }}\frac{\left\langle {{{\dot{n}}}^{A}}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,n}^{A}\left( 0 \right){{T}_{B}}} \right] \right| \\ {} & {} & +O(1/{{T}^{2}}). \\ \end{array}\end{eqnarray} \tag{ 39 }$$

It is then clear that if we take $\theta =\pi /4$ and T_A=T_B then every transition will be suppressed from $O(1/T)$ to $O(1/{{T}^{2}})$ under these assumptions.

The question remaining is, when can we make these conditions hold? A natural case that covers a wide range of adiabatic protocols is the case where the eigenvalue gap is symmetric. That is to say that ${{\gamma }_{g,n}}(s)={{\gamma }_{g,n}}(1-s)$ for all s and $n\ne g$. It is difficult to find a second adiabatic path that satisfies the conditions in section 3.3, for the choice ${{f}_{A}}=f,$ because the anti-symmetry required by (37), equation (37) necessitates the use of adiabatic paths similar to those in section 3.2. Such paths will typically violate (36) because including the reversal near s = 0 and s = 1 will change the gap integral.

A better approach is to modify both paths. It is easy to see by substitution that if we let f_B be given by (29) and (30) and then choose f_A to be the time reversed version of this path (i.e. ${{f}_{A}}(s)=1-{{f}_{B}}(1-s)$) then

$$\begin{eqnarray}\begin{array}{rcl} H({{f}_{A}}(s)) & = & (1-{{f}_{A}}(s)){{H}_{0}}+{{f}_{A}}(s){{H}_{1}} \\ {} & = & {{f}_{B}}(1-s){{H}_{0}}+(1-{{f}_{B}}(1-s)){{H}_{1}}. \\ \end{array}\end{eqnarray} \tag{ 40 }$$

The assumption that ${{\gamma }_{g,n}}(s)={{\gamma }_{g,n}}(1-s)$ then directly implies that ${{\gamma }_{g,n}}({{f}_{A}}(s))={{\gamma }_{g,n}}({{f}_{B}}(1-s)),$ which gives us the desired result of

$$\begin{eqnarray}&&\int _{0}^{1}{{\gamma }_{g,n}}({{f}_{B}}(s)){\rm d}s=\int _{0}^{1}{{\gamma }_{g,n}}({{f}_{A}}(s)){\rm d}s.\end{eqnarray} \tag{ 41 }$$

This fact becomes immediately obvious in light of figure 5, where the spectrum for a Search Hamiltonian with f_A and f_B chosen to be time reverses of each other is given. After substituting (41) into (39) and using the assumption that ${{\gamma }_{g,n}}(1)={{\gamma }_{g,n}}(0),$ we see that

$$\begin{eqnarray}&&\left| n(1) \right\rangle \left\langle n(1) \right|\left| \phi \right\rangle =O(1/{{T}^{2}}).\end{eqnarray} \tag{ 42 }$$

Thus for any adiabatic path parameterized by f(s) and any evolution time T, we can always choose two paths f_A and f_B that both incur diabatic errors that cancel to $O(1/{{T}^{2}})$. This fact is demonstrated numerically for a Hamiltonian that satisfies these requirements in figure 6. In contrast, without these assumptions, only the dominant transition is suppressed to $O(1/{{T}^{2}})$ which implies that the diabatic errors scale as $O(1/T)$ for sufficiently long evolutions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This result is much stronger than that of [32], which leads to suppression of all diabatic errors only if the gap integrals for each transition are rational multiples of each other and, even then, will only work at specially chosen values of T. This precludes the techniqueʼs use for almost all Hamiltonians. In contrast, coherent control of the adiabatic path allows all of the transitions to be suppressed for a wide class of adiabatic protocols and this result holds for any T. This clearly demonstrates that coherently controlled adiabatic evolution allows us to circumvent the limitations of existing adiabatic optimization schemes. We will also see this below, where we show that these methods can be used in concert with boundary cancellation methods.

In principle, our approach for improving the scaling of adiabatic evolution can be improved from $O(1/{{T}^{2}})$ to $O(1/{{T}^{3}})$ or higher by interfering more adiabatic paths. A challenge facing this approach is that the expression for the error in the adiabatic approximation becomes much more complex as we transition from first to second order in powers of $1/T$. As noted in [10, 13], these expressions become much simpler if we take several of the derivatives to be zero at s = 0 and s = 1. Here we use this simplification to show that our strategy can be used in concert with boundary cancellation methods. This allows the error scaling to be improved by one order without increasing the number of zero derivatives on the boundary. This is important because the time needed for the system to transition to the adiabatic regime increases with the number of derivatives that are set to zero at s = 0 and s = 1.

If the first m derivatives of H(s) are set to zero on the boundary then the remaining error in the adiabatic approximation is [32]:

$$\begin{eqnarray}&&\begin{array}{l} \left( \mathbb{I}-\left| g(1) \right\rangle \left\langle g(1) \right| \right)U\left( 1,0 \right)\left| g\left( s \right) \right\rangle \\ \quad =\left. \mathop{\sum }\limits_{n\ne g}{{{\rm e}}^{-{\rm i}\int _{0}^{1}{{E}_{n}}\left( \xi \right)T{\rm d}\xi }}\frac{\left\langle n(s) \right|(\partial _{s}^{m+1}H(s))\left| g(s) \right\rangle {{{\rm e}}^{{\rm i}\int _{0}^{s}{{\gamma }_{g,n}}(\xi ){\rm d}\xi T}}}{-{\rm i}\gamma _{g,n}^{m+2}(s){{T}^{m+1}}} \right|_{s=0}^{1} \\ \qquad \times \left| n\left( 1 \right) \right\rangle +O\left( \frac{1}{{{T}^{m+2}}} \right) \\ \end{array}\end{eqnarray} \tag{ 43 }$$

This expression is analogous to (9), as can be seen by substituting

$$\begin{eqnarray}\begin{array}{rcl} \left\langle \dot{n}(s)|m(s) \right\rangle & \to & \frac{\left\langle n(s) \right|(\partial _{s}^{m+1}H(s))\left| m(s) \right\rangle }{\gamma _{g,n}^{m+1}(s)}, \\ T & \to & {{T}^{m+1}} \\ \end{array}\end{eqnarray} \tag{ 44 }$$

into (9). Equation (43) implicitly assumes that $\left\langle \dot{n}(s)|n(s) \right\rangle =0$, which can always be assumed to be true because the phases of $\left| n(s) \right\rangle$ are arbitrary. All of the previous methods can then be used after making these substitutions and using a higher-order polynomial to perform the interpolation.

For example, we can generalize the method of section 5 with the following modification to the conditions required for both f_A and f_B:

$$\begin{eqnarray}\begin{array}{rcl} f_{B}^{(m+1)}(0) & = & -f_{A}^{(m+1)}(0) \\ f_{B}^{(m+1)}(1) & = & -f_{A}^{(m+1)}(1) \\ f_{B}^{(j)}(\Delta ) & = & f_{A}^{(j)}(\Delta )\;\forall \;j\in [0,m+2] \\ f_{B}^{(j)}(1-\Delta ) & = & f_{A}^{(j)}(1-\Delta )\;\forall \;j\in [0,m+2] \\ f_{B}^{(j)}(1) & = & 0\;\forall \;j\in (0,m], \\ \end{array}\end{eqnarray} \tag{ 45 }$$

and picking T_A = T_B with $\theta =\pi /4$. This enables exponentially accurate adiabatic approximations if m is chosen as a function of T.

Alternatively, cancellation of the leading order transition to $O(1/{{T}^{m+1}})$ can be obtained by using the exact same ideas within the methods of sections 3.1 and 3.2 after using the substitutions in (44) and conditions similar to (45) for the polynomial interpolations. It should also be noted that in systems where the adiabatically transported subspace is a one- or two-dimensional subspace, such approaches also can be used to make the overall error scaling $O(1/{{T}^{m+1}})$. This raises the possibility that controlled adiabatic evolution can be combined with boundary cancellation methods to substantially reduce the cost of performing a high-accuracy adiabatic state preparation.

There are two types of costly resources in coherently controlled adiabatic evolutions. The first is the cost of evolving the adiabatic register, denoted $\left| \psi \right\rangle$ in figure 1. The second is the cost of performing the required rotations on the control register. In this section we will provide a complete cost analysis of this model under the assumption that it is being simulated using a circuit-based quantum computer that is further equipped with oracles that compute the necessary properties of the Hamiltonian. We will then conclude that a coherently controlled adiabatic using sparse, row-computable Hamiltonian evolution is polynomially equivalent to the circuit model. Other appropriate cost models, such as bounding the energy and time required to implement the controlled Hamiltonians using k-local Hamiltonians, will not be discussed here. Since we show that coherently controlled adiabatic evolution is polynomially equivalent to the circuit model, it will immediately follow that it is also polynomially equivalent to adiabatic quantum computation using local Hamiltonians.

The first important result that we need to show this is an upper bound on the number of oracle queries needed to simulate a time-dependent Hamiltonian within fixed error on a quantum computer. We will use this result to upper-bound the query complexity of performing the controlled adiabatic evolutions. In order to understand the theorem, we will define a smoothness classification for Hamiltonians:

Definition 1. The set of operators $\{{{H}_{j}}\;:j=1,...,m\}$ is $\Lambda -p$-smooth on $\mathcal{I}\subseteq \mathbb{R}$ if $\Lambda \geqslant {{\left( \sum _{j=1}^{m}\left\|\partial _{t}^{p}{{H}_{j}}(t)\right\| \right)}^{1/(p+1)}}$, for all $t\in \mathcal{I}$ and $p\in \{0,1,\cdots ,P\}$.

Now with definition 1 we can state the following corollary, which gives the query complexity of the simulaton algorithm. The number of one-and two-qubits needed for the simulation is at most proportional to the number of queries made.

Corollary 1. (Cor. 6 of [36]) Let $\{{{H}_{\mu }}\;:\mathbb{R}\mapsto {{\mathbb{C}}^{{{2}^{n}}\times {{2}^{n}}}};\mu =1,\ldots ,M\}$ be a set of time-dependent Hermitian operators that is $\Lambda -2k$-smooth on $\mathcal{I}=({{t}_{0}},{{t}_{0}}+\Delta t)\backslash \{{{t}_{1}},\ldots ,{{t}_{L}}\}$, where ${{t}_{0}}\lt {{t}_{1}}\lt \cdots \lt {{t}_{L}}\lt {{t}_{0}}+\Delta t$, with the additional conditions

• 1.  
  $\exists \;{{H}_{{\rm max} }}\in \mathbb{R}\;:{{H}_{{\rm max} }}\geqslant {{{\rm max} }_{t\in [{{t}_{0}},{{t}_{0}}+\Delta t]}}\parallel H(t)\parallel$,
• 2.  
  $0\lt \epsilon \leqslant {\rm min} \{1,27{{(5/3)}^{k-1}}{{{\rm d}}^{2}}\Lambda \Delta t\}$,
• 3.  
  ${{N}_{T}}$ satisfies ${{N}_{T}}\geqslant \left \lceil {{{\rm log} }_{2}}\left( \frac{({{{\rm max} }_{t\in \mathcal{I},\mu }}\parallel {{\partial }_{t}}{{H}_{\mu }}(t)\parallel )(32kM{{d}^{2}}){{(5/3)}^{k-1}}\Delta {{t}^{2}}}{\epsilon } \right)\right \rceil$,
• 4.  
  ${{n}_{H}}$ satisfies ${{n}_{H}}\geqslant 2\left \lceil {{{\rm log} }_{2}}\left( \frac{32\;kM{{d}^{2}}{{(5/3)}^{k-1}}\Lambda \Delta t}{\epsilon } \right)\right \rceil +6$ and
• 5.  
  $\Delta t/{{2}^{{{N}_{T}}}}\lt {{{\rm min} }_{\ell =0,\ldots ,L}}({{t}_{\ell +1}}-{{t}_{\ell }})$ with ${{t}_{L+1}}:={{t}_{0}}+\Delta t$,

where ${{N}_{T}}$ is the number of bits used to represent t and ${{n}_{H}}$ is the number of qubits used to encode the matrix elements of H. Then the query complexity for simulating evolution generated by $H(t)={{\sum }_{\mu }}T_{\mu }^{\dagger }{{H}_{\mu }}(t){{T}_{\mu }}$, for a fixed set of unitary basis changing operators ${{T}_{\mu }}$, within an error of $\epsilon$ using time-ordered Trotter–Suzuki formulas with error $O(\Delta {{t}^{2k+1}})$ is

$$\begin{eqnarray}&&{{N}_{{\rm queries}}}\leqslant 12CM{{d}^{2}}{{5}^{k-1}}\left[ (L+1)+24k{{d}^{2}}\Lambda \Delta t{{\left( \frac{5}{3} \right)}^{k}}{{\left( \frac{6{{d}^{2}}\Lambda \Delta t}{(\epsilon /3)} \right)}^{1/2k}} \right],\end{eqnarray} \tag{ 46 }$$

where $C$ is the number of oracle calls needed to simulate a one-sparse Hamiltonian, and the number of basis change operations is at most ${{N}_{{\rm queries}}}/(3C{{d}^{2}})$.

Note that we need to use a result for simulating piecewise smooth Hamiltonians because the interpolations used in our method will typically cause $H({{f}_{j}}(s))$ to be non-analytic at either one or two points. The result of corollary 1 is then useful because it provides the cost of performing such a simulation despite such complications. Note that in the cases we consider L = 1 for partially anti-symmetric boundary conditions and L = 2 for completely anti-symmetric boundary conditions. Also, for simplicity, we cite a method that does not use adaptively chosen timesteps. Such adaptive methods are given in [36] and lead to similar scaling where $\Lambda$ is replaced by the time average of the instantaneous values of $\Lambda$.

There are two types of oracles that are required by this corollary. First, an oracle is required that outputs the location of the jth non-zero matrix element in a given row, where $j\leqslant d$ if H is d-sparse.

$$\begin{eqnarray}&&{{O}_{1}}(j)\;:\;\left| i \right\rangle \left| 0 \right\rangle \to \left| i \right\rangle \left| {{L}_{j}}(i) \right\rangle \end{eqnarray} \tag{ 47 }$$

where L_j(i) gives the jth non-zero element in row i. We cost a single query to O₁ as n queries to a single qubit oracle because each query to this oracle yields n bits, and it is more realistic to cost the algorithm by the number of qubits output if the dimension of the Hilbert space is large. The corollary also requires an oracle for matrix elements of H(s)

$$\begin{eqnarray}&&{{O}_{2}}\;:\ \left| i \right\rangle \left| k \right\rangle \left| s \right\rangle \left| 0 \right\rangle \to \left| i \right\rangle \left| k \right\rangle \left| s \right\rangle \left| H{{(s)}_{i,k}} \right\rangle .\end{eqnarray} \tag{ 48 }$$

We also use (for convenience) a new oracle, O_f, whose role is to prepare a quantum state encoding the particular value of f_A(s) or f_B(s) that is needed in a given timestep. In general, if we wish to find the value of f_p(s) we use the oracle in the following way:

$$\begin{eqnarray}&&{{O}_{f}}(s)\left| p \right\rangle \left| 0 \right\rangle =\left| j \right\rangle \left| {{f}_{p}}(s) \right\rangle .\end{eqnarray} \tag{ 49 }$$

This oracle is crucial to our approach because it allows us to remove the multiple controls used in figure 1. For example,

$$\begin{eqnarray*}&&{{O}_{2}}{{O}_{1}}{{O}_{f}}(s)\mathop{\sum }\limits_{p=1}^{N}\mathop{\sum }\limits_{x=1}^{{{2}^{n}}}{{b}_{p}}{{a}_{j}}\left| p \right\rangle \left| x \right\rangle \left| 0 \right\rangle \left| 0 \right\rangle \left| 0 \right\rangle \end{eqnarray*}$$

$$\begin{eqnarray}&&=\mathop{\sum }\limits_{p=1}^{N}\mathop{\sum }\limits_{x=1}^{{{2}^{n}}}{{b}_{p}}{{a}_{j}}\left| p \right\rangle \left| x \right\rangle \left| {{L}_{\ell }}(x) \right\rangle \left| {{f}_{p}}(s) \right\rangle \left| H{{({{f}_{p}}(s))}_{x,L(x,i)}} \right\rangle .\end{eqnarray} \tag{ 50 }$$

Our cost analysis of the controlled adiabatic evolution follows by converting the controlled evolution in figure 1 into the evolution of a single larger Hamiltonian. This larger Hamiltonian can then be simulated by conventional means (such as a Trotter–Suzuki decomposition as per corollary 1).

Theorem 2. Assume that we wish to simulate a coherently controlled adiabatic evolution that uses the controlled evolutions $\{\mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\,H({{f}_{1}}(s)){{T}_{1}}}},\ldots ,\mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\,H({{f}_{p}}(s)){{T}_{p}}}}\}$ such that H(s) is a Hamiltonian satisfying $\parallel H^{\prime} (s)\parallel \leqslant \Gamma$ and each $H({{f}_{j}})$ is Λ – $2k$-smooth and all remaining assumptions of corollary 1 are held for $\Delta t={{{\rm max} }_{j}}{{T}_{j}}$. The query complexity of performing the simulation within error at most using kth-order time-ordered Trotter–Suzuki formulas and oracles that yield one bit per query obeys

$$\begin{eqnarray}&&{{N}_{{\rm queries}}}\leqslant 12CM{{d}^{2}}{{5}^{k-1}}\left[ (L+1)+24k{{d}^{2}}\Lambda {{{\rm max} }_{j}}{{T}_{j}}{{\left( \frac{5}{3} \right)}^{k}}{{\left( \frac{6{{d}^{2}}\Lambda {{{\rm max} }_{j}}{{T}_{j}}}{(\epsilon /6)} \right)}^{1/2k}} \right],\end{eqnarray} \tag{ 51 }$$

where

$$\begin{eqnarray}&&C\leqslant 4n({{z}_{n}}+2)+3{{n}_{H}}+\left \lceil 2{{{\rm log} }_{2}}\left( \frac{6\Gamma {{{\rm max} }_{j}}{{T}_{j}}}{\epsilon } \right)\right \rceil ,\end{eqnarray} \tag{ 52 }$$

z_n is the number of times that $n\mapsto \lceil 2{{{\rm log} }_{2}}\;n\rceil$ must be iterated before achieving a value that is less than or equal to 6 and ${{N}_{f(s)}}\leqslant \lceil {{{\rm log} }_{2}}(\Gamma )\rceil$.

Proof. To see this, note that the controlled unitary evolutions in figure 1 produce a time-evolution operator of the following block-diagonal form:

$$\begin{eqnarray*}\left[ \begin{array}{ccccccccccccccc} {{U}_{1}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {\mathbb{1}} \\ \end{array} \right]\times \cdots \times \left[ \begin{array}{ccccccccccccccc} {\mathbb{1}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {{U}_{p}} \\ \end{array} \right]=\left[ \begin{array}{ccccccccccccccc} {{U}_{1}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {{U}_{p}} \\ \end{array} \right],\end{eqnarray*}$$

where ${{U}_{1}},\ldots ,{{U}_{p}}$ are the p controlled adiabatic evolutions. By expanding out the unitaries as time-ordered operator exponentials, we see that the ideal time evolution operator is of the form

$$\begin{eqnarray}\left[ \begin{array}{ccccccccccccccc} {{U}_{1}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {{U}_{p}} \\ \end{array} \right]=\left[ \begin{array}{ccccccccccccccc} \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\,H({{f}_{1}}(s)){\rm d}s{{T}_{1}}}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\,H({{f}_{p}}(s)){\rm d}s{{T}_{p}}}} \\ \end{array} \right].\end{eqnarray} \tag{ 53 }$$

Consider the Hamiltonian $H={{\sum }_{j}}\left| j \right\rangle \left\langle \;j \right|\otimes {{H}_{j}}$. It is easy to see using Taylor expansion and the fact that each of the terms in H commute that

$$\begin{eqnarray}\begin{array}{rcl} {{{\rm e}}^{-{\rm i}Ht}} & = & {{{\rm e}}^{-{\rm i}{{\sum }_{j}}\left| \,j \right\rangle \left\langle \,j \right|\otimes {{H}_{j}}t}} \\ {} & = & \mathop{\prod }\limits_{j=1}^{p}{{{\rm e}}^{-{\rm i}\left| \,j \right\rangle {{\left\langle \,j \right|}^{{}}}}}\otimes {{H}_{j}}{{t}^{{}}} \\ {} & = & \mathop{\prod }\limits_{j=1}^{p}(\left| j \right\rangle {{\left\langle \;j \right|}^{{}}}\otimes {{{\rm e}}^{-{\rm i}{{H}_{j}}t}}+({\mathbb{1}}-\left| \;j \right\rangle \left\langle \;j \right|)\otimes {\mathbb{1}}) \\ {} & = & \mathop{\oplus }\limits_{j=1}^{p}\left| j \right\rangle {{\left\langle \;j \right|}^{{}}}\otimes {{{\rm e}}^{-{\rm i}{{H}_{j}}t}}. \\ \end{array}\end{eqnarray} \tag{ 54 }$$

Here $\oplus$ represents the direct sum operation. Thus we have that

$$\begin{eqnarray}{{{\rm e}}^{-{\rm i}Ht}}=\left[ \begin{array}{ccccccccccccccc} {{{\rm e}}^{-{\rm i}{{H}_{1}}t}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {{{\rm e}}^{-{\rm i}{{H}_{p}}t}} \\ \end{array} \right].\end{eqnarray} \tag{ 55 }$$

Now let ${\bf H}(s)={{\sum }_{j}}\left| j \right\rangle \left\langle \;j \right|\otimes {{H}_{j}}(s)$. It then follows from the definition of the ordered-operator exponential and the block-diagonal structure of (55) that

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}H(s){\rm d}sT}} & = & \mathop{{\rm lim} }\limits_{r\to \infty }\mathop{\prod }\limits_{j=1}^{r}\left[ \begin{array}{ccccccccccccccc} {{{\rm e}}^{-{\rm i}{{H}_{1}}(j/r)T/r}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & {{{\rm e}}^{-{\rm i}{{H}_{p}}(j/r)T/r}} \\ \end{array} \right] \\ {} & = & \left[ \begin{array}{ccccccccccccccc} \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}{{H}_{1}}(s){\rm d}sT}} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}{{H}_{p}}(s){\rm d}sT}} \\ \end{array} \right]. \\ \end{array}\end{eqnarray} \tag{ 56 }$$

It then follows that the controlled evolutions in (53) can be expressed as a simulation of a single time-dependent Hamiltonian by taking ${{H}_{j}}(s)\to H({{f}_{j}}(s)){{T}_{j}}/T$ in (56). For simplicity, let us take $T={{{\rm max} }_{j}}{{T}_{j}}$.

Next we need to find the properties of the dilated Hamiltonian H(s) that describes the controlled evolutions in the controlled adiabatic evolution. Firstly, assuming each H_j is the sum of M Hamiltonians that can be efficiently transformed into d-sparse matrices, it follows that ${\bf H}$ can be expressed as a similar sum. Similarly, since $H(s)={{\sum }_{j}}\left| \;j \right\rangle \left\langle \;j \right|\otimes H({{f}_{j}}(s)){{T}_{j}}/T$, it follows from the fact that ${\bf H}(s)$ has a direct product structure and the assumption that each $H({{f}_{j}}(s))$ is $\Lambda$ – $2k$-smooth that for any non-negative integer $q\leqslant 2k$

$$\begin{eqnarray}&&\parallel \partial _{s}^{q}{\bf H}(s)\parallel =\mathop{{\rm max} }\limits_{j}\parallel \partial _{s}^{H}({{f}_{j}}(s))\parallel \frac{{{T}_{j}}}{T}\leqslant {{\Lambda }^{q+1}}.\end{eqnarray} \tag{ 57 }$$

Hence, for $T={{{\rm max} }_{j}}{{T}_{j}}$, H(s) is at most $\Lambda -2k$-smooth.

We then have from (46) that the cost of simulating the effective Hamiltonian ${\bf H}(s)$ using kth-order Trotter–Suzuki formulas is at most

$$\begin{eqnarray}&&\begin{array}{l} {{N}_{{\rm queries}}}\leqslant 12CM{{d}^{2}}{{5}^{k-1}} \\ \quad \times \left[ (L+1)+24k{{d}^{2}}\Lambda \mathop{{\rm max} }\limits_{j}{{T}_{j}}{{\left( \frac{5}{3} \right)}^{k}}{{\left( \frac{6{{d}^{2}}\Lambda \mathop{{\rm max} }\limits_{j}{{T}_{j}}}{(\epsilon /3)} \right)}^{1/2k}} \right]. \\ \end{array}\end{eqnarray} \tag{ 58 }$$

The remaining issue is the calculation of C. In order to compute C we need to first show that we can simulate a query to the Hamiltonian oracles for ${\bf H}$ using those for H. We specifically require two oracles: one that computes the locations of the ith (potentially) non-zero matrix element in any row of ${\bf H}$ and another that evaluates that matrix element at a fixed value of s.

The oracle for finding the column index for a specified element in row x of ${\bf H}$ can be constructed as follows. The oracle O₁ has the property that

$$\begin{eqnarray}&&{{O}_{1}}(q)\left| x \right\rangle \left| y(q) \right\rangle ,\end{eqnarray} \tag{ 59 }$$

where y(q) is the column index of the qth element in row x. Then for any j we can construct the corresponding oracle by exploiting the block diagonal structure of ${\bf H}$ via

$$\begin{eqnarray}&&{{{\bf O}}_{1}}(q)\left| j \right\rangle \left| x \right\rangle \left| 0 \right\rangle =\left| j \right\rangle \left| x \right\rangle \left| y(q)+{{2}^{n}}(j-1) \right\rangle =\left| j \right\rangle \left| x \right\rangle \left| j \right\rangle \left| y(q) \right\rangle .\end{eqnarray} \tag{ 60 }$$

The oracle ${{{\bf O}}_{1}}(q)$ can therefore be enacted using one query to O₁ and a polynomial size arithmetic circuit.

The second oracle ${{{\bf O}}_{2}}(q)$ gives for a specific value of s that is specified, the value of H(s). Specifically, after taking into account the block diagonal structure of ${\bf H}$, we need the oracle to be of the form

$$\begin{eqnarray}&&{{{\bf O}}_{2}}\left| j \right\rangle \left| x \right\rangle \left| y \right\rangle \left| s \right\rangle \left| 0 \right\rangle =\left| j \right\rangle \left| x \right\rangle \left| y \right\rangle \left| {{[H({{f}_{j}}(s))]}_{x,y}} \right\rangle .\end{eqnarray} \tag{ 61 }$$

This oracle can be implemented using one query to O_f and one query to O₂.

In [36], it is assumed that the time is provided to the oracles via classical control. Here, we assume that the time is provided via a quantum register so we must add the cost of preparing this register to the cost, C, of simulating a one-sparse matrix. Lemma 9 of [36] gives us that the query complexity (costed at 1/per bit of output yielded by ${{{\bf O}}_{1}}$ and ${{{\bf O}}_{2}}$) is

$$\begin{eqnarray}&&C\leqslant 4n({{z}_{n}}+2)+3{{n}_{H}}.\end{eqnarray} \tag{ 62 }$$

For each one-sparse Hamiltonian that appears in the Trotter–Suzuki decomposition, the time register must be initialized once [36]. This causes an additional source of error and if we are to fit it within our error budget, we must reduce the error in other parts of the simulation algorithm. There are three sources of error in the simulation algorithm: Trotter–Suzuki error, error due to finite n_H and error due to finite n_T (we have neglected errors in synthesizing single qubit operations etc). Each of these three sources of error is chosen to be at most $\epsilon /3$ in [36]. Therefore, if we reduce the error tolerance in the Trotter–Suzuki approximation to $\epsilon /6$ and allow an error tolerance of $\epsilon /6$ for approximating f_j(s) then the overall error will remain at most $\epsilon$. Thus the overall complexity becomes

$$\begin{eqnarray}\begin{array}{rcl} {{N}_{{\rm queries}}} & \leqslant & 12CM{{d}^{2}}{{5}^{k-1}}\left[ (L+1)+24k{{d}^{2}}\Lambda \mathop{{\rm max} }\limits_{j}{{T}_{j}}{{\left( \frac{5}{3} \right)}^{k}}{{\left( \frac{6{{d}^{2}}\Lambda \mathop{{\rm max} }\limits_{j}{{T}_{j}}}{(\epsilon /6)} \right)}^{1/2k}} \right]. \\ \end{array}\end{eqnarray} \tag{ 63 }$$

The error in ${{{\rm e}}^{-{\rm i}{\bf H}(s)T}}$ is at most $\parallel \Delta {\bf H}(s)\parallel T$ [37], where $\Delta {\bf H}(s)$ is the error in implementing the Hamiltonian. By Taylorʼs theorem this is at most $\Gamma {{{\rm max} }_{j}}|\Delta {{f}_{j}}(s)|{{{\rm max} }_{j}}{{T}_{j}}$, where $\Delta {{f}_{j}}(s)$ is the error incurred by approximating f_j(s) to a finite number of digits. Let us define the number of digits used to express f_j as ${{n}_{{{f}_{j}}}}$. Then the error in f_j is $\Delta {{f}_{j}}\leqslant {{2}^{-{{n}_{{{f}_{j}}}}}}$. Hence it suffices to choose

$$\begin{eqnarray}&&{{2}^{-{{n}_{{{f}_{j}}}}}}\Gamma \mathop{{\rm max} }\limits_{j}{{T}_{j}}=\epsilon /6.\end{eqnarray} \tag{ 64 }$$

Thus since we have to both compute the value of f_j(s) to ${{n}_{{{f}_{j}}}}$ bits of precision using queries to O₂ and then uncompute it, equations (64) and (62) give

$$\begin{eqnarray}&&C\leqslant 4n({{z}_{n}}+2)+3{{n}_{H}}+2\left \lceil {{{\rm log} }_{2}}\left( \frac{6\Gamma {{{\rm max} }_{j}}{{T}_{j}}}{\epsilon } \right)\right \rceil ,\end{eqnarray} \tag{ 65 }$$

as claimed.

We therefore see from theorem 2 that this model of adiabatic computation can be efficiently simulated using the posited oracles under reasonable smoothness assumptions. This naturally leads to the following corollary:

Corollary 3. Let ${{f}_{j}}(s):j=1,\ldots ,p$ efficiently computable functions, $H(s)=\sum _{\mu =1}^{M}{{T}_{\mu }}{{H}_{\mu }}(s)T_{\mu }^{\dagger }$ where each ${{H}_{\mu }}({{f}_{j}}(s))$ is a d-sparse row-computable matrix for all s and $p\in O({\rm poly}(n))$. If the conditions of theorem 2 are satisfied for the adiabatic paths $\{{{f}_{1}},\ldots ,{{f}_{p}}\}$ then controlled adiabatic evolution under $\{H({{f}_{1}}(s)),\ldots ,H({{f}_{p}}(s))\}$ is polynomially equivalent to both the circuit model and in turn adiabatic quantum computation.

Proof. We know that a circuit simulation of the controlled adiabatic evolution is efficient under the assumptions of theorem 2 given access to the oracles O₁, O₂ and O_f. If H(s) is row computable, then it implies that there exist efficient algorithms to find the locations and values of each non-zero matrix element of H(s). Thus O₁ and O₂ can be implemented efficiently by the definition of row computability.

O_f can be efficiently computed for each j by assumption, and hence for any fixed j and s the state $|{{f}_{j}}(s)\rangle$ can be prepared efficiently. Furthermore, because $p\in O({\rm poly}(n))$, it follows that the state $\sum _{j=1}^{p}{{a}_{j}}\left| j \right\rangle |{{f}_{j}}(s)\rangle$ can be efficiently prepared. Thus O_f can be efficiently simulated in the circuit model as well. This implies that quantum computers can efficiently simulate this class of coherently controlled adiabatic evolutions.

Local Hamiltonians are a subset of d-sparse Hamiltonians. Therefore the class of adiabatic evolutions considered includes a set of Hamiltonians that generate a family of evolutions that are polynomially equivalent to the circuit model [38]. Thus if we ignore the control register, then the controlled adiabatic evolution can be reduced to a universal adiabatic quantum computer. Thus our model of computation is polynomially equivalent to both the circuit model and adiabatic quantum computation.

We now see that controlled adiabatic quantum computation using piecewise smooth, sparse, bounded, row-computable Hamiltonians is not an exponentially more powerful model of computation than traditional adiabatic computation. Apart from showing that this is a reasonable model of quantum computation, it also shows that the maximum evolution time used is a reasonable metric for the cost of the evolution (once made dimensionless by multiplying by a characteristic energy of the system). For most of the adiabatic paths considered, the contribution of the derivatives of the Hamiltonian to $\Lambda$ is negligible and thus in practice it suffices to ignore their contributions. Also, because this algorithm scales near-linearly with the evolution time, this analysis clearly shows that our model can only potentially provide sub-polynomial speedups over circuit-based quantum computation for fixed d and M.

Our methods improve the performance of adiabatic state preparation by utilizing knowledge about the Hamiltonian. Such knowledge is never perfect, in practice, because of noise and experimental imperfections. The question remains whether coherently controlled adiabatic evolutions inherit the robustness of traditional adiabatic approaches or if the extra knowledge the approach requires causes them to be innately more fragile. For simplicity, we restrict our attention to adiabatic dynamics controlled with one qubit. The generalization to many qubits is straightforward.

Let us first consider errors in the Hamiltonians H_A and H_B applied to the adiabatic register. The effective undisturbed time evolution operator for the system (after absorbing the individual evolution times T_A and T_B into the Hamiltonians) is

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\left| 0 \right\rangle \left\langle 0 \right|\otimes {{H}_{A}}(s){\rm d}sT}}\mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\left| 1 \right\rangle \left\langle 1 \right|\otimes {{H}_{B}}(s){\rm d}sT}} \\ \quad =\mathcal{T}{{{\rm e}}^{-{\rm i}\int _{0}^{1}\left| 0 \right\rangle \left\langle 0 \right|\otimes {{H}_{A}}(s)+\left| 1 \right\rangle \left\langle 1 \right|\otimes {{H}_{B}}(s){\rm d}sT}}. \\ \end{array}\end{eqnarray} \tag{ 66 }$$

Thus the system evolution can be described by a single time-dependent Hamiltonian with eigenvalue gap $\Gamma :={{{\rm min} }_{\nu ,s}}\{\gamma _{\nu ,g}^{A}(s),\gamma _{\nu ,g}^{B}(s)\}$. In particular, let us assume that a perturbation Hamiltonian $\lambda {{H}_{P}}$ is added to the effective Hamiltonian, where ${{H}_{P}}=\left| 0 \right\rangle \left\langle 0 \right|\otimes H_{0}^{P}+\left| 1 \right\rangle \left\langle 1 \right|\otimes H_{1}^{P}$. We take this form to ensure that the imperfection on the adiabatic register does not modify the state of the control register. As the total effective system is gapped, we anticipate that the resultant adiabatic evolution will be robust to these errors if the magnitude of the perturbation is small relative to $\Gamma$.

Our intuition can be validated using perturbation theory to leading order. The shift in the instantaneous eigenvalues for the Hamiltonian is

$$\begin{eqnarray}&&\Delta {{E}_{\nu ,g}}(s)=\lambda \left\langle \nu (s) \right|{{H}_{P}}(s)\left| \nu (s) \right\rangle -\lambda \left\langle g(s) \right|{{H}_{P}}(s)\left| g(s) \right\rangle \leqslant 2\lambda \parallel {{H}_{P}}(s)\parallel +O({{\lambda }^{2}}).\end{eqnarray} \tag{ 67 }$$

Similarly, the shift in the instantaneous eigenvector $\left| g(s) \right\rangle$ is

$$\begin{eqnarray}&&\begin{array}{l} \lambda \mathop{\sum }\limits_{\nu \ne g}\left( \frac{\left\langle 0 \right|\left\langle \nu (s) \right|{{H}^{P}}(s)\left| 0 \right\rangle \left| g(s) \right\rangle }{\gamma _{\nu ,g}^{A}(s)}\left| 0 \right\rangle +\frac{\left\langle 0 \right|\left\langle \nu (s) \right|{{H}^{P}}(s)\left| 1 \right\rangle \left| g(s) \right\rangle }{\gamma _{\nu ,g}^{B}(s)}\left| 1 \right\rangle \right)\left| \nu (s) \right\rangle \\ \quad +O({{\lambda }^{2}}). \\ \end{array}\end{eqnarray} \tag{ 68 }$$

Equation (68) implies that, for fixed dimension, the error in the final state will be small if $\lambda {{{\rm max} }_{s}}\parallel {{H}^{P}}(s)\parallel \ll \Gamma$. This suggests that the gap protects the ground state against this sort of noise exactly in the same manner as in traditional adiabatic evolution. On the other hand, the gap provides no corresponding protection against the shifts in the phase that H^P incurs. We see from (67) that such phases have a negligible impact if

$$\begin{eqnarray}&&\lambda \int _{0}^{1}\parallel {{H}^{P}}(s)\parallel {\rm d}sT\ll 1.\end{eqnarray} \tag{ 69 }$$

Thus controlled adiabatic evolution will function well under such imperfections if $T\ll 1/\lambda$. This suggests that these errors impose a time cutoff, after which the errors in the phases in the two adiabatic paths differ enough for phase cancellation to fail. In such cases, the system will tend to revert to performance that is comparable to uncontrolled adiabatic evolution because (68) shows that the ground state will still be protected against such noise by the gap. This behavior is illustrated in figure 7.

Zoom In Zoom Out Reset image size

Another way to look at this is to ask the question of how much control is needed in order to guarantee a particular magnitude of error given a fixed evolution time. To ensure that errors in the Hamiltonian do not degrade the asymptotic scaling of coherently controlled adiabatic evolution, λ must be of order $O\left( \frac{1}{T} \right)$. We find a similar bound if we do not assume that H_P has a direct sum structure i.e. ${{H}_{P}}\ne \left| 0 \right\rangle \left\langle 0 \right|\otimes H_{0}^{P}+\left| 1 \right\rangle \left\langle 1 \right|\otimes H_{1}^{P}$ using results from [37].

Timing errors are also an important contributor to the errors in controlled adiabatic evolutions such as these . Suppose that T_A and T_B are shifted by ${{\Delta }_{A}}$ and ${{\Delta }_{B}}$ respectively from the values that satisfy (19). Similarly as for an eigenvalue shift, imperfect timing will cause dephasing errors between the evolutions. Expanding (18) in powers of ${{\Delta }_{A}}$ and ${{\Delta }_{B}}$, we find that the timing errors are of order $O\left( \frac{1}{{{T}^{2}}} \right)$ when ${{\Delta }_{a}},{{\Delta }_{B}}\in O\left( \frac{1}{T} \right)$.

Finally there is the issue of state preparation errors for the control register. In order to ensure that the overall diabatic error is $O(\frac{1}{{{T}^{2}}})$, it suffices to ensure that the state preparation error is at most $O(\frac{1}{T})$. This is because the error for each of the individual adiabatic paths is already $O(\frac{1}{T})$, which means that $O(\frac{1}{T})$ errors in the operations on the control register result in errors of $O(\frac{1}{{{T}^{2}}})$ in the ground state. This means that the control accuracy needed in order to achieve this improved scaling is modest.

These results show that our ideas will work well for performing adiabatic state preparation on a quantum computer as well as for sufficiently short evolutions on a controlled adiabatic device even under certain imperfections. Such results are not surprising, as our techniques exploit knowledge about the dynamics to gain an advantage over traditional approaches. We find, however, that the method is as robust to noise as the result of [32], which does not use controlled adiabatic evolution. This shows that the level of control needed does not preclude the use of this approach in non-fault-tolerant applications. It is also interesting to note that these dephasing errors in the instantaneous eigenbasis of the Hamiltonian can actually be corrected by using the control registers to add an appropriate phase to one of the computational braches, whereas such errors can only be corrected by changing the evolution time in [32].

Our work shows that traditional approaches to adiabatic state preparation or quantum computing that only utilize classical control over the adiabatic subsystem are unnecessarily restrictive. In particular, we show that much greater control over the adiabatic process can be achieved by relaxing these conditions to allow an adiabatic system to perform coherent control of a polynomially large number of adiabatic evolutions. Furthermore, by introducing measurements on these registers we gain the ability to apply non-unitary operations on the system. This is significant because it not only allows dissipative dynamics to be naturally included into the state preparation process but also allows post-selection and feedback to be used to generate gates that would be otherwise difficult to engineer. Although this approach requires more complex control, it retains what is perhaps the most significant feature of adiabatic state preparation algorithms: protection against certain types of noise by the spectral gap. Also, the small control register may require far less quantum error correction than would be required to construct a circuit-based quantum computer that could simulate the entire protocol.

Apart from providing a richer paradigm for performing state preparation, we show that coherently controlled adiabatic evolution can provide performance advantages over existing adiabatic strategies for state preparation. Specifically, it allows us to combine the best features of local adiabatic evolution and boundary cancellation methods, which are two optimization strategies that are traditionally at odds with each other. These combined strategies provide better error scaling than any known adiabatic optimization technique for evolution times that are close to the transition between the Landau–Zener regime and the adiabatic regime. We have also shown that using controlled adiabatic evolutions to prepare initial states is also robust to errors. This means that the technique can be realistically applied in state preparation problems outside of quantum computation. We finally provide an explicit quantum simulation algorithm for simulating our protocol, and in turn traditional adiabatic algorithms, that explicitly gives the cost of simulating the controlled adiabatic evolution using a quantum computer and find that this cost scales near-linearly in the evolution time.

There are several natural applications of our method to quantum computing. First, adiabatic state preparation is widely used in quantum chemistry simulation algorithms to prepare an approximation to the ground state of the molecule in cases where naïve approximations, such as the Hartree–Fock approximation, fail [4]. Traditional circuit synthesis methods for preparing the FCI groundstate do not apply in such cases because the ground state is not known a priori. Thus adiabatic state preparation provides one of the few known methods that can provide states that have high fidelity with the FCI ground state, which is needed to simulate chemically important expectation values such as dipole moments. Controlled adiabatic evolution may substantially accelerate the preparation of high-fidelity simulation of such quantities.

Another application of our methods involves simulating quantum field theories. Adiabatic evolution is used as an important step in the algorithm for simulating ${{\phi }^{4}}$ theories by Jordan, Lee and Preskill [39]. A key step in their method is to use adiabatic evolution to transform the basis from that of the free theory to the interacting theory to enable state preparation and measurement. The number of operations needed to simulate the adiabatic evolution is a major contribution to the complexity of the algorithm and so methods like ours could lead to substantially improved methods for simulating field theories.

These results only begin to scratch the surface of what is possible within this paradigm. First, apart from suppressing diabatic errors, our approach removes several important restrictions that make developing computationally useful adiabatic algorithms challenging. In particular, the use of measurement and coherent control provides an important generalization beyond the approach of [40], which used transitions between the ground state and the excited states to enable an exponential speedup that is not known to be achievable within the standard adiabatic paradigm. Second, our approach explicitly uses linear combinations of unitary operations that are nearly unitary. This raises an interesting question of whether truly non-unitary processes will be of use in optimizing adiabatic passage. Progress towards this goal has already been reported in [41]. Also, techniques similar to ours may be of value in phase randomization protocols similar to [42]. Ultimately, these ideas may even lead to more natural ways of performing error correction or suppression in a coherently controlled adiabatic quantum computer. These are just a few examples of the many avenues of research that are opened by this work.

This research was carried out while MK was visiting the Institute for Quantum Computing at the University of Waterloo and MK is grateful for their kind hospitality. We thank M Mosca, I Hen and D Nagaj for valuable comments and feedback. NW acknowledges funding from USARO-DTO, NSERC and CIFAR. MK acknowledges support from QIMABOS, APVV-0646-10 (COQI), IQC, Industry Canada and NSERC.

The approaches used to cancel the first-order transitions for a two-level system can be generalized for larger systems as well. In a case of a three-level system, we must ensure that transitions to both the first and second excited states are $O(1/T)$. This occurs if the following conditions are met:

$$\begin{eqnarray}\begin{array}{rcl} 0 & = & {{{\rm cos} }^{2}}(\theta )\left[ \frac{\left\langle {{{\dot{e}}}_{1}}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{1}}}^{A}\left( 1 \right){{T}_{A}}}-{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}}}\frac{\left\langle \mathop{{{e}_{1}}}\limits^{.}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{1}}}^{A}\left( 0 \right){{T}_{A}}} \right], \\ {} & {} & +{{{\rm sin} }^{2}}(\theta )\left[ \frac{\left\langle \dot{e}_{1}^{\prime }(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{1}}}^{B}\left( 1 \right){{T}_{B}}}-{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi {{T}_{B}}}}\frac{\left\langle \mathop{{{e}_{1}}}\limits^{.}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{1}}}^{B}\left( 0 \right){{T}_{B}}} \right] \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} 0 & = & {{{\rm cos} }^{2}}(\theta )\left[ \frac{\left\langle \mathop{{{e}_{2}}}\limits^{.}(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{2}}}^{A}\left( 1 \right){{T}_{A}}}-{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}}}\frac{\left\langle \mathop{{{e}_{2}}}\limits^{.}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{1}}}^{B}\left( 0 \right){{T}_{A}}} \right] \\ {} & {} & +{{{\rm sin} }^{2}}(\theta )\left[ \frac{\left\langle \dot{e}_{2}^{\prime }(1)|g(1) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{2}}}^{B}\left( 1 \right){{T}_{B}}}-{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi {{T}_{B}}}}\frac{\left\langle \dot{e}_{2}^{\prime }(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{2}}}^{B}\left( 0 \right){{T}_{B}}} \right], \\ \end{array}.\end{eqnarray} \tag{ A.2 }$$

where the the first evolution corresponds to a Hamiltonian $H({{f}^{A}}(s))$ and the second one to $H({{f}^{B}}(s))$ with its states denoted by primes, parameterizing the Hamiltonian with a single function as in the last paragraph. Moreover, we assume that $H({{f}^{A}}(s))$ and $H({{f}^{B}}(s))$ are equal at the beginning and the end of evolution (but their derivatives with respect to s differ).

It is straightforward to cancel transitions at certain (discrete) times using our knowledge from the 2-level case when we realize

$$\begin{eqnarray}&&\left\langle {{{\dot{e}}}_{1}}(0)|g(0) \right\rangle =\left\langle {{e}_{1}}(0)\left| \dot{f}(s){{|}_{s=0}}{{H}_{1}} \right|g(0) \right\rangle \end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&\left\langle {{{\dot{e}}}_{2}}(0)|g(0) \right\rangle =\left\langle {{e}_{2}}(0)\left| \dot{f}(s){{|}_{s=0}}{{H}_{1}} \right|g(0) \right\rangle \end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&\left\langle {{{\dot{e}}}_{1}}(1)|g(1) \right\rangle =\left\langle {{e}_{1}}(1)\left| -\dot{f}(s){{|}_{s=1}}{{H}_{0}} \right|g(1) \right\rangle \end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&\left\langle {{{\dot{e}}}_{2}}(1)|g(1) \right\rangle =\left\langle {{e}_{2}}(1)\left| -\dot{f}(s){{|}_{s=1}}{{H}_{0}} \right|g(1) \right\rangle .\end{eqnarray} \tag{ A.6 }$$

Hence, by choosing ${{f}^{A}}(s)$ and ${{f}^{B}}(s)$ as in section 3.1 and using (25), we get rid of terms containing derivatives at the end for both levels. This approach trivially generalizes to higher-dimensional systems.

Now we need to fix T_A and T_B in order to remove the contributions from the boundary by requiring that evolutions gain opposite phases. We can rewrite already simplified (A.1), (A.2) as

$$\begin{eqnarray}\left( \begin{array}{ccccccccccccccc} \int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}(\xi ){\rm d}\xi & \int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi \\ \int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi & \int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{T}_{A}} \\ {{T}_{B}} \\ {} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} (2n+1)\pi \\ (2m+1)\pi \\ \end{array} \right)\end{eqnarray} \tag{ A.7 }$$

This system of equation has a solution, unless the determinant of the matrix equals zero. Note that with this approach we get only a discrete set of times T_A and T_B for which the error vanishes, in contrast to many of our prior methods.

Error suppression can also be achieved for arbitrary time if we use more than two evolutions. A 2-level inspired solution uses 4 unitaries, U_A, U_B, U_C and U_D, where U_A and U_C are given by Hamiltonian $H({{f}^{A}}(s))$ and U_B and U_D by $H({{f}^{B}}(s))$. We pick the functions ${{f}^{A}},{{f}^{B}}$ based on section 3.1. The goal is then to find times T_A–T_D and weights a–d such that the following equations are satisfied:

$$\begin{eqnarray}\begin{array}{rcl} 0 & = & \frac{\left\langle {{{\dot{e}}}_{1}}(1)|g(1) \right\rangle }{\gamma _{g,{{e}_{1}}}^{A}\left( 1 \right)}\left[ \frac{a}{{{T}_{A}}}-\frac{b}{{{T}_{B}}}+\frac{c}{{{T}_{C}}}-\frac{d}{{{T}_{D}}} \right] \\ {} & {} & +\frac{\left\langle \mathop{{{e}_{1}}}\limits^{.}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{1}}}^{A}\left( 0 \right)}\left[ \frac{a}{{{T}_{A}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}}}+\frac{b}{{{T}_{B}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi {{T}_{B}}}} \right. \\ {} & {} & +\left. \frac{c}{{{T}_{C}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}(\xi ){\rm d}\xi {{T}_{C}}}}+\frac{d}{{{T}_{D}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi {{T}_{D}}}} \right] \\ \end{array}\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}\begin{array}{rcl} 0 & = & \frac{\left\langle {{{\dot{e}}}_{2}}(1)|g(1) \right\rangle }{\gamma _{g,{{e}_{2}}}^{A}\left( 1 \right)}\left[ \frac{a}{{{T}_{A}}}-\frac{b}{{{T}_{B}}}+\frac{c}{{{T}_{C}}}-\frac{d}{{{T}_{D}}} \right] \\ {} & {} & +\frac{\left\langle \mathop{{{e}_{2}}}\limits^{.}(0)|g(0) \right\rangle }{-{\rm i}\gamma _{g,{{e}_{2}}}^{A}\left( 0 \right)}\left[ \frac{a}{{{T}_{A}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}}}+\frac{b}{{{T}_{B}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi {{T}_{B}}}} \right. \\ {} & {} & +\left. \frac{c}{{{T}_{C}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{C}}}}+\frac{d}{{{T}_{D}}}{{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi {{T}_{D}}}} \right] \\ \end{array}\end{eqnarray} \tag{ A.9 }$$

First, the normalization condition

$$\begin{eqnarray}&&a+b+c+d=1\end{eqnarray} \tag{ A.10 }$$

must hold. Second, we choose b and T_B such that they cancel the error on the first level from evolution by U_A. This is exactly the same problem we solved for a 2-level system, hence the proper b and T_B are

$$\begin{eqnarray}&&{{T}_{B}}=\frac{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi {{T}_{A}}+\left( 2n+1 \right)\pi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}\left( \xi \right){\rm d}\xi }\end{eqnarray} \tag{ A.11 }$$

$$\begin{eqnarray}&&b=\frac{a{{T}_{B}}}{{{T}_{A}}}.\end{eqnarray} \tag{ A.12 }$$

The same can be done for U_C and U_D:

$$\begin{eqnarray}&&{{T}_{D}}=\frac{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi {{T}_{C}}+\left( 2n+1 \right)\pi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}\left( \xi \right){\rm d}\xi }\end{eqnarray} \tag{ A.13 }$$

$$\begin{eqnarray}&&d=\frac{c{{T}_{D}}}{{{T}_{C}}}.\end{eqnarray} \tag{ A.14 }$$

This suppresses the first transition out of $\left| g \right\rangle$ (typically the transition to the first excited state) to $O(1/{{T}^{2}})$. In addition, the errrors from the derivatives at the end on the second excited state cancel as well.

Finally, we are left with

$$\begin{eqnarray}\begin{array}{rcl} 0 & = & \frac{a}{{{T}_{A}}}\left[ {{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}}}+{{{\rm e}}^{\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }}}(\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}+(2n+1)\pi ) \right] \\ {} & {} & +\frac{c}{{{T}_{C}}}\left[ {{{\rm e}}^{{\rm i}\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{C}}}}+{{{\rm e}}^{\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }(\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}(\xi ){\rm d}\xi {{T}_{C}}+(2m+1)\pi )}} \right]. \\ \end{array}\end{eqnarray} \tag{ A.15 }$$

Therefore we can set the value of T_C and ratio of a and c and we are still free to choose arbitrary T_A. After some algebra we can rewrite (A.15)

$$\begin{eqnarray}\begin{array}{l} {} & {} & \frac{a}{{{T}_{A}}}{{{\rm e}}^{\frac{{\rm i}}{2}\left( \int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{A}}+\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }\left( \int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi {{T}_{A}}+\left( 2n+1 \right)\pi \right) \right)}} \\ {} & {} & \qquad \times {\rm cos} \left[ \frac{1}{2}\left( \int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{A}} \right. \right. \\ {} & {} & \quad +\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }\left. \left. \left( \int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi {{T}_{A}}+\left( 2n+1 \right)\pi \right) \right) \right] \\ {} & {} & =-\frac{c}{{{T}_{C}}}{{{\rm e}}^{\frac{{\rm i}}{2}\left( \int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{C}}+\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }\left( \int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi {{T}_{C}}+\left( 2\,m+1 \right)\pi \right) \right)}} \\ {} & {} & \qquad \times {\rm cos} \left[ \frac{1}{2}\left( \int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi {{T}_{C}} \right. \right. \\ {} & {} & \quad +\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }\left. \left. \left( \int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi {{T}_{C}}+\left( 2\;m+1 \right)\pi \right) \right) \right]. \\ \end{array}\end{eqnarray} \tag{ A.16 }$$

We can ensure that both sides of the equation pick the same phases (up to $2k\pi$) by setting T_C and then we only need a value of k for which the cosines would have the same sign. That is possible, unless

$$\begin{eqnarray}&&\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{A}(\xi ){\rm d}\xi +\frac{\int _{0}^{1}\gamma _{g,{{e}_{2}}}^{B}(\xi ){\rm d}\xi }{\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{B}(\xi ){\rm d}\xi }\int _{0}^{1}\gamma _{g,{{e}_{1}}}^{A}\left( \xi \right){\rm d}\xi \equiv 0\;{\rm mod} \;2\pi .\end{eqnarray} \tag{ A.17 }$$

This procedure can also be used to find paths that cancel multiple transitions in higher-dimensional systems, but a closed form may not necessarily exist, unlike the present case.