Maze Solving by a Quantum Walk with Sinks and Self-Loops: Numerical Analysis

: Maze-solving by natural phenomena is a symbolic result of the autonomous optimization induced by a natural system. We present a method for ﬁnding the shortest path on a maze consisting of a bipartite graph using a discrete-time quantum walk, which is a toy model of many kinds of quantum systems. By evolving the amplitude distribution according to the quantum walk on a kind of network with sinks, which is the exit of the amplitude, the amplitude distribution remains eternally on the paths between two self-loops indicating the start and the goal of the maze. We performed a numerical analysis of some simple cases and found that the shortest paths were detected by the chain of the maximum trapped densities in most cases of bipartite graphs. The counterintuitive dependence of the convergence steps on the size of the structure of the network was observed in some cases, implying that the asymmetry of the network accelerates or decelerates the convergence process. The relation between the amplitude remaining and distance of the path is also discussed brieﬂy.


Introduction
Maze-solving methods are important because they have practical applications and provide insight into the invisible intelligence that underlies them. Maze-solving problems can be regarded as a subset of the shortest path problem [1], which is a practical problem in daily life. To solve the maze problem, a maze can be expressed as a network and then solved by an algorithm, such as the depth-first search or the breadth-first search algorithm [2]. There are also maze-solving methods that exploit natural phenomena.
Such methods have been studied experimentally using the Belousov-Zhabotinsky reaction mixtures [3], amoeboid organisms [4], gas discharge [5], and photons in a waveguide array [6]. In these experiments, the result of maze-solving has a symbolic aspect in that it represents the autonomous optimization of the natural system. In this way, the pursuit and modeling of the optimization process in maze solving by a natural phenomenon, can provide a path to a deeper understanding of that phenomenon. The quantum walk model, which has been studied as a quantum counterpart of random walk, has been applied to describe various transportation phenomena in nature [7]. It was first studied as the time-evolution of probability distribution, mainly on a onedimensional network. In the discrete-time quantum walk (DTQW) model, each node has a state vector of complex amplitudes whose dimension corresponds to the number of neighboring nodes. Each evolution is composed of a coin operation and a shift operation; after multiplying the unitary matrix (coin operation), the complex amplitude is transferred into an element of the state vector of a neighboring node (shift operation). By considering time-dependent or site-dependent unitary matrices, the quantum walk can express many kinds of transport dynamics.
The study of quantum walks was extended to arbitrarily connected networks from an early stage [8] because the quantum search on graphs by quantum walks was proposed [9][10][11] as an alternative to Grover's search algorithm [12]. When dealing with discrete-time quantum walks on an arbitrarily connected network, the concept of scattering quantum walks (SQWs) can simplify the model [13].
In an SQW, the state vectors are placed on the edges rather than the nodes. The dimension of each state vector is two: this corresponds to the two directions of an edge between two nodes. Moreover, each node has a scattering matrix that corresponds to the unitary matrix in the coin operation. The time evolution is composed of an intrusion in the node, the multiplication of the scattering matrix, and an escape from the node. The dynamics of an SQW are equivalent to those of a DTQW except for the location of the state vectors.
Recently, the concepts of consecutive injection and corresponding emission into and from the system were incorporated into quantum walks on arbitrarily connected networks [14,15]. For quantum walks on a network with entrances and/or exits, the steadystate [14], trapped-state [15], analogy to an electrical circuit [16], and relationship to the dressed photon phenomenon [17] have been discussed. In particular, the emergence of a trapped state between two self-loops on a network with an exit sink [15] directly motivated the present study, which applies this concept to maze-solving.
Maze-solving using quantum walks has been studied by Hillery, Koch, and Reitzner on an N-tree maze [18] and a chain of stars [19,20]. Their works are the extension of their studies on quantum search and finding structural anomalies in networks [21][22][23][24][25][26][27]. They characterized the start and goal nodes in the maze by reflection with phase inversion, which can be regarded as a pair of structural anomalies.
In this paper, we numerically examine a maze-solving method that uses a quantum walk on a network. The presented method is an application of the emergence of a trapped eigenstate on a network with sinks, and it provides an alternative to previously reported methods [18][19][20]. Although the mathematical foundation of this method was given by Konno, Segawa, and Štefaňák [28], the results presented here are non-trivial because the interaction among multiple trapped eigenstates and the initial condition is generally difficult to characterize as of now.
We show the effectiveness of the method for some examples of the maze with and without cycles and also show the undesirable cases for which this method does not work. The dependence of the number of steps for convergence (convergence steps) on the size of the network structure was also investigated and found to be counterintuitive in certain cases. We also make a tentative discussion about the amount of amplitude remaining on a path and its relative amount among the multiple paths from the numerical results.

