# A Comparison of Dijkstra’s Algorithm Using Fibonacci Heaps, Binary Heaps,
and Self-Balancing Binary Trees

Rhyd Lewis  School of Mathematics,
Cardiff University, Cardiff, Wales.
<EMAIL_ADDRESS>http://www.rhydlewis.eu

###### Abstract

This paper describes the shortest path problem in weighted graphs and examines
the differences in efficiency that occur when using Dijkstra’s algorithm with
a Fibonacci heap, binary heap, and self-balancing binary tree. Using C++
implementations of these algorithm variants, we find that the fastest method
is not always the one that has the lowest asymptotic complexity. Reasons for
this are discussed and backed with empirical evidence.

##  1 Introduction

Dijkstra’s algorithm is an efficient, exact method for finding shortest paths
between vertices in edge- and arc-weighted graphs. It is particularly useful
in transportation problems when we want to determine the shortest (or fastest)
route between two geographic locations on a road network [7, 10, 14]. It is
also applicable in areas such as telecommunication, social network analysis,
arbitrage, and currency exchange [15, 18].

In this paper, we examine the changes in computational complexity and
computing times that occur when using either a self-balancing binary tree,
binary heap, or Fibonacci heap within Dijkstra’s algorithm. As part of this
work, we give an efficient C++ implementation of the algorithm and of
Fibonacci heaps. Many programming languages, including C++, contain versions
of self-balancing binary trees and binary heaps as part of their libraries;
however, implementations of Fibonacci heaps are less common. Existing C++
implementations of Fibonacci heaps are also buggy, inefficient, and/or
difficult to use. This is not the case for the custom implementation used
here, which has been fully tested and evaluated.

The next section formally defines the shortest path problem and surveys
several algorithms for solving it. Section 3 gives a detailed description of
Dijkstra’s algorithm, while Section 4 shows how the efficiency of this method
can be improved through the use of priority queues. In Section 5 we describe
implementations of four variants of Dijkstra’s algorithm, which are then
compared and evaluated in Section 6. Conclusions are drawn in Section 7.

##  2 Problem Definition and Existing Algorithms

Let $G=(V,A)$ be an arc-weighted, directed graph in which $V$ is a set of $n$
vertices, and $A$ is a set of $m$ arcs (directed edges). In addition, let
$\Gamma(u)$ denote the set of vertices that are _neighbours_ of a vertex $u$.
That is, $\Gamma(u)=\\{v:(u,v)\in A\\}$. Finally, we also define a nonnegative
_weight_ (or length) $w(u,v)$ for each arc $(u,v)\in A$. The weight (or
length) of a path is defined by the sum of the weights of its arcs.

According to Cormen et al. [9], three problems involving shortest paths on
arc-weighted graphs can be distinguished:

The single-source single-target shortest path problem.

This involves finding the shortest path between a particular _source_ vertex
$s$ and _target_ vertex $t$. In other words, we want to identify the
$s$-$t$-path in $G$ whose length (weight) is minimal among all possible
$s$-$t$-paths.

The single-source shortest path problem.

This involves determining the shortest path from a source $s$ to all other
reachable vertices in $G$. In this sense, we are seeking a “shortest path tree
rooted at $s$”. An example of such a tree is shown in Figure 1.

The all-pairs shortest path problem.

This involves finding the shortest path between every pair of vertices in $G$.
That is, we are seeking the shortest $u$-$v$-paths for all $u,v\in V$.

Figure 1: An example graph with $n=100$ vertices (the black circles). In this
case, arcs are drawn with straight lines, and arc weights correspond to the
lengths of these lines. This particular graph is also symmetric in that
$(u,v)\in A$ if and only if $(v,u)\in A$, with $w(u,v)=w(v,u)$. The
highlighted arcs show a shortest path tree rooted at the vertex $s$. Because
this graph is a strongly connected component, the shortest path tree is also a
spanning tree.

In this work, we will use Dijkstra’s algorithm to solve the second problem in
the above list. It can, however, also be used for the other two variants. To
do this with the single-source single-target shortest path problem, we simply
need to halt Dijkstra’s algorithm as soon as the target vertex becomes
“distinguished” (see Section 3). For the all-pairs shortest path problem,
meanwhile, it is sufficient to execute Dijkstra’s algorithm $n$ times, using
each vertex $u\in V$ as the source in turn.

Several other algorithms also exist for the above three problems. In cases
where graphs feature negative arc weights, a more suitable alternative is the
$O(nm)$-time Bellman-Ford algorithm [9]. Although this algorithm has a higher
growth rate than that of Dijkstra’s, it has the added advantage of being able
to detect if a particular instance of the shortest path problem is “ill-
defined”, in that it contains negative cycles. Bellman-Ford can also be
augmented with additional data structures to form Moore’s algorithm [16]
which, though still featuring a complexity of $O(nm)$, usually features faster
run times than Bellman-Ford.

For the single-source single-target shortest path problem, other specialised
algorithms exist, though none of these is known to run asymptotically faster
than Dijkstra’s algorithm. One well-known alternative is the A* algorithm of
Hart et al. [11, 12]. This is a heuristic-based variant of Dijkstra’s
algorithm and usually gives much faster run times in applications involving
transportation networks. Algorithms for variants of the single-source single-
target shortest path problem have also been proposed by Yen [19] and Bhandari
[8]. Yen’s algorithm is used to find the $k$ shortest $s$-$t$-paths, where $k$
is a user-defined parameter. The methods of Bhandari, meanwhile, are used to
produce a pair of shortest $s$-$t$-paths that are either edge-disjoint and/or
vertex-disjoint.

Finally, alternative algorithms are also available for the all-pairs shortest
path problem. One well-known approach is the $O(n^{3})$ Floyd-Warshall
algorithm, which is also able to handle graphs containing negative arc
weights. Another option in the presence of negative weights is the algorithm
of Johnson [13]. This operates by transforming the input graph into a second
graph that has no negative weights but which maintains the same shortest-paths
structure as the original. This can be achieved through a single application
of the Bellman-Ford algorithm. After this, $n$ applications of Dijkstra’s
algorithm can then be performed.

##  3 Dijkstra’s Algorithm

In this section, we give a more detailed description of Dijkstra’s algorithm.
We also describe its underlying data structures and derive its complexity.

To produce a shortest-path tree rooted at $s$, Dijkstra’s algorithm operates
by maintaining a set $D$ of so-called “distinguished vertices”. Initially,
only the source vertex $s$ is considered distinguished. During execution,
further vertices are then added to $D$, one at a time, until all reachable
vertices have been inserted. Two other data structures are also maintained.
First, a “label” $L(u)$ is stored for each vertex $u\in V$ in the graph.
During execution, $L(u)$ stores the length of the shortest $s$-$u$-path that
uses distinguished vertices only. Consequently, on termination of the
algorithm, $L(u)$ gives the length of the shortest $s$-$u$-path in the graph.
If a vertex $u$ has a label $L(u)=\infty$, then no $s$-$u$-path is possible.
Finally, a “predecessor” $P(u)$ is also stored for each vertex $u\in V$.
During execution, $P(u)$ stores the vertex that occurs before $u$ in the
shortest $s$-$u$ path (of length $L(u)$) that uses distinguished vertices
only. If a vertex $u$ has no predecessor, then $P(u)=\textsc{Null}$. On
termination, these predecessor values can be used to construct the shortest
paths from $s$ to all reachable vertices.

In its most basic form, Dijkstra’s algorithm can now be described by just
three steps. These are given in Algorithm 1. In these steps, note that one
vertex is inserted into $D$ at each iteration. This gives $O(n)$ iterations of
the algorithm in total. Furthermore, within each iteration, we need to
identify the vertex $u\notin D$ with the minimum label (an $O(n)$ operation)
and then examine (and possibly update) the labels of all vertices
$v\in\Gamma(u)$. This leads to an overall complexity of $O(m+n^{2})$. Since
the upper bound of $m=O(n^{2})$, this complexity can be simplified to
$O(n^{2})$. Output from an example run of this algorithm is shown in Figure 2.

input : An arc-weighted graph $G=(V,A)$ and source vertex $s\in V$

output : A populated label array $L$ and predecessor array $P$

1 Set $L(u)=\infty$ and $P(u)=\textsc{Null}$ for all $u\in V$. Also let
$D=\emptyset$, and set $L(s)=0$.

2 Choose a vertex $u\in V$ such that: (a) its value for $L(u)$ is minimal; (b)
$L(u)$ is less than $\infty$; and (c) $u$ is not in $D$. If no such vertex
exists, then end; otherwise insert $u$ into $D$ and go to Step 3.

3 For all neighbours $v\in\Gamma(u)$ that are not in $D$, if
$L(u)+w(u,v)<L(v)$ then set $L(v)=L(u)+w(u,v)$ and set $P(v)=u$. Now return to
Step 2.

Algorithm 1 Dijkstra’s Algorithm (Basic Form) Figure 2: Example output from
Algorithm 1 using the indicated eight-vertex graph and source vertex
$s=v_{1}$. In this case, the graph is undirected. Consequently, an arc $(u,v)$
exists if and only if the arc $(v,u)$ exists. In all cases, $w(u,v)=w(v,u)$,
as shown. The shortest path tree rooted at $s$ (defined by $P$) is shown by
the bold lines in the graph.

On completion of Dijkstra’s algorithm, the sequence of vertices that occurs in
each shortest path starting at $s$ is stored in $P$. The shortest $s$-$u$-path
(for all $u\in V$) can be constructed using the Get-Path procedure shown in
Algorithm 2. As shown, this operates by starting at $u$, and taking each
preceding vertex until the source $s$ is encountered. The $s$-$u$-path is then
the reverse of this sequence. For example, the shortest $s$-$v_{8}$-path from
Figure 2 is written $\pi=(s,v_{2},v_{3},v_{8})$. Note that, because all arc
weights are assumed to be nonnegative, the paths returned by Dijkstra’s
algorithm will always be “simple”. That is, they will never contain the same
vertex more than once.

input : The arc weighted graph $G$, predecessors $P$, source vertex $s$, and
an arbitrary vertex $u$

output : A vertex sequence $\pi$ corresponding to the shortest $s$-$u$-path in
$G$

1 Let $\pi=()$, and let $v=u$

2 if  _$P(u)\neq\textsc{Null}$_ then

3 while  _$v\neq s$_ do

4 Append $v$ to $\pi$ and set $v=P(v)$

