# Hierarchical deposition and scale-free networks: a visibility algorithm
approach

Jonas Berx  Institute for Theoretical Physics, KU Leuven, B-3001 Leuven,
Belgium

(September 3, 2024)

###### Abstract

The growth of an interface formed by the hierarchical deposition of particles
of unequal size is studied in the framework of a dynamical network generated
by a horizontal visibility algorithm. For a deterministic model of the
deposition process, the resulting network is scale-free with dominant degree
exponent $\gamma_{e}=\ln{3}/\ln{2}$ and transient exponent $\gamma_{o}=1$. An
exact calculation of the network diameter and clustering coefficient reveals
that the network is scale invariant and inherits the modular hierarchical
nature of the deposition process. For the random process, the network remains
scale free, where the degree exponent asymptotically converges to $\gamma=3$,
independent of the system parameters. This result shows that the model is in
the class of fractional Gaussian noise (fGn) through the relation between the
degree exponent and the series’ Hurst exponent $H$. Finally, we show through
the degree-dependent clustering coefficient $C(k)$ that the modularity remains
present in the system.

##  I Introduction

The deposition of particles onto a substrate is a subject of considerable
practical importance, and it has found applications in various fields Barabási
and Stanley (1995). A plethora of deposition processes exists within the area
of statistical physics, which often reduce to variations of either the
_random_ or the _ballistic_ deposition process Family and Vicsek (1985);
Family (1986); Meakin  _et al._ (1986). These surface growth models are
understood to be the microscopic discrete versions of the continuous models of
respectively Edwards and Wilkinson Edwards and Wilkinson (1982); Vvedensky
_et al._ (1993); Buceta  _et al._ (2014) (if surface relaxation is taken into
account), and Kardar, Parisi and Zhang Kardar  _et al._ (1986); Bertini and
Giacomin (1997). Through experimental and theoretical studies performed over
the past decades, it is known that such surfaces are rough and generally
display fractal behaviour Gomes-Filho  _et al._ (2021). One assumption that is
often made in the discrete variations is that the particles are geometrically
identical in every aspect, and that the deposition steps occur _sequentially_
, i.e. one-by-one. These assumptions are not generally true for real-world
applications. Here, we will consider the hierarchical deposition model (HDM)
that was introduced in Ref. Indekeu and Fleerackers (1998), which assumes that
particles are deposited according to a power law as a function of their size
in a _synchronous_ fashion, where particles of the same size are all deposited
simultaneously in “generations”, analogous to a finite density aggregation
process where the incoming particle flux can be controlled Cheng  _et al._
(1987); Baiod  _et al._ (1988); Krug and Meakin (1991). Similar studies on
systems with power-law particle distributions investigated sequential
deposition by including a power-law distributed noise, resulting in rare-event
dominated fluctuations Zhang (1990); Hosseinabadi and Masoudi (2019). Natural
processes that show a power-law size distribution of the particles can range
from the large-scale, e.g. the deposition of pyroclastic material from
volcanic eruptions Pioli  _et al._ (2019), or the mass distribution of
fragments of an asteroid breaking up in the Earth’s atmosphere Brykina and
Egorova (2021), to small laboratory scales such as the debris of nuclei
multifragmentation Finn  _et al._ (1982). Generally, a panoply of solid
materials (rock, sand, concrete, etc.), possess a power-law size distribution
resulting from the breaking of larger objects Turcotte (1986). The HDM has
been extended to include different scaling of particle size and deposition
probability Posazhennikova and Indekeu (2000, 2001); Indekeu  _et al._ (2000),
and has been used to model the deposition of particles with spins
Posazhennikova and Indekeu (2014). Furthermore, it has found applications in
the description of bacterial biofilms Indekeu and Giuraniuc (2004); Indekeu
and Sznajd-Weron (2003) and has been used to describe coastline formation and
fractal percolation Berx  _et al._ (2021).

In this work, we map the topology of the surface formed by the HDM into a
complex network by means of a _visibility algorithm_ (VA). The essence of the
VA is to create a network from a set of data by assigning a node to each datum
and assigning edges based on the mutual visibility between two data, i.e., if
a line of visibility is not “intersected” by any intermediate data. This
algorithm was originally developed Lacasa  _et al._ (2008) to uncover
structures in time series data, uch as finding signatures of self-organised
criticality (SOC) in avalanche-based data Kaki  _et al._ (2022), and it has
found applications in astrophysics Yu  _et al._ (2012), medicine Ahmadlou  _et
al._ (2010), fluid mechanics Juniper and Sujith (2018) and several other
fields Paiva  _et al._ (2022). Generally, the VA comes in two types: the
Natural Visibility Algorithm (NVA) and the Horizontal Visibility Algorithm
(HVA). The former considers the visibility line as a direct connection between
two data and hence it is inclined with respect to the chosen common baseline.
The latter is based on whether or not two data can have mutual horizontal
visibility (illustrated in the middle panel of Fig. 1) and has no inclination
with respect to the baseline. It is not difficult to see that the graph
generated by the HVA is always a subgraph of the NVA, as horizontal visibility
implies natural visibility between two data. An extensive comparison between
the two choices is performed in Ref. Lacasa  _et al._ (2008). The visibility
algorithm has been extended to include, e.g., directed networks, which can be
used to determine time series (ir)reversibility Lacasa  _et al._ (2012),
limited penetrable (i.e., see-through) visibility networks Wang  _et al._
(2018), and multilayer networks to study multivariate data Lacasa  _et al._
(2015). A comprehensive review of VAs can be found in Ref. Nuñez  _et al._
(2012).

The setup of this work is as follows. In section II, we introduce a
deterministic version of the hierarchical deposition model and calculate
various topological quantities: the degree distribution, graph diameter,
clustering coefficient and adjacency matrix spectrum. In section III, we
extend the model to incorporate random deposition and erosion with
probabilities $P$ and $Q$, respectively. We revisit the topological quantities
and calculate them numerically. We analytically argue that the degree exponent
is $\gamma=3$ through a connection with the Hurst exponent $H$ and the
Hausdorff dimension $D_{f}$. Finally, in section IV, we present conclusions
and a future outlook.

##  II Hierarchical deterministic deposition

Before studying the random version of the HDM, let us start by introducing the
hierarchical deterministic deposition model (HDDM). This model will serve to
illustrate the main concepts in an analytical and mechanistic fashion.
Consider the deposition of rigid particles on a one-dimensional substrate,
where we assume that deposited particles are squares with different sizes. If
the deposition occurs in a viscous medium such as water, the larger particles
are deposited first, followed by the smaller particles. Additionally, we
assume that the number of particles $N(s)$ of size $s$ follows a hyperbolic
distribution

$N(s)=\lambda^{-1}N(s/\lambda)\,,$ (1)

where $\lambda\in\mathbb{N}$ and $\lambda>1$. Particles are deposited
simultaneously instead of sequentially, in contrast with, e.g., the solid-on-
solid or the ballistic deposition models. For simplicity, let us only consider
$\lambda=3$ and divide the unit interval $[0,1]$ into segments of length
$1/3$. Deposit a square with side length $1/3$ on the middle segment while
leaving the others empty. Repeating this procedure of segmentation and
subsequent deposition of particles results in a “castle” landscape as shown in
Fig. 1(a), where one step of segmentation and deposition corresponds to a
single generation $n$.

Figure 1: (a) Iterative construction of the fractal landscape for the first
three generations of the HDDM. (b) Height profile for the fractal landscape.
The shortest path length between the leftmost and central nodes is indicated
by dashed orange lines. (c) The complex network based on the HVA as defined by
equation (4). The different stages of each figure are indicated by colors:
generation $n=1$ (black), $n=2$ (dark blue), and $n=3$ (light blue)

.

In Ref. Indekeu and Fleerackers (1998), it was shown that this construction
results in a logarithmic fractal, whereby the surface length increment quickly
saturates to the constant $2/3$. This logarithmic fractal character remains
present for a random model with deposition probability $P$.

From the above construction, it is clear that not every value for the
interface height, measured with respect to the baseline, is accessible. In
particular, the set of height values (in ascending order) in a generation $n$
is the subset of the first $2^{n}$ values of the set $3^{-n}S(0,1)$, i.e.,
$H^{(n)}=3^{-n}S^{(n)}(0,1)\subset 3^{-n}S(0,1)$ where
$S(0,1)=\\{0,1,3,4,9,10,12,13,\dots\\}$ is the _Stanley sequence_ (OEIS
A005836). The height distribution function, supported on $H^{(n)}$, is then

$P^{(n)}(y=h_{i}\in H^{(n)})=\left(\frac{2}{3}\right)^{n}2^{-s(i)}\,,$ (2)

where $i$ is the index of $h$ in the set $H^{(n)}$, starting from $i=0$, and
the sequence $s(i)$ indicates the number of $1$’s in a binary expansion of
$i$, defined through the recurrence

$\displaystyle s(0)$ $\displaystyle=0\,,$ $\displaystyle s(2i)=s(i)\,,$
$\displaystyle s(2i+1)=s(i)+1\,.$ (3)

By construction, the height at the central point of the interface grows to the
maximal value $h_{max}=\sum_{i=1}^{n}3^{-i}$, which is bounded by $1/2$ for
large $n$.

We now show that the HDDM can be mapped to a scale-free network, the
horizontal visibility graph (HVG), by means of the HVA Lacasa  _et al._
(2008).

Two data points $x_{i}$ and $x_{j}$ possess mutual horizontal visibility if
the following criterion holds:

$\displaystyle\min{\\{x_{i}\,,x_{j}\\}}>x_{k}$ $\displaystyle\text{for all
$k$}\quad\text{such that}\,i<k<j\,.$ (4)

