# The three way Dirac operator and
dynamical Turing and Dirac induced patterns on nodes and links

Riccardo Muolo  Department of Systems and Control Engineering, Tokyo Institute
of Technology,
2 Chome-12-1 Ookayama, Meguro-ku, Tokyo, 152-8552, Japan  Department of
Mathematics & naXys, Namur Institute for Complex Systems, Université de Namur,
Rue Grafé 2, Namur, 5000, Belgium  Timoteo Carletti  Department of Mathematics
& naXys, Namur Institute for Complex Systems, Université de Namur,
Rue Grafé 2, Namur, 5000, Belgium  Ginestra Bianconi
<EMAIL_ADDRESS>Centre for Complexity Science, School of
Mathematical Sciences, Queen Mary University of London,
Mile End Road, London, E1 4NS, UK  The Alan Turing Institute, The British
Library, 96 Euston Road, London, NW1 2DB, UK

###### Abstract

Topological signals are dynamical variables not only defined on nodes but also
on links of a network that are gaining significant attention in non-linear
dynamics and topology and have important applications in brain dynamics. Here
we show that topological signals on nodes and links of a network can generate
dynamical patterns when coupled together. In particular, dynamical patterns
require at least three topological signals, here taken to be two node signals
and one link signal. In order to couple these signals, we formulate the 3-way
topological Dirac operator that generalizes previous definitions of the 2-way
and 4-way topological Dirac operators. We characterize the spectral properties
of the 3-way Dirac operator and we investigate the dynamical properties of the
resulting Turing and Dirac induced patterns. Here we emphasize the distinct
dynamical properties of the Dirac induced patterns which involve topological
signals only coupled by the 3-way topological Dirac operator in absence of the
Hodge-Laplacian coupling. While the observed Turing patterns generalize the
Turing patterns typically investigated on networks, the Dirac induced patterns
have no equivalence within the framework of node based Turing patterns. These
results open new scenarios in the study of Turing patterns with possible
application to neuroscience and more generally to the study of emergent
patterns in complex systems.

###### keywords:

Higher-order networks , Turing patterns , Topological signals , Dirac operator

##  1 Introduction

Networks [1, 2] represent the discrete architecture of complex systems by
encoding the set of interactions (links) existing among their constituent
elements (nodes). Network science is based on the fundamental assumption that
the network topology encodes important information about the network dynamics
and hence about the function of the underlying complex interacting system. The
statistical and combinatorial properties of the network structure have been
shown to be key in shaping the phase diagram of dynamical processes ranging
from epidemics to percolation [3, 4]. Yet, we are still far from a complete
understanding of the interplay between network topology and dynamics that
would be necessary, for instance, to transform our understanding of brain
dynamics.

Until recently, the dynamical state of a network has been exclusively
described following a node centered approach, where the dynamical variables
are exclusively defined on the nodes of the networks. Nowadays, it is
increasingly recognized that this approach has important limitations and that
in a number of cases the description of network dynamics should include
higher-order topological signals, i.e., dynamical variables not only defined
on nodes but also on links. Examples of link signals include synaptic signals
between neurons or brain link signals at the level of brain regions [5, 6] and
in general currents in biological transportation networks. Moreover, on
higher-order networks one can define topological signals also on triangles,
tetrahedra, and so on. This shift of perspective significantly enhances the
topological characterization of their structure[7, 8, 9, 10, 11, 12] and the
node-centric investigation of their dynamics [13, 14, 15, 16, 17].

The study of topological signals requires to combine topology with non-linear
dynamics and it is emerging as a field that can transform our understanding of
the interplay between structure and dynamics on simple and higher-order
networks [18, 19]. Interestingly, it has been found that topological signals
can undergo collective critical phenomena such as topological higher-order
Kuramoto model [20, 21, 22, 23], global topological synchronization [24] and
higher-order diffusion [25, 26, 27]. Moreover, real-world topological signals
can be predicted and processed with topological machine learning algorithms
[28, 29, 30, 31, 32].

On a simplicial complex, topological signals are modeled, treated and
processed by combining non-linear dynamics with algebraic topology and
discrete calculus inherited by the simplicial support (which reduces to a
network if the simplicial complex is $1$-dimensional). In particular,
topological signals of a given dimension are often treated by using the Hodge
Laplacian [18, 33] that describes diffusion from $n$-dimensional simplices to
$n$-dimensional simplices.

Recently, the discrete topological Dirac operator [34] has been shown to be
the most suitable operator in situations in which topological signals of
different dimensions interact and cross-talk; indeed the Dirac operator allows
to project the signal on $n$-dimensional signal to one dimension up or down.
For instance in a network the Dirac operator can be used to project node
signals to link signals and vice versa.

The Dirac operator has been originally defined in non-commutative geometry and
in quantum graphs [35], but only recently it has been show to lead to
topological field theories [34, 36, 37] and to be a key operator in complex
systems. Indeed, it is key to treat dynamics of coupled topological signals
[38, 39, 40, 41, 28, 42] on top of having important application to topological
analysis of real data [43, 44].

In this work, we combine topology and non-linear dynamics to study dynamical
pattern formation of topological signals of nodes and links, induced by the
Dirac operator.

On continuous domains, reaction-diffusion equations are partial differential
equations that can be used to describe chemical reaction systems, ecosystems,
neuronal dynamics and fluid-dynamics, and are a core subject of non-linear
dynamics [45]. Turing patterns are spatial patterns emerging from a diffusion-
driven instability of the homogeneous equilibrium state of a reaction-
diffusion system and are pivotal to describe pattern formation in the natural
world. In a nutshell, a system of two interacting species is perturbed about
its spatially homogeneous stable equilibrium, which becomes unstable due to
diffusion, giving rise to the celebrated Turing instability from which
patterns may originate [46]. The original framework in which Turing conceived
this pattern formation mechanism was morphogenesis, but, nowadays, it finds
applications in many different field, including biology [47], neuroscience
[48] and even quantum mechanics [49] and nano-materials [50]. Turing theory
has been extended on regular lattices in the 70s by Othmer and Scriven [51,
52], but it is only recently that increasing attention has been devoted to
characterize Turing patterns on networks, starting from the seminal work of
Nakao and Mikhailov [53]. Turing patterns have been thoroughly studied on
different networks topologies, such as directed [54], non-normal [55] and
geometric [56] ones, with applications ranging from ecology [57] to control
[58, 59, 60], to name a few. Moreover, they have also been observed on
multilayer [61, 62, 63], temporal [64, 65, 66] and higher-order [67, 68, 69]
networks.

In all these approaches it is assumed that the dynamics is only localized on
the nodes of the network. In order to address this limitation, recently
Giambagli et al. [40] have characterized the formation of topological Turing
patterns on nodes and links of networks. In this setting, the dynamical state
of the network is encoded by one topological signal (species) on the nodes and
one topological signal (species) on the links of the network coupled together
by the Dirac operator. The topological Turing patterns that emerge in this
scenario are inhomogeneous on nodes and links, however they are stationary,
i.e., they characterize a heterogeneous static asymptotic state of the
network.

In this work, we build a mathematical framework that is able to generate
dynamical Turing patterns on nodes and links. The key element of the model is
the assumption that the system state is encoded in three different types of
topological signals (species): two anchored to nodes and one on links. In
order to treat these topological signals, here we define the 3-way Dirac
operator and its associated gamma matrix. The Dirac operator projects the two
nodes signals into the links and the single link signal into the nodes, while
the gamma matrix is used to compress the projected nodes signals into a single
link signal and to expand the single projected link signal into two nodes
signals. The 3-way Dirac operator couples nodes and links topological signals
and is shown here to induce dynamical Turing patterns of topological nodes and
links signals displaying very rich and non-trivial dynamical properties.

##  2 Introduction to discrete exterior calculus and to the topological Dirac
operator

We consider a network $G=(V,E)$ comprising a set $V$ of $N_{0}=|V|$ nodes, and
a set $E$ of $N_{1}=|E|$ links. According to the usual algebraic topology
setting links are undirected, yet they are characterized by an orientation
from one endnode to the other. Assuming that the network dynamics is captured
by a topological signal on the nodes and a topological signal on the links,
the full dynamical state of the network is encoded by the topological spinor
$\bm{\Phi}\in C^{0}\oplus C^{1}$ given by

$\displaystyle{\mathbf{\Phi}}=\binom{u}{w}\,,$ (1)

where $u\in C^{1}$ is $0$-cochain here encoded by a vector taking real values
on each node of the network and $v\in C^{1}$ is a $1$-cochain here encoded by
a vector taking real values on each link of the network, i.e.

$\displaystyle u=\left(\begin{array}[]{c}u_{1}\\\ u_{2}\\\ \vdots\\\
u_{N_{0}}\end{array}\right)\text{ and }w=\left(\begin{array}[]{c}w_{1}\\\
w_{2}\\\ \vdots\\\ w_{N_{1}}\end{array}\right).$ (10)

The discrete exterior calculus is a branch of mathematics that allows to
define the discrete gradient and the discrete divergence of these topological
signals. In particular the discrete gradient $\delta_{1}:C^{0}\to C^{1}$ is
the linear operator that acts on topological nodes signals (e.g., the signal
$u$) and provides a link topological signal (e.g., $g$ with $g=\delta_{1}u$);
in the considered example we have

$\displaystyle g_{rs}=u_{s}-u_{r}\,,\text{for all link $[r,s]$.}$ (11)

Moreover the discrete divergence $\delta_{1}^{\star}:C^{1}\to C^{2}$ is the
linear operator that acts on link topological signals (e.g., $w$) and provides
a node topological signal (e.g., $f$ with $f=\delta_{1}^{\star}w$) where in
the considered case we have for all node $r$

$\displaystyle f_{r}=\sum_{\ell\in E^{+}_{r}}w_{\ell}-\sum_{\ell\in
E^{-}_{r}}w_{\ell}\,,$ (12)

where $E^{+}_{r}$ is the set of links oriented toward node $r$ while
$E^{-}_{r}$ is the set of links oriented from node $r$ toward their other
endnode. Therefore the discrete divergence associates to a node $r$ the
difference between the outward and inward flow from/to the node. The two
linear operators can be represented by the $N_{0}\times N_{1}$ boundary matrix
${\bf B}$ whose elements are given by

$\displaystyle{\bf B}_{r\ell}=\begin{cases}-1~{}\mbox{ if }~{}\ell=[r,s],\\\
~{}1~{}\mbox{ if }~{}\ell=[s,r].\end{cases}$

In particular we have ${f}={\bf B}w$ and $g={\bf B}^{\top}u$. On a network the
Hodge Laplacians ${\bf L}_{[0]}$ and ${\bf L}_{[1]}$ describe diffusion from
nodes to nodes passing through links and from links to links passing through
nodes. Therefore the Hodge Laplacians of a network can be build by contracting
the boundary matrix in the two possible ways, leading to ${\bf L}_{[0]}={\bf
B}{\bf B}^{\top}$ and ${\bf L}_{[1]}={\bf B}^{\top}{\bf B}$. The most simple
definition of the Hodge-Dirac operator [35, 34] on a network is the operator
$\partial:C^{0}\oplus C^{1}\to C^{0}\oplus C^{1}$ that maps topological
spinors to topological spinors and is defined as
$\partial=\delta_{1}^{\star}\oplus\delta_{1}$. In particular if
$\bm{\Phi}=u\oplus w$ we have

$\displaystyle\partial(u\oplus w)=\delta_{1}^{\star}w\oplus\delta_{1}u\,,$
(13)

or in matrix form we have

$\displaystyle\bm{\partial}=\left(\begin{array}[]{cc}0&\mathbf{B}\\\
\mathbf{B}^{\top}&0\end{array}\right)\,.$ (16)

The Hodge-Dirac operator can be considered as the “square root” of the higher-
order Laplacian $\boldsymbol{\mathcal{L}}$. Indeed we have

$\displaystyle\bm{\partial}^{2}=\boldsymbol{\mathcal{L}}=\left(\begin{array}[]{cc}\mathbf{L}_{[0]}&0\\\
0&\mathbf{L}_{[1]}\end{array}\right)\,.$ (19)

