The topological Dirac equation of networks and simplicial complexes

We define the topological Dirac equation describing the evolution of a topological wave function on networks or on simplicial complexes. On networks, the topological wave function describes the dynamics of topological signals or cochains, i.e. dynamical signals defined both on nodes and on links. On simplicial complexes the wave function is also defined on higher-dimensional simplices. Therefore the topological wave function satisfies a relaxed condition of locality as it acquires the same value along simplices of dimension larger than zero. The topological Dirac equation defines eigenstates whose dispersion relation is determined by the spectral properties of the Dirac (or chiral) operator defined on networks and generalized network structures including simplicial complexes and multiplex networks. On simplicial complexes the Dirac equation leads to multiple energy bands. On multiplex networks the topological Dirac equation can be generalized to distinguish between different mutlilinks leading to a natural definition of rotations of the topological spinor. The topological Dirac equation is here initially formulated on a spatial network or simplicial complex for describing the evolution of the topological wave function in continuous time. This framework is also extended to treat the topological Dirac equation on $1+d$ spaces describing a discrete space-time with one temporal dimension and $d$ spatial dimensions with $d\in \{1,2,3\}$. This work includes also the discussion of numerical results obtained by implementing the topological Dirac equation on simplicial complex models and on real simple and multiplex network data.


Introduction
Topology [1][2][3][4][5] is emerging as a very powerful tool to reveal important structural properties of interacting networks and simplicial complexes. Indeed persistence homology [6,7] can reveal large scale-properties of data that cannot be captured by more traditional network measures that characterize mostly the combinatorial properties of networks. New results [8][9][10][11][12][13][14][15][16][17][18][19] are showing that topology can also greatly enrich our understanding of dynamics. In fact networks can sustain topological signals that are dynamical variables not only defined on the nodes but also on the links. For instance topological signals defined on links, or 1-cochains, can model fluxes associated to links. On simplicial complexes, which are discrete interacting structures not only formed by nodes and links but also by higher-dimensional simplices, topological signals can be defined on simplices of every dimension including nodes, links, triangles, tetrahedra and so on. By combining algebraic topology and dynamics, important progress has been made in capturing the higher-order synchronization dynamics and the higherorder diffusion dynamics of topological signals [8][9][10]. Moreover progress as been made in formulating signal processing techniques to analyse data in the form of topological signals [14][15][16][17] Topology is receiving a surge of interest also in theoretical physics [20] and is increasingly used in high energy physics, quantum gravity and in condensed matter as well. In particular topology reveals important properties of new states of matter such as topological insulators [21] and higher-order superconductors [22]. Interestingly quantum algorithms [23] have also been proposed to speed up Topological Data Analysis.
An important question is whether topology can also be explored to investigate the properties of quantum complex networks [24].
Despite quantum complex network is a term used to indicate very different approaches, all approaches that go under this name combine network and graph theory with quantum mechanics in order to capture non-trivial effects induced by the architecture of complex networks on quantum dynamics. We can broadly classify the approaches proposed so far in two major classes. The first class of models of quantum complex networks interpret network as "physical systems". This approach includes the study of quantum critical dynamics defined on complex networks instead of the traditionally investigated lattices [25,26], or quantum dynamics on non-commutative geometry of graphs [27], and approaches that associate quantum-oscillators to the nodes of a network and couple them to a single probe whose frequency can be modulated freely [28][29][30]. The second class of approaches uses networks as "abstract computational objects" that can represent quantum states or the entanglement present in many-body systems [31][32][33][34][35][36][37][38][39]. This has lead to the definition of Von Neumann entropy of networks [33], and its generalizations based on the spectral properties of networks [35]. This class of models includes also models of networks [40] and simplicial complexes [34,39] whose statistical properties are described by quantum statistics. Finally the formulation of complex networks built from weighted matrices of mutual information characterizing interacting many-body systems [37,38] follows also in this class of models.
In this work our aim is to use algebraic topology to investigate the interplay between topology and the evolution of wave functions. We formulate the topological Dirac equation describing the dynamics of a topological wave function. On networks, the topological Dirac equation acts on a discretized wave function which is defined on each node and each link, i.e. is the combination of different topological signals (a 0-cochain and a 1-cochain) associated to the network. This implies that the discretized wave function includes dynamical variables that are defined also on links, i.e. we are ready to sacrifice slightly the notion of locality of the wave function. We have therefore that the whole wave function of the network is defined by a vector here called topological spinor which is formed by two blocks: one characterizing the value of the wave function on the nodes one characterizing the value of the wave function on the links of the network. The topological Dirac equation is then simply written by appropriately using the discrete Dirac operator (also known as chiral operator) associated to the network. We note that while in discrete topology and discrete geometry the name of Dirac operator is used to indicate different operators such as the one used in Ref. [41], here we leverage on the definition adopted in Ref. [23] to propose quantum algorithms for speeding up Topological Data Analysis.
The topological Dirac equation determines a class of quantum complex networks that uses networks both as a physical system indicating the substrate of the quantum dynamics (as the wave function is defined on the network) and as computational representation of the wave function (as the wave function is a combination of 0-cochains and 1-cochains defined on networks). In some sense the topological Dirac equation goes in the opposite direction of quantum graphs [42][43][44]: while quantum graphs describe piecewise continuous wave functions obeying the linear Schröedinger equation on each link, the topological Dirac equation describes a completely discretized wave-function taking a given (constant) value on each node and on each link.
Our aim is to explore the relation between the spectral properties of the network and the eigenstates of the topological Dirac equation.
The spectral properties of networks [45][46][47] are known to be fundamental to investigate diffusion on networks and simplicial complexes [9,10,13,15] and they provide a natural link between topology and dynamics. Characterizing the spectral properties of graphs and networks constitute and important branch of graph theory [45], network theory [46,47] and machine learning [48]. Interestingly the spectral properties of networks have recently attracted also increasing attention in the quantum gravity community and in the field of quantum networks [30,[49][50][51].
Here we reveal the relation between the topological Dirac equation and the spectral properties of simple networks, simplicial complexes and multiplex networks.
While networks are only objects formed by nodes and links generalized network structures including simplicial complexes and multilayer networks are able to encode more complex interacting patterns. Simplicial complexes [11,19] capture the many body interactions of a complex systems including not only pairwise interactions (like networks) but also higher-order interactions such as three-body interactions (triangles) and fourbody interactions (tetrehedra) and so on. Multilayer networks [52] are able to capture network of networks describing complex systems in which a given set of nodes is related by different types of interactions. As the vast majority of complex systems from the brain to technological and transportation network is formed by elements (nodes) connected via interactions that have different nature and connotation, multilayer networks are ubiquitous when describing complex systems.
The spectral properties of simplicial complexes [9, 10, 15] and multilayer networks [53] are attracting great interest in recent years. Here we show how their spectral properties affects the topological Dirac equation. On networks the topological Dirac equation describes the dynamics of wave functions with relativistic dispersion relation in which the role of the momentum is played by the eigenvalue of the graph Laplacian. On simplicial complexes the topological Dirac equation has a dispersion relation that includes more than a single energy band with each energy band determined by the spectral properties of the nth order up-Laplacian. The eigenstates corresponding to the eigenvalues in the n-th energy band are also localized exclusively on n and n + 1 simplices.
The topological Dirac equation defined on networks can be compared with the Dirac equation in one dimension, and in a lattice does not distinguishes links according to their directionality. In order to distinguish between x-links, y-links and z-links in lattices of dimension d = 2 and d = 3 we introduce here the directional topological Dirac equation. While several approaches to quantum gravity ranging from spin networks to spin foams [54][55][56][57][58] associate spins to links of the network here the approach is developed in the simple case in which the directions of the links are orthogonal leading to a very simplified treatment of the rotations induced by the different directions of the links.
Inspired by the Dirac equation in two and three dimension we can formulate a directional topological Dirac equation for duplex networks, (i.e. multiplex networks formed by two layers). This directional topological Dirac equation is able to discriminate between the different types of interactions connecting any two pair of nodes encoded by the so called multilinks [59]. This implies that the directional topological Dirac equation will treat differently multilinks of types (1, 0), (0, 1) and (1, 1) characterizing the interaction of pair of nodes connected only in layer 1, only in layer 2 or in both layers respectively. This approach leads to a natural definition of a topological spinor formed by two 0-cochains and two 1-cochains and naturally associates a different rotation of the topological spinor to different types of multilinks.
All the above results are obtained for the topological Dirac equation or its directional generalizations which are acting on a topoloigcal spinor defined on a network describing the spatial topology, while time is treated as a continuous variable. To address the interesting problem whether the topological Dirac equation can be defined on a discretized space-time. Here we treat the illustrative case of a 1 + d dimensional lattice having d ∈ {1, 2, 3} spatial dimension and one temporal dimension. Interestingly we can treat the space as a simple Euclidean space in R d+1 dimension while the directional Dirac operator can be defined differently for different directions w. In particular choosing the spatial directional Dirac operators as Hermitian and the temporal directional Dirac operator as an suitably defined anti-Hermitian matrix allows us to obtain eigenstates whose component defined on the links is an eigenvector of the discrete d'Alembert operator.
The paper is defined as follows: In Sec. In this section we introduce the topological Dirac equation defined on networks. To this end we consider a network G = (V, E) formed by a set V of N nodes (or vertices) {1, 2, . . . , N } and a set E of M links (or edges) { 1 , 2 , . . . , M } connecting the nodes by pairwise interactions. Here we assume that the network G has an arbitrary topology that might include chains and square lattices but can be used to treat also complex network topologies capturing the complexity of interacting systems ranging from the Internet to brain networks.
The topology of the network [1] can be captured by the incidence matrix B [1] representing the boundary operator of the network mapping each link of the network to a linear combination of its two endnodes. In particular the incidence matrix B [1] of a network (indicated here for ease of notation simply by B) is a rectangular matrix of size N × M having elements given by The incidence matrix B is known to define the graph Laplacian L [0] describing diffusion from node to node through links and the 1-down-Laplacian L down [1] describing diffusion from link to link though nodes: L down The Dirac operator (also called the chiral operator) [23] of a network is defined starting from the incidence matrix of the network and can be represented by a square (N + M ) × (N + M ) matrix given by where b ∈ C has absolute value |b| = 1 and where b is its complex conjugate. Typically we can choose b equal to the real unit (b = b = 1) or b equal to the imaginary unit (b = −b = i). Interestingly the square of the Dirac operator is a topological Laplacian represented by a block diagonal matrix formed by two blocks: one acting on the nodes degree of freedom and reducing to the graph Laplacian and one acting on the link degree of freedom and reducing to the 1-up-Laplacian, i.e. [1] .
In this work we want to use this Dirac operator to define a topological Dirac equation defined on a network. Quantum wave equations including the Dirac equations and the Schröedinger equations are typically defined on continuous space-time. When these equations are put on a lattice typically it is assumed that the wave function is defined on the nodes of the lattice (network). Here we assume instead that the wave function of a network is not only defined by a set of variables φ defined on the nodes of the network but includes also a set of variable χ defined on the links of the network. Therefore we relax a bit the notion of point-like definition of the wave function and consider also degree of freedoms associated to links. Therefore we consider the wave function described by the topological spinor with φ indicating a 0-cochain and χ indicating a 1-cochain associated to the network, i.e.
where a 0-cochain and a 1-cochain indicate functions defined on nodes {1, 2, . . . N } and links { 1 , 2 , . . . , M } respectively. While in Sec. 5 we will touch on the problem of modelling a discrete Lorentzian space-time, to start with we consider that the network, although of arbitrary topology only captures the discrete space, while we keep time as a continuous variable. The topological Dirac equation is defined as the wave-equation where the Hamiltonian of the topological Dirac equation is taken to be Here m 0 ≥ 0 is a parameter of the dynamics we call the mass and the matrix β is a (N + L) × (N + L) block diagonal matrix of elements where I P indicates the identity matrix of dimension P . Interestingly we note that the anti-commutator between D and β vanishes, i.e.
and that Let us now consider the eigenstates at constant energy, satisfying Let us write this equation explicitly by considering the part of the wave function φ defined on nodes and the part of the wave function χ defined on links. In this way we obtain, or equivalently From this equation it can be shown that φ should be an eigenvector of the graph Laplacian L [0] and χ should be an eigenvector of the 1-down-Laplacian L down [1] . In fact with simple manipulations of Eqs. (15) we obtain By indicating with λ the generic eigenvalue of B and with λ the corresponding eigenvalue of B † this equation implies the energy E satisfy the relativistic dispersion relation This means that there are positive and negative energy eigenstates corresponding to energy values E = ± |λ| 2 + m 2 . As long as the mass is strictly positive, i.e. m 0 > 0 the eigenfunctions associated to the positive and negative solutions are given by ψ λ,+ and ψ λ,− respectively. These eigenstates are given by where w R λ and w L λ are the right and left eigenvectors of the incidence matrix B associated to the eigenvalue λ, and where we adopted the normalization condition ψ λ,± |ψ λ,± >= ψ λ,± |β|ψ λ,± >= ±1, similar to the one used for the massive Dirac equation [60]. In Figure 1 and Figure 2 we represent examples of eigenstates of the topological Dirac equation corresponding to positive and negative energy values E obtained starting from real network datasets. As it is clear from the figures and from the analytical expression of the eigenstates, these eigenstates are strongly depending on the topology of the underlying network and on its spectral properties. Let us now consider the density of states associated to the  dispersion relation Eq. (18). We observe that µ = |λ| 2 indicates the eigenvalue of the graph Laplacian L [0] . Indicating by ρ(µ) the density of positive eigenvalue of the graph Laplacian we can obtain the density g(E) of positive energy eigenvalues E > m 0 from the relation Using E = µ + m 2 0 we therefore obtain In particular this implies that if the network has a spectral dimension d S [46,47], i.e. if ρ(µ) scales like then g(E) scales as Since the non-zero eigenvalues of the graph Laplacian L [0] are the same eigenvalues of the 1-down Laplacian L down [1] (and with the same degeneracy) the eigenstates at negative energy E = − |λ| 2 + m 2 0 have the same degeneracy of the states of positive energy E = |λ| 2 + m 2 0 as long as |E| > m 0 . Therefore we have that the density of negative energy eigenstates is given by g(−E) as long as |E| > m 0 . In Figure 3 we display the cumulative density of eigenvalues G(E) indicating the fraction of positive eigenvalues with energy greater than m 0 and smaller than E for the real networks investigated in Figure 1 and Figure 2. The density of states at energy E = m 0 and energy E = −m 0 deserves some particular attention. These energy states correspond to the eigenvalue µ = 0 of the graph Laplacian L [0] and of the 1-down Laplacian L down [1] . These two matrices have the same non-zero eigenvalues. In particular the positive eigenvalues coincide together with their degeneracy. However the matrix L [0] is a N × N matrix and the matrix L [1] is a M × M matrix therefore the zero eigenvalue has a different degeneracy for the two Laplacians. If we assume, without loss of generality that the network is connected, then the number of links M satisfies M ≥ N − 1 where M = N − 1 implies that the network is a tree and L = N implies that the network is a ring (or chain with periodic boundary conditions). An important result of spectral graph theory is that for a connected network the graph Laplacian L [0] has a zero eigenvalue µ = 0 with multiplicity 1 (corresponding to the fact that the network has a single connected component). The L down [1] Laplacian of the same network will have N − 1 non-zero eigenvalues while the multiplicity of the zero eigenvalue will be M −(N −1). In other words the multiplicity of the zero eigenvalue of the 1-down-Laplacian L down [1] is equal to the number of independent cycles in the network.

Topological Dirac equation on a one dimensional chain
The topological Dirac equation can be defined on a network of arbitrary topology, however it is instructive to investigate its properties on a 1-dimensional chain. For a 1-dimensional chain of size N with periodic boundary conditions we have that the number of links is equal to the number of nodes, (i.e. M = N ) and that the eigenvalues of the graph Laplacian µ can be expressed as where k x = 2πn x /N indicates the wave-number and 0 ≤n x ≤ N − 1. For k x 1 we can approximate the expression for µ with The density of eigenvalues ρ(µ) ∝ k −1/2 x and the degeneracy of the eigenstates with positive energy E = m 0 and negative energy E = −m 0 is the same. The eigenstates of the topological Dirac equation are formed by a 0-cochain φ and a 1-cochain χ that are proportional to the left and right eigenvector of the incidence matrix B forming a 1 dimensional Fourier basis with wave-number k x . In this case the topological Dirac equation is most similar to the Dirac equation defined in one dimension [21] with the most notable difference that the spinor is given a topological interpretation. However the topological Dirac equation applied to higher dimensional lattices does not distinguish between the the discrete gradient (performed by the application of the incidence matrix B) performed in different orthogonal directions. In this respect the topological Dirac equation differs significantly from the Dirac equation in d = 2 and d = 3 dimensions. In Sec. 3 we will investigate how the topological Dirac equation can be modified to get closer to the Dirac equation in dimensions d = 2 and d = 3 taking into account the different directions of the links. Before we explore this generalization, however, we devote the next section to treat the topological Dirac equation on simplicial complexes.

Topological Dirac equation of simplicial complexes
Simplicial complexes are generalized network structures which explicitly indicate the existence of higher-order interactions. Indeed simplicial complexes do not only include nodes and links but also include higher-order simplices such as triangles, tetrahedra and so on. A n-dimensional simplex α is a set of n + 1 nodes A simplicial complex K is a set of simplices closed under the inclusion of subsets. Therefore if a simplex α belongs to the simplicial complex, i.e. α ∈ K any simplex α ⊂ α including a subset of the nodes of α also belongs to the simplicial complex, i.e. α ∈ K. Here an in the following we indicate with N [n] the number of n-dimensional simplices of the considered simplicial complex. Therefore N [0] = N indicates the number of nodes of the simplicial complex and N [1] = M indicates the number of links of the simplicial complex. Simplicial complexes describe discrete topological spaces and are the fundamental objects studied by algebraic topology. Interestingly, when we associate a length to the links of the simplices, simplicial complexes can also be used to study arbitrary discrete geometries, and as such they have been extensively used in quantum gravity to describe discrete (or discretized) space-time. Here we build on algebraic topology to formulate a higher-order version of the topological Dirac equation on simplicial complexes. In algebraic topology [1], simplices are also given an orientation, typically chosen to be induced by the node labels, so that for instance a link [i, j] has positive orientation if j > i and negative otherwise. Oriented n-dimensional simplices can be chosen as the base of the space of n-chains C n whose elements are linear combinations of ndimensional simplices with coefficients in Z (or R). The definition of n-chains allows us to explore the topology of the simplicial complex algebraically. An important operator in algebraic topology is the boundary operator ∂ n : C n → C n−1 that maps each n-simplex to the (n − 1)-chain representing the (n − 1)-dimensional simplices at its boundary. The boundary operator ∂ n indeed acts on any simplex α given by Eq. (28) as It follows that for n = 1 every link is mapped to its two endnodes, i.e.
and for n = 2 every triangle is mapped to the three links at its boundary with the correct orientation, i.e.
Let us indicate with B [n] the matrix that represents the n-boundary operator. We notice immediately that B [1] = B indicates the incidence matrix of a network as defined in Eq.
(1) where this definition applies also to higher dimensional simplicial complexes. The incidence matrices B [n] have the notable property that their subsequent concatenation is null, i.e. for any n > 1 This algebraic relation reflects the important topological property that the boundary of the boundary is null. The higher-order incidence matrices can be used to define the higher-order Laplacians L [n] . For n = 0 the higher-order Laplacian reduces to the graph Laplacian L [0] defined in Eq. (2), for n > 0 the higher-order Laplacian is given by with L down and L [n] commute and therefore can be simultaneously diagonalized. Moreover thanks to the Hodge decomposition [1] any non-zero eigenvector of L [n] is either a non-zero eigenvector of L up [n] or a non-zero eigenvector of L down [n] . Given this important insight for algebraic topology we can generalize the topological Dirac equation to higher-order topologies. In particular we can define the n-order Dirac operator with n ≥ 1 as an operator represented by a (N [n−1] + N [n] ) square matrix given by where b [n] ∈ C has absolute value |b [n] | = 1 and where b [n] is its complex conjugate. Therefore we have that for n = 1 the n-order Dirac operator reduces to the Dirac operator D defined in the previous paragraph. The square of the n-order Dirac operator leads to the Laplacian operator acting on (n − 1)-dimensional simplices and n-dimensional simplices and describing diffusion among (n − 1)-simplices through nsimplices and diffusion among n simplices through (n − 1)-simplices, i.e.
Using the n-order Dirac operator we can define the n-order topological Dirac equation of a simplicial complex. This equation is defined on a topological spinor ψ [n] including a (n − 1)-cochain χ [n−1] and a n-cochain χ [n] , i.e. the topological spinor ψ [n] is a (N [n−1] + N [n] )-dimensional column vector given by According to the definition used in this section χ [0] is the 0-cochain that can be identified with the 0-cochain indicated in the previous section by φ. Similarly χ [1] can be identified with χ used in the previous section. The n-order topological Dirac equation on simplicial complex, describes the evolution of the n-order topological spinor ψ [n] according to the equation where the Hamiltonian H [n] , given by is defined in terms of the n-order Dirac operator D [n] and the matrix β [n] given by By indicating with λ n the generic eigenvalue of the boundary operator B [n] and proceeding in analizing the eigenstate of the Hamiltonian using steps similar to the ones used in the previous paragraph, one can easily see that the energy E [n] associated to the eigenstates of H [n] satisfies the dispersion relation As long as the mass m 0 is strictly positive the eigenstates associated the the positive energy states (ψ λ,+ ) and the ones associated with the negative energy states (ψ λ,− ) can be normalized as the eigenstates of the topological Dirac equation on networks and they are given by where w R [n]λ and w L [n]λ are the right and left eigenvectors of the boundary operator B [n] associated to the eigenvalue λ. Interestingly, it is also possible to consider a higher-order topological Dirac equation defined on a spinor including all possible n-cochains defined on the simplicial complex. Therefore, if the simplicial complex is d dimensional, i.e. if the simplicial complex includes simplices of dimensions up to n = d, we define the topological spinor [1] . . .
and we consider the higher-order topological Dirac equation given by with HamiltonianH [d] given bỹ where now we have lifted the dimension of D [n] and the dimension of β [n] to the dimension of the spinor Ψ [d] . This leads to the matrix representation for D [n] and β [n] given by Note that from Eq. (32) it follows that the n-order Dirac operators obey for any n = n . Therefore we have with L being the block diagonal matrix with diagonal blocks given by the n-order Laplacian, i.e.
Interestingly this latter formulation of the higher-order Dirac equation on simplicial complexes admits eigenstates that are exclusively localized on (n − 1)-dimensional and n-dimensional simplices for 0 < n ≤ d where d is the dimension of the simplicial complex. To be concrete we can consider this higher-order topological Dirac equation on a simplicial complex of dimension d = 3. In this case the HamiltonianH [3] reads The eigenstate Ψ [3] = (χ [0] , χ [1] , χ [2] , χ [3] ) corresponding to the energy E satisfies the system of equations This system of equations can be easily generalized to arbitrary dimension d. By investigating this equation it is possible to see that this system of equation admits d positive and d negative energy bands E [n] with the nth (positive and negative) energy bands depending only on the eigenvalues of the n-incidence matrices λ n and the mass m 0 . We have in particular that the energy dispersion relation of the different bands is where for small |λ 1 | m 0 and small |λ d | m 0 we have Interestingly for m 0 = 0 all the energy bands have relativistic dispersion relation For each band, the density of states with |E In particular for the case in which the mass is null, m 0 = 0, we have for all 0 < n ≤ d. It has been recently pointed out that geometrical simplicial complexes can admit differebr higher-order spectral dimensions d  In this case we obtain that the higher-order Dirac equation admits d positive and d negative energy bands with density of states g [n] obeying In Figure 4 we show the multiband spectrum of the simplicial complex model of emergent hyperbolic geometry called "Network Geometry with Flavor" [34,39] showing that the different bands can be associated to different spectral dimensions. The eigenstates associated to the n-th energy band are localized on (n − 1) and n-dimensional cochains. Possibly by introducing opportune interaction terms one could induce transitions among states in different bands and even hybridizations of these different eigenstates. We finish this section by noticing that the topological Dirac equation can be also extended to treat hypegraphs in which hyperedges are directed, i.e. the hyperedges have input and output nodes such as in chemical reaction networks. Indeed for such hypegraphs the incidence matrix has been recently defined in Ref. [13] which can easily allow the generalization of the topological Dirac operator and Dirac equation to hypegraphs. This is a very interesting extension of the topological Dirac equation, however due to space limitation we omit this discussion in the present paper.

Directional topological Dirac equation in d = 2-dimensional lattices
In general network topologies describing discrete spaces with non trivial local fluctuations of the curvature it might be non trivial to define the equivalent to the orthogonal directions of lattices. For this reason the topological Dirac equation introduced in Sec. 1 only includes the β matrix but not the equivalent of all the γ matrices present in the d = 3 dimensional Dirac equation. However for 2-dimensional and 3-dimensional lattices the orthogonal spatial directions are naturally defined. Therefore we can explore whether it is possible to define a variation of the topological Dirac equation which distinguishes between the different spatial directions.
To this end we consider a square portion of a square lattice with sides of length √ N with periodic boundary condition (a torus) and we classify links of the d = 2 square lattice according to their direction (see Figure 5). For instance we consider x-type links the links connecting adjacent nodes along the x-direction and y-type links connecting adjacent nodes along the y-direction. Therefore a x-type link will connect a node i of coordinates r i = (x i , y i ) to the nodes j of coordinates r j = r i ± e x modulo periodic boundary conditions, where e x = (1, 0). Similarly a y-type link will connect a node i of coordinates r i = (x i , y i ) to the nodes j of coordinates r j = r i ± e y modulo periodic boundary conditions, where e y = (0, 1). Given this classification of different types of links, distinguished according to their direction, the incidence matrix B of a d = 2 dimensional lattice can be written as where B (w) with w ∈ {x, y} describes the incidence matrix between nodes and links of type w, i.e.
, and is a w-type link 0 otherwise.
Introducing this decomposition of the incidence matrix in contributions coming from links of orthogonal directions opens the possibility to decompose also the Dirac operation accordingly and to introduce some additional rotation in the topological space of nodes and links. To this end, for lattices of dimension d = 2 we can define the directional Dirac operatorD as with D (w) describing the contribution of w-type links to the directional Dirac operator. We can proceed by considering rotations due to the different directions of the links. Therefore, inspired by the standard Dirac equation in dimension d = 2 [21] we define D (x) and D (y) as With this choice of the directional Dirac operator, we consider always the topological spinor ψ = (φ, χ) of dimension N + M defined on Eq. (6) and we consider a modified topological Dirac equation given by with HamiltonianH taken to bē whereD is the directional Dirac operator in two dimension given by Eq. (60) and the matrix β is given by Eq. (10). Note however a very notable consequence of this definition of directional topological Dirac equation aimed at treating differently orthogonal space directions: while the anti-commutator between D (w) and the matrix β is zeros for both w ∈ {x, y}, i.e.
the anti-commutator between D (x) and D (y) is given by the not zero matrix This is a notable difference with respect to the Dirac equation where the corresponding anti-commutator vanishes as the derivative respect to x and to y commute. Despite this difference it can be shown that the directional Dirac operator has eigenfunctions that satisfy the relativistic dispersion relation where λ = (λ x , λ y ) with λ x and λ y representing the eigenvalues of the directional incidence matrices B (x) and B (y) respectively. In particular we have that λ x and λ y can be expressed in terms of the wave-number k = (k x , k y ) as where k x = 2πn x / √ N and k y = 2πn y / √ N withn x ,n y integers between zero and N 1/2 − 1. In order to show this results we can decompose the link component χ of the topological spinor in two blocks, χ = (χ (x) , χ (y) ) where χ (w) indicates the wave function calculated on links of type w. The wave function at constant energy E then satisfies Eq. (13) which can be rewritten as Simple algebraic operations show that φ should be an eigenvector of the graph Laplacian, with the graph Laplacian L [0] given by where We note that L (x) and L (y) commute, i.e.
therefore the eigenstates of the directional Dirac equation have a dispersion relation where µ is the generic eigenvalue of the graph Laplacian L [0] that can be written as It follows that for |k| 1, i.e. µ 1 the density of eigenvalues ρ(µ) of the graph Laplacian follows Eq. (24)  x + λ 2 y + m 2 0 has elements where k = (k x , k y ) and where N is a normalization constant, while the components χ (x) and χ (y) satisfy Here we note that since L [0] commutes with L (w) the component χ of the eigenstates of the directional topological Dirac equation is also eigenvector of the L down [1],(x) and the L down [1],(y) Laplacians given by L down

Directional topological Dirac equation in d = 3 dimensional lattices
In this section we consider a cubic portion of a three dimensional lattice with sides of length N 1/3 and periodic boundary conditions. We distinguish between links of type x, y and z according to their direction (see Figure 5). In particular two nodes i and j of coordinates r i = (x i , y i , z i ) and r j = (x j , y j , z j ) are connected by a w-type link if and only if r j = r i ± e (w) modulo periodic boundary conditions, where e (x) = (1, 0, 0), e (y) = (0, 1, 0), and e (z) = (0, 0, 1). The incidence matrix B of the three dimensional lattice can be written as where B (w) indicates the directional incidence matrix between w-type links and nodes that admits the same definition as in two dimensions (Eq. (105)). It turns out that for lattices in dimension d = 3 the topological spinor of dimension N + M is not sufficient to distinguish between the three orthogonal spatial directions. Indeed, if we want to define the directional topological Dirac operator that distinguishes between the three different types of links we need to consider the topological spinorΨ of dimension 2(N + M ) given bỹ where Φ is a vector of two 0-cochains φ (r) with r ∈ {1, 2} defined on the nodes of the network and X is a vector of two 1-cochains χ (r) with r ∈ {1, 2} defined on the links of the network, i.e.
Having defined the directional incidence matrix and the topological spinor, the natural next step is to define a directional Dirac operator given bȳ where the w-directional Dirac operators D (w) define independent rotations of the topological spinor. The w-directional Dirac operator D (w) has a block structure given by where B (w) is a matrix of dimension 2N × 2M having a different definition depending on the type of the link. If we define the "Pauli" matrices σ(F) as where F is a generic matrix of dimension N × M , we can express B (w) for w ∈ {x, y, z} as Having defined the directional Dirac operatorD in dimension d = 3 (given by Eq.(81)), we can formulate the directional topological Dirac equation in d = 3 as with HamiltonianH given bȳ where β is now the matrix We notice that β anti-commutes with each D (w) , i.e.
for w ∈ {x, y, z} and the anti-commutators between the directional d = 3-dimensional Dirac operators associated to different types of links w = w are given by the matrices where [B (w) ] † B (w ) + [B (w ) ] † B (w) can be expressed as with ijs indicating the Levi-Civita symbol. Similarly to what we have discussed for the directional topological Dirac equation in dimension d = 2, also the eigenstates of the directional topological Dirac equation in dimension d = 3 obey the dispersion relation where µ is the generic eigenvalue of the graph Laplacian L [0] given by Here the w-type directional graph Laplacian L (w) is given by with each pair of directional graph Laplacian commuting with each other, i.e.
By indicating with λ w the generic eigenvalue of the incidence matrix B (w) we have with where k = (k x , k y , k z ) indicates the wave-number and where due to the periodic boundary conditions k w takes the values k w = 2πn w /N 1/3 withn w integer between 0 and N 1/3 − 1. It follows that for |k| 1, i.e. µ 1 the density of eigenvalues ρ(µ) of the graph Laplacian follows Eq. (24) and the density g(E) of eigenstates of positive energy E follows Eq. (25) with d S = d = 3. In order to show that the eigenstates of the directional topological Dirac equation in d = 3 follow the dispersion relation given by Eq. (91), we indicate the component X of the eigenstates as X = (X (x) , X (y) , X (z) ) where X (w) is the component of X calculated on links of type w. With this choice of notation, the directional topological Dirac equations has eigenstates associated to the energy E which satisfy From this equation it follows that Φ satisfies where L [0] is given by and each directional L (w) can be written as with r ∈ {1, 2}) form a d = 3-dimensional Fourier basis with associated wave-number k = (k x , k y , k z ) and for each value of the energy E = ± |λ x | 2 + |λ y | 2 + |λ z | 2 + m 2 0 have elements where N is a normalization constant. The component can be evaluated as with w ∈ {x, y, z}. Interestingly since L [0] commutes with every directional L (w) , we obtain that X (r) obtained in this way are also eigenvectors of the L down [1],(w) Laplacians given by for w ∈ {x, y, z}.

