Quasi-adiabatic Grover search via the WKB approximation

TL;DR

This paper introduces a quasi-adiabatic WKB approximation for quantum dynamics at intermediate times, effectively capturing subtle features missed by traditional adiabatic methods, especially in Grover search problems.

Contribution

It develops a novel WKB-based quasi-adiabatic approximation applicable to intermediate evolution times in quantum systems, improving analysis of quantum search algorithms.

Findings

The method captures subtle dynamical features missed by adiabatic approximation.

It reproduces the quadratic Grover speedup with the optimal schedule.

Sensitivity to the interpolation schedule can lead to nonsensical results.

Abstract

In various applications one is interested in quantum dynamics at intermediate evolution times, for which the adiabatic approximation is inadequate. Here we develop a quasi-adiabatic approximation based on the WKB method, designed to work for such intermediate evolution times. We apply it to the problem of a single qubit in a time-varying magnetic field, and to the Hamiltonian Grover search problem, and show that already at first order, the quasi-adiabatic WKB captures subtle features of the dynamics that are missed by the adiabatic approximation. However, we also find that the method is sensitive to the type of interpolation schedule used in the Grover problem, and can give rise to nonsensical results for the wrong schedule, due to loss of normalization. Conversely, it reproduces the quadratic Grover speedup when the well-known optimal schedule is used.

Figures

Equations

$$y(r)\sim e^{\theta (r)/ϵ}[y_0(r)+ϵy_1(r)+{ϵ}^2{y}_2(r)+\dots ],$$

$$y$$

$$y^{\prime }$$

$$y^{\prime \prime }$$

$$H[r(s)]=[1-r(s)]H_{\mathrm{initial}}+r(s)H_{\mathrm{final}}.$$

$$i\frac{d}{dt}|\chi (t)⟩=\mu H(r[s(t)])|\chi (t)⟩,$$

$$i\frac{d}{ds}|\chi (s)⟩=\mu t_{f}H[r(s)]|\chi (s)⟩.$$

$$iϵ\frac{d}{dr}|\chi (r)⟩=g(r)H(r)|\chi (r)⟩,$$

$$ϵ\equiv \frac{1}{\mu t_f},$$

$$i\frac{\hslash }{t_f}\frac{d}{ds}|\chi ⟩=\left(-\frac{{\hslash }^2}{2m}{\partial }_x^2+V(x,s{t}_f)\right)|\chi ⟩,$$

$$H_{\mathrm{Grover}}[r(s)]=[1-r(s)](I-|u⟩⟨u|)+r(s)(I-|m⟩⟨m|),$$

$$|u⟩\equiv \frac{1}{\sqrt{2^n}}\sum _{x\in \{0,1{\}}^n}|x⟩,$$

$$|m^{\perp }⟩\equiv \frac{1}{\sqrt{K}}\sum _{\begin{array}{c}x\in \{0,1{\}}^n\\ x\ne m\end{array}}|x⟩,$$

$$H(s)=\left(\begin{array}{cc}[1-r(s)]\frac{K}{K+1}& -[1-r(s)]\frac{\sqrt{K}}{K+1}\\ -[1-r(s)]\frac{\sqrt{K}}{K+1}& 1-[1-r(s)]\frac{K}{K+1}\end{array}\right).$$

$$|\chi (s)⟩=\psi (s)|m⟩+\varphi (s)|m^{\perp }⟩,$$

$$iϵ{\psi }^{\prime }$$

$$iϵ{\varphi }^{\prime }$$

$$\psi (0)=\frac{1}{\sqrt{K+1}},\varphi (0)=\sqrt{\frac{K}{K+1}},$$

$${ϵ}^2(1-r){\psi }^{\prime \prime }+ϵ(ϵa_{1,1}+a_{1,2}){\psi }^{\prime }+a_0\psi =0$$

$${ϵ}^2(1-r){\varphi }^{\prime \prime }+ϵ(ϵa_{1,1}+a_{1,2}){\varphi }^{\prime }+(a_0+igϵ)\varphi =0,$$

$$a_0$$

$$a_{1,1}$$

$$r^{\prime }(s)=c_{\alpha }\Delta (r{)}^{\alpha },$$

$$\Delta (r)=\sqrt{1-\frac{4Kr(1-r)}{K+1}}.$$

$$g_0(r)$$

$$g_1(r)$$

$$g_2(r)$$

$$g_3(r)$$

$$(1-r)z_j^{(2)}+a_{1,1}{z}_{j-1}^{(1)}+a_{1,2}{z}_j^{(1)}+a_0{y}_j=0$$

$$(1-r)z_j^{(2)}+a_{1,1}{z}_{j-1}^{(1)}+a_{1,2}{z}_j^{(1)}+a_0{y}_j+ig{y}_{j-1}=0,$$

$$({\theta }^{\prime }{)}^2+ig{\theta }^{\prime }-g^2\frac{Kr(1-r)}{K+1}=0,$$

$${\theta }_{\pm }^{\prime }=-\frac{ig}{2}\left[1\pm \Delta (r)\right].$$

$$\frac{{\psi }_0^{\prime }}{{\psi }_0}$$

