Nonlinear quantum search using the Gross–Pitaevskii equation

We solve the unstructured search problem in constant time by computing with a physically motivated nonlinearity of the Gross–Pitaevskii type. This speedup comes, however, at the novel expense of increasing the time-measurement precision. Jointly optimizing these resource requirements results in an overall scaling of N1/4. This is a significant, but not unreasonable, improvement over the N1/2 scaling of Grover's algorithm. Since the Gross–Pitaevskii equation approximates the multi-particle (linear) Schrödinger equation, for which Grover's algorithm is optimal, our result leads to a quantum information-theoretic lower bound on the number of particles needed for this approximation to hold, asymptotically.


Introduction
Abrams and Lloyd [1] argued that a nonlinear quantum theory could result in unreasonable computational advantages by giving two examples of nonlinear algorithms that solve NP-complete and #P problems in polynomial time.Both of their algorithms can be implemented by a nonlinear Schrödinger-type evolution in which the time derivatives of the state components depend upon their hyperbolic tangents [2,3].The derivative of tanh x at x = 0 is 1, so this is a strongly nonlinear system in which 0 is an unstable fixed point.The strength of the nonlinearity provides a large computational advantage, but it also makes the system highly susceptible to noise [1,2,3].
An obvious question is whether a modest, physically motivated nonlinearity can still produce a computational advantage.In particular, consider Bose-Einstein condensates (BECs).In 1924-25, Bose and Einstein predicted that cooling a dilute gas of bosons near absolute zero would cause the atoms to occupy their lowest quantum state, forming a new state of matter where quantum effects would be macroscopically apparent [4,5,6].It took seventy years, however, for these BECs to be experimentally produced [7,8,9].In general, describing such many-body systems is difficult because of the many interaction terms.But under certain conditions, one can assume that only two-body contact interactions contribute and the s-wave scattering length a is much less than the interparticle spacing.Then using mean field theory, one finds that the system is approximately described by a nonlinear Schrödinger equation with a cubic nonlinearity: where H 0 includes the kinetic energy and trapping potential, m is the mass of the condensate atom, and N 0 is the number of condensate atoms.‡ The validity of this celebrated Gross-Pitaevskii equation [10,11] requires that N 0 be much greater than 1-but how much greater?
We provide a quantum information-theoretic solution to this question in the context of solving the unstructured quantum search problem [12] in continuous-time [13] using the Gross-Pitaevskii equation as the governing equation.The cubic nonlinearity in (1) has zero derivative at zero, making it softer than those considered by Abrams and Lloyd [1].We quantify the computational advantage that such a nonlinearity provides for the unstructured search problem compared to standard quantum computation; this requires considering time-measurement precision as a physical resource.Since this advantage cannot persist when the Gross-Pitaevskii equation is recognized as an approximation to an underlying multi-particle Schrödinger equation, for which Grover's algorithm is optimal, we arrive at a quantum information-theoretic lower bound on the number of particles, N 0 , needed for this approximation to hold, asymptotically.

Setup
The system evolves in a N -dimensional Hilbert space with computational basis {|0 , . . ., |N − 1 }.The initial state |ψ(0) is an equal superposition |s of all these basis states: The goal is to find a particular "marked" basis state, which we label |w .Let's first review the linear solution.We use Childs and Goldstone's [14] notation for Farhi and Gutmann's [13] Hamiltonian: where γ is a parameter, inversely proportional to mass.The system evolves in the twodimensional subspace spanned by |w and |s .One might (correctly) reason that the success of the algorithm depends on the value of γ.This can be seen in figure 1, which shows the difference in eigenvalues of −H 0 and the overlaps of its nontrivial eigenvectors with |s and |w .When γ takes a critical value of γ c = 1/N , the eigenstates of −H 0 are proportional to ±|w + |s , and the corresponding energy gap is 2/ √ N .So the Schrödinger evolution rotates the state from |s to |w in time π √ N /2.‡ More generally, the cubic nonlinear Schrödinger equation is the equation of motion for φ 4 field theory.In the nonlinear regime, we include an additional nonlinear "self-potential" V (t) so that the system evolves according to the Gross-Pitaevskii equation (1): ]ψ(r, t), where g > 0. This corresponds to a BEC with attractive interactions, and thus a negative scattering length [15,16].Heuristically, as probability accumulates at the marked state due to the |w w| term in H 0 , the self-potential attracts more probability, speeding up the search.Thus we expect larger g to result in a faster algorithm.
In the computational basis, the self-potential is Even with this nonlinearity, the system remains in the subspace spanned by {|w , |s } throughout its evolution.We define a vector which is orthonormal to |w .Then the state of the system |ψ(t) can be written as Writing the Gross-Pitaevskii equation in this {|w , |r } basis, we get d dt where we've defined

