††thanks: Formerly Dana Vaknin

# Percolation in heterogeneous spatial networks with long-range interactions

Guy Amit<EMAIL_ADDRESS>Department of Physics, Bar-Ilan University, Ramat
Gan, Israel  Dana Ben Porath  Department of Physics, Bar-Ilan University,
Ramat Gan, Israel  Sergey V. Buldyrev  Department of Physics, Yeshiva
University, New York, USA  Amir Bashan<EMAIL_ADDRESS>Department of
Physics, Bar-Ilan University, Ramat Gan, Israel

###### Abstract

We study the emergence of a giant component in a spatial network where the
distribution of the metric distances between the nodes is scale-invariant, and
the interaction between the nodes has a long-range power-law behavior. The
nodes are positioned in the metric space using a Lèvy flight procedure, with
an associated scale-invariant step probability density function, and is then
followed by a process of connecting each pair of nodes with a probability
function that depends on the distance between them. A natural way to analyze
the system is to consider the total probability for an edge between steps in
term of their indexes, by summing over their possible positions. By doing so,
a correspondence is found between this model and a model of percolation in a
one-dimensional lattice with long-range interactions, which allows the
identification of the conditions for which a percolation transition is
possible. We find that the emergence of a giant component and percolation
transitions is determined by a complicated phase diagram, that exhibits a
transition from weak long-range interactions to strong long-range
interactions.

Introduction.–In this work we ask and answer the following question: If the
nodes of a spatial network are embedded in a metric space in a scale-invariant
manner, and the interactions between the nodes are governed by a power-law
relationship, under what conditions will a giant component emerge?

Spatial networks Barthélemy (2011) are a class of networks where the nodes are
embedded in a metric space, meaning that the nodes have associated coordinates
that allow to define distances between them. Typically, the edges between the
nodes are constructed such that a pair of nodes with a small distance between
them have a larger probability of being connected with an edge compared with
nodes which are far apart. Many real-world systems can be modeled as spatial
networks, including transportation systems Zeng et al. (2019, 2020); Cogoni
and Busonera (2021) power-grids Yang et al. (2017), communications Lambiotte
et al. (2008), urban streets Wu et al. (2004); Zhang et al. (2019)
infrastructure Rinaldi et al. (2001); Chang (2009); Hines et al. (2010) and
brain activity Fornito et al. (2016); Breakspear (2017). It is of particular
importance to know whether a network contains a Giant Component (GC), i.e., a
connected cluster of nodes whose size, $M$, remains greater than $cN$, where
$0<c\leq 1$ is some constant, for any $N\to\infty$. The existence of a GC in a
network signifies the ability of the network to function properly Molloy and
Reed (1998); Strogatz (2001); Newman et al. (2001); Newman (2002).

Here, we consider scale-invariant spatial networks in term of the spatiality,
i.e., the metric distance between the nodes follows a power-law distribution.
To mitigate confusion, we emphasize that when we say scale-invariant we refer
to the metric distance between the nodes, and we reserve the term scale-free
to the degree distribution Aldous and Ganesan (2013). We model the spatial
placement of the nodes in the network using a Lèvy flight process, a
generalization of Brownian motion, that is characterized by infinite variance,
which was used to model many physical phenomena Shlesinger et al. (1995), and
in particular, foraging problems Viswanathan et al. (2000); Boyer et al.
(2006); Viswanathan et al. (2008); Codling et al. (2008); Sims et al. (2008);
James et al. (2011).

In addition to the scale-invariant placement of the nodes, the interactions
are Long-Ranged (LR), i.e., the bond probability is taken from a slowly
decaying power-law function of the distance between them. Examples of long-
range interactions in physical systems can be found in a wide range of
subjects, such as astrophysics, plasma physics, hydrodynamics, atomic physics,
and nuclear physics Campa et al. (2009, 2014); Defenu et al. (2021). Long-
range interactions in complex networks were also used to model disease
spreading Solomon et al. (2000); Newman (2002), social aggregation Solomon et
al. (2000); Frasco et al. (2014) and financial transactions Stauffer (2001).

The percolation problem of geometric networks with Short-Range (SR)
interactions has been the subject of many recent studies Gilbert (1961); Dall
and Christensen (2002); Penrose (2003); Xie and Rogers (2016); Dettmann and
Georgiou (2016); Murphy et al. (2018); Plaszczynski et al. (2022); Gross and
Havlin (2022). In addition, the seminal problem of percolation on a lattice
with long-range interactions has also been extensively studied Kac and
Thompson (1969); Dyson (1969); Stephen and Aharony (1981); Schulman (1983);
Zhang et al. (1983); Imbrie and Newman (1988); Aizenman and Barsky (1987);
Newman and Schulman (1986); Coppersmith et al. (2002); Moukarzel (2006);
d’Auriac and Iglói (2016); Gori et al. (2017). However, real-world spatial
networks are typically heterogeneous (in marked contrast with lattices) and
can have long-range interactions. Therefore, it is important to understand the
emergence of a GC in scenarios that involve both of these conditions.