If two points $x_{i},x_{j}$ have mutual visibility, an edge exists between the
two nodes $i,j$ in the network created by the HVA. This algorithm has been
used to characterise time series Luque  _et al._ (2009) by converting them to
complex networks, enabling one to study the series combinatorically. To the
author’s knowledge, VAs have not received much attention in the context of
deposition processes, with only one work Brú  _et al._ (2014) that deduces
growth exponents for the Edwards-Wilkinson (EW), Kardar-Parisi-Zhang (KPZ) and
Molecular Beam Epitaxy (MBE) equations, as well as for the discrete Random
Deposition (RD), Random Deposition with Surface Relaxation (RDSR) and Eden
models. Another related work Kartha (2017) studies the random deposition of
patchy particles in the context of VAs. We will now study different
topological aspects of the network associated with the HDDM.

###  II.1 The degree distribution

When the HVA is applied to the HDDM, the node degree $k_{j}(n)$ can be found
for every individual node $j=1,...,3^{n}$ in a generation $n$ by summing all
entries in the $j-$th column of the adjacency matrix $A^{(n)}$, i.e.,
$k_{j}(n)=\sum_{i=1}^{3^{n}}A^{(n)}_{ij}$. After some algebra, which is
outlined in appendix A, and by making use of the hierarchical structure of
$A^{(n)}$, the $j$th node degree can be found as follows

$\begin{split}k_{j}(n)&=\left(1+\delta_{1,\alpha_{1}}\right)+\sum\limits_{k=1}^{n-1}\Bigg{\\{}\delta_{0,\alpha_{k+1}}\cdot\sum\limits_{i=0}^{2^{k-1}-1}\left(\delta_{j-\sigma_{nk},\,\Xi_{ki}+1/2}+\delta_{j-\sigma_{nk},\,\Xi_{ki}+3/2}\right)\\\
&+\delta_{2,\alpha_{k+1}}\cdot\sum\limits_{i=0}^{2^{k-1}-1}\left(\delta_{3^{k}+1-j+\sigma_{nk},\,\Xi_{ki}+1/2}+\delta_{3^{k}+1-j+\sigma_{nk},\,\Xi_{ki}+3/2}\right)+\delta_{1,\alpha_{k+1}}\cdot
2^{k}\left(\delta_{j-\sigma_{nk},\,1}+\delta_{j-\sigma_{nk},\,3^{k}}\right)\Bigg{\\}}\end{split}$
(5)

where $\sigma_{nk}$ and $\Xi_{ki}$ are defined as:

$\displaystyle\sigma_{nk}$
$\displaystyle=\sum\limits_{l=k}^{n}\alpha_{l+1}3^{l}\,,\qquad\Xi_{ki}=\frac{3^{k}}{2}+\frac{3i}{2}+\sum\limits_{m=1}^{i}\frac{3^{a(m)}}{2}\,.$
(6)

The $\alpha_{i}$’s, where $i\in\\{1,...,n\\}$ are the digits in the base 3
representation of the node index $j$. For example, for $j=17$ in generation
$n=4$, the number $j$ can be written as $17=2\cdot 3^{0}+2\cdot 3^{1}+1\cdot
3^{2}+0\cdot 3^{3}$, hence,
$\vec{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})=(2,2,1,0)$.
Furthermore, the function $a(m)$ is the so-called ruler function (OEIS
A001511) Allouche and Shallit (2003), where the $n$th term in this integer
sequence is the highest power of 2 that divides $2n$. It is clear from
equation (5) that not every possible integer value of $k$ occurs in every
generation. In particular, the number of distinct values for the degree grows
as $2n+1$ for $n>1$. In a generation $n>1$, $k$ can only take values in the
set $K^{(n)}$, where
$K^{(n)}=\left\\{1,2,\Lambda^{(n)}_{\text{odd}},\Lambda^{(n)}_{\text{even}}\right\\}$,
and where the sets of odd and even degrees are defined, respectively, as
$\Lambda^{(n)}_{\text{odd}}=\left\\{1+2^{\frac{j-1}{2}}\right\\}_{j\,\text{odd}}$
and
$\Lambda^{(n)}_{\text{even}}=\left\\{2+2^{\frac{j}{2}-1}\right\\}_{j\,\text{even}}$,
with $j=3,4,...,2n-1$. For the first generation, the set of possible degrees
is $K^{(1)}=\left\\{1,2\right\\}$.

After some algebra one can deduce the $n$th generation degree distribution
$P^{(n)}(k)$ from the node degree as described by equation (5), i.e.,

$\displaystyle P^{(n)}(1)$ $\displaystyle=\frac{2}{3^{n}}\,,$ $\displaystyle
P^{(n)}(2)=\frac{1}{3}\,,$ $\displaystyle
P^{(n)}(3)=\left(\frac{2}{3}\right)^{n-1}-\frac{2}{3^{n}}\,,$ (7)
$\displaystyle P^{(n)}(4)$
$\displaystyle=\frac{5}{9}-\left(\frac{2}{3}\right)^{n-1}+\frac{1}{3^{n}}\,,$
$\displaystyle
P^{(n)}\left(\Lambda^{(n)}_{\text{odd}}\right)=2^{\frac{1-j}{2}}\left(\frac{2}{3}\right)^{n}\,,$
$\displaystyle
P^{(n)}\left(\Lambda^{(n)}_{\text{even}}\right)=2\left(3^{-\frac{j}{2}}-2^{-\frac{j}{2}}\left(\frac{2}{3}\right)^{n}\right)\,.$

In Fig. 2, the numerically determined degree distribution is shown for $n=10$
(red crosses), together with the theoretically predicted even (odd) degree
distributions $P_{e}$ ($P_{o}$), for values of $n=10$ (black dots) and $n=20$
(blue squares). For visualisation purposes, lines are shown for the even
(full) and odd (dashed) distributions.

Note that the degree distribution behaves as a power law, i.e., $P_{e}\sim
k^{-\gamma_{e}}$, and $P_{o}\sim k^{-\gamma_{o}}$, with exponents $\gamma_{e}$
and $\gamma_{o}$. The exact values of $\gamma_{e}$ and $\gamma_{o}$ can be
calculated in the large $n$ limit to be $\gamma_{e}=\log{3}/\log{2}\approx
1.585$ and $\gamma_{o}=1$. Both values were confirmed numerically, as
indicated in the inset in Fig. 2. The value of
$\gamma_{e}=\log{(3)}/\log{(2)}$ has been found previously for the
deterministic network introduced by Barabási, Ravasz, and Vicsek Barabási  _et
al._ (2001); Iguchi and Yamada (2005). Eliminating $j$ from
$k=1+2^{\frac{j-1}{2}}$ and
$P^{(n)}\left(k\in\Lambda^{(n)}_{\text{odd}}\right)$, one finds that for large
$k$, the odd degree distribution is
$P^{(n)}\left(k\in\Lambda^{(n)}_{\text{odd}}\right)\sim(2/3)^{n}\,k^{-1}$. It
follows that this vanishes for large $n$. A similar calculation reveals that
no such behaviour is present for the even degrees, which will consequently
dominate $P(k)$ for large $n$.

Figure 2: The numerical degree distribution (red crosses) for the $n=10$ HVG.
$P_{e}(k)$ and $P_{o}(k)$ (7) are shown (respectively black full and dashed
lines) for $n=10$ (black dots) and $n=20$ (blue squares). Inset: exponents
$\gamma_{e}$ (dots) and $\gamma_{o}$ (stars) as a function of $n$.

Equipped with the exact expression for the degree distribution (7), one can
calculate a number of properties. The moments of the distribution can be found
as

$\langle k^{m}(n)\rangle=\sum_{k}k^{m}P^{(n)}(k)\,.$ (8)

Note that this expression only converges for $m<\gamma-1$, so only the mean
degree is finite for the HVG associated with the HDDM. The mean degree
$\langle k\rangle$ (where we suppressed $n$) can be determined exactly from
equation (8) for $m=1$ as $\langle
k\rangle=4\left(1-\left(2/3\right)^{n}\right)$. This expression converges to
$\langle k^{\infty}\rangle=4$, indicating that the network becomes sparse for
large $n$, a property that could also be determined from the increasing
sparsity of the adjacency matrix of the HVG. The mean degree can be used as a
measure of the irregularity of the of the underlying data, since for an
infinite periodic series of period T, it is given by

$\langle k\rangle=4\left(1-\frac{1}{2T}\right)\,,$ (9)

which becomes $\langle k\rangle=2$ for a constant series and $\langle
k\rangle=4$ for a completely aperiodic series where $T\rightarrow\infty$ Nuñez
_et al._ (2012). The surface height profile hence becomes completely irregular
in the long time limit; any periodicity in the model is washed out.

Through a similar calculation we find that the second moment $\langle
k^{2}\rangle$ diverges for large generations $n$ as $\langle k^{2}\rangle\sim
2(4/3)^{n}$, as expected for scale-free networks. Since scale-free networks
arise naturally from fractal series in the context of VAs, one needs another
measure to determine whether the underlying series is, e.g. a self-affine or
self-similar fractal Lacasa  _et al._ (2008). One such measure is the network
diameter $D$, which we now proceed to calculate exactly.

###  II.2 Network diameter and clustering coefficient

The diameter $D$ of a graph is defined as $D=\max_{ij}{l_{ij}}$, where
$l_{ij}$ is the shortest path distance between nodes $i$ and $j$. For the HVG,
the diameter $D(n)$ in generation $n$ is

$D(n)=2(2^{n}-1)\,.$ (10)

This can easily be shown by considering that in the fractal landscape the
outermost nodes $i=1$ and $j=3^{n}$ are the two most distant nodes in the
corresponding HVG. The total path distance between these nodes can hence be
split into two identical distances by virtue of the reflection symmetry about
the central node. Therefore, we only consider the left side of the landscape
for the calculation of the diameter. Assume now that we know the shortest path
distance between the outermost left node and the central node in generation
$n-1$. Then, again by virtue of the symmetry, we only need to consider the
left side of the different copies of the $n-1$ landscape, because the leftmost
node of such a copy is already visible from the central node of the copy to
the left.