Critical Gamma
Before proceeding with further analytical calculations, we build some intuition by examining two plots.For constant γ and g, the success probability as a function of time, |α(t)| 2 , is plotted in figure 2 along with the linear result.The nonlinear algorithm underperforms the linear one in this case.This is true in general for constant γ and g, and it can be understood by examining the time-dependence of the critical value of γ, which is the value of γ that ensures that the eigenstates of A are in the form ±|w + |s .Initially, γ c = 1/N .Then, as shown in figure 3, it shifts to a larger value.If γ is constant, it will not follow this shift, we will no longer have the desired eigenstates, and the algorithm will perform poorly.
To determine how γ c varies with time, we find the eigenvectors of A and choose γ so that they have the desired form ±|w + |s .To eliminate fractions in the subsequent algebra, we rescale the nonlinearity coefficient g by defining Solving the characteristic equation gives the eigenvalues of A: where the gap between them is and we've defined The corresponding eigenvectors of A are The critical value of γ ensures that these eigenvectors have the form ±|w + |s .That is, Solving this yields: Note that in the linear limit (G = 0), this reduces to γ c = 1/N , as expected.Importantly, since δ varies with time, (3) implies γ c also varies with time, in agreement with our previous discussion about figures 2 and 3.

Runtime
For the remainder of the paper, we choose time-varying γ = γ c according to (3).Before analytically determining the consequences of this, let's again consider a plot.Figure 4 shows the success probability as a function of time.There are several observations.First, the success probability reaches 1, which occurs because we constructed the eigenstates to make this happen.Second, as N increases, the runtime remains constant.Third, the success probability is periodic.Finally, the peak in success probability becomes increasingly narrow for large N .Let's now analytically prove the second, third, and fourth observations.To begin, we explicitly write out (2) to get two coupled, first-order ordinary differential equations for α(t) and β(t): We can decouple these equations for the success probability The details of this decoupling procedure are given in Appendix A. To solve this uncoupled equation, we use separation of variables and integrate from t = 0 to t and x = 1/N to x, which yields Solving for x, the success probability as a function of time is From this, the success probability reaches 1 when the tangent term is zero, which first occurs at time This runtime is exactly constant for G = 1.Also, when G = Θ(1), the runtime for large N is π/2 √ G, and thus asymptotically constant (and arbitrarily small!).From (8), we also see that the success probability is periodic with a period of 2t * .Now let's prove that the peak in success probability is narrow by finding its width, thus proving all our observations about figure 4. Using (7), the difference in time at which the success probability reaches a height of 1 The tan −1 makes it difficult to determine the scaling with N , so we Taylor expand it: When G = N κ , the first term scales as Θ(N 1/2 ) when κ ≤ −1 and Θ(N −1/2−κ ) when κ > −1, for large N .To determine whether keeping this first term alone is sufficient, we use Taylor's remainder theorem to bound the error which has the same scaling for large N as the first term in the Taylor series for ∆t.Thus it suffices to keep only the first term.
For constant G, the width in success probability is Θ(1/ √ N ), which agrees with our observation from figure 4 that the peak in success probability is increasingly narrow as N increases.Thus we must measure the system with increasing time precision.This behavior is opposite the linear case.That is, when G = 0 the width is Θ( √ N ), so the time at which we measure the result can be increasingly imprecise as N increases.

