Bifurcation Curves of Two-Dimensional Quantum Walks

TL;DR

This paper analyzes the asymptotic behavior of two-dimensional quantum walks, showing that regions of polynomial decay are bounded by explicitly computable algebraic bifurcation curves, enhancing understanding of their long-term dynamics.

Contribution

It introduces a method to explicitly compute bifurcation curves bounding polynomial decay regions in two-dimensional quantum walks.

Findings

Regions of polynomial decay are bounded by algebraic curves.

Explicit examples of bifurcation curves are provided for various quantum walks.

The results improve understanding of the asymptotic structure of quantum walk dynamics.

Abstract

The quantum walk differs fundamentally from the classical random walk in a number of ways, including its linear spreading and initial condition dependent asymmetries. Using stationary phase approximations, precise asymptotics have been derived for one-dimensional two-state quantum walks, one-dimensional three-state Grover walks, and two-dimensional four-state Grover walks. Other papers have investigated asymptotic behavior of a much larger set of two-dimensional quantum walks and it has been shown that in special cases the regions of polynomial decay can be parameterized. In this paper, we show that these regions of polynomial decay are bounded by algebraic curves which can be explicitly computed. We give examples of these bifurcation curves for a number of two-dimensional quantum walks.

Figures

Equations

$$Q\Psi =\left[\begin{array}{ccc}{T}_{\_}{\sigma }_{\_}1& & \\ & \ddots & \\ & & T_{\_}{\sigma }_{\_}n\end{array}\right]U\Psi$$

$$\mathcal{F}[T_{\_}\sigma \psi ](\theta )=e^{i\sigma \cdot \theta }\mathcal{F}[\psi ](\theta )$$

$$\mathcal{F}[Q\Psi ](\theta )=\left[\begin{array}{ccc}{e}^{i{\sigma }_{\_}1\cdot \theta }& & \\ & \ddots & \\ & & e^{i{\sigma }_{\_}n\cdot \theta }\end{array}\right]U\mathcal{F}[\Psi ]=M(\theta )\mathcal{F}[\Psi ](\theta )$$

$$M(\theta {)}^t\mathcal{F}[{\Psi }_{\_}0](\theta )=\sum _{\_}{m=1}^n{\lambda }_{\_}j(\theta {)}^t{E}_{\_}m(\theta )$$

$$Q^t{\Psi }_{\_}0=\frac{1}{(2\pi {)}^d}\sum _{\_}{m=1}^n{\int }_{\_}\Vert \theta {\Vert }_{\_}\infty \le \pi E_{\_}m(\theta )\exp \left[i(t{H}_{\_}m(\theta )-x\cdot \theta )\right]d\theta$$

$$Q^t{\Psi }_{\_}0=\frac{1}{(2\pi {)}^d}\sum _{\_}{m=1}^n{\int }_{\_}\Vert \theta {\Vert }_{\_}\infty \le \pi E_{\_}m(\theta )\exp \left[it(H_{\_}m(\theta )-X\cdot \theta )\right]d\theta$$

$$I(t)={\int }_{\_}{\mathbb{R}}^{d}g(x)e^{itf(x)}dx.$$

$$S=\left\{x\in {\mathbb{R}}^d:\nabla f(x)=0,\det \left(\left[\frac{{\partial }^{2}f}{\partial x_{\_}j\partial x_{\_}k}\right]\right)\ne 0\right\}$$

$$\det \left(\left[\frac{{\partial }^2{H}_{\_}m}{\partial {\theta }_{\_}j\partial {\theta }_{\_}k}\right](\theta )\right)=0$$

$$\nabla H_{\_}m(\theta )=X$$

$$p\left(\lambda (\theta );x_{\_}k(\theta )\right)=0$$

$$\frac{\partial \lambda }{\partial {\theta }_{\_}j}=i\frac{\partial H}{\partial {\theta }_{\_}j}{e}^{iH(\theta )}=i\frac{\partial H}{\partial {\theta }_{\_}j}\lambda (\theta )$$

$$H_{\_}j=\frac{{\lambda }_{\_}j}{i\lambda }$$

$${\lambda }_{\_}j-i\lambda X_{\_}j=0$$

$${\lambda }_{\_}jk=i(H_{\_}jk\lambda +H_{\_}j{\lambda }_{\_}k)$$

$$H_{\_}jk=\frac{i}{{\lambda }^2}\left({\lambda }_{\_}j{\lambda }_{\_}k-{\lambda }_{\_}jk\lambda \right)$$

$$\det \left(\left[{\lambda }_{\_}j{\lambda }_{\_}k-{\lambda }_{\_}jk\lambda \right]\right)=0$$

$$\frac{\partial p}{\partial \lambda }\frac{\partial \lambda }{\partial {\theta }_{\_}j}+\frac{\partial p}{\partial x_{\_}j}\frac{\partial x_{\_}j}{\partial {\theta }_{\_}j}=0$$

$${\lambda }_{\_}j=-\frac{i{p}_{\_}j{x}_{\_}j}{p_{\_}\lambda }$$

$$p_{\_}j{x}_{\_}j+\lambda p_{\_}\lambda X_{\_}j=0$$