This is shown graphically with blue lines in Fig. 1. The resulting diameter of
the left side of the landscape is twice the left diameter from the previous
generation plus one from the connection between the two copies. From this
construction, it is clear that $D(1)=2$, $D(2)=6$, $D(3)=14$, … Hence, the
recurrence for the total diameter is $D(n+1)=2D(n)+2$. Solving for $D(n)$
results in equation (10). It is clear that for large $n$, the diameter
increases as a power law of the number of nodes $N(n)=3^{n}$, i.e., $D(n)\sim
N(n)^{\epsilon}$, making the HVG a self-similar, scale invariant fractal
network Lacasa  _et al._ (2008). The exponent $\epsilon$ can easily be found
to be the inverse of the degree exponent for the even degree nodes
$\epsilon=\gamma_{e}^{-1}$. Such a power-law growth of the diameter is a
signature of hub repulsion in fractal networks Lacasa  _et al._ (2008); Song
_et al._ (2006), where the hubs of the network are the nodes with the highest
connectivity. In our deterministic model, these hubs are located in the
locally monotonic regions of the surface and are effectively repulsed from one
another by the intermittent nodes corresponding to the highest data.
Conversely, for stochastic self-affine series such as Brownian motion $B(t)$,
where $B(t)=\lambda^{1/2}B(t/\lambda)$, the associated VG evidences the small-
world property. Hence, the network diameter or associated quantities such as
the average path length can be used to differentiate between different kinds
of fractality in the underlying height profile.

By inspecting the clustering coefficient $C(k)$ for the nodes with degree $k$,
it becomes clear that there is a one-to-one correspondence between the
clustering coefficient of a vertex and its degree: $C(k)$ is proportional to
$k^{-1}$ and stays stationary, i.e., nodes with degree $k$ always have the
same clustering coefficient $C(k)$, independent of the generation. A simple
geometrical argument reveals that the clustering coefficient is the ratio of
the number of triangles that a node $i$ with degree $k$ is part of, normalised
by the maximal number of triangles, i.e., ${k\choose 2}$. So for e.g. $k=2$,
the number of triangles is one (except for the central node of the network,
which is not a part of any triangle). Hence, $C(k=2)=1$. For $k=3$, one of the
connections will be one of the neighbour data, while the other two form one
triangle. Hence, $C(k=3)=1/3$. In Table 1, the clustering coefficients are
listed as a function of the node degree $k$ and generation $n$.

$k$ | $C^{(1)}(k)$ | $C^{(2)}(k)$ | $C^{(3)}(k)$ | $C^{(4)}(k)$ | $C^{(5)}(k)$
---|---|---|---|---|---
2 | 0 | 1 | 1 | 1 | 1
3 | - | $1/3$ | $1/3$ | $1/3$ | $1/3$
4 | - | - | $1/2$ | $1/2$ | $1/2$
5 | - | - | $3/10$ | $3/10$ | $3/10$
6 | - | - | - | $1/3$ | $1/3$
9 | - | - | - | $7/36$ | $7/36$
10 | - | - | - | - | $1/5$
17 | - | - | - | - | $15/136$
Table 1: Clustering coefficient $C^{(n)}(k)$ for the first five generations of
the visibility graph.

We find that the nodes with even and odd degree have the following form,
respectively

$\displaystyle C\left(k\in\Lambda^{(n)}_{\text{even}}\right)$
$\displaystyle=\frac{2}{k},$ $\displaystyle
C\left(k\in\Lambda^{(n)}_{\text{odd}}\right)=\frac{2(k-2)}{k(k-1)},$ (11)

which are independent of the generation $n$. By construction, we also have
$C(k=1)=0$, as nodes with only a single connection (edge nodes) can never form
triangles. Hence, for large $k$, $C(k)\sim 2k^{-1}$, as was previously found
to indicate hierarchical clustering Ravasz  _et al._ (2002); Dorogovtsev  _et
al._ (2002). The scaling law indicates that a hierarchy of nodes with
different degrees of modularity coexists in the network. Since the local
clustering coefficient is based on a three-point relationship, i.e., it is a
function of three elements in the adjacency matrix, it expresses a local
convexity property of the hierarchical surface Donner and Donges (2012). As a
consequence, local minima of the surface tend to possess the highest
clustering coefficient, as can for example be seen for internal nodes where
$h=0$, which always have $k=2$ and thus $C=1$, since the nearest neighbours of
these nodes must be mutually visible for the node to be a local minimum. The
clustering properties of the visibility graph can consequently be used as a
measure for the local convexity of the underlying surface.

The mean clustering coefficient $\langle C^{(n)}\rangle$ can be found by

$\begin{split}\langle C^{(n)}\rangle&=\sum_{k}C(k)P^{(n)}(k)\,,\end{split}$
(12)

For large $n$, the mean clustering coefficient saturates to a constant value
$\langle C^{\infty}\rangle$. It can easily be shown that the contributions
from the odd degrees vanish for large $n$, as well as the majority of factors
from the even degrees. The terms that remain are

$\begin{split}\langle
C^{\infty}\rangle&=\lim\limits_{n\rightarrow\infty}\langle
C^{(n)}\rangle=\frac{11}{18}+8\sum\limits_{l=3}^{\infty}\frac{3^{-l}}{4+2^{l}}\end{split}$
(13)

where in the sum over the even degrees we have made the change of variables
$l=j/2$. It can be checked numerically that $\langle C^{n}\rangle$ converges
to the value $\langle C^{\infty}\rangle\approx 0.642$. This is shown in Fig. 3
together with the relations (11).

Figure 3: Numerical (black) and exact (orange) mean clustering coefficient
$\langle C^{(n)}\rangle$ as a function of $n$. The black dashed line indicates
$\langle C^{\infty}\rangle\approx 0.64185$. Inset: $C(k)\sim 2k^{-1}$ for even
(full line) and odd degrees (dashed line).

###  II.3 Adjacency matrix and the eigenvalue spectrum

The adjacency matrix $A^{(n)}$ is the matrix with values $a_{ij}=1$ if the
nodes $i$ and $j$ have mutual visibility, and value $a_{ij}=0$ otherwise. The
$A^{(n+1)}$ matrix is a block tridiagonal matrix with on the diagonal the
adjacency matrix $A^{(n)}$ from the previous generation, repeated three times.
It has the general form:

$A^{(n+1)}=\begin{pmatrix}A^{(n)}&C^{\intercal}&0\\\
C&A^{(n)}&B^{\intercal}\\\ 0&B&A^{(n)}\end{pmatrix}$ (14)

The matrix $B$ is related to $C$ by $B=JC^{\intercal}J$, where
$J=antidiag(1,...,1)$ is the $3^{n}\times 3^{n}$ exchange matrix.

For the adjacency matrix spectrum, we follow the representation used in
Flanagan  _et al._ (2019) and rescale the eigenvalue index for every
generation $n$ such that the smallest eigenvalue corresponds to 0 and the
largest corresponds to 1. This way, we can show that the spectrum converges to
a fixed shape with a hierarchical structure. In Fig. 4(a), we show the
spectrum in the rescaled coordinates for the first seven generations, where
the point size is inversely related to the generation number. Note that the
shape of the spectrum appears to be fractal, which we can surmise to be a
consequence of Cauchy’s interlacing theorem, as $A^{(n)}$ is always a proper
subgraph of $A^{(n+1)}$. However, the exact shape of the curve eludes any
analytic description so far and will be left for future work. The largest
eigenvalue $\lambda_{max}$ of the adjacency matrix can be numerically
determined to depend on the system size $N=3^{n}$ as $\lambda_{max}\sim
N^{\frac{1}{4}}\sim 3^{\frac{n}{4}}$. In Fig. 4(b), the eigenvalue
distribution $\rho(\lambda)$ is drawn for $n=8$, showing that the distribution
possesses power-law tails where $\rho(\lambda)\sim|\lambda|^{4.04}$. The
scaling of the largest eigenvalue and the tail of the eigenvalue spectrum is
very similar to the results previously found for the Barabási-Albert scale-
free network Goh  _et al._ (2001).

Figure 4: (a) Eigenvalue spectrum $\\{\lambda_{i}\\}$ as a function of the
rescaled coordinate $i$ for generations $n=1,...,7$. The spectrum converges to
a fixed shape for large $n$. Inset: zoomed image around the central part of
the spectrum. The fine structure is clearly visible. (b) Eigenvalue
distribution $\rho(\lambda)$ for $n=8$ in semilog scale.

The largest eigenvalue of a visibility graph’s adjacency matrix has been
proposed as a measure to characterise the complexity of the associated
sequence by means of the so-called Graph Index Complexity (GIC) Kim and
Wilhelm (2008), which is a rescaled version of $\lambda_{max}$. It has been
used for example as a measure to discriminate between randomness and chaos in
the underlying data of the graph Fioriti  _et al._ (2012); Flanagan  _et al._
(2019). In this work, we will not elaborate on this subject further but leave
a detailed analysis of the spectral properties for future investigation.

We now move on from the HDDM and introduce an element of randomness into the
process. We repeat our study of the graph-theoretical quantities.

##  III Hierarchical random deposition

Let us generalise the model introduced in the previous section by introducing
an element of randomness. In a generation $n$, a number $\lambda^{n}$ of
disjoint sites are available. These sites are either located on the substrate,
if no deposition event has occurred on this site prior to generation $n$, or
on previously deposited particles. A square particle with linear size
$s=\lambda^{-n}$ is deposited on every site with a probability $P\in[0,1]$.
Additionally, we allow for the erosion of a block of similar size that is
located on the interface with probability $Q\in[0,1]$, where the total
probability of either deposition or erosion is bounded, i.e., $0<P+Q<1$.
Erosion events thus remove material that has either been deposited in prior
generations, or is part of the substrate. Lateral attachment of incoming
particles, i.e., ballistic deposition, is not allowed and consequently
overhangs cannot be formed. This process has been named the hierarchical
random deposition model (HRDM). The parameters $P$ and $Q$ allow us to tune
the amount of noise in the growth process. For $Q=0$ and $P\rightarrow 0$
(i.e., the dilute limit), classical sequential deposition is recovered Krug
and Meakin (1991).

