# The Quadratic Wasserstein Metric With Squaring Scaling For Seismic Velocity
Inversion

Zhengyang Li  Department of Mathematical Sciences, Tsinghua University,
Beijing, China 100084.  Yijia Tang  School of Mathematical Sciences, Shanghai
Jiao Tong University, Shanghai, China 200240.  Jing Chen  Hao Wu

###### Abstract

The quadratic Wasserstein metric has shown its power in measuring the
difference between probability densities, which benefits optimization
objective function with better convexity and is insensitive to data noise.
Nevertheless, it is always an important question to make the seismic signals
suitable for comparison using the quadratic Wasserstein metric. The squaring
scaling is worth exploring since it guarantees the convexity caused by data
shift. However, as mentioned in [Commun. Inf. Syst., 2019, 19:95-145], the
squaring scaling may lose uniqueness and result in more local minima to the
misfit function. In our previous work [J. Comput. Phys., 2018, 373:188-209],
the quadratic Wasserstein metric with squaring scaling was successfully
applied to the earthquake location problem. But it only discussed the inverse
problem with few degrees of freedom. In this work, we will present a more in-
depth study on the combination of squaring scaling technique and the quadratic
Wasserstein metric. By discarding some inapplicable data, picking seismic
phases, and developing a new normalization method, we successfully invert the
seismic velocity structure based on the squaring scaling technique and the
quadratic Wasserstein metric. The numerical experiments suggest that this
newly proposed method is an efficient approach to obtain more accurate
inversion results.

Keywords: Optimal Transport, Wasserstein metric, Waveform inversion, Seismic
velocity inversion, Squaring Scaling.

††∗ Corresponding author.

##  1 Introduction

Full waveform inversion (FWI) has been receiving wide attention in recent
years [9, 14, 22, 32, 36, 37] due to its high-resolution imaging in
geophysical properties. Generally, it can be formulated as a PDE constrained
optimization problem in mathematics, which consists of two parts [31]: the
forward modeling of seismic wavefield, and the optimization problem searching
for suitable model parameters to minimize the mismatch between the predicted
and observed seismic signals. In previous decades, limited by the computing
power, most tomography methods were based on the ray theory, which ignores
finite frequency phenomena such as wave-front healing and scattering [15], and
thus results in low-resolution inverion results. With the rapid development of
computing power and the forward modeling method, more accurate synthetic
signals could be computed by directly simulating seismic wave propagation.
This makes it possible to obtain high-resolution results by FWI, which could
provide important information for seismic hazard assessment [28] and
exploration geophysics [31].

The $L^{2}$ metric-based model is the simplest and most common FWI. However,
it suffers from the well-known cycle skipping problem [31] that the solution
may be trapped in the local minima during the iteration, leading to incorrect
inversion results. The quadratic Wasserstein metric ($W_{2}$) from the Optimal
transport (OT) theory [29, 30] seems to be a solution to the above problem. It
measures the difference between two probability distributions by minimizing
the transport cost from one distribution to the other, which is insensitive to
the data noise and keeps convexity to the data shift, dilation, and partial
amplitude change [10, 11]. The number of local minima of the FWI model based
on this metric is therefore significantly reduced. Thus, it is favored by
researchers and has been widely applied to the earthquake location and seismic
tomography [5, 10, 11, 12, 13, 36, 37, 38]. In applying the quadratic
Wasserstein metric to the seismic inverse problem, there is a critical
problem. The quadratic Wasserstein metric compares the normalized and
nonnegative data while the seismic signal does not meet this requirement.
Thus, various techniques are developed to deal with this problem, e.g., linear
scaling [36], squaring scaling [5], and exponential scaling [26]. Among all
these methods, squaring scaling is considered to maintain the convexity of the
optimization objective function. But this method seems to lose uniqueness and
result in additional minima. This may be the reason why we haven’t seen the
application of squaring scaling and quadratic Wasserstein metric to the
velocity inversion problem. Moreover, there are also some other metrics based
on the OT theory, e.g., the WFR metric and the KR norm, which have been
successfully applied to the seismic inverse problem [23, 24, 38].

In our previous work [5], the quadratic Wasserstein metric with squaring
scaling is successfully applied to the earthquake location problem. The
squaring scaling ensures the differentiability and nice convexity property,
leading to a large convergent domain and accurate inversion results. However,
it is still a challenging problem for velocity inversion with a large number
of degrees of freedom since the squaring scaling may lose uniqueness and
result in additional local minima to the misfit function [12]. In this work,
we would like to provide a comprehensive approach to the seismic velocity
inversion based on squaring scaling and the quadratic Wasserstein metric. The
key ingredient of this work consists of two parts. First, for seismic velocity
inversion, the fundamental geophysical characteristic of seismic signals
should be taken into account. For example, certain erroneous seismic signals
and multi-arrival seismic signals, which have destructive effects on the
inverse process, should be deleted in the preprocessing stage. Moreover, a
more accurate optimal transport map can be obtained by picking appropriate
seismic phases. Secondly, a new normalization method is developed to obtain a
more accurate optimal transport map for the squared seismic signals. From
this, we can calculate better sensitivity kernels, which are more consistent
with physical intuition.

The rest of the paper is organized as follows. In Section 2, we briefly review
the mathematical formula of seismic velocity inversion and the basics of the
quadratic Wasserstein metric. We discuss important issues in the inversion and
present detailed implementations in Section 3. Meanwhile, we illustrate the
necessity of our method by some toy models. In Section 4, the numerical
experiments are provided to demonstrate the effectiveness and efficiency of
our method. Finally, we conclude the paper in Section 5.

##  2 The quadratic Wasserstein metric and seismic velocity inversion

We review the full waveform seismic tomography and the adjoint state method in
this section. The mathematical formulation of seismic velocity inversion can
be written as the PDE constrained optimization problem,

$c_{T}(\boldsymbol{x})=\operatorname*{argmin}_{c(\boldsymbol{x})}\Xi(c(\boldsymbol{x})),\quad\Xi(c(\boldsymbol{x}))=\sum_{i=1}^{N}\sum_{j=1}^{M}\chi_{ij}(c(\boldsymbol{x})),$
(2.1)

where index $(i,j)$ indicates the source-receiver pair. We used $N$ seismic
events, and considered $M$ seismic signals for each event. Correspondingly,
the misfit function $\chi_{ij}$ is defined as

$\chi_{ij}(c(\boldsymbol{x}))=\mathcal{D}(s_{ij}(t;c(\boldsymbol{x})),d_{ij}(t)).$
(2.2)

Here, $\mathcal{D}$ is the distance function that measures the difference
between the real seismic signal $d_{ij}(t)$ and the synthetic signal
$s_{ij}(t;c(\boldsymbol{x}))$, which can be regarded as the solution

$d_{ij}(t)=u_{i}(\boldsymbol{\eta}_{j},t;c_{T}(\boldsymbol{x})),\quad
s_{ij}(t;c(\boldsymbol{x}))=u_{i}(\boldsymbol{\eta}_{j},t;c(\boldsymbol{x})),$
(2.3)

of the following acoustic wave equation with the initial boundary condition

$\displaystyle\frac{\partial^{2}u_{i}(\boldsymbol{x},t;c(\boldsymbol{x}))}{\partial
t^{2}}=\nabla\cdot\left(c^{2}(\boldsymbol{x})\nabla
u_{i}(\boldsymbol{x},t;c(\boldsymbol{x}))\right)+R(t-\tau_{i})\delta(\boldsymbol{x}-\boldsymbol{\xi}_{i}),\quad\boldsymbol{x}\in\Omega,t>0,$
(2.4) $\displaystyle u_{i}(\boldsymbol{x},0;c(\boldsymbol{x}))=\frac{\partial
u_{i}(\boldsymbol{x},0;c(\boldsymbol{x}))}{\partial
t}=0,\quad\boldsymbol{x}\in\Omega,$ (2.5)
$\displaystyle\boldsymbol{n}\cdot\left(c^{2}(\boldsymbol{x})\nabla
u_{i}(\boldsymbol{x},t;c(\boldsymbol{x}))\right)=0,\quad\boldsymbol{x}\in\partial\Omega,t>0.$
(2.6)

Here, the locations of the earthquake and receiver station are
$\boldsymbol{\xi}_{i}$ and $\boldsymbol{\eta}_{j}$, the origin time of the
earthquake is $\tau_{i}$. The seismic rupture is modeled by the point source
$\delta(\boldsymbol{x}-\boldsymbol{\xi})$ since its scale is much smaller
compared to the scale of seismic wave propagation [1, 20]. And the source time
function is simplified as the Ricker wavelet

$R(t)=A\left(1-2\pi^{2}f_{0}^{2}t^{2}\right)e^{-\pi^{2}f_{0}^{2}t^{2}},$ (2.7)

where $f_{0}$ denotes the dominant frequency, and $A$ is the normalization
factor. The outward unit normal vector to the simulation domain boundary
$\partial\Omega$ is $\boldsymbol{n}$. In practice, the perfectly matched layer
absorbing boundary condition [17] is used to deal with the propagation of
waves outside the area. In this section, we use the reflection boundary
condition to simplify the derivation.

######  Remark 1.

Here, we consider the trace by trace strategy [36] to apply the 1-D quadratic
Wasserstein metric to the waveform inversion. Considering the fact that
receiver stations are located far from each other on the geological scale,
this approach is more in line with physical reality and also easier in
mathematics.

###  2.1 The adjoint method

Below, we briefly review the adjoint method [11, 25] for solving the
optimization problems (2.1)-(2.7). For small perturbation of seismic velocity
structure $\delta c$, it causes the perturbation of the wavefield

$\delta
u_{i}(\boldsymbol{x},t;c(\boldsymbol{x}))=u_{i}(\boldsymbol{x},t;c+\delta
c)-u_{i}(\boldsymbol{x},t;c).$ (2.8)

For the sake of brevity, we will omit the parameter $c(\boldsymbol{x})$ of the
wavefield and the signals in the following. The perturbation $\delta
u_{i}(\boldsymbol{x},t)$ satisfies the equations

$\displaystyle\frac{\partial^{2}\delta u_{i}(\boldsymbol{x},t)}{\partial
t^{2}}=\nabla\cdot\left(c^{2}(\boldsymbol{x})\nabla\delta
u_{i}(\boldsymbol{x},t)\right)$ (2.9) $\displaystyle\quad\quad\quad\quad\quad\
\ +\nabla\cdot\left(\left(2c(\boldsymbol{x})+\delta
c(\boldsymbol{x})\right)\delta c(\boldsymbol{x})\nabla(u_{i}+\delta
u_{i})(\boldsymbol{x},t)\right),\quad\boldsymbol{x}\in\Omega,$
$\displaystyle\delta u_{i}(\boldsymbol{x},0)=\frac{\partial\delta
u_{i}(\boldsymbol{x},0)}{\partial t}=0,\quad\boldsymbol{x}\in\Omega,$ (2.10)
$\displaystyle\boldsymbol{n}\cdot\left(c^{2}(\boldsymbol{x})\nabla\delta
u_{i}(\boldsymbol{x},t)+\left(2c(\boldsymbol{x})+\delta
c(\boldsymbol{x})\right)\delta c(\boldsymbol{x})\nabla(u_{i}+\delta
u_{i})(\boldsymbol{x},t)\right)=0,\quad\boldsymbol{x}\in\partial\Omega.$
(2.11)

Multiply test function $w_{i}(\boldsymbol{x},t)$ on equation (2.9) and
integrate it on $\Omega\times[0,t_{f}]$ for sufficient large time $t_{f}$.
Using integration by parts yields

$\int_{0}^{t_{f}}\int_{\Omega}\frac{\partial^{2}w_{i}}{\partial t^{2}}\delta
u_{i}\mathrm{d}\boldsymbol{x}\mathrm{d}t-\int_{\Omega}\left.\frac{\partial
w_{i}}{\partial t}\delta
u_{i}\right|_{t=t_{f}}\mathrm{d}\boldsymbol{x}+\int_{\Omega}\left.w_{i}\frac{\partial\delta
u_{i}}{\partial t}\right|_{t=t_{f}}\mathrm{d}\boldsymbol{x}\\\
=\int_{0}^{t_{f}}\int_{\Omega}\nabla\cdot(c^{2}\nabla w_{i})\delta
u_{i}\mathrm{d}\boldsymbol{x}\mathrm{d}t-\int_{0}^{t_{f}}\int_{\partial\Omega}\boldsymbol{n}\cdot(c^{2}\nabla
w_{i})\delta
u_{i}\mathrm{d}\zeta\mathrm{d}t-\int_{0}^{t_{f}}\int_{\Omega}\left(2c+\delta
c\right)\delta c\nabla w_{i}\cdot\nabla(u_{i}+\delta
u_{i})\mathrm{d}\boldsymbol{x}\mathrm{d}t\\\
\approx\int_{0}^{t_{f}}\int_{\Omega}\nabla\cdot(c^{2}\nabla w_{i})\delta
u_{i}\mathrm{d}\boldsymbol{x}\mathrm{d}t-\int_{0}^{t_{f}}\int_{\partial\Omega}\boldsymbol{n}\cdot(c^{2}\nabla
w_{i})\delta
u_{i}\mathrm{d}\zeta\mathrm{d}t-\int_{0}^{t_{f}}\int_{\Omega}2c\delta c\nabla
w_{i}\cdot\nabla u_{i}\mathrm{d}\boldsymbol{x}\mathrm{d}t,$ (2.12)

