<figure><img src="image_1.png" /></figure>

HAL
open science

# Searching via Nonlinear Quantum Walk on the 2D-Grid

Di Molfetta Giuseppe, Basile Herzog

► To cite this version:

Di Molfetta Giuseppe, Basile Herzog. Searching via Nonlinear Quantum Walk on the 2D-Grid. Algorithms, 2020, 13 (11), pp.305. 10.3390/a13110305. hal-03594725

HAL Id: hal-03594725

https://amu.hal.science/hal-03594725

Submitted on 26 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Searching via nonlinear quantum walk on the 2D-grid

Giuseppe Di Molfetta <fnref n="1" />,* and Basile Herzog <fnref n="1" />,<fnref n="2" />

Received: date; Accepted: date; Published: date

<figure><img src="image_4.png" /></figure>

check for
updates

**Abstract:** We provide numerical evidence that the nonlinear searching algorithm introduced by Wong and Meyer [1], rephrased in terms of quantum walks with effective nonlinear phase, can be extended to the finite 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 [2]. 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. We also proved that there exists an optimal choice of the walker parameters to avoid that the time measurement precision affects the complexity searching time of the algorithm.

## 1. Introduction

Searching 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 [3]. The algorithm comes in many variants and has been rephrased in many ways, including quantum walks [4]. Quantum walks (QW) are essentially local unitary gates that drive the evolution of a particle on a graph [5], 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 unify and encompass both systems [6,7]. They have been used as a mathematical framework to express many quantum algorithms, e.g. [8–10], but also many quantum simulation schemes e.g. [11,12]. In particular, it has been shown that many of these QW admit, as their continuum limit, the Dirac equation [13,14], 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 [15] 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 [15], it has been proved 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 [16], showing that such systems are capable of universal quantum computation. Moreover, there exists many physical systems which may be described by a nonlinear *effective* equation, such as

<footnote n="1">¹ Aix-Marseille Université, CNRS, LIS, Marseille, France; e-mail@e-mail.com</footnote>

<footnote n="2">² Université de Lorraine, LPCT, Nancy, France; e-mail@e-mail.com</footnote>

<footnote>* giuseppe.dimolfetta@lis-lab.fr</footnote>
Bose-Einstein condensates (BEC)[17,18]. Quite interestingly, the experimental setup proposed by Alberti and Windemberg in [19], showed that a spinor BEC can simulate, in quasi-momentum space, a 1D QW with weak non-linearities. Notice that such results, does 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 comes from the weakly interacting BEC. In the following, in order to avoid confusion, we will introduce nonlinearities in the very same spirit of [20]. Abrams and Lloyd consider a hypothetical nonlinear quantum mechanics, 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 have been considered in several recent studies and may appear either under the form of nonlinear phases, (e.g. in Kerr medium [21] and [22]) or via a feed-forward quantum diffusion operator [23]. 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 [24], 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.

## 2. Model

### 2.1. The linear algorithm

Let us consider a quantum walk in continuous time, over a two dimensional grid of size $L \times L$, where $N = L^2$ is the number of vertices $v \in V = \{1, ..., N\}$. The Hilbert space of the quantum walk is spanned by the set of vertices, such that:

$$ \mathcal{H} := \operatorname{span}\{|v\rangle : v \in V\}. \qquad (1) $$

A general state of the walker $|\psi(t)\rangle$ is a vector, whose complex components $\psi_v = \langle v|\psi$, evolve according the Schrödinger equation:

$$ i \frac{d\psi_v(t)}{dt} = \sum_w H_{vw} \psi_w(t), \qquad (2) $$

where $H_{vw}$ are the coefficients of the Hamiltonian. We consider that $H$ is local, or in other terms that $H_{vw} \neq 0$ if and only if $v$ and $w$ are adjacent.

In this framework, the search problem corresponds to find a marked vertex $|\bar{v}\rangle$, given an initial state $|\psi(t=0)\rangle$, usually constructed as a uniform superposition over all vertices :

