Conserved Quantities in Linear and Nonlinear Quantum Search

Conserved quantities play important roles in physics: energy, momentum, and angular momentum are conserved in systems that are time, translation, and rotation invariant, respectively. In quantum mechanics, the wave function ψ of a system with Hamiltonian H₀ evolves according to Schrödinger’s equation, 1iℏ∂ψ∂t=H0ψ,i\hbar {{\partial \psi } \over {\partial t}} = {H_0}\psi , where we will work in units where the reduced Planck constant ħ = 1 throughout this paper. The Hamiltonian H₀ must be Hermitian (i.e., self-adjoint), which causes the time evolution to be unitary. Then, |ψ|² = 〈ψ|ψ〉 is a conserved quantity, meaning a normalized wave function stays normalized, which is necessary for the probabilistic interpretation of the wave function à la the Born rule. This may not be the only conserved quantity, however. For example, if the Hamiltonian is time-independent, then the expected value of the Hamiltonian 〈H₀ = 〈ψ|H₀|ψ〉, i.e., the average energy of an ensemble of identically prepared systems, is also conserved [1]. The time-independence of the Hamiltonian and the conservation of its expected value hold for many quantum algorithms based on continuoustime quantum walks [2–4], where the particle walks on the vertices of a graph, and the kinetic energy of the system is proportional to the discrete Laplacian. For overviews of quantum walks, see [5–8]. In Section 2, we will give an example of this by reviewing how a continuous-time quantum walk searches the complete graph for a marked vertex [3,9], which is the combinatorial formulation of the unstructured search problem of Grover’s algorithm [10]. In the process, we will give a new derivation of the evolution of the algorithm for arbitrary jumping rates, as well as give a new derivation—using only elementary calculus—of the critical jumping rate that causes the algorithm to succeed with certainty. We will explicitly show that for the Hamiltonian used in the search algorithm, 〈H₀〉 is conserved.

Some many-body quantum systems are approximately described by an effectively nonlinear Schrödinger-type equation with a cubic nonlinearity proportional to |ψ|²ψ, i.e., 2iℏ∂ψ∂t=(H0+λ|ψ|2)ψ,i\hbar {{\partial \psi } \over {\partial t}} = \left( {{H_0} + \lambda |\psi {|^2}} \right)\psi , where the effective Hamiltonian 3H(t)=H0+λ|ψ|2H(t) = {H_0} + \lambda |\psi {|^2} depends on time because it depends on the state ψ, which evolves with time. A Bose-Einstein condensate [11–13] is an example of such a system, with the Gross-Pitaevskii equation [14,15] describing it in the mean-field limit and taking the form of the cubic nonlinear Schrödinger equation (2); for reviews, see [16,17]. As before, |ψ|² = 1 is conserved, but since H(t) is time-dependent, 〈H(t)〉 = 〈ψ|H(t)|ψ〉 is not conserved. When H₀ is a constant, however, then 〈H0+12λ|ψ|2〉=〈ψ|H0+12λ|ψ|2|ψ〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle = \langle \psi |{H_0} + {1 \over 2}\lambda |\psi {|^2}|\psi \rangle is conserved [18]. In Section 3, we will give an example of this by reviewing how a nonlinear quantum walk with λ > 0, which corresponds to a Bose-Einstein condensate with repulsive interactions, searches the complete graph for a marked vertex [19], and we will explicitly show that 〈H0+12λ|ψ|2〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle is conserved for the algorithm.

If H₀ depends on time, however, then 〈H0+12λ|ψ|2〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle is generally not conserved. In Section 4, we give an example of this by reviewing how a nonlinear quantum walk with λ < 0, which corresponds to a Bose-Einstein condensate with attractive interactions, searches a complete graph of N vertices for a marked vertex [20]. In this algorithm, the jumping rate of the quantum walk is a critical, time-dependent function γ_c(t), which causes H₀ to be time-dependent as well. While 〈H0+12λ|ψ|2〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle is not conserved in this case, we show that the critical jumping rate is such that 〈H₀/[γ_c(t)N]〉 is conserved instead.

Thus, in this tutorial, we will review three continuous-time quantum walk algorithms for searching the complete graph: a linear algorithm in Section 2, a nonlinear algorithm corresponding to a Bose-Einstein condensate with repulsive interactions in Section 3, and a nonlinear algorithm corresponding to a Bose-Einstein condensate with attractive interactions in Section 4. For each, we will demonstrate a conserved quantity related to the (effective) Hamiltonian of the system, thus making connections between quantum algorithms, conservation laws, and manybody quantum systems. We conclude in Section 5.

In quantum mechanics, the Hamiltonian H₀ is the operator corresponding to the total energy of the system, which consists of a kinetic energy term that is proportional to Laplace’s operator ▽² and a potential energy term V [1]: H0=−ℏ22m∇2+V.{H_0} = - {{{\hbar ^2}} \over {2m}}{\nabla ^2} + V.

For a quantum walk on a graph of N vertices, the vertices label computational basis states |1〉, |2〉, …, |N〉, and in general, the state is a superposition or complex linear combination over the vertices. Then, Laplace’s operator in the continuum should be replaced by the discrete Laplacian L: H0=−γL+V,{H_0} = - \gamma L + V, where L is an N × N matrix with L_ij = 1 if vertices i and j are adjacent, L_ii = – deg(i) is the negative of the number of neighbors of vertex i, and L_ij = 0 otherwise. For the complete graph of N vertices, an example of which is shown in Figure 1, the discrete Laplacian has 1 for each off-diagonal element and –(N – 1) on the diagonal. We have also written the coefficient of L as a parameter γ, which corresponds to the jumping rate of the quantum walk. To construct a quantum search algorithm [3], an oracle “marks” a particular vertex |a〉 that we want to find, which manifests as a potential energy term: 4H0=−γL−|a〉〈a|.{H_0} = - \gamma L - |a\rangle \langle a|.

Figure 1.

A complete graph of N = 6 vertices. A vertex is marked, as indicated by a double circle. In the search algorithm, vertices that evolve identically have the same label and color.

Initially, we have no knowledge of where the marked vertex might be, so we guess each vertex equally by beginning in a uniform superposition over all N vertices: 5|ψ(0)〉=1N∑i=1N|i〉.|\psi (0)\rangle = {1 \over {\sqrt N }}\sum\limits_{i = 1}^N | i\rangle .

Then, the system evolves by Schrödinger’s equation (1), and since H₀ is time-independent, 6|ψ(t)〉=e−iH0t|ψ(0)〉,|\psi (t)\rangle = {e^{ - i{H_0}t}}|\psi (0)\rangle , where we have taken ħ = 1. The success probability at time t is the probability that measuring the position of the particle at time t results in it being found at the marked vertex, so the success probability is 7p(t)=|〈a∣ψ(t)〉|2.p(t) = \;|\langle a\mid \psi (t)\rangle {|^2}.

The evolution of this success probability is plotted in Figure 2 for search on a complete graph of N = 100 vertices and with various values of the jumping rate γ. In Figure 2a, the jumping rate starts off small as γ = 0.001 in the solid black curve, and the success probability is relatively unchanged from its initial value of 1/100 = 0.01. As the jumping rate increases, however, the success probability evolves with a higher and higher peak, eventually reaching a peak success probability of 1 when γ = 0.01. In Figure 2b, γ is increased further, and the peak success probability decreases, with larger γ resulting in the success probability staying near its initial value.