where the higher-order terms are ignored in the last step since we can
naturally assume that $\left\|\delta u_{i}\right\|\ll\left\|u_{i}\right\|$ and
$\left\|\delta c(\boldsymbol{x})\right\|\ll\left\|c(\boldsymbol{x})\right\|$.

On the one hand, the perturbation of misfit $\delta\chi_{ij}$ results from the
wave speed perturbation $\delta c(\boldsymbol{x})$, which writes

$\delta\chi_{ij}(c)=\mathcal{D}\big{(}s_{ij}(t)+\delta
s_{ij}(t),d_{ij}(t)\big{)}-\mathcal{D}\big{(}s_{ij}(t),d_{ij}(t)\big{)}\\\
\approx\langle Q_{ij}(t),\;\delta
s_{ij}(t)\rangle=\int_{0}^{t_{f}}Q_{ij}(t)\delta s_{ij}(t)\mathrm{d}t.$

Here, $Q_{ij}(t)$ indicates the Fréchet gradient of the distance $\mathcal{D}$
with respect to the synthetic data $s_{ij}(t)$:

$Q_{ij}(t)=\nabla_{s}\mathcal{D}(s,d)\big{|}_{s=s_{ij}(t),d=d_{ij}(t)},$
(2.13)

which will be specified later. Let $w_{i}(\boldsymbol{x},t)$ satisfy the
adjoint equation

$\displaystyle\frac{\partial^{2}w_{i}(\boldsymbol{x},t)}{\partial
t^{2}}=\nabla\cdot\left(c^{2}(\boldsymbol{x})\nabla
w_{i}(\boldsymbol{x},t)\right)+\sum_{j=1}^{M}Q_{ij}(t)\delta(\boldsymbol{x}-\boldsymbol{\eta}_{j}),\quad\boldsymbol{x}\in\Omega,$
(2.14) $\displaystyle w_{i}(\boldsymbol{x},t_{f})=\frac{\partial
w_{i}(\boldsymbol{x},t_{f})}{\partial t}=0,\quad\boldsymbol{x}\in\Omega,$
(2.15) $\displaystyle\boldsymbol{n}\cdot\left(c^{2}(\boldsymbol{x})\nabla
w_{i}(\boldsymbol{x},t)\right)=0,\quad\boldsymbol{x}\in\partial\Omega.$ (2.16)

Multiply $\delta u_{i}(\boldsymbol{x},t)$ on equation (2.14), integrate it on
$\Omega\times[0,t_{f}]$ and subtract (2.12) to obtain

$\sum_{j=1}^{M}\int_{0}^{t_{f}}Q_{ij}(t)\delta
s_{ij}(t)\mathrm{d}t=\sum_{j=1}^{M}\int_{0}^{t_{f}}\int_{\Omega}Q_{ij}(t)\delta(\boldsymbol{x}-\boldsymbol{\eta}_{j})\delta
u_{i}(\boldsymbol{x},t)\mathrm{d}t\\\
=-\int_{0}^{t_{f}}\int_{\Omega}2c(\boldsymbol{x})\delta
c(\boldsymbol{x})\nabla w_{i}(\boldsymbol{x},t)\cdot\nabla
u_{i}(\boldsymbol{x},t)\mathrm{d}\boldsymbol{x}\mathrm{d}t.$

The linear relationship between $\delta\Xi$ and $\delta c(\boldsymbol{x})$ is
established as

$\delta\Xi(c)=\sum_{i=1}^{N}\sum_{j=1}^{M}\delta\chi_{ij}(c)=\sum_{i=1}^{N}\int_{\Omega}K_{i}(\boldsymbol{x})\delta
c(\boldsymbol{x})\mathrm{d}\boldsymbol{x},$ (2.17)

where the sensitivity kernel of the $i$-th source for $c(\boldsymbol{x})$ is
defined as

$K_{i}(\boldsymbol{x})=-\int_{0}^{t_{f}}2c(\boldsymbol{x})\nabla
w_{i}(\boldsymbol{x},t)\cdot\nabla u_{i}(\boldsymbol{x},t)\mathrm{d}t.$ (2.18)

###  2.2 The quadratic Wasserstein metric

As we discussed at the beginning of this section, the synthetic signal
$s_{ij}(t)$ and real seismic signal $d_{ij}(t)$ are time series. As we know,
the quadratic Wasserstein metric between the 1-D probability density functions
has an analytic form [5, 29, 30, 36], i.e.,

$W_{2}^{2}(f,g)=\int_{0}^{t_{f}}\left|t-T(t)\right|^{2}f(t)\mathrm{d}t,\quad
T(t)=G^{-1}\left(F(t)\right).$ (2.19)

Here $f(t),\;g(t)$ are probability density functions defined on $[0,t_{f}]$
and $F(t),\;G(t)$ are cumulative density functions defined on $[0,t_{f}]$,

$F(t)=\int_{0}^{t_{f}}f(\tau)\mathrm{d}\tau,\quad
G(t)=\int_{0}^{t_{f}}g(\tau)\mathrm{d}\tau.$

Note that the seismic signals are not probability density functions. We need
to transform them into nonnegative and normalized functions for the quadratic
Wasserstein metric comparison. In other words, the misfit function defined in
(2.2) can be written as

$\chi_{ij}=\mathcal{D}(s_{ij}(t),d_{ij}(t))=W^{2}_{2}(\mathcal{P}(s_{ij}(t)),\mathcal{P}(d_{ij}(t))).$
(2.20)

The operator $\mathcal{P}$ converts the seismic signals into probability
density functions, including processing them into nonnegative and normalized
time series. In the later part, we will discuss this in detail. Thus, we can
obtain the expression of the Fréchet gradient [5, 36] mentioned in (2.13),

$\nabla_{s}\mathcal{D}(s,d)=\nabla_{f}W^{2}_{2}(f,g)|_{f=\mathcal{P}(s),g=\mathcal{P}(d)}\cdot\nabla_{s}\mathcal{P}(s)=\left\langle
2\int_{0}^{t}\tau-T(\tau)\mathrm{d}\tau,\nabla_{s}\mathcal{P}(s)\right\rangle.$
(2.21)

##  3 Data preprocessing and new normalization

In this section, we discuss two important issues when carrying out seismic
velocity inversion. First of all, when using real data for inversion, we do
not use all the data in each iteration. Some data, such as the case where the
direct wave and the reflected wave arrive simultaneously, are difficult to use
and can be ignored. In order to avoid the mismatch between different types of
seismic phases, we only retain the direct waves in the real seismic signals
and the synthetic signals. This processing procedure ensures reasonable
optimal transport maps and accurate sensitivity kernels. Secondly, we will
carefully design the operator $\mathcal{P}$ to get a better OT map $T$. In the
following, we will present detailed implementations and discussions.

###  3.1 Selecting source-receiver pairs and picking seismic phases

The complex subsurface structures, such as the velocity discontinuity
interfaces, may lead to different types of seismic phases, including the
direct wave and the reflected wave. These seismic waves propagate along
different wave paths and carry distinct underground structure information.
Sometimes, the direct wave and the reflected wave arrive simultaneously and
can not be distinguished, called the multipath phenomenon [27]. It is not
trivial to extract robust information from this kind of constraint. In
practice, these source-receiver pairs are always manually excluded to avoid
interference caused by unreliable constraints [3, 16]. We will also use this
strategy in this study.

From the perspective of signal processing, different phases of the real
seismic signal and the synthetic signal should be matched separately. If there
is a matching error, for example, part of the direct wave of the synthetic
signal is matched with part of the reflected wave of the real seismic signal,
it would lead to the optimal transport map being inconsistent with basic
seismic knowledge and further result in the artifacts in the sensitivity
kernel [10]. In particular, for the squaring scaling and quadratic Wasserstein
metric based seismic velocity inversion, this problem is more prominent. The
reason is that the quadratic Wasserstein metric requires mass conservation and
global match. When the masses of the real seismic signal and the synthetic
signal are unbalanced in the same phase, the mass transportation between
different phases will occur, causing the inconsistency between the OT map with
seismic reality. Moreover, the squaring scaling could further magnify the
problem. The idea of solving the above problems is also easy. By picking the
phases, we only match the same phases of the real seismic signals and the
synthetic signals. This is a common strategy in seismic inversion [21, 6], and
it can be achieved simply by calculating the arrival time of the direct phase
and the reflected phase [7, 33].

Figure 1: Illustration of the two-layer model. Left: the real seismic velocity
model with a high-velocity anomaly; Right: the initial velocity model. The
green inverted triangles indicate the receiver stations and the white stars
indicate the earthquakes. The specific source-receiver pair is highlighted by
the black star and inverted triangle. The cyan and tan dashed lines are the
direct wave path and the reflected wave path, respectively.

Next, we explain the necessity of the above-mentioned data preprocessing
method. The initial and real seismic velocity models are shown in Figure 1,
and the parameter settings can be found in Section 4.1. The main goal is to
detect the high-velocity anomaly above the Moho discontinuity.

Whether initial or real seismic velocity models, there are at least two paths
from the earthquake hypocenter to the receiver station: the direct wave (cyan
dashed lines) and the reflected wave (tan dashed lines). In the real seismic
velocity model, the wave amplitude of the direct wave signal is slightly
smaller since it partially reflects when passing through the high-velocity
anomaly. On the other hand, the reflected wave signal should be the same since
the velocity structure on the reflected wave path is the same in the initial
and real seismic velocity models, see Figure 1 for illustration.

Figure 2: Illustration of the Optimal Transport map between the real seismic
signal and synthetic signal (left) and the sensitivity kernel (right). The
mass transportation from the direct wave of the synthetic signal to the
reflected wave of the real seismic signal (within the green box of the upper
left subgraph) will cause artifacts in the sensitivity kernel, which arise
around the reflected wave path (the blue dashed lines of the upper right
subgraph). In the lower subgraphs, we can obtain the satisfactory OT map and
sensitivity kernel since only direct waves are picked.

The above difference between the real seismic signal and the synthetic signal
is further magnified by the squaring scaling. It leads to unreasonable mass
transportation from the direct wave of the synthetic signal to the reflected
wave of the real seismic signal (upper left subgraph of Figure 2). Therefore,
there will be artifacts in the sensitivity kernel $K_{i}(\boldsymbol{x})$, as
we illustrate in the upper right subgraph of Figure 2. On the other hand, if
we only consider the direct waves for inversion, the above-mentioned
difficulties will be easily solved, as we illustrate in the lower subgraphs of
Figure 2.

######  Remark 2.

In fact, the reflected wave signals are also important to constrain the
underground velocity structures [16]. The reflection phases can also be
similarly picked, processed, and used for inversion by our approach. However,
the utilization of the reflected wave is not trivial, and more technical
details are required in practice [2, 35, 39]. Thus, we will not discuss the
issues of the reflected wave in the following sections.

###  3.2 New normalization method

As it is well known, the quadratic Wasserstein metric measures the difference
between two probability density functions, which is not directly suitable for
seismic signals. Thus, some processing procedures, i.e., choosing an
appropriate operator $\mathcal{P}$ in (2.20) are required to convert the
seismic signals into probability density functions. Several different
approaches, e.g., linear scaling [36], squaring scaling [5], and exponential
scaling [26], have been proposed to address this issue. Among these methods,
the squaring scaling maintains convexity very well, and it is worthy of more
discussions.

The normalization operator with squaring scaling consists of two ingredients:
squaring seismic signal to ensure non-negativity and normalization to
guarantee the same mass. A natural approach is

$\mathcal{P}_{1}(s(t))=\frac{s^{2}(t)}{\left\|s^{2}(t)\right\|},$ (3.1)

in which

$\left\|s(t)\right\|=\int_{0}^{t_{f}}s(t)\mathrm{d}t.$

Substitute the above formula into equation (2.20), the form of the misfit
function is given by

$\chi=\mathcal{D}(s(t),d(t))=W^{2}_{2}\left(\frac{s^{2}(t)}{\left\|s^{2}(t)\right\|},\frac{d^{2}(t)}{\left\|d^{2}(t)\right\|}\right).$

Here the subscript indices $i$ and $j$ are dropped for simplicity. According
to the discussions in Section 2.2, we need to compute the inverse of the
following cumulative distribution function

$G(t)=\int_{0}^{t_{f}}\frac{d^{2}(t)}{\left\|d^{2}(t)\right\|}\mathrm{d}t.$

However, $G^{-1}(t)$ is not well defined when the real seismic signal $d(t)=0$
in certain interval. Correspondingly, there will be difficulties in the
computation of the misfit function.

In order to avoid the above-mentioned problem, we can make a slight upward
shift on the squared signal before the normalization, i.e.,

$\mathcal{P}_{2}(s(t))=\frac{s^{2}(t)+\varepsilon}{\left\|s^{2}(t)+\varepsilon\right\|}.$
(3.2)

Here $\varepsilon>0$ is a small parameter. However, the misfit function in
(2.20) with this normalization operator

$\chi=\mathcal{D}(s(t),d(t))=W^{2}_{2}\left(\frac{s^{2}(t)+\varepsilon}{\left\|s^{2}(t)+\varepsilon\right\|},\frac{d^{2}(t)+\varepsilon}{\left\|d^{2}(t)+\varepsilon\right\|}\right)$

still leads to unreasonable mass transportation (green box in the upper left
subgraph of Figure 3) since the additional mass does not equal

$\frac{\varepsilon}{\left\|s^{2}(t)+\varepsilon\right\|}\neq\frac{\varepsilon}{\left\|d^{2}(t)+\varepsilon\right\|}.$

This again leads to artifacts in the sensitivity kernel
$K_{i}(\boldsymbol{x})$ (upper right subgraph of Figure 3).