5

6 Append $s$ to $\pi$ and then reverse $\pi$

7

Algorithm 2 Get-Path

As mentioned earlier, Dijkstra’s algorithm is also exact, meaning that it is
guaranteed to determine the shortest $s$-$u$-path for all reachable vertices
$u\in V$. Short proofs of this correctness can be found in several well-known
textbooks such as [9] and [17].

##  4 Using Priority Queues

Although the complexity of Dijkstra’s algorithm is $O(n^{2})$, for sparse
graphs its run times can be significantly improved by making use of a priority
queue. During execution, this priority queue is used to hold the labels of all
vertices that have been considered by the algorithm but that are not yet
marked as distinguished. It should also allow us to quickly identify the
undistinguished vertex that has the minimum label value. This “improved”
version of Dijkstra’s algorithm is expressed in Algorithm 3.

input : An arc weighted graph $G=(V,A)$, and source vertex $s\in V$

output : A populated label array $L$ and predecessor array $P$

1 For all $u\in V$, set $L(u)=\infty$, set $D(u)=$ false, and set $P(u)=$ Null

2 Set $L(s)=0$ and insert the ordered pair $(L(s),s)$ into $Q$

3 while  _$Q$ is not empty _ do

4 Let $(L(u),u)$ be the element in $Q$ with the minimum value for $L(u)$

5 Remove the element $(L(u),u)$ from $Q$

6 Set $D(u)=$ true

7 foreach  _$v\in\Gamma(u)$ such that $D(v)=$ false_ do

8 if  _$L(u)+w(u,v) <L(v)$_ then

9 if  _$L(v) <\infty$_ then

10 Decrease the key of $(L(v),v)$ to $L(u)+w(u,v)$. That is, replace the
element $(L(v),v)$ in $Q$ with the element $(L(u)+w(u,v),v)$

11

12 else

13 Insert the element $(L(u)+w(u,v),v)$ into $Q$

14

15 Set $L(v)=L(u)+w(u,v)$ and set $P(v)=u$

16

17

Algorithm 3 Dijkstra

As shown, the Dijkstra procedure in Algorithm 3 uses four data structures,
$D$, $L$, $P$ and $Q$. The first three of these contain $n$ elements and
should allow direct access (e.g., by using arrays). $D$ is used to mark the
distinguished vertices, while $L$ and $P$ hold the labels and predecessors of
each vertex as before. In this pseudocode, the priority queue is denoted by
$Q$. At each iteration, $Q$ is used to identify the element $(L(u),u)$,
representing the undistinguished vertex $u$ with the minimal label value. In
the remaining instructions, this element is removed from $Q$, $u$ is marked as
distinguished and, if necessary, adjustments are made to the labels of
undistinguished neighbours of $u$ and the corresponding entries in $Q$.

The running time of Dijkstra now depends on the data structure used for the
priority queue $Q$. Our first option here is to use a self-balancing binary
tree. These are a class of binary trees that automatically keep their height
logarithmic to the number of elements they contain. In C++, implementations of
a self-balancing binary tree are provided by the std::set container, usually
using a red-black tree [2].

A second option for $Q$ is to use a binary heap. Binary heaps are data
structures that take the form of _complete_ binary trees. Because of this
restriction, unlike self-balancing binary trees they can be implemented using
an array. This allows them to be stored in contiguous memory and also means
that the parent and children of any node can be determined using arithmetic on
array indices, as opposed to pointers. In C++, binary heaps are provided by
the std::priority_queue container [1].

Our final option for $Q$ is to use a Fibonacci heap. This data structure
operates by maintaining several heap-ordered trees. It also features better
amortized running times than the previous two options for some of the
operations used by Dijkstra’s algorithm. Note, however, that C++ does not
contain a Fibonacci heap in its standard library. Instead, a custom class is
required.

| Identify-Minimum | Remove-Minimum | Insert | Decrease-Key
---|---|---|---|---
Self-balancing binary tree | $O(1)$ | $O(1)^{\dagger}$ | $O(\lg n)$ | $O(\lg n)$
Binary heap | $O(1)$ | $O(\lg n)$ | $O(\lg n)$ | $O(\lg n)$
Fibonacci heap | $O(1)$ | $O(\lg n)^{\dagger}$ | $O(1)$ | $O(1)^{\dagger}$
Table 1: Complexities of several operators using self-balancing binary trees,
binary heaps, and Fibonacci heaps. Here, $n$ represents the number of elements
in the data structure. All entries are worst-case run times except for those
marked by †, which are amortised run times. Note that the Decrease-Key
operation is not available in the binary heap implementation provided by
std::priority_queue in C++ [1].

Table 1 considers these three alternatives for $Q$ and uses big O notation to
summarise the complexities of operations relevant to Dijkstra’s algorithm.
Further information on how these data structures work “under the hood” can be
found in [9]. For self-balancing binary trees, observe that the removal of the
minimum element $(L(u),u)$ (Line 5 of Dijkstra) takes constant amortised time.
At Line 10 of the algorithm, Decrease-Key operations are then carried out by
finding and removing the element $(L(v),v)$ in $Q$ and inserting the new
element $(L(u)+w(u,v),v)$. This process has a complexity $O(\lg n)$.
Similarly, the Insert operation on Line 12 also has a complexity of $O(\lg
n)$. Using a self-balancing binary tree for $Q$, therefore, leads to an
overall complexity for Dijkstra of $O(m\lg n)$.

As shown in Table 1, the operations with binary heaps show similar
complexities to self-balancing binary trees, though the removal of the minimum
element is now $O(\lg n$) as opposed to constant amortised time. Note,
however, that the std::priority_queue container in C++ does not feature the
functionality for performing Decrease-Key operations or for removing arbitrary
elements in $Q$. The use of this data structure, therefore, requires
modifications to the Dijkstra procedure. Specifically, at Line 10, instead of
replacing the element $(L(v),v)$ in $Q$ with $(L(u)+w(u,v),v)$, the latter
element is now simply inserted into $Q$ alongside the former. This means that,
unlike previously, a vertex $v$ can occur in several elements of $Q$. Because
of this, an additional check is now required between Lines 5 and 6 of Dijkstra
to determine if the selected vertex $u$ is already distinguished (that is, if
$D(u)=$ true). If this is the case, then the remaining steps in the while-loop
should not be considered, and the process should return to Line 3. Note that
this modification increases the overall complexity of Dijkstra to $O(m\lg m$).
On the other hand, binary heaps are usually seen to be faster than self-
balancing binary trees because they use contiguous memory, require fewer
memory allocations, and therefore tend to feature lower constant factors in
their operators. The effects of this trade-off are considered in Section 6.

Finally, Table 1 also shows the complexities of these operations using
Fibonacci heaps. For this data structure, Insert and Decrease-Key operations
both take place in constant amortised time. In particular, unlike self-
balancing binary trees, the Decrease-Key operation does not involve a removal
followed by an insertion; instead, it operates by identifying the correct
element $(L(v),v)$ in the heap, and then lowering the first value of the
element to $L(u)+w(u,v)$, modifying its position in the heap as applicable.
The use of a Fibonacci heap for $Q$, therefore, leads to an overall complexity
for Dijkstra of $O(m+n\lg n)$. This is better than both of the previous
options. Despite this, however, Fibonacci heaps are often considered to be
slow in practice due to their larger memory consumption and the high constant
factors contained in their operators. Indeed, it is noted by Cormen at al. [9]
that:

> “the constant factors and programming complexity of Fibonacci heaps makes
> them less desirable than ordinary binary (of $k$-ary) heaps for most
> applications. Thus Fibonacci heaps are predominantly of theoretical
> interest.”

This claim will also be investigated further in the next section. A summary of
the complexities of these algorithm variants is provided in Table 2.

Variant | Complexity | Comments
---|---|---
Basic form (Section 3) | $O(n^{2})$ | Minimum label $L(u)$ found using linear search. Optimal
|  | bound for dense graphs.
Self-balancing binary tree | $O(m\lg n)$ | More efficient than the basic form with sparse graphs.
Binary heap | $O(m\lg n)$ | Same complexity as previous. Note that the C++ implementation of
|  | binary heaps (std::priority_queue) does not allow the removal of
|  | arbitrary elements, so the size of the heap is $O(m)$, leading to a
|  | complexity of $O(m\lg m)$.
Fibonacci heap | $O(m+n\lg n)$ | Lower complexity than the previous variants.
Table 2: Complexities of the different variants of Dijkstra’s algorithm
considered in this paper. Recall that $n$ gives the number of vertices in the
graph, and $m$ is the number of arcs.

##  5 An Implementation

In this section, we consider four C++ implementations of Dijkstra’s algorithm.
The first three use a self-balancing binary tree (std::set), a binary heap
(std:priority_queue), and a Fibonacci heap, as described in the previous
section. The fourth implements the basic form of Dijkstra’s algorithm seen in
Section 3.

A complete listing of our code is shown in Appendix A and can be downloaded at
[3]. The bespoke FibonacciHeap class is defined on Lines 79 to 296. Lines 299
to 448 then give our four implementations of Dijkstra’s algorithm:
dijkstraFibonacci(...), dijkstraTree(...), dijkstraHeap(...), and
dijkstraBasic(...). The main() function on Lines 464 onwards applies these
methods to a small toy graph. The output is shown at the end of Appendix A.

The following features in this code should be noted.

  * •

Here, graphs are defined using the custom Graph class. Graph objects are
directed and are stored using the adjacency list representation for weighted
graphs [6]. It is assumed that the vertices are labelled from $0$ to $n-1$
(where $n$ is the number of vertices). The algorithm will also work perfectly
well on undirected graphs by ensuring that if an edge $\\{u,v\\}$ is present
in the graph, then both of the arcs $(u,v)$ and $(v,u)$ are present in the
adjacency list.

  * •

All arc weights in the graph are assumed to be nonnegative integers. The code
uses the inbuilt C++ constant INT_MAX to represent infinity values. Exceeding
this value will result in overflow and incorrect behaviour.

  * •

Each version of Dijkstra’s algorithm returns two arrays (vectors): the label
vector $L$ and the predecessor vector $P$.

  * •

Lines 450 to 462 of the code also give a function for extracting a path from
$P$. This corresponds to Algorithm 2 above.

In the following, all trials were performed on a 64-bit Windows 10 machine
with a 3.3 GHz Pro Intel Core i5-4590 CPU and 8 GB of RAM. In our case, the
code was compiled using Microsoft Visual Studio 2019 under release mode.

