Mohamed Maghenem
email: mohamed.maghenem@cnrs.fr

Elena Panteley
email: elena.panteley@cnrs.fr

Antonio Loría
email: antonio.loria@cnrs.fr

Singular-Perturbations-Based Analysis of Dynamic Consensus in Directed Networks of Heterogeneous Nonlinear Systems

Keywords: Consensus, multi-agent systems, singular perturbations, network systems, synchronization Mathematics Subject Classification-34D15, 34D05, 34C25, 93C10, 34C28

We analyze networked heterogeneous nonlinear systems, with diffusive coupling and interconnected over a generic static directed graph. Due to the network's hetereogeneity, complete synchronization is impossible, in general, but an emergent dynamics arises. This may be characterized by two dynamical systems evolving in two time-scales. The first, "slow", corresponds to the dynamics of the network on the synchronization manifold. The second, "fast", corresponds to that of the synchronization errors. We present a framework to analyse the emergent dynamics based on the behavior of the slow dynamics. Firstly, we give conditions under which if the slow dynamics admits a globally asymptotically stable equilibrium, so does the networked systems. Secondly, we give conditions under which, if the slow dynamics admits an asymptotically stable orbit and a single unstable equilibrium point, there exists a unique periodic orbit that is almost-globally asymptotically stable. The emergent behavior is thus clear, the systems asymptotically synchronize in frequency and, in the limit, as the coupling strength grows, the emergent dynamics approaches that of the slow system. Our analysis is established using singularperturbations theory. In that regard, we contribute with original statements on stability of disconnected invariant sets and limit cycles.

Introduction and motivation

Networks of nonlinear heterogeneous systems are both, ubiquitous and complex. Their ubiquity motivates their study across numerous research disciplines, as varied as Engineering Systems theory [START_REF] Chow | Time-Scale Modeling of Dynamic Networks with Applications to Power Systems[END_REF], Complexity theory [START_REF] Heylighen | Self-organization, emergence and the architecture of complexity[END_REF] or, even, Philosophy of Science [START_REF] Silberstein | Reduction, emergence and explanation[END_REF]. Their complexity is motor for two apparently antagonistic trains of thought that attempt to explain the collective behavior of networked systems in a broad sense: reductionism and emergentism. The first asserts that any whole can be reduced to its constituent parts-as in the case of networked linear systems [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF], while the tenet of emergentism is that a new behavior appears as a consequence of the interaction of the said parts [START_REF] Silberstein | Reduction, emergence and explanation[END_REF]-as in networks of heterogeneous nonlinear systems. What is more, one of the accepted definitions of Complexity is that it corresponds to the difference between the network as a whole and the sum of its parts and, in that regard, nonlinearity is a necessary condition for complexity to appear [START_REF] Arecchi | A Critical Approach to Complexity and Self Organization, ch. in Mathematical Undecidability[END_REF].

In this paper, we show that, to some extent, both emergentism and reductionism are not necessarily mutually exclusive, but their respective underlying postulates are useful to assess the behavior of networks of heterogeneous nonlinear systems. We focus on systems with dynamics given by ẋi = f i (x i ) + u i , i ∈ {1, 2, . . . , N },

x i ∈ R n , (1)
where i ∈ {1, 2, . . . , N }, x i ∈ R n is the state of the ith system and u i ∈ R n is the decentralized control input to each system, defined as the consensus control law

u i := -σ l i1 (x i -x 1 ) + . . . + l iN (x i -x N ) , (2)
where l ij are different non-negative real numbers denoting the individual interconnection weights and the scalar parameter σ > 0 is the common coupling strength. The control law [START_REF] Heylighen | Self-organization, emergence and the architecture of complexity[END_REF] is reminiscent of that commonly used in the literature on consensus control, in which, the coupling strength σ = 1. This is specifically the case for networks of linear systems, in which case complexity hardly appears and the focus turns towards relaxing the various conditions pertaining to the nature of the interconnections. These may be linear [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF]; nonlinear [START_REF] Arcak | Passivity as a design tool for group coordination[END_REF][START_REF] Isidori | Robust output synchronization of a network of heterogeneous nonlinear agents via nonlinear regulation theory[END_REF]; time-varying [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF][START_REF] Chowdhury | Persistence-based convergence rate analysis of consensus protocols for dynamic graph networks[END_REF]; switching [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF][START_REF] Adhikari | An emerging dynamics approach for synchronization of linear heterogeneous agents interconnected over switching topologies[END_REF] or even state dependent [START_REF] Teel | Smooth time-varying stabilization of driftless systems over communication channels[END_REF][START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF]. Other topology aspects, such as whether the graph is directed [START_REF] Martin | Time scale modeling for consensus in sparse directed networks with time-varying topologies[END_REF] or the interconnections are signed [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF], may also alter consensus. For networks of heterogeneous nonlinear systems, however, the coupling strength plays a central role. Different kinds of emergent behavior may arise depending on whether σ is "weak" [START_REF] Pogromsky | Pattern generation in diffusive networks: How do those brainless centipedes walk?[END_REF][START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] or "strong" [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF].

Our main interest in this paper is to assess the behavior of the corresponding closed-loop system, specifically, in the case that the coupling gain σ is larger than a certain threshold, but we restrict our analysis to networks with an underlying static directed graph. Akin to [START_REF] Panteley | A stability-theory perspective to synchronisation of heterogeneous networks[END_REF] and [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], we analyze the closed-loop networked system via a change of coordinates-introduced in [START_REF] Maghenem | Singular-perturbations-based analysis of synchronization in heterogeneous networks: a case-study[END_REF]-that exhibits an intrinsic dichotomous structure composed of two dynamics defined in orthogonal spaces. On one hand, one has a reduced-order dynamics with state x m ∈ R n (miscalled emergent dynamics in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]) and, on the other, the dynamics of synchronization errors, denoted by e i := x i -x m . This characterization of the networked system is driven by the objective of characterizing synchronization phenomena that may appear (or not) in view of the systems' interconnections.

In [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], we say that the systems reach dynamic consensus if, for all i ≤ N , the synchronization errors e i converge to zero asymptotically. The dynamic consensus paradigm generalizes the more common equilibrium consensus, in which case, the reduced-order dynamics is null, i.e. ẋm = 0, because the collective behavior is static and the state of the reduced-order dynamics satisfies x m (t) ≡ x m (0), where x m (0) is a weighted average of the nodes' states' initial values. In the case of a network of oscillators, the reduced-order dynamics may admit an asymptotically stable equilibrium (an example is provided in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]) or an asymptotically stable attractor [START_REF] Panteley | Practical dynamic consensus of Stuart-Landau oscillators over heterogeneous networks[END_REF]. However, in general, asymptotic dynamic consensus is unreachable due to the heterogeneity [START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF]. An exception is that of systems that admit an internal model [START_REF] Wieland | An internal model principle is necessary and sufficient for linear output synchronization[END_REF][START_REF] De Persis | On the internal model principle in the coordination of nonlinear systems[END_REF][START_REF] Kim | Output consensus of heterogeneous uncertain linear multi-agent systems[END_REF]. Otherwise, for general nonlinear heterogeneous systems, dynamic consensus may be guaranteed only in a practical sense [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF].

In this paper, we establish two statements, each addressing one of two possible cases: one in which the (slow) reduced-order dynamics admits a globally asymptotically stable equilibrium and another one, in which, it admits a periodic solution. In the first case, we give sufficient conditions, under which, the origin for the networked system is globally asymptotically stable. In the second case, we prove that the synchronization errors converge to a unique attractive periodic orbit, so the systems synchronize in frequency. Moreover, for "large" values of the coupling strength σ, this orbit is "close" to that generated by periodic solutions of the reduced dynamics. Thus, the emergent dynamics approaches that of the reduced-order system, as the coupling gain grows.

The analysis is based on the recognized premise that in self-organized complex systems, emergence is multi-level, as occurring in multiple timescales [START_REF] Heylighen | Self-organization, emergence and the architecture of complexity[END_REF][START_REF] Haken | Synergetics[END_REF]. Here, we only consider two, one that is slow and pertains to the reduced-order system and another that is fast and pertains to the synchronization errors. We rely on singular perturbation theory, on statements found in [START_REF] Anosov | Limit cycles of systems of differential equations with small parameters in the highest derivatives[END_REF][START_REF] Kokotović | Singular perturbations and order reduction in control theory-an overview[END_REF][START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF][START_REF] Khalil | Nonlinear systems[END_REF], but also on original refinements of some statements from [START_REF] Anosov | Limit cycles of systems of differential equations with small parameters in the highest derivatives[END_REF] for systems admitting disconnected sets composed of equilibria and periodic orbits.

The model-reduction-and-multi-time-scale perspective is certainly not new, neither in systems theory [START_REF] Chow | Time-Scale Modeling of Dynamic Networks with Applications to Power Systems[END_REF] nor in other disciplines. In the seminal work [START_REF] Chow | Time scale modeling of sparse dynamic networks[END_REF], which follows up on [START_REF] Chow | Time-Scale Modeling of Dynamic Networks with Applications to Power Systems[END_REF], the authors consider a modular network composed of sparsely connected clusters of densely interconnected dynamical systems modeled by simple integrators-the paradigm is motivated by that of large electrical networks. Using classical singularperturbation theory [START_REF] Anosov | Limit cycles of systems of differential equations with small parameters in the highest derivatives[END_REF][START_REF] Kokotović | Singular perturbations and order reduction in control theory-an overview[END_REF], it is showed that such networks achieve synchronization at two levels, within and among the clusters. The analysis is based on relating the network's sparsity to a singular-perturbation parameter. These concepts have been revisited in many succeeding works, such as [START_REF] Bıyık | Area aggregation and time-scale modeling for sparse nonlinear networks[END_REF] and [START_REF] Martin | Time scale modeling for consensus in sparse directed networks with time-varying topologies[END_REF]. In the former, for networks of simple integrators through sector nonlinearities, and in the latter, for linear homogeneous systems interconnected through time-varying persistently-exciting gains a la [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF]. On the other hand, in [START_REF] Tognetti | Synchronization via output feedback for multi-agent singularly perturbed systems with guaranteed cost[END_REF][START_REF] Rejeb | Control design with guaranteed cost for synchronization in networks of linear singularly perturbed systems[END_REF], networks of linear homogeneous singularly-perturbed systems are considered. Thus, in all of the above, the setting is fundamentally different from the one adopted here.

In [START_REF] Maghenem | Singular-perturbations-based analysis of synchronization in heterogeneous networks: a case-study[END_REF], for a particular case-study of networked Andronov-Hopf oscillators, we use a coordinate transformation to exhibit the presence of the two-timescale emergent dynamics and singular-perturbation theory to analyze the collective behavior under the premise that the reduced-order system admits an asymptotically stable orbit. Based on the coordinate transformation introduced in [START_REF] Maghenem | Singular-perturbations-based analysis of synchronization in heterogeneous networks: a case-study[END_REF], singular-perturbation theory is used in [START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF] on a wider class of nonlinear systems with rank-deficient coupling to establish synchronization in the practical sense. Thus, there are several articles in the literature that explicitly use reduction and singular-perturbation theory, even in a multi-agent context. Yet, we are not aware of any such work whose scope covers generic nonlinear heterogeneous systems interconnected over directed graphs and characterize the collective behavior with higher precision. Conceptually, in phase with the emergentism posit, we exhibit the emergence of a complex (in the sense of [START_REF] Arecchi | A Critical Approach to Complexity and Self Organization, ch. in Mathematical Undecidability[END_REF]) dynamic behavior, as a result of the systems' interactions. At the same time, we give a more precise characterization (well beyond practical asymptotic stability of the synchronization manifold) of the collective behavior of networked systems based on that of a reduced-order model.

The rest of the paper is organized as follows. In Section 2, we exhibit the network's reduced-order and synchronization dynamics, under an invertible coordinate transformation. In Sections 3 and 4, we present our main results for the two cases described above, respectively. In Section 5, we revisit the case-study of Stuart-Landau oscillators, and we conclude with some remarks in Section 6. The paper is completed with technical appendices.

Notation. Given a nonempty set K ⊂ R n , |x| K := inf y∈K |x -y|, where |s| denotes the Euclidean norm of s, defines the distance between x and the set K. For a nonempty set O ⊂ R n , K\O denotes the subset of elements of K that are not in O. For a matrix A ∈ R n×n , A -1 denotes its inverse, A denotes its transpose, and |A| denotes its norm. For a matrix Γ ∈ R n×n , λ min (Γ) and λ max (Γ) denote the smallest and the largest eigenvalues of Γ, respectively. By 1 N ∈ R N , we denote the vector whose entries are equal to 1. For a sequence