With a simple trick, we can solve the problem of unequal additional masses by
modifying the normalization operator as

$\mathcal{P}_{3}(s(t))=\frac{\frac{s^{2}(t)}{\left\|s^{2}(t)\right\|}+\varepsilon}{1+t_{f}\varepsilon}.$
(3.3)

We can clearly see that regardless of the values of $s(t)$ and $d(t)$, the
additional mass is $\frac{\varepsilon}{1+t_{f}\varepsilon}$. As a result, we
can avoid all the mentioned troubles. Both the OT map and the sensitivity
kernel are satisfactory, as we illustrate in the lower subgraphs of Figure 3.

######  Remark 3.

In the squaring scaling, a parameter $\varepsilon$ is added to avoid the
singularity. It is noted that large $\varepsilon$ could destroy the convexity
property. On the other hand, there will still be numerical singularities when
$\varepsilon$ is small. In practice, $\varepsilon$ is feasible in a relatively
large range, e.g., $10^{-4}\sim 10^{-2}$. In the following numerical
experiments, we select $\varepsilon=10^{-3}$.

Figure 3: Illustration of the Optimal Transport map between the real seismic
signal and synthetic signal (left) and the sensitivity kernel (right). In the
upper subgraphs, the newly created mass by the operator $\mathcal{P}_{2}$
could not be balanced, which leads to unreasonable mass transportation (upper
left) and artifacts in the sensitivity kernel (upper right). In the lower
subgraphs, we can obtain the satisfactory OT map and sensitivity kernel since
a new operator $\mathcal{P}_{3}$ is used.

##  4 Numerical Experiments

In this section, we present two numerical experiments to investigate the
validity of our inversion method based on the quadratic Wasserstein metric
with squaring scaling. We use the finite difference method to solve the
acoustic wave equation [8, 19, 36]. The perfectly matched layer boundary
condition [17] is applied to absorb the outgoing wave. The delta source
function is discretized by piecewise polynomial given in [34]

