11institutetext: Ana Paula Millán 22institutetext: Amsterdam UMC, Vrije
Universiteit Amsterdam, Department of Clinical Neurophysiology and MEG Center,
Amsterdam Neuroscience, De Boelelaan 1117, Amsterdam, The Netherlands.
22email<EMAIL_ADDRESS>33institutetext: Juan G. Restrepo
44institutetext: Department of Applied Mathematics, University of Colorado at
Boulder, Boulder, CO 80309, USA44email<EMAIL_ADDRESS>Joaquín J. Torres 66institutetext: Departamento de Electromagnetismo y Física
de la Materia and Instituto Carlos I de Física Teórica y Computacional,
Universidad de Granada, 18071, Granada, Spain 66email:
<EMAIL_ADDRESS>Ginestra Bianconi 88institutetext:
School of Mathematical Sciences, Queen Mary University of
London99institutetext: The Alan Turing Institute, The British Library, London,
UK 99email<EMAIL_ADDRESS>
# Geometry, Topology and Simplicial Synchronization

Ana Paula Millán  Juan G. Restrepo  Joaquín J. Torres and Ginestra Bianconi

###### Abstract

Simplicial synchronization reveals the role that topology and geometry have in
determining the dynamical properties of simplicial complexes. Simplicial
network geometry and topology are naturally encoded in the spectral properties
of the graph Laplacian and of the higher-order Laplacians of simplicial
complexes. Here we show how the geometry of simplicial complexes induces
spectral dimensions of the simplicial complex Laplacians that are responsible
for changing the phase diagram of the Kuramoto model. In particular,
simplicial complexes displaying a non-trivial simplicial network geometry
cannot sustain a synchronized state in the infinite network limit if their
spectral dimension is smaller or equal to four. This theoretical result is
here verified on the Network Geometry with Flavor simplicial complex
generative model displaying emergent hyperbolic geometry. On its turn
simplicial topology is shown to determine the dynamical properties of the
higher-order Kuramoto model. The higher-order Kuramoto model describes
synchronization of topological signals, i.e., phases not only associated with
the nodes of a simplicial complexes but associated also to higher-order
simplices, including links, triangles and so on. This model displays
discontinuous synchronization transitions when topological signals of
different dimension and/or their solenoidal and irrotational projections are
coupled in an adaptive way.

[scale=.25]figure1.pdf

Figure 1:  Simplicial complexes encode the rich simplicial network topology
and geometry of data and models, which strongly affects the higher-order
dynamics. In this chapter we will see how this interplay between structure and
dynamics can enrich our understanding of synchronization dynamics defined on
simplicial complexes.

##  1 Introduction

The interplay between structure and dynamics of complex networks
barabasi2016network ; Doro_crit ; Barrat  has been at the forefront of network
theory since the beginning of the field. In this context it has been found
that the combinatorial and statistical properties of complex networks have
unexpected effects on dynamics. For instance, a scale-free degree distribution
changes the phase diagram of a wide range of dynamical processes including
percolation, epidemic spreading, and the Ising model. The recent interest on
higher-order networks battiston2020networks ; bianconi2021higher ;
bianconi2015interdisciplinary ; eliassirad  provides an opportunity to bring a
fresh perspective to this subject. Indeed, higher-order networks, and in
particular simplicial complexes, constitute the ideal mathematical framework
to capture the simplicial network topology and geometry of data. Here we
reveal that the network topology and geometry of simplicial complexes can be
crucial to define higher-order dynamics. The dynamical process considered in
this chapter is synchronization strogatz2000kuramoto ;
arenas2008synchronization ; boccaletti2018synchronization , captured by the
Kuramoto model kuramoto1975  and the recently introduced higher-order Kuramoto
model millan2020explosive . The dynamical properties of these dynamical
processes will be shown to be highly dependent on the spectral properties
torres2020simplicial ; Barbarossa  of the simplicial complexes
courtney2016generalized ; bianconi2016network ; bianconi2017emergent ;
mulder2018network ; courtney2017weighted  on which they are defined. The main
message of this chapter is summarized in Figure 1, which highlights the role
of network topology and network geometry in shaping higher-order network
dynamics. In particular, in this chapter we will disclose how the spectral
properties of simplicial complexes are foundational to reveal the relation
between higher-order network geometry, topology and dynamics. Note that while
our approach to simplicial synchronization is based on simplicial network
geometry and topology, other approaches based on a combinatorial definition of
higher-order interactions have been pursued in the literature
skardal2019abrupt ; skardal2019higher ; skardal2020memory , as covered by the
Skardal and Arenas chapter of this book.

[scale=.20]figure2.pdf

Figure 2:  Schematic representation of the two classes of simplicial complex
models considered in this work. The Network Geometry with Flavor is a model of
growing simplicial complexes describing emergent hyperbolic network
geometries. The left panel shows a realization of the NGF simplicial complex
of dimension $d=2$ and flavor $s=-1$. The configuration model of simplicial
complexes (right panel) is a maximum entropy model enforcing a given sequence
of generalized degrees of the nodes. Right panel reprinted with permission
from Ref. courtney2016generalized . ©Copyright (2016) by the American Physical
Society.

##  2 Simplicial complex models

Simplicial network models are ideal to test, in a well-controlled setting, the
interplay between simplicial network geometry, topology and dynamics. Here we
focus in particular on two large classes of simplicial complex models with
very distinct structural properties: the Network Geometry with Flavor (NGF)
bianconi2016network ; bianconi2017emergent ; mulder2018network ;
courtney2017weighted  and the configuration model of simplical complexes
courtney2016generalized  (see schematic illustrations of the two models in
Figure 2). These models are implemented in codes available at the repository
GitHubrepository .

###  2.1 The Network Geometry with Flavor (NGF)

The Network Geometry with Flavor (NGF) bianconi2016network ;
bianconi2017emergent ; mulder2018network ; courtney2017weighted  is a general
mathematical framework for growing simplicial complexes that displays emergent
hyperbolic network geometry. In other words, the NGF model generates
simplicial complexes with hyperbolic geometry that evolve following purely
combinatorial and stochastic rules that do not make any use of the natural
hyperbolic embedding of the simplicial complexes.

The NGFs are simplicial complexes characterized by two main parameters: the
dimension of the simplicial complex $d$ and the flavor $s$ which is a
parameter that takes values $s\in\\{-1,0,1\\}$. The NGFs are generated by a
dynamical process which starting at time $t=1$ from a single $d$-dimensional
simplex proceeds at each time $t>1$ by adding a new $d$-dimensional simplex to
the simplicial complex. The new $d$-dimensional simplex includes one new node
and is glued to a $(d-1)$-dimensional face $\alpha$ of the existing simplicial
complex chosen according to the probability

$\displaystyle\Pi_{\alpha}=\frac{1-s+sk_{d,d-1}(\alpha)}{\sum_{\alpha^{\prime}}1-s+sk_{d,d-1}(\alpha^{\prime})},$
(1)

where $k_{d,d-1}(\alpha)$ indicates the generalized degree
courtney2016generalized  of a $(d-1)$-dimensional face $\alpha$, i.e., the
number of $d$-dimensional simplices incident to the $(d-1)$-dimensional face
$\alpha$. This model generates emergent hyperbolic simplicial complexes which
satisfy Gromov criteria Gromov  of hyperbolicity and are $\delta$-hyperbolic
with $\delta=1$ for every value of the flavor $s$ millan2021local . Moreover,
for flavor $s=-1$ the generated simplicial complexes form $d$-dimensional
hyperbolic manifolds bianconi2017emergent . The network skeleton of the NGFs
are small world, display hierarchical community structure and are scale-free
for $d\geq 2-s$ bianconi2016network ; bianconi2017emergent ; mulder2018network
. Interestingly, in the case $d=1$ and $s=1$ the NGF reduces to the Barabási-
Albert model, and for $d=1$ and $s=-1$ the NGF reduces to random Apollonian
networks.

This model can be generalized in different directions. Instead of considering
simplicial complexes, one can use a similar model to generate cell complexes
by gluing together convex regular polytopes mulder2018network . Another
possibility is to consider weighted simplicial complexes or to allow any new
node to be incident to more than one $d$-dimensional simplex
courtney2017weighted . Finally, the faces can be assigned a fitness that can
be used to modulate the attachment probability $\Pi_{\alpha}$ causing
topological phase transitions for certain fitness distributions
bianconi2016network ; bianconi2017emergent .

###  2.2 Configuration model of simplicial complexes

The configuration model of simplicial complexes courtney2016generalized  is a
maximum entropy model of pure $d$-dimensional simplicial complexes. A pure
$d$-dimensional simplicial complex $\mathcal{K}$ can be fully encoded in a
$(d+1)$-dimensional adjacency tensor indicating the presence of each
$d$-dimensional facet of the simplicial complex. In particular adjacency
tensor ${\bf a}$ has elements $a_{\alpha}=1$ if the $d$-dimensional simplex
$\alpha$ is present in the simplicial complex, otherwise $a_{\alpha}=0$.

The configuration model of simplicial complexes courtney2016generalized  is
the least biased ensemble of simplicial complexes with a given generalized
degree sequence of the nodes

$\displaystyle{\bf k}_{d,0}=\\{k_{d,0}(1),k_{d,0}(2),\ldots k_{d,0}(N)\\}$ (2)

where $k_{d,0}(i)$ indicates the generalized degree of the node $i$, i.e., the
number of $d$-dimensional complexes of the node $i\in\\{1,2,\ldots,N\\}$.

The configuration model of simplicial complexes is fully characterized by the
probability $P(\mathcal{K})$ assigned to each pure $d$-dimensional simplicial
complex $\mathcal{K}$ of $N$ nodes. The probability $P(\mathcal{K})$ maximizes
the entropy $S$ of the simplicial complex ensemble

$\displaystyle S=-\sum_{\mathcal{K}}P(\mathcal{K})\ln P(\mathcal{K}),$ (3)

given the constrain that each node $i$ has generalized degree $k_{d,0}(i)$,
i.e.,

$\displaystyle\sum_{\alpha\in Q_{d}(N)|i\subset\alpha}a_{\alpha}=k_{d,0}(i),$
(4)

where $Q_{d}(N)$ indicates the set of all possible $d$-dimensional simplices
of a simplicial complex formed by $N$ nodes. Therefore the configuration model
of simplicial complexes is characterized by the uniform distribution

$\displaystyle
P(\mathcal{K})=\frac{1}{\mathcal{N}}\prod_{i=1}^{N}\delta\left(\sum_{\alpha\in
Q_{d}(N)|i\subset\alpha}a_{\alpha},k_{d,0}(i)\right),$ (5)

where here and in the following $\delta(a,b)$ indicates the Kronecker delta.
In Ref. courtney2016generalized  the Authors proposed an algorithm for
sampling simplicial complexes from this ensemble. This algorithm
GitHubrepository  uses the mapping of simplicial complexes to factor graphs.
The configuration model of simplicial complexes is a very valuable null model
of simplicial complexes, and provides an ideal benchmark to study dynamical
processes on higher-order networks.

##  3 Laplacians

###  3.1 Graph Laplacian

The graph Laplacian describes linear diffusion on a network and it is an
important operator that encodes the interplay between network structure and
dynamics Barrat ; burioni2005random . The graph Laplacian can also capture the
underlying network geometry of the skeleton of a simplicial complex, i.e., of
the network obtained from the simplicial complex by retaining only its nodes
and links.

