Nonlinear diffusion on networks: Perturbations and consensus dynamics

1.   Introduction

Diffusion is certainly among the most renowned subjects in dynamical systems. The ordinary differential equation $\dot{\mathit{\boldsymbol{x}}} = - L \mathit{\boldsymbol{x}}$, where $L$ is the graph Laplacian, describes a (linear) diffusion process on a network. In fact, the Laplacian matrix plays the role of a discrete counterpart of the Laplacian differential operator, $\Delta$, present for example in the heat equation, $\dot u = \Delta u$. A lot of effort has been made to precisely analyze the behavior of Laplacian systems, especially for the linear case ^[48,57,58]. More recently, nonlinear extensions of the Laplacian dynamics have been explored ^{[1,3,31,32,55]}.

In the present paper, we aim to analyze a class of nonlinear Laplacian systems defined by

where $\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{x}}(t)\in\mathbb{R}^n$, $F(\mathit{\boldsymbol{x}})$ is a nonlinear vector field, $0 < \epsilon \ll 1$ is a small parameter, and $H(\mathit{\boldsymbol{x}}, \Lambda)$ is a perturbation with $\Lambda\in\mathbb R^p$ an extra set of parameters; more precise definitions follow in the next sections. Eq (1.1) describes a nonlinear diffusion process on a network. In the literature, such systems are also known as Absolute Laplacian Flows (ALFs) ^[3,55]. Similarities with the continuum case can be found by looking at drift-diffusion equations ^[11], in particular with the porous-medium equation, $\dot u = \Delta u^m$ ^[56]. The porous-medium equation has several potential applications, for example: flow of gases ^[38,41], theoretical biology ^[54], and particle-based models ^[28] (where the power considered is $m = 2$). A network model for the porous-medium equation has been proposed in ^[19]; similar systems have also been studied in ^[9,47]. Let us moreover notice that by considering the rewriting $F(\mathit{\boldsymbol{x}}) = : K(\mathit{\boldsymbol{x}}) \mathit{\boldsymbol{x}}$ one can interpret the matrix $K(\mathit{\boldsymbol{x}})$ as a state-dependent diffusion coefficient, leading, in turn, to another interpretation in terms of the network structure. In fact, we can define a state-dependent Laplacian, $L_{ \text{K}}(\mathit{\boldsymbol{x}}): = L K(\mathit{\boldsymbol{x}})$, which can be seen either as the Laplacian of a directed graph under the definition of out-degree Laplacian ^[2,46], or as a particular case of the Laplacian tensor for consensus models on hypergraphs ^[52]. Nonlinear diffusion processes have also been considered in discrete time systems, in particular for lattice networks ^[45].

Within the context of Laplacian dynamics the diffusion of 'information' among the agents can lead the system to reach a common state, referred to as consensus. As such, Laplacian systems find a remarkable application in consensus problems ^[15,48,50], which have been extensively studied for their importance in several branches of science, e.g., synchronization, opinion dynamics, coordination of robots, rendezvous problems, among many others.

One significant challenge in studying networked dynamical systems is the high dimensionality of the problem. Furthermore, the nonlinearities inherent in such models contribute to an even greater degree of complexity. The contributions of this paper address such challenges in particular contexts where the model we consider is a nonlinear version of the well-known consensus dynamics. For complete graphs, under appropriate conditions accounting for the symmetries of the network system, we provide a reduction to a planar system, see Lemma 4. Hence, in this particular setup, a large dimensional nonlinear system can be studied with low-dimensional techniques of dynamical systems. For example, by exploiting the aforementioned reduction, we prove that, generically, the singularities along the consensus space are transcritical; see Proposition 8. Later, for nondegenerate transcritical singularities, we provide conditions for the existence of ''canard solutions'' in Proposition 10. These canards are relevant in the context of perturbation problems since they are solutions that are not expected from zeroth-order approximations. Although the previous results are for complete graphs, we also notice that, under some conditions, the canard condition can be provided for more general graph structures; see Proposition 11. Finally, the insights that we gain from our analysis are used in some numerical experiments showing some intriguing phenomena that we associate with the canards; see Section 6.2.

The paper is structured as follows: In Section 2, we define the (unperturbed) ALF and exploit some general properties. Section 3 is dedicated to understanding the role of symmetries in ALFs; for a brief introduction to the main concepts related to symmetries see Appendix 8. In Section 4, we introduce a perturbation term leading to Eq (1.1). As we will see, the perturbations, in general, give rise to a slow-fast behavior of the system. A brief summary of the main ideas of singular perturbation theory useful for the main analysis can be found in Appendix 9. The slow-fast consensus dynamics of ALFs with complete graph structure is studied in Section 5, where we also provide conditions under which systems with arbitrary graph structure admit a maximal canard. In Section 6, we explore numerically different network structures, spatio-temporal configurations, and we describe the effects of canard conditions in different setups. A brief discussion and the conclusions are drawn in Section 7.

2.   Definitions and preliminaries

From a mathematical point of view, an undirected network structure can be described by a graph, $\mathcal{G} = \{ \mathcal{V}, \mathcal{E} \}$, which consists of a finite set of nodes (or vertices) $\mathcal{V} = \{ 1, ..., n \}$, and a set of edges $e_{ij} \in \mathcal{E}$, where each $e_{ij}$ represents a connection between the node $j$ and $i$ ^[24]. Furthermore, if we assign a set of weights $\mathcal{W}$ to the graph $\mathcal{G}$, then we obtain a weighted graph, $\mathcal{G}_ \text{w} = \{ \mathcal{V}, \mathcal{E}, \mathcal{W} \}$. We assume (unless otherwise stated) that the weights, $w_{ij} \in \mathcal{W}$, are positive reals. Such a choice is natural both for the algebraic properties that follow from it ^[5], and for the modeling perspective, since often the weights represent some measure of distance or strengths of the coupling ^[42]. In other contexts, negative weights are also considered ^[53], for example in models representing inhibition or antagonism; in this paper we do not consider such situations. The edges and the weights are in one-to-one correspondence, in fact, we assign to each edge $e_{ij}$ a weight $w_{ij}$. Throughout this paper, we consider only undirected graphs. Therefore, we have that the weights are symmetric, that is $w_{ij} = w_{ji}$.

Definition 1. A simple graph is an undirected graph without self-loops or multiple edges.

Unless otherwise stated, we consider simple graphs. For weighted graphs, we assume a simple unweighted structure, i.e., $\mathcal{G}_ \text{w} \setminus \mathcal{W}$ is a simple graph. Moreover, we make no nomenclature distinction between graphs and weighted graphs with unit weights, because the algebraic structures in such cases are exactly the same.

We now recall a few well-known concepts of graph theory ^[7,24,42]. A path is a sequence of edges that joins a sequence of distinct vertices. A connected component of a graph is a sub-graph such that every pair of nodes is connected by a path. A relevant class of graphs is given by the complete graphs.

Definition 2. A complete graph with $n$ nodes, denoted by $K_n$, is a graph such that for all pairs of nodes there is an edge connecting them, i.e., $\forall i, j \in \mathcal{V}$, with $i \neq j$, $\exists \ e_{ij} \in \mathcal{E}$.

Given a graph it is possible to define some algebraic structures on it. We recall here the ones that are useful for the purpose of the paper.

● The adjacency matrix, $A$, of a graph is the $n \times n$ matrix with components $A_{ij} = w_{ij}$ if $e_{ij} \in \mathcal{E}$, and $A_{ij} = 0$ if $e_{ij} \notin \mathcal{E}$.

● The degree matrix, $\Delta$, of a graph is the $n \times n$ diagonal matrix with (diagonal) components $\Delta_{ii} = \sum_j w_{ij}$.

● The Laplacian matrix, $L$, of a graph is the $n \times n$ matrix defined by $L : = \Delta - A$.

For simple graphs, let us recall that the Laplacian matrix, $L$, is symmetric, degenerate and semi-positive definite; moreover, the following proposition holds ^[5].

Proposition 1. The algebraic multiplicity of the $0$ eigenvalue of $L$, $\mu_a(0)$, is equal to the number of connected components of the graph. Furthermore, the vector $\mathbf{1} = (1, ..., 1)^\intercal$ is an eigenvector of $L$ with eigenvalue $0$, i.e., $\mathbf{1} \in \text{ker}(L)$.

We now introduce the class of systems we are going to investigate in the present paper. Let us assign to each node of a graph, $i \in \mathcal{V}$, a state $x_i \in \mathbb{R}$ and a smooth function

which we call the response function. So, given the set of nodes $\mathcal{V}$, we obtain a state vector $\mathit{\boldsymbol{x}} = (x_1, ..., x_n)^\intercal$, and a response vector field $F(\mathit{\boldsymbol{x}}) = \left(f(x_1), \dots, f(x_n) \right)^\intercal$. An Absolute Laplacian Flow ^[3,55] is defined by the equation