The HVA can now be applied in exactly the same manner as was outlined in
section II. In Fig. 5, we illustrate the landscape formed by random deposition
and the horizontal visibility lines between the data, as well as the
associated complex network. Due to the randomness introduced in the HRDM,
analytical progress is challenging so we will proceed numerically. However, we
show analytically that the fractal nature of the set of admissible height
values remains fractal (for large $n$) even when randomness is introduced into
the model.

###  III.1 The degree distribution

Based on the seminal work on HVGs performed in Ref. Luque  _et al._ (2009),
one would expect that the degree distribution becomes exponential due to the
random nature of the deposition. There are, however, some caveats to such a
claim. In the following, we show that the degree distribution is scale-free
with exponent $\gamma=3$.

Figure 5: The landscape formed by hierarchical random deposition with $P=0.5$,
$Q=0.1$, $n=3$ and $\lambda=3$ (left) and the associated visibility lines
(right). The inset in the right figure is the complex network formed by
applying the horizontal visibility criterion (4).

In their initial work Luque  _et al._ (2009), the authors proved that the
degree distribution associated with a bi-infinite (time) series generated from
an uncorrelated random variable $X$ with an absolutely continuous probability
density function $f(x)$ with $x\in[0,1]$ has the following exponential form:

$P(k)=\frac{1}{3}\left(\frac{2}{3}\right)^{k-2}=\frac{1}{3}\mathrm{e}^{-\mu_{c}(k-2)}\,,$
(15)

where $\mu_{c}=\ln{(3/2)}$. Note that for finite $n$, the set of admissible
heights in the hierarchical deposition model is also finite, and the height
distribution function is hence discrete. Therefore, one cannot a priori expect
an exponential degree distribution for $n<\infty$.

However, for uncorrelated random variables taking discrete integer values up
to a number $m$, an exact formula for $P(k)$ can be found, which does not
possess a generic algebraic closed form Lacasa (2016). For
$m\rightarrow\infty$, the degree distribution reduces to the exponential form
that was found for the continuous case. The premise of the authors’ results
was based on the equispacedness of the support of the distribution, i.e., for
the field $\mathbb{F}=(h_{1},h_{2},\dots,h_{n})\subset\mathbb{N}$, it holds
that $h_{i}=h_{i-1}+c$, $c\in\mathbb{N}^{+}$, for $m<\infty$, since the HVG is
invariant under monotonic transformations in the series. For asymptotically
unbounded support, i.e., the continuous case, one can interpret an
unequispaced subset of the natural numbers simply as the set of natural
numbers (equispaced) where some elements have probability zero. Since in the
continuous case the marginal distribution does not matter, one recovers the
exponential distribution (15) even when the support seems to be unequispaced.
For finite $m$, however, it is not clear how the unequispacedness affects the
shape of the degree distribution.

In the HRDM, however, the height of a site can take discrete values that, when
scaled with $\lambda^{n}$, are elements of a growing subset of $\mathbb{Z}$,
where the increments between elements are unequispaced. In appendix B, we
calculate that in the long-time limit, $n\rightarrow\infty$, the cumulative
height distribution
$\mathcal{F}(h)=\lim_{n\rightarrow\infty}\mathcal{P}^{(n)}(H\leq h)$ is given
by

$\mathcal{F}(h)=\frac{1}{2\pi
i}\int_{\mathbb{R}}\mathrm{d}t\,\frac{\mathrm{e}^{it/(\lambda-1)}-\mathrm{e}^{-ith}}{t}\prod_{j=1}^{\infty}\left[P\mathrm{e}^{it\lambda^{-j}}+Q\mathrm{e}^{-it\lambda^{-j}}+(1-P-Q)\right]\,.$
(16)

Note that although $\mathcal{F}$ is continuous, the height distribution
function $\mathcal{P}(H=h)=\lim_{n\rightarrow\infty}\mathcal{P}^{(n)}(H=h)$ is
a generalisation of the Cantor distribution and is therefore not absolutely
continuous with respect to Lebesgue measure, making it a singular
distribution. The set of final accessible heights is a fractal due to the
retrodictive nature of the deposition process, i.e., every final height can be
reached by a unique path de la Torre  _et al._ (2000).

These properties constitute the major differences between our results and the
results in Ref. Lacasa (2016). Since the set of admissible (scaled) height
values is a finite and non-uniformly spaced subset of $\mathbb{Z}$ for
$n<\infty$, and their distribution depends on $P$, $Q$ and $\lambda$, the
degree distribution of the HVG associated with the hierarchical random
deposition model cannot generally be mapped one-to-one to the one found in
Ref. Lacasa (2016). Asymptotically, the height distribution is singular, so
the degree distribution cannot be expected to be exponential. Since we have
established that the set of admissible height values is a fractal for
$n\rightarrow\infty$ and $P+Q\neq 0$, the degree distribution is expected to
converge to a power law $P(k)\sim k^{-\gamma}$.

For $n\rightarrow\infty$, we can find the degree exponent of this power law
distribution exactly. In Ref. Lacasa  _et al._ (2009), the authors showed
empirically that the relation between the exponent $\gamma$ of the visibility
graph degree distribution and the Hurst exponent $H$ (alternatively, the local
roughness exponent) of the associated data series for fractional Brownian
motion (fBm) and fractional Gaussian noise (fGn) is given by respectively

$\displaystyle\gamma$ $\displaystyle=3-2H\,,$ for fBm (17)
$\displaystyle\gamma$ $\displaystyle=5-2H\,,$ $\displaystyle\text{for fGn}\,.$
(18)

For self-affine processes, the Hurst exponent is related to the fractal
dimension $D_{f}$ through $D_{f}=2-H$, leading to the relations between
$\gamma$ and $D_{f}$, i.e.,

$\displaystyle\gamma$ $\displaystyle=2D_{f}-1\,,$ for fBm (19)
$\displaystyle\gamma$ $\displaystyle=2D_{f}+1\,,$ $\displaystyle\text{for
fGn}\,.$ (20)

For the HRDM, the fractal dimension is equal to the Euclidean dimension
Indekeu and Fleerackers (1998), so the degree exponent can be either
$\gamma=1$ (fBm) or $\gamma=3$ (fGn), since $D_{f}=1$. However, the mean
degree for HVGs is in fact bounded from above by $\langle k\rangle\leq 4$
Luque  _et al._ (2009); Gutin  _et al._ (2011). With $\gamma=1$, the mean
degree would become unbounded. Hence, we can conclude that in the long time
limit, the degree exponent converges to $\gamma=3$. Alternatively, a simple
detrended fluctuation analysis (DFA) reveals that indeed $H=1$ for both the
deterministic and random HDM.

Note that since the degree exponent is directly related to the topological
dimension for the HRDM, it is independent of $P$, $Q$ and $\lambda$. The
degree exponent thus presents an alternative approach for the calculation of
the local roughness exponent and for the characterisation of the underlying
process as either fBm or fGn . Note that the concomitant global roughness
exponent cannot be determined in this manner, since surfaces can exhibit
anomalous scaling Das Sarma  _et al._ (1996); López  _et al._ (1997) where the
global and local roughness exponents are unequal. Moreover, it should be noted
that for the HRDM the roughness is a power law of only the length already in
the first generation. Consequently, there is no initial growth that can be
found as a power law of time Berx  _et al._ (2021). Hence, the usual scaling
exponents $\beta$ or $z$ are nonexistent in this model. A detailed treatment
of this problem can be found in Ref. Giuraniuc (2006).

To numerically find the degree exponent $\gamma$ of the scale-free networks
generated by the HVA, we solve the following transcendental equation for the
maximum likelihood estimator (MLE) $\hat{\gamma}$ Clauset  _et al._ (2009):

$\frac{\zeta^{\prime}(\hat{\gamma},k_{\text{min}})}{\zeta(\hat{\gamma},k_{\text{min}})}=-\frac{1}{m}\sum\limits_{i=1}^{m}\ln{k_{i}}\,,$
(21)

where $\zeta$ is the Hurwitz zeta function Olver  _et al._ (2010), $m$ is the
total number of observed values of the degrees $k$, and $k_{\text{min}}$ is
the cutoff degree above which the distribution is a discrete power law. The
prime denotes differentiation with respect to the first argument. To remove
boundary effects for small system size, we impose periodic boundary
conditions, but only on the two outer nodes, which now neighbour each other.
This fixes the minimal possible degree to be $k=2$. We combine our analysis
with a Kolmogorov-Smirnov test Clauset  _et al._ (2009) on the cumulative
degree distribution to find the optimal $\hat{\gamma}$ and $k_{\text{min}}$
that fit the simulation data. In appendix C, we list the exponents for
increasing $n$, together with the numerical errors.

In Fig. 6(a), the degree distribution is shown for increasing generations, for
$P=Q=0.25$ and $\lambda=3$, averaged over 5000 realisations, indicating that
the network is scale-free. In the inset, the exponent $\gamma$ is shown as a
function of the generation $n$. An exponential fit yields the expression
$\gamma(n)=3.00213+12.2109\,\mathrm{e}^{-0.34159n}$, which supports the
theoretical assertion that the exponent equals $\gamma=3$ for
$n\rightarrow\infty$. The dependence of $\gamma$ on $P$ and $Q$ is shown in
Fig. 7(a). While the degree exponent converges to $\gamma=3$, the speed of
this convergence depends on the deposition and erosion probabilities.