The graph Laplacian $\bm{L}_{[0]}$ is the discrete version of the Laplacian
operator defined on continuous functions. The graph Laplacian of a network
with $N_{[0]}$ nodes is a $N_{[0]}\times N_{[0]}$ matrix of elements

$\left[L_{[0]}\right]_{ij}=k_{i}\delta({i,j})-a_{ij},$ (6)

where $k_{i}$ is the degree of node $i$ and $a_{ij}$ are the elements of the
adjacency matrix of the network. In some cases it is useful to consider a
generalization of this operator called the normalized Laplacian
$\bm{\hat{L}}_{[0]}$ that has elements

$\left[\hat{L}_{[0]}\right]_{ij}=\delta({i,j})-\frac{a_{ij}}{k_{i}}.$ (7)

For instance, the normalized Laplacian is commonly employed to describe random
walk dynamics on a network.

Both the standard and the normalized Laplacians have real eigenvalues
$0=\lambda_{1}\leq\lambda_{2}\leq\ ...\ \leq\lambda_{N_{[0]}}$. The density of
eigenvalues is described by the spectral density,

$\rho(\lambda)=\frac{1}{N_{[0]}}\sum_{i=1}^{N_{[0]}}\tilde{\delta}(\lambda-\lambda_{i}),$
(8)

where $\tilde{\delta}\left(x\right)$ indicates the delta function. The
eigenvectors of the normalized Laplacian also encode relevant properties of
the underlying network.

###  3.2 Spectral dimension

For many complex networks, the smallest non-zero eigenvalue (also called
Fiedler eigenvalue) $\lambda_{2}$ remains finite as the network size
increases. In this case, the network is said to have a _spectral gap_. On the
contrary, if $\lambda_{2}\to 0$ as $N\to\infty$, and the density of
eigenvalues $\rho(\lambda)$ scales as

$\rho(\lambda)\propto\lambda^{d_{S}^{[0]}/2-1},$ (9)

for $\lambda\ll 1$, the network is said to present a _spectral dimension_
$d_{S}^{[0]}$ rammal1983random ; burioni2005random ; burioni1996universal .
The spectral dimension can be interpreted as the perceived dimension of the
network by diffusion processes, and it is a notable feature of networks with
an underlying geometrical nature. This is a definition of dimension that is
alternative to the Hausdorff dimension $d_{H}$ characterizing the scaling of
the diameter $D$ of a network with the network size $N$, i.e., $D\propto
N^{1/d_{H}}$.

For Euclidean lattices of dimension $d$, the spectral dimension coincides with
the Hausdorff dimension and we have $d_{S}^{[0]}=d_{H}=d$. However, in general
the spectral dimension of a network skeleton of a simplicial complex does not
have to coincide with the topological dimension of the simplicial complex, nor
with its Hausdorff dimension wu2015emergent ; rammal1983random ;
millan2020complex ; millan2021local .

While Euclidean lattices display a finite spectral dimension, random graphs
and the configuration model of networks are instead characterized by a finite
spectral gap. The presence of a finite spectral gap indicates the mean-field
nature of the network interactions, and the absence of a clear notion of
locality for these networks. Interestingly, also the configuration model of
simplicial complexes displays a finite spectral gap.

Remarkably, the emergent network geometry of NGF reveals itself on their
significant spectral properties. Indeed the NGF network skeletons, together
with other small-world of simplicial complexes tadic1 , have a finite spectral
dimension whose value can be tuned according to the different control
parameters mulder2018network  although the NGFs are small world, i.e., they
have an infinite Hausdorff dimension $d_{H}=\infty$.

Figure 3: Spectral dimension of NGF networks. Panel (a) The cumulative
distribution $\rho_{c}(\lambda)$ of eigenvalues of the NGF with flavor $s=-1$
is shown for dimension $d=2,3,4,5$. The power-law scaling of
$\rho_{c}(\lambda)$ observed for small values of $\lambda$ indicates that the
skeleton of NGF has a finite spectral dimension. Panels (b-d) The fitted
spectral dimension $d_{S}^{[0]}$ of the skeleton of NGF simplicial and cell
complexes being formed by simplices (panel (b)), hypercubes (panel (c)) and
orthoplexes (panel (c)), is shown for values of the flavor $s\in\\{-1,0,1\\}$
as indicated in the legend of panel (d). Lines indicate best fit of the
$d_{S}^{[0]}$ versus $d$ dependence using parabolic (panel (b)) and
exponential (panels (c), (d)) functional forms. Data from Ref.
millan2019synchronization  and Ref. mulder2018network . Details of the fits
are described in Ref. mulder2018network .

.

For $s=-1$, NGF networks formed purely by $d$-dimensional simplices have
$d_{S}^{[0]}\sim d$ for $d\in\\{2,3,4\\}$ as shown in Figure 3
millan2018complex . More generally, for NGFs formed by regular polytopes,
$d_{S}^{[0]}$ increases with the dimension $d$ of the polytopes, and it
saturates for hypercubes and orthoplexes ($d_{S}^{[0]}\leq 3$)
millan2019synchronization . It was shown in Ref. mulder2018network  that a
similar trend is observed for different flavours of the NGF networks:
$d_{S}^{[0]}$ grows proportional to $d^{2}$ for simplicial NGF networks,
whereas it saturates at a value $d_{S}^{[0]}=\bar{d}_{S}$ with
$2\leq\bar{d}_{S}\leq 3$ for hypercube and orthoplex NGF networks. Moreover,
it was shown recently in Ref. millan2021local  that, in a generalization of
the NGF model in which different polytopes are glued together in the same
higher-order network, the spectral dimension of the network skeleton can be
continuously tuned as a function of the fraction of simplexes in the cell
complex.

Thus, not only the dimension of the building blocks shapes the spectral
dimension of the networks, but the specific nature and symmetry of these
building blocks also play a role in the emerging spectral dimension of the
network skeleton (see Figure 3).

###  3.3 Higher-order Laplacians

Important topological aspects of simplicial complexes are reflected in the
spectral properties of the higher-order Laplacians that generalize the graph
Laplacian to describe diffusion that occurs among higher-order simplices. The
graph Laplacian describes diffusion occurring among nodes connected by links.
Similarly the $n$-th order up-Laplacian describes the diffusion occurring
among $n$-dimensional simplices connected by $(n+1)$-dimensional simplices and
the $n$-th order down-Laplacian describes the diffusion occurring among
$n$-dimensional simplices connected by $(n-1)$-dimensional simplices. The
higher-order Laplacians capture the topology of simplicial complexes. For
instance their spectrum encodes the Betti numbers, i.e., the number of
$n$-dimensional cavities of the simplicial complex.

The higher-order Laplacians are defined in terms of the incidence matrices of
the simplicial complex which represent the boundary operators playing a
fundamental role in algebraic topology.

Here we provide a brief introduction to the necessary elements of algebraic
topology needed to define higher-order Laplacians.

We consider a $d$-dimensional simplicial complex formed by $N_{[n]}$ simplices
of dimension $n$. The simplices of the simplicial complexes are associated
with an orientation induced by the labelling of the nodes so that the link
$[i,j]$ has a positive orientation if $j>i$ and so on (see Figure 4).

We consider algebraic entities called $n$-chains that are linear combinations
of $n$-dimensional simplices with coefficients in $\mathbb{Z}$. In a less
informal definition $n$-chains are the elements of a free abelian group
${\mathcal{C}}_{n}$ with basis on the $n$-simplices of the simplicial complex.
The boundary map is a linear map
$\partial_{n}:{\mathcal{C}}_{n}\to{\mathcal{C}}_{n-1}$ defined by its action
on each simplex $\alpha=[i_{0},i_{1},i_{2}\ldots,i_{n}]$. In particular the
boundary map associates to every $n$-dimensional simplex
$\alpha=[i_{0},i_{1},i_{2}\ldots,i_{n}]$ a linear combination of the
$(n-1)$-dimensional oriented faces at its boundary, given by

$\displaystyle\partial_{n}[i_{0},i_{1}\ldots,i_{n}]=\sum_{p=0}^{n}(-1)^{p}[i_{0},i_{1},\dots,i_{p-1},i_{p+1},\dots,i_{n}].$
(10)

It follows that the image of the boundary operator $\partial_{n}$ are the
$(n-1)$-chains that are at the boundary of $n$-chains. For instance we have

$\displaystyle\partial_{2}([1,2,3])=[1,2]+[2,3]-[1,3],$ (11)

i.e., the image of a triangle is the linear combination of the links at its
boundary with the correct orientation. Additionally, from this definition it
is also easy to see that a cyclic $n$-chain is in the kernel of the boundary
map $\partial_{n}$ independently of whether the cyclic $n$-chain is the
boundary of a $(n+1)$-chain. For instance we have

$\displaystyle\partial_{1}([1,2]+[2,3]-[1,3])=[2]-[1]+[3]-[2]-[3]+[1]=0,$ (12)

whether the simplex $[123]$ belongs to the simplicial complex or not. One
important topological property of the boundary operator is that the “boundary
of a boundary is null” which implies

$\displaystyle\mbox{im}\ \partial_{n+1}\subseteq\mbox{ker}\ \partial_{n}$ (13)

or equivalently

$\displaystyle\partial_{n}\partial_{n+1}=0.$ (14)

For instance we have

$\displaystyle\partial_{1}\partial_{2}[1,2,3]=\partial_{1}([1,2]+[2,3]-[1,3])=0.$
(15)

The boundary map $\partial_{n}$ can be represented by a $N_{[n-1]}\times
N_{[n]}$ incidence matrix ${\bf B}_{[n]}$ if we adopt as a basis of the space
${\mathcal{C}}_{n}$ an ordered set of the $n$-dimensional simplices $\alpha$,
and as a basis of the space ${\mathcal{C}}_{n-1}$ an ordered set of the
$(n-1)$-dimensional simplices $\hat{\alpha}$.

If the basis of $n$-chains ${\mathcal{C}}_{n}$ is given by the $n$-simplices
$\\{\alpha_{1},\alpha_{2},\ldots\alpha_{s}\ldots\\}$ and the basis of
$(n-1)$-chains ${\mathcal{C}}_{n-1}$ is given by the $(n-1)$-simplices
$\\{\hat{\alpha}_{1},\hat{\alpha}_{2},\ldots\hat{\alpha}_{r}\ldots\\}$ the
action of the boundary map over any arbitrary $n$-dimensional simplex
$\alpha_{s}=[i_{0},i_{1}\ldots,i_{n}]$ given by Eq. (10) can be expressed as

$\displaystyle\partial_{n}\alpha_{s}=\sum_{r=1}^{N_{n-1}}[B_{[n]}]_{rs}\hat{\alpha}_{r}.$
(16)

This equation fully determines the incidence matrices ${\bf B}_{[n]}$. Since
we have seen that the “boundary of a boundary is null” then the incidence
matrices follow ${\bf B}_{[n]}{\bf B}_{[n+1]}={\bf 0}$ and also ${\bf
B}_{[n+1]}^{\top}{\bf B}_{[n]}^{\top}={\bf 0}$.

Figure 4: Example of oriented simplicial complex. This illustrates an example
of a $2$-dimensional oriented simplicial complex having associated incidence
matrices given by Eqs. (27). Reprinted figure with permission from Ref.
millan2020explosive . ©Copyright (2020) by the American Physical Society.

