Understanding and controlling N-dimensional quantum walks via dispersion relations: application to the two-dimensional and three-dimensional Grover walks—diabolical points and more

Margarida Hinarejos; Armando Pérez; Eugenio Roldán; Alejandro Romanelli; Germán J de Valcárcel

doi:10.1088/1367-2630/15/7/073041

1. Introduction

Quantum walks (QWs) are becoming widespread in physics. Originally introduced as quantum versions of random classical processes [1] and quantum cellular automata [2], they have been intensively investigated, especially in connection with quantum information science [3–9], and have been recently shown to constitute a universal model for quantum computation [10, 11]. Remarkably, in the last decade this simple quantum diffusion model has found connections with very diverse phenomena such as optical diffraction and other wave phenomena [12, 13], quantum games [14], Anderson localization and quantum chaos [15–18], and has even been identified as a generalized quantum measurement device [19]. Moreover, through proposals for and actual implementations, QWs have found connections with the dynamics of quite different physical systems such as linear optical systems and optical cavities [20–28], optical lattices [29–32], trapped ions [33–35] and Bose–Einstein condensates [36], among others. In recent years, experimental capabilities in performing QWs have increased steadily and today there are devices that, e.g., implement two-dimensional (2D) QWs with a large number of propagation steps [37]. For recent reviews on the field, see [38, 39].

Recently, some of us have been stressing on the interest of analyzing the dynamics of coined QWs from the perspective of their dispersion relation [13, 40], which has also been used for purposes different from ours (see e.g. [41, 42]). We think this is an interesting approach that is particularly suitable for non-localized initial conditions, i.e. when the initial state is described by a probability distribution extending to a finite region in the lattice. When this region is not too small, the evolution of the QW can be derived, to a good degree of approximation, from familiar linear differential wave equations. Certainly the use of extended initial conditions has been rare up to now [44], but we think it is worth stressing in that they allow one to reach final probability distributions which, interestingly, can be tailored to some extent. It is also worth mentioning that extended initial conditions make a connection with multiparticle QWs for non-interacting particles [43].

In [40], we have recently applied the above viewpoint in studying alternate QWs [45] in N dimensions. Alternate QWs are a simpler version of multidimensional QWs than the standard one [46] as they only require a single qubit as a coin (used N times per time step) [45] instead of using a 2N-qubit as the coin (used once per time step), as originally defined in [46]. Here we shall carry out, from the dispersion relation viewpoint, the study of multidimensional QWs in their standard form. Below we consider a general treatment, valid for N dimensions, and illustrate our results with the special cases of the Grover 2D and three-dimensional (3D) QWs. We note that the 2D Grover QW has received most of the theoretical attention paid to multidimensional QWs [46–58].

We shall derive continuous wave equations, whose form depends on the probability distribution at time t = 0. To leading order, these equations are partial differential equations that can be written in the form

$\begin{equation} \mathrm{i}\frac{\partial A^{\left(s\right)}({\bf X},t)}{\partial t}=-\frac{1}{2}\sum_{i,j=1}^{N}\varpi_{ij}^{\left(s\right)}\frac{\partial^{2}A^{\left(s\right)}}{\partial X_{i}\partial X_{j}}+\cdots \end{equation} \tag{ 1 }$

with $A^{\left (s\right )}$ a continuous amplitude probability, coefficients $\varpi _{ij}^{\left (s\right )}$ determined by the QW dispersion relation properties and X_i the spatial coordinates in a reference frame moving with group velocity (see below for full details). There is a lot of accumulated knowledge regarding solutions of the above continuous equations and this knowledge allows, to some extent, for qualitative knowledge of the long-term probability distribution of the QW for particular initial conditions. This qualitative knowledge is particularly appealing to us, as we are interested, not in obtaining approximate continuous solutions to the discrete QW, but in getting a quick intuition of the QW evolution for different initial probability distributions. This allows one, in particular, to reach a desired asymptotic distribution by suitably tailoring the initial (extended) state [13]. Of course, this qualitative knowledge will be only approximate (as the continuous solutions are strictly valid only for infinitely extended initial distributions) but, as we show below, even for relatively narrow initial distributions the qualitative knowledge turns out to be quite accurate. In each case, we will present an exact numerical simulation, obtained from the discrete map of the QW, that illustrates the agreement with the qualitative analysis described above. As far as we know, this is the first time that continuous approximations to the discrete QW are derived for a multidimensional QW.

The continuous equations we derive, however, cannot apply near degeneracies in the dispersion relations, as in their vicinity the eigensolutions of the QW vary wildly from point to point. Hence, close to degeneracies a specific treatment is necessary and we develop such a treatment for a particular type of degeneracy, namely conical intersections. We apply this treatment to the 2D Grover walk that exhibits such conical intersections in the dispersion relation, a fact that seems to have been noticed only recently [40] and that we study here in mathematical detail. These conical intersections, called diabolical points [59–61], determine a type of dynamics similar to that of massless Dirac fermions or of electrons in graphene [62–66] and appear in very different systems: quantum triangular billiards [59], conical refraction in crystal optics [60, 61], the already mentioned graphene electrons [62–66], the spectra of polyatomic molecules [67, 68], optical lattices [69] or acoustic surface waves [70], just to mention a few. In all cases, the remarkable physical properties exhibited by these systems are directly associated with the existence of diabolical points. In the 3D Grover walk, we also find degeneracies, but these are different from the diabolical points. Diabolical points on the dispersion relations also play a role in the analysis of topological properties of the QW. Configurations with different topologies or Chern numbers can be manufactured, for example, by introducing a spatial variation of one of the coin angles in alternate QWs, so that these regions possess different topologies. At the boundaries delimiting those regions the gap in the dispersion relation closes, and one can observe the formation of bound states or unidirectionally propagating modes, which are protected by topological invariants and by the existence of the gap away from the boundary [71].

The rest of the paper is organized as follows. In section 2, we introduce the N-dimensional QW and solve it formally. In section 3, we formally derive continuous equations in the general N-dimensional case and discuss how one can treat conical intersections. In sections 4 and 5, we apply our results to the 2D and 3D Grover QWs, respectively. Finally, in section 6 we present our main conclusions.

2. N-dimensional discrete quantum walks: generalities

In the N-dimensional QW, the walker moves at discrete time steps $t\in \mathbb {N}$ across an N-dimensional lattice of sites $\mathbf {x}\equiv \left (x_{1},\ldots ,x_{N}\right )\in \mathbb {Z}^{N}$ . The walker is endowed with a 2N-dimensional coin that, after convenient unitary transformation, determines the direction of displacement. The Hilbert space of the whole system (walker + coin) then has the form

$\begin{equation} \mathcal{H}=\mathcal{H}_{\mathrm{P}}\otimes\mathcal{H}_{\mathrm{C}}, \end{equation} \tag{ 2 }$

where the position space, $\mathcal {H}_{\mathrm {P}}$ , is spanned by $\left \{ \left \vert \mathbf {x}\right \rangle \equiv \left \vert x_{1},\ldots ,x_{N}\right \rangle :x_{\alpha }\in \mathbb {Z};\,\,\alpha =1,\ldots ,N\right \}$ ( $\left \langle \mathbf {x}\right .\left \vert \mathbf {x}^{\prime }\right \rangle =\delta _{\mathbf {x,x}^{\prime }}$ ), and the coin space, $\mathcal {H}_{\mathrm {C}}$ , is spanned by 2N orthonormal quantum states $\left \{ \left \vert \alpha _{\eta }\right \rangle :\alpha =1,\ldots ,N;\eta =\pm \right \}$ . Note that α is associated with the axis and η with the direction. For example, in the popular one-dimensional QW (N = 1) we would have just $\left \vert 1_{+}\right \rangle$ and $\left \vert 1_{-}\right \rangle$ , which are equivalent to the $\left \vert R\right \rangle$ and $\left \vert L\right \rangle$ (for right and left) states commonly used in the literature. This notation is introduced in order to easily deal with an arbitrary number of dimensions.

The state of the total system at time t is represented by the ket $\left \vert \psi _{t}\right \rangle$ , which can be expressed in the form

$\begin{equation} \vert \psi_{t}\rangle =\sum_{\mathbf{x}}\sum_{\alpha=1}^{N}\sum_{\eta=\pm}\psi_{\mathbf{x},t}^{\alpha,\eta} \vert \mathbf{x}\rangle \otimes\vert \alpha_{\eta}\rangle , \end{equation} \tag{ 3 }$

where the projections

$\begin{equation} \psi_{\mathbf{x},t}^{\alpha,\eta}=(\langle \alpha_{\eta}\vert \otimes\langle \mathbf{x}\vert )\vert \Psi_{t} \end{equation} \tag{ 4 }$

are wavefunctions on the lattice. We find it convenient to define, at each point x, the following ket:

$\begin{equation} \vert \psi_{\mathbf{x},t}\rangle =\left\langle \mathbf{x}\right.\vert \psi_{t}\rangle =\sum_{\alpha=1}^{N}\sum_{\eta=\pm}\psi_{\mathbf{x},t}^{\alpha,\eta}\vert \alpha_{\eta}\rangle, \end{equation} \tag{ 5 }$

which is an (unnormalized) coin state, so that ψ^α,η_x,t = 〈α_η|ψ_x,t〉. As $\vert \psi _{\mathbf {x},t}^{\alpha ,\eta }\vert ^{2}=\vert(\langle \alpha _{\eta }\vert \otimes\langle \mathbf {x}\vert)\vert \psi _{t} \rangle \vert ^{2}$ is the probability of finding the walker at $\left (\mathbf {x},t\right )$ and the coin in state $\left \vert \alpha _{\eta }\right \rangle$ , the probability of finding the walker at $\left (\mathbf {x},t\right )$ irrespective of the coin state is then

$\begin{equation} P_{\mathbf{x},t}=\sum_{\alpha=1}^{N}\sum_{\eta=\pm}\vert \psi_{\mathbf{x},t}^{\alpha,\eta}\vert ^{2}=\langle \psi_{\mathbf{x},t}\vert \psi_{\mathbf{x},t}\rangle , \end{equation} \tag{ 6 }$

where we used the fact that $\sum _{\alpha =1}^{N}\sum _{\eta =\pm } \vert \alpha _{\eta } \rangle \langle \alpha _{\eta }\vert$ is the identity in $\mathcal {H}_{\mathrm {C}}$ . Clearly, $\sum _{\mathbf {x}}P_{\mathbf {x},t}=1$ because $\sum _{\mathbf {x}}\left \vert \mathbf {x}\right \rangle \left \langle \mathbf {x}\right \vert$ is the identity in $\mathcal {H}_{\mathrm {P}}$ .

The dynamical evolution of the system is ruled by

$\begin{equation} \vert \psi_{t+1}\rangle ={\hat{U}}\vert \psi_{t}\rangle , \end{equation} \tag{ 7 }$

where the unitary operator

$\begin{equation} \hat{U}=\hat{D}\circ(\hat{I}\otimes\hat{C}) \end{equation} \tag{ 8 }$

is given in terms of the identity operator in $\mathcal {H}_{\mathrm {P}}$ , $\hat {I}$ , and two more unitary operators. On the one hand, $\hat {C}$ is the so-called coin operator (an operator in $\mathcal {H}_{\mathrm {C}}$ ), which can be written in its more general form as

$\begin{equation} \hat{C}=\sum_{\alpha,\alpha^{\prime}=1}^{N}\,\,\sum_{\eta,\eta^{\prime}=\pm}C_{\alpha^{\prime},\eta^{\prime}}^{\alpha,\eta}\vert \alpha_{\eta}\rangle \langle \alpha_{\eta^{\prime}}^{\prime}\vert , \end{equation} \tag{ 9 }$

where the matrix elements $C_{\alpha ^{\prime },\eta ^{\prime }}^{\alpha ,\eta }\equiv \langle \alpha _{\eta } \vert \hat {C} \vert \alpha _{\eta ^{\prime }}^{\prime }\rangle$ can be arranged as a 2N × 2N unitary square matrix C. On the other hand, $\hat {D}$ is the conditional displacement operator in $\mathcal {H}$ ,

$\begin{equation} \hat{D}=\sum_{\mathbf{x}}\sum_{\alpha=1}^{N}\sum_{\eta=\pm}\left\vert \mathbf{x}+\eta\mathbf{u}_{\alpha}\right\rangle \left\langle \mathbf{x}\right\vert \otimes\vert \alpha_{\eta}\rangle \langle \alpha_{\eta}\vert , \end{equation} \tag{ 10 }$

where u_α is the unit vector along direction x_α; note that, depending on the coin state $\vert \alpha _{\eta }\rangle$ , the walker moves one site to the positive or negative direction of x_α if η = + or −, respectively.

Projecting (7) onto $\left \langle \mathbf {x}\right \vert$ and using (4) and (8)–(10), we obtain straightforwardly

$\begin{equation} \vert \psi_{\mathbf{x},t+1}\rangle =\sum_{\alpha=1}^{N}\sum_{\eta=\pm}\vert \alpha_{\eta}\rangle \langle \alpha_{\eta}\vert \hat{C}\vert \psi_{\mathbf{x}-\eta\mathbf{u}_{\alpha},t}\rangle , \end{equation} \tag{ 11 }$