$$ |\psi(0)\rangle = \frac{1}{\sqrt{N}} \sum_{v \in V} |v\rangle. \qquad (3) $$

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) = |\langle \bar{v} | e^{-iHt} | \psi(0) \rangle|^2$ is maximised. In the following, we refer to $T$, as the searching time and $\bar{p} = p(T)$, 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:

$$ H_L = H_0 + H_{\text{oracle}}, \qquad (4) $$
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 [25], 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 [26] that by introducing a second register, the coin state, as in the discrete time quantum walk, the quantum speedup $O(\sqrt{N} \log N)$ is achieved even for a two dimensional lattice. The reason has been investigated in [4], where again Childs and Goldston, inspired by Ambainis, Kempe and Rivosh [27], 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 [28], 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 [2], 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 Fig.1. In order to compare this model with the coined quantum walk, similarly to the staggered fermion model [29], 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 possible marked item, keeping the advantage of the Hamiltonian over the cells. Let us now review more formally the model, as introduced in [2].

Consider the lattice partitioned in $n := N/4$ cells, as in Fig.1, with even $L$ and periodic boundary conditions. The free lattice Hamiltonian,

$$H_0 |v\rangle = (-1)^{v_x} (|v + e_x\rangle - |v - e_x\rangle) \\ + (-1)^{v_x+v_y} (|v + e_y\rangle - |v - e_y\rangle) \qquad (5)$$

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 :

$$[H_0, T_i] = 0 \quad T_i |v\rangle = |v + 2e_i\rangle. \qquad (6)$$

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 \in [l]^2$ with $l = L/2$, and each of the four internal vertices will be labeled by $\sigma \in \mathbb{Z}_2^2$. In particular, we can write the following transformations, which map the vertices coordinates into cell ones :

$$r_i = \lfloor \frac{v_i}{2} \rfloor, \qquad (7)$$

$$\sigma_i = v_i - 2r_i \qquad (8)$$

where $\lfloor \cdot \rfloor$ denotes the floor function. Now, the Hilbert space is spanned by $|r, \sigma\rangle$ and the Hamiltonian reads:

$$H_0 |r, \sigma\rangle = \sum_{i=x,y} (-1)^{s_i(\sigma)} (|r + \sigma_i e_i, \sigma + e_i\rangle - |r - \bar{\sigma}_i e_i, \sigma + e_i\rangle) \qquad (9)$$
where $s_x(\sigma) = \sigma_x, s_y(\sigma) = \sigma_y + \sigma_y$ and $\bar{\sigma}_i = 1 - \sigma_i$.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** The 2d-grid and its factorisation into cells (grey squares). The marked vertex (circled) is disconnected from its neighbors.</figcaption></figure>

In this new setting, the searched vertex is:

$$|\bar{\sigma}\rangle = |w, \alpha\rangle \qquad (10)$$

with $w \in [l]^2$ and $\alpha \in \mathbb{Z}_2^2$.

We now define the oracle Hamiltonian, previously introduced in [2], by:

$$H_{\text{oracle}} = -|w, \alpha\rangle \langle w, \alpha| H_0 - H_0 |w, \alpha\rangle \langle w, \alpha| \qquad (11)$$

Notice that the above oracle Hamiltonian, $H_{\text{oracle}}$, added to $H_0$, has the effect of disconnecting the marked vertex from its neighborhood, as in Fig. 1. Indeed, by the definition, for any $v$, we have $\langle v| H_0 |v\rangle = \langle w, \alpha| H_0 |w, \alpha\rangle = 0$, and it follows that

$$H_L |w, \alpha\rangle = (H_{\text{oracle}} + H_0) |w, \alpha\rangle = -|w, \alpha\rangle \langle w, \alpha| \quad H_0 \quad |w, \alpha\rangle - H_0 |w, \alpha\rangle + H_0 |w, \alpha\rangle = 0. \qquad (12)$$

As consequence, one can not reach the neighbors of $|w, \alpha\rangle$ from itself. Moreover, one can not reach $|w, \alpha\rangle$ 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 [25], $H_{\text{oracle}} \propto |w, \alpha\rangle \langle w, \alpha|$. Indeed, it has been proven in [2,30,31], that this redefinition is justified by the 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 [15], in the discrete time setting.