Figure 2.

Success probability versus time for (linear) quantum search on the complete graph with N = 100 vertices and various jumping rates. In (a), the solid black curve is γ = 0.001, the dashed red curve is γ = 0.005, the dotted green curve is γ = 0.008, the dot-dashed blue curve is γ = 0.009, and the dot-dot-dashed orange curve is γ = 0.01. In (b), the solid black curve is γ = 0.011, the dashed red curve is γ = 0.012, the dotted green curve is γ = 0.015, the dot-dashed blue curve is γ = 0.02, and the dot-dot-dashed orange curve is γ = 0.03.

To prove this behavior, we begin by noting that due to the symmetry of the problem, there are only two types of vertices: the marked vertex and the unmarked vertices, labeled in Figure 1 as a and b. Vertices of the same type evolve identically, so the system evolves in a two-dimensional (2D) subspace spanned by the marked vertex |a〉 and the uniform superposition of unmarked vertices, which we will call |b〉: |b〉=1N−1∑i unmarked |i〉.|b\rangle = {1 \over {\sqrt {N - 1} }}\mathop \sum \limits_{iunmarked} |i\rangle .

In this basis, the initial state (5) is 8|ψ(0)〉=1N|a〉+N−1N|b〉=1N(1N−1),|\psi (0)\rangle = {1 \over {\sqrt N }}|a\rangle + \sqrt {{{N - 1} \over N}} |b\rangle = {1 \over {\sqrt N }}\left( \matrix{ 1 \cr \sqrt {N - 1} \cr} \right), and the Hamiltonian (4) is 9=(γ(N−1)−1−γN−1−γN−1γ).{H_0} = \left( {\matrix{ {\gamma (N - 1) - 1} & { - \gamma \sqrt {N - 1} } \cr { - \gamma \sqrt {N - 1} } & \gamma \cr } } \right).

The (unnormalized) eigenvectors and eigenvalues of H₀ are ψ1=(γN−2γ−1−ΔE)|a〉−2γN−1|b〉,E1=γN−1−ΔE2,ψ2=(γN−2γ−1+ΔE)|a〉−2γN−1|b〉,E2=γN−1+ΔE2,\matrix{ {{\psi _1} = (\gamma N - 2\gamma - 1 - \Delta E)|a\rangle - 2\gamma \sqrt {N - 1} |b\rangle ,} \hfill & {{E_1} = {{\gamma N - 1 - \Delta E} \over 2},} \hfill \cr {{\psi _2} = (\gamma N - 2\gamma - 1 + \Delta E)|a\rangle - 2\gamma \sqrt {N - 1} |b\rangle ,} \hfill & {{E_2} = {{\gamma N - 1 + \Delta E} \over 2},} \hfill \cr } where 10ΔE=E2−E1=γ2N2−2γN+4γ+1\Delta E = {E_2} - {E_1} = \sqrt {{\gamma ^2}{N^2} - 2\gamma N + 4\gamma + 1} is the energy gap. Expressing the initial state (8) as a linear combination of these eigenvectors, |ψ(0)〉=−γN+1−ΔE4γNΔEψ1+γN−1−ΔE4γNΔEψ2.|\psi (0)\rangle = {{ - \gamma N + 1 - \Delta E} \over {4\gamma \sqrt N \Delta E}}{\psi _1} + {{\gamma N - 1 - \Delta E} \over {4\gamma \sqrt N \Delta E}}{\psi _2}.

Then, using (6), the state of the system at time t is 11|ψ(t)〉=e−iH0t(−γN+1−ΔE4γNΔEψ1+γN−1−ΔE4γNΔEψ2)=−γN+1−ΔE4γNΔEe−iE1tψ1+γN−1−ΔE4γNΔEe−iE2tψ2=e−i(γN−1)t/24γNΔE[ (−γN+1−ΔE)eiΔEt/2ψ1+(γN−1−ΔE)e−iΔEt/2ψ2 ]=e−i(γN−1)t/24γNΔE{ (−γN+1−ΔE)eiΔEt/2[ (γN−2γ−1−ΔE)|a〉−2γN−1|b〉 ]+(γN−1−ΔE)e−iΔEt/2[ (γN−2γ−1+ΔE)|a〉−2γN−1|b〉 ] } =e−i(γN−1)t/24γNΔE{ [ (eiΔEt/2−e−iΔEt/2)((−γN+1)(γN−2γ−1)+(ΔE)2) +(eiΔEt/2+e−iΔEt/2)(−ΔE(γN−2γ−1)+(γN−1)ΔE) ]|a〉−2γN−1 [ (eiΔEt/2−e−iΔEt/2)(−γN+1) +(eiΔEt/2+e−iΔEt/2)(−ΔE) ]|b〉 }=e−i(γN−1)t/24γNΔE{ [ 2isin(ΔEt2)[ (−γN+1)(γN−2γ−1)+(ΔE)2 ] +2cos(ΔEt2)2γΔE ]|a〉−2γN−1[ 2isin(ΔEt2)(−γN+1)+2cos(ΔEt2)(−ΔE) ]|b〉 }.\matrix{ {|\psi (t)\rangle } \hfill & { = {e^{ - i{H_0}t}}\left( {{{ - \gamma N + 1 - \Delta E} \over {4\gamma \sqrt N \Delta E}}{\psi _1} + {{\gamma N - 1 - \Delta E} \over {4\gamma \sqrt N \Delta E}}{\psi _2}} \right)} \hfill \cr {} \hfill & { = {{ - \gamma N + 1 - \Delta E} \over {4\gamma \sqrt N \Delta E}}{e^{ - i{E_1}t}}{\psi _1} + {{\gamma N - 1 - \Delta E} \over {4\gamma \sqrt N \Delta E}}{e^{ - i{E_2}t}}{\psi _2}} \hfill \cr {} \hfill & { = {{{e^{ - i(\gamma N - 1)t/2}}} \over {4\gamma \sqrt N \Delta E}}\left[ {( - \gamma N + 1 - \Delta E){e^{i\Delta Et/2}}{\psi _1} + (\gamma N - 1 - \Delta E){e^{ - i\Delta Et/2}}{\psi _2}} \right]} \hfill \cr {} \hfill & { = {{{e^{ - i(\gamma N - 1)t/2}}} \over {4\gamma \sqrt N \Delta E}}\left\{ {( - \gamma N + 1 - \Delta E){e^{i\Delta Et/2}}\left[ {(\gamma N - 2\gamma - 1 - \Delta E)|a\rangle - 2\gamma \sqrt {N - 1} |b\rangle } \right]\left. { + \;(\gamma N - 1 - \Delta E){e^{ - i\Delta Et/2}}\left[ {(\gamma N - 2\gamma - 1 + \Delta E)|a\rangle - 2\gamma \sqrt {N - 1} |b\rangle } \right]} \right\}} \right.} \hfill \cr {} \hfill & { = {{{e^{ - i(\gamma N - 1)t/2}}} \over {4\gamma \sqrt N \Delta E}}\left\{ {\left[ {\left( {{e^{i\Delta Et/2}} - {e^{ - i\Delta Et/2}}} \right)\left( {( - \gamma N + 1)(\gamma N - 2\gamma - 1) + {{(\Delta E)}^2}} \right)} \right.} \right.\left. { + \left( {{e^{i\Delta Et/2}} + {e^{ - i\Delta Et/2}}} \right)( - \Delta E(\gamma N - 2\gamma - 1) + (\gamma N - 1)\Delta E)} \right]|a\rangle - 2\gamma \sqrt {N - 1} \left[ {\left( {{e^{i\Delta Et/2}} - {e^{ - i\Delta Et/2}}} \right)( - \gamma N + 1)} \right.\left. {\left. { + \left( {{e^{i\Delta Et/2}} + {e^{ - i\Delta Et/2}}} \right)( - \Delta E)} \right]|b\rangle } \right\}} \hfill \cr {} \hfill & { = {{{e^{ - i(\gamma N - 1)t/2}}} \over {4\gamma \sqrt N \Delta E}}\left\{ {\left[ {2i\sin \left( {{{\Delta Et} \over 2}} \right)\left[ {( - \gamma N + 1)(\gamma N - 2\gamma - 1) + {{(\Delta E)}^2}} \right]} \right.\left. { + 2\cos \left( {{{\Delta Et} \over 2}} \right)2\gamma \Delta E} \right]|a\rangle - \left. {2\gamma \sqrt {N - 1} \left[ {2i\sin \left( {{{\Delta Et} \over 2}} \right)( - \gamma N + 1) + 2\cos \left( {{{\Delta Et} \over 2}} \right)( - \Delta E)} \right]|b\rangle } \right\}.} \right.} \hfill \cr }

