Searching via Nonlinear Quantum Walk on the 2D-Grid

: We provide numerical evidence that the nonlinear searching algorithm introduced by Wong and Meyer , rephrased in terms of quantum walks with effective nonlinear phase, can be extended to the ﬁnite 2-dimensional grid, keeping the same computational advantage with respect to the classical algorithms. For this purpose, we have considered the free lattice Hamiltonian, with linear dispersion relation introduced by Childs and Ge The numerical simulations showed that the walker ﬁnds the marked vertex in O ( N 1/4 log 3/4 N ) steps, with probability O ( 1/ log N ) , for an overall complexity of O ( N 1/4 log 5/4 N ) , using amplitude ampliﬁcation. We also proved that there exists an optimal choice of the walker parameters to avoid the time measurement precision affecting the complexity searching time of the algorithm.


Introduction
Searching for an element among an unstructured database of size N takes O(N) iterations, resulting in a linear complexity time. In 1996, L. Grover came up with a quantum algorithm that speeds up this brute force O(N) problem into a O( √ N) problem [1]. The algorithm comes in many variants and has been rephrased in many ways, including quantum walks [2]. Quantum walks (QW) are essentially local unitary gates that drive the evolution of a particle on a graph [3], and although they may appear defined in a discrete and in a continuous time setting, it has been recently shown that a new family of "plastic" QW unifies and encompasses both systems [4,5]. They have been used as a mathematical framework to express many quantum algorithms, e.g., [6][7][8], but also many quantum simulation schemes e.g., [9,10]. In particular, it has been shown that many of these QWs admit, as their continuum limit, the Dirac equation [11,12], providing 'quantum simulation schemes', for the future quantum computers, to simulate all free spin-1/2 fermions. More interestingly, it has been recently proven by one of the authors, that the Grover algorithm is indeed a naturally occurring phenomenon, i.e., spontaneously implemented by some kind of particles in nature [13] over arbitrary surfaces with topological defects. From a theoretical perspective, a Grover search on a graph, rephrased in terms of QW, is an alternation of a diffusion step and an oracle step. The nodes of the graph represent elements of the configuration space of a problem, and whose edges represent the existence of a local transformation between two configurations. So far, the QW search has only been used to look for 'marked nodes', i.e., good configurations within the configuration space, as specified by an oracle. In [13], it has been proven that instead of using them to look for 'good' solutions within the configuration space of a problem, we could use them to look for topological properties of the entire configuration space.
The generalisation to an interacting multi-walkers scenario has already been explored in quantum algorithmics [14], showing that such systems are capable of universal quantum computation. Moreover, there exist many physical systems which may be described by a nonlinear effective equation, such as Bose-Einstein condensates (BEC) [15,16]. Quite interestingly, the experimental setup proposed by Alberti and Windemberg in [17] showed that a spinor BEC can simulate, in quasi-momentum space, a 1D QW with weak non-linearities. Notice that such results do not aim to suggest a scheme for implementing nonlinear quantum mechanics, which we know is linear. However, it paves the way to simulate, via a mean-field approach, one particle nonlinear dynamics, where nonlinearities come from the weakly interacting BEC. In the following, in order to avoid confusion, we will introduce nonlinearities in the very same spirit of [18]. Abrams and Lloyd consider a hypothetical nonlinear quantum mechanics equation, i.e., a nonlinear Schrödinger equation for closed single quantum mechanical systems. In particular, they show that nonlinearity in quantum computation could make quantum systems solve NP-complete problems in polynomial time.
Nonlinearity in quantum walks has been considered in several recent studies and may appear either under the form of nonlinear phases (e.g., in Kerr medium [19,20]) or via a feed-forward quantum diffusion operator [21]. In this manuscript, we will give numerical and analytical evidence that nonlinearity leads to a clear computational advantage on the two dimensional grid with respect to the linear case, consistently with previous results on complete graph [22], paving the way to extend our nonlinear scheme to higher dimensional physical dimensional grids.
The manuscript is organized as follows: In Section 2, we will introduce the QW in continuous time and discrete space; in Section 3, we will present our numerical simulations for the linear and nonlinear algorithm, and in Section 4, we will derive analytically the searching time. Finally, in Section 5, we discuss the results and conclude.

