Understanding and controlling N-dimensional quantum walks via dispersion relations. Application to the 2D and 3D Grover walks: Diabolical points and more

The discrete quantum walk in N dimensions is analyzed from the perspective of its dispersion relations. This allows understanding known properties, as well as designing new ones when spatially extended initial conditions are considered. This is done by deriving wave equations in the continuum, which are generically of the Schr\"odinger type, and allow devising interesting behaviors, such as ballistic propagation without deformation, or the generation of almost flat probability distributions, what is corroborated numerically. There are however special points where the energy surfaces display intersections and, near them, the dynamics is entirely different. Applications to the two- and three-dimensional Grover walks are presented.

Recently, some of us have been insisting in the interest of analyzing the dynamics of coined QWs from the perspective of their dispersion relation [13,40], which has been also used for purposes different to 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 on 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 insisting in that they allow to reach final probability distributions which, interestingly, can be tailored to some extent. It is also worth mentioning that extended initial conditions also make a connection with multiparticle quantum walks for noninteracting 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 -qudit 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 make a general treatment, valid for N dimensions, and illustrate our results with the special cases of the Grover two-and three-dimensional QWs. We notice that the 2D Grover-QW has received most of the theoretical attention paid to multidimensional QWs [46][47][48][49][50][51][52][53][54][55][56][57][58].
We shall derive continuous wave equations, whose form depend on the probability distribution at time t = 0. To the leading order, these equations are partial differential equations that can be written in the form with A (s) a continuous amplitude probability, coefficients (s) ij determined by the QW dispersion relation properties, and X i are the spatial coordinates in a reference frame moving with the group velocity (see below for full details). There is a lot of accumulated knowledge regarding the solutions of the above continuous equations and this knowledge allows, to some extent, for a 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, 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 degeneracies, namely conical intersections. We apply this treatment to the twodimensional 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][60][61], determine a type of dynamics similar to that of massless Dirac fermions or of electrons in graphene [62][63][64][65][66] and appear in very different systems: quantum triangular billiards [59], conical refraction in crystal optics [60,61], the already mentioned graphene electrons [62][63][64][65][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 the diabolical points. In the three-dimensional Grover walk we also find degeneracies, but these are different to the diabolical points. Diabolical points on the dispersion relations play also a role on 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 QW's, such 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 the next section we introduce the N -dimensional QW and solve it formally. In Section III we formally derive continuous equations in the general N -dimensional case, and we discuss how to treat conical intersections. In Sections IV and V we apply our results to the two-and threedimensional Grover QWs, respectively. Finally, in Sect. VI we give our main conclusions.

II. N-DIMENSIONAL DISCRETE QUANTUM WALKS. GENERALITIES
In the N -dimensional QW the walker moves at discrete time steps t ∈ N across an N -dimensional lattice of sites x ≡ (x 1 , . . . , x N ) ∈ Z N . The walker is endowed with a 2N -dimensional coin which, after a convenient unitary transformation, determines the direction of displacement. The Hilbert space of the whole system (walker+coin) has then the form where the position space, H P , is spanned by , and the coin space, H C , is spanned by 2N orthonormal quantum states {|α η : α = 1, . . . , N ; η = ±}. 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 |1 + and |1 − , which are the equivalent to the |R and |L (for right and left) states commonly used in the literature. This notation is introduced in order to deal easily with an arbitrary number of dimensions.
The state of the total system at time t is represented by the ket |ψ t , which can be expressed in the form where the projections are wave functions on the lattice. We find it convenient to define, at each point x, the following ket which is an (unnormalized) coin state, so that ψ α,η x,t = α η | ψ x,t . As ψ α,η x,t 2 = |( α η | ⊗ x|) |ψ t | 2 is the probability of finding the walker at (x, t), and the coin in state |α η , the probability of finding the walker at (x, t) irrespectively of the coin state is, then, where we used the fact that N α=1 The dynamical evolution of the system is ruled by where the unitary operator is given in terms of the identity operator in H P ,Î, and two more unitary operators. On the one handĈ is the so-called coin operator (an operator in H C ), which can be written in its more general form aŝ where the matrix elements C α,η α ,η ≡ α η |Ĉ α η can be arranged as a 2N × 2N unitary square matrix C. On the other handD is the conditional displacement operator in HD where u α is the unit vector along direction x α ; note that, depending on the coin state |α η , the walker moves one site to the positive or negative direction of x α if η = + or η = −, respectively. Projecting (7) onto x| and using (4) and (8)-(10) we get straightforwardly which further projected onto α η | leads to Equation (11), or equivalently (12), is the NDQW 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 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, for example, [72]). First we define the DFT pair where k = (k 1 , . . . , k N ) and k α ∈ [−π, π] is the (quasi-)momentum vector [73]. Applying the previous definitions to the map (11) we readily get where we defined a coin operator in the quasi-momentum spaceĈ k α = k · u α , whose matrix elements read Projection of (15) onto α η | and use of (16,17) leads tõ Hence the nonlocal maps (11,12) become local in the momentum representation (15,18). This allows solving formally the QW dynamics very easily because map (15) implies and hence the eigensystem ofĈ k (or of C k in matrix form) is most useful for solving the problem, as we do next.
As the operatorĈ k is unitary, its eigenvalues have all the form λ  . Once the eigensystem ofĈ k is known, implementing (19) is trivial: Given an initial distribution of the walker in position representation |ψ x,0 we compute its DFT ψ k,0 via (13), as well as the projections so that ψ k,0 = sf . Now recalling (19) we arrive to where we used λ , while in position representation we get, using (14), Hence the QW is formally solved: all we need is to compute the eigensystem ofĈ k and the initial state in reciprocal space ψ k,0 , which determines the weight functions f (s) k through (20). Equation (22) shows that the QW dynamics corresponds to the superposition of 2N independent walks, labeled by s. According to (23) the ω (s) k '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, iff k0 , which is an unnormalizable plane wave and thus unphysical.
In order to avoid possible confusions, 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 |1 + , |1 − , . . . |N + , |N − .