{A i } N i=1 ⊂ Π N i=1 R ni×ni , blkdiag i∈{1,2,...,N }
{A i } is the block-diagonal matrix whose i-th diagonal block corresponds to the matrix A i . By ⊗, we denote the Kronecker product. For a complex number λ ∈ C, e(λ) denotes the real part of λ and (λ) denotes the imaginary part of λ.

2 On strongly-coupled connected networks

The model and standing assumptions

Consider a group of N nonlinear systems as in (1) driven by the distributed control inputs in [START_REF] Heylighen | Self-organization, emergence and the architecture of complexity[END_REF], where each l ij ≥ 0 is constant but not necessarily equal to l ji . In particular, when l ij is strictly positive, then there exists an interconnection from the jth node to the ith node, but l ji may be null, in which case, the interconnection is unidirectional. More precisely, we pose the following hypothesis.

Assumption 1 (connected di-graph) The network's graph is connected.

Remark 1 We stress that under Assumption 1, any kind of directed graph containing a rooted spanningtree (with or without cycles) is considered. Moreover, under Assumption 1 the Laplacian L has exactly one eigenvalue (say, λ 1 ) that equals zero, while the others have positive real part, i.e., 0 = λ 1 < e {λ 2 } ≤ . . . ≤ e {λ n }. Furthermore, the right eigenvector corresponding to the simple eigenvalue λ 1 = 0, is v r = 1 N ∈ R N , while the left eigenvector, denoted v l , contains only non-negative elements [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF] and satisfies 1 N v l = 1.

•

In addition, for each unit, we impose the following regularity and structural hypotheses:

Assumption 2 (equilibrium) The functions f i are continuously differentiable and, without loss of generality, we assume that f i (0) = 0 for all i ∈ {1, 2, . . . , N }.

Assumption 3 (Semi passive units) Each agent is input-to-state strictly semi-passive; namely, for each i ∈ {1, 2, . . . , N }, there exists a continuously differentiable and radially unbounded storage function V i : R n → R ≥0 , a positive constant ρ, a continuous function H i : R n → R, and a continuous function

ψ i : R n → R ≥0 such that Vi (x i ) ≤ x T i u i -H i (x i )
and

H i (x) ≥ ψ i (|x|) for all |x| ≥ ρ.
Assumption 3 is useful to assess the boundedness of solutions for system (1) in closed loop with (2) for a sufficiently large coupling strength σ. More precisely, we have the following result.

Lemma 1 (Global ultimate boundedness [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]) Consider the systems in (1) in closed loop with the control inputs in (2) and let Assumptions 1-3 be satisfied. Then, the closed-loop system is (uniformly in σ) ultimately bounded; namely, there exist σ * > 0 and r > 0 such that, for any R ≥ 0, there exists τ R ≥ 0 such that, for each solution x(t) starting from x o ∈ R nN , we have

|x o | ≤ R =⇒ |x(t)| ≤ r ∀t ≥ τ R , ∀σ ≥ σ * ,
where x ∈ R nN denotes the network's state, i.e., x = [x 1 , . . . , x N ] .

Under the assumptions listed above, we investigate the problem of assessing the behavior of the networked closed-loop system (1)- [START_REF] Heylighen | Self-organization, emergence and the architecture of complexity[END_REF]. To this end, as it is customary, let us collect the individual interconnection coefficients l ij into the Laplacian matrix

L := [ ij ] ∈ R N ×N , where ij = k∈Ni a ik i = j -a ij i = j.
Then, replacing (2) in (1), we see that the overall network dynamics takes the form

ẋ = F (x) -σ[L ⊗ I n ]x, (3)
where the function F : R nN → R nN is given by

F (x) := [f 1 (x 1 ) • • • f N (x N )] .
As in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] and [START_REF] Panteley | A stability-theory perspective to synchronisation of heterogeneous networks[END_REF], to analyze the behavior of the network system (3), we acknowledge its dichotomous nature. In these references as well as in many others-see, e.g., [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF][START_REF] Liu | Output synchronization of dynamical networks with incrementallydissipative nodes and switching topology[END_REF][START_REF] Delellis | On quad, lipschitz, and contracting vector fields for consensus and synchronization of networks[END_REF], synchronization is defined as the property of the trajectories of each individual system following the trajectories of an "averaged" unit with state

x m := [v l ⊗ I n ]x. ( 4
)
The quotes in "averaged" are superfluous in the case of undirected networks, in which case v l = 1 N , so

x m = 1 N N i=1
x i , but for directed connected networks the state x m is more generically defined as a weighted average of the respective systems' states since v li ≥ 0 for all i ∈ {1, 2, ..., N } and v l 1 N = 1.

In either case, a sensible way to define the synchronization errors e is as the difference between the units' states and x m , that is,

e := x -[1 N ⊗ I n ]x m . (5)
Thus, in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] and [START_REF] Panteley | A stability-theory perspective to synchronisation of heterogeneous networks[END_REF] the collective behavior of network systems is studied in function of the dynamics of the "averaged" unit x m and that of the synchronization errors e.

In the next section, we introduce another change of coordinates to rewrite system (3) in an equivalent form that exhibits two motions; one that is generated by the averaged dynamics and another by a projection of the synchronization errors e on a certain subspace. This coordinate transformation is not a simple artifice for analysis, it exhibits two time-scales that are inherent to networked systems satisfying Assumption 1 and subject to a sufficiently large coupling σ.

Intrinsic two-time-scales decomposition

After Assumption 1 and Remark 1, because λ 1 = 0 has multiplicity one, the Laplacian admits the following Jordan-block decomposition of over R N ×N :

L = U 0 0 0 Λ U -1 , (6)
where Λ ∈ R (N -1)×(N -1) is composed by the Jordan blocks corresponding to the N -1 non-zero eigenvalues.

Remark 2 Note that even though a Jordan decomposition does not necessarily exist with a real matrix U , it is always possible to use the spectral-invariant-subspace decomposition as in Lemma 13 to generate a real matrix U .

•

The convertible matrix U is constituted, column-wise, of the right eigenvector of the Laplacian, 1 N , and a left-invertible matrix V ∈ R N ×(N -1) , which consists of the eigenvectors corresponding to the nonzero eigenvalues of L. That is,

U = [1 N V ], U -1 = v l V † , (7)
where V † ∈ R (N -1)×N , and

v l V = 0, V † V = I N -1 . (8)
So, using ( 7) and ( 8), we also have the useful identity

V V † = I N -1 N v l .
Now, using U -1 , we define the new coordinates

x := [U -1 ⊗ I n ]x (9)
and their inverse transformation

x := [U ⊗ I n ]x. ( 10
)
The interest of the coordinate x is that it consists in the familiar "averaged" states x m and a projection of the synchronization errors e defined in (5) onto the subspace that is generated by V † , which is orthogonal to the right eigenvector 1 N . To better see this, note that such projection yields

[V † ⊗ I n ]e = [V † ⊗ I n ] x -[1 N ⊗ I n ]x m
In the sequel, we refer to the left-hand side of the latter equation as the projected synchronization errors,

e v := [V † ⊗ I n ]x. (11)
Hence, in view of ( 4), ( 7), [START_REF] Chowdhury | Persistence-based convergence rate analysis of consensus protocols for dynamic graph networks[END_REF], and (11), we have

x = x m e v = [v l ⊗ I n ]x [V † ⊗ I n ]x . ( 12
)
In the new coordinates, the network system (3) is equivalently written as

ẋ = [U -1 ⊗ I n ] F (x) -σ[L ⊗ I n ]x ,
which consists in two interconnected dynamics, that of the "averaged" state x m and that of the projected synchronization errors e v . Therefore, the behavior of the trajectories of (3) may be assessed via the behavior of the latter dynamics. To this end, we use x = [x m e v ] and U = [1 N V ] in [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF] to write

x = [1 N ⊗ I n ]x m + [V ⊗ I n ]e v . (13)
So, differentiating on both sides of (4), using (3), [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF], and the fact that v l L = 0, we obtain

ẋm = F m (x m ) + G m (x m , e v ), (14)
where

F m (x m ) := [v l ⊗ I n ]F [1 N ⊗ I n ]x m and G m (x m , e v ) := [v l ⊗ I n ] F [1 N ⊗ I n ]x m + [V ⊗ I n ]e v -F [1 N ⊗ I n ]x m .
Note that F m (x m ) effectively corresponds to an "averaged" dynamics of the systems in [START_REF] Chow | Time-Scale Modeling of Dynamic Networks with Applications to Power Systems[END_REF],

F m (x m ) = N i=1 v li f i (x m ),
G m (x m , 0) = 0 and, under Assumption 2, all these functions are smooth and there exists a continuous function h :

R nN → R ≥0 such that G m (x m , e v ) ≤ h(x m , e v )|e v | ∀(x m , e v ) ∈ R nN . ( 15
)
On the other hand, by differentiating on both sides of e v = [V † ⊗ I n ]x and using ( 6), [START_REF] Isidori | Robust output synchronization of a network of heterogeneous nonlinear agents via nonlinear regulation theory[END_REF], and (13), we obtain the synchronization-errors dynamics

ėv = -σ[Λ ⊗ I n ]e v + G e (x m , e v ), (16)
where

G e (x m , e v ) := [V † ⊗ I n ]F [V ⊗ I n ]e v + [1 N ⊗ I n ]x m .
The complete collective behavior of the networked control system (3), up to the globally invertible coordinate transformation in [START_REF] Chowdhury | Persistence-based convergence rate analysis of consensus protocols for dynamic graph networks[END_REF], may be assessed by analyzing that of the interconnected systems ( 14) and [START_REF] Pogromsky | Pattern generation in diffusive networks: How do those brainless centipedes walk?[END_REF]. We see that the systems in (1) under the action of the control laws in (2) synchronize if the errors e v tend asymptotically to zero. However, the characterization of the networked systems' behavior would be incomplete unless one can ascertain what the individual systems do when they synchronize. Indeed, a priori, not even boundedness of solutions is guaranteed (whence Assumption 3). To assess any kind of stable behavior, we analyze the network system (3) on the synchronization subspace corresponding to e v = 0. On such a subspace, we have the reduced-order dynamics

ẋm = F m (x m ). (17)
So, it is clear that the motion of the synchronized systems is fully determined by that of the reduced-order dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]. In this regard, it is important to underline that [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF], as well as the "averaged" dynamics [START_REF] Martin | Time scale modeling for consensus in sparse directed networks with time-varying topologies[END_REF] are independent of the coupling gain σ. This dynamics is inherent to the network and appears simply as a consequence of the graph's connectivity imposed by Assumption 1. The synchronization dynamics [START_REF] Pogromsky | Pattern generation in diffusive networks: How do those brainless centipedes walk?[END_REF], on the other hand, clearly depends on the coupling strength σ. In this paper we are interested in investigating the synchronization behavior for 'large' values of the coupling strength. More precision about the meaning of 'large' is given farther below.

We consider two scenarii of major interest. The first pertains to the case in which the reduced-order dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] admits the origin as a globally asymptotically stable equilibrium point. This case covers, in particular, the classical problem of consensus for heterogeneous nonlinear systems interconnected over generic connected graphs, since in this case ẋm ≡ 0. Our main statement in this case (Theorem 1) is that not only the networked system achieves dynamic consensus but the origin {x = 0} is GAS for [START_REF] Silberstein | Reduction, emergence and explanation[END_REF]. The second scenario pertains to the case in which the reduced-order dynamics admits an unstable equilibrium and a stable periodic orbit. Our main statement in this case (Theorem 2) establishes sufficient conditions for almost global asymptotic stability. For instance, for a network of periodic heterogeneous nonlinear oscillators, it is possible to assess the conditions under which they synchronize and to characterize the resulting collective oscillatory behavior.

The analysis of (3) is carried out using singular-perturbations theory. For this, we write the network system (3) in the familiar singular-perturbation form [START_REF] Kokotović | Singular perturbation methods in control: analysis and design[END_REF][START_REF] Khalil | Nonlinear systems[END_REF]

ẋm = F m (x m ) + G m (x m , e v ) (18a) ε ėv = -(Λ ⊗ I n ) e v + εG e (x m , e v ), ε := 1/σ, (18b)
in which we recognize two time scales, "slow" and "fast", corresponding, respectively, to the dynamics of the averaged-unit states x m and the projected synchronization errors e v . Now, in accordance with singularperturbation theory; see [30, p. 358], the behavior of ( 18) is ineluctably determined by that of the slow dynamics, obtained by setting ε = 0, which clearly corresponds to the reduced-order model [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]. Thus, the rest of the paper is devoted to the analysis of ( 18) in the two cases evoked above.