which further projected onto $\langle \alpha _{\eta }\vert$ leads to

$\begin{equation} \psi_{\mathbf{x},t+1}^{\alpha,\eta}=\sum_{\alpha^{\prime}=1}^{N}\sum_{\eta^{\prime}=\pm}C_{\alpha^{\prime},\eta^{\prime}}^{\alpha,\eta}\psi_{\mathbf{x}-\eta\mathbf{u}_{\alpha},t}^{\alpha^{\prime},\eta^{\prime}}.\end{equation} \tag{ 12 }$

Equation (11), or equivalently (12), is the N-dimensional QW map in position representation; it shows that, at each (discrete) time, the wavefunctions at each point are coherent linear superpositions of wavefunctions at neighboring points at a previous time, the weights of the superposition being given by the coin operator matrix elements C^α,η_α',η'. Next we proceed to derive the solution of map (12).

Given the linearity of the map and the fact that it is space invariant (C^α,η_α',η' do not depend on space), a useful technique here is the spatial discrete Fourier transform (DFT), which has been used many times in QW studies (see e.g. [72]). First we define the DFT pair

$\begin{equation} \vert \tilde{\psi}_{\mathbf{k,}t}\rangle \equiv\sum_{\mathbf{x}}\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\vert \psi_{\mathbf{x},t}\rangle , \end{equation} \tag{ 13 }$

$\begin{equation} \vert \psi_{\mathbf{x},t}\rangle \equiv\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{x}}\vert \tilde{\psi}_{\mathbf{k},t}\rangle , \end{equation} \tag{ 14 }$

where $\mathbf {k}=\left (k_{1},\ldots ,k_{N}\right )$ and $k_{\alpha }\in \left [-\pi ,\pi \right ]$ is the (quasi-)momentum vector [73]⁵. Applying the previous definitions to the map (11), we readily obtain

$\begin{equation} \vert \tilde{\psi}_{\mathbf{k},t+1}\rangle =\hat{C}_{\mathbf{k}}\vert \tilde{\psi}_{\mathbf{k},t}\rangle , \end{equation} \tag{ 15 }$

where we defined a coin operator in quasi-momentum space

$\begin{equation} \hat{C}_{\mathbf{k}}\equiv\sum_{\alpha=1}^{N}\sum_{\eta=\pm}\left\vert \alpha_{\eta}\right\rangle \left\langle \alpha_{\eta}\right\vert \hat{C}\mathrm{e}^{-\mathrm{i}\eta k_{\alpha}}, \end{equation} \tag{ 16 }$

k_α = k·u_α, whose matrix elements read

$\begin{equation} \langle \alpha_{\eta}\vert \hat{C}_{\mathbf{k}}\vert \alpha_{\eta^{\prime}}^{\prime}\rangle \equiv\left(C_{\mathbf{k}}\right)_{\alpha^{\prime},\eta^{\prime}}^{\alpha,\eta}=\mathrm{e}^{-\mathrm{i}\eta k_{\alpha}}C_{\alpha^{\prime},\eta^{\prime}}^{\alpha,\eta}. \end{equation} \tag{ 17 }$

Projection of (15) onto $\langle \alpha _{\eta } \vert$ and use of (16) and (17) leads to

$\begin{equation} \tilde{\psi}_{\mathbf{k},t+1}^{\alpha,\eta}=\sum_{\alpha^{\prime}=1}^{N}\sum_{\eta^{\prime}=\pm}\mathrm{e}^{-\mathrm{i}\eta\mathbf{k}\cdot\mathbf{u}_{\alpha}}C_{\alpha^{\prime},\eta^{\prime}}^{\alpha,\eta}\tilde{\psi}_{\mathbf{k},t}^{\alpha^{\prime},\eta^{\prime}}. \end{equation} \tag{ 18 }$

Hence the non-local maps (11) and (12) become local in the momentum representation (15) and (18). This allows solving formally the QW dynamics very easily because map (15) implies

$\begin{equation} \vert \tilde{\psi}_{\mathbf{k},t}\rangle =(\hat{C}_{\mathbf{k}})^{t}\vert \tilde{\psi}_{\mathbf{k},0}\rangle, \end{equation} \tag{ 19 }$

and hence the eigensystem of $\hat {C}_{\mathbf {k}}$ (or of C_k in matrix form) is most useful for solving the problem, as we do next.

As the operator $\hat {C}_{\mathbf {k}}\mathcal {\ }$ is unitary, its eigenvalues all have the form $\lambda _{\mathbf {k}}^{\left (s\right )}=\exp (-\mathrm {i}\omega _{\mathbf {k}}^{\left (s\right )})$ , s = 1,...,N, with $\omega _{\mathbf {k}}^{\left (s\right )}$ real. We will need to know the $\hat {C}_{\mathbf {k}}$ eigenstates too, $\{\vert \phi _{\mathbf {k}}^{\left (s\right )} \rangle\} _{s=1}^{2N}$ . Once the eigensystem of $\hat {C}_{\mathbf {k}}$ is known, implementing (19) is trivial. Given an initial distribution of the walker in position representation $\vert \psi _{\mathbf {x},0}\rangle$ , we compute its DFT $\vert \tilde {\psi }_{\mathbf {k},0}\rangle$ via (13), as well as the projections

$\begin{equation} \skew3\tilde{f}_{\mathbf{k}}^{\left(s\right)}=\langle \phi_{\mathbf{k}}^{\left(s\right)}\vert \tilde{\psi}_{\mathbf{k},0}\rangle , \end{equation} \tag{ 20 }$

so that $\vert \tilde {\psi }_{\mathbf {k},0}\rangle =$ $\sum _{s}\tilde {f}_{\mathbf {k}}^{\left (s\right )}\vert \phi _{\mathbf {k}}^{\left (s\right )}\rangle$ . Now recalling (19), we arrive at

$\begin{equation} \vert \tilde{\psi}_{\mathbf{k},t}\rangle =\sum_{s=1}^{2N}\mathrm{e}^{-\mathrm{i}\omega_{\mathbf{k}}^{\left(s\right)}t}\skew3\tilde{f}_{\mathbf{k}}^{\left(s\right)}\vert \phi_{\mathbf{k}}^{\left(s\right)}\rangle , \end{equation} \tag{ 21 }$

where we used $\lambda _{\mathbf {k}}^{\left (s\right )}=\exp(-\mathrm {i}\omega _{\mathbf {k}}^{\left (s\right )})$ , while in position representation we obtain, using (14),

$\begin{equation} \vert \psi_{\mathbf{x},t}\rangle =\sum_{s=1}^{2N}\vert \psi_{\mathbf{x},t}^{\left(s\right)}\rangle , \end{equation} \tag{ 22 }$

$\begin{equation} \vert \psi_{\mathbf{x},t}^{\left(s\right)}\rangle =\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\mathrm{e}^{\mathrm{i}\left(\mathbf{k}\cdot\mathbf{x-}\omega_{\mathbf{k}}^{\left(s\right)}t\right)}\skew3\tilde{f}_{\mathbf{k}}^{\left(s\right)}\vert \phi_{\mathbf{k}}^{\left(s\right)}\rangle . \end{equation} \tag{ 23 }$

Hence the QW is formally solved: all we need is to compute the eigensystem of $\hat {C}_{\mathbf {k}}$ and the initial state in reciprocal space $\vert \tilde {\psi }_{\mathbf {k},0}\rangle$ , which determines the weight functions $\skew3\tilde {f}_{\mathbf {k}}^{\left (s\right )}$ through (20).

Equation (22) shows that the QW dynamics corresponds to the superposition of 2N independent walks, labeled by s. According to (23) the $\omega _{\mathbf {k}}^{\left (s\right )}$ 's are the frequencies of the map, each of which defines a dispersion relation in the system (2N in total). Note as well that what we have done in the end is to decompose the QW dynamics in terms of plane waves. In particular, if $\skew3\tilde {f}_{\mathbf {k}}^{\left (s\right )}=\delta ^{(N)}\left (\mathbf {k}-\mathbf {k}_{0}\right )$ , this means that $\vert \tilde {\psi }_{\mathbf {k},0}\rangle$ is different from zero only for k = k₀, $\vert \psi _{\mathbf {x},t}^{\left (s\right )}\rangle =\left (2\pi \right )^{-N}\exp [\mathrm {i}(\mathbf {k}_{0}\cdot \mathbf {x-}\omega _{\mathbf {k}_{0}}^{\left (s\right )}t)]\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )} \rangle$ , which is an unnormalizable plane wave and thus unphysical.

In order to avoid possible confusion, to conclude this initial part we state that the ordering of the coin base elements we will be using in the matrix representations of operators and kets is $\left \vert 1_{+}\right \rangle ,\left \vert 1_{-}\right \rangle ,\ldots ,\left \vert N_{+}\right \rangle ,\left \vert N_{-}\right \rangle$ .

3. Continuous wave equations for spatially extended initial conditions

In this section, we describe the evolution of spatially extended initial conditions that are close to the plane waves we have introduced. We define such spatially extended states as those having a width that is appreciably larger than the lattice spacing (taken as unity in this work). These types of initial states are wavepackets that are easily expressed in reciprocal space;

$\begin{equation} \vert \tilde{\psi}_{\mathbf{k},0}\rangle =\sum_{s=1}^{2N}\tilde{F}_{\mathbf{k}-\mathbf{k}_{0}}^{\left(s\right)}\vert \phi_{\mathbf{k}_{0}}^{\left(s\right)}\rangle , \end{equation} \tag{ 24 }$

where k₀ is a reference (carrier) wavevector, chosen at will, $\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )} \rangle$ are the associated eigenvectors and $\tilde {F}_{\mathbf {k}}^{\left (s\right )}$ is a narrow function of k, centered at k = 0 ( $\tilde {F}_{\mathbf {k}-\mathbf {k}_{0}}^{\left (s\right )}$ is centered at k = k₀) and having a very small width $\Delta k^{\left (s\right )}{\ll }\pi.$ ⁶ Note that we have chosen an initial coin state that is independent of k. From here one must distinguish between regular points, where eigenvectors $\vert \phi _{\mathbf {k}}^{\left (s\right )} \rangle$ have a smooth dependence on k close to k₀, and degeneracy points, where eigenvectors have wild variations around them, as we will see later.

According to (14) the initial condition (24) reads in position representation

$\begin{equation} \vert \psi_{\mathbf{x},0}\rangle =\mathrm{e}^{\mathrm{i}\mathbf{k}_{0}\cdot\mathbf{x}}\sum_{s=1}^{2N}F_{\mathbf{x},0}^{\left(s\right)}\vert \phi_{\mathbf{k}_{0}}^{\left(s\right)}\rangle , \end{equation} \tag{ 25 }$

where

$\begin{equation} F_{\mathbf{x},0}^{\left(s\right)}=\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\mathrm{e}^{\mathrm{i}\left(\mathbf{k}-\mathbf{k}_{0}\right)\cdot\mathbf{x}}\tilde{F}_{\mathbf{k}-\mathbf{k}_{0}}^{\left(s\right)} \end{equation} \tag{ 26 }$

is a wide and smooth function of x because $\tilde {F}_{\mathbf {k}}^{\left (s\right )}$ is concentrated around k = 0. Hence in real space our initial condition consists of a coin state $\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )}\rangle$ equal at all points, multiplied by the carrier $\exp \left [\mathrm {i}\left (\mathbf {k}_{0}\cdot \mathbf {x}\right )\right ]$ and by a wide and smooth function of space, $F_{\mathbf {x},0}^{\left (s\right )}$ . We see then that the types of initial conditions we are dealing with are very close to the plane waves of momentum k₀ of the QW, given by $\mathrm {e}^{\mathrm {i}\mathbf {k}_{0}\cdot \mathbf {x}}\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )} \rangle$ .

3.1. Regular points

For regular points (i.e. far from degeneracies) the initial condition (24) determines the weight functions (20) as

$\begin{equation} \tilde{f}_{\mathbf{k}}^{\left(s\right)}=\sum_{s^{\prime}}\tilde{F}_{\mathbf{k}-\mathbf{k}_{0}}^{\left(s^{\prime}\right)}\langle\phi_{\mathbf{k}}^{\left(s\right)}\mid\phi_{\mathbf{k}_{0}}^{\left(s^{\prime}\right)}\rangle=\tilde{F}_{\mathbf{k}-\mathbf{k}_{0}}^{\left(s\right)}+O(\Delta k^{\left(s\right)}), \end{equation} \tag{ 27 }$

where we take into account that the eigenvectors $\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )} \rangle$ vary smoothly around k = k₀ and that only for k ≈ k₀ is the function $\tilde{F}_{\mathbf {k}-\mathbf {k}_{0}}^{(s)}$ non-vanishing. Note that we cannot make such an approximation in the case of degeneracy points because of the strong variations of the eigenvectors around them.

Now, the system will evolve according to (23). Approximating the eigenvector $\vert \phi _{\mathbf {k}}^{\left (s\right )}\rangle$ appearing in (23) by $\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )}\rangle$ using the same arguments as before, the partial waves $\vert \psi _{\mathbf {x},t}^{\left (s\right )}\rangle$ can be written as