Therefore the non-zero eigenvalues of the Dirac operator are given by plus or
minus the square root of the non-zero eigenvalues of ${\bf L}_{[0]}$ (which is
by the way isospectral to ${\bf L}_{[1]}$ on a generic network). The Hodge-
Dirac operator admits eigenvectors that are either harmonic (i.e., associated
to the zero eigenvalue of the Hodge-Dirac operator) or chiral. Here by chiral
eigenvectors we refer to the relations of eigenvectors associated with non-
zero eigenvalues with the same absolute value. Indeed it can be easily proved
that $\bm{\partial}$ anticommutes (i.e.,
$\\{\bm{\partial},\bm{\gamma}_{0}\\}=0$) with $\bm{\gamma}_{0}$, where

$\displaystyle\bm{\gamma}_{0}=\left(\begin{array}[]{cc}{\bf I}_{N_{0}}&{0}\\\
{0}&-{\bf I}_{N_{1}}\end{array}\right)\,,$ (22)

where ${\bf I}_{X}$ indicates the $X\times X$ identity matrix. This result
implies that if $(u,w)^{\top}$ is an eigenvector of $\bm{\partial}$ with
eigenvalue $\hat{\lambda}$, then $(u,-w)^{\top}$ is an eigenvector of
$\bm{\partial}$ with eigenvalue $-\hat{\lambda}$. This relation between the
non-harmonic eigenvectors of the Hodge-Dirac operator is called chirality.
This Hodge-Dirac operator has been used in Ref.[40] to model Turing patterns
on nodes and links. However, by using a topological spinor formed by a single
topological nodes signal and a single topological links signal, can only
account for stationary Turing patterns.

Interestingly the Hodge-Dirac operator can be coupled with group operations
enforced by the so called gamma matrices [34]. For instance on a lattice the
gamma matrices allows to distinguish between gradient and divergence performed
along different directions. In particular on a three dimensional lattice in
order to distinguish between the $x$, $y$, and $z$ directions one needs to
consider two topological signals on the nodes and two topological signals on
the links and express the gamma matrices in terms of the Pauli matrices acting
on the two dimensional node (or link) signal. However no existing approach is
able to jointly treat and process an odd number of topological signals such as
two signal on the nodes and one signal on the links or viceversa two signals
on the links and on signal on the nodes.

Figure 1: Schematic representation of the dynamical state of the network for
dynamical Turing patterns. We consider a network $G=(V,E)$ whose dynamical
states includes two node topological signals (species densities represented
with yellow triangles and yellow stars respectively on the enlarged node,
although they are defined on each node) and one link topological signal
schematically represent with green arrows associated to links.

##  3 The 3-way Dirac operator: 2 species on the nodes and 1 species on the
links

In this work we consider a network dynamical state captured by two topological
signals on the nodes and one topological signal on the links encoded by the
topological spinor ${\Phi}\in C^{0}\oplus C^{0}\oplus C^{1}$ (see Fig.1). Note
that this approach can be readily generalized to the case in which the
dynamical state of the network is encoded into a single topological node
signal and two topological link signals.

The topological spinor encodes the 3-way topological signals of the network
and has block structure

$\displaystyle{\mathbf{\Phi}}=\left(\begin{array}[]{c}\bm{\chi}\\\
\bm{\psi}\end{array}\right)\,,$ (25)

where $\bm{\chi}\in C^{0}\oplus C^{0}$ is defined on nodes and $\bm{\psi}\in
C^{1}$ is defined on links. In particular we will use the notation

$\displaystyle\bm{\chi}=\left(\begin{array}[]{c}{u}\\\
{v}\end{array}\right),\quad\bm{\psi}={w}\,,$ (28)

with $u\in C^{0}$, $v\in C^{0}$ indicating the two node signals and $w\in
C^{1}$ indicating the link signal. For ease of notation here and in the
following we indicate with $\mathcal{N}=2N_{0}+N_{1}$ and with
$\mathcal{M}=N_{0}+2N_{1}$. We have therefore that $\bm{\Phi}$ is the
$\mathcal{N}$ dimensional column vector, given by

$\displaystyle\bm{\Phi}=\left(\begin{array}[]{c}{u}\\\ {v}\\\
{w}\end{array}\right)\,.$ (32)

The 3-way Hodge-Dirac operator $\bm{\partial}:C^{0}\oplus C^{0}\oplus C^{1}\to
C^{0}\oplus C^{1}\oplus C^{1}$ acts on the topological spinor $\Phi$
projecting independently the two node signals encoded by $\bm{\chi}$ into the
links and the single link signal $\bm{\psi}$ into the nodes. We define the
3-way Hodge-Dirac operator $\bm{\partial}$ as the
$\mathcal{M}\times\mathcal{N}$ matrix given by

$\displaystyle\bm{\partial}=\left(\begin{array}[]{ccc}{\bf 0}&{\bf
I}_{1}\otimes\mathbf{B}\\\ {\bf I}_{2}\otimes\mathbf{B}^{\top}&{\bf
0}\end{array}\right)\,,$ (35)

or alternatively as

$\displaystyle\bm{\partial}=\left(\begin{array}[]{ccc}0&0&\mathbf{B}\\\
\mathbf{B}^{\top}&0&0\\\ 0&\mathbf{B}^{\top}&0\end{array}\right)\,.$ (39)

In contrast with the Hodge-Dirac operator, the Dirac operator
${{\bm{\mathcal{D}}}\hskip-6.25958pt/}:C^{0}\oplus C^{0}\oplus C^{1}\to
C^{0}\oplus C^{0}\oplus C^{1}$ allows the node and link signals to cross talk
as it maps topological spinors into topological spinors. The Dirac operator is
obtained by multiplying the Hodge-Dirac operator by the
$\mathcal{N}\times\mathcal{M}$ gamma matrix $\bm{\gamma}$, i.e.

$\displaystyle{\bm{\mathcal{D}}}\hskip-6.25958pt/=\bm{\gamma}{\bm{\partial}}\,.$
(40)

Therefore the gamma matrix $\bm{\gamma}$ defines the way in which the two
nodes signals projected on the links are compressed and how the single link
signal projected on the nodes is expanded into the two nodes signals of the
topological spinor. The gamma matrix $\bm{\gamma}$ is defined in terms of
2-dimensional column vector ${\bm{\alpha}}=(\alpha_{u},\alpha_{v})^{\top}$ and
a 2-dimensional row vector ${\bm{\beta}}=(\beta_{u},\beta_{v})$, as

$\displaystyle\bm{\gamma}=\left(\begin{array}[]{ccc}{\bm{\alpha}}\otimes{\bf
I}_{N_{0}}&{\bf 0}\\\ {\bf 0}&{\bm{\beta}}\otimes{\bf
I}_{N_{1}}\end{array}\right)\,,$ (43)

or alternatively

$\displaystyle\bm{\gamma}=\left(\begin{array}[]{ccc}{\alpha}_{u}{\bf
I}_{N_{0}}&0&0\\\ {\alpha}_{v}{\bf I}_{N_{0}}&0&0\\\ {0}&{\beta}_{u}{\bf
I}_{N_{1}}&{\beta}_{v}{\bf I}_{N_{1}}\end{array}\right)\,.$ (47)

It follows that the 3-way Dirac operator can be also expressed as

${\bm{\mathcal{D}}}\hskip-6.25958pt/=\bm{\gamma}\bm{\partial}=\left(\begin{array}[]{ccc}0&0&\alpha_{u}\mathbf{B}\\\
0&0&\alpha_{v}\mathbf{B}\\\
\beta_{u}\mathbf{B}^{\top}&\beta_{v}\mathbf{B}^{\top}&0\end{array}\right)\,.$
(48)

As previously shown in Eq. (19) the Dirac operator defined on a topological
spinors having the same number of topological signals on nodes and links, can
be interpreted as the “square-root” of the higher-order Laplacian. Therefore
it is natural to investigate whether a similar condition holds true for 3-way
Dirac operator. One can straightforwardly obtain

${\bm{\mathcal{D}}}\hskip-6.25958pt/^{2}=\bm{\mathcal{L}}=\left(\begin{array}[]{ccc}D_{uu}\mathbf{L}_{[0]}&D_{uv}\mathbf{L}_{[0]}&0\\\
D_{vu}\mathbf{L}_{[0]}&D_{vv}\mathbf{L}_{[0]}&0\\\
0&0&D_{ww}\mathbf{L}_{[1]}\end{array}\right)\,,$ (49)

with
$D_{uu}=\alpha_{u}\beta_{u},D_{uv}=\alpha_{u}\beta_{v},D_{vu}=\alpha_{v}\beta_{u},D_{vv}=\alpha_{v}\beta_{v}$
and $D_{ww}=\alpha_{u}\beta_{u}+\alpha_{v}\beta_{v}$, assumed to be non-
negative. It is instructive to consider the case in which
$\bm{\alpha}=(1,1)^{\top}$ and $\bm{\beta}=(1,1)$. In this case we have
${\bm{\mathcal{D}}}\hskip-6.25958pt/={\bm{\mathcal{D}}}\hskip-6.25958pt/_{\
0}$ with

$\displaystyle{\bm{\mathcal{D}}}\hskip-6.25958pt/_{\
0}=\left(\begin{array}[]{ccc}0&0&\mathbf{B}\\\ 0&0&\mathbf{B}\\\
\mathbf{B}^{\top}&\mathbf{B}^{\top}&0\end{array}\right)\text{ and
}{\bm{\mathcal{D}}}\hskip-6.25958pt/_{\
0}^{2}=\bm{\mathcal{L}}_{0}=\left(\begin{array}[]{ccc}\mathbf{L}_{[0]}&\mathbf{L}_{[0]}&0\\\
\mathbf{L}_{[0]}&\mathbf{L}_{[0]}&0\\\
0&0&2\mathbf{L}_{[1]}\end{array}\right).$ (56)

On a side note we observe that the Dirac operator defined in Eq.(56) is
reminiscent of the topological operator used for studying duplex networks with
link overlap in Ref. [70].

Given the diagonal block structure of $\bm{\mathcal{L}}$, its eigenvalues are
either the eigenvalues of $\bm{\mathcal{L}}_{uv}$ or the ones of
$\bm{\mathcal{L}}_{w}$ with

$\displaystyle\bm{\mathcal{L}}_{uv}=\left(\begin{array}[]{cc}D_{uu}\mathbf{L}_{[0]}&D_{uv}\mathbf{L}_{[0]}\\\
D_{vu}\mathbf{L}_{[0]}&D_{vv}\mathbf{L}_{[0]}\end{array}\right),\quad\bm{\mathcal{L}}_{w}=D_{ww}\mathbf{L}_{[1]}\,.$
(59)

From this property it follows that the eigenvalues of the Dirac operator are
the square roots of the eigenvalues of $\mathcal{L}$ taken both with positive
and negative sign. Therefore, while the Hodge Laplacians are semi-definite
positive the Dirac operator is not semi-definite positive and, as we will
discuss in the following, might also display complex eigenvalues.

We observe that the three way Dirac operator
${\bm{\mathcal{D}}}\hskip-6.25958pt/~{}$ commutes with the matrix
$\bm{\gamma}_{0}$ given by

$\bm{\gamma}_{0}=\left(\begin{array}[]{ccc}\mathbf{I}_{N_{0}}&0&0\\\
0&\mathbf{I}_{N_{0}}&0\\\ 0&0&-\mathbf{I}_{N_{1}}\end{array}\right)\,.$ (60)

Hence if a spinor $\boldsymbol{\Phi}=(u,v,w)^{\top}$ is an eigenvector of the
three way Dirac operator with eigenvalue ${\hat{\lambda}}\neq 0$, i.e. if

${\bm{\mathcal{D}}}\hskip-6.25958pt/\bm{\Phi}={\hat{\lambda}}\bm{\Phi}\,,$
(61)

then $\bm{\gamma}_{0}\mathbf{\Phi}=(u,v,-w)^{\top}$ is an eigenvector of the
three way Dirac operator with eigenvalue $-{\hat{\lambda}}$. Indeed we have

${\bm{\mathcal{D}}}\hskip-6.25958pt/\bm{\gamma}_{0}\bm{\Phi}=-\bm{\gamma}_{0}{\bm{\mathcal{D}}}\hskip-6.25958pt/\bm{\Phi}=-{\hat{\lambda}}\bm{\gamma}_{0}\bm{\Phi}\,.$
(62)

From these results it follows that the harmonic eigenvectors of the three way
Dirac operator are

