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

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.


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 [28,[38][39][40][41][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 nanomaterials [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 where u ∈ C 1 is 0-cochain here encoded by a vector taking real values on each node of the network and v ∈ C 1 is a 1-cochain here encoded by a vector taking real values on each link of the network, i.e.
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 δ 1 : C 0 → 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 = δ 1 u); in the considered example we have Moreover the discrete divergence δ ⋆ 1 : C 1 → 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 = δ ⋆ 1 w) where in the considered case we have for all node r 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 × N 1 boundary matrix B whose elements are given by In particular we have f = Bw and g = B ⊤ u.On a network the Hodge Laplacians L [0] and 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 The most simple definition of the Hodge-Dirac operator [34,35] on a network is the operator ∂ : C 0 ⊕ C 1 → C 0 ⊕ C 1 that maps topological spinors to topological spinors and is defined as or in matrix form we have The Hodge-Dirac operator can be considered as the "square root" of the higher-order Laplacian L. Indeed we have Therefore the non-zero eigenvalues of the Dirac operator are given by plus or minus the square root of the non-zero eigenvalues of L [0] (which is by the way isospectral to 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 ∂ anticommutes (i.e., {∂, γ 0 } = 0) with γ 0 , where where I X indicates the X × X identity matrix.This result implies that if (u, w) ⊤ is an eigenvector of ∂ with eigenvalue λ, then (u, −w) ⊤ is an eigenvector of ∂ with eigenvalue − λ.This relation between the nonharmonic 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.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 Φ ∈ C 0 ⊕ C 0 ⊕ 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 where χ ∈ C 0 ⊕ C 0 is defined on nodes and ψ ∈ C 1 is defined on links.In particular we will use the notation with u ∈ C 0 , v ∈ C 0 indicating the two node signals and w ∈ C 1 indicating the link signal.For ease of notation here and in the following we indicate with N = 2N 0 + N 1 and with M = N 0 + 2N 1 .We have therefore that Φ is the N dimensional column vector, given by The 3-way Hodge-Dirac operator ∂ : or alternatively as In contrast with the Hodge-Dirac operator, the Dirac operator D / : Therefore the gamma matrix γ 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 γ is defined in terms of 2-dimensional column vector α = (α u , α v ) ⊤ and a 2-dimensional row vector β = (β u , β v ), as or alternatively It follows that the 3-way Dirac operator can be also expressed as As previously shown in Eq. ( 7) 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 with It is instructive to consider the case in which α = (1, 1) ⊤ and β = (1, 1).In this case we have On a side note we observe that the Dirac operator defined in Eq.( 19) is reminiscent of the topological operator used for studying duplex networks with link overlap in Ref. [70].
Given the diagonal block structure of L, its eigenvalues are either the eigenvalues of L uv or the ones of From this property it follows that the eigenvalues of the Dirac operator are the square roots of the eigenvalues of 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 D / commutes with the matrix γ 0 given by Hence if a spinor Φ = (u, v, w) ⊤ is an eigenvector of the three way Dirac operator with eigenvalue λ ̸ = 0, i.e. if then γ 0 Φ = (u, v, −w) ⊤ is an eigenvector of the three way Dirac operator with eigenvalue − λ.Indeed we have From these results it follows that the harmonic eigenvectors of the three way Dirac operator are where χ harm indicates the generic harmonic eigenvector of L uv while ψ harm the generic harmonic eigenvector of L w .
Note that the harmonic eigenvectors χ harm read Where χ0 are the harmonic right singular vectors of the boundary matrix, and where χ is an arbitrary N 0 dimensional vector.The non zero eigenvalues λk of the 3-way Dirac operator are given by where µ k are the singular eigenvalues of the boundary operator.Note that with a general choice of the gamma matrix γ the eigenvalues of the Dirac operator are real only if and becomes purely imaginary (i.e., appears in complex conjugate pairs) for ∆ < 0 while they are identically zero in the trivial case ∆ = 0.
The eigenvectors ϕ + , ϕ − associated to the non-zero eigenvalues λk with the same absolute value are related by chirality that for the three-way Dirac operator takes the form where χµ are the right singular vectors and ψ µ are the left singular vectors corresponding to a non-zero singular value µ = µ k of the boundary operator.
Let us now discuss the spectral properties of the Dirac operator D / = D / 0 defined in Eq. ( 19) in the case in which α = (1, 1) ⊤ and β = (1, 1).As long as the network is connected the zero eigenvalue of L has multiplicity N 1 + 2 while the non-zero eigenvalues are given by twice the non-zero eigenvalues µ of L [0] (or ). Therefore the eigenvalues of the Dirac operator D / 0 in this case are real and given by λk = 0 with multiplicity

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 Φ follows the dynamics where F(Φ, ∂Φ) 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 ∂Φ vanishes when there are no interactions between simplices of different dimensions.Nonetheless, this formalism allows us to consider also nonlinear terms in B and B ⊤ , 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 Φ = (χ, ψ) ⊤ , where χ = (u, v) ⊤ , encodes for the two node signals and ψ = 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 ∂Φ is instead given by In other words χ is the one dimensional projected signal on the nodes and ψ is the two dimensional projected Accordingly, our choice of the reaction term F(Φ, ∂Φ) is 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 with c 1 ∈ R, c 2 ∈ 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. ( 34) explicitly, obtaining the following 3-species reaction-diffusion system 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 B and B ⊤ 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 F(Φ) 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 1 is in the kernel of the Dirac operator [24,40].
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 ∂Φ ⋆ = (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 hence, we can study the local stability of this homogeneous state considering nodes and links separately.
Let us assume Φ * = (u * , v * , w * ) ⊤ to be a fixed point for the above system, meaning that w * is a fixed point of h(0, 0, w) and (u * , v * ) ⊤ 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 From the assumption that (u * , v * ) ⊤ 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 * ) ⊤ is a stable equilibrium for system (36) provided that Let us now perturb system (35) with a spatially inhomogeneous perturbation (δu, δv, δw) ⊤ 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 where, for sake of notation, we have considered the terms in c 1 of Eq. ( 35) as part of the reaction functions, We can proceed as in [40] and use the eigenvectors of L [0] and L [1] to perform the singular value decomposition of B. Let us recall that L [0] and L [1] are, in this case, isospectral and their eigenvalues µ 2 k are the square of the singular values µ k of B. By projecting the perturbations δu and δv on the eigenbasis of L [0] , and the perturbation δw on the eigenbasis of L [1] , Eq. ( 38) becomes where δ ûk (resp.δv k , δ ŵk ) is the perturbation expressed in the new basis.To study the linear stability of the perturbation, we hypothesize that δ ûk (t), δv k (t), δ ŵk (t) ∼ e λ k t .To ensure the existence of a nontrivial solution, the linear growth rate λ k must satisfy the following condition where |J k | denotes the determinant of the matrix .
We can rewrite this condition as a 3-rd order polynomial in the variable λ k of the form whose coefficients are where π 1 and π 2 are polynomials in µ 2 k , whose explicit expressions can be found in Appendix A. By means of the Routh-Hurwitz criterion [71,72], we can study the stability of polynomial (41).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 λ the maximum real part of the λ k , and its associate imaginary part ϱ.Here λ and ϱ are the real and imaginary part of the dispersion relation, respectively, and they are both functions of the continuous parameter µ k .The condition for having Turing instability is thus the existence of a µ k for which λ(µ k ) > 0. Let us observe that being the support discrete, there can be finite size effects, namely λ(µ k ) < 0 for all k while the continuous function λ 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 λ(µ k ) > 0 and the corresponding imaginary part is non-zero, i.e., ϱ(µ k ) ̸ = 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 ϱ implies also the presence of a solution with −ϱ.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 ϱ [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.

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] 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) ⊗ (1 N0 , 1 N0 , 1 N1 ) 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 (42) 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 L [0] admits (1 N0 , 1 N0 ) 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 2D lattice with periodic boundary condition of 36 nodes, i.e., 6 × 6.
The conditions (37) ensuring the stability of the homogeneous equilibrium become in the present case As shown in Appendix 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 × 6 square lattice with periodic boundary conditions.In Fig. 2 In both bifurcation diagrams the red curves, i.e., σ 1 +σ 2 = 0 and σ 1 σ 2 = ξ 1 ξ 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 ϱ, i.e., the imaginary part associated to unstable mode λ, 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)).We now compare the dynamical patterns observed in presence of diffusion (c 2 ̸ = 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 µ k = 0 is stable (indeed λ(0) < 0), while there exist some modes that are unstable (λ(µ 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 ϱ 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 network 1 .
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 µ 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 µ k can be found in Appendices Appendix A and Appendix B. By considering the obtained dispersion relation we note that also in this case, as in the previous one, its imaginary part ϱ 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 2D lattice 2 .Remarkably, Dirac-induced patterns present several distinct dynamical properties that are not present in Turing pattern observed for c 2 ̸ = 0.
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 1 Such dynamics can be better appreciated in the video available as Supplementary Material of the online version. 2 The dynamics can be better appreciated in the video available as Supplementary Material of the online version.patterns exclusively driven by the Dirac operator.Also in this case, we consider a set of parameters such that the stability of the homogeneous equilibrium is lost once perturbed, as revealed by the real part of the dispersion relation in Fig. 4a).We can also remark that the latter is divergent.When working on continuous support, where all instability modes are present, this is a problem because the maximum critical mode as an infinite wave number, hence the solution is not a physical one.However, on discrete support not every possible instability mode is present, and on any finite network the spectrum of the Laplacian is bound.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 the two dynamical patterns, Laplace-driven and Dirac-induced, respectively, let us project the dynamics in the 3D 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, ℓ, we consider one of its end-node, i, and we build the orbit (u i (t), v i (t), w ℓ (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 3D 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.

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 ̸ = 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 1 in the kernel of the Dirac operator.
In order to investigate the dynamical properties of the model (42) 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 dispersion relation is computed as for the case with diffusion (Fig. 2, 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 (λ > 0) and they are oscillatory (the corresponding ϱ ̸ = 0), as shown by the dynamics of the 3 interacting species, in panel c).
The network is a 2D lattice with periodic boundary conditions of 36 nodes and the parameters are 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 with that of Turing patterns, which are highly sensible to even small variations on the initial conditions.

Phenomenology on arandom graph
We further investigate the nature and the phenomenology of the patterns by analyzing the case of a random Eulerian network, shown in Fig.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 (λ > 0) and they are oscillatory (the corresponding ϱ ̸ = 0), as shown by the dynamics of the 3 interacting species, in panel c).
The network is a 2D lattice with periodic boundary conditions of 36 nodes and the parameters are the initial perturbation is the order of 10 −2 .
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)) 3 .
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   path and we orient the links according to the node label, then it is immediate to show that the network admits a constant eigenvector 1 in the kernel of the Dirac operator.
In order to investigate the dynamical properties of the model (42) on this random topology, we distinguish between Turing patterns and Dirac-induced patterns.In both cases we consider parameter values for which the homogeneous equilibrium becomes unstable and (in principle) oscillatory.In Fig. 9  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 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.
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.

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  affect jointly nodes and links signals with possible applications to neuroscience.whose coefficients are where π 1 and π 2 are polynomials in µ 2 k .
The first polynomial has degree 2 and has the following coefficients 1 (µ while the second one has degree 3 and coefficients are π

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.
whose coefficients are , where π 1 and π 2 are polynomials in µ 2 k .
The first polynomial has degree 2 and has the following coefficients while the second one has degree 3 and coefficients are π Let us observe that The Routh-Hurwitz criterion [71,72] gives us the following necessary and sufficient conditions for stability 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 π 2 (µ 2 k ) = Fµ 4 k + Gµ 2 k + H can take negative values.We can observe that 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 F < 0 (so that we are sure that π 2 takes some negative values, independently of G) or G < 0 (π 2 may be negative, independently of 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 ).Then, it follows that λ 0 satisfies A(λ 0 ) = 0 and the second term λ 1 can be obtained by solving A ′ (λ 0 )λ 1 + B(λ 0 ) = 0.Because A is a linear function, its root can be straightforwardly determined This implies that for sufficiently large µ k the roots of (A1) will be close to λ 0 within a term of order 1/µ 2 k .

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: where and ϱ 1 and ϱ 2 are now linear curves in µ 2 k .
The first one has the following form ϱ 1 (µ while the second one is ϱ 2 (µ 2 k ) = Zµ 2 k + K, explicitly: The Routh-Hurwitz criterion [71,72] gives us the following necessary and sufficient conditions for stability a > 0, b > 0, bc − ad > 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 ϱ 1 < 0. Since the coefficient K is always positive, the condition reduces to imposing Z negative 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.
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 Φ ∈ C 0 ⊕ C 1 given by

Figure 1 :
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.
acts on the topological spinor Φ projecting independently the two node signals encoded by χ into the links and the single link signal ψ into the nodes.We define the 3-way Hodge-Dirac operator ∂ as the M × N matrix given by 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 N × M gamma matrix γ, i.e.D / = γ∂ .
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 χ can interact with each other and with the projected signals χ and similarly on the links the projected signals on the links ψ can interact with each other and with the link signal ψ.
, we report the bifurcation diagram of the model in the parameter space (σ 1 , σ 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).
(a-b), we plot the dispersion relation corresponding to our choice of parameter values and, in panel c), we display the resulting 17 a) b)

µ k ≫ 1 . 4 k
The first step is to rewrite the latter as followsp(λ) = µ 4 k A(λ) + µ 2 k B(λ) + C(λ) ,where A, B and C are polynomials in λ of degree 1, 2 and 3.By dividing the previous expression by µ we obtain a second polynomial in λ with the same roots as the previous one.Let us first consider µ k ≫ 1 and write the root of the polynomial as λ(µ k ) = λ 0 + λ 1 /µ 2 k + O(µ −4 k ) 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 ϱ 3 (µ 2 k ) = J µ 2 k + L, where now      J = bX − Z = (f u + g
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 timevarying 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