$$\frac{{\varphi }_0^{\prime }}{{\varphi }_0}$$

$$\frac{{\psi }_0^{\prime }}{{\psi }_0}$$

Full text

I Introduction

There exist only a handful of Hamiltonian-based quantum algorithms Farhi and Gutmann (1998); Farhi et al. (2000), designed to run on analog quantum computers Feynman (1985); Peres (1985); Margolus (1990), that exhibit a provable quantum speedup Albash and Lidar (2016). The adiabatic version of the Grover search problem is one such example Roland and Cerf (2002). The existence of this speedup is proven using the adiabatic theorem Jansen et al. (2007), i.e., it is based on an asymptotic analysis in the total evolution time. This is in contrast to the circuit model version of the Grover problem Grover (1996), for which a closed-form analytical solution is known for arbitrary evolution times and arbitrary initial amplitude distributions Biham et al. (2000, 1999). No such closed form analytical solution of the Hamiltonian version of Grover’s algorithm is known as of yet.

The Wentzel-Kramers-Brillouin (WKB) method is a famous technique for approximating differential equations which has found applications in many domains of physics and mathematics, including optics, acoustics, astrophysics, elasticity, and quantum mechanics (see Ref. Schlissel (1977) for a mathematical history of the WKB method). In this work, we adopt the WKB method to provide approximate analytical solutions to the Hamiltonian Grover problem. The WKB method we use is quasi-adiabatic (as opposed to semiclassical A. Messiah (1999)): the small parameter is the inverse of the total evolution time (not $\hbar$, which we set to $1$). We choose to focus on the Grover problem since this problem is well studied and understood, but the WKB method is widely applicable and easily generalizable to other Hamiltonian-based quantum algorithms. We thus expect it to be a useful tool in analyzing such algorithms beyond the adiabatic approximation.

We compare the results of the WKB approximation with a numerically exact solution. Strikingly, we find that the quality of WKB results depends strongly on the interpolation schedule from the initial to the final Hamiltonian. The WKB approximation is reliable already at low order for the schedule that generates a quantum speedup for the Grover problem Roland and Cerf (2002), but fails for the other schedules we tested. These other schedules are characterized by a different dependence on the power of the inverse spectral gap.

The structure of the paper is as follows. We briefly review the quasi-adiabatic WKB method in Sec. II. The method is applied to the Grover problem in Sec. III, and the WKB solutions are derived in Sec. IV. The results are discussed and analyzed in Sec. V, where we perform a comparison with the numerically exact solution. We conclude in Sec. VI.

We remark that there are other tools available to study quasi-adiabatic dynamics: Adiabatic perturbation theory is a popular method Teufel (2003). In Appendix A we study a particular variant of adiabatic perturbation theory from Ref. Hagedorn and Joye (2002) and compare it to our method. Our method is not a variant of adiabatic perturbation theory because at the lowest order we do not recover the adiabatic evolution.

II Quasi-adiabatic WKB for interpolating Hamiltonians

We start by briefly reviewing the asymptotic WKB expansion technique (for background see, e.g., Ref. Holmes (2012)), and connect it to interpolating Hamiltonians of the type used in adiabatic quantum computing.

II.1 WKB as an asymptotic expansion

The WKB expansion

[TABLE]

is an ansatz used for the solution of ordinary differential equations in $y(r)$ that contain a small parameter, $\epsilon$, multiplying the highest derivative. This ansatz is an asymptotic expansion in $\epsilon$, i.e., there is no guarantee that it will provide a unique or even a convergent solution. In fact, the asymptotic series for $y(r)$ is usually divergent; the general term $\epsilon^{n}y_{n}(r)$ starts to increase after a certain value $n=n_{\max}$, which can be estimated for second order differential equations of the form $\epsilon^{2}y^{\prime\prime}(r)=Q(r)y(r)$, if $Q(r)$ is analytic Winitzki (2005). The number $n_{\max}$ can be interpreted as the number of oscillations between $r_{0}$ [the point at which $y(r)$ needs to be evaluated] and the turning point $r_{\star}$ [i.e., where $Q(r_{\star})=0]$ closest to $r_{0}$. In this work we will only be concerned with the expansion up to order $\epsilon$ for a second order differential equation. For later convenience, we list the expressions for the derivatives of the ansatz:

[TABLE]

with $y_{k}\equiv 0$ if $k<0$, and where the number of primes denotes the order of the derivative.

II.2 Interpolating Hamiltonians

We consider interpolating Hamiltonians of the form

[TABLE]

which depend on time only via the dimensionless time $s\equiv t/t_{f}$. Here $t_{f}$ denotes the total evolution time and is the only timescale in the problem. The “interpolation schedule” $r(s)$ is strictly increasing, differentiable, and satisfies the boundary conditions $r(0)=0$ and $r(1)=1$. The derivative of the inverse of $r(s)$, viz. $s^{\prime}(r)$, is therefore also strictly positive. This allows us to divide by $s^{\prime}$ when we need to.

Consider now the Schrödinger equation for this evolution

[TABLE]

where $\mu$ is an energy scale, and $H(\cdot)$ is dimensionless, e.g., a linear combination of Pauli matrices. Writing everything in terms of $s$, we get