The Linear Algorithm
Let us consider a quantum walk in continuous time, over a two dimensional grid of size L × L, where N = L 2 is the number of vertices v ∈ V = {1, ..., N}. The Hilbert space of the quantum walk is spanned by the set of vertices, such that: A general state of the walker |ψ(t) is a vector, whose complex components ψ v = v ψ evolve according to the Schrödinger equation: where H vw are the coefficients of the Hamiltonian. We consider that H is local, or in other terms that H vw = 0, if and only if, v and w are adjacent. In this framework, the search problem corresponds to find a marked vertex |v , given an initial state |ψ(t = 0) , usually constructed as a uniform superposition over all vertices : The main idea is to let the quantum walker evolve on the grid for such a time that a measure of his state, on the computational basis, can be used to find the marked element with constant probability. Thus, we need to find the shortest time T, such that the probability p(t) = | v| e −iHt |ψ(0) | 2 is maximised. In the following, we refer to T as the searching time andp = p(T) as the success probability. In a continuous time setting, in order to search a single marked node, we let the walker evolve under the action of the following Hamiltonian: where H 0 is the lattice free Hamiltonian of the walker and H oracle is the oracle Hamiltonian which is used to perturb the walker free Hamiltonian upon the marked vertex. In the original work by Childs and Goldstone [2], it has been proven that, for a two-dimensional lattice, the quadratic speed up with respect to the classical search algorithm is lost. Later, the same authors proved in [23] that by introducing a second register, the coin state, as in the discrete time quantum walk, the quantum speedup O( √ N log N) is achieved even for a two-dimensional lattice. The reason has been investigated in [2], where again Childs and Goldston, inspired by Ambainis, Kempe and Rivosh [24] pointed out that the optimality loss in two spatial dimension for a continuous time quantum walk came from the quadratic dispersion relation of the free Hamiltonian. Since that, several other results confirmed this interpretation, as for example in [25], where a coinless walker with linear dispersion relation over a honeycomb lattice led to the same quantum speedup. In this paper, we will use the same model introduced by Childs and Ge [26], where a linear dispersion law is achieved introducing periodic inhomogeneities into the lattice Hamiltonian: the idea is to group, periodically, multiple neighbouring vertices into cells, as in Figure 1. In order to compare this model with the coined quantum walk, similarly to the staggered fermion model [27], we can notice that the coin states may be seen as embedded into the cells. The dimension of the coinspace is mimicked by the number of nodes grouped in each cell. Although this construct may suggest that now the searching algorithm will look for a marked cell and no longer for a specific vertex, we will show that an opportune choice of the oracle H oracle allows us to look for any single vertex as a possible marked item, keeping the advantage of the Hamiltonian over the cells. Let us now review more formally the model, as introduced in [26].
Consider the lattice partitioned in n := N/4 cells, as in Figure 1, with even L and periodic boundary conditions. The free lattice Hamiltonian, with e i , the unit vector in the i-direction, is translation invariant of length 2, and, due to the periodic conditions, the above Hamiltonian is well defined on the border of the lattice. More formally: Because, the 2D-grid is now factorised into cells, each composed by four vertices, it is convenient to introduce a new coordinates system. The cell will be labeled by r, such that r ∈ [l] 2 with l = L/2, and each of the four internal vertices will be labeled by σ ∈ Z 2 2 . In particular, we can write the following transformations, which map the vertices coordinates into cell ones : where . denoting the floor function. Now, the Hilbert space is spanned by |r, σ and the Hamiltonian reads: where In this new setting, the searched vertex is : with w ∈ [l] 2 and α ∈ Z 2 2 . We now define the oracle Hamiltonian, previously introduced in [26], by: Notice that the above oracle Hamiltonian, H oracle , added to H 0 , has the effect of disconnecting the marked vertex from its neighborhood, as in Figure 1. Indeed, by the definition, for any v, we have v| H 0 |v = w, α| H 0 |w, α = 0, and it follows that As a consequence, one cannot reach the neighbors of |w, α from itself. Moreover, one cannot reach |w, α from its neighbors, neither because the total Hamiltonian H L is real and then it has to be symmetric.
Note that this oracle Hamiltonian is different from the one usually used, for example in [2], H oracle ∝ |w, α w, α|. Indeed, it has been proven in [25,26,28] that this redefinition is justified by the 2 2 cell-partition of the lattice and the consequent linear behavior of the dispersion relation near the ground state of H 0 . Moreover, the disconnection of the marked vertex from the rest of the lattice and the on-site potentials at the neighbours of the marked item is reminiscent of the natural occurring searching algorithm studied in [13], in the discrete time setting. Now, following Equation (12), looking at the marked vertex |w, α , coincides to find a neighbour of |w, α with constant probability, or with a sufficiently large probability that can be amplified with reasonable computational overhead. In other terms, we have to search the following state: Childs and Ge [26] proved that the overlap of |Γ with the evolved state e −iHT |s is Ω(1/ log(N)), with T = O( N log(N)), where the initial state |s reads: Let us shortly discuss this result, redirecting the interested reader to the detailed calculations in [26]. In order to analyse the algorithm, we have to compute the overlap | Γ| e −iHT |s | and the searching time T for which it is maximised. By analysing the spectrum of the Hamiltonian H, we can prove that H has two reals eigenvalues: with the corresponding eigenstates |ψ ± : where √ R ± := w, α| H 0 |ψ ± = 0 and I 2 is the identity matrix of dimensions 2. Now, taking the overlap of each of the eigenstates, with the initial state, we obtain: where the last approximation holds for large n and has been proven rigorously in [26]. From the above equation, we derive the approximated expression for the initial state: and letting the system evolve for T = π 2|E ± | : Finally, the overlap shows that the system will be sufficiently close to the state |Γ , in a time T = π/2|E ± | = O ( N log N).
In other terms, the Equation (20) shows us that the system evolves in the two dimensional space spanned by |s and |Γ , and after a searching time T = π 2|E ± | , the unitary e −iHT rotates the state from |s to |Γ .
We can now summarize the above in the following Algorithm 1: Algorithm 1: Searching algorithm for the 2d-grid Initialization: |ψ = |s ; In order to study numerically the algorithm, we proceed as follows: (i) Prepare, the initial state as in Equation (14); (ii) Let the walker evolve with time; (iii) Quantify the searching time T(N) before the walker reaches its success probabilityp(N) of being localized in a ball of radius 1 around the marked element. Then, estimate this success probability, at fixed N ; (iv) Characterize T(N) andp(N), i.e., the way the success probability and the searching time depend upon the total number grid points.