Time-Measurement Precision
This time-measurement precision requirement of the nonlinear algorithm requires additional resources.In particular, time and frequency standards are currently defined by atomic clocks, such as NIST-F1 in the United States [17].An atomic clock with n clock ions used as atomic oscillators has a time-measurement precision of 1/ √ n clock when the ions are acted upon independently.This can be improved using quantum entanglement, reducing the time-measurement precision to 1/n clock [18,19].Even with this improvement, our constant-time nonlinear search algorithm would require O( √ N ) ions in an atomic clock to have sufficiently high time-measurement precision to measure the peak in success probability.So, although our nonlinear algorithm runs in constant time, the total resource requirement is still O( √ N ), the same as the linear algorithm.This raises the possibility that nonlinear quantum mechanics may not provide efficient solutions to NP-complete and #P problems when all the resource requirements are taken into consideration [1].
In our case, however, we can settle for a smaller improvement in runtime and reduce the time-measurement precision and total resource requirement.If we let G decrease as N κ for κ ≤ 0, then the runtime is t * = Θ(N −κ/2 ), and the time-measurement precision is ∆t = Θ(N −1/2−κ ), where we've assumed for both that κ > −1, since for κ ≤ −1, ∆t = Θ(N 1/2 ), independently of G.This time-measurement precision requires O(N 1/2+κ ) ions in an atomic clock that utilizes entanglement.We assume, as in the setup for Grover's algorithm, that log N qubits can be used to encode the N -dimensional Hilbert space; these should also be included in the required "space" resources.Multiplying the time and "space" requirements together, which preserves the time-space tradeoff inherent in naïve parallelization, the resulting total resource requirement takes a minimum value of O(N 1/4 log N ) when κ = −1/2 (so that the runtime is N 1/4 and the time-measurement precision is constant).The success probability as a function of time at this jointly optimized value of G is plotted in figure 5; note that the peak width is independent of N .This significant-but not unreasonable-improvement over the Θ( √ N log N ) timespace resource requirements of the linear quantum search algorithm is consistent with our expectation that a modest nonlinearity should result in a modest speedup.

Repulsive Interactions
Our nonlinear search algorithm was based on the intuition that attractive interactions speed up the accumulation of success probability.By the same intuition, repulsive interactions, where G < 0, should yield a worse runtime.Our derivation of (3) for the critical value of γ is unchanged if we flip the sign of G, so (A.9) and ( 6) are still valid for repulsive interactions.These equations yield critical points x * = 1/N , 1, and (G − 1)/N G, corresponding to minima, maxima, and stationary points, respectively.
When G > −1/(N − 1), the success probability is unhindered by the stationary point and reaches a maximum value of 1, as shown in the dashed curve of figure 6.When G < −1/(N − 1), however, reaching this maximum is precluded by the presence of a stationary point, as shown in the dashed and dot-dashed curves of figure 6.
We can explicitly prove that repulsive interations (G < 0) will underperform the linear (G = 0) algorithm.From (6), So when G < 0, the magnitude of dx/dt at a particular value of x is less than when G = 0. Then success probability will increase more slowly for repulsive interactions  3).The solid line is the linear (g = 0) case, the dashed line is the nonlinear g = −0.5 case, the dotted line is the nonlinear g = −1 case, and the dot-dashed line is the nonlinear g = −1.5 case.
(except initially, where they increase at the same rate).Thus it will underperform the linear algorithm.