Even though the analysis relies on the study of the system in singularly-perturbed form, [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], we remark that our two main statements are formulated for system (3), which remains the main subject study in this paper. Therefore, we finish this section by re-expressing the properties of (3) in Assumptions 1-3 in terms of [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], in the form of the following, rather evident, statement that is extensively used in the sequel. Lemma 2 Consider system (3) such that Assumptions 1-3 hold. Then, the resulting system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], with states defined in [START_REF] Teel | Smooth time-varying stabilization of driftless systems over communication channels[END_REF], enjoys the following properties:

(i) the functions F m , G m , and G e are continuously differentiable;

(ii) the origin {(x m , e v ) = (0, 0)} as an isolated equilibrium point;

(iii) the solutions to (18) are globally uniformly (in σ) ultimately bounded;

(iv) the matrix Λ is Hurwitz.

Case I: Global asymptotic stability

Consider the network system (3) in its equivalent singular-perturbation representation [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. In the case that the reduced-order dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] admits the origin as a globally asymptotically stable equilibrium, under Assumptions 1-3, global asymptotic stability of the origin for [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] follows under a sufficiently small perturbation ε > 0. According to [START_REF] Teel | Semi-global practical asymptotic stability and averaging[END_REF], we recall the practical-stability notions of a closed subset A ⊂ R n for a general nonlinear system of the form

ẋ = f (x, ε) x ∈ R n , ε ∈ [0, 1]. (19)
• The set A is globally practically attractive if, for each β > 0, there exists ε > 0 such that, for each ε ≤ ε , every solution x(t) to [START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF], there exists T > 0 such that ||x(t o + T )|| ≤ β.

• The set A is globally practically stable if there exists κ ∈ K such that, for each β > 0, there exists ε > 0 such that, for each ε ≥ ε , we have

|x(t)| ≤ κ (|x(t o )|) + β ∀t ≥ t o .
• The set A is globally asymptotically practically stable if it is globally practically attractive and globally practically stable.

Our first statement is the following.

Theorem 1 (GAS) Consider system (3) under Assumptions 1-3. In addition, assume that for system [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF], the origin {x m = 0} is globally asymptotically stable. Then, (i) the origin for ( 3) is globally asymptotically practically stable.

Furthermore, assume in addition that there exists a continuously differentiable Lyapunov function

V m : R n → R ≥0 and a class K ∞ function α such that ∂V m (x m ) ∂x m F m (x m ) ≤ -α (|x m |) 2 (20)
and there exists c r > 0 such that, for all x ∈ B r -see [START_REF] Teel | Smooth time-varying stabilization of driftless systems over communication channels[END_REF],

max x∈Br G e (x m , e v ) -G e (0, e v ) , ∂V (x m ) ∂x m ≤ c r α (|x m |) , (21)
where r > 0 is the global ultimate bound established in Lemma 1. Then, (ii) there exists σ > 0 such that, for all σ ≥ σ , the origin for ( 3) is globally asymptotically stable.

Proof. The statement in Item (i) follows by establishing global asymptotic practical stability of the origin {x = 0} for [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. To this end, we first remark that, after Assumptions 2-3 and Lemma 2, system (18) is globally uniformly (in σ) ultimately bounded. That is, there exist σ * > 0 and r > 0 such that, for any R ≥ 0, there exists τ R ≥ 0 such that

|x o | ≤ R ⇒ |x(t)| ≤ r ∀t ≥ τ R , ∀σ ≥ σ * .
Hence, it suffices to show (i) only inside the compact set B r . To do so, we show that, for every positive constant r f < r, there exists σ f > 0 such that, for each σ ≥ σ f , the compact set B r f is globally asymptotically stable for (18) on B r [START_REF] Teel | Semi-global practical asymptotic stability and averaging[END_REF]. We first rewrite the perturbed system (18) in the following form

ẋm =F m (x m ) + G m (x m , e v ) ε ėv = -(Λ ⊗ I n ) e v + εG e (0, e v ) + ε [G e (x m , e v ) -G e (0, e v )] . (22)
Next, using Assumption 2, we conclude the existence of a constant d r > 0 such that, for each x ∈ B r ,

max x∈Br {|G m (x m , e v )| , |G e (0, e v )|} ≤ d r |e v | , (23)
|G e (x m , e v ) -G e (0, e v )| ≤ d r |x m | ∀x ∈ B r . (24)
Furthermore, after Item (iv) in Lemma 2, the linear system

ėv = -(Λ ⊗ I n ) e v
is exponentially stable. As a result, there exists P ∈ R (N -1)×(N -1) symmetric and positive definite such that

P Λ + Λ P ≤ -I N .
Finally, we use the converse Lyapunov theorem in [START_REF] Khalil | Nonlinear systems[END_REF]Theorem 3.14] to conclude that, since the origin for (17) is globally asymptotically stable, there exists a function V m : R n → R ≥0 and positive constants α rr f and β rr f such that, for each x ∈ B r \B r f , we have

∂V m (x m ) ∂x m ≤ β rr f |x m | , (25)
∂V m (x m ) ∂x m F m (x m ) ≤ -α rr f |x m | 2 . ( 26
)
Therefore, under (23)-( 26), the total derivative of the Lyapunov function candidate

V (x) := V m (x m ) + V e (e v ), (27)
where V e (e v ) := e v (P ⊗ I n )e v , satisfies

V (x) ≤ -α rr f |x m | 2 + d r λ max (P ) + β rr f |x m | |e v | -[σ -d r λ max (P )] |e v | 2 ∀ x ∈ B r \B r f .
Thus, for σ f := d 2 r λ max (P ) + β rr f 2 + 2d r λ max (P ), we conclude that, for each σ ≥ σ f , we have

V (e v , x m ) ≤ - 1 2 α rr f |x m | 2 - 1 2 σ |e v | 2 ∀ x ∈ B r \B r f .
The statement in Item (ii) follows by establishing global asymptotic stability of the origin {x = 0} for (18). However, since (18) is globally uniformly (in σ) ultimately bounded, it is enough to prove that the origin is asymptotically stable on for all solutions contained in B r . So in the rest of the proof, we restrict the analysis on B r . Let us reconsider the Lyapunov function V introduced in [START_REF] Anosov | Limit cycles of systems of differential equations with small parameters in the highest derivatives[END_REF]. Note that the time derivative of V along the solutions to (18) satisfies

V (x) = -α(|x m |) 2 + ∂V m ∂x m G m (x m , e v ) -σ |e v | 2 + 2e v P [G e (0, e v ) + (G e (x m , e v ) -G e (0, e v ))] . (28)
Now, after (23) and ( 21), we conclude that

∂V (x) ∂x m G m (x m , e v ) ≤ c r d r α (|x m |) |e v | ∀ x ∈ B r , (29)
and, consequently, from (28), [START_REF] Wieland | An internal model principle is necessary and sufficient for linear output synchronization[END_REF], and (29), it follows that

V (x) ≤ -α(|x m |) 2 + c r [2λ max (P ) + d r ]α(|x m |)|e v | -[σ -2d r λ max (P )]|e v | 2 ∀x ∈ B r .
Therefore, for σ := c 2 r [λ max (P ) + d r ] 2 + 4d r λ max (P ), we conclude that, for each σ ≥ σ ,

V (x) ≤ - 1 2 α(|x m |) 2 - 1 2 σ |e v | 2 ∀ x ∈ B r .
Remark 3 The first statement in Theorem 1 can be deduced using [18, Corollaries 1 and 3], a detailed proof is included in this paper for the sake of self containedness.

•

Remark 4

The regularity conditions in (20)-( 21) are required to ensure negativity of the time derivative of the Lyapunov function V along the solutions to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. Even though they may appear conservative, since they are certainly not necessary -see [START_REF] Khalil | Nonlinear systems[END_REF]Exercise 9.24], it is important to stress that the origin is not necessarily globally asymptotically stable if these conditions do not hold. In Section 5 we provide an example that illustrates this claim. Furthermore, the inequalities (20)-( 21) are required to hold only in a compact set containing the origin. •

Case II: Periodic behavior

In this section, we present our second and main statement, which pertains to the case where (17) admits a periodic orbit that is attractive from almost all initial conditions. Under this condition, Theorem 2 below establishes that for σ > 0 sufficiently large, the network system (3) also admits a unique periodic orbit and it is globally attractive from almost all initial conditions. In particular, frequency synchronization is achieved and the synchronization errors can be made arbitrary small by choosing σ sufficiently large. It is important to stress that our main statement establishes a precise periodic behavior for the network system (3) rather that just approaching the periodic solution to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF].

For completeness and clarity, we start by recalling some notions and tools related to the stability of periodic solutions to nonlinear systems.

Consider the system ẋ = f (x)

x ∈ R n , (30)
where f : R n → R n is locally Lipschitz.

Definition 1 (Periodic solution and periodic orbit) A solution t → φ(t), or simply φ(t), to (30) is said to be α-periodic if there exists α > 0 (the period) such that, for each t ≥ 0,

φ(t + α) = φ(t) and φ(t + s) = φ(t) ∀s ∈ (0, α).
Moreover, if the system (30) admits a periodic solution φ, we say that it admits a (closed) periodic orbit γ ⊂ R n generated by the image of φ.

Then, according with Lyapunov theory, we may single out the following desired properties for periodic solutions.

Definition 2 (Orbital Stability) Let γ be a periodic orbit for [START_REF] Khalil | Nonlinear systems[END_REF].

• The orbit γ is orbitally stable if, for each ε > 0, there exist δ > 0 and T ≥ 0, such that, for each initial condition

x o satisfying |x o | γ ≤ δ, the solution φ starting from x o satisfies |φ(t)| γ ≤ ε for all t ≥ T .
• The orbit γ is orbitally asymptotically stable, if it is orbitally stable and attractive; i.e., if there exists R ∈ (0, +∞] such that, for each

x o satisfying |x o | γ ≤ R, the solution φ starting from x o satisfies lim t→∞ |φ(t)| γ = 0.
• The orbit γ is globally orbitally asymptotically stable if it is orbitally asymptotically stable with R = +∞ and almost globally asymptotically stable if it is orbitally asymptotically stable for all x o ∈ R n \D, where D ⊂ R n has a null Lebesgue measure.

Finally, we recall some orbital stability criteria in terms of the so-called characteristics multipliers [39, Section III.7] which, for linear periodic systems, are the counterpart of eigenvalues for linear autonomous systems. To see this, we assume that f is continuously differentiable and we consider the α-periodic matrix A(t) := ∂f ∂x (φ(t)), where φ(t) is the α-periodic solution to [START_REF] Khalil | Nonlinear systems[END_REF] generating the orbit γ. After Floquet theory-see e.g., [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF] and [START_REF] Perko | Differential Equations and Dynamical Systems[END_REF], there exist an α-periodic non-singular matrix P : [t o , +∞] → R n×n and a constant matrix B ∈ R n×n such that the transition matrix associated to the linear time-varying system

ẋ = A(t)x (31)
is given by X(t) := P (t)e Bt and the non-singular change of coordinates y := P (t) -x transforms the linear system (31) into ẏ = By.

Definition 3 (Characteristic multipliers)

The characteristic multipliers of the α-periodic matrix A(t) are the eigenvalues of the matrix e Bα .

Definition 4 (Non-singular periodic orbit) The periodic orbit γ generated by the periodic solution φ(t) is non-singular if the matrix A(t) := ∂f ∂x (φ(t)) admits a simple characteristic multiplier equal to 1.

Lemma 3 (Theorem 2.1, Section VI.2. [START_REF] Hale | Ordinary Differential equations[END_REF]) Consider system (30) with f continuously differentiable and let φ be a non-trivial α-periodic solution generating the orbit γ. Assume that the matrix A(t) := ∂f ∂x (φ(t)) is non-singular and all the characteristic multipliers, except one, have modulus strictly less than 1. Then, the resulting periodic orbit γ is asymptotically orbitally stable.

Sufficient conditions for orbital stability