$\phi_{0}^{\text{harm }}=\left(\begin{array}[]{c}\bm{\chi}^{\text{harm}}\\\
0\end{array}\right)\quad\phi_{1}^{\text{harm }}=\left(\begin{array}[]{c}0\\\
{\bm{\psi}}^{\text{harm }}\\\ \end{array}\right)\,,$ (63)

where $\bm{\chi}^{\text{harm}}$ indicates the generic harmonic eigenvector of
$\bm{\mathcal{L}}_{uv}$ while ${\bm{\psi}}^{\text{harm }}$ the generic
harmonic eigenvector of $\bm{\mathcal{L}}_{w}$.

Note that the harmonic eigenvectors $\bm{\chi}^{\text{harm}}$ read

$\displaystyle\left(\begin{array}[]{c}\hat{\bm{\chi}}_{0}\\\ 0\\\
\end{array}\right),\quad\left(\begin{array}[]{c}0\\\
\hat{\bf\chi}_{0}\end{array}\right),\quad\left(\begin{array}[]{c}-\beta_{v}\hat{\bm{\chi}}\\\
\beta_{u}\hat{\bm{\chi}}\end{array}\right)\,.$ (70)

Where $\hat{\bm{\chi}}_{0}$ are the harmonic right singular vectors of the
boundary matrix, and where $\hat{\bm{\chi}}$ is an arbitrary $N_{0}$
dimensional vector. The non zero eigenvalues $\hat{\lambda}_{k}$ of the
$3$-way Dirac operator are given by

$\displaystyle\hat{\lambda}_{k}=\pm\sqrt{\alpha_{u}\beta_{u}+\alpha_{v}\beta_{v}}\mu_{k}$
(71)

where $\mu_{k}$ are the singular eigenvalues of the boundary operator. Note
that with a general choice of the gamma matrix $\bm{\gamma}$ the eigenvalues
of the Dirac operator are real only if

$\displaystyle\Delta=\alpha_{u}\beta_{u}+\alpha_{v}\beta_{v}>0,$ (72)

and becomes purely imaginary (i.e., appears in complex conjugate pairs) for
$\Delta<0$ while they are identically zero in the trivial case $\Delta=0$.

The eigenvectors $\boldsymbol{\phi}^{+}$, $\boldsymbol{\phi}^{-}$ associated
to the non-zero eigenvalues $\hat{\lambda}_{k}$ with the same absolute value
are related by chirality that for the three-way Dirac operator takes the form

$\phi^{+}=\mathcal{N}\left(\begin{array}[]{c}\alpha_{u}\hat{\bm{\chi}}_{\mu}\\\
\alpha_{v}\hat{\bm{\chi}}_{\mu}\\\
\sqrt{\Delta}{\bm{\psi}}_{\mu}\end{array}\right)\quad\phi^{-}=\mathcal{N}\left(\begin{array}[]{c}\alpha_{u}\hat{\bm{\chi}}_{\mu}\\\
\alpha_{v}\hat{\bm{\chi}}_{\mu}\\\
-\sqrt{\Delta}{\bm{\psi}}_{\mu}\end{array}\right)\,,$ (73)

where $\hat{\chi}_{\mu}$ are the right singular vectors and $\bm{\psi}_{\mu}$
are the left singular vectors corresponding to a non-zero singular value
$\mu=\mu_{k}$ of the boundary operator.

Let us now discuss the spectral properties of the Dirac operator
${\bm{\mathcal{D}}}\hskip-6.25958pt/={\bm{\mathcal{D}}}\hskip-6.25958pt/_{\
0}$ defined in Eq. (56) in the case in which $\bm{\alpha}=(1,1)^{\top}$ and
$\bm{\beta}=(1,1)$. As long as the network is connected the zero eigenvalue of
${\bm{\mathcal{L}}}$ has multiplicity $N_{1}+2$ while the non-zero eigenvalues
are given by twice the non-zero eigenvalues $\mu$ of $L_{[0]}$ (or
equivalently of $\mathbf{L}_{[1]}$ since on a network $\mathbf{L}_{[0]}$ is
isospectral to $\mathbf{L}_{[1]}$). Therefore the eigenvalues of the Dirac
operator ${\bm{\mathcal{D}}}\hskip-6.25958pt/_{\ 0}$ in this case are real and
given by

$\displaystyle\hat{\lambda}_{k}=0\quad\text{ with multiplicity }N_{1}+2,$ (74)
$\displaystyle\hat{\lambda}_{k}=\pm\sqrt{2}\mu_{k}\,.$

##  4 Theory of topological 3-way Turing patterns

Our goal will be now to define a dynamical system for the $3$-way topological
spinor defined in the previous section, that can display dynamical Turing
patterns on nodes and links. To this end, we assume that in absence of
interactions, the dynamical state of the network, captured by the topological
spinor $\bm{\Phi}$ follows the dynamics

$\displaystyle\dot{\bm{\Phi}}={\bf F}(\bm{\Phi},\bm{\partial}\bm{\Phi}),$ (75)

where ${\bf F}(\bm{\Phi},\bm{\partial}\bm{\Phi})$ determines the reaction
dynamics involving the three signals $u,v,w$. Let us observe that, as we will
see further in this section, the term $\bm{\partial}\bm{\Phi}$ vanishes when
there are no interactions between simplices of different dimensions.
Nonetheless, this formalism allows us to consider also nonlinear terms in
$\mathbf{B}$ and $\mathbf{B}^{\top}$, making the latter more general. Let us
now specify more in detail the assumed structure of this reaction term. As
previously stated the topological spinor
$\bm{\Phi}=(\bm{\chi},\bm{\psi})^{\top}$, where $\bm{\chi}=(u,v)^{\top}$,
encodes for the two node signals and $\bm{\psi}=w$ encodes for the single link
signal. On each node the signal $u$ and $v$ can interact, while they cannot
interact with the link signal. The projected topological spinor
$\bm{\partial}\bm{\Phi}$ is instead given by

$\displaystyle\bm{\partial}\bm{\Phi}=\left(\begin{array}[]{c}\bm{\hat{}}{\bm{\chi}}\\\
\hat{\bm{\psi}}\end{array}\right),$ (78)

with

$\displaystyle\hat{\bm{\chi}}={\bf
B}w,\quad\hat{\bm{\psi}}=\left(\begin{array}[]{c}{\bf B}^{\top}u\\\ {\bf
B}^{\top}v\end{array}\right).$ (81)

In other words $\hat{\bm{\chi}}$ is the one dimensional projected signal on
the nodes and $\hat{\bm{\psi}}$ is the two dimensional projected signals on
the links. Since on nodes only node signals can interact, and on links only
link signals can interact, it follows that on nodes all the signals encoded in
$\bm{\chi}$ can interact with each other and with the projected signals
$\hat{\bm{\chi}}$ and similarly on the links the projected signals on the
links $\hat{\bm{\psi}}$ can interact with each other and with the link signal
${\bm{\psi}}$.

Accordingly, our choice of the reaction term ${\bf
F}(\bm{\Phi},\bm{\partial}\bm{\Phi})$ is

$\displaystyle{\bf
F}(\bm{\Phi},\bm{\partial}\bm{\Phi})=\left(\begin{array}[]{c}f_{1}({\bm{\chi}},\hat{\bm{\chi}})\\\
f_{2}({\bm{\chi}},\hat{\bm{\chi}})\\\
f_{3}({\bm{\psi}},\hat{\bm{\psi}})\end{array}\right)=\left(\begin{array}[]{c}f(u,v,{\bf
B}w)\\\ g(u,v,{\bf B}w)\\\ h({\bf B}^{\top}u,{\bf
B}^{\top}v,w)\end{array}\right),$ (88)

where $f,g,h$ act on their arguments element-wise.

Here we want to show that topological dynamical Turing patterns emerge when we
introduce interactions between the topological signals of nodes and nearby
links captured by the Laplacian and Dirac operators. We therefore consider the
reaction-diffusion equations

$\dot{\bm{\Phi}}={\bf
F}(\bm{\Phi},\bm{\partial}\bm{\Phi})-c_{1}{\bm{\mathcal{D}}}\hskip-6.25958pt/\bm{\Phi}-c_{2}\bm{\mathcal{L}}\bm{\Phi},$
(89)

with $c_{1}\in\mathbb{R}$, $c_{2}\in\mathbb{R}^{+}$. Note that when $c_{1}=0$
the interactions will be only diffusive, i.e., driven by the Hodge Laplacians
while when $c_{2}=0$ the interaction is exclusively driven by the Dirac
operator. The possible patterns emerging in the former case will thus named
Turing patterns, while we will deal with Dirac patterns, in the latter
scenario. Let us observe that interestingly, this latter scenarios that we
introduce here, leads to very non trivial dynamical patterns as we will show
in the next sections. Note that these Dirac induced patterns have no
equivalent in the framework of node centered theory of Turing pattern
formation.

We can rewrite Eq. (89) explicitly, obtaining the following $3$-species
reaction-diffusion system

$\begin{cases}\dot{u}=f(u,v,\mathbf{B}w)-c_{2}\mathbf{L}_{[0]}(D_{uu}u+D_{uv}v)-c_{1}\alpha_{u}\mathbf{B}w,\\\
\dot{v}=g(u,v,\mathbf{B}w)-c_{2}\mathbf{L}_{[0]}(D_{vu}u+D_{vv}v)-c_{1}\alpha_{v}\mathbf{B}w,\\\
\dot{w}=h(\mathbf{B}^{\top}u,\mathbf{B}^{\top}v,w)-c_{2}D_{ww}\mathbf{L}_{[1]}-c_{1}\mathbf{B}^{\top}(\beta_{u}u+\beta_{v}v)\,.\end{cases}$
(90)

Let us observe that cross-diffusion emerges naturally from the formalism
developed in the previous section.

#### Remark

Let us observe that, if we take functions $f$, $g$ and $h$ such that the terms
in $\mathbf{B}$ and $\mathbf{B}^{\top}$ are linear (as it will be the case for
the example that we will hereafter present), we can map the problem into a new
one with a local reaction term ${\bf F}(\bm{\Phi})$ and a linear Dirac part.
The relevance of this interpretation will become clear in the next Sections,
when we will show that an oscillatory instability can emerge solely due to the
Dirac coupling.

In order to study the emergence of Turing and Dirac patterns, we need the
system to exhibit a stable homogeneous equilibrium and that turns unstable
once subjected to a spatially inhomogeneous perturbation. Before studying the
stability of the homogeneous state, let us recall that such equilibrium is a
solution for the interconnected system provided the constant eigenvector ${\bf
1}$ is in the kernel of the Dirac operator [40, 24]. For the case of networks,
this implies that the divergence of the constant link signal is null, i.e.,
that there is an orientation of the simplex such that for each node $r$ a
equal number of links point toward it and point outward it [40]. This relation
implies that the network must be Eulerian, i.e., to have only nodes with even
degree.

When the coupling is not active, we have
$\bm{\partial}\bm{\Phi}^{\star}=(0,0,0)$, hence, the system reduces to $N_{0}$
isolated systems of two interacting species on the nodes and $N_{1}$ systems
of one species on the links, of the following form

$\begin{cases}\dot{u}=f(u,v,0),\\\ \dot{v}=g(u,v,0),\\\
\dot{w}=h(0,0,w)\,,\end{cases}$ (91)

hence, we can study the local stability of this homogeneous state considering
nodes and links separately.
Let us assume $\bm{\Phi}^{*}=(u^{*},v^{*},w^{*})^{\top}$ to be a fixed point
for the above system, meaning that $w^{*}$ is a fixed point of $h(0,0,w)$ and
$(u^{*},v^{*})^{\top}$ for the system $f(u,v,0),~{}g(u,v,0)$. For the equation
on the links, the linear stability condition is trivially $h_{w}<0$, where
with $h_{w}$ we indicate the derivative of function $h$ with respect to
variable $w$ evaluated at the homogeneous state $w^{*}$. For what concerns the
dynamics on the nodes, it reduces to studying the stability of the following
$2$ dimensional system

$\begin{cases}\dot{u}=f(u,v,0)\\\ \dot{v}=g(u,v,0)\,.\end{cases}$

From the assumption that $(u^{*},v^{*})^{\top}$ is a fixed point for the above
system, we can linearize around it and obtain the Jacobian matrix, whose
stability conditions are $f_{u}+g_{v}<0$ and $f_{u}g_{v}-f_{v}g_{u}>0$. Hence,
the homogeneous vector $(u^{*},v^{*},w^{*})^{\top}$ is a stable equilibrium
for system (91) provided that