III. 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 which is appreciably larger than the lattice spacing (taken as unity in this work). This type of initial states are wavepackets which are easily expressed in reciprocal space, where k 0 is a reference (carrier) wavevector, chosen at will, φ k−k0 is centered at k = k 0 ) and having a very small width ∆k (s) π [74]. Notice that we have chosen an initial coin state which is independent of k. From here one must distinguish between regular points, where eigenvectors φ (s) k have a smooth dependence on k close to k 0 , 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 where F (s) is a wide and smooth function of x becauseF (s) k is concentrated around k = 0. Hence in real space our initial condition consists of a coin state φ x,0 . We see then that the type of initial conditions we are dealing with are very close to the plane waves of momentum k 0 of the QW, given by e ik0·x φ (s) k0 .

A. Regular points
For regular points (i.e., far from degeneracies) the initial condition (24) determines the weight functions (20) asf where we took into account that the eigenvectors φ (s) k0 vary smoothly around k = k 0 and that only for k ≈ k 0 the functionF (s) k−k0 is 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). Ap- using the same arguments as before, the partial waves ψ (s) x,t can be written as where we have defined new functions F , in agreement with (25).
are determined. Instead of doing so by brute force, i.e. by integrating -maybe numerically-Eq. (29), what we do now is to look for a continuous wave equation that, by comparison with other known cases, sheds light onto 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 Eq. (12) without any approximation.
To avoid confusion, we introduce a function F (s) (x, t) of continuous real arguments exactly in the same way as in (29) as there is nothing forbidding such a definition. We have then F (s) In order to simplify the derivation we rewrite Eq. (29) by making the variable change k − k 0 → 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 Let us obtain the (approximate) wave equation. First we take the time derivative of (30) and get Then we Taylor expand the function ω k0 (except in the exponent: otherwise large errors at long times would be introduced) around k = 0, etc. In the latter equation, v is the group velocity at the point k = k 0 (see discussion below). Now it is trivial to transform the so-obtained right hand side of Eq. (31) as a sum of spatial derivatives of F (s) (x, t), leading to the result which is the sought continuous wave equation. One should understand that, in general, few terms are necessary on the right hand side of Eq.
is a slowly varying function of space, as commented.
The first term on the right hand side of the wave equation is an advection term, implying that the initial condition F to the leading order. In fact we can get rid of this term by defining a moving reference frame X such that which is an N -dimensional Schrödinger-like equation, to the leading order. This means that the evolution of the wave packet consists of the advection of the initial wave packet 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 article as it governs the nontrivial 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 0 , the local variations of ω around k 0 determine the type of wave equation controlling the QW dynamics. Just two additional remarks: First, if the initial condition only projects onto one of the sheets of the dispersion relation (hence the sum in s in (24) reduces to a single element), all F (s) (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 the leading order. Second, if the initial state projects onto several eigenvectors the probability of finding the walker at (x, t), irrespectively of the coin state, follows from (6), (22) and (28), and reads to the leading order, where we used φ = δ s,s . Hence P x,t is just the sum of partial probabilities: there is not interference among the different sub-QW's.

B. 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 Eq. (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 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 by 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 ω (s) k are constant around the diabolical point.
In order to alleviate the notation we will denote by ω D the value of the degenerate frequencies at Notice that there can be other sheets with ω kD in the 2D Grover map). We consider an initial wave packet defined in reciprocal space by where |Ξ is an eigenstate of the coin operatorĈ k at k = k D , andF k−kD is again a narrow function centered at k = k D . In the position representation this initial condition reads, see (14), where Hence, as in the regular case, |Ξ represents the state of the coin at the initial condition. Clearly |Ξ is not unique because of the degeneracy. We do not use the notation φ (s) kD but instead |Ξ 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) asf We will assume that |Ξ 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 but 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 [77].
A main difference with the regular case is seen at this stage because there is not any choice of the initial coin state |Ξ making that only one value of s be populated (remind that φ (s) k vary strongly around k = k D ). This has consequences as we see next.
Now the system will evolve according to (23). The partial waves ψ (s) 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 ω (s =1,2) kD − ω D = 0 according to assumption (ii) above. As in the regular case, the evolution has a fast part, given by the carrier wave e i(kD·x−ω D t) , and a slow part (both in time and in space) given by F 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 φ (s) kD+k | Ξ = 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 F k0 , see (28). We will not try deriving such a wave equation but conform ourselves with trying to understand the evolution of the system under the assumed conditions.
What we can say is that the frequency offsets ω (s=1,2) (44) are conical for small k (i.e. in the vicinity of the diabolical point, which is the region selected byF k ); in other words, ω where c is the group speed (the modulus of the group velocity) and k = |k|; see the 2D Grover walk case in Eq.
For the rest of sheets we have assumed ω (s =1,2) kD − ω D = 0. Hence we can write, from (44) while the rest of partial waves verify 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 wave packet. We thus consider in the following that φ (s =1,2) kD+k | Ξ = 0. Regarding the sheets s = 1, 2 we cannot progress unless we assume some property of the eigenvectors φ (s) kD+k . We will assume that φ (s=1,2) kD+k (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 whichF k only depends on k (spherical symmetry), i.e.F k =F k . Hence we write the integral in (45) in (N -dimensional) spherical coordinates as where we wrote d N k = k N −1 dkd N Ω and d N Ω is the Ndimensional solid angle element. 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.

IV. APPLICATION TO THE 2D GROVER QW
So far we have derived general expressions for the Ndimensional QW. In this section we consider the special case of the two-dimensional QW using the Grover coin operator.
A. Diagonalization of the 2D 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 1 , k 2 ), has the form whose diagonalization yields the eigenvalues λ where 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 λ For the sake of later use we note that the frequency Ω k reads, close to the diabolical point k = 0, where k = k (cos θ, sin θ). 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 ω 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 Ref. [49]. All this can be put more formally in terms of the group velocity v g (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 wave packet, centered in k-space around some value k 0 , moves (there are other effects affecting extended wavepackets which we analyze below). In our case ω (s=3,4) k are constant, hence their gradient is null: these eigenvalues entail no motion. On the contrary ω (s=1,2) k depend on k (50,51) and then define non null group velocities: where the plus and minus signs stand for s = 1, 2, respectively. Using (51) these velocities become Hence the maximum group velocity modulus in the 2D Grover walk is 1 √ 2 , as can be easily checked and is obtained, in particular, along the diagonals k 2 = ±k 1 . This is the maximum speed at which any feature of the QW can propagate. Figure 2 displays a vector plot of v (1,2) g (k) where a number of interesting points (in k space) is revealed: On the one hand there are points of null group velocity, located at k = (0, ±π) and at k = (±π, 0), in correspondence with the saddle points observed in Fig. 1. On the other hand there are five singular points at k = (0, 0) 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 Fig. 1.
We find it worth mentioning that in two-dimensional QWs with coin operators other that 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 a 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).