As mentioned above, generally speaking, the standing assumption in this section is that the reduced-order dynamics (17) admits an orbitally asymptotically stable periodic solution. However, we remark that some nonlinear systems defined on compact and convex sets and that admit a limit cycle, also admit at least one equilibrium point [START_REF] Basener | The brouwer fixed point theorem applied to rumour transmission[END_REF]. This imposes particular richness to the network's collective behavior and considerable difficulty to analyze it since it translates into studying the stability of a disconnected invariant set. In that light, we pose the following, more precise hypothesis.

Assumption 4

The reduced-order dynamics (17) admits a unique compact invariant subset ω ⊂ R n that is globally attractive; namely, for each x mo ∈ R n , the solution x m (t) starting from x mo satisfies lim sup

t→+∞ |x m (t)| ω = 0 (32)
Furthermore, the set ω is compose a non-singular periodic orbit γ o , of period α o , that is orbitally asymptotically stable and

• the origin {x m = 0}, when the latter is repulsive.

• the homoclinic orbit γ 1 := W u o (0) ∩ W s o (0)
when the origin is hyperbolic 1 , where W s o (0) and W u o (0) are, respectively, the global stable and unstable manifolds of the origin.

Remark 5 Note that the global attractivity property in [START_REF] Bıyık | Area aggregation and time-scale modeling for sparse nonlinear networks[END_REF] plus the structure of the invariant set ω imply the existence of a Lyapunov function enjoying useful properties along the solutions to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]; see Lemma 11 in Appendix I.

• Lemma 4 Under Assumption 4 and Item (i) in Lemma 2, the periodic orbit γ o is almost globally orbitally asymptotically stable for [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF].

We are ready to present our main statement.

Theorem 2 (Almost global orbital asymptotic stability) Consider the network system (3) under the Assumptions 1-3 and such that the reduced-order dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] satisfies Assumption 4. Then, there exists σ f > 0, such that, for all σ ≥ σ f , the networked system (3) admits a unique nontrivial periodic orbit O 1/σ , of period α 1/σ and that is almost globally orbitally asymptotically stable. Moreover, O 1/σ → O o , where

O o := {x ∈ R nN : x m ∈ γ o and e v = 0}
-see [START_REF] Teel | Smooth time-varying stabilization of driftless systems over communication channels[END_REF], and

α 1/σ → α o as σ → ∞.
The proof of Theorem 2, which is provided farther below, follows a sequence of logical steps to assess the existence, uniqueness, and almost global orbital asymptotic stability of an orbit for [START_REF] Silberstein | Reduction, emergence and explanation[END_REF]. The analysis relies on studying the singularly-perturbed system ( 18), but we emphasize that the available literature on stability (of the origin or a compact set) for singularly-perturbed systems [START_REF] Tikhonov | Systems of differential equations containing small parameters in the derivatives[END_REF], [START_REF] Kokotović | Singular perturbations and order reduction in control theory-an overview[END_REF], [START_REF] Khalil | Nonlinear systems[END_REF] does not apply to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], when (17) admits a limit cycle and an isolated equilibrium point. Therefore, the proof of Theorem 2 relies on technical lemmata that are presented next, but, for clarity of exposition, the proofs of these lemmata are given in the Appendix. • Lemma 6 establishes that, given a torus sufficiently tight around O o , for each coupling gain σ sufficiently large, system (17) admits a unique periodic orbit O 1/σ contained in the torus;

• Lemma 7 establishes that each such orbit O 1/σ is (locally) asymptotically stable and admits the aforementioned torus as a basin of attraction;

• Lemmata 8 and 18 (in the Appendix) provide a local analysis around the origin, to establish that it attracts only the solutions starting from a null-measure set.

Technical Lemmata

We start by introducing the following notations. Correspondingly to γ o ⊂ R n and γ 1 ⊂ R n , which denote, respectively, the closed periodic and homoclinic orbits for system [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]-see Assumption 4, we introduce their "lifting" Γ o ⊂ R nN and Γ 1 ⊂ R nN in the space of system (18), as

Γ o := (x m , e v ) ∈ R nN : x m ∈ γ o and e v = 0 , (33)
Γ 1 := (x m , e v ) ∈ R nN : x m ∈ γ 1 and e v = 0 . (34)
Furthermore, we denote by T ρ the torus defined as

T ρ := (e v , x m ) ∈ R N (n-1) × R n : |(x m , e v )| Γo ≤ ρ , (35)
and we use Γ ε ⊂ R nN to denote a closed orbit generated by a periodic solution to system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], if it exists.

That is, Γ ε is a subset in the space of (x m , e v ) that consists in the image points generated by the parameterized solutions of (18), t → (x m (t), e v (t)), that are periodic with period α ε . The first technical lemma provides a statement for system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] on global practical asymptotic (GApS) stability of the set

Ω := Γ o ∪ Γ 1 ,
where Γ o and Γ 1 are introduced in ( 33)- [START_REF] Rejeb | Control design with guaranteed cost for synchronization in networks of linear singularly perturbed systems[END_REF].

Lemma 5 (GApS of Ω) Consider system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] such that Items (i)-(iii) of Lemma 2 hold and let Assumption 4 be satisfied for the corresponding reduced-order system [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]. Then, the set Ω := Γ o ∪ Γ 1 is GApS for [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. In particular, for any ρ > 0, there exists ε 1 (ρ) > 0, such that, for each ε ≤ ε 1 (ρ) and for each initial condition x o = (e vo , z mo ) ∈ R nN , the solution (x m (t), e v (t)) satisfies lim t→∞ |x m (t), e v (t)| Ω ≤ ρ.

The next lemma establishes that, for all sufficiently small values of ρ > 0, there exist sufficiently small values of ε, such that there exists a unique periodic orbit Γ ε ⊂ T ρ generated by a solution to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] of period

α ε ≈ α o .
Lemma 6 (Existence of Γ ε ) Consider system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] such that Items (i) and (iii) of Lemma 2 hold and let Assumption 4 hold for the reduced-order dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]. Then, there exist ρ o > 0 and a class K function ε o such that, for each ρ ∈ (0, ρ o ] and for each ε ≤ ε o (ρ), system (18) has a unique nontrivial periodic orbit Γ ε , which is strictly contained in T ρ . Moreover, the period α ε of the solution to (18) generating the orbit Γ ε tends to α o , which is the period of the solution to (17) generating the orbit Γ o .

Remark 6

The existence result in Lemma 6 follows from a direct application of Anosov Theorem-see Lemma 10 in the Appendix.

•

The next lemma establishes local asymptotic orbital stability of all periodic orbits Γ ε lying inside the torus T ρ for sufficiently small values of ε and ρ. Moreover, we show that the corresponding domain of attraction is uniform in ε.

Lemma 7 (Stability of Γ ε ) Let system (18) satisfy Items (i) and (iii) of Lemma 2 and let Assumption 4 be satisfied for the corresponding reduced-order dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]. Then, there exist ε * * > 0 and ρ * * > 0 such that, for each ε ≤ ε * * , each periodic orbit Γ ε ⊂ T ρ * * generated by an α ε -periodic solution to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], with α ε sufficiently close to α o , is asymptotically orbitally stable with a domain of attraction that contains T ρ * * .

Remark 7 Lemma 7 is reminiscent of a statement established by Anosov-[27, Theorem 5]-that pertains to the case in which the periodic orbit γ o for (17) is only non-singular (or hyperbolic). Although it is claimed in [START_REF] Anosov | Limit cycles of systems of differential equations with small parameters in the highest derivatives[END_REF] that the proof therein translates directly to the case where γ o is non-singular and asymptotically stable, in this paper, we provide an original proof for the latter case using the theory of perturbed matrices [START_REF] Moro | On the lidskii-vishik-lyusternik perturbation theory for eigenvalues of matrices with arbitrary jordan structure[END_REF][START_REF] Stewart | Matrix Perturbation Theory[END_REF].

•

The next statement links those from Lemmata 5-7. It establishes that if ε is sufficiently small, then the periodic behavior of the reduced order system ( 17) is preserved for system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] as well as its stability properties.

Proposition 1 Consider the dynamical system (18) under the assumption that Items (i)-(iv) of Lemma 2 hold and assume further that the reduced-order dynamics (17) satisfies Assumption 4.

Then, there exists ρ o > 0 such that, for each ρ ∈ (0, ρ o ], there exists ε 2 (ρ) > 0 such that, for each ε ∈ (0, ε 2 (ρ)], (i) system (18) admits a unique orbit Γ ε ⊂ T ρ generated by a non-trivial (α ε )-periodic solution, with

Γ ε → Γ o and α ε → α o as ε → 0;
(ii) Γ ε is (locally) asymptotically stable;

(iii) for any initial condition xo ∈ R nN the corresponding solution x(t) to (18) either converges to

Γ ε or to a ρ-neighborhood of Γ 1 , that is, lim sup t→∞ |x(t)| Γ1 ≤ ρ. (36)
Proof. Items (i)-(iv) in Lemma 2 and Assumption 4 imply that the statements of Lemmata 5-7 hold. Then, let Lemma 7 generate ε * * > 0 and ρ * * > 0. Furthermore, let Lemma 6 generate

(ρ o , ε o (•)), for each ρ ∈ (0, ρ o ], let Lemma 5 generate (ε 1 (ρ), ε 1 (ρ * * )), and let ε ≤ ε 2 (ρ) := min{ε * * , ε o (ρ), ε o (ρ * * ), ε 1 (ρ), ε 1 (ρ * * )}
be arbitrarily fixed.

After Lemma 6, there exists a unique periodic orbit Γ ε ⊂ T ρ ∩ T ρ * * = T min{ρ,ρ * * } generated by a solution to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. Now, given any sequence { i } ∞ i=1 that converges to zero and such that i ≤ ε 2 (ρ) for all i ∈ {1, 2, ...}, from the above, we know that, for i large enough, the unique orbit Γ i satisfies Γ i ⊂ T ρi , where

ρ i := ε -1 o ( i )- note that ε -1
o exists and is of class K because so it ε o . Item [i] of the proposition follows since T ρi converges to Γ o -see [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF].

Next, after Lemma 7, we conclude that Γ ε is orbitally asymptotically stable and T ρ * * is inside the domain of attraction of Γ ε . This establishes Item (ii).

Finally, from Lemma 5 we conclude that each solution to (18) either converges to T min{ρ,ρ * * } ⊂ T ρ * * , so it also converges to Γ ε , or it converges to a min{ρ, ρ * * }-neighborhood of Γ 1 . This establishes Item (iii).

The last technical lemma provides a local stability analysis around the origin of [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. It states that the origin is a hyperbolic equilibrium point, for all sufficiently-small values of ε. Furthermore, inspired by the Stable Manifold Theorem [START_REF] Coddington | Theory of ordinary differential equations[END_REF]Theorem 13.4.1], we show that the stable and unstable manifolds around the origin are uniquely defined on a neighborhood whose size does not shrink with ε.

Lemma 8 (Local behavior around the origin) Consider system (18) and let Items (i)-(ii) of Lemma 2 hold. Assume further that the corresponding reduced-order dynamics (17) satisfies Assumption 4. Then, there exist ρ * > 0, ε * > 0, a neighborhood of the origin denoted U ⊂ R n , and r > 0 such that, for each ε ∈ (0, ε * ), (i) system (18) admits a unique unstable and stable manifolds (W u ε (0), W s ε (0)) defined on U ; (ii) for each x(t) bounded solution to (18) starting from xo ∈ U \W s ε (0), there exists T 1 > 0 such that |x(t)| ≥ r for all t ≥ T 1 ;

(iii) for each x(t) solution to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] such that |x(0)| Γ1 ≤ ρ * , there exists T 2 > 0 such that |x(T 2 )| < r.

Proof of Theorem 2

Under Assumptions 1-3, Items (i)-(iv) of Lemma 2 hold for system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. This and Assumption 4 imply that the statements of Proposition 1 and Lemmata 5-8 hold. Then, let Lemma 8 generate ρ * > 0 and ε * > 0 and let Proposition 1 generate (ρ o , ε 2 (min{ρ * , ρ o }/4)). We show that the statement of Theorem 2 holds with σ f := 1/ min {ε * , ε 2 (min{ρ * , ρ o }/4)} in the following four steps:

1. First, for any σ ≥ σ f or, equivalently, any ε = 1/σ satisfying ε ≤ min {ε * , ε 2 (min{ρ * , ρ o }/4)}, we use Item (i) of Proposition 1 to conclude the existence of a unique nontrivial periodic orbit Γ ε generated by a periodic solution to (18) of period α ε . From Item (ii) of the same Proposition it follows that Γ ε is locally asymptotically stable. In addition, from Item (iii) of Proposition 1 it follows that each solution x(t) to system (18) either converges to the orbit Γ ε , otherwise, it converges to a (min{ρ * , ρ o }/4)neighborhood of Γ 1 ; that is, [START_REF] Liu | Output synchronization of dynamical networks with incrementallydissipative nodes and switching topology[END_REF] holds with ρ = min{ρ * , ρ o }/4 and, consequently, there exists T < ∞ such that