$\delta_{h}(x)=\left\\{\begin{array}[]{ll}\frac{1}{h}\left(1-\frac{5}{4}\left|\frac{x}{h}\right|^{2}-\frac{35}{12}\left|\frac{x}{h}\right|^{3}+\frac{21}{4}\left|\frac{x}{h}\right|^{4}-\frac{25}{12}\left|\frac{x}{h}\right|^{5}\right),&\left|x\right|\leq
h,\\\
\frac{1}{h}\left(-4+\frac{75}{4}\left|\frac{x}{h}\right|-\frac{245}{8}\left|\frac{x}{h}\right|^{2}+\frac{545}{24}\left|\frac{x}{h}\right|^{3}-\frac{63}{8}\left|\frac{x}{h}\right|^{4}+\frac{25}{24}\left|\frac{x}{h}\right|^{5}\right),&h<\left|x\right|\leq
2h,\\\
\frac{1}{h}\left(18-\frac{153}{4}\left|\frac{x}{h}\right|+\frac{255}{8}\left|\frac{x}{h}\right|^{2}-\frac{313}{24}\left|\frac{x}{h}\right|^{3}+\frac{21}{8}\left|\frac{x}{h}\right|^{4}-\frac{5}{24}\left|\frac{x}{h}\right|^{5}\right),&2h<\left|x\right|\leq
3h,\\\ 0,&\left|x\right|>3h.\end{array}\right.$

Here $h$ is related to the mesh size.

###  4.1 The Two-Layer Model

Consider the two-layer model in a bounded domain
$\Omega=[0,80\;km]\times[0,60\;km]$, which consists of the crust, the
uppermost mantle, and the Moho discontinuity at a depth of $30\;km$, see
Figure 1 for illustration. The real seismic velocity model includes a $+15\%$
high-velocity anomaly in the crust, given by

$c_{T}(x,z)=\left\\{\begin{array}[]{ll}6.67\ km/s,&(x,z)\in[35\ km,45\
km]\times[10\ km,20\ km],\\\ 8.1\ km/s,&z>30\ km,\\\ 5.8\
km/s,&others.\end{array}\right.$

Our goal is to perform the seismic velocity inversion to detect this high-
velocity anomaly. Correspondingly, the initial velocity model without high-
velocity anomaly is as follows

$c_{0}(x,z)=\left\\{\begin{array}[]{ll}5.8\ km/s,&z\leq 30\ km,\\\ 8.1\
km/s,&z>30\ km.\end{array}\right.$

The computational time interval is $[0\;s,21\;s]$. The inversion grid step is
$2\;km$ and the number of degrees of freedom amounts to $1200$. The space and
time steps in the forward simulation are $0.2\ km$ and $0.01\ s$,
respectively. The dominant frequency of the earthquakes in (2.7) is
$f_{0}=2\;Hz$. We randomly choose $25$ receiver stations deployed on the
surface and $80$ earthquakes distributed in the study region.

We then perform the seismic velocity inversion by using the quadratic
Wasserstein metric with squaring scaling. As a comparison, the inversion is
also performed with the traditional $L^{2}$ metric.

To quantitatively compare the results of different methods, we also compute
the relative model error

$RME=\frac{\int_{\Omega}|c_{k}(\boldsymbol{x})-c_{T}(\boldsymbol{x})|^{2}\mathrm{d}\boldsymbol{x}}{\int_{\Omega}|c_{0}(\boldsymbol{x})-c_{T}(\boldsymbol{x})|^{2}\mathrm{d}\boldsymbol{x}},$

and the relative misfit function

$RMF=\frac{\Xi(c_{k}(\boldsymbol{x}))}{\Xi(c_{0}(\boldsymbol{x}))},$

where $c_{k}(\boldsymbol{x})$ indicates the velocity model in the $k$-th
iteration.

Figure 4: The inversion results of the two-layer model. Upper subgraphs: the
result for $L^{2}$ metric after 20 steps (upper left); the convergent
trajectories of the relative model error (upper middle); the convergent
trajectories of the relative misfit function (upper right). In the middle and
the lower subgraphs, we present the results for the $W_{2}$ metric with the
operators $\mathcal{P}_{2}$ and $\mathcal{P}_{3}$, respectively. From left to
right, the inversion iteration steps are $20$, $40$, and $80$. All the results
are shown in the same color bar.

In Figure 4, we present the inversion results of $L^{2}$ metric and $W_{2}$
metric. Obviously, the $L^{2}$-based inversion could not capture the $+15\%$
high-velocity anomaly (upper left subgraph of Figure 4). Although the misfit
function decreases in the iteration (upper middle subgraph of Figure 4), the
model error increases (upper right subgraph of Figure 4).

In Figure 4 and Table 1, we also compare the inversion results of the
quadratic Wasserstein metric with different operators $\mathcal{P}_{2}$ and
$\mathcal{P}_{3}$. From the convergent trajectories (upper middle and upper
right subgraphs of Figure 4), we can see the relative model error and the
relative misfit function of the operator $\mathcal{P}_{3}$ both have a faster
descent rate than those of the operator $\mathcal{P}_{2}$. Quantitatively, we
can see from Table 1 that the operator $\mathcal{P}_{3}$ only needs half of
the iteration steps of the operator $\mathcal{P}_{2}$ to achieve almost the
same relative model error and relative misfit function. This significantly
saves the expensive computational cost of the seismic velocity inversion
problem. Finally, it can be seen from the middle and lower subgraphs of Figure
4, the velocity inversion results of the operator $\mathcal{P}_{3}$ are
significantly better than those of the operator $\mathcal{P}_{2}$ under the
same iteration steps. The above discussions show that our approach has higher
efficiency and better inversion results.

Table 1: The two-layer model. Relative Model Error and Relative Misfit Function of $W_{2}$ with the operators $\mathcal{P}_{2}$ and $\mathcal{P}_{3}$ in $20$, $40$ and $80$ iteration steps, respectively. Iteration Steps | Relative Model Error | Relative Misfit Function
---|---|---
$W_{2}$ with $P_{2}$ |  $W_{2}$ with $P_{3}$ |  $W_{2}$ with $P_{2}$ |  $W_{2}$ with $P_{3}$
$20$ | $3.69\times 10^{-1}$ | $2.15\times 10^{-1}$ | $4.99\times 10^{-3}$ | $6.90\times 10^{-4}$
$40$ | $2.23\times 10^{-1}$ | $1.04\times 10^{-1}$ | $3.61\times 10^{-4}$ | $9.41\times 10^{-5}$
$80$ | $8.35\times 10^{-2}$ | $3.04\times 10^{-2}$ | $1.25\times 10^{-5}$ | $2.75\times 10^{-6}$

###  4.2 The Crustal Root Model

Let us consider the crustal root model, a kind of subsurface structure usually
found along the orogen. This model consists of the two-layered crust divided
by the Conrad discontinuity. A dipping and discontinuous Moho interface
separates the crust and the mantle. The depiction of these tectonic features
helps us better understand the forming of the old mountains. In mathematics,
we consider this three-layer model in the bounded domain
$\Omega=[0,80\;km]\times[0,80\;km]$. Three layers are divided by the Conrad
discontinuity at $20\;km$ depth and the Moho discontinuity whose location
$(x,L(x))$ is formulated with a quadratic function is given by

$L(x)=\left\\{\begin{array}[]{ll}36+\frac{25}{1600}x^{2}\ km,&0\ km\leq x\leq
40\ km,\\\ 36\ km,&40\ km<x\leq 80\ km.\end{array}\right.$

The seismic wave speed at each layer refers to the AK135 model [18],
generating the real seismic velocity model (Figure 5, left)

$c_{T}(x,z)=\left\\{\begin{array}[]{ll}5.8\ km/s,&z\leq 20\ km,\\\ 6.5\
km/s,&20\ km<z\leq L(x),\\\ 8.04\ km/s,&others.\end{array}\right.$

Our goal is to perform the seismic velocity inversion to detect this crustal
root. Correspondingly, the initial velocity model (Figure 5, right) without
crustal root anomaly is as follows

$c_{0}(x,z)=\left\\{\begin{array}[]{ll}5.8\ km/s,&z\leq 20\ km,\\\ 6.5\
km/s,&20\ km<z\leq 36\ km\\\ 8.04\ km/s,&others.\end{array}\right.$

The computational time interval is $[0\;s,21\;s]$. The inversion grid step is
$2\;km$ and the number of degrees of freedom amounts to $1600$. The space and
time steps in the forward simulation are $0.2\ km$ and $0.01\ s$,
respectively. The dominant frequency of the earthquakes in (2.7) is
$f_{0}=2\;Hz$. We randomly choose $40$ receiver stations deployed on the
surface and $80$ earthquakes distributed in the study region.

Figure 5: Illustration of the crustal root model. Left: the real seismic
velocity model. Right: the initial velocity model. The green inverted
triangles and the white stars indicate the receiver stations and the
earthquakes, respectively.

Similar to subsection 4.1, we present the inversion results of $L^{2}$ metric
and $W_{2}$ metric with the operators $\mathcal{P}_{2}$ and $\mathcal{P}_{3}$
in Figure 6. Obviously, the $L^{2}$-based inversion could not capture the
crustal root structure. The relative model error and the relative misfit
function with respect to different normalization operators are given in Table
2. Correspondingly, the convergent trajectories are output in the upper middle
and upper right subgraphs of Figure 6. In the middle and lower subgraphs of
Figure 6, the inversion results are also presented. From which, we can draw
the same conclusions as those in subsection 4.1.

Figure 6: The inversion results of the crustal root model. Upper subgraphs: the result for $L^{2}$ metric after 40 steps (upper left); the convergent trajectories of the relative model error (upper middle); the convergent trajectories of the relative misfit function (upper right). In the middle and the lower subgraphs, we present the results for the $W_{2}$ metric with the operators $\mathcal{P}_{2}$ and $\mathcal{P}_{3}$, respectively. From left to right, the inversion iteration steps are $40$, $80$, and $160$. All the results are shown in the same color bar. Table 2: The crustal root model. Relative Model Error and Relative Misfit Function of $W_{2}$ with the operators $\mathcal{P}_{2}$ and $\mathcal{P}_{3}$ in $40$, $80$ and $160$ iteration steps, respectively. Iteration Steps | Relative Model Error | Relative Misfit Function
---|---|---
$W_{2}$ with $P_{2}$ |  $W_{2}$ with $P_{3}$ |  $W_{2}$ with $P_{2}$ |  $W_{2}$ with $P_{3}$
$40$ | $6.43\times 10^{-1}$ | $5.59\times 10^{-1}$ | $5.47\times 10^{-3}$ | $6.35\times 10^{-4}$
$80$ | $5.37\times 10^{-1}$ | $4.68\times 10^{-1}$ | $7.83\times 10^{-4}$ | $1.74\times 10^{-4}$
$160$ | $4.32\times 10^{-1}$ | $3.99\times 10^{-1}$ | $1.33\times 10^{-4}$ | $6.11\times 10^{-5}$

##  5 Conclusion

What we have seen from the above is the solution to the problem that the
seismic velocity inversion based on squaring scaling and the quadratic
Wasserstein metric is difficult, as mentioned in [Commun. Inf. Syst., 2019,
19:95-145] and [Meth. Appl. Anal., 2019, 2:133-148]. We can not only solve the
seismic velocity inversion with a large number of degrees of freedom. By
introducing a better normalization operator, the convergence efficiency is
significantly improved. We would like to combine the above techniques with the
double-difference traveltime adjoint tomography [4], which has significant
advantages in real seismic data. This may result in a more robust and reliable
seismic velocity inversion method. We are currently investigating this
interesting topic and hope to report this in an independent publication.

## Acknowledgments

This work was supported by National Natural Science Foundation of China (Grant
No. 11871297) and Tsinghua University Initiative Scientific Research Program.

## References

  * [1] K. Aki and P.G. Richards, Quantitative Seismology: Theory and Methods volume II, W.H. Freeman & Co (Sd), 1980.
  * [2] R. Brossier, S. Operto and J. Virieux, Velocity model building from seismic reflection data by full-waveform inversion, Geophysical Prospecting, 63(2), 354-367, 2015.
  * [3] S.-J. Chang and C.-E. Baag, Crustal Structure in Southern Korea from Joint Analysis of Regional Broadband Waveforms and Travel Times, Bulletin of the Seismological Society of America, 96(3), 856–870, 2006.
  * [4] J. Chen, G.X. Chen, H. Wu, J.Y. Yao and P. Tong, Adjoint tomography of NE Japan revealed by common-source double-difference traveltime data, preprint.
  * [5] J. Chen, Y.F. Chen, H. Wu and D.H. Yang, The quadratic Wasserstein metric for Earthquake Location, J. Comput. Phys., 373, 188-209, 2018.
  * [6] Y. Chen, J. Hill, W. Lei, M. Lefebvre, J. Tromp, E. Bozdag and D. Komatitsch, Automated time-window selection based on machine learning for full-waveform inversion, SEG Technical Program Expanded Abstracts, 1604-1609, 2017.
  * [7] R. Chu, S. Ni, A. Pitarka and D.V. Helmberger, Inversion of Source Parameters for Moderate Earthquakes Using Short-Period Teleseismic P Waves, Pure and Applied Geophysics, 171(7), 1329–1341, 2014.
  * [8] M.A. Dablain, The application of high-order differencing to the scalar wave equation, Geophysics, 51(1), 54-66, 1986.
  * [9] M. Dunlop and Y.N. Yang, New likelihood functions and level-set prior for Bayesian full-waveform inversion, In SEG Technical Program Expanded Abstracts 2020, pages 825-829. Society of Exploration Geophysicists, 2020.
  * [10] B. Engquist and B.D. Froese, Application of the Wasserstein metric to seismic signals, Commun. Math. Sci., 12(5), 979-988, 2014.
  * [11] B. Engquist, B.D. Froese and Y.N. Yang, Optimal transport for seismic full waveform inversion, Commun. Math. Sci., 14(8), 2309-2330, 2016.
  * [12] B. Engquist and Y.N. Yang, Seismic imaging and optimal transport, Communications in Information and Systems, 19(2), 95-145, 2019.
  * [13] B. Engquist and Y.N. Yang, Seismic inversion and the data normalization for optimal transport, Methods and Applications of Analysis, 26(2), 133-148, 2019.
  * [14] B. Engquist and Y.N. Yang, Optimal Transport Based Seismic Inversion:Beyond Cycle Skipping, Comm. Pure Appl. Math.. https://doi.org/10.1002/cpa.21990, 2021.
  * [15] S.-H. Hung, F.A.Dahlen and G. Nolet, Wavefront healing: a banana–doughnut perspective, Geophys. J. Int., 146(2), 289-312, 2001.
  * [16] X. Huang, D. Yang, P. Tong, J. Badal and Q. Liu, Wave equation-based reflection tomography of the 1992 Landers earthquake area, Geophysical Research Letters, 43(5), 1884-1892, 2016.
  * [17] D. Komatitsch and J. Tromp, A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation, Geophys. J. Int., 154(1), 146-153, 2003.
  * [18] B.L.N. Kennett, E.R. Engdahl and R. Buland, Constraints on seismic velocities in the Earth from traveltimes, Geophys. J. Int., 122(1), 108-124, 1995.
  * [19] J.S. Li, D.H. Yang, H. Wu and X. Ma, A low-dispersive method using the high-order stereo-modelling operator for solving 2-D wave equations, Geophys. J. Int., 210(3), 1938-1964, 2017.
  * [20] R. Madariaga, Seismic Source Theory, in Treatise on Geophysics (Second Edition), pp. 51-71, ed. Gerald, S., Elsevier B.V., 2015.
  * [21] A. Maggi, C. Tape, M. Chen, D. Chao and J. Tromp, An automated time-window selection algorithm for seismic tomography, Geophys. J. Int., 178(1), 257–281, 2009.
  * [22] L. Métivier, R. Brossier, J. Virieux and S. Operto, Full waveform inversion and the truncated Newton method, SIAM Journal on Scientific Computing, 35(2), B401-B437, 2013.
  * [23] L. Metivier, R. Brossier, Q. Merigot, E. Oudet and J. Virieux, Measuring the misfit between seismograms using an optimal transport distance: application to full waveform inversion, Geophys. J. Int., 205(1), 345–377, 2016.
  * [24] L. Metivier, R. Brossier, Q. Merigot, E. Oudet and J. Virieux, An optimal transport approach for seismic tomography: application to 3D full waveform inversion, Inverse Problems, 32(11), 115008, 2016.
  * [25] R.-E. Plessix, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, Geophys. J. Int., 167(2), 495–503, 2006.
  * [26] L. Qiu, J. Ramos-Martínez, A. Valenciano, Y. Yang and B. Engquist, Full-waveform inversion with an exponentially encoded optimal-transport norm, In SEG Technical Program Expanded Abstracts 2017, pages 1286–1290. Society of Exploration Geophysicists, 2017.
  * [27] N. Rawlinson, M. Sambridge and J. Hauser, Multipathing, reciprocal traveltime fields and raylets, Geophys. J. Int., 181(2), 1077-1092, 2010.
  * [28] C. Tape, Q.Y. Liu, A. Maggi and J. Tromp, Seismic tomography of the southern California crust based on spectral-element and adjoint methods, Geophys. J. Int., 180(1), 433-462, 2010.
  * [29] C. Villani, Topics in Optimal transportation, Graduate Studies in Mathematics, American Mathematical Society, 2003.
  * [30] C. Villani, Optimal Transport: Old and New, Springer Science $\&$ Business Media, 2008.
  * [31] J. Virieux and S. Operto, An overview of full-waveform inversion in exploration geophysics, Geophysics, 74(6), WCC1-WCC26, 2009.
  * [32] J. Wang, D.H. Yang, H. Jing and H. Wu, Full waveform inversion based on the ensemble Kalman filter method using uniform sampling without replacement, Science Bulletin, 64(5), 321-330, 2019.
  * [33] S. Wang and H. Tkalčić, Seismic event coda-correlation: Toward global coda-correlation tomography, Journal of Geophysical Research: Solid Earth, 125, e2019JB018848, 2020.
  * [34] X. Wen, High Order Numerical Quadratures to One Dimensional Delta Function Integrals, SIAM J. Sci. Comput., 30(4), 1825-1846, 2008.
  * [35] S. Xu, D.L. Wang, F. Chen, G. Lambaré and Y. Zhang, Inversion on Reflected Seismic Wave, SEG Technical Program Expanded Abstracts, 1-7, 2012.
  * [36] Y.N. Yang, B. Engquist, J.Z. Sun and B.D. Froese, Application of Optimal Transport and the Quadratic Wasserstein Metric to Full-Waveform Inversion, Geophysics, 83(1), R43-R62, 2018.
  * [37] Y.N. Yang and B. Engquist, Analysis of optimal transport and related misfit functions in full-waveform inversion, Geophysics, 83(1), A7-A12, 2018.
  * [38] D. T. Zhou, J. Chen, H. Wu, D.H. Yang and L.Y. Qiu, The Wasserstein-Fisher-Rao metric for waveform based earthquake location, J. Comp. Math., accepted.
  * [39] W. Zhou, R. Brossier, S. Operto and J. Virieux, Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation, Geophys. J. Int., 202(3), 1535-1554, 2015.


# Delocalization-localization dynamical phase transition of random walks on
graphs

Giorgio Carugno<EMAIL_ADDRESS>Department of Mathematics, King’s
College London, Strand, London WC2R 2LS, UK  Pierpaolo Vivo
<EMAIL_ADDRESS>Department of Mathematics, King’s College London,
Strand, London WC2R 2LS, UK  Francesco Coghi<EMAIL_ADDRESS>Nordita,
KTH Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg
12, SE-106 91 Stockholm, Sweden

###### Abstract

We consider random walks evolving on two models of connected and undirected
graphs and study the exact large deviations of a local dynamical observable.
We prove, in the thermodynamic limit, that this observable undergoes a first-
order dynamical phase transition (DPT). This is interpreted as a ‘co-
existence’ of paths in the fluctuations that visit the highly connected bulk
of the graph (delocalization) and paths that visit the boundary
(localization). The methods we used also allow us to characterize analytically
the scaling function that describes the finite size crossover between the
localized and delocalized regimes. Remarkably, we also show that the DPT is
robust with respect to a change in the graph topology, which only plays a role
in the crossover regime. All results support the view that a first-order DPT
may also appear in random walks on infinite-size random graphs.

##  I Introduction

Random walks on graphs are versatile tools to model real-world noisy dynamical
processes embedded in spatial structures Hughes1995 ; Noh2004 ; Barrat2008 ;
Newman2010 ; Latora2017 ; Masuda2017 . These processes describe both natural
and man-made phenomena such as the spreading of infectious diseases Barrat2008
; Pastor-Satorras2015 , the transport of vesicles in cell cytoskeletons
Julicher1997 , the propagation of information in communication networks
Castellano2009 ; Liu2014 , and the robustness of networks to random failures
Barrat2008  to name just a few examples. Often, the focus in these
applications is towards time-averaged quantities including stationary
distributions, and energy and particle currents. Indeed, these are observables
commonly used in applications to gather information on the average state
occupation and mobility in network structures Barrat2008 ; Masuda2017 . On the
other hand, fluctuations are also fundamental to understand the behavior of
physical systems living in unstable environments as rare events are often
responsible for the evolution dynamics Albeverio2006 ; Kishore2011 . However,
much less is known about them and in the last decades many research efforts
have been deployed towards the development of theoretical frameworks that
allow for their study, e.g., large deviation theory DenHollander2000 ;
Touchette2009 ; Dembo2010 ; Chetrite2015 ; Jack2020a ; Carugno2022 .

Recently, signatures of a dynamical phase transition (DPT), viz. a transition
between different fluctuation mechanisms, has been identified in the study of
the mean degree (connectivity) visited by unbiased random walks evolving on
sparse random graphs DeBacco2016 ; Coghi2019 ; Gutierrez2021 . There are good
grounds to consider it as a first-order DPT where we observe the coexistence
of two ‘phases’ characterized by random walk paths that visit the whole graph,
and paths localized in dangling chains, i.e., lowly connected structures of
the graph. However, a rigorous proof for ensembles of random graphs is still
lacking and, in fact, the community still debates on the real nature and
interpretation of DPTs Whitelam2018 ; Whitelam2021 .

In this paper, we contribute to the debate by analyzing two exactly solvable
models where the transition appears to be first-order and characterized by an
absorbing dynamics. This sees, on the one hand, the random walk fully
localized in dangling structures, and on the other hand, the random walk fully
absorbed by the bulk of the graph, which acts as an entropic basin and allows
the random walk to be fully delocalized. We make use of a theoretical
framework for the calculation of large deviations that we developed in
Carugno2022  and that allows us to: (i) consider general time-additive
observables, (ii) analytically characterize the behavior of random walks on
finite-size graphs, and (iii) rigorously study the scaling (with respect to
the size of the graph) of fluctuations around the critical value of the DPT.
Remarkably, in agreement with Whitelam2018 ; Whitelam2021  we notice that an
important ingredient for the appearance of a first order DPT in both models is
the presence of absorbing dynamics, generated by different scalings of the
hopping probabilities in the graph. Furthermore, we notice that although the
first order DPT appears in both the models we investigated, the scaling of the
fluctuations around the transition is different and we argue that it is both
function of the dynamical process and of the inherent topology of the network.

A brief outline of the paper follows. In Section II we set up a general model
of an URW hopping on a graph, discuss the general form of observables that we
consider in this manuscript, and introduce the theory of large deviations in
this setting. In Section III we collect our results related to two exactly
solvable models. In Section IV we conclude the paper by summarizing the
results obtained and briefly discussing open questions.

##  II Setting and large deviations

We consider an unbiased discrete-time random walk (URW)
$X=\left(X_{\ell}\right)_{\ell=1}^{n}=(X_{1},X_{2},\dots,X_{n})$ evolving on a
finite connected graph $G=(V,E)$, with $V$ denoting the set of $N$ vertices
(or nodes) and $E$ the set of edges (or links). The topology of the graph is
encoded in the symmetric adjacency matrix $A$, which has components

$A_{ij}=\begin{cases}1&i\in\partial j\\\ 0&\text{otherwise}\,,\end{cases}$ (1)

where $\partial j$ denotes the set of neighbors of node $j$. Notice that we
choose to consider an unweighted symmetric graph for simplicity, but our
methods can be easily generalized to more structured cases. The dynamics of
the random walk is defined by the transition matrix $\Pi$ having components

$\Pi_{ij}=\frac{A_{ij}}{k_{i}}\,,$ (2)

where $k_{i}=\sum_{j\in V}A_{ij}$ is the degree of node $i$, viz. the number
of edges in which node $i$ participates. The matrix $\Pi$ characterizes the
uniform probability of going from a vertex $X_{\ell}=i$ at time $\ell$ to a
vertex $X_{\ell+1}=j$ at time $\ell+1$ – that is, the probability of
transitioning from $i$ to $j\in\partial i$ does not depend on $j$.
Furthermore, for simplicity we restrict the random walk to be ergodic, viz.
$\Pi$ is irreducible and aperiodic. In the rest of the manuscript we use the
index $\ell$ to refer to time and the indices $i$ and $j$ to refer to nodes of
the graph.

The long-time behavior of the URW is well understood. Thanks to ergodicity,
the random walk has a unique stationary distribution

$\rho_{i}=\frac{k_{i}}{\sum_{j\in V}k_{j}}\,,$ (3)

which is found to be proportional to the degree of each node. Furthermore, the
URW is also reversible, viz. it is an equilibrium process, as it satisfies the
detailed balance condition

$\rho_{i}\Pi_{ij}=\rho_{j}\Pi_{ji}\,$ (4)

for each pair of nodes in $V$.

In this setting, we assume that the URW $X$ accumulates a cost in time given
by

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}f(X_{\ell})\,,$ (5)

where $f$ is any function of the vertex state. In nonequilibrium statistical
mechanics, this cost is also called a dynamical observable Touchette2009  and,
depending on $f$, it may represent interesting physical quantities, such as
occupation times Chetrite2015 , internal energy Sekimoto2010 , chemical
concentrations Dykman1998 , activities Gutierrez2021 , and entropy production
rates Coghi2019 . Because of the ergodicity of the URW, in the long-time limit
the observable $C_{n}$ converges with probability $1$ to the ergodic average

$\sum_{i\in V}\rho_{i}f(i)\eqqcolon c^{*}\,.$ (6)

This convergence property is often used to estimate properties of large graphs
such as degree distributions or centrality measures, by running random walks
(or, generally speaking, agents) on the graphs for long times Newman2010 .

Following the introduction, here we study fluctuations of $C_{n}$ around the
typical value $c^{*}$ by calculating its probability distribution
$\mathbb{P}(C_{n}=c)$ in the $n\rightarrow\infty$ limit. The probabilistic
theory of large deviations DenHollander2000 ; Touchette2009 ; Dembo2010  tells
us that this distribution has an exponentially decaying form

$\mathbb{P}(C_{n}=c)=e^{-nI(c)+o(n)}\,,$ (7)

described by the rate function $I$ given by the following limit

$I(c)=-\lim_{n\rightarrow\infty}\frac{1}{n}\log\mathbb{P}(C_{n}=c)\,.$ (8)

The rate function $I$ is a pivotal object in the theory of large deviations as
it characterizes the fluctuations of $C_{n}$ to leading order in $n$; it is a
non-negative function and it is equal to $0$ for ergodic random walks only at
$c^{*}$ (where the probability concentrates exponentially fast with time).

Much effort is drawn towards the development of methods that allow one to
calculate $I$ in (8) efficiently Touchette2009 . Spectral and variational
techniques can both be implemented and, depending on the particular model
studied, it may well be that some work better than others Coghi2021PhD .
Spectral techniques based on moment generating functions have the merit to
reformulate the problem in a different setting—similarly to a microcanonical-
canonical change of ensemble—whereas variational techniques based on the
contraction principle DenHollander2000 ; Touchette2009  are useful to find
probabilistic bounds Hoppenau2016 . In the following, we will base our large
deviation study on the techniques discussed in Carugno2022  which try to merge
the pros of spectral and variational methods.

In line with Carugno2022  and previous works Whittle1955 ; Dawson1957 ;
Goodman1958 ; Billingsley1961 ; CsiszaR1987 ; DenHollander2000 ; Polettini2015
; Dembo2010 , in order to calculate the rate function $I$ associated with the
observable $C_{n}$ in (5) we move the focus on to the study of the higher-
dimensional pair-empirical occupation measure

$L^{(2)}_{n}(i,j)=\frac{1}{n}\sum_{\ell=1}^{n}\delta_{X_{\ell},i}\delta_{X_{\ell+1},j}=\nu_{ij}\hskip
28.45274pt\forall i,j\in V\,,$ (9)

which counts the fraction of jumps $\nu_{ij}$ that the URW makes between each
couple of nodes in the graph – see Carugno2022 . Remarkably, the value of
$C_{n}$ can be deduced via the formula

$C_{n}=\sum_{i,j\in V}f(i)L_{n}^{(2)}(i,j)\,.$ (10)

We can calculate the rate function (8) by means of the Gärtner–Ellis theorem
Touchette2009 ; Dembo2010 . To do so we need to introduce the scaled cumulant
generating function (SCGF) of $C_{n}$, which is defined as

$\lambda_{s,N}[\nu^{*}]=\lim_{n\rightarrow\infty}\frac{1}{n}\mathbb{E}\left[e^{nsC_{n}}\right]\
,$ (11)

and calculate its Legendre–Fenchel transform, i.e.,

$I(c)=\sup_{s\in\mathbb{R}}\\{sc-\lambda_{s,N}[\nu^{*}]\\}\,.$ (12)

In Carugno2022 , we showed that $\lambda_{s,N}[\nu^{*}]$ can be obtained
minimizing the following action

$\displaystyle\lambda_{s,N}$
$\displaystyle[\nu]=\lambda_{1,N}[\nu]+\lambda_{2,N}[\nu]+\lambda_{3,N}[\nu]+\lambda_{4,N}[\nu]$
(13) $\displaystyle\lambda_{1,N}$
$\displaystyle[\nu]=\sum_{i=1}^{N}\sum_{j=1}^{N}\nu_{ij}\left(\log\left(\sum_{k=1}^{N}\nu_{ik}\right)-\log(\nu_{ij})\right)$
(14) $\displaystyle\lambda_{2,N}$
$\displaystyle[\nu]=\sum_{i=1}^{N}\sum_{j=1}^{N}\log(\Pi_{ij})\ \nu_{ij}$ (15)
$\displaystyle\lambda_{3,N}$
$\displaystyle[\nu]=s\sum_{i=1}^{N}f(i)\sum_{j=1}^{N}\ \nu_{ij}$ (16)
$\displaystyle\lambda_{4,N}$
$\displaystyle[\nu]=\epsilon\left(\sum_{i=1}^{N}\sum_{j=1}^{N}\nu_{ij}-1\right)+\sum_{i=1}^{N}\eta_{i}\left(\sum_{j=1}^{N}\nu_{ij}-\sum_{j=1}^{N}\nu_{ji}\right)\
,$ (17)

with respect to $\nu_{ij}$, $\epsilon$ and $\eta_{i}$, which are respectively
the fraction of jumps from node $i$ to node $j$, and the Lagrange multipliers
fixing the normalization constraint and the global balance. We remark that
these formulas are valid for any finite $N$. In our setting the rate function
$I$ in (12) reduces to

$I(c)=-\lambda_{1,N}[\nu^{*}]-\lambda_{2,N}[\nu^{*}]-\lambda_{4,N}[\nu^{*}]\,,$
(18)

where the dependence on the fluctuation $c$ enters through the minimizer
$\nu^{*}$, which depends on the optimized tilting parameter $s^{*}(c)$, i.e.,
$\nu^{*}\equiv\nu^{*}(s^{*}(c))$. For further details on the methods we used
to derive (13)-(17) and on a useful physical characterization of the action
(13) we refer the reader to Carugno2022  and related bibliography.

Although the equations for the minimum of (13) may be complicated to solve
analytically for complex models, the proposed approach has several advantages
that will also be highlighted in the next sections when studying simplified
scenarios. Firstly, the explicit form of the action (13) for any finite $N$
allows us to study directly the scaling of the fluctuations with varying graph
size. As we will see in the following, this is an important feature that will
help us in characterizing fluctuations around critical points. Furthermore,
the tilting parameter $s$ is responsible for biasing the dynamics of the URW
via (16) to realize a fluctuation $c$ for the observable $C_{n}$ fixed by the
Legendre duality equation

$c=\frac{d\lambda_{s,N}[\nu]}{ds}\,.$ (19)

Therefore, the minimizer $\nu^{*}$ of the action $\lambda_{s,N}$ characterizes
the typical configuration of jumps on the graph $G$ that give rise to the
fluctuation $c$ defined by (19) (similarly to Gutierrez2021 ).
It is natural to introduce a biased dynamics for which $C_{n}=c$ is realized
in the typical state: this biased dynamics has been thoroughly characterized
in its general form in Chetrite2015 ; Chetrite2015a  and also in the setting
of random walks on graphs in Coghi2019 . The process characterized by the
biased dynamics is known as driven (or effective/auxiliary) process and in
this context is a locally-biased version of the URW whose transition
probability matrix has been described in Coghi2019 . Within our approach, we
can also fully define the driven process and its transition matrix, which
reads

$\left(\Pi_{s}\right)_{ij}=\frac{\nu^{*}_{ij}(s)}{\sum_{k\in
V}\nu^{*}_{ik}(s)}\,.$ (20)

In other words, the driven process is the effective biased random walk that
explains how a fluctuation $C_{n}=c$ is created up to time $n$.

We remark that the method used here to calculate the SCGF
$\lambda_{s,N}[\nu^{*}]$ via the minimization of (13) is equivalent to
spectral methods based on the so-called tilted matrix Touchette2009 ;
Dembo2010 ; Touchette2018 ; Coghi2021PhD . In particular, in Carugno2022  we
show that the Euler–Lagrange equations for the minima of (13) are a useful re-
writing of the dominant eigenvalue problem associated with the tilted matrix.
The main pros of the method reviewed in Carugno2022  are (i) to give a clear
physical interpretation of all the terms in the action and SCGF and (ii) to
express the driven process in terms of the minimizers of (13), which are the
optimal jumps that create a fluctuation $C_{n}=c$.

##  III Delocalization-localization dynamical phase transition

Recently, it has been pointed out DeBacco2016 ; Coghi2019 ; Gutierrez2021
that an URW that accumulates a cost proportional to the degree of each visited
node, e.g., $f(X_{l})=k_{X_{l}}$ in (5), and that runs on the largest
connected component of an Erdős–Rényi random graph seems to undergo a DPT.
This transition is localized in the fluctuations of the mean degree visited
when this is lower than the mean connectivity of the graph and is interpreted
as a ‘co-existence’ of paths that visit nodes with low degree and paths that
visit the whole graph DeBacco2016 ; Coghi2019 . [Noticeably, another DPT may
arise when the random walk visits more often highly connected regions of the
graph and localizes around the highest degree node of the graph. Although as
interesting, in this manuscript we will not focus on this behavior.] The
Erdős–Rényi random graph in question is picked from a ‘canonical’ ensemble of
random graphs having a fixed number of nodes $N$ and edges randomly placed
with a small probability $p$ between each pair of nodes such that $Np=\bar{c}$
is fixed to be the mean degree of the graph. Noticeably, signatures of the DPT
disappear when $\bar{c}$ is large, revealing that a fundamental topological
ingredient for the appearance of the transition is the presence of lowly
connected structures in the graph (such as trees and dangling chains of
nodes). These structures carry strong spatial correlations—a node of degree
one is likely to be connected with a node of degree two in a dangling chain
and these correlations are responsible for the dynamics of the random walk
when visiting low-degree nodes. The overall picture is that of a random walk
whose behavior fluctuates between two distinct phases characterized by (i)
being localized in lowly-connected regions of the graph and (ii) being spread
over the bulk (most connected region) of the graph which acts as an ergodic
basin absorbing the dynamics (see model in Appendix A of Coghi2019  and also
Whitelam2018 ).

As far as the current state of the art is concerned, it is not clear whether
such a DPT appears in the infinite-size limit of ensemble of random graphs.
However, as mentioned in the previous paragraph, various numerical studies
indicate an abrupt change in the mechanisms that generate fluctuations,
endorsing the idea of a dynamical phase transition DeBacco2016 ; Coghi2019 ;
Gutierrez2021 .

In the following, by applying the theory discussed in Section II, we
analytically characterize the DPT in two models, which catch what we think are
the most relevant physical features of this phenomenon. We believe that these
characteristics are shared by the dynamics of URWs on Erdős–Rényi random
graphs: in particular, the heterogeneity of the scaling of the degree and the
presence of lowly connected regions such as dangling chains. We show that (i)
the DPT is first order – that is, $\lambda_{s,N}[\nu^{*}]$ as a function of
$s$ has a non-differentiable point $s_{c}$ in which the first derivative is
discontinuous; (ii) the behavior around $s_{c}$ is characterized by a scaling
function which is not universal, and depends on both the dynamics of the URW
and the topology of the graph; (iii) the driven process can be fully
characterized, allowing us to understand how fluctuations arise in time. These
results give further evidence of the presence of a first-order DPT in random
walks exploring random graphs.

###  III.1 Bulk-dangling model

The first model we look at is based on an URW with transition matrix (2)
collecting a cost (5) of the form

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}\frac{k_{X_{\ell}}}{\bar{k}}\,,$ (21)

by visiting a graph of $N$ nodes composed by a fully connected bulk of $N-2$
nodes and a single dangling chain of $2$ nodes (see Fig. 1). The structure of
this graph incorporates two relevant features: (i) the presence of a spatially
correlated dangling chain and (ii) a fully connected bulk that allows the URW
to uniformly spread over the network. Given this structure, the mean degree of
the graph is $\bar{k}=((N-3)(N-3)+(N-2)+2+1)/N$, which evidently scales
linearly with $N$ for large-size graphs. This feature allows us to deduce the
behavior of the observable $C_{n}$ in two opposite scenarios in the
$N\rightarrow\infty$ limit. Evidently, if the random walk is uniformly spread
over the bulk $C_{n}\sim 1$, whereas if the random walk is localized in the
chain $C_{n}\sim 0$. We argue that this behavior does not depend on the length
$L$ of the dangling chain – the choice $L=2$ is made to simplify calculations.

Figure 1: Bulk-dangling model graph for $N=11$. Node $1$ has degree $k_{1}=1$,
node $2$ has degree $k_{2}=2$ and node $3$ has degree $k_{3}=N-2=9$: the first
two nodes represent the dangling chain, while node $3$ bridges the chain with
the bulk, being part of the latter. All the remaining $N-3=8$ nodes are of
type $4$, having degree $k_{4}=N-3=8$. Together with node $3$, they form the
fully connected bulk.

Following Section II we calculate the action $\lambda_{s,N}[\nu]$ in (13).
Because of the bulk of the graph being highly symmetric, i.e., every link in
the bulk is equivalent, and the global balance imposed on the dynamics, i.e.,
incoming and outgoing flow of a node being equal, we are only left with four
degrees of freedom (variables) $\nu_{ij}$ that determine the action. We name
the fraction of jumps

$\displaystyle\nu_{12}$ for both directions: $1\rightarrow 2$ and
$2\rightarrow 1$ (22) $\displaystyle\nu_{23}$ for both directions:
$2\rightarrow 3$ and $3\rightarrow 2$ (23) $\displaystyle\nu_{34}$ for both
directions of each link: $3\rightarrow 4$ and $4\rightarrow 3$ (24)
$\displaystyle\nu_{44}$ $\displaystyle\;\;\text{for both directions of each
link in the bulk}\,.$ (25)

Notice that if at this stage we also imposed the normalization constraint,
i.e., $\sum_{i,j\in V}\nu_{ij}=1$, we would be left with only three degrees of
freedom. However, for practical reasons in the calculation of the minimum of
the action $\lambda_{s,N}[\nu]$, we leave this last constraint as an implicit
parametrization with a Lagrange multiplier $\epsilon$ entering the action.

The action can be explicitly written as

$\lambda_{s,N}[\nu]=h(s,N,\nu_{12},\nu_{23},\nu_{34},\nu_{44})+\epsilon(1-2\nu_{12}-2\nu_{23}-2(N-3)\nu_{34}-(N-3)(N-4)\nu_{44})\,,$
(26)

with the long form of the function $h$ postponed to the Appendix A. The
minimum and minimizers of the action (26) can be found by solving the saddle-
point equations that can be cast in the following linear system:

$\displaystyle\nu_{12}=\nu_{23}\frac{a(s,N,\epsilon)}{1-a(s,N,\epsilon)}$ (27)
$\displaystyle\nu_{23}=\nu_{34}\frac{b(s,N,\epsilon)}{1-b(s,N,\epsilon)}$ (28)
$\displaystyle\nu_{34}=\nu_{23}\frac{c(s,N,\epsilon)}{1-(N-3)c(s,N,\epsilon)}$
(29)
$\displaystyle\nu_{44}=\nu_{34}\frac{d(s,N,\epsilon)}{1-(N-4)d(s,N,\epsilon)}$
(30) $\displaystyle
2\nu_{12}+2\nu_{23}+2(N-3)\nu_{34}+(N-3)(N-4)\nu_{44}=1\,,$ (31)

with

$\displaystyle a(s,N,\epsilon)=e^{3\frac{s}{\bar{k}}-\log 2+2\epsilon}$ (32)
$\displaystyle b(s,N,\epsilon)=\frac{e^{N\frac{s}{\bar{k}}-\log
2-\log(N-2)+2\epsilon}}{1-e^{3\frac{s}{\bar{k}}-\log 2+2\epsilon}}$ (33)
$\displaystyle
c(s,N,\epsilon)=\frac{e^{(2N-5)\frac{s}{\bar{k}}-\log(N-2)-\log(N-3)+2\epsilon}}{1-e^{(N-3)\frac{s}{\bar{k}}-\log(N-3)+\epsilon}}$
(34) $\displaystyle
d(s,N,\epsilon)=e^{(N-3)\frac{s}{\bar{k}}-\log(N-3)+\epsilon}\,.$ (35)

We report the form of the minimizers
$\nu^{*}=(\nu^{*}_{12},\nu^{*}_{23},\nu^{*}_{34},\nu^{*}_{44})$ in the
Appendix A. It is important to remark that these minimizers are not yet in a
fully explicit form as they depend on the Lagrange parameter $\epsilon$
(which, in Carugno2022 , has also been proved to be the negative SCGF
$\lambda_{s,N}[\nu^{*}]$). However, the Lagrange parameter $\epsilon$—hence,
also the SCGF we are after—can be determined by solving the normalization
constraint in (31) after having replaced the form of the minimizers
$\nu^{*}=(\nu^{*}_{12},\nu^{*}_{23},\nu^{*}_{34},\nu^{*}_{44})$. The equation
reads

$\begin{split}2(N-2)(N-3)-2&\tau(N-2)(N-4)e^{(N-3)\frac{s}{\bar{k}}}-\tau^{2}(N-3)\left((N-2)e^{3\frac{s}{\bar{k}}}+e^{N\frac{s}{\bar{k}}}+2e^{(2N-5)\frac{s}{\bar{k}}}\right)+\\\
&\tau^{3}(N-4)\left((N-2)e^{N\frac{s}{\bar{k}}}+e^{(2N-3)\frac{s}{\bar{k}}}\right)+\tau^{4}(N-3)e^{(2N-2)\frac{s}{\bar{k}}}=0\,,\end{split}$
(36)

with $\tau=e^{\epsilon}$. Noticeably, (36) can also be derived by imposing
that the matrix of the coefficients of the four linear equations (27)–(30) has
a nullspace of dimension greater than zero—that is, when its determinant is
zero. Equation (36) is fourth order in $\tau$ and hence it admits four
solutions of which only one is physical. This can be selected by noticing
that, since $\nu^{*}$ must be positive, the right hand side of the four
equations (27)–(30) is also positive (we postpone the exact form of the
inequality constraints to the Appendix A). In this way, we obtained the SCGF
$\lambda_{s,N}[\nu^{*}]$ valid for any finite-size graph.

(a) SCGF $\lambda_{s,N}$

(b) Derivative of SCGF $\lambda^{\prime}_{s,N}$

(c) Rate function $I(c)$

Figure 2: Large deviation study for the bulk-dangling model. In all three
figures, different colors correspond to a different number $N$ of nodes: i)
light blue is $N=15$; ii) orange is $N=25$; iii) green is $N=100$; iv) black
is $N\to\infty$. All finite $N$ curves where obtained by solving (36)
numerically, while analytical expressions for the black curves are presented
in (37) for figure (a), (38) for figure (b) and (39) for figure (c).