Taking the norm-square of the amplitude of |a〉, the success probability (7) at time t is 12p(t)=116γ2N(ΔE)2{ 4sin2(ΔEt2)[ (−γN+1)(γN−2γ−1)+(ΔE)2 ]2 +16γ2(ΔE)2cos2(ΔEt2) }=((−γN+1)(γN−2γ−1)+(ΔE)22γNΔE)2sin2(ΔEt2)+1Ncos2(ΔEt2).\matrix{ {p(t)} \hfill & { = {1 \over {16{\gamma ^2}N{{(\Delta E)}^2}}}\left\{ {4{{\sin }^2}\left( {{{\Delta Et} \over 2}} \right){{\left[ {( - \gamma N + 1)(\gamma N - 2\gamma - 1) + {{(\Delta E)}^2}} \right]}^2}} \right.\left. { + 16{\gamma ^2}{{(\Delta E)}^2}{{\cos }^2}\left( {{{\Delta Et} \over 2}} \right)} \right\}} \hfill \cr {} \hfill & { = {{\left( {{{( - \gamma N + 1)(\gamma N - 2\gamma - 1) + {{(\Delta E)}^2}} \over {2\gamma \sqrt N \Delta E}}} \right)}^2}{{\sin }^2}\left( {{{\Delta Et} \over 2}} \right) + {1 \over N}{{\cos }^2}\left( {{{\Delta Et} \over 2}} \right).} \hfill \cr }

This agrees with the curves in Figure 2 with N = 100 and various values of γ. The success probability reaches a maximum value when the sine is 1 and the cosine is 0, i.e., at time 13t*=πΔE,{t_*} = {\pi \over {\Delta E}}, at which the success probability peaks at 14p*=((−γN+1)(γN−2γ−1)+(ΔE)22γNΔE)2.{p_*} = {\left( {{{( - \gamma N + 1)(\gamma N - 2\gamma - 1) + {{(\Delta E)}^2}} \over {2\gamma \sqrt N \Delta E}}} \right)^2}.

To the best of our knowledge, the above results for |ψ(t)〉 (11), p(t) (12), t_* (13), and p_* (14) for arbitrary γ are new, as previous literature [9] focused on the specific case when γ = 1/N, motivated by minimizing the energy gap [3] or by degenerate perturbation theory [21]. Now, we will provide another way to arrive at the γ = 1/N case, namely, using elementary calculus to prove that it is the value of γ that maximizes the peak success probability p_* (14). Substituting the energy gap ΔE (10) into p_* (14) and differentiating the result with respect to γ, we get dp*dγ=4(N−1)(γ2N2−1)N(γ2N2−2γ(N−1)+1)2.{{d{p_*}} \over {d\gamma }} = {{4(N - 1)\left( {{\gamma ^2}{N^2} - 1} \right)} \over {N{{\left( {{\gamma ^2}{N^2} - 2\gamma (N - 1) + 1} \right)}^2}}}.

Setting this equal to 0 and solving for γ, the peak success probability p_* is maximized when γ takes a critical value of γc=1N,{\gamma _c} = {1 \over N}, at which the energy gap (10) is ΔEc=2N,\Delta {E_c} = {2 \over {\sqrt N }}, and the peak success probability (14) is p*c=1.{p_{*c}} = 1.

That is, γ = γ_c = 1/N makes the algorithm deterministic, i.e., the marked vertex is found with certainty. Using (13), this occurs at time t*c=π2N,{t_{*c}} = {\pi \over 2}\sqrt N , which is the expected O(N)O(\sqrt N ) scaling of Grover’s algorithm [10]. Furthermore, using (12), the success probability evolves with time as pc(t)=sin2(tN)+1Ncos2(tN),{p_c}(t) = {\sin ^2}\left( {{t \over {\sqrt N }}} \right) + {1 \over N}{\cos ^2}\left( {{t \over {\sqrt N }}} \right), which agrees with the dot-dot-dashed orange curve in Figure 2a, where N = 100 and γ = 1/100 = 0.01.

Moving on to conserved quantities, since the Hamiltonian (4) is time-independent, its expected value is conserved, i.e., 〈H〉 = 〈ψ(t)|H|ψ(t)〉 is constant in time [1]. We will explicitly show this in the 2D subspace of our algorithm for pedagogical purposes, as the calculation for the nonlinear algorithm in the next section will use similar methods. To begin, we write the state of the system in the {|a〉, |b〉} basis as 15|ψ(t)〉=α(t)|a〉+β(t)|b〉=(α(t)β(t))|\psi (t)\rangle = \alpha (t)|a\rangle + \beta (t)|b\rangle = \left( \matrix{ \alpha (t) \cr \beta (t) \cr} \right) for complex α and β, where for convenience, we have stopped writing the explicit time-dependence of α and β. The Hamiltonian (9) takes the form 16H0=(abb*d),{H_0} = \left( {\matrix{ a & b \cr {{b^*}} & d \cr } } \right), where 17a=γ(N−1)−1,b=−γN−1, and d=γa = \gamma (N - 1) - 1,\quad b = - \gamma \sqrt {N - 1} ,\quad {\rm{ and }}\quad d = \gamma are real and time-independent, although we allow for b to be complex for greater generality. Then, the expected value of H₀ is 18〈 H0 〉=〈ψ|H0|ψ〉=(α*β*)(abb*d)(αβ)=(α*β*)(aα+bβb*α+dβ)=aα*α+bα*β+b*β*α+dβ*β.\matrix{ {\left\langle {{H_0}} \right\rangle } \hfill & { = \langle \psi |{H_0}|\psi \rangle } \hfill \cr {} \hfill & { = \left( {\matrix{ {{\alpha ^*}} \hfill & {{\beta ^*}} \hfill \cr } } \right)\left( {\matrix{ a & b \cr {{b^*}} & d \cr } } \right)\left( \matrix{ \alpha \cr \beta \cr} \right)} \hfill \cr {} \hfill & { = \left( {\matrix{ {{\alpha ^*}} \hfill & {{\beta ^*}} \hfill \cr } } \right)\left( \matrix{ a\alpha + b\beta \cr {b^*}\alpha + d\beta \cr} \right)} \hfill \cr {} \hfill & { = a{\alpha ^*}\alpha + b{\alpha ^*}\beta + {b^*}{\beta ^*}\alpha + d{\beta ^*}\beta .} \hfill \cr }