$$\left(\frac{\partial p_{\_}\lambda }{\partial \lambda }\frac{\partial \lambda }{\partial {\theta }_{\_}k}+\frac{\partial p_{\_}\lambda }{\partial x_{\_}k}\frac{\partial x_{\_}k}{\partial {\theta }_{\_}k}\right)\frac{\partial \lambda }{\partial {\theta }_{\_}j}+\frac{\partial p}{\partial \lambda }\frac{{\partial }^2\lambda }{\partial {\theta }_{\_}j\partial {\theta }_{\_}k}+\left(\frac{\partial p_{\_}j}{\partial \lambda }\frac{\partial \lambda }{\partial {\theta }_{\_}k}+\frac{\partial p_{\_}j}{\partial x_{\_}k}\frac{\partial x_{\_}k}{\partial {\theta }_{\_}k}\right)\frac{\partial x_{\_}j}{\partial {\theta }_{\_}j}=0$$

$${\lambda }_{\_}jk=\frac{x_{\_}j{x}_{\_}k}{p_{\_}{\lambda }^3}\left[p_{\_}\lambda \lambda p_{\_}j{p}_{\_}k-p_{\_}\lambda k{p}_{\_}j{p}_{\_}\lambda -p_{\_}{\lambda }_{\_}j{p}_{\_}k{p}_{\_}\lambda +p_{\_}jk{p}_{\_}{\lambda }^2\right]$$

$$\left(\frac{\partial p_{\_}\lambda }{\partial \lambda }\frac{\partial \lambda }{\partial {\theta }_{\_}j}+\frac{\partial p_{\_}\lambda }{\partial x_{\_}j}\frac{\partial x_{\_}j}{\partial {\theta }_{\_}j}\right)\frac{\partial \lambda }{\partial {\theta }_{\_}j}+\frac{\partial p}{\partial \lambda }\frac{{\partial }^2\lambda }{\partial {\theta }_{\_}{j}^2}+\left(\frac{\partial p_{\_}j}{\partial \lambda }\frac{\partial \lambda }{\partial {\theta }_{\_}j}+\frac{\partial p_{\_}j}{\partial x_{\_}j}\frac{\partial x_{\_}j}{\partial {\theta }_{\_}j}\right)\frac{\partial x_{\_}j}{\partial {\theta }_{\_}j}+\frac{\partial p}{\partial x_{\_}j}\frac{{\partial }^2{x}_{\_}j}{\partial {\theta }_{\_}{j}^2}=0$$

$${\lambda }_{\_}jj=\frac{x_{\_}j}{p_{\_}{\lambda }^3}\left[p_{\_}\lambda \lambda p_{\_}{j}^2{x}_{\_}j-2p_{\_}\lambda j{p}_{\_}j{p}_{\_}\lambda x_{\_}j+p_{\_}jj{p}_{\_}{\lambda }^2{x}_{\_}j+p_{\_}j{p}_{\_}{\lambda }^2\right]$$

$$\text{Res}(f,g;x)=F(a_{\_}k,b_{\_}k)=\det \left[\begin{array}{cccccccc}{a}_{\_}n& a_{\_}n-1& \dots & a_{\_}0& 0& & & \\ b_{\_}n& b_{\_}n-1& \dots & b_{\_}0& 0& \ddots & & \\ 0& a_{\_}n& \dots & a_{\_}1& a_{\_}0& \ddots & & \\ 0& b_{\_}n& \dots & b_{\_}1& b_{\_}0& \ddots & & \\ & & \ddots & \ddots & \ddots & \ddots & & \\ & & & & & a_{\_}n& \dots & a_{\_}0\\ & & & & & b_{\_}n& \dots & b_{\_}0\end{array}\right]$$

$$p_{\_}0(\lambda ;x,y)={\lambda }^4+x(\theta ){\lambda }^3+y(\theta ){\lambda }^2+x(\theta )\lambda +1$$

$$p_{\_}0(\lambda ;x,y)=({\lambda }^2+a(\theta )\lambda +1)({\lambda }^2+b(\theta )\lambda +1)$$

$$\nabla H(\theta )=\frac{a(\theta )\nabla x(\theta )-\nabla y(\theta )}{(2a(\theta )-x(\theta ))\sqrt{4-a(\theta {)}^2}}$$

$$G_{\_}4=\frac{1}{2}\left[\begin{array}{cccc}-1& 1& 1& 1\\ 1& -1& 1& 1\\ 1& 1& -1& 1\\ 1& 1& 1& -1\end{array}\right]$$

$$M({\theta }_{\_}1,{\theta }_{\_}2)=\frac{1}{2}\left[\begin{array}{cccc}-e^{i{\theta }_{\_}1}& e^{i{\theta }_{\_}1}& e^{i{\theta }_{\_}1}& e^{i{\theta }_{\_}1}\\ e^{-i{\theta }_{\_}1}& -e^{-i{\theta }_{\_}1}& e^{-i{\theta }_{\_}1}& e^{-i{\theta }_{\_}1}\\ e^{i{\theta }_{\_}2}& e^{i{\theta }_{\_}2}& -e^{i{\theta }_{\_}2}& e^{i{\theta }_{\_}2}\\ e^{-i{\theta }_{\_}2}& e^{-i{\theta }_{\_}2}& e^{-i{\theta }_{\_}2}& -e^{-i{\theta }_{\_}2}\end{array}\right]$$

