# Boosting Graph Search with Attention Network for Solving the General
Orienteering Problem

Jing Xu State Grid Corporation of China<EMAIL_ADDRESS>Jintao Su Zhejiang
University<EMAIL_ADDRESS>Tao Xiao State Grid Corporation of China
<EMAIL_ADDRESS>Zongtao Liu1 Jing Xu2 Jintao Su1 Tao Xiao2&Yang Yang1
1Zhejiang University
2State Grid Corporation of China
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Recently, several studies have explored the use of neural network to solve
different routing problems, which is an auspicious direction. These studies
usually design an encoder-decoder based framework that uses encoder embeddings
of nodes and the problem-specific context to produce node sequence(path), and
further optimize the produced result on top by beam search. However, existing
models can only support node coordinates as input, ignore the self-referential
property of the studied routing problems, and lack the consideration about the
low reliability in the initial stage of node selection, thus are hard to be
applied in real-world.

In this paper, we take the orienteering problem as an example to tackle these
limitations. We propose a novel combination of a variant beam search algorithm
and a learned heuristic for solving the general orienteering problem. We
acquire the heuristic with an attention network that takes the distances among
nodes as input, and learn it via a reinforcement learning framework. The
empirical studies show that our method can surpass a wide range of baselines
and achieve results close to the optimal or highly specialized approach. Also,
our proposed framework can be easily applied to other routing problems. Our
code is publicly
available111https://anonymous.4open.science/repository/7cb20ede-b50e-4a9a-99a5-e3c3626bf1a1/.

##  1 Introduction

The orienteering problem(OP) is an important routing problem that originates
from the sports game of orienteering. In this game, each competitor starts at
a specific control point, visits control points in the tour, and will arrive
at a destination point (usually the same as the start point). Each control
point has an associated prize, and the travel between control points involves
a certain cost. The goal is to select a sequence of points such that the total
prize is maximized within the cost constraint. This problem relates to several
practical applications, such as resource delivery, tourist trip guide and
single-ring design in telecommunication networks(Golden et al. (1987);
Souffriau et al. (2008); Thomadsen and Stidsen (2003)). Golden et al. (1987)
have shown the OP to be an NP-hard problem, which motivates vigorous research
into the design of heuristic solvers and approximation algorithms.

A recent trend in the machine learning community for solving routing problems
is using deep neural networks to learn heuristic algorithms. With the help of
other meta-algorithm such as beam search, some of these methods can achieve
near start-of-the-art performance on several tasks (Li et al. (2018); Kool et
al. (2018); Vinyals et al. (2015); Khalil et al. (2017); Nazari et al.
(2018)). These tasks include the traveling salesman problem(TSP) and the
orienteering problem(OP), and the vehicle routing problem(VRP). These studies
deem solving such problems as a sequence generation process from the
coordinates of nodes, and use recurrent neural networks(RNNs) or graph neural
networks(GNNs) to build a bridge between node patterns and the policy. The OP
belongs to the routing problem family and thus is fit for this methodology.

However, current learning-based algorithms have several limitations. Firstly,
existing methods do not support direct input of travel cost information. They
implicitly assume the travel cost, saying travel distance here, between each
pair of points is their Euclidean distance, and take the coordinate of each
position as network input to indirectly acquire the distance information.
However, the distance metric used in real-world is not limited to the
Euclidean distance, thus such kind of methods has limited practicability. For
the term convenience, we here define the OP which only supports coordinates as
input as the coordinate orienteering problem, and define the OP that supports
direct input of costs among locations as the general orienteering problem.
Secondly, the routing problems considered in the previous works have self-
referential property, but these works lack insight into it. For instance, in
the orienteering problem, each step of point selection can be viewed as the
initial step in another distinct orienteering problem, in which the competitor
starts from the last selected points and will choose from the rest points.
Lacking consideration about such property in the existing studies might
degrade the utilization of samples when training.

Besides, most existing studies usually rely on a step-by-step beam search to
acquire a better solution (Kaempfer and Wolf (2018); Vinyals et al. (2015);
Nowak et al. (2017)); however, it might not work well on the routing problem
with relative large problem size. The step-by-step beam search at each step
prunes non-promising partial solutions based on a heuristic function and
stores a fix-sized set (also called beam) of best alternatives. However, as
the difficulty in the early stage of a routing problem is much more than that
in the later stage, the estimated heuristic score (probability or action
value) to perform beam search in the initial step might also be less credible.
To overcome this shortcoming, the stored partial solutions in the early stage
of a search procedure should be more than those in the later stage (since the
less credible heuristic scores might mislead the pruning). How to design such
a search algorithm is an important issue.

In this paper, we therefore propose a novel combination of a variant beam
search algorithm and the learned heuristic for solving the general
orienteering problem. We acquire the heuristic with an attention network that
takes both node attribute (prize) and edge attribute (cost) as input, and
learn it via a reinforcement learning framework. Rather than formulating such
routing tasks as sequence generation, we instead consider each step of point
selection separately; that is,we view each state in the decision process as
another different OP. This insight helps provide more direct reinforcement in
each state and improves the utilization of training instances. The
experimental result shows that our method can surpass a wide range of
baselines and achieve results close to the optimal or highly specialized
approach.

##  2 Related Work

Orienteering problem.  Over the years, many challenging applications in
logistics, tourism, and other fields were modeled as OP. Meanwhile, quite a
few studies for orienteering problems have been conducted. Golden et al.
(1987) prove that the OP is NP-hard, i.e., no polynomial-time algorithm has
been designed or is expected to be developed to solve this problem to
optimality. Thus the exact solution algorithms are very time consuming, and
for practical applications, heuristics will be necessary. Exact solution
methods (using branch-and-bound and branch-and-cut ideas) have been presented
by Laporte and Martello (1990); Ramesh et al. (1992). Heuristic approaches,
which use traditional operation research techniques, have been developed by
Tsiligirides (1984); Gendreau et al. (1998); Tang and Miller-Hooks (2005).

Applications of neural networks in routing problems.  Using neural networks
(NNs) for optimizing decisions in routing problems is a long-standing
direction, which can date back to the early work (Hopfield and Tank (1985);
Wang et al. (1995)). These studies usually design an NN and learn a solution
for each instance in an online fashion. Recent studies focus on using (D)NNs
to learn about an entire class of problem instances. The critical point of
these studies is to model permutation-invariant and variable-sized input
appropriately.

Earlier researches leverage sequence-to-sequence models and the attention
mechanism to produce path, and update network parameters via supervised or
reinforcement learning (Vinyals et al. (2015); Bello et al. (2016); Nazari et
al. (2018)). More recently, several studies leverage the power of GNNs that
can handle variable-sized and order-independent input to tackle the routing
problems (Khalil et al. (2017); Nowak et al. (2017); Kool et al. (2018)).
However, these studies only focus on coordinate routing problem; thus the
practicability is limited in real-word. Besides, current methods usually use a
set-to-sequence/node-to-sequence model to perform routing, which ignores the
self-referential property of many routing problems (including the OP).

##  3 Problem Definition

Formally, we define each orienteering problem instance $s$ as a tuple with six
elements $\langle V,C,R,T,v_{start},v_{end}\rangle$. More specifically,
$V=\\{v_{1},...,v_{M}\\}$ refers to the set of nodes (control points)
containing the start node $v_{start}$, the prized nodes, and the end node
$v_{end}$, where $M$ denotes the number of nodes. $C$ refers to the cost map
where $C_{v_{i},v_{j}}$ represents the travel cost from $v_{i}$ to $v_{j}$.
$R$ is the prize map where $R_{v}$ indicates the prize collected when $v$ is
visited. Therefore, the problem is to find a path from $v_{start}\in V$ to
$v_{end}\in V$ such that the total cost of the path is less than the
prescribed cost limit $T$, and the overall prize collected from the nodes
visited is maximized. The orienteering problem can also be formalized as a
mixture integer problem, which can be referenced in Gunawan et al. (2016).

For notation convenience, here we introduce some operators or functions among
these elements. $cost(P)$, $prize(P)$, and $last(P)$ represent the total cost,
the total prize, and the last selected point of the (partial) selected path
$P$, respectively. $(P+v_{i})$ represents a selected path obtained by adding a
new node $v_{i}$ to the original path $P$. $(V-P)$ represents the node set
that excludes the nodes in $P$ from $V$.

##  4 Method

###  4.1 Cost-Level Beam Search with the Learned Heuristic

Before introducing the search method proposed in this study, we first present
an exact search method. We divide the maximal cost $T$ into $\lceil
T/\tau\rceil$ intervals with the length of $\tau$. For each range, we maintain
a queue to save all partial paths(solutions) with the total cost in that
interval . Starting from interval $t=0$, we iteratively retrieve the
incomplete solution $P$ from the front of the queue and scan all the nodes
that are unselected. For each node $v$ scanned, if $cost(P+v+v_{end})$ is no
greater than the cost limit $T$, we then add the updated partial path $(P+v)$
to the queue in interval $(t+\lceil C_{last(P),v}/\tau\rceil)$. Finally, from
all stored paths with the end node, the path with the highest total prize is
the optimal path.

This method cannot finish computation in polynomial time. As the value of $t$
increases, the number of stored partial paths per window increases
exponentially. To reduce the space and time occupied by the search, at each
step of the iteration, it is practical to prune some partial paths based on a
predefined heuristic score $f(P,s)$ with the input of the current partial
solution $P$ and the problem instance $s$. To follow up this idea, we apply
the idea of beam search and replace the queue maintained in each window $t$ by
a priority queue sized $K$, which is used to save the paths with the K highest
heuristic scores in each window. We introduce a data-driven heuristic score
function $f$ to estimate the total prize that can be reached along the partial
path $P$:

$f(P,s)=prize(P)+e(s^{\prime};\theta)$ (1)

where $f$ consists of two terms, in which the term $prize(P)$ computes the
total prize of the given partial path, and the term $e(s^{\prime};\theta)$ is
an estimation of the subsequent prize along the current path under $\pi$ till
the end of the problem. $e(\cdot~{};\theta)$ is an estimation function
parameterized $\theta$ that computes the total potential prize given by an OP
under the policy $\pi_{\theta}$ and $s^{\prime}$ represents the subproblem
$\langle V-P,C,S,T-cost(P),last(P),v_{end}\rangle$ of the original OP.

A question that naturally arises here is how to obtain a reliable function
$e(\cdot~{};\theta)$ to estimate the total prize of an orienteering problem.
In the next section, we will present our solution with the attention-based
neural networks learned by Q-learning.

Algorithm 1 Cost-Level Beam Search

1: initialize a list of priority queues PQ[$\lceil T/\tau\rceil$]

2: initialize the empty path $P_{g}$ and insert it to PQ[0]

3: for t=0 to $\lceil T/\tau\rceil$ do

4: while PQ[t] is not empty do

5: pop partial path $P_{c}$ from PQ[t]

6: for v in $(V-P_{c})$ do

7: pq = PQ[$\lceil cost(P_{c}+v)/\tau\rceil$]

8: if $cost(P_{c}+v+v_{end})>=T$ then continue

9: insert $(P_{c}+v)$ to pq

10: if $prize(P_{c}+v)>prize(P_{g})$ then $P_{g}=(P_{c}+p)$

11: if the size of pq is larger than K then

12: pop $argmin_{P}f(P,s)$ from pq

13: end if

14: end for

15: end while

16: end for

17: return path $P_{g}$

###  4.2 Prize Estimation via Attention Network

Kool et al. (2018) propose a network architecture based on self-attention to
model several coordinate routing problems, including TSP, OP, and several VRP
variants. With the benefit of reinforcement learning, the learned heuristic
achieves better performance over other learning-based heuristic methods.
Following the idea of using the attention mechanism to model the node
interactions with permutation-invariant property, we propose an attention-
based network, named Attention Network (AN). Our proposed structure is not
limited in handling coordinate OP, but can directly take all the node
attributes (prize of each node), edge attributes (costs among nodes) and
global attribute (the remaining cost of the agent) as input.

Figure 1: An illustration of our Attention Network.

Edge-Aware Graph Attention Networks Layer.  The graph attention
networks(GAT)Veličković et al. (2017) utilize the attention mechanism to
aggregate the attributes of node neighbors adaptively, and extract high-level
representation feature vector of each node. We modify the original structure
of GAT; thus it can aggregate information from the connected edges, as well as
the node neighbors and global information:

$\displaystyle\alpha_{vk}=\frac{\exp\left(\text{LeakyReLU}\left(\mathbf{a}^{T}\left[{\mathbf{W}}[\mathbf{x}_{v}\|\mathbf{x}_{k}\|\mathbf{u}_{vk}\|\mathbf{g}]\right]\right)\right)}{\sum_{j\in\mathcal{N}_{v}}\exp\left(\text{LeakyReLU}\left(\mathbf{a}^{T}\left[{\mathbf{W}}[\mathbf{x}_{v}\|\mathbf{x}_{j}\|\mathbf{u}_{vj}\|\mathbf{g}]\right]\right)\right)}$
(2)
$\displaystyle\mathbf{h}_{\mathcal{N}_{v}}=\sigma\left(\sum_{k\in\mathcal{N}_{v}}\alpha_{vk}\left[{\mathbf{W}}[\mathbf{x}_{v}\|\mathbf{x}_{k}\|\mathbf{u}_{vk}\|\mathbf{g}]\right]\right)$
(3) $\displaystyle
MHA:\mathbf{h}_{\mathcal{N}_{v}}=\sigma\left(\frac{1}{M}\sum_{m=1}^{M}\sum_{k\in\mathcal{N}_{v}}\alpha_{vk}^{m}\left[{\mathbf{W}^{m}}[\mathbf{x}_{v}\|\mathbf{x}_{k}\|\mathbf{u}_{vk}\|\mathbf{g}]\right]\right)$
(4)

where $x_{v}$ denotes $v$’s node attribute (prize $R_{v}$), $u_{vw}$ denotes
the edge attribute (travel cost $C_{v,w}$) between node $v$ and $w$, and $g$
indicates the problem-level feature (the maximal cost $T$). After aggregating
information from connected edges and node neighbors, we obtain the
intermediate representation of each node $H=\\{h_{v_{1}},h_{v_{i}},...,\\\
h_{v_{M}}\\}$, where $h_{v}$ is the representation of $v$.

Transformer Encoding Layer(TEL).  Recently, Kool et al. (2018) also
demonstrates the effectiveness of the Transformer to model many routing
problems dealing with coordinates of points. Following this idea, we use the
encoding layer of the Tranformer in this study to extract the node features
for the intermediate representation set
$H=\\{h_{v_{1}},h_{v_{2}},...,h_{v_{M}}\\}$. Since the resulting node
representations are invariant to the input order, we do not use a positional
encoding. Feeding the immediate embeddings $H$ into $N$ stacked Tranformer
encoding layers, we obtain the updated embeddings
$H^{\prime}=\\{h^{\prime}_{v_{1}},h^{\prime}_{v_{2}},...,h^{\prime}_{v_{M}}\\}$.

Linear Projection Layer.  Finally, we use the linear combination of $v_{i}$’s
embedding $h^{\prime}_{i}$ to estimate the potential total prize $r_{i}$ when
selecting $v_{i}$ in the next step. where $U\in\mathbf{R}^{|H|}$ is a
learnable weight vector. Thus, we obtain the prize estimation
$e(s;\theta)=max_{v\in V}q(s,v;\theta)$. Although it would be more convenient
here to estimate the total prize of the given OP directly, we adopt this way
for the sake of facilitating Q-learning to train the networks.

###  4.3 Training via Q-learning

We use a fitted double Q-learning algorithm (Van Hasselt et al. (2016)) to
obtain the potential future prize of each selection. We formulate the states,
actions, rewards, transitions, and policy in the reinforcement learning
framework as follows:

  * •

States: We define the state $s$ as the current subproblem of the original OP,
i.e., $\langle V-P,C,S,T-cost(P),last(P),v_{end}\rangle$.

  * •

Actions: The action $v$ is a node selection of $V$ that is in the part of the
current state $s$. Both the definition of states and actions is applicable
across node set with various size.

  * •

Rewards: We define the reward $r$ of each action $v$ as the prize collected
when the node is visited.

  * •

Transitions: A transition is a tuple $\langle
s_{t},v_{t},r_{t+1},s_{t+1}\rangle$ which represents the agent takes action
$v_{t}$ in state $s_{t}$ and then observes the immediate reward $r_{t+1}$ and
resulting state $s_{t+1}$.

  * •

Policy: Given the current state $s$, we can compute the q-value
$q(s,v;\theta)$ of each unselected node $v$ and apply a deterministic greedy
policy $\pi(v|s)=\operatorname*{\arg\\!\max}_{v\in|V|}q(s,v;\theta)$.

We use double Q-learning (Van Hasselt et al. (2016)) to learn a greedy policy
$\pi_{\theta}$ parametrized by the attention network and obtain the expected
total prizes $q(s,v;\theta)$ of each action $v$ simultaneously. The network
parameters can be updated by minimizing the following loss function of the
transitions picked up from the replayed memory. More details can be referenced
in supplementary part.

###  4.4 Discussions

Comparison with the step-by-step beam search.  The learned heuristic score
function $f$ of (1) can also be applied to the step-by-step beam search.
However, as the routing in the early stage is much complicated than that in
the later stage, the estimated heuristic score in the early stage might also
be less reliable. In the initial stage of step-by-step beam search, the
proportion of reliable selection could be very low and thus degrade the
routing performance. Compared with the step-by-step beam search, the cost-
level beam search store the early partial paths as many as possible.

Running times.  As the computation of the estimated action values in each time
window is parallelizable, the presented beam search algorithm can be
significantly accelerated. The basic idea is to compute the action value on
the GPU or TPU.

Extensibility to other routing problems.  By simple modification, our
framework can be easy to be extended to other routing problems. Here we take
the classical TSP as an example. Similar to OP, we set the learned heuristic
score function as $f(P,s)=-cost(P)-e(s^{\prime},\theta)$, which represents the
negative of the expected total cost of the whole tour, while the edge
attribute is placed the same as those in the OP (the node and global
attributes are not needed). The reward in each selection is defined as the
negative of the increased cost. The cost constraint $T$ can be set as the
total travel cost given by another simple heuristic method.

##  5 Experimental Results

Table 1: Performance comparison. The gap % is w.r.t. the best value across all methods. |  | $n=20,T=2$ | $n=50,T=3$ | $n=100,T=4$
---|---|---|---|---
Prize | Method | Obj. | Gap | Time | Obj. | Gap | Time | Obj. | Gap | Time
Uniform | Gurobi | 5.85 | $0.00\%$ | (1s) | - | -
Compass | $5.84$ | $0.17\%$ | (0s) | 16.46 | $0.00\%$ | (2s) | 33.30 | $0.00\%$ | (6s)
Random | $2.13$ | $63.59\%$ | (0s) | $3.33$ | $79.77\%$ | (0s) | $4.37$ | $86.88\%$ | (0s)
Tsili (greedy) | $4.85$ | $17.09\%$ | (0s) | $12.46$ | $22.94\%$ | (0s) | $25.48$ | $23.48\%$ | (0s)
AN-dqn (greedy) | ${5.15}$ | ${11.97\%}$ | (0s) | ${13.91}$ | ${15.49\%}$ | (0s) | ${26.51}$ | ${20.39\%}$ | (1s)
AN-a2c (greedy) | ${4.97}$ | ${15.04\%}$ | (0s) | ${14.13}$ | ${14.16\%}$ | (0s) | ${24.31}$ | ${27.00\%}$ | (1s)
AN-at (greedy) | ${3.98}$ | ${31.97\%}$ | (0s) | ${8.22}$ | ${50.06\%}$ | (0s) | ${9.75}$ | ${70.72\%}$ | (1s)
AN-pn (greedy) | ${4.46}$ | ${23.76\%}$ | (0s) | ${11.10}$ | ${32.56\%}$ | (1s) | ${10.06}$ | ${69.79\%}$ | (2s)
Tsili (beam) | $5.54$ | $5.30\%$ | (0s) | $13.54$ | $15.63\%$ | (2s) | $25.91$ | $22.19\%$ | (9s)
AN-dqn (beam) | ${5.74}$ | ${18.80\%}$ | (1s) | ${15.76}$ | ${4.25\%}$ | (10s) | ${30.61}$ | ${8.07\%}$ | (40s)
AN-a2c (beam) | ${5.57}$ | ${4.79\%}$ | (1s) | ${15.13}$ | ${8.08\%}$ | (5s) | ${27.40}$ | ${17.72\%}$ | (23s)
AN-at (beam) | ${4.45}$ | ${23.93\%}$ | (1s) | ${9.45}$ | ${16.46\%}$ | (4s) | ${10.55}$ | ${68.32\%}$ | (20s)
AN-pn (beam) | ${4.84}$ | ${17.26\%}$ | (1s) | ${13.75}$ | ${42.59\%}$ | (4s) | ${12.08}$ | ${63.72\%}$ | (24s)
CS (20, 0.05) | ${5.82}$ | ${0.51\%}$ | (1s) | ${16.13}$ | ${2.00\%}$ | (7s) | ${30.93}$ | ${7.12\%}$ | (40s)
Distance | Gurobi | 5.39 | $0.00\%$ | (4s) | - | -
Compass | $5.37$ | $0.37\%$ | (0s) | 16.17 | $0.00\%$ | (3s) | 33.19 | $0.00\%$ | (8s)
Random | $1.89$ | $64.94\%$ | (0s) | $2.90$ | $82.07\%$ | (0s) | $3.81$ | $88.52\%$ | (0s)
Tsili (greedy) | $4.08$ | $24.30\%$ | (0s) | $12.46$ | $22.94\%$ | (0s) | $25.69$ | $22.60\%$ | (0s)
AN-dqn (greedy) | ${4.77}$ | ${11.50\%}$ | (0s) | ${13.70}$ | ${15.28\%}$ | (0s) | ${26.60}$ | ${19.86\%}$ | (1s)
AN-a2c (greedy) | ${4.57}$ | ${15.21\%}$ | (0s) | ${13.70}$ | ${15.28\%}$ | (0s) | ${22.32}$ | ${32.75\%}$ | (1s)
AN-at (greedy) | ${3.62}$ | ${32.84\%}$ | (0s) | ${9.29}$ | ${42.55\%}$ | (1s) | ${10.52}$ | ${68.30\%}$ | (2s)
AN-pn (greedy) | ${3.75}$ | ${30.43\%}$ | (0s) | ${11.37}$ | ${29.68\%}$ | (1s) | ${10.77}$ | ${67.55\%}$ | (2s)
Tsili (beam) | $5.26$ | $1.68\%$ | (0s) | $15.50$ | $14.50\%$ | (2s) | $30.53$ | $8.04\%$ | (10s)
AN-dqn (beam) | ${5.27}$ | ${11.5\%}$ | (1s) | ${15.82}$ | ${4.14\%}$ | (9s) | ${30.28}$ | ${8.77\%}$ | (35s)
AN-a2c (beam) | ${5.02}$ | ${6.86\%}$ | (1s) | ${15.13}$ | ${6.43\%}$ | (6s) | ${27.83}$ | ${16.15\%}$ | (24s)
AN-at (beam) | ${4.42}$ | ${18.00\%}$ | (1s) | ${10.68}$ | ${33.95\%}$ | (4s) | ${11.10}$ | ${66.56\%}$ | (21s)
AN-pn (beam) | ${4.24}$ | ${21.34\%}$ | (1s) | ${12.15}$ | ${24.86\%}$ | (4s) | ${12.48}$ | ${62.40\%}$ | (17s)
CS (20, 0.05) | ${5.35}$ | ${0.74\%}$ | (1s) | ${15.99}$ | ${1.11\%}$ | (6s) | ${31.16}$ | ${6.12\%}$ | (42s)

###  5.1 Dataset Creation

To evaluate our proposed method against other algorithms, we generate problem
instances with different settings. We sample each node $v$’s coordinates
$x_{v}$ uniformly at random in the unit square and then compute the travel
cost from node $v$ to node $w$ by their Euclidean distance $dis(v,w)$, where
$dis(v,w)$ denotes the Euclidean distance between node $v$ and $w$. We set the
distance metrics as Euclidean distance because although many OP algorithms can
handle the OP with different distance metrics, most released implementations
can only support planar coordinates as input for the sake of convenience. One
more to mention is that selecting the Euclidean distance metric does not mean
that we will compare the performance of our approach with those specialized
methods designed for coordinate OP like Khalil et al. (2017); Kool et al.
(2018), because such comparison is unfair.

For the prize distribution, we consider three different settings described by
Fischetti et al. (1998); Kool et al. (2018). Besides, we normalize the prizes
to make the normalized value fall between 0 and 1. Uniform.
$r^{\prime}_{v}\thicksim DiscreteUniform(1,100);r_{v}=r^{\prime}_{v}/100$. The
prize of each node is discretized uniformed.