$\begin{equation} \fl \vert \psi_{\mathbf{x},t}^{\left(s\right)}t\rangle =F_{\mathbf{x},t}^{\left(s\right)}\exp[\mathrm{i}(\mathbf{k}_{0}\cdot\mathbf{x-}\omega_{\mathbf{k}_{0}}^{\left(s\right)}t)]\vert \phi_{\mathbf{k}_{0}}^{\left(s\right)}\rangle +O(\Delta k^{\left(s\right)}), \end{equation} \tag{ 28 }$

$\begin{equation} \fl F_{\mathbf{x},t}^{\left(s\right)} =\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\exp\left(\mathrm{i}\left[\left(\mathbf{k}-\mathbf{k}_{0}\right)\cdot\mathbf{x-}\left(\omega_{\mathbf{k}}^{\left(s\right)}-\omega_{\mathbf{k}_{0}}^{\left(s\right)}\right)t\right]\right) \tilde{F}_{\mathbf{k}-\mathbf{k}_{0}}^{\left(s\right)}, \end{equation} \tag{ 29 }$

where we have defined new functions $\{ F_{\mathbf {x},t}^{\left (s\right )} \} _{s=1}^{2N}$ . Note that at t = 0, $\vert \psi _{\mathbf {x},0}^{\left (s\right )}\rangle =F_{\mathbf {x},0}^{\left (s\right )}\exp \left [\mathrm {i}\left (\mathbf {k}_{0}\cdot \mathbf {x}\right )\right ]\vert \phi _{\mathbf {k}_{0}}^{\left (s\right )}\rangle$ , in agreement with (25).

The problem is then solved if functions $\{ F_{\mathbf {x},t}^{\left (s\right )}\} _{s=1}^{2N}$ are determined. Instead of doing so by brute force, i.e. by integrating—maybe numerically—equation (29), what we do now is to look for a continuous wave equation that, by comparison with other known cases, sheds light on the expected behavior of the QW. In order to test the quality of our analysis based on our continuous equations and the knowledge of the dispersion relations, we will present some plots that are obtained directly by iteration of equation (12) without any approximation.

To avoid confusion, we introduce a function $F^{(s)}({\bf x},t)$ of continuous real arguments in exactly the same way as in (29), as there is nothing forbidding such a definition. We then have $F_{\mathbf {x},t}^{\left (s\right )}=F^{(s)}({\bf x},t)$ for $\mathbf {x\in \mathbb {Z}}^{N}$ and $t\in \mathbb {N}$ . In order to simplify the derivation, we rewrite equation (29) by making the variable change k − k₀ → k. Finally, as for the limits of integration (originally from −π to +π for each dimension), we extend them from −∞ to +∞ in agreement with the continuous limit. We then have

$\begin{equation} F^{\left(s\right)}({\bf x},t) \equiv\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\exp\left[\mathrm{i}\left(\mathbf{k}\cdot\mathbf{x-}\left[\omega_{\mathbf{k}_{0}+\mathbf{k}}^{\left(s\right)}-\omega_{\mathbf{k}_{0}}^{\left(s\right)}\right]t\right)\right] \tilde{F}_{\mathbf{k}}^{\left(s\right)}. \end{equation} \tag{ 30 }$

Let us obtain the (approximate) wave equation. First we take the time derivative of (30) and obtain

$\begin{equation} \fl \mathrm{i}\partial_{t}F^{\left(s\right)}({\bf x},t) =\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\exp\left[\mathrm{i}\left(\mathbf{k}\cdot\mathbf{x-}\left[\omega_{\mathbf{k}_{0}+\mathbf{k}}^{\left(s\right)}-\omega_{\mathbf{k}_{0}}^{\left(s\right)}\right]t\right)\right] \left[\omega_{\mathbf{k}_{0}+\mathbf{k}}^{\left(s\right)}-\omega_{\mathbf{k}_{0}}^{\left(s\right)}\right]\tilde{F}_{\mathbf{k}}^{\left(s\right)}. \end{equation} \tag{ 31 }$

Then we Taylor expand the function $[\omega _{\mathbf {k}_{0}+\mathbf {k}}^{\left (s\right )}-\omega _{\mathbf {k}_{0}}^{\left (s\right )}]$ (except in the exponent: otherwise large errors at long times would be introduced) around k = 0,

$\begin{equation} \omega_{\mathbf{k}_{0}+\mathbf{k}}^{\left(s\right)}-\omega_{\mathbf{k}_{0}}^{\left(s\right)}=\sum_{i=1}^{N}\varpi_{i}^{\left(s\right)}k_{i}+\frac{1}{2!}\sum_{i,j=1}^{N}\varpi_{ij}^{\left(s\right)}k_{i}k_{j}+\cdots, \end{equation} \tag{ 32 }$

where

$\begin{eqnarray} \varpi_{i}^{(s)}&=&\!\left.\partial\omega_{\mathbf{k}_{0}+\mathbf{k}}^{(s)}/\partial k_{i}\right\vert _{\mathbf{k}=0}\nonumber\\ \\\nonumber &=&\!\left.\partial\omega_{\mathbf{k}}^{(s)}/\partial k_{i}\right\vert _{\mathbf{k}=\mathbf{k}_{0}}=\left[\mathbf{v}_{\mathrm{g}}^{(s)}(\mathbf{k}_{0})\right]_{i},\nonumber \end{eqnarray} \tag{ 33 }$

$\begin{eqnarray}\varpi_{ij}^{(s)}&=&\!\left.\partial^{2}\omega_{\mathbf{k}_{0}+\mathbf{k}}^{(s)}/\partial k_{i}\partial k_{j}\right\vert _{\mathbf{k}=0}\nonumber\\ \\\nonumber &=&\!\left.\partial^{2}\omega_{\mathbf{k}}^{(s)}/\partial k_{i}\partial k_{j}\right\vert _{\mathbf{k}=\mathbf{k}_{0}},\nonumber \end{eqnarray} \tag{ 34 }$

etc. In the latter equation, $\mathbf {v}_{\mathrm {g}}^{(s)}(\mathbf {k}_{0})=\left .\nabla _{\mathbf {k}}\omega _{\mathbf {k}}^{(s)}\right \vert _{\mathbf {k}=\mathbf {k}_{0}}$ is the group velocity at the point k = k₀ (see the discussion below). Now it is trivial to transform the so-obtained right-hand side of equation (31) as a sum of spatial derivatives of $F^{(s)}({\bf x},t)$ , leading to the result

$\begin{equation} \fl \mathrm{i}\frac{\partial F^{(s)}({\bf x},t)}{\partial t} =-\mathrm{i}\mathbf{v}_{\mathrm{g}}^{(s)}({\bf k}_0)\cdot\nabla F^{(s)} -\frac{1}{2!}\sum_{i,j=1}^{N}\varpi_{ij}^{(s)}\frac{\partial^{2}F^{(s)}}{\partial x_{i}\partial x_{j}}+\cdots, \end{equation} \tag{ 35 }$

which is the sought continuous wave equation. One should understand that, in general, few terms are necessary on the right-hand side of equation (35) because $F^{(s)}({\bf x},t)$ is a slowly varying function of space, as commented on.

The first term on the right-hand side of the wave equation is an advection term, implying that the initial condition $F_{0}^{(s)}({\bf x})\equiv F^{(s)}({\bf x},0)$ is shifted in time as $F^{(s)}({\bf x},t)=F_{0}^{(s)}(\mathbf {x}-\mathbf {v}_{\mathrm {g}}^{(s)}(\mathbf {k}_{0})t)$ without distortion, to leading order. In fact, we can get rid of this term by defining a moving reference frame X such that

$\begin{equation} \mathbf{X}=\mathbf{x}-\mathbf{v}_{\mathrm{g}}^{(s)}({\bf k}_0)t,\ \ \ A^{(s)}({\bf X},t)=F^{(s)}({\bf x},t) \end{equation} \tag{ 36 }$

and then equation (35) becomes

$\begin{equation} \mathrm{i}\frac{\partial A^{(s)}({\bf X},t)}{\partial t}=-\frac{1}{2}\sum_{i,j=1}^{N}\varpi_{ij}^{(s)}\frac{\partial^{2}A^{(s)}}{\partial X_{i}\partial X_{j}}+\cdots, \end{equation}{} \tag{ 37 }$

which is an N-dimensional Schrödinger-like equation, to leading order. This means that the evolution of the wavepacket consists of the advection of the initial wavepacket at the corresponding group velocity and, on top of that, the wavefunction itself evolves according to (37). Equation (37) is a main result of this paper as it governs the non-trivial dynamics of the QW when spatially extended initial conditions, as we have defined them, are considered. It evidences the role played by the dispersion relations as anticipated: for distributions whose DFT is centered around some k₀, the local variations of ω around k₀ determine the type of wave equation controlling the QW dynamics.

A few additional remarks. Firstly, if the initial condition only projects onto one of the sheets of the dispersion relation (hence the sum s in (24) reduces to a single element), all $F^{(s)}({\bf X},t)$ will be zero, except the one corresponding to the chosen eigenvector. This means that the coin state is preserved along the evolution, to leading order. Secondly, if the initial state projects onto several eigenvectors the probability of finding the walker at $({\bf x},t)$ , irrespectively of the coin state, follows from (6), (22) and (28), and reads

$\begin{equation} P_{\mathbf{x},t}=\sum_{s}\left\vert F_{\mathbf{x},t}^{(s)}\right\vert ^{2}=\sum_{s}\left\vert A^{(s)}({\bf X},t)\right\vert ^{2} \end{equation} \tag{ 38 }$

to leading order, where we used $\langle \phi _{\mathbf {k}_{0}}^{(s)}\mid \phi _{\mathbf {k}_{0}}^{\left (s^{\prime }\right )}\rangle =\delta _{s,s^{\prime }}$ . Hence P_x,t is just the sum of partial probabilities: there is no interference among the different sub-QWs.

Finally, one might consider an initial condition consisting in a linear combination of initial terms of the form (25), each having a different k₀. It is easy to see that we can generalize the above treatment to this case, each term evolving in accordance with the corresponding k₀. Of course, this situation may show more complicated behavior arising from interference effects among the different terms. We have not considered this richer evolution here, and we have restricted ourselves to a simpler study, where the evolution can be traced back to the analysis of the dispersion relation around a single k₀ point.

3.2. Degeneracy points

This case requires special care because there are eigenvectors of the QW that vary strongly around the degeneracy, forbidding the very initial approximation taken in the case of regular points, namely equation (27). Given the singular nature of the problem we will try to give a rather general (but not fully general) theory here. We will present a treatment that covers the case of (hyper-) conical intersections. These appear in the 2D Grover walk (as we will see in the next section), in the 3D alternate QW [40] (that will not be treated here) and, very likely, in some other cases. The basic assumptions are: (i) there are just two (hyper-) sheets in the dispersion relation that display a conical intersection at k = k_D (diabolical point) and we label them as s = 1,2 for the sake of definiteness; and (ii) there are other sheets, degenerate with the diabolical ones at k = k_D, whose associated frequencies $\omega _{\mathbf {k}}^{(s)}$ are constant around the diabolical point.

In order to alleviate the notation, we will denote by ω_{_D} the value of the degenerate frequencies at k = k_D:

$\begin{equation} \omega_{_{\mathrm{D}}}=\omega_{\mathbf{k}_{\mathrm{D}}}^{\left(1\right)}=\omega_{\mathbf{k}_{\mathrm{D}}}^{\left(2\right)}. \end{equation} \tag{ 39 }$

Note that there can be other sheets with $\omega _{\mathbf {k}_{\mathrm {D}}}^{\left (s\neq 1,2\right )}=\omega _{_{\mathrm {D}}}$ (as it happens with $\omega _{\mathbf {k}_{\mathrm {D}}}^{\left (3\right )}$ in the 2D Grover map).

We consider an initial wavepacket defined in reciprocal space by

$\begin{equation} \vert \tilde{\psi}_{\mathbf{k},0}\rangle =\tilde{F}_{\mathbf{k}-\mathbf{k}_{\mathrm{D}}}\left\vert \Xi\right\rangle , \end{equation} \tag{ 40 }$

where $\left \vert \Xi \right \rangle$ is an eigenstate of the coin operator $\hat {C}_{\mathbf {k}}$ at k = k_D, and $\tilde {F}_{\mathbf {k}-\mathbf {k}_{\mathrm {D}}}$ is again a narrow function centered at k = k_D. In position representation this initial condition reads, see (14),

$\begin{equation} |\psi_{\mathbf{x},0}\rangle=\mathrm{e}^{\mathrm{i}\mathbf{k}_{\mathrm{D}}\cdot\mathbf{x}}F_{\mathbf{x},0}\left\vert \Xi\right\rangle ,\end{equation} \tag{ 41 }$