Now, following Eq. (12), looking at the marked vertex $|w, \alpha\rangle$, coincides to find a neighbour of $|w, \alpha\rangle$ 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:

$$|\Gamma\rangle = \frac{1}{2} H_0 |w, \alpha\rangle. \qquad (13)$$
Childs and Ge [2] proved that the overlap of $|\Gamma\rangle$ with the evolved state $e^{-iHT}|s\rangle$ is $\Omega(1/\sqrt{\log(N)}),$ with $T = O(\sqrt{N \log(N)})$, where the initial state $|s\rangle$ reads :

$$s := \sqrt{\frac{4}{N}} \sum_{r} |r, \alpha\rangle. \qquad (14)$$

Let us shortly discuss this result, redirecting the interested reader to the detailed calculations in [2]. In order to analyse the algorithm, we have to compute the overlap $|\langle\Gamma|e^{-iHT}|s\rangle|$ 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:

$$E_{\pm} \approx \pm \sqrt{\frac{16\pi}{N \log N}}, \qquad (15)$$

with the corresponding eigenstates $|\psi_{\pm}\rangle$:

$$|\psi_{\pm}\rangle = \sqrt{R_{\pm}}(H_0 - \mathbb{I}_2 E_{\pm})^{-1} |w, \alpha\rangle \qquad (16)$$

where $\sqrt{R_{\pm}} := \langle w, \alpha | H_0 | \psi_{\pm} \rangle \neq 0$ and $\mathbb{I}_2$ is the Identity matrix of dimensions 2. Now, taking the overlap of each of the eigenstates, with the initial state we get :

$$\langle \psi_{\pm} | |s\rangle = - \frac{\sqrt{R_{\pm}}}{E_{\pm} \sqrt{N/4}} \approx \mp \frac{1}{\sqrt{2}} \qquad (17)$$

where the last approximation holds for large $n$ and has been proven rigorously in [2]. From the above equation, we derive the approximated expression for the initial state:

$$|s\rangle \approx \frac{1}{\sqrt{2}} (|\psi_{-}\rangle - |\psi_{+}\rangle) \qquad (18)$$

and letting the system evolve for $T = \frac{\pi}{2|E_{\pm}|}$:

$$e^{-iHT} |s\rangle \approx \frac{1}{\sqrt{2}} (|\psi_{-}\rangle + |\psi_{+}\rangle). \qquad (19)$$

Finally the overlap

$$|\langle\Gamma|e^{-iHt}|s\rangle|^2 \approx \frac{\pi}{\log N} \sin^2(2|E_{\pm}|t), \qquad (20)$$

shows that the system will be sufficiently close to the state $|\Gamma\rangle$, in a time $T = \pi/2|E_{\pm}| = O(\sqrt{N \log N})$. In other terms, the Eq. (20) shows us that the system evolves in the two dimensional space spanned by $|s\rangle$
and $|\Gamma\rangle$, and after a searching time $T = \frac{\pi}{2|E_{\pm}|}$, the unitary $e^{-iHT}$ rotates the state from $|s\rangle$ to $|\Gamma\rangle$. We can now summarize the above in the following algorithm:

<table><caption>**Algorithm 1:** Searching algorithm for the 2d-grid</caption><tr><td>Initialization:</td><td></td></tr><tr><td></td><td>|ψ⟩ = |s⟩;</td></tr><tr><td></td><td>t = 0;</td></tr><tr><td>while</td><td>t &lt; T = O(√N log N) do</td></tr><tr><td></td><td>|ψ⟩ = exp(-iH<sub>L</sub>Δt) |ψ⟩;</td></tr><tr><td></td><td>t = t + Δt;</td></tr><tr><td>end</td><td></td></tr><tr><td>Measure</td><td>|ψ⟩</td></tr></table>