Model.–We begin by detailing the process of building the network. A graphical
description is presented in Fig. 1. The $N$ nodes of the network are embedded
in a $d$-dimensional space with a radially-symmetric Lèvy flight process, such
that in each step the direction is picked from a uniform distribution, and the
step length $r$ is chosen from a distribution that, for large $r$, has a
power-law behavior $\sim\frac{1}{r^{1+\alpha}}$. The probability density
function of the position of the steps is therefore given by Wolfgang and
Baschnagel (2013)

$f(\mathbf{x})\sim\frac{1}{\lvert\mathbf{x}\rvert^{d+\alpha}},\quad\frac{\lvert\mathbf{x}\rvert}{\gamma}\gg
1,$ (1)

where $\mathbf{x}\in\mathbb{R}^{d}$, $\alpha>0$ is the stability parameter,
and the probability density function has an associated scaling factor $\gamma$
such that $\lvert\mathbf{x}\rvert/\gamma$ is dimensionless. When $\alpha<2$
the distances between the nodes is scale-invariant and a mean distance is
undefined.

Figure 1: Graphical description of the model. (a) A Lèvy flight process
positions the nodes in the $d$-dimensional metric space. The step-lengths
(black dashed line) are taken from a scale-invariant probability density
function $f(\mathbf{x})$, with an associated stability parameter $\alpha$. (b)
The nodes are then connected with edges according to a probability function
$p(r)$ which depends on the distance between them and has an associated long-
range interaction parameter $\mu$. The different colored nodes represent
different separated clusters. The red cluster is the largest connected
component (c) Visualization of the function $f(\mathbf{x})$ with its long-tail
behavior that depends on the dimensionality $d$ and on the stability parameter
$\alpha$. (d) Visualization of the function $p(r)$ with its associated length
$L$ and the long-range interaction parameter $\mu$. (e) We calculate the total
probability that step $i$ is connected to step $j$ by summing over all
possible locations of step $j$. We find a correspondence between the existence
of a GC in this problem and the problem of percolation in a one-dimensional
lattice with long-range interactions.

Then, each pair of nodes, $i$ and $j$, are connected with an edge with a
probability that depends on the distance between them
$p(\lvert\mathbf{x}_{i}-\mathbf{x}_{j}\rvert)$. We choose the following long-
range probability function,

$p(r_{ij})=\begin{cases}\rho,&r_{ij}\leq L\\\
\rho\left(\frac{r_{ij}}{L}\right)^{\mu-1},&r_{ij}>L\end{cases},\quad\mu\leq 1$
(2)

where $r_{ij}=\lvert\mathbf{x}_{i}-\mathbf{x}_{j}\rvert$, $L>0$ and the
parameter $0<\rho\leq 1$ is the probability that two nodes, which are at a
distance $r_{ij}\leq L$ apart, are connected with an edge. The parameter $\mu$
is responsible for long-range interactions (see Fig. 1(d)). When $\mu=1$ the
probability for an edge between any two nodes is equal to $\rho$ independently
of the distance between them, i.e., an Erdős–Rènyi network. When $\mu=-\infty$
the function $p(r_{ij})$ is given by

$p(r_{ij})=\begin{cases}\rho,&r_{ij}\leq L\\\
0,&r_{ij}>L\end{cases},\quad\mu=-\infty.$ (3)

In this case, the nodes are connected with a probability $\rho$ if the
distance between them is less than or equal to $L$, and 0 otherwise, i.e., a
soft geometric network Gilbert (1961); Dettmann and Georgiou (2016); Penrose
(2016); Wilsher et al. (2020); Erba et al. (2020); Plaszczynski et al. (2022).

Results.–The probability for an edge between node $i$ and node $j$ can be
calculated explicitly. We denote the total probability that step $i$ is
connected with an edge to step $j$ after $\lvert i-j\rvert$ steps with
$\mathcal{P}_{ij}^{\text{SE}}$ (where SE stands for Spatial Embedding). To
find $\mathcal{P}_{ij}^{\text{SE}}$, we sum over all possible positions of
node $j$, $\mathbf{x}_{j}$, whilst treating the position of node $i$,
$\mathbf{x}_{i}$, as a constant (assuming $j>i$ without loss of generality),
i.e.,

$\displaystyle\mathcal{P}_{ij}^{\text{SE}}(\mathbf{x}_{i})=\int\dots\int
p(\lvert\mathbf{x}_{i}-\mathbf{x}_{j}\rvert)f(\mathbf{x}_{j}-\mathbf{x}_{j-1})\dots
f(\mathbf{x}_{i+1}-\mathbf{x}_{i})\differential[d]{\mathbf{x}_{j}}\dots\differential[d]{\mathbf{x}_{i+1}}=(p\,{\textstyle*d}\,f^{{*d}\lvert
i-j\rvert})(\mathbf{x}_{i}),$ (4)

where ${*d}$ is the $d$-dimensional convolution, $f^{{*d}\lvert i-j\rvert}$ is
a $\lvert i-j\rvert$-times $d$-dimensional convolution power of
$f(\mathbf{x})$, and the integrals are over the entire volume of the
$d$-dimensional space.