##  6 Empirical Evaluation

To assess the performance of the four variants of Dijkstra’s algorithm, timed
tests were carried out on two graph topologies: dense planar graphs and random
graphs. Planar graphs are a type of graph that can be drawn on a plane so that
no arcs intersect. Here, they were formed by randomly placing $n$ vertices
into a $(10,000\times 10,000)$-unit square before generating a random Delaunay
triangulation to give a graph with approximately (but not exceeding) $6n-12$
arcs. The weight of each arc was then set to the Euclidean distance between
its two endpoints, giving $w(u,v)=w(v,u)$ for all $(u,v)\in A$. An example is
shown in Figure 3. These planar graphs can be considered similar to road
networks which, as noted, are an important application area of shortest path
algorithms. Our random graphs, meanwhile, were generated by creating $n$
vertices and then, for each ordered pair of vertices $u,v$, adding the arc
$(u,v)$ with a probability $p$. This generation process leads to graphs with
approximately $pn(n-1)$ arcs. Here, each arc was assigned a weight between $1$
and $10,000$, selected at random.

Figure 3: Example of the dense planar graphs used in our trials. This
particular instance has $n=1000$ vertices and $m=5984$ arcs. Figure 4:
Execution times of the four algorithm variants with dense planar graphs. Each
point on the graph is the mean taken from 100 runs. Error bars show one
standard deviation on either side of the mean.

Execution times of the four algorithm variants with our dense planar graphs
are summarised in Figure 4 using a range of values for $n$. Across this range,
we see that the use of a self-balancing binary tree gives faster run times
than Fibonacci heaps, though these differences are very marginal. On the other
hand, the use of a binary heap gives a noticeable improvement over these two
variants, with run times dropping by approximately half. That said, each of
these three algorithms can compute shortest path trees of up to a million
vertices in well under a second with these graphs. As noted, the maximum
number of arcs in a directed planar graph is $6n-12$ meaning that, in these
cases, $m=O(n)$. Consequently, we can consider these three variants of
Dijkstra’s algorithm to have a complexity of $O(n\lg n)$ here. Given that
binary heaps involve fewer computational overheads than self-balancing binary
trees and Fibonacci heaps, this helps to explain the superior performance of
the binary heap variant in these cases. Finally, observe that the basic form
of Dijkstra’s algorithm shows both higher means and variances in run times
compared to the other three options. This is to be expected due to the
sparsity of these graphs. The gap between the basic version and the other
three also widens with increases in $n$, highlighting the quadratic growth
rate of the former.

Figure 5: Execution times of the four algorithm variants with random graphs
using $p=0.1$, $0.5$ and $0.9$ respectively. Each point on the graph is the
mean taken from 100 runs. Error bars are not shown here because all standard
deviations were seen to be less than 0.006 seconds.

The charts in Figure 5 show the results of the same experiments using random
graphs with $p=0.1$, $0.5$ and $0.9$ respectively. Note that, in these cases,
the number of arcs is much higher than the number of vertices; consequently,
the value of $m$ is now the dominant factor in the variants using priority
queues. Since $m\approx pn(n-1)$, increases to $p$ and/or $n$ therefore result
in longer run times.

For these graphs, we see that the relative performance of the Fibonacci heap
variant improves, with its results being almost indistinguishable from those
of the binary heap. The reasons for this are that, with these denser graphs,
the number of neighbours per vertex is larger. This brings a higher number of
Decrease-Key and Insert operations during execution, which are more efficient
with Fibonacci heaps. Despite this, however, the additional overheads required
by Fibonacci heaps seem to prevent the algorithm from improving on the binary
heap’s run times. In cases where $n$ and/or $p$ are high, the basic version of
Dijkstra’s algorithm also sometimes outperforms the variant using self-
balancing binary trees. This is particularly the case for graphs with $p=0.9$
where the number of arcs $m$ in these cases is close to $n^{2}$.

Finally, note that the largest instances considered here involve $n=10,000$
vertices, density $p=0.9$, and therefore approximately $90$ million arcs. In
our runs, such graphs were seen to occupy around $900$ MB of memory, but the
run times of the three variants using priority queues were still well below
half a second in all cases. The times taken to load these graphs into RAM are
not included in the above timings, however.

##  7 Conclusions

This paper has described the shortest path problem and shown how the
performance of Dijkstra’s algorithm can be affected by the choice of data
structure used for its priority queue. Using a C++ implementation tested over
a large range of problem instances, we have seen that the best-performing
algorithm does not always have the lowest complexity. Indeed, on the whole,
the best performance has been seen when using Dijkstra’s algorithm with a
binary heap, even though its complexity of $O(m\lg m)$ is higher than the
other variants. For dense graphs, however, the variant using a Fibonacci heap
features very similar run times to the binary heap version. Interestingly,
this binary-heap variant is also the chosen method of implementation in
several open-source libraries including NetworkX [5] and GraphHopper [4].

As noted, our current implementation operates by loading the entire graph into
RAM before execution. It also assumes that vertices are labelled with indices
from $0$ to $n-1$. In cases where these conditions are not possible, our code
will need to be modified to make use of associative arrays instead of vectors.
In C++ these are provided by the std::map and std::unordered_map containers.

## References

  * [1] C++ priority queue documentation.  https://en.cppreference.com/w/cpp/container/priority˙queue.  Accessed 2023-03-15.
  * [2] C++ set documentation.  https://en.cppreference.com/w/cpp/container/set.  Accessed 2023-03-15.
  * [3] C++ source code and results datasets.  https://doi.org/10.5281/zenodo.7741249.  Accessed 2023-03-15.
  * [4] Graphhopper implementation of Dijkstra’s algorithm (in java).  https://github.com/graphhopper/graphhopper/blob/master/core/src/main/java/com/graphhopper/routing/Dijkstra.java.  Accessed 2023-03-15.
  * [5] Networkx implementation of Dijkstra’s algorithm (in python).  https://networkx.org/documentation/stable/˙modules/networkx/algorithms/shortest˙paths/weighted.html#single˙source˙dijkstra.  Accessed 2023-03-15.
  * [6] Weighted graph representation.  https://www.tutorialspoint.com/weighted-graph-representation-in-data-structure.  Accessed 2023-03-15.
  * [7] H. Bast, D. Delling, A. Goldberg, M. Müller-Hannemann, T. Pajor, P. Sanders, D. Wagner, and R. Werneck.  Route Planning in Transportation Networks, pages 19–80.  Springer International Publishing, Cham, 2016.
  * [8] R. Bhandari.  Survivable Networks: Algorithms for Diverse Routing.  Kluwer Academic Publishers, 1999.
  * [9] T. Cormen, C. Leiserson, and R. Rivest.  Introduction to Algorithms.  The MIT Press, 1st edition, 2000.
  * [10] L. Fu, D. Sun, and L. Rilett.  Heuristic shortest path algorithms for transportation applications: State of the art.  Computers and Operations Research, 33(11):3324–3343, 2006.
  * [11] P. Hart, N. Nilsson, and B. Raphael.  A formal basis for the heuristic determination of minimum cost paths.  IEEE Transactions on Systems Science and Cybernetics, 4(2):100–107, 1968.
  * [12] P. Hart, N. Nilsson, and B. Raphael.  Correction to ‘A formal basis for the heuristic determination of minimum cost paths’.  ACM SIGART Bulletin, 37:28–29, 1972.
  * [13] D. Johnson.  Efficient algorithms for shortest paths in sparse networks.  Journal of the ACM, 24(1):1–13, 1977.
  * [14] R. Lewis.  Algorithms for finding shortest paths in networks with vertex transfer penalties.  Algorithms, 13(11), 2020.
  * [15] R. Lewis.  Who is the centre of the movie universe? Using python and networkx to analyse the social network of movie stars.  arXiv, 2002.11103, 2020.  https://arxiv.org/pdf/2002.11103.pdf.
  * [16] F. Moore.  The shortest path through a maze.  Technical report, Bell Telephone System, 1959.  vol. 3523.
  * [17] K. Rosen.  Discrete Mathematics and its Applications.  Mcgraw Hill, 8th edition, 2018.
  * [18] R. Sedgewick.  Algorithms in Java, Part 5: Graph Algorithms.  AddisonWesley Professional, 3rd edition, 2003.
  * [19] J. Yen.  Finding the $k$ shortest loopless paths in a network.  Management Science, 17(11):712–716, 1971.

##  Appendix A Code Listing and Example Run

The following code, together with this paper’s experimental data can be be
downloaded at [3]

⬇

1#include <iostream>

2#include <climits>

3#include <algorithm>

4#include <vector>

5#include <tuple>

6#include <set>

7#include <queue>

8#include <time.h>

9

10using namespace std;

11

12const int infty = INT_MAX;

13

14//Code for printing a vector

15template<typename T>

16ostream& operator<<(ostream& s, vector<T> t) {

17 s << "[";

18 for (size_t i = 0; i < t.size(); i++) {

19 s << t[i] << (i == t.size() - 1 ? "" : ",");

20 }

21 return s << "] ";

22}

23

24//Struct used for each element of the adjacency list.

25struct Neighbour {

26 int vertex;

27 int weight;

28};

29//Graph class (uses adjacency list)

30class Graph {

31 public:

32 int n; //Num. vertices

33 int m; //Num. arcs

34 vector<vector<Neighbour> > adj;

35 Graph(int n) {

36 this->n = n;

37 this->m = 0;

38 this->adj.resize(n, vector<Neighbour>());

39 }

40 ~Graph() {

41 this->n = 0;

42 this->m = 0;

43 this->adj.clear();

44 }

45 void addArc(int u, int v, int w) {

46 this->adj[u].push_back(Neighbour{ v, w });

47 this->m++;

48 }

49};

50

51//Struct and comparison operators used with std::set std::priority_queue)

52struct QueueItem {

53 int label;

54 int vertex;

55};

56struct minQueueItem {

57 bool operator() (const QueueItem& lhs, const QueueItem& rhs) const {

58 return tie(lhs.label, lhs.vertex) < tie(rhs.label, rhs.vertex);

59 }

60};

61struct maxQueueItem {

62 bool operator() (const QueueItem& lhs, const QueueItem& rhs) const {

63 return tie(lhs.label, lhs.vertex) > tie(rhs.label, rhs.vertex);

64 }

65};

66