Full text

1 Introduction

The quantum walk is a discrete quantum mechanical system which serves as an analogue to the classical random walk. While the probability distribution of a classical random walk approximates a normal distribution with a linearly growing variance in large time, the quantum walk has more involved asymptotic behavior. A paper by Ambainis et. al. [3] provides a detailed asymptotic description of the one-dimensional Hadamard walk (a quantum walk governed by a $2\times 2$ Hadamard matrix) for large time. The large time behavior of the quantum walk contrasts with that of the classical random walk in both the shape of the probability density and its rate of spread. Ambainis et. al. showed that the standard deviation of the Hadamard walk grows at $O(t)$ and in particular the distribution at time $t$ is almost entirely contained in the interval $\left[-\frac{t}{\sqrt{2}},\frac{t}{\sqrt{2}}\right]$. In this interval, the distribution contains highly oscillatory peaks at the boundaries and is approximately uniform near the origin.

These results were among the first in a line of papers describing asymptotic behavior of more general quantum walks. A few years later, Konno et. al. [14] extended these results to two state quantum walks governed by general unitary matrices. Inui et. al. [12] provided asymptotic results on the three-state one dimensional Grover walk with a focus on localization phenomena. Analysis of two-dimensional quantum walks started with a brief survey by Mackay et. al. [19], and continued with a description of localization in the two-dimensional Grover walk by Inui et. al. [11], and a differential geometric interpretation of two-dimensional quantum walks by Baryshnikov et. al. [5]. The asymptotic behavior of the one dimensional Hadamard walk described by Ambainis et. al. has served as a common baseline through these papers. The standard deviation of these walks grows as $O(t)$, and the probability distributions are almost entirely contained in a linearly expanding subset of the domain, which we term as the region of polynomial decay. In Kuklinski [17], these regions of polynomial decay are investigated further. In particular, it was shown that in certain low state two-dimensional quantum walks, the regions of polynomial decay can be explicitly parameterized. However, these parameterizations are often unwieldy and cannot easily give a description of the bounds on the region.

In this paper, we show that for any quantum walk defined on $\mathbb{Z}^{d}$, one can write down a collection of algebraic curves which bound the regions of polynomial decay. This procedure builds upon previous methods to analyze the asymptotic behavior of the quantum walk. First, we conduct an eigenvalue decomposition of the quantum walk operator via a Fourier transform. The state of the quantum walk particle at time $t$ becomes a sum of integrals whose asymptotic behavior we analyze using the method of stationary phase. Solving for points of stationary phase in these integrals leads to the previously discussed parametric representation of the regions of polynomial decay. We will take this a step further and compute bifurcation curves corresponding to these integrals by computing a Hessian determinant of the phase of the integrand. These bifurcation curves divide the space $\mathbb{R}^{d}$ into a discrete collection of subsets, a finite number of which are found to belong to the region of polynomial decay. The bifurcation curves are found to be solutions to a system of multivariate polynomial equations. We use a Gröbner basis computation to derive an implicit algebraic representation of these curves. However, this system of multivariate polynomial equations is quite large which puts strain on our algorithm, thus only in the simplest cases can we derive the bifurcation curves. We present a non-rigorous ad-hoc method for computing bifurcation curves for more complicated quantum walks.

The remainder of this paper is organized as follows. The definition of the quantum walk as well as a discussion of stationary phase approximations is given in section 2. In section 3 we explicitly write the multivariate polynomial system corresponding to the bifurcation curves and discuss solution techniques. In section 4 we compute bifurcation curves for several examples of two-dimensional quantum walks.

2 Definitions and Method of Stationary Phase

We begin by defining the quantum walk on a group, as first introduced by Acevedo et. al. [2]:

Definition 2.1

Let $(G,\cdot)$ be a group, let $\Sigma\subset G$ with $|\Sigma|=n$, and let $U$ be an $n\times n$ unitary matrix. The quantum walk operator $Q:\ell^{2}(G\times\Sigma)\rightarrow\ell^{2}(G\times\Sigma)$ corresponding to the triple $(G,\Sigma,U)$ may be written as the composition $Q=T(I\otimes U)$ where for $g\in G$ and $\sigma\in\Sigma$, $T:|g\rangle|\sigma\rangle\mapsto|\sigma\cdot g\rangle|\sigma\rangle$. We denote this correspondence as $Q\leftrightarrow(G,\Sigma,U)$.

The ordered pair $(G,\Sigma)$ can be thought of as an undirected Cayley graph which admits loops [8]. In this paper, we will primarily consider $G=\mathbb{Z}^{2}$. Let $C_{\_}d\subset\mathbb{Z}^{d}$ be the set of unit directional vectors and let $\tilde{C_{\_}d}=C_{\_}d\cup\{0\}$. Two of the unitary matrices which we will use in this paper are the Grover matrix and the Hadamard matrix. If $I_{\_}n$ is the $n\times n$ identity matrix and $\text{\bf 1}_{\_}n$ is the $n\times n$ matrix filled with ones, then we define the $n\times n$ Grover matrix as $G_{\_}n=\frac{2}{n}\text{\bf 1}_{\_}n-I_{\_}n$. The Hadamard matrix is defined as $H=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\ 1&-1\end{bmatrix}$.