.

As an example we can consider the simplicial complex shown in Figure 4 whose
incidence matrices are given by

$\displaystyle{\bf
B}_{[1]}=\begin{array}[]{c|cccc}&$[1,2]$&$[1,3]$&$[2,3]$&$[3,4]$\\\
\hline\cr$[1]$&-1&-1&0&0\\\ $[2]$&1&0&-1&0\\\ $[3]$&0&1&1&-1\\\
$[4]$&0&0&0&1\\\ \end{array},\ {\bf B}_{[2]}=\begin{array}[]{c|c}&$[1,2,3]$\\\
\hline\cr$[1,2]$&1\\\ $[1,3]$&-1\\\ $[2,3]$&1\\\ $[3,4]$&0\end{array}.$ (27)

The graph Laplacian can be expressed in terms of the incidence matrix ${\bf
B}_{[1]}$ of the graph as

$\displaystyle{\bf L}_{[0]}={\bf B}_{[1]}{\bf B}^{\top}_{[1]}.$ (28)

Similarly, in a simplicial complex the higher-order Laplacian ${\bf L}_{[n]}$
(with $n>0$) simplices2 ; thesis ; Barbarossa  is the $N_{[n]}\times N_{[n]}$
matrix defined as

$\displaystyle{\bf L}_{[n]}={\bf L}^{[\text{down}]}_{[n]}+{\bf
L}^{[\text{up}]}_{[n]},$ (29)

where

$\displaystyle{\bf L}^{[\text{down}]}_{[n]}={\bf B}^{\top}_{[n]}{\bf
B}_{[n]},\quad{\bf L}^{[\text{up}]}_{[n]}={\bf B}_{[n+1]}{\bf
B}^{\top}_{[n+1]}.$ (30)

The $n$-Laplacian is positive semi-definite and, therefore, it has $N_{[n]}$
non negative eigenvalues
$0\leq\lambda_{1}\leq\lambda_{2}\leq\ldots\lambda_{r}\leq\ldots\leq\lambda_{N_{[n]}}.$
A notable property of the spectrum of the $n$-th order Laplacian is that the
degeneracy of zero eigenvalues is given by the $n$-th Betti number. Despite
the fact that the construction of the higher-order Laplacians described above
seems to rely on the choice of the orientations adopted for the simplices of
the simplicial complex, it is possible to show that the higher-order
Laplacians are independent on the orientation of the simplices as long as such
orientation is induced by the labeling of the nodes.

We note that also the up and the down Laplacians are positive semi-definite
and that from the definition of these matrices it follows immediately than the
non-zero eigenvalues in the spectrum of ${\bf L}_{[n]}^{[\text{up}]}$ are the
same as the non-zero eigenvalues in the spectrum of ${\bf
L}_{[n]}^{[\text{down}]}$.

By considering the property of the incidence matrices such that ${\bf
B}_{[n]}{\bf B}_{[n+1]}=0,$ it is possible to derive the Hodge decomposition
of the space of $n$-chains, which reads

$\displaystyle\mathbb{R}^{D_{n}}=\mbox{img}({\bf
B}_{[n]}^{\top})\oplus\mbox{ker}({\bf L}_{[n]})\oplus\mbox{img}({\bf
B}_{[n+1]}).$ (31)

This implies that the higher-order Laplacian ${\bf L}_{[n]}$ can be
simultaneously diagonalized with the $n$-th order up and $n$-th order down
Laplacian and that the non-zero eigenvectors of ${\bf L}_{[n]}$ are either
non-zero eigenvector of ${\bf L}_{[n]}^{[\text{down}]}$ or non-zero
eigenvectors of ${\bf L}_{[n]}^{[\text{up}]}$. Therefore there is a basis in
which ${\bf L}_{[n]}$, ${\bf L}_{[n]}^{[\text{down}]}$ and ${\bf
L}_{[n]}^{[\text{up}]}$ have diagonal form given by

$\displaystyle{\bf U}^{-1}{\bf L}_{[n]}{\bf U}=\left(\begin{matrix}{\bf
D}^{[\text{down}]}_{[n]}&{\bf 0}&{\bf 0}\\\ {\bf 0}&{\bf 0}&{\bf 0}\\\ {\bf
0}&{\bf 0}&{\bf D}^{[\text{up}]}_{[n]}\end{matrix}\right),\ {\bf U}^{-1}{\bf
L}_{[n]}^{[\text{down}]}{\bf U}=\left(\begin{matrix}{\bf
D}^{[\text{down}]}_{[n]}&{\bf 0}&{\bf 0}\\\ {\bf 0}&{\bf 0}&{\bf 0}\\\ {\bf
0}&{\bf 0}&{\bf 0}\end{matrix}\right),\ {\bf U}^{-1}{\bf
L}_{[n]}^{[\text{up}]}{\bf U}=\left(\begin{matrix}{\bf 0}&{\bf 0}&{\bf 0}\\\
{\bf 0}&{\bf 0}&{\bf 0}\\\ {\bf 0}&{\bf 0}&{\bf
D}^{[\text{up}]}_{[n]}\end{matrix}\right),$

where ${\bf D}_{[n]}^{[\text{up}]}$ and ${\bf D}_{[n]}^{[\text{down}]}$ are
diagonal matrices having positive diagonal elements.

###  3.4 Higher order spectral dimension

The notion of spectral dimension (see section 3.2) can be generalized to
$n$-order up-Laplacians with important consequences for higher-order
simplicial complex dynamics. Here, we will focus on the model of NGF
simplicial complexes introduced in section 2.1. In section 3.2, we have shown
that the graph Laplacian of NGFs displays a finite spectral dimension
$d_{S}^{[0]}$ mulder2018network ; millan2018complex ;
millan2019synchronization . Interestingly, higher-order up-Laplacians ${\bf
L}_{[n]}^{[\text{up}]}$ and the higher-order down-Laplacians ${\bf
L}_{[n]}^{[\text{down}]}$ of NGFs also display a finite spectral dimension.

Figure 5:  The cumulative density of non-zero eigenvalues
$\rho_{c}^{up}(\lambda)$ of the ${\bf L}^{[\text{up}]}_{[n]}$ for NGF of
dimension $d=3$ and flavor $s=-1$ (panel (a)), $s=0$ (panel (b)), and $s=1$
(panel (c)) for $n=0$ (blue solid lines), $n=1$ (red dashed lines), $n=2$
(yellow dotted lines) and $n=3$ (purple dot-dashed lines). Considered NGFs
sizes are $N_{[0]}=2000$ nodes, $N_{[1]}=5994$ links, $N_{[2]}=5992$
triangles, and $N_{[3]}=1997$ tetrahedra. Adapted figure from
torres2020simplicial .

In particular, the higher-order up-Laplacians of NGFs display a finite
spectral dimension $d_{S}^{[n]}$ depending on the order $n$, the dimension of
the simplicial complex $d$ and the flavor parameter $s$ torres2020simplicial .
Therefore, we can define different spectral dimensions for $0<n<d-1$. In order
to show this remarkable geometrical property of NGFs in Figure 5 we provide
numerical evidence of the scaling of the cumulative density of non-zero
eigenvalues $\rho_{c}^{\text{up}}(\lambda)$ of the ${\bf
L}^{[\text{up}]}_{[n]}$ with $\lambda$ for $\lambda\ll 1$, given by

$\displaystyle\rho_{c}^{[\text{up}]}(\lambda)\propto\lambda^{d_{S}^{[n]}/2},$
(32)

for different value of the order $n$ of the up-Laplacian, and the flavor $s$
of the NGF. It follows that an NGF model is not characterized by a single
spectral dimension, i.e., the spectral dimension $d_{S}^{[0]}$ of the graph
Laplacian, rather the NGF simplicial complexes have a higher-order network
geometry encoded in a vector of spectral dimensions

$\displaystyle{\bf
d}_{S}=\left(d_{S}^{[0]},d_{S}^{[1]},d_{S}^{[2]}\ldots,d_{S}^{[d-1]}\right).$
(33)

Consequently, the diffusion dynamics defined on simplices of different order
$n$ of the same NGF simplicial complex can be significantly different
torres2015brain . We finally note that for deterministic Apollonian and
pseudo-fractal simplicial complexes that constitute the deterministic
counterpart of NGF simplicial complexes, the higher-order spectral dimension
can be predicted analytically by the real-space renormalization group
reitz2020higher  showing that the higher-order spectral dimension of these
structures depends on the order $n$ and remains finite as long as $0\leq
n<d-1$.

##  4 Simplicial synchronization

Synchronization is a fundamental dynamical state observed in a wide variety of
complex systems and capturing among other phenomena important aspects of brain
dynamics and circadian rhythms. The Kuramoto model kuramoto1975 ;
strogatz2000kuramoto ; acebron2005kuramoto ; arenas2008synchronization ;
boccaletti2018synchronization  is a stylized model that explains how coupled
oscillators, that in absence of interactions would have different intrinsic
frequencies, can start to follow a collective coherent motion when their
coupling constant $\sigma$, measuring the strength of their interaction, is
larger than a critical value $\sigma_{c}$ also called synchronization
threshold.

In order to model the coupling between the oscillators the Kuramoto model
considers a network of $N_{[0]}$ nodes and associates a phase $\theta_{i}$ to
each node $i\in\\{1,2,\ldots,N_{[0]}\\}$ of the network. Therefore in the
Kuramoto model, the dynamical state of the network is determined by the vector
$\bm{\theta}$ of phases associated with its nodes given by

$\displaystyle\bm{\theta}=(\theta_{1},\theta_{2},\ldots,\theta_{N_{[0]}})^{\top}.$
(34)

Each phase $\theta_{i}$ describes the dynamical state of an oscillator that in
absence of interactions oscillates at an intrinsic frequency $\omega_{i}$
drawn independently from a distribution $g(\omega)$. Common choices for
$g(\omega)$ are the unimodal Gaussian or Lorentzian distributions.

The equations determining the dynamics of the phases associated with the nodes
include an important contribution indicating the coupling among the phases of
neighbour nodes. This coupling term has a strength modulated by the coupling
constant $\sigma$. In , it is assumed that this contribution expresses the
tendency of the phase of any given node to oscillate together with the phases
of its neighbour nodes. The resulting standard Kuramoto dynamics is captured
by the differential equations

$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sigma\sum_{j=1}^{N_{[0]}}a_{ij}\sin\left(\theta_{j}-\theta_{i}\right),$
(35)

valid for every node $i$ of the network, where $a_{ij}$ is the generic element
of the adjacency matrix of the network. The level of synchronization in the
system is measured by the Kuramoto order parameter,

$Z_{0}=R_{0}\textrm{e}^{\mathrm{i}\Theta}=\frac{1}{N}\sum_{j=1}^{N_{[0]}}\textrm{e}^{\mathrm{i}\theta_{j}},$
(36)

where $R_{0}$ and $\Theta$ are both real and where $0\leq R_{0}\leq 1$
measures the overall coherence and $\Theta=\Theta(t)$ is the phase of global
oscillations.

The relation between the Kuramoto model and the graph Laplacian of the
underlying network is revealed when the Kuramoto model is linearized for
$|\theta_{i}-\theta_{j}|\ll 1$ for every pair of neighbour nodes $(i,j)$. In
this limit the Kuramoto model can be shown to be described by the system of
equations