67//Struct used for each Fibonacci heap node

68struct FibonacciNode {

69 int degree;

70 FibonacciNode* parent;

71 FibonacciNode* child;

72 FibonacciNode* left;

73 FibonacciNode* right;

74 bool mark;

75 int key;

76 int nodeIndex;

77};

78//Fibonacci heap class

79class FibonacciHeap {

80 private:

81 FibonacciNode* minNode;

82 int numNodes;

83 vector<FibonacciNode*> degTable;

84 vector<FibonacciNode*> nodePtrs;

85 public:

86 FibonacciHeap(int n) {

87 //Constructor function

88 this->numNodes = 0;

89 this->minNode = NULL;

90 this->degTable = {};

91 this->nodePtrs.resize(n);

92 }

93 ~FibonacciHeap() {

94 //Destructor function

95 this->numNodes = 0;

96 this->minNode = NULL;

97 this->degTable.clear();

98 this->nodePtrs.clear();

99 }

100 int size() {

101 //Number of nodes in the heap

102 return this->numNodes;

103 }

104 bool empty() {

105 //Is the heap empty?

106 if (this->numNodes > 0) return false;

107 else return true;

108 }

109 void insert(int u, int key) {

110 //Insert the vertex u with the specified key (value for L(u)) into the
Fibonacci heap. O(1) operation

111 this->nodePtrs[u] = new FibonacciNode;

112 this->nodePtrs[u]->nodeIndex = u;

113 FibonacciNode* node = this->nodePtrs[u];

114 node->key = key;

115 node->degree = 0;

116 node->parent = NULL;

117 node->child = NULL;

118 node->left = node;

119 node->right = node;

120 node->mark = false;

121 FibonacciNode* minN = this->minNode;

122 if (minN != NULL) {

123 FibonacciNode* minLeft = minN->left;

124 minN->left = node;

125 node->right = minN;

126 node->left = minLeft;

127 minLeft->right = node;

128 }

129 if (minN == NULL || minN->key > node->key) {

130 this->minNode = node;

131 }

132 this->numNodes++;

133 }

134 FibonacciNode* extractMin() {

135 //Extract the node with the minimum key from the heap. O(log n) operation,
where n is the number of nodes in the heap

136 FibonacciNode* minN = this->minNode;

137 if (minN != NULL) {

138 int deg = minN->degree;

139 FibonacciNode* currChild = minN->child;

140 FibonacciNode* remChild;

141 for (int i = 0; i < deg; i++) {

142 remChild = currChild;

143 currChild = currChild->right;

144 _existingToRoot(remChild);

145 }

146 _removeNodeFromRoot(minN);

147 this->numNodes\--;

148 if (this->numNodes == 0) {

149 this->minNode = NULL;

150 }

151 else {

152 this->minNode = minN->right;

153 FibonacciNode* minNLeft = minN->left;

154 this->minNode->left = minNLeft;

155 minNLeft->right = this->minNode;

156 _consolidate();

157 }

158 }

159 return minN;

160 }

161 void decreaseKey(int u, int newKey) {

162 //Decrease the key of the node in the Fibonacci heap that has index u.
O(1) operation

163 FibonacciNode* node = this->nodePtrs[u];

164 if (newKey > node->key) return;

165 node->key = newKey;

166 if (node->parent != NULL) {

167 if (node->key < node->parent->key) {

168 FibonacciNode* parentNode = node->parent;

169 _cut(node);

170 _cascadingCut(parentNode);

171 }

172 }

173 if (node->key < this->minNode->key) {

174 this->minNode = node;

175 }

176 }

177 private:

178 //The following are private functions used by the public methods above

179 void _existingToRoot(FibonacciNode* newNode) {

180 FibonacciNode* minN = this->minNode;

181 newNode->parent = NULL;

182 newNode->mark = false;

183 if (minN != NULL) {

184 FibonacciNode* minLeft = minN->left;

185 minN->left = newNode;

186 newNode->right = minN;

187 newNode->left = minLeft;

188 minLeft->right = newNode;

189 if (minN->key > newNode->key) {

190 this->minNode = newNode;

191 }

192 }

193 else {

194 this->minNode = newNode;

195 newNode->right = newNode;

196 newNode->left = newNode;

197 }

198 }

199 void _removeNodeFromRoot(FibonacciNode* node) {

200 if (node->right != node) {

201 node->right->left = node->left;

202 node->left->right = node->right;

203 }

204 if (node->parent != NULL) {

205 if (node->parent->degree == 1) {

206 node->parent->child = NULL;

207 }

208 else {

209 node->parent->child = node->right;

210 }

211 node->parent->degree\--;

212 }

213 }

214 void _cut(FibonacciNode* node) {

215 _removeNodeFromRoot(node);

216 _existingToRoot(node);

217 }

218 void _addChild(FibonacciNode* parentNode, FibonacciNode* newChildNode) {

219 if (parentNode->degree == 0) {

220 parentNode->child = newChildNode;

221 newChildNode->right = newChildNode;

222 newChildNode->left = newChildNode;

223 newChildNode->parent = parentNode;

224 }

225 else {

226 FibonacciNode* child1 = parentNode->child;

227 FibonacciNode* child1Left = child1->left;

228 child1->left = newChildNode;

229 newChildNode->right = child1;

230 newChildNode->left = child1Left;

231 child1Left->right = newChildNode;

232 }

233 newChildNode->parent = parentNode;

234 parentNode->degree++;

235 }

236 void _cascadingCut(FibonacciNode* node) {

237 FibonacciNode* parentNode = node->parent;

238 if (parentNode != NULL) {

239 if (node->mark == false) {

240 node->mark = true;

241 }

242 else {

243 _cut(node);

244 _cascadingCut(parentNode);

245 }

246 }

247 }

248 void _link(FibonacciNode* highNode, FibonacciNode* lowNode) {

249 _removeNodeFromRoot(highNode);

250 _addChild(lowNode, highNode);

251 highNode->mark = false;

252 }

253 void _consolidate() {

254 int deg, rootCnt = 0;

255 if (this->numNodes > 1) {

256 this->degTable.clear();

257 FibonacciNode* currNode = this->minNode;

258 FibonacciNode* currDeg, * currConsolNode;

259 FibonacciNode* temp = this->minNode, * itNode = this->minNode;

260 do {

261 rootCnt++;

262 itNode = itNode->right;

263 } while (itNode != temp);

264 for (int cnt = 0; cnt < rootCnt; cnt++) {

265 currConsolNode = currNode;

266 currNode = currNode->right;

267 deg = currConsolNode->degree;

268 while (true) {

269 while (deg >= int(this->degTable.size())) {

270 this->degTable.push_back(NULL);

271 }

272 if (this->degTable[deg] == NULL) {

273 this->degTable[deg] = currConsolNode;

274 break;

275 }

276 else {

277 currDeg = this->degTable[deg];

278 if (currConsolNode->key > currDeg->key) {

279 swap(currConsolNode, currDeg);

280 }

281 if (currDeg == currConsolNode) break;

282 _link(currDeg, currConsolNode);

283 this->degTable[deg] = NULL;

284 deg++;

285 }

286 }

287 }

288 this->minNode = NULL;

289 for (size_t i = 0; i < this->degTable.size(); i++) {

290 if (this->degTable[i] != NULL) {

291 _existingToRoot(this->degTable[i]);

292 }

293 }

294 }

295 }

296};

297//End of FibonacciHeap class

298

299tuple<vector<int>, vector<int>> dijkstraFibonacci(Graph& G, int s) {

300 //Dijkstra’s algorithm using a Fibonacci heap object

301 int u, v, w;

302 FibonacciHeap Q(G.n);

303 vector<int> L(G.n), P(G.n);

304 vector<bool> D(G.n);

305 for (int u = 0; u < G.n; u++) {

306 D[u] = false;

307 L[u] = infty;

308 P[u] = -1;

309 }

310 L[s] = 0;

311 Q.insert(s, 0);

312 while (!Q.empty()) {

313 u = Q.extractMin()->nodeIndex;

314 D[u] = true;

315 for (auto& neighbour : G.adj[u]) {

316 v = neighbour.vertex;

317 w = neighbour.weight;

318 if (D[v] == false) {

319 if (L[u] + w < L[v]) {

320 if (L[v] == infty) {

321 Q.insert(v, L[u] + w);

322 }

323 else {

324 Q.decreaseKey(v, L[u] + w);

325 }

326 L[v] = L[u] + w;

327 P[v] = u;

328 }

329 }

330 }

331 }

332 return make_tuple(L, P);

333}

334

335tuple<vector<int>, vector<int>> dijkstraTree(Graph& G, int s) {

336 //Dijkstra’s algorithm using a self-balancing binary tree (C++ set)

337 int u, v, w;

338 set<QueueItem, minQueueItem> Q;

339 vector<int> L(G.n), P(G.n);

340 vector<bool> D(G.n);

341 for (u = 0; u < G.n; u++) {

342 D[u] = false;

343 L[u] = infty;

344 P[u] = -1;

345 }

346 L[s] = 0;

347 Q.emplace(QueueItem{ 0,s });

348 while (!Q.empty()) {

349 u = (*Q.begin()).vertex;

350 Q.erase(*Q.begin());

351 D[u] = true;

352 for (auto& neighbour : G.adj[u]) {

353 v = neighbour.vertex;

354 w = neighbour.weight;

355 if (D[v] == false) {

356 if (L[u] + w < L[v]) {

357 if (L[v] == infty) {

358 Q.emplace(QueueItem{ L[u] + w, v });

359 }

360 else {

361 Q.erase({ L[v], v });

362 Q.emplace(QueueItem{ L[u] + w, v });

363 }

364 L[v] = L[u] + w;

365 P[v] = u;

366 }

367 }

368 }

369 }

370 return make_tuple(L, P);

371}

372

373tuple<vector<int>, vector<int>> dijkstraHeap(Graph& G, int s) {

374 //Dijkstra’s algorithm using a binary heap (C++ priority_queue)

375 int u, v, w;

376 priority_queue<QueueItem, vector<QueueItem>, maxQueueItem> Q;

377 vector<int> L(G.n), P(G.n);

378 vector<bool> D(G.n);

379 for (u = 0; u < G.n; u++) {

380 D[u] = false;

381 L[u] = infty;

382 P[u] = -1;

383 }

384 L[s] = 0;

385 Q.emplace(QueueItem{ 0,s });

386 while (!Q.empty()) {

387 u = Q.top().vertex;

388 Q.pop();

389 if (D[u] != true) {

390 D[u] = true;

391 for (auto& neighbour : G.adj[u]) {

392 v = neighbour.vertex;

393 w = neighbour.weight;

394 if (D[v] == false) {

395 if (L[u] + w < L[v]) {

396 Q.emplace(QueueItem{ L[u] + w, v });

397 L[v] = L[u] + w;

398 P[v] = u;

399 }

400 }

401 }

402 }

403 }

404 return make_tuple(L, P);

405}