Model and Method
In this study, the maze is composed of nodes and edges that connect pairs of nodes. The number of nodes is finite, but pairs of nodes can be connected arbitrarily without limit. The distance between two nodes is given by the smallest number of edges connecting them. Therefore, only distances expressed in positive integers are considered. The start and goal can be placed at any node in the network, even at nodes that are not dead ends. To run the quantum walk, scattering matrices and state vectors are placed on nodes and edges, respectively, as in previous studies on SQWs [13].
The state vectors consist of two complex amplitudes, which express the two directions of the quantum walkers on the edges. As in quantum mechanics, the density of the walkers on an edge is given by the square of the complex amplitude. At each evolution of time, the vector of the incoming component is multiplied by the scattering matrix, generating the vector of the outgoing component. The scattering matrix of the d-dimensional unitary matrix is placed at each node, where d is the number of edges connected to that node. Specifically, we use the scattering matrix of the Grover walk, which, in concrete form, is given by  where a i is the incoming complex amplitude from the i-th edge, and b i is the outgoing complex amplitude to the i-th edge. An example with d = 3 is given in Figure 1a.
To implement maze-solving, two self-loops and a sink are introduced, and the conceptual diagram is shown in Figure 1b. Self-loops are the same as edges except that they are only attached to a single node. As a result, a self-loop has a one-dimensional state vector, where the outgoing amplitude from the node becomes the next incoming amplitude without being modified. For this method, one self-loop is attached to the start node, and the other is attached to the goal node to which a sink node is also attached. The sink node has only one edge, which is connected to the goal node, and its scattering matrix is a zero matrix. The sink serves as the exit from the network for complex amplitudes. The initial amplitude "1" is placed at the self-loop of the start node. To discover the correct path, the initial amplitude should be placed on the path between the start and goal, and it should be kept at a distance from the sink node. In this method, placing the localized amplitude at the self-loop of the start node is the best initial condition for solving the maze correctly without requiring any prior knowledge of the structure of the maze. Finally, note that all the amplitudes in the system are denoted by the real numbers even though the quantum walks are defined using complex amplitudes.
Maze solving was studied for simple structures only because of the large amount of computational time involved by the current code (by the current code on our standard personal computer, the calculation of 10 5 steps took several hours because of unoptimized function-call overhead). For the maze without a cycle, a tree-like structure and a single line with branches were investigated. For the maze with a cycle, independent multiple paths and a ladder-like structure were investigated.
Two undesirable cases, namely, a maze with odd cycles and a maze showing eternal vibration are also investigated. For each kind of structure, the dependence of convergence steps on the size of the structure was investigated. The convergence was judged by the stability of the second decimal place for all the amplitudes in the network, and the convergence steps were expressed with an accuracy of one (or two) significant digits.
For discussion regarding the amount of amplitude remaining, about five digits after the decimal point were considered. The numerical error estimated from the squared sum of the amplitudes was of nearly the same order as the double-precision real number error computed using code written in Python. The source code is available in the repository [29].