$\displaystyle\dot{\bm{\theta}}=\bm{\omega}-\sigma{\bf L}_{[0]}\bm{\theta},$
(37)

where $\bm{\omega}$ indicates the vector of elements $\omega_{i}$ with
$i\in\\{1,2,\ldots,N_{[0]}\\}$.

The Kuramoto model has been analytically solved only on a fully connected
network, although important progress has been made in understanding the
Kuramoto model in random complex networks arenas2008synchronization ;
restrepo2005onset ; chavez2005synchronization ,

In the fully connected network and in random networks the Kuramoto model
displays a second order phase transition at the synchronization threshold
$\sigma=\sigma_{c}$ when the number of nodes goes to infinity, i.e.,
$N_{[0]}\to\infty$. For $\sigma<\sigma_{c}$ the Kuramoto model is in an
incoherent state characterized by having a zero order parameter $R_{0}=0$. For
$\sigma>\sigma_{c}$ the Kuramoto model is in a coherent state characterized by
a non-zero order parameter $R_{0}>0$ acebron2005kuramoto ;
strogatz2000kuramoto ; arenas2008synchronization ;
boccaletti2018synchronization .

In this chapter we show how the network geometry and topology of simplicial
complexes, directly acting on the spectral properties of the graph Laplacian
and the higher-order Laplacians, can dramatically change the dynamical
properties of the synchronization process on higher-order networks.

Figure 6:  Schematic representation of the Kuramoto model and the higher-order
Kuramoto model capturing dynamics of topological signals. The Kuramoto model
(panel (a)) captures the emergence of a synchronized state among coupled
oscillators described by phases associated with the nodes of a network. The
higher-order Kuramoto model (panel (b)) reveals the synchronization of
topological signals on simplicial complexes, i.e., oscillators associated not
only to the nodes of a simplicial complex, but also to higher-dimensional
simplices such as links or triangles. Interestingly, topological signals of
different dimension can co-exist and co-evolve and can be non-trivially
coupled leading to simultaneous explosive transitions.

  * •

Simplicial network geometry and the Kuramoto model. First we will illustrate
how the phase diagram of the Kuramoto model changes when the model is defined
on the skeleton of simplicial complexes with distinct higher-order network
geometry. We have previously introduced the finite spectral dimension
$d_{S}^{[0]}$ of the graph Laplacian as a fundamental observable of the
higher-order network geometry of networks and simplicial complexes. In the
following we will discuss how the dynamics of the Kuramoto model depends on
$d_{S}^{[0]}$ revealing that network geometries have a more rich dynamics than
random networks.

  * •

Simplicial network topology and the higher-order Kuramoto model. Secondly we
will reveal how the Kuramoto model can be considered as a limiting case of a
much wider class of higher-order Kuramoto models describing coupling of
topological signals. Topological signals are phases of oscillators associated
not only to the nodes but also to higher-order simplices of a simplicial
complex. For instance topological signals can be associated with both nodes
and links of a simplicial complex as schematically described by Figure 6. In
this case the dynamical state of a simplicial complex is captured by the
vector $\bm{\theta}$ of the phases associated with the nodes (defined in Eq.
(34)) and by the vector $\bm{\phi}$ of the phases associated with the links of
the simplicial complex given by

$\displaystyle\bm{\phi}=\left(\phi_{\ell_{1}},\phi_{\ell_{2}},\ldots,\phi_{\ell_{N_{[1]}}}\right)^{\top}.$
(38)

The phases associated with the links of the simplicial complexes are
topological signals that have the potential to capture the dynamics of fluxes
in brain networks deville_brain  and biological transportation networks
katifori1 . The newly formulated higher-order Kuramoto model opens new
scenarios for characterizing how topology affects dynamics on higher-order
networks and simplicial complexes. As we will discuss in Sec. V the higher-
order Kuramoto model has a linearized dynamics described by the higher-order
Laplacians of the simplicial complex, and can display the simultaneous
explosive synchronization transition of the soleinodal and the irrotational
component of topological signals and even the simultaneous explosive
synchronization transition of topological signals of different dimensions.

##  5 Kuramoto model on simplicial network geometry

###  5.1 Synchronization on simplicial network skeletons with finite spectral
dimension

In order to investigate the role that simplicial network geometry has on the
Kuramoto model we explore the phase diagram of the normalized Kuramoto model
on networks with finite spectral dimension millan2019synchronization . Our
theoretical results are then validated by simulations performed over the
skeleton of the NGF with flavor $s=-1$. Indeed NGF with flavor $s=-1$ provide
a very suitable benchmark to test our theoretical results as they display a
spectral dimension $d_{S}^{[0]}$ that can be changed by tuning the dimension
$d$ of the simplicial complex millan2018complex .

The normalized Kuramoto model determines the dynamics of the phases
$\bm{\theta}$ associated with the nodes of a network. The only difference with
the standard Kuramoto model is that the coupling between the phase of a given
node $i$ and the phases of its neighbour nodes is normalized with the node
degree $k_{i}$. Therefore, the normalized Kuramoto model is dictated by the
differential equations

$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sigma\sum_{j=1}^{N_{[0]}}\frac{a_{ij}}{k_{i}}\sin\left(\theta_{j}-\theta_{i}\right),$
(39)

where here and in the following we consider internal frequencies $\omega_{i}$
drawn independently from a normal distribution, i.e.,
$\omega_{i}\sim\mathcal{N}(0,1).$ The normalization of the coupling term by
the degree of the node $i$ is a very efficient way to screen out the effects
of the heterogeneity of the degrees of the nodes and single out only the
effects due to the geometrical nature of the network of their interactions.

The linearized equation of the normalized Kuramoto model is therefore
determined by the normalized Laplacian ${\hat{\bf L}}_{[0]}$ instead of the
graph Laplacian ${\bf L}_{[0]}$, i.e.,

$\displaystyle\dot{\bm{\theta}}=\bm{\omega}-\sigma{\bf\hat{L}}_{[0]}\bm{\theta}.$
(40)

The analytical investigation of the stability of the synchronized phase
indicates that the spectral dimension $d_{S}^{[0]}$ of the (normalized) graph
Laplacian of the network plays a fundamental role for determining the phase
diagram of the normalized Kuramoto model in the limit of infinite network size
$N_{[0]}\to\infty$ millan2019synchronization . Interestingly, we note that the
spectral dimension of the graph Laplacian and the normalized graph Laplacian
take the same value under very general regularity conditions of the network
burioni1996universal . For the very heterogeneous NGF, numerical results show
that the two spectral dimensions differ by a small amount as long as the
topological dimension $d$ is small. Depending on the value of the spectral
dimension the phase diagram of the normalized Kuramoto model defined in the
infinite network limit changes drastically millan2019synchronization :

  * (1)

For networks with finite spectral dimension $d_{S}^{[0]}\leq 2,$ the Kuramoto
model cannot synchronize and is found in the incoherent state for every value
of the coupling constant $\sigma$.

  * (2)

For networks with spectral dimension $2<d_{S}^{[0]}\leq 4,$ global
synchronization is not achievable in the infinite network limit but an
entrained state can be observed. Therefore it is possible the Kuramoto model
can have a transition between an incoherent state and an entrained phase.

  * (3)

Only for networks with spectral dimension $d_{S}^{[0]}>4$ it is possible to
see a synchronized phase.

These results reveal how the dynamics of the Kuramoto model depends on the
simplicial network geometry on which is defined, and extend previous results
valid on regular lattices of dimension $d$ hong2005collective ;
hong2007entrainment .

###  5.2 Frustrated synchronization on Network Geometry with Flavor

The NGF constitutes a perfect model to investigate numerically the role that
simplicial network geometry has on the dynamics of the Kuramoto model
millan2018complex ; millan2019synchronization . Indeed, as we discussed
previously, the NGF displays an emergent hyperbolic network geometry and for
flavor $s=-1$ generates random hyperbolic manifolds. The simplicial network
geometry of NGF is also reflected on their spectral properties. Specifically,
NGFs have a finite spectral dimension $d_{S}^{[0]}$ which for flavor $s=-1$
and $d\in\\{2,3,4\\}$ can be approximated by $d_{S}^{[0]}\simeq d$. It follows
that by changing the dimension $d$ of the NGF with flavor $s=-1$ we can
explore the dynamics of the Kuramoto model when the global synchronization
state is not stable in the infinite network limit.

A computational finite size analysis of the Kuramoto model millan2018complex ;
millan2019synchronization  reveals, in agreement with the theoretical
expectations, that for $d=2$ and $d_{S}^{[0]}\simeq 2$, global synchronization
is never achieved for large network sizes. On the other hand, synchronization
in NGFs with $d=3$ and $d=4$ and $d_{S}^{[0]}\simeq d$ is only possible for
finite networks. In fact its onset occurs for higher couplings when the system
size increases, revealing, in agreement with the theoretical expectation, that
in the limit $N_{[0]}\to\infty$ this state is never achieved.

Interestingly, we observe that for NGF with $d=3$ and $d=4$ and
$d_{S}^{[0]}\simeq d$ the Kuramoto model exhibits a phase with entrained
synchronization that we call a phase of frustrated synchronization for a wide
range of coupling values. In the frustrated synchronization regime the order
parameter $R_{0}$ displays strong temporal fluctuations. This phase is
observed on finite NGF between the incoherent state and the globally
synchronized state, and for $d=3$ a much broader regime of large fluctuations
is observed than for the $d=4$ case.

Interestingly, the frustrated synchronization phase of NGFs is not only
characterized by strong temporal fluctuations of the global order parameter
but displays also strong spatial fluctuations induced by the non-trivial
hyperbolic network geometry of the NGF. As such, the frustrated
synchronization phase can capture an important mechanism for inducing spatio-
temporal fluctuations in brain dynamics.

The hyperbolic network geometry of NGF has a strong hierarchical nature that
is responsible for the emergence of a relevant community structure. In order
to study how the dynamics of the Kuramoto model is affected by the community
structure of the NGFs, we define mesoscopic synchronization order parameters
$Z_{\textrm{mod}}$ that characterize the dynamical state of each community:

$Z_{\textrm{mod}}=R_{\textrm{mod}}\textrm{e}^{\mathrm{i}\Theta_{\textrm{mod}}}=\frac{1}{\left|\mathcal{C}\right|}\sum_{j\in\mathcal{C}}\textrm{e}^{\mathrm{i}\theta_{j}},$
(41)

where $\mathcal{C}$ is the set of nodes in the community and
$\left|\mathcal{C}\right|$ the total number of nodes in said community. Figure
7b displays the trajectory of $Z_{\textrm{mod}}=Z_{\textrm{mod}}(t)$ in the
complex plane for some exemplary modules of an NGF with flavour $s=-1$ in
$d=3$, for the coupling that leads to the largest fluctuations of $R_{0}$ as a
function of time. As shown in the figure, different modules display different
synchronization regimes and may oscillate at different frequencies. Due to the
underlying geometrical structure of NGFs, these modules correspond to
spatially localized regions.

We note here that the frustrated synchronization observed in NGF can be
related with analogous phases observed in other hierarchical models
moretti2013griffiths ; villegas2014frustrated  where temporal fluctuations of
the synchronization order parameter are observed. However, the combination of
both temporal and spatial fluctuations is a specific property of the
frustrated synchronization in NGF due to their rich simplicial network
geometry.