where $F_{\mathbf {x},0}=\int \frac {\mathrm {d}^{N}\mathbf {k}}{\left (2\pi \right )^{N}}\mathrm {e}^{\mathrm {i}\left (\mathbf {k-k}_{\mathrm {D}}\right )\cdot \mathbf {x}}\tilde {F}_{\mathbf {k}-\mathbf {k}_{\mathrm {D}}}$ is the DFT of $\tilde {F}_{\mathbf {k}}$ . Hence, as in the regular case, $\left \vert \Xi \right \rangle$ represents the state of the coin at the initial condition. Clearly $\left \vert \Xi \right \rangle$ is not unique because of the degeneracy. We do not use the notation $\vert \phi _{\mathbf {k}_{\mathrm {D}}}^{(s)}\rangle$ but instead $\left \vert \Xi \right \rangle$ because the former are ill defined, as we have seen in the example of the 2D Grover walk.

The initial condition (40) determines the weight functions (20) as

$\begin{equation} \skew3\tilde{f}_{\mathbf{k}}^{(s)}=\skew3\tilde{F}_{\mathbf{k}-\mathbf{k}_{\mathrm{D}}}\langle\phi_{\mathbf{k}}^{(s)}\mid\Xi\rangle. \end{equation} \tag{ 42 }$

We will assume that $\left \vert \Xi \right \rangle$ does not project on regular sheets (those not degenerate with the conical intersection); otherwise those projections will evolve as described in the previous case of regular points. This is equivalent to saying that, in the remainder of this subsection, all sums on s will be restricted to sheets that are degenerate at the conical intersection.

Definition (42) poses no problem except for k = k_D. However, as this is a single point (a set of null measure) its influence on the final result, given by integrals in k, is null⁷.

A main difference with the regular case is seen at this stage because there is no choice of the initial coin state $\left \vert \Xi \right \rangle$ ensuring that only one value of s is populated (recall that $\vert \phi _{\mathbf {k}}^{(s)} \rangle$ vary strongly around k = k_D). This has consequences, as we see next.

Now the system will evolve according to (23). The partial waves $\vert \psi _{\mathbf {x},t}^{(s)} \rangle$ now read

$\begin{equation} \vert \psi_{\mathbf{x},t}^{(s)}\rangle =\mathrm{e}^{\mathrm{i}\left(\mathbf{k}_{\mathrm{D}}\cdot\mathbf{x-}\omega_{_{\mathrm{D}}}t\right)}\left\vert \mathbf{F}_{\mathbf{x},t}^{(s)}\right\rangle , \end{equation} \tag{ 43 }$

$\begin{equation} \fl \left\vert \mathbf{F}_{\mathbf{x},t}^{(s)}\right\rangle =\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\exp[\mathrm{i}(\mathbf{k}\cdot\mathbf{x-}[\omega_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}-\omega_{_{\mathrm{D}}}]t)] \tilde{F}_{\mathbf{k}}\,\,\langle\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\mid\Xi\rangle\mid\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\rangle, \end{equation} \tag{ 44 }$

where we made the variable change k − k_D → k and extended the limits of integration to infinity as corresponding to the continuous limit. Note that $\omega _{\mathbf {k}_{\mathrm {D}}}^{\left (s\neq 1,2\right )}-\omega _{_{\mathrm {D}}}=0$ according to assumption (ii) above. As in the regular case, the evolution has a fast part, given by the carrier wave $\mathrm {e}^{\mathrm {i}\left (\mathbf {k}_{\mathrm {D}}\cdot \mathbf {x-}\omega _{_{\mathrm {D}}}t\right )}$ , and a slow part (both in time and in space) given by $\vert \mathbf {F}_{\mathbf {x},t}^{(s)}\rangle$ .

Note that the previous approximations reproduce the correct result at t = 0, equation (41), upon using that $\sum _{s}\vert \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{(s)}\rangle \langle \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{(s)}\vert$ is very approximately the unit operator in the degenerate subspace, to which $\left \vert \Xi \right \rangle$ belongs.

It is important to understand that the vector (coin) part of the state evolves in time in this diabolical point situation, because more than one of the sheets that become degenerate at the conical intersection become necessarily populated, i.e. there are at least two s values for which $\langle \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{(s)}\mid \Xi \rangle \neq 0$ , unlike the regular case.

Finding a wave equation in this case is by far more complicated than in the regular case, as the vector part of $\vert \mathbf {F}_{\mathbf {x},t}^{(s)}\rangle$ depends on time (because $\vert \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{(s)} \rangle$ cannot be approximated by $\vert \phi _{\mathbf {k}_{\mathrm {D}}}^{(s)}\rangle$ ), unlike the regular case where it can be approximated by the constant vector $\vert \phi _{\mathbf {k}_{0}}^{(s)}\rangle$ ; see (28). We will not try deriving such a wave equation but confine ourselves to trying to understand the evolution of the system under the assumed conditions.

We can say that the frequency offsets $[\omega _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{\left (s=1,2\right )}-\omega _{_{\mathrm {D}}}]$ in (44) are conical for small k (i.e. in the vicinity of the diabolical point, which is the region selected by $\tilde {F}_{\mathbf {k}}$ ); in other words, $[\omega _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{\left (s=1,2\right )}-\omega _{_{\mathrm {D}}}]\approx \pm ck$ , where c is the group speed (the modulus of the group velocity) and $k=\left \vert \mathbf {k}\right \vert$ ; see the 2D Grover walk case in equation (52), where $c=1/\sqrt {2}$ . For the rest of the sheets, we have assumed $\omega _{\mathbf {k}_{\mathrm {D}}}^{\left (s\neq 1,2\right )}-\omega _{_{\mathrm {D}}}=0$ . Hence we can write, from (44),

$\begin{equation} \fl \left\vert \mathbf{F}_{\mathbf{x},t}^{\left(s=1,2\right)}\right\rangle =\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\exp\left[\mathrm{i}\left(\mathbf{k}\cdot\mathbf{x}\mp ckt\right)\right]\langle\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\mid\Xi\rangle\mid\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\rangle \end{equation} \tag{ 45 }$

while the rest of the partial waves verify

$\begin{equation} \fl \left\vert \mathbf{F}_{\mathbf{x},t}^{\left(s\neq1,2\right)}\right\rangle =\int\frac{\mathrm{d}^{N}\mathbf{k}}{\left(2\pi\right)^{N}}\exp\left[\mathrm{i}\left(\mathbf{k}\cdot\mathbf{x}\right)\right] \tilde{F}_{\mathbf{k}}\langle\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\mid\Xi\rangle\mid\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\rangle=\left\vert \mathbf{F}_{\mathbf{x},0}^{\left(s\neq1,2\right)}\right\rangle , \end{equation} \tag{ 46 }$

which is a constant, equal to its initial value. Hence, the projection of the initial state onto the sheets of constant frequency does not evolve, as expected, leading to a possible localization of a part of the initial wavepacket. We thus consider in the following that $\langle \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{\left (s\neq 1,2\right )}\mid \Xi \rangle =0$ .

Regarding the sheets s = 1,2, we cannot progress unless we assume some property of the eigenvectors $\vert \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{(s)}\rangle$ . We will assume that $\vert \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{\left (s=1,2\right )}\rangle$ (for k close to zero, i.e. around the diabolical point) depend just on the angular part of k (let us call it Ω) but not on its modulus k. We also restrict ourselves to cases in which $\skew3\tilde {F}_{\mathbf {k}}$ depends only on k (spherical symmetry), i.e. $\tilde {F}_{\mathbf {k}}=\tilde {F}_{k}$ . Hence we write the integral in (45) in (N-dimensional) spherical coordinates as

$\begin{equation} \left\vert \mathbf{F}_{\mathbf{x},t}^{(s=1,2)}\right\rangle =\left(2\pi\right)^{-N}\int\mathrm{d}k\exp\left(\mp\mathrm{i}ckt\right)k^{N-1}\tilde{F}_{k}\left\vert \mathbf{E}_{\mathbf{x}}^{(s)}(k)\right\rangle , \end{equation} \tag{ 47 }$

$\begin{equation} \left\vert \mathbf{E}_{\mathbf{x}}^{(s)}(k)\right\rangle =\int\mathrm{d}^{N}\Omega\exp\left(\mathrm{i}\mathbf{k}\cdot\mathbf{x}\right)\langle\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\mid\Xi\rangle\mid\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\rangle, \end{equation} \tag{ 48 }$

where we write d^Nk = k^N−1dk d^NΩ and d^NΩ is the N-dimensional solid angle element.

As mentioned at the end of the previous subsection, one might also study more complicated initial situations, by combining initial terms, each consisting of a different k₀. As an example, assume that k₀ can take two different values k₀₁ and k₀₂. We can face different situations. If both values are far from degeneracies, the discussion of section 3.1 applies. If one of them, say k₀₁, shows a degeneracy, the treatment developed in the present subsection has to be taken into account for that particular term in the initial state. Finally, one can have a situation in which both values are close to each other and lie close to a degeneracy point k₀₁ ≃ k₀₂ ≃ k_D. This case might lead to very complicated behavior, which is beyond the scope of this paper (in fact, as we show later, even an initial condition consisting of a single value k₀ close to a degeneracy can give rise to very rich phenomena).

Up to here the theory is somehow general. However, we cannot progress unless we particularize to some special case. We shall do this below for the 2D Grover walk.

4. Application to the two-dimensional Grover quantum walk

So far we have derived general expressions for the N-dimensional QW. In this section, we consider the special case of the 2D QW using the Grover coin operator.

4.1. Diagonalization of the two-dimensional Grover map: dispersion relations and diabolical points

The so-called Grover coin of dimension 2N has matrix elements C^α,η_α',η' = 1/N − δ_α,α'δ_η,η'. In the present case (N = 2), the corresponding matrix C_k (17), with k = (k₁,k₂), has the form

$\begin{equation} C_{\mathbf{k}}=\frac{1}{2}\left(\begin{array}{@{}cccc@{}}-\mathrm{e}^{-\mathrm{i}k_{1}} \mathrm{e}^{-\mathrm{i}k_{1}} & \mathrm{e}^{-\mathrm{i}k_{1}} & \mathrm{e}^{-\mathrm{i}k_{1}}\\ \mathrm{e}^{\mathrm{i}k_{1}} & -\mathrm{e}^{\mathrm{i}k_{1}} & \mathrm{e}^{\mathrm{i}k_{1}} & \mathrm{e}^{\mathrm{i}k_{1}}\\ \mathrm{e}^{-\mathrm{i}k_{2}} & \mathrm{e}^{-\mathrm{i}k_{2}} & -\mathrm{e}^{-\mathrm{i}k_{2}} & \mathrm{e}^{-\mathrm{i}k_{2}}\\ \mathrm{e}^{\mathrm{i}k_{2}} & \mathrm{e}^{\mathrm{i}k_{2}} & \mathrm{e}^{\mathrm{i}k_{2}} & -\mathrm{e}^{\mathrm{i}k_{2}} \end{array},\right) \end{equation} \tag{ 49 }$

whose diagonalization yields the eigenvalues $\lambda _{\mathbf {k}}^{(s)}=\exp (-\mathrm {i}\omega _{\mathbf {k}}^{(s)})$ with

$\begin{equation} \omega_{\mathbf{k}}^{(1)}=\pi+\Omega_{\mathbf{k}},\qquad\quad\omega_{\mathbf{k}}^{(2)}=\pi-\Omega_{\mathbf{k}},\qquad\quad \omega_{\mathbf{k}}^{(3)}=\pi,\ \omega_{\mathbf{k}}^{\left(4\right)}=0, \end{equation} \tag{ 50 }$

where

$\begin{equation} \Omega_{\mathbf{k}}=\arccos\left[\frac{1}{2}\left(\cos{k_{1}}+\cos{k_{2}}\right)\right]\in\left[0,\pi\right]. \end{equation} \tag{ 51 }$

Note that adding a multiple integer of 2π to any of the ω's does not change anything because time is discrete and runs in steps of 1; see (23).

We see, according to (50), that the last two eigenvalues $\lambda _{\mathbf {k}}^{\left (3,4\right )}=\exp(-\mathrm {i}\omega _{\mathbf {k}}^{\left (3,4\right )})=-1,1$ do not depend on k, and the first two eigenvalues $\lambda _{\mathbf {k}}^{\left (1,2\right )}=\exp (-\mathrm {i}\omega _{\mathbf {k}}^{\left (1,2\right )})=-\exp (\mp \mathrm {i}\Omega _{\mathbf {k}})$ are complex conjugates of each other; hence we can choose $\Omega _{\mathbf {k}}\in \left [0,\pi \right ]$ without loss of generality. A plot of the dispersion relations (50) is given in figure 1, where five degeneracy points (the origin and the four corners; see also figure 2) are observed with three-fold degeneracies. Moreover, the frequencies ±Ω_k have a conical form (diabolo-like) at those degeneracies and this is why we call them diabolical points, in analogy to other diabolos found in different systems [59–63, 67–70]. At the central diabolical point the degeneracy is between $\omega _{\mathbf {k}}^{\left (s=1,2,4\right )}$ , while at the corners it is between $\omega _{\mathbf {k}}^{\left (s=1,2,3\right )}$ .