In order to study numerically the algorithm, we proceed as follows: (i) Prepare, the initial state as in Eq. (14); (ii) Let the walker evolve with time; (iii) Quantify the searching time $T(N)$ before the walker reaches its success probability $\bar{p}(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)$ and $\bar{p}(N)$, i.e. the way the success probability and the searching time depend upon the total number grid points.

## 2.2. Adding a nonlinearity

Here we will consider a hypothetical nonlinear quantum mechanics in the same spirit of [20]. 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 to the nonlinear Schrödinger equation:

$$ \partial_t \psi(\mathbf{r}, t) = (H_0 - g|\psi(\mathbf{r}, t)|^2) \psi(\mathbf{r}, t), \quad (21) $$

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 [32,33], acting as a self-potential. One concrete example of this family of quantum systems is the Bose-Einstein condensate (BEC) [34,35], effectively described by Eq. (21) in the limit of large particle number [36]. In the recent years, BECs have been realised in several experimental setups [37-39] and some of them inspired experimental realisation of quantum walks [19,40]. 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 [24,41] 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 respect to the classical case. In fact, Wong and Meyer used the same linear Hamiltonian introduced by Farhi and Gutmann [42], which is not optimal on the square grid, as rigorously proved in [25]. 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 Eq. (4), including a self-potential of the kind:

$$ H_{NL} = -g \sum_{\mathbf{r},\sigma} |\langle \mathbf{r}, \sigma \rangle \psi(t)|^2 |\mathbf{r}, \sigma \rangle \langle \mathbf{r}, \sigma | \quad (22) $$

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:

$$ H = H_L + H_{NL}. \quad (23) $$
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 :

$$
|\psi(t)\rangle = a(t) |\Gamma\rangle + b(t) |s\rangle, \qquad (24)
$$

where now, the probability amplitudes are time-dependent. Moreover, Γ and s are eigenstates of H<sub>NL</sub>, in fact :

$$
\begin{align}
H_{NL} |s\rangle &= -g \sum_{\mathbf{r},\sigma} |\langle\mathbf{r}, \sigma\rangle s|^2 |\mathbf{r}, \sigma\rangle \langle\mathbf{r}, \sigma\rangle s \nonumber \\
&= -g \sum_{\mathbf{r},\sigma} \frac{4}{N} |\mathbf{r}, \sigma\rangle \sqrt{\frac{4}{N}} = -g \frac{4}{N} |s\rangle \tag{25}
\end{align}
$$

and,

$$
\begin{align}
H_{NL} |\Gamma\rangle &= -g \sum_{\mathbf{r},\sigma} |\langle\mathbf{r},\sigma\rangle \Gamma|^2 |\mathbf{r},\sigma\rangle \langle\mathbf{r},\sigma\rangle \Gamma \nonumber \\
&= -g \sum_{\mathbf{r},\sigma} \frac{1}{4} |\mathbf{r},\sigma\rangle \langle\mathbf{r},\sigma| \frac{H_0}{2} |\mathbf{w},\sigma\rangle = -g \frac{1}{4} |\Gamma\rangle. \tag{26}
\end{align}
$$

In this subspace, the nonlinear Schrödinger equation reads:

$$
\partial_t \begin{pmatrix} a \\ b \end{pmatrix} = -iH \begin{pmatrix} a \\ b \end{pmatrix} \qquad (27)
$$

where we have dropped *t* parameter from the amplitudes to lighten the notation.

The total Hamiltonian is represented by the following matrix:

$$
H = \begin{pmatrix} \frac{-g|a|^2}{4} & E \\ E & \frac{-4g|b|^2}{N} \end{pmatrix}, \qquad (28)
$$

where $E = |E_{\pm}|$.