Tree-like Structure
We first examine maze-solving for the tree-like structures. This structure has no cycles, and there is only one path from the start to the goal. Figure 2a-c show the results of the amplitude distribution and the number of steps after convergence for the tree-like structures of 2 N leaves for N = 1, 2, and 3. From the results, we observed that only the shortest path emerges as a chain of the eternally remaining densities, whereas the densities on the dead ends vanish during the evolution. The number of convergence steps seems to increase by digits according to the increase of N in these cases. The case of N = 4 was also examined. However, the distribution did not converge even after 10 6 steps that took three days. Figure 2d shows the time profiles of the densities on selected edges, where the label 0-3, for example, denotes the edge between nodes 0 and 3. The densities fluctuate strongly at first and then converge to zero or to positive values. The speed of convergence varies according to the position of the edge; the greater the distance to the sink node is, the slower the speed of convergence.
To consider the influence of the extra branching at dead ends on the convergence steps, the cases of decreased and increased extra branching based on Figure 2c were examined. For the case with decreased branching as shown in Figure 2e, the convergence steps decreased, which was an intuitive result. However, the decrease in convergence steps was more for the case with increased extra branching as shown in Figure 2f. This counterintuitive dependence is difficult to explain for the present. However, it can be suggested that the extent of asymmetry in the network accelerated the convergence.
For the cases of Figure 2c,e,f, the absolute values of the converged amplitudes on the correct paths, including self-loops, were all 0.08. That value seems to have been determined by the distance between the start and the goal nodes for the case of the network without cycles. Table 1 lists the relation between the distance between the start and the goal and the absolute value of the amplitude remaining on an edge for each case. Edges indicates the number of edges on a path including self-loops. (Edges = 2 × Distance + 2).
A rational expression approximating the amount of amplitude was attached for each case. For these cases, the amplitudes can be expressed by the inverse of the number of edges included in the path. Namely, the sum of amplitudes along the path is "1.0" for all the cases. However, note that the "1.0" does not indicate all the amplitudes injected into the system because that is not the square sum of the amplitudes.

A Line with Branches
To investigate the dependence of the convergence steps on the placement of the branches, the maze-solving for various patterns of a line with shallow dead ends was examined. Figure 3a shows the result for a simple line constructed based on the correct path of Figure 2c. The number of convergence steps decreased by two orders of magnitude from the case shown in Figure 2c. Figure 3b shows the result for a line with four shallow dead-ends. Nearly the same result as Figure 2e was obtained, as the difference between them was only the length of the dead ends. Figure 3c-e shows the results for patterns of placement of one to three dead ends, respectively. The distribution of the amplitudes is omitted, but ±0.08, which is the same as Figure 2c,e,f, is on the correct path, and 0.00 is on the dead-end edges in all cases. The convergence steps varied not only by the number of dead ends but also by the positions. As the trend shows, the convergence steps became larger with increasing dead ends; however, exceptions were observed depending on the positions of the dead ends. The convergence steps became larger for the case where dead ends were attached closed to the goal node. This seems counterintuitive considering the quick convergence near the sink, which was observed in Figure 2d.
The number of convergence steps for the case of Figure 3e was much larger than for Figure 3b or Figure 2c. For these cases, the asymmetry significantly decelerated the convergence speed, which is in contrast to the acceleration due to asymmetry observed in Figure 2e,f. A maze without cycles can be solved by this method; however, the dependence of the convergence steps on the network structure is difficult to predict intuitively.