**Figure 1.** (Main figure) Dispersion relation sketch around the diabolical point. The value of ω/π is plotted for the different sheets (i.e. different eigenvalues). The middle (flat) surface corresponds to ω⁽³⁾_k = π; the surface with values ranging from 1 to 2 corresponds to ω⁽¹⁾_k. The bottom plane is associated with ω⁽⁴⁾_k and the surface with values ranging from 0 to 1 corresponds to ω⁽²⁾_k. (Inset) Detail of ω⁽¹⁾_k and ω⁽²⁾_k around the origin (orange and blue curves, respectively) for k₂ = 0 as a function of k₁, showing the intersection at the diabolical point.
Download figure:
Standard image High-resolution image

**Figure 2.** Vector plot of the group velocity (shown as arrows), superimposed on a density plot of the velocity modulus equation (51).
Download figure:
Standard image High-resolution image

For the sake of later use, we note that the frequency Ω_k reads, close to the diabolical point k = 0,

$\begin{equation} \Omega_{\mathbf{k}}\simeq\frac{k}{\sqrt{2}}-\frac{k^{3}\cos\left(2\theta\right)}{48\sqrt{2}}+O(k^{5}), \end{equation} \tag{ 52 }$

where $\mathbf {k}=k\left (\cos \theta ,\sin \theta \right )$ . This means that, very close to the diabolical point, the frequency Ω_k actually has a conical dependence on the wavenumber. Later we will understand the consequences of the existence of these points.

The fact that $\omega _{\mathbf {k}}^{\left (s=3,4\right )}$ are constant results in $\vert \psi _{\mathbf {x},t}^{\left (s=3,4\right )} \rangle$ , see (23), not having wave character; only $\vert \psi _{\mathbf {x},t}^{\left (s=1,2\right )}\rangle$ are true waves. Hence propagation in the 2D Grover walk occurs only if the initial condition projects onto subspaces 1 or 2; otherwise the walker remains localized. This explains why strong localization in the 2D Grover map is usually observed, as first noted in [49]. All this can be put more formally in terms of the group velocity $\mathbf {v}_{\mathrm {g}}({\bf k})$ of the waves in the system, given by the gradient (in k) of the wave frequency. The group velocity, as in any linear wave system, has the meaning of velocity at which an extended wavepacket, centered in k-space around some value k₀, moves (there are other effects affecting extended wavepackets, which we analyze below). In our case, $\omega _{\mathbf {k}}^{\left (s=3,4\right )}$ are constant; hence their gradient is null: these eigenvalues entail no motion. In contrast, $\omega _{\mathbf {k}}^{\left (s=1,2\right )}$ depend on k ((50) and (51)) and then define non-null group velocities

$\begin{equation} \mathbf{v}_{\mathrm{g}}^{(1,2)}({\bf k})=\pm\nabla_{\mathbf{k}}\Omega_{\mathbf{k}}, \end{equation} \tag{ 53 }$

where the plus and minus signs stand for s = 1,2, respectively. Using (51) these velocities become

$\begin{equation} \mathbf{v}_{\mathrm{g}}^{\left(1,2\right)}({\bf k})=\pm\frac{(\sin{k_{1}},\sin{k_{2}})}{\sqrt{4-(\cos{k_{1}}+\cos{k_{2}})^{2}}}. \end{equation} \tag{ 54 }$

Hence the maximum group velocity modulus in the 2D Grover walk is $\frac {1}{\sqrt {2}}$ , as can be easily checked, and is obtained, in particular, along the diagonals k₂ = ± k₁. This is the maximum speed at which any feature of the QW can propagate. Figure 2 displays a vector plot of $\mathbf {v}_{\mathrm {g}}^{(1,2)}({\bf k})$ where a number of interesting points (in k space) are revealed. On the one hand, there are points of null group velocity, located at $\mathbf {k}=\left (0,\pm \pi \right )$ and at $\mathbf {k}=\left (\pm \pi ,0\right )$ , in correspondence with the saddle points observed in figure 1. On the other hand, there are five singular points at $\mathbf {k}=\left (0,0\right )$ as well as at the four corners. These points are singular because the group velocity is undefined at them (they look like sources or sinks). We note that these points coincide with the diabolical points in figure 1.

We find it worth mentioning that in 2D QWs with coin operators other than the Grover, diabolical points are also found, as we have checked by using the DFT coin [46] that displays several conical intersections. In this case of the DFT coin, none of the four leaves of the dispersion relation is constant. We mention too that behavior very similar to that of the 2D Grover walk is found in the alternate QW [40], except for the fact that in this last case there are no constant surfaces in the dispersion relation (and hence no localization phenomena).

4.2. Eigenvectors

The eigenvectors of the 2D Grover walk matrix (49) are given by Watabe et al [56]. Use of these eigenvalues in (23) allows us to finally compute the state of the system at any time.

Whenever k is not close to a diabolical point, these eigenvectors vary smoothly around k. Here we just want to study the behavior of the eigenvectors close to the diabolical point at $\mathbf {k}=\mathbf {k}_{\mathrm {D}}\equiv \left (0,0\right )$ —at the corners, the conclusions are similar. We find it convenient to use polar coordinates $\left (k_{1},k_{2}\right )=\left (k\cos \theta ,k\sin \theta \right )$ . Performing the limit for k → 0, we find

$\begin{eqnarray} \phi_{\mathbf{k}}^{\left(1\right)} & =&\frac{1}{2\sqrt{2}}\left(\begin{array}{@{}c@{}} 1+\sqrt{2}\cos\theta\\ 1-\sqrt{2}\cos\theta\\ -1-\sqrt{2}\sin\theta\\ -1+\sqrt{2}\sin\theta \end{array}\right),\nonumber \\ \phi_{\mathbf{k}}^{\left(2\right)} & =&\frac{1}{2\sqrt{2}}\left(\begin{array}{@{}c@{}} 1-\sqrt{2}\cos\theta\\ 1+\sqrt{2}\cos\theta\\ -1+\sqrt{2}\sin\theta\\ -1-\sqrt{2}\sin\theta \end{array}\right),\nonumber \\ \phi_{\mathbf{k}}^{\left(3\right)} & =&\frac{1}{\sqrt{2}}\left(\begin{array}{@{}c@{}} \sin\theta\\ -\sin\theta\\ \cos\theta\\ -\cos\theta \end{array}\right),\quad \phi_{\mathbf{k}}^{\left(4\right)}=\frac{1}{2}\left(\begin{array}{@{}c@{}} 1\\ 1\\ 1\\ 1 \end{array}\right). \end{eqnarray} \tag{ 55 }$

Hence, close to the diabolical point k_D, there are three eigenvectors, corresponding to s = 1,2,3, displaying a strong azimuthal dependence, while $\phi _{\mathbf {k}}^{\left (4\right )}$ is constant around k_D (its variations are smooth and hence tend to zero as k → k_D). Thus, even if only $\phi _{\mathbf {k}}^{\left (1\right )}$ and $\phi _{\mathbf {k}}^{\left (2\right )}$ participate in the conical intersection and define the diabolical point, $\phi _{\mathbf {k}}^{\left (3\right )}$ is also affected by this feature because it is associated with eigenvalue $\lambda _{\mathbf {k}}^{\left (3\right )}=-1$ , which is degenerate with $\lambda _{\mathbf {k}_{\mathrm {D}}}^{\left (s=1,2\right )}$ . The fact that the diabolical point is singular is easy to understand: depending on the direction we approach it, the eigenvectors of the problem are different.

Just at the diabolical point, $\lambda _{\mathbf {k}}^{\left (s=1,2,3\right )}=-1$ and hence there is a 3D eigensubspace formed by all vectors orthogonal to the fourth eigenvector, namely $\phi _{\mathbf {k}_{\mathrm {D}}}^{\left (4\right )}=\frac {1}{2}\mathrm {col}\left (1,1,1,1\right )$ , which is equal to $\phi _{\mathbf {k}}^{\left (4\right )}$ ; see (55). A sensible choice for these eigenvectors follows from noting that $\phi _{\mathbf {k}}^{\left (1\right )}+\phi _{\mathbf {k}}^{\left (2\right )}$ in (55) is independent of the azimuth and is orthogonal to $\phi _{\mathbf {k}}^{\left (4\right )}$ . Hence one of the basis elements for this 3D degenerate subspace can be chosen as

$\begin{equation} \phi_{_{\mathrm{D}}}=\frac{\phi_{\mathbf{k}}^{\left(1\right)}+\phi_{\mathbf{k}}^{\left(2\right)}}{\sqrt{2}}=\frac{1}{2}\,\mathrm{col}\left(1,1,-1,-1\right), \end{equation} \tag{ 56 }$

which has the outstanding property of being the single vector that projects, close to the diabolical point k_D, only onto $\phi _{\mathbf {k}}^{\left (s=1,2\right )}$ , i.e. onto the eigenspaces of the diabolo. As for the other two eigenvectors, we just impose their orthogonality with ϕ_{_D} and $\phi _{\mathbf {k}_{\mathrm {D}}}^{\left (4\right )}$ and we choose them arbitrarily as

$\begin{equation} \phi_{_{\mathrm{D}}}^{\prime}=\frac{1}{\sqrt{2}}\,\mathrm{col}\left(1,-1,0,0\right),\qquad\quad \phi_{_{\mathrm{D}}}^{\prime\prime}=\frac{1}{\sqrt{2}}\,\mathrm{col}\left(0,0,-1,1\right). \end{equation} \tag{ 57 }$

Note that we are using a 'primed' notation instead of keeping the notation with the label s, except for s = 4, because none of these eigenvectors keeps being an eigenvector for k ≠ k_D (only for special directions θ) and hence there is no sense in attaching any of these eigenvectors to a particular sheet of the dispersion relation. To conclude we list the projections of these eigenvectors onto the eigenvectors (55) in the neighborhood of k_D:

$\begin{equation} \begin{array}{@{}ll@{}} \left\langle\phi_{_{\mathrm{D}}}\mid\phi_{\mathbf{k}}^{\left(1\right)}\right\rangle=\left\langle\phi_{_{\mathrm{D}}}\mid\phi_{\mathbf{k}}^{\left(2\right)}\right\rangle=\frac{1}{\sqrt{2}}, & \left\langle\phi_{_{\mathrm{D}}}\mid\phi_{\mathbf{k}}^{\left(3\right)}\right\rangle=0,\\ \left\langle\phi_{_{\mathrm{D}}}^{\prime}\mid\phi_{\mathbf{k}}^{\left(1\right)}\right\rangle=-\left\langle\phi_{_{\mathrm{D}}}^{\prime}\mid\phi_{\mathbf{k}}^{\left(2\right)}\right\rangle=\frac{\cos\theta}{\sqrt{2}},\quad & \left\langle\phi_{_{\mathrm{D}}}^{\prime}\mid\phi_{\mathbf{k}}^{\left(3\right)}\right\rangle=-\sin\theta,\\ \left\langle\phi_{_{\mathrm{D}}}^{\prime\prime}\mid\phi_{\mathbf{k}}^{\left(1\right)}\right\rangle=-\left\langle\phi_{_{\mathrm{D}}}^{\prime\prime}\mid\phi_{\mathbf{k}}^{\left(2\right)}\right\rangle=\frac{\sin\theta}{\sqrt{2}}, & \left\langle\phi_{_{\mathrm{D}}}^{\prime\prime}\mid\phi_{\mathbf{k}}^{\left(3\right)}\right\rangle=\cos\theta \end{array} \end{equation} \tag{ 58 }$

and $\langle \phi _{_{\mathrm {D}}}\vert \phi _{\mathbf {k}}^{ (4)}\rangle =\langle \phi _{_{\mathrm {D}}}^{\prime }\vert \phi _{\mathbf {k}}^{(4)} \rangle =\langle \phi _{_{\mathrm {D}}}^{\prime \prime }\vert \phi _{\mathbf {k}}^{(4)}\rangle =0$ . All this has strong consequences on the QW dynamics near k_D, as we will show below.

4.3. Continuous wave equations

In order to illustrate how the evolution of the QW can be controlled via an appropriate choice of the initial conditions, we will show some results obtained numerically from the evolution of the 2D quantum walk with the Grover operator. Unless otherwise specified, the initial profile in position space is a Gaussian centered at the origin with cylindrical symmetry (independent of s):

$\begin{equation} \mid\!\!\psi_{\mathbf{x},0}\rangle=\mathrm{e}^{\mathrm{i}\mathbf{k}_{\mathrm{0}}\cdot\mathbf{x}}F_{\mathbf{x},0}\left\vert \phi\right\rangle , \end{equation} \tag{ 59 }$

where $\left \vert \phi \right \rangle$ represents the initial state of the coin, and

$\begin{equation} F_{\mathbf{x},0}=\mathcal{N}\mathrm{e}^{\mathrm{i}\mathbf{k}_{0}\cdot\mathbf{x}}\mathrm{e}^{-\frac{x_{1}^{2}+x_{2}^{2}}{2\sigma^{2}}} \end{equation} \tag{ 60 }$

with x = (x₁,x₂), which implies that in the continuous limit we also have a Gaussian, in momentum space, centered at k₀ = (k₁,k₂). Here, $\mathcal {N}$ is a constant that guarantees the normalization of the state. In the numerics, we have taken a sufficiently large value σ = 50 for the Gaussian, so as to make possible the connection with the continuous limit, although we will also discuss the situation for smaller values of σ in order to show the robustness of the results.