To evaluate (4), we apply a Fourier transform Wolfgang and Baschnagel (2013).
We define $\tilde{f}(\boldsymbol{\nu})$ and $\tilde{p}(\boldsymbol{\nu})$ as
the Fourier transform of $f(\mathbf{x})$ and $p(\lvert\mathbf{x}\rvert)$,
respectively.

$\displaystyle\begin{split}&\tilde{f}(\boldsymbol{\nu})=\int_{\mathbb{R}^{d}}f(\mathbf{x})e^{-i\boldsymbol{\nu}\cdot\mathbf{x}}\differential[d]{\mathbf{x}},\\\
&\tilde{p}(\boldsymbol{\nu})=\int_{\mathbb{R}^{d}}p(\mathbf{x})e^{-i\boldsymbol{\nu}\cdot\mathbf{x}}\differential[d]{\mathbf{x}},\end{split}$
(5)

where $\boldsymbol{\nu}\in\mathbb{R}^{d}$. With these definitions at hand, the
Fourier transform of $\mathcal{P}_{ij}^{\text{SE}}(\mathbf{x}_{i})$ is

$\tilde{\mathcal{P}}_{ij}^{\text{SE}}(\boldsymbol{\nu})=\tilde{p}(\boldsymbol{\nu})(\tilde{f}(\boldsymbol{\nu}))^{\lvert
i-j\rvert}.$ (6)

Since $f(\mathbf{x})$ and $p(\lvert\mathbf{x}\rvert)$ are independent of the
step index, when calculating the probability $\mathcal{P}_{ij}^{\text{SE}}$ we
can set $\mathbf{x}_{i}=0$ at the origin without loss of generality, i.e.,
$\mathcal{P}_{ij}^{\text{SE}}=\mathcal{P}_{ij}^{\text{SE}}(\mathbf{x}_{i}=0)$.
Finally, $\mathcal{P}_{ij}^{\text{SE}}$ can be calculated via the inverse
Fourier transform

$\mathcal{P}_{ij}^{\text{SE}}=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\tilde{p}(\boldsymbol{\nu})(\tilde{f}(\boldsymbol{\nu}))^{\lvert
i-j\rvert}\differential[d]{\boldsymbol{\nu}}.$ (7)

Using (7) we will find the conditions for which percolation is possible.
Before continuing, we will do a short interlude and talk about percolation on
a one-dimensional lattice with long-range interactions.

Consider a one-dimensional lattice, where nodes are embedded on the lattice
points, and each pair of nodes $i$ and $j$ are connected with an edge with a
probability $\mathcal{P}_{ij}^{\text{1DL}}$ (see Fig. 1(e); 1DL stands for
1-Dimensional Lattice). A well studied scenario of this problem involves long-
range interactions of the form

$\mathcal{P}_{ij}^{\text{1DL}}=\frac{C}{\lvert i-j\rvert^{1+\sigma}},\quad
0<C\leq 1.$ (8)

With this form of interactions, GC exists under the following conditions Gori
et al. (2017). For $\sigma\geq 1$, a GC exists if and only if $C=1$. This
means that a GC exists only if nearest-neighboring nodes are guaranteed to be
connected with an edge. For $0<\sigma<1$, a GC exists if $C$ is larger than
some critical value $C_{c}$. Schulman Schulman (1983) showed that $C_{c}\geq
1/2\zeta(1+\sigma)$ where $\zeta(x)$ is the Riemann zeta function,
$\zeta(x)=\sum_{n=1}^{\infty}\frac{1}{n^{x}}$, while Gori et al Gori et al.
(2017) found precise numerical estimates of $C_{c}$. For $\sigma\leq 0$ a GC
exists for arbitrarily small $C$.

Coming back to our random walk model, although the nodes are embedded in a
$d$-dimensional space, Eq. (7) defines the probability that two nodes are
connected in term of their indexes, not their position. If the long-term
behavior of (7) is similar to (8), we can find the conditions under-which a GC
exists, by corresponding the random walk problem in $d$ dimensions to the
percolation threshold of a one-dimensional lattice. Indeed, with the use of
asymptotic expansion, we find that the SE model has a similar power-law
probability of $\mathcal{P}_{ij}^{\text{1DL}}$, with $\sigma$ being a function
$\alpha$, $\mu$ and $d$, and $C$ is proportional to $\rho$,

$\mathcal{P}_{ij}^{\text{SE}}\Longleftrightarrow\mathcal{P}_{ij}^{\text{1DL}},\quad\text{for
}\lvert i-j\rvert\gg 1.$ (9)

In other words, although the nodes are embedded in a $d$-dimensional space,
the probability for an edge between any two nodes depends simply on the order
in which the nodes appeared in the Lèvy flight path, which can be thought of
as a one dimensional lattice such that each lattice point represents an index
of the path. In this correspondence, the factor $\rho$ from the random walk
model (7) takes the role of the constant $C$ from the 1D lattice model (8). We
expect that for certain choices of $\mu$, $\alpha$ and $d$, a critical value,
$0<\rho_{c}<1$, can exist such that a GC will emerge as $\rho$ is increased
from $\rho<\rho_{c}$ to $\rho>\rho_{c}$.