$\begin{cases}f_{u}+g_{v}<0\\\ f_{u}g_{v}-f_{v}g_{u}>0\\\ h_{w}<0\
.\end{cases}$ (92)

Let us now perturb system (90) with a spatially inhomogeneous perturbation
$(\delta u,\delta v,\delta w)^{\top}$ and linearize it around the homogeneous
equilibrium point that we have proved to be a solution also of the coupled
system. We hence obtain the following system for the dynamics of the
perturbation

$\begin{cases}\delta\dot{u}=f_{u}\delta u+f_{v}\delta v+f_{w}\mathbf{B}\delta
w-c_{2}\mathbf{L}_{[0]}(D_{uu}\delta u+D_{uv}\delta v)\\\
\delta\dot{v}=g_{u}\delta u+g_{v}\delta v+g_{w}\mathbf{B}\delta
w-c_{2}\mathbf{L}_{[0]}(D_{vu}\delta u+D_{vv}\delta v)\\\
\delta\dot{w}=h_{u}\mathbf{B}^{\top}\delta u+h_{v}\mathbf{B}^{\top}\delta
v+h_{w}\delta w-c_{2}D_{ww}\mathbf{L}_{[1]}\delta w\,,\end{cases}$ (93)

where, for sake of notation, we have considered the terms in $c_{1}$ of Eq.
(90) as part of the reaction functions, namely
$f_{w}=\partial_{\mathbf{B}w}f-c_{1}\alpha_{u}$,
$g_{w}=\partial_{\mathbf{B}w}g-c_{1}\alpha_{v}$,
$h_{u}=\partial_{\mathbf{B}^{\top}u}h-c_{1}\beta_{u}$ and
$h_{v}=\partial_{\mathbf{B}^{\top}v}h-c_{1}\beta_{v}$.

We can proceed as in [40] and use the eigenvectors of $\mathbf{L}_{[0]}$ and
$\mathbf{L}_{[1]}$ to perform the singular value decomposition of
$\mathbf{B}$. Let us recall that $\mathbf{L}_{[0]}$ and $\mathbf{L}_{[1]}$
are, in this case, isospectral and their eigenvalues $\mu_{k}^{2}$ are the
square of the singular values $\mu_{k}$ of $\mathbf{B}$. By projecting the
perturbations $\delta u$ and $\delta v$ on the eigenbasis of
$\mathbf{L}_{[0]}$, and the perturbation $\delta w$ on the eigenbasis of
$\mathbf{L}_{[1]}$, Eq. (93) becomes

$\begin{cases}\delta\dot{\hat{u}}_{k}=f_{u}\delta\hat{u}_{k}+f_{v}\delta\hat{v}_{k}+f_{w}\mu_{k}\delta\hat{w}_{k}-\mu_{k}^{2}c_{2}(D_{uu}\delta\hat{u}_{k}+D_{uv}\delta\hat{v}_{k}),\\\
\delta\dot{\hat{v}}_{k}=g_{u}\delta\hat{u}_{k}+g_{v}\delta\hat{v}_{k}+g_{w}\mu_{k}\delta\hat{w}_{k}-\mu_{k}^{2}c_{2}(D_{vu}\delta\hat{u}_{k}+D_{vv}\delta\hat{v}_{k}),\\\
\delta\dot{\hat{w}}_{k}=h_{u}\mu_{k}\delta\hat{u}_{k}+h_{v}\mu_{k}\delta\hat{v}_{k}+h_{w}\delta\hat{w}_{k}-\mu_{k}^{2}c_{2}D_{ww}\delta\hat{w}_{k},\end{cases}$
(94)

where $\delta\hat{u}_{k}$ (resp. $\delta\hat{v}_{k},\delta\hat{w}_{k}$) is the
perturbation expressed in the new basis. To study the linear stability of the
perturbation, we hypothesize that
$\delta\hat{u}_{k}(t),\delta\hat{v}_{k}(t),\delta\hat{w}_{k}(t)\sim
e^{\lambda_{k}t}$. To ensure the existence of a nontrivial solution, the
linear growth rate $\lambda_{k}$ must satisfy the following condition

$\left|J_{k}\right|=0\,,$ (95)

where $|J_{k}|$ denotes the determinant of the matrix

$J_{k}=\left(\begin{array}[]{ccc}f_{u}-\mu_{k}^{2}c_{2}D_{uu}-\lambda_{k}&f_{v}-\mu_{k}^{2}c_{2}D_{uv}&\mu_{k}f_{w}\\\
g_{u}-\mu_{k}^{2}c_{2}D_{vu}&g_{v}-\mu_{k}^{2}c_{2}D_{vv}-\lambda_{k}&\mu_{k}g_{w}\\\
\mu_{k}h_{u}&\mu_{k}h_{v}&h_{w}-\mu_{k}^{2}c_{2}D_{ww}-\lambda_{k}\end{array}\right)\,.$

We can rewrite this condition as a $3$-rd order polynomial in the variable
$\lambda_{k}$ of the form

$a\lambda_{k}^{3}+b\lambda_{k}^{2}+c\lambda_{k}+d=0,$ (96)

whose coefficients are

$\begin{cases}a=1\\\
b=c_{2}\mu_{k}^{2}(D_{uu}+D_{vv}+D_{ww})-(f_{u}+g_{v}+h_{w})\\\
c=A\mu_{k}^{4}+B\mu_{k}^{2}+C\equiv\pi_{1}(\mu_{k}^{2})\\\
d=\mathcal{A}\mu_{k}^{4}+\mathcal{B}\mu_{k}^{2}+\mathcal{C}\equiv\pi_{2}(\mu_{k}^{2})\,,\end{cases}$

where $\pi_{1}$ and $\pi_{2}$ are polynomials in $\mu_{k}^{2}$, whose explicit
expressions can be found in A. By means of the Routh-Hurwitz criterion [71,
72], we can study the stability of polynomial (96). When it is unstable, we
have Turing patterns. The explicit conditions are rather cumbersome, but it is
easy to check that Turing patterns can be obtained.

We will call $\lambda$ the maximum real part of the $\lambda_{k}$, and its
associate imaginary part $\varrho$. Here $\lambda$ and $\varrho$ are the real
and imaginary part of the dispersion relation, respectively, and they are both
functions of the continuous parameter $\mu_{k}$. The condition for having
Turing instability is thus the existence of a $\mu_{k}$ for which
$\lambda(\mu_{k})>0$. Let us observe that being the support discrete, there
can be finite size effects, namely $\lambda(\mu_{k})<0$ for all $k$ while the
continuous function $\lambda$ can assume positive values. This means that the
networked system may not exhibit Turing patterns, while the same exact system
defined on a continuous support will develop patterns.

Moreover, if $\lambda(\mu_{k})>0$ and the corresponding imaginary part is non-
zero, i.e., $\varrho(\mu_{k})\neq 0$, then the system can exhibit dynamical
(oscillatory) patterns [73]. Note that since we are considering case for which
$J_{k}$ is a real valued matrix, any non-zero value of $\varrho$ implies also
the presence of a solution with $-\varrho$. Let us finally remark that a non-
zero imaginary part of a critical mode is not a sufficient condition for the
existence of wave patterns, as there are cases in which the patterns are
stationary despite a non-zero $\varrho$ [74]. In fact, one cannot know a
priori which kind of pattern is obtained by solely using the information
contained in the linear stability analysis, and indeed the problem of pattern
prediction is still open.

##  5 Emergence of dynamical Turing patterns on nodes and links

In order to show evidence of the emergence of dynamical topological Turing
pattern, we consider the following dynamical system inspired by the excitable
dynamics of the FitzHugh-Nagumo neuronal model [75, 76]

$\begin{cases}\dot{u}=\sigma_{1}u-\eta_{1}u^{3}+\xi_{1}v+\zeta_{1}\mathbf{B}w-c_{2}\mathbf{L}_{[0]}(D_{uu}u+D_{uv}v)-c_{1}\alpha_{u}\mathbf{B}w\\\
\dot{v}=\sigma_{2}v+\xi_{2}u+\zeta_{2}\mathbf{B}w-c_{2}\mathbf{L}_{[0]}(D_{vu}u+D_{vv}v)-c_{1}\alpha_{v}\mathbf{B}w\\\
\dot{w}=\sigma_{3}w+\zeta_{3}\mathbf{B}^{\top}u+\zeta_{4}\mathbf{B}^{\top}v-c_{2}D_{ww}\mathbf{L}_{[1]}-c_{1}\mathbf{B}^{\top}(\beta_{u}u+\beta_{v}v)\,,\end{cases}$
(97)

whose homogeneous equilibrium state is given by

$(u^{*},v^{*},w^{*})=(0,0,0)\,.$

######  Remark 1.

Let us observe that in this particular case, the homogeneous vector
$(0,0,0)\otimes(\mathbf{1}_{N_{0}},\mathbf{1}_{N_{0}},\mathbf{1}_{N_{1}})$
trivially belongs to the kernel of the Dirac operator for any network.
However, for sake of continuity with the general framework above presented, we
decided to apply the model (97) to an Eulerian graph as support.

Let us observe that a similar claim could be applied any time the homogeneous
equilibrium is of the form $(u^{*},v^{*},0)$, because the network Laplace
matrix $\mathbf{L}_{[0]}$ admits $(\mathbf{1}_{N_{0}},\mathbf{1}_{N_{0}})$ as
eigenvector associated to the $0$ eigenvalue for any connected network.

In order for the homogeneous equilibrium on the nodes to exist, we need the
network to be connected. The same condition for the variables on the links is
attained if the network is Eulerian, i.e., every node has an even degree, as
shown in [40]. Let us observe that the equilibrium of the chosen model would
allow to relax the latter condition, nonetheless we want our theory to be
general. The support on which the dynamics take place is a $2$D lattice with
periodic boundary condition of $36$ nodes, i.e., $6\times 6$.

The conditions (92) ensuring the stability of the homogeneous equilibrium
become in the present case

$\begin{cases}\sigma_{2}<-\sigma_{1}\\\
\displaystyle\sigma_{2}>\frac{\xi_{2}\xi_{1}}{\sigma_{1}}\\\
\sigma_{3}<0\,.\end{cases}$ (98)

As shown in A, although obtaining the explicit form of the conditions for
Turing instability is challenging, relying on numerical results is possible.

Here we show results of dynamical Turing and Dirac pattern on two different
Eulerian network structures: a $2D$ square lattice tessellating a torus and
random graph with given degree sequence.

###  5.1 Phenomenology on a torus (2D square lattice with periodic boundary
conditions)

We consider a $2D$ square lattice tessellating the torus. Specifically we
assume to deal with a $6\times 6$ square lattice with periodic boundary
conditions. In Fig. 2, we report the bifurcation diagram of the model in the
parameter space $(\sigma_{1},\sigma_{2})$ computed from the dispersion
relation for the cases in which both $c_{1}$ and $c_{1}$ are non-zero (panel
a) and in which the interaction is exclusively driven by the Dirac operator,
i.e., $c_{2}=0$ (panel b).

In both bifurcation diagrams the red curves, i.e., $\sigma_{1}+\sigma_{2}=0$
and $\sigma_{1}\sigma_{2}=\xi_{1}\xi_{2}$, determine the boundary of the
stability region of the homogeneous equilibrium, which is unstable outside the
region (white region). The system exhibits patterns for parameters in the
green and blue regions. In the latter, patterns are dynamics, while in the
former they are stationary, because $\varrho$, i.e., the imaginary part
associated to unstable mode $\lambda$, is zero. The yellow region determines
parameters values for which the homogeneous equilibrium remains stable, even
after an inhomogeneous perturbation; hence, no patterns are observed. Note
that, for this choice of parameters, the latter region is observed in presence
of diffusion, i.e., $c_{2}>0$, while it disappears in the case of the
exclusive Dirac coupling, i.e. $c_{2}=0$ (see panel b)).