Remark 1. Due to Proposition 1, the Laplacian decomposes into a direct sum $L = L^{(1)} \oplus \dots \oplus L^{(\mu_a(0))}$, and, therefore, an ALF decomposes in $\mu_a(0)$ independent systems. So, in general, we consider connected graphs.

Note that if we choose the response vector field to be the identity, i.e., $F(\mathit{\boldsymbol{x}}) = \mathit{\boldsymbol{x}}$, then Eq (2.1) becomes the widely studied (linear) Laplacian dynamics, $\dot{\mathit{\boldsymbol{x}}} = - L \mathit{\boldsymbol{x}}$, ^[48,57,58]. Indeed, the ALF is a direct nonlinear generalization of the linear Laplacian dynamics ^[1,3,31,55]. ALFs have received considerably less attention than the linear Laplacian dynamics, although they have strong relations with consensus problems. For example, ALFs have been considered as models for nonlinear communications protocols in ^[43,61].

An important observation is that, in contrast to linear Laplacian dynamics and depending on the response function, ALFs may exhibit bifurcations of the consensus state, see (2.3). Therefore, our objective is to characterize the dynamic behavior near consensus for such systems under small perturbations.

2.1.   General properties

ALFs have some generic properties derived from the symmetric and degenerate structure of the Laplacian matrix. First, notice that since $\mathbf{1}\in\ker L$, then for any $h\in \mathcal{C}^m(\mathbb{R}^n, \mathbb{R})$, $m\geq3$, one has that $L(F(\mathit{\boldsymbol{x}})+h(\mathit{\boldsymbol{x}})\mathit{\boldsymbol{1}}) = LF(\mathit{\boldsymbol{x}})$, meaning that ALFs are invariant under response field transformations $F(\mathit{\boldsymbol{x}}) \mapsto F(\mathit{\boldsymbol{x}}) + h(\mathit{\boldsymbol{x}}) \mathbf{1}$. As a consequence, we have that the response vector field is not unique: in fact, there exists an infinite set of vector fields leading to the same dynamics. This freedom can be used in practice to choose the "simplest" response field. Such a property reminds in some way of a gauge freedom ^[37], where one can choose among a class of response fields $F^{(h)}(\mathit{\boldsymbol{x}}) = F(\mathit{\boldsymbol{x}}) + h(\mathit{\boldsymbol{x}}) \mathbf{1}$ without affecting the evolution equation (2.1). Another important property of ALFs is the existence of a constant of motion.

Lemma 1. Given an ALF (2.1), there exists a constant of motion $k$,

where $\langle \cdot, \cdot \rangle$ is the usual Euclidean scalar product in $\mathbb{R}^n$.

Proof. The time evolution of $k$ is given by $\dot{k} = \langle \mathbf{1}, \dot{\mathit{\boldsymbol{x}}} \rangle$. Then we have $\langle \mathbf{1}, \dot{\mathit{\boldsymbol{x}}} \rangle = - \mathbf{1}^\intercal L F(\mathit{\boldsymbol{x}}) = - (L \mathbf{1})^\intercal F(\mathit{\boldsymbol{x}}) = \; 0$. Thus, $k$ is a constant (in time) value.

It is well known that the constant of motion $k$ is related to the arithmetic mean of the states, $k = n \langle \mathit{\boldsymbol{x}} \rangle$, where $\langle \mathit{\boldsymbol{x}} \rangle : = \langle \mathbf{1}, \mathit{\boldsymbol{x}} \rangle /n$. Thus, we can interpret the presence of the constant of motion $k$ as a conservation law for $\langle \mathit{\boldsymbol{x}} \rangle$. More generally, every quantity of the form $c k$, with $c$ a constant, is a constant of motion for an ALF.

Lemma 2. The set of equilibria of an ALF (2.1) is given by

Proof. In order to have an equilibrium for Eq (2.1) we need $F(\mathit{\boldsymbol{x}}) \in \ker(L)$, which means $f(x_1) = \dots = f(x_n)$, implying in turn that $\mathit{\boldsymbol{x}}$ is an equilibrium if and only if $\mathit{\boldsymbol{x}} \in E$.

Corollary 1. The set

is a subset of equilibria of system (2.1), i.e., $C \subseteq E$.

The set $C$ is the so-called consensus space ^[32,48,50]. As we already mentioned, one of our goals is to characterize the behavior near the consensus space (2.3). Our characterization aims to be qualitative, in the sense of topological equivalence ^[27,62].

Lemma 3. The consensus space can be parametrized by the constant k, i.e.,

Proof. Follows from straightforward computations.

From the previous lemma, we denote a parametrized point in $C$ by $\mathit{\boldsymbol{x}}^*_k = (x_k^*, \dots, x_k^*)^\intercal$, where

Notice that, since the system has a constant of motion, the evolution is actually restricted to a family of $(n-1)$-dimensional hyperplanes

It is known that for linear Laplacian dynamics, systems associated with weighted graphs are topologically equivalent to the corresponding unweighted ones ^[8]. An important observation is that ALFs on weighted and unweighted graphs are equivalent near consensus on $P_k$.

Definition 3. Two vector fields, $X, \; Y$, are said to be topologically equivalent if there exists a homeomorphism which takes orbits of $X$ to orbits of $Y$, preserving directions but not necessarily time parametrisation.

Proposition 2. For any $f$ and $k$ fixed, in a small neighborhood of $\mathit{\boldsymbol{x}}^*_k \in C$ and restricted to $P_k$, the weighted, $\dot{\mathit{\boldsymbol{x}}} = -L_ \text{w}F(\mathit{\boldsymbol{x}})$, and the unweighted, $\dot{\mathit{\boldsymbol{x}}} = -LF(\mathit{\boldsymbol{x}})$, ALFs are topologically equivalent.

Proof. Let $L_ \text{w}$ be the Laplacian of a weighted graph $\mathcal{G}_ \text{w}$. Then the Jacobian of the associated ALF restricted to the consensus space is given by $- L_ \text{w} \text{D}F(\mathit{\boldsymbol{x}}^*_k)$, where $\text{D}F(\mathit{\boldsymbol{x}}^*_k)$ is the Jacobian of $F(\mathit{\boldsymbol{x}})$ evaluated on the consensus set parametrized by $k$. On the other hand, by setting all the weights to one on $\mathcal{G}_ \text{w}$ we obtain the simple graph $\mathcal{G}$. Thus, for the ALF with unweighted graph $\mathcal{G}$ the Jacobian is $- L \text{D}F(\mathit{\boldsymbol{x}}^*_k)$, where $L$ is the Laplacian matrix associated to $\mathcal{G}$. We notice that $\text{D}F(\mathit{\boldsymbol{x}}^*_k) = \text{Id } \text{d}_xf(x^*_k)$, where $\text{d}_x f$ is a short-hand notation for the derivative of $f$ with respect to $x$, and $\text{Id}$ is the $n\times n$ identity matrix. So, the two matrices to compare are $- \text{d}_xf(x^*_k) L_ \text{w}$ and $- \text{d}_xf(x^*_k) L$. We remark that $- \text{d}_xf(x^*_k)$ is a common scalar factor and that the matrices $L_ \text{w}$ and $L$ do not depend on $x^*_k$.

Let us recall that the Laplacian matrix is a semi-positive definite matrix ^[5], which means that its eigenvalues are all greater than or equal to zero. Such a property holds for both simple graphs and positive-weighted graphs. Moreover, since in our context $L_ \text{w}$ and $L$ are symmetric, we can transform (by a similarity transformation) the Jacobians to the diagonal forms

where $\{0, \mu_1, \dots, \mu_{n-1}\}$ are the eigenvalues of $L_ \text{w}$, and $\{0, \lambda_1, \dots, \lambda_{n-1}\}$ are the eigenvalues of $L$; in both cases the non-zero eigenvalues are positive. Moreover, since the $0$ eigenvalue is associated to the one-dimensional set of equilibria $C$ (Proposition 1 and Corollary 1), we can reduce the analysis to the hyperplane $P_k$, which leads to the $(n-1)$-dimensional diagonal matrices

where $x^*_k$ is now a parameter. Whenever the common term $\text{d}_xf(x^*_k)$ is nonzero, we have that both matrices have eigenvalues with the same sign. Therefore, topological equivalence follows from the Hartman-Grobman Theorem. Actually, it is possible to explicitly compute the homeomorphism between the two systems described by Eq (2.5). Since the matrices are diagonal, effectively, we compute $n-1$ homeomorphisms between one-dimensional systems. The flows of the $i$-th components are respectively

Thus, the homeomorphism, $z:\mathbb{R} \to \mathbb{R}$, such that $z(x_{\mu_i}(t, x_0)) = x_{\lambda_i}(t, z(x_0))$, is given by

It is straightforward to check that Eq (2.6) is a homeomorphism of the real line, as it is also proven in ^[29]. Notice that the transformation (2.6) does not depend on the common factor of the systems (2.5), and therefore, it holds also when $- \text{d}_xf(x^*_k) = 0$.

Remark 2. Having in mind Proposition 2, from now on we assume that ALFs have all weights equal to one, i.e., they are defined by a simple (unweighted) graph $\mathcal{G}$.