|x(t)| Γ1 = |(x m (t), e v (t))| Γ1 ≤ min{ρ * , ρ o }/2 ∀t ≥ T. (37)
2. Let us now introduce the backward propagation of the stable manifold W s ε (0) introduced for (18) in Lemma 8. That is, we introduce set R(W s ε (0)) := {x(t) : t ≤ 0, x(0) ∈ W s ε (0)} and prove by contradiction that the solution x to (18) satisfying (37) must start from the set R(W s ε (0)). Indeed, assume that the opposite holds. Then, using Item (iii) in Lemma 8, we conclude that the solution x must enter B r at some T * ≥ T . In particular, x(T * ) ∈ U and x(T * ) / ∈ W s ε (0). Now, using Item (ii) in Lemma 8, we conclude that there exists T 1 > 0 such that

|x(T * + t)| ≥ r ∀t ≥ T 1 .
However, since |x(T * + T 1 )| Γ1 ≤ min{ρ * , ρ o }/2, it follows that x must enter B r again under Item (iii) of Lemma 8, which contradicts Item (ii).

3. Next, we show that the set R(W s ε (0)) is a null measure set using contradiction. That is, let S o ⊂ R(W s ε (0)) such that µ(S o ) = 0. Assume without loss of generality that for some T < 0, we have

S o ⊂ {x(t) : t ∈ [-T, 0], x(0) ∈ W s ε (0), x solution to (18)}.
Note that R b (T, S o ) := {x(T ) : x(0) ∈ S o } ⊂ W s ε (0) with µ(W s ε (0)) = 0. However, using Lemma 18 from the Appendix, we conclude that µ(S o ) = 0.

4. Finally, using the inverse transformation [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF], it follows that the orbit

O ε := {x ∈ R nN : (x e , e v ) ∈ Γ ε }
is almost GAS for (3). The second statement follows from Lemma 6.

Case study: a network of Andronov-Hopf oscillators

For illustration of our main theoretical findings, we address, as a case-study, the analysis of a network of Andronov-Hopf oscillators. The equation of such oscillator represents a normal form of the bifurcation carrying the same name and is given by ż = -ν|z| 2 z + µz,

where z ∈ C denotes the state of the oscillator and ν, µ ∈ C are constant parameters; ν := ν R + jν I and µ := µ R + jµ I . The analysis of ( 38) is well documented in the literature. For instance, via Lyapunovexponents-based methods, as in [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] and [START_REF] Perko | Differential Equations and Dynamical Systems[END_REF], or using Lyapunov's direct method, as in [START_REF] Matthews | Phase diagram for the collective behavior of limit-cycle oscillators[END_REF] and [START_REF] Pham | Stable concurrent synchronization in dynamic system networks[END_REF]. The behavior on the phase plane is illustrated in Fig. 1.

Figure 1: Trajectories of the Andronov-Hopf oscillator on the complex plane when µ R > 0. In this case, the origin is an unstable equilibrium and all trajectories tend to the stable limit-cycle of radius √ µ R and period 2π µ I . Now, for the purpose of this paper, we consider N forced Andronov-Hopf oscillators for which, for simplicity, we assume that ν = 1, the general case where ν = 1 will be considered elsewhere. That is, consider the systems

żi = -|z i | 2 z i + µ i z i + u i i ∈ {1, 2, . . . , N } (39)
where z i ∈ C denotes the i-th oscillator's state, u i ∈ C is its control input, and µ i := µ Ri + iµ Ii is a constant complex parameter. Now we design each input u i as in (2) while replacing (x i , x j ) therein by (z i , z j ). Then, in compact form, the overall network dynamics in closed loop takes the form

ż = F (z) -σLz, (40)
where

z := [z 1 • • • z N ] and the function F : C N → C N is given by F (z) := [f 1 (z 1 ) f 2 (z 2 ) • • • f N (z N )] , f i (z i ) := -|z i | 2 z i + µ i z i .
Next, using [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF] and the coordinate transformation defined in [START_REF] Teel | Smooth time-varying stabilization of driftless systems over communication channels[END_REF], we can rewrite the network dynamics in the singular perturbation form, where the reduced-order dynamics defined by [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]

has the form żm = f m (z m ) := -|z m | 2 z m + µ m z m , (41)
where z m = v l z, the parameter µ m ∈ C corresponds to the weighted average of the µ i s, i.e.,

µ m := µ mR + iµ mI , µ mR := N j=1 v lj µ Rj , µ mI := N j=1 v lj µ Ij .
Hence, the reduced-order dynamics ( 41) is an Andronov-Hopf oscillator, but we note that µ mR and µ mI do not correspond to simple averages of the individual parameters µ Ri and µ Ii ; through the Laplacian's eigenvector v l , they depend on the network's topology.

Two cases that are of interest are studied next.

Case in which µ mR = 0

The parameter µ mR , which depends both on the individual parameters of the oscillators and on the structure of the network, is particularly important. Indeed, as we shall show, for σ sufficiently large, the asymptotic collective behavior of the network is uniquely defined by the sign of µ mR .

Proposition 2 Consider the network of Andronov-Hopf oscillators defined in [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF] interconnected over a connected directed graph. Then, the following statements hold true:

(i) If µ mR < 0 then there exists σ f > 0 such that, for all σ ≥ σ f , the origin {z = 0} is GAS.

(ii) If µ mR > 0 then there exists σ f > 0 such that, for each σ ≥ σ f , system (40) has a unique nontrivial periodic orbit which is almost globally asymptotically stable.

Proof.

The network model [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF] satisfies Assumptions 1-3 regardless of the sign of µ mR . Indeed, Assumption 1 is explicit in the Proposition. On the other hand, each single unit has a polynomial vector field, so Assumption 2 holds. To see that Assumption 3 on semi-passivity of each unit in (39) holds, we evaluate the total derivative of the Lyapunov function candidate

V (z i ) := 1 2 z * i z i ,
where z * i is the complex conjugate of z i , to obtain

V (z i ) = -|z i | 4 + µ Ri |z i | 2 + z * i u i .
Now, to complete the proof, we note that if µ mR < 0, it follows that V (z m ) is negative definite and the origin for the reduced-order dynamics (41) is globally asymptotically stable. Hence, after Theorem 1, Item (i) of the proposition follows.

In the opposite case, if µ mR > 0, the total derivative of V (z m ) along the solutions to (41) satisfies

V (z m ) = -|z m | 4 + µ mR |z m | 2 . ( 42
)
Hence, system (41) admits a compact invariant set that is composed of two disjoint invariant subsets, that is

ω = {z m ∈ C : |z m | = √ µ mR } ∪ z m = 0 .
Moreover, the invariant orbit

γ o := {z m ∈ C : |z m | = √ µ mR }
is almost globally asymptotically stable and the origin {z m = 0} is anti-stable -see [START_REF] Panteley | On the stability and robustness of Stuart-Landau oscillators[END_REF]. As a result, Item (ii) of the proposition follows using Theorem 2.

Case in which µ mR = 0

When µ mR = 0, we use [START_REF] Basener | The brouwer fixed point theorem applied to rumour transmission[END_REF] to conclude that the reduced-order dynamics (41) admits the origin as a GAS equilibrium. Hence, using Theorem 1, we conclude GApS of the origin for [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF]. Now, we would like to answer the following two questions:

1. Are there cases where, although being globally practically asymptotically stable, the origin for ( 40) is locally unstable?

2. Under which conditions, the origin for ( 40) is locally exponentially stable?

3. Are we able, to establish global exponential stability?

Consider the linearization of (40) around the origin, which is given by

ż = A σ z, A σ := A o -σL, A o := blkdiag {µ 1 , µ 2 , ..., µ N } . ( 43
)
To answer the first question, in the following example, we provide a case where the matrix A σ has an eigenvalue with a strictly positive real part for all σ > 0 sufficiently large. Hence, we provide a positive answer to the question in 1).

Example 1 Consider a network of two Andronov-Hopf oscillators, N = 2, and let µ R1 = -µ R2 = µ, and µ I1 = µ I2 = 0. Clearly, in this case, µ mR = 0. Note that in this particular case

A σ = µ 0 0 -µ -σL.
Next, we let for simplicity µ := 1 and L := 1 -1 -1 1 , so

A σ = -σ + 1 σ σ -σ -1 .
Direct calculations show that the eigenvalues of the matrix A ε are given by

λ 1,2 (A σ ) = -σ ∓ σ 2 + 1
and therefore the matrix A σ has 2 positive eigenvalues for all σ > 0. Hence, linearization of the network dynamics is unstable and therefore (see Theorem 4.7, Khalil, 3rd) the origin is unstable for the network of Andronov-Hopf oscillators.

Next, to answer the second question, we consider the following assumption.

Assumption 5 There exists σ > 0 such that, for all σ ≥ σ , the smallest (in norm) eigenvalue of A σ , denoted by λ 1 (A σ ), has a strictly negative real part.

When Assumption 5 holds, we show that the origin for ( 40) is locally exponentially stable.

Proposition 3 Consider the network of Andronov-Hopf oscillators in [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF], interconnected over a connected directed graph, and such that Assumption 5 holds and µ mR = 0. Then, there exists σ f > 0 such that, for all σ ≥ σ f , the origin {z = 0} is locally exponentially stable.

In the following example, we modify Example 1 so that Assumption 5 holds.

Example 2 Consider a network of two Andronov-Hopf oscillators, N = 2, and let µ R1 = -µ R2 = 1, and

µ I1 = -µ I2 = 2.
Clearly, in this case, µ mR = 0. Note that in this particular case

A σ = 1 -2i 0 0 -1 + 2i -σL.
Next, we let L := 1 -1 -1 1 to obtain

A σ = -σ + 1 -2i σ σ -σ -1 + 2i .
Direct calculations show that the characteristic polynomial of A σ is given by

P (λ) := λ 2 + 2σλ + (3 + 4i).
The roots of the polynomial P are the eigenvalues of the complex matrix A σ and are given by

λ 1,2 := -σ ± σ 2 -3 -4i.
As a result, the real parts of λ 1,2 are given by

(λ 1,2 ) := -σ ± (σ 2 -3) + (σ 2 -3) 2 + 16 √ 2 .
Clearly, (λ 1,2 ) are strictly negative for σ > 0 sufficiently large. Hence, linearization of the network dynamics is exponentially stable.

The third question is an open question to the best of our knowledge. Indeed, since we already know that the origin for ( 40) is GApS, one way to give a positive answer for the third question is by showing that the origin for [START_REF] Floquet | Sur les quations diffŕentielles lináires à coefficients périodiques[END_REF], which is already locally exponentially stable under Assumption 5, has a basin of attraction that does not shrink as σ get bigger or shrinks at a specific rate. However, we are not, so far, able to establish this property.

Conclusion

We presented a framework to analyse networks of heterogeneous nonlinear systems. Our approach allows to qualitatevely characterize the collective behavior for "large" values of the coupling gains. The proposed approach, however, does not give much information about the quantitative properties of such behavior. In particular, we do not provide explicit values of the coupling strength under which the networked system exhibit the established behavior. Characterizing the emergent behavior, such as orbital asymptotic stability, both qualitatively and quantitatively (in terms of the coupling gain) is still an open problem. Furthermore, beyond the analysis problems solved in this paper, the control design problem is widely open. To find conditions under which a network of heterogeneous systems may be controlled so that it admits a desired reduced-order dynamics. Finally, we believe that extending the proposed framework for general classes of nonlinear systems such as hybrid systems is an interesting perspective as well.

has unique solutions on [0, t 1 ] when starting from B r . Furthermore, there exists r 1 ≥ r > 0 such that, for each initial condition zo ∈ B r , the solution z(t) to (45) starting from zo satisfies

|z(t)| ≤ r 1 ∀t ∈ [0, t 1 ].
• the origin of the boundary-layer system

y = g(z, y + h(z), 0) (46)
is exponentially stable, uniformly in z.

Then, there exist µ * > 0, ε * ∈ (0, ε o ], and M > 0 such that, for each ε ∈ (0, ε * ), the solution (z(t), e(t)) to (44) starting from (z o , e o ) ∈ B r × B ρ , the solution z(t) to (45) starting from z o , and the solution y(t) to [START_REF] Coddington | Theory of ordinary differential equations[END_REF]