B. Eigenvectors
The eigenvectors of the 2D Grover walk matrix (49) are given by [56]. Use of these eigenvalues in (23) allows finally to 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 k = k D ≡ (0, 0) -at the corners the conclusions are similar. We find it convenient to use polar coordinates (k 1 , k 2 ) = (k cos θ, k sin θ). Performing the limit for k → 0 we find Hence, close to the diabolical point k D , there are three eigenvectors, corresponding to s = 1, 2, 3, displaying a strong azimuthal dependence, while φ (4) k is constant around k D (its variations are smooth and hence tend to zero as k → k D ). Thus, even if only φ (1) k and φ (2) k participate in the conical intersection and define the diabolical point, φ k is also affected by this feature because it is associated with eigenvalue λ . The fact that the diabolical point is singular is easy to understand: depending on the direction we approach to it the eigenvectors of the problem are different.
Just at the diabolical point λ k : see (55). A sensible choice for these eigenvectors follows from noting that φ (55) is independent of the azimuth and is orthogonal to φ (4) k . Hence one of the basis elements for this 3D degenerate subspace can chosen as which has the outstanding property of being the single vector that projects, close to the diabolical point k D , only onto φ (s=1,2) k , i.e. onto the eigenspaces of the diabolo. As for the other two eigenvectors we just impose their orthogonality with φ D and φ (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 : and φ D φ All this has strong consequences on the QW dynamics near k D as we will show below.

C. 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 α or η ): with x = (x 1 , x 2 ), which implies that in the continuous limit we also have a Gaussian, in momentum space, centered at k 0 = (k 1 , k 2 ). Here, 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 it 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 1 = k 2 = π/2, with an initial coin state 1 1, 0, −1) that projects on the s = 1 and s = 2 branches with equal weights. The group velocity is given by v (s) g (k 0 ) = ±( 1 2 , 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 Fig. 3. It must be noticed that, since the second derivatives vanish at this point ( ij = 0), the Gaussian distribution moves with no appreciable distortion. On the right bottom panel of this figure, we show the result for σ = 5 and we observe that, in this case, the results obtained for large σ are very robust even for such a low value of σ.  However, if we choose an initial value k 0 close to the origin, then the second derivatives take a large value, with the consequence that now the two pieces experience a considerable distortion along the line perpendicular to the line of motion, as can be seen from Fig. 4, where k 1 = k 2 = 0.01π, implying a curvature ij 7.96 in the perpendicular direction.
In this case, at variance with the case considered in Fig. 3, we have to face a peculiar situation, as the distribution is now peaked around a value close to the diabolical point k D = (0, 0). As far 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 diabolic point (to be studied in detail in the next subsection). This transition is clearly observed in Fig. 5, where the resulting evolution is depicted as σ is increased from σ = 5 to σ = 30, the latter subpanel already showing a close resemblance with the result in Fig.  4. Distributions initially centered around non-singular points, however, are much more robust against a smaller σ, as previously shown in Fig. 3.
A completely different situation arises when k 0 corresponds to one of the saddle points (such as the one located at k sp = (0, π)). As the group velocity vanishes at that point, i.e., v (1,2) g (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 Eq. (37) reduces to the paraxial wave equation of optical diffraction [75]. In our case, by considering s = 1 as the initial coin state, one finds 11 = −1/2, 12 = 21 = 0, 22 = 1/2. Therefore, Eq. (37) becomes This wave equation is similar to that of paraxial optical diffraction [75], differing from it in that the sign of spatial derivatives are different for the two spatial direc- k 0 ) and k1 = k2 = 0.01π.
tions. Hence this hyperbolic equation describes a situation in which there is diffraction in one direction and antidiffraction in the other direction, as it happens in certain optical systems [76]. This causes that the cylindrical symmetry proper of diffraction is replaced by a rectangular symmetry. Basing 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,75] with sinc (x) = sin (πx) /πx and σ 0 a constant that accounts for the width of the distribution. In Fig. 6 we observe the evolution of the Grover walk with the above initial condition and the initial coin φ , 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 σ/σ 0 has been chosen to optimize the result [13]). This is the expected result, but Fig. 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 −5 ) as compared to the plateau (around 4 · 10 −7 ). The existence of that central peak is a consequence of the common initial coin |φ k0,0 . This coin state guarantees that the initial distribution projects just onto the relevant branch of the dispersion relation at its central Figure 6: Probability distribution as a function of the dimensionless (x1, x2) position at time t=9000, for initial conditions at the saddle point k0 = (0, π). The coin initial state is | φ , and the starting space distribution as given by Eq. (61). We have taken a value σ0 = 15 for this simulation. The central peak has been cut for a better display of the rest of the features. spatial frequency but, as the width of the distribution is finite, for the larger values of |k − k 0 | it will provide a non-negligible projection onto the static branches of the dispersion relation and the localized part of these projections are what we see as a central peak in Fig. 6.