Figure 7: Frustrated synchronization on NGF characterized by spatio-temporal
fluctuations of the order parameter. Panel (a) The synchronization order
parameter $R_{0}(T)$ calculated at time $T$ is plotted versus the coupling
constant $\sigma$ revealing the regime of frustrated synchronization (top
panel). The time series $R_{0}(t)$ are shown (bottom panel) for the values of
the coupling constant $\sigma$ indicated by arrows in the top panel. These
time series reveal the temporal fluctuations of the order parameter in the
frustrated synchronization regime. Panel (b) The spatial fluctuations of the
order parameter in the frustrated synchronization regime are revealed by the
local order parameter $Z_{\textrm{mod}}$ of four different communities of the
NGF calculated for $\sigma=5$ (top panel). In this representation, a circular
trajectory describes a situation of global oscillations of the nodes in the
community, with constant $R_{\textrm{mod}}(t)$. Random trajectories around the
origin describe unsynchronized communities. Partially synchronized
communities, on the other hand, may describe more complex trajectories. The
bottom panel shows the corresponding time series of $R_{\textrm{mod}}(t)$.
Panel (c) Synchronization transition as given by the mean order parameter
$\bar{R}_{0}$ averaged over time and its variance $\sigma_{R_{0}}$, as
functions of $\sigma$, for different network sizes $N=100$ (blue), $200$,
$400$, …, $3200$ (black). The apparent onset of the synchronized regime is
retarded to large values of the coupling constant $\sigma$ for large network
sizes, revealing that the NGF of dimension $d=3$ cannot sustain a synchronized
phase in the limit of infinite network sizes. All simulations reported in the
figures are obtained for NGF of flavor $s=-1$, dimension $d=3$ and number of
nodes $N_{[0]}=1600$. In panel (a) $T=500$. Figure adapted from
millan2018complex .

##  6 Higher-order Kuramoto model: a topological approach to synchronization

###  6.1 Synchronization of topological signals

Simplicial complexes are formed by nodes and higher-order simplices including
links, triangles, tetrahedra, and so on. As such, simplicial complexes have
the ability to sustain topological signals, i.e., dynamic variables not only
associated with the nodes of their network skeleton but also associated with
links, triangles, and so on millan2020explosive ; ghorbanchian2020higher ;
Barbarossa . Topological signals have the ability to capture dynamics
associated with links, such as fluxes in brain networks deville_brain  and
biological transportation networks katifori1 ; katifori2 .

Here we will present the higher-order Kuramoto model millan2020explosive  that
reveals how topological signals can undergo continuous and discontinuous
synchronization transitions. Interestingly, we will observe that this
synchronization transition can be detected only if the signals are filtered
with the appropriate topological operators. Therefore, the model not only
captures a new topological critical phenomenon but also prescribes a way to
process real data in order to investigate whether this topological
synchronization phenomenon can be observed in real systems such as the brain
or biological transportation networks.

###  6.2 Higher-order Kuramoto dynamics

The higher-order Kuramoto model millan2020complex  describes the
synchronization of topological signals defined on simplices of dimension $n$.
For instance, one can consider topological signals defined on the links of a
simplicial complex (case $n=1$) or alternatively one can consider signals
defined on the triangles of a simplicial complex (case $n=2$). The higher-
order Kuramoto model is the most natural extension of the Kuramoto model to
capture the synchronization of higher-order topological signals.

The standard Kuramoto dynamics describes the dynamics of the phases
$\bm{\theta}$ associated with the nodes of the network. This dynamics, defined
by Eq. (35), can be expressed in terms of the incidence matrix $\bm{B}_{[1]}$
(see Appendix) as

$\displaystyle\dot{\bm{\theta}}=\bm{\omega}-\sigma\bm{B}_{[1]}\sin\bm{B}_{[1]}^{T}\bm{\theta}.$
(42)

The higher-order Kuramoto dynamics describes instead the dynamics of
topological signals $\bm{\phi}$ with $\phi_{\alpha}$ indicating the phase
associated with the simplex $\alpha$ of dimension $n>0$. Therefore, using the
insights coming from algebraic topology, the natural definition of the simple
higher-order Kuramoto model is

$\displaystyle\dot{\bm{\phi}}=\bm{\hat{\omega}}-\sigma\bm{B}_{[n+1]}\sin\bm{B}_{[n+1]}^{T}\bm{\phi}-\sigma\bm{B}_{[n]}^{T}\sin\bm{B}_{[n]}\bm{\phi},$
(43)

where $\bm{\hat{\omega}}$ is the vector of intrinsic frequencies
$\hat{\omega}_{\alpha}$ associated with each $n$-dimensional simplex $\alpha$
drawn independently from a normal distribution, i.e.,
$\hat{\omega}_{\alpha}\sim\mathcal{N}(\Omega_{1},1/\tau_{1})$. As the standard
Kuramoto model can be related to the graph Laplacian via linearization (see
Eq. (37)), the higher-order Kuramoto model can be related to the higher-order
Laplacian upon linearization, leading to

$\displaystyle\dot{\bm{\phi}}=\bm{\hat{\omega}}-\sigma{\bf L}_{[n]}\bm{\phi}.$
(44)

The higher-order Kuramoto dynamics defined on topological signals $\bm{\phi}$
associated with $n$ dimensional simplices can be projected on $(n+1)$ and
$(n-1)$ dimensional simplices by applying to the signals the incidence
matrices. For instance a dynamics defined on topological signals associated
with links can be projected on nodes or on triangles. Specifically, we have
that the projected dynamics $\bm{\phi}^{[-]}$ on $(n-1)$-dimensional simplices
and the projected dynamics $\bm{\phi}^{[+]}$ on $(n+1)$-dimensional simplices
is given by

$\displaystyle\bm{\phi}^{[-]}$ $\displaystyle=$ $\displaystyle{\bf
B}_{[n]}\bm{\phi},$ $\displaystyle\bm{\phi}^{[+]}$ $\displaystyle=$
$\displaystyle{\bf B}_{[n+1]}^{\top}\bm{\phi},$ (45)

where, for $n=1$, ${\bf B}_{[n]}$ indicates the discrete divergence and ${\bf
B}_{[n+1]}^{\top}$ indicates the discrete curl. Therefore $\bm{\phi}^{[-]}$
indicates the irrotational component of $\bm{\phi}$ while $\bm{\phi}^{[+]}$
indicates the solenoidal component of $\bm{\phi}$. For the simple higher-order
Kuramoto model defined in Eq. (43), by recalling that ${\bf
B}_{[n+1]}^{\top}{\bf B}_{[n]}^{\top}={\bf 0}$ and that ${\bf B}_{[n]}{\bf
B}_{[n+1]}={\bf 0},$ it follows that the projected topological signals
$\bm{\phi}^{[-]}$ and $\bm{\phi}^{[+]}$ obey the uncoupled system of equations

$\displaystyle\dot{\bm{\phi}}^{[-]}$ $\displaystyle=$
$\displaystyle\bm{B}_{[n]}\bm{\hat{\omega}}-\sigma\bm{L}_{[n-1]}^{[\textrm{up}]}\sin\left(\bm{\phi}^{[-]}\right),$
$\displaystyle\dot{\bm{\phi}}^{[+]}$ $\displaystyle=$
$\displaystyle\bm{B}^{T}_{[n+1]}\bm{\hat{\omega}}-\sigma\bm{L}_{[n+1]}^{[\textrm{down}]}\sin\left(\bm{\phi}^{[+]}\right).$
(46)

Therefore, the solenoidal and the irrotational components of the topological
signals are decoupled for the simple higher-order Kuramoto model.

The higher-order Kuramoto dynamics is remarkable from two perspectives:

  * (1)

First of all it uses topology to naturally define the higher-order
interactions between the topological signals. Indeed the incidence matrices
define higher-order interactions with a clear prescription indicating the
coupled variables and the sign of their interactions. For example, the higher-
order Kuramoto dynamics for $n=1$ dimensional simplices of the simplicial
complex shown in Figure 4 reads

$\displaystyle\dot{\theta}_{[1,2]}$ $\displaystyle=$
$\displaystyle\hat{\omega}_{[1,2]}-\sigma\sin(\theta_{[2,3]}-\theta_{[1,3]}+\theta_{[1,2]})-\sigma\left[\sin(\theta_{[1,2]}-\theta_{[2,3]})+\sin(\theta_{[1,3]}+\theta_{[1,2]})\right],$
$\displaystyle\dot{\theta}_{[1,3]}$ $\displaystyle=$
$\displaystyle\hat{\omega}_{[1,3]}+\sigma\sin(\theta_{[2,3]}-\theta_{[1,3]}+\theta_{[1,2]})-\sigma\left[\sin(\theta_{[1,3]}+\theta_{[1,2]})+\sin(\hat{\theta}_{[3]})\right],$
$\displaystyle\dot{\theta}_{[2,3]}$ $\displaystyle=$
$\displaystyle\hat{\omega}_{[2,3]}-\sigma\sin(\theta_{[2,3]}-\theta_{[1,3]}+\theta_{[1,2]})-\sigma\left[\sin(\theta_{[2,3]}-\theta_{[1,2]})+\sin(\hat{\theta}_{[3]})\right],$
$\displaystyle\dot{\theta}_{[3,4]}$ $\displaystyle=$
$\displaystyle\hat{\omega}_{[3,4]}-\sigma\left[\sin(\theta_{[3,4]})-\sin(\hat{\theta}_{[3]})\right],$
(47)

with $\hat{\theta}_{[3]}$ indicating the three-body interaction

$\displaystyle\hat{\theta}_{[3]}=\theta_{[13]}+\theta_{[23]}-\theta_{[34]}.$
(48)

Therefore, the choice of the higher-order interactions in the higher-order
Kuramoto model is naturally dictated by topology.

  * (2)

The synchronization of the higher-order Kuramoto model is only detectable if
the right topological filtering of the data is performed. Indeed the naïve
order parameter

$\displaystyle
R_{n}=\frac{1}{N_{[n]}}\left|\sum_{\alpha=1}^{N_{[n]}}e^{\mathrm{i}\phi_{\alpha}}\right|$
(49)

associated with the unfiltered topological signal $\bm{\phi}$ does not detect
any synchronization transition (see Figure 8). Instead, the order parameter
associated with the solenoidal and the irrotational components of the
topological signal do detect the synchronization transition of the topological
signals (see Figure 8). These order parameters are given by

$\displaystyle
R_{n}^{[-]}=\frac{1}{N_{[n-1]}}\left|\sum_{\alpha=1}^{N_{[n-1]}}e^{\mathrm{i}\phi_{\alpha}^{[-]}}\right|,\
\ \
R_{n}^{[+]}=\frac{1}{N_{[n+1]}}\left|\sum_{\alpha=1}^{N_{[n+1]}}e^{\mathrm{i}\phi_{\alpha}^{[+]}}\right|,$
(50)

or alternatively by

$\displaystyle
R_{n}^{\downarrow}=\frac{1}{N_{[n]}}\left|\sum_{\alpha=1}^{N_{[n]}}e^{\mathrm{i}\phi_{\alpha}^{\downarrow}}\right|,\
\ \
R_{n}^{\uparrow}=\frac{1}{N_{[n]}}\left|\sum_{\alpha=1}^{N_{[n]}}e^{\mathrm{i}\phi_{\alpha}^{\uparrow}}\right|,$
(51)