Figure 6: (a) The HVG degree distribution for increasing $n$, with $P=Q=0.25$,
and $\lambda=3$. Inset: Evolution of $\gamma(n)$. Symbols indicate numerical
results, the full line is an exponential fit. (b) The degree dependence of the
clustering coefficient $C(k)\sim k^{-\delta}$ for increasing $n$ with
$P=Q=0.3$, shifted vertically for visibility purposes. Results are averaged
over 5000 realisations.

The mean degree $\langle k\rangle$, degree variance $\sigma^{2}_{k}$, and mean
clustering coefficient $\langle C\rangle$ are shown in Fig. 7(b)-(d). It can
be seen that for $P=0$ and $P=1$, the mean degree is $\langle k\rangle=2$ and
the variance and mean clustering are $\sigma_{k}^{2}=\langle C\rangle=0$,
since the surface height profile is flat and the associated HVG becomes
regular. For these points, the degree distribution is trivially
$P(k)=\delta_{k,2}$ and hence the degree exponent $\gamma$ diverges to
infinity. For a scale-free network with exponent $\gamma$, the mean degree and
degree variance are

$\displaystyle\langle k\rangle$
$\displaystyle=\frac{\zeta(\gamma-1)-1}{\zeta(\gamma)-1}$ (22)
$\displaystyle\sigma_{k}^{2}$
$\displaystyle=\frac{\zeta(\gamma-2)\left[\zeta(\gamma)-1\right]-\zeta(\gamma-1)\left[\zeta(\gamma-1)-2\right]-\zeta(\gamma)}{\left(\zeta(\gamma)-1\right)^{2}}\,,$
(23)

respectively. For $\gamma=3$, this evaluates to $\langle k\rangle\approx 3.2$
and $\sigma_{k}^{2}\rightarrow\infty$. For the clustering coefficient $C$ no
such general analytical expression is known, since it depends on the specific
wiring on the level of vertex neighbourhoods, and hence on the mechanism that
formed the network. When we again consider the relationship between the local
clustering coefficients $C_{i}$, $i=1,2,\dots,\lambda^{n}$ and associated
vertex degrees, we find that for one degree $k$, multiple values of $C_{i}$
are possible. Hence, to study $C(k)$, we average the $C_{i}$ that correspond
to one value of $k$. One can find that $C$ exhibits a power-law relation,
i.e., $C(k)\sim k^{-\delta}$, with $\delta>0$. For increasing $n$, the
exponent becomes $\delta=1$. This is shown in Fig. 6(b). The hierarchical
clustering that was present for the HDDM persists even when randomness is
introduced.

It can be shown numerically (see appendix C), by means of a Floyd-Warshall
algorithm Cormen  _et al._ (2001), that the network diameter is again a power
law of the number of vertices, i.e., $D(N)\sim N^{\epsilon}$. This indicates
that the underlying height profile is self-similar and not self-affine, the
same result that was found in Section II for the HDDM. The scaling exponent
$\epsilon$ now depends on the system parameters but always takes a value
$0<\epsilon\leq 1$. The maximal value of $\epsilon=1$ is reached for flat
profiles where the diameter is maximal, since it is the distance between the
edge nodes of the system, i.e., $D(N)=N-1$.

From Fig 7, we find that for a quantity $X(P,Q)$, with
$X\in\left\\{\gamma,\langle k\rangle,\sigma^{2}_{k},\langle
C\rangle\right\\}$, the following invariance properties hold,

$\displaystyle X(1-a,a)=X(1-a,0)=X(0,a)$ (24)

with $0\leq a\leq 1$. Since this invariance holds for $\gamma$ (and trivially
for $\langle k\rangle$ and $\sigma_{k}^{2}$), we can conclude that the degree
distribution is identical for parameters choices that obey the invariance,
even for finite generation $n$.

Figure 7: (a) The degree exponent $\gamma$ as a function of $P$ for increasing
$Q$ for $n=8$. (b)-(d) Mean clustering coefficient $\langle C\rangle$, mean
degree $\langle k\rangle$ and variance $\sigma^{2}_{k}$ as a function of $P$
for increasing $Q$, with $n=5$. Results are averaged over 5000 realisations.

The different network quantities shown in Fig. 7 can be used to study the
structure of the underlying surface profile. Let us for example inspect the
mean degree $\langle k\rangle$ as a function of $S=P+Q$, i.e., the total
probability that a deposition or erosion event occurs at each site. Since $S$
is degenerate for different combinations of $P$ and $Q$, we average quantities
with $S$ fixed. One can immediately see (appendix C) that the highest value of
$\langle k\rangle$ occurs for $S\approx 0.6$. Similarly, the highest mean
clustering coefficient, highest degree variance, and lowest degree exponent
also occur for $S\approx 0.6$. Intuitively, one can expect that the mean
degree and clustering coefficient are maximal at the same value of $S$, since
the higher average convexity of the surface, characterised by a high value of
$\langle C\rangle$, implies that it is more irregular, which leads to a high
value of $\langle k\rangle$. Moreover, since the mean degree of a power-law
distribution only depends on the degree exponent $\gamma$, we can see that
when the mean degree grows closer to the asymptotic value $\langle
k\rangle\approx 3.2$ the degree exponent decreases to $\gamma\rightarrow 3$,
explaining the lower value of $\gamma$ for $S\approx 0.6$.

Before concluding, one remark is in order. While for $\lambda>2$ the
deposition history is unique, this is not the case anymore for $\lambda=2$.
Two different deposition histories can lead to the same landscape after two
generations by means of partial levelling of vertical segments Indekeu and
Fleerackers (1998), where subsequent deposition and erosion of neighbouring
blocks annihilates the vertical segment in between; the model thus becomes
non-retrodictive. This degeneracy can possibly impact the properties of the
associated HVG and its network topology. We defer to future research for a
more detailed study.

##  IV Conclusion and outlook

A mapping between the HDM for surface growth and complex networks was
introduced, based on the HVA. We have shown that the resulting network is
sparse and scale-free through the calculation of the degree distribution and
exponents, and by calculating the first two moments of this distribution. We
first performed these calculations for a deterministic version of the HDM and
subsequently extended them to random deposition with the addition of a erosion
probability. For the deterministic version, we then calculated the graph
diameter and derived an exact expression for the clustering coefficient, which
indicate that the network is fractal, self-similar, and exhibits hierarchical
clustering. The spectrum of the adjacency matrix was determined numerically
and it was shown that it converges to a yet undetermined fractal shape and its
spectral density possesses power-law tails. For the random deposition model,
we studied the resulting HVG numerically and showed analytically that
$\gamma\rightarrow 3$ asymptotically, indicating that the HRDM belongs to the
fractional Gaussian noise stochastic processes. We showed that the
hierarchical clustering and modularity found for the deterministic model is
robust when introducing randomness into the system, and that it remains
fractal through the scaling of the network diameter.

The connection between surface growth phenomena and complex networks is still
largely unexplored. The results from this work indicate that the structure of
the growing surface is encoded in the topology of the associated visibility
network, and that it can be discovered through careful analysis of different
network measures. Subsequent research could extend the above results to more
complex situations, such as the HRDM with lateral attachment, resembling a
hierarchical variation of the discrete model underlying the KPZ equation
(i.e., ballistic deposition). Such extensions present an interesting problem.
Since a point on the surface can possibly “see” both another point and an
accompanying overhang, the associated HVG can be a multigraph.

More general, one can construct a classification of “network universality
classes” for surface growth phenomena based on topological properties and
exponents determined from visibility graphs. This can complement or refine the
existing theory of universality classes for which the exponents are well-
known. In particular, future research could extend the VA formalism to growth
processes (or stochastic processes in general) that exhibit anomalous scaling
Mandelbrot  _et al._ (2002); Chen  _et al._ (2017).

Another subject of future study can include the search for the other scaling
exponents usually associated with dynamical scaling, i.e., the growth exponent
$\beta$ and dynamic exponent $z$, and the corresponding relations between
them. Since both exponents are concerned with the dynamical properties of the
surface, it may be possible to find explicit time-dependent network quantities
in the HVG and relate specific network growth exponents to $\beta$ and $z$.
Mapping the dynamical scaling relations back to these network exponents can
reveal hitherto unknown associations between them.

While the visibility algorithm approach to interface growth presents an
interesting avenue of research, classical methods to calculate, e.g., growth
exponents are still relatively easy to implement and do not usually
necessitate the use of other methods. Hence, more future research is needed to
identify the areas where VAs can really make a difference, either by
facilitating calculations or by decreasing computational load.

###### Acknowledgements.

The author is grateful for the spirited discussions with J. O. Indekeu and R.
Tielemans, whose insights were very enlightening, and for the suggestions
received from L. Lacasa. He also acknowledges the support of the G-Research
grant.

##  Appendix A Derivation of the individual node degree for the HDDM

Consider the adjacency matrix of the visibility graph in generation $n$, i.e.,
$A^{(n)}$ with dimensions $3^{n}\times 3^{n}$, see e.g., Fig. 8. There are
$3^{n}$ nodes with degree $k_{j}$, where $j=1,2,...,3^{n}$ denotes the node
index. Let us for convenience shift this index $j$ by one, $j^{\prime}=j-1$
and expand the result into base 3, writing the digits to a vector
$\vec{\alpha}_{j}$. For example, for the eighteenth node, $j=18$, hence
$j^{\prime}=17$ in generation $n=4$, and the number $j^{\prime}$ can be
written as $17=2\cdot 3^{0}+2\cdot 3^{1}+1\cdot 3^{2}+0\cdot 3^{3}$, hence,
$\vec{\alpha}_{18}=(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})=(2,2,1,0)$,
where every $\alpha_{i}\in\\{0,1,2\\}$ for $i=1,2,...,n$. Thus, every degree
$k_{j}$ acquires a unique vector $\vec{\alpha}_{j}$, which is associated with
its position in the adjacency matrix. We now consider the lowest-order
adjacency matrix $A^{(1)}$, which is given by