D. Dynamics at the diabolical points
In this case N = 2, d N Ω = dθ, with θ the azimuth, and We remind that here |Ξ is assumed to project just onto φ (s=1,2) kD+k . This means that |Ξ coincides with |φ D , see (56). According to the projections (58) and to the vectors (55) we can write, in matrix form, where the upper (lower) sign in ± or ∓ corresponds to s = 1(2), respectively, and we remind that k = (k cos θ, k sin θ). Upon writing k · x = kx cos (θ − ϕ), where x = x (cos ϕ, sin ϕ) and performing the integral we arrive to where J n are the Bessel functions of the first kind of order n. Using this result in (47) we get, in matrix form, Let us compute the probability of finding the walker at some point, regardless the coin state, as given by (6), (22), (43) and (65). After simple algebra we obtain where p (0) are amplitude probabilities. We see that the spatial dependence of the probability amplitude on the azimuth ϕ in (65,64) 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 when the initial coin |Ξ coincides with |φ D , see (56), we predict (i) no localization and, (ii) when the initial wave packet F x,0 is radially symmetric (it only depends on x = |x|), an evolution of the probability P x,t that keeps the radial symmetry along time, according to (66,67,68).
Equations (67) and (68) describe the evolution of the probability of an initial wave packet 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 consider next 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 (67) and (68) become p (1) wheref κ =F k=κ/σ is non null only for κ 1. Hence when ct σ, cos ct σ κ and sin ct σ κ are strongly oscillating functions of κ, around zero, what would make p (0,1) x,t to vanish. However in the integrals defining such amplitude probabilities other oscillating functions, J 0 x σ κ and J 1 x σ κ , appear. Thus it can be expected that when the oscillations of the latter and the ones of the former are partially in phase, a non null value of the integrals is got. According to the asymptotic behavior of J n (z) at large z, J n (z) ≈ 2/πz cos (z − nπ/2 − π/4) for z n 2 − 1/4 [78], 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 (69) and (70), after expressing the products of trigonometric functions as sums, become to the leading order, where the approximation follows from disregarding highly oscillating terms. As an application of this result we consider next the Gaussian casẽ where σ is the standard deviation of the initial probability in real space, Eq. (59). We readily obtain p (0) x,t = p where ξ = x−ct σ 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 (72) in (66), 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 Fig. 8, where two unequal spikes, whose maxima are located at ξ = −1.74623 and at ξ = 0.550855, separated by a zero at ξ = −0.765951, are apparent. We remind 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 Fig. 7 is to be compared to our Fig. 8.
As an example that illustrates the above features, we have studied the evolution of a state of the form Eq. (59) for the diabolic point k D = (0, 0), and a coin state as defined in Eq. (56). As seen from Fig. 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 Fig. 8, which shows a radial cut of the previous figure.
We also plot for comparison our analytical result, as obtained form Eqs. (66) and (72) , which shows a 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 Fig. 9 the evolution after 200 time steps, of the probability distribution when we start from a state, also centered around the diabolic point, but the vector coin is now φ (2) kD with θ = π/2. Most of the probability is directed along the positive x 2 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.