Let us first consider a simple case having k₁ = k₂ = π/2, with an initial coin state $\frac {1}{\sqrt {2}} (\mid\! \phi _{\mathbf {k}_{0}}^{(1)}\rangle +\mid\! \phi _{\mathbf {k}_{0}}^{(2)}\rangle )$ , where $\mid\! \phi _{\mathbf {k}_{0}}^{(1)}\rangle =\frac {1}{\sqrt {2}}\,\mathrm {col}\left (1,0,-1,0\right )$ and $\mid\! \phi _{\mathbf {k}_{0}}^{(2)}\rangle =\frac {1}{\sqrt {2}}\,\mathrm {col}\left (0,1,0,-1\right )$ that projects on the s = 1 and 2 branches with equal weights. The group velocity is given by $\mathbf {v}_{\mathrm {g}}^{(s)}(\mathbf {k}_{0})=\pm (\frac {1}{2},\frac {1}{2})$ , respectively, for s = 1,2, implying that the original Gaussian will split into two pieces that will move along the diagonal of the lattice in opposite directions, as shown in figure 3. It must be noted that, since the second derivatives vanish at this point (ϖ_ij = 0), the Gaussian distribution moves with no appreciable distortion. In the right panel of this figure, we show the result for σ = 5 and observe that, in this case, the results obtained for large σ are very robust even for such a low value of σ.

**Figure 3.** Probability distribution, as a function of the dimensionless (x₁,x₂) position and for the three times t = 0, 400 and 800. The initial condition (59) is given by the coin state $\left \vert \phi \right \rangle =\frac {1}{\sqrt {2}}(\mid \phi _{\mathbf {k}_{0}}^{(1)}\rangle +\mid \phi _{\mathbf {k}_{0}}^{(2)}\rangle )$ with k₁ = k₂ = π/2. (Left panel) 3D plot; the central peak corresponds to t = 0, and produces two symmetric peaks that move away from the origin: each pair of symmetric peaks corresponds to a different value of the time t. (Middle and right panels) Top view of the previous figure, for both σ = 50 (middle) and σ = 5 (right). The numbers indicate, in each case, the value of time for different snapshots.
Download figure:
Standard image High-resolution image

**Figure 3.** Probability distribution, as a function of the dimensionless (x₁,x₂) position and for the three times t = 0, 400 and 800. The initial condition (59) is given by the coin state $\left \vert \phi \right \rangle =\frac {1}{\sqrt {2}}(\mid \phi _{\mathbf {k}_{0}}^{(1)}\rangle +\mid \phi _{\mathbf {k}_{0}}^{(2)}\rangle )$ with k₁ = k₂ = π/2. (Left panel) 3D plot; the central peak corresponds to t = 0, and produces two symmetric peaks that move away from the origin: each pair of symmetric peaks corresponds to a different value of the time t. (Middle and right panels) Top view of the previous figure, for both σ = 50 (middle) and σ = 5 (right). The numbers indicate, in each case, the value of time for different snapshots.
Download figure:
Standard image High-resolution image

However, if we choose an initial value k₀ close to the origin, then the second derivatives take a large value, with the consequence that now the two pieces experience considerable distortion along the line perpendicular to the line of motion, as can be seen from figure 4, where k₁ = k₂ = 0.01π, implying a curvature ϖ_ij ≃ 7.96 in the perpendicular direction.

**Figure 4.** Probability distribution, as a function of the dimensionless (x₁,x₂) position and for the three times t =0, 400 and 800. The initial condition (59) is given by the coin state $\left \vert \phi \right \rangle =\frac {1}{\sqrt {2}}(\mid \phi _{\mathbf {k}_{0}}^{(1)}\rangle +\mid \phi _{\mathbf {k}_{0}}^{(2)}\rangle )$ with k₁ = k₂ = 0.01π. (Left panel) 3D plot; as in figure 3, the central peak corresponds to t = 0, and produces two symmetric distributions that move away from the origin: each pair of symmetric distributions corresponds to a different value of the time t. (Right panel) Top view of the previous figure, showing the value of the corresponding times.
Download figure:
Standard image High-resolution image

In this case, at variance with the case considered in figure 3, we have to face a peculiar situation, as the distribution is now peaked around a value close to the diabolical point $\mathbf {k}_{\mathrm {D}}=\left (0,0\right )$ . As long as the initial distribution does not contain this point with a significant probability (i.e. for a large value of σ), the behavior will be as described above. However, if σ decreases below a certain limit, the evolution is dominated by the singular behavior of the diabolical point (to be studied in detail in the next subsection). This transition is clearly observed in figure 5, where the resulting evolution is depicted as σ is increased from σ = 5 to 30, the latter subpanel already showing close resemblance to the result in figure 4. Distributions initially centered around non-singular points, however, are much more robust against a smaller σ, as previously shown in figure 3.

**Figure 5.** Same as figure 4, corresponding to four different values of σ. From left to right and from top to bottom: σ = 5, 10, 20 and 30, respectively. The initial condition (59) is given by the coin state $\left \vert \phi \right \rangle =\frac {1}{\sqrt {2}}(\mid \phi _{\mathbf {k}_{0}}^{(1)}\rangle +\mid \phi _{\mathbf {k}_{0}}^{(2)}\rangle )$ and k₁ = k₂ = 0.01π. The distribution is shown for three times t = 0, 400 and 800. The central peak corresponds to t = 0, and expands from the origin in a manner that depends on the value of σ. For σ = 5 and 10, the probability distribution has an almost cylindrical shape for each value of time. As σ is increased, the probability distribution becomes similar to the one found (for each value of t) in figure 4.
Download figure:
Standard image High-resolution image

A completely different situation arises when k₀ corresponds to one of the saddle points (such as the one located at $\mathbf {k}_{\mathrm {sp}}=\left (0,\pi \right )$ ). As the group velocity vanishes at that point, i.e. $\mathbf {v}_{\mathrm {g}}^{(1,2)}(\mathbf {k}_{0})=(0,0)$ , the center of the probability distribution is expected to remain at rest but will spread in time analogously to diffracting optical waves, as equation (37) reduces to the paraxial wave equation of optical diffraction [74]. In our case, by considering s = 1 as the initial coin state, one finds ϖ₁₁ = −1/2,ϖ₁₂ = ϖ₂₁ = 0,ϖ₂₂ = 1/2. Therefore, equation (37) becomes

$\begin{equation} \mathrm{i}\frac{\partial A^{(1)}({\bf X},t)}{\partial t}=-\frac{1}{4}\frac{\partial^{2}A^{(1)}}{\partial X_{1}^{2}}+\frac{1}{4}\frac{\partial^{2}A^{(1)}}{\partial X_{2}^{2}}. \end{equation} \tag{ 61 }$

This wave equation is similar to that of paraxial optical diffraction [74], differing from it in that the sign of spatial derivatives is different for the two spatial directions. Hence this hyperbolic equation describes a situation in which there is diffraction in one direction and anti-diffraction in the other direction, as occurs in certain optical systems [75]. This causes the cylindrical symmetry proper of diffraction to be replaced by a rectangular symmetry. Based on the analogy with diffraction, it is possible to tailor the initial distribution in order to reach a desired asymptotic distribution as we already demonstrated in the one-dimensional QW [13]. Hence if we want to reach a final homogeneous (top hat) distribution, we must use an initial condition of the form [13, 74]

$\begin{equation} F_{\mathbf{x},0}=\mathcal{N}\mathrm{e}^{\mathrm{i}\mathbf{k}_{0}\cdot\mathbf{x}}\mathrm{e}^{-\frac{x_{1}^{2}+x_{2}^{2}}{2\sigma^{2}}}\mathrm{sinc}\left(x_{1}/\sigma_{0}\right)\mathrm{sinc}\left(x_{2}/\sigma_{0}\right), \end{equation} \tag{ 62 }$

with $\mathrm {sinc}\left (x\right )=\sin \left (\pi x\right )/\pi x$ and σ₀ a constant that accounts for the width of the distribution. In figure 6, we observe the evolution of the Grover walk with the above initial condition and the initial coin $\vert \phi _{\mathbf {k_{0}},0}^{(1)} \rangle =\frac {1}{2}\mathrm {col}\left (1+\mathrm {i},1+\mathrm {i},1-\mathrm {i},1-\mathrm {i}\right )$ , equal for all populated sites. The figure shows a final distribution that is quite homogeneous along most of its support, except in its outer borders, where it shows a smooth but rapid fall-out (the ratio σ/σ₀ has been chosen to optimize the result [13]). This is the expected result, but figure 6 also shows a central peak that has been cut for a better display, the height of this central peak being quite large (around 1.5 × 10⁻⁵) as compared to the plateau (around 4 × 10⁻⁷). The existence of that central peak is a consequence of the common initial coin $\vert \phi _{\mathbf {k_{0}},0}\rangle$ . This coin state guarantees that the initial distribution projects just onto the relevant branch of the dispersion relation at its central spatial frequency but, as the width of the distribution is finite, for the larger values of $\left \vert \mathbf {k}-\mathbf {k}_{0}\right \vert$ it will provide a non-negligible projection onto the static branches of the dispersion relation and the localized part of these projections is what we see as a central peak in figure 6.

**Figure 6.** Probability distribution as a function of the dimensionless (x₁,x₂) position at time t = 9000, for initial conditions at the saddle point k₀ = (0,π). The coin initial state is $|\phi _{\mathbf {k}_{0}}^{(1)}\rangle =\frac {1}{2}{\rm col}(1+{{\rm i}},1+{{\rm i}},1-{{\rm i}},1-{{\rm i}})$ , and the starting space distribution is given by equation (62). We have taken a value σ₀ = 15 for this simulation. The central peak has been cut for a better display of the rest of the features.
Download figure:
Standard image High-resolution image

4.4. Dynamics at the diabolical points

In this case, N = 2, d^NΩ = dθ, with θ the azimuth, and

$\begin{equation} \left\vert \mathbf{E}_{\mathbf{x}}^{(s)}(k)\right\rangle =\int_{0}^{2\pi}\mathrm{d}\theta\exp\left(\mathrm{i}\mathbf{k}\cdot\mathbf{x}\right)\langle\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\mid\Xi\rangle\mid\phi_{\mathbf{k}_{\mathrm{D}}+\mathbf{k}}^{(s)}\rangle. \end{equation} \tag{ 63 }$

We recall that here $\vert \Xi \rangle$ is assumed to project just onto $\vert \phi _{\mathbf {k}_{\mathrm {D}}+\mathbf {k}}^{\left (s=1,2\right )}\rangle$ . This means that $\left \vert \Xi \right \rangle$ coincides with $\left \vert \phi _{\mathrm {D}}\right \rangle$ ; see (56). According to projections (58) and to vectors (55), we can write, in matrix form,

$\begin{equation} \mathbf{E}_{\mathbf{x}}^{(s)}(k)=\frac{1}{4}\int_{0}^{2\pi}\mathrm{d}\theta\exp\left(\mathrm{i}\mathbf{k}\cdot\mathbf{x}\right)\left(\begin{array}{@{}c@{}} 1\pm\sqrt{2}\cos\theta\\ 1\mp\sqrt{2}\cos\theta\\ -1\mp\sqrt{2}\sin\theta\\ -1\pm\sqrt{2}\sin\theta \end{array}\right), \end{equation} \tag{ 64 }$

where the upper (lower) sign in ± or ∓ corresponds to s = 1 (2), respectively, and we recall that $\mathbf {k}=(k\cos \theta ,k\sin \theta)$ . Upon writing $\mathbf {k}\cdot \mathbf {x}=kx\cos \left (\theta -\varphi \right )$ , where $\mathbf {x}=x\left (\cos \varphi ,\sin \varphi \right )$ , and performing, the integral, we arrive at

$\begin{equation} \mathbf{E}_{\mathbf{x}}^{(s)}(k)=\frac{\pi}{2}\left(\begin{array}{@{}c@{}} J_{0}(kx)\pm\mathrm{i}\sqrt{2}J_{1}\kern-1pt(kx)\cos\varphi\\ J_{0}(kx)\mp\mathrm{i}\sqrt{2}J_{1}\kern-1pt(kx)\cos\varphi\\ -J_{0}(kx)\mp\mathrm{i}\sqrt{2}J_{1}\kern-1pt(kx)\sin\varphi\\ -J_{0}(kx)\pm\mathrm{i}\sqrt{2}J_{1}\kern-1pt(kx)\sin\varphi \end{array}\right), \end{equation} \tag{ 65 }$

where J_n are Bessel functions of the first kind of order n. Using this result in (47), we obtain, in matrix form,

$\begin{equation} \mathbf{F}_{\mathbf{x},t}^{(s)}=(2\pi)^{-2}\int\mathrm{d}k\exp\left(\mp\mathrm{i}kct\right)k\tilde{F}_{k}\mathbf{E}_{\mathbf{x}}^{(s)}(k). \end{equation} \tag{ 66 }$

Let us compute the probability of finding the walker at some point, regardless of the coin state, as given by (6), (22), (43) and (66). After simple algebra, we obtain