We want to show that this is constant in time. Using the chain rule, its time-derivative is 19d〈 H0 〉dt=∂〈 H0 〉∂αdαdt+∂〈 H0 〉∂βdβdt+∂〈 H0 〉∂α*dα*dt+∂〈 H0 〉∂β*dβ*dt=∂〈 H0 〉∂αα˙+∂〈 H0 〉∂ββ˙+∂〈 H0 〉∂α*α˙*+∂〈 H0 〉∂β*β˙*,\matrix{ {{{d\left\langle {{H_0}} \right\rangle } \over {dt}}} \hfill & { = {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \alpha }}{{d\alpha } \over {dt}} + {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \beta }}{{d\beta } \over {dt}} + {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\alpha ^*}}}{{d{\alpha ^*}} \over {dt}} + {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\beta ^*}}}{{d{\beta ^*}} \over {dt}}} \hfill \cr {} \hfill & { = {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \alpha }}\dot \alpha + {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \beta }}\dot \beta + {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\alpha ^*}}}{{\dot \alpha }^*} + {{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\beta ^*}}}{{\dot \beta }^*},} \hfill \cr } where we used a dot to indicate a time derivative. Let us find the partial derivatives of 〈H₀〉 and then the time derivatives of α, β, and their conjugates. First, we differentiate (18) to get 20∂〈 H0 〉∂α=aα*+b*β*,∂〈 H0 〉∂β=bα*+dβ*,∂〈 H0 〉∂α*=aα+bβ,∂〈 H0 〉∂β*=b*α+dβ.\matrix{ {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \alpha }}} \hfill & { = a{\alpha ^*} + {b^*}{\beta ^*},} \hfill \cr {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \beta }}} \hfill & { = b{\alpha ^*} + d{\beta ^*},} \hfill \cr {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\alpha ^*}}}} \hfill & { = a\alpha + b\beta ,} \hfill \cr {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\beta ^*}}}} \hfill & { = {b^*}\alpha + d\beta .} \hfill \cr }

Then, to find α˙,β˙,α˙*\dot \alpha ,\dot \beta ,{{\dot \alpha }^*}, and β˙*{{\dot \beta }^*} in (19), we use Schrödinger’s equation (1): i(α˙β˙)=(abb*d)(αβ),i\left( \matrix{ {\dot \alpha } \cr {\dot \beta } \cr} \right) = \left( {\matrix{ a & b \cr {{b^*}} & d \cr } } \right)\left( \matrix{ \alpha \cr \beta \cr} \right), which multiplied out is iα˙=aα+bβ,iβ˙=b*α+dβ.\matrix{ {i\dot \alpha = a\alpha + b\beta ,} \hfill \cr {i\dot \beta = {b^*}\alpha + d\beta .} \hfill \cr }

Solving for α˙{\dot \alpha } and β˙{\dot \beta } and taking the complex conjugate of each equation, 21α˙=−i(aα+bβ),β˙=−i(b*α+dβ),α˙*=i(aα*+b*β*),β˙*=i(bα*+dβ*).\matrix{ {\dot \alpha = - i(a\alpha + b\beta ),} \hfill \cr {\dot \beta = - i\left( {{b^*}\alpha + d\beta } \right),} \hfill \cr {{{\dot \alpha }^*} = i\left( {a{\alpha ^*} + {b^*}{\beta ^*}} \right),} \hfill \cr {{{\dot \beta }^*} = i\left( {b{\alpha ^*} + d{\beta ^*}} \right).} \hfill \cr }

Plugging (20) and (21) into (19), we get 22d〈 H0 〉dt=(aα*+b*β*)[−i(aα+bβ)]+(bα*+dβ*)[ −i(b*α+dβ) ]+(aα+bβ)[ i(aα*+b*β*) ]+(b*α+dβ)[ i(bα*+dβ*) ]=0,\matrix{ {{{d\left\langle {{H_0}} \right\rangle } \over {dt}}} \hfill & { = \;\left( {a{\alpha ^*} + {b^*}{\beta ^*}} \right)[ - i(a\alpha + b\beta )] + \left( {b{\alpha ^*} + d{\beta ^*}} \right)\left[ { - i\left( {{b^*}\alpha + d\beta } \right)} \right]} \hfill \cr {} \hfill & { + (a\alpha + b\beta )\left[ {i\left( {a{\alpha ^*} + {b^*}{\beta ^*}} \right)} \right] + \left( {{b^*}\alpha + d\beta } \right)\left[ {i\left( {b{\alpha ^*} + d{\beta ^*}} \right)} \right]} \hfill \cr {} \hfill & { = 0,} \hfill \cr } so we have proved that the expected value of the Hamiltonian is constant.

Our derivation also shows that there is a relationship between the derivatives of the expected value of the Hamiltonian and the time derivatives of the amplitudes. Comparing (20) and (21), 23∂〈 H0 〉∂α=−iα˙*,∂〈 H0 〉∂β=−iβ˙*,∂〈 H0 〉∂α*=iα˙,∂〈 H0 〉∂β*=iβ˙.\matrix{ {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \alpha }} = - i{{\dot \alpha }^*},} \hfill \cr {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial \beta }} = - i{{\dot \beta }^*},} \hfill \cr {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\alpha ^*}}} = i\dot \alpha ,} \hfill \cr {{{\partial \left\langle {{H_0}} \right\rangle } \over {\partial {\beta ^*}}} = i\dot \beta .} \hfill \cr }

Plugging these into (19), we get a cleaner proof that the expected value of the Hamiltonian is conserved: d〈 H0 〉dt=−iα˙*α˙−iβ˙*β˙+iα˙α˙*+iβ˙β˙*=−i|α|2−i|β|2+i|α|2+i|β|2=0.\matrix{ {{{d\left\langle {{H_0}} \right\rangle } \over {dt}}} \hfill & { = - i{{\dot \alpha }^*}\dot \alpha - i{{\dot \beta }^*}\dot \beta + i\dot \alpha {{\dot \alpha }^*} + i\dot \beta {{\dot \beta }^*}} \hfill \cr {} \hfill & { = - i|\alpha {|^2} - i|\beta {|^2} + i|\alpha {|^2} + i|\beta {|^2}} \hfill \cr {} \hfill & { = 0.} \hfill \cr }

Of course, the cleanest proof that the expected value of this time-independent Hamiltonian is conserved follows from the fact that the time evolution multiplies any eigenvector ψ_j of H₀ by the phase e^iE_jt. This leaves the norm-squared of the projection of the state onto each eigenspace invariant, and also, a fortiori, 〈H₀〉. This proof fails, however, for the time-dependent evolutions we will consider in the next two sections. This is why we gave the preceding argument, which can be translated into the nonlinear settings.