Topological Dirac equation of multiplex networks
Interestingly the directional topological Dirac operator developed to treat the different types of links in d = 3 dimensional lattices, can be adapted to investigate the topological Dirac equation on multiplex networks [52] as well. Here we consider the case of a duplex network, i.e. a multiplex network formed by two layers. Each layer of the network describes a different network of interactions, between the same set of nodes. For instance in the brain network of C. elegans the two layer can represent electrical gap junctions or chemical synapses between the neurons. In a single network each pair of nodes can be either connected or not connected by a link. In a duplex network instead each pair of nodes can be connected in multiple ways [59]: they can be connected only in layer 1, they can be connected only in layer 2 or they can be connected in both layers. If they are only connected in layer 1, we say that they are connected by a multilink m = (1, 0), if they are only connected in layer 2 we say that they are connected by a multilink m = (0, 1), if they are connected in both layers we say they are connected by a multilink (1, 1) (see Figure 5). Note that each pair of nodes can be connected only by a single type of multilink. This decomposition of multilinks of different types can be compared with the treatment of different types of links of the d = 3 dimensional lattice discussed in the previous section. This allows us to extend to duplex networks the formalism developed to treat the directional topological Dirac equation for d = 3 dimensional lattices. Let us indicate with N the total number of nodes of the duplex network and with M the total number of multilinks of any type m ∈ {(1, 0), (0, 1), (1, 1)}. In order to treat the directional topological Dirac equation on duplex networks, we decompose the incidence matrix B of size N × M in three terms depending on the three types of multilinks where here B m with m ∈ {(1, 0), (0, 1), (1, 1)} describes the incidence matrix between nodes and multilinks of type m, i.e.
This decomposition allows us to consider the directional Dirac operatorD where the m directional Dirac operator is defined as where µ is the eigenvalue of the graph Laplacian L [0] given by with L m given by However for duplex networks, we have that a very notable difference with respect to the d = 3 dimensional lattices. Indeed the Laplacians corresponding to distinct multilinks in the vast majority of cases do not commute, i.e.
[L m , L m ] = 0, if m = m. Therefore the eigenstates (Φ (r) , X (r) ) with r ∈ {1, 2}, corresponding to the energy E = ± µ + m 2 0 have the Φ (r) component with elements  [63] with the first layer including collaborations or works published in the arxiv reporsitory physics.soc-ph while the second layer includes collaborations of work published in the arxiv repository cond-mat.dis-nn. The total number of nodes is N = 4537, the total number of pairs of nodes connected in both layers is O = 650.
(1,0) , X However since the graph Laplacians associated to different types of multilinks do not commute, i.e. Eq. (112) applies, the component X of the eigenstates of the directional topological Dirac equation on duplex networks is not an eigenvector of the Laplacians for any type of multilink m ∈ {(1, 0), (0, 1), (1, 1)}. In Table 1, Table 2 and Table 3 we evaluate to what extent the multi-Laplacians L m associated to different multilinks of real multiplex networks do not commute. The metric used to evaluate the deviation from zero of the commutator [L m , L m ] is the Euclidean metric given by the square-root of the normalized sum of the squares of the elements of the matrix. The data that we consider includes the collaboration network of scientists working in Network Science, the C.elegans duplex connectome, the duplex network of the Arabidopsis genetic layers. The Euclidean metric measured on the commutators [L m , L m ] display a wide range of values with the C.elegans duplex network displaying the largest values indicating the most significant deviation from the commutating scenario.  Table 2. The Euclidean metric of the commutator between the Laplacian L m and the Laplacian L m for a duplex networks of gap-junction and synaptic connections in C. elegans.. The duplex network is constructed from the dataset published in [64] and curated in [65] Table 3. The Euclidean metric of the commutator between the Laplacian L m and the Laplacian L m for a duplex networks capturing Arabidopsis genetic layers. The duplex network is constructed from the dataset published in [66] and curated in [67] with the first layer including direct interactions the second layer physical association. The total number of nodes is N = 6903, the total number of pairs of nodes connected in both layers is O = 628.
(1, 0) 0 0.620 0.120 (0, 1) 0.620 0 0.076 (1,1) 0.120 0.076 0 which also time is discretized. To start investigating this question we consider a 1 + d dimensional lattice with d ∈ {1, 2} representing the entire space-time with d dimensions indicating space and an additional dimension indicating time, so that the links will be distinguished between x-type, y-type links and t-type links. For simplicity we assume the space is isotropic portion of the lattice in R d+1 with sides of length N 1/(d+1) and we consider period boundary conditions. For 1 + 1 lattice incidence matrix B can be decomposed as for 1 + 2 lattice the incidence matrix can be decomposed as where B (w) with w ∈ {x, y, t} indicates the directional incidence matrix defined by Eq. (105) . Starting from these directional incidence matrices we can define the directional Dirac operatorD as where for 1 + 1 lattices the sum extends to w ∈ {x, y} and for 1 + 2 lattices it extends to w ∈ {x, y, t}. Here the w-directional Dirac operators D (w) corresponding to the spatial directions w ∈ {x, y} have the same defintion as for d dimensional spatial lattices, and are given by Eq. (65) with the directional Dirac operator D (t) corresponding to the temporal direction is defined as We notice that while D (x) and D (y) are Hermitian, we have chosen D (t) to be anti-Hermitian. With this choice of the t-directional Dirac operator D (t) , adopting for β the definition given by Eq. (10), one observes that the anti-commutator between β and D (w) vanishes, i.e.
for w ∈ {x, y, t}. However the anti-commutator between D (w) with w ∈ {x, y} and D (t) is given by .
The directional topological Dirac equation on this 1 + d space-time is given by the eigenvalue problem D + m 0 β ψ = 0 (122) with the topological spinor ψ = (φ, χ) defined in Eq. (6). Let us study this equation in the case of a 1 + 1 dimensional lattice (as the case 1 + 2 is a straightforward generalization). By writing this eigenvalue problem for the component φ and the component χ separately, where the component χ is decomposed according to the type of link in which the 1-cochain is defined, χ = (χ (x) , χ (t) ) we obtain From this system of equations it is possible to observe that φ must be an eigenvector of the discrete D'Alembert operator with eigenvalue m 2 0 . Therefore the component φ of the topological spinor calculated on a node i of coordinates (t i , x i ) is given by where N is a normalization constant. Here ω = 2πn ω /N 1/2 , k x = 2πn x / √ N with 0 ≤n w ≤ N 1/2 − 1. Let us indicate with λ x and E t the eigenvalues of the directional incidence matrices B (x) and B (t) with with for w ∈ {x, y}. In order to satisfy the directional topological Dirac equation on 1 + 1 space-time E t and λ x must satisfy the dispersion relation which, given the discrete nature of the spectrum can impose constraints on the possible value of the mass m 0 and eigenvalues of B t and B x that might be considered. If this relation can be satisfied the eigenstate exists, and the topological Dirac equation on discrete space time has a solution with the components χ (x) and χ (y) given by By following similar steps it is possible to show that for the 1 + 2 topological Dirac equation, the component φ is an eigenvector of the operator with eigenvalue m 2 0 and its element φ i calculated on a node i of coordinates (t i , x i , y i ) is given by where N is a normalization constant and k w = 2πn w /N 1/3 with 0 ≤n w ≤ N 1/3 − 1.
Indicating with λ x , λ y and E t the eigenvalues of the directional incidence matrices B (x) , B (y) and B (t) respectively we see that for 1 + 2 dimensional space-time this eigenvalues must satisfy the dispersion relation Therefore, due to the discrete nature of the spectrum, solving for the spectrum of the directional topological Dirac equation reduces to a problem connected to number theory.