We solve Eq. (7) with asymptotic expansion, i.e., for large $\lvert
i-j\rvert\gg 1$. The full derivation is presented in the Supplementary
Material (SM). A key step in the analysis is that, due the Generalized Central
Limit Theorem Gnedenko and Kolmogorov (1954), the term
$(\tilde{f}(\boldsymbol{\nu}))^{\lvert i-j\rvert}$ on the r.h.s. of (7), for
$\lvert i-j\rvert\gg 1$, is given by

$\displaystyle\begin{split}(\tilde{f}(\boldsymbol{\nu}))^{\lvert
i-j\rvert}&=\exp\left\\{-\lvert\gamma\nu\rvert^{\alpha}\lvert
i-j\rvert\right\\},\quad\text{for }0<\alpha<2,\\\
(\tilde{f}(\boldsymbol{\nu}))^{\lvert
i-j\rvert}&=\exp\left\\{-\lvert\gamma\nu\rvert^{2}\lvert
i-j\rvert\right\\},\quad\text{for }\alpha\geq 2,\end{split}$ (10)

where $\nu=\lvert\boldsymbol{\nu}\rvert$.

For the geometric case ($\mu=-\infty$) the first order term of
$\mathcal{P}_{ij}^{\mathrm{SE}}$ is given by

$\mathcal{P}_{ij}^{\text{SE}}(\mu=-\infty)\sim\frac{d}{2^{d}\alpha}\frac{\left(\frac{L}{\gamma}\right)^{d}\Gamma\left(\frac{d}{\alpha}\right)}{(\Gamma\left(\frac{d}{2}+1\right))^{2}}\times\frac{\rho}{\lvert
i-j\rvert^{d/\alpha}}.$ (11)

This relation already hints to the connection with the 1D lattice case, as we
can observe that the factor $1+\sigma$ in (8) is replaced here with
$d/\alpha$. Indeed, even in the geometric case, which is characterized by
short-range interactions in the $\mathbb{R}^{d}$ metric space, there are long-
range interactions in the index-space. However, for this case we found no
numerical evidence of percolation transition or a critical value
$0<\rho_{c}<1$. In general $\rho_{c}>1$, unless $\alpha\geq d$, in-which-case
$\rho_{c}=0$ (which coincides with a return probability of 1 Wolfgang and
Baschnagel (2013)). For further discussion see the SM.

In the long-range case ($\mu>-\infty$), we find the following form of
$\mathcal{P}_{ij}^{\mathrm{SE}}$

$\begin{cases}\begin{aligned}
&\mathcal{P}_{ij}^{\text{SE}}\sim\frac{\rho}{\lvert
i-j\rvert^{\frac{1-\mu}{\alpha}}},\quad&&\text{for }1-d<\mu<1,\\\
&\mathcal{P}_{ij}^{\text{SE}}\sim\frac{\rho}{\lvert
i-j\rvert^{d/\alpha}},\quad&&\text{for }\mu<1-d,\end{aligned}\end{cases}$ (12)

where $\alpha$ is replaced by $2$ if $\alpha>2$ due to (10). The derivations
are presented in the SM.

Let us consider what Eq. (12) means. The system behaves qualitatively
differently as $\mu$ is changed from $\mu>\mu^{*}$ to $\mu<\mu^{*}$, where

$\mu^{*}=1-d.$ (13)

When $\mu>\mu^{*}$ the dominating power of the long-range interaction is of
the form $\frac{1-\mu}{\alpha}$. When $\mu<\mu^{*}$ the dominating power is of
the form $d/\alpha$, exactly like in the geometric, short-range, case (11).
The case of $\mu>\mu^{*}$ is typically called strong long-range interactions,
and of $\mu<\mu^{*}$ weak long-range interactions Fisher et al. (1972); Defenu
et al. (2020). The transition from strong LR to weak LR in the case of $d=2$
is presented in Fig. 2. Corresponding examples of networks, which exhibit
cases where phase transitions are possible and impossible, are presented in
Fig. 3. A full picture of the phase diagram in the $d=2$ is presented in Fig.
4, and for other values of $d$ in the SM.

Figure 2: Transition from weak LR to strong LR interactions. (a) The
probability for an edge, $\mathcal{P}_{ij}^{\mathrm{SE}}$, as a function of
$\lvert i-j\rvert$ for different values of $\mu$. Here $\alpha=1$ and $d=2$.
In this case, from Eq. (13), $\mu^{*}=-1$. For $\mu>\mu^{*}$ the power
behavior of the probability for an edge depends on $\mu$ and is equal to
$(1-\mu)/\alpha$ (the strong LR case). For $\mu<\mu^{*}$ the power behavior is
independent of $\mu$ and is equal to $d/\alpha$ (the weak LR case). Each value
of $\mu$ is averaged over 500 realizations. (b) Comparison between the
numerically calculated power of $\mathcal{P}_{ij}^{\mathrm{SE}}$ and the
theory. For each line in (a) we calculated a simple linear fit. The bar-lines
show the confidence intervals. Here $N=10^{4}$. Figure 3: Examples of networks
with and without phase transitions. (a) A case of $\alpha=1.5$ and
$\mu=-0.25$. These parameters correspond to the mean field universal behavior
region (green area in Fig. 4), which has $\rho_{c}=0$ (i.e., no phase
transitions). Here we chose $\rho=0.05$. (b) A case of $\alpha=0.8$ and
$\mu=-0.8$. In this region a GC is guaranteed not to exist ($\rho_{c}>1$, red
area in Fig. 4). Here we chose $\rho=1$. (c) and (d) Cases of $\alpha=1.2$ and
$\mu=-0.8$. This is a region where phase transitions are possible, and in the
strong LR case (grey area in Fig. 4). For (c) $\rho=0.5$ and for (d)
$\rho=0.05$, the system exhibits a phase transition (i.e.,
$0.05<\rho_{c}<0.5$). In all cases, the bright red cluster is the largest
component.