Many-body quantum systems can evolve by effective nonlinearities. Bose-Einstein condensates are a quintessential example, and for pairwise contact interactions in the mean field limit, they evolve by the Gross-Pitaevskii equation, which is a Schrödinger-type equation with a cubic nonlinearity (2). A nonlinear quantum search algorithm [20] can be constructed from this by taking H₀ = –γL – |a〉〈a| (4) so that the effective Hamiltonian (3) becomes 24H(t)=−γL−|a〉〈a|+λ∑i=1N|〈i∣ψ(t)〉|2|i〉〈i|.H(t) = - \gamma L - |a\rangle \langle a| + \lambda \sum\limits_{i = 1}^N | \langle i\mid \psi (t)\rangle {|^2}|i\rangle \langle i|.

When the nonlinearity coefficient λ is positive, this physically corresponds to bosons with repulsive interactions, and when λ is negative, the interactions are attractive. The sign and strength of λ can be tuned, for example, by Feshbach resonance [22]. In this section, we consider λ > 0, and in the next section, we consider λ < 0.

As shown in [19], for large N, when γ takes a critical value of γc=2−λ2N{\gamma _c} = {{2 - \lambda } \over {2N}} and the nonlinearity coefficient λ lies in the range 0 < λ < λ_c, where λc=42+N′{\lambda _c} = {4 \over {2 + {{\sqrt N }^\prime }}} then the system evolves from the initial uniform superposition (5) to the marked vertex |a〉 in time roughly t*=π3N.{t_*} = {\pi \over {\sqrt 3 }}\sqrt N .

This runtime is slower than the linear algorithm’s πN/2\pi \sqrt N /2 (13) by a factor of 2/3≈1.1552/\sqrt 3 \approx 1.155, i.e., it is slower by 15.5%. As discussed in [19], although this algorithm searches in the same O(N)O(\sqrt N ) time as the linear quantum algorithm, but with a worse coefficient, it nonetheless shows that it is theoretically possible for a Bose-Einstein condensate with repulsive interactions evolving by an effective nonlinearity to implement a quantum search algorithm with an O(N)O(\sqrt N ) scaling.

Note from [19] that λ_c is not a tight upper bound, i.e., there may exist some nonlinearity coefficient λ > λ_c for which the success probability reaches 1, but eventually as λ increases, the success probability stops reaching 1. This is shown in Figure 3 for search on a complete graph of N = 100 vertices, for which λc=4/(2+100)=4/12=1/3≈0.333{\lambda _c} = 4/(2 + \sqrt {100} ) = 4/12 = 1/3 \approx 0.333. The solid black curve is λ = 0.2, since 0 < 0.2 < 0.333, the success probability must evolve to 1, as it does. The dashed red curve is λ = 0.4, and while it is greater than λ_c = 0.333, the success probability still reaches 1. The dotted green curve is λ = 0.611, and while the success probability struggles to grow, it does eventually reach 1 off the right of the graph. In contrast, the dot-dashed blue curve is λ = 0.612, and the success probability does not reach 1 anymore. A larger nonlinearity coefficient reaches a smaller and smaller success probability, such as the dot-dot-dashed orange curve with λ = 0.8. A tight upper bound on λ such that the success probability reaches 1 remains an open question.

Figure 3.

Success probability versus time for nonlinear quantum search on the complete graph with N = 100 vertices, λ = (2 – λ)/2N, and various values of λ. The solid black curve is γ = 0.2, the dashed red curve is λ = 0.6, the dotted green curve is λ = 0.611, the dot-dashed blue curve is λ = 0.612, and the dot-dot-dashed orange curve is λ = 0.8.

Since the effective Hamiltonian of the search algorithm depends on time, 〈H(t)〉 is generally not conserved. Instead, the expected value of H₀ plus half the nonlinear self-potential is conserved, i.e., 〈 H0+12λ|ψ|2 〉\left\langle {{H_0} + {1 \over 2}\lambda |\psi {|^2}} \right\rangle is conserved, as this is the quantity that corresponds to the total energy [18]. Conserved quantities for the Gross-Pitaevskii equation have been explored in some other contexts, including as a classical Hamiltonian [19], the focusing cubic and quintic Gross-Pitaevskii hierarchies [23], cubic Gross-Pitaevskii hierarchy on ℝ [24], and the one-dimensional Gross-Pitaevskii equation [25,26], but here we concentrate on conserved quantities related to the expected value of the effective Hamiltonian.

For instructional purposes, let us use the nonlinear search algorithm to explicitly show why 〈H(t)〉 = 〈H₀ + λ|ψ|²〉 is not conserved, while 〈 H0+12λ|ψ|2 〉\left\langle {{H_0} + {1 \over 2}\lambda |\psi {|^2}} \right\rangle is conserved. Even with the nonlinearity, as before, the system evolves in the same 2D subspace spanned by |a and |b〉, so the state |ψ(t)〉 can be written as (15). Now, the effective Hamiltonian (24) is 25H(t)=(γ(N−1)−1+λ|α|2−γN−1−γN−1γ+λN−1|β|2).H(t) = \left( {\matrix{ {\gamma (N - 1) - 1 + \lambda |\alpha {|^2}} & { - \gamma \sqrt {N - 1} } \cr { - \gamma \sqrt {N - 1} } & {\gamma + {\lambda \over {N - 1}}|\beta {|^2}} \cr } } \right).

This takes the form H(t)=(a+f|α|2bb*d+g|β|2),H(t) = \left( {\matrix{ {a + f|\alpha {|^2}} & b \cr {{b^*}} & {d + g|\beta {|^2}} \cr } } \right), where a, b, and d are defined in (17), and 26f=λ and g=λN−1f = \lambda \quad {\rm{ and }}\quad g = {\lambda \over {N - 1}} are real and time independent. Now, the expected value of the effective Hamiltonian is 27〈H(t)〉=〈ψ|H(t)|ψ〉=(α*β*)(a+f|α|2bb*d+g|β|2)(αβ)=(α*β*)(aα+bβ+f|α|2αb*α+dβ+g|β|2β)=aα*α+bα*β+f|α|4+b*β*α+dβ*β+g|β|4.\matrix{ {\langle H(t)\rangle } \hfill & { = \langle \psi |H(t)|\psi \rangle } \hfill \cr {} \hfill & { = \left( {\matrix{ {{\alpha ^*}} \hfill & {{\beta ^*}} \hfill \cr } } \right)\left( {\matrix{ {a + f|\alpha {|^2}} & b \cr {{b^*}} & {d + g|\beta {|^2}} \cr } } \right)\left( \matrix{ \alpha \cr \beta \cr} \right)} \hfill \cr {} \hfill & { = \left( {\matrix{ {{\alpha ^*}} \hfill & {{\beta ^*}} \hfill \cr } } \right)\left( \matrix{ a\alpha + b\beta + f|\alpha {|^2}\alpha \cr {b^*}\alpha + d\beta + g|\beta {|^2}\beta \cr} \right)} \hfill \cr {} \hfill & { = a{\alpha ^*}\alpha + b{\alpha ^*}\beta + f|\alpha {|^4} + {b^*}{\beta ^*}\alpha + d{\beta ^*}\beta + g|\beta {|^4}.} \hfill \cr }