One way to view this quantum walk operator is as a linear combination of translations. Let $\psi\in\ell^{2}(G)$ and let $T_{\_}\sigma:\ell^{2}(G)\rightarrow\ell^{2}(G)$ be a translation operator which acts as $T_{\_}\sigma(\psi(g))=\psi(\sigma^{-1}\cdot g)$ with $\sigma\in\Sigma$. Then we can visualize $Q$ acting on the vector $\Psi=\left[\psi_{\_}{\sigma_{\_}1},...,\psi_{\_}{\sigma_{\_}n}\right]$ as follows:

[TABLE]

When $G=\mathbb{Z}^{d}$, we can gain a better understanding of $Q$ through the application of a discrete Fourier transform. If $\theta\in\mathbb{R}^{d}$, the multi-dimensional Fourier transform acts on a translation as:

[TABLE]

Applying a Fourier transform to equation (1) and using equation (2), we find:

[TABLE]

We refer to $M(\theta)$ as the multiplier matrix of the quantum walk operator.

After $t$ steps of the quantum walk, the Fourier transform of the state becomes

$\mathcal{F}[Q^{t}\Psi](\theta)=M(\theta)^{t}\mathcal{F}[\Psi](\theta)$. Thus, if we wish to study long term behavior of the quantum walk, we must take large powers of the multiplier matrix. If $\{\lambda_{\_}m(\theta)\}_{\_}{m=1}^{n}$ is the set of eigenvalues of $M(\theta)$, then we can write the transform of a quantum walk with initial condition $\Psi_{\_}0$ as:

[TABLE]

Here, the $E_{\_}m(\theta)$ are scaled eigenvectors of $M(\theta)$ whose weights depend on the initial condition. Since $M(\theta)$ is unitary, we can write $\lambda_{\_}m(\theta)=e^{iH_{\_}m(\theta)}$ where $H_{\_}m:[-\pi,\pi]^{n}\rightarrow\mathbb{R}$. To return to the original domain, we conduct an inverse Fourier transform on equation (4):

[TABLE]

Here, $\lVert\cdot\rVert_{\_}\infty$ refers to the $\ell^{\infty}$ norm [16]. We let $x=Xt$ such that our spatial variable of interest is now $X\in\mathbb{R}^{d}$. This scales position space such that we will no longer be observing a linearly expanding spatial region but a stationary one. Substituting this into equation (5), we have:

[TABLE]

We use the method of stationary phase [6] to asymtptotically evaluate the integrals in equation (6). Consider the following $d$-dimensional oscillatory integral:

[TABLE]

where $g$ and $f$ are smooth functions with compact support. Consider the set:

[TABLE]

whose members we refer to as points of stationary phase or nondegenerate critical points. If $\nabla f(x)=0$ and the determinant of the Hessian matrix also vanishes, we say that $x$ is a degenerate critical point. We state two results from Stein [22] relating this set to an approximation of $I(t)$ in equation (7):

Proposition 2.1

Suppose $f$ has no critical points in the support of $g$. Then $I(t)=O(t^{-n})$ for every $n\geq 0$.

This proposition says that if the set of critical points is empty, then the integral $I(t)$ in equation (7) decays superpolynomially (often this superpolynomial decay can be shown to be exponential). If the set of critical points is nonempty, then we can use the following proposition:

Proposition 2.2

If $f$ has a nondegenerate critical point at $x_{\_}0$ and $g$ is supported in a sufficiently small neighborhood of $x_{\_}0$, then $I(t)=O(t^{-d/2})$.

Thus, if the set of nondegenerate critical points is nonempty, then the integral $I(t)$ decays polynomially, slower than the decay dictated by proposition 1.

We apply these propositions to the integrals in equation (6). Let $f_{\_}j(\theta)=H_{\_}j(\theta)-X\cdot\theta$ correspond to the expression in the exponential for the $j^{\text{th}}$ integrand. This integral has a critical point if $\nabla f_{\_}j(\theta)=0$, or if $\nabla H_{\_}j(\theta)=X$ for some $\theta\in[-\pi,\pi]^{d}$. By letting $\theta$ vary in this range we can trace out a region in $\mathbb{R}^{d}$ on which the amplitudes decay polynomially. In this way, $\nabla H_{\_}j$ maps the region of integration $[-\pi,\pi]^{d}$ into the spatial domain on which the quantum walk resides. This representation requires that we can analytically solve for the eigenvalues $\lambda_{\_}m(\theta)$ of $M(\theta)$ which is not always possible, especially for quantum walks with many states. Moreover, this parametric representation is often complicated and it is not immediately apparent how to connect this representation to the more general structure of the region of polynomial decay, or even to a mathematical description of its boundaries.