Notice that the topological equivalence described by Proposition 2, defines a local equivalence close to the consensus space. However, in what follows the assumption of a simple graph structure is global. This is motivated by the later analysis that is dedicated to the consensus space, which is the focus of the paper. The reader should keep in mind that away of the consensus space the equivalence is not proven, and, therefore, does not apply in general. The results we stated until this point rely mainly on the algebraic properties of the Laplacian. For simple graphs, together with Proposition 2, it is particularly useful to exploit the symmetry properties induced by the graph, which in turn can provide insightful information for the ALF.

3.   Symmetries

In the following sections, we are going to look at how the symmetry properties of the response function and of the graph affect the equilibria and the dynamics of ALFs. Some fundamental notions of group theory and its applications are contained in Appendix A.

3.1.   Equilibria of ALFs

Corollary 1 shows that the consensus space is always a subset of equilibria for (2.1). Now, we show that, starting from the consensus space, the group transformations leaving the response function invariant induce new equilibria.

Proposition 3. Let $f$ be a $\Gamma$-invariant response function. Then, the group $\Gamma_n : = \underbrace{\Gamma\times\cdots\times\Gamma}_{n \mathit{\text{-times}}}$ maps the equilibria coinciding with the consensus space, $C$, to equilibria, i.e.,

Proof. Let us consider $-L F(\Gamma_n \mathit{\boldsymbol{x}}^*)$, where $\mathit{\boldsymbol{x}}^* \in C$. We recall that by definition every component of the vector field $F$ has the same response function $f$. Also, each response function is invariant under $\Gamma$, i.e., $f(\Gamma x) = f(x)$. So we have that $F(\Gamma_n \mathit{\boldsymbol{x}}^*) = (f(\Gamma x_1^*), \dots, f(\Gamma x_n^*))^\intercal = (f(x_1^*), \dots, f(x_n^*))^\intercal = F(\mathit{\boldsymbol{x}}^*)$. Therefore $-L F(\Gamma_n \mathit{\boldsymbol{x}}^*) = -L F(\mathit{\boldsymbol{x}}^*) = 0$.

Remark 3. From Proposition 3 it follows that any subgroup of $\Gamma_n$ maps $C$ to equilibria.

The invariance of the response function leads to a richer structure of the equilibria, similarly to what happens for equivariant bifurcations ^[12,26]. Such a richer structure can give rise, for example, to bipartite consensus ^[3,39,64] where the agents converge to two separate clusters of consensus. A simple exemplification of the appearance of bipartite consensus in ALFs is given by considering the response function $f(x) = x^2$. It is clear that the symmetry group for such a function is $\mathbb{Z}_2$, in particular we are considering the representation $\{ 1, -1 \}$. From Proposition 3 we have that the set of equilibria consists of all the possible combinations $\pm x_1 = \pm x_2 = \dots = \pm x_n$, and, therefore, we have bipartite consensus among the potential equilibria of the system. Let us notice that the linear set-up for bipartite consensus requires negative weights on the network structure in order to obtain such a result, while for ALFs bipartite consensus follows from nonlinearity, and symmetry properties of the response function. It is interesting to notice that the lines of bipartite consensus, $\pm x_1 = \pm x_2 = \dots = \pm x_n$, arising from generic even response functions are saddles; this property follows from straightforward computations of the Jacobian.

3.2.   Equivariant ALFs

Since we are considering systems with a network structure, the symmetry group associated with the equivariance properties, i.e., the automorphism group of the graph, denoted by $\text{Aut} (\mathcal{G})$, is a subgroup of the symmetric group over $n$ elements, denoted by $\mathfrak{S}_n$. Moreover, if an ALF is $\mathfrak{S}_n \supseteq \Gamma$-equivariant, it follows that it is also equivariant under any subgroup $\Sigma \subset \Gamma$. Therefore, $\text{Fix}(\Sigma)$ is also a flow-invariant set, and $\text{Fix}(\Gamma)\subseteq \text{Fix}(\Sigma)$. From these arguments, if we are considering a $\Gamma$-equivariant ALF, the analysis of the dynamics can be reduced to the fixed-point spaces $\text{Fix}(\Sigma)$, for all subgroups $\Sigma$ of $\Gamma$. Considering that the phase space of an ALF is $\mathbb{R}^n$, when acting on a state, $\mathit{\boldsymbol{x}}$, we use the standard/permutation representation of $\mathfrak{S}_n$ in $\mathbb{R}^n$. An element of the standard representation is called a permutation matrix.

Definition 4. A permutation matrix $\sigma \in \mathfrak{S}_n$ is an $n \times n$ matrix with components

where the expression $\sigma(i) = j$ is the abstract action of the group element $\sigma$ on an element of the vertex set $\mathcal{V} = \{1, \dots, n\}$, meaning that the component $\sigma_{ij}$ is $1$ if the vertex $i$ is permuted with the vertex $j$, and is zero otherwise.

From now on, if not otherwise explicitly mentioned, when we write the action of a permutation group, we implicitly consider the standard representation.

Proposition 4. $\mathcal{G}$-ALFs are $\text{Aut}(\mathcal{G})$-equivariant, where $\text{Aut}(\mathcal{G})$ is a subgroup of $\mathfrak{S}_n$ in the standard representation.

Proof. Let $\sigma \in \text{Aut}(\mathcal{G})$ be an $n \times n$ permutation matrix. We start by applying a permutation to the vector field, $\sigma L F(\mathit{\boldsymbol{x}})$. Let us notice that $[L, \sigma] = 0$ for all $\sigma \in \text{Aut}(\mathcal{G})$, where $[\cdot, \cdot]$ is the commutator, as it follows from the definition of symmetry for a graph (Definition 9, Appendix A). Thus, we have that $\sigma L F(\mathit{\boldsymbol{x}}) = L \sigma F(\mathit{\boldsymbol{x}})$. Notice that a permutation applied to the vector field $F(\mathit{\boldsymbol{x}}) = (f(x_1), \dots, f(x_n))^\intercal$ exchanges the order of the components, but since the functions are indistinguishable we have that $\sigma F(\mathit{\boldsymbol{x}}) = F(\sigma \mathit{\boldsymbol{x}})$. Therefore, $\sigma L F(\mathit{\boldsymbol{x}}) = L F(\sigma \mathit{\boldsymbol{x}})$, which proves the statement.

Remark 4. Note that the action of $\text{Aut}(\mathcal{G})\subseteq\mathfrak{S}_n$ on the (nonlinear) vector field $F(\mathit{\boldsymbol{x}})$ is linear.

Proposition 5. Let $\text{Aut}(\mathcal{G})$ be an irreducible subgroup of $\mathfrak{S}_n$ in the standard representation. If $\text{Ord}(\text{Aut}(\mathcal{G})) \geq n$ then $\text{Fix}(\text{Aut}(\mathcal{G})) = C$.

Proof. Let $\text{Ord}(\text{Aut}(\mathcal{G})) = k \geq n$. The fixed-point space is obtained by considering the equations

where $\sigma^{(s)} \in \text{Aut}(\mathcal{G})$, $s = 1, \dots, k$. Suppose, by contradiction, that Eq (3.1) does not fix all the components to the same value, but $m < n$ components are fixed to one common value and the other $n-m$ components are fixed to a different common value, i.e., Eq (3.1) leads to

As a consequence, we have that the permutation matrices $\sigma^{(s)}$ can be written as follows:

which is a block form. This would imply that the representation is decomposable and therefore not irreducible, which is a contradiction. The case where the components are fixed to a number of distinct values $l > 2$ is analogous, and would give rise to block matrices with $l$ blocks. In turn, such block matrices would be associated to a decomposable representation, leading again to a contradiction.

We have established a link between ALFs and equivariant systems. Let us notice that, if we have a dynamical system equivariant under a finite group, such dynamical system can be seen as a cell-network ^[25]. A cell-network is a system of differential equations symmetric under a group(oid) associated to a network that could have directed and multiple edges. So, there is a strict correspondence between cell-networks and the ALF framework. For example, a $K_n$-ALF is a cell-network with complete graph structure, $K_n$. This alternative point of view provides an extra tool in the analysis of ALFs. Indeed, from Proposition 12 it follows that $\text{Fix}(\text{Aut}(\mathcal{G}))$ is flow-invariant, and so, in turn, the consensus space is an invariant space for a $K_n$-ALF, or in the language of cell-networks it is a pattern of synchrony, which is already known from Corollary 1. Actually, Corollary 1 not only holds for all ALFs, but also tells us that the consensus space consists of equilibria of the system, and is not just an invariant set. However, Proposition 5 does not rely on the algebraic structure of the Laplacian, it is a general result for a class of equivariant systems. Let us notice that the results following from equivariance are, in general, more robust than the algebraic ones. Therefore, the symmetry perspective could be very useful when some transformation is applied to the system.