Figure 2: Instability regions in the $(\sigma_{1},\sigma_{2})$ space The red
curves are the bifurcation curves determining the stability of the homogeneous
equilibrium, namely, the first two equations of system (98). Panel a) shows
the case with diffusion (Turing patterns), while the case of Dirac-induced
patterns is shown in panel b). Having fixed $\sigma_{3}<0$, the homogeneous
equilibrium is unstable in the white region, while it is stable in the colored
area. The blue shaded area is where the patterns are oscillatory, while in the
green one they are stationary. The yellow region is where the homogeneous
equilibrium remains stable even after the spatially inhomogeneous
perturbation, which never occurs for Dirac-induced patterns with this set of
parameters. The colored regions are obtained by computing the real dispersion
relation and, if it is positive, checking if the imaginary part of the
positive modes is zero (green region) or nonzero (blue region). The parameters
are
$\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
for panel a) we have $\sigma_{3}=-2,~{}c_{2}=1$, while for panel b)
$\sigma_{3}=-10,~{}c_{2}=0$.

We now compare the dynamical patterns observed in presence of diffusion
($c_{2}\neq 0$) and in the case when it is absent ($c_{2}=0$). We chose
parameters values such that, in principle, dynamical (oscillatory) Turing
patterns are obtained, i.e., parameters in blue region. Let us first discuss
the phenomenology of the dynamical Turing patterns observed in presence of
diffusion. The dynamical nature of these topological patterns is revealed by
the dispersion relation displayed in Fig. 3a), where we can observe that the
mode associated to $\mu_{k}=0$ is stable (indeed $\lambda(0)<0$), while there
exist some modes that are unstable ($\lambda(\mu_{k})>0$ for some $k$) driving
the formation of Turing patterns. In panel b) of the same figure, we can
observe that the imaginary part of the dispersion relation $\varrho$ is non-
zero for corresponding unstable modes. For this reason, the obtained Turing
patterns on nodes and links are dynamical, i.e., oscillatory, as it is
apparent from Fig. 3c), where we display the time series for the two species
on the nodes, i.e., $u$ and $v$, and for the species on the links, i.e., $w$.
Let us point out that the form of Turing patterns is perturbation dependent
and different perturbation yield different (oscillatory) patterns. Lastly, in
Fig. 4, we display a snapshot of the patterns for species $u$ on the nodes and
$w$ on the links visualized on the network111Such dynamics can be better
appreciated in the video available as Supplementary Material of the online
version..

Figure 3: Dispersion relation and dynamical Turing patterns Panel a) shows the
real part of the dispersion relation $\lambda$, while panel b) its respective
imaginary part $\varrho$. The blue curves are computed by considering a
continuous parameter $\mu_{k}$, while the cyan dots indicate the eigenvalues
of the Hodge Laplacians. The latter determine the effective dispersion
relation and the emergence or not of Turing patterns. The chosen setting is
such that Turing patterns emerge ($\lambda>0$) and they are oscillatory (the
corresponding $\varrho\neq 0$), as shown by the dynamics of the $3$
interacting species, in panel c). The network is a $2$D lattice with periodic
boundary conditions of $36$ nodes and the parameters are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\sigma_{3}=-2,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$. Figure 4: Turing patterns
on the 2D lattice with periodic boundary conditions Snapshot, at $60$ time
units, of the patterns for species $u$ (nodes) and $w$ (links) visualized on
the network. Such dynamics can be better appreciated in the video (see the
Supplementary Material of the online version of this paper). The periodicity
of the lattice can be visualized thanks to the node labeling. The parameters
are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\sigma_{3}=-2,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$.

We now compare the phenomenology of the Turing patterns obtained in presence
of diffusion, i.e., $c_{2}>0$, with the ones arising once solely using the
Dirac operator, i.e., $c_{2}=0$. Also in this case, we consider a set of
parameters for which the stability of the homogeneous equilibrium is lost once
perturbed, as revealed by the real part of the dispersion relation in Fig.
5a). We can also remark that, for this choice of parameters, the latter is
divergent for large $\mu_{k}$. When working on continuous support, where all
instability modes are present, this is a problem because the maximum critical
mode has an infinite wave number, hence the solution is not a physical one
(long wave instability). However, on discrete support the finite size of the
latter induces an upper bound on the largest possible Laplace eigenvalue. In
particular it is well known that on networks the maximum eigenvalue of the
$0$-Laplacian (and the isomorphic $1$-Laplacian) is bounded for networks with
bounded degrees. Finally, even for networks in with unbounded degrees it is
possible to consider a bounded dispersion relation by adopting the normalized
Dirac [77] and Laplacian operators. A further discussion of the dispersion
relation for large $\mu_{k}$ can be found in Appendices A and B. By
considering the obtained dispersion relation we note that also in this case,
as in the previous one, its imaginary part $\varrho$ is non-zero for the
unstable modes (see Figure 5, panel b)); the resulting dynamical patterns can
be visualized in panel c). Lastly, in Fig. 6, we display a snapshot of the
patterns for species $u$ on the nodes and $w$ of the links visualized on the
$2$D lattice222The dynamics can be better appreciated in the video available
as Supplementary Material of the online version.. Remarkably, Dirac-induced
patterns present several distinct dynamical properties that are not present in
Turing pattern observed for $c_{2}\neq 0$.

Figure 5: Dispersion relation and Dirac-induced patterns Panel a) shows the
real part of the dispersion relation $\lambda$, while panel b) its respective
imaginary part $\varrho$. The dispersion relation is computed as for the case
with diffusion (Fig. 3, panels a and b). We can observe that the real part of
the dispersion relation is divergent, however this does not pose a problem
given that our support is discrete and finite, and so is its spectrum. The
chosen setting is such that Dirac-induced patterns emerge ($\lambda>0$) and
they are oscillatory (the corresponding $\varrho\neq 0$), as shown by the
dynamics of the $3$ interacting species, in panel c). The network is a $2$D
lattice with periodic boundary conditions of $36$ nodes and the parameters are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\sigma_{3}=-10,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=0,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$. Figure 6: Dirac-induced
patterns on the 2D lattice with periodic boundary conditions Snapshot, at $60$
time units, of the patterns for species $u$ (nodes) and $w$ (links) visualized
on the network. Such dynamics can be better appreciated in the video (see the
Supplementary Material of the online version of this paper). The parameters
are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\sigma_{3}=-10,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=0,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$.

In particular we observe that nodes and links are phase-synchronized and
divided in two clusters (see Fig. 5c)), which recalls a phenomenon called
cluster synchronization, often occurring when the network possesses symmetries
[78], as it is for our case. Moreover, such patterns do not change when
varying the initial conditions, at contrast with Turing patterns. To better
visualize the difference in the behavior of the two dynamical patterns,
Laplace-driven and Dirac-induced, respectively, let us project the dynamics in
the $3$D phase space $(u,v,w)$. The results are shown in Fig. 7, with panel a)
depicting the “attractor” for Laplace-driven patterns, while panel b) for
Dirac-induced patterns. More precisely, we select one link, $\ell$, we
consider one of its end-node, $i$, and we build the orbit
$(u_{i}(t),v_{i}(t),w_{\ell}(t))$. Let us observe that, while a representation
in the $(u,v)$ phase space is unique, the shown representation is not, because
there are more links thank nodes and there is an ambiguity in the selection of
the end-node for each link. For what concerns Dirac-induced patterns, however,
even by reshuffling the $w$ in the triplets, the attractor does not change.
This is different for Turing patterns and, indeed, the $3$D attractor as the
sole scope of showing the difference between the two cases, rather than an
accurate representation of the former. An important fact is that the Dirac
attractor is robust with respect to the initial condition and all the
trajectories converge to it, no matter the size of the initial perturbation.
Such behavior is remarkably different from the one arising with Turing
patterns, which are highly sensible to even small variations on the initial
conditions.

Figure 7: Projected attractors a) Turing patterns, b) Dirac-induced patterns.
The network is a $2$D lattice with periodic boundary conditions of $36$ nodes.
The parameters are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
for panel a) we have $\sigma_{3}=-2,~{}c_{2}=1$, while for panel b)
$\sigma_{3}=-10,~{}c_{2}=0$; the initial perturbation is the order of
$10^{-2}$.

###  5.2 Phenomenology on a random graph

We further investigate the nature and the phenomenology of the patterns by
analyzing the case of a random Eulerian network, shown in Fig. 8. It is
important to deal with a random structure, so we can better understand what
causes the regularity of the Dirac-induced pattern on the lattice.
Specifically we consider a random network with $N_{0}=20$ nodes and degree
distribution $P(2)=8/20$, $P(4)=9/20$, $P(6)=3/20$ and $P(k)=0$ for all $k\neq
2,4,6$. Specifically this network has all the nodes with even degree, hence it
is Eulerian. This implies that if we label the nodes subsequently along any of
its Eulerian paths and we orient the links according to the node label, then
it is immediate to show that the network admits a constant eigenvector ${\bf
1}$ in the kernel of the Dirac operator.

Figure 8: Random Eulerian Network The considered random Eulerian network of
$N_{0}=20$ nodes with degree distribution $P(2)=8/20$, $P(4)=9/20$,
$P(6)=3/20$ and $P(k)=0$ for all $k\neq 2,4,6$.

In order to investigate the dynamical properties of the model (97) on this
random topology, we distinguish between Turing patterns, i.e., once the
Laplace matrix is present, and Dirac-induced patterns, once the coupling is
solely realized by using the Dirac operator. In both cases we consider
parameter values for which the homogeneous equilibrium becomes unstable and
(in principle) oscillatory. In Fig. 10(a-b), we plot the dispersion relation
corresponding to this choice of parameter values and, in panel c), we display
the resulting dynamical Turing patterns. These are qualitatively similar to
those obtained on the lattice network, with the only difference that the
transient is longer. The latter is due to the fact that, for this choice of
parameters, the critical mode is closer to zero (compare Fig. 9a) with
3a))333By comparing the Figures, one can observe that the spectral gap of the
random Eulerian network is larger, hence motivating a slightly different
choice of parameters. Let us recall that, on networks, the conditions for
instability are necessary but not sufficient..

Figure 9: Dispersion relation and Turing patterns Panel a) shows the
dispersion relation $\lambda$, while panel b) its associated imaginary part
$\varrho$. As for the case of $2$D lattice, the chosen setting is such that
Dirac-induced patterns emerge ($\lambda>0$) and they are oscillatory (the
corresponding $\varrho\neq 0$), as shown by the dynamics of the $3$
interacting species, in panel c). The network is that of Fig. 8 and the
parameters are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\sigma_{3}=-2,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$.

Let us now repeat the numerical analysis for the case without diffusion, i.e.,
$c_{2}=0$. Again, the chosen setting yields oscillatory patterns, i.e.,
parameters in the blue region of Fig. 2b), as can be appreciated by looking at
the dispersion relation, depicted in Fig. 10(a-b); the resulting Dirac-induced
patterns can be visualized in panel c). We can observe that they are not
clustered as for the case on lattice, but we can somehow spot some kind of
regularity. We corroborate this qualitative observation by plotting the
trajectories in the phase space, shown in Fig. 11, for Turing patterns (panel
a) and Dirac-induced patterns (panel b). There, we can see that Turing
patterns are qualitatively equivalent to those obtained on the lattice
network; moreover, as expected, they are not robust with respect to the
initial perturbation. On the contrary, Dirac-induced patterns are, again,
robust to the initial conditions and their regularity can be visualized
through the closed curves in the phase space.

Figure 10: Dispersion relation and Turing patterns Panel a) shows the real
part of the dispersion relation $\lambda$, while panel b) its respective
imaginary part $\varrho$. As for the case of $2$D lattice, the chosen setting
is such that Dirac-induced patterns emerge ($\lambda>0$) and they are
oscillatory (the corresponding $\varrho\neq 0$), as shown by the dynamics of
the $3$ interacting species, in panel c). The network is that of Fig. 8 and
the parameters are
$\sigma_{1}=1,~{}\sigma_{2}=-1.4,~{}\sigma_{3}=-10,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=0,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$. Figure 11: Projected
attractors a) Turing patterns, b) Dirac-induced patterns. The network is the
Random Eulerian Network of Fig. 8. a) The parameters are
$\sigma_{1}=0.3,~{}\sigma_{2}=-2,~{}\sigma_{3}=-2,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$.
b) The parameters are
$\sigma_{1}=0.3,~{}\sigma_{2}=-2,~{}\sigma_{3}=-10,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}c_{1}=-1,~{}\alpha_{u}=1.1,~{}\alpha_{v}=7,~{}\beta_{u}=7,~{}\beta_{v}=1.1$;
the initial perturbation is the order of $10^{-2}$.

These observations allow us to conclude that Dirac-induced patterns are robust
with respect to the initial conditions and, might be sensible to the
symmetries of the underlying network topology.