Distance.  $r_{v}=(1+\lfloor 99\cdot\frac{dis(x_{v_{start}},x_{v})}{max_{v\in
V}dis(x_{v_{start}},x_{v})}\rfloor)/100$. Each node has a (discretized) prize
that is proportional to the distance to the start node, which is a challenging
task as the node with the largest prize are the furthest away from the start
nodeFischetti et al. (1998).

We choose the maximal travel cost $T$ as a value close to half of the average
optimal length of the corresponding TSP tour. Setting approximately the half
of the nodes can be visited results in the most challenging problem instances
(Vansteenwegen et al. (2011)). Finally, we set the value of $T$ for the OP
sized 20, 50, 100 as 2, 3, and 4, respectively.

###  5.2 Performance Comparison

Comparable Methods. We evaluate the performance of our method under the above-
mentioned instance generation settings . Firstly, we compare our proposed
method with an exact method Gurobi (Optimization (2014)) and a state-of-the-
art heuristic method Compass (Kobeaga et al. (2018)).We further compare with
several policy-based methods, which compute the node selection probability
$p(v|s)$ or action value $q(s,v)$ in state $s$.These methods include (1)Random
selection (Random), (2)Tsili (Tsiligirides (1984)), (3)policy learned by deep
Q-learning algorithm (AN-dqn), (4)policy learned by advantage actor-critic
algorithm (AN-a2c). We report their results based on two meta-algorithm:
(1)greedy: the node with the highest probability or action value is selected
in each state. (2) beam: we apply a step-by-step beam search for some of the
below methods to obtain the target path. The beam size is set as 100.

(a) With/without MHA in GAT. (b) Varing #layer of TLE. (c) Varing hidden size.

Figure 2: Model parameter/structure analysis. The x-axis is the number of
episodes. The y-axis denotes the average total prize on 200 test instances
with greedy policy using the estimated action values from AN. This experiment
is conducted five times, and we show the average performance and the
corresponding standard error.

In addition to the above methods as benchmark approaches, the comparison
between $AN-a2c$, $AN-at$, and $AN-pn$ also serves to demonstrate the
effectiveness of our network architecture compared with the encoder-decoder
schemes used in the previous works. This is because $AN-a2c$ can be view as a
policy gradient version of our proposed DQN framework proposed in Section 4.3,
which considers each step of node selection as a distinct OP. In contract,
$AN-at$ and $AN-pn$ view each selection is related to the previous ones (via
RNN or self-attention mechanism). Lastly, we report the performance of our
cost-level beam search algorithm CS(K, $\tau$). $K$ and $\tau$ represent the
beam size and cost interval, respectively.

Comparison Results.  Table 1 reports the performance and average running time
of each approach on 10K test instances. Presenting running time is not for the
comparison between the efficiency of different methods. The results of
different approaches might not be comparable due to the differences in their
hardware usages (e.g., CPU or GPU), details of implementation (e.g.,
C++/Python with/without some optimizations), the settings of hyperparameters
(e.g., beam size). We report the time spent mainly to show that in general,
our method and implementation are time-affordable with the help of GPU
acceleration.

We note that our search method with the learned heuristic ($CS$) surpasses all
comparative methods except the exact approach Gurobi and the specialized
genetic algorithm Compass (the gap is relatively small compared with those of
other benchmarks), which demonstrates the effectiveness of our proposed
methods. AN-a2c is a competitive method concerning other baselines with
attention network as a shared bottom. However, since the step-by-step beam
search retains relative fewer partial paths at the start of the selection
process, our search method significantly outperforms this variant.

The comparison between $AN-dqn(beam)$ and $CS$ demonstrates the effectiveness
of the cost-level beam search.

By comparison of the result of AN-a2c, AN-at, and AN-pn in a greedy fashion or
a beam search fashion, we can find which scheme of network design works better
for the OP. We find that An-a2c significantly outperforms the rest two
approaches. This shows it is necessary to take into account the self-
referential property of the OP and use a non-sequential network to model this
problem.

###  5.3 Parameter/Structure Analysis

We conduct sensitive analysis in this section. We present the performance of
AN-dqn (greedy policy using DQN) on 200 instances in the learning process.
Each instance has 50 prized nodes, and the prize of each node is proportional
to its distance to the start node (distance). We note that by introducing
multi-head, the performance increases, which shows that the multi-head helps
the GAT better aggregate edge information. Following the GAT, the Transformer
encoding layer(TEL) further extracts more high-level representation for nodes.
From Figure 2(b), we observe that the performance improves as the number of
layer increases. Lastly, we explore the relationship between the hidden size
$\mathbb{R}^{H}$ and the quality of the learned model.

###  5.4 Visualization

Finally, we visualize the output paths generated by $AN-dqn$, $CS-p$, and
$CS-l$ with the problem size of 50 and the prize type of distance in Figure 3.
We observe that compared with the results of the greedy selection with action
values ($AN-dqn$) or the cost-level beam search with a simple heuristic
function ($CS-p$), the combination of our cost-level beam search and the
learned heuristic ($CS-l$) can conduct more look-forwarding selection, thus
have better performance.

(a) result of AN-dqn on Instance 1 (b) result of AN-dqn on Instance 2

Figure 3: Visualization of the output paths by AN-dqn, and CS. The star marker
denotes the start/end node, and the circle marker represents the prized node,
the area of which is proportional to its prize.

##  6 Conclusion and Future Work

In this paper, we solve the NP-hard problem of the general orienteering
problem by conducting a novel combination of a cost-level beam search method
and the learned heuristic score function. To estimate the potential prize of
the current partial solution, we design an attention-based network. It first
uses a variant GAT to aggregate both the edge attribute(travel cost) and the
node attribute(prize), then applies stacked Transformer encoding layers to
produce high-level node representations. Rather than modeling the OP in an
encoder-decoder structure, we use the extracted node representations to
directly estimate the target action values, which is a better scheme for the
OP with the self-referential property. We further conduct sufficient
experiments to show that our method surpasses most of the benchmark methods
and get results close to those of the exact method and the highly optimized
and specialized method. We also note that the performance of our method has
been improved by introducing the deep neural network. Besides, with a simple
modification, our proposed framework can also be straightforward to apply to
other routing problems.

## References

  * Bello et al. [2016] Irwan Bello, Hieu Pham, Quoc V Le, Mohammad Norouzi, and Samy Bengio.  Neural combinatorial optimization with reinforcement learning.  arXiv preprint arXiv:1611.09940, 2016.
  * Fischetti et al. [1998] Matteo Fischetti, Juan Jose Salazar Gonzalez, and Paolo Toth.  Solving the orienteering problem through branch-and-cut.  INFORMS Journal on Computing, 10(2):133–148, 1998.
  * Gendreau et al. [1998] Michel Gendreau, Gilbert Laporte, and Frédéric Semet.  A tabu search heuristic for the undirected selective travelling salesman problem.  European Journal of Operational Research, 106(2-3):539–545, 1998\.
  * Golden et al. [1987] Bruce L Golden, Larry Levy, and Rakesh Vohra.  The orienteering problem.  Naval Research Logistics (NRL), 34(3):307–318, 1987.
  * Gunawan et al. [2016] Aldy Gunawan, Hoong Chuin Lau, and Pieter Vansteenwegen.  Orienteering problem: A survey of recent variants, solution approaches and applications.  European Journal of Operational Research, 255(2):315–332, 2016\.
  * Hopfield and Tank [1985] John J Hopfield and David W Tank.  “neural” computation of decisions in optimization problems.  Biological cybernetics, 52(3):141–152, 1985.
  * Kaempfer and Wolf [2018] Yoav Kaempfer and Lior Wolf.  Learning the multiple traveling salesmen problem with permutation invariant pooling networks.  arXiv preprint arXiv:1803.09621, 2018.
  * Khalil et al. [2017] Elias Khalil, Hanjun Dai, Yuyu Zhang, Bistra Dilkina, and Le Song.  Learning combinatorial optimization algorithms over graphs.  In NIPS, pages 6348–6358, 2017.
  * Kobeaga et al. [2018] Gorka Kobeaga, María Merino, and Jose A Lozano.  An efficient evolutionary algorithm for the orienteering problem.  Computers & Operations Research, 90:42–59, 2018.
  * Kool et al. [2018] Wouter Kool, Herke van Hoof, and Max Welling.  Attention, learn to solve routing problems!  arXiv preprint arXiv:1803.08475, 2018.
  * Laporte and Martello [1990] Gilbert Laporte and Silvano Martello.  The selective travelling salesman problem.  Discrete applied mathematics, 26(2-3):193–207, 1990.
  * Li et al. [2018] Zhuwen Li, Qifeng Chen, and Vladlen Koltun.  Combinatorial optimization with graph convolutional networks and guided tree search.  In NIPS, pages 539–548, 2018.
  * Nazari et al. [2018] Mohammadreza Nazari, Afshin Oroojlooy, Lawrence Snyder, and Martin Takác.  Reinforcement learning for solving the vehicle routing problem.  In NIPS, pages 9839–9849, 2018.
  * Nowak et al. [2017] Alex Nowak, Soledad Villar, Afonso S Bandeira, and Joan Bruna.  A note on learning algorithms for quadratic assignment with graph neural networks.  stat, 1050:22, 2017.
  * Optimization [2014] Gurobi Optimization.  Inc.,“gurobi optimizer reference manual,” 2015, 2014.
  * Ramesh et al. [1992] R Ramesh, Yong-Seok Yoon, and Mark H Karwan.  An optimal algorithm for the orienteering tour problem.  ORSA Journal on Computing, 4(2):155–165, 1992.
  * Souffriau et al. [2008] Wouter Souffriau, Pieter Vansteenwegen, Joris Vertommen, Greet Vanden Berghe, and Dirk Van Oudheusden.  A personalized tourist trip design algorithm for mobile tourist guides.  Applied Artificial Intelligence, 22(10):964–985, 2008.
  * Tang and Miller-Hooks [2005] Hao Tang and Elise Miller-Hooks.  A tabu search heuristic for the team orienteering problem.  Computers & Operations Research, 32(6):1379–1407, 2005.
  * Thomadsen and Stidsen [2003] Tommy Thomadsen and Thomas K Stidsen.  The quadratic selective travelling salesman problem.  Informatics and Mathematical Modeling Technical Report, 2003.
  * Tsiligirides [1984] Theodore Tsiligirides.  Heuristic methods applied to orienteering.  Journal of the Operational Research Society, 35(9):797–809, 1984\.
  * Van Hasselt et al. [2016] Hado Van Hasselt, Arthur Guez, and David Silver.  Deep reinforcement learning with double q-learning.  In AAAI, 2016.
  * Vansteenwegen et al. [2011] Pieter Vansteenwegen, Wouter Souffriau, and Dirk Van Oudheusden.  The orienteering problem: A survey.  European Journal of Operational Research, 209(1):1–10, 2011.
  * Veličković et al. [2017] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * Vinyals et al. [2015] Oriol Vinyals, Meire Fortunato, and Navdeep Jaitly.  Pointer networks.  In NIPS, pages 2692–2700, 2015.
  * Wang et al. [1995] Qiwen Wang, Xiaoyun Sun, Bruce L Golden, and Jiyou Jia.  Using artificial neural networks to solve the orienteering problem.  Annals of Operations Research, 61(1):111–120, 1995.


# Analysis of Many-body Localization Landscapes and Fock Space Morphology via
Persistent Homology

Gregory A. Hamilton<EMAIL_ADDRESS>Bryan K. Clark  Institute for Condensed
Matter Theory and IQUIST and Department of Physics, University of Illinois
Urbana-Champaign, Urbana, IL 61801, USA

###### Abstract

We analyze functionals that characterize the distribution of eigenstates in
Fock space through a tool derived from algebraic topology: persistent
homology. Drawing on recent generalizations of the localization landscape
applicable to mid-spectrum eigenstates, we introduce several novel persistent
homology observables in the context of many-body localization that exhibit
transitional behavior near the critical point. We demonstrate that the
persistent homology approach to localization landscapes and, in general,
functionals on the Fock space lattice offer insights into the structure of
eigenstates unobtainable by traditional means.

##  I Introduction