We mention that a critical condition for the stationary phase propositions is the smoothness imposed on $f$ and $g$. In the context of proposition 2.1, a discontinuity in an $n^{\text{th}}$ derivative of one of these functions will lead to a slower polynomial rate of decay than the superpolynomial decay guaranteed by the proposition for smooth amplitude and phase functions. It is not trivial to show that the eigenvalues and eigenvectors of the multiplier matrix are smooth. The following proposition was proved by Rainer [20]:

Proposition 2.3

Let $A(t)=[A_{\_}{ij}(t)]_{\_}{1\leq i,j\leq n}$ be a $\mathcal{C}$-curve of normal complex matrices, i.e., the entries $A_{\_}{ij}$ belong to $\mathcal{C}(\mathbb{R},\mathbb{C})$, such that $P_{\_}A$ is normally nonflat. Then there exists a global $\mathcal{C}$-parameterization of the eigenvalues and the eigenprojections of $A$.

Here, $\mathcal{C}$ is a subalgebra of $C^{\infty}$ (Rainer notes that it may be true that $\mathcal{C}=C^{\infty}$), which the entries of $M(\theta)$ can be shown to belong to. However, this proposition only applies to families of normal matrices whose entries are smothly parameterized by a single variable. Families of normal matrices parameterized by more than one variable may not admit a smooth selection of eigenvalues (e.x. consider the family of Hermitian matrices $\begin{bmatrix}x&y\\ y&-x\end{bmatrix}$). This means that for one-dimensional quantum walks, we can provably distinguish between a region of polynomial decay and a region of superpolynomial decay. This distinction is not guaranteed for higher dimensional quantum walks, although in the examples that follow a smooth selection of eigenvalues and eigenvectors can be demonstrated. In any case, for general higher dimensional quantum walks a lack of smoothness in the eigenvalues and eigenvectors of the multiplier matrix would not negate existence of a distinct region of polynomial decay bounded by algebraic curves, it simply does not guarantee qualitatively faster rates of decay outside this region.

3 Multivariate Polynomial System

Instead of focusing our attention on the locus of critical points, we illustrate a method to construct algebraic surfaces representing the bifurcation curves of the quantum walk. These bifurcation curves will more appropriately describe the structure of the region of polynomial decay than does the aforementioned critical point parameterization. The degenerate critical points $\theta$ of the integral in equation (6) must simultaneously satisfy:

[TABLE]

[TABLE]

Using implicit differentiation, one can represent equations (8) and (9) as multivariate polynomial equations in terms of derivatives of the characteristic polynomial. Let $p_{\_}0(\lambda;\theta_{\_}k)$ be the characteristic polynomial of $M(\theta)$ in $\lambda$. Since the entries of $M(\theta)$ are linear combinations of terms of the form $e^{i\theta_{\_}k}$, we can write a multivariate polynomial $p(\lambda;x_{\_}k)$ such that the roots in $\lambda$ of $p(\lambda;e^{i\theta_{\_}k})$ coincide with the eigenvalues of $M(\theta)$ (i.e. the roots of $p_{\_}0(\lambda;\theta_{\_}k)$ in $\lambda$). Thus, if $\lambda(\theta)$ is an eigenvalue of $M(\theta)$, we have:

[TABLE]

Recall that $\lambda(\theta)=e^{iH(\theta)}$, and that we are searching for derivatives of $H(\theta)$ to use in equations (8) and (9). As such, let us take a derivative of this equation with respect to $\theta_{\_}j$:

[TABLE]

From here, let us refer to partial derivatives via subscripts (not to be confused with the indexing of eigenvalues in the previous section) and suppress mention of $\theta$. Rearranging terms, we find:

[TABLE]

If $X=(X_{\_}1,...,X_{\_}n)$, we can substitute equation (12) into equation (9) to find:

[TABLE]

If we take a second derivative of equation (11) with respect to $\theta_{\_}k$, we have:

[TABLE]

By isolating $H_{\_}{jk}$ and substituting the expression in equation (12) for $H_{\_}j$, we may write:

[TABLE]

We can similarly substitute equation (15) into equation (8) such that:

[TABLE]

Both equations (13) and (16) may be used to describe the bifurcation curves, however these equations are dependent on derivatives of $\lambda(\theta)$. We can solve for these derivatives in terms of the characteristic polynomial by taking derivatives of equation (10) with respect to $\theta_{\_}j$. Let us take a first derivative with respect to $\theta_{\_}j$:

[TABLE]

Noticing that $\frac{\partial x_{\_}j}{\partial\theta_{\_}j}=ix_{\_}j$, we can rearrange terms to write:

[TABLE]

This can be substituted into equation (13) to find:

[TABLE]

We take an additional derivative of equation (17) with respect to $\theta_{\_}k$. If $j\neq k$, then we can write:

[TABLE]

If we rearrange terms and substitute equation (18), we have:

[TABLE]

Meanwhile, if $j=k$ then we must account for an additional term:

[TABLE]

Expanding this expression and making similar substitutions, we find:

[TABLE]

We can substitute equations (21) and (23) into equation (16) and eliminate the $p_{\_}\lambda$ denominator factors to arrive at a polynomial equation in $\lambda$ and $\{x_{\_}k\}$, which we term the exponential Hessian determinant.