##  6 Conclusions

The topological approach to network dynamics implies that the dynamical state
of a network is encoded in the direct sum of topological signals on nodes
($0$-cochains) and topological signals on links ($1$-cochains). This approach
significantly changes our description of dynamical processes on networks.
Topological signals on nodes and links can be coupled by the Dirac operator
that represents the “square-root” of the Laplacian. Here we have shown that
dynamical Turing and Dirac-induced patterns can set up on nodes and links of
Eulerian networks as a result of the instability of the homogeneous steady
state solution. This approach extends previous results on node-based Turing
patterns defined on networks, simplicial complexes and time-varying networks.
Let us observe that the setting where a single dynamical variable is
associated to nodes and a single dynamical variable is associated to links
cannot account for dynamical Turing and Dirac patterns [40]. Here we have
shown that taking three topological signals (two on the nodes and one on the
links) can instead account for dynamical patterns on nodes and links. In order
to describe their dynamics we have defined the $3$-way Hodge-Dirac operator
that is projecting the node signals into the links and, vice versa, is
projecting the link signals into the nodes. This operator is then combined to
an appropriate gamma matrix that has the ability to compress the projected
node signals into a single link signal and expand the projected link signal
into two distinct node signals. In the absence of coupling, the topological
signals of each dimension obey the same dynamical equation, while the coupling
can be induced by either the higher-order Laplacian matrix or the $3$-way
Dirac operator, or a combination of both operators. We called Dirac-induced
patterns the ones induced exclusively by the Dirac operator, while we call all
the other patterns Turing patterns. The formulated theory of Turing and Dirac
induced patterns leads to the instability of the uniform steady state solution
and the onset of a dynamical behavior that displays remarkable differences
observed for a significant region of the parameter space. The difference
between the observed dynamical behavior of Turing and Dirac induced patterns
is here highlighted by considering two very distinct topologies: a
$2$-dimensional torus tessellated by a square lattice and a random (Eulerian)
graph. In both cases, it emerges that the Dirac induced patterns have distinct
dynamical properties displaying a clustering of the dynamics characterized by
an increased robustness to noise with respect to Turing patterns.

We believe that this work can open new scenarios in the treatment of
topological Turing instabilities affect jointly nodes and links signals with
possible applications to neuroscience.

## Author contribution statement

Riccardo Muolo: Conceptualization, Formal analysis, Methodology, Software,
Investigation, Writing - Original draft preparation
Timoteo Carletti: Formal analysis, Investigation, Writing - Original draft
preparation
Ginestra Bianconi: Conceptualization, Formal analysis, Methodology,
Investigation, Writing - Original draft preparation

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Data availability

Data will be made available upon reasonable request to the corresponding
author.

## Funding

R.M. acknowledges JSPS KAKENHI JP22K11919, JP22H00516, and JST CREST JP-
MJCR1913 for financial support. During part of this work, R.M. was supported
by a FRIA-FNRS Fellowship, funded by the Walloon Region, Grant FC 33443.

## Acknowledgements

The authors are grateful to Benedetta Franceschiello, Lorenzo Giambagli and
Lucille Calmon for interesting discussions.

##  Appendix A Conditions for the emergence of Turing patterns

We can determine the conditions to obtain Turing patterns by studying the
stability of the following polynomial

$a\lambda_{k}^{3}+b\lambda_{k}^{2}+c\lambda_{k}+d=0,$ (A1)

whose coefficients are

$\begin{cases}a=1,\\\
b=c_{2}\mu_{k}^{2}(D_{uu}+D_{vv}+D_{ww})-(f_{u}+g_{v}+h_{w}),\\\
c=\pi_{1}(\mu_{k}^{2}),\\\ d=\pi_{2}(\mu_{k}^{2}),\end{cases}$

where $\pi_{1}$ and $\pi_{2}$ are polynomials in $\mu_{k}^{2}$.

The first polynomial has degree $2$ and has the following coefficients
$\pi_{1}(\mu_{k}^{2})=\mathcal{A}\mu_{k}^{4}+\mathcal{B}\mu_{k}^{2}+\mathcal{C}$,
explicitly

$\begin{cases}\mathcal{A}=(D_{vv}D_{ww}+D_{uu}D_{ww}+D_{uu}D_{vv}-D_{uv}D_{vu})c_{2}^{2},\\\
\mathcal{B}=[D_{uv}g_{u}+D_{vu}f_{v}-D_{uu}(h_{w}+g_{v})+\\\
~{}~{}~{}~{}~{}~{}~{}~{}-D_{vv}(h_{w}+f_{u})-D_{ww}(f_{u}+g_{v})]c_{2}-(f_{w}h_{u}+g_{w}h_{v}),\\\
\mathcal{C}=f_{u}g_{v}-f_{v}g_{u}+h_{w}(f_{u}+g_{v}),\end{cases}$

while the second one has degree $3$ and coefficients are
$\pi_{2}(\mu_{k}^{2})=\mathcal{E}\mu_{k}^{6}+\mathcal{F}\mu_{k}^{4}+\mathcal{G}\mu_{k}^{2}+\mathcal{H}$,
explicitly

$\begin{cases}\mathcal{E}=D_{ww}(D_{uu}D_{vv}-D_{uv}D_{vu}){c_{2}}^{3},\\\
\mathcal{F}=[f_{w}(h_{v}D_{vu}-h_{u}D_{vv})+g_{w}(h_{u}D_{uv}-h_{v}D_{uu})]c_{2}+\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+[D_{ww}(f_{v}D_{vu}-f_{u}D_{vv}+g_{u}D_{uv}-g_{v}D_{uu})+h_{w}(D_{uu}D_{vv}+D_{uv}D_{vu})]c_{2}^{2},\\\
\mathcal{G}=g_{w}(f_{u}h_{v}-f_{v}h_{u})+f_{v}(g_{v}h_{u}-g_{u}h_{v})+[D_{ww}(f_{u}g_{v}-f_{v}g_{u})+h_{w}(f_{u}D_{vv}+g_{v}D_{uu}-f_{v}D_{vu}-g_{u}D_{uv}]c_{2}.\\\
\mathcal{H}=h_{w}(f_{v}g_{u}-f_{u}g_{v})\end{cases}$

Let us observe that

$D_{uu}D_{vv}-D_{uv}D_{vu}=\alpha_{u}\beta_{u}\alpha_{v}\beta_{v}-\alpha_{v}\beta_{v}\alpha_{v}\beta_{u}=0,$

hence $\mathcal{E}=0$ and $\pi_{2}(\mu_{k}^{2})$ is de facto a polynomial of
degree $2$, i.e.,
$\pi_{2}(\mu_{k}^{2})=\mathcal{F}\mu_{k}^{4}+\mathcal{G}\mu_{k}^{2}+\mathcal{H}$.

The Routh-Hurwitz criterion [71, 72] gives us the following necessary and
sufficient conditions for stability

$\displaystyle a>$ $\displaystyle~{}0,$ (A2) $\displaystyle b>$
$\displaystyle~{}0,$ $\displaystyle bc-ad>$ $\displaystyle~{}0,$
$\displaystyle d>$ $\displaystyle~{}0.$

We are interested in the conditions to obtain Turing patterns, i.e., at least
one of the roots of the polynomial has a positive real part. This means that
it is enough that one of the (A2) is violated.

The first and second conditions cannot be violated: $a>0$ is trivial, while
$b>0$ needs to be satisfied due to the stability of the homogeneous
equilibrium, which gives us that $f_{u}+g_{v}<0$ and $h_{w}<0$. Let us hence
study the fourth condition $d>0$. Namely, we need check if the polynomial
$\pi_{2}(\mu_{k}^{2})=\mathcal{F}\mu_{k}^{4}+\mathcal{G}\mu_{k}^{2}+\mathcal{H}$
can take negative values. We can observe that
$\mathcal{H}=h_{w}(f_{v}g_{u}-f_{u}g_{v})>0$, again for the conditions on the
homogeneous equilibrium. Then we need that either $\mathcal{F}<0$ (so that we
are sure that $\pi_{2}$ takes some negative values, independently of
$\mathcal{G}$) or $\mathcal{G}<0$ ($\pi_{2}$ may be negative, independently of
$\mathcal{F}$, but not for every set of parameters). The same principles apply
to imposing the negativity of the third of the Routh-Hurwitz conditions (A2),
i.e., $bc-ad<0$, which can be obtained in an analogous way. In the main text
we put ourselves in a setting in which Eqs. (A2) are violated and hence the
system exhibits Turing patterns, which are, as expected, dynamics.

Let us conclude this section by studying the asymptotic behavior of the roots
of the polynomial (A1) for $\mu_{k}\gg 1$. The first step is to rewrite the
latter as follows

$p(\lambda)=\mu_{k}^{4}A(\lambda)+\mu_{k}^{2}B(\lambda)+C(\lambda)\,,$

where $A$, $B$ and $C$ are polynomials in $\lambda$ of degree $1$, $2$ and
$3$. By dividing the previous expression by $\mu_{k}^{4}$ we obtain a second
polynomial in $\lambda$ with the same roots as the previous one. Let us first
consider $\mu_{k}\gg 1$ and write the root of the polynomial as
$\lambda(\mu_{k})=\lambda_{0}+\lambda_{1}/\mu_{k}^{2}+\mathcal{O}(\mu_{k}^{-4})$.
Then, it follows that $\lambda_{0}$ satisfies $A(\lambda_{0})=0$ and the
second term $\lambda_{1}$ can be obtained by solving
$A^{\prime}(\lambda_{0})\lambda_{1}+B(\lambda_{0})=0$. Because $A$ is a linear
function, its root can be straightforwardly determined

$\displaystyle\lambda_{0}=-\frac{1}{c_{2}(\alpha_{u}\beta_{u}+\alpha_{v}\beta_{v})^{2}}$
$\displaystyle\Big{[}$
$\displaystyle\alpha_{v}\beta_{v}f_{w}h_{u}-\alpha_{v}\beta_{u}f_{w}h_{v}+\alpha_{v}^{2}\beta_{v}^{2}c_{2}f_{u}-\alpha_{v}^{2}\beta_{u}\beta_{v}c_{2}f_{v}-\alpha_{u}\beta_{v}g_{w}h_{u}$
$\displaystyle+$
$\displaystyle\alpha_{u}\beta_{u}g_{w}h_{v}-\alpha_{u}\alpha_{v}\beta_{v}^{2}c_{2}g_{u}+\alpha_{u}\alpha_{v}\beta_{u}\beta_{v}c_{2}g_{v}+\alpha_{u}\alpha_{v}\beta_{u}\beta_{v}c_{2}f_{u}$
$\displaystyle-$
$\displaystyle\alpha_{u}\alpha_{v}\beta_{u}^{2}c_{2}f_{v}-\alpha_{u}^{2}\beta_{u}\beta_{v}c_{2}g_{u}+\alpha_{u}^{2}\beta_{u}^{2}c_{2}g_{v}\Big{]}.$

This implies that for sufficiently large $\mu_{k}$ the roots of (A1) will be
close to $\lambda_{0}$ within a term of order $1/\mu_{k}^{2}$.

##  Appendix B Conditions for the emergence of Dirac-induced patterns

To compute the conditions yielding the emergence of Dirac-induced patterns,
one can proceed in the same way as before. Moreover, let us remark that the
conditions for the stability of the homogeneous equilibrium are the same of
the previous case, of the uncoupled dynamics has not changed.

The fact that $c_{2}=0$ simplifies the coefficients of the polynomial whose
stability needs to be determined. In fact, now we are dealing with the
following:

$a\lambda_{k}^{3}+b\lambda_{k}^{2}+c\lambda_{k}+d=0,$ (B1)

where

$\begin{cases}a=1\\\ b=-(f_{u}+g_{v}+h_{w}),\\\ c=\varrho_{1}(\mu_{k}^{2}),\\\
d=\varrho_{2}(\mu_{k}^{2}),\end{cases}$

and $\varrho_{1}$ and $\varrho_{2}$ are now linear curves in $\mu_{k}^{2}$.

The first one has the following form
$\varrho_{1}(\mu_{k}^{2})=\mathcal{X}\mu_{k}^{2}+\mathcal{Y}$, explicitly

$\begin{cases}\mathcal{X}=-(f_{w}h_{u}+g_{w}h_{v}),\\\
\mathcal{Y}=f_{u}g_{v}-f_{v}g_{u}+h_{w}(f_{u}+g_{v}),\end{cases}$

while the second one is
$\varrho_{2}(\mu_{k}^{2})=\mathcal{Z}\mu_{k}^{2}+\mathcal{K}$, explicitly:

$\begin{cases}\mathcal{Z}=g_{w}(f_{u}h_{v}-f_{v}h_{u})+f_{v}(g_{v}h_{u}-g_{u}h_{v}),\\\
\mathcal{K}=h_{w}(f_{v}g_{u}-f_{u}g_{v}).\end{cases}$

The Routh-Hurwitz criterion [71, 72] gives us the following necessary and
sufficient conditions for stability

$\displaystyle a>$ $\displaystyle~{}0,$ (B2) $\displaystyle b>$
$\displaystyle~{}0,$ $\displaystyle bc-ad>$ $\displaystyle~{}0,$
$\displaystyle d>$ $\displaystyle~{}0.$

As before, $a$ and $b$ are always positive and such conditions can never be
violated. Let us hence study the condition $d<0$, i.e., we need to impose that
$\varrho_{1}<0$. Since the coefficient $\mathcal{K}$ is always positive, the
condition reduces to imposing $\mathcal{Z}$ negative

$g_{w}(f_{u}h_{v}-f_{v}h_{u})+f_{v}(g_{v}h_{u}-g_{u}h_{v})<0,$ (B3)

which is our first condition for Dirac-induced patterns. Let us then proceed
in studying the Routh-Hurwitz condition $bc-ad<0$, which translates into a new
linear curve $\varrho_{3}(\mu_{k}^{2})=\mathcal{J}\mu_{k}^{2}+\mathcal{L}$,
where now

$\begin{cases}\mathcal{J}=b\mathcal{X}-\mathcal{Z}=(f_{u}+g_{v}+h_{w})(f_{w}h_{u}+g_{w}h_{v})-g_{w}(f_{u}h_{v}-f_{v}h_{u})-f_{v}(g_{v}h_{u}-g_{u}h_{v}),\\\
\mathcal{L}=b\mathcal{Y}-\mathcal{K}=-(f_{u}+g_{v}+h_{w})[f_{u}g_{v}-f_{v}g_{u}+h_{w}(f_{u}+g_{v})]-h_{w}(f_{v}g_{u}-f_{u}g_{v}).\end{cases}$

Hence, patterns are obtained if $\mathcal{J}<0$ or if $\mathcal{L}<0$, namely
if at least one the two condition

$\begin{cases}(f_{u}+g_{v}+h_{w})(f_{w}h_{u}+g_{w}h_{v})<g_{w}(f_{u}h_{v}-f_{v}h_{u})+f_{v}(g_{v}h_{u}-g_{u}h_{v}),\\\
(f_{u}+g_{v}+h_{w})[f_{u}g_{v}-f_{v}g_{u}+h_{w}(f_{u}+g_{v})]>h_{w}(f_{v}g_{u}-f_{u}g_{v}),\end{cases}$
(B4)

is satisfied. Let us observe that the second condition is never satisfied
because of the conditions for the homogeneous equilibrium. In fact, the
expression $(f_{u}+g_{v}+h_{w})[f_{u}g_{v}-f_{v}g_{u}+h_{w}(f_{u}+g_{v})]$ is
always negative, while $h_{w}(f_{v}g_{u}-f_{u}g_{v})$ is always positive.

Summing up, the conditions for Dirac-induced patterns are the following

$\begin{cases}g_{w}(f_{u}h_{v}-f_{v}h_{u})+f_{v}(g_{v}h_{u}-g_{u}h_{v})<0,\\\
(f_{u}+g_{v}+h_{w})(f_{w}h_{u}+g_{w}h_{v})<g_{w}(f_{u}h_{v}-f_{v}h_{u})+f_{v}(g_{v}h_{u}-g_{u}h_{v}).\end{cases}$
(B5)

Also in this case, we can study the asymptotic behavior of the roots of the
polynomial (B1) for $\mu_{k}\gg 1$. Let us thus rewrite the latter polynomial
as follows

$p(\lambda)=\mu_{k}^{2}\tilde{B(}\lambda)+\tilde{C}(\lambda)\,,$

where $\tilde{B}$ and $\tilde{C}$ are polynomials in $\lambda$ of degree $1$
and $3$. By dividing the previous expression by $\mu_{k}^{2}$, we obtain again
a polynomial in $\lambda$ sharing the same roots as the previous one. Let us
assume $\mu_{k}\gg 1$ and write the root of the polynomial as
$\lambda(\mu_{k})=\lambda_{0}+\lambda_{1}/\mu_{k}^{2}+\mathcal{O}(\mu_{k}^{-4})$.
Then, $\lambda_{0}$ satisfies $\tilde{B}(\lambda_{0})=0$ and the second term
$\lambda_{1}$ can be obtained by solving
$\tilde{B}^{\prime}(\lambda_{0})\lambda_{1}+\tilde{C}(\lambda_{0})=0$. Being
$\tilde{B}$ a first degree polynomial, its root can be straightforwardly
determined

$\lambda_{0}=-\frac{f_{w}g_{v}h_{u}-f_{w}g_{u}h_{v}-f_{v}g_{w}h_{u}+f_{u}g_{w}h_{v}}{g_{w}h_{v}+f_{w}h_{u}}\,.$

This implies that for sufficiently large $\mu_{k}$ the roots of (B1) will be
close to $\lambda_{0}$ within a correction of order $1/\mu_{k}^{2}$.

##  Appendix C Stationary patterns

Throughout this work, we have always dealt with dynamical patterns, which are
made possible by the formalism hereby introduced. However, as we have
mentioned in the Main Text, there are cases in which the patterns are
stationary, even when one deals with allowing for oscillatory instabilities,
such as $3$ interacting species [73] or hyperbolic reaction-diffusion systems
[79].
For sake of completeness, in Fig. C1, we show, for both the $2$D lattice and
the Random Eulerian Network, the dispersion relation and the patterns in a
setting in which the latter are stationary, i.e., the parameters are chosen
from the green region of Fig. 2a). We can see that, when the real part of the
dispersion relation is positive and the corresponding imaginary part is zero,
while the modes with non-zero imaginary part have a negative real part. Hence,
the pattern is stationary.
Let us conclude by stressing that the fact the imaginary part of the
dispersion relation is non-zero in correspondence of the unstable modes does
not automatically guarantee a dynamical pattern. In fact, for hyperbolic
reaction-diffusion systems on networks, there have been found settings in
which the modes responsible for the instability have zero imaginary part,
which is non-zero for negative modes, but patterns are oscillatory [74]. Also
the opposite case can verify, hence, one may have the unstable modes with zero
imaginary part, but the patterns can be oscillatory. The latter cases are not
common and in general the imaginary part gives a good indication of the
resulting pattern, even though one always needs to perform the numerical
experiment before claiming the nature of the pattern.