Adding a Nonlinearity
Here, we will consider the hypothetical nonlinear quantum mechanics in the same spirit of [18]. However, let us notice that there are physical systems, in which multiple interacting quantum particles, under particular conditions, can be effectively described as a single particle obeying the nonlinear Schrödinger equation: where H 0 is the free Hamiltonian of the single particle and g is the coupling constant of the non linear systems. The nonlinear term corresponds to a Kerr law [29,30], acting as a self-potential.
One concrete example of this family of quantum systems is the Bose-Einstein condensate (BEC) [31,32], effectively described by Equation (21) in the limit of large particle number [33]. In the recent years, BECs have been realised in several experimental setups [34][35][36] and some of them inspired the experimental realisation of quantum walks [17,37]. In the following, we will consider a BEC with attractive interaction, i.e., g > 0, in the context of quantum algorithms. Wang and Meyer have previously investigated such systems applied to spatial search [22,38] on complete graphs. Although they found a reasonable improvement in time complexity over the linear case, their algorithm is not directly reducible to the 2-dimensional grid, keeping the same quantum speedup with respect to the classical case. In fact, Wong and Meyer used the same linear Hamiltonian introduced by Farhi and Gutmann [39], which is not optimal on the square grid, as rigorously proven in [2]. In the following, we will show that the same nonlinear spatial search introduced by Wang and Meyer can be optimal, even in 2D, considering the free Hamiltonian (9). The basic idea is to modify Equation (4), including a self-potential of the kind: where the coupling positive constant g amplifies the accumulation of the probability at the marked state due to H oracle , speeding up the search. The overall time-dependent Hamiltonian now reads: Similarly to the linear algorithm, here we have to compute the success probability p, or equivalently we have to find a searching time for which the overlap between Γ and the evolved state ψ(t) is maximised. Before doing so, we must ensure that the non-linear contribution is as reasonable as possible. In other words, we wonder whether it is possible for the system to continue to evolve in the subspace generated by |Γ and |s . We will see that this is possible at the cost of rescaling the spectrum at fixed N.
Let us look into the subspace spanned by the two orthogonal vectors {|Γ , |s }, in which |ψ(t) can be decomposed as follows: where now, the probability amplitudes are time-dependent. Moreover, Γ and s are eigenstates of H NL , in fact: and, In this subspace, the nonlinear Schrödinger equation reads: where we have dropped t parameter from the amplitudes to lighten the notation. The total Hamiltonian is represented by the following matrix: where E = |E ± |.
Remark that, by adding g|a| 2 4 to the diagonal coefficients, we only introduce a global phase, which does not change the dynamics. Thus, the rescaled total Hamiltonian is : with The eingenvalues of the above matrix are : with eigenstates, respectively: Notice that, if the perturbation gδ is sufficiently small and we rescale the spectrum, we can recover the eigenstates of the form (±1, 1) and thus, force the system to oscillate between them, similarly to the linear case. Indeed, let us multiply the Hamiltonian H L by a term (1 + c(N)gδ), with c(N) being a real function of N. The eigenvalue E now transforms as follows: In order to make the eigenstates (32) approximately of the form (±1, 1), we have to impose 2E(1 + c(N)gδ) gδ or equivalently Moreover, an upper bound of c(N) can be obtained as follows: since |ψ(t = 0) = |s and |β(t)| 2 = | s ψ(t)| 2 . Thus, Finally, keeping only the leading term 1/2E in Equation (34), we can bound c(N) as follows Notice that the above inequality implies that g cannot scale faster than √ N. Finally, we can summarise the nonlinear Algorithm 2 as follows: Algorithm 2: Searching with a nonlinear algorithm for the 2d-grid Initialization: |ψ = |s ; Contrary to the linear Algorithm 1, an exact analytical treatment to explicitly calculate T , for the Algorithm 2, is difficult. In the following, first, we will provide strong numerical evidence that such a search scheme allows a clear temporal advantage over the linear case, deriving numerically T and second, we will give an approximate analytical proof for it. In conclusion, we will discuss the time measurement precision and how it affects the choice of the parameters. Figure 2 the searching time, T as a function of c at fixed N and, in Figure 3, the corresponding success probability for several values of c. In all cases, we have chosen g = log N/π, which ensures that gδ ∼ 1. Notice that, for c = 0, the searching time T is the same as for the linear algorithm, T, but the success probability is lower, because the initial eigenstate |s is not rotated to the state |Γ , as in the linear case. In Figure 2, we can also remark that the probability peak becomes narrower, the more c increases, which would require more attempts to estimate the maximum of the probability curve, making the algorithm sub-optimal [22], as we discuss later. These considerations suggest to choose a value for c which, while ensuring the good oscillatory behaviour of the system, leaves the peaks wide enough. Looking at Figure 3, a reasonable choice is