Validity of the Gross-Pitaevskii Equation
Of course, the cubic nonlinearity we've exploited is not fundamental, but rather occurs in an effective description of an interacting multi-particle quantum system (e.g., a BEC).So we must include the number of particles N 0 in our resource accounting.Each particle interacts with the potential at the marked site, so in the framework of Zalka's optimality proof for Grover's algorithm [20] (generalized to continuous time [21]), there are N 0 oracles, each responding to a log N bit query.Zalka showed that the product of the space requirements and the square of the time requirements is lower bounded by N , i.e., (N 0 log N )(N 1/4 ) 2 = Ω(N ).Solving for the number of particles, N 0 = Ω(N 1/2 / log N ).This is a quantum information-theoretic lower bound on the number of particles necessary for the Gross-Pitaevskii equation to be the correct asymptotic description of the multi-particle (linear) quantum dynamics.
Notice that once we account for the scaling of N 0 in the space requirements, the product of the time and space requirements is O(N 3/4 ), worse than the O(N 1/2 log N ) of Grover's algorithm.In fact, if we calculate for the general case G = N κ , where κ need not be chosen to optimize the product of the time and space (ignoring N 0 ) resources, Zalka's bound implies N 0 = Ω(max{1, N 1+κ / log N }), so the total time-space requirements are O(N 1+κ/2 ) for κ > −1, and O(N 1/2 log N ) when κ = −1.This is optimized for κ = −1, i.e., by Grover's algorithm.On the other hand, Zalka's bound is strongest when κ = 0, in which case it implies that N 0 = Ω(N/ log N ).That is, the existence of the constant time nonlinear algorithm we found in section 4 implies this stronger lower bound on N 0 , despite the O(N 1/2 ) number of clock ions required.To our knowledge, this is the first lower bound derived on the scaling of N 0 required for the Gross-Pitaevskii equation be a good asymptotic approximation.This bound also is significantly stronger than the bound implied by the physically plausible requirement that the volume of the multi-particle condensate, and thus N 0 , be of at least the order of the volume of space in which the N possible discrete locations are defined.Were we working in any fixed, finite dimension, e.g., on a cubic lattice, the volume would be proportional to N , implying N 0 = Ω(N ).But we are not; the complete graph with equal pairwise transition rates is realized by the vertices and edges of an equilateral (N − 1)-dimensional simplex.With edges of length 1, this has volume N/2 N −1 /(N − 1)!, which is much smaller than N , and also much smaller than our bound of N/ log N .

Critical Gamma is a Continuous Rescaling of Time
We previously derived the critical value of γ so that the eigenstates of the Hamiltonian are proportional to ±|w + |s .Now we examine what the critical value of γ does from another perspective.Recall the "Hamiltonian" we've been using is Recall γ = γ c is chosen according to (3): Then the Hamiltonian becomes The last term continuously redefines the "zero" of energy, so we can drop it.That is, it only changes the overall phase of the system, which has no measurable effect.Then the Hamiltonian is Importantly, H FG = −|s s| − |w w| is the Hamiltonian from Farhi and Gutmann's "analog analogue" of Grover's algorithm [13], and it is optimal.Our nonlinear algorithm has a factor of γN , so it effectively follows their optimal algorithm, but with a continuously rescaled time.That is, the system evolves according to i dψ γN dt = H FG ψ.
Let's call the rescaled time τ (t) so that dτ = γN dt.Then This has success probability given by ( 11) of [13]: Plugging in for τ , Since γ c N = 1 + Gδ = 1 − G + GN x, we get This integral transcendental equation gives x(t).While the difficulty of solving this equation makes it less useful in practice, it does reveal our nonlinear algorithm's relationship with the linear, optimal algorithm.In particular, a different control policy for γ will cause the system to evolve along a different, slower path.While not a proof, this is an argument for the optimality of our algorithm.

Multiple Marked States
Our analysis naturally extends to the case of k marked states.Let M be the set of marked basis states.As before, the system evolves in a two-dimensional subspace:

Figure 2 . 2 Figure 3 .
Figure 2. Success probability as a function of time for N = 1024 and γ = 1/N constant.The solid curve is the linear (g = 0) case, and the dashed curve is the nonlinear g = 1 case.

Figure 4 .
Figure 4. Success probability as a function of time for nonlinear search with G = 1 and γ = γ c as defined in (3).The solid line is N = 100 and the dashed line is N = 1000.

Figure 5 .
Figure 5. Success probability as a function of time for nonlinear search with G = N −1/2 and γ = γ c as defined in (3).The solid line is N = 100 and the dashed line is N = 1000.The peaks have same width, independent of N .

Figure 6 .
Figure 6.Success probability as a function of time for N = 1024 and γ = γ c as defined in (3).The solid line is the linear (g = 0) case, the dashed line is the nonlinear g = −0.5 case, the dotted line is the nonlinear g = −1 case, and the dot-dashed line is the nonlinear g = −1.5 case.

τ
= γN dt, and the equation of motion becomes i dψ dτ = H FG ψ.