$A^{(1)}=\begin{pmatrix}0&1&0\\\ 1&0&1\\\ 0&1&0\end{pmatrix}\,.$ (25)

For every subsequent generation $n$, the adjacency matrix $A^{(1)}$ represents
a “unit cell”, which is repeated along the diagonal of every $A^{(n)}$. From a
physical standpoint, this can be attributed to the mutual visibility of
neighbors, where every node can always see its direct neighbors.

Figure 8: (a)-(c) Adjacency matrices for generations $n=2,3,4$, respectively.
Black squares indicate matrix elements $a_{ij}=1$, while white regions
indicate that $a_{ij}=0$. (d) A general depiction of the hierarchical
structure of the adjacency matrix $A^{(n+1)}$.

The first element in the vector $\vec{\alpha}$ indicates the column of the
node index within this unit cell, i.e., $\alpha_{1}=0,1,2$ indicates the left,
middle and right columns, respectively. The contribution of the position in
this unit cell to the total node degree is the sum of elements in the column
$\alpha_{1}$, and can be written consisely as $1+\delta_{\alpha_{1},1}$, where
$\delta_{ij}$ is the Kronecker delta, i.e., left and right columns contribute
$1$ to the degree, while the middle column contributes $2$, which can easily
be deduced from equation (25).

We now look at higher “levels” $k$ in the adjacency matrix through the
subsequent values in the $\vec{\alpha}$ vector; when $\alpha_{k}=0,1,2$, the
index is located in the left, middle or right column of the $k$th level of the
adjacency matrix. Let us look at the following example for $n=2$:

$\begin{split}A^{(2)}&=\left(\begin{array}[]{ccc}A^{(1)}&(C^{(1)})^{t}&0\\\
C^{(1)}&A^{(1)}&J^{(1)}C^{(1)}J^{(1)}\\\
0&J^{(1)}(C^{(1)})^{t}J^{(1)}&A^{(1)}\end{array}\right)\\\
&=\left(\begin{array}[]{ccccccccc}0&1&0&0&0&0&0&0&0\\\ 1&0&1&1&0&0&0&0&0\\\
0&1&0&1&0&0&0&0&0\\\ 0&1&1&0&1&0&0&0&0\\\ 0&0&0&1&0&1&0&0&0\\\
0&0&0&0&1&0&1&1&0\\\ 0&0&0&0&0&1&0&1&0\\\ 0&0&0&0&0&1&1&0&1\\\
0&0&0&0&0&0&0&1&0\end{array}\right)\,,\end{split}$ (26)

where the first level $(k=1)$ is made up of the $3\times 3$ matrices on the
main diagonal, and the second level $(k=2)$ is the full matrix. This simple
example of the second generation of the network already illustrates the nested
structure of the adjacency matrix.

Consider now for example the column with index $j=8$, which corresponds to
$j^{\prime}=j-1=7$, and which has the following representation in our base 3
notation: $\vec{\alpha}_{8}=(\alpha_{1},\alpha_{2})=(1,2)$. From this, one can
read off that the index $j=8$ is located in the right column of the $k=2$
level, which is the uppermost level, and in the middle column of the $k=1$
level, or unit cell. From the latter, the node degree $k_{8}$ gains a
contribution of $1+\delta_{\alpha_{1},1}=2$. The remaining factor of one
needed to obtain the exact result $k_{8}=3$ will be discussed now.

For levels $k\geq 2$, the block matrices located on the sub- and
superdiagonals are nonzero and are given by transposing the matrix $C^{(n)}$
either along the main diagonal, or along the antidiagonal. For $n=2$, the
block matrices are given by

$\begin{split}C^{(1)}&=\begin{pmatrix}0&1&1\\\ 0&0&0\\\ 0&0&0\\\
\end{pmatrix}\\\ (C^{(1)})^{t}&=\begin{pmatrix}0&0&0\\\ 1&0&0\\\ 1&0&0\\\
\end{pmatrix}\\\ J^{(1)}C^{(1)}J^{(1)}&=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
1&1&0\\\ \end{pmatrix}\\\ J^{(1)}(C^{(1)})^{t}J^{(1)}&=\begin{pmatrix}0&0&1\\\
0&0&1\\\ 0&0&0\\\ \end{pmatrix}\,,\end{split}$ (27)

where $J^{(n)}$ is the $3^{n}\times 3^{n}$ exchange matrix, i.e., for $n=1$,
this is

$J^{(1)}=\begin{pmatrix}0&0&1\\\ 0&1&0\\\ 1&0&0\\\ \end{pmatrix}\,.$ (28)

The first row in the $C^{(n)}$ matrix will be denoted as the vector
$\vec{v}^{(n)}=(\vec{0},\vec{v}^{(n-1)},\vec{v}^{(n-1)})$, where
$\vec{0}=(0,...,0)$ is the zero vector with $3^{n-1}$ elements. We list here
the first three $\vec{v}^{(n)}$:

$\begin{split}\vec{v}^{(0)}&=(1)\\\ \vec{v}^{(1)}&=(0,1,1)\\\
\vec{v}^{(2)}&=(0,0,0,0,1,1,0,1,1)\\\ \end{split}$ (29)

From the structure of (29), it can be seen that for $n\geq 2$ the first
$(3^{n}-1)/2$ elements in the vector $\vec{v}^{(n)}$ are identically zero. The
first two nonzero elements are then located at positions $j=(3^{n}+1)/2$ and
$j=(3^{n}+3)/2$. Subsequent nonzero elements are separated by gaps of unequal
size.

With moderate effort, one can deduce that subsequent nonzero elements are
located at indices that are shifted by a factor of
$(3i+\sum_{m=1}^{i}3^{a(m)})/2$ with respect to the positions $j=(3^{n}+1)/2$
and $j=(3^{n}+3)/2$. Herein, $a(m)$ is the so-called ruler function (OEIS
A001511), the elements of which are the exponent of the largest power of $2$
which divides a given number $2m$ Allouche and Shallit (2003). This sequence
can be characterized by the following recurrence relations:

$\displaystyle a(2m+1)$ $\displaystyle=1$ $\displaystyle a(2m)=1+a(m)\,.$ (30)

A plot of the magnitude of the elements of the ruler function $a(m)$ is shown
for $m=1,2,\dots,255$ in Fig. 9. One can now easily see why it is called as
such, since the rising and falling of the sequence as a function of the index
resembles the markings on a classical ruler.

Figure 9: Graphical representation of the ruler sequence $a(m)$ generated by
the recurrence relations (30) for $m=1,2,...255$. The pattern resembles the
markings on a classical ruler.

After some simple algebra, the elements $v_{j}^{(n)}$, with $j=1,...,3^{n}$
can be fully determined by

$v_{j}^{(n)}=\sum\limits_{i=0}^{2^{n-1}-1}\left(\delta_{j,\,\Xi_{ni}+\frac{1}{2}}+\delta_{j,\,\Xi_{ni}+\frac{3}{2}}\right)$
(31)

with

$\Xi_{ni}=\frac{1}{2}[3^{n}+3i+\sum\limits_{m=1}^{i}3^{a(m)}]\,.$ (32)

The number of nonzero elements in $v^{(n)}$ can be deduced by noticing that in
a generation $n+1$, this number has doubled with respect to generation $n$,
and the first generation possesses two nonzero elements. Hence, the total
number of nonzero elements is equal to $2^{n}$.

When the index $j$ is located in either the left or right column of the $k$th
level (i.e., when $\alpha_{k}=0$ or $\alpha_{k}=2$, respectively), it picks up
an extra factor of $+1$ for the node degree $k_{j}$ based on its position
within the vectors $\vec{v}^{(k-1)}$ for $\alpha_{k}=0$ or
$\vec{v}^{(k-1)}\cdot J^{(k-1)}$ for $\alpha_{k}=2$. Moreover, when $j$ is
located in the middle column $(\alpha_{k}=1)$ at positions $j=3^{k-1}+1$ or
$j=2\cdot 3^{k-1}$, the node degree $k_{j}$ picks up an extra $2^{k-1}$.

Summing over levels $k$ and shifting the indices for simplicity, we arrive
after some cumbersome algebra at the final result for the node degree
$k_{j}(n)$, yielding the expression in equation (5).

##  Appendix B Height distribution in the HRDM

We calculate the height probability density function (PDF)
$\mathcal{P}^{(n)}(H=h)$ or, alternatively, the cumulative height distribution
function (CDF) $\mathcal{F}^{n}(H\leq h)$ for the random hierarchical
deposition model. Let us consider the height as a geometric random walk and
let $\lambda^{-n}$ be the step size at time $n$, and $X_{n}$ a random variable
that can assume the values $-\lambda^{-n},0$ or $\lambda^{n}$, according to
whether the $n$th step is to the “left”, “middle” or “right”, which correspond
to particle deposition, no deposition, or erosion, respectively. The random
walk is then the countable sequence

$X_{1},\,X_{1}+X_{2},\,X_{1}+X_{2}+X_{3},\ldots$ (33)

and the height $H^{(n)}$ in generation $n\geq 1$ is consequently

$H^{(n)}=\sum_{j=1}^{n}X_{j}\,,$ (34)

which, for the geometric random walk with $\lambda\geq 2$ equals

$H^{(n)}=\sum_{j=1}^{n}a_{j}\lambda^{-j}\,,\,a_{j}\in\\{-1,0,1\\}$ (35)