Let us notice that, when the order of $\text{Aut}(\mathcal{G})$ is less than $n$, we have $\dim \text{Fix}(\text{Aut}(\mathcal{G})) > 1$, so the fixed-point space is related to clustering, and $C \subset \text{Fix}(\text{Aut}(\mathcal{G}))$; we illustrate such a case in the following example. We consider the path graph $P_n$ with $n > 2$, which has $\mathbb{Z}_2$ as group of symmetry and $\text{Ord}(\mathbb{Z}_2) = 2$. We can see that the equalities imposed by $\mathbb{Z}_2 \mathit{\boldsymbol{x}} = \mathit{\boldsymbol{x}}$ are

which do not lead to consensus in general, but to clustering $\{ \mathit{\boldsymbol{x}} \in \mathbb{R}^n \ | \ x_1 = x_n, \dots, x_n = x_1 \}$.

In summary, in this section we outlined a different perspective on the consensus space. Let us recall that a weighted ALF, near the consensus space, is equivalent to an ALF with simple graph structure (Proposition 2). Thus, the weighted ALF acquires the property of equivariance under the automorphism group of the (simple) graph; the symmetry group is a subgroup of the symmetric group. Therefore, an ALF possesses also a cell-network structure. At this point, forgetting momentarily the specific properties of ALFs, we can study the connections between the consensus space and the symmetry properties. It turns out that the consensus space arises as an invariant space for a class of subgroups of $\mathfrak{S}_n$ (Proposition 5).

We conclude this section by recalling a theorem highlighting the fundamental role of the graph structure in the definition of a dynamic network. Although, as we have seen, graphs and groups are intimately related, if we consider the logic flow of definitions for dynamic networks the graph structure should be defined first.

Theorem 1 (Frucht ^[22]). For any finite group $\Gamma$, there exists a finite graph $\mathcal{G}$ such that $\text{Aut}(\mathcal{G}) \cong \Gamma$.

Theorem 1 ensures that for any given finite group we can construct a graph which is invariant under the group. However, uniqueness is not provided, au contraire, in general there are infinitely many graphs with such a property ^[7]. For this reason, the definition of the graph structure of a dynamic network precedes the group properties.

4.   Perturbations

A perturbation is a small 'change' of the differential equations governing the system. The order of magnitude can be quantified by a parameter $\epsilon$, $0 < \epsilon \ll 1$. A perturbation is called singular if the solutions of the perturbed ODEs have a different qualitative behavior with respect to the solutions of the unperturbed system ^[60]. On the other hand, regular perturbations are the ones that preserve the qualitative structure. For this reason, singular perturbations are the most interesting if we aim to characterize a system qualitatively. In Appendix B, we briefly review some of the standard terminology of singular perturbations.

4.1.   Perturbed ALFs

As mentioned in the introduction, we now consider perturbations of ALFs given by

where $0 < \epsilon \ll 1$, $H(\mathit{\boldsymbol{x}}, \Lambda) = (h_1(\mathit{\boldsymbol{x}}, \Lambda), \dots, h_n(\mathit{\boldsymbol{x}}, \Lambda))^\intercal$, and $\Lambda$ is a set of parameters. The consensus space $C$ is always a set of equilibrium points for $\epsilon = 0$ (Corollary 1). Moreover, $C$ is one-dimensional and hence generic perturbations are singular ^[60]. The system as written in (4.1) is in non-standard form. In many cases, it is useful to transform a singularly perturbed system into a standard form, one reason among all is that a large amount of results in singular perturbation theory are stated for standard forms.

Proposition 6. Let $l \in \mathcal{V}$, and $\tilde{\mathit{\boldsymbol{x}}} = (\tilde{x}_1, \dots, \tilde{x}_{n-1})^\intercal$ be an $(n-1)$-dimensional vector. Let $\tilde{L}$ be the $(n-1)\times(n-1)$ matrix obtained from the Laplacian $L$ by removing the $l$-th column and row, i.e.,

and let $v^{(l)}$ be the $(n-1)$-dimensional vector defined as follows:

The coordinate transformation

where $k$ is the constant of motion (2.2), transforms system (4.1) into the standard form

where $g_i (\tilde{\mathit{\boldsymbol{x}}}, k, \Lambda) : = h_i \left(\tilde{x}_1, \dots, x_l \left(k, \tilde{\mathit{\boldsymbol{x}}}\right), \dots, \tilde{x}_n, \Lambda \right)$.

Proof. Substituting the coordinate change Eq (4.2) into Eq (4.1) we obtain the set of equations for $\tilde{x}_i$. Then, the equation for $k$ is retrieved by evaluating the scalar product $\langle \mathbf{1}, \dot{\mathit{\boldsymbol{x}}} \rangle$.

Let us notice that the perturbation transforms the constant of motion into a slow variable of the system. Since the constant $k$ is related to the arithmetic mean of the states, $\langle \mathit{\boldsymbol{x}} \rangle$, perturbing the system, in general, means that the states' density is changing, e.g., the concentration of the system is increasing/decreasing. In other words, the perturbation acts as a drift/dissipation in the undergoing diffusive evolution. Notice that if the perturbation $H$ is 'perpendicular' to the consensus space, $\langle \mathbf{1}, H(\mathit{\boldsymbol{x}}, \Lambda) \rangle = 0$, then $\dot k = 0$, which in turn means that the perturbation problem is regular.

As mentioned, we are interested in singular perturbations for two major reasons. First of all, in the setting of our problem where the unperturbed system has a set of equilibria of dimension greater than zero, generic perturbations are singular. The second reason is that singular perturbations are the ones that could lead to a different qualitative behavior of the system. In the following section, we look at the near consensus behavior under a singular perturbation, and we characterize the dynamics for a selected class of systems.

5.   Consensus dynamics

In this section, we further specialize our analysis to complete networks, $K_n$. We will show how this assumption allows a reduction of the equations, and a detailed description of the consensus behavior will be given.

5.1.   The layer problem

Let us consider a perturbed $K_n$-ALF in the standard form (4.3). The layer problem, obtained by setting $\epsilon$ to zero, reads as

The system described by Eq (5.1) is $(n-1)$-dimensional with one parameter, $k$. It is worth noting that although the graph structure in (5.1) seems to be lost, some information is readily available. On the one hand, the term $\sum_{j} \tilde{L}_{ij}f(\tilde{x}_j)$ corresponds to a $K_{(n-1)}$-ALF, i.e., a complete graph with $(n-1)$ nodes. On the other hand, the scalar term $v^{(l)}_i f\left(x_l \left(k, \tilde{\mathit{\boldsymbol{x}}} \right) \right)$, can be interpreted as a higher-order interaction ^[6] naturally arising from the reduction performed.

Furthermore, and for our specific purposes, there is a more relevant property that survives the transformation: The symmetry. In fact, starting with an $\mathfrak{S}_n$-equivariant unperturbed system, after the transformation into the standard form we get a layer problem which is $\mathfrak{S}_{n-1}$-equivariant. This is a first example of robustness of the symmetries. Since we are interested in the consensus dynamics, we consider Eq (5.1) on the invariant space given by $\text{Fix}(\mathfrak{S}_{n-1})$. Let us notice that, in general, the invariant space $\text{Fix}(\mathfrak{S}_{n-1})$ is not globally stable. For example, in a neighborhood of an attracting branch of the consensus space it is stable. Instead, in a neighborhood of a repelling part of the consensus space it is unstable. This means that the restriction to $\text{Fix}(\mathfrak{S}_{n-1})$ provides a robust description of the system only in the basin of attraction of a stable region of $\text{Fix}(\mathfrak{S}_{n-1})$; in the unstable regions the restriction works only if the initial conditions lie exactly in the invariant space and if the perturbation preserves $\text{Fix}(\mathfrak{S}_{n-1})$. Then, restricting Eq (5.1) to $\text{Fix}(\mathfrak{S}_{n-1}) = \{ \tilde{\mathit{\boldsymbol{x}}} \in \mathbb{R}^{n-1} \ | \ \tilde{x}_1 = \dots = \tilde{x}_{n-1} \}$ we obtain

where $x: = \tilde{x}_1 = \dots = \tilde{x}_{n-1}$, and we used the properties of the Laplacian of a complete graph. So, in the invariant space the system reduces to a one dimensional system (with one parameter). The former consensus space, $C = \text{Fix}(\mathfrak{S}_n)$, gives the equilibria of Eq (5.2). This is a consequence of the fact that the layer problem is a coordinate transformation of the unperturbed non-standard form Eq (2.1), for which we know the consensus space is a set of equilibria.

Remark 5. The parametrized consensus states (2.4) given by $x^*_k = k/n$ are equilibria of Eq (5.2).

Since the equilibria of the layer equation constitute the critical manifold for the slow-fast system, then the consensus space is a (part of the) critical manifold of the system.

Proposition 7. The stability of the consensus set (2.4) is determined by the derivative of the response function. In particular, we can decompose the consensus into three components,

where $C^a$ is attracting (stable), $C^r$ is repelling (unstable), and $C^s$ is singular.