starting from e o -h(z o ) with |e o -h(z o )| ≤ µ * satisfy, for each t ∈ [0, t 1 ], |z(t) -z(t)| ≤ M ε e(t) -h(z(t)) -y t ε ≤ M ε.
Remark 8 Note that system (18) verifies all smoothness conditions in the first item of Lemma 9 after Item (i) in Lemma 2. We also know that the trajectories of both ( 18) and (17) exist and are uniquely defined on the interval [0, ∞]. Furthermore, the boundary-layer model [START_REF] Coddington | Theory of ordinary differential equations[END_REF], in the case of [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], is given by the linear system ėv = -(Λ ⊗ I n )e v , which is exponentially stable under Item (iv) in Lemma 2. Finally, under lemma 1, we know the existence of r ≥ 0 such that the set

B r := (x m , e v ) ∈ R n × R n(N -1) : |(x m , e v )| ≤ r
is globally attractive. Thus, all the Items in Lemma 9 are satisfied. •

Existence of periodic orbits for singularly perturbed systems

The following result, which is a consequence of the main statements in [START_REF] Anosov | Limit cycles of systems of differential equations with small parameters in the highest derivatives[END_REF], establishes the existence of periodic solutions for singularly-perturbed systems.

Lemma 10 Consider the singularly perturbed system (44) such that the following properties hold:

1. the functions f and g are continuous with respect to (z, e, ε) and differentiable with respect to z and e. Moreover, the derivatives of f and g with respect to z and e depend continuously on (z, e, ε).

2. there is a unique function h : R mz → R me such that g(z, h(z), 0) = 0.

3. the equilibrium state y = 0 of (46) is hyperbolic uniformly in z.

4. the unperturbed system (45) has a nontrivial nonsingular periodic orbit γ o ⊂ R mz .

Then, there exists ρ o > 0 and a class K function ε o such that for each ρ ∈ (0, ρ o ] and for each ε ≤ ε o (ρ), the system (44) has a unique nontrivial periodic orbit Γ ε , which is strictly contained in the ρ-neighborhood of Γ o , where

Γ o := {(z, e) ∈ R mz × R me : z ∈ γ o and e = h(z)}.
Moreover, the period α ε of the periodic solution to (44) generating the orbit Γ ε tends to α o the period of the solution to (45) generating the orbit Γ o .

Eigenvalues of linearly perturbed matrices

In this section, we recall the general result due to Lidskii-Vishik-Lyusternik [44, Theorem 2.1], where the first order expansion of the eigenvalues of the matrix (A + εB), with respect to ε, and the explicit formulas for the leading coefficients are given. Before recalling the result, we start introducing few notions. Given a complex matrix A ∈ C n×n , we assume without loss of generality that A has a single multiple eigenvalue λ ∈ C. Hence, the Jordan form of A is given by J := blkdiag J 11 , ..., J 1r1 , ......, J q1 , ..., J qrq .

Note that, for each j ∈ {1, 2, ..., q}, the blocks {J ji } rj i=1 have the same size denoted n j . We assume without loss of generality that n 1 < n 2 < ... < n q . Furthermore, for each J ji , i ∈ {1, 2, ..., r j } and j ∈ {1, 2, ...., q}, distinct right and left eigenvectors are associated, which we denoted by p ji ∈ C n and q ji ∈ C n , respectively. To know which right and left eigenvectors are associated to the block J ji , we use the nonsingular matrices T , U ∈ R n such that J = U AT , and use the decompositions

T = [T 11 , ..., T 1r1 , ......, T q1 , ..., T qrq ], U = [U 11 , ..., U 1r1 , ......, U q1 , ..., U qrq ].
As a result, we pick p ji as the first column of T ji and q ji as the last row U ji . Now, Collecting the aforementioned vectors allows us to introduce the following matrices

P j := [p j1 , ..., p jrj ], Q j :=    q j1 . . . p jrj    j ∈ {1, 2, ..., q} and
P s := [P 1 , ..., P s ], Q s :=    Q 1 . . . Q s    s ∈ {1, 2, ..., q}.
Finally, we let E 1 := I and E s := 0 0 0 I for all s ∈ {2, 3, ..., q}.

Lemma 12 (Lidskii Theorem) Given j ∈ {1, 2, ..., q}, if the matrix Q j-1 BP j-1 is nonsingular. Then, there are r j n j eigenvalues of A + εB admitting a first order expansion

λ l jk (A + εB) = λ + ζ 1 n j jk ε 1 n j + o ε 1 n j k ∈ {1, 2, ..., r j }, l ∈ {1, 2, ..., n j },
where {ζ jk } rj k=1 is the root of the equation det (Q j BP j -ζE j ) = 0.

The invariant subspaces of perturbed matrices

The set X ⊂ R n×n is an invariant subspace for the square matrix

A ∈ R n×n if Ax ∈ X ∀x ∈ X .
Furthermore, we let a matrix X ∈ R n×k , k ∈ {1, 2, ..., n}, whose columns form a basis for X . It is well established that there is a unique matrix L ∈ R k×k such that AX = XL. The matrix L is the representation of A on X with respect to the basis X, and the eigenvalues of L are eigenvalues of A. Now, we let the columns of X form an orthonormal basis for X . Let the matrix [X Y ] be orthonormal. Consequently, we obtain that

[X Y ] A[X Y ] = L X AY 0 L c , L c := Y AY. (49)
Note that the set X is said to be a simple invariant subspace if the intersection between the eigenvalues of L and L c is empty. The following result is recalled from [45, Theorem 1.5. Page 224].

Lemma 13 (Spectral resolution) Let X be a simple invariant subspace having the form (49) with respect to the matrix (X, Y ). Then, there exist matrices

X c ∈ R n×(n-k) and Y c ∈ R n×k such that [X X c ] -= [Y Y c ] and [X X c ] -A[X X c ] = L 0 0 L c .
Next, we let X be a simple invariant subspace of A and let a matrix B ∈ R n×n . We show that when A -B is sufficiently small, then there is a unique invariant subspace Y of B that approaches X as |A -B| goes to zero. The following result is recalled from [52, Theorem 15.5.1].

Lemma 14 (Continuity of invariant subspaces) For a matrix A ∈ R n×n admitting an invariant subspace X , the following properties are equivalent:

• The invariant subspace X is simple.

• For each ε > 0, there exists δ > 0 such that, each B ∈ R n×n , with |A -B| ≤ δ, has a unique invariant subspace Y satisfying dist(X , Y) ≤ ε, where dist(•, •) is the Grassmann distance between linear subspaces [START_REF] Ye | Schubert varieties and distances between subspaces of different dimensions[END_REF].

Next, we analyze the smooth variation of the invariant subspace X when the corresponding matrix A depends smoothly on its parameters. The following lemma can be found in [START_REF] Gohberg | Invariant Subspaces of Matrices with Applications[END_REF]Theorem 18.7.2].

Lemma 15 (Analytic families of invariant subspaces) Let A : C → C n×n be analytic. Given z o ∈ C, we let Γ be a contour in the complex plan such that the intersection between Γ and the eigenvalues of A(z o ) is empty. Then, the invariant subspace X (z) of A(z) corresponding to the eigenvalues of A(z) lying inside Γ is analytic on a neighborhood of z o .

The stable manifold theorem

Consider the dynamical system ẋ = Ax + g(x)

x ∈ R n ,

where A has n eigenvalues {λ 1 (A), λ 2 (A), ..., λ n (A)} with (λ j (A)) = 0 ∀j ∈ {1, 2, ..., n},

i.e., the origin x = 0 is a hyperbolic equilibrium for [START_REF] Panteley | On the stability and robustness of Stuart-Landau oscillators[END_REF]. Moreover, the function g : R n → R n is continuously differentiable, and there exists κ ∈ K such that

|g(x)| ≤ κ(|x|)|x|. (52)
Next, we let (E s (A), E u (A)) be the spectral invariant subspaces associated with the stable and the unstable eigenvalues of A, respectively. Note that there exist positive constants (c s (A), c u (A)) (called, respectively, forward and backward overshoots) such that

|e tA x s | ≤ c s (A)e -rs(A)t |x s | ∀t ≥ 0, ∀x s ∈ E s (A), |e tA x u | ≤ c u (A)e ru(A)t |x u | ∀t ≤ 0, ∀x u ∈ E u (A), (53)
The following result and its proof can be found in [START_REF] Birnir | Dynamical systems theory[END_REF]Theorem 4.1].

Lemma 16 (The stable manifold theorem) Consider system (50) such that g is continuously differentiable and (51)-( 52) hold. Then, on the neighborhood of the origin B γ , where γ > 0 is chosen so that, for some ∆ > 0, we have

c s (A)γ + κ(γ) c s (A) r s (A) -µ(A) + c u (A) r u (A) + µ(A) ∆ ≤ ∆,
with µ(A) := min{r u (A), r s (A)}/2, there exists a continuously differentiable function h s : Proj E s (A) (B γ ) → E u (A) such that the stable manifold W s (0) is the graph of h s and the following properties hold:

1. h s (0) = 0 and ∂h s ∂xs (0) = 0, 2. x(0) ∈ W s (0) =⇒ x(t) ∈ W s (0) for all t ≥ 0, 3. when x(0) ∈ B γ \W s (0) and x(t) is bounded, it follows that there exists t 1 > 0 such that

|x(t)| ≥ γ/2 ∀t ≥ t 1 .
Appendix II: Auxiliary Lemmata Given a function of two scalar variables that is smooth in one and only continuous in the other, the following original lemma shows the existence of a smooth approximation to any given 'nonuniform' degree of precision.

Lemma 17 Consider a function T : [0, 1] × [0, α o ] → R n×n such that τ → T (ε, τ ) is continuous, ε → T (ε, τ )
is continuously differentiable, and τ → T (0, τ ) is continuously differentiable. Then, for each ρ > 0, there exists T :

[0, 1] × [0, α o ] → R n×n continuously differentiable such that T (0, τ ) = T (0, τ ) ∀τ ∈ [0, α o ], (55)
| T (ε, τ ) -T (ε, τ )| ∞ ≤ ρε + o(ε) ∀τ ∈ [0, α o ], (56)
and

lim ε→0 Ṫ (ε, τ ) = lim ε→0 ∂ T (ε, τ ) ∂τ = Ṫ (0, τ ). ( 57
)
Proof. Since the matrix T is continuously differentiable in ε and continuous in τ , then it admits a first-order Taylor expansion of the form

T (ε, τ ) = T (0, τ ) + a(τ )ε + g(ε, τ ),
where a : [0, α o ] → R n×n is continuous and g : [0, 1] × [0, α o ) → R n×n enjoys the same continuity and smoothness properties as T . Furthermore, there exists M > 0 such that, for each τ ∈ [0, α o ], we have

|g(ε, τ )| ≤ M ε 2 ∀ε ∈ [0, 1].
Now, we choose the matrix T as T (ε, τ ) = T (0, τ ) + â(τ )ε,

where â : [0, α o ] → R n×n is a continuously differentiable approximation of a on [0, α o ] satisfying sup τ ∈[0,αo] {|a(τ ) -â(τ )|} ≤ ρ.
Note that to obtain the latter inequality, we used Stone-Weierstrass theorem stating that every continuous function defined on a closed interval [0, α o ] can be uniformly approximated as closely as desired by a polynomial function [START_REF] Stone | The generalized weierstrass approximation theorem[END_REF].

However, since the right-hand side in [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] is continuously differentiable, we the continuous dependence of the solutions on the initial data [56, Theorem V.2.1] to conclude that, for each > 0, there exists i ∈ {1, 2, ...} such that, for each i ≥ i , we have

|y mi (t) -x m (t + t i )| ≤ ∀t ∈ [0, T ].
The latter implies that x m (t i ) must lie inside B r , which contradicts (59).

Proof of Lemma 5

We start noting, by definition of the sets Ω and ω, that (x m , e v ) ∈ Ω if and only if e v = 0 and x m ∈ ω. As a result,

|(x m , e v )| Ω ≤ |e v | + |x m | ω .
Then, to prove the lemma, we show that ( ) there exists ε 1 ∈ K such that, for each ρ > 0, for each ε ≤ ε 1 (ρ), and for each initial condition (x mo , e vo ) ∈ R n × R n(N -1) , the solution (x m (t), e v (t)) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] satisfies