The distribution $\mathcal{P}$ of $H$ is a measure whose cumulative
distribution function, defined as $\mathcal{F}(x)=\mathcal{P}([-1/2,x))$ is
continuous and non-decreasing. Note that the measure is singular with respect
to the Lebesgue measure, and hence we will from now on only work with the
cumulative distribution function $\mathcal{F}$. We show that $\mathcal{F}$ can
be found exactly for the hierarchical deposition model for
$n\rightarrow\infty$, albeit in an integral form that must be solved
numerically.

######  Theorem 1.

Consider the limiting ($n\rightarrow\infty$) hierarchical random deposition
process where $H=\sum_{j=1}^{\infty}X_{j}$ and $X_{j}$ are independent random
variables that have values $-\lambda^{-j}$, $0$ or $\lambda^{j}$ with
probabilities $Q$, $1-P-Q$ and $P$, respectively. The cumulative distribution
function $\mathcal{F}(h)=\mathcal{P}(H\leq h)$ is the following:

$\mathcal{F}(h)=\frac{1}{2\pi
i}\int_{\mathbb{R}}\mathrm{d}t\,\frac{\mathrm{e}^{it/(\lambda-1)}-\mathrm{e}^{-ith}}{t}\prod_{j=1}^{\infty}\left[P\mathrm{e}^{it\lambda^{-j}}+Q\mathrm{e}^{-it\lambda^{-j}}+(1-P-Q)\right]$
(36)

###### Proof.

The characteristic function $\varphi_{H}(t)$ for the height for
$n\rightarrow\infty$ can be written as

$\begin{split}\varphi_{H}(t)&=E\left[\mathrm{e}^{itH}\right]=E\left[\mathrm{e}^{it\sum_{j}X_{j}}\right]\\\
&=\prod_{j=1}^{\infty}\left[P\mathrm{e}^{it\lambda^{-j}}+Q\mathrm{e}^{-it\lambda^{-j}}+(1-P-Q)\right]\,,\end{split}$
(37)

where in the final equality we have used the fact that the $X_{j}$ are
independent random variables and the expected value can be factored.
Furthermore, we can define the characteristic function as
$\varphi_{H}(t)=\int_{\mathbb{R}}\mathrm{e}^{itx}\mathrm{d}\mathcal{F}(x)$. To
proceed, we define the indicator function $\chi_{h}(x)$ as

$\chi_{h}(x)=\begin{cases}1\,,&\qquad if\,x\in(-\frac{1}{\lambda-1},h]\\\
0\,,&\qquad else\,,\end{cases}$ (38)

which can also be defined through its Fourier transform $\hat{\chi}_{h}(t)$,
i.e.,

$\chi_{h}(x)=\int_{\mathbb{R}}\frac{\mathrm{e}^{itx}}{2\pi
it}\left(\mathrm{e}^{it/(\lambda-1)}-\mathrm{e}^{-ith}\right)\,\mathrm{d}t\,.$
(39)

The cumulative distribution function can now be written as

$\begin{split}\mathcal{F}(h)&=\int_{\mathbb{R}}\chi_{h}(x)\mathrm{d}\mathcal{F}(x)\\\
&=\frac{1}{2\pi
i}\int_{\mathbb{R}}\frac{\mathrm{e}^{it/(\lambda-1)}-\mathrm{e}^{-ith}}{t}\varphi_{H}(t)\mathrm{d}t\,,\end{split}$
(40)

which, upon inserting the expression for $\varphi_{H}(t)$ proves the theorem.
∎

The cumulative distribution function as defined by equation (36) cannot be
calculated analytically as far as we know, but numerical approximations by
truncating the infinite product can yield high-accuracy results. In Fig. 10,
we show the PDF and the CDF for the ternary Cantor set and a set with $P=1/3$
and $Q=1/4$, and compare results obtained by exact counting of all possible
heights with the infinite-generation approximation in equation (36). The
results coincide quite well.

Figure 10: (a)-(b) Height PDF and CDF, respectively, for $P=Q=1/2$, with the
ternary Cantor set as support. (c)-(d) Height PDF and CDF, respectively, for
$P=1/3$ and $Q=1/4$. Orange lines indicate exact distributions and black dots
are numerical approximations obtained from Eq. (36). Note that the same
pattern repeats at multiple scales, indicating the fractal behaviour of the
PDF. For both parameter choices, $n=5$ and $\lambda=3$ were used to find the
exact PDF/CDF, and the product in equation (36) is truncated to $j=100$.

Note that from the characteristic function all of the moments of
$\mathcal{P}(h)$ can be calculated by differentiation, i.e., $\langle
h^{m}\rangle=i^{-m}\varphi^{(m)}_{H}(0)$, where the superscript denotes
$m$-fold differentiation. After some simple algebra, the mean and variance are
then found to equal

$\displaystyle\langle h\rangle$ $\displaystyle=\frac{P-Q}{\lambda-1}$ (41)
$\displaystyle\sigma_{h}^{2}$ $\displaystyle=\langle h^{2}\rangle-\langle
h\rangle^{2}=\frac{(P+Q)-(P-Q)^{2}}{\lambda^{2}-1}\,.$ (42)

For $\lambda=3$, $P=Q=1/2$, the mean and variance equal $\langle h\rangle=0$
and $\sigma_{h}^{2}=1/8$, respectively. For $P=1/3$ and $Q=1/4$, they equal
$\langle h\rangle=1/24$ and $\sigma_{h}^{2}=83/1152\approx 0.072$,
respectively. This is confirmed by numerical calculation.

##  Appendix C Numerical degree exponents, moments, diameter and clustering
for the HRDM

To numerically find the degree exponent $\gamma$ of the scale-free networks
generated by the HVA, we solve the transcendental equation for the maximum
likelihood estimator (MLE) $\hat{\gamma}$ Clauset  _et al._ (2009):

$\frac{\zeta^{\prime}(\hat{\gamma},k_{\text{min}})}{\zeta(\hat{\gamma},k_{\text{min}})}=-\frac{1}{m}\sum\limits_{i=1}^{m}\ln{k_{i}}\,,$
(43)

where $\zeta$ is the Hurwitz zeta function Olver  _et al._ (2010), $m$ is the
total number of observed values of the degrees $k$, and $k_{\text{min}}$ is
the cutoff degree above which the distribution is a discrete power law. The
prime denotes differentiation with respect to the first argument. The standard
error $\sigma$ associated with the MLE is

$\sigma=\left[m\left(\frac{\zeta^{\prime\prime}(\hat{\gamma},k_{\text{min}})}{\zeta(\hat{\gamma},k_{\text{min}})}\right)-\left(\frac{\zeta^{\prime}(\hat{\gamma},k_{\text{min}})}{\zeta(\hat{\gamma},k_{\text{min}})}\right)^{2}\right]^{-\frac{1}{2}}\,.$
(44)

Alternatively, the exponent can be found approximately, whereby the true
power-law distributed integers are approximated as continuous reals rounded to
the nearest integer, i.e.,

$\begin{split}\hat{\gamma}^{\prime}&\simeq
1+m\left[\sum\limits_{i=1}^{m}\ln\frac{k_{i}}{k_{\text{min}}-\frac{1}{2}}\right]^{-1}\,.\end{split}$
(45)

The error can be found in the same manner as before. We complement our
analysis with a Kolmogorov-Smirnov (KS) test Clauset  _et al._ (2009) on the
cumulative degree distribution to find the optimal $\hat{\gamma}$ (or
$\hat{\gamma}^{\prime}$) and $k_{\text{min}}$ that fit the simulation data. We
show the MLE results for the estimators $\hat{\gamma}$ and
$\hat{\gamma}^{\prime}$ in Table 2.

Table 2: Evolution of $\hat{\gamma}$ and $\hat{\gamma}^{\prime}$ exponents with associated standard errors $\sigma$ and $\sigma^{\prime}$. Parameters are $P=Q=0.25$ and $\lambda=3$ and results are averaged over 5000 realizations. n | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---
$\hat{\gamma}$ | 6.476 | 5.809 | 4.888 | 4.437 | 4.045 | 3.810
$\sigma$ | 0.066 | 0.040 | 0.016 | 0.013 | 0.008 | 0.012
$\hat{\gamma}^{\prime}$ | 7.379 | 6.179 | 5.081 | 4.507 | 4.083 | 3.830
$\sigma^{\prime}$ | 0.017 | 0.008 | 0.004 | 0.004 | 0.003 | 0.006

As a function of the total event probability $S=P+Q$, the mean clustering
coefficient $\langle C\rangle$, the mean degree $\langle k\rangle$, the degree
variance $\sigma_{k}^{2}$ and the degree exponent $\gamma$ are shown in Fig.
11. The lines connecting data points are only for visualisation purposes.

Figure 11: (a) - (c) Mean clustering coefficient, mean degree and degree
variance as a function of $S=P+Q$ for $n=5$. (d) Degree exponent $\gamma$ as a
function of $S=P+Q$ for $n=8$. Results are averaged over 5000 realisations.

The network diameter $D$ is shown in Fig. 12(a) as a function of the number of
vertices $N=3^{n}$ for fixed deposition probability $P=0.2$. Full lines are
power-law fits to the numerical data. The power-law behaviour of $D(N)\sim
N^{\epsilon}$ is clear from Fig. 12(b). Lines are only shown for visualisation
purposes.

Figure 12: (a) The HVG diameter $D$ as a function of $N=3^{n}$ for fixed
$P=0.2$ and $0\leq Q\leq P$ in steps of $\Delta Q=0.1$, and (b) the scaling
exponent $\epsilon$ as a function $P$ for fixed $n=6$ and $0\leq Q\leq 1$ n
steps of $\Delta Q=0.1$. Results are averaged over 5000 realisations.