406

407tuple<vector<int>, vector<int>> dijkstraBasic(Graph& G, int s) {

408 //Basic Dijkstra’s algorithm (O(n^2) complexity)

409 int u, v, w, minL;

410 size_t i, uPos;

411 vector<int> L(G.n), P(G.n), Candidates;

412 vector<bool> D(G.n);

413 for (u = 0; u < G.n; u++) {

414 D[u] = false;

415 L[u] = infty;

416 P[u] = -1;

417 }

418 L[s] = 0;

419 Candidates.push_back(s);

420 while (!Candidates.empty()) {

421 uPos = 0;

422 minL = L[Candidates[0]];

423 for (i = 1; i < Candidates.size(); i++) {

424 if (L[Candidates[i]] < minL) {

425 minL = L[Candidates[i]];

426 uPos = i;

427 }

428 }

429 u = Candidates[uPos];

430 swap(Candidates[uPos], Candidates.back());

431 Candidates.pop_back();

432 D[u] = true;

433 for (auto& neighbour : G.adj[u]) {

434 v = neighbour.vertex;

435 w = neighbour.weight;

436 if (D[v] == false) {

437 if (L[u] + w < L[v]) {

438 if (L[v] == infty) {

439 Candidates.push_back(v);

440 }

441 L[v] = L[u] + w;

442 P[v] = u;

443 }

444 }

445 }

446 }

447 return make_tuple(L, P);

448}

449

450vector<int> getPath(int u, int v, vector<int>& P) {

451 //Get the u-v-path specified by the predecessor vector P

452 vector<int> path;

453 int x = v;

454 if (P[x] == -1) return path;

455 while (x != u) {

456 path.push_back(x);

457 x = P[x];

458 }

459 path.push_back(u);

460 reverse(path.begin(), path.end());

461 return path;

462}

463

464int main() {

465 //Construct a small example graph. (A directed cycle on 5 vertices here.
All arcs have weight 10)

466 Graph G(5);

467 G.addArc(0, 1, 10);

468 G.addArc(1, 2, 10);

469 G.addArc(2, 3, 10);

470 G.addArc(3, 4, 10);

471 G.addArc(4, 0, 10);

472

473 //Set the source vertex and declare some variables

474 int s = 0;

475 vector<int> L, P;

476

477 //Execute Dijkstra’s algorithm using a Fibonacci heap

478 clock_t start = clock();

479 tie(L, P) = dijkstraFibonacci(G, s);

480 double duration1 = ((double)clock() - start) / CLOCKS_PER_SEC;

481

482 //Execute Dijkstra’s algorithm using a self-balancing binary tree

483 start = clock();

484 tie(L, P) = dijkstraTree(G, s);

485 double duration2 = ((double)clock() - start) / CLOCKS_PER_SEC;

486

487 //Execute Dijkstra’s algorithm using a binary heap

488 start = clock();

489 tie(L, P) = dijkstraHeap(G, s);

490 double duration3 = ((double)clock() - start) / CLOCKS_PER_SEC;

491

492 //Execute basic version of Dijkstra’s algorithm

493 start = clock();

494 tie(L, P) = dijkstraBasic(G, s);

495 double duration4 = ((double)clock() - start) / CLOCKS_PER_SEC;

496

497 //Output some information

498 cout << "Input graph has " << G.n << " vertices and " << G.m << " arcs\n";

499 cout << "Dijkstra with Fibonacci heap took " << duration1 << " sec.\n";

500 cout << "Dijkstra with self-balancing tree took " << duration2 << "
sec.\n";

501 cout << "Dijkstra with binary heap took " << duration3 << " sec.\n";

502 cout << "Dijkstra (basic form) took " << duration4 << " sec.\n";

503 cout << "Shortest paths from source to each vertex are as follows:\n";

504 for (int u = 0; u < G.n; u++) {

505 cout << "v-" << s << " to v-" << u << ",\t";

506 if (L[u] == infty) {

507 cout << "len = infinity. No path exists\n";

508 }

509 else {

510 cout << "len = " << L[u] << "\tpath = " << getPath(s, u, P) << "\n";

511 }

512 }

513}

Running this code produces the following output

⬇

1Input graph has 5 vertices and 5 arcs

2Dijkstra with Fibonacci heap took 1.9e-05 sec.

3Dijkstra with self-balancing tree took 1.2e-05 sec.

4Dijkstra with binary heap took 1.1e-05 sec.

5Dijkstra (basic form) took 1.5e-05 sec.

6Shortest paths from source to each vertex are as follows:

7v-0 to v-0, len = 0 path = []

8v-0 to v-1, len = 10 path = [0,1]

9v-0 to v-2, len = 20 path = [0,1,2]

10v-0 to v-3, len = 30 path = [0,1,2,3]

11v-0 to v-4, len = 40 path = [0,1,2,3,4]


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Observation of non-Hermitian edge burst in quantum dynamics

Lei Xiao  Beijing Computational Science Research Center, Beijing 100084, China
Wen-Tan Xue  Institute for Advanced Study, Tsinghua University, Beijing,
100084, China  Fei Song  Institute for Advanced Study, Tsinghua University,
Beijing, 100084, China  Yu-Min Hu  Institute for Advanced Study, Tsinghua
University, Beijing, 100084, China  Wei Yi<EMAIL_ADDRESS>CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei 230026, China  CAS Center For Excellence in Quantum Information
and Quantum Physics, Hefei 230026, China  Zhong Wang
<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua
University, Beijing, 100084, China  Peng Xue<EMAIL_ADDRESS>Beijing
Computational Science Research Center, Beijing 100084, China

###### Abstract

The non-Hermitian skin effect, by which the eigenstates of Hamiltonian are
predominantly localized at the boundary, has revealed a strong sensitivity of
non-Hermitian systems to the boundary condition. Here we experimentally
observe a striking boundary-induced dynamical phenomenon known as the non-
Hermitian edge burst, which is characterized by a sharp boundary accumulation
of loss in non-Hermitian time evolutions. In contrast to the eigenstate
localization, the edge burst represents a generic non-Hermitian dynamical
phenomenon that occurs in real time. Our experiment, based on photonic quantum
walks, not only confirms the prediction of the phenomenon, but also unveils
its complete space-time dynamics. Our observation of edge burst paves the way
for studying the rich real-time dynamics in non-Hermitian topological systems.

Figure 1: Experimental implementation. (a) The domain-wall geometry of the
non-Hermitian quantum walk. The operations of $S,R,L$ contained in $U$ are
pictorially shown. (b) Experimental setup. Photon pairs are created by the
spontaneous parametric down conversion process in a type-II cut PPKTP crystal.
One of the photon is injected into the quantum-walk interferometric network,
and the other is used as the trigger. The walker photon passes the polarizing
beam splitter (PBS) and the half-wave plate (HWP), so that its polarization is
prepared in the coin state $\ket{0}$. It then undertakes the quantum walk
through the network containing partially polarizing beam splitters (PPBSs),
HWPs, beam displacers (BDs). Finally, avalanche photodiodes (APDs) are used to
detect the walker photons that coincide with the trigger photons. Figure 2:
Edge burst in non-Hermitian quantum walks. The fixed parameters are
$\theta_{1,2}^{L}=0.85\pi$, $\theta_{1}^{R}=0.12\pi$ and $\gamma=0.8$. (a)(b)
Brillouin zone (BZ) and generalized Brillouin zone (GBZ) for
$\theta_{2}^{R}=0.12\pi$ and $\theta_{2}^{R}=0.48\pi$. (c)(d) Energy spectra
(for the right region in which the walker is initialized) under the periodic
boundary condition (PBC) for two indicated values of $\theta_{2}^{R}$. (e)(f)
Experimentally measured $P(x)$ of a $14$-step non-Hermitian quantum walk with
the initial state $\ket{x_{0}=10}\otimes\ket{0}$. (g)(h) The space-time-
resolved loss probability $p(x,t)$ for the two values of $\theta_{2}^{R}$.
Error bars represent the statistical uncertainty under the assumption of
Poissonian statistics.  Figure 3: (a)(b) Experimentally measured time-
dependent total loss probability $P(t)$ for different starting positions
$x_{0}=5,6,7,8,9,10$, respectively. For (a), $\theta_{2}^{R}=0.12\pi$; for
(b), $\theta_{2}^{R}=0.48\pi$. Other parameters are the same as those in Fig.
2. (c) The measured survival probability at the domain wall,
$|\psi(x=-1,t)|^{2}$, for different starting positions.
$\theta_{2}^{R}=0.12\pi$ (upper panel) and $0.48\pi$ (lower panel). (d) (e)
Experimentally measured $P(x)$ ($T=14$) for $x_{0}=6,8,10$;
$\theta_{2}^{R}=0.12\pi$ for (d) and $0.48\pi$ for (e). (f) The measured
relative height $P_{\text{domain}}/P_{\text{min}}$ versus $x_{0}$. Hollow and
solid symbols represent numerically evaluated results and experimental data,
respectively. Dashed lines are obtained by numerical fitting experimental
data.

Non-Hermitian physics has attracted increasing attention in a vast variety of
contexts ranging from classical waves to open quantum systems Ashida ; BBK21 .
Intriguingly, the spatial boundary plays a much more dramatic role in non-
Hermitian systems than in Hermitian ones. In particular, for certain non-
Hermitian systems, the eigenstates concentrate predominantly at the boundary,
which is known as the non-Hermitian skin effect (NHSE) YW18 ; YSW18 ; KEB+18 ;
AVT18 ; LT19 ; SYW19 ; OKS+20 ; XDW+20 ; HHI+20 ; Ma ; GBW+20 ; ZYF22 . Among
many other consequences, it implies a fundamental revision of the principle of
bulk-boundary correspondence XDW+20 ; HHI+20 .