lim t→∞ (|e v (t)| + |x m (t)| ω ) ≤ ρ.
In turn, to prove ( ), we first show that ( ) there exists ε v ∈ K such that, for each ρ v > 0, for each ε ≤ ε v (ρ v ), and for each initial condition (x mo , e vo ) ∈ R n × R n(N -1) , the solution (x m (t), e v (t)) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] satisfies

lim t→∞ |e v (t)| ≤ ρ v . (60)
To prove ( ), we use the following four properties:

• The dynamics of e v in (18b) can be expressed as in [START_REF] Panteley | Practical dynamic consensus of Stuart-Landau oscillators over heterogeneous networks[END_REF]; namely, in the following form

ε ėv = -(Λ ⊗ I n ) e v + εG e (0, e v )
+ ε [G e (x m , e v ) -G e (0, e v )] .

• By assumption, system (18) is globally (uniformly in σ := 1/ε) ultimately bounded. That is, there exists ε * > 0 and r > 0 such that, for any R ≥ 0, there exists T R ≥ 0 such that, for each ε ≤ ε * and for each solution (x m (t), e v (t)) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] starting from (x mo , e vo ) ∈ R n × R n(N -1) , we have

|(x mo , e vo )| ≤ R =⇒ |(x m (t), e v (t))| ≤ r ∀t ≥ T R .
• We recall from ( 23)-( 24) the existence of a constant d r > 0 such that, for each (x m , e v ) ∈ B r , we have

|G e (0, e v )| ≤ d r |e v | , |G e (x m , e v ) -G e (0, e v )| ≤ d r |x m | . (61)
• Using Item (iv) in Lemma 2, we conclude that the system ėv = -(Λ ⊗ I n ) e v is exponentially stable and there exists P ∈ R (N -1)×(N -1) symmetric and positive definite such that

P Λ + Λ P ≤ -I N -1 .
As a result, combining the aforementioned four properties, we conclude that the time derivative of the Lyapunov function candidate V e (e v ) := e v (P ⊗ I n )e v , under (61) and along the solution (x m (t), e v (t)) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]

λ min (Q) = d r λ max (P ) 2 rε [λ min (Q) -(1 + r)d r λ max (P )λ min (Q)ε]
.

Hence, by introducing the class K function

ε v (ρ v ) := min ε * , 1 (1 + r)d r λ max (P ) × ρ v drλmax(P ) 2 r (1+r)drλmax(P )λmin(Q) + ρ v ,
we conclude that, for each ρ v > 0 and for each ε ≤ ε v (ρ v ), each solution (x m (t), e v (t)) to ( 18) satisfies (60), which proves ( ).

To complete the proof of ( ), we rewrite the dynamics of x m in (18b) as in [START_REF] Martin | Time scale modeling for consensus in sparse directed networks with time-varying topologies[END_REF]; namely, we have

ẋm = F m (x m ) + G m (x m , e v ),
where

G m (x m , e v ) := F m (x m , e v ) -F m (x m ).
Next, we recall from (15) the existence of h :

R n × R n(N -1) → R ≥0 continuous such that G m (x m , e v ) ≤ h(x m , e v )|e v | ∀(x m , e v ) ∈ R n × R n(N -1) .
Hence, for each ρ v > 0, for each ε ≤ ε v (ρ v ), and for each solution (x m (t), e v (t)) to ( 18) starting from B R , under the global ultimate boundedness of ( 18) and (60), we conclude the existence of

T R1 > T R > 0 and c r > 0 such that |G m (x m (t), e v (t))| ≤ c r ρ v ∀t ≥ T R1 .
Next, we use Assumption 4 to show that, according to the notation in Appendix II.C, the decomposition ω := γ o ∪ γ 1 has no cycles and we can only have γ o < γ 1 . Indeed, since the periodic orbit γ o is asymptotically orbitally stable, we conclude that R(γ o ) = γ o ; hence, R(γ o ) ∩ U (γ 1 ) = γ o ∩ U (γ 1 ) and, since γ o is invariant, compact, and γ o ∩ γ 1 = ∅, it follows that R(γ o ) ∩ U (γ 1 ) = ∅. Thus, we cannot have γ 1 < γ o . The latter also implies that we cannot have a 2-cycle. Next, we exclude the possibility of having a 1-cycle using asymptotic orbital stability of γ o and by definition of the set γ 1 . Finally, the compact disconnected subset ω := γ o ∪ {0} is globally attractive for the unperturbed dynamics [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF]. Hence, using Lemma 11 in Appendix I, we conclude the existence of a continuously differentiable Lyapunov function V : R n → R ≥0 , class K ∞ functions α, ᾱ, α, and a positive constant c ≥ 0, such that

α (|x m | ω ) ≤ V (x m ) ≤ ᾱ (|x m | ω + c) and ∂V ∂x m (x m )F m (x m ) ≤ -α (|x m | ω ) ∀x m ∈ R n .
Hence, along the solution (x m (t), e v (t)) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] starting from B R , the time-derivative of the Lyapunov function V satisfies the following upper-bound

V (x m (t)) ≤ -α (|x m (t)| ω ) + ∂V ∂x m (x m (t)) c r ρ v ∀t ≥ T R1 .
Next, using the global ultimate boundedness of ( 18) and the continuity of ∂V ∂x m

, we conclude the existence of b r and

T R2 > T R > 0 such that ∂V ∂x m (x m (t)) ≤ b r ∀t ≥ T R2 .
As a result, for all t ≥ sup {T R1 , T R2 }, we have

V (x m (t)) ≤ -α α -(V (x m (t))) + b r c r ρ v .
The last step invokes the comparison lemma [30, Lemma 2.5], that establishes the existence of µ

∈ K ∞ such that lim t→∞ V (x m (t)) ≤ µ (b r c r ρ v ) ,
and consequently, lim

t→∞ |x m (t)| ω = lim t→∞ |(x m (t), 0)| Ω ≤ α -(µ (b r c r ρ v )) .
Finally, to complete the proof, it is enough to take ε 1 (ρ) := ε v (χ -1 (ρ)), where

χ(ρ v ) := α -(µ (b r c r ρ v )) + ρ v .

Proof of Lemma 6

The proof is based on a direct application of Lemma 10. The first item in Lemma 10 is holds under Assumption 2, the second item holds with e v = h(x m ) = 0, the third item is satisfied since the boundarylayer model of ( 18) is the linear system ėv = -(Λ ⊗ I n )e v with Λ Hurwitz since the graph is connected; see Lemma 2, and the last item in Lemma 10 holds under Assumption 4.

Proof of Lemma 7

Consider a periodic orbit Γ ε ⊂ T ρ , for some ρ ∈ (0, 1], generated by an α ε -periodic solution to (18) denoted by xε (t) := (x mε (t), e vε (t)). Also, we let the α o -periodic solution to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF], denoted by x mo (t), generating the periodic orbit γ o introduced in Assumption 4. We also let α := α ε -α o . Finally, we introduce the following error coordinate x = x -xε .

In the coordinates x, we re-express system (18) as follows

ẋ = A(x ε (t))x + g(x ε (t), x), (62)
where

A(x ε (t)) :=     ∂Fm ∂x m ∂Gm ∂e v ∂Ge ∂x m - 1 ε (Λ ⊗ I n ) + ∂Ge ∂e v     xm = xmε(t) ev = evε(t), g(x ε (t), x) := F (x ε (t) + x) -F (x ε (t)) -A(x ε (t))x.
Note that g(x ε (t), x) is continuously differentiable in x, continuous in xε (t), and we can find κ

∈ K (inde- pendent of ε since Γ ε ⊂ T ρ ⊂ T ρ=1 ) such that |g(x ε (t), x)| ≤ κ(x)|x|. (63)
To prove the lemma, we follow the following steps:

1. We start re-describing (62) using the new time scale τ := αo αε t, which gives us

x := dx dτ = α ε α o A xε α ε α o τ x := A 1 (x 1 (τ )) x + α ε α o g (x 1 (τ ), xε ) ,
where

A 1 (•) := α ε α o A(•) and x1 (τ ) := xε α ε α o τ .
2. Next, we re-express A 1 (x 1 (τ )) as

A 1 (x 1 (τ )) = A (x o (τ )) + ∆ 1 (x 1 (τ ), xo (τ ), ε, ρ, α), (64)
where xo (τ ) := (x mo (τ ), 0) and

∆ 1 (x 1 (τ ), xo (τ ), ε, ρ) := A 1 (x 1 (τ )) -A 1 (x o (τ )) + α α o A(x o (τ )).
Furthermore, we show that lim (ρ,ε, α)→(0,0,0)

|∆ 1 (x 1 (τ ), xo (τ ), ε, ρ, α)| = 0 ∀τ ∈ [0, α o ] (65)
by following two steps.

• We first apply Tikhonov Theorem (see Lemma 9 and Remark 8) on the singularly perturbed system [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], to conclude the existence of ε * > 0 and M > 0 such that, for each ε ∈ [0, ε * ], the solution xε (t) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] and the signal xo (t), with x mo (t) solution to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF], satisfy

|x ε (0) -xo (0)| ≤ ρ =⇒ |x ε (t) -xo (t)| ≤ M ρ ∀t ∈ [0, α ε ]. (66)
• Now, using (66), we conclude that, for each τ ∈ [0, α o ], we have

xε α ε α o τ = xo α ε α o τ + O(ρ) = xo (τ ) -∆ xo (τ, ε, α) + O(ρ), where ∆ xo (τ, ε, α) := xo α ε α o τ - α α o τ -xo α ε α o τ .
Under Item (i) in Lemma 2, there exists M ρ > 0 such that, for each τ ∈ [0, α o ], we have

|∆ xo (τ, ε, α)| ≤ sup 0≤t≤αε | ẋo (t)| |α| ≤ sup 0≤t≤αε F (x o (t)) |α| ≤ M ρ |α| ,
where, for xo (t) := (x m (t), 0), we have

F (x o (t)) := F m (x m (t)) + G m (x m (t), 0) G e (x m (t), 0) .
As a result, since both x1 (τ ) and xo (τ ) are α o -periodic, we conclude that lim (ρ,ε, α)→(0,0,0)

|x 1 (τ ) -xo (τ )| = 0 ∀τ ≥ 0. ( 67
)
As a result, using (67) and Item (i) in Lemma 2, we conclude that (65) holds.

3. Now, we decompose the matrix A(x o (τ )) as

A(x o (τ )) := ∂Fm ∂x m (x mo (τ )) 0 0 0 + 1 ε A ε (x o (τ )),
where

A ε := 0 0 0 -Λ ⊗ I n + ε    0 ∂Gm ∂e v ∂Ge ∂x m ∂Ge ∂e v    xm = xmo(τ ) ev = 0.
Under Lemma 12, we know that A ε (x o (τ )) admits n eigenvalues of the form

λ j (A ε (x o (τ ))) = o(ε) ∀j = {1, 2, ..., n}.
Furthermore, we consider the non-singular matrix T (ε, τ ) transforming A ε (x o (τ )), as in Lemma 13, into the block-diagonal form 1) is the representation of A ε (x o (τ )) on the invariant spectral subspace corresponding to the remaining eigenvalues. Now, for ε > 0 sufficiently small, we use Lemma 14 to conclude that the transformation matrix T (ε, τ ) is continuous. Furthermore, using Lemma 15, we conclude that T (ε, τ ) is analytic (smooth) in ε. As a result, for each τ ∈ [0, α o ), we have lim ε→0 T (ε, τ ) = I nN and lim Hence, we can find ρ * * > 0 and ε * * > 0 such that, for each ε ∈ (0, ε * * ], the origin y = 0 for (69) is uniformly exponentially stable on the set {y ∈ R nN : |y| ≤ ρ * * }.