Many-body localization (MBL) is a paradigmatic example of a dynamical phase
transition whereby eigenstates fail to thermalize [1; 2; 3]. MBL is
characterized by a multitude of observables including vanishing conductivity
[4; 5], logarithmic entanglement growth [6], area law entanglement [7], level
statistics [8], quasilocal integrals of motion [9; 10; 2], and a litany of
others [11; 1]. A variety of models feature in the MBL literature, including
the disordered model considered here, quasiperiodic potentials [12], and
quadratic Stark potentials [13]. While microscopic mechanisms for the phase
transition have been proposed and elucidated upon [14; 15; 16; 17], much
controversy still surrounds the stability of the MBL phase in the
thermodynamic limit [18; 19] as well as the structure of eigenstates near the
transition [20].

In response to this challenge a large body of recent work probes the
correlational structure of eigenstates in Fock space (here taken to be the
configurational $\sigma_{z}$ basis). The guiding intuition for this approach
is, as the many-body Hamiltonian maps to a tight binding model (TBM) on a Fock
space lattice (or graph), eigenstate spatial correlations on the lattice
dictate the behavior of autocorrelation functions that diagnose the phase
transition [21].

Fock space eigenstate structure has been a fundamental question since the
early days of Anderson localization (1BL), as the inverse participation ratio
(IPR) and participation entropies $S_{i}^{P}$ serve as indicators for the
Anderson transition [11; 22; 23; 3]. From the participation entropy stems a
notion of multifractality with multifractal exponents identifying the
universality class; an analogous approach to the participation entropy exists
in the many-body case.

In pursuit of a more nuanced understanding of the transition and MBL
stability, recent approaches look beyond the participation entropies to
quantify the “anatomy” of Fock space correlations [24; 21; 25; 17]. From this
TBM picture several phenomenological models for the MBL transition have been
proposed, including so-called “avalanche” mechanisms [26; 27; 28], and,
pertinent to this work, percolation models [24; 29; 30; 21; 25; 31; 32; 33].
We approach this Fock space viewpoint from a morphological lens by leveraging
a new functional: the localization landscape.

New insights into the 1BL problem hinge on a functional called the
localization landscape (LL) [34], which bounds remarkably well the spatial
extent of low-energy eigenstates [34; 35; 36]. The LL is the solution to a
simple differential equation:

$\displaystyle H\bm{u}=\bm{1},$ (1)

for Hamiltonian $H$ that satisfies mild constraints ($\bm{1}$ is the vector of
all ones). The LL (or rather, its inverse $V_{u}=1/\bm{u}$) operates as an
“effective potential”, and can be used in place of the physical potential to
well approximate both the integrated density of states via the classical Weyl
law [36; 37; 38] and localization lengths for low-energy eigenstates [38]. The
effective potential, being numerically tractable, has already been used
extensively in the context of quantum transport in disordered semiconductors
[39; 37; 36].

Connections between percolation theory and transport extend back several
decades [40], and the effective potential $V_{u}$ lends itself to an
interpretation wherein a transition from localized to extended eigenstates is
percolative, as has been recently hypothesized [39; 41]. This core idea, that
the morphological shift in an effective potential (or more generally, a
potential in a configuration or phase space) might serve as an indicator of a
phase transition, has been well-developed in the context of classical phase
transitions and serves as the impetus for our work here [42; 43; 44; 45].

As we describe in Sec. II, the LL is perfectly well-defined for the many-body
case by the TBM mapping [41], and variants of the LL equally apply to interior
eigenstates [46; 47]. The localization landscape variant we use here
(discussed in detail in Sec. II) is termed the $L_{2}$ landscape, and takes
the form

$\displaystyle(u^{(2)})^{2}=\text{diag}(M^{-1}),\,M:=H^{\dagger}H,$ (2)

where $H$ is the Hamiltonian [46]. As we discuss in Sec. II, with a proper
choice of energy $E_{0}$ the $L_{2}$ landscape bounds mid-spectrum eigenstates
just as the LL bounds low-energy eigenstates. In what follows we analyze the
$L_{2}$ landscape of a 1D spinless fermion model with open boundary
conditions, given by

$\displaystyle H=H_{t}+H_{W}+H_{V}$ (3) $\displaystyle H_{t}=t\sum_{\langle
ij\rangle}\left(c_{i}^{\dagger}c_{j}+\text{ H.c.}\right),$ $\displaystyle
H_{W}=\sum_{i}h_{i}\left(n_{i}-\frac{1}{2}\right),$ $\displaystyle
H_{V}=V\sum_{\langle
ij\rangle}\left(n_{i}-\frac{1}{2}\right)\left(n_{j}-\frac{1}{2}\right),$

with $P(h_{i})\in[-W,W]$ uniform and $\langle ij\rangle$ denoting nearest
neighbor (NN) pairs. The dimensionless parameter $V/t$ sets the strength of NN
interactions; for a choice of $(t,V)=(1/2,1)$ this model reduces to the
standard XXZ model. The $L_{2}$ landscape, while simple to define, is a
complex functional on the Fock space lattice. To ascertain what the $L_{2}$
landscape can tell us about the MBL transition, we leverage a new tool from
computational geometry and topological data analysis: persistent homology.

Persistent homology offers a computable handle on the small and large-scale
structure of a point cloud or underlying functional by building a filtration
of topological objects, usually simplicial complexes. Persistent homology is
largely parameter-free and yields novel observables fundamentally tied to
discrete Morse theory and algebraic topology [48]. While persistent homology
has recently been used to identify classical and quantum phase transitions
[49; 50; 51], we leverage it here to gain a deeper understanding of whether
the $L_{2}$ landscape exhibits signs of a percolation transition.

To this end we generate a large set of observables that quantify the
morphology structure of the LL with respect to the high-dimensional Fock
space. The observables and topological summaries we garner from this approach
paint a complex picture. While some results are easy to interpret in the
context of standard observables, still others defy simple interpretations.
Part of this difficulty lies in the incoincidence of persistent homology with
standard statistical techniques, an impediment we confront several times in
this work.

Nonetheless, we do find indications that persistent homology accesses the
fundamental structure of mid-spectrum eigenstates and coincides with critical
points found in previous works. We provide novel ways to probe the “anatomy”
of Fock space at the morphological level. In Sec. IV we describe many more
ways future research can extend and refine our approach.

To briefly state the novel contributions of our work:

  * •

We describe a new tool, persistent homology, to probe the morphology of Fock
space. In contrast to recent investigations of Fock space correlations [33],
this technique is parameter-free and probes clustering at all scales.

  * •

We identify the proliferation of extremal points on the $L_{2}$ as an
indicator for the MBL transition.

  * •

We demonstrate several topological observable that probe the clustering
properties of the $L_{2}$ landscape scale with the Hilbert space dimension.
The scaling exponents exhibit cross-over behavior near the MBL transition.

  * •

We give a new, persistent homology-based definition for fractal dimensions
applicable to both many-body eigenstates and functionals on Fock space.

The outline of this exploratory work is as follows. In Sec. II we describe the
localization landscape, its extension to the many-body case, and the variant
on the LL utilized in this work. We also give an introduction to persistent
homology and its primary output, the persistence diagram. In Sec. III we apply
our persistent homology pipeline to the Hamiltonian in Eq. 3. We give a
phenomenological understanding of several observables extracted from the
persistent homology pipeline, and assess the implications for percolative MBL
phase transition hypothesis. Finally, in Sec. IV we summarize our work and
suggest future directions of study.

Figure 1: Pictorial representation of a superlevel set filtration for a one-
dimensional function. The dashed line denote different level sets, while the
black lines underneath each box denote the intervals in the superlevel set. As
the level decreases, intervals emerge and merge until a final connected
components is formed. The bottom plot shows the corresponding persistence
diagram. The rightmost point $(a,-\infty)$ is displayed but is ignored with
respect to topological summaries of the persistence diagram. The green
dendrogram denotes the merge tree derived from the subset filtration. External
vertices correspond to the maxima, while internal vertices correspond to the
minima.

##  II Theory

In this section we overview the LL and $L_{2}$ landscape. We motivate and
describe how to apply persistent homology to Fock space functionals before
describing topological summaries of Fock space morphology extracted from
persistence diagrams and hierarchical clustering.

###  II.1 Localization landscapes

The original LL construction takes as input a single-body Hamiltonian $H$ with
eigenspectrum $\\{(E^{\beta},\ket{\phi^{\beta}})\\}$, expressed in orthonormal
basis $\\{\ket{I}\\}$ such that $H$ is positive-definite and $H_{IJ}<0$ for
all $I\neq J$. For Eq. 3 these conditions are easily achieved by adding an
offset to $H\mapsto-H+\Gamma I$ to ensure positivity (though care should be
taken regarding the choice of $\Gamma$ [52]). Then the solution $u$ to the
equation

$\displaystyle Hu=1$ (4)

is component-wise positive. Alternatively, we can express $u$ as [46]

$\displaystyle
u_{J}=\sum_{\beta}\frac{\phi_{J}^{\beta}}{E^{\beta}}\sum_{I}\phi_{I}^{\beta},$
(5)

where $\phi_{J}^{\beta}=\braket{J}{\phi^{\beta}}$.

Remarkably, $u$ bounds eigenstate amplitudes, as we have

$\displaystyle|\phi_{J}^{\beta}|\leq|E^{\beta}|||\phi^{\beta}||_{\infty}u_{J},$
(6)

where $||\phi||^{p}:=\left(\sum_{I}|\phi_{I}|^{p}\right)^{1/p}$ denotes the
$\ell_{p}$ norm. The inverse $V_{u}:=1/u$ serves as an effective potential,
the valleys of which indicate where low-energy eigenstates tend to localize.
The peak ranges, or “valley network” of $V_{u}$ (the unstable manifold of $u$
in the context of discrete Morse theory) demarcate the domains of eigenstate
support [38; 46].

To apply the LL to the many-body case, we map the Hamiltonian in Eq. 3 to a
TBM on the Fock space lattice/graph $G_{\mathcal{F}}=(V,E)$. Here we closely
follow the notation and setup given in Ref. [31], in which several properties
of the Fock space lattice and on-site energies are described.

The physical lattice has $N$ sites with $N_{e}=\lceil\nu N\rceil$ fermions at
filling fraction $\nu$, set to $\lfloor N/2\rfloor$ in this work. The
dimension of the Fock space graph is then
$|V|=\,^{N}C_{N_{e}}=\binom{N}{N_{e}}$, with nodes $V$ indexed by
$\\{\ket{I}\\}$, the eigenstates of $H_{0}=H_{W}+H_{V}$. This is in contrast
to the Anderson basis which diagonalizes $H_{t}+H_{W}$. Each $I$ can be
represented by a tuple of $N$ occupation numbers $n_{i}^{(I)}=0$ or $1$ for
each physical site $i$ such that $\sum_{i}n_{i}^{(I)}=N_{e}$. We note that
$N_{\mathcal{H}}\propto\exp(s_{\infty}N)$ where

$\displaystyle s_{\infty}=-\left(\nu\ln\nu+(1-\nu)\ln(1-\nu)\right)$ (7)

is the per-site configurational entropy [31].

The edges of $G_{\mathcal{F}}$ correspond to states (nodes) $I,J$ connected by
$H_{t}$, i.e.,

$\displaystyle H=\sum_{I}\mathcal{E}_{I}\ket{I}\bra{I}+\sum_{\langle
IJ\rangle}t\ket{I}\bra{J}.$ (8)

We depict in the bottom right of Fig. 2 $G_{\mathcal{F}}$ for $N=8=2N_{e}$.

We denote by $N(I)$ the neighborhood of $I$, such that the coordination number
$Z_{I}=|N(I)|$. The average coordination number over $G_{\mathcal{F}}$ is
given by [31]

$\displaystyle\overline{Z}=2(1-\nu)N_{e}.$ (9)

The competition of the extensive coordination number $\overline{Z}$ (favoring
delocalization), with the strong correlations of neighboring on-site energies
(discussed in Sec. III) is fundamental to the enigma of whether MBL survives
in the thermodynamic limit [31].

Turning back to the LL, the TBM in Eq. 8 $u_{2}$ is defined so long as $H\geq
0$, which is always achievable by adding a positive offset. In fact, any
stoquastic Hamiltonian yields a positive $u$ [52]. From $V_{u}$ a derived,
degenerate Agmon metric on Fock space bounds the extent of low-energy
eigenstate spread into classically disallowed regions [52].

However, one significant shortcoming of the LL is its ineffectiveness in
bounding eigenstates other those close to the ground state. This follows from
the definition of $u$ in Eq. 5: low-energy eigenstates contribute more to the
LL. An analogue dual formulation of the LL can be made for eigenstates at the
top of the spectrum [53], but both parts of the spectrum localize at zero
disorder strength in the thermodynamic limit for the Hamiltonian in Eq. 8.
Thus, the traditional LL is not well-equipped to probe mid-spectrum eigenstate
Fock space anatomy.