The equations (10), (16), and (19) make up a system of $n+2$ multivariate polynomial equations in $2n+1$ variables; these are $\lambda$, $\{x_{\_}k\}$, and $\{X_{\_}k\}$. Using a Gröbner basis calculation [7], we can reduce this system to a single equation of $n$ spatial variables $\{X_{\_}k\}$. Unfortunately the exponential Hessian determinant is often prohibitively large and the system requires significant computational resources to solve. However, we present a more feasible naïve method of bifurcation curve computation which, while not rigorously supported, generates curves that bear striking visual resemblance to the quantum walk boundaries. Consider the polynomial system $f(x)=\sum_{\_}{k=0}^{n}a_{\_}kx^{k}=0$ and $g(x)=\sum_{\_}{k=0}^{n}b_{\_}kx^{k}=0$. We wish to find a resultant multivariate polynomial $F(a_{\_}k,b_{\_}k)$ such that selections of coefficients $\{a_{\_}k,b_{\_}k\}$ which admit simultaneous solutions of the polynomial system in $x$ also satisfy the equation $F(a_{\_}k,b_{\_}k)=0$. Such a resultant may be computed using the determinant of a $2n\times 2n$ Sylvester matrix [23]:

[TABLE]

In this notation, $x$ is the variable being cancelled. There are two shortcomings with this formula. First, this resultant will often overrepresent solutions in the system in the sense that solutions of $F(a_{\_}k,b_{\_}k)=0$ in $\{a_{\_}k,b_{\_}k\}$ may not admit solutions in the corresponding polynomial system. For example, consider the system $f(x)=ax+b=0$ and $g(x)=cx+d=0$ such that the resultant satisfies $F(a,b,c,d)=ad-bc=0$. If we let $a=c=0$, then the resultant equation is trivially satisfied, but any nonzero choice of $b$ or $d$ leads to a polynomial system with no solutions. The second shortcoming of this procedure is that the resultant will not take into account any a priori restrictions on the cancelled variable. For instance in our current quantum walk example, we will require that the cancelled variables satisfy $|\lambda|=|x_{\_}k|=1$. It is often difficult to discern which portions of the generated bifurcation curves satisfy these conditions, as these variables are absent from the resultant. Though this method overrepresents solutions of the polynomial system, it will not miss any of the solutions and we use our non-rigorous judgment to hypothesize which ones truly exist in the system.

We use a simple extension of this method to solve a larger system of polynomial equations. Suppose we have a system of $n$ polynomial equations in $n$ variables with a set of variable coefficients. Let us write these polynomials as $p_{\_}{0,k}(x_{\_}1,...,x_{\_}n)$ where $1\leq k\leq n$. Using the Sylvester matrix determinant, let $p_{\_}{1,k}(x_{\_}1,...,x_{\_}{n-1})=\text{Res}(p_{\_}{0,k},p_{\_}{0,n};x_{\_}n)$ for $1\leq k\leq n-1$, in other words we choose a base polynomial $p_{\_}{0,n}$ and compute resultants with the remaining polynomials in the system by cancelling $x_{\_}n$. The new system $\{p_{\_}{1,k}\}$ has $n-1$ equations and $n-1$ variables. We can continue inductively to arrive at a total multivariate resultant polynomial in the coefficient variables. As with the Gröbner basis calculation, even small systems in multiple variables can lead to extremely large resultant polynomials. To combat this, we factor the intermediate polynomials in the system to create a tree of possible solutions, these putative solutions being more feasible to derive. Also, depending on the structure of the system, different choices of base polynomials as well as different orders of variable cancellations can often lead to significant changes in computation time.

Up to this point, the naïve method we have described for computing bifurcation curves is legitimate, we need only take care to ensure that we judiciously select correct bifurcation curves from the overrepresentation provided. However, the exponential Hessian determinant still provides a massive roadblock and renders this naïve method just as intractible as the Gröbner basis method. It has been observed that replacing the exponential Hessian determinant with the far simpler equation $p_{\_}\lambda=0$ results in bifurcation curve solutions which visually bound the regions of polynomial decay, though it is not clear how these curves can be rigorously justified. If we let $p_{\_}\lambda=0$, then there are no restrictions on the spatial variables $X_{\_}j$ in equation (20), and these are the variables of interest. A Gröbner basis calculation would fail to render spatial bifurcation curves in this case while the naïve method generates outputs. We have found that letting $p_{\_}\lambda=0$ be the first base polynomial, and cancelling the variables $\{x_{\_}k\}$ before cancelling $\lambda$ leads to more digestible factors for the algorithm. In the subsequent section, we will clearly state when the displayed bifurcation curves result from the Gröbner basis calculation and when they are derived from the naïve method.

4 Examples

In this section, we compute parameterizations of regions of polynomial decay for five different two-dimensional quantum walks and compute bifurcation curves where possible. In the following cases, characteristic polynomials of the multiplier matrix will often be symmetric quartic polynomials. Suppose we have a characteristic polynomial

[TABLE]

such that $x$ and $y$ are functions of a vector $\theta$. This polynomial can be factored as:

[TABLE]