J ε (τ ) := Λ o (A ε (x o (τ ))) 0 n×n(N -1) 0 n(N -1)×n Λ 1 (A ε (x o (τ ))) , where Λ o (A ε (x o (τ ))) ∈ R n×n is the representation of A ε (x o (τ )) on the invariant spectral subspace corresponding to {λ j (A ε (x o (τ )))} n j=1 . Similarly, Λ 1 (A ε (x o (τ ))) ∈ R n(N -1)×n(N -

Proof of Lemma 8

To establish the proof, we start analyzing the linearization of (18) around the origin, which is given by ẋ = (A ε /ε)x, where Furthermore, we consider the non-singular matrix T (ε) transforming A ε , as in Lemma 13, into the blockdiagonal form

A ε :=          0 0 0 -(Λ ⊗ I n )   + ε     ∂Fm ∂x m (0)
J ε :=   Λ u (A ε ) 0 0 Λ s (A ε ) 0 n×n(N -1) 0 n(N -1)×n Λ 1 (A ε )   ,
where Λ u (A ε ) ∈ R k×k is the representation of A ε on the spectral subspace corresponding to the unstable eigenvalues in {λ j (A ε )} n j=1 , Λ s (A ε ) ∈ R (n-k)×(n-k) is the representation of A ε on the spectral subspace corresponding to the stable eigenvalues in {λ j (A ε )} n j=1 . Similarly, Λ 1 (A ε ) ∈ R n(N -1)×n(N -1) is the representation of A ε on the spectral subspace corresponding to the remaining eigenvalues of A ε .

For ε > 0 sufficiently small, we use Lemma 15 to conclude that T (ε) is analytic (smooth) in ε. As a result, we have is a representation ∂Fm ∂x m (0) on its stable subspace.

Λ u (A ε ) = ε ∂F m ∂x m (0) u + o(ε), Λ s (A ε ) = ε ∂F m ∂x m (0) s + o(ε), Λ 1 (A ε ) = -(Λ ⊗ I n ) + O(ε).
According to the aforementioned properties, we conclude the existence of ε > 0 such that, for each ε ∈ (0, ε ], the origin is a hyperbolic equilibrium for [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. On the other hand, following the notation of Lemma 16, we conclude that (71)

Moreover, the forward and the backward overshoots of A ε /ε satisfy

c s (A ε /ε) = c s (A ε ) = c s (J ε ) = c s (blkdiag{Λ s (A ε ), Λ 1 (A ε )}) , c u (A ε /ε) = c u (A ε ) = c u (J ε ) = c u (Λ u (A ε )). (72)
Using (70), we conclude that for ε * small enough, the forward and backward overshoots (c u (Λ u (A ε )), c s (Λ s (A ε ))) can be chosen as

As a result, under (71), (72), and (73), we conclude that, for ε > 0 sufficiently small, we can find γ > 0 and ∆ > 0 such that, for each ε ∈ (0, ε ], we have

c s (A ε /ε)γ + κ(γ) c s (A ε /ε) r s (A ε /ε) -µ(A ε /ε) + c u (A ε /ε) r u (A ε /ε) + µ(A ε /ε) ∆ ≤ ∆.
Hence, using Lemma 16, Items (i) and (ii) in Lemma 8 follow with U := B γ and r := γ/2.

To prove Item (iii), we apply Tikhonov Theorem-see Lemma 9 and Remark 8-to conclude that, given T > 0, we can find ρ * > 0 and ε * > 0 sufficiently small and M > 0 such that, for each ε ∈ [0, ε * ], the solution x(t) to [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] and the signal xo (t) := (x mo (t), 0), with x mo (t) solution to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] Now, we pick T > 0 as the largest time a solution x m (t) to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] starting from Γ 1 \B r takes to enter the ball B r/2 . Such a T > 0 always exists by definition of the set Γ 1 and is finite since the set Γ 1 is compact. By taking ρ * = r 3M , we conclude that |x(T )| must be inside B r .

•

  Lemma 5 establishes global asymptotic practical stability of the {(x m , e v ) ∈ R nN : x m ∈ γ o ∪ γ 1 and e v = 0}; 1 i.e. the eigenvalues of A := ∂fm(xm) ∂xm xm=0 has 0 < k < n eigenvalues with positive real part and n -k with negative real part.

  where r s (A) := {min | (λ j (A))| : (λ j (A)) < 0} r u (A) := {min | (λ j (A))| : (λ j (A)) > 0}.

ε→0 1 εT 1 ε Λ 1 ( 1 εFinally 5 ..

 11115 Λ o (A ε (x o (τ ))) = 0, lim ε→0 [Λ 1 (A ε (x o (τ ))) + (Λ ⊗ I n )] = 0.(68) 4. At this point, we introduce the change of coordinates xε := T (ε, τ )x, where T is a continuously differentiable non-singular approximation of T chosen according to Lemma 17. In the new coordinates,xε + T A ε (x o (τ )) Txε + T ∆ 1 T -+ Ṫ Txε + α ε α o T g x1 (τ ), Txε .The latter can be further expressed asx ε = ∂Fm ∂x m (x mo (τ )-I nN ]x ε + T ∆ 1 T -+ Ṫ Txε + [ T -T ]A(x o (τ ))[ T --T -]x ε + α ε α o T g x1 (τ ), Txε , A ε (x o (τ ))) xε + Λ o (A ε (x o (τ )-I nN ]x ε + T ∆ 1 T -+ Ṫ Txε + [ T -T ]A(x o (τ ))[ T --T -]x ε + α ε α o T g x1 (τ ), Txε . -I nN ] + T ∆ 1 T -+ Ṫ T -+ [ T -T ]A(x o (τ ))[ T --T -], h ε (τ, xε ) := α ε α o T g x1 (τ ), Txε .As a last step, applying Floquet Theory, we conclude the existence of a non-singular α o -periodic matrix P o : R ≥0 → R n×n and B o ∈ R n×n , such that Ṗo (τ ) + P o (τ ) -∂F m (x mo (τ )) ∂x m P o (τ ) = B o . Now, using the change of coordinates y := P(τ ) -xε , with P(τ ) := blkdiag P o (τ ), I n(N -1) ,we obtainy = B(ε, τ )y + Ā3 (τ, ρ, ε, α)y + P(τ ) -h ε (τ, P(τ )y) ,(69)where Ā3 (τ, ρ, ε, α) := P(τ ) -A 3 (τ, ρ, ε, α)P(τ ),B(ε, τ ) := B o 0 0 (Λ⊗In)-[Λ1(Aε(xo(τ )))+(Λ⊗In)] -εNote that Ā3 is continuous in its arguments and, using (65), (68), and Lemma 17, we conclude that lim(ρ,ε, α)→0 | Ā3 (τ, ε, ρ, α)| = 0 ∀τ ∈ [0, α o ].Furthermore, for the block-diagonal matrix B(ε, τ ), we use Assumption 4 and Lemma 3, to conclude that the upper block B o is Hurwitz. Moreover, the lower block (Λ⊗In)-[Λ1(Aε(xo(τ )))+(Λ⊗In)] -ε, for sufficiently small ε > 0, has all its characteristic multipliers inside the unit circle since Λ is Hurwitz andlim ε→0 [Λ 1 (A ε (x o (τ ))) + (Λ ⊗ I n )] = 0 ∀τ ∈ [0, α o ].Finally, h ε (τ, xε ) is continuous in τ and continuously differentiable in xε , and under (63), we can find κ h ∈ K such that P(τ ) -h ε (τ, P(τ )y) ≤ κ h (|y|)|y| ∀τ ∈ [0, α o ].

.

  Using Lemma 12, we know that A ε admits n eigenvalues of the formλ j (A ε ) = λ j ∂F m ∂x m (0) ε + o(ε) ∀j = {1,2, ..., n}.

∈

  R (n-k)×(n-k)

  r s (A ε /ε) := {min | (λ j (A ε /ε))| : (λ j (A ε /ε)) < 0} = r s ∂F m ∂x m (0) + O(ε), r u (A ε /ε) := {min | (λ j (A ε /ε))| : (λ j (A ε /ε)) > 0} = r u ∂F m ∂x m (0) + O(ε).

c

  s (Λ s (A ε )) = c s ∂F m ∂x m (0) s + c s (Λ) , c u (Λ u (A ε )) = c u ∂F m ∂x m (0) u ∀ε ∈ (0, ε * ]. Now, if we let g(x) := F (x) -(A ε /ε) x,we conclude that g does not depend on ε and we can find κ ∈ K such that |g(x)| ≤ κ(|x|)|x|.

  , satisfy |x(0) -xo (0)| ≤ ρ * =⇒ |x(t) -xo (t)| ≤ M ρ * ∀t ∈ [0, T ].

  starting from (x mo , e vo ) ∈ B R , satisfies Ve (e v (t)) ≤ d r λ max (P )r |e v (t)| -1 ε -d r λ max (P ) |e v (t)|

				2
			∀t ≥ T R .
	Now, for ε ≤	1 drλmax(P ) , we obtain
		Ve (e v (t)) ≤ d r λ max (P )r -	[1 -(1 + r)d r λ max (P )ε] λ max (P )ε	V e (e v (t))
			∀t ≥ T R .
	As a result,
		lim t→+∞	|e v (t)| 2 ≤ lim

t→+∞

V e (e v (t))

This work was supported by the French ANR via project HANDY, contract number ANR-18-CE40-0010.

Appendix I: Background Distortion of the solutions to singularly perturbed systems

We recall a version of the well-known Tikhonov Theorem [START_REF] Khalil | Nonlinear systems[END_REF]Theorem 9.1], originally published in [START_REF] Tikhonov | Systems of differential equations containing small parameters in the derivatives[END_REF], which allows to approximate the trajectories of a singularly perturbed system by the trajectories of its unperturbed form, along a finite interval of time.

Lemma 9 (Tikhonov Theorem) Consider the singularly perturbed system ż = f (z, e, ε)

Assume that there exist t 1 , r, ρ, and ε o > 0 such that, for each (z, e, ε)

we have

• the functions f , g, and their first partial derivatives with respect to (z, e, ε) are continuous.

• the function e = h(z), solution to the equation 0 = g(z, e, 0), and the Jacobian matrix ∂g(z,e,0) ∂e have continuous first partial derivatives with respect to their arguments.

Input-to-state stability of compact disconnected subsets

In this section, we recall a version of [START_REF] Angeli | Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF]Theorem 1] on the characterization of input-to-state stability of compact disconnected subsets for nonlinear systems of the form

where D ⊂ R m d is a closed subset and f : R n × D → R n is continuously differentiable. Consider a compact subset ω ⊂ R n that is invariant for the unperturbed system

The attracting and repulsing subsets of ω are given by

Furthermore, we consider the following assumption.

Assumption 6

The set ω is the union of disconnected and nonempty subsets

Definition 5 We say that ω i < ω j , for some i, j ∈ {1, 2, ..., k}, if U (ω i ) ∩ R(ω j ) = ∅, which implies the existence of a solution to [START_REF] Matthews | Phase diagram for the collective behavior of limit-cycle oscillators[END_REF] relating ω i to ω j .

Next, we introduce the following properties:

• A filtration ordering is a numbering of the ω i s so that ω i < ω j implies i ≤ j.

The result in [51, Theorem 1] is based on the following assumption:

Assumption 7 The decomposition of ω into ω := k i=1 ω i has no cycles and forms a filtration ordering of ω.

Lemma 11 (Angeli-Effimov Theorem) Consider the nonlinear system in (47) and let ω be a compact subset that is invariant for the unperturbed system ẋ = f (x, 0) such that Assumptions 6-7 hold. Then, the following two properties are equivalent.

• The system enjoys the asymptotic gain property; namely, there exists a class K ∞ function η such that each solution x to (47) satisfies lim sup

• The system admits an input-to-state Lyapunov function; namely, there exists a C 1 function V : R n → R ≥0 , class K ∞ functions α, ᾱ, α, γ, and a positive constant c ≥ 0, such that

As a result, (55) holds. Furthermore, we note that

which implies that (56) also holds. Finally, (57) holds under (55) and the continuous differentiability of T .

Consider the dynamical system

where f : R n → R n is continuously differentiable and the origin x = 0 is a hyperbolic equilibrium point.

The following lemma allows us to show that the propagation of the local stable manifold W s (0) using the backward solutions to (58) is a null-measure set.

Lemma 18 Consider the dynamical system (58) and let S o ⊂ R n and T > 0 such that, for each x o ∈ S o , the solution x(t) is well defined on [0, T ]. Furthermore, for each t ∈ [0, T ], we define the reachable set

Appendix III: Proofs

Proof of Lemma 4

By Assumption 4, γ o is globally orbitally asymptotically stable. Also, ω is globally attractive. Since γ o and γ 1 are disconnected, the solutions must converge to either one. Then, to show the statement of the Lemma, we show that the solutions to [START_REF] Tumash | Synchronization patterns in Stuart-Landau networks: a reduced system approach[END_REF] converging to γ 1 must start from a null measure set which is the global stable manifold W s o (0). Indeed, to find a contradiction, we let x mo / ∈ W s o (0) such that the solution x m (t) to (17) starting from x mo converges γ 1 . Now, using Lemma 16, we conclude the existence of r > 0 and t 1 > 0 such that

As a result, the solution x m (t) must converge to γ 1 \B r . Now, we let a strictly increasing sequence of times {t 1 , t