By carefully taking the limit $N\rightarrow\infty$ in the polynomial equation
(36) we can also explicitly obtain the SCGF in the infinite-size limit of the
graph, which is

$\lambda_{s,\infty}=\begin{cases}-\frac{\log 2}{2}&s\leq-\frac{\log 2}{2}\\\
s&s>-\frac{\log 2}{2}\end{cases}$ (37)

and highlights a non-differentiable point at $s_{c}=-\log 2/2$. The derivative
of $\lambda_{s,\infty}$, according to (19), is

$\frac{d\lambda_{s,\infty}}{ds}=\begin{cases}0&s\leq-\frac{\log 2}{2}\\\
1&s>-\frac{\log 2}{2}\,,\end{cases}$ (38)

which explicitly describes the fluctuation $C_{n}=c$ happening with varying
tilting parameter $s$. This confirms our expectations: on the left of the
critical point $s_{c}$ the random walk is localized in the dangling chain—the
only region of the graph where the cost accumulated $C_{n}$ (see (5)) does not
scale with the size $N$—whereas on the right of $s_{c}$ the random walk is
spread in the bulk where it accumulates a cost that scales linearly with $N$.
Furthermore, the value of $s_{c}$—and with it all the left branch of
$\lambda_{s,\infty}$ in (37)—can easily be interpreted as the mean entropy of
the random walk that, localized in the dangling chain, keeps going back and
forth from the node of degree $1$ to the node of degree $2$ (see also
Carugno2022 ). Eventually, the rate function $I$ can easily be obtained by
Legendre transforming the two analytical branches of (37) and by connecting
them with a linear section or by implementing directly (8); in both cases we
obtain