## References

  * Barabási and Stanley (1995) A.-L. Barabási and H. E. Stanley,  _Fractal Concepts in Surface Growth_ (Cambridge University Press, 1995).
  * Family and Vicsek (1985) F. Family and T. Vicsek, J. Phys. A: Math. Gen. 18, L75 (1985).
  * Family (1986) F. Family, J. Phys. A: Math. Gen. 19, L441 (1986).
  * Meakin  _et al._ (1986) P. Meakin, P. Ramanlal, L. M. Sander, and R. C. Ball, Phys. Rev. A 34, 5091 (1986).
  * Edwards and Wilkinson (1982) S. F. Edwards and D. R. Wilkinson, Proc. R. Soc. A 381, 17 (1982).
  * Vvedensky  _et al._ (1993) D. D. Vvedensky, A. Zangwill, C. N. Luse, and M. R. Wilby, Phys. Rev. E 48, 852 (1993).
  * Buceta  _et al._ (2014) R. C. Buceta, D. Hansmann, and B. von Haeften, J. Stat. Mech.: Theory Exp. 2014 (12), P12028.
  * Kardar  _et al._ (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. 56, 889 (1986).
  * Bertini and Giacomin (1997) L. Bertini and G. Giacomin, Commun. Math. Phys. 183, 571 (1997).
  * Gomes-Filho  _et al._ (2021) M. S. Gomes-Filho, A. L. Penna, and F. A. Oliveira, Results Phys. 26, 104435 (2021).
  * Indekeu and Fleerackers (1998) J. O. Indekeu and G. Fleerackers, Physica A 261, 294 (1998).
  * Cheng  _et al._ (1987) Z. Cheng, L. Jacobs, D. Kessler, and R. Savit, J. Phys. A: Math. Gen. 20, L1095 (1987).
  * Baiod  _et al._ (1988) R. Baiod, D. Kessler, P. Ramanlal, L. Sander, and R. Savit, Phys. Rev. A 38, 3672 (1988).
  * Krug and Meakin (1991) J. Krug and P. Meakin, Phys. Rev. A 43, 900 (1991).
  * Zhang (1990) Y.-C. Zhang, Physica A 170, 1 (1990).
  * Hosseinabadi and Masoudi (2019) S. Hosseinabadi and A. A. Masoudi, Phys. Rev. E 99, 012130 (2019).
  * Pioli  _et al._ (2019) L. Pioli, C. Bonadonna, and M. Pistolesi, Sci. Rep. 9, 10006 (2019).
  * Brykina and Egorova (2021) I. G. Brykina and L. A. Egorova, Adv. Astron. 2021, 9914717 (2021).
  * Finn  _et al._ (1982) J. E. Finn, S. Agarwal, A. Bujak, J. Chuang, L. J. Gutay, A. S. Hirsch, R. W. Minich, N. T. Porile, R. P. Scharenberg, B. C. Stringfellow, and F. Turkot, Phys. Rev. Lett. 49, 1321 (1982).
  * Turcotte (1986) D. L. Turcotte, Geophys. Res. Solid Earth 91, 1921 (1986).
  * Posazhennikova and Indekeu (2000) A. Posazhennikova and J. O. Indekeu, Mod. Phys. Lett. B 14, 119 (2000).
  * Posazhennikova and Indekeu (2001) A. Posazhennikova and J. O. Indekeu, Int. J. Thermophys. 22, 1123 (2001).
  * Indekeu  _et al._ (2000) J. Indekeu, G. Fleerackers, A. Posazhennikova, and E. Bervoets, Physica A 285, 135 (2000).
  * Posazhennikova and Indekeu (2014) A. Posazhennikova and J. O. Indekeu, Physica A 414, 240 (2014).
  * Indekeu and Giuraniuc (2004) J. Indekeu and C. Giuraniuc, Physica A 336, 14 (2004).
  * Indekeu and Sznajd-Weron (2003) J. O. Indekeu and K. Sznajd-Weron, Phys. Rev. E 68, 061904 (2003).
  * Berx  _et al._ (2021) J. Berx, E. Bervoets, C. V. Giuraniuc, and J. O. Indekeu, Physica A 574, 125998 (2021).
  * Lacasa  _et al._ (2008) L. Lacasa, B. Luque, F. Ballesteros, J. Luque, and J. C. Nuño, PNAS 105, 4972 (2008).
  * Kaki  _et al._ (2022) B. Kaki, N. Farhang, and H. Safari, Sci. Rep. 12, 16835 (2022).
  * Yu  _et al._ (2012) Z. G. Yu, V. Anh, R. Eastes, and D.-L. Wang, Nonlinear Process. Geophys. 19, 657 (2012).
  * Ahmadlou  _et al._ (2010) M. Ahmadlou, H. Adeli, and A. Adeli, J. Neural Transm. 117, 1099 (2010).
  * Juniper and Sujith (2018) M. Juniper and R. Sujith, Annu. Rev. Fluid Mech. 50, 661 (2018).
  * Paiva  _et al._ (2022) L. R. Paiva, S. G. Alves, L. Lacasa, O. DeSouza, and O. Miramontes, arXiv:2204.11690 (2022).
  * Lacasa  _et al._ (2012) L. Lacasa, A. Nuñez, É. Roldán, J. M. R. Parrondo, and B. Luque, EPJ B 85, 217 (2012).
  * Wang  _et al._ (2018) M. Wang, A. L. M. Vilela, R. Du, L. Zhao, G. Dong, L. Tian, and H. E. Stanley, Sci. Rep. 8, 5130 (2018).
  * Lacasa  _et al._ (2015) L. Lacasa, V. Nicosia, and V. Latora, Sci. Rep. 5, 15508 (2015).
  * Nuñez  _et al._ (2012) A. M. Nuñez, L. Lacasa, J. P. Gomez, and B. Luque, in  _New Frontiers in Graph Theory_, edited by Y. Zhang (IntechOpen, Rijeka, 2012) Chap. 6.
  * Luque  _et al._ (2009) B. Luque, L. Lacasa, F. Ballesteros, and J. Luque, Phys. Rev. E 80, 046103 (2009).
  * Brú  _et al._ (2014) A. Brú, E. Alós, J. C. Nuño, and M. F. de Dios, Sci. Rep. 4, 7550 (2014).
  * Kartha (2017) M. J. Kartha, Phys. Lett. A 381, 556 (2017).
  * Allouche and Shallit (2003) J.-P. Allouche and J. Shallit,  _Automatic Sequences: Theory, Applications, Generalizations_ (Cambridge University Press, 2003).
  * Barabási  _et al._ (2001) A.-L. Barabási, E. Ravasz, and T. Vicsek, Physica A 299, 559 (2001).
  * Iguchi and Yamada (2005) K. Iguchi and H. Yamada, Phys. Rev. E 71, 036144 (2005).
  * Song  _et al._ (2006) C. Song, S. Havlin, and H. A. Makse, Nat. Phys. 2, 275 (2006).
  * Ravasz  _et al._ (2002) E. Ravasz, A. L. Somera, D. A. Mongru, Z. N. Oltvai, and A.-L. Barabási, Science 297, 1551 (2002).
  * Dorogovtsev  _et al._ (2002) S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. E 65, 066122 (2002).
  * Donner and Donges (2012) R. V. Donner and J. F. Donges, Acta Geophys. 60, 589 (2012).
  * Flanagan  _et al._ (2019) R. Flanagan, L. Lacasa, and V. Nicosia, J. Phys. A: Math. Theor. 53, 025702 (2019).
  * Goh  _et al._ (2001) K.-I. Goh, B. Kahng, and D. Kim, Phys. Rev. E 64, 051903 (2001).
  * Kim and Wilhelm (2008) J. Kim and T. Wilhelm, Physica A 387, 2637 (2008).
  * Fioriti  _et al._ (2012) V. Fioriti, A. Tofani, and A. D. Pietro, Complex Syst. 21 (2012).
  * Lacasa (2016) L. Lacasa, J. Phys. A: Math. Theor. 49, 35LT01 (2016).
  * de la Torre  _et al._ (2000) A. C. de la Torre, A. Maltz, H. O. Mártin, P. Catuogno, and I. García-Mata, Phys. Rev. E 62, 7748 (2000).
  * Lacasa  _et al._ (2009) L. Lacasa, B. Luque, J. Luque, and J. C. Nuño, EPL 86, 30001 (2009).
  * Gutin  _et al._ (2011) G. Gutin, T. Mansour, and S. Severini, Physica A 390, 2421 (2011).
  * Das Sarma  _et al._ (1996) S. Das Sarma, C. J. Lanczycki, R. Kotlyar, and S. V. Ghaisas, Phys. Rev. E 53, 359 (1996).
  * López  _et al._ (1997) J. M. López, M. A. Rodríguez, and R. Cuerno, Phys. Rev. E 56, 3993 (1997).
  * Giuraniuc (2006) C. V. Giuraniuc,  _Deposition and network models in statistical and condensed matter physics with interdisciplinary applications_ , Ph.D. thesis, KU Leuven (2006).
  * Clauset  _et al._ (2009) A. Clauset, C. R. Shalizi, and M. E. J. Newman, SIAM Review 51, 661 (2009).
  * Olver  _et al._ (2010) F. Olver, D. Lozier, R. Boisvert, and C. Clark,  _The NIST Handbook of Mathematical Functions_ (Cambridge University Press, New York, NY, 2010).
  * Cormen  _et al._ (2001) T. Cormen, C. Leiserson, M. Press, R. Rivest, C. Stein, and M.-H. P. Company,  _Introduction To Algorithms_ , Introduction to Algorithms (MIT Press, 2001).
  * Mandelbrot  _et al._ (2002) B. Mandelbrot, F. Damerau, M. Frame, K. McCamy, J. Van Ness, and J. Wallis,  _Gaussian Self-Affinity and Fractals: Globality, The Earth, 1/f Noise, and R/S_, Selecta (Old or New), Volume H (Springer, 2002).
  * Chen  _et al._ (2017) L. Chen, K. E. Bassler, J. L. McCauley, and G. H. Gunaratne, Phys. Rev. E 95, 042141 (2017).