Proof. By straightforward computations we obtain the Jacobian of the system evaluated at the consensus, $- \text{d}_xf\left(k/n\right) n$, which implies the proposition.

For $0 < \epsilon\ll1$, thanks to Theorem 2 (Fenichel), Appendix B, compact submanifolds of $C^{a/r}$ perturb to $O(\epsilon)$-close locally invariant slow manifolds of the singularly perturbed system (4.3) (denoted correspondingly by $C^{a/r}_\epsilon$). For the singularities $C^s$ we need further analyses. We are going to consider the case of transcritical singularities which, as we will show, turn out to be generic due to the geometric properties of the problem.

5.2.   The slow-fast dynamics around the consensus space

We start by considering a symmetry condition on the perturbation. While restrictive, we impose it to guarantee the invariance of $\text{Fix}(\mathfrak{S}_{n-1})$ for the perturbed system (4.3).

Lemma 4. Let $H(\mathit{\boldsymbol{x}}, \Lambda)$ be a perturbation such that $h_i = h$, $\forall i \neq l$, and $h_l = \tilde{h}$. Then, $\text{Fix}(\mathfrak{S}_{n-1})$ is an invariant space for the perturbed system (4.3).

Proof. Let us consider the equations in standard form (4.3) with a perturbation satisfying the conditions $h_i = h$, $\forall i \neq l$, and $h_l = \tilde{h}$. So, if we restrict the equations to

we obtain

where $g(x, k, \Lambda) : = h(x, \dots, k - (n-1) x, \dots, x, \Lambda)$, and $\tilde{g}(x, k, \Lambda) : = \tilde{h}(x, \dots, k - (n-1) x, \dots, x, \Lambda)$. Then $\text{Fix}(\mathfrak{S}_{n-1})$ is an invariant space.

Assuming a perturbation satisfying Lemma 4, we are allowed to continue our analysis on the invariant space (5.3), which we will refer to as $(x, k)$-plane. We remark once again that, in a neighborhood of a stable region of $\text{Fix}(\mathfrak{S}_{n-1})$, the reduction is robust; while in unstable regions the reduction is preserved only by fine tuning of initial conditions and perturbation.

5.2.1.   Genericity of transcritical intersections

The analysis to be performed in this section is local. Thus, let us now consider that the response function can be written as

where $a_j$ are real constant coefficients, $N \in \mathbb{N}$, and we assume that at least one of the coefficients $a_j$, $j\neq1$, is different from zero. We know that the consensus space, for an ALF, does not depend on the particular choice of the response function. On the other hand, the remaining equilibrium sets, that together with $C$ give rise to the full critical manifold, depend on the form of the response function. By detailing the local form of $f$, we are able to further describe the local shape of the other branches of the critical manifold crossing the consensus space. We focus our analysis on a neighborhood of the consensus space. The following proposition shows that if there is a singularity on the consensus space, implying that another branch of equilibria intersects with the consensus space, then we can derive the linear approximation of such intersecting branches of equilibria.

Proposition 8. Let the response function $f$ be of the form (5.5) with $N\geq2$. In the $(x, k)$-plane, if a branch of the critical manifold intersects the consensus space, then such a branch has the local form

where $(x_k^s, \ldots, x_k^s)\in C^s$.

Proof. Let us consider the function $\tilde{\phi}(x, y) : = - [f(x)-f(y + x)]$ and the auxiliary system

where $y$ acts as a parameter, with $y$ close to zero. We are going to evaluate the slope of the tangents of the curves of zeros of $\tilde{\phi}(x, y)$ intersecting the line $y = 0$. Considering the response function (5.5), we expand $f(x + y)$ by using the binomial theorem,

where, for brevity, we collect in $\mathcal O(y^ux^v)$ all higher order terms of such order for all combinations $u+v = N$. So, we can rewrite $\tilde{\phi}(x, y)$ as

where we factored out $y$, and set $l: = k-1$. We regularize $\tilde{\phi}(x, y)$ by removing the line $y = 0$ from the set of zeros, i.e. $\tilde{\phi}(x, y) = :y \tilde{\phi}^ \text{reg}(x, y)$, where

We continue by disregarding the higher-order terms as they do not play a role in our forthcoming arguments. The regularized auxiliary system reads

Let us compute the Jacobian, $J$, for the regularized auxiliary system. We get

where

Since we are interested in what happens at the intersection with $y = 0$, we evaluate the Jacobian on that line, namely

The eigenvalues of $J_0$ are determined by the equation $\det(J_0 - \mu\; \text{Id}) = 0$, which in our case reads

and so the eigenvalues are

Let us notice that, in general, for the regularized auxiliary system the Jacobian $J_0$ cannot be interpreted as the linearisation of the system, because $y = 0$ is no more an equilibrium. However, there could be other equilibria on the line $y = 0$ at the intersections with the other curves of equilibria. At such points the eigenvector with zero eigenvalue gives the direction of the tangent of the intersecting curve of zeros. The eigenvector $\tilde{v}^{(\mu_1)} = (\tilde{v}^{(\mu_1)}_x, \tilde{v}^{(\mu_1)}_y)^\intercal$ must satisfy the equation $J_0\tilde{v}^{(\mu_1)} = 0$, that is

Therefore, since for $j\neq1$ at least one $a_j\neq0$, $\tilde{v}^{(\mu_1)}_x = - \frac{1}{2} \tilde{v}^{(\mu_1)}_y$, and so $\tilde{v}^{(\mu_1)} \propto (1, -2)^\intercal$, where $\propto$ is the proportionality symbol. Notice that $\tilde{v}^{(\mu_1)}$ does not depend on $x$, which means that wherever the intersection takes place, the tangent will have the direction of $\tilde{v}^{(\mu_1)}$. At this point, we need to obtain the equivalent vector for the system of our interest (5.2). Note that the equation $y = k-nx$ puts in relation the auxiliary system with the layer problem (5.2). Thanks to this relation we can construct the transformation connecting the tangent vectors of the two systems, in particular we have

and then the inverse transformation reads

Thus, the tangents of the critical manifold at the intersections with the consensus space, in the original system, have direction given by $v_{\mu_1} = \Pi^{-1} \tilde{v}_{\mu_1} \propto (1, n-2)$, which proves that the slope of the tangent in the $(x, k)$-plane is $n-2$. Therefore, the tangent lines have an equation of the form $\kappa(x) = \tilde{\kappa} + (n-2)x$. In order to find $\tilde{\kappa}$, we use the fact that $\kappa(x^s_k) = n x^s_k$, obtaining $\tilde{\kappa} = 2 x^s_k$; the statement is proven. A sketch of the result of this proposition is shown in Figure 1.

Proposition 8 tells us that, up to linear approximation, all other branches of equilibria of the critical manifold crossing the consensus space are parallel to each other and are in fact transverse to the consensus space. Moreover, for $n > 2$ such branches of the critical manifold locally are also transverse to the fast-flow. As a consequence, we have that pitchfork singularities are not possible on the consensus space for the class of ALFs under examination, this statement is also confirmed by checking the conditions for pitchfork singularities ^[34]. By a similar reasoning, fold singularities on the consensus space are also excluded. Finally, Hopf singularities do not appear either by the fact that the layer system has always real eigenvalues. Thus, one can conclude that for $K_n$-ALFs transcritical singularities are generic on the consensus space, in the sense of codimension one bifurcations; while on other branches of the critical manifold it is possible to have other singularities, see some examples in Section 6.

Remark 6. Although in this section we have only considered complete networks, we conjecture that other topologies also present, generically, transcritical singularities along the consensus space. The main argument is that, no matter the network structure, as long as the consensus space exists, it is a one dimensional linear subspace of the phase-space. Singularities of this linear subspace appear as intersections with other branches of the critical manifold induced by the nonlinear response function $F$. Locally, and generically, one would expect that these intersections are transversal, leading to transcritical singularities (possibly in higher dimensions). In particular, recalling the result of Proposition 3, we know that the symmetries induced by the invariance of the response function generates linear branches of the critical manifold intersecting transversely the consensus set. The great advantage of considering complete networks is the reduction to a planar problem. This reduction requires extra work for other network topologies. Nevertheless, our conjecture will be supported by numerical simulations in Section 6.

5.2.2.   Nondegenerate transcritical singularities

We now state the conditions for (nondegenerate) transcritical singularities.

Proposition 9. Let $\mathit{\boldsymbol{x}}^s_k \in C^s$ be a singular consensus point. Let $f$ and $h_i$, $i = 1, \dots, n$, be $\mathcal{C}^m$-functions with $m \geq 3$, and assume that Lemma 4 holds. If

then $\mathit{\boldsymbol{x}}^s_k$ is a nondegenerate transcritical singularity of Eq (5.4).

Proof. We check the conditions for transcritical singularities on planar systems provided by Krupa and Szmolyan ^[34]. Given a planar system of the form (B.1) then the origin is a transcritical singularity if

We now prove that these conditions applied to our system lead to the conditions of the statement. Conditions Eqs (5.6) and (5.7) are satisfied because $x^s_k$ is a singular consensus point. For the non-degeneracy condition (5.8), we obtain