Whereas the NHSE has revealed intriguing static properties such as novel
behaviors of eigenstates and energy spectra, in this work we unveil a striking
dynamic boundary effect in non-Hermitian systems. We experimentally observe
that in a class of lossy quantum walk of single photons, the loss rate is
drastically enhanced at the boundary. Specifically, for a lossy particle
initially located at a position far from the boundary of a lattice system, the
space-resolved loss has a surprisingly high boundary peak, in sharp contrast
to the common expectation that the particle loss should decay away from the
initial position. Remarkably, the relative height of the edge peak even grows
as the distance between the initial position and boundary increases. This
striking phenomenon, dubbed non-Hermitian edge burst, has been predicted in
recent theories XHS+22 ; WLZ21 .

Since both the NHSE and edge burst involve boundary localization, it is
tempting to attribute the latter to the former. However, it turns out that
NHSE does not guarantee the emergence of edge burst. Closing the gap of the
imaginary part of energy spectrum (i.e., the imaginary gap or dissipative gap)
is the other necessary condition, which highlights the rich implication of
spectral profile and topology in non-Hermitian systems ZYF22 ; OKS+20 . At a
deeper level, a novel dynamic bulk-edge scaling relation has been suggested as
the origin of edge burst XHS+22 . Thus, the edge burst signifies an
unprecedented interplay between non-Hermitian topological physics and non-
Hermitian dynamical phenomena.

Lossy quantum walk.—To study the non-Hermitian edge burst, we design a one-
dimensional quantum walk ADZ93 ; FG98 ; W01 ; C09  with the Floquet operator

$U=R\left(\frac{\theta_{2}}{2}\right)SR\left(\frac{\theta_{1}}{2}\right)L(\gamma).$
(1)

The shift operator $S=\sum_{x}|x-1\rangle\langle x|\otimes|0\rangle\langle
0|+|x+1\rangle\langle x|\otimes|1\rangle\langle 1|$, so that the walker’s
position is shifted from the site $x$ to $x-1$ or $x+1$ according to the coin
state $\ket{0}$ or $\ket{1}$. The coin state is rotated along the $y$ axis by
$R(\theta)=\mathds{1}_{\text{w}}\otimes e^{-i\theta\sigma_{y}}$, where
$\mathds{1}_{\text{w}}=\sum_{x}\ket{x}\bra{x}$ is the identity operator. The
operator $L(\gamma)=\mathds{1}_{\text{w}}\otimes\begin{pmatrix}1&0\\\
0&e^{-2\gamma}\end{pmatrix}$ generates a state-selective loss. For our
photonic platform, it is more convenient to create a domain wall instead of an
open boundary [see Fig. 1(a)]. The left (L) and right (R) regions are
characterized by coin parameters $\theta_{1,2}^{L}$ and $\theta_{1,2}^{R}$,
respectively. The dynamics of the non-Hermitian quantum walk follows

$|\psi(t)\rangle=U^{t}|\psi(0)\rangle,$ (2)

where $\ket{\psi(0)}$ is the initial state and $t$ is the integer discrete
time. One can also define an effective non-Hermitian Hamiltonian
$H_{\text{eff}}$ by $U=\exp(-iH_{\text{eff}})$, which shares the same
eigenstates as $U$.

The Floquet operator $U$ defined in Eq. (1) and the associated
$H_{\text{eff}}$ exhibit the NHSE, which originates from the state-dependent
directional hoppings built in the model (akin to Refs. XDW+20 ; XDW+21 ). In
the presence of a domain wall [Fig. 1(a)], all the eigenstates of $U$ exhibit
localization at the domain wall when the non-Hermiticity is nonzero, i.e.,
$\gamma\neq 0$. Accordingly, the generalized Brillouin zone (GBZ) deviates
from the unit circle [see Figs. 2(a) and (b)] YW18 ; Yokomizo2019 ; DY19 .
Here, we focus on two sets of parameters, $\theta_{2}^{R}=0.12\pi$ and
$\theta_{2}^{R}=0.48\pi$, with other parameters fixed as
$\theta_{1,2}^{L}=0.85\pi$, $\theta_{1}^{R}=0.12\pi$, and $\gamma=0.8$. In
Figs. 2(c) and (d), we show the energy spectrum of $H_{\text{eff}}$, which
clearly indicates that the imaginary gap (the gap between $0$ and the maximum
imaginary part of the spectrum) is zero for $\theta_{2}^{R}=0.12\pi$ but
nonzero for $\theta_{2}^{R}=0.48\pi$. In fact, the imaginary gap vanishes
along the lines $\theta_{1}=2\pi n\pm\theta_{2}$ ($n\in\mathbb{Z}$) (see
Supplementary Information).

Observation of edge burst.—In our experiment, a walker is initialized at a
site $x_{0}$, which evolves under Eq. (2) in discrete time steps. The key
quantity for edge burst is the probability $P(x)$ that the walker escapes from
the position $x$. In practice, one can measure the space-time-resolved loss
$p(x,t)$ from $t=1$ to $t=T$, with $T$ being a large integer so that the loss
is almost complete. The sum over $t$ then gives

$P(x)=\sum_{t=1}^{T}p(x,t).$ (3)

According to the specific form of loss adopted here, we have

$\displaystyle
p(x,t)=(1-e^{-4\gamma})|\bra{1}\otimes\bra{x}\psi(t-1)\rangle|^{2}.$ (4)

It may also be written as
$p(x,t)=|\bra{1}\otimes\bra{x}M|\psi(t-1)\rangle|^{2}$ with
$M=\mathds{1}_{\text{w}}\otimes\begin{pmatrix}0&0\\\
0&\sqrt{1-e^{-4\gamma}}\end{pmatrix}$, which can be implemented by a partial
measurement via the PPBS [see Fig. 1(a)] at the time step $t$. We also define
a time-dependent total loss probability

$P(t)=\sum_{t^{\prime}=1}^{t}\sum_{x}p(x,t^{\prime}),$ (5)

so that the survival probability after a $t$-step evolution is $1-P(t)$. In
our quantum-walk platform, $p(x,t)$ can be readily extracted from photon-
number measurements (see Methods), and $P(x)$, $P(t)$ can be obtained from
Eqs. (3)(5).

We implement a $14$-step ($T=14$) quantum walk with initial walker location
$x_{0}=10$. The space-resolved loss probability $P(x)$ is shown in Figs. 2(e)
and (f) for the aforementioned two sets of parameters. In both (e)
($\theta_{2}^{R}=0.12\pi$) and (f)($\theta_{2}^{R}=0.48\pi$), we observe that
the loss probability initially decays away from $x_{0}$. Moreover, the $P(x)$
profile is asymmetric around $x_{0}$, which can be naturally attributed to the
NHSE.

The surprising feature is an exceptionally high peak emerging at the domain
wall in Fig. 2(e). Intuitively, one may resort to the NHSE to explain this
edge burst. However, the NHSE is also strong for the parameters of Fig. 2(f),
yet the edge burst is not seen there. Therefore, the origin of edge burst
cannot be explained by the NHSE alone. In fact, the imaginary gap plays an
essential role here XHS+22 . The corresponding imaginary gap, shown in Figs.
2(c) and (d), is zero and nonzero for Figs. 2(e) and (f), respectively.

To unveil the space-time profile of walker’s loss, we plot $p(x,t)$ for the
above two sets of parameters. Figs. 2(g) and (h) show that the walker
propagates almost ballistically with concurrent loss along the trajectory. In
the case of edge burst, a large loss peak in $p(x,t)$ emerges when the walker
hits the domain wall. It also indicates that the burst occurs around a
particular time, before which it is indiscernible.

Furthermore, we vary the initial position $x_{0}=5,6,7,8,9,10$ and measure the
time-dependent loss probability $P(t)$. As shown in Fig. 3(a), for
$\theta_{2}^{R}=0.12\pi$ (with edge burst), $P(t)$ suddenly increases near the
domain wall. In contrast, in Fig. 3(b), for $\theta_{2}^{R}=0.48\pi$ (without
edge burst), $P(t)$ increases steadily with $t$ without sudden change.
Similarly, the space-resolved survival probability $|\psi(x=-1,t)|^{2}$ at the
domain wall at each step $t$ behaves differently with and without the edge
burst [see Fig. 3(c)]. The value of $|\psi(x=-1,t)|^{2}$ is significantly
larger in the presence of edge burst. In Fig. 3(d), we show that the edge
burst remains robust when the starting position varies. In contrast, when the
edge burst is absent, $P(x)$ decays rapidly as $x_{0}$ moves away from the
domain wall [see Fig. 3(e)].

To further characterize the edge burst, we measure the relative height
$P_{\text{domain}}/P_{\text{min}}$, where $P_{\text{domain}}\equiv P(x=-1)$ is
the probability that the photon escapes from the domain wall $x=-1$, and
$P_{\text{min}}\equiv\min_{x=-1,\cdots,x_{0}}\\{P(x)\\}$ is the minimum of
$P(x)$ in the interval between the initial location $x_{0}$ and the domain
wall location $x=-1$. The edge burst is characterized by
$P_{\text{domain}}/P_{\text{min}}\gg 1$, while its absence means that
$P_{\text{domain}}/P_{\text{min}}$ is on the order of unity. As shown in Fig.
3(f), for $\theta_{2}^{R}=0.48\pi$, the measured relative height remains close
to $1$ as $x_{0}$ increases. In stark contrast, for $\theta_{2}^{R}=0.12\pi$,
the relative height increases with $x_{0}$ and fits well with a linear
relation $P_{\text{domain}}/P_{\text{min}}\sim x_{0}$. Thus, the relative
height grows as the initial walker position moves away from the domain wall.
While counterintuitive, this behavior is a consequence of a novel bulk-edge
scaling relation XHS+22 .

Discussions.—We present the first experimental observation of the non-
Hermitian edge burst by using discrete-time non-Hermitian quantum walk of
photons. Our experiment not only demonstrates that edge burst originates from
the intriguing interplay between two unique non-Hermitian concepts, the NHSE
and imaginary gap, but also unveils the real-time dynamics of this phenomenon.
The observation of non-Hermitian edge burst paves the way for investigating
the real-time dynamics in non-Hermitian topological systems, which remains
largely unexplored. From a practical perspective, the edge burst may offer a
promising non-Hermitian approach for the on-demand harvesting of light or
particles at a prescribed position.

Methods