where $\bm{\phi}^{\downarrow}={\bf L}_{[n]}^{[\text{down}]}\bm{\phi}$ and
$\bm{\phi}^{\uparrow}={\bf L}_{[n]}^{[\text{up}]}\bm{\phi}$.

The synchronization transition described by the simple higher-order Kuramoto
model leads to a continuous transition occurring at zero coupling, i.e., the
synchronization threshold is $\sigma_{c}=0$ as long as $n>0$. However, the
higher-order Kuramoto model admits a formulation called explosive higher order
Kuramoto model that displays instead a discontinuous transition at a non zero
coupling $\sigma_{c}>0$.

Figure 8:  Higher-order Kuramoto dynamics. The synchronization of higher-order
topological signals is captured by the simple and explosive higher-order
Kuramoto models. Panel (a) The order parameters $R_{1}^{[+]}$ and
$R_{1}^{[-]}$ reveal the synchronization transition of topological signals
defined on links ($n=1$) of the configuration model of simplicial complexes.
Panel (b) The naïve order parameter $R_{1}$ does not reveal the
synchronization transition of topological signals defined on the links of the
configuration model of simplicial complex, while the order parameters
$R_{1}^{\uparrow}$ and $R_{1}^{\downarrow}$, sensible to the irrotational and
solenoidal decomposition of the signal, do reveal the transition (as well as
$R_{1}^{[+]}$ and $R_{1}^{[-]}$ shown in panel (a)). The underlying network is
the same for both panels and has $N_{[0]}=1000$ nodes, $N_{[1]}=5299$ links
and $N_{[2]}=4147$ triangles. The generalized degree of the nodes is power-law
distributed with power-law exponent $\gamma=2.8$. Panel (c) shows the
theoretical expectations provided by an effective phenomenological model
treated in Ref. millan2020explosive  for the simple (top) and explosive
(bottom) higher-order Kuramoto models. Reprinted figure with permission from
Ref. millan2020explosive . ©Copyright (2020) by the American Physical Society.

The explosive higher-order Kuramoto dynamics millan2020explosive  implements
an adaptive coupling of the projected dynamics of $\bm{\phi}^{[+]}$ and
$\bm{\phi}^{[-]}$ through their global order parameters. The adopted adaptive
coupling is inspired by analogous couplings previously applied to multilayer
and simple networks zhang2015explosive . The explosive higher-order Kuramoto
model millan2020explosive  is defined by the system of equations

$\dot{\bm{\phi}}=\bm{\hat{\omega}}-\sigma
R^{[-]}\bm{B}_{[n+1]}\sin\bm{B}^{T}_{[n+1]}\phi-\sigma
R^{[+]}\bm{B}^{T}_{[n]}\sin\bm{B}_{[n]}\phi.$ (52)

It follows that the dynamics projected on the $(n+1)$ and $(n-1)$-dimensional
simplices now obeys the coupled system of equations

$\displaystyle\dot{\bm{\phi}}^{[+]}$ $\displaystyle=$
$\displaystyle\bm{B}^{T}_{[n+1]}\bm{\hat{\omega}}-\sigma
R^{[-]}\bm{L}_{[n+1]}^{[\textrm{down}]}\sin\left(\bm{\phi}^{[+]}\right),$
$\displaystyle\dot{\bm{\phi}}^{[-]}$ $\displaystyle=$
$\displaystyle\bm{B}_{[n]}\bm{\hat{\omega}}-\sigma
R^{[+]}\bm{L}_{[n-1]}^{[\textrm{up}]}\sin\left(\bm{\phi}^{[-]}\right).$ (53)

Numerical simulations on the configuration model of simplicial complexes with
power-law distribution of generalized degrees reveal that the explosive
higher-order Kuramoto model displays a discontinuous phase transition. The
nature of the transition confirms the theoretical expectations obtained with
an approximate phenomenological approach. This transition is clearly detected
by a discontinuity in $R^{[+]}_{n}$ and $R^{[-]}_{n}$ and in
$R^{\downarrow}_{n}$ and $R^{\uparrow}_{n}$ as well, but is not captured by
the naïve order parameter $R_{n}$ (see Fig. 8).

In Ref. millan2020explosive  it has been shown that the nature of the phase
transition does not change if the generalized degree distribution is more
uniform or if the simplicial complex has a non trivial network geometry.
Interestingly, simple and explosive higher-order Kuramoto models can be
investigated on simplicial complexes constructed from real connectomes leading
to continuous (for the simple higher-order Kuramoto model) and discontinuous
(for the explosive higher-order Kuramoto model) synchronization transitions.

###  6.3 Coupled topological signals

So far we have considered the synchronization of topological signals defined
on simplices of dimension $n>0$. However, topological signals associated with
simplices of different dimension can co-exist and co-evolve. For instance,
phases associated with the nodes of a simplicial complex can be coupled with
phases associated with its links. In this section we will show how different
topological signals, i.e., phases defined on simplicial complexes of different
dimensions, can be coupled to each other leading to simultaneous explosive
synchronization transitions. For simplicity of presentation we will focus on
phases defined on nodes and links, but we emphasize that our formalism allows
one to consider the interaction of more general topological signals.

We start by considering Model 1, an explosive higher-order Kuramoto model of
coupled signals of nodes and links. This model differs from the explosive
higher-order Kuramoto model defined in the previous section as it includes an
additional adaptive coupling among the topological signals of nodes and links,
and obeys the equations

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=$
$\displaystyle\bm{\omega}-\sigma R_{1}^{[-]}{\bf B}_{[1]}\sin({\bf
B}^{\top}_{[1]}\bm{\theta}),$ $\displaystyle\dot{\bm{\phi}}$ $\displaystyle=$
$\displaystyle\hat{\bm{\omega}}-\sigma R_{0}R_{1}^{[+]}{\bf
B}^{\top}_{[1]}\sin({\bf B}_{[1]}\bm{\phi})-R_{1}^{[-]}\sigma{\bf
B}_{[2]}\sin({\bf B}_{[2]}^{\top}\bm{\phi}).$ (54)

This model simulated in a wide variety of simplicial complexes including the
NGF, the configuration model of simplicial complexes, and the clique complex
of real connectomes displays a simultaneous explosive (i.e., discontinuous)
synchronization of the topological signals defined on nodes and of the
soleinodal and irrotational component of the topological signals defined on
links ghorbanchian2020higher . Indeed, at a critical threshold
$\sigma=\sigma_{c}$ we observe a discontinuity in the three order parameters
$R_{0},R_{1}^{[-]}$ and $R_{1}^{[+]}$. In Figure 9 we present numerical
evidence of this discontinuous transition by displaying the corresponding
hysteresis loop in the order parameters. In particular, instead of plotting
the order parameters obtained at each value of $\sigma$ starting from random
initial conditions as in Figure 8 here we display the order parameters along
the forward and backward synchronization transitions obtained by first
adiabatically increasing and then decreasing the coupling constant $\sigma$.

t]

Figure 9:  Simultaneous explosive synchronization of coupled topological
signals. Model 1 displays the simultaneous explosive (discontinuous)
transition captured by the order parameters $R_{0}$, $R_{1}^{[-]}$ and
$R^{[+]}_{1}$. Panels (a), (b), and (c) provide numerical evidence of this
discontinuous transition on different simplicial network topologies: the NGF
with flavor $s=-1$ and dimension $d$ (panel (a)), the configuration model of
pure $d=3$ simplicial complex with power-law generalized degree distribution
with power-law exponent $\gamma=2.8$ (panel (b)) and the clique complex of the
Caenorhabditis elegans (C.elegans) connectome coming from Ref. celegans
(panel (c)). Figure adapted from ghorbanchian2020higher .

Interestingly, this model (Model 1) admits a variation (Model 2) that
describes the simultaneous synchronization of the topological signals
associated with the nodes and the irrotational component of the topological
signal associated with the links. This simpler version of the explosive
higher-order Kuramoto model of coupled topological signals can be also defined
on pairwise networks and is amenable to an exact analytical treatment on fully
connected networks and to an accurate annealed approximation solution on
random networks with a given degree distribution. More specifically, Model 2
only couples the signal of the nodes with the signal of the links projected to
the nodes, as described by the differential equations

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=$
$\displaystyle\bm{\omega}-\sigma R_{1}^{[-]}{\bf B}_{[1]}\sin({\bf
B}^{\top}_{[1]}\bm{\theta}),$ (55) $\displaystyle\dot{\bm{\phi}}$
$\displaystyle=$ $\displaystyle\hat{\bm{\omega}}-\sigma R_{0}{\bf
B}^{\top}_{[1]}\sin({\bf B}_{[1]}\bm{\phi})-\sigma{\bf B}_{[2]}\sin({\bf
B}_{[2]}^{\top}\bm{\phi}).$ (56)

Projecting the dynamics of the link phases down to 0-simplices as in the
previous section, we introduce $\bm{\psi}$ for simplicity of notation with

$\displaystyle{\bm{\psi}}\equiv{\bm{\phi}}^{[-]}={\bf B}_{[1]}{\bm{\phi}}.$
(57)

By left multiplying Eq. (56) by ${\bf B}_{[1]}$, we obtain the closed system
of equations for ${\bm{\theta}}$ and ${\bm{\psi}}$

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=$
$\displaystyle\bm{\omega}-\sigma R_{1}^{[-]}{\bf B}_{[1]}\sin({\bf
B}^{\top}_{[1]}\bm{\theta}),$ $\displaystyle\dot{\bm{\psi}}$ $\displaystyle=$
$\displaystyle\tilde{\bm{\omega}}-\sigma R_{0}{\bf L}_{[0]}\sin({\bm{\psi}}),$
(58)

where $\tilde{\bm{\omega}}={\bf B}_{[1]}\hat{\bm{\omega}}$. Here we assume
$\omega\sim\mathcal{N}(0,1)$ and
$\hat{\omega}_{\alpha}\sim\mathcal{N}(\Omega_{1},1/\tau_{1})$. With these
hypotheses the internal frequencies of the links projected on the nodes
$\\{\tilde{\omega}_{i}\\}_{i=1,2,\ldots,N_{[0]}}$ are Gaussian correlated
variables with average

$\displaystyle\left\langle{\tilde{\omega}_{i}}\right\rangle=\left[\sum_{j<i}a_{ij}-\sum_{j>i}a_{ij}\right]\Omega_{1}$
(59)

and with correlation matrix ${\bf C}$ of elements
$C_{ij}=\left\langle{\tilde{\omega}_{i}\tilde{\omega}_{j}}\right\rangle-\left\langle{\tilde{\omega}_{i}}\right\rangle\left\langle{\tilde{\omega}_{j}}\right\rangle$
given by

$\displaystyle{\bf C}={\bf L}_{[0]}\frac{1}{\tau_{1}^{2}}.$ (60)

To understand the nature of the synchronization transition analytically when
Model 2 is defined on an uncorrelated random graph, in the following we
discuss the solution of the model in the annealed approximation. The annealed
approximation is a widely used approximation to study dynamical processes on
random uncorrelated networks which consists in substituting the adjacency
matrix entries of the network $a_{ij}$ by their average values in an
uncorrelated random network with given degree sequence, i.e.,

$\displaystyle a_{ij}\to\frac{k_{i}k_{j}}{\langle{k}\rangle N}.$ (61)

Using this approximation, Eqs. (58) can be recast into the differential
equations