Numerical simulations show in
(whose value is 5.52 for N = 900 in Figure 3). For the following, we will keep the above choice for our analysis. Now, we run both algorithms, the linear and the nonlinear ones, and we derive from both, a numerical characterization for the searching times T(N) (linear), T (N) (nonlinear), and the success probabilitiesp(N). The first peak probability in the linear and in the nonlinear case behaves as 1/ log(N), as is shown in Figure 4, approximately with the same pre-factor ∼4.3.

Scale Analysis
Although an exact analytical description is prohibitive, in this section, we will provide some elements of analysis that will convince the reader of the robustness of the numerical results obtained in the previous section. Let's start from the eigenvalues given in Equation (31), and rescale them using Equation (33): Setting g(N) = log N/π and c(N) = N log N 16π , as in our numerical simulation, we can clearly distinguish two different behaviours for the probability |a(t)| 2 : • A strictly linear regime, when δ(t) = 0, which appears periodically with a period • A nonlinear regime, when δ(t) is maximum. In this case, (1 + cgδ(t)) ∼ cgδ(t) = O( N log N) and consequently the period is constant: Now, we have to characterize how the system goes from the first regime to the second one and under which conditions. With no loss of generality, and because we are solely interested in the first probability peak behaviour, we focus on the first half period of the dynamics. We propose the following ansatz: where A, C 0 and C 1 are three real constants. Notice that the above self-consistent equation describes both regimes for the probability |a(t)| 2 : indeed, when |a(t = 0)| 2 = 0, we recover the Equation (20). By choosing C 0 = 1 and A = 1, and as soon as |a| 2 increases, the argument in the sine increases as well, speeding up the dynamics, as observed in Figure 3. As T 1 T 0 , in the limit of large N, there exists a critical time t s , at which the Equation (43) coincides with |a(t)| 2 = Aπ log N sin 2 C 1 |a(t)| 2 t T 1 .
We may argue that t s scales in the same way as the searching time T with N. Numerical simulations show that this is approximately true and it is shown in Figure 6 for several values of c. Along the first half period, the ansatz works surprisingly well. Moreover, if we compute the characteristic time t s at which the system transits from the first behavior with period T 0 to the second one with period T 1 , we derive exactly the same scaling laws, we have for T and T . Indeed, at t s , we have |α(t s )| 2 t s /T 1 ∼ t s /T 0 and it follows: In conclusion, by keeping only the first term in the development of arcsin x: and using Equations (41) and (42), we finally recover which is consistent with the fits realized in the previous section and confirms the time complexity of the nonlinear algorithm that we numerically assessed in Section 3.