Directly in response to this limitation, a new localization landscape was
introduced in Ref. [46], termed the $L_{2}$ landscape. The function $u$ is
modified to take the form

$\displaystyle u_{I}^{(2)}:=\sqrt{M^{-1}_{II}},\,M=H^{\dagger}H.$ (10)

Similar to Eq. 6, this new landscape satisfies

$\displaystyle|\phi_{I}^{\beta}|\leq|E^{\beta}|||\phi^{\beta}||_{2}\sqrt{M^{-1}_{II}},$
(11)

where by orthonormality the $\ell_{2}$ norm equals unity (hence the name
“$L_{2}$ landscape”). Since $M>0$ for any $H$, we shift $H\to H-(E_{0}+i\eta)$
and probe the eigenspectrum near $E_{0}$, where $\eta\in\mathbb{R}^{+}$ is an
infinitesimal offset. From this shift the $L_{2}$ landscape can be written as
a function of the resolvent $G^{\pm}(E_{0}):=(E_{0}\pm i\eta-H)^{-1}$:

$\displaystyle(u^{(2)}_{J}(E_{0}))^{2}=\left(G^{+}(E_{0})G^{-}(E_{0})\right)_{JJ},$
(12)

which for finite systems $u^{(2)}$ is exactly given by

$\displaystyle(u_{J}^{(2)}(E_{0}))^{2}=\sum_{\beta}\frac{|\phi_{J}^{\beta}|^{2}}{|E_{0}+i\eta-
E_{\beta}|^{2}}.$ (13)

Closely related is the local density of states

$\displaystyle\rho_{J}(E_{0}):=\lim_{\eta\to 0^{+}}\frac{1}{2\pi
i}\left(G^{+}_{JJ}(E_{0})-G^{-}_{JJ}(E_{0})\right).$ (14)

It therefore follows [46]

$\displaystyle\rho_{J}(E_{0})=\lim_{\eta\to 0^{+}}\eta(u_{J}^{(2)})^{2}.$ (15)

The behavior of the resolvent $G^{\pm}$ is the central consideration of
several recent works [15], which we discuss further in Sec. III.

The $L_{2}$ landscape incorporates information from the full eigenspectrum,
and has a nonlinear dependence on the offset $\eta$, the consequences of which
we address in Sec. III. Due to the connectivity of $G_{\mathcal{F}}$, it is
not immediately obvious what observables can be extracted from the $L_{2}$
landscape. However, the authors in Ref. [46] where the $L_{2}$ landscape was
first defined posed an open question of how fractal dimensions like those
extracted from the participation entropies can be generalized to the $L_{2}$
landscape. An analogue of the fractal dimension for $\rho(E_{0})$ was recently
proposed [54], but we explore in Sec. III a novel construction by leveraging a
tool originating in Morse theory and rooted in algebraic topology: persistent
homology.

###  II.2 Persistent Homology: motivation and theory

Persistent homology characterizes the small and large-scale topological
structure of an underlying dataset by computing how the homology of a
filtration of topological spaces changes as a function of scale. This
fundamental idea underlies the many generalizations of persistent homology;
for the interested reader there exist many reviews such as Ref. [55].

Rather than focus on the technical aspects of persistent homology, here we
give an illustrative example quite similar in spirit to the construction we
later present in Sec. III. Consider a function $f:[0,1]\to\mathbb{R}$ depicted
in Fig. 1. Our objective is to characterize the structure of $f$; i.e., how
the extrema relate to one another. To quantify this structure, we consider the
superlevel set $L_{c}^{+}(f)=\\{x\in[0,1]|f(x)\geq c\\}$. For any value of $c$
(the dashed lines in Fig. 1, known as the filtration value) the superlevel set
is a set of intervals in $[0,1]$, indicated by the black bars at the bottom of
each panel. The number of intervals is the number of connected components
(clusters); in topological terms, the Betti number $\beta_{0}$. As we decrease
$c$, the number of connected clusters fluctuates until, when
$c=c_{min}:=\min_{x}f(x)$, we have $L_{c_{max}}^{+}=[0,1]$, the domain. Along
the filtration we can monitor which clusters merge into one another.