V. APPLICATION TO THE 3D GROVER QW
In this section we briefly analyze the three-dimensional Grover QW in order to illustrate similarities and differences with the two-dimensional case. In 3D, the Grovercoin operator reads and the diagonalization of the corresponding matrix C k (17), with k = (k 1 , k 2 , k 3 ), yields six eigenvalues λ cos Ω Remember 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 (74) for particular values of k 3 and for particular values of ω are given in Figs. 10 and 11, respectively.
A large variety of propagation properties can be expected depending on the particular region in k-space that the initial distribution occupies. Notice the existence of  two constant dispersion relations, ω (5,6) k , 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 5-fold degeneracy at k = 0, where ω (j) k = π for j = 1, 2, 3, 4, 6. There is also a 3-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 R 3 (i.e., k || has two null components and an arbitrary third component), where ω (j) k = π for j = 1, 2, 6. For ω = 0 these degeneracies migrate to the corners and borders, respectively, of the first Brillouin zone, see Fig. 11. The frequencies Ω (±) k , close to the degeneracies for ω = π, read cos Ω In the case that k takes the form k = k || ≡ ke i , we obtain cos Ω cos Ω which shows the above-mentioned degeneracy. We see that the point degeneracy occurring at k = 0, cos Ω (+) k = −1, is not a true diabolical point because, for fixed frequency, the dispersion relation has not spherical symmetry (hence the dispersion relation surface is not a hyperdiabolo). This implies that the theory developed in Sec. III.B cannot be applied to this case. (We notice, however, that it can be applied when coins different to the Grover one are used as it occurs, for example, in the Alternate 3D QW, where true 3D diabolical points exist [40].) Its not difficult to see that the group velocities v (1,2,3,4) and v (5,6) g (k) = 0. We see that the extrema of the dispersion relations for branches 1,2,3 and 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 k 0 = (0, 0, ±π) at which Ω (+) k presents degeneracies but not Ω (−) k , which takes the value Ω (−) k = arccos 1/3. In order to particularize the wave equation (1) to this case we must calculate the nine coefficients which turn out to be 33 /4 = ∓1/4 √ 2 the rest being zero. Hence the continuous description in this case is given by i ∂A (3,4) ∂t which exhibits an anisotropic diffraction (diffraction in the (X 1 , X 2 ) plane and anti-diffraction, with a different coefficient, in the X 3 direction). For k 0 = (π, 0, 0) and k 0 = (0, π, 0) the result is the same after making the changes X 1 → X 3 and X 2 → X 3 , respectively. Eq. (83) 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 (except for an overall constant), where is any of the coefficients ij that appear in Eq. (83), and x represents one of the coordinates X i with i = 1, 2, 3. This implies a characteristic time scale t ∼ σ 2 /| | 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 with N an appropriate constant defining the normalization, and |Ξ 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 0 . We have not explored those cases where k 0 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 a 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 = (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 v  Fig. 12 we show two different views corresponding to the evolved probability distribution. Within the timescale displayed on this Figure, the initial Gaussian has not appreciably broaden, and we can instead observe the motion of the center of the wave packet according to what is expected from the group velocity. A somehow opposite behavior is observed if we choose a point giving a zero velocity, such as k 0 = (0, 0, π). As described above, we expect no motion of the center and a broadening which depends on the axis we choose. This broadening can be approximately described by Eq. (83) and, for the values of σ we are assuming here, one needs to follow the time evolution during a considerably larger time scale. As a consequence, higher order effects neglected in deriving Eq. (83) may accumulate. In Fig. 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 Eq. (83) 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 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 k0 = (0, 0, π) and |Ξ the eigenvector of the coin operator corresponding to s = 3. The figure also shows for comparison the analytical result as obtained from Eq. (83). The curve with red symbols corresponds to a cut along the x3 axis, to be compared with the analytical approximation (yellow curve). The curve with blue symbols is the plot for a cut along the x3 axis, to be compared with the analytical approximation (green).

VI. CONCLUSIONS
In this article 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 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 insist here in that our goal was not to obtain approximate solutions to the QW, thus quantifying their degree of accuracy, but rather to get 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, specially far from degeneracies.
We have also shown that the two dimensional 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 threedimensional 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 a future work.
We want to insist in 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 Ref. [28,37]. An exception to this general rule are those distributions which are peaked close to the diabolic 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 wave-packet. 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.
in any direction or, more formally, as the largest eigenvalue of the covariance matrix σ, such that σ 2 ij =´d N k ki −ki kj −kj ρ