Implementation.—For the experimental implementation, we adopt the scheme of
single-photon discrete-time quantum walks illustrated in Fig. 1(b). Photon
pairs are created by spontaneous parametric down conversion, where a $20$mm
type-II periodically poled potassium titanyl phosphate (PPKTP) crystal is
pumped by a $405$nm continuous wave diode laser with the power of $1$mW. One
photon serves as a trigger, and the other as a heralded single photon
undertaking the quantum walk. The photon polarizations are adopted as the coin
state. The walker photon is initialized in the spatial mode $|x_{0}\rangle$
with the internal state $|0\rangle$, i.e.
$|\psi(0)\rangle=|x_{0}\rangle\otimes|0\rangle$. The localized initial state
is prepared by passing the walker photons through a half-wave plate (HWP) and
a polarizing beam splitter (PBS).

For the quantum-walk dynamics, the shift operator $S$ is implemented by a beam
displacer (BD) whose optical axis is cut in the way so that the vertically
polarized photons are directly transmitted and the horizontally polarized
photons are laterally displaced into a neighboring mode. The coin rotation
$R(\frac{\theta_{1(2)}}{2})$ is realized by two HWPs at $0$ and
$\frac{\theta_{1(2)}}{4}$, respectively. The loss operator $L(\gamma)$ is
realized by a partially polarizing beam splitter (PPBS), which completely
transmits the coin state $\ket{0}$ but reflects the coin state $\ket{1}$ with
a probability $e^{-4\gamma}$. At last, avalanche photodiodes (APDs) are used
to detect the walker photons coinciding with the trigger photons. The total
number of coincidences is approximately $23000$.

The measurements are based on photon-number counting. The space-time-resolved
probability $p(x,t)$ can be calculated from the photon number through

$p(x,t)=\frac{N(x,t)}{\sum_{x^{\prime}}N^{\prime}(x^{\prime},t)+\sum_{t^{\prime}=1}^{t}\sum_{x^{\prime}}N(x^{\prime},t^{\prime})},$
(6)

where $N(x,t)$ is the number of photons escaping from the position $x$ at the
time step $t$, and $N^{\prime}(x,t)$ is the number of remaining photons at $x$
after a $t$-step evolution.

Finally, the space-resolved survival probability at $x$ can be calculated as

$|\psi(x,t)|^{2}=\frac{N^{\prime}(x,t)}{\sum_{x^{\prime}}N^{\prime}(x^{\prime},t)+\sum_{t^{\prime}=1}^{t}\sum_{x^{\prime}}N(x^{\prime},t^{\prime})}.$
(7)

Note. After completing this work, we learned of a related experiment by a team
at Southern University of Science and Technology.

Acknowledgments

This work has been supported by the National Natural Science Foundation of
China (Grant Nos. 92265209, 12025401, 12125405, 11974331 and 12104036.)

## References

  * (1) Ashida, Y., Gong, Z. & Ueda, M. Non-Hermitian physics. Adv. Phys. 69, 249-435 (2020).
  * (2) Bergholtz, E. J., Budich, J. C. & Kunst, F. K. Exceptional topology of non-Hermitian systems. Rev. Mod. Phys. 93, 015005 (2021).
  * (3) Yao, S. & Wang, Z. Edge states and topological invariants of non-Hermitian systems. Phys. Rev. Lett. 121, 086803 (2018).
  * (4) Yao, S., Song, F. & Wang, Z. Non-Hermitian chern bands. Phys. Rev. Lett. 121, 136802 (2018).
  * (5) Kunst, F. K., Edvardsson, E., Budich, J. C. & Bergholtz, E. J. Biorthogonal bulk-boundary correspondence in non-Hermitian systems. Phys. Rev. Lett. 121, 026808 (2018).
  * (6) Martinez Alvarez, V. M., Barrios Vargas, J. E. & Foa Torres, L. E. F. Non-Hermitian robust edge states in one dimension: Anomalous localization and eigenspace condensation at exceptional points. Phys. Rev. B 97, 121401(R) (2018).
  * (7) Lee, C. H. & Thomale, R. Anatomy of skin modes and topology in non-Hermitian systems. Phys. Rev. B 99, 201103(R) (2019).
  * (8) Song, F., Yao, S. & Wang, Z. Non-Hermitian skin effect and chiral damping in open quantum systems. Phys. Rev. Lett. 123, 170401 (2019).
  * (9) Zhang, K., Yang, Z. & Fang, C. Correspondence between winding numbers and skin modes in non-Hermitian systems. Phys. Rev. Lett. 125, 126402 (2020).
  * (10) Okuma, N., Kawabata, K., Shiozaki, K. & Sato, M. Topological origin of non-Hermitian skin effects. Phys. Rev. Lett. 124, 086801 (2020).
  * (11) Xiao, L., Deng, T., Wang, K., Zhu, G., Wang, Z., Yi W. & Xue, P. Non-Hermitian bulk-boundary correspondence in quantum dynamics. Nat. Phys. 16, 761-766 (2020).
  * (12) Helbig, T., Hofmann, T., Imhof, S., Abdelghany, M., Kiessling, T., Molenkamp, L. W., Lee, C. H., Szameit, A., Greiter, M. & Thomale, R. Generalized bulk-boundary correspondence in non-Hermitian topolectrical circuits. Nat. Phys. 16, 747-750 (2020).
  * (13) Wang, W., Wang, X. & Ma, G. Non-Hermitian morphing of topological modes. Nature 608, 50-55 (2022).
  * (14) Ghatak, A., Brandenbourger, M., Van Wezel, J. & Coulais, C. Observation of non-Hermitian topology and its bulk-edge correspondence in an active mechanical metamaterial. Proc. Natl. Acad. Sci. U.S.A. 117, 29561-29568 (2020).
  * (15) Xue, W. T., Hu, Y. M., Song, F. & Wang, Z. Non-Hermitian edge burst. Phys. Rev. Lett. 128, 120401 (2022).
  * (16) Wang, L., Liu, Q. & Zhang, Y. Quantum dynamics on a lossy non-Hermitian lattice, Chin. Phys. B 30, 020506 (2021).
  * (17) Aharonov, Y., Davidovich, L. & Zagury, N. Quantum random walks. Phys. Rev. A 48, 1687 (1993).
  * (18) Farhi, E. & Gutmann, S. Quantum computation and decision trees. Phys. Rev. A 58, 915 (1998).
  * (19) Watrous, J. Quantum simulations of classical random walks and undirected graph connectivity. J. Comput. Syst. Sci. 62, 376-391 (2001).
  * (20) Childs, A. M. Universal computation by quantum walk. Phys. Rev. Lett. 102, 180501 (2009).
  * (21) Xiao, L., Deng, T., Wang, K., Wang, Z., Yi, W. & Xue, P. Observation of non-Bloch parity-time symmetry and exceptional points. Phys. Rev. Lett. 126, 230402 (2021).
  * (22) Yokomizo, K. & Murakami, S. Non-Bloch band theory of non-Hermitian systems. Phys. Rev. Lett. 123, 066404 (2019).
  * (23) Deng, T. S. & Yi, W. Non-Bloch topological invariants in a non-Hermitian domain wall system. Phys. Rev. B 100, 035102 (2019).

Supplemental Material for “Observation of non-Hermitian edge burst in quantum
dynamics”

##  Appendix A Effective Hamiltonian and generalized Brillouin zone

In this section, we derive an expression for the effective Hamiltonian
$H_{\text{eff}}$ in momentum space. First, we transform the real-space
nonunitary Floquet operator $U$ and its conjugate transpose $U^{\dagger}$ into
the momentum-space $U_{k}$ and $U_{k}^{\dagger}$:

$\displaystyle
U_{k}=d_{0}\sigma_{0}-id_{1}\sigma_{1}-id_{2}\sigma_{2}-id_{3}\sigma_{3},$
$\displaystyle
U_{k}^{\dagger}=d_{0}^{\ast}\sigma_{0}+id_{1}^{\ast}\sigma_{1}+id_{2}^{\ast}\sigma_{2}+id_{3}^{\ast}\sigma_{3}$
(S1)

where $\sigma_{1,2,3}$ are the Pauli matrices and $\sigma_{0}$ is the identity
matrix, and

$\displaystyle d_{0}=e^{-\gamma}(\cosh\gamma\cos
k\cos\frac{\theta_{1}+\theta_{2}}{2}+i\sinh\gamma\sin
k\cos\frac{\theta_{1}-\theta_{2}}{2}),$ $\displaystyle
d_{1}=e^{-\gamma}(\cosh\gamma\sin
k\sin\frac{\theta_{1}-\theta_{2}}{2}+i\sinh\gamma\cos
k\sin\frac{\theta_{1}+\theta_{2}}{2}),$ $\displaystyle
d_{2}=e^{-\gamma}(\cosh\gamma\cos
k\sin\frac{\theta_{1}+\theta_{2}}{2}-i\sinh\gamma\sin
k\sin\frac{\theta_{1}-\theta_{2}}{2}),$ $\displaystyle
d_{3}=e^{-\gamma}(-\cosh\gamma\sin
k\cos\frac{\theta_{1}-\theta_{2}}{2}+i\sinh\gamma\cos
k\cos\frac{\theta_{1}+\theta_{2}}{2}).$ (S2)

Note that the relation
$\sqrt{d_{0}^{2}+d_{1}^{2}+d_{2}^{2}+d_{3}^{2}}=e^{-\gamma}$ is satisfied. The
eigenvalue and eigenvector can be derived from

$\displaystyle
U_{k}|\psi_{\pm}\rangle=\lambda_{\pm}|\psi_{\pm}\rangle,~{}U_{k}^{\dagger}|\chi_{\pm}\rangle=\lambda^{\ast}_{\pm}|\chi_{\pm}\rangle.$
(S3)

Straightforward calculations lead to

$\lambda_{\pm}=d_{0}\pm it_{0},~{}\lambda^{\ast}_{\pm}=d_{0}^{\ast}\mp
it_{0}^{\ast},$ (S4) $\displaystyle|\psi_{\pm}\rangle$
$\displaystyle=\frac{1}{d_{1}+id_{2}}\left(\begin{array}[]{c}d_{3}\mp t_{0}\\\
d_{1}+id_{2}\end{array}\right),$ (S7) $\displaystyle\langle\chi_{\pm}|$
$\displaystyle=\frac{1}{d_{1}-id_{2}}(d_{3}\mp t_{0},d_{1}-id_{2}),$ (S8)