$I(c)=\begin{cases}\frac{\log 2}{2}-c\frac{\log 2}{2}&0\leq c\leq 1\\\
\infty&\text{otherwise}\,,\end{cases}$ (39)

and we remark that the non-differentiable point $s_{c}$ for the SCGF
$\lambda_{s,\infty}$ is mapped onto the linear section characterizing the rate
function. We graphically show in Fig. 2 the SCGF, its derivative, and the rate
function for finite-size graphs and in the infinite-size limit.

The non-differentiability of the SCGF can be physically related to a first-
order DPT that is interpreted here as a coexistence of paths that either visit
predominantly the bulk of the graph ($C_{n}\sim 1$) or are localized in the
dangling chain ($C_{n}\sim 0$). A further characterization of this DPT is
given by identifying the mechanisms that give rise to the fluctuations around
the critical point $s_{c}$.

As it appears from formula (38) and Fig. 2(b), the normalized mean-degree
visited by the URW (21) computed from (19) is a piece-wise constant function
of the tilting parameter $s$ in the infinite-size limit of the graph. We refer
to the region $s<s_{c}$ ($s>s_{c}$) corresponding to the localized
(delocalized) phase as $s^{-}$ ($s^{+}$) and we study in these two regions the
scaling with $N$ of the transition probabilities of the driven process (20).
The calculation can be done by properly taking the $N\rightarrow\infty$ limit
of the minimizer
$\nu^{*}=(\nu^{*}_{12},\nu^{*}_{23},\nu^{*}_{34},\nu^{*}_{44})$ and inserting
the result in (20). We get the following two transition matrices that
characterize the probability of stepping from a node to another one in the
graph of Fig. 1:

$\Pi_{s^{-}}=\left(\begin{array}[]{cccccc}0&1&0&\cdots&0&\cdots\\\
1+O(N^{-1})&0&O(N^{-1})&\cdots&0&\cdots\\\
0&-\frac{e^{s}}{3s}+O(N^{-1})&0&O(N^{-1})&\cdots&\cdots\\\
\vdots&&\ddots&&\ddots&\cdots\\\
0&0&1-\sqrt{2}e^{s}&0&O(N^{-1})&\cdots\end{array}\right)$ (40)
$\Pi_{s^{+}}=\left(\begin{array}[]{cccccc}0&1&0&\cdots&0&\cdots\\\
\frac{e^{-2s}}{2}+O(N^{-1})&0&\left(1-\frac{e^{-2s}}{2}\right)+O(N^{-1})&\cdots&0&\cdots\\\
0&O(N^{-1})&0&O(N^{-1})&\cdots&\cdots\\\ \vdots&&\ddots&&\ddots&\cdots\\\
0&0&O(N^{-1})&0&O(N^{-1})&\cdots\end{array}\right)\,.$ (41)

Evidently, for fluctuations obtained by fixing $s<s_{c}$ the random walk is
biased towards localizing in the dangling chain, e.g., if the random walk is
on node two, the probability of hopping onto node one is one order of
magnitude (with respect to the system size) bigger than moving towards the
fully connected bulk. For $s>s_{c}$ instead, the bulk behaves as an entropic
basin absorbing the random walk and allowing it to be fully spread over the
graph.

Figure 3: Crossover regime of the SCGF of the bulk-dangling model around
$s_{c}$ as a function of the scaling variable $\tilde{s}$ for different values
of $N$. As $N$ increases, the colored curves collapse into the limiting curve
predicted theoretically (45).

We conclude the study of the bulk-dangling model showing how fluctuations
scale with the graph size locally around the critical point $s_{c}$ (in
analogy with the study carried out in Whitelam2018 ). This can be done by
centering and rescaling the tilting variable $s$ as

$s=-\frac{\log 2}{2}+\frac{\tilde{s}}{N}\,,$ (42)

and the Lagrange parameter (or negative SCGF)

$\epsilon=\frac{\log 2}{2}-\frac{\tilde{\epsilon}_{\tilde{s}}}{N}\,,$ (43)

in the polynomial equation (36). Using this scaling and expanding the
polynomial to leading order in $N$ we obtain

$\lambda_{\tilde{s},N}\approx-\frac{\log
2}{2}+\frac{\tilde{\epsilon}_{\tilde{s}}}{N}\,,$ (44)

with

$\tilde{\epsilon}_{\tilde{s}}=\frac{1}{8}\left(\sqrt{2}+4\tilde{s}-7\log
2+\sqrt{2+16\sqrt{2}+2\sqrt{2}\log 2-\log^{2}2-8\sqrt{2}\tilde{s}-8\log
2\tilde{s}+16\tilde{s}^{2}}\right)\,,$ (45)

which explains how the SCGF locally scales as a function of the graph size
around $s_{c}$. We report in Fig. 6 the function
$\tilde{\epsilon}_{\tilde{s}}$ (translated to be centered in $(0,0)$ and not
in $(-\log 2/2,-\log 2/2$). Evidently, the function continuously joins the two
branches of fluctuations separated by the critical point $s_{c}$ in Fig. 2(a):
on the left, for $\tilde{s}\ll 0$, $\tilde{\epsilon}_{\tilde{s}}$ tends to $0$
(hence, $-\log 2/2$)), on the right, for $\tilde{s}\gg 0$, it behaves linearly
with respect to $\tilde{s}$.