The transversality condition (5.9) gives,

which implies $n \neq 2$. Finally, the last condition (5.10) leads to $\langle \mathbf{1}, H(\mathit{\boldsymbol{x}}^s_k, \Lambda) \rangle \neq 0$.

We consider the case where all the singularities on the consensus space are nondegenerate transcritical, i.e., the conditions of Proposition 9 are met on $C^s$. So, in our current setting, the consensus space consists of a cascade of transcritical singularities changing the stability from attracting to repelling, and vice-versa. If at $x^s_k$ the stability transition, in the direction of the slow flow, is from attracting to repelling, we call $x^s_k$ a type-1 transcritical point. Type-1 transcritical points are characterized by the sign ratio

If at $x^s_k$ the stability transition, in the direction of the slow flow, is from repelling to attracting, we call $x^s_k$ a type-2 transcritical point. Type-2 transcritical points are characterized by the sign ratio $\rho = 1$.

Under perturbations the consensus space, being a component of the critical manifold, extends to a slow manifold as predicted by Fenichel, Theorem 2. More specifically, compact normally hyperbolic submanifolds of the consensus perturb to slow manifolds with the same stability properties. The continuation of such slow manifolds close to a singularity depends on the perturbation. Suppose we follow a trajectory starting close to an attracting section of the consensus space. When such trajectory reaches a neighborhood of a singularity its behavior depends specifically on the perturbation and on the geometry of the aforementioned slow manifolds. In particular, a trajectory reaching a neighborhood of a transcritical point could proceed in three different manners.

Exchange of stability The trajectory exchanges to another branch, which intersects the consensus space transversely, of the critical manifold. In other words, the trajectory does not follow anymore the consensus space, but another stable branch of the critical manifold that results from the nonlinear response function.

Fast escape The trajectory, once traversing a small neighborhood of the singularity, continues along the fast flow.

Bifurcation delay The trajectory continues along an unstable slow manifold that is perturbed from a normally hyperbolic repelling submanifold of the consensus space for a $O(1)$ time.

The third characteristic behavior is called canard, and acts as a separating case between the other two ^[35,60]. Usually, the presence of a canard in a planar system is associated to a critical value of the parameters of the system. The description in terms of the parameters of the system has the advantage of enabling the unfolding of the singularity to be studied. Here, however, we adopt a geometric perspective on the problem. Near fold singularities a geometric perspective on canard solutions in the plane has been treated extensively in ^[13]. In our case of study we show that the requirement for a canard solution can be stated as an algebraic condition on the perturbation.

Proposition 10. Let $\mathit{\boldsymbol{x}}^s_k \in C^s$ be a singular consensus point satisfying the conditions of Proposition 9. If the perturbation satisfies the property

where $\propto$ is the proportionality sign, then for $\epsilon$ sufficiently small, system (4.3) admits canard solutions.

When a perturbation satisfies the condition of Proposition 10, we refer to it as critical perturbation, and we use the notation $H^ \text{crit}(\mathit{\boldsymbol{x}}, \Lambda)$.

Proof. We use again the results from ^[34] to check that our statement holds. Considering a system of the form (B.1), satisfying conditions (5.6)–(5.10), then there exists a parameter $\lambda$ controlling the behavior at the transcritical singularity (which we remind to be at the origin for this benchmark system). The parameter $\lambda$ is given by

where

So, the parameter $\lambda$ for the system (5.4), at a transcritical point $x_k^s$, reads

For type-1 transcritical points, we have $\rho = -1$, and the condition for canards is $\lambda = \lambda^ \text{crit} : = 1$. So, we obtain the equation

Since the case $n = 2$ is excluded by Proposition 9, we have $\tilde{h}(\mathit{\boldsymbol{x}}^s_k, \Lambda) = h(\mathit{\boldsymbol{x}}^s_k, \Lambda)$. Therefore the critical perturbation $H^ \text{crit}(\mathit{\boldsymbol{x}}, \Lambda)$ at $\mathit{\boldsymbol{x}}^s_k$ reads

which agrees with the statement. In order to complete the proof, we show that for a type-2 transcritical point a critical perturbation admits faux-canards. Let us recall that faux-canards for planar systems are generic, and appear for $\lambda < \lambda^ \text{crit}$ ^[34]. Using the fact that for a type-2 transcritical point $\rho = 1$, we get $\lambda = -1$.

In Figure 2, we sketch the blow-up of the transcritical points under a critical perturbation ^[33,34,40]. The blow-up picture makes clear the different behavior of the type-$1/2$ transcritical points and illustrates how the continuation of the consensus space occurs.

Note that we can interpret the result of Proposition 10 as follows: If the perturbation is tangent to the consensus space at a singularity, then canard solutions exist. As a matter of fact, assuming a perturbation that is (everywhere) tangent to the consensus space, it is possible to prove a result that holds for every network structure.

Proposition 11. Given an ALF with arbitrary graph structure $\mathcal{G}$, and $H(\mathit{\boldsymbol{x}}, \Lambda)$ a perturbation satisfying $H(\mathit{\boldsymbol{x}}^*, \Lambda) \propto \mathbf{1}$ and $H(\mathit{\boldsymbol{x}}^*, \Lambda) \neq 0$, $\forall \mathit{\boldsymbol{x}}^* \in C$. Then the consensus space is a trajectory of the perturbed system.

Proof. From Corollary 1 we know that the consensus space is an equilibrium of the unperturbed system; moreover, since $H(\mathit{\boldsymbol{x}}^*, \Lambda) \propto \mathbf{1}$, $\forall \mathit{\boldsymbol{x}}^* \in C$, we have that $C$ is an invariant space for the perturbed vector field. So, projecting Eq (4.1) onto $C$ we obtain

$\mathit{\boldsymbol{x}}^* \in C$. Given the assumption that the perturbation is non-zero on the consensus space, which is a one-dimensional space, we obtain that $C$ is a trajectory of the perturbed system.

In essence, Proposition 11 gives some conditions under which a perturbed equivariant ALF (4.1) exhibits a maximal canard. Let us also notice that Proposition 11 holds for all simple graphs, and it does not make any particular assumption on the singularities of the consensus space. In other words, paying the price of making more stringent assumptions on the perturbation we can cover the cases where the slow-fast theory would be difficult, or not even possible to apply.

6.   Numerical simulations

In this section, we present several examples of ALFs, exploiting generic and canard behavior from different perspectives. The simulations we show are performed in Mathematica ^[63], where the precision of the numerical simulations can be easily set. We show phase portraits, time series and spatio-temporal plots highlighting the nonlinear diffusion and drift mechanism on non-complete graphs. In view of Proposition 2, our simulations are restricted to simple graphs but we emphasize that for weighted graphs one would obtain equivalent results.

6.1.   $K_3$ networks: Generic and canard behavior

We start by looking at a $K_3$-ALF considering two nonlinear response functions, one that is even $f(x) = (x-1)^2(x+1)^2$ (Figure 3), and one odd $f(x) = (x-1)(x+1)^2$ (Figure 4). The two response functions induce different critical manifolds, Figures 3a and 4c. Nevertheless, the typical behavior is similar, see Figures 3b and 4a.

It is worth recalling now a few issues with numerical simulations of transcritical points, see more details in ^[18]. Roughly speaking, canard solutions, as shown in Figure 3c, contract first towards the critical manifold, cross a singularity, and then follow the unstable branch of the critical manifold enough time to compensate the aforementioned contraction before leaving an $\mathcal O(\epsilon)$-neighborhood of the critical manifold. For the particular case of a $K_3$-ALF with response function $f(x) = (x-1)^2(x+1)^2$, one can check that the transverse eigenvalue along the critical manifold is $-\frac{4}{9}k(9-k^2)$ implying that the contraction and expansion rate is symmetric with respect to the origin (this balance is more precisely given by the entry-exit relation ^[18,35]). Such a time symmetry should be visible as a spatial symmetry in correspondence with a transcritical point on the consensus manifold. As one can see in Figure 3c, the symmetry described is not very well reproduced by the numerical simulation, although the computations are performed with a working precision of $50$ digits. There are several factors that induce such a problem. First, we have that the symmetry is exact only in the limit for $\epsilon \to 0$, so we expect that a numerical simulation will just approximate such symmetry. The simulation displayed in Figure 3c is performed with $\epsilon = 1/10$, which is quite large to clearly see this effect. However, by lowering $\epsilon$ one would need to further increase the working precision, leading to problems with the time performance of the CPU. Since our present work has a theoretical focus, we do not dive deeper in the numerical analysis of these delicate systems, leaving the problem of highly accurate simulations for future works.