$\begin{equation} P_{\mathbf{x},t}=\sum_{s,s^{\prime}}\left[\mathbf{F}_{\mathbf{x},t}^{(s)}\right]^{\dagger}\cdot\mathbf{F}_{\mathbf{x},t}^{\left(s^{\prime}\right)}=\left[p_{\mathbf{x},t}^{\left(0\right)}\right]^{2}+\left[p_{\mathbf{x},t}^{\left(1\right)}\right]^{2}, \end{equation} \tag{ 67 }$

where

$\begin{equation} p_{\mathbf{x},t}^{\left(0\right)} =(2\pi)^{-1}\int_{0}^{\infty}\mathrm{d}k\cos(kct)J_{0}(kx)k\tilde{F}_{k}, \end{equation} \tag{ 68 }$

$\begin{equation} p_{\mathbf{x},t}^{(1)} =(2\pi)^{-1}\int_{0}^{\infty}\mathrm{d}k\sin(kct)J_{1}\kern-1pt(kx)k\tilde{F}_{k}\end{equation} \tag{ 69 }$

are amplitude probabilities. We see that the spatial dependence of the probability amplitude on the azimuth φ in (66) and (65) disappears when looking at the full probability P_x,t. This means that a coin-sensitive measurement of the probability should display angular features.

As a summary of what we know up to here about the neighborhood of the diabolical point, we can say that (i) when the initial coin $\left \vert \Xi \right \rangle$ coincides with $\left \vert \phi _{\mathrm {D}}\right \rangle$ , see (56), we predict no localization and (ii) when the initial wavepacket F_x,0 is radially symmetric (it only depends on $x=\left \vert \mathbf {x}\right \vert$ ), we predict an evolution of the probability P_x,t that keeps the radial symmetry along time, according to (67)–(69).

Equations (68) and (69) describe the evolution of the probability of an initial wavepacket centered at a diabolical point and with radial symmetry. The actual probability P_x,t can be computed from them numerically, but little can be said in general. In order to gain some insight, we next consider the relevant long-time limit, which turns out to be analytically tractable. First of all, it is convenient to scale the radial wavenumber k to the (loosely defined) width of the initial state in real space, F_x,0, which we denote by σ. Accordingly, we define κ = σk so that (68) and (69) become

$\begin{equation} p_{\mathbf{x},t}^{\left(0\right)} =\left(2\pi\sigma^{2}\right)^{-1}\int_{0}^{\infty}\mathrm{d}\kappa\cos\left(\frac{ct}{\sigma}\kappa\right)J_{0}\left(\frac{x}{\sigma}\kappa\right)\kappa\tilde{f}_{\kappa}, \end{equation} \tag{ 70 }$

$\begin{equation} p_{\mathbf{x},t}^{\left(1\right)} =\left(2\pi\sigma^{2}\right)^{-1}\int_{0}^{\infty}\mathrm{d}\kappa\sin\left(\frac{ct}{\sigma}\kappa\right)J_{1}\left(\frac{x}{\sigma}\kappa\right)\kappa\tilde{f}_{\kappa}, \end{equation} \tag{ 71 }$

where $\tilde {f}_{\kappa }=\tilde {F}_{k=\kappa /\sigma }$ is non-null only for κ ≲ 1. Hence when ct ≫ σ, $\cos \left (\frac {ct}{\sigma }\kappa \right )$ and $\sin \left (\frac {ct}{\sigma }\kappa \right )$ are strongly oscillating functions of κ, around zero, which would make $p_{\mathbf {x},t}^{\left (0,1\right )}$ vanish. However in the integrals defining such amplitude probabilities, other oscillating functions, $J_{0}\left (\frac {x}{\sigma }\kappa \right )$ and $J_{1}\left (\frac {x}{\sigma }\kappa \right )$ , appear. Thus it can be expected that when the oscillations of the latter and those of the former are partially in phase, a non-null value of the integrals is obtained. According to the asymptotic behavior of $J_{n}\left (z\right )$ at large z, $J_{n}\left (z\right )\approx \sqrt {2/\pi z}\cos \left (z-n\pi /2-\pi /4\right )$ for $z\gg \left \vert n^{2}-1/4\right \vert$ [76], this partial phase matching occurs when x ≈ ct, which is expected from physical considerations. Hence we are led to consider the limit x = ct + σξ, with ct ≫ σ, where ξ is the scaled radial offset with respect to ct, in which case (70) and (71), after expressing the products of trigonometric functions as sums, become

$\begin{eqnarray} p_{\mathbf{x},t}^{\left(0\right)} & =&p_{\mathbf{x},t}^{\left(1\right)}\nonumber \\ & \approx&\left(2\pi\sigma^{2}\right)^{-1}\sqrt{\frac{\sigma}{2\pi ct}}\int_{0}^{\infty}\mathrm{d}\kappa\sin\left(\kappa\xi+\pi/4\right)\sqrt{\kappa}\tilde{f}_{\kappa}\nonumber \\ & =&\left(2\pi\right)^{-1}/\sqrt{2\pi ct}\int_{0}^{\infty}\mathrm{d}k\sin\left(k\sigma\xi+\pi/4\right)\sqrt{k}\tilde{F}_{k} \end{eqnarray} \tag{ 72 }$

to leading order, where the approximation follows from disregarding highly oscillating terms. As an application of this result, we next consider the Gaussian case $\tilde {F}_{k}=2\sigma \sqrt {\pi }\exp \left (-\frac {1}{2}\sigma ^{2}k^{2}\right )$ , where σ is the standard deviation of the initial probability in real space, equation (60). We readily obtain

$\begin{eqnarray} &&\fl p_{\mathbf{x},t}^{(0)} =p_{\mathbf{x},t}^{(1)}=\frac{1}{8\sqrt{\sigma ct}}\mathrm{e}^{-\xi^{2}/4} \left\{\vphantom{\frac{\xi^{2}I_{5/4}(\xi^{2}/4)-(\xi^{2}-2)I_{1/4}(\xi^{2}/4)}{\sqrt{\left\vert \xi\right\vert }}} \xi\sqrt{\left\vert \xi\right\vert }\left[I_{-1/4}(\xi^{2}/4)-I_{3/4}(\xi^{2}/4)\right]\right.\nonumber \\ &&+\left.\frac{\xi^{2}I_{5/4}(\xi^{2}/4)-(\xi^{2}-2)I_{1/4}(\xi^{2}/4)}{\sqrt{\left\vert \xi\right\vert }}\right\} , \end{eqnarray} \tag{ 73 }$

where $\xi =\frac {x-ct}{\sigma }$ and I_n is the modified Bessel function of the first kind of order n. The total probability P_x,t ≡ P(ξ) follows from using (73) in (67), from which we observe that the initial width σ acts only as a scale factor for the height and shape of the probability at long times (ct ≫ σ). A plot of the probability can be seen in figure 8, where two unequal spikes, whose maxima are located at ξ = −1.746 23 and at 0.550 855, separated by a zero at ξ = −0.765 951, are apparent. We recall that x = ct + σξ is a radial coordinate; hence these maxima correspond to two concentric rings separated by a dark ring. The situation is fully analogous to the so-called 'Poggendorff rings' appearing in the paraxial propagation of a beam incident along an optic axis of a biaxial crystal, a situation in which a conical intersection (a diabolical point) is present. Remarkably, the result we have obtained for P_x,t fully coincides with that described in [60] whose figure 7 is to be compared with our figure 8.

**Figure 7.** (Left panel) Probability distribution as a function of the dimensionless (x₁,x₂) position, at time t = 400. The initial condition is given by equation (56) around the diabolical point k_D = (0,0), with σ = 50. (Middle and right panels) Top view of the previous figure with σ = 50 (middle) and σ = 5 (right), but now at time t = 100.
Download figure:
Standard image High-resolution image

**Figure 8.** Radial cut, with x₂ = 0, for the top panel of figure 7 (blue, dotted curve), showing the Poggendorff rings (see the text for an explanation). Also shown for comparison is our analytical result, as obtained from equations (67) and (73) (solid, red curve).
Download figure:
Standard image High-resolution image

As an example that illustrates the above features, we have studied the evolution of a state of the form of equation (60) for the diabolical point k_D = (0,0), and a coin state as defined in equation (56). As seen from figure 7, the probability distribution keeps its cylindrical symmetry while it expands in position. It can also be seen from this figure that the general features remain valid for a wider distribution (in k space) corresponding to σ = 5, even if the distribution in spatial coordinates is narrower for a given time step. The Poggendorff rings can be appreciated in figure 8, which shows a radial cut of the previous figure. We also plot for comparison our analytical result, as obtained from equations (67) and (73), which shows good agreement with the (exact) numerical result, as can be seen from this figure.

All these features crucially depend on the choice of the coin initial conditions, as clearly observed when these conditions are chosen differently. As an example, we show in figure 9 the evolution, after 200 time steps, of the probability distribution when we start from a state, also centered around the diabolical point, but the vector coin is now $\vert \phi _{\mathbf {k}_{\mathrm {D}}}^{\left (2\right )}\rangle$ with θ = π/2. Most of the probability is directed along the positive x₂-axis, which corresponds to the choice of θ. In fact, by changing the value of this parameter one can, at will, obtain a chosen direction for propagation.

**Figure 9.** Probability distribution as a function of the dimensionless (x₁,x₂) position, at time t = 200. The initial condition is given by $\vert \phi _{\mathbf {k}_{\mathrm {D}}}^{(2)}\rangle$ around the diabolical point with θ = π/2.
Download figure:
Standard image High-resolution image

$\vert \phi _{\mathbf {k}_{\mathrm {D}}}^{(2)}\rangle $ — **Figure 9.** Probability distribution as a function of the dimensionless (x₁,x₂) position, at time t = 200. The initial condition is given by $\vert \phi _{\mathbf {k}_{\mathrm {D}}}^{(2)}\rangle$ around the diabolical point with θ = π/2.
Download figure:
Standard image High-resolution image

5. Application to the three-dimensional Grover quantum walk

In this section, we briefly analyze the 3D Grover QW in order to illustrate similarities and differences with the 2D case. In 3D, the Grover-coin operator reads

$\begin{equation} C=\frac{1}{3}\left(\begin{array}{@{}cccccc@{}} -2 & 1 & 1 & 1 & 1 & 1\\ 1 & -2 & 1 & 1 & 1 & 1\\ 1 & 1 & -2 & 1 & 1 & 1\\ 1 & 1 & 1 & -2 & 1 & 1\\ 1 & 1 & 1 & 1 & -2 & 1\\ 1 & 1 & 1 & 1 & 1 & -2 \end{array}\right) \end{equation} \tag{ 74 }$

and the diagonalization of the corresponding matrix C_k (17), with k = (k₁,k₂,k₃), yields six eigenvalues $\lambda _{\mathbf {k}}^{(s)}=\exp \left (-\mathrm {i}\omega _{\mathbf {k}}^{(s)}\right )$ with

$\begin{equation} \omega_{\mathbf{k}}^{(1,2)}=\pm\Omega_{\mathbf{k}}^{(+)},\qquad\quad\omega_{\mathbf{k}}^{(3,4)}=\pm\Omega_{\mathbf{k}}^{(-)},\qquad\quad \omega_{\mathbf{k}}^{(5)}=0,\ \omega_{\mathbf{k}}^{(6)}=\pi, \end{equation} \tag{ 75 }$

where

$\begin{equation} \fl \cos\Omega_{\mathbf{k}}^{\left(\pm\right)} =-\frac{1}{3}\left[\sum_{i=1,2,3}\cos{k_{i}}\right.\pm\left.\sqrt{\sum_{i=1,2,3}^{j>i}\left(\cos^{2}{k_{i}}-\cos{k_{i}\cos}k_{j}\right)}\right].\end{equation} \tag{ 76 }$

Recall that adding a multiple integer of 2π to any of the ω's does not change anything because time is discrete and runs in steps of 1. This implies that ω = −π and π represent the same frequency.

The graphical representation is more complicated in the 3D case. Plots of the dispersion relations (75) for particular values of k₃ and for particular values of ω are given in figures 10 and 11, respectively.

**Figure 10.** Representations of the Grover 3D six dispersion relation surfaces for k₃ = 0 (left), π/2 (center) and π (right).
Download figure:
Standard image High-resolution image

**Figure 11.** Representations of the Grover 3D isofrequency surfaces, as obtained from $\Omega _{\mathbf {k}}^{\left (+\right )}$ in the first octant of the Brillouin zone for ω = 0.95π (a), 0.8π (b), 0.6π (c), 0.4π (d), 0.2π (e) and 0.05π (f). For ω = π the two dispersion relations coincide with the cube axes and the cube origin, while for ω = 0 they have migrated to the opposite side of the cube.
Download figure:
Standard image High-resolution image

**Figure 11.** Representations of the Grover 3D isofrequency surfaces, as obtained from $\Omega _{\mathbf {k}}^{\left (+\right )}$ in the first octant of the Brillouin zone for ω = 0.95π (a), 0.8π (b), 0.6π (c), 0.4π (d), 0.2π (e) and 0.05π (f). For ω = π the two dispersion relations coincide with the cube axes and the cube origin, while for ω = 0 they have migrated to the opposite side of the cube.
Download figure:
Standard image High-resolution image