In conclusion, the critical point $s_{c}$ marks a first-order DPT for the
observable $C_{n}$ in (21), however, thanks to a proper rescaling showed in
(42) and (43) we can get more precise information on how fluctuations scale
with the system size around the critical point $s_{c}$.

###  III.2 Two-state Markov chain

In this Subsection we analyze another model which takes its cue from the
findings in the previous model and is also inspired by the works of
Whitelam2018 ; Coghi2019 . The model is a two-state Markov chain as
represented in Fig. 4. If the Markov chain is found on the state on the left,
namely $1$, at a certain time $\ell$, it collects a unitary reward $1$,
whereas if it is on the right, namely $b$, it collects a reward $b\geq 1$
(eventually $b\rightarrow\infty$). Further, although the probability of moving
from the left to the right is totally unbiased, the probability of moving from
the right to the left inversely scales with the reward $b$.

Figure 4: Sketch of the two-state Markov chain model, composed by state $1$
and state $b$. The transition probabilities – depicted in the figure above the
arrows – read explicitly: $p_{11}=1/2$, $p_{1b}=1/2$, $p_{b1}=1/b$,
$p_{bb}=1-1/b$.

The observable we focus on has the general form in (5) and in this particular
scenario it reduces to

$C_{n}=\frac{1}{n}\sum_{\ell=1}^{n}\frac{X_{\ell}}{b}\,,$ (46)

which is the mean reward collected over time renormalized by the maximum
reward $b$.

This model tries to catch once again the most relevant physical ingredients
that may lead to a delocalization-localization first-order DPT. In doing so,
however, we take a further simplification: we try to rule out as much as we
can the graph topology, replacing bulk and dangling contributions with two
single states which respectively give a reward of $b$ and $1$ to the observed
cost in (46). In comparison with the previous model, the reward $b$ should be
analogous to the graph size $N$—a random walk lost in the bulk of a graph
observes nodes with a degree that scales with $N$ in (21)—whereas the reward
$1$ should mimic the observed degree in the dangling chain. Furthermore, the
inverse scaling with $b$ of the probability of moving from the right state to
the left one, should give rise to an absorbing dynamics in the right state for
$b\rightarrow\infty$. As we will see in the following, the topology of the
graph does not play a pivotal role in the appearance of the first-order DPT,
but it may play an important role in determining the exact fluctuations in the
crossover regime around the critical point.

Once again, following Section II we calculate the action $\lambda_{s,b}[\nu]$
in (13). Because the global balance is imposed on the dynamics, we only have
to deal with three degrees of freedom (variables) $\nu_{ij}$ that determine
the action. These are

$\displaystyle\nu_{11}$ fraction of jumps from $1$ to $1$ (47)
$\displaystyle\nu_{1b}$ for both directions: $1\rightarrow b$ and
$b\rightarrow 1$ (48) $\displaystyle\nu_{bb}$ $\displaystyle\;\;\text{fraction
of jumps from $b$ to $b$}\,.$ (49)

Notice that if we also imposed the normalization constraint, i.e.,
$\sum_{i,j\in\left\\{1,b\right\\}}\nu_{ij}=1$, we would be left with only two
degrees of freedom. However, analogously to the previous Subsection, we leave
this last constraint as an implicit parametrization with a Lagrange multiplier
$\epsilon$ entering the action.

The action can be explicitly written as

$\begin{split}\lambda_{s,b}[\nu]&=-\nu_{11}\log\nu_{11}-\nu_{bb}\log\nu_{bb}-2\nu_{1b}\log\nu_{1b}+(\nu_{1b}+\nu_{bb})\log(\nu_{1b}+\nu_{bb})+(\nu_{1b}+\nu_{11})\log(\nu_{1b}+\nu_{11})+\\\
&\hskip 28.45274pt+\frac{s}{b}(\nu_{11}+\nu_{1b})+s(\nu_{1b}+\nu_{bb})-\log
2\nu_{11}+\log\frac{(b-1)}{b}\nu_{bb}+\nu_{1b}(-\log 2-\log
b)+\epsilon(2\nu_{1b}+\nu_{11}+\nu_{bb}-1)\,.\end{split}$ (50)

The minimum and minimizers of the action (50) can be found by solving the
saddle-point equations and imposing the normalization constraint. We get

$\displaystyle\nu_{11}=\frac{e^{\frac{s}{b}+\epsilon}(-b+(b-1)e^{s+\epsilon})}{-4b+2(b-1)e^{s+\epsilon}+be^{\frac{s}{b}+\epsilon}}$
(51)
$\displaystyle\nu_{bb}=\frac{(b-1)e^{s+\epsilon}(-2+e^{\frac{s}{b}+\epsilon})}{-4b+2(b-1)e^{s+\epsilon}+be^{\frac{s}{b}+\epsilon}}$
(52)
$\displaystyle\nu_{1b}=\frac{1}{\frac{b}{b-(b-1)e^{s+\epsilon}}-\frac{2}{-2+e^{\frac{s}{b}+\epsilon}}}\,,$
(53)

as still functions of the Lagrange multiplier $\epsilon$. This last can be
determined by, for instance, using the equation for the minimum of
$\lambda_{s,b}[\nu]$ w.r.t $\nu_{1b}$ and by replacing the values of
$\nu_{11}$ and $\nu_{bb}$ with those in (51) and (52). We find that the SCGF
$\lambda_{s,b}[\nu^{*}]$ is analytically given by

$\lambda_{s,b}=-\epsilon=\log\left[\frac{1}{4b}\left(be^{\frac{s}{b}}+2(b-1)e^{s}+\sqrt{4(b-1)^{2}e^{2s}+b^{2}e^{\frac{2s}{b}}-4(b-3)be^{\frac{bs+s}{b}}}\right)\right]\,.$
(54)

By carefully taking the limit $b\rightarrow\infty$ of (54) we explicitly
obtain the SCGF in the infinite size limit of the reward, which reads

$\lambda_{s,\infty}=\begin{cases}-\log 2&s\leq-\log 2\\\ s&s>-\log
2\end{cases}$ (55)

and highlights the appearance of a non-differentiable point at $s_{c}=-\log 2$
(this value can always be interpreted as the mean entropy of the random walk
localized in $1$). The derivative of $\lambda_{s,\infty}$, as in (19),
explicitly describes the fluctuation $C_{n}=c$ happening with varying tilting
parameter $s$ and analogously to (37) we obtain

$\frac{d\lambda_{s,\infty}}{ds}=\begin{cases}0&s\leq-\log 2\\\ 1&s>-\log
2\,.\end{cases}$ (56)

This says that on the left of the critical point $s_{c}$ the Markov chain is
localized in the left state where it accumulates a cost that does not scale
with the reward $b$, whereas on the right of $s_{c}$ the Markov chain is
localized in the right state where it accumulates a cost $b$ at every step. We
remark that at finite $N$ the critical point is absent, replaced by a
crossover region where the Markov chain visits both nodes for a finite
fraction of time.

The rate function $I$ can also be easily obtained as explained in the previous
Subsection and reads

$I(c)=\begin{cases}\log 2-c\log 2&0\leq c\leq 1\\\
\infty&\text{otherwise}\,.\end{cases}$ (57)

We graphically show in Fig. 5 the SCGF, its derivative, and the rate function
for the finite reward case and in the infinite-reward $b$ limit. Noticeably,
the rate function obtained in (39) is exactly half the rate function obtained
above here. This is consequence of the mean entropy
$\lambda_{1,N}+\lambda_{2,N}$ (14), (15) that the random walk has in the
localized state: in the bulk-dangling model is half with respect to the two-
state model presented here.

(a) SCGF $\lambda_{s,b}$

(b) Derivative of SCGF $\lambda^{\prime}_{s,b}$

(c) Rate function $I(c)$

Figure 5: Large deviation study for the two-state model. In all three figures,
different colors correspond to a different value $b$ of the reward: i) light
blue is $b=15$; ii) orange is $b=50$; iii) green is $b=250$; iv) black is
$b\to\infty$. All finite $b$ curves where obtained from (54), while analytical
expressions for the black curves are presented in (55) for figure (a), (56)
for figure (b) and (57) for figure (c).

The non-differentiability of the SCGF can be physically related to a first-
order DPT also in this case. Once again, this is interpreted as a coexistence
of paths that are either absorbed by the state $b$ ($C_{n}\sim 1$) or are
localized in the state $1$ ($C_{n}\sim 0$). We can further characterize this
DPT by writing the driven process (20) that leads to fluctuations for
$s^{-}\equiv s<s_{c}$ or for $s^{+}\equiv s>s_{c}$. This can be done by
properly taking the $b\rightarrow\infty$ limit of the minimizer
$\nu^{*}=(\nu^{*}_{11},\nu^{*}_{1b},\nu^{*}_{bb})$ and inserting the results
in (20). We get the following two transition matrices:

$\Pi_{s^{-}}=\left(\begin{array}[]{cc}1+O(b^{-1})&O(b^{-1})\\\
1-2e^{s}+O(b^{-1})&2e^{s}+O(b^{-1})\end{array}\right)$ (58)

$\Pi_{s^{+}}=\left(\begin{array}[]{cc}\frac{1}{2e^{s}}+O(b^{-1})&1-\frac{1}{2e^{s}}+O(b^{-1})\\\
O(b^{-1})&1+O(b^{-1})\end{array}\right)\,.$ (59)

For $s<s_{c}$ the Markov chain is biased towards localizing in the state $1$,
whereas for $s>s_{c}$, the state $b$ absorbs the Markov chain. This is very
similar to what we have seen in the bulk-dangling model, with the only
difference that now the role of the topology has been replaced by different
rewards on the two states of the chain.

Although this structural change in the model does not seem to affect the
appearance of a first-order DPT, we notice that fluctuations scale differently
around the critical point $s_{c}=-\log 2$. This is made evident by rescaling
the tilting parameter $s$ and the SCGF $\lambda_{s,b}$ similarly to the
previous Subsection, we obtain

$\lambda_{\tilde{s},b}\approx-\log
2+\frac{\tilde{\epsilon}_{\tilde{s}}}{2\sqrt{b}}\,,$ (60)

with

$\tilde{\epsilon}_{\tilde{s}}=\tilde{s}+\sqrt{4+\tilde{s}^{2}}\,.$ (61)

Eq. (60) describes fluctuations locally around $s_{c}$ for large (but finite)
reward $b$. Also in this case, we plot in Fig. 6 the function
$\tilde{\epsilon}_{\tilde{s}}$ along with $b$-finite scalings.

Figure 6: Crossover regime of the SCGF of the two-state Markov chain model
around $s_{c}$ as a function of the scaling variable $\tilde{s}$ for different
values of $b$. As $b$ increases, the colored curves collapse into the limiting
curve predicted theoretically (61).

The scaling of fluctuations is much different from what we found in (44) for
the bulk-dangling model. We argue that the exact form of the scaling is not
only determined by the dynamics of the model, but also by the topology of the
graph considered. Indeed, the two-state Markov chain is only composed by two
nodes playing the role of bulk and dangling chain, whereas in the bulk-
dangling model, as previously mentioned, we count four key nodes (see Fig. 2)
and among these, differently from the two-state Markov chain, the orange and
red one play the role of a gate between the bulk and the yellow node of degree
one.

To further corroborate this argument we also studied generalizations of the
two-state Markov chain analyzed so far. These are obtained by considering as a
probability to escape the rightmost state the value $b^{-\gamma}$, with
$\gamma\geq 1$, and by rescaling, or not, the reward $b$ by $b^{-\gamma}$. In
all cases considered (not shown here), the scaling of fluctuations around the
critical point for the DPT are different from the case of the bulk-dangling
model. These results support the argument that changing the dynamics does not
make up for having different graph topologies.

To investigate the robustness of the aforementioned DPT, we investigated also
other variants of the two-state model, where the reward of the two nodes are
set to $1$ and $k$ respectively. In this version of the model $k$ does not
depend on $b$, and $1/b$ is only the transition probability to remain in state
$k$. Interestingly, the DPT appears also in this case when $b$ goes to
infinity, but both the critical tilting parameter $s_{c}$ and the behavior for
$s<s_{c}$ depend on the value of $k$. This is in line with previous works on
two state models Whitelam2018 , although with a remarkable difference. In
these works, the author considers models where the transition matrix is
symmetric, so that both nodes become absorbing in the appropriate limit. As a
consequence, $s_{c}$ in those models is exactly $0$. In our work, instead, we
are naturally led to consider non-symmetric transition matrices, such that
only one of the two states becomes absorbing. To win this asymmetric absorbing
dynamics, an infinitesimal $s$ is not sufficient, hence $s_{c}\neq 0$.

##  IV Conclusions

In this manuscript, we have shown the appearance of a first-order dynamical
phase transition in two models that catch the relevant physical aspects of the
dynamics of random walks on random graphs, for which the dynamical phase
transition has hitherto only been argued. In both models, the random walk
collects a cost—with general form given in (5)—which scales differently in
different regions of the graph. In the bulk-dangling model, very much
similarly to a random walk on a random graph, the cost scales proportionally
to the size of the graph in the bulk, whereas it gives only a constant
contribution in the dangling chain. In the two-state Markov chain instead, we
greatly simplified the topology of the graph and made the cost scale with a
reward rather than keeping it linked to the graph structure. As a consequence,
to keep the analogy with random walks on random graphs, we also suitably
rescaled the transition probabilities inversely with the reward. We analyzed
both models by applying a large deviation framework Carugno2022  that allowed
us to carry out analytical results.