Differentiating this with respect to time using the chain rule, 28d〈H(t)〉dt=∂〈H(t)〉∂αdαdt+∂〈H(t)〉∂βdβdt+∂〈H(t)〉∂α*dα*dt+∂〈H(t)〉∂β*dβ*dt=∂〈H(t)〉∂αα˙+∂〈H(t)〉∂ββ˙+∂〈H(t)〉∂α*α˙*+∂〈H(t)〉∂β*β˙*.\matrix{ {{{d\langle H(t)\rangle } \over {dt}}} \hfill & { = {{\partial \langle H(t)\rangle } \over {\partial \alpha }}{{d\alpha } \over {dt}} + {{\partial \langle H(t)\rangle } \over {\partial \beta }}{{d\beta } \over {dt}} + {{\partial \langle H(t)\rangle } \over {\partial {\alpha ^*}}}{{d{\alpha ^*}} \over {dt}} + {{\partial \langle H(t)\rangle } \over {\partial {\beta ^*}}}{{d{\beta ^*}} \over {dt}}} \hfill \cr {} \hfill & { = {{\partial \langle H(t)\rangle } \over {\partial \alpha }}\dot \alpha + {{\partial \langle H(t)\rangle } \over {\partial \beta }}\dot \beta + {{\partial \langle H(t)\rangle } \over {\partial {\alpha ^*}}}{{\dot \alpha }^*} + {{\partial \langle H(t)\rangle } \over {\partial {\beta ^*}}}{{\dot \beta }^*}.} \hfill \cr }

The partial derivatives of the expected value of the effective Hamiltonian can be found by differentiating (27), noting that |α|⁴ = (α*)²α² and |β|⁴ = (β*)²β²: 29∂〈H(t)〉∂α=aα*+b*β*+2f|α|2α*,∂〈H(t)〉∂β=bα*+dβ*+2g|β|2β*,∂〈H(t)〉∂α*=aα+bβ+2f|α|2α,∂〈H(t)〉∂β*=b*α+dβ+2g|β|2β.\matrix{ {{{\partial \langle H(t)\rangle } \over {\partial \alpha }} = a{\alpha ^*} + {b^*}{\beta ^*} + 2f|\alpha {|^2}{\alpha ^*},} \hfill \cr {{{\partial \langle H(t)\rangle } \over {\partial \beta }} = b{\alpha ^*} + d{\beta ^*} + 2g|\beta {|^2}{\beta ^*},} \hfill \cr {{{\partial \langle H(t)\rangle } \over {\partial {\alpha ^*}}} = a\alpha + b\beta + 2f|\alpha {|^2}\alpha ,} \hfill \cr {{{\partial \langle H(t)\rangle } \over {\partial {\beta ^*}}} = {b^*}\alpha + d\beta + 2g|\beta {|^2}\beta .} \hfill \cr }

To find α˙,β˙,α˙*\dot \alpha ,\dot \beta ,{{\dot \alpha }^*}, and β˙*{{\dot \beta }^*} in (28), we use the cubic nonlinear Schrödinger equation (2), i(α˙β˙)=(a+f|α|2bb*d+g|β|2)(αβ),i\left( \matrix{ {\dot \alpha } \cr {\dot \beta } \cr} \right) = \left( {\matrix{ {a + f|\alpha {|^2}} & b \cr {{b^*}} & {d + g|\beta {|^2}} \cr } } \right)\left( \matrix{ \alpha \cr \beta \cr} \right), which multiplied out is iα˙=aα+bβ+f|α|2α,iβ˙=b*α+dβ+g|β|2β.\matrix{ {i\dot \alpha = a\alpha + b\beta + f|\alpha {|^2}\alpha ,} \hfill \cr {i\dot \beta = {b^*}\alpha + d\beta + g|\beta {|^2}\beta .} \hfill \cr }

Solving for α˙{\dot \alpha } and β˙{\dot \beta } and taking the complex conjugate of each equation, 30α˙=−i(aα+bβ+f|α|2α),β˙=−i(b*α+dβ+g|β|2β),α˙*=i(aα*+b*β*+f|α|2α*),β˙*=i(bα*+dβ*+g|β|2β*).\matrix{ {\dot \alpha = - i\left( {a\alpha + b\beta + f|\alpha {|^2}\alpha } \right),} \hfill \cr {\dot \beta = - i\left( {{b^*}\alpha + d\beta + g|\beta {|^2}\beta } \right),} \hfill \cr {{{\dot \alpha }^*} = i\left( {a{\alpha ^*} + {b^*}{\beta ^*} + f|\alpha {|^2}{\alpha ^*}} \right),} \hfill \cr {{{\dot \beta }^*} = i\left( {b{\alpha ^*} + d{\beta ^*} + g|\beta {|^2}{\beta ^*}} \right).} \hfill \cr }

Comparing (29) and (30), we see that these are no longer proportional to partial derivatives of the expected value of the Hamiltonian, i.e., (23) no longer holds. Then, (28) is not equal to zero, so the expected value of the effective Hamiltonian is not conserved.

The reason why (29) is not proportional to (30) is the extra factor of 2 in the last term. We can eliminate this factor of 2 by dividing the nonlinear self-potential by 2. That is, we consider H0+12λ|ψ|2=(γ(N−1)−1+12λ|α|2−γN−1−γN−1γ+λ2(N−1)|β|2),{H_0} + {1 \over 2}\lambda |\psi {|^2} = \left( {\matrix{ {\gamma (N - 1) - 1 + {1 \over 2}\lambda |\alpha {|^2}} & { - \gamma \sqrt {N - 1} } \cr { - \gamma \sqrt {N - 1} } & {\gamma + {\lambda \over {2(N - 1)}}|\beta {|^2}} \cr } } \right), which takes the form 31H0+12λ|ψ|2=(a+12f|α|2bb*d+12g|β|2),\matrix{ {{H_0} + {1 \over 2}\lambda |\psi {|^2} = \left( {\matrix{ {a + {1 \over 2}f|\alpha {|^2}} & b \cr {{b^*}} & {d + {1 \over 2}g|\beta {|^2}} \cr } } \right)} \cr } , where a, b, and d are defined in (17), and f and g are defined in (26). The expected value of this is 32〈〈H0+12λ|ψ|2=〈 ψ|H0+12λ|ψ|2|ψ 〉=(α*β*)(a+12f|α|2bb*d+12g|β|2)(αβ)=(α*β*)(aα+bβ+12f|α|2αb*α+dβ+12g|β|2β)=aα*α+bα*β+12f|α|4+bβ*α+dβ*β+12g|β|4,\matrix{ {\left\langle {{{\left. {\langle {H_0} + {1 \over 2}\lambda |\psi } \right|}^2}} \right.} \hfill & { = \left\langle {\psi |{H_0} + {1 \over 2}\lambda |\psi {|^2}|\psi } \right\rangle } \hfill \cr {} \hfill & { = \left( {\matrix{ {{\alpha ^*}} \hfill & {{\beta ^*}} \hfill \cr } } \right)\left( {\matrix{ {a + {1 \over 2}f|\alpha {|^2}} & b \cr {{b^*}} & {d + {1 \over 2}g|\beta {|^2}} \cr } } \right)\left( \matrix{ \alpha \cr \beta \cr} \right)} \hfill \cr {} \hfill & { = \left( {\matrix{ {{\alpha ^*}} \hfill & {{\beta ^*}} \hfill \cr } } \right)\left( \matrix{ a\alpha + b\beta + {1 \over 2}f|\alpha {|^2}\alpha \cr {b^*}\alpha + d\beta + {1 \over 2}g|\beta {|^2}\beta \cr} \right)} \hfill \cr {} \hfill & { = a{\alpha ^*}\alpha + b{\alpha ^*}\beta + {1 \over 2}f|\alpha {|^4} + b{\beta ^*}\alpha + d{\beta ^*}\beta + {1 \over 2}g|\beta {|^4},} \hfill \cr } which has a time derivative of 33d〈H012λ|ψ|2〉dt=∂〈H0+12λ|ψ|2〉∂αdαdt+∂〈H0+12λ|ψ|2〉∂βdβdt+∂〈H0+12λ|ψ|2〉∂α*dα*dt+∂〈H0+12λ|ψ|2〉∂β*dβ*dt=∂〈H0+12λ|ψ|2〉∂αα˙+∂〈H0+12λ|ψ|2〉∂ββ˙+∂〈H0+12λ|ψ|2〉∂α*α˙*+∂〈H0+12λ|ψ|2〉∂β*β˙*.\matrix{ {{{d\langle {H_0}{1 \over 2}\lambda |\psi {|^2}\rangle } \over {dt}}} \hfill & { = {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \alpha }}{{d\alpha } \over {dt}} + {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \beta }}{{d\beta } \over {dt}}} \hfill \cr {} \hfill & { + {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\alpha ^*}}}{{d{\alpha ^*}} \over {dt}} + {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\beta ^*}}}{{d{\beta ^*}} \over {dt}}} \hfill \cr {} \hfill & { = {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \alpha }}\dot \alpha + {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \beta }}\dot \beta } \hfill \cr {} \hfill & { + {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\alpha ^*}}}{{\dot \alpha }^*} + {{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\beta ^*}}}{{\dot \beta }^*}.} \hfill \cr }