Independent Multiple Paths
Next, we examined maze-solving on multiple independent paths of different lengths. This structure includes cycles, which makes maze-solving difficult even in classical schemes. Figure 4a-c,e,f shows the numerical results for the networks with M paths, where the length of the Mth path is 2M. After convergence, in all the examples shown here, the densities remain on all the paths between the start and the goal; however, the maximum densities are only observed on the shortest path, while smaller densities are observed farther from the shortest path. By regarding the path of the maximum densities as the correct path, the maze-solving was successful for these examples. Figure 4d shows the time profiles of the edges on each path. The speed of convergence was higher than in the case of other structures of a similar scale. The reason for this is unclear, but a lack of branching on the paths may be responsible for the high speed. Among the three paths, the speeds of convergence did not differ significantly, and they did not depend on the distance from the sink unlike in the tree-like structure. In general, the convergence steps increased by the addition of other paths. However, a counter-intuitive decrease of the convergence steps was observed in Figure 4e,f.
The absolute values of amplitudes, after convergence, decrease as the length of the path becomes longer; however, they are not constant for the length of paths because a slight decrease was observed by additional paths. Table 2 lists the relation between the distances of paths and amplitude remaining on an edge. For Figure 4a, the amplitude is the inverse of the number of edges included in the path, which is the same as given in Table 1. However, the rule looks broken in the case of multiple paths. Table 2. The relation between the distance and remaining amplitude for Figure 4 (The waypoint of a path, the distance between the start and goal on the path, the number of edges in the path, the amplitude remaining on an edge on the path, and an approximate rational expression of the amplitude. Only the relative ratios are shown for Figure 4f because an appropriate rational number was not found).    Step 4-7 5-8 6-9 Most of the amplitudes were assigned rational expressions; however, the rule determining the absolute value (or a positive integer of the denominator) is not clear. However, the relative amounts of amplitudes among paths in each case were found to be in inverse proportion to the distance between the start and goal exactly for these cases. The relative amounts of amplitude were determined not by the number of edges but by the distances.

Ladder-like Structure
As in other small examples of mazes with cycles, the ladder-like structures with L paths were examined. The difference from the previous subsection is that the edges are shared among the different paths. Figure 5a,b,d shows the results of L = 1, 3, and 4, respectively. For the cases shown here, the shortest paths are indicated by the chain of maximum densities, while smaller densities are observed farther from the shortest path, meaning that the maze-solving was successful. The absolute values of amplitude after convergence seem to correspond to the distance of each path by considering Figure 5b,d. However, this is not the case in Figure 5a.
For the cases of L = 2 and 5, undesirable eternal vibrations were observed and mazesolving could not work, which will be discussed in a later subsection. For the case of L = 6, the convergence was difficult to realize owing to the limitation of the computational times. However, it was not an eternal vibration judging from the actual calculations. Figure 5c shows the time profiles of the edges in each path of Figure 5b. For the three paths, the speeds of convergence did not differ significantly, and they did not depend on the distance from the sink unlike what was observed in the tree-like structure. The convergence is faster than for the tree-like structure but slower than for the independent multiple paths. The number of convergence steps in Figure 5d is smaller than that in Figure 5b, exhibiting the difficulty faced in predicting the convergence speed from the structure of the maze. Table 3 lists the relation between the distances of paths and amplitude remaining on an edge in Figure 5. The amplitudes can be expressed by rational numbers; however, the meanings of the denominator numbers are not clear as in Figure 4. The absolute values of the amplitudes in Figure 5 are generally smaller than those of the same scale case in Figure 4.
The ratio among the paths seems to have meaning; however, the reason has not been determined except for the relation between the longest and second-longest paths. The relative ratio of the longest path and the second-longest path is thought to be in inverse proportion to the ratio of the length of the non-shared part of each path. This hypothesis is complemented in the next subsection. Table 3. The relation between the distance and remaining amplitude for Figure 5 (The waypoint of a path, the distance between the start and goal on the path, the number of edges in the path, the amplitude remaining on an edge on the path, and an approximate rational expression of the amplitude).

Small Maze
To demonstrate slightly complicated cases, solutions of small mazes by the method are presented. The maze includes some dead ends and two paths to the goal as shown in Figure 6a. As in the other related cases shown above, the maximum density remains on the shortest path, and the densities of the dead-end paths vanish. The maze-solving worked correctly for a small maze with both dead-ends and cycles. Figure 6b shows the time profiles of the densities on selected edges. The convergence speeds were not so different from each other, as in the case of ladder-like structures.
The result of another maze that is slightly modified from Figure 6a is shown in Figure 6c. The convergence step was 22,000 for this example; however, a drastic increase of convergence steps was often observed by another slight modification of the structure. It might be in rare cases that the complex maze could be solved in a permissible computational time. Table 4 lists the relations between the distances of paths and amplitude remaining on an edge in Figure 6. The hypothesis that the relative ratio of the longest path and the second-longest path is in inverse proportion to the ratio of the distances of the non-shared part of each path was also confirmed for these cases.  Table 4. The relation between the distance and remaining amplitude for Figure 6 (The waypoint of a path, the distance between the start and goal on the path, the amplitude remaining on an edge on the path, and an approximate rational expression of the amplitude).