[TABLE]

One can also write the problem in terms of $r$. This yields:

[TABLE]

where $g(r)\equiv s^{\prime}(r)$, $s(r):[0,1]\mapsto[0,1]$, and where

[TABLE]

is the dimensionless small parameter for our WKB expansion. Since $\epsilon$ is small for large $t_{f}$, we call our method “quasi-adiabatic WKB”.111We remark that the quasi-adiabatic WKB approximation should not be confused with the traditional WKB approximation associated with the $\hbar\to 0$ limit. The latter is typically used as a semiclassical approximation in one-dimensional position-momentum quantum mechanics, involving a potential barrier (see, e.g., Ref. A. Messiah (1999)). The quasi-adiabatic and semiclassical WKB approximations are not interchangeable. This can be seen from the Schrödinger equation for a particle in a one-dimensional potential:

$i\frac{\hbar}{t_{f}}\frac{d}{ds}\ket{\chi}=\left(-\frac{\hbar^{2}}{2m}\partial_{x}^{2}+V(x,st_{f})\right)\ket{\chi},$

where again $s=t/t_{f}$ and $V(x,t)$ is a space- and time-dependent potential energy function. It is evident that there is no way to trade both $\hbar$ and $1/t_{f}$ for a single small parameter, since they appear together as the product $\hbar t_{f}$.

III The Grover problem via the quasi-adiabatic WKB approximation

Recall that the Grover problem can be formulated as finding an item in an unsorted list of $N=2^{n}$ items, in the smallest number of queries Grover (1997). This admits a quadratic quantum speedup, as was first shown by Grover in the circuit model Grover (1996). It is also one of the few instances where an adiabatic algorithm was discovered which recovers the quantum speedup. The crucial insight, which eluded the first attempt Farhi et al. (2000), was that the speedup obtained in the circuit model could be recovered in the adiabatic model provided the right interpolation schedule $r(s)$ is chosen, namely, a schedule that drives the system more slowly when the gap is smaller Roland and Cerf (2002) (see also Refs. Jansen et al. (2007); Rezakhani et al. (2010)).

In the Hamiltonian Grover algorithm one uses the $n$-qubit interpolating Hamiltonian

[TABLE]

where $m\in\{0,1\}^{n}$ is the marked state and

[TABLE]

is the uniform superposition state. The system is initialized in the state $\ket{u}$. It can be easily checked that the dynamics described by this Hamiltonian is restricted to $\mathcal{S}\equiv\mathrm{Span}\{\ket{u},\ket{m}\}$. Let $K\equiv 2^{n}-1$ and define

[TABLE]

so that $\ket{u}=(\ket{m}+\sqrt{K}\ket{m^{\perp}})/\sqrt{K+1}$. Note that $\{\ket{m},\ket{m^{\perp}}\}$ is an orthonormal basis for $\mathcal{S}$. In this basis, the Hamiltonian is

[TABLE]

Henceforth we restrict our analysis to this two-dimensional Hamiltonian and will not return to the high-dimensional Hamiltonian that gave rise to it.

Let

[TABLE]

i.e., henceforth $\psi(s)$ is the amplitude of the marked state (ground state of the final Hamiltonian), and $\phi(s)$ is the amplitude of the unmarked component (the excited state of the final Hamiltonian).

From Eq. (6), the Schrödinger equation for a general interpolation becomes

[TABLE]

The boundary conditions are

[TABLE]

so it follows from Eqs. (13) that $\psi^{\prime}(0)=\phi^{\prime}(0)=0$.

We now turn the above coupled first order system into two decoupled second order differential equations:

[TABLE]

where

[TABLE]

The function $g(r)=s^{\prime}(r)$ uniquely determines the schedule $r(s)$. We shall consider four different schedules corresponding to choices $\alpha\in\{0,1,2,3\}$ in

[TABLE]

where $c_{\alpha}$ is a constant that depends on $K$ (see Refs. Jansen et al. (2007); Roland and Cerf (2002)) and $\Delta(r)$ is the eigenvalue gap of the Hamiltonian in Eq. (11), given by:

[TABLE]

Equation (17) forces the schedule to become slower (faster) when the gap is smaller (larger).

The linear schedule [$r(s)=s$] corresponds to the choice $\alpha=0$, and the schedule discovered by Roland and Cerf Roland and Cerf (2002) corresponds to $\alpha=2$. We also analyze schedules corresponding to $\alpha=1,3$. To find the constant $c_{\alpha}$, we integrate Eq. (17) and use the boundary condition $s(1)=1$. The expressions for the schedules thus obtained, expressed in terms of the corresponding $g_{\alpha}(r)$ functions [recall that $g(r)\equiv s^{\prime}(r)$], are as follows:

[TABLE]

We now turn to the construction of the WKB solutions for both amplitudes $\psi$ and $\phi$ for each of the schedules.

IV Constructing the WKB solutions

To derive the WKB solutions, we substitute the WKB ansatz [Eqs. (2)] into Eqs. (15a) and (15b).222A Mathematica® notebook containing code for obtaining the WKB expressions used in our analysis is provided at https://tinyurl.com/WKB-notebook. Then, we set the terms multiplying different orders $\epsilon^{j}$ to zero, which yields the following recursive set of equations for $j\geq 1$:

[TABLE]

where $\psi$ [Eq. (15a)] is reconstructed from Eq. (20a), and $\phi$ [Eq. (15b)] is reconstructed from Eq. (20b). We consider only the lowest three orders in $\epsilon$ below.

First, isolating the $\epsilon^{0}$ term [i.e., setting $j=0$ in both Eqs. (20a) and (20b)], we obtain the eikonal equation:

[TABLE]

which is a quadratic equation in $\theta^{\prime}$, so that:

[TABLE]

Turning to the $\epsilon^{1}$ term, we obtain the transport equations:

[TABLE]

Here $\psi_{0},\phi_{0}$ are the parts of the WKB approximant that correspond to $y_{0}$, which was a general placeholder for the lowest order term. Further, Eq. (23a) is obtained from Eq. (20a), and Eq. (23b) is obtained from Eq. (20b), both after setting $j=1$ and using the eikonal equation (21) to eliminate the $y_{1}$ term. Let $\Theta_{\pm}\equiv\theta_{\pm}^{\prime}/g(r)$. It is easy to check that the transport equations then become:

[TABLE]

Since, by Eq. (22), $\Theta_{\pm}=-\frac{i}{2}\left[1\pm\Delta(r)\right]$ is independent of $g$, it follows that $\psi_{0}$ and $\phi_{0}$ do not depend on the interpolation $g$.

Further, using $\Theta_{+}=-(i+\Theta_{-})$ and $\Theta^{\prime}_{\pm}=\mp\frac{i}{2}\Delta^{\prime}$, it is straightforward to show that the r.h.s. of Eq. (24a) for $\psi_{0}^{\pm}$ is identical to the r.h.s. of Eq. (24b) for $\phi_{0}^{\mp}$. Thus, after integration we have $\psi_{0}^{\pm}(r)=c_{0}^{\pm}\phi_{0}^{\mp}(r)$, where $c_{0}^{\pm}$ is (the exponential of) an integration constant.

Moreover, the r.h.s. of Eq. (24a) corresponding to $\Theta_{\pm}$ is easily seen to be equal to $-\frac{1}{2}\left[\frac{\Delta^{\prime}(r)}{\Delta(r)}+\frac{1}{1-r}\pm\frac{1}{(1-r)\Delta(r)}\right]$. Hence, integrating Eqs. (24) yields:

[TABLE]

where $\tilde{c}_{0}^{\pm}$ and $\tilde{d}_{0}^{\pm}$ are integration constants. Or, using the explicit form for the gap given in Eq. (18):

[TABLE]

where $c_{0}^{\pm}=e^{\tilde{c}_{0}^{\pm}}$ and $d_{0}^{\pm}=e^{\tilde{d}_{0}^{\pm}}$.

Finally, turning to the $\epsilon^{2}$ term [i.e., setting $j=2$ in Eqs. (20a) and (20b)] yields:

[TABLE]

where $w\equiv\frac{y_{1}}{y_{0}}$. Here $y$ represents both $\psi$ and $\phi$. We used the eikonal equation to eliminate the $y_{2}$ term, and the transport equations to obtain $y_{0}^{\prime}/y_{0}$ in $w^{\prime}$. Solving Eq. (27) yields $y_{1}$. Note that here we cannot remove the dependence of $y_{1}$ on the interpolation $g$.

We can now assemble the different functions into a solution. Given the interpolation $g$, we can integrate Eq. (22) to find $\theta_{\pm}$, resulting in two solutions $\psi^{\pm}$ and $\phi^{\pm}$. This means that we have to consider linear combinations of these two solutions. Thus

[TABLE]

where the constants $A_{\psi,\phi},B_{\psi,\phi}$ are determined using the boundary conditions $\psi(0)=\frac{1}{\sqrt{K+1}}$, $\phi(0)=\sqrt{\frac{K}{K+1}}$, and $\psi^{\prime}(0)=\phi^{\prime}(0)=0$. Note that despite the fact that $\psi_{0}$ and $\phi_{0}$ do not depend on $g$, the parameter $\theta$ does, via $\theta_{\pm}=\int g\Theta_{\pm}dr$. Therefore even at the lowest order, the approximate solution retains a dependence on the interpolation $g$.

The only constraints our solutions must satisfy are the differential equations (20) and the boundary conditions. Thus we are free to choose the integration constants ($c_{0}^{\pm},d_{0}^{\pm}$, and others that would arise at higher orders $j\geq 2$), and henceforth we choose them to be equal at all orders, such that only $A,B,C,D$ are undetermined until we use the boundary conditions.

It is important to remember that the WKB approximation method does not enforce normalization. Hence, generically, the WKB approximation to a quantum state is unnormalized, resulting in approximations to probabilities that may be greater than $1$.333For convenience, we will abuse terminology somewhat and refer to the WKB approximations to physical probabilities as “probabilities” even though they may not be normalized. It should be clear from the context whether we are referring to approximated probabilities or to actual probabilities. Thus, care must be taken when applying this approximation technique to estimate physical quantities, and in particular one must check that normalization holds. For some of the examples we study here, such nonsensical probabilities indeed arise. In Sec. V.2, we study whether the norm of the WKB approximation is an indicator of approximation quality and also whether renormalization can be improve the WKB approximation.