Now, when we take the partial derivatives of (32), the factors of 2 are eliminated: ∂〈H0+12λ|ψ|2〉∂α=aα*+bβ*+f|α|2α*,∂〈H0+12λ|ψ|2〉∂β=bα*+dβ*+g|β|2β*,∂〈H0+12λ|ψ|2〉∂α*=aα+bβ+f|α|2α,∂〈H0+12λ|ψ|2〉∂β*=bα+dβ+g|β|2β.\matrix{ {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \alpha }}} \hfill & { = a{\alpha ^*} + b{\beta ^*} + f|\alpha {|^2}{\alpha ^*},} \hfill \cr {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \beta }}} \hfill & { = b{\alpha ^*} + d{\beta ^*} + g|\beta {|^2}{\beta ^*},} \hfill \cr {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\alpha ^*}}}} \hfill & { = a\alpha + b\beta + f|\alpha {|^2}\alpha ,} \hfill \cr {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\beta ^*}}}} \hfill & { = b\alpha + d\beta + g|\beta {|^2}\beta .} \hfill \cr }

Comparing these with (30), these are exactly proportional to the time derivatives of α*, β*, α, and β, i.e., ∂〈H0+12λ|ψ|2〉∂α=−iα˙*,∂〈H0+12λ|ψ|2〉∂β=−iβ˙*,∂〈H0+12λ|ψ|2〉∂α*=iα˙,∂〈H0+12λ|ψ|2〉∂β*=iβ˙.\matrix{ {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \alpha }}} \hfill & { = - i{{\dot \alpha }^*},} \hfill \cr {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial \beta }}} \hfill & { = - i{{\dot \beta }^*},} \hfill \cr {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\alpha ^*}}}} \hfill & { = i\dot \alpha ,} \hfill \cr {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {\partial {\beta ^*}}}} \hfill & { = i\dot \beta .} \hfill \cr }

Then, substituting into (33), we get ∂〈H0+12λ|ψ|2〉dt=−iα˙*α˙−iβ˙*β˙+iα˙α˙*+iβ˙β˙*=0.\matrix{ {{{\partial \langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle } \over {dt}}} \hfill & { = - i{{\dot \alpha }^*}\dot \alpha - i{{\dot \beta }^*}\dot \beta + i\dot \alpha {{\dot \alpha }^*} + i\dot \beta {{\dot \beta }^*}} \hfill \cr {} \hfill & { = 0.} \hfill \cr }

Thus, 〈H0+12λ|ψ|2〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle is conserved in the nonlinear quantum search algorithm with repulsive interactions.

In this section, we consider the nonlinear quantum search algorithm with negative values of λ, which corresponds to a Bose-Einstein condensate with attractive interactions. Intuitively, as the amplitude at the marked vertex builds up, the attractive interaction should draw additional amplitude, speeding up the algorithm. Indeed, as shown in [20], when the jumping rate is the time-varying critical function 34γc(t)=1N[ 1−λ(|α|2−|β|2N−1) ]{\gamma _c}(t) = {1 \over N}\left[ {1 - \lambda \left( {|\alpha {|^2} - {{|\beta {|^2}} \over {N - 1}}} \right)} \right] [which is determined independently of the marked vertex |a〉, and whose time dependence comes from the amplitudes α and β (15)], then the algorithm evolves from the uniform superposition over all N vertices to the marked vertex in time t*=11−λπ2N.{t_*} = {1 \over {\sqrt {1 - \lambda } }}{\pi \over 2}\sqrt N .

This is demonstrated in Figure 4, where we plot the evolution of the success probability for search on the complete graph of N = 100 vertices. The solid black, dashed red, dotted green, and dot-dashed blue curves correspond to λ = 0, –1, –2, and –3, and they reach success probabilities of 1 at time t*=(1/1−λ)(π/2)N≈15.708{t_*} = (1/\sqrt {1 - \lambda } )(\pi /2)\sqrt N \approx 15.708, 11.107, 9.069, and 7.854, respectively, and so the more negative the nonlinearity coefficient, the faster the algorithm.

Figure 4.

Success probability versus time for nonlinear quantum search on the complete graph with N = 100 vertices, γ = γ_c(t), and various values of λ. The solid black curve is λ = 0, the dashed red curve is λ = –1, the dotted green curve is λ = –2, and the dot-dashed blue curve is λ = –3.

This speedup comes, however, at the expense of increasing the necessary time-measurement precision. In Figure 4, as the magnitude of the nonlinearity coefficient increases, the peak in success probability becomes more narrow, which necessitates greater time-measurement precision to catch the peak. This was quantified in [20], which showed that the width of the peak scales as O(N/(1−λ))O(\sqrt N /(1 - \lambda )), and so the time-measurement precision is constant when λ=Θ(N)\lambda = \Theta (\sqrt N ), at which the runtime t_* = O(N^1/4), which is a quadratic improvement over the linear and nonlinear repulsive algorithms.

Since the effective Hamiltonian H(t) = H₀ + λ|ψ|² depends on time, its expected value is not conserved. H0+12λ|ψ|2{{\rm{H}}_0}{\rm{ + }}{1 \over 2}\lambda |\psi {|^2} takes the form given in (31), but now a, b, and d depend on time because γ varies with time. Then, the time-derivative of 〈H0+12λ|ψ|2〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle from (33) will include additional terms involving time-derivatives of a, b, and d, and 〈H0+12λ|ψ|2〉\langle {H_0} + {1 \over 2}\lambda |\psi {|^2}\rangle will not be conserved. Instead, let us show that another quantity is conserved.