There is another detail to take into account: In Figure 3c, we have canard solutions that cross three transcritical points on the consensus space without leaving the repelling intermediate section (light blue and pink curves). Here the question is: Where should we expect the canards to leave the consensus manifold? As we just discussed, we expect some symmetry, but here we need to be more careful stating with respect to which transcritical point. We argue that, for these particular solutions, the intermediate section has a null effect, and, therefore, the solution should leave the consensus after the third transcritical point (counting from top to bottom) after a distance (approximately) equal to the distance traveled close to the attracting section before the first transcritical point. In order to prove the aforementioned statement, we use the slow-divergence integral ^[13,14,17]. So, let us consider the slow-divergence integral along the consensus manifold for the system under analysis. Since the perturbation is constant the divergence of the vector field in the $(x, k)$-plane, evaluated on the consensus space, is simply $- n \text{d}_x f (k/n)$. Therefore the slow-divergence integral becomes

We now notice two facts. First, the intermediate section is divided into two spatially equal sections by a transcritical point. Second, on the first subsection $\text{d}_x f (k/n) > 0$, and on the second subsection we have $\text{d}_x f (k/n) < 0$. Therefore the slow-divergence integral on the intermediate section of the consensus space is null, so on the intermediate section the attracting and repelling effect compensate each other. The considerations we have done hold for the dynamics near consensus under a constant critical perturbation. It is worth mentioning that the nonlocal interaction of transcritical singularities via hysteresis processes could give rise to an enhanced delay for the canard ^[21,59].

6.2.   Spatio-temporal configurations

In this section, we consider the spatio-temporal configurations induced by different network structures both under generic and critical perturbations, the latter giving rise to canards. The simulations are not subject to any reduction, and the evolution of the full system is computed. Given a graph we consider a planar embedding, and we assume the obtained coordinates as representative of a portion of a square. The spatio-temporal configurations are then obtained by evaluating a heat map where the colors, from blue to yellow, are in correspondence with the state's values of the system. Each frame represents a unit time step in the sampled interval, and the time flow is represented as follows:

We considered the following network structures: random networks (Figures 5 and 6), a complete graph with $10$ nodes (Figure 7), and a lattice with $30 \times 30$ nodes (Figures 8–10). Throughout this section we adopt the response function $f(x) = (x-1)^2(x+1)^2$.

Spatio-temporal patterns are usually considered in the realm of partial differential equations ^[4,16,44] where a wide range of equations, e.g., reaction-diffusion systems, can give rise to complex patterns observable in nature. Here, the continuous space is replaced by a network and consequently in place of partial differential equations we have a set of ordinary differential equations. Given the analogy with diffusion and drift mechanisms that perturbed ALFs undergo it is quite natural to explore the potential spatio-temporal patterns appearing under different conditions. The most remarkable observation we can draw from the simulations is that under generic perturbations no recognizable pattern appear, while under a critical perturbation we observe transient patterns. In other words, our simulations suggest that spatio-temporal patterns are induced by canards. Such an effect is related to the delayed loss of stability produced by a canard. In fact, we generically have that the individual states, once they reach the singularity, leave consensus one after the other, with no particular order. On the other hand, under a critical perturbation a significant amount of states transition at the same time and such a configuration persists for several units of time. Such phenomenon is more clearly visible when the number of nodes is large, to see this compare the plots for a complete graph with $10$ nodes (Figure 7c), a random graph with $30$ nodes (Figure 6a) and a lattice with $900$ nodes (Figure 10).

It is worth mentioning that spatio-temporal patterns are commonly studied and understood for reaction-diffusion equations where the two species interaction could give rise to patterns. Here, we are studying equations that resemble a drift-diffusion, and the (transient) pattern structure arises from the nonlinearity of the diffusion mechanism and the geometry of the drift.

In general, we can notice that the nonlinearity given by the response function induces a different diffusion mechanism. By looking at the generic cases, illustrated in Figures 5a, 7a, and 9, we see that the transition between the initial consensus state to the final one exhibits jumps instead of a homogeneous process. Such an effect is due to the fact that the nonlinearity induces an intricate structure of the critical manifold, as we have already seen in the previous examples. We can partially observe the homogeneous process of diffusion when the canard conditions apply, see Figures 6a and 7c (for the lattice in Figure 10 such an effect is suppressed by the image quality compression).

As a final remark, we observe that the transition time between two consensus sections is visibly shorter in the canard cases.

7.   Conclusions and discussion

This paper has been dedicated to the analysis of a class of dynamical systems describing nonlinear diffusion on a simple graph structure. The analysis starts by looking at (unperturbed) ALFs on graphs with positive weights. We emphasize the importance of the fiber-wise topological equivalence proved in Proposition 2. Indeed, a weighted graph, in general, loses the symmetry properties that its own unweighted counterpart has. The fact that a positive-weighted graph induces an evolution of the system that is topologically equivalent to the one induced by the unweighted counterpart, allows us to consider unweighted simple graphs. An ALF with a simple graph structure inherits the symmetry properties of the underlying graph, enabling us to use the tools from equivariant dynamical systems theory, such as fixed-point spaces and invariant properties of the response function. The structure of the equations allows for a detailed analysis of the equilibria, that in turn led us to study the local structure of singularities on the consensus space for complete graphs (Propositions 8 and 10), and to prove the existence of a maximal canard for generic simple graphs (Proposition 11).

Several new directions arise from the results of our paper. For example, one could consider response functions $f(\mathit{\boldsymbol{x}})$ depending also on "neighboring" nodes; in such a case the network structure also appears inside the response vector field, and not only through the Laplacian matrix. Moreover, one may be interested in studying the case of heterogeneous response functions, i.e., different nodes might contribute differently to the mutual interaction. Regarding the graph structure, one may take into account negative weights or directed graphs, which, even for the linear case, give rise to much richer dynamics than the positive-weights case. Hypergraphs are also becoming more and more relevant. Incorporating a hypergraph structure in our model could take different directions. The simplest choice would be to study which class of hypergraph models can be associated with an ALF. Of course, such a choice would lead unavoidably to a restricted family of hypergraphs. On the other hand, one could extend the model to hypergraphs, by replacing the graph Laplacian with the hypergraph Laplacian, and the other quantities accordingly. In this case, the system will differ dramatically from the one we studied. However, our analysis could serve as a guideline: Exploiting symmetries and algebraic properties could lead to relevant simplifications for the study. In a similar line of thought, reduction techniques similar to those presented here, but for other network topologies than complete networks, seems worth pursuing. Another possibility is to study Laplacian systems where the states have dimensions larger than one, for example, one might consider Laplacian systems where the states evolve in the real plane, or even in a two-dimensional manifold. The dynamics of planar systems entail interesting phenomena, such as periodic orbits, that could potentially give rise to new phenomena unseen in the one-dimensional framework. For similar reasons, states of dimension three might be relevant as well, since it is well-known that it is for three-dimensional systems that chaos, at the level of the nodes, may appear. In all the aforementioned cases, a natural question is to investigate, and classify the possible singularities that arise and the influence of the network on them (for example, in this paper, we showed that complete ALFs can only have transcritical singularities along the consensus manifold). Moreover, one would then be interested in the behavior under small perturbations, which, as seen here, is also closely related to the network structure. Finally, given its relevance for network dynamics, it would be interesting to further investigate which perturbations lead to persistence, or on the other hand destroy, the singular canards of Proposition 10. One possibility is already given in Proposition 11, but its generalization seems rather promising.

Author contributions

The principal author is RB, who conducted the research and wrote the article. HJK contributed to the research phase through active discussion and undertook some revisions during the writing process.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors thank the anonymous reviewers whose thorough comments helped to improve the manuscript. All datasets are available upon request to the corresponding author.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A.   Groups

We recall some basic notions about group theory ^[23,49].

Definition 5. Let $\Gamma$ be a set and $\star: \Gamma \times \Gamma \to \Gamma$ a binary operation. A group is a couple $(\Gamma, \star)$ satisfying the following properties:

1. Associativity: $(\gamma_i \star \gamma_j) \star \gamma_k = \gamma_i \star (\gamma_j \star \gamma_k)$, $\forall \gamma_i, \gamma_j, \gamma_k \in \Gamma$.

2. Identity: $\exists \ e \in \Gamma$ such that $e \star \gamma = \gamma \star e = \gamma$, $\forall \gamma \in \Gamma$.

3. $Inverse:$ $\forall \gamma \in \Gamma$, $\exists \ \gamma^{-1} \in \Gamma$ such that $\gamma^{-1} \star \gamma = \gamma \star \gamma^{-1} = e$.

The operation $\star$ is called group operation. Note that, closure under the group operation follows from the definition of the operation itself. Unless strictly necessary, we refer to a group just by the set, e.g., $\Gamma$, implicitly considering it equipped with a group operation; we also omit the group operation symbol, e.g., $\gamma_i \gamma_j = \gamma_i \star \gamma_j$.

Definition 6. A group is said to be finite if the set $\Gamma$ has a finite number of elements. The order of a finite group, $\text{Ord}(\Gamma)$, is the number of elements in the group.