Remarkably, regardless of the precise details of the model, a first-order
dynamical phase transition in the cost accumulated by the random walk always
appears (see Fig. 2 and 5). This is interpreted as a coexistence of paths that
visit regions of the graph where the cost scales proportionally with the
relevant physical parameter of the model (size $N$ or reward $b$) and paths
that visit regions that only contribute to the cost with constant increments.
We gave further evidence for this interpretation by also calculating the
relevant driven process, which explains how fluctuations arise in time (see
Eqs. (40) and (41), and (58) and (59)). Furthermore, by zooming around the
critical value for the transition, we exactly determined how fluctuations
scale either with the system size $N$ or the reward $b$. Since the scaling
turns out to be different in the two models investigated, we argue that
although the dynamical phase transition is robust to topological changes in
the model in the thermodynamic limit, the exact structure of the graph plays a
role—along with the dynamics—for finite systems.

These results support the idea that also random walks on sparse random graphs
undergo first-order phase transitions in the fluctuations of the mean-degree
visited DeBacco2016 ; Coghi2019  for infinite-size graphs. However, a full
proof has yet to be advanced. We believe that by implementing the large
deviation framework discussed in Carugno2022  and in this manuscript one
should be able to average the relevant action over the infinite realizations
of the random graph ensemble, obtaining the scaled cumulant generating
function in the thermodynamic-size limit. Furthermore, it would also be
interesting to study transient regimes in time—and not only the asymptotics
given by the large deviation theory—of the URW fluctuations. This could be
done in principle by following ideas presented in Polettini2015 ; Causer2022 ;
Carugno2022 .

Evidently, there is much scope for future work, both theoretical, as just
mentioned, and applied. Related to the latter, for instance, by appropriately
tuning the tilting parameter one could exploit the driven processes to
generate optimal explorers of networks, a topic that has recently gained much
attention Adam2019 ; Carletti2020 .

## References

  * (1) B. D. Hughes, Random walks and random environments: Random walks.  Oxford University Press, 1995.
  * (2) J. D. Noh and H. Rieger, “Random walks on complex networks,” Physical Review Letters, vol. 92, no. 11, p. 118701, 2004.
  * (3) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical processes on complex networks.  Cambridge University Press, 2008.
  * (4) M. E. J. Newman, Networks: An Introduction.  Oxford University Press, 2010.
  * (5) V. Latora, V. Nicosia, and G. Russo, Complex Networks: principles, methods and applications.  Cambridge University Press, 2017.
  * (6) N. Masuda, M. A. Porter, and R. Lambiotte, “Random walks and diffusion on networks,” Physics Reports, vol. 716-717, pp. 1–58, 2017.
  * (7) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, “Epidemic processes in complex networks,” Reviews of Modern Physics, vol. 87, no. 3, p. 925, 2015.
  * (8) F. Jülicher, A. Ajdari, and J. Prost, “Modeling molecular motors,” Reviews of Modern Physics, vol. 69, no. 4, pp. 1269–1281, 1997.
  * (9) C. Castellano, S. Fortunato, and V. Loreto, “Statistical physics of social dynamics,” Reviews of Modern Physics, vol. 81, no. 2, pp. 591–646, 2009\.
  * (10) C. Liu and Z. K. Zhang, “Information spreading on dynamic social networks,” Communications in Nonlinear Science and Numerical Simulation, vol. 19, no. 4, pp. 896–904, 2014.
  * (11) S. Albeverio, V. Jentsch, and H. Kantz, eds., Extreme Events in Nature and Society.  The Frontiers Collection, Berlin, Heidelberg: Springer Berlin Heidelberg, 2006.
  * (12) V. Kishore, M. S. Santhanam, and R. E. Amritkar, “Extreme events on complex networks,” Physical Review Letters, vol. 106, no. 13, p. 188701, 2011\.
  * (13) F. den Hollander, Large Deviations.  American Mathematical Society, 2000.
  * (14) H. Touchette, “The large deviation approach to statistical mechanics,” Physics Reports, vol. 478, no. 1-3, pp. 1–69, 2009.
  * (15) A. Dembo and O. Zeitouni, Large Deviations Techniques and Applications, vol. 38 of Stochastic Modelling and Applied Probability.  Springer Berlin Heidelberg, 2010.
  * (16) R. Chetrite and H. Touchette, “Nonequilibrium Markov processes conditioned on large deviations,” Annales Henri Poincare, vol. 16, no. 9, pp. 2005–2057, 2015.
  * (17) R. L. Jack, “Ergodicity and large deviations in physical systems with stochastic dynamics,” European Physical Journal B, vol. 93, no. 4, p. 74, 2020.
  * (18) G. Carugno, P. Vivo, and F. Coghi, “Graph-combinatorial approach for large deviations of Markov chains,” arXiv:2201.00582, 2022.
  * (19) C. De Bacco, A. Guggiola, R. Kühn, and P. Paga, “Rare events statistics of random walks on networks: localisation and other dynamical phase transitions,” Journal of Physics A: Mathematical and Theoretical, vol. 49, no. 18, p. 184003, 2016.
  * (20) F. Coghi, J. Morand, and H. Touchette, “Large deviations of random walks on random graphs,” Physical Review E, vol. 99, no. 2, p. 022137, 2019.
  * (21) R. Gutierrez and C. Perez-Espigares, “Generalized optimal paths and weight distributions revealed through the large deviations of random walks on networks,” Physical Review E, vol. 103, no. 2, p. 022319, 2021.
  * (22) S. Whitelam, “Large deviations in the presence of cooperativity and slow dynamics,” Physical Review E, vol. 97, p. 62109, 2018.
  * (23) S. Whitelam and D. Jacobson, “Varied phenomenology of models displaying dynamical large-deviation singularities,” Physical Review E, vol. 103, no. 3, p. 032152, 2021.
  * (24) K. Sekimoto, Stochastic Energetics, vol. 799 of Lecture Notes in Physics.  Springer Berlin Heidelberg, 2010.
  * (25) M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, “Large fluctuations and optimal paths in chemical kinetics,” The Journal of Chemical Physics, vol. 100, no. 8, p. 5735, 1998.
  * (26) F. Coghi, Large deviation methods for the study of nonequilibrium systems: variational and spectral approaches.  PhD thesis, Queen Mary, University of London, 2021.
  * (27) J. Hoppenau, D. Nickelsen, and A. Engel, “Level 2 and level 2.5 large deviation functionals for systems with and without detailed balance,” New Journal of Physics, vol. 18, no. 8, p. 083010, 2016.
  * (28) P. Whittle, “Some distribution and moment formulae for the Markov chain,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 17, no. 2, pp. 235–242, 1955.
  * (29) R. Dawson and I. J. Good, “Exact Markov probabilities from oriented linear graphs,” The Annales of Mathematical Statistics, vol. 28, no. 4, pp. 946–956, 1957.
  * (30) L. A. Goodman, “Exact probabilities and asymptotic relationships for some statistics from $m$-th order Markov chains,” The Annals of Mathematical Statistics, vol. 29, no. 2, pp. 476–490, 1958.
  * (31) P. Billingsley, “Statistical methods in Markov chains,” The Annales of Mathematical Statistics, vol. 32, no. 1, pp. 12–40, 1961.
  * (32) I. Csiszár, T. M. Cover, and B. S. Choi, “Conditional Limit theorems under Markov conditioning,” IEEE Transactions on Information Theory, vol. 33, no. 6, pp. 788–801, 1987.
  * (33) M. Polettini, “BEST statistics of Markovian fluxes: a tale of Eulerian tours and Fermionic ghosts,” Journal of Physics A: Mathematical and Theoretical, vol. 48, no. 36, p. 365005, 2015.
  * (34) R. Chetrite and H. Touchette, “Variational and optimal control representations of conditioned and driven processes,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2015, no. 12, p. P12001, 2015\.
  * (35) H. Touchette, “Introduction to dynamical large deviations of Markov processes,” Physica A: Statistical Mechanics and its Applications, vol. 504, pp. 5–19, 2018.
  * (36) L. Causer, M. C. Bañuls, and J. P. Garrahan, “Finite time large deviations via matrix product states,” Physical Review Letters, vol. 128, no. 9, p. 090605, 2022.
  * (37) I. Adam, D. Fanelli, T. Carletti, and G. Innocenti, “Reactive explorers to unravel network topology,” The European Physical Journal B 2019 92:5, vol. 92, no. 5, pp. 1–8, 2019.
  * (38) T. Carletti, M. Asllani, D. Fanelli, and V. Latora, “Nonlinear walkers and efficient exploration of congested networks,” Physical Review Research, vol. 2, no. 3, p. 033012, 2020.

##  Appendix A Details on ‘Bulk-dangling model’

The exact form of the function $h$ appearing in (26) reads

$\begin{split}h(s,N,\nu_{12}&,\nu_{23},\nu_{34},\nu_{44})=\nu_{12}\left(\frac{3s}{\bar{k}}-\log
2\right)+\nu_{23}\left(\frac{Ns}{\bar{k}}-\log(2(N-2))\right)+\nu_{34}(N-3)\left(\frac{(2N-5)s}{\bar{k}}-\log((N-3)(N-2))\right)+\\\
&+\nu_{44}(N-3)(N-4)\left(\frac{(N-3)s}{\bar{k}}-\log(N_{3})\right)+\nu_{12}\log\left(\frac{\nu_{12}+\nu_{23}}{\nu_{12}}\right)+\nu_{23}\left(\log\left(\frac{\nu_{12}+\nu_{23}}{\nu_{23}}\right)+\log\left(\frac{\nu_{23}+(N-3)\nu_{34}}{\nu_{23}}\right)\right)+\\\
&+\nu_{34}(N-3)\left(\log\left(\frac{\nu_{23}+(N-3)\nu_{34}}{\nu_{34}}\right)+\log\left(\frac{\nu_{34}+(N-4)\nu_{44}}{\nu_{34}}\right)\right)+\nu_{44}(N-3)(N-4)\log\left(\frac{\nu_{34}+(N-4)\nu_{44}}{\nu_{44}}\right)+\\\
&+\epsilon\left(1-2\nu_{12}-2\nu_{23}-2(N-3)\nu_{34}-(N-3)(N-4)\nu_{44}\right)\,.\end{split}$
(62)

The minimizers $\nu^{*}=(\nu^{*}_{12},\nu^{*}_{23},\nu^{*}_{34},\nu^{*}_{44})$
of the action (26) are explicitly given by

$\displaystyle\begin{split}\nu_{12}(s,N,\epsilon)&=\frac{a(s,N,\epsilon)}{(a(s,N,\epsilon)-1)}(N-3)\frac{b(s,N,\epsilon)}{(b(s,N,\epsilon)-1)}\times\\\
&\times\frac{((-1+a(s,N,\epsilon))(-1+b(s,N,\epsilon))(-1+c(s,N,\epsilon)(-4+N)))}{((-2+a(s,N,\epsilon)(-1+b(s,N,\epsilon))(-2+c(s,N,\epsilon)(-4+N))+(1+b(s,N,\epsilon))c(s,N,\epsilon)(-4+N))(-3+N))}\end{split}$
(63)
$\displaystyle\begin{split}\nu_{23}(s,N,\epsilon)&=(N-3)\frac{b(s,N,\epsilon)}{(b(s,N,\epsilon)-1)}\times\\\
&\times\frac{((-1+a(s,N,\epsilon))(-1+b(s,N,\epsilon))(-1+c(s,N,\epsilon)(-4+N)))}{((-2+a(s,N,\epsilon)(-1+b(s,N,\epsilon))(-2+c(s,N,\epsilon)(-4+N))+(1+b(s,N,\epsilon))c(s,N,\epsilon)(-4+N))(-3+N))}\end{split}$
(64)
$\displaystyle\begin{split}\nu_{34}(s,N,\epsilon)&=\frac{((-1+a(s,N,\epsilon))(-1+b(s,N,\epsilon))(-1+c(s,N,\epsilon)(-4+N)))}{((-2+a(s,N,\epsilon)(-1+b(s,N,\epsilon))(-2+c(s,N,\epsilon)(-4+N))+(1+b(s,N,\epsilon))c(s,N,\epsilon)(-4+N))(-3+N))}\end{split}$
(65)
$\displaystyle\begin{split}\nu_{44}(s,N,\epsilon)&=\frac{c(s,N,\epsilon)}{((N-4)c(s,N,\epsilon)-1)}\times\\\
&\times\frac{((-1+a(s,N,\epsilon))(-1+b(s,N,\epsilon))(-1+c(s,N,\epsilon)(-4+N)))}{((-2+a(s,N,\epsilon)(-1+b(s,N,\epsilon))(-2+c(s,N,\epsilon)(-4+N))+(1+b(s,N,\epsilon))c(s,N,\epsilon)(-4+N))(-3+N))}\end{split}$
(66)

The inequality constraints that select the physical solution of (36) are

$\displaystyle\epsilon$ $\displaystyle>\frac{3s}{2\bar{k}}-\frac{\log 2}{2}$
(67) $\displaystyle\epsilon$
$\displaystyle>-\frac{1}{2}\log\left(\frac{2(N-2)}{(N-2)e^{3\frac{s}{\bar{k}}}+e^{N\frac{s}{\bar{k}}}}\right)$
(68) $\displaystyle 0$
$\displaystyle>\tau^{2}(N-3)+\tau(N-2)(N-4)e^{-(N-2)\frac{s}{\bar{k}}}\tau-(N-2)(N-3)e^{-(2N-5)\frac{s}{\bar{k}}}$
(69) $\displaystyle\epsilon$
$\displaystyle>(N-3)\frac{s}{\bar{k}}-\log(N-3)+\log(N-4).$ (70)