Remark that, by adding $\frac{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 :

$$
H' = \begin{pmatrix} 0 & E \\ E & g\delta \end{pmatrix}, \qquad (29)
$$

with

$$
\delta(t) = \left( \frac{|a|^2}{4} - \frac{4|\beta|^2}{N} \right). \qquad (30)
$$
The eigenvalues of the above matrix are:

$$E'_{\pm} = \frac{1}{2}g\delta \pm \frac{1}{2}\sqrt{(g\delta)^2 + 4E^2} \quad (31)$$

with eigenstates, respectively:

$$|\psi'_{\pm}\rangle = \left( \frac{g\delta}{2E} \pm \frac{\sqrt{(g\delta)^2 + 4E^2}}{2E}, 1 \right). \quad (32)$$

Notice that, if the perturbation $g\delta$ is sufficiently small and we rescale the spectrum, we can recover the eigenstates of the form $(\pm 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\delta)$, with $c(N)$ a real function of $N$. The eigenvalue $E$, now transforms as follows:

$$E \rightarrow E(1 + c(N)g\delta). \quad (33)$$

In order to make the eigenstates (32) approximately of the form $(\pm 1, 1)$, we have to impose $2E(1 + c(N)g\delta) \gg g\delta$ or equivalently

$$c(N) \gg \frac{1}{2E} - \frac{1}{g\delta}. \quad (34)$$

Moreover, an upper bound of $c(N)$ can be obtained as follows:

$$\delta(0) = \frac{|a(0)|^2}{4} - \frac{4|\beta(0)|^2}{N} = -4/N \quad (35)$$

since $|\psi(t=0)\rangle = |s\rangle$ and $|\beta(t)|^2 = |\langle s \rangle \psi(t)|^2$. Thus,

$$c(N) < c_{max} = \frac{N}{4g}. \quad (36)$$

Finally, keeping only the leading term $1/2E$ in Eq. (34), we can bound $c(N)$ as follows:

$$\frac{N}{4g} > c(N) \gg \frac{1}{2} \sqrt{\frac{N \log N}{16\pi}} = \frac{1}{2E}. \quad (37)$$

Notice that, the above inequality implies that $g$ cannot scale faster than $\sqrt{N}$.

Finally, we can summarise the nonlinear algorithm as follows:

<caption>**Algorithm 2:** Searching with a nonlinear algorithm for the 2d-grid</caption>

Initialization:

$|\psi\rangle = |s\rangle$;

$t = 0$;

while $t < T'$ do

$|\psi\rangle = \exp(-i\Delta t ((1 + gc\delta) H_L + H_{NL})) |\psi\rangle$; 	$\left| t = t + \Delta t; \right|$

end

Measure $|\psi\rangle$

Contrary to the linear algorithm (1), an exact analytical treatment to explicitly calculate $T'$ 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.

### 3. Numerical results

Numerical simulations show in Fig.2 the searching time, $T'$ as a function of $c$ at fixed $N$ and, in Fig.3, the corresponding success probability for several values of $c$. In all cases, we have chosen $g = \log N/\pi$, which ensures that $g\delta \sim 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\rangle$ is not rotated to the state $|\Gamma\rangle$, as in the linear case. In Fig.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 [24], 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 Fig.3, a reasonable choice is

$$c = 1/E = \sqrt{\frac{N \log N}{16\pi}}, \qquad (38)$$

(whose value is 5.52 for $N = 900$ in Fig. 3). For the following, we will keep the above choice for our analysis.

<figure><img src="image_6.png" /><figcaption>Figure 2. Searching time of the algorithm as a function of $c$, $N = 900$.</figcaption></figure>

Now we run both algorithms, the linear and the nonlinear one and we derive from both, a numerical characterization for the searching times $T(N)$ (linear), $T'(N)$ (nonlinear), and the success probabilities $\bar{p}(N)$. The first peak probability in the linear and in the nonlinear case, behaves as $1/\log(N)$, as it is shown in Fig.4, approximately with the same pre-factor $\sim 4.3$. The Fig.5 shows the searching times, $T$, and $T'$. As expected, in the linear model, we recover $T = O(\sqrt{N \log N})$. The nonlinear case, instead, shows a significant advantage respect to the linear case, with

$$T' = O(N^{1/4} \log^{3/4} N), \qquad (39)$$

which yields an overall complexity algorithm of $O(N^{1/4} \log^{5/4} N)$, using amplitude amplification.
<figure><img src="image_7.png" /><figcaption>Figure 3. The probability $p(t) = |\langle \Gamma | e^{-iHt} | s \rangle|^2$ for different values of $c - N = 900$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. Success probability as a function of N for the nonlinear and linear algorithm (inset)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5. Searching times in the nonlinear and the linear algorithm (inset).</figcaption></figure>
## 4. 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 Eq. (31), and rescale them using Eq. (33):

$$E'_{\pm} = \frac{1}{2} g\delta \pm \frac{1}{2} \sqrt{(g\delta)^2 + 4E^2(1 + c g\delta)^2} \quad (40)$$

Setting $g(N) = \log N/\pi$ and $c(N) = \sqrt{N \log N}/16\pi$, as in our numerical simulation, we can clearly distinguish two different behaviours for the probability $|a(t)|^2$:

* A strictly linear regime, when $\delta(t) = 0$, which appears periodically with a period

$$T_0 \sim 1/\Delta\lambda \sim \frac{1}{E} = O(\sqrt{N \log N}) \quad (41)$$

* A nonlinear regime, when $\delta(t)$ is maximum. In this case, $(1 + c g\delta(t)) \sim c g\delta(t) = O(\sqrt{N \log N})$ and consequently the period is constant:

$$T_1 \sim 1/\Delta\lambda = O(1). \quad (42)$$

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 to the first probability peak behaviour, we focus on the first half period of the dynamics. We propose the following ansatz:

$$|a(t)|^2 = \frac{A\pi}{\log N} \sin^2 \left( \frac{C_0 t}{T_0} + C_1 |a(t)|^2 \frac{t}{T_1} \right) \quad (43)$$

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 Eq. (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 on Fig. 3. As $T_1 \ll T_0$, in the limit of large $N$, there exists a critical time $t_s$, at which the Eq. (43) coincides with $|a(t)|^2 = \frac{A\pi}{\log N} \sin^2 \left( C_1 |a(t)|^2 \frac{t}{T_1} \right)$.

We may argue that $t_s$ scales in the same way than the searching time $T'$ with $N$. Numerical simulations shows that this is approximately true and it is shown in Fig. 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 $|a(t_s)|^2 t_s / T_1 \sim t_s / T_0$ and it follows:

$$\frac{1}{\log N} \sin^2(t_s/T_0) \sim T_1/T_0. \quad (44)$$

In conclusion, by keeping only the first term in the development of $\arcsin x$:

$$t_s = O(\sqrt{T_1 T_0 \log N}) \quad (45)$$
<figure><img src="image_10.png" /><figcaption>Figure 6. N = 900 - Probability over |Γ⟩ for different values of c and numerical resolution of eq. (43).</figcaption></figure>

and using Eqs. (41) and (42), we finally recover

$$t_s = O(N^{1/4} \log^{3/4} N), \qquad (46)$$

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 Sec. 3.

## 5. Discussion

We provide numerical evidence that the nonlinear searching introduced by Wong and Meyer [1] 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 proved in [2]. 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 proved by Wong and Meyer in [24] 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 [24]. 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 choose both c and g in $O(\sqrt{N})$, we would end up
with a runtime of $t_s = O(\log^{5/4} N)$. However, the period $T_1$ of the Eq. 42 would no longer be a constant, it would rather decrease with $N$:

$$T_2 = O\left(\frac{\log N}{\sqrt{N}}\right). \qquad (47)$$

The peak width at half maximum, shown in Fig. 7, is given by:

$$\sin^2(\Delta t/2T_2) = \frac{1}{2}, \qquad (48)$$

which yields to

$$\Delta t \sim T_2. \qquad (49)$$

Now, considering that $n$ ions are required in an atomic clock to obtain a time precision of $O(1/n)$ [43],

<figure><img src="image_11.png" /><figcaption>Figure 7. Width of the probability peak - Insight: Searching time</figcaption></figure>

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.
6. Acknowledgements

The authors acknowledge inspiring conversations with Thomas Wong. This work has been funded by the Pépinière d'Excellence 2018, AMIDEX fondation, project DiTiQuS and the ID 60609 grant from the John Templeton Foundation, as part of the "The Quantum Information Structure of Spacetime (QISS)" Project.
1. Meyer, D.A.; Wong, T.G. Nonlinear quantum search using the Gross-Pitaevskii equation. New Journal of Physics **2013**, *15*, 063014.

2. Childs, A.M.; Ge, Y. Spatial search by continuous-time quantum walks on crystal lattices. Physical Review A **2014**, *89*. doi:10.1103/physreva.89.052337.

3. Grover, L.K. A fast quantum mechanical algorithm for database search, 1996, [arXiv:quant-ph/quant-ph/9605043].

4. Childs, A.M.; Goldstone, J. Spatial search by quantum walk. Physical Review A **2004**, *70*, 022314.

5. Venegas-Andraca, S.E. Quantum walks: a comprehensive review. Quantum Information Processing **2012**, *11*, 1015–1106.

6. Di Molfetta, G.; Arrighi, P. A quantum walk with both a continuous-time limit and a continuous-spacetime limit. Quantum Information Processing **2020**, *19*, 47.

7. Manighalam, M.; Di Molfetta, G. Continuous Time Limit of the DTQW in 2D+ 1 and Plasticity. arXiv preprint arXiv:2007.01425 **2020**.

8. Ambainis, A. Quantum walk algorithm for element distinctness. SIAM Journal on Computing **2007**, *37*, 210–239.

9. Childs, A.M.; Cleve, R.; Deotto, E.; Farhi, E.; Gutmann, S.; Spielman, D.A. Exponential algorithmic speedup by a quantum walk. Proceedings of the thirty-fifth annual ACM symposium on Theory of computing, 2003, pp. 59–68.

10. Childs, A.M. Universal computation by quantum walk. Physical review letters **2009**, *102*, 180501.

11. Hatifi, M.; Di Molfetta, G.; Debbasch, F.; Brachet, M. Quantum walk hydrodynamics. Scientific reports **2019**, *9*, 1–7.

12. Di Molfetta, G.; Pérez, A. Quantum walks as simulators of neutrino oscillations in a vacuum and matter. New Journal of Physics **2016**, *18*, 103038.

13. Di Molfetta, G.; Brachet, M.; Debbasch, F. Quantum walks as massless Dirac fermions in curved space-time. Physical Review A **2013**, *88*, 042301.

14. Arrighi, P.; Di Molfetta, G.; Márquez-Martín, I.; Pérez, A. Dirac equation as a quantum walk over the honeycomb and triangular lattices. Physical Review A **2018**, *97*, 062111.

15. Roget, M.; Guillet, S.; Arrighi, P.; Di Molfetta, G. Grover Search as a Naturally Occurring Phenomenon. Physical Review Letters **2020**, *124*, 180501.

16. Childs, A.M.; Gosset, D.; Webb, Z. Universal computation by multiparticle quantum walk. Science **2013**, *339*, 791–794.

17. Ebrahimi Kahou, M.; Feder, D.L. Quantum search with interacting Bose-Einstein condensates. Phys. Rev. A **2013**, *88*, 032310. doi:10.1103/PhysRevA.88.032310.

18. Kevrekidis, P.G.; Frantzeskakis, D.J.; Carretero-González, R. Emergent nonlinear phenomena in Bose-Einstein condensates: theory and experiment; Vol. 45, Springer Science & Business Media, 2007.

19. Alberti, A.; Wimberger, S. Quantum walk of a Bose-Einstein condensate in the Brillouin zone. Physical Review A **2017**, *96*, 023620.

20. Abrams, D.S.; Lloyd, S. Nonlinear Quantum Mechanics Implies Polynomial-Time Solution for NP-Complete and P-problems. Physical Review Letters **1998**, *81*, 3992–3995. doi:10.1103/physrevlett.81.3992.

21. Molfetta, G.D.; Debbasch, F.; Brachet, M. Nonlinear Optical Galton Board: thermalization and continuous limit, 2015, [arXiv:quant-ph/1506.04323].

22. Navarrete-Benlloch, C.; Pérez, A.; Roldán, E. Nonlinear optical Galton board. Physical Review A **2007**, *75*. doi:10.1103/physreva.75.062333.

23. Shikano, Y.; Wada, T.; Horikawa, J. Discrete-time quantum walk with feed-forward quantum coin. Scientific reports **2014**, *4*, 4427.

24. Meyer, D.A.; Wong, T.G. Nonlinear quantum search using the Gross-Pitaevskii equation. New Journal of Physics **2013**, *15*, 063014. doi:10.1088/1367-2630/15/6/063014.
25. Childs, A.M.; Goldstone, J. Spatial search by quantum walk. Physical Review A **2004**, *70*, doi:10.1103/physreva.70.022314.

26. Childs, A.M.; Goldstone, J. Spatial search and the Dirac equation. Physical Review A **2004**, *70*, 042312.

27. Ambainis, A.; Kempe, J.; Rivosh, A. Coins Make Quantum Walks Faster, 2004, [arXiv:quant-ph/quant-ph/0402107].

28. Foulger, I.; Gnutzmann, S.; Tanner, G. Quantum search on graphene lattices. Physical review letters **2014**, *112*, 070504.

29. Susskind, L. Lattice fermions. Physical Review D **1977**, *16*, 3031.

30. Foulger, I.; Gnutzmann, S.; Tanner, G. Quantum Search on Graphene Lattices. Phys. Rev. Lett. **2014**, *112*, 070504. doi:10.1103/PhysRevLett.112.070504.

31. Chakraborty, S.; Novo, L.; Roland, J. Finding a marked node on any graph via continuous-time quantum walks. Physical Review A **2020**, *102*, 022227. doi:10.1103/PhysRevA.102.022227.

32. Kohl, R.; Biswas, A.; Milovic, D.; Zerrad, E. Optical soliton perturbation in a non-Kerr law media. Optics & Laser Technology **2008**, *40*, 647–662.

33. Zhang, Z.y.; Liu, Z.h.; Miao, X.j.; Chen, Y.z. New exact solutions to the perturbed nonlinear Schrödinger's equation with Kerr law nonlinearity. Applied Mathematics and Computation **2010**, *216*, 3064–3072.

34. Einstein, A. Quantentheorie des einatomigen idealen Gases. SB Preuss. Akad. Wiss. phys.-math. Klasse **1924**.

35. Bose, S.N. Plancks gesetz und lichtquantenhypothese **1924**.

36. Erdős, L.; Schlein, B.; Yau, H.T. Derivation of the Gross-Pitaevskii equation for the dynamics of Bose-Einstein condensate. Annals of Mathematics **2010**, pp. 291–370.

37. Cornell, E.A.; Wieman, C.E. Nobel Lecture: Bose-Einstein condensation in a dilute gas, the first 70 years and some recent experiments. Reviews of Modern Physics **2002**, *74*, 875.

38. Nikuni, T.; Oshikawa, M.; Oosawa, A.; Tanaka, H. Bose-Einstein condensation of dilute magnons in TlCuCl 3. Physical review letters **2000**, *84*, 5868.

39. O'Callaghan, J. Quantum clouds in orbit, 2020.

40. Xie, D.; Deng, T.S.; Xiao, T.; Gou, W.; Chen, T.; Yi, W.; Yan, B. Topological Quantum Walks in Momentum Space with a Bose-Einstein Condensate. Phys. Rev. Lett. **2020**, *124*, 050502. doi:10.1103/PhysRevLett.124.050502.

41. Meyer, D.A.; Wong, T.G. Quantum search with general nonlinearities. Physical Review A **2014**, *89*, 012312.

42. Farhi, E.; Gutmann, S. Analog analogue of a digital quantum computation. Physical Review A **1998**, *57*, 2403.

43. Bollinger, J.J.; Itano, W.M.; Wineland, D.J.; Heinzen, D.J. Optimal frequency measurements with maximally correlated states. Phys. Rev. A **1996**, *54*, R4649–R4652. doi:10.1103/PhysRevA.54.R4649.

<figure><img src="image_12.png" /></figure>

<footnote>**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.</footnote>