Bifurcation Curves of Two-Dimensional Quantum Walks

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.

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

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_{\lx@text@underscore}d\subset\mathbb{Z}^{d}$ be the set of unit directional vectors and let $\tilde{C_{\lx@text@underscore}d}=C_{\lx@text@underscore}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_{\lx@text@underscore}n$ is the $n\times n$ identity matrix and $\text{\bf 1}_{\lx@text@underscore}n$ is the $n\times n$ matrix filled with ones, then we define the $n\times n$ Grover matrix as $G_{\lx@text@underscore}n=\frac{2}{n}\text{\bf 1}_{\lx@text@underscore}n-I_{\lx@text@underscore}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_{\lx@text@underscore}\sigma:\ell^{2}(G)\rightarrow\ell^{2}(G)$ be a translation operator which acts as $T_{\lx@text@underscore}\sigma(\psi(g))=\psi(\sigma^{-1}\cdot g)$ with $\sigma\in\Sigma$. Then we can visualize $Q$ acting on the vector $\Psi=\left[\psi_{\lx@text@underscore}{\sigma_{\lx@text@underscore}1},...,\psi_{\lx@text@underscore}{\sigma_{\lx@text@underscore}n}\right]$ as follows:

$$Q\Psi=\begin{bmatrix}T_{\lx@text@underscore}{\sigma_{\lx@text@underscore}1}&~&~\\ ~&\ddots&~\\ ~&~&T_{\lx@text@underscore}{\sigma_{\lx@text@underscore}n}\end{bmatrix}U\Psi$$

(1)

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:

$$\mathcal{F}[T_{\lx@text@underscore}\sigma\psi](\theta)=e^{i\sigma\cdot\theta}\mathcal{F}[\psi](\theta)$$

(2)

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

$$\mathcal{F}[Q\Psi](\theta)=\begin{bmatrix}e^{i\sigma_{\lx@text@underscore}1\cdot\theta}&~&~\\ ~&\ddots&~\\ ~&~&e^{i\sigma_{\lx@text@underscore}n\cdot\theta}\end{bmatrix}U\mathcal{F}[\Psi]=M(\theta)\mathcal{F}[\Psi](\theta)$$

(3)

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_{\lx@text@underscore}m(\theta)\}_{\lx@text@underscore}{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_{\lx@text@underscore}0$ as:

$$M(\theta)^{t}\mathcal{F}[\Psi_{\lx@text@underscore}0](\theta)=\sum_{\lx@text@underscore}{m=1}^{n}\lambda_{\lx@text@underscore}j(\theta)^{t}E_{\lx@text@underscore}m(\theta)$$

(4)

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

$$Q^{t}\Psi_{\lx@text@underscore}0=\frac{1}{(2\pi)^{d}}\sum_{\lx@text@underscore}{m=1}^{n}\int_{\lx@text@underscore}{\lVert\theta\rVert_{\lx@text@underscore}\infty\leq\pi}E_{\lx@text@underscore}m(\theta)\exp\left[i(tH_{\lx@text@underscore}m(\theta)-x\cdot\theta)\right]d\theta$$

(5)

Here, $\lVert\cdot\rVert_{\lx@text@underscore}\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:

$$Q^{t}\Psi_{\lx@text@underscore}0=\frac{1}{(2\pi)^{d}}\sum_{\lx@text@underscore}{m=1}^{n}\int_{\lx@text@underscore}{\lVert\theta\rVert_{\lx@text@underscore}\infty\leq\pi}E_{\lx@text@underscore}m(\theta)\exp\left[it(H_{\lx@text@underscore}m(\theta)-X\cdot\theta)\right]d\theta$$

(6)

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

$$I(t)=\int_{\lx@text@underscore}{\mathbb{R}^{d}}g(x)e^{itf(x)}dx.$$

(7)

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

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

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):

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:

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_{\lx@text@underscore}j(\theta)=H_{\lx@text@underscore}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_{\lx@text@underscore}j(\theta)=0$, or if $\nabla H_{\lx@text@underscore}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_{\lx@text@underscore}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_{\lx@text@underscore}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]:

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.

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:

$$\det\left(\left[\frac{\partial^{2}H_{\lx@text@underscore}m}{\partial\theta_{\lx@text@underscore}j\partial\theta_{\lx@text@underscore}k}\right](\theta)\right)=0$$

(8)

$$\nabla H_{\lx@text@underscore}m(\theta)=X$$

(9)

Using implicit differentiation, one can represent equations (8) and (9) as multivariate polynomial equations in terms of derivatives of the characteristic polynomial. Let $p_{\lx@text@underscore}0(\lambda;\theta_{\lx@text@underscore}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_{\lx@text@underscore}k}$, we can write a multivariate polynomial $p(\lambda;x_{\lx@text@underscore}k)$ such that the roots in $\lambda$ of $p(\lambda;e^{i\theta_{\lx@text@underscore}k})$ coincide with the eigenvalues of $M(\theta)$ (i.e. the roots of $p_{\lx@text@underscore}0(\lambda;\theta_{\lx@text@underscore}k)$ in $\lambda$). Thus, if $\lambda(\theta)$ is an eigenvalue of $M(\theta)$, we have:

$$p\left(\lambda(\theta);x_{\lx@text@underscore}k(\theta)\right)=0$$

(10)

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_{\lx@text@underscore}j$:

$$\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}j}=i\frac{\partial H}{\partial\theta_{\lx@text@underscore}j}e^{iH(\theta)}=i\frac{\partial H}{\partial\theta_{\lx@text@underscore}j}\lambda(\theta)$$

(11)

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:

$$H_{\lx@text@underscore}j=\frac{\lambda_{\lx@text@underscore}j}{i\lambda}$$

(12)

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

$$\lambda_{\lx@text@underscore}j-i\lambda X_{\lx@text@underscore}j=0$$

(13)

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

$$\lambda_{\lx@text@underscore}{jk}=i(H_{\lx@text@underscore}{jk}\lambda+H_{\lx@text@underscore}j\lambda_{\lx@text@underscore}k)$$

(14)

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

$$H_{\lx@text@underscore}{jk}=\frac{i}{\lambda^{2}}\left(\lambda_{\lx@text@underscore}j\lambda_{\lx@text@underscore}k-\lambda_{\lx@text@underscore}{jk}\lambda\right)$$

(15)

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

$$\det\left(\left[\lambda_{\lx@text@underscore}j\lambda_{\lx@text@underscore}k-\lambda_{\lx@text@underscore}{jk}\lambda\right]\right)=0$$

(16)

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_{\lx@text@underscore}j$. Let us take a first derivative with respect to $\theta_{\lx@text@underscore}j$:

$$\frac{\partial p}{\partial\lambda}\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}j}+\frac{\partial p}{\partial x_{\lx@text@underscore}j}\frac{\partial x_{\lx@text@underscore}j}{\partial\theta_{\lx@text@underscore}j}=0$$

(17)

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

$$\lambda_{\lx@text@underscore}j=-\frac{ip_{\lx@text@underscore}jx_{\lx@text@underscore}j}{p_{\lx@text@underscore}\lambda}$$

(18)

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

$$p_{\lx@text@underscore}jx_{\lx@text@underscore}j+\lambda p_{\lx@text@underscore}\lambda X_{\lx@text@underscore}j=0$$

(19)

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

$$\left(\frac{\partial p_{\lx@text@underscore}\lambda}{\partial\lambda}\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}k}+\frac{\partial p_{\lx@text@underscore}\lambda}{\partial x_{\lx@text@underscore}k}\frac{\partial x_{\lx@text@underscore}k}{\partial\theta_{\lx@text@underscore}k}\right)\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}j}+\frac{\partial p}{\partial\lambda}\frac{\partial^{2}\lambda}{\partial\theta_{\lx@text@underscore}j\partial\theta_{\lx@text@underscore}k}+\left(\frac{\partial p_{\lx@text@underscore}j}{\partial\lambda}\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}k}+\frac{\partial p_{\lx@text@underscore}j}{\partial x_{\lx@text@underscore}k}\frac{\partial x_{\lx@text@underscore}k}{\partial\theta_{\lx@text@underscore}k}\right)\frac{\partial x_{\lx@text@underscore}j}{\partial\theta_{\lx@text@underscore}j}=0$$

(20)

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

$$\lambda_{\lx@text@underscore}{jk}=\frac{x_{\lx@text@underscore}jx_{\lx@text@underscore}k}{p_{\lx@text@underscore}\lambda^{3}}\left[p_{\lx@text@underscore}{\lambda\lambda}p_{\lx@text@underscore}jp_{\lx@text@underscore}k-p_{\lx@text@underscore}{\lambda k}p_{\lx@text@underscore}jp_{\lx@text@underscore}\lambda-p_{\lx@text@underscore}{\lambda_{\lx@text@underscore}j}p_{\lx@text@underscore}kp_{\lx@text@underscore}\lambda+p_{\lx@text@underscore}{jk}p_{\lx@text@underscore}\lambda^{2}\right]$$

(21)

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

$$\left(\frac{\partial p_{\lx@text@underscore}\lambda}{\partial\lambda}\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}j}+\frac{\partial p_{\lx@text@underscore}\lambda}{\partial x_{\lx@text@underscore}j}\frac{\partial x_{\lx@text@underscore}j}{\partial\theta_{\lx@text@underscore}j}\right)\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}j}+\frac{\partial p}{\partial\lambda}\frac{\partial^{2}\lambda}{\partial\theta_{\lx@text@underscore}j^{2}}+\left(\frac{\partial p_{\lx@text@underscore}j}{\partial\lambda}\frac{\partial\lambda}{\partial\theta_{\lx@text@underscore}j}+\frac{\partial p_{\lx@text@underscore}j}{\partial x_{\lx@text@underscore}j}\frac{\partial x_{\lx@text@underscore}j}{\partial\theta_{\lx@text@underscore}j}\right)\frac{\partial x_{\lx@text@underscore}j}{\partial\theta_{\lx@text@underscore}j}+\frac{\partial p}{\partial x_{\lx@text@underscore}j}\frac{\partial^{2}x_{\lx@text@underscore}j}{\partial\theta_{\lx@text@underscore}j^{2}}=0$$