where $a+b=x$ and $ab+2=y$, otherwise $a(\theta)=\frac{1}{2}\left[x(\theta)\pm\sqrt{x(\theta)^{2}-4y(\theta)+8}\right]$ and $b(\theta)$ takes the opposite sign. These quadratic factors allow for an explicit representation of the eigenvalues, but recall that $\lambda=e^{iH}$ and we are searching for $\nabla H$. Using implicit differentiation and previous equations, we find that the following holds:

[TABLE]

This formula will grant us a parametric representation for the region of polynomial decay in these examples.

4.1 Four-State Grover Walk

We first consider the two-dimensional four-state Grover walk $(\mathbb{Z}^{2},C_{\_}2,G_{\_}4)$. Recall that the $4\times 4$ Grover matrix is written as:

[TABLE]

We write the corresponding multiplier matrix $M(\theta_{\_}1,\theta_{\_}2)$ as:

[TABLE]

The characteristic polynomial of the multiplier matrix thus satisfies:

[TABLE]

We pause to note that the solution $\lambda=\pm 1$ causes a breakdown in the stationary phase approximation in that the only critical point that exists corresponding to this eigenvalue is at $(X_{\_}1,X_{\_}2)=(0,0)$ and is degenerate. This phenomenon is known as localization and has been explored by several authors [11] [12] [15] [21]; we will not elaborate any further on it in this paper. We will term constant solutions to the characteristic polynomial as trivial.

As the non-trivial portion of the characteristic polynomial is quadratic, we can explicitly solve for the non-trivial eigenvalues and use a reduced version of equation (24) to construct a parameterization of the region of polynomial decay:

[TABLE]

It was shown in Kuklinski [17] that this parameterization traces the circle $2X_{\_}1^{2}+2X_{\_}2^{2}\leq 1$ in $\mathbb{R}^{2}$, albeit in an atypical way.

In this example, we can in fact solve for the bifurcation curve. By letting $x_{\_}1=e^{i\theta_{\_}1}$ and $x_{\_}2=e^{i\theta_{\_}2}$, we can write the characteristic equation in a different way:

[TABLE]

The corresponding exponential Hessian determinant is small enough that the system may be efficiently reduced via Gröbner basis computation. The result is as we expect:

[TABLE]

4.2 Five-State Grover Walk

We explore a variant of the four-state Grover walk [4] with the operator $Q\leftrightarrow(\mathbb{Z}^{2},\tilde{C_{\_}2},G_{\_}5)$ The $5\times 5$ Grover matrix is written as:

[TABLE]

We write the corresponding multiplier matrix $M(\theta_{\_}1,\theta_{\_}2)$ as:

[TABLE]

The characteristic polynomial of this multiplier matrix satisfies

[TABLE]

where $c_{\_}1=\cos\theta_{\_}1$ and $c_{\_}2=\cos\theta_{\_}2$. Notice that the eigenvalue $\lambda=1$ leads to localization in this walk as well.

Using equation (24), we can write an explicit parameterization of the region of polynomial decay:

[TABLE]

Here, $s_{\_}1=\sin\theta_{\_}1$, $s_{\_}2=\sin\theta_{\_}2$, and $a=\frac{1}{5}\left(3c_{\_}1+3c_{\_}2+4\pm\sqrt{9c_{\_}1^{2}+9c_{\_}2^{2}-2c_{\_}1c_{\_}2-16c_{\_}1-16c_{\_}2+16}\right)$. Notice that in this case, the choice of plus/minus in $a(\theta_{\_}1,\theta_{\_}2)$ results in different regions of polynomial decay.

To find the bifurcation curves of this quantum walk, we rewrite the characteristic polynomial:

[TABLE]

In this example, the exponential Hessian determinant is several pages long, so a Gröbner basis calculation is computationally infeasable. We choose to illustrate the naïve method for this system. Let $P_{\_}1=p$, $P_{\_}2=p_{\_}1x_{\_}1+\lambda p_{\_}\lambda X_{\_}1$, $P_{\_}3=p_{\_}2x_{\_}2+\lambda p_{\_}\lambda X_{\_}2$, and $P_{\_}4=p_{\_}\lambda$. We first cancel $x_{\_}1$ from this system by letting $P_{\_}4$ be the base polynomial:

[TABLE]

[TABLE]

[TABLE]

The $Q_{\_}k$ polynomials are irreducible. We can omit factors of $\lambda$ in the second two equations as $|\lambda|=1$. If we substitute the factor $\lambda=\pm 1$ from one of these equations, substitute into the remaining two, and then cancel $x_{\_}2$ from the resulting system, the generated curves do not visually fit the system, so we ignore this option and proceed with solving the system $\{Q_{\_}1,Q_{\_}2,Q_{\_}3\}$. By choosing $Q_{\_}2$ to be the base polynomial and eliminating $x_{\_}2$ from the system, we find:

[TABLE]

[TABLE]

If we substitute $\lambda=-1$ from the second equation into $R_{\_}1$, then we find:

[TABLE]