Discussion
We provide numerical evidence that the nonlinear searching introduced by Wong and Meyer [22] can be extended to the finite 2-dimensional grid, keeping the same computational advantage with respect to the classical algorithm. For this purpose, we have considered a different Hamiltonian, admitting a linear dispersion relation, as proven in [26]. The numerical simulations showed that the walker finds the marked vertex in O(N 1/4 log 3/4 N) steps with probability O(1/ log N), for an overall complexity of O(N 1/4 log 5/4 N), using amplitude amplification. These results are consistent with those proven by Wong and Meyer in [22] using a nonlinear Schrödinger equation on a complete graph. However, the optimality of the above algorithm strongly depends on the time measurement precision. As we mentioned earlier, the width of the probability peaks may affect the complexity of the algorithm, as already discussed in [22]. The narrower the peak, the less efficient the algorithm will be. The optimality of the nonlinear algorithm depends on the choice of g and c. For instance, if we had chosen both c and g in O( √ N), we would end up with a runtime of t s = O(log 5/4 N). However, the period T 1 of the Equation (42) would no longer be a constant, it would rather decrease with N: The peak width at half maximum, shown in Figure 7, is given by : which yields to ∆t ∼ T 2 .
Now, considering that n ions are required in an atomic clock to obtain a time precision of O(1/n) [40], we have to require an extra O(1/T 2 ) temporal resource, which dramatically affects the total complexity, making us lose any kind of advantage over the linear algorithm.
In conclusion, our work presents several advantages with respect to the previous ones: (i) quantum walks are easily implementable in our labs by several physical systems; (ii) graphs having sets of vertices of constant degree, are more natural and pave the way to a 3-dimensional generalisation. An extension of the above results to other kind of nonlinearities will deserve future investigation.