The end result of capturing the progression of homology groups is a
persistence diagram showing the filtration value at which a connected
component appears in the filtration (known as the birth time) on the ordinate,
and when a cluster merges with another cluster (the death time on the
abscissa. Each persistence point has a distance from the diagonal: the farther
from the diagonal, the longer-lived the cluster was in the filtration. Fig. 1
depicts the persistence diagram with notation indicating the minima and maxima
of the functional that identify each cluster. Note that the point
$(a,-\infty)$ indicates the largest component, the domain, which technically
has infinite lifetime. Next to the persistence diagram in Fig. 1 we depict the
dendrogram, or merge tree corresponding to the superlevel set filtration:
external (internal) nodes denote the maxima (minima).

One remarkable property of persistent homology is its stability with respect
to perturbations of the underlying functional. In particular, a series of
algebraic stability results imply the Wasserstein metric of two persistence
diagrams is continuous with respect to $\ell_{p}$ norms between the underlying
filtration functionals [56; 57]. Put in simpler terms, small perturbations to
the superlevel sets result in small perturbations to the persistence diagrams.

Figure 2: Example of the Fock space $G_{\mathcal{F}}$; colors of nodes denote
the value of a random potential. For each value of
$\epsilon\in\\{0.7,0.5,0.2,0.0\\}$ the superlevel set (subgraph)
$G[V_{\epsilon}]$ is displayed.

The unique idea of persistent homology is that the aforementioned pipeline is
readily generalizable to higher dimensions. The filtration of intervals
generalizes to filtrations of simplicial complexes, and the notion of clusters
generalizes to homology groups of a sequence of simplicial complexes. A
filtration need not only be superlevel set filtrations of an underlying
functional: any filtered sequence of topological spaces (given a total
ordering) has a unique representation via a persistence diagram [58].

The toy example described above informs our study of the morphology of the
$L_{2}$ landscape. Instead of a functional on a continuous domain, we perform
a superlevel set filtration on $G_{\mathcal{F}}$ given the function
$u^{(2)}(E_{0})$ on $V(G)$. For each filtration value $c$, the analogue to the
superlevel set is the induced subgraph $G[V_{c}]$, where $V_{c}:=\\{I\in
V(G_{\mathcal{F}})|u_{I}^{(2)}\geq c\\}$. Note that the induced subgraph
$G[V_{c}]$ includes both the vertex set $V_{c}$ as well as all of the edges
$E$ that have both endpoints in $V_{c}$.

While perhaps a less intuitive than a sublevel set filtration, we chose to
perform a superlevel filtration for one major reason. The peaks of the $L_{2}$
landscape indicate where eigenstates localize, while the minima demarcate and
separate the regions eigenstate near $E_{0}$ localize. If we performed a
sublevel filtration, in the extended regime a large fraction of Fock space
sites $I$ would be connected into one cluster well before the maxima
indicating the localization centers would appear. Upon appearing these maxima
would immediately merge into a connected component and therefore not be
recorded in the persistence diagram. By performing a superlevel filtration we
ensure the localization center feature in the diagrams.

In Fig. 2 we depict for visual reference a superlevel set filtration on
$G_{\mathcal{F}}$ for $N=8=2N_{e}$ in 1D and with open boundary conditions.
Due to the absence of odd-length cycles (the graph is bipartite), our
filtration fails to produce any interesting homology in dimensions greater
than zero, as any loops generated persistent throughout the rest of the
filtration, and no higher order simplices are generated. However, a different
interaction term $H_{t}$ will generate higher dimensional homology, as we
discuss in Sec. IV. As noted above, our pipeline for 0D persistent homology
corresponds to single-linkage hierarchical clustering (depicted via the merge
tree in Fig. 1).

There are a multitude of other persistent homology construction that we could
have attempted in this work. Of particular interest would be filtrations based
on the distance function derived from the Agmon metric [52; 41], which has
deep theoretical connections to localization bounds. While an Agmon metric for
the $L_{2}$ does not currently exist, it is an exciting prospect and the
subject of future work.

##  III Results

With the theoretical underpinnings of the many-body localization landscape and
the persistent homology pipeline established, we turn our methodology to the
prototypical MBL: the disordered Heisenberg Hamiltonian given in Eq. 3. Our
objective again is to assess whether the $L_{2}$ landscape shows indications
of a morphological transition upon approaching the MBL phase. To build these
indicators we leverage novel topological summaries and observables generated
from persistence diagrams and merge trees. Broadly speaking, we find these
observables best discriminate the transition when compared across system
sizes.

For the Hamiltonian in Eq. 3 we assume open boundary conditions and performed
$1000$ disorder realizations for the $N=[10,14]$ system sizes and over $500$
for the $N=15$ system size. We sampled $W\in[0.1,8.0]$ in intervals of $\delta
W=.2$ as well as sampled $W\in\\{10,12,14,16,18,20\\}$ deep in the MBL phase
for comparison. We use exact diagonalization to compute the eigenspectrum for
system sizes $N\in[10,15]$.

The critical disorder strength for the phase transition in Eq. 3 is typically
given in the range $W_{c}\in[3.5,4]$, though some works place $W_{c}$ at much
higher values [59; 28; 60]. Whether or not the phase transition survives the
thermodynamic limit remains an open question, with recent work posited a
strong distinction between the MBL phase and the finite-size MBL regime that
can be probed with numerics [28]. Given our limited system sizes we take no
position on the stability issue or the absence of mobility edges in the
thermodynamic limit in this work, choosing instead to see if morphological
indicators coincide with critical values claimed in the literature.

Before moving to the persistent homology topological summaries, we briefly
digress to discuss the statistics of eigenenergies in Fock space. Given the
$L_{2}$ landscape’s close correspondence to the local density of states and
$G^{\pm}$, the topological summaries explored below probe similar statistics,
greatly expounded upon in recent works we now refer to [31].

Note that both the on-site energies and eigenenergies $\mathcal{E},E$ share
the same $\overline{E}$, as $H_{t}+H_{V}$ is disorder independent. Due to
parity and our choice of $V=1$ we have $\overline{E}=-1/4$ for even system
sizes. Meanwhile the variance $\mu_{E}^{2}$ goes as $\mathcal{O}(N)$.

We depict in Fig. 3 the eigenenergy probability distribution (normalized
histogram) for $W\in[1.0,5.0,12.0]$ at $N=15$, as well as the Gaussian
approximation to the density of states

$\displaystyle
D(\omega)=(2\pi\mu_{E}^{2})^{-1}\exp\left(-(\omega-\overline{E})^{2}/2\mu_{E}^{2}\right),$
(16)

shown as the dashed line for comparison.

Figure 3: Normalized histogram for $L=15$ of the eigenenergies for
$W=(1,5,12)$, the dashed lines show the Gaussian approximation given in Eq. 16

To evaluate $u^{(2)}$, we must select $E_{0}$ and a finite value for $\eta$
which should be smaller than the level spacing near $E_{0}$, but not so small
as to impact numerical stability [46]. The choice of $\eta$ dictates the
nonlinear extent to which eigenstates far from $E_{0}$ contributes to the
$L_{2}$ landscape. To first order Eq. 13 shifts under $\eta\to\eta+d\eta$ as

$\displaystyle(u_{J}^{(2)})^{2}\to(u_{J}^{(2)})^{2}-2\eta
d\eta\sum_{\beta}\frac{|\phi_{J}^{\beta}|^{4}}{|E_{0}+i\eta-E_{\beta}|^{4}}.$
(17)

One natural choice for $\eta$ is the eigenvalue spacing near $E_{0}$, given by
$(N_{\mathcal{H}}D(E_{0}))^{-1}$. Here $D(E_{0})$ is the density of states
given in Eq. 16, such that $D(E_{0}=\overline{E})=(\sqrt{2\pi}\mu_{E})$. This
then implies $\eta=\sqrt{2\pi}\mu_{E}/N_{\mathcal{H}}$.

However, as noted above, $\mu_{E}$ depends both on $N$ and $W$. Moreover, the
level statistics of the model in Eq. 3 has long been known to shift under the
phase transition from GOE (Wegner-Dyson) to Poissonian statistics [61]. The
precise nature of the level statistics near the transition is still very much
controversial [62], and so it is unclear if a more canonical choice for $\eta$
can be made. We make the choice $\eta=\sqrt{2\pi}\mu_{E}/N_{\mathcal{H}}$ for
two main reasons.

First, the spectral width $\mu_{E}$ is non-parametric in the sense that we do
not need to average over some number $k$ level spacing near $\overline{E}$.
The number of level spacing necessary for estimation could very well depend
upon the disorder strength, rendering this technique possibly biased.

Second, the spectral width $\mu_{E}$ is fundamental to the analysis given in
Refs. [31], on which we have based our notation and exposition for Fock space
spectral statistics. There the authors analyze the stability of the MBL phase
in the thermodynamic limit by considering functionals of the resolvent
$G_{IJ}^{\pm}$ (particularly, the Feenburg self-energy) under a rescaling
$\tilde{E}_{\beta}\to(E_{\beta}-\overline{E})/\mu_{E}$; i.e., the
eigenspectrum is normalized to have vanishing mean and a standard deviation of
one. Under this rescaling the offset becomes
$\tilde{\eta}=\eta/\mu_{E}=\sqrt{2\pi}/N_{\mathcal{H}}$. The rescaling proves
central to a mean-field theory approach to the self-consistent determination
of the typical self-energy [31]. Note that $u^{(2)}(\overline{E})$ would be
rescaled to $\tilde{u}^{(2)}(0)$ with offset $\tilde{\eta}$, and so we have
$\tilde{u}^{(2)}(0)=\mu_{E}u^{(2)}(\overline{E})$. Choosing $\eta$ is a
nuanced issue and should be the study of future work.

Figure 4: Representative persistence diagrams for disorders strengths
$W=\\{1.0,3.0,5.0,8.0\\}$, plot on a log-log scale. Note that $d<b$ due to the
superlevel set filtration. Inset shows the same persistence diagram on a
linear scale.

Persistence Diagram –  Turning to the output of the persistent homology
pipeline, Fig. 4 depicts representative persistence diagrams for the disorder
strengths $W=\\{2,8,14\\}$. Recall that the abscissa and ordinate are the
birth and death times, respectively; due to the superlevel set filtration, the
death times are smaller than the birth times. We have plotted the persistence
diagram in both log-log and linear scale (inset) for comparison. The
persistence pairs for the extended regime $W=2$ cluster at high filtration
values with a broad spread of death times. In contrast, in the localized
regime the birth and death times are much smaller (corresponding to a small
value of $u^{(2)}$).

A perhaps naive interpretation of this behavior is that
$|\phi_{J}^{\beta}|^{2}\sim(N_{\mathcal{H}})^{-1}$ on
$\mathcal{O}(N_{\mathcal{H}})$ many sites in the extended regime. In contrast,
in the localized regime $|\phi_{J}^{\beta}|^{2}\sim N_{\mathcal{H}}^{-\alpha}$
for $\alpha<1$ on $\mathcal{O}(N_{\mathcal{H}}^{\alpha})$ sites. We note that
the number of persistence pairs for each diagram corresponds to the total
number of maxima and is therefore not constant, a fact that proves
consequential when we consider a generalization of the fractal dimension.

Figure 5: Statistics on the birth, death, and lifetime distributions on the
persistence diagrams. (a-b) depict the mean (solid), median (dashed) and
standard deviation of the birth values as a function of disorder strength $W$.
(c-d) and (e-f) show the corresponding statistics for the death and lifetimes
values, respectively.

$(b,d,l)$ statistics –  Many topological summaries center on properties of the
distributions of the set birth, death, and lifetime values, which we write as
vectors $b,\,d,\,l$, respectively. Note that we exclude the persistence pair
corresponding to the infinitely long-lived component ($(a,-\infty)$ in Fig.
1). As we explore in Sec. III.0.1, norms of these vectors correspond to
observables relevant to percolation theory.

Fig. 5(a) and (b) depicts the arithmetic mean $\overline{b}$, geometric mean
$l_{typ}$ (dashed lines), and standard deviation $\mu_{b}$ of the birth
distribution against $W$. Fig. 5(c-f) show equivalent data for the death and
lifetime distributions. All statistics are calculated within for each disorder
realization before taking averaging over realizations.

Plots (a) and (c) show the mean birth and death times decay monotonically with
$W$, visually reflecting the change in spread of points in Fig. 4. The
geometric mean (or typical value) of the birth/death times deviate strongly
from the arithmetic mean towards the localized regime, with the deviation more
pronounced with larger system sizes and particularly for the birth times. This
implies large outliers in the birth times, coincident with the onset of
localization centers in the localized regime. In contrast, the difference
between arithmetic and geometric mean death times is less pronounced, as the
death times (minima of $u^{(2)}$) are lower bounded by zero.

While the arithmetic mean of $l$ is the difference in the mean birth and death
times, the same cannot be said for the typical lifetime $l_{typ}$, shown in
(e). Here the distinction between arithmetic and geometric mean is even more
pronounced, and a peak is observed at $W\approx 1.5$ for $l_{typ}$ and
$W\approx 2$ for $\overline{l}$ for the $N=15$ system size.

Plot (b), (d), and (f) depict $\mu_{X}$ for $X=(b,d,l)$, respectively, wherein
we see peaks generally obey $\mu_{l}^{max}\geq\mu_{b}^{max}>\mu_{d}^{max}$.
This follows from the fact that $\mu_{l}$ should add in quadrature, i.e.,
$\mu_{l}^{2}\approx\mu_{b}^{2}+\mu_{d}^{2}$. The asymmetry in the peak
$\mu_{b},\,\mu_{d}$ values again stems from the fact that both $u^{(2)}\geq 0$
and $d_{i}<b_{i}$ for any persistence pair $(b_{i},d_{i})$.

Figure 6: Fractal scaling $\nu_{p}$ for $p=\\{0,0.5,1.0,1.5\\}$, corresponding
to plots (a-d), respectively. Figure 7: (a) Slope $\tilde{\beta}_{p}$ of
$\log||l||_{p}^{p}$ against $\log N_{\mathcal{H}}$ for
$p=(0.0,0.5,1.0,1.5,2.0)$. (b) Fractal dimension $\dim_{p}(PD)$ for the case
$p\to 0$. (c) The fractal dimension plotted on a semilog scale.

####  III.0.1 $\ell_{p}$-norms of the persistence diagrams

We now assess to what extent $\ell_{p}$ norms of the $(b,d,l)$ vectors are
sensitive to the phase transition. As we outline below, several of these norms
correspond to meaningful quantities in the context of percolation theory, and
thus should help shed light on the morphological structure of Fock space under
the $L_{2}$ landscape.

Fractal dimensions –  The first observable we consider is how the $\ell_{p}$
norms of $l$ scale with $(N,W)$. We define the function

$\displaystyle\nu_{p}(N,W)=\frac{\log\langle||l||_{p}^{p}\rangle}{\log
N_{\mathcal{H}}},$ (18)

where $\langle\cdot\rangle$ denotes averaging over disorder realizations
holding $W$ fixed. For $p<1$ $\nu_{p}$ is a semi-norm that emphasizes small-
scale homological features. While persistent homology is often used to
identify large-scale topological structure, the small-scale structure is
relevant in the context of fractality and estimating curvature [63; 64].
Indeed, $\nu_{p}$ was first defined to give a notion of “persistent fractal
dimension” derived from persistent homology and comparable to more traditional
measures of fractality like the correlation and box-counting dimensions [63].
In its original formulation the persistent fractal dimension is applicable to
filtrations built with respect to a finite metric space, such a Vietoris-Rips
or Cěch filtration [63]. In our case we apply the same methodology to the
$u^{(2)}$ landscape. We interpret $\nu_{p}$ as quantifying the relative weight
of small-scale homological features (short-lived clusters) as a function of
$(N,W)$. Note that $||l||_{0}$ counts the number of $u^{(2)}$ maxima.

We define a persistent fractal dimension as [63]

$\displaystyle\text{dim}_{p}(PD)=\frac{p}{1-\tilde{\beta}_{p}},$ (19)
$\displaystyle\tilde{\beta}_{p}=\lim_{n\to\infty}\sup\nu_{p}.$ (20)

Here $\tilde{\beta}_{p}$ amounts to the slope of $\log||l||_{p}^{p}$ against
$\log N_{\mathcal{H}}$. Note that the tilde on $\beta$ is to distinguish
$\tilde{\beta}_{0}$ from the Betti number $\beta_{0}$, and does not refer to
any rescaling such as $u^{(2)}\to\tilde{u}^{(2)}$.

Fig. 6(a-d) depicts $\nu_{p}$ for $p\in\\{0,0.5,1.0,1.5\\}$. For $p<1$ we
observe a noticeable “kink” in $\nu_{p}$ as $W$ increases, close to
$W_{c}\approx 3.5$. For $p>0$ $\nu_{p}$ decreases monotonically with $W$,
which upon considering Fig. 5(e) seems contradictory. However, given that
$\overline{l}=||l||_{1}/||l||_{0}$, the peak in Fig. 5(e) stems from the
relative decays rates of $\nu_{1}$ versus $\nu_{0}$.

For each $p\in\\{0,.5,1.0,1.5,2.0\\}$ we calculated $\tilde{\beta}_{p}$,
depicted in Fig. 7(a). While $\tilde{\beta}_{p}$ visually appears to scale as
$\log W$ in the extended regime, there is a small sub-logarithmic factor. Note
that for $p\geq 1/2$ the fractal dimension $\text{dim}_{p}(PD)$ is either
divergent at some value of $W$ or negative; however, $\tilde{\beta}_{p}<0$ for
$p=0$. Fig. 7(b) and (c) depict $\text{dim}_{p}(PD)/p$ for $p=0^{+}$. The
inset (c) shows the same plot as (b) but on a log scale. For $W\geq W_{c}$
$\text{dim}_{0+}(PD)/0^{+}\approx\log W^{a}+b$ with $(a,b)=(.79,3.0)$.

There are several interesting remarks to make concerning the case of $p=0$.
First, note that $||l||_{0}$, which we also write as $|PD|$, equates to the
total number of persistence pairs in a persistence diagram (hence $|PD|$
denotes cardinality). This quantity is invariant under any monotonically
increasing transformation of the underlying filtration; in particular, $|PD|$
is the same whether we use $u^{(2)}$ or $\tilde{u}^{(2)}$. Second, $|PD|$
equates to the number of maxima in the $L_{2}$ landscape. These peaks
qualitatively correspond with Fock space sites where eigenstates near
$\overline{E}$ also peak. Third, we can compare the scaling behavior of
$\nu_{0}$ against the scaling of the Hilbert space dimension $N_{\mathcal{H}}$
against $N$. As noted above, $N_{\mathcal{H}}\propto\exp^{s_{\infty}N}$, where
$s_{\infty}$ is the configurational entropy; therefore, $\log
N_{\mathcal{H}}\propto s_{\infty}N$. Thus, we could equivalently have computed
$\tilde{\beta}_{0}$ by fitting $\log||l||_{0}$ against $N$, the system size.
Then $\tilde{\beta}_{0}\to\tilde{\beta}_{0}/s_{\infty}>1$ which implies the
fractal dimension would become negative.

To try and understand this behavior a bit more analytically, we return to
$\nu_{0}$ and consider the behavior of the local density of states
$D_{I}(\overline{E}_{0})$. Deep in the ergodic phase the eigenstates can be
well-approximated by Gaussian random vectors with vanishing mean and a
standard deviation that scales as $N_{\mathcal{H}}^{-1/2}$ [20]. By time-
reversal symmetry the eigenstates are real, and thus the contribution of each
eigenstate to the $L_{2}$ landscape is effectively random apart from the
denominator, which can also be well-approximated by a Gaussian, see Eq. 16.
The lack of correlations for $u^{(2)}$ serves as a proxy for a proliferation
of local maxima, as neighboring Fock space sites are effectively uncorrelated.
As the disorder becomes non-vanishing and the model moves from integrability
the eigenstates drift from uncorrelated Gaussians, thus lowering the number of
maxima. A more simple picture is that, as the eigenstates begin to localize in
regions (clusters of Fock space), the appearance of large amplitude maxima
necessitates a decrease in the total number of maxima.

In the other extreme, deep in the MBL phase, we expect eigenstates near
$\overline{E}_{0}$ to be tightly clustered around localization centers. Away
from these centers, the $L_{2}$ landscape is very small (see Fig. 4) and again
effectively random, which drives up the number of local maxima.

This qualitative description holds as well for the non-interacting $(V=0)$
case, as deep in the MBL phase the model is perturbatively connected to the
non-interacting limit [20]. We can leverage the perturbative connection to the
non-interacting limit to partly explain the $\log W$ behavior in $\nu_{p}$. As
the non-interacting eigenstates are Slater determinants of the Anderson
orbitals, the localization length for an eigenstate scales as $(\log W)^{-1}$
[3]. Thus the proliferation of uncorrelated regions of the $L_{2}$ landscape
should rise in proportion to the decrease in the localization length, which by
proxy implies the number of local maxima scale with $\log W$ as observed.
Because the number of maxima decrease in the ergodic regime and increase in
the MBL regime, there is a natural point between where the number of maxima is
minimal, near $W=3.0$.

Connectivity threshold –  In the context of percolation theory, one figure of
merit is the probability of obtaining a giant connected component on Erdős-
Rényi random graphs, wherein a probability is assigned to each edge [65]. For
these systems of random graphs a large body of literature focuses on analogues
of the MBL localization/delocalization transition with respect to the
adjacency matrix [66; 67]. In the context of our study of the $L_{2}$
landscape morphology, an analogous figure of merit is played by the
connectivity threshold $d_{f}$, which we define as the filtration value at
which only one connected component remains. By construction, this corresponds
to the smallest value of $u^{(2)}$ (global minimum), or equivalently, as
$||d||_{-\infty}$. Fig. 8 depicts the connectivity threshold $d_{f}$ against
the $W\in[2.0,3.0]$ wherein we observe the different system sizes cross at
various points. The inset shows $d_{f}$ on a log scale across the entire set
of $W$ probed. In the localized regime the power law decay goes as
$W^{-\alpha}$ with $\alpha\approx 1.15$, which we depict with the dashed black
line. The crossing points in the main plot of Fig. 8 are again substantially
lower than $W_{c}$; it is unclear if the reduced finite size effects at larger
values of $N$ would resolve these crossing points.

Figure 8: Log of the connectivity threshold $d_{f}$, the last superlevel set
wherein at most two connected clusters remain, against $W$. The main plot
depict the connectivity threshold near where different system size curves
coincide, the inset depicts the curves from $W\in[0,20]$. the dashed black
line depicts an the approximate power law scaling $d_{f}\propto W^{-\alpha}$
for $\alpha=1.15$.

Maximum spanning tree –  Also related to percolation theory is the notion of
maximum (MaxST) and minimum spanning trees (MinST). On a graph $G=(V,E)$ with
edge weight function $w:E\to\mathbb{R}$, the MaxST (MinST) (if it exists, else
a spanning forest) is an acyclic connected subgraph such that the sum of the
edge weights is maximal (minimal). The weight of the MaxST corresponds to the
sum of the death times $||d||_{1}$ of the persistent diagram and is therefore
equivalent to the integral of the reduced Betti number with respect to the
filtration [68]. The reduced Betti number at filtration value $c$ is
specifically the rank of the reduced homology group $H_{0}(G[V_{c}],I_{max})$,
where $I_{max}:=\text{argmax}_{V(G)}u_{I}^{(2)}$. Put more plainly, the
reduced Betti number is generally one less than the regular Betti number.

Fig. 9 (a) depicts $w(\text{MaxST})/N_{\mathcal{H}}$ against $W$ on a log-log
scale, whereby the power law scaling in the localized regime is manifest. The
number of terms in $w(\text{MaxST})$ scales with $|PD|\propto\log W$, which
implies the average death time scales as $W^{\gamma}/\log W$ for some scaling
exponent $\gamma$ that varies slightly for different system sizes.

Figure 9: (a) the weight of the maximal spanning tree $w(\text{MaxST})$
against $W$, normalized by $N_{\mathcal{H}}$. (b) The persistent entropy
$S(PD)$ against $W$. Figure 10: The exponent $\alpha_{S(PD)}$ from fitting
$S(PD)\propto N^{\alpha}$. Inset shows the range $W\in[2,7]$.

Persistent entropy – We turn now to an observable that closely mirrors the
participation entropy in accounting both for cardinality and magnitude of
persistence pairs: the persistent entropy. By normalizing $l$ to a probability
vector $\tilde{l}:=l/||l||_{1}$, we define the persistent entropy as

$\displaystyle S(PD):=-S(\tilde{l}),$ (21)

where $S(\cdot)$ as before is the entropy function [69]. We can equate $S(PD)$
to a functional of the $\ell_{p}$ norms:

$\displaystyle S(PD)=\frac{S(l)}{||l||_{1}}-\sum_{i}\ln l_{i}.$ (22)

The persistent entropy has been extensively used in the context of detecting
both classical and quantum phase transitions [70]. Qualitatively, for a fixed
value $|PD|$ a large $S(PD)$ implies the lifetimes are largely constant. In
contrast, a small $S(PD)$ implies the prevalence of vanishing few clusters
with lifetime disproportionately large. Clearly $S(PD)\leq\log|PD|$, just as
the Schmidt rank bounds the entanglement entropy and $\log N_{\mathcal{H}}$
bounds the participation entropies.

Fig. 9(b) depicts the persistent entropy $S(PD)$ as a function of $W$. In
contrast to the observables examined above, $S(PD)$ flattens out to nearly
constant, with a very slow, sublogarithmic scaling in $W$. Qualitatively, near
constant $S(PD)$ is analogous to the area law entanglement entropy exhibited
by localized states. Heuristically, the lifetimes of the Fock site clusters
indicating eigenstates near $\overline{E}$ roughly correspond to the large
Schmidt coefficients of a localized eigenstate under a bipartition. We expect
a relatively large contribution to the $L^{2}$ landscape from these
localization centers while the rest of the landscape is vanishingly small.
There is a competition between the number of local maxima (increasing as $\log
W$) versus the relative persistence of these maxima. Fig. 9 indicates that the
maxima proliferation outweighs the disproportional lifetime of these maxima in
slowly driving up the persistent entropy.

Just as $S(\rho_{A})$ scales with $N$ in the extended regime, so too does
$S(PD)$; in Fig. 10 we depict the coefficient $\alpha$ from a fit
$S(PD)\propto N^{\alpha}$, with inset zoomed into the region $W\in[2.0,7.0]$.
The coefficient $\alpha$ plateaus in the region $W\in[4,6]$ before slowly
increasing with $W$. This behavior is very similar to Fig. 7 (b) which
depicted the fractal dimension, though the minimum value for $\alpha$ seems be
at a higher $W$ and the scaling is not quite logarithmic.

Figure 11: The peak value of $\beta_{0}$ during the filtration (equivalently,
the maximum number of independent clusters at any point during the
filtration), normalized by $N_{\mathcal{H}}$. Note that $\beta_{0}^{max}$ goes
as $\log W$ in the localized regime. Figure 12:  Scaling exponents from
fitting $\beta_{0}^{max}\propto N_{\mathcal{H}}^{a}\left(\log
N_{\mathcal{H}}\right)^{b}$. The horizontal dashed lines at $0$ and $.5$
roughly correspond to the maximum (minimum) value of $b$ ($a$), which occurs
at $W\approx 3.0$. Inset depicts representative fitting of $\beta_{0}^{max}$
against $N_{\mathcal{H}}$ at $W=3.0$.

Maximum Betti number – The final observable we consider is the maximum value
of the Betti number $\beta_{0}^{max}$ encountered during the superlevel set
filtration. This observable has been leveraged in recent works as an indicator
for quantum phase transitions [51]. Fig. 11 depicts
$\beta_{0}^{max}/N_{\mathcal{H}}$ against $W$ wherein we see that, quite
similar to $|PD|$, a logarithmic growth in the extended regime. Appealing to
the heuristical argument given above, in the extended regime $\beta_{0}^{max}$
occurs very late in the filtration, growing in proportional to the extent of
decorrelated Fock space. Fig. 12 depicts the scaling of
$\beta_{0}^{max}\propto N_{\mathcal{H}}^{a}\left(\log
N_{\mathcal{H}}\right)^{b}$ with exponents $(a,b)$ disorder-dependent.
Intriguingly, at $W\approx 3.0$ $(a,b)\approx(.5,0)$, which implies
$\beta_{0}^{max}\propto\sqrt{N_{\mathcal{H}}}$. This is roughly the same
disorder value for which the minimum in $\text{dim}_{0+}(PD)$ occurred. Though
it is at present unclear if the extremal behavior of the scaling exponents
$(a,b)$ would shift for larger values of $N$, the simplicity of the scaling at
$W\approx 3.0$ indicates deeper behavior that merits exploration in future
work.

##  IV Conclusion

This work takes an exploratory approach to leveraging persistent homology as a
descriptor for the complex correlational structure of eigenstates in Fock
space. We utilize a new construction, the $L_{2}$ landscape, to procure
topological summaries of the mid-spectrum Fock space structure.

Overall, we found that several observables, including the cardinality $|PD|$
of the persistence diagrams (Fig. 7), the persistent entropy (Fig. 10) and the
peak Betti number (Fig. 11) all exhibit scaling in either $N$ or
$N_{\mathcal{H}}$ near $W_{c}$ indicative of a transition. This is especially
evident for $|PD|$ which is perhaps the most easily interpreted observable
presented here, as it corresponds to the number of maxima of an approximate
local density of states.

Along with these observables we proposed a novel notion of a fractal dimension
readily applicable to other functionals defined on Fock space, including but
not limited to self-energies and the local density of states. While a recent
study expanded the participation entropy definition to a normalized version of
the local density of states [54], our approach is fundamentally different. The
notion of fractal dimension presented here incorporates not only the Fock
space lattice but also the complex morphology of level sets.

Still other observables painted a less intuitive picture of the MBL phase
transition, and it should be the aim of future work to discriminate
observables that inaccurately pinpoint the phase transition due to finite-size
effects from observables that do not probe the transition at all. Two possible
indicators for the transition, the connectivity threshold $d_{f}$ (Fig. 8) and
the $(b,d,l)$ statistics (Fig. 5), show substantial drift in the crossing of
different system size curves, indicative of the finite-size effects at the
relatively small system sizes probed here.

Throughout this work we have posed a heuristical understanding that is largely
based upon the local density of states which, by Eq. 15, the $L_{2}$ is
connected to. However, lacking a deeper, more analytic understanding of how
$u^{(2)}$ evolves with $W$ limits the scope to which a topological data
analysis lens can shed light on the transition. Much of this limitation stems
from the difficulty in interpreting persistent homology in terms of standard
statistical observables. While persistent homology is strongly connected to
discrete Morse theory [71], to our knowledge there exists a limited amount of
work in the topological data analysis literature to bridge homological and
statistical observables. We conclude with several proposed future directions
to take this exploratory analysis, as the use of persistent homology in the
context of quantum information and many-body systems is still nascent.

Persistent Homology for functionals –  The PH pipeline explored here can
equally be applied to individual eigenstate as performed in Ref. [49] for the
single-body Aubry-André model. In fact, the persistent homology pipeline could
be readily applied to any functional represented on $G_{\mathcal{F}}$,
including the local density of states and the self-energies. The edges
$E(G_{\mathcal{F}})$ can have filtration values that differ from their
boundary nodes.

As noted in Sec. II, the original $LL$ gives rise to a degenerate Agmon metric
on $G_{\mathcal{F}}$. From this metric a series of localization theorems limit
the extent eigenstate support extended into “classically disallowed regions”
[52]. This Agmon metric determines a filtration on the complete graph $G$ with
vertex set $V(G)$ the Fock space sites. However, a complete graph would become
computationally prohibitive with large $N$, at which point approximate methods
such as witness complexes could be used [72].

Higher-dimensional Homology –  The analysis given in this work is restricted
to 0D homology due to the interaction graph $G_{\mathcal{F}}$: the
bipartiteness of the graph precludes any $2$-simplices from forming. This
obstruction can, of course, be overcome by simply choosing a different
Hamiltonian; e.g., including a $\sigma_{x}$ term that breaks total
magnetization conservation. Under this inclusion persistent homology could be
used to study the occurrence of higher-dimensional homology features like
cycles and voids.

Approximations for larger system sizes –  For larger system sizes exact
diagonalization is computationally expensive. Thus, approximate forms of the
$L_{2}$ landscape could be computed via recursive approximations to the
diagonal of the Green function [46]. The original LL equation in Eq. 4 can be
written as a Rayleigh quotient and is therefore amenable to DMRG and tensor
network approaches. Recent tensor network representations of Green functions
could be leveraged to scale the $u^{(2)}$ calculation [73].

###### Acknowledgements.

GH acknowledges useful conversations with Felix Leditzky, Di Luo, Yuliy
Baryshnikov, Shmuel Weinberger, Giuseppe De Tomasi, and Henry Adams. This work
made use of the Illinois Campus Cluster, a computing resource that is operated
by the Illinois Campus Cluster Program (ICCP) in conjunction with the National
Center for Supercomputing Applications (NCSA) and which is supported by funds
from the University of Illinois at Urbana-Champaign. We acknowledge support
from the Department of Energy grant DOE DESC0020165.

## References

  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Reviews of Modern Physics 91, 021001 (2019).
  * [2] S. J. Thomson and M. Schiro, SciPost Physics , 39.
  * Scardicchio and Thiery [2017] A. Scardicchio and T. Thiery, arXiv:1710.01234 [cond-mat] (2017), arXiv: 1710.01234.
  * Alet and Laflorencie [2018a] F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018a).
  * Gopalakrishnan  _et al._ [2015] S. Gopalakrishnan, M. Müller, V. Khemani, M. Knap, E. Demler, and D. A. Huse, Physical Review B 92, 104202 (2015).
  * Kim  _et al._ [2014] I. H. Kim, A. Chandran, and D. A. Abanin, arXiv preprint arXiv:1412.3073 (2014).
  * Pekker and Clark [2017] D. Pekker and B. K. Clark, Physical Review B 95, 035116 (2017).
  * Alet and Laflorencie [2018b] F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498–525 (2018b), arXiv:1711.03145 [cond-mat].
  * De Roeck and Imbrie [2017] W. De Roeck and J. Z. Imbrie, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20160422 (2017).
  * Imbrie  _et al._ [2017] J. Z. Imbrie, V. Ros, and A. Scardicchio, Annalen der Physik 529, 1600278 (2017), arXiv: 1609.08076.
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Physical Review B 91, 081103 (2015), arXiv: 1411.0660.
  * Cookmeyer  _et al._ [2020] T. Cookmeyer, J. Motruk, and J. E. Moore, Physical Review B 101, 174203 (2020).
  * Schulz  _et al._ [2019] M. Schulz, C. A. Hooley, R. Moessner, and F. Pollmann, Physical Review Letters 122, 040606 (2019), arXiv: 1808.01250.
  * Morningstar  _et al._ [2020] A. Morningstar, D. A. Huse, and J. Z. Imbrie, Physical Review B 102, 125134 (2020).
  * [15] F. Monteiro, T. Micklitz, M. Tezuka, and A. Altland, , 21.
  * Vidmar  _et al._ [2021] L. Vidmar, B. Krajewski, J. Bonca, and M. Mierzejewski, Physical Review Letters 127, 230603 (2021).
  * Roy and Logan [2019] S. Roy and D. E. Logan, SciPost Physics 7, 042 (2019).
  * Imbrie [2016] J. Z. Imbrie, Physical review letters 117, 027201 (2016).
  * Chandran  _et al._ [2015] A. Chandran, C. R. Laumann, and V. Oganesyan,  (2015), arXiv:1509.04285 [cond-mat].
  * Roy [2022] S. Roy, Physical Review B 106, L140204 (2022).
  * Roy and Logan [2021a] S. Roy and D. E. Logan, arXiv:2106.09036 [cond-mat, physics:quant-ph] (2021a), arXiv: 2106.09036.
  * Macé  _et al._ [2019] N. Macé, F. Alet, and N. Laflorencie, Physical Review Letters 123, 180601 (2019).
  * Tikhonov and Mirlin [2021] K. S. Tikhonov and A. D. Mirlin, arXiv:2102.05930 [cond-mat] (2021), arXiv: 2102.05930.
  * Roy  _et al._ [2019a] S. Roy, J. T. Chalker, and D. E. Logan, Physical Review B 99, 104206 (2019a), arXiv: 1812.06101.
  * Roy and Logan [2021b] S. Roy and D. E. Logan, Physical Review B 104, 174201 (2021b).
  * De Tomasi  _et al._ [2020] G. De Tomasi, I. M. Khaymovich, F. Pollmann, and S. Warzel, arXiv:2011.03048 [cond-mat, physics:quant-ph] (2020), arXiv: 2011.03048.
  * Szoldra  _et al._ [2021] T. Szoldra, P. Sierant, K. Kottmann, M. Lewenstein, and J. Zakrzewski, arXiv:2106.01811 [cond-mat, physics:nlin, physics:physics, physics:quant-ph] (2021), arXiv: 2106.01811.
  * Morningstar  _et al._ [2022] A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Physical Review B 105, 174205 (2022).
  * Roy  _et al._ [2019b] S. Roy, D. E. Logan, and J. T. Chalker, Physical Review B 99, 220201 (2019b).
  * Goblot  _et al._ [2020] V. Goblot, A. Štrkalj, N. Pernet, J. L. Lado, C. Dorow, A. Lemaître, L. L. Gratiet, A. Harouri, I. Sagnes, S. Ravets, A. Amo, J. Bloch, and O. Zilberberg, Nature Physics 16, 832–836 (2020), arXiv: 1911.07809.
  * Logan and Welsh [2019] D. E. Logan and S. Welsh, Physical Review B 99, 045131 (2019), arXiv: 1806.01688.
  * Thiery  _et al._ [2018] T. Thiery, F. Huveneers, M. Müller, and W. De Roeck, Physical Review Letters 121, 140601 (2018), arXiv: 1706.09338.
  * Prelovšek  _et al._ [2018] P. Prelovšek, O. S. Barišić, and M. Mierzejewski, Physical Review B 97, 035104 (2018).
  * Filoche and Mayboroda [2012] M. Filoche and S. Mayboroda, Proceedings of the National Academy of Sciences 109, 14761 (2012).
  * Arnold  _et al._ [2018] D. Arnold, G. David, M. Filoche, D. Jerison, and S. Mayboroda, arXiv:1711.04888 [math] (2018), arXiv: 1711.04888.
  * Arnold  _et al._ [2015] D. L. Arnold, G. David, D. Jerison, S. Mayboroda, and M. Filoche, “The effective confining potential of quantum states in disordered media,”  (2015).
  * Pelletier  _et al._ [2021] P. Pelletier, D. Delande, V. Josse, A. Aspect, S. Mayboroda, D. Arnold, and M. Filoche, arXiv:2111.13155 [quant-ph] (2021), arXiv: 2111.13155.
  * Shamailov  _et al._ [2021] S. Shamailov, D. J. Brown, T. A. Haase, and M. Hoogerland, SciPost Physics Core 4, 017 (2021).
  * Li  _et al._ [2017] C.-K. Li, M. Piccardo, L.-S. Lu, S. Mayboroda, L. Martinelli, J. Peretti, J. S. Speck, C. Weisbuch, M. Filoche, and Y.-R. Wu, Physical Review B 95, 144206 (2017).
  * Kirkpatrick [1973] S. Kirkpatrick, Reviews of modern physics 45, 574 (1973).
  * Balasubramanian  _et al._ [2020] S. Balasubramanian, Y. Liao, and V. Galitski, Physical Review B 101, 014201 (2020).
  * Bel-Hadj-Aissa  _et al._ [2021] G. Bel-Hadj-Aissa, M. Gori, R. Franzosi, and M. Pettini, Journal of Statistical Mechanics: Theory and Experiment 2021, 023206 (2021).
  * Gori  _et al._ [2022] M. Gori, R. Franzosi, G. Pettini, and M. Pettini, Journal of Physics A: Mathematical and Theoretical (2022), 10.1088/1751-8121/ac7f09, arXiv:2207.05077 [cond-mat, physics:math-ph, physics:physics].
  * Sale  _et al._ [2022] N. Sale, J. Giansiracusa, and B. Lucini, Physical Review E 105, 024121 (2022).
  * Donato  _et al._ [2016] I. Donato, M. Gori, M. Pettini, G. Petri, S. De Nigris, R. Franzosi, and F. Vaccarino, Physical Review E 93, 052138 (2016).
  * Herviou and Bardarson [2020] L. Herviou and J. H. Bardarson, Physical Review B 101, 220201 (2020), arXiv: 2004.02903.
  * Lemut  _et al._ [2020] G. Lemut, M. Pacholski, O. Ovdat, A. Grabsch, J. Tworzydło, and C. Beenakker, Physical Review B 101, 081405 (2020).
  * Zomorodian and Carlsson [2005] A. Zomorodian and G. Carlsson, Discrete I& Computational Geometry 33, 249–274 (2005).
  * He  _et al._ [2022] Y. He, S. Xia, D. G. Angelakis, D. Song, Z. Chen, and D. Leykam, arXiv:2204.13276 [cond-mat, physics:physics, physics:quant-ph] (2022), arXiv: 2204.13276.
  * Cole  _et al._ [2021] A. Cole, G. J. Loges, and G. Shiu, Physical Review B 104, 104426 (2021).
  * Olsthoorn and Balatsky [2021] B. Olsthoorn and A. V. Balatsky, arXiv:2110.10214 [cond-mat, physics:quant-ph] (2021), arXiv: 2110.10214.
  * Filoche  _et al._ [2021] M. Filoche, S. Mayboroda, and T. Tao, Journal of Mathematical Physics 62, 041902 (2021).
  * Lyra  _et al._ [2015] M. L. Lyra, S. Mayboroda, and M. Filoche, EPL (Europhysics Letters) 109, 47001 (2015).
  * Murphy  _et al._ [2011] N. C. Murphy, R. Wortis, and W. A. Atkinson, Physical Review B 83, 184206 (2011), arXiv:1011.0659 [cond-mat].
  * Otter  _et al._ [2017] N. Otter, M. A. Porter, U. Tillmann, P. Grindrod, and H. A. Harrington, EPJ Data Science 6, 1 (2017).
  * Cohen-Steiner  _et al._ [2007] D. Cohen-Steiner, H. Edelsbrunner, and J. Harer, Discrete I& Computational Geometry 37, 103–120 (2007).
  * Bauer and Lesnick [2015] U. Bauer and M. Lesnick, Journal of Computational Geometry , 162 (2015), arXiv: 1311.3681.
  * Carlsson  _et al._ [2009] G. Carlsson, V. De Silva, and D. Morozov, in  _Proceedings of the twenty-fifth annual symposium on Computational geometry_ (2009) pp. 247–256.
  * Devakul and Singh [2015] T. Devakul and R. R. Singh, Physical Review Letters 115, 187201 (2015).
  * Doggen  _et al._ [2018] E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Physical Review B 98, 174202 (2018).
  * Serbyn and Moore [2016] M. Serbyn and J. E. Moore, Physical Review B 93, 041424 (2016), arXiv: 1508.07293.
  * Sierant and Zakrzewski [2019] P. Sierant and J. Zakrzewski, Physical Review B 99, 104205 (2019), arXiv: 1808.02795.
  * Jaquette and Schweinhart [2020] J. Jaquette and B. Schweinhart, Communications in Nonlinear Science and Numerical Simulation 84, 105163 (2020).
  * Bubenik  _et al._ [2020] P. Bubenik, M. Hull, D. Patel, and B. Whittle, Inverse Problems 36, 025008 (2020), arXiv: 1905.13196.
  * Bollobás [1998] B. Bollobás, in  _Modern graph theory_ (Springer, 1998) pp. 215–252.
  * Alt  _et al._ [2021a] J. Alt, R. Ducatez, and A. Knowles, Communications in Mathematical Physics 388, 507–579 (2021a).
  * Alt  _et al._ [2021b] J. Alt, R. Ducatez, and A. Knowles,  (2021b), arXiv:2109.03227 [math-ph].
  * Skraba  _et al._ [2020] P. Skraba, G. Thoppe, and D. Yogeshwaran,  (2020), arXiv:1701.00239 [math].
  * Myers  _et al._ [2019] A. Myers, E. Munch, and F. A. Khasawneh, Physical Review E 100, 022314 (2019).
  * Tran  _et al._ [2021] Q. H. Tran, M. Chen, and Y. Hasegawa, Physical Review E 103, 052127 (2021).
  * Saucan [2020] E. Saucan, arXiv preprint arXiv:2003.03844 (2020).
  * De Silva and Carlsson [2004] V. De Silva and G. E. Carlsson, in  _PBG_ (2004) pp. 157–166.
  * Shinaoka  _et al._ [2020] H. Shinaoka, D. Geffroy, M. Wallerberger, J. Otsuki, K. Yoshimi, E. Gull, and J. Kuneš, SciPost Physics 8, 012 (2020).