One final general comment is in order. It turns out that the differential equations (15a) and (15b) have the following unfortunate property: substituting the WKB approximation to $\psi$ into Eq. (13a) and solving for $\phi$ does not yield a good approximation to $\phi$. On other hand, the WKB approximation to $\phi$ does yield a good approximation to $\phi$. This is why we need to perform the WKB approximation separately for each of the amplitudes.

V Results

In this section we analyze the quality of the approximate solutions by comparing them with the solutions obtained via numerical integration of the Schrödinger equation. The results obtained by numerical integration are sufficiently accurate that we can take the numerical solution to be a good proxy to the exact solution. We denote the numerically obtained solution by $\ket{\chi_{\mathrm{Num}}}$ and the solution obtained from the WKB approximation by $\ket{\chi_{\mathrm{WKB}}}$.

V.1 Single Qubit in a magnetic field

As a simple test, we first apply the formalism developed in Sec. II to the case $K=n=1$, which models a qubit in a time-varying magnetic field that changes from the $x$-direction to the $z$-direction, with a linear interpolation $r(s)=s$:

[TABLE]

where $\sigma^{x}\equiv\ket{m}\bra{m^{\perp}}+\ket{m^{\perp}}\bra{m}$ and $\sigma^{z}\equiv\ket{m^{\perp}}\bra{m^{\perp}}-\ket{m}\bra{m}$. Thus, the eikonal equation (22) becomes

[TABLE]

where $\Delta(r)\equiv\sqrt{1-2r(1-r)}$. Therefore the two energy levels of this problems are $E_{\pm}(s)=-i\theta_{\pm}^{\prime}$. Similarly, the solutions [Eqs. (26)] of the transport equations yield, after setting $K=1$,

[TABLE]

where we have chosen the integration constants to remove overall numerical factors.

Next, we may use these solutions to obtain the first-order correction. For this we obtain from Eq. (27):

[TABLE]

From these expressions and the boundary conditions we construct two solutions: $\ket{\chi_{\mathrm{WKB}}^{(0)}}$ (using $\psi_{0},\phi_{0}$) and $\ket{\chi_{\mathrm{WKB}}^{(1)}}$ (using $\psi_{0},\phi_{0}$ and $\psi_{1},\phi_{1}$). We expect $\ket{\chi_{\mathrm{WKB}}^{(1)}}$ to be a better approximation to the exact solution than $\ket{\chi_{\mathrm{WKB}}^{(0)}}$ and we expect the quality of approximation to improve with increasing $t_{f}$, i.e., with decreasing $\epsilon$. We also consider the naive adiabatic approximation, which we define as the instantaneous ground state of $H(r)$.