The remaining factors from the second polynomial do not satisfy $|\lambda|=1$ so we ignore these. Although this is not proof, equation (28) visually replicates the boundary of both regions of polynomial decay. Notice that although we have found two distinct regions of polynomial decay in the parametric representation, this irreducible bifurcation curve traces the boundaries of both regions. We note that the outer curve has a maximum distance of $\sqrt{\frac{3}{5}}$ and a minimum distance of $\frac{1}{\sqrt{2}}$ from the origin, while the inner curve has a maximum distance of $\frac{1}{\sqrt{3}}$ and a minimum distance of $\frac{1}{\sqrt{10}}$ from the origin. These maxima are attained on the cardinal axes, and the minima are attained on the axes $y=\pm x$.

4.3 Triangular Quantum Walk

We now analyze the triangular Grover walk with operator $Q\leftrightarrow(\mathbb{Z}^{2},\Sigma,G_{\_}3)$ where

$\Sigma=\{(0,1),(-1,-1),(1,-1)\}$. The multiplier matrix of this operator may be written as:

[TABLE]

Letting $x_{\_}k=e^{i\theta_{\_}k}$, the eigenvalues of this matrix satisfy the following equation:

[TABLE]

The structure of this polynomial does not lend itself to an analytic solution of $\lambda(\theta)$ without invoking the cubic formula [10]. We do not write the formula here, but we may graphically display this parameterization.

Furthermore, the Hessian determinant is too large to facilitate a Gröbner basis calculation for the bifurcation curves. However, using the naïve elimination procedure outputs the bifurcation curve:

[TABLE]

This equation represents an ellipse centered at $(X_{\_}1,X_{\_}2)=(0,-\frac{1}{3})$ with vertical major axis length $\frac{4}{3}$ and horizontal minor axis length $\frac{2}{\sqrt{3}}$.

4.4 Hexagonal Quantum Walk

We now illustrate an example of a quantum walk which traverses a hexagonal lattice on $\mathbb{Z}^{2}$. The hexagonal or honeycomb lattice has been the subject of a few quantum walk investigations [13] [18]. Let us define the set:

[TABLE]

Consider the quantum walk operator $Q\leftrightarrow(\mathbb{Z}^{2},\Sigma,X\otimes G_{\_}3)$ where $X$ is the Pauli-$X$ gate [9] $X=\begin{bmatrix}0&1\\ 1&0\end{bmatrix}$. In this case, amplitudes will travel on two separate hexagonal lattices. To be clear, this hexagonal Grover walk takes place on a subset of $\mathbb{Z}^{2}$ spanned by the elements of $\Sigma$. The multiplier matrix of this quantum walk takes the form:

[TABLE]

The characteristic polynomial of the multiplier matrix of this quantum walk operator is written as:

[TABLE]

The factor $\lambda^{2}-1$ indicates that localization is present in this walk. The remaining eigenvalues satisfy the following equation:

[TABLE]

Since the characteristic polynomial is quadratic in $\lambda^{2}$, we can efficiently parameterize the region of polynomial decay:

[TABLE]

Again, the characteristic polynomial is too large for a Gröbner basis calculation, but the naïve method leads to two possible bifurcation curves:

[TABLE]

These equations represent two concentric ellipses with major axis of length $\sqrt{3}$ times the length of the minor axis. It seems likely that the larger ellipse in equation (31) is a bifurcation curve of the system, but it is uncertain whether the smaller ellipse in equation (32) is a legitimate bifurcation curve.

4.5 Four-State Hadamard Walk

In the final example we consider a quantum walk governed by a different unitary matrix. Let $Q\leftrightarrow(\mathbb{Z}^{2},C_{\_}2,H\otimes H)$ such that the corresponding multiplier matrix becomes:

[TABLE]

The characteristic polynomial of this matrix becomes:

[TABLE]

This is not a symmetric quartic polynomial, but the parametrization of the region of polynomial decay may stille be solved using a modification of equation (24):

[TABLE]

Here, we let $s_{\_}{1+2}=\sin(\theta_{\_}1+\theta_{\_}2)$, $c_{\_}{1+2}=\cos(\theta_{\_}1+\theta_{\_}2)$, and $a=-\frac{1}{2}\left[c_{\_}1+c_{\_}2\pm\sqrt{(c_{\_}1+c_{\_}2)^{2}-4c_{\_}{1+2}+4}\right]$.

By letting $x_{\_}1=e^{i\theta_{\_}1}$ and $e^{i\theta_{\_}2}$, we can rewrite the characteristic polynomial:

[TABLE]

This characteristic polynomial is too large for a Gröbner basis calculation, and even the naïve method cannot generate a complete set of outputs. However, this algorithm is capable of generating the equations for the two main ellipses in the region of polynomial decay:

[TABLE]

The major axes of these ellipses are parallel to the line $x=y$ and have length 1 while the minor axes have length $\frac{\sqrt{2}}{2}$. Though it did not appear in the calculation, we also predict that the bifurcation curve set also includes a rhombus and a $16^{\text{th}}$ order algebraic curve.

5 Conclusion

In this paper we have detailed a process to compute bifurcation curves of two-dimensional quantum walks, as well as describe a non-rigorous algorithm to trace bifurcation curves for more complicated examples which are difficult to solve analytically. In addition, we have provided parameterizations of the regions of polynomial decay. These methods are not unique to two-dimensional quantum walks, and with sufficient computational resources could potentially be extended to computing bifurcation surfaces for higher-dimensional quantum walks.