Directional Topological Dirac equation on 1 + 3 dimensional space-time lattice
In this section we further extend the directional topological Dirac equation to 1 + 3 dimensions. To this end we consider the topological spinor Ψ = (Φ, X) defined in Eq. (79). The incidence matrix B can be written as a sum of directional incidence matrix B (w) with w ∈ {x, y, z, t}. The directional Dirac operators D (w) with w ∈ {x, y, z, } have the same definition as for d = 3 spatial lattices given by Eq. (82). The directional Dirac operator in the time direction is defined instead as with B (t) given by We notice that if we define β according to Eq. (87), we have Moreover the anti-commutators between the spatial Dirac operators and the temporal Dirac operator are given by can be expressed as where σ(F) are given by Eq. (83). The directional topological Dirac equation on a 1 + 3 dimensional space-time can therefore we written as where the directional Dirac operatorD is given bȳ The study of this equation shows that the topological spinor Ψ = (Φ, X) has the component Φ which is an eigenvector of the discrete D'Alembert operator with eigenvalue m 2 0 , where L (w) are given by for w ∈ {t, x, y, z}. The eigenvectors of are eigenvector of the spatial B (w) with eigenvalues λ w and of the temporal B (t) with eigenvalue E t where E t and λ w which must satisfy where λ = (λ x , λ y , λ z ). Given that we have assumed to have a finite 1 + 3 space-time with periodic boundary conditions, λ w and E t can only take discrete values fixed by the pair (ω, k) with with k w and ω quantized and given by k w = 2πn w /N 1/4 , ω = 2πn ω /N 1/4 withn w and n ω taking integer values between zero and N 1/4 −1. The component Φ (r) with r ∈ {1, 2} of the eigenstates of energy E t has elements Φ (1) where we have indicated with k = (k x , k y , k z ) and components X (r) = X (r) (x) , X (r) (y) , X (r) (z) , X (r) (t) given by for the spatial directions w ∈ {x, y, z} and for the temporal direction t, X (t) are:

Conclusions
In this work we have formulated the topological Dirac equation and its directional variations and we have studied the properties of this wave equation on networks (simple and multilayer) and on simplicial complexes. The main difference between the Dirac equation and the topological Dirac equation is that the wave function of the topological Dirac equation is discretized and defined on nodes, links and higher dimensional simplices. Therefore the topological spinor has a strict topological interpretation and relaxes the requirement of the locality of the wave function to the extent that it includes for instance 1-cochains defined on the links of the network. We have shown that the topological Dirac equation obeys a relativistic dispersion relation in which the eigenvalue of the Dirac operator plays the role of the momentum. The topological Dirac operator has been extended to treat also simplicial complexes. In this context we observe that the eigenfunctions of the higher-order topological Dirac operator can be localized only on n and n − 1 dimensional simplices and that correspondingly the dispersion relation is described by independent energy bands associated to the spectrum of the n-order up Laplacian. While the topological Dirac equation is defined on an arbitrary topology, we have considered a variation of this equation, the directional Dirac equation on d = 2 and d = 3 dimensional lattices and multiplex networks. The directional equation treats interactions of different types differently introducing a rotation of the spinor that is a function of the type of the interaction. While for lattices the type of interaction distinguishes between the different orthogonal directions of the links on multiplex networks the interactions of different types indicate different types of multilinks. Therefore we see that if we want to distinguish between different multilinks of multiplex networks the directional topological Dirac equation naturally leads to rotations of the topological spinor.
Finally we have discussed the extension of this work to treat a discrete network describing the entire space-time topology over which the directional topological Dirac equation can be defined. In particular we formulate the directional topological Dirac equation defined on 1 + d lattices with d ∈ {1, 2, 3} spatial dimensions. This sheds some light on the notoriously difficult problem of defining a Lorentzian metric for a discrete geometry. Interestingly here the lattice is treated as an Euclidean lattice in R 1+d and the Loretzian nature of the lattice is only due to the opportune definition of the directional spatial and temporal Dirac operators.
This work can be expanded in different directions. From a Network Science perspective, it is interesting to study also the dissipative dynamics described by the Dirac operator, this would entail making a Wick rotation of the time in the topological Dirac equation. Research in this direction would be important to investigate whether the Dirac operator could be used for signal processing of different topological signals simultaneously and whether it could be explored as an important operator to define convolution of simplicial neural networks.
Many other research questions remain open if one desires to push forward the mapping between this mathematical framework and physics systems. One important question that arises is whether one can observe symmetry breaking and the emergence of a renormalized mass from the study of non-linear topological Dirac equations similarly to what happens for the non-linear Dirac equation captured by the Nambu Jona-Lasinio model [68]. In this case the renormalized mass might acquire a different value depending on the spectral properties of the network or of the simplicial complex expanding our understanding of the relation between discrete topology, geometry and dynamics. Another important direction that could be explored is the introduction of appropriate topological fields that might induce interaction terms able to describe transitions between different eigenstates of the topological Dirac operator. These terms could eventually lead to transition across different bands of the topological Dirac operator defined on simplicial complexes. Finally the directional topological Dirac equation defined on discrete-space time could be extended to general discrete networks and simplicial complex geometries.
We hope that these results will inspire further applied and theoretical work along these lines revealing the rich interplay between the topology of networks and simplicial complexes and the evolution of wave functions and topological signals defined on them.

Data availability.
The fungi data published in [61] is freely available at the repository [62]. The multiplex network datasets used in this work are all freely available at the repository [69].

Code availability.
The code to generate the simplicial complex "Network Geometry with Flavor" [34,39] is freely available at the repository [70]. All other codes used in this work are available upon request.