First, let us consider the weak LR region ($\mu<\mu^{*}$). In Fig. 4 it is the
region below the dotted line of $\mu=-1$. In this case, we compare our model
with the 1D lattice model (8), by matching the powers of $\lvert i-j\rvert$
and find that

$1+\sigma=d/\alpha,\quad\text{for }\mu<\mu^{*},$ (14)

where $\alpha$ is replaced by $2$ if $\alpha>2$. This equation will define for
us the conditions for which a phase transition is possible or impossible, for
any value of $\alpha$, $\mu$, and $d$, by using the conditions for $\sigma$ in
the 1DL model Gori et al. (2017). If $\alpha\geq d$ ($\sigma\leq 0$, the green
area in Fig. 4 below the dotted $\mu=-1$ line) a GC is guaranteed to exist for
arbitrarily small value of $\rho$. In other words, $\rho_{c}=0$ and there are
no phase transitions. For $d=1$ and $d=2$ this region coincides with a return
probability of $1$ Wolfgang and Baschnagel (2013).

For $d/2<\alpha<d$ ($0<\sigma<1$, the grey zone in Fig. 4 below the dotted
$\mu=-1$ line) a phase transition is possible, i.e., there is a possibility of
a critical value $\rho_{c}$ which is between 0 and 1. A phase transition will
be apparent in this case as $\rho$ is changed from $\rho<\rho_{c}$ to
$\rho>\rho_{c}$. Finding an exact expression for the value of $\rho_{c}$ as a
function of $\alpha$, $\mu$, and $d$ is an extremely difficult task which we
will not attempt here. We note that a phase transition is possible in this
region, but not guaranteed. For some values of $\alpha$ and $\mu$ we find
$\rho_{c}=0$ whilst for others $\rho_{c}>1$. Fig. 4 shows an example where
phase transition was detected. For $\alpha\leq d/2$ ($\sigma\geq 1$, the red
zone in Fig. 4 below the dotted $\mu=-1$ line), a GC is guaranteed not to
exists, i.e., $\rho_{c}>1$ 111To be more precise, for a GC to exist in this
region each pair of consecutive steps much be connected with an edge
(corresponding to $C=1$ in the 1D lattice case). This will only be true if
$f(\mathbf{x})$ is bounded such that $f(\lvert\mathbf{x}\rvert>L)=0$ (meaning
$\alpha=\infty$) and $\rho=1$..

To summarize, in the weak LR region, where $\sigma+1=d/\alpha$ (14),

$\mu<\mu^{*}\begin{cases}\begin{aligned} &\rho_{c}=0\text{ (GC
exists)},\quad&&\text{for }\alpha\geq d\\\ &0<\rho_{c}<1\text{
possible},\quad&&\text{for }d/2<\alpha<d\\\ &\rho_{c}>1\text{ (No
GC)},\quad&&\text{for }\alpha\leq d/2.\end{aligned}\end{cases}\\\ $ (15)

In the strong LR region (Fig. 4, above the dotted $\mu=-1$ line), the
correspondence with the 1D lattice case is summarized by

$1+\sigma=\frac{1-\mu}{\alpha},\quad\text{for }\mu>\mu^{*}.$ (16)

In complete analogy to the weak LR case, we find that

$\mu>\mu^{*}\begin{cases}\begin{aligned} &\rho_{c}=0\text{ (GC
exists)},\quad&&\text{for }\alpha\geq 1-\mu\\\ &0<\rho_{c}<1\text{
possible},\quad&&\text{for }\frac{1-\mu}{2}<\alpha<1-\mu\\\ &\rho_{c}>1\text{
(No GC)},\quad&&\text{for }\alpha\leq\frac{1-\mu}{2}.\end{aligned}\end{cases}$
(17)

In Fig. 4 these regions are separated by the black diagonal lines
222Furthermore, without proving here, we note that the case of $\mu=\mu^{*}$
is characterized by interactions of the form $\frac{\log(\lvert
i-j\rvert)}{\lvert i-j\rvert^{\beta}}$, for some power $\beta$..

Along with the theoretical phase diagram, Fig. 4 also presents a simulated
construction of the phase diagram. There, we detect the existence of phase
transitions without explicitly calculating the value of $\rho_{c}$ (a
numerically expensive task), by finding the general region where the Second
Largest Component (SLC) is maximal Margolina et al. (1982). We calculate the
SLC for only three values of $\rho$, 0.05, 0.5, and 0.95. If the maximal SLC
is $<0.05$ we colored the corresponding squares green. If the maximal SLC is
$>0.95$, we colored the corresponding squares red. If the maximal SLC is found
in $\rho=0.5$, we colored the corresponding squared gold. These are the
regions where phase transitions were detected.