$\displaystyle\dot{\bm{\theta}}$ $\displaystyle=$
$\displaystyle\bm{\omega}-\sigma R_{1}^{[-]}\hat{R}_{0}{\bf
k}\cdot\sin(\bm{\theta}-\hat{\Theta}),$ $\displaystyle\dot{\bm{\psi}}$
$\displaystyle=$ $\displaystyle\tilde{\bm{\omega}}+\sigma
R_{0}\hat{R}_{1}^{[-]}{\bf k}\sin\hat{\Psi}-\sigma R_{0}{\bf
k}\odot\sin\bm{\psi},$ (62)

where $\odot$ indicates the Hadamard product (element by element
multiplication) and where two auxiliary complex order parameters are defined
as

$\displaystyle\hat{R}_{0}e^{\mathrm{i}\hat{\Theta}}=\sum_{i=1}^{N_{[0]}}\frac{k_{i}}{\langle{k}\rangle{N_{[0]}}}e^{\mathrm{i}\theta_{i}},\quad\hat{R}_{1}^{[-]}e^{\mathrm{i}\hat{\Psi}}=\sum_{i=1}^{N_{[0]}}\frac{k_{i}}{\langle{k}\rangle{N_{[0]}}}e^{\mathrm{i}\psi_{i}},$
(63)

with $\hat{R}_{0},\hat{\Theta},\hat{R}_{1}^{{[-]}}$ and $\hat{\Psi}$ being
real. Let us indicate with $g(\omega)$ the probability distribution of the
internal frequencies of the nodes and with $G_{i}(\tilde{\omega})$ the
marginal probability distribution of the internal frequencies of the links
projected on node $i$, i.e., the marginal probability that
$\tilde{\omega}_{i}=\tilde{\omega}$. With this notation it is possible to
derive the analytic solution of Eqs. (62) which gives the following expression
for the order parameters

$\displaystyle R_{0}=\frac{1}{{N_{[0]}}}\sum_{i=1}^{N_{[0]}}r_{0}(i),$
$\displaystyle\hat{R}_{n}=\sum_{i=1}^{N_{[0]}}\frac{k_{i}}{\langle{k}\rangle{N_{[0]}}}r_{0}(i),$
$\displaystyle
R_{1}^{[-]}=\frac{1}{{N_{[0]}}}\sum_{i=1}^{N_{[0]}}r_{1}^{[-]}(i),$
$\displaystyle\hat{R}_{1}^{[-]}=\sum_{i=1}^{N_{[0]}}\frac{k_{i}}{\langle{k}\rangle{N_{[0]}}}r_{1}^{[-]}(i),$

with ${r}_{0}(i)$ and ${r}_{1}^{[-]}(i)$ given by

$\displaystyle{r}_{0}(i)$ $\displaystyle=$
$\displaystyle\int_{|\hat{c}_{i}|<1}d\omega
g(\omega)\sqrt{1-\left(\frac{\omega-\Omega_{0}}{\sigma
k_{i}\hat{R}_{0}R_{1}^{[-]}}\right)^{2}},$ $\displaystyle{r}_{1}^{[-]}(i)$
$\displaystyle=$
$\displaystyle\int_{|\hat{d}_{i}|<1}d\tilde{\omega}G_{i}(\tilde{\omega})\sqrt{1-\left(\frac{\tilde{\omega}}{\sigma
k_{i}\hat{R}_{0}}\right)^{2}},$

and $\hat{c}_{i}$ and $\hat{d}_{i}$ indicating

$\displaystyle\hat{c}_{i}=\frac{\omega-\Omega_{0}}{\sigma
k_{i}\hat{R}_{0}R_{1}^{[-]}},\quad\hat{d}_{i}=\frac{\tilde{\omega}}{\sigma
k_{i}\hat{R}_{0}}.$ (66)

Figure 10 shows excellent agreement between the simulation results of Model 2
and the analytical prediction obtained in the annealed approximation for a
Poisson network with average degree $c=12$ (panel (a)) and for an uncorrelated
scale-free network with minimum degree $m=6$ and power-law exponent
$\gamma=2.5$ (panel (b)).

t]

Figure 10: Theoretical prediction for Model 2 of explosive higher-order
Kuramoto model of coupled dynamical signals applied to random and fully
connected networks. The hysteresis loop for Model 2 is shown for different
types of networks. The backward transition is fully captured by the
theoretical expectations (solid black lines). Panels (a), (b) and (c) show the
order parameters $R_{1}$ and $R_{1}^{[-]}$ versus the coupling constant
$\sigma$ for: A Poisson network with average degree $c=12$ and $N_{[0]}=1600$
nodes (panel (a)), a scale-free network with minimum degree $m=6$, power-law
exponent $\gamma=2.5$ and $N_{[0]}=1600$ nodes (panel (b)) and a fully
connected network of different network sizes $N_{[0]}=500$ (cyan symbols),
$N_{[0]}=1000$ (green symbols)and $N_{[0]}=2000$ (purple symbols) (panel (c)).
From panel (c) it is evident that the forward transition occurs at higher
values of $\sigma$ for larger network size, confirming the theoretical
prediction indicating that the transition is driven by finite size effects and
it is absent in the infinite network limit. Figures adapted from
ghorbanchian2020higher .

Model 2 of explosive higher-order synchronization of coupled topological
signals of nodes and links can be also solved exactly on a fully connected
network. Before discussing these theoretical results let us highlight that
when treating Model 2 on a fully connected network the model parameters need
to be rescaled appropriately to give a well defined transition in the large
network limit. In particular the coupling constant $\sigma$ and $\tau_{1}$ are
rescaled to

$\displaystyle\sigma$ $\displaystyle\to$ $\displaystyle\frac{\sigma}{N},$
$\displaystyle\tau_{1}$ $\displaystyle\to$ $\displaystyle\tau_{1}\sqrt{N}.$
(67)

Moreover, for simplicity we set $\Omega_{1}=0$. With these hypotheses the
marginal distribution $G_{i}(\tilde{\omega})=G(\tilde{\omega})$ for every node
$i$ of the network can be derived to be equal to

$\displaystyle G(\tilde{\omega})$ $\displaystyle=$
$\displaystyle\frac{\tau_{1}}{\sqrt{2\pi/\bar{c}}}\exp\left[-\tau_{1}^{2}\bar{c}\frac{\tilde{\omega}^{2}}{2}\right],$
(68)

with $\bar{c}=N/(N-1)$. The self-consistent equations for the order parameters
$R_{0}=\hat{R}_{0}$ and $R_{1}^{[-]}=\hat{R}_{1}^{[-]}$ are given by

$\displaystyle 1=\sigma
R_{1}^{[-]}h\left(\sigma^{2}R_{0}^{2}(R_{1}^{{[-]}})^{2}\right),$
$\displaystyle R_{1}^{[-]}=\sigma
R_{0}\tau_{1}\sqrt{\bar{c}}h\left(\sigma^{2}\tau_{1}^{2}R_{0}^{2}\right),$
(69)

where the scaling function $h(x)$ is given by

$\displaystyle
h(x)=\sqrt{\frac{\pi}{2}}e^{-x/4}\left[I_{0}\left(\frac{x}{4}\right)+I_{1}\left(\frac{x}{4}\right)\right],$
(70)

with $I_{0}$ and $I_{1}$ indicating the modified Bessel functions.

These equations agree perfectly well with direct simulation of Model 2
dynamical Eqs. (56) on a fully connected network as it can be appreciated from
Figure 10(c). Moreover a closer look to these equations reveals an important
aspect of these transitions. While the backward transition has a well defined
limit as $N\to\infty$, the forward transition occurs at larger value of the
coupling constant for large network size $N$, and is only determined by finite
size fluctuations, therefore the transition disappears in the limit
$N\to\infty$. Interestingly, this lack of a proper hysteresis loop can be also
predicted for Model 2 defined on uncorrelated random graphs with finite second
moment of the degree distribution, starting from their annealed approximation
solution.

##  7 Conclusions

In this chapter our goal has been to provide evidence that the interplay
between simplicial complex structure and dynamics is mediated by simplicial
geometry and topology. The spectral properties of the graph Laplacian and the
higher-order Laplacian have been used here to reveal how simplicial
synchronization is shaped by topological and geometry of the simplicial
complex. In particular, we investigated how simplicial network geometry
changes the phase diagram of the Kuramoto model defined on the network
skeleton of simplicial complexes with notable geometrical properties and
characterized by a finite spectral dimension. We have shown that a spectral
dimension smaller or equal than four but larger than two can lead to a regime
of frustrated synchronization characterized by large spatio-temporal
fluctuations of the order parameter, while a spectral dimension smaller or
equal than two leads to a Kuramoto model in the incoherent state for every
finite value of the coupling constant. These theoretical results have been
shown to apply to the simplicial complexes generated by the modelling
framework called Network Geometry with Flavor (NGF) that is able to generate
simplicial complexes with tunable spectral dimension of the graph Laplacian.
Interestingly, the NGF are characterized also by displaying higher-order
spectral dimension of the higher-order up Laplacian that describe higher-order
diffusion processes.

This chapter introduces also a set of models for capturing synchronization of
topological signals, i.e., phases not only associated with the nodes of a
simplicial complex but also to the higher-order simplices such as links,
triangles, and so on. This higher-order synchronization reveals itself in the
order parameter of the irrotational and solenoidal projection of the
topological signals. In the simple higher-order Kuramoto model the
irrotational and solenoidal projection of the topological signal are uncoupled
and undergo a sychronization transition at $\sigma_{c}=0$. However, when these
two projections are coupled to each other by an adaptive global coupling the
synchronization becomes explosive, i.e., discontinuous, and occurs at a non-
zero value of coupling constant.

The higher-order Kuramoto model can be further extended to capture coupled
topological signals of different dimension, for instance coupling phases
associated with nodes and to links of a network or of a simplicial complex.
This generalized higher-order Kuramoto model can lead to an explosive phase
transition affecting simultaneously the phases associated with the nodes and
the irrotational and solenoidal projection of the phases associated with the
links. Interestingly, the higher-order Kuramoto model of coupled topological
signals defined on nodes and links can be treated analytically using the
annealed approximation when it is defined on a random uncorrelated network and
can be solved exactly on a fully connected network. This solution confirms the
discontinuous nature of the transition of the explosive higher-order Kuramoto
model and sheds light on the stability of the hysteresis loop associated with
the transition on finite networks. The mathematical framework that we have
proposed here can be explored and modified in different directions and we
believe that an in-depth analysis of the model and its variations will provide
important insights on the interplay between topology and dynamics. For
instance, we note that the higher-order Kuramoto model has been recently
modified lee  to investigate also the properties of a consensus model finding
interesting results.

In conclusion, this chapter aims to provide an overview of the relation
between network geometry topology and dynamics. We believe this topic will
flourish in the incoming years and will transform our understanding of the
relation between structure and dynamics of higher-order networks. Therefore
our expectation is that this research line will play a relevant role for
providing new insights in a variety of applications including brain dynamics
and biological transportation networks.