Undesirable Cases 1: Odd Cycle
Here, we show some examples of undesirable cases where maze-solving did not work. First, this method cannot be applicable for a maze that includes odd cycles. Figure 7a shows a network with a single odd cycle whose length is 5. The cycle that consists of nodes 1, 2, 5, 6, and 3 is the odd cycle. When the amplitude distribution converged, the absolute value of the amplitude between the exit of the cycle and the goal (edges between nodes 4 and 5) became small. The correct path was not indicated by the maximum densities, meaning that the maze-solving went wrong.  Figure 7b shows an attempt of solving for the network with two sequential odd-cycles. The effects of the two odd-cycles were not canceled out, and only a small amplitude reached the goal. The solving method presented cannot apply to the network with odd-cycles.

Undesirable Cases 2: Eternal Vibration
Even though the odd cycle was not involved in the network, undesirable eternal vibration was observed in some cases. Figure 8a shows the network of the ladder-like structure for L = 2, which is exhibiting eternal vibration. In this case, only the edges between nodes 5 and 6, were stabilized. The amplitudes of other edges, from the start to the cycle, exhibit a constant vibration pattern eternally. Figure 8b shows the time profiles of the densities of some edges. The constant amplitude vibrations seem to continue eternally and not converge. The inset shows the details of the vibrational behavior. The same patterns are seen to be repeating. The eternal vibration was observed only for the ladderlike structure of L = 2 and 5.
We found that the eternal vibration was suppressed by the addition of an extra dead end. Figure 8c shows the network in which one dead-end is attached to Figure 8a. The amplitude distribution converged, and the shortest path was indicated by the maximum densities. Figure 8d shows the time profile of the density for some selected edges in Figure 8c. The reason for the stabilization is unclear at present; however, a small perturbation of the network may have a significant influence on the behavior of quantum walks.