where $t_{0}=\sqrt{e^{-2\gamma}-d_{0}^{2}}$. Since the effective Hamiltonian
$H_{\text{eff}}(k)$ is related to $U_{k}$ through
$U_{k}=e^{-iH_{\text{eff}}}$, the quasienergy spectrum of $H_{\text{eff}}(k)$
is

$E_{\pm}(k)=i\ln\lambda_{\pm}(k)=\pm\arccos(\cosh\gamma\cos
k\cos\frac{\theta_{1}+\theta_{2}}{2}+i\sinh\gamma\sin
k\cos\frac{\theta_{1}-\theta_{2}}{2})-i\gamma.$ (S9)

Specifically, for $\theta_{1}=\theta_{2}$, we have

$E_{-}(k=\pi/2)=-\text{arccos}(i\sinh\gamma)-i\gamma=-\text{arccos}(\sin
i\gamma)-i\gamma=-\pi/2,$ (S10)

so that $\text{Im}[E_{-}(k=\pi/2)]=0$, i.e., the imaginary gap closes at
$k=\pi/2$. For $\theta_{1}=-\theta_{2}$, the imaginary gap closes at $k=0$
because

$E_{+}(k=0)=\text{arccos}(\cosh\gamma)-i\gamma=\text{arccos}(\cos(i\gamma))-i\gamma=0.$
(S11)

Thus, the imaginary gap closes when $\theta_{1}=2\pi n\pm\theta_{2}$
($n\in\mathbb{Z}$). While the eigenvectors in Eq. (S8) are not orthogonal, one
can derive a set of bi-orthonormal eigenvectors
$\\{|\tilde{\psi}_{\pm}\rangle,|\tilde{\chi}_{\pm}\rangle\\}$:

$\displaystyle|\tilde{\psi}_{\pm}\rangle=$
$\displaystyle\frac{|\psi_{\pm}\rangle}{\sqrt{\langle\chi_{\pm}|\psi_{\pm}\rangle}}=\frac{1}{\sqrt{2t_{0}(t_{0}\mp
d_{3})}}\left(\begin{array}[]{c}d_{3}\mp t_{0}\\\
d_{1}+id_{2}\end{array}\right),$ (S14)
$\displaystyle\langle\tilde{\chi}_{\pm}|=$
$\displaystyle\frac{\langle\chi_{\pm}|}{\sqrt{\langle\chi_{\pm}|\psi_{\pm}\rangle}}=\frac{1}{\sqrt{2t_{0}(t_{0}\mp
d_{3})}}(d_{3}\mp t_{0},d_{1}-id_{2}),$ (S15)

which satisfy

$\langle\tilde{\chi}_{m}|\tilde{\psi}_{n}\rangle=\delta_{mn},~{}\sum_{m=+,-}|\tilde{\psi}_{m}\rangle\langle\tilde{\chi}_{m}|=\mathds{1}.$
(S16)

It follows that

$U_{k}=\lambda_{+}|\tilde{\psi}_{+}\rangle\langle\tilde{\chi}_{+}|+\lambda_{-}|\tilde{\psi}_{-}\rangle\langle\tilde{\chi}_{-}|,$
(S17)

and the effective Hamiltonian $H_{\text{eff}}(k)$ can be written as

$H_{\text{eff}}=i\ln{\lambda_{+}}|\tilde{\psi}_{+}\rangle\langle\tilde{\chi}_{+}|+i\ln{\lambda_{-}}|\tilde{\psi}_{-}\rangle\langle\tilde{\chi}_{-}|.$
(S18)

To derive the generalized Brillouin zone (GBZ) YW18 ; Yokomizo2019 , we
rewrite the Floquet operator $U$ as

$\displaystyle U=\sum_{x}|x-1\rangle\langle x|\otimes A_{0}+|x+1\rangle\langle
x|\otimes A_{1},$ (S19)

where

$\displaystyle A_{0}=$ $\displaystyle
R_{c}(\frac{\theta_{2}}{2})P_{0}R_{c}(\frac{\theta_{1}}{2})L_{c}(\gamma),$
$\displaystyle A_{1}=$ $\displaystyle
R_{c}(\frac{\theta_{2}}{2})P_{1}R_{c}(\frac{\theta_{1}}{2})L_{c}(\gamma),$
(S20)

with $L_{c}(\gamma)=\left(\begin{array}[]{@{\,}cc@{\,}}1&0\\\
0&e^{-2\gamma}\\\ \end{array}\right)$, $R_{c}(\theta)=e^{-i\theta\sigma_{y}}$,
$P_{0}=|0\rangle\langle 0|$ and $P_{1}=|1\rangle\langle 1|$. In view of the
translational symmetry inside the bulk, the eigenstate $|\varphi\rangle$ of
$U$ can be expressed as

$|\varphi\rangle=\sum_{x,j}\beta_{j}^{x}|x\rangle\otimes|\phi_{j}\rangle_{c},$
(S21)

where $|\phi_{j}\rangle_{c}$ is the coin state and $\beta_{j}$ is the spatial-
mode function. Inserting Eq. (S21) into eigen-equation
$U|\varphi\rangle=\lambda|\varphi\rangle$, we obtain

$(A_{0}\beta+\frac{A_{1}}{\beta}-\lambda)|\phi\rangle_{c}=0,$ (S22)

which has nontrivial solutions only when

$\det[A_{0}\beta+\frac{A_{1}}{\beta}-\lambda]=0.$ (S23)

In an explicit form, Eq. (S23) is a quadratic equation of $\beta$:

$[\sin(\frac{\theta_{1}}{2})\sin(\frac{\theta_{2}}{2})-e^{2\gamma}\cos(\frac{\theta_{1}}{2})\cos(\frac{\theta_{2}}{2})]\beta^{2}+(\frac{1}{\lambda}+e^{2\gamma}\lambda)\beta+e^{2\gamma}\sin(\frac{\theta_{1}}{2})\sin(\frac{\theta_{2}}{2})-\cos(\frac{\theta_{1}}{2})\cos(\frac{\theta_{2}}{2})=0.$
(S24)

In the thermodynamic limit, the GBZ equation is determined by
$|\beta_{1}(\lambda)|=|\beta_{2}(\lambda)|$ YW18 ; Yokomizo2019 . Thus, we
obtain

$|\beta_{1}|=|\beta_{2}|=\sqrt{|\frac{e^{2\gamma}\sin(\frac{\theta_{1}}{2})\sin(\frac{\theta_{2}}{2})-\cos(\frac{\theta_{1}}{2})\cos(\frac{\theta_{2}}{2})}{\sin(\frac{\theta_{1}}{2})\sin(\frac{\theta_{2}}{2})-e^{2\gamma}\cos(\frac{\theta_{1}}{2})\cos(\frac{\theta_{2}}{2})}|}=\sqrt{|\frac{\cosh\gamma\cos\frac{\theta_{1}+\theta_{2}}{2}-\sinh\gamma\cos\frac{\theta_{2}-\theta_{1}}{2}}{\cosh\gamma\cos\frac{\theta_{1}+\theta_{2}}{2}+\sinh\gamma\cos\frac{\theta_{2}-\theta_{1}}{2}}|}.$
(S25)

Therefore, the GBZ is a circle in the complex plane, as shown in Fig. 2 in the
main text. When $|\beta|<1(|\beta|>1)$, the skin modes are localized at the
left (right) edge. According to Eq. (S25), when
$\cos\frac{\theta_{1}+\theta_{2}}{2}\cos\frac{\theta_{2}-\theta_{1}}{2}>0$
$(\cos\frac{\theta_{1}+\theta_{2}}{2}\cos\frac{\theta_{2}-\theta_{1}}{2}<0)$,
the skin modes are localized at the left (right) edge. For the two sets of
parameters used in the main article, the skin modes are localized at the
domain wall.

##  Appendix B Numerical fitting for larger time steps

In the experiment, we have found that the relative height can be well fitted
by $P_{\text{domain}}/P_{\text{min}}\sim x_{0}$. Thus, the relative height
grows as $x_{0}$ increases. In this section, we add numerical simulations with
more steps to further demonstrate this behavior.

As illustrated in Fig. S1, we fit the loss probability $P(x=-1)$ at the domain
wall, $P(x)$ in bulk, and the relative height
$P_{\text{domain}}/P_{\text{min}}$. The results show that when the edge burst
exists [Fig. S1(a,c,e)], both $P(x=-1)$ and $P(x)$ follow power laws:
$P(x=-1)\sim x_{0}^{-\alpha_{d}}$ and $P(x)\sim(x_{0}-x)^{-\alpha_{b}}$, with
certain $\alpha_{d}$ and $\alpha_{b}$. The fitting for the relative height is
$P_{\text{domain}}/P_{\text{min}}\sim x_{0}^{1.0802}$, which is close to the
$P_{\text{domain}}/P_{\text{min}}\sim x_{0}$ behavior predicted by theory and
supported by our experiment. Notably, the fitting for $\alpha_{d,b}$ are
$\alpha_{d}=0.4717$ and $\alpha_{b}=1.4751$, so that
$\alpha_{b}-\alpha_{d}=1.0034$, which agrees well with the predicted bulk-edge
scaling relation in Ref. XHS+22 .

When the edge burst is absent [Fig. S1(b,d,f)], the fitting turns out to be
exponential: $P(x=-1)\sim\beta_{d}^{x_{0}-x}$ and
$P(x)\sim\beta_{b}^{x_{0}-x}$, with certain $\beta_{b}$ and $\beta_{d}$ that
are approximately equal. The relative height
$P_{\text{domain}}/P_{\text{min}}$ is almost constant as $x_{0}$ varies.

Figure S1:  Numerical simulations for the loss probabilities $P(-1)$ (for the
domain wall) and $P(x)$ (for the bulk), and the relative height
$P_{\text{domain}}/P_{\text{min}}$ . The coin parameters are fixed as
$\theta_{1,2}^{L}=0.85\pi$ and $\theta_{1}^{R}=0.12\pi$. For the upper row
(red), $\theta_{2}^{R}=0.12\pi$, and the edge burst is present; for the lower
row (blue), $\theta_{2}^{R}=0.48\pi$, and the edge burst is absent. (a)(b) The
loss probability $P(x=-1)$ versus $x_{0}$. (c)(d) $P(x)$ versus $x_{0}-x$.
(e)(f) The relative height $P_{\text{domain}}/P_{\text{min}}$. The dots are
from numerical simulations, and the black solid lines are the fitting results.