Given a vector space, a representation, $\psi$, is a homomorphism from the group to the general linear group of the vector space. Since the states, $\mathit{\boldsymbol{x}}$, of an ALF belong to $\mathbb{R}^n$, we are interested in how $\Gamma$ acts on $\mathbb{R}^n$, therefore, we consider representations $\psi : \Gamma \to GL(\mathbb{R}^n)$. To avoid new unnecessary notation, we omit to explicitly write the representation. It is implicit that when we write a group element applied to a vector we are considering the representation on the appropriate vector space, e.g., $\gamma \mathit{\boldsymbol{x}} = \psi(\gamma) \mathit{\boldsymbol{x}}$.

Definition 7. Given two groups $\Gamma^{(1)}$, $\Gamma^{(2)}$, the direct product $\Gamma^{(1)} \times \Gamma^{(2)}$ is a group with set of elements given by the Cartesian product

and operation defined as follows

$\forall \gamma_{i}^{(1)}, \gamma_{k}^{(1)} \in \Gamma^{(1)}$ and $\forall \gamma_{j}^{(2)}, \gamma_{l}^{(2)} \in \Gamma^{(2)}$.

It follows that the direct product can be extended to an arbitrary number of groups, $\Gamma^{(1)} \times \Gamma^{(2)} \times \dots \times \Gamma^{(n)}$. Another important concept is the one of fixed-point space, i.e., the set of points that are kept fixed by the action of the group.

Definition 8. Let $\Sigma \subseteq \Gamma$ be a subgroup of $\Gamma$. The fixed-point space (under the action in $\mathbb{R}^n$) of the group $\Sigma$ is denoted by $\text{Fix}(\Sigma)$ and it is given by

A direct interplay between graph theory and group theory is given by the graph's automorphism group. For graphs the automorphism group is a subgroup of the (finite dimensional) symmetric group $\mathfrak{S}_n$ ^[51], i.e., the group of permutations of $n$ symbols.

Definition 9. Given a (simple) graph $\mathcal{G}$, the automorphism group, $\text{Aut}(\mathcal{G})$, is the group of permutations which preserves the adjacency structure of the vertices of $\mathcal{G}$.

Remark 7. From an algebraic point of view, preservation of adjacency structure correspond to the fact that $\sigma \in \text{Aut}(\mathcal{G})$ if and only if $[A, \sigma] = 0$, where $[\cdot, \cdot]$ represent the commutator of matrices ^[7]. Notice that the same statement holds also if we substitute the adjacency matrix with the Laplacian matrix.

A.1.   Invariance and equivariance

Definition 10. Let $\Gamma$ be a group^*. A function $f$ is called $\Gamma$-invariant if

*Notice that this is not related to $\text{Aut}(\mathcal{G})$, we use the same symbol $\Gamma$ for different groups in different context.

where $\circ$ denotes the composition of maps.

The concept of symmetry, in dynamical systems theory, refers to the fact that there exists a group mapping solutions to solutions. When a vector field is given, such a concept can be translated in terms of equivariance of the vector field.

Definition 11. Let $\Gamma$ be a group, and $X$ a vector field. Then $X$ is called $\Gamma$-equivariant if

An element of the group, $\gamma \in \Gamma$, is called a symmetry of the vector field.

Remark 8. It is straightforward to see that if a system is $\Gamma$-equivariant, then $\Gamma$ maps solutions to solutions.

A solution of a dynamical system is a simple example of an invariant set under the flow. Invariant sets, in general, are extremely relevant in the study of ODEs, since they allow to reduce the analysis to potentially simpler equations. We recall here the general definition of invariant set under the flow induced by a vector field. Later we state the interplay between symmetries and invariant sets.

Definition 12. Let $\phi$ be the flow generated by a vector field on $\mathbb{R}^n$. A set $\mathcal{M} \subset \mathbb{R}^n$ is said to be invariant (respectively positively/negatively invariant) under the flow if for every $x\in\mathcal M$, $\phi_t(x)\in\mathcal M$, $\forall t \in \mathbb{R}$ (respectively $\forall t \in \mathbb{R}_{+} / \mathbb{R}_{-}$).

In order to keep a compact notation, we call a set invariant under the flow a $\phi$-invariant set. From equivariance of the vector field we can immediately infer the existence of a $\phi$-invariant set.

Proposition 12. The fixed-point space, $\text{Fix}(\Gamma)$, is a $\phi$-invariant set for $\Gamma$-equivariant systems ^[26].

Appendix B.   Singular perturbations

In general, the symbols used in this section are not related to those in the main text.

A singularly perturbed system of ODEs, in the standard form, is also called a slow-fast system, and it has the form

where $0 < \epsilon \ll 1$, $x \in \mathbb{R}^n$, $y \in \mathbb{R}^m$, and the time-parametrisation $t$ is called fast-time. Eq (B.1) is said to be in the standard form because the time separation between the fast-variables, $x$, and the slow-variables, $y$, is explicit. In the non-standard form a singularly perturbed system has the form

where $z \in \mathbb{R}^{n+m}$, and the vector field $F(z, 0)$ has at least one set of equilibria of dimension greater than zero ^[60].

Defining the slow-time, $\tau: = \epsilon t$, Eq (B.1) can be rewritten as

Setting $\epsilon$ to zero in Eq (B.2), we obtain the reduced system

which is a set of constrained differential equations. Instead, setting $\epsilon$ to zero in Eq (B.1) we obtain the layer problem

where the slow-variables, $y$, play the role of parameters.

Definition 13. The critical manifold is defined as the set

Remark 9. Regardless of the commonly accepted name, the critical manifold is not necessarily a manifold.

The critical manifold is the constraint space for the reduced system (B.3), or equivalently, it is the set of equilibria of the layer problem (B.4).

Definition 14. A point $x^* \in \mathcal{C}_0$ is called normally hyperbolic if the Jacobian matrix $\text{D}_x f(x^*, y, 0)$, associated to the layer problem (B.4), has no eigenvalues on the imaginary axis.

Definition 15. A submanifold of the critical manifold $\mathcal{M}_0 \subseteq \mathcal{C}_0$ is called normally hyperbolic if all points in $\mathcal{M}_0$ are normally hyperbolic. Moreover, invariant manifolds that are normally hyperbolic are called normally hyperbolic invariant manifold (NHIM) ^[30].

Theorem 2 (Fenichel). For $\epsilon$ sufficiently small, compact normally hyperbolic submanifolds of the critical manifold, $\mathcal{M}_0 \subseteq \mathcal{C}_0$, persists as $O(\epsilon)$-close locally invariant slow manifolds $\mathcal{M}_\epsilon$ diffeomorphic to $\mathcal{M}_0$. Moreover, $\mathcal{M}_\epsilon$ is normally hyperbolic and has the same stability properties with respect to the fast variables as $\mathcal{M}_0$.

Remark 10. Theorem 2 is restricted to the results relevant for the present paper. An exhaustive version can be found in ^[20,35]. The terminology $O(\epsilon)$-close means that, as $\epsilon \to 0$, $\mathcal{M}_\epsilon$ has Hausdorff distance $O(\epsilon)$ from $\mathcal{M}_0$.

For nonhyperbolic points, $x^* \in \mathcal{C}_0$, also known as singularities, Fenichel's theorem does not provide information about the perturbed system in a neighborhood of $x^*$. The process of analysis of the system in a neighborhood of a singularity is called desingularisation. We briefly introduce the blow-up (desingularisation) for nilpotent singularities ^[33].

Definition 16. A point $x^* \in \mathcal{C}_0$ is a nilpotent singularity if all the eigenvalues of the Jacobian, at $x^*$, are zero.

When a singular point has a mixture of zero eigenvalues, and non-zero (real part) eigenvalues it is possible to reduce the problem to a nilpotent singularity via a center manifold reduction ^[10].

Definition 17. Let $x^*$ be a nilpotent singularity. The spherical blow-up is a local transformation

where $\mathbb{S}^{n+m}_{x^*}$ is the $n+m$-dimensional sphere, embedded in $\mathbb{R}^{n+m+1}$, centred at $x^*$, $(\bar{x}, \bar{y}, \bar{\epsilon}) \in \mathbb{S}^{n+m}_{x^*}$, $I \subseteq \mathbb{R}_{\geq 0}$ is an interval such that $0 \in I$, and $r \in I$.

Let us notice that, for $r > 0$ the transformation (B.5) is a diffeomorphism, while for $r = 0$ it is not invertible. Roughly speaking, the blow-up maps the point $x^*$ into a $(n+m)$-dimensional sphere. The core idea of the blow-up is to analyze the induced flow on the sphere and deduct the behavior of the singularity. However, the blow-up transformation by itself does not desingularise the system. After the blow-up it is necessary to perform a conformal transformation (time reparametrization) in order to regularize the system, i.e., retrieve hyperbolicity. We described the spherical blow-up, but several other blow-up transformations are possible, e.g., quasihomogeneous, directional, or even to other manifolds, not necessarily spheres ^[36]. The choice of the appropriate blow-up should be done by considering the regularization, and the feasibility of the analysis. Since this work is not concerned with a direct application of the blow-up, but rather with the implementation of results involving the blow-up technique, we presented the spherical set-up which provides the most clear geometrical visualization.