First, we review from [20] that the critical function γ_c (t) for the jumping rate causes H(t) to be proportional to H₀. To begin, H(t)=−γc(t)L−|a〉〈a|+λ∑i=1N|〈i∣ψ(t)〉|2|i〉〈i|=−γc(t)L−|a〉〈a|+λ|〈a∣ψ(t)〉|2|a〉〈a|+λ∑i≠a|〈i∣ψ(t)〉|2|i〉〈i|=−γc(t)L−|a〉〈a|+λ|α|2|a〉〈a|+λ|β|2N−1∑i≠a|i〉〈i|=−γc(t)L−(1−λ|α|2)|a〉〈a|+λ|β|2N−1∑i≠a|i〉〈i|.\matrix{ {H(t)} \hfill & { = - {\gamma _c}(t)L - |a\rangle \langle a| + \lambda \sum\limits_{i = 1}^N | \langle i\mid \psi (t)\rangle {|^2}|i\rangle \langle i|} \hfill \cr {} \hfill & { = - {\gamma _c}(t)L - |a\rangle \langle a| + \lambda |\langle a\mid \psi (t)\rangle {|^2}|a\rangle \langle a| + \lambda \sum\limits_{i \ne a} | \langle i\mid \psi (t)\rangle {|^2}|i\rangle \langle i|} \hfill \cr {} \hfill & { = - {\gamma _c}(t)L - |a\rangle \langle a| + \lambda |\alpha {|^2}|a\rangle \langle a| + \lambda {{|\beta {|^2}} \over {N - 1}}\sum\limits_{i \ne a} | i\rangle \langle i|} \hfill \cr {} \hfill & { = - {\gamma _c}(t)L - \left( {1 - \lambda |\alpha {|^2}} \right)|a\rangle \langle a| + \lambda {{|\beta {|^2}} \over {N - 1}}\sum\limits_{i \ne a} | i\rangle \langle i|.} \hfill \cr }

Now, recall that the jumping rate satisfies (34). Multiplying both sides of (34) by N and distributing λ yields γc(t)N=1−λ|α|2−λ|β|2N−1.{\gamma _c}(t)N = 1 - \lambda |\alpha {|^2} - \lambda {{|\beta {|^2}} \over {N - 1}}.

Rearranging, 1−λ|α|2=γc(t)N+λ|β|2N−1.1 - \lambda |\alpha {|^2} = {\gamma _c}(t)N + \lambda {{|\beta {|^2}} \over {N - 1}}.

Substituting this into our previous expression for H(t), H(t)=−γc(t)L−(γc(t)N+λ|β|2N−1)|a〉〈a|+λ|β|2N−1∑i≠a|i〉〈i|=−γc(t)L−γc(t)N|a〉〈a|+λ|β|2N−1∑i=1N|i〉〈i|=−γc(t)L−γc(t)N|a〉〈a|+λ|β|2N−1𝕀\matrix{ {H(t)} \hfill & { = - {\gamma _c}(t)L - \left( {{\gamma _c}(t)N + \lambda {{|\beta {|^2}} \over {N - 1}}} \right)|a\rangle \langle a| + \lambda {{|\beta {|^2}} \over {N - 1}}\sum\limits_{i \ne a} | i\rangle \langle i|} \hfill \cr {} \hfill & { = - {\gamma _c}(t)L - {\gamma _c}(t)N|a\rangle \langle a| + \lambda {{|\beta {|^2}} \over {N - 1}}\sum\limits_{i = 1}^N | i\rangle \langle i|} \hfill \cr {} \hfill & { = - {\gamma _c}(t)L - {\gamma _c}(t)N|a\rangle \langle a| + \lambda {{|\beta {|^2}} \over {N - 1}}{\rm{I}}} \hfill \cr } where 𝕀 is the identity matrix, which can be dropped since it corresponds to a rescaling of energy or a global, unobservable phase. Then, H(t)=−γc(t)L−γc(t)N|a〉〈a|=γc(t)N(−1NL−|a〉〈a|)=γc(t)NH0|γc.\matrix{ {H(t)} \hfill & { = - {\gamma _c}(t)L - {\gamma _c}(t)N|a\rangle \langle a|} \hfill \cr {} \hfill & { = {\gamma _c}(t)N\left( { - {1 \over N}L - |a\rangle \langle a|} \right)} \hfill \cr {} \hfill & { = {{\left. {{\gamma _c}(t)N{H_0}} \right|}_{{\gamma _c}}}.} \hfill \cr }

That is, the effective Hamiltonian for the nonlinear search algorithm with attractive interactions is a time-varying rescaling of the linear search algorithm’s Hamiltonian evaluated at its critical jumping rate of γ_c = 1/N. Solving for the linear algorithm’s Hamiltonian at its critical jumping rate, H0|γc=H(t)γc(t)N.{\left. {{H_0}} \right|_{{\gamma _c}}} = {{H(t)} \over {{\gamma _c}(t)N}}.

Since 〈H₀〉 is conserved for arbitrary constant γ, as we proved in (22), we have that 〈 H(t)γc(t)N 〉=0,\left\langle {{{H(t)} \over {{\gamma _c}(t)N}}} \right\rangle = 0, and we have proved a conserved quantity for the nonlinear search algorithm with attractive interactions.

We have reviewed three quantum search algorithms governed by continuous-time quantum walks. The first used a linear quantum walk, and we derived the evolution of the system for arbitrary jumping rates. We used elementary calculus to prove that there is a critical jumping rate γ_c that causes the success probability to reach 1 in time πN/2\pi \sqrt N /2. We showed that the expected value of the Hamiltonian H₀ is conserved, which is an example of the fact that all time-independent Hamiltonians have constant expected values. Second, we considered a nonlinear quantum walk with repulsive interactions, whose effective Hamiltonian was H(t) = H₀ + λ|ψ|² with λ > 0. For a range of possible critical jumping rates, the success probability reaches 1 in time πN/3\pi \sqrt {N/3} , which is slower than the linear algorithm by a constant factor, and we showed that the expected value of H0+12λ|ψ|2{H_0} + {1 \over 2}\lambda |\psi {|^2} is conserved. Third, for the nonlinear quantum walk with attractive interactions, the effective Hamiltonian has λ < 0, and when the jumping rate takes a critical function γ_c (t), the success probability reaches 1 in time πN/(1−λ)/2\pi \sqrt {N/(1 - \lambda )} /2, which is faster than the other algorithms for all λ < 0. The jumping rate varies such that the effective Hamiltonian H(t) is equal to a time-dependent rescaling of the linear algorithm’s Hamiltonian H₀ evaluated at its critical jumping rate γ_c, namely, H(t) = γ_c(t)NH₀|_γc, and since the expected value of the time-independent linear Hamiltonian H₀ is conserved for an arbitrary constant jumping rate, the expected value of H(t)/[γ_c (t)N] is conserved for the third algorithm. These results provide a perspective on quantum walk search algorithms through the lens of conservation laws and many-body quantum theory.