## References

  * (1) Albert László Barabási.  Network Science.  Cambridge University press, 2016.
  * (2) Sergey N Dorogovtsev, Alexander V Goltsev, and José FF Mendes.  Critical phenomena in complex networks.  Reviews of Modern Physics, 80(4):1275, 2008.
  * (3) Alain Barrat, Marc Barthelemy, and Alessandro Vespignani.  Dynamical processes on complex networks.  Cambridge University Press, 2008.
  * (4) Federico Battiston, Giulia Cencetti, Iacopo Iacopini, Vito Latora, Maxime Lucas, Alice Patania, Jean-Gabriel Young, and Giovanni Petri.  Networks beyond pairwise interactions: structure and dynamics.  Physics Reports, 2020.
  * (5) Ginestra Bianconi.  Higher-order networks: An introduction to simplicial complexes.  Cambridge University Press, 2021.
  * (6) Ginestra Bianconi.  Interdisciplinary and physics challenges of network theory.  EPL (Europhysics Letters), 111(5):56001, 2015.
  * (7) Leo Torres, Ann S Blevins, Danielle S Bassett, and Tina Eliassi-Rad.  The why, how, and when of representations for complex systems.  arXiv preprint arXiv:2006.02870, 2020.
  * (8) Steven H Strogatz.  From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 143(1-4):1–20, 2000.
  * (9) Alex Arenas, Albert Díaz-Guilera, Jurgen Kurths, Yamir Moreno, and Changsong Zhou.  Synchronization in complex networks.  Physics Reports, 469(3):93–153, 2008.
  * (10) Stefano Boccaletti, Alexander N Pisarchik, Charo I Del Genio, and Andreas Amann.  Synchronization: from coupled systems to complex networks.  Cambridge University Press, 2018.
  * (11) Yoshiki Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In Huzihiro Araki, editor, International Symposium on Mathematical Problems in Theoretical Physics, pages 420–422, Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.
  * (12) Ana P Millán, Joaquín J Torres, and Ginestra Bianconi.  Explosive higher-order kuramoto dynamics on simplicial complexes.  Physical Review Letters, 124(21):218301, 2020.
  * (13) Joaquín J Torres and Ginestra Bianconi.  Simplicial complexes: higher-order spectral dimension and dynamics.  Journal of Physics: Complexity, 1(1):015002, 2020.
  * (14) Sergio Barbarossa and Stefania Sardellitti.  Topological signal processing over simplicial complexes.  IEEE Transactions on Signal Processing, 68:2992–3007, 2020.
  * (15) Owen T Courtney and Ginestra Bianconi.  Generalized network structures: the configuration model and the canonical ensemble of simplicial complexes.  Physical Review E, 93(6):062311, 2016.
  * (16) Ginestra Bianconi and Christoph Rahmede.  Network geometry with flavor: from complexity to quantum geometry.  Physical Review E, 93(3):032315, 2016.
  * (17) Ginestra Bianconi and Christoph Rahmede.  Emergent hyperbolic network geometry.  Scientific Reports, 7:41974, 2017.
  * (18) Daan Mulder and Ginestra Bianconi.  Network geometry and complexity.  Journal of Statistical Physics, 173(3-4):783–805, 2018.
  * (19) Owen T. Courtney and Ginestra Bianconi.  Weighted growing simplicial complexes.  Physical Review E, 95(6):062301, 2017.
  * (20) Per Sebastian Skardal and Alex Arenas.  Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes.  Physical review letters, 122(24):248301, 2019.
  * (21) Per Sebastian Skardal and Alex Arenas.  Higher-order interactions in complex networks of phase oscillators promote abrupt synchronization switching.  arXiv preprint arXiv:1909.08057, 2019.
  * (22) Per Sebastian Skardal and Alex Arenas.  Memory selection and information switching in oscillator networks with higher-order interactions.  Journal of Physics: Complexity, 2020.
  * (23) Repository for higher-order network codes.
  * (24) Mikhael Gromov.  Hyperbolic groups.  In Essays in group theory, pages 75–263. Springer, 1987.
  * (25) Ana P Millán, Reza Ghorbanchian, Nicolò Defenu, Federico Battiston, and Ginestra Bianconi.  Local topological moves determine global diffusion properties of hyperbolic higher-order networks.  arXiv preprint arXiv:2102.12885, 2021.
  * (26) Raffaella Burioni and Davide Cassi.  Random walks on graphs: ideas, techniques and results.  Journal of Physics A: Mathematical and General, 38(8):R45, 2005\.
  * (27) Rammal Rammal and Gérard Toulouse.  Random walks on fractal structures and percolation clusters.  Journal de Physique Lettres, 44(1):13–22, 1983.
  * (28) Raffaella Burioni and Davide Cassi.  Universal properties of spectral dimension.  Physical Review Letters, 76(7):1091, 1996.
  * (29) Zhihao Wu, Giulia Menichetti, Christoph Rahmede, and Ginestra Bianconi.  Emergent complex network geometry.  Scientific Reports, 5:10073, 2015.
  * (30) Ana P Millán, Giacomo Gori, Federico Battiston, Tilman Enss, and Nicolò Defenu.  Complex networks with tuneable dimensions as a universality playground.  arXiv preprint arXiv:2006.10421, 2020.
  * (31) Marija Mitrović Dankulov, Bosiljka Tadić, and Roderick Melnik.  Spectral properties of hyperbolic nanonetworks with tunable aggregation of simplexes.  Physical Review E, 100(1):012309, 2019.
  * (32) Ana P Millán, Joaquín J Torres, and Ginestra Bianconi.  Synchronization in network geometries with finite spectral dimension.  Phys. Rev. E, 99(2):022307, 2019.
  * (33) A P Millán, J J Torres, and G Bianconi.  Complex network geometry and frustrated synchronization.  Scientific Reports, 8(9910), 2018.
  * (34) Abubakr Muhammad and Magnus Egerstedt.  Control using higher order laplacians in network topologies.  In Proc. of 17th International Symposium on Mathematical Theory of Networks and Systems, pages 1024–1038. Citeseer, 2006.
  * (35) Timothy E Goldberg.  Combinatorial laplacians of simplicial complexes.  Senior Thesis, Bard College, 2002.
  * (36) Joaquín J Torres and J Marro.  Brain performance versus phase transitions.  Scientific Reports, 5:12216, 2015.
  * (37) Marcus Reitz and Ginestra Bianconi.  The higher-order spectrum of simplicial complexes: a renormalization group approach.  Journal of Physics A: Mathematical and Theoretical, 2020.
  * (38) Juan A Acebrón, Luis L Bonilla, Conrad J Pérez Vicente, Félix Ritort, and Renato Spigler.  The kuramoto model: a simple paradigm for synchronization phenomena.  Reviews of Modern Physics, 77(1):137, 2005.
  * (39) Juan G Restrepo, Edward Ott, and Brian R Hunt.  Onset of synchronization in large networks of coupled oscillators.  Physical Review E, 71(3):036151, 2005.
  * (40) M Chavez, D-U Hwang, A Amann, HGE Hentschel, and S Boccaletti.  Synchronization is enhanced in weighted complex networks.  Physical Review Letters, 94(21):218701, 2005.
  * (41) Weiyu Huang, Thomas AW Bolton, John D Medaglia, Danielle S Bassett, Alejandro Ribeiro, and Dimitri Van De Ville.  A graph signal processing perspective on functional brain imaging.  Proceedings of the IEEE, 106(5):868–885, 2018.
  * (42) Miguel Ruiz-Garcia and Eleni Katifori.  Topologically controlled emergent dynamics in flow networks.  arXiv preprint arXiv:2001.01811, 2020.
  * (43) H Hong, Hyunggyu Park, and M Y Choi.  Collective synchronization in spatially extended systems of coupled oscillators with random frequencies.  Physical Review E, 72(3):036217, 2005.
  * (44) Hyunsuk Hong, Hugues Chaté, Hyunggyu Park, and Lei Han Tang.  Entrainment transition in populations of random frequency oscillators.  Physical Review Letters, 99(18):184101, 2007.
  * (45) Paolo Moretti and Miguel A Muñoz.  Griffiths phases and the stretching of criticality in brain networks.  Nature Communications, 4:2521, 2013.
  * (46) Pablo Villegas, Paolo Moretti, and Miguel A Muñoz.  Frustrated hierarchical synchronization and emergent complexity in the human connectome network.  Scientific Reports, 4:5990, 2014.
  * (47) Reza Ghorbanchian, Juan G Restrepo, Joaquín J Torres, and Ginestra Bianconi.  Higher-order simplicial synchronization of coupled topological signals.  arXiv preprint arXiv:2011.00897, 2020.
  * (48) Jason W Rocks, Andrea J Liu, and Eleni Katifori.  Revealing structure-function relationships in functional flow networks via persistent homology.  Physical Review Research, 2(3):033234, 2020.
  * (49) Xiyun Zhang, Stefano Boccaletti, Shuguang Guan, and Zonghua Liu.  Explosive synchronization in adaptive and multilayer networks.  Physical review letters, 114(3):038701, 2015.
  * (50) Lav R Varshney, Beth L Chen, Eric Paniagua, David H Hall, and Dmitri B Chklovskii.  Structural properties of the caenorhabditis elegans neuronal network.  PLoS Comput Biol, 7(2):e1001066, 2011.
  * (51) Lee DeVille.  Consensus on simplicial complexes, or: The nonlinear simplicial laplacian.  arXiv preprint arXiv:2010.07421, 2020.

## Appendix: Kuramoto dynamics expressed in terms of the incidence matrix

In this Appendix our aim is to show that Eq. (35) that we rewrite here for
convenience,

$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sigma\sum_{j=1}^{N_{[0]}}a_{ij}\sin\left(\theta_{j}-\theta_{i}\right),$
(71)

is equivalent to Eq. (42) given by

$\displaystyle\dot{\bm{\theta}}=\bm{\omega}-\sigma\bm{B}_{[1]}\sin\bm{B}_{[1]}^{T}\bm{\theta}.$
(72)

In order to show this let us observe that the incidence matrix ${\bf B}_{[1]}$
has elements given by

$\displaystyle[B_{[1]}]_{i\ell}=\left\\{\begin{array}[]{cc}-1&\mbox{if}\
\ell=[i,j],\\\ 1&\mbox{if}\ \ell=[j,i],\\\
0&\mbox{otherwise}.\end{array}\right.$ (76)

To show the equivalence between Eq. (71) and Eq. (72) let us start by
rewriting Eq. (72) element by element, getting

$\displaystyle\dot{\theta}_{i}=\omega_{i}-\sigma\sum_{\ell\in
S_{1}}[B_{[1]}]_{i\ell}\sin\left(\sum_{j\in S_{0}}[B_{[1]}^{\top}]_{\ell
j}\theta_{j}\right),$ (77)

where we indicate with $S_{1}$ the set of all links present in the simplicial
complex or network under consideration.

Let us consider the particular link $\ell=[i,j]$ in this case we have

$\displaystyle[B_{[1]}]_{i\ell}\sin\left(\sum_{j\in
S_{0}}[B_{[1]}^{\top}]_{\ell
j}\theta_{j}\right)=-a_{ij}\sin(\theta_{j}-\theta_{i}).$ (78)

Equivalently, if we consider the same link with opposite orientation
$\ell=[j,i]$ we get

$\displaystyle[B_{[1]}]_{i\ell}\sin\left(\sum_{j\in
S_{0}}[B_{[1]}^{\top}]_{\ell
j}\theta_{j}\right)=a_{ij}\sin(\theta_{i}-\theta_{j})=-a_{ij}\sin(\theta_{j}-\theta_{i}).$
(79)

Since the incidence matrix ${\bf B}_{[1]}$ has non-zero elements only among
nodes and links incident to each other, it follows that Eq. (71) is equivalent
to Eq. (72).