Figure C1: Dispersion relation and stationary Turing patterns In panels a) and
b) we show the real and imaginary part of the dispersion relation,
respectively, for the case of Turing patterns (i.e., with diffusion) on the
$2$D lattice. We can see that, when the real part of the dispersion relation
is positive and the corresponding imaginary part is zero, while the modes with
non-zero imaginary part have a negative real part. The dynamics of the
stationary patterns for the $3$ interacting species $(u,v,w)$ are shown in
panel c). Panels d), e) and f) show the analogous setting but for the Random
Eulerian Networks. The parameters for both cases are
$\sigma_{1}=0.2,~{}\sigma_{2}=-7.8,~{}\sigma_{3}=-2,~{}\eta_{1}=-3,~{}\xi_{1}=-0.5,~{}\xi_{2}=5,~{}\zeta_{1}=0.1,~{}\zeta_{2}=6,~{}\zeta_{3}=6,~{}\zeta_{4}=1,~{}c_{1}=-1,~{}c_{2}=1,~{}\alpha_{u}=1,~{}\alpha_{v}=1,~{}\beta_{u}=1,~{}\beta_{v}=0.1$;
the initial perturbation is the order of $10^{-2}$.

## References

  * Barabasi [2016] Albert-Laszlo Barabasi,  _Network Science_ (Cambridge University Press, 2016).
  * Newman [2010] Mark Newman,  _Networks_ (Oxford University Press, 2010).
  * Barrat  _et al._ [2008] Alain Barrat, Marc Barthelemy, and Alessandro Vespignani,  _Dynamical processes on complex networks_ (Cambridge University Press, 2008).
  * Dorogovtsev  _et al._ [2008] Sergey N Dorogovtsev, Alexander V Goltsev, and José FF Mendes, “Critical phenomena in complex networks,” Reviews of Modern Physics 80, 1275 (2008).
  * Santoro  _et al._ [2022] A. Santoro, F. Battiston, G. Petri, and E. Amico, “Unveiling the higher-order organization of multivariate time series,” arXiv preprint arXiv:2203.10702 (2022).
  * Faskowitz  _et al._ [2022] J. Faskowitz, R.F. Betzel, and O. Sporns, “Edges in brain networks: Contributions to models of structure and function,” Netw. Neurosci. 6, 1–28 (2022).
  * Giusti  _et al._ [2015] Chad Giusti, Eva Pastalkova, Carina Curto, and Vladimir Itskov, “Clique topology reveals intrinsic geometric structure in neural correlations,” Proceedings of the National Academy of Sciences 112, 13455–13460 (2015).
  * Reimann  _et al._ [2017] Michael W Reimann, Max Nolte, Martina Scolamiero, Katharine Turner, Rodrigo Perin, Giuseppe Chindemi, Paweł Dłotko, Ran Levi, Kathryn Hess, and Henry Markram, “Cliques of neurons bound into cavities provide a missing link between structure and function,” Frontiers in computational neuroscience 11, 48 (2017).
  * Petri  _et al._ [2013] Giovanni Petri, Martina Scolamiero, Irene Donato, and Francesco Vaccarino, “Topological strata of weighted complex networks,” PLoS One 8, e66506 (2013).
  * Lee  _et al._ [2021] Y. Lee, J. Lee, S.M. Oh, D. Lee, and B. Kahng, “Homological percolation transitions in growing simplicial complexes,” Chaos 31, 041102 (2021).
  * Bobrowski and Skraba [2022] Omer Bobrowski and Primoz Skraba, “Homological percolation: The formation of giant k-cycles,” International Mathematics Research Notices 2022, 6186–6213 (2022).
  * Santos  _et al._ [2019] Fernando AN Santos, Ernesto P Raposo, Maurício D Coutinho-Filho, Mauro Copelli, Cornelis J Stam, and Linda Douw, “Topological phase transitions in functional brain networks,” Physical Review E 100, 032414 (2019).
  * Mulas  _et al._ [2020] R. Mulas, C. Kuehn, and J. Jost, “Coupled dynamics on hypergraphs: Master stability of steady states and synchronization,” Phys. Rev. E 101, 062313 (2020).
  * Skardal and Arenas [2020] Per Sebastian Skardal and Alex Arenas, “Higher order interactions in complex networks of phase oscillators promote abrupt synchronization switching,” Communications Physics 3, 218 (2020).
  * Gambuzza  _et al._ [2018] Lucia Valentina Gambuzza, Mattia Frasca, and Vito Latora, “Distributed control of synchronization of a group of network nodes,” IEEE Transactions on Automatic Control 64, 365–372 (2018).
  * Gambuzza  _et al._ [2021] L.V Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, “Stability of synchronization in simplicial complexes,” Nat. Comm. 12, 1–13 (2021).
  * de Arruda  _et al._ [2021] G.F. de Arruda, M. Tizzani, and Y. Moreno, “Phase transitions and stability of dynamical processes on hypergraphs,” Comm. Phys. 4, 1–9 (2021).
  * Bianconi [2021a] G. Bianconi,  _Higher-Order Networks:An introduction to simplicial complexes_ (Cambridge University Press, 2021).
  * F. and et al. [2021] Battiston F. and et al., “The physics of higher-order interactions in complex systems,” Nat. Phys. 17, 1093–1098 (2021).
  * Millán  _et al._ [2020] A.P. Millán, J.J. Torres, and G. Bianconi, “Explosive higher-order Kuramoto dynamics on simplicial complexes,” Phys. Rev. Lett. 124, 218301 (2020).
  * Ghorbanchian  _et al._ [2021] Reza Ghorbanchian, Juan G Restrepo, Joaquín J Torres, and Ginestra Bianconi, “Higher-order simplicial synchronization of coupled topological signals,” Communications Physics 4, 120 (2021).
  * DeVille [2021] L. DeVille, “Consensus on simplicial complexes: Results on stability and synchronization,” Chaos 31, 023137 (2021).
  * Arnaudon  _et al._ [2022] A. Arnaudon, R.L. Peach, G. Petri, and P. Expert, “Connecting Hodge and Sakaguchi-Kuramoto through a mathematical framework for coupled oscillators on simplicial complexes,” Comm. Phys. 5, 1–12 (2022).
  * Carletti  _et al._ [2023] T. Carletti, L. Giambagli, and G. Bianconi, “Global topological synchronization on simplicial and cell complexes,” Phys. Rev. Lett. 130, 187401 (2023).
  * Torres and Bianconi [2020] J.J. Torres and G. Bianconi, “Simplicial complexes: higher-order spectral dimension and dynamics,” J. phys. Complex. 1, 015002 (2020).
  * Ziegler  _et al._ [2022] C. Ziegler, P.S. Skardal, H. Dutta, and D. Taylor, “Balanced Hodge Laplacians optimize consensus dynamics over simplicial complexes,” Chaos 32, 023128 (2022).
  * Schaub  _et al._ [2020] M.T. Schaub, A.R. Benson, P. Horn, G. Lippner, and A. Jadbabaie, “Random walks on simplicial complexes and the normalized Hodge 1-Laplacian,” SIAM Rev. 62, 353–391 (2020).
  * Battiloro  _et al._ [2023] Claudio Battiloro, Lucia Testa, Lorenzo Giusti, Stefania Sardellitti, Paolo Di Lorenzo, and Sergio Barbarossa, “Generalized simplicial attention neural networks,” arXiv preprint arXiv:2309.02138 (2023).
  * Barbarossa and Sardellitti [2020] Sergio Barbarossa and Stefania Sardellitti, “Topological signal processing over simplicial complexes,” IEEE Transactions on Signal Processing 68, 2992–3007 (2020).
  * Schaub  _et al._ [2021] M.T. Schaub, Y. Zhu, J.-B. Seby, T.M. Roddenberry, and S. Segarra, “Signal processing on higher-order networks: Livin’on the edge… and beyond,” Signal Process. 187, 108149 (2021).
  * Ebli  _et al._ [2020] S. Ebli, M. Defferrard, and G. Spreemann, “Simplicial neural networks,” arXiv preprint https://arxiv.org/abs/2010.03633 (2020).
  * Hensel  _et al._ [2021] Felix Hensel, Michael Moor, and Bastian Rieck, “A survey of topological machine learning methods,” Frontiers in Artificial Intelligence 4, 681108 (2021).
  * Horak and Jost [2013] Danijela Horak and Jürgen Jost, “Spectra of combinatorial laplace operators on simplicial complexes,” Advances in Mathematics 244, 303–336 (2013).
  * Bianconi [2021b] G. Bianconi, “The topological dirac equation of networks and simplicial complexes,” J. phys. Complex. 2, 035022 (2021b).
  * Post [2009] O. Post, “First order approach and index theorems for discrete and metric graphs,” in  _Annales Henri Poincaré_ , Vol. 10 (Springer, 2009) pp. 823–866.
  * Bianconi [2022] Ginestra Bianconi, “Dirac gauge theory for topological spinors in 3+ 1 dimensional networks,” Journal of Physics A: Mathematical and Theoretical 56, 275001 (2022).
  * Bianconi [2023] Ginestra Bianconi, “The mass of simple and higher-order networks,” arXiv preprint arXiv:2309.07851 (2023).
  * Calmon  _et al._ [2022] Lucille Calmon, Juan G Restrepo, Joaquín J Torres, and Ginestra Bianconi, “Dirac synchronization is rhythmic and explosive,” Communications Physics 5, 253 (2022).
  * Calmon  _et al._ [2023a] Lucille Calmon, Sanjukta Krishnagopal, and Ginestra Bianconi, “Local Dirac synchronization on networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science 33 (2023a).
  * Giambagli  _et al._ [2022] L. Giambagli, L. Calmon, R. Muolo, T. Carletti, and G. Bianconi, “Diffusion-driven instability of topological signals coupled by the dirac operator,” Phys. Rev. E 106 (2022).
  * Calmon  _et al._ [2023b] Lucille Calmon, Michael T Schaub, and Ginestra Bianconi, “Dirac signal processing of higher-order topological signals,” New Journal of Physics 25 (2023b).
  * Nurisso  _et al._ [2023] Marco Nurisso, Alexis Arnaudon, Maxime Lucas, Robert L Peach, Paul Expert, Francesco Vaccarino, and Giovanni Petri, “A unified framework for simplicial kuramoto models,” arXiv preprint arXiv:2305.17977 (2023).
  * Lloyd  _et al._ [2016] S. Lloyd, S. Garnerone, and P. Zanardi, “Quantum algorithms for topological and geometric analysis of data,” Nat. Comm. 7, 1–7 (2016).
  * Wee  _et al._ [2023] JunJie Wee, Ginestra Bianconi, and Kelin Xia, “Persistent dirac for molecular representation,” Scientific Reports 13, 11183 (2023).
  * Murray [2001] J.D. Murray,  _Mathematical biology II: Spatial models and biomedical applications_ (Springer-Verlag, 2001).
  * Turing [1952] A.M. Turing, “The chemical basis of morphogenesis,” Phil. Trans. R. Soc. Lond. B 237, 37 (1952).
  * Nakamasu  _et al._ [2009] A. Nakamasu, G. Takahashi, A. Kanbe, and S. Kondo, “Interactions between zebrafish pigment cells responsible for the generation of turing patterns,” Proc. Natl. Acad. Sci. U.S.A. 106, 8429 (2009).
  * Byrne  _et al._ [2019] Áine Byrne, Daniele Avitabile, and Stephen Coombes, “Next-generation neural field model: The evolution of synchrony within patterns and waves,” Physical Review E 99, 012313 (2019).
  * Kato and Nakao [2022] Y. Kato and H. Nakao, “Turing instability in quantum activator–inhibitor systems,” Sci. Rep. 12 (1) (2022).
  * Fuseya  _et al._ [2021] Y. Fuseya, H. Katsuno, K. Behnia, and A. Kapitulnik, “Nanoscale turing patterns in a bismuth monolayer,” Nat. Phys. 17 (2021).
  * Othmer and Scriven [1971] H.G. Othmer and L.E. Scriven, “Instability and dynamic pattern in cellular networks,” J. Theor. Biol. 32, 507 (1971).
  * Othmer and Scriven [1974] H.G. Othmer and L.E. Scriven, “Non-linear aspects of dynamic pattern in cellular networks,” J. Theor. Biol. 43, 83–112 (1974).
  * Nakao and Mikhailov [2010] H. Nakao and A.S. Mikhailov, “Turing patterns in network-organized activator-inhibitor systems,” Nat. Phys. 6, 544 (2010).
  * Asllani  _et al._ [2014a] M. Asllani, J.D. Challenger, F.S. Pavone, L. Sacconi, and D. Fanelli, “The theory of pattern formation on directed networks,” Nature Communication 5 (2014a).
  * Muolo  _et al._ [2019] R. Muolo, M. Asllani, D. Fanelli, P.K. Maini, and T. Carletti, “Patterns of non-normality in networked systems,” J. Theor. Biol. 480, 81 (2019).
  * van der Kolk  _et al._ [2023] J. van der Kolk, G. García-Pérez, N.E. Kouvaris, M.A. Serrano, and M. Boguñá, “Emergence of geometric turing patterns in complex networks,” Phys. Rev. X 13 (2023).
  * Hata  _et al._ [2014a] S. Hata, H. Nakao, and A.S. Mikhailov, “Dispersal-induced destabilization of metapopulations and oscillatory turing patterns in ecological networks,” Sci. Rep. 4 (2014a).
  * Hata  _et al._ [2012] S. Hata, H. Nakao, and A.S. Mikhailov, “Global feedback control of turing patterns in network-organized activator-inhibitor systems,” Eur. Phys. Lett. 98, 64004 (2012).
  * Gao  _et al._ [2022] S. Gao, L. Chang, I. Romić, Z. Wang, M. Jusup, and P. Holme, “Optimal control of networked reaction–diffusion systems,” J. R. Soc. Interface 19, 20210739 (2022).
  * Buscarino  _et al._ [2019] A. Buscarino, C. Corradino, L. Fortuna, and M. Frasca, “Turing patterns via pinning control in the simplest memristive cellular nonlinear networks,” Chaos 29, 103145 (2019).
  * Asllani  _et al._ [2014b] M. Asllani, D.M. Busiello, T. Carletti, D. Fanelli, and G. Planchon, “Turing patterns in multiplex networks,” Phys. Rev. E 90, 042814 (2014b).
  * Asllani  _et al._ [2016] M. Asllani, T. Carletti, and D. Fanelli, “Tune the topology to create or destroy patterns,” Eur. Phys. J. B , 89 (2016).
  * Kouvaris  _et al._ [2015] N.E. Kouvaris, S. Hata, and A. Diaz-Guilera, “Pattern formation in multiplex networks,” Scientific Reports 5, 10840 (2015).
  * Petit  _et al._ [2017] J. Petit, B. Lauwens, D. Fanelli, and T. Carletti, “Theory of turing patterns on time varying networks,” Phys. Rev. Lett. 119, 148301 (2017).
  * Van Gorder [2021] R.A. Van Gorder, “A theory of pattern formation for reaction-diffusion systems on temporal networks,” Proc. R. Soc. A 477, 20200753 (2021).
  * Carletti and Fanelli [2022] T. Carletti and D. Fanelli, “Theory of synchronisation and pattern formation on time varying networks,” Chaos Solit. Fractals 156, 112180 (2022).
  * Carletti  _et al._ [2020] T. Carletti, D. Fanelli, and S. Nicoletti, “Dynamical systems on hypergraphs,” J. phys. Complex. 1, 035006 (2020).
  * Muolo  _et al._ [2023] R. Muolo, L. Gallo, V. Latora, M. Frasca, and T. Carletti, “Turing patterns in systems with high-order interaction,” Chaos Solit. Fractals 166, 112912 (2023).
  * Gao  _et al._ [2023] S. Gao, L. Chang, M. Perc, and Z. Wang, “Turing patterns in simplicial complexes,” Phys. Rev. E 107, 014216 (2023).
  * Ghorbanchian  _et al._ [2022] Reza Ghorbanchian, Vito Latora, and Ginestra Bianconi, “Hyper-diffusion on multiplex networks,” Journal of Physics: Complexity 3 (2022).
  * Routh [1877] E.J. Routh,  _Stability of a given state of motion_ (MacMillan and Co, London, 1877).
  * Hurwitz [1895] A. Hurwitz, “Ueber die bedingungen, unter welchen eine gleichung nur wurzeln mit negativen reellen theilen besitzt,” Math. Ann. 46, 273 (1895).
  * Hata  _et al._ [2014b] S. Hata, H. Nakao, and A.S. Mikhailov, “Sufficient conditions for wave instability in three-component reaction–diffusion systems,” Prog. Theor. Exp. Phys. 1 (2014b).
  * Carletti and Muolo [2021] T. Carletti and R. Muolo, “Finite propagation enhances Turing patterns in reaction–diffusion networked systems,” J. phys. Complex. 2, 045004 (2021).
  * FitzHugh [1961] R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophys. J. 1, 445 (1961).
  * Nagumo  _et al._ [1962] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” Proc. IRE 50, 2061 (1962).
  * Baccini  _et al._ [2022] F. Baccini, F. Geraci, and G. Bianconi, “Weighted simplicial complexes and their representation power of higher-order network data and topology,” Phsy. Rev. E 106, 034319 (2022).
  * Pecora  _et al._ [2014] L.M. Pecora, F. Sorrentino, A.M. Hagerstrom, T.E. Murphy, and R. Roy, “Cluster synchronization and isolated desynchronization in complex networks with symmetries,” Nat. Comm. 5, 1–8 (2014).
  * Zemskov and Horsthemke [2016] E.P. Zemskov and W. Horsthemke, “Diffusive instabilities in hyperbolic reaction-diffusion equations,” Phys. Rev. E 93, 032211 (2016).