Figure 1 shows that the WKB approximation is able to capture the correct population dynamics.444Figures were made with the help of the MaTeX package for Mathematica® by Szabolcs Horvát (see url: http://szhorvat.net/pelican/latex-typesetting-in-mathematica.html). In more detail, Fig. 1 shows that the approximation captures oscillations not present in a naive adiabatic approximation, and Fig. 1 shows that the quality of the approximation improves from the lowest order to the next order of the WKB approximation.

Next, consider the final ground state probability, $p_{\mathrm{GS}}(t_{f})$. In Fig. 2, we see that $\ket{\chi_{\mathrm{WKB}}^{(0)}}$ is already sufficient to capture the asymptotic scaling of $p_{\mathrm{GS}}$ with $t_{f}$. Further, $\ket{\chi_{\mathrm{WKB}}^{(1)}}$ captures the oscillations in $p_{\mathrm{GS}}(t_{f})$, with an accuracy that grows with increasing $t_{f}$. Performing a series expansion of $\left|\braket{1}{\chi_{\mathrm{WKB}}^{(0)}(1)}\right|^{2}$ in powers of $\frac{1}{t_{f}}$, we obtain the leading order term to be $\frac{1}{4t_{f}^{2}}$. As we see in Fig. 2, this asymptotic prediction is close to the asymptotic scaling of the numerical solution.

Finally, consider the time-averaged trace-norm distance between two time-evolving states $\ket{\chi_{1}(t)}$ and $\ket{\chi_{2}(t)}$:

[TABLE]

The results comparing the numerical solution to the naive adiabatic approximation and the two lowest WKB approximation orders are shown in Fig. 3. As expected, the naive adiabatic approximation becomes better as $t_{f}$ increases, and the same is true for the WKB approximations, which are both more accurate than the adiabatic approximation. Moreover, the first-order WKB approximation is better than the adiabatic approximation according to the time-averaged trace-norm distance metric.

Taken together, the results for the $n=1$ case show that both the zeroth-order WKB approximation and the first-order WKB approximation consistently improve upon the naive adiabatic approximation, and the first-order WKB approximation can be used to pick out more subtle features of the quantum evolution.

V.2 The $n$-qubit Grover problem

We next turn to a study of the Grover problem as a function of problem size $n$, with $n>1$. The quantity of interest to us is how long we need to run the adiabatic algorithm before a certain threshold probability of success $p_{\mathrm{Th}}$ is exceeded. The associated threshold timescale is defined as:

[TABLE]

Here, $p_{\mathrm{GS}}(t)$ represents the probability of finding the ground state at the end an adiabatic evolution of time $t$. We choose $p_{\mathrm{Th}}=0.95$ (we have checked that the results are insensitive to changing $p_{\mathrm{Th}}$).

First, in Fig. 4 we show how $t_{f}^{\mathrm{Th}}(n)$ scales for the numerical solution, under the four different schedules defined in Eqs. (19). It appears as though the scaling for the $g_{3}(r)$ schedule is better than the theoretically optimal scaling of $2^{n/2}$ Bennett et al. (1997), but this is a small $n$ effect as shown in Fig. 5.

Next, we examine how well the WKB approximation does in predicting these scalings. In Fig. 6 we plot the scaling of $t_{f}^{\mathrm{Th}}(n)$ for the same four schedules, under the lowest order WKB approximation. Only the $g_{2}(r)$ schedule (which slows as the inverse-square of the gap) yields the correct scaling of $t_{f}^{\mathrm{Th}}(n)$. This is also the schedule which yields the smallest instantaneous and time-averaged trace-norm distance, as shown in Figs. 7(a) and 7(b), respectively. For the other schedules, Fig. 6 shows that the WKB approximation gives answers that are dramatically different from the exact solution. Furthermore, for the $g_{0}$ and $g_{1}$ schedules, the scaling with $n$ of $t_{f}^{\mathrm{Th}}(n)$ violates the query complexity bound Bennett et al. (1997).

Why do the approximations for the $g_{0}$, $g_{1}$, and $g_{3}$ schedules give us the wrong scalings, while the approximation for the $g_{2}$ schedule gives us the correct scaling? A possible answer lies in the steepness of the final-time approximate success probability curves for the different schedules. In Fig. 8 we show the $p_{\mathrm{GS}}(t_{f})$ curves for all four schedules for $n=4$ ($K=15$) as predicted by the first-order WKB approximation (the highest order at which we are able to obtain analytic expressions). For the $g_{0}$ and $g_{1}$ schedules, we see that the final ground state probability rises very sharply and exceeds unity (very slightly for $g_{0}$), and thereby becomes nonsensical [see Fig. 8], while for the $g_{2}$ and $g_{3}$ schedules, $p_{\mathrm{GS}}(t_{f})\leq 1$ [see Fig. 8]. Further, we observe that the curves are ordered from steepest to shallowest rise as $g_{0}$, $g_{1}$, $g_{2}$, $g_{3}$. We conjecture that this rise in $p_{\mathrm{GS}}$ with $t_{f}$ continues to slow down with increasing $\alpha$. Thus the $g_{2}$ schedule captures the right scaling [in Fig. 6] by capturing the right steepness: for $\alpha<2$ the rise is too steep, and for $\alpha>2$ the rise is too shallow. A full explanation of this phenomenon is left to future work, but we speculate that the $g_{0}$ and $g_{1}$ schedules correspond to effective Hamiltonians that no longer represent the Grover problem.

Given that the WKB approximation gives consistent results only for the $g_{2}(r)$ schedule, we focus on this schedule and examine where the WKB approximation performs better than the naive adiabatic approximation. As can be seen in Fig. 9, for the $g_{2}(r)$ schedule, the WKB approximation always has an advantage over the naive adiabatic approximation. Further, it is clear that the advantage is bigger for smaller evolution times $t_{f}$ and for larger problem sizes $n$.

As we have indicated above, the WKB approximants are generically not normalized: they can be sub-normalized or super-normalized. Two questions arise: (1) Is the degree of non-normalization a good indicator of approximation quality of the WKB aproximation? (2) Does renormalization by fiat improve the quality of the approximation?

First, in Fig. 10 we plot the norm of the WKB approximation at the lowest order for the case of $n=6$ and $t_{f}=60$ as a function of the anneal parameter $r$ for all four schedules. In this regime, the WKB approximation is sub-normalized for all schedules. Arranging the schedules from farthest from normalization to closest to normalization, we have: $g_{3},g_{2},g_{1},g_{0}$, with $g_{0}$ and $g_{1}$ closest to being normalized. On the other hand, we have seen [Fig. 7(a) and Figs. 6,7(b)] that the best approximation was obtained for the $g_{2}$ schedule. Thus, we conclude that the degree of non-normalization is not a good indicator of the quality of approximation.

Second, we consider renormalization of the WKB approximation, by which we mean:

[TABLE]

where we have denoted the renormalized WKB approximants as $\ket{\chi_{\mathrm{rWKB}}}$. Note that it is somewhat ad hoc to normalize the approximant: a renormalization step is not part of the standard WKB approximation procedure. With this caveat, we now analyze the behavior of the WKB approximants after renormalization.

First we consider the time-averaged trace-norm distance, $\overline{\mathcal{D}}$ [Eq. (33)]. As shown in Fig. 11, renormalizing the WKB approximation does improve the approximation.

However, the situation changes when we consider the threshold timescale $t_{f}^{\mathrm{Th}}(n)$ [Eq. (34)], shown in Fig. 12. We see that renormalizing the WKB approximation yields highly unphysical results for the $g_{0},g_{1}$, and $g_{2}$ schedules. In particular, for the $g_{0}$ and $g_{1}$ schedules, we see that the renormalized WKB approximation predicts a scaling for $t_{f}^{\mathrm{Th}}(n)$ that decreases with problem size $n$. For the $g_{2}$ schedule we see a scaling of $\mathcal{O}(1)$. So, while the unnormalized WKB predicted the correct scaling for the $g_{2}$ schedule, the renormalized WKB does not retain that feature. On the other hand, for the $g_{3}$ schedule, we see that the renormalized WKB predicts the correct scaling of $\mathcal{O}(2^{n/2})$, fixing the incorrect scaling of the unnormalized WKB approximation for that schedule.

VI Summary and Conclusions

We have presented a straightforward technique to obtain an analytic asymptotic approximation to slowly evolving $2$-level systems by adapting the WKB method. We have applied it to a problem that is motivated by adiabatic quantum computation: the Hamiltonian Grover search problem. This problem has a physical Hilbert space of dimension $2^{n}$, but is effectively constrained to a $2$-dimensional subspace. We have seen that in this case when $n=1$, the WKB method provides good approximations, especially to the population dynamics. We saw that the WKB approximation can capture fluctuations in the population that are absent in the purely adiabatic (ground state) evolution. Thus, the WKB is quasi-adiabatic. For completeness, in Appendix A, we compare our WKB approximation to the asymptotic expansion method of Hagedorn and Joye Hagedorn and Joye (2002), and show that the latter misses the oscillations that are captured by the quasi-adiabatic WKB expansion.

Turning to the Grover problem with $n>1$ and with different interpolation schedules, we observed that the WKB approximation yields meaningful results only for the schedule which slows quadratically with the ground state gap. For this $g_{2}(r)$ schedule, the WKB approximation is able to capture the scaling with $n$ of $t_{f}^{\mathrm{Th}}$, and hence recovers the quantum speedup, even at the lowest approximation order. On the other hand, for the schedules that slow down more slowly than quadratically with the gap, the WKB approximation violates normalization and predicts an impossible faster-than-quadratic quantum speedup for the Grover problem.

We also saw that, using the time-averaged trace-norm distance, for the $g_{2}(r)$ schedule, the WKB approximation always does better than a naive adiabatic approximation, and the advantage becomes more pronounced for larger system sizes and shorter evolution times.

Turning to the question of the whether the norm of the WKB approximation is a good signal of approximation quality, we saw that this is not the case. Further, we saw that enforcing renormalization by fiat gives mixed results. On the one hand, it lead to an improvement in the time-averaged trace-distance and also gave the right predictions for scaling of the threshold timescale for the $g_{3}$ schedule. On the other hand, it gave incorrect predictions for the scaling of the threshold timescales for the $g_{0},g_{1},$ and $g_{2}$ schedules, especially degrading the prediction for the $g_{2}$ schedule compared to its unnormalized counterpart.

An interesting problem for future work is to provide more rigorous justifications and explanations for when and where the WKB approximation provides good approximations. With a better understanding of the regions where the WKB approximation performs well, and if the approximation errors are better controlled, the method could be used in the design of quantum control protocols to implement quantum gates A.P. Peirce, M.A. Dahleh, H. Rabitz (1988); D’alessandro and Dahleh (2001); Grace et al. (2007); Barnes and Das Sarma (2012); Palittapongarnpim et al. (2017).

Acknowledgements.

We are grateful to an anonymous referee who provided several constructive suggestions. We also thank Tameem Albash and Milad Marvian for helpful discussions and comments. This work was supported under ARO grant number W911NF-12-1-0523 and NSF grant number INSPIRE-1551064.

Appendix A Comparison with the method of Hagedorn and Joye

In this section, we recap the asymptotic expansion of Hagedorn and Joye Hagedorn and Joye (2002), which is a powerful tool for proving adiabatic theorems. In particular, the Hagedorn and Joye method can be used to prove bounds on the error incurred due to their asymptotic expansion; in fact, the main goal of Ref. Hagedorn and Joye (2002) was to show that the adiabatic approximation can provide exponentially small errors if the Hamiltonian is analytic in the time-variable (see also Refs. Nenciu (1993); Lidar et al. (2009); Ge et al. (2016)). Here, we analyze its utility as a computational tool.

Hagedorn and Joye (HJ) propose the following method to obtain asymptotic approximations to the time-dependent Schrödinger equation

[TABLE]

Note that the above equation is of the form of Eq. (6), with $\epsilon\equiv\frac{1}{\mu t_{f}}$ and $\mathcal{H}(r)\equiv s^{\prime}(r)H(r)$.

They obtain a theorem which states that for any value of the small parameter $\epsilon$, one can write down an approximation for $\ket{\chi(r)}$ which takes the form of a power series in $\epsilon$. The quality of the approximation (as measured by the 2-norm) scales as $e^{-\frac{1}{\epsilon}}$ provided that the number of terms in the series scales as $1/\epsilon$. More precisely:

Theorem 1** (Hagedorn and Joye (2002)).**

Assume reasonable smoothness and gap conditions on the Hamiltonian. We can then recursively obtain an asymptotic expansion of the form

[TABLE]

such that for any $r$, there exist positive $G$, $C(g)$, and $\Gamma(g)$ such that for all $g\in(0,G)$, the vector $\ket{\chi^{(\lfloor{g/\epsilon}\rfloor)}_{\mathrm{HJ}}(r,\epsilon)}$ satisfies

[TABLE]

for all $\epsilon\leq 1$. Here, $\ket{\chi(r,\epsilon)}$ is the Schrödinger evolved wavefunction starting from the initial condition $\ket{\chi(0,\epsilon)}=\ket{\chi^{(\lfloor{g/\epsilon}\rfloor)}_{\mathrm{HJ}}(0,\epsilon)}$.

We explore the usefulness of this asymptotic expansion as an approximation tool and thus we do not estimate the number of terms that are necessary to provide an exponentially small error. Instead, we develop the approximation for two orders and compare the resulting asymptotic expansion with the WKB method.

Let us develop the terms in the HJ expansion (as given in Ref. Hagedorn and Joye (2002)). We substitute the asymptotic ansatz

[TABLE]

into the Schrödinger equation, and equate the terms multiplying the same order of $\epsilon$, which results in the following expression for the $O(\epsilon^{j})$ term

[TABLE]

Here $\ket{\Phi(r)}$ is the eigenstate being (approximately) followed, and the other components of the above formula are obtained recursively by using:

[TABLE]

where, in going from Eq. (41b) to Eq. (41c), we integrated by parts and used $\braket{\Phi(q)}{\psi_{j-1}^{\perp}(q)}=0$. Also, $P_{\perp}(r)\equiv I-\ket{\Phi(r)}\bra{\Phi(r)}$ is the instantaneous projector on to the complement of $\ket{\Phi(r)}$; $E(r)$ is the eigenvalue being quasi-adiabatically followed; and $[H(r)-E(r)]_{R}^{-1}$ is the reduced resolvent, i.e., the inverse of $[H(r)-E(r)]$ restricted to the complement of $\ket{\Phi(r)}$.

In order to compare the HJ expansion with the WKB approximation, we will compare the $N$-th order expansion provided by both methods. Note that the $N$-th order of the HJ expansion includes terms up to $\mathcal{O}(\epsilon^{N+1})$. This means that we will be comparing the zeroth order of WKB (i.e., $\ket{\chi_{\mathrm{WKB}}^{(0)}}$) with

[TABLE]

and the first order WKB (i.e., $\ket{\chi_{\mathrm{WKB}}^{(1)}}$) with

[TABLE]

For two-level systems such as the one that we are concerned with, we obtain the following simplified expressions, where “$\mathrm{GS}$” and “$\mathrm{Exc}$” denote the ground and excited states respectively and $\Delta$ represents the spectral gap:

[TABLE]

We have assumed that the ground state is being followed and hence set $\ket{\Phi}=\ket{\chi_{\mathrm{GS}}}$. We have also used the fact that for real-valued Hamiltonians in two dimensions $\braket{\chi_{\mathrm{GS}}^{\prime}}{\chi_{\mathrm{GS}}}=0$. (Note that this does not mean $\ket{\chi_{\mathrm{GS}}^{\prime}}=\ket{\chi_{\mathrm{Exc}}}$ because $\ket{\chi_{\mathrm{GS}}^{\prime}}$ is generally not normalized and carries a non-trivial phase.)

We now restrict to the case of a qubit in a magnetic field.

First, consider $\ket{\chi_{\mathrm{HJ}}^{(0)}}$ (which includes terms up to order $\epsilon$). Figure 13(a) shows that $\ket{\chi_{\mathrm{HJ}}^{(0)}}$ provides an approximation that is ‘too adiabatic’. In particular, it fails to capture the oscillations that are captured by the WKB approximation, as seen in Fig. 1. Furthermore, from the form of $\ket{\chi_{\mathrm{HJ}}^{(0)}}$ it is clear that this approximation will predict $p_{\mathrm{GS}}(t_{f})=1$ always:

[TABLE]

Next, consider $\ket{\chi_{\mathrm{HJ}}^{(1)}}$ (which includes terms up to order $\epsilon^{2}$). Figure 13(b) shows that this too provides an approximation which fails to capture the oscillations that are present in the numerical solution and also in the lowest order WKB solution. Thus, we conclude that the WKB method is more suitable for developing analytic approximations.

While we pointed out some of the disadvantages of the HJ method as an approximation technique, we remark that the method is particularly useful to prove scaling results. For example, consider,

[TABLE]

In the first line, we used the fact $\ket{\chi_{\mathrm{GS}}}$ is orthogonal to any (unnormalized) state that carries the $\perp$ symbol. In the second line, we used the fact that

[TABLE]

is purely imaginary. Thus the HJ expansion captures the $1-\mathcal{O}(\frac{1}{t_{f}^{2}})$ scaling of the final ground state probability.