(22)

Expanding this expression and making similar substitutions, we find:

$$\lambda_{\lx@text@underscore}{jj}=\frac{x_{\lx@text@underscore}j}{p_{\lx@text@underscore}\lambda^{3}}\left[p_{\lx@text@underscore}{\lambda\lambda}p_{\lx@text@underscore}j^{2}x_{\lx@text@underscore}j-2p_{\lx@text@underscore}{\lambda j}p_{\lx@text@underscore}jp_{\lx@text@underscore}\lambda x_{\lx@text@underscore}j+p_{\lx@text@underscore}{jj}p_{\lx@text@underscore}\lambda^{2}x_{\lx@text@underscore}j+p_{\lx@text@underscore}jp_{\lx@text@underscore}\lambda^{2}\right]$$

(23)

We can substitute equations (21) and (23) into equation (16) and eliminate the $p_{\lx@text@underscore}\lambda$ denominator factors to arrive at a polynomial equation in $\lambda$ and $\{x_{\lx@text@underscore}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_{\lx@text@underscore}k\}$, and $\{X_{\lx@text@underscore}k\}$. Using a Gröbner basis calculation [7], we can reduce this system to a single equation of $n$ spatial variables $\{X_{\lx@text@underscore}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_{\lx@text@underscore}{k=0}^{n}a_{\lx@text@underscore}kx^{k}=0$ and $g(x)=\sum_{\lx@text@underscore}{k=0}^{n}b_{\lx@text@underscore}kx^{k}=0$. We wish to find a resultant multivariate polynomial $F(a_{\lx@text@underscore}k,b_{\lx@text@underscore}k)$ such that selections of coefficients $\{a_{\lx@text@underscore}k,b_{\lx@text@underscore}k\}$ which admit simultaneous solutions of the polynomial system in $x$ also satisfy the equation $F(a_{\lx@text@underscore}k,b_{\lx@text@underscore}k)=0$. Such a resultant may be computed using the determinant of a $2n\times 2n$ Sylvester matrix [23]:

$$\text{Res}(f,g;x)=F(a_{\lx@text@underscore}k,b_{\lx@text@underscore}k)=\det\begin{bmatrix}a_{\lx@text@underscore}n&a_{\lx@text@underscore}{n-1}&\ldots&a_{\lx@text@underscore}0&0&~&~&~\\ b_{\lx@text@underscore}n&b_{\lx@text@underscore}{n-1}&\ldots&b_{\lx@text@underscore}0&0&\ddots&~&~\\ 0&a_{\lx@text@underscore}n&\ldots&a_{\lx@text@underscore}1&a_{\lx@text@underscore}0&\ddots&~&~\\ 0&b_{\lx@text@underscore}n&\ldots&b_{\lx@text@underscore}1&b_{\lx@text@underscore}0&\ddots&~&~\\ ~&~&\ddots&\ddots&\ddots&\ddots&~&~\\ ~&~&~&~&~&a_{\lx@text@underscore}n&\ldots&a_{\lx@text@underscore}0\\ ~&~&~&~&~&b_{\lx@text@underscore}n&\ldots&b_{\lx@text@underscore}0\end{bmatrix}$$

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_{\lx@text@underscore}k,b_{\lx@text@underscore}k)=0$ in $\{a_{\lx@text@underscore}k,b_{\lx@text@underscore}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_{\lx@text@underscore}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_{\lx@text@underscore}{0,k}(x_{\lx@text@underscore}1,...,x_{\lx@text@underscore}n)$ where $1\leq k\leq n$. Using the Sylvester matrix determinant, let $p_{\lx@text@underscore}{1,k}(x_{\lx@text@underscore}1,...,x_{\lx@text@underscore}{n-1})=\text{Res}(p_{\lx@text@underscore}{0,k},p_{\lx@text@underscore}{0,n};x_{\lx@text@underscore}n)$ for $1\leq k\leq n-1$, in other words we choose a base polynomial $p_{\lx@text@underscore}{0,n}$ and compute resultants with the remaining polynomials in the system by cancelling $x_{\lx@text@underscore}n$. The new system $\{p_{\lx@text@underscore}{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_{\lx@text@underscore}\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_{\lx@text@underscore}\lambda=0$, then there are no restrictions on the spatial variables $X_{\lx@text@underscore}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_{\lx@text@underscore}\lambda=0$ be the first base polynomial, and cancelling the variables $\{x_{\lx@text@underscore}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.

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

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

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

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

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:

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

(24)

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