A large variety of propagation properties can be expected depending on the particular region in k-space that the initial distribution occupies. Note the existence of two constant dispersion relations, $\omega _{\mathbf {k}}^{\left (5,6\right )}$ , which will give rise to localization phenomena, as already discussed in the 2D case.

There exist several degeneracies in the above dispersion relations, appearing in the crossings of several dispersion curves. For ω = π, there is a five-fold degeneracy at k = 0, where $\omega _{\mathbf {k}}^{\left (j\right )}=\pi$ for j = 1,2,3,4,6. There is also a three-fold degeneracy occurring when k takes the form k = k_|| ≡ ke_i, where e_i is one of the vectors of the canonical basis in $\mathbb {R}^{3}$ (i.e. k_|| has two null components and an arbitrary third component), where $\omega _{\mathbf {k}}^{\left (j\right )}=\pi$ for j = 1,2,6. For ω = 0 these degeneracies migrate to the corners and borders, respectively, of the first Brillouin zone; see figure 11. The frequencies $\Omega _{\mathbf {k}}^{\left (\pm \right )}$ , close to the degeneracies for ω = π, read

$\begin{eqnarray}\fl\left.\cos\Omega_{\mathbf{k}}^{\left(\pm\right)}\right\vert _{\mathbf{k}=0} &\simeq&-1+\frac{1}{6}\left\vert \mathbf{k}\right\vert ^{2}\mp\frac{1}{6}\sqrt{\left\vert \mathbf{k}\right\vert ^{4}-3\left(k_{x}^{2}k_{y}^{2}+k_{x}^{2}k_{z}^{2}+k_{y}^{2}k_{z}^{2}\right)}+O\left(k^{4}\right). \end{eqnarray} \tag{ 77 }$

In the case that k takes the form k = k_|| ≡ ke_i, we obtain

$\begin{equation} \cos\Omega_{\mathbf{k}}^{\left(+\right)}=-1, \end{equation} \tag{ 78 }$

$\begin{equation} \cos\Omega_{\mathbf{k}}^{\left(-\right)}=-\frac{1}{3}(1+2\cos k), \end{equation} \tag{ 79 }$

which shows the above-mentioned degeneracy.

We see that the point degeneracy occurring at k = 0, $\cos \Omega _{\mathbf {k}}^{\left (+\right )}=-1$ , is not a true diabolical point because, for fixed frequency, the dispersion relation has no spherical symmetry (hence the dispersion relation surface is not a hyperdiabolo). This implies that the theory developed in section 3.2 cannot be applied to this case. (We note, however, that it can be applied when coins different from the Grover one are used, as occurs e.g. in the alternate 3D QW, where true 3D diabolical points exist [40].)

It is not difficult to see that the group velocities

$\begin{equation} \mathbf{v}_{\mathrm{g}}^{(1,2,3,4)}({\bf k})=\pm\nabla_{\mathbf{k}}\Omega_{\mathbf{k}}^{(\pm)} \end{equation} \tag{ 80 }$

now have the expression

$\begin{equation} \left(\mathbf{v}_{\mathrm{g}}^{(1,2)}({\bf k})\right)_{i} =\mp\frac{(2\cos\Omega_{\mathbf{k}}^{\left(+\right)}+\sum_{j\neq i}\cos k_{j})\sin k_{i}}{2\left(3\cos\Omega_{\mathbf{k}}^{\left(+\right)}+\sum_{j}\cos k_{j}\right)\sin\Omega_{\mathbf{k}}^{\left(+\right)}}, \end{equation} \tag{ 81 }$

$\begin{equation} \left(\mathbf{v}_{\mathrm{g}}^{(3,4)}({\bf k})\right) _{i}=\mp\frac{(2\cos\Omega_{\mathbf{k}}^{\left(-\right)}+\sum_{j\neq i}\cos k_{j})\sin k_{i}}{2\left(3\cos\Omega_{\mathbf{k}}^{\left(-\right)}+\sum_{j}\cos k_{j}\right)\sin\Omega_{\mathbf{k}}^{\left(-\right)}} \end{equation} \tag{ 82 }$

and $\mathbf {v}_{\mathrm {g}}^{(5,6)}({\bf k})=0$ .

We see that the extrema of the dispersion relations for branches 1–4 giving points of null velocity (those for which a Schrödinger-type equation is a priori well suited for extended initial conditions) occur, in particular, when all sin k_i = 0. At some of these points, we find the just discussed degeneracies; hence at them a Schrödinger-type equation is not appropriate. However, some of the points of null velocity do not correspond to degeneracies, in particular those k with two null components and the third one equal to ±π. Consider, for example, the point $\mathbf {k}_{0}=\left (0,0,\pm \pi \right )$ at which $\Omega _{\mathbf {k}}^{\left (+\right )}$ presents degeneracies but not $\Omega _{\mathbf {k}}^{\left (-\right )}$ , which takes the value $\Omega _{\mathbf {k}}^{\left (-\right )}=\arccos 1/3$ . In order to particularize the wave equation (1) to this case, we must calculate the nine coefficients

$\begin{equation} \ \varpi_{ij}^{\left(3,4\right)}=\pm\left.\partial^{2}\Omega_{\mathbf{k}}^{\left(-\right)}/\partial k_{i}\partial k_{j}\right\vert _{\mathbf{k}=\mathbf{k}_{0}} \end{equation} \tag{ 83 }$

which turn out to be $\varpi _{11}^{\left (3,4\right )}=\varpi _{22}^{\left (3,4\right )}=-\varpi _{33}^{\left (3,4\right )}/4=\mp 1/4\sqrt {2}$ , the rest being zero. Hence the continuous description in this case is given by

$\begin{equation} \mathrm{i}\frac{\partial A^{\left(3,4\right)}}{\partial t}=\pm\frac{1}{8\sqrt{2}}\left(\frac{\partial^{2}}{\partial X_{1}^{2}}+\frac{\partial^{2}}{\partial X_{2}^{2}}-4\frac{\partial^{2}}{\partial X_{3}^{2}}\right)A^{\left(3,4\right)}, \end{equation} \tag{ 84 }$

which exhibits anisotropic diffraction (diffraction in the $\left (X_{1},X_{2}\right )$ plane and anti-diffraction, with a different coefficient, in the X₃ direction). For $\mathbf {k}_{0}=\left (\pi ,0,0\right )$ and $\mathbf {k}_{0}=\left (0,\pi ,0\right )$ , the result is the same after making the changes X₁ → X₃ and X₂ → X₃, respectively. Equation (84) can be solved via a Fourier transform method. For a starting Gaussian profile as given below, the solution will factorize in three one-dimensional functions of the form

$\begin{equation} \frac{\mathrm{e}^{-\frac{\mathrm{i}x^{2}}{2t\varpi+2\mathrm{i}\sigma^{2}}}}{\sqrt{\sigma^{2}-\mathrm{i}t\varpi}} \end{equation} \tag{ 85 }$

(except for an overall constant), where ϖ is any of the coefficients ϖ_ij that appear in equation (84), and x represents one of the coordinates X_i with i = 1,2,3. This implies a characteristic time scale t ∼ σ²/|ϖ| for the Gaussian to broaden.

We now show two examples that illustrate the above results. As in the 2D cases, we start with an initial condition (41) defined by a Gaussian shape

$\begin{equation} F_{\mathbf{x},0}=\mathcal{N}\mathrm{e}^{\mathrm{i}\mathbf{k}_{0}\cdot\mathbf{x}}\mathrm{e}^{-\frac{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}{2\sigma^{2}}} \end{equation} \tag{ 86 }$

with $\mathcal {N}$ an appropriate constant defining the normalization and $\left \vert \Xi \right \rangle$ a constant (position-independent) vector in the coin space, which is chosen to be one of the eigenvectors of the matrix C_k (17) at the point k = k₀. We have not explored those cases where k₀ lies at, or very close to, any of the degeneracies discussed above. The problem at these points seems to be much more involved than in the 2D case. For example, the extension of the eigenvectors obtained in [56] to this case encounters singular behavior at the axis (which are degeneracy lines). One needs to carefully take care of the appropriate limit as one approaches these lines, and perform a systematic study in order to find the relevant linear combinations providing a sensible time evolution. Such a study is beyond the scope of this paper.

We first start with a value k₀ = (0.1,0.2,0.3)π that lies far from any degeneracy or points with zero group velocity. For s = 3, the group velocity is $\mathbf {v}_{\mathrm {g}}^{(3)}(\mathbf {k}_{0})=(-0.028,-0.232,-0.704)$ . In figure 12, we show two different views corresponding to the evolved probability distribution. Within the time scale displayed in this figure, the initial Gaussian has not appreciably broadened, and we can instead observe the motion of the center of the wavepacket according to what is expected from the group velocity. Somehow opposite behavior is observed if we choose a point giving zero velocity, such as k₀ = (0,0,π). As described above, we expect no motion of the center and a broadening that depends on the axis we choose. This broadening can be approximately described by equation (84) and, for the values of σ we are assuming here, one needs to follow the time evolution on a considerably larger time scale. As a consequence, higher-order effects neglected in deriving equation (84) may accumulate. In figure 13, we show two different radial cuts of the 3D probability distribution. As observed from this figure, the accumulated discrepancy between the exact numerical evolution and the approximated analytical equation (84) may depend on the spatial direction. However, given all these considerations, we see that our derived analytical result gives a reasonable description for the main features of the time evolution, even in 3D.

**Figure 12.** Top view of the probability distribution in position space on the x₂ = 0 plane (left) and on the x₁ = 0 plane (right) for a time step t = 125. The initial distribution is given by equations (41) and (86), with k₀ = (0.1,0.2,0.3)π and $\left \vert \Xi \right \rangle$ the eigenvector of the coin operator corresponding to s = 3.
Download figure:
Standard image High-resolution image

$\left \vert \Xi \right \rangle $ — **Figure 12.** Top view of the probability distribution in position space on the x₂ = 0 plane (left) and on the x₁ = 0 plane (right) for a time step t = 125. The initial distribution is given by equations (41) and (86), with k₀ = (0.1,0.2,0.3)π and $\left \vert \Xi \right \rangle$ the eigenvector of the coin operator corresponding to s = 3.
Download figure:
Standard image High-resolution image

**Figure 13.** Two different radial cuts of the evolved probability distribution, represented by P(r), where r is the coordinate along the cut, at t = 1000 when the initial distribution is a spherically symmetric Gaussian function with σ = 20, with k₀ = (0,0,π) and $\left \vert \Xi \right \rangle$ the eigenvector of the coin operator corresponding to s = 3. The figure also shows for comparison the analytical result as obtained from equation (84). The curve with red symbols corresponds to a cut along the x₃-axis, to be compared with the analytical approximation (yellow curve). The curve with blue symbols is the plot for a cut along the x₁ axis, to be compared with the analytical approximation (green).
Download figure:
Standard image High-resolution image

6. Conclusions

In this paper, we have formally solved the standard multidimensional QW, taking into account extended initial conditions of arbitrary shape and width. In the limit of large width, the continuous limit is particularly well suited and clearly reveals the wave essence of the walk. The use of the dispersion relations is the central argument behind this view as its analysis allows us to get much insight into the propagation properties of the walk, even allowing for the tailoring of the initial distribution in order to reach a desired asymptotic distribution, as we have demonstrated. We must stress here that our goal was not to obtain approximate solutions to the QW, thus quantifying their degree of accuracy, but rather to obtain qualitative insight into the long-term solutions. The equations we have derived will be quite accurate for wide initial conditions, but the interesting point is that they provide a good qualitative description for relatively narrow initial distributions, especially far from degeneracies.

We have also shown that the 2D Grover QW exhibits a diabolical point in its dispersion relations, and have analyzed in detail the dynamics around this point, closely connecting the Grover walk with the phenomenon of conical refraction [60, 61]. In the 3D Grover walk, we have found other types of degeneracies whose influence on the dynamics we have not illustrated. The detailed study of the construction of the eigenvectors providing a clear qualitative interpretation turns out to be much more complicated in the 3D than in the 2D case, and we leave this extensive study for future work.

We stress that the asymptotic distributions found in the continuous limit, valid for large initial distributions, can be qualitatively reached for not very broad distributions, say σ ∼ 5, and this should be easy to implement in systems such as the optical interferometers of [28, 37]. An exception to this general rule are those distributions that are peaked close to the diabolical point, since in this case a broad distribution will be dominated by the singular nature of that point.

As a final comment, we would like to say, on the one hand, that the multidimensional QW can be viewed as a simulator of a large variety of linear differential equations depending on the particular region of the dispersion relation that governs the evolution of the initial wavepacket. On the other hand, we stress that continuous approximations to the QW must follow the dispersion relation if they are to be taken as good approximations.

Acknowledgments

This work was supported by projects FIS2011-26960, FPA2011-23897 of the Spanish Government, Generalitat Valenciana grant PROMETEO/2009/128 and FEDER. AR acknowledges financial support from PEDECIBA, ANII, Universitat de València and Generalitat Valenciana.