Discussion
In applying the proposed method to mazes without odd cycles, we verified that the paths between the start and goal emerge as trapped states of the quantum walk, and the density on the shortest path was maximized autonomously. As the network without odd cycles is regarded as a bipartite graph, we concluded that the method can be applied to the bipartite graph except for the case where the eternal vibration emerges. The condition for the occurrence of the eternal vibration is not clear as of now as only a few examples were considered.
The key features of the proposed method are the self-loops at the start and the goal and the sink node attached to the goal. In previous studies, the start and goal were marked by reflection with phase inversion placed at the dead ends [18][19][20]. The correct path was then judged by the transient profile of the probabilities. Our method partially improves upon past works by incorporating self-loops, which can be placed anywhere in the maze, and by determining the correct path according to the eternally remaining densities.
We now consider the remaining densities on the correct path in terms of knowledge that has been proven mathematically. The eigenstate of the time evolution operator of the quantum walk with sinks was constructed on the path between two-self loops [28]. This eigenstate is called the trapped state, and it is not absorbed by the sink. In the Grover walk, the eigenstates are constructed between two self-loops and also around the cycles [28].
To generate a trapped state, the initial amplitude should be placed on the edge that is to be included in the trapped state. This was the reason why the initial amplitude was placed on the self-loop of the start node. If the initial amplitude is placed randomly at the edge, the trapped state on the correct path does not always emerge because the initial edge may not be included on the path between the start and the goal. Even if the initial state is on the correct path, a trapped state may also emerge in the cyclic structure that includes the initial edge. In this case, the shortest path may not have the maximum density. When placing the initial amplitude on a self-loop, the initial edge is not included in any cyclic structure in the network, and only the paths between the start and the goal emerge.
The role of the sink should also be considered. The dynamics of this type of Grover walk can be separated into an electric current component that propagates rapidly, and a random-walk component that propagates slowly [16]. The emergence of the trapped state results from the electric current component; hence, to observe the trapped state, the random-walk component must be eliminated by the sink.
Even without the mathematical knowledge above, the amplitude distribution after convergence can be interpreted by the simple rules observed in the numerical results.
The key rule for determining amplitude distribution is that the sum of incoming/outgoing amplitudes to/from a node must be zero, separately. This rule is mathematically and numerically exact at all the nodes in the examples that converged. Figure 9a shows an example of amplitude distribution around a node on the correct path in the maze.
As the sum of incoming and outgoing amplitudes should be zero separately, amplitudes of plus and minus emerge alternately on the line. Figure 9b shows an example of an amplitude distribution around a dead-end node. As only one amplitude is incoming to the dead-end node, that should be zero to make the sum of the incoming amplitude zero. This is why the amplitudes vanish on the path to the dead-end. Figure 9c shows an example of amplitude distribution around a self-loop. In this case, the amplitude on the self-loop acts as both incoming and outgoing amplitudes to keep the sum zero for both. This is the reason why the signs of the amplitudes are the same on the edges connected to the node involving the self-loop. These facts fit all the nodes included in the numerical results after convergence.
The sum rule above can be used to also explain the amplitude distribution around the even cycle and odd cycle. Figure 9d shows an example of the amplitude distribution around an even cycle. When the large positive amplitude enters the cycle, two small negative amplitudes are generated at the first branching node. Both amplitudes move on the cycle by changing the sign alternately and meet again on the join node. If the cycle is an even cycle, two small negative amplitudes make a large positive amplitude to the outside of the cycle to keep the sum rule. For the case of an odd cycle (Figure 9e), two amplitudes meet at the join node with different signs. To maintain the sum rule, only the smaller amplitude, which is nearly zero, generates the output. This is the reason why the maze, including the odd cycle, cannot be solved by this method.
The maze-solving speed of this method is clearly considerably slower than that of other known algorithms. Although the examples were limited, the convergence steps were difficult to predict by intuition in observing the structure of the network. At present, the intuitive unified parameter that connects the network structure and convergence speed has not been determined mathematically. Further analysis, considering some other aspect, such as the symmetry of the graph, may be required.
The general reason that the maximum densities emerge on the shortest path remains unclear at present; however, some tentative rules were observed numerically. In many cases, the absolute values of amplitude that remained could be approximated to a rational number composed of integers. When there is only one path to the goal, the absolute values of amplitudes on an edge become the inverse number of the number of edges included in the path. This is observed in Figures 2, 3, 4a and 5a. The preserved amount is not the square of the amplitude but the absolute value of amplitude. It is the same as the sum rule discussed above.
When there are multiple paths to the goal and they do not share edges mutually, the ratio of the absolute values of amplitudes is in inverse proportion to the ratio of the distance of the paths. This is observed in Figure 4b-f. When there are two paths to the goal and they share some edges, the ratio of the absolute values of amplitudes is in inverse proportion to the ratio of the distance of the non-shared part of the paths. This is observed in Figure 6a,c. When there are more than two paths to the goal and they share some edges, the dependence of the relative amount of amplitudes on the distance is unclear; however, certain rules clearly exist. This was seen in Figure 5b,d.
To apply the method presented for actual problems of the shortest path finding, the lengths of all the paths should be expressed by positive integers. The odd-loops must not be included; however, the odd-loops would be eliminated by a slight modification of the distance in the process of discretization of the network. The eternal vibration is still the obstacle of the path-finding problem. This should be analyzed more both mathematically and numerically. Additionally, this method cannot involve the negative distance that is considered in some classical algorithms.
Studying this maze-solving method may not appear to be of much use from the viewpoint of computational algorithms; however, it may help to understand the mechanisms of autonomous features that can be observed in a natural system because the quantum walk is a toy model that can be applied to the energy transportation in quantum fields, such as dressed photon phenomena [30].
While the emergence of the shortest path or some other optimized structure in a natural phenomenon may seem mysterious at first glance, they may have an analogy in maze-solving using the quantum walks. Moreover, the implicit existence of the sink node may play an important role in such systems. Funding: This research received no external funding.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.