An accurate way to estimate $\rho_{c}$ explicitly is to use the measure
$Q_{G}$ proposed in Gori et al. (2017)

$Q_{G}=\left\langle\frac{\sum_{\mathcal{C}}\mathcal{N}_{\mathcal{C}}^{4}}{(\sum_{\mathcal{C}}\mathcal{N}_{\mathcal{C}}^{2})^{2}}\right\rangle,$
(18)

where $\mathcal{N}_{\mathcal{C}}$ is the number of nodes in cluster
$\mathcal{C}$ and the averaging is done over the different realizations.
$Q_{G}$ measures the spread of the cluster size, and acts as a Binder cumulant
of the system Binder (1981). When a GC is present, $Q_{G}\rightarrow 1$,
otherwise $Q_{G}\rightarrow 0$. The curves of $Q_{G}$ as a function of $\rho$
for a different number of nodes, $N$, will cross near the critical $\rho_{c}$
(if one exists). Examples of such plots for various choices of $\alpha$ and
$\mu$ are presented in the SM. Calculation of the critical $\rho_{c}$ using
$Q_{G}$ for various values of $\alpha$ and $\mu$ is also presented in the SM.

As we mentioned, in the SM we show the phase diagram of the system for $d=1$
to $d=9$. Note that for $d>3$, a GC can not emerge in the weak LR region,
i.e., $\rho_{c}(d>4,\mu<\mu^{*})>1$.

Figure 4: Phase diagram for the case of $d=2$. (a) Theoretical prediction of
the different regions. The diagram shows the different regions where a GC can
and cannot exist (green and red respectively), and where percolation
transition is possible (grey), as well as the strong and weak LR regions
(above and below the line $\mu=\mu^{*}$, respectively). (b) Simulated phase
diagram, colored according to detection of phase transitions with the use of
the maximal values of the second largest component. Squares are colored green
if $\rho_{c}<0.05$, red if $\rho_{c}>0.95$ and gold if $0.05<\rho_{c}<0.95$.
Here $N=100$ and the number of realizations is 100. (c) Simulations of the
size of the GC for three cases of $\mu$ and $\alpha$ marked by x in the phase
diagram plot, for various number of nodes $N$, averaged over 1000 realizations
(red: $\alpha=0.6$, $\mu=-0.25$; orange: $\alpha=0.8$, $\mu=-0.25$; green:
$\alpha=1.8$, $\mu=-0.5$). The black arrow shows the crossing of the size of
the which is near the critical value $\rho_{c}$ (inset: zoom in on the red
curves). (d) Simulations of the size of the second largest component. The
arrow marks the maximal values, which is near the critical value $\rho_{c}$.

Discussion.–In this work, we present a first and necessary step towards
answering the following question: How do structures form in spatially scale-
invariant networks with long-range interactions? This question is important as
many real-world systems can display both of these properties. For example, bus
stations in a country can have a scale-invariant distribution of distances,
but some long-range bus lines between distance cities do often exist
Shanmukhappa et al. (2018); Wang et al. (2020). For a fully connected
transport system, we need to understand how the proclivity of long-range
lines, together with the scale-invariant distribution of distances between the
stops, compete with each other to form a large structure.

Our solution relies on the way we positioned the nodes in the metric space.
Since we used a Lévy flight procedure, which is a Markovian process that
places the nodes sequentially, we were able to calculate the probability that
any two steps are connected, and compare it to the known model of 1D
percolation. This in turn produced a complicated and rich behavior of the
system. First, we were able to find conditions for which a GC will always
exist or will never exist. Then, we found conditions for which a phase
transition is possible, and we detected these phase transitions using
simulations. This combines techniques and behaviors from both long-range
interaction theory and random walk theory in beautiful harmony. It is
important to note, however, that we have no particular reason to suspect that
real-world systems are produced by a similar Lèvy flight procedure, and a more
complete theory that can handle heterogeneous scale-invariant placement of
nodes that is not procedural or Markovian is needed.

Other than the conditions for which a GC emerges, we also found a transition
from weak to strong long-range interactions. This is in line with a known
behavior of Ising model with long-range interactions. A critical value
$\sigma^{*}$ exists that separates the system into two cases, one where the
critical behavior is the same as the short-range analog, and one where the LR
interactions are relevant and the system has a phase transition with a
peculiar LR critical behavior. In this model a third region, one where mean
field universal behavior is expected, also exists Fisher et al. (1972); Dutta
and Bhattacharjee (2001); Defenu et al. (2017). This corresponds to our green
regions in the phase-diagram plots and so it is expected that mean-field
theory will be applicable there. Note that a crucial difference between our
model and the Ising model is that the parameter $\alpha$ assumes the role of
the parameter $d$ in the Ising case, making the phase-diagram plot continuous.

Many other properties of the network are crying out to be analyzed. For
example, although we detect phase transitions in the predicated region of the
phase diagram, we do not have a mathematical theory to explain how $\rho_{c}$
changes in this region, or why parts of this region do not exhibit a phase
transition. We plan to study in the future other basic properties of this
network, such as connectivity, clustering coefficients, centrality, and
betweenness measures, and to analyze the critical exponents of these phase
transitions.

We thank Shlomo Havlin, Eli Barkai, and Stas Burov for fruitful discussions
and helpful suggestions. This work was supported by the United States-Israel
Binational Science Foundation (BSF), Jerusalem, Israel (Grant No. 2020255).

## References

  * Barthélemy (2011) M. Barthélemy, Physics reports 499, 1 (2011).
  * Zeng et al. (2019) G. Zeng, D. Li, S. Guo, L. Gao, Z. Gao, H. E. Stanley, and S. Havlin, Proceedings of the National Academy of Sciences 116, 23 (2019).
  * Zeng et al. (2020) G. Zeng, J. Gao, L. Shekhtman, S. Guo, W. Lv, J. Wu, H. Liu, O. Levy, D. Li, Z. Gao, et al., Proceedings of the National Academy of Sciences 117, 17528 (2020).
  * Cogoni and Busonera (2021) M. Cogoni and G. Busonera, Physical Review E 104, L012301 (2021).
  * Yang et al. (2017) Y. Yang, T. Nishikawa, and A. E. Motter, Science 358, eaan3184 (2017).
  * Lambiotte et al. (2008) R. Lambiotte, V. D. Blondel, C. De Kerchove, E. Huens, C. Prieur, Z. Smoreda, and P. Van Dooren, Physica A: Statistical Mechanics and its Applications 387, 5317 (2008).
  * Wu et al. (2004) J. Wu, Z. Gao, H. Sun, and H. Huang, Modern Physics Letters B 18, 1043 (2004).
  * Zhang et al. (2019) L. Zhang, G. Zeng, D. Li, H.-J. Huang, H. E. Stanley, and S. Havlin, Proceedings of the National Academy of Sciences 116, 8673 (2019).
  * Rinaldi et al. (2001) S. M. Rinaldi, J. P. Peerenboom, and T. K. Kelly, IEEE control systems magazine 21, 11 (2001).
  * Chang (2009) S. E. Chang, Frontiers of Engineering 39, 36 (2009).
  * Hines et al. (2010) P. Hines, S. Blumsack, E. C. Sanchez, and C. Barrows, in _2010 43rd Hawaii International Conference on System Sciences_ (IEEE, 2010), pp. 1–10.
  * Fornito et al. (2016) A. Fornito, A. Zalesky, and E. Bullmore, _Fundamentals of brain network analysis_ (Academic Press, 2016).
  * Breakspear (2017) M. Breakspear, Nature neuroscience 20, 340 (2017).
  * Molloy and Reed (1998) M. Molloy and B. Reed, Combinatorics, probability and computing 7, 295 (1998).
  * Strogatz (2001) S. H. Strogatz, nature 410, 268 (2001).
  * Newman et al. (2001) M. E. Newman, S. H. Strogatz, and D. J. Watts, Physical review E 64, 026118 (2001).
  * Newman (2002) M. E. Newman, Physical review E 66, 016128 (2002).
  * Aldous and Ganesan (2013) D. Aldous and K. Ganesan, Proceedings of the National Academy of Sciences 110, 8782 (2013).
  * Shlesinger et al. (1995) M. F. Shlesinger, G. M. Zaslavsky, and U. Frisch (1995).
  * Viswanathan et al. (2000) G. Viswanathan, V. Afanasyev, S. V. Buldyrev, S. Havlin, M. Da Luz, E. Raposo, and H. E. Stanley, Physica A: Statistical Mechanics and its Applications 282, 1 (2000).
  * Boyer et al. (2006) D. Boyer, G. Ramos-Fernández, O. Miramontes, J. L. Mateos, G. Cocho, H. Larralde, H. Ramos, and F. Rojas, Proceedings of the Royal Society B: Biological Sciences 273, 1743 (2006).
  * Viswanathan et al. (2008) G. M. Viswanathan, E. Raposo, and M. Da Luz, Physics of Life Reviews 5, 133 (2008).
  * Codling et al. (2008) E. A. Codling, M. J. Plank, and S. Benhamou, Journal of the Royal society interface 5, 813 (2008).
  * Sims et al. (2008) D. W. Sims, E. J. Southall, N. E. Humphries, G. C. Hays, C. J. Bradshaw, J. W. Pitchford, A. James, M. Z. Ahmed, A. S. Brierley, M. A. Hindell, et al., Nature 451, 1098 (2008).
  * James et al. (2011) A. James, M. J. Plank, and A. M. Edwards, Journal of the Royal Society Interface 8, 1233 (2011).
  * Campa et al. (2009) A. Campa, T. Dauxois, and S. Ruffo, Physics Reports 480, 57 (2009).
  * Campa et al. (2014) A. Campa, T. Dauxois, D. Fanelli, and S. Ruffo, _Physics of long-range interacting systems_ (OUP Oxford, 2014).
  * Defenu et al. (2021) N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni (2021), eprint arXiv:2109.01063.
  * Solomon et al. (2000) S. Solomon, G. Weisbuch, L. de Arcangelis, N. Jan, and D. Stauffer, Physica A: Statistical Mechanics and its Applications 277, 239 (2000).
  * Frasco et al. (2014) G. F. Frasco, J. Sun, H. D. Rozenfeld, and D. Ben-Avraham, Physical Review X 4, 011008 (2014).
  * Stauffer (2001) D. Stauffer, Advances in Complex Systems 4, 19 (2001).
  * Gilbert (1961) E. N. Gilbert, Journal of the society for industrial and applied mathematics 9, 533 (1961).
  * Dall and Christensen (2002) J. Dall and M. Christensen, Physical review E 66, 016121 (2002).
  * Penrose (2003) M. Penrose, _Random geometric graphs_ , vol. 5 (OUP Oxford, 2003).
  * Xie and Rogers (2016) Z. Xie and T. Rogers, Physical Review E 93, 032310 (2016).
  * Dettmann and Georgiou (2016) C. P. Dettmann and O. Georgiou, Physical Review E 93, 032313 (2016).
  * Murphy et al. (2018) C. Murphy, A. Allard, E. Laurence, G. St-Onge, and L. J. Dubé, Physical Review E 97, 032309 (2018).
  * Plaszczynski et al. (2022) S. Plaszczynski, G. Nakamura, C. Deroulers, B. Grammaticos, and M. Badoual, Physical Review E 105, 054151 (2022).
  * Gross and Havlin (2022) B. Gross and S. Havlin, _Percolation in Spatial Networks: Spatial Network Models Beyond Nearest Neighbours Structures_ (Cambridge University Press, 2022).
  * Kac and Thompson (1969) M. Kac and C. J. Thompson, Journal of mathematical physics 10, 1373 (1969).
  * Dyson (1969) F. J. Dyson, Communications in Mathematical Physics 12, 91 (1969).
  * Stephen and Aharony (1981) M. Stephen and A. Aharony, Journal of Physics C: Solid State Physics 14, 1665 (1981).
  * Schulman (1983) L. S. Schulman, Journal of Physics A: Mathematical and General 16, L639 (1983).
  * Zhang et al. (1983) Z. Zhang, F. Pu, and B. Li, Journal of Physics A: Mathematical and General 16, L85 (1983).
  * Imbrie and Newman (1988) J. Z. Imbrie and C. M. Newman, Communications in Mathematical Physics 118, 303 (1988).
  * Aizenman and Barsky (1987) M. Aizenman and D. J. Barsky, Communications in Mathematical Physics 108, 489 (1987).
  * Newman and Schulman (1986) C. M. Newman and L. S. Schulman, Communications in Mathematical Physics 104, 547 (1986).
  * Coppersmith et al. (2002) D. Coppersmith, D. Gamarnik, and M. Sviridenko, in _Mathematics and computer science II_ (Springer, 2002), pp. 147–159.
  * Moukarzel (2006) C. F. Moukarzel, Physica A: Statistical Mechanics and its Applications 372, 340 (2006).
  * d’Auriac and Iglói (2016) J.-C. A. d’Auriac and F. Iglói, Physical Review E 94, 062126 (2016).
  * Gori et al. (2017) G. Gori, M. Michelangeli, N. Defenu, and A. Trombettoni, Physical Review E 96, 012108 (2017).
  * Wolfgang and Baschnagel (2013) P. Wolfgang and J. Baschnagel, _Stochastic Processes_ (Springer Cham, 2013).
  * Penrose (2016) M. D. Penrose, The Annals of Applied Probability 26, 986 (2016).
  * Wilsher et al. (2020) M. Wilsher, C. P. Dettmann, and A. Ganesh, Physical Review E 102, 062312 (2020).
  * Erba et al. (2020) V. Erba, S. Ariosto, M. Gherardi, and P. Rotondo, Physical Review E 102, 012306 (2020).
  * Gnedenko and Kolmogorov (1954) B. Gnedenko and A. Kolmogorov, Am. J. Math 105 (1954).
  * Fisher et al. (1972) M. E. Fisher, S.-k. Ma, and B. Nickel, Physical Review Letters 29, 917 (1972).
  * Defenu et al. (2020) N. Defenu, A. Codello, S. Ruffo, and A. Trombettoni, Journal of Physics A: Mathematical and Theoretical 53, 143001 (2020).
  * Margolina et al. (1982) A. Margolina, H. Herrmann, and D. Stauffer, Physics Letters A 93, 73 (1982).
  * Binder (1981) K. Binder, Zeitschrift für Physik B Condensed Matter 43, 119 (1981).
  * Shanmukhappa et al. (2018) T. Shanmukhappa, I. W.-H. Ho, and C. K. Tse, Physica A: Statistical Mechanics and its Applications 502, 295 (2018).
  * Wang et al. (2020) Y. Wang, Y. Deng, F. Ren, R. Zhu, P. Wang, T. Du, and Q. Du, Cities 96, 102406 (2020).
  * Dutta and Bhattacharjee (2001) A. Dutta and J. Bhattacharjee, Physical Review B 64, 184106 (2001).
  * Defenu et al. (2017) N. Defenu, A. Trombettoni, and S. Ruffo, Physical Review B 96, 104432 (2017).