Master Stability Function: the formalism that separates dynamics from structure in networks of coupled systems. Theory, applications, and a generalization fordiffusive coupling

Abstract

The goal of this paper is to introduce the formalism of the Pecora and Carroll’s MSF (acronym of Master Stability Function), emphasizing its educational potential for physics students (and related subjects) and the meticulous reconstruction of the theory, accompanied by some classic applications. We begin with a detailed description of the classical MSF formalism as a tool for assessing synchronizability in dynamical networks. Next, we apply the formalism to three typical nonlinear dynamical systems under different coupling schemes: the Lorenz system (strange attractor), the Duffing forced oscillator (forced nonlinearity), and the Van der Pol forced oscillator (relaxational oscillator). The choice of these systems is due to their historical relevance and the contribution of their authors to nonlinear dynamics and the study of chaos. Based on the numerical results, we discuss the behavior of MSFs in general coupled nonlinear dynamical systems, not limited to phase-reduced or weakly coupled models, as already reported in the literature. Finally, we explore a generalization of the classical formalism where the specific algebraic form of the diffusive coupling is relaxed. This broadens the scope of MSF formalism, while maintaining its analytical structure, allowing us to formally treat cases such as the Kuramoto model and other systems that do not fit into the traditional formulation, thereby highlighting the versatility and elegance of the formalism. This paper is primarily intended as an educational contribution, while also offering a concise methodological synthesis and a formal extension for general diffusive coupling.

Keywords:
Master Stability Function; Synchronizability; Dynamical networks; Nonlinear systems.

1. By Way of Introduction

Synchronization is an emerging collective phenomenon: it is not found in any isolated element, but arises spontaneously from the interaction among them. It is present in various areas of science, from physics to biology, chemistry, engineering and social life [¹, ²]. It is something we see in countless contexts, from metronomes on an oscillating platform¹[³, ⁴], that synchronize their frequencies and phases even if they start out in a disorderly fashion, to the collective glow of male fireflies [⁵, ⁶, ⁷], where individuals light up randomly and, spontaneously, the synchronization emerges as a collective phenomenon. There is no one “conducting the orchestra” or anything like that, but the flash of each individual stimulates the flashing of others, and suddenly patterns of synchronization begin to form. There are other examples, including pendulums that oscillate in unison [⁸], the coherence of neurons that fire collectively [⁹, ¹⁰], or chaotic systems [¹, ², ¹¹] that come to share the same trajectory on the attractor. All these examples illustrate the surprising or remarkable tendency of certain dynamical systems to synchronize their behavior.

From a theoretical point of view, however, understanding when and how this synchronization occurs in networks of $N$ nonlinear dynamical systems, each described by $n$ equations, is a considerably complicated task. In principle, the problem involves dealing with a set of $nN$ ordinary differential equations describing the dynamics of the coupled systems. The explicit form of these equations that will be considered in this paper is introduced in Section 2, where the full model, notation and underlying assumptions are discussed in detail. In most cases, their analytical integration cannot be done using elementary methods [¹²]. In such situations, one must resort to numerical methods, which, in many scenarios, can also become impractical, even with the technological advances and computing power currently available.

It was in this context that Louis M. Pecora and Thomas L. Carroll introduced the formalism of the Master Stability Function (MSF) [¹³, ¹⁴]. This formalism is one of the most powerful tools for studying synchronization in highly complex dynamical systems (e.g., chaotic systems) when they are coupled through pairs of nodes in a network [²]. The great genius of this method lies in the fact that it reduces a high-dimensional problem (i.e., the stability of synchronization in coupled dynamical networks) to a single parametric linear equation, whose explicit form is given in equation (12) and will be derived step by step in Subsection 2.3. This equation elegantly separates the influence of the individual dynamics from that of the network structure. Thus, without needing to know the network topology in detail, it is possible to determine, based solely on local dynamics and the coupling scheme, which networks can or cannot achieve a stable synchronized state. In fact, despite its wide use in research, the core logic of the MSF can be presented in a remarkably self-contained way, which makes it especially suitable for this educational exposition.

Beyond the purely mathematical formulation, it is useful to keep in mind a simple physical intuition of what is happening in a network of coupled dynamical systems. Synchronization results from the interplay between the intrinsic dynamics of each node and the mutual coupling between them, which usually has a diffusive character [², ¹⁵, ¹⁶]. Each individual system tends to follow its own intrinsic dynamics, which may be periodic or chaotic, while the diffusive coupling term acts as an interaction that continuously compares the states of different nodes and tries to reduce their differences [¹³]. Synchronization emerges precisely from the competition between these two effects [¹⁷, ¹⁸]: when the interaction is sufficiently strong, the coupling is able to suppress the growth of discrepancies between dynamical systems and promotes a synchronized collective evolution. In this sense, the MSF formalism tells us the conditions under which such differences are attenuated rather than amplified, thereby indicating when a synchronous state may be expected to be physically possible.

From a methodological standpoint, the MSF provides the largest Lyapunov exponent associated with the directions transverse to the synchronization subspace [¹³, ¹⁹]. The necessary condition for the synchronous state (or synchronization manifold) to be linearly stable is precisely that this exponent is negative [¹³].

The purpose of this paper is to provide an accessible, self-contained and pedagogical introduction to this formalism to physics students and related disciplines, reconstructing step by step the theoretical basis of the classical MSF paradigm and highlighting its conceptual elegance through the analysis of synchronization in typical chaotic dynamical systems. To this end, we have chosen three classical nonlinear systems in the history of dynamics and chaos: the Lorenz system, the forced Duffing oscillator, and the Van der Pol oscillator. Through these examples, we will investigate the form of the MSF for each of them, under different linear coupling schemes between their components, and discuss the different qualitative scenarios that the MSF can exhibit, which should be understood as generic phenomenological profiles reported in the literature for a broader class of coupled nonlinear systems [²⁰].

Finally, as a methodological extension, we explore a natural generalization of the classical formalism, in which we consider a generalized form of diffusive coupling. This extension preserves the analytical structure of the theory but broadens its scope, allowing for the incorporation of more general interactions and thus making it applicable to cases where the coupling does not fit the specific form assumed in the traditional formulation (e.g., the Kuramoto model [²¹, ²²]). We hope that this will not only provide a robust tool for studying synchronization, but also highlight the beauty and versatility of a formalism that continues to inspire research in the area of dynamical systems. Our goal is therefore to guide the reader through the MSF framework in a self-contained way, and then show how its logic naturally extends beyond the classical setting in a form that remains simple and usable, without departing from its original spirit and emphasizing understanding over novelty. In the conclusions we summarize the results and limitations of the MSF method. We also comment on recent extensions of the method not covered by this paper.

2. Determining the Synchronizabilityof a Network Through the MSF Approach

Consider a network of $N$ identical coupled dynamical systems. Each dynamical system is described by a state vector $x_i\in {ℝ}^n$, with $i=\mathrm{1,\; 2},\dots ,N$. Then, in the absence of coupling, each system evolves according to the same differential equation,

(1) $${\dot{x}}_i=f\left(x_i\right),$$

where $f:{ℝ}^n\to {ℝ}^n$. For simplicity, our derivation, will consider only continuous-time systems. Nevertheless, the MSF approach also applies to discrete-time systems [²³], including chaotic maps [²⁴].

When coupling is introduced, the systems interact through a given coupling function, $h:{ℝ}^n\to {ℝ}^n$ such that,

$${\dot{x}}_i=f\left(x_i\right)+\sigma \sum _{j=1}^N{a}_{ij}{\omega }_{ij}\left(h\left(x_j\right)-h\left(x_i\right)\right),$$

where $\sigma \in ℝ$ is called the coupling strength, ${\omega }_{ij}>0$ are the connection weights, and $a_{ij}$ are the entries of the adjacency matrix²$A=\left(a_{ij}\right)$, given by:

$$a_{ij}=\{\begin{array}{cc}1\mathit{\text{if nodes i and j are connected}},\hfill & \\ 0\mathit{\text{otherwise}},\hfill & \end{array}$$

such that $\text{dim}A=N\times N$.

In Appendix A we provide a quick and concise review of some basic concepts from graph theory. This review is not just included for completeness: throughout this section we will make extensive use of those notions. Therefore, the appendix provides the minimum background needed to follow the developments made here without constantly referring back to specialized literature.

Although the purpose of this paper is to present the MSF formalism to students of physics (and related areas) with little experience in dynamical systems and network theory, we assume that the reader has a basic knowledge of linear algebra and ordinary differential equations. Standard textbooks suitable for physics students include, e.g., [²⁵, ²⁶, ²⁷, ²⁸] for linear algebra and [²⁹, ³⁰, ³¹, ³²] for ordinary differential equations and nonlinear dynamics.

2.1. The dynamical model

2.1.1. Assumptions

• The systems that make up the network are identical. This is an ad hoc assumption of the formalism.

• The dynamics of the $i$-th system is governed by the decoupled vector field $f\left(x_i\right)$ and by the diffusive coupling term, i.e., $h\left(x_j\right)-h\left(x_i\right)$, for every connected pair of nodes $i$ and $j$.

• Simple network: there are no self-connections, i.e., $a_{ii}=0$ for all $i$.

• Undirected network: the adjacency matrix is symmetric, $a_{ij}=a_{ji}$.

• Unweighted network: uniform coupling strength for all connections, i.e., ${\omega }_{ij}=1$.

• Weak coupling regime: $\sigma$ sufficiently small so that coupling does not dominate the intrinsic dynamics [³³].

2.1.2. The model

The coupled dynamics of the network is governed by the following system of $n\times N$ ordinary differential equations:

(2) $${\dot{x}}_i=f\left(x_i\right)+\sigma \sum _{j=1}^N{a}_{ij}\left(h\left(x_j\right)-h\left(x_i\right)\right),$$

where $i,j=\mathrm{1,\; 2},\dots ,N$, $x_i\in {ℝ}^n$, $f:{ℝ}^n\to {ℝ}^n$ and $h:{ℝ}^n\to {ℝ}^n$. Throughout this work, we assume that the vector field $f$ and the coupling function $h$ are locally Lipschitz continuous (and, in practice, at least continuously differentiable $C^1$), which ensures the local existence and uniqueness of solutions to the associated initial value problems³. An accessible presentation of these existence and uniqueness conditions, oriented to physics students can be found in the book Nonlinear Dynamics and Chaos by Steven Strogatz [³¹], see § 2.5.

2.1.3. The Laplacian matrix

The combinatorial Laplacian matrix, $L=\left(l_{ij}\right)$, of dimension $N\times N$, is defined by:

$$l_{ij}=\{\begin{array}{cc}{k}_i\mathit{\text{if i = j}},\hfill & \\ -1\mathit{\text{if nodes i and j are connected}},\hfill & \\ 0\mathit{\text{otherwise}},\hfill & \end{array}$$

where

$$k_i=\sum _{j=1}^N{a}_{ij}$$

is the degree of the $i$-th node in the network. Equivalently, in matrix form one can write $L=D-A$, where $D=\text{diag}(k_1,k_2,\dots ,k_N)$ is the degree matrix. Thus, an alternative expression is

$$l_{ij}={\delta }_{ij}{k}_i-a_{ij},$$

where ${\delta }_{ij}$ denotes the Kronecker delta.

Therefore, due to the diffusive nature of the interaction, it is possible to write the coupling of equation (2) in terms of the Laplacian matrix. Indeed,

$$\begin{array}{cc}\hfill {\dot{x}}_i& =f\left(x_i\right)+\sigma \sum _{j=1}^N{a}_{ij}\left(h\left(x_j\right)-h\left(x_i\right)\right)\hfill \\ & =f\left(x_i\right)+\sigma \left\{\sum _{j=1}^N{a}_{ij}h\left(x_j\right)-\left(\sum _{j=1}^N{a}_{ij}\right)h\left(x_i\right)\right\}\hfill \\ & =f\left(x_i\right)+\sigma \left(\sum _{j=1}^N{a}_{ij}h\left(x_j\right)-k_{i}h\left(x_i\right)\right)\hfill \\ & =f\left(x_i\right)+\sigma \left(\sum _{j=1}^N{a}_{ij}h\left(x_j\right)-\sum _{j=1}^N{\delta }_{ij}{k}_{i}h\left(x_j\right)\right)\hfill \\ & =f\left(x_i\right)+\sigma \sum _{j=1}^N\left(a_{ij}-{\delta }_{ij}{k}_i\right)h\left(x_j\right)\hfill \\ & =f\left(x_i\right)-\sigma \sum _{j=1}^N\left({\delta }_{ij}{k}_i-a_{ij}\right)h\left(x_j\right)\hfill \\ & =f\left(x_i\right)-\sigma \sum _{j=1}^N{l}_{ij}h\left(x_j\right).\hfill \end{array}$$

2.1.4. Diffusive form

We can rewrite the model in terms of the Laplacian matrix,

(3) $${\dot{x}}_i=f\left(x_i\right)-\sigma \sum _{j=1}^N{l}_{ij}h\left(x_j\right),$$

which makes explicit the diffusive nature of the coupling mentioned earlier. Moreover, the coupling function $h\left(x_j\right)$ can be reinterpreted as an output function, which generates an output from the state $x_j$ and sends it to the other systems in the network. Equation (3) corresponds to the standard configuration used in the literature to describe the dynamics of a network of $N$ coupled identical systems under diffusive interactions [², ¹⁴, ²⁰, ³⁸, ³⁹, ²³, ⁴⁰].

It is worth emphasizing that, although equation (2) is written explicitly as a difference $h\left(x_j\right)-h\left(x_i\right)$ and the adjacency matrix satisfies $a_{ii}=0$ according to the assumptions of the model (i.e., there are no self-connections), the equivalent Laplacian form in equation (3) naturally contains diagonal terms $l_{ii}=k_i\ne 0$. These diagonal contributions do not represent self-interactions, but arise solely from the algebraic reorganization of the same diffusive coupling when expressed in Laplacian form.

2.2. The synchronization manifold

2.2.1. Existence

Note that the Laplacian is a symmetric and zero-row sum matrix, i.e., for a fixed $i$,

$$\begin{array}{l}\sum _{j=1}^N{l}_{ij}=\sum _{j=1}^N{l}_{ji}=\sum _{j=1}^N\left(a_{ij}-{\delta }_{ij}{k}_j\right)\\ =\sum _{j=1}^N{a}_{ij}-k_i=k_i-k_i=0,\end{array}$$

i.e.,

$$\sum _{j=1}^N\left\{\begin{array}{c}\mathit{\text{entries}}\\ \mathit{\text{from the i-th}}\\ \mathit{\text{row or column}}\end{array}\right\}=0.$$

This implies that, if $s\left(t\right)$ is a solution of equation (3) for $\sigma =0$, then the constraint

(4) $$s=x_1=x_2=\cdots =x_N$$

defines a solution of equation (3) for any value of $\sigma$, for any coupling function $h$ and regardless of the network topology [⁴¹], i.e., there exists a completely synchronized state in which all systems in the network behave in the same way.

Actually, given the form of the coupling in equation (2), which is the same as equation (3), it becomes clearer that the state $x_i\left(t\right)=s\left(t\right)$ (in which all systems evolve synchronously following the solution of the isolated or decoupled system, i.e., equation (1)), with $\dot{s}=f\left(s\right)$, is always a solution of equation (3). These two perspectives are fully equivalent and simply reflect different but algebraically identical ways of writing the same diffusive coupling.

We emphasize that the equality in equation (4) is understood as a vector equality in ${ℝ}^n$, so that complete synchronization means that all subsystems follow the same trajectory in phase space, with all state variables synchronized componentwise. Moreover, for a given initial condition, the synchronous trajectory $s\left(t\right)$ is uniquely defined, consistently with the assumed local Lipschitz continuity of the vector field.

In this sense, the equation (4) represents a subspace of the vector space of the solutions of the equation (3), in which all dynamical systems behave as an isolated system following the dynamical equation (1). This subspace is called the synchronization manifold [²⁰, ⁴², ⁴³, ¹³].

2.2.2. Underlying variational problem

Let us now study the stability of the synchronized solution $x_i=s$, since it is not enough for the solution to exist, but it also needs to be stable if we want to have a chance to observe it in practice.

Since all $N$ dynamical systems lie in the synchronization subspace given by equation (4), they will evolve synchronously. The central question, then, is: is the state $s$ stable in the presence of small perturbations? In other words, for which values of $\sigma$ is the synchronization manifold $s$ linearly stable? Also, what types of networks “favor” synchronization and, more generally, what is the role of network topology? These are the questions answered by the MSF!

2.2.3. Linearization around $s$

Let us consider

(5) $$x_i=s+\delta x_i,$$

so that, since $s$ evolves according to the isolated system,

$${\dot{x}}_i=\dot{s}+{\dot{\delta x}}_i=f\left(s\right)+{\dot{\delta x}}_i,$$

and, substituting $x_i$ for the equation (3),

$${\dot{\delta x}}_i={\dot{x}}_i-f\left(s\right)=f\left(x_i\right)-\sigma \sum _{j=1}^N{l}_{ij}h\left(x_j\right)-f\left(s\right).$$

Now, we linearize around $s$,

$$\begin{array}{cc}\hfill {\dot{\delta x}}_i& ={f\left(x_i\right)|}_{x_i=s}+{Df\left(x_i\right)|}_{x_i=s}\hfill \\ & \cdot \left(x_i-s\right)+𝒪\left({\Vert x_i-s\Vert }^2\right)\hfill \\ & -\sigma \sum _{j=1}^N{l}_{ij}\{h\left(x_j\right){|}_{x_j=s}+Dh\left(x_j\right){|}_{x_j=s}\hfill \\ & \cdot (x_j-s)+𝒪(\parallel x_j-s{\parallel }^2)\}-f(s)\hfill \\ & \approx \cancel{f\left(s\right)}+Df\left(s\right)\delta x_i\hfill \\ & -\sigma \sum _{j=1}^N{l}_{ij}\left(h\left(s\right)+Dh\left(s\right)\delta x_j\right)-\cancel{f\left(s\right)}\hfill \\ & =Df\left(s\right)\delta x_i-\sigma \left(\sum _{j=1}^N{l}_{ij}h\left(s\right)+\sum _{j=1}^N{l}_{ij}Dh\left(s\right)\delta x_j\right),\hfill \end{array}$$

where $Df\left(s\right)$ and $Dh\left(s\right)$ are the Jacobian matrices of the functions $f$ and $h$ evaluated at $s\left(t\right)$, respectively. The assumption that $f$ and $h$ are locally Lipschitz continuous (and, in practice, at least $C^1$) ensures that these Jacobians are well defined along the synchronous trajectory and that the linearization procedure is valid. In the approximation above, we retain only the terms that are linear in the perturbations and neglect terms of quadratic or higher order, i.e., $𝒪\left({\Vert \delta x_i\Vert }^2\right)$.

Moreover, recalling that the Laplacian matrix has zero-row sum property, i.e., ${\sum }_{j=1}^N{l}_{ij}=0$, we obtain

(6) $${\dot{\delta x}}_i=Df\left(s\right)\delta x_i-\sigma \sum _{j=1}^N{l}_{ij}Dh\left(s\right)\delta x_j.$$

Here, we require that all perturbations $\delta x_i$ tend to zero, since the asymptotic behavior $\delta x_i\to 0$ guarantees that the solution $s$ is stable: small perturbations do not grow but instead dissipate along the system’s evolution [¹³].

2.3. Linear stability analysis of the synchronized solution

We now proceed to the manipulation of the perturbation equation, equation (6), until we isolate the structure that gives rise to the MSF.

2.3.1. Perturbation equations in Kronecker representation

The $N$ coupled equations given by equation (6), each of dimension $n$, can be rewritten compactly using the Kronecker tensor product, denoted by $\otimes$. This allows us to write the equations as a single block, i.e., to express them as a single matrix system of dimension $nN$. In Appendix B, we provide a brief review of tensor products with explicit examples. Readers not familiar with this concept are encouraged to read the appendix at this point, before continuing with the perturbation analysis.

Let $X$ be the $nN$-dimensional column vector obtained by stacking the $N$ subsystem state vectors $x_i\in {ℝ}^n$ into a single column vector, i.e.,

$$\begin{array}{l}X={(x_1^T,x_2^T,\dots ,x_N^T)}^T\\ =(x_1^1,x_1^2,\dots ,x_1^n,x_2^1,x_2^2,\dots ,x_2^n,\dots ,x_N^1,x_N^2,\dots ,x_N^n){}^T,\end{array}$$

where ${(\cdot )}^T$ represents the transpose operation. Note our index convention. Throughout the paper, the subscript $i=\mathrm{1,\; 2},\dots ,N$ labels the node (subsystem), whereas the superscript $k=\mathrm{1,\; 2},\dots ,n$ labels the component of the corresponding state vector $x_i\in {ℝ}^n$.

Similarly, we consider our perturbation vector $\delta X$, also of dimension $nN\times 1$, as

(7) $$\delta X=\left(\begin{array}{c}\delta x_1\\ \delta x_2\\ ⋮\\ \delta x_N\end{array}\right)=(\delta x_1^1,\delta x_1^2,\dots ,\delta x_1^n,\delta x_2^1,\delta x_2^2,\dots ,\delta x_2^n,\dots ,\delta x_N^1,\delta x_N^2,\dots ,\delta x_N^n){}^T,$$

where each $\delta x_i\in {ℝ}^n$ represents a perturbation of the state of node $i$. By stacking all equations of the form given in equation (6), we obtain the compact matrix equation:

(8) $$\dot{\delta X}=\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)-\sigma L\otimes Dh\left(s\right)\right)\delta X.$$

A straightforward way to verify this result is to compute the $i$-th block of dimension $n$ from equation (8) and compare it with equation (6), corresponding to the $i$-th node, as shown step by step in Appendix C.

2.3.2. Decoupling into Laplacian transverse eigenmodes

The main idea in this algebraic step is to project $\delta X$ onto the eigenspace spanned by the eigenvectors of the Laplacian matrix. This projection can be performed in block form, preserving the internal block structure. By doing so, the equation (8) can be block-diagonalized into $N$ eigenmodes, transverse to the synchronization manifold and mutually decoupled.

Recall that, since we are assuming a simple and undirected network, the Laplacian is a symmetric matrix and therefore diagonalizable, with a set of real eigenvalues and eigenvectors that form an orthonormal basis for ${ℝ}^N$. Let us also remember that the Laplacian is zero-row sum, so the smallest eigenvalue is ${\lambda }_1=0$. Also, let us assume that the network is connected, so that there is only one zero eigenvalue, while the remaining ones are strictly positive. Thus, the spectrum can be ordered as $0={\lambda }_1<{\lambda }_2\le \cdots \le {\lambda }_N$, showing that $L$ is positive semidefinite. See Appendix A (Theorem 1), where we describe some important spectral properties of the Laplacian.

Furthermore, the eigenvector (already normalized) associated with eigenvalue ${\lambda }_1=0$ is

$$v_1=\frac{1}{\sqrt{N}}\left(\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right),$$

whose components are all equal, indicating that it points along the direction of the synchronous state. Indeed, $v_1$ determines the synchronization manifold, since this manifold is geometrically represented by the hyperdiagonal in phase space, defined by the constraint $x_1=x_2=\cdots =x_N=s$, which is precisely the subspace where $v_1$ lies. On the other hand, the other eigenvectors $v_i$ of the Laplacian, corresponding to the eigenvalues ${\lambda }_i$ (for $i\ge 2$), form a basis for the tangent space (or transverse space) of the synchronization manifold [²³], i.e., they provide the perpendicular directions to the hyperdiagonal.

Within the MSF framework, we aim to ensure that the synchronization manifold (geometrically, the hyperdiagonal in phase space) is attractive along all transverse directions, thereby guaranteeing synchronization stability.

This idea is illustrated schematically in Figure 1, where the synchronization manifold appears as the hyperdiagonal in phase space, aligned with the eigenvector $v_1$, while the transverse directions (associated with the other eigenvectors of the Laplacian) are shown as directions along which perturbations must decay to ensure stability. For synchronization to be stable, the hyperdiagonal must be attractive in all these transverse directions, ensuring that any perturbation decays over time.

Figure 1
Schematic representation of the synchronization manifold as the hyperdiagonal in phase space, aligned with the eigenvector $v_1$. The transverse directions, associated with the other Laplacian eigenvectors, indicate directions along which the hyperdiagonal must be attractive for synchronization to be stable.

Since the Laplacian eigenvectors form an orthonormal basis of ${ℝ}^N$, the perturbation vector $\delta X$ of dimension $nN$, given in equation (7), can be written as a linear combination of this basis using, in addition, the Kronecker product:

(9) $$\delta X=\sum _{k=1}^N{v}_k\otimes {\xi }_k,$$

where ${\xi }_k\in {ℝ}^n$ represents the projection of the perturbation $\delta X$ onto the direction of the eigenvector $v_k$. Differentiating, we obtain

$$\dot{\delta X}=\sum _{k=1}^N{v}_k\otimes {\dot{\xi }}_k.$$

We calculate the two contributions on the right side of equation (8) separately, by substituting equation (9) into it and making use of the mixed-product property of the Kronecker product⁴,

$$\begin{array}{l}\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)\right)\delta X=\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)\right)\left({\displaystyle \sum _{k=1}^N}{v}_k\otimes {\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)\right)\left(v_k\otimes {\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}\left({\mathrm{𝟙}}_N{v}_k\right)\otimes \left(Df\left(s\right){\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}{v}_k\otimes \left(Df\left(s\right){\xi }_k\right),\end{array}$$

$$\begin{array}{l}\left(L\otimes Dh\left(s\right)\right)\delta X=\left(L\otimes Dh\left(s\right)\right)\left({\displaystyle \sum _{k=1}^N}{v}_k\otimes {\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}\left(L\otimes Dh\left(s\right)\right)\left(v_k\otimes {\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}\underset{{\lambda }_k{v}_k}{\underbrace{\left(L{v}_k\right)}}\otimes \left(Dh\left(s\right){\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}{\lambda }_k{v}_k\otimes \left(Dh\left(s\right){\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}{v}_k\otimes \left({\lambda }_{k}Dh\left(s\right){\xi }_k\right).\end{array}$$

Now, let’s substitute these expressions into equation (8),

$$\begin{array}{l}{\displaystyle \sum _{k=1}^N}{v}_k\otimes {\dot{\xi }}_k=\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)\right)\delta X-\sigma \left(L\otimes Dh\left(s\right)\right)\delta X\\ ={\displaystyle \sum _{k=1}^N}{v}_k\otimes \left(Df\left(s\right){\xi }_k\right)\\ -\sigma {\displaystyle \sum _{k=1}^N}{v}_k\otimes \left({\lambda }_{k}Dh\left(s\right){\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}{v}_k\otimes \left(Df\left(s\right){\xi }_k-\sigma {\lambda }_{k}Dh\left(s\right){\xi }_k\right)\\ ={\displaystyle \sum _{k=1}^N}{v}_k\otimes \left(Df\left(s\right)-\sigma {\lambda }_{k}Dh\left(s\right)\right){\xi }_k.\end{array}$$

Next, we multiply both sides on the left by $v_i^T\otimes {\mathrm{𝟙}}_n$ and, using again the mixed-product property of the Kronecker product together with the orthonormality of the eigenvectors, i.e., $v_i^T{v}_k=⟨v_i|v_k⟩={\delta }_{ik}$, we can decouple the transverse eigenmodes,

(10) $$\begin{array}{cc}& \left(v_i^T\otimes {\mathrm{𝟙}}_n\right)\left(\sum _{k=1}^N{v}_k\otimes {\dot{\xi }}_k\right)\hfill \\ & =\left(v_i^T\otimes 1_n\right)\left(\sum _{k=1}^N{v}_k\otimes \left(Df\left(s\right)-\sigma {\lambda }_{k}Dh\left(s\right)\right){\xi }_k\right)\hfill \\ & \sum _{k=1}^N\left(v_i^T\otimes {\mathrm{𝟙}}_n\right)\left(v_k\otimes {\dot{\xi }}_k\right)\hfill \\ & =\sum _{k=1}^N\left(v_i^T\otimes {\mathrm{𝟙}}_n\right)\left(v_k\otimes \left(Df\left(s\right)-\sigma {\lambda }_{k}Dh\left(s\right)\right){\xi }_k\right)\hfill \\ & \sum _{k=1}^N\left(v_i^T{v}_k\right)\otimes \left({\mathrm{𝟙}}_n{\dot{\xi }}_k\right)\hfill \\ & =\sum _{k=1}^N\left(v_i^T{v}_k\right)\otimes \left({\mathrm{𝟙}}_n\left(Df\left(s\right)-\sigma {\lambda }_{k}Dh\left(s\right)\right){\xi }_k\right)\hfill \\ & \sum _{k=1}^N{\delta }_{ik}\otimes {\dot{\xi }}_k=\sum _{k=1}^N{\delta }_{ik}\otimes \left(Df\left(s\right)-\sigma {\lambda }_{k}Dh\left(s\right)\right){\xi }_k\hfill \\ & \sum _{k=1}^N{\delta }_{ik}{\dot{\xi }}_k=\sum _{k=1}^N{\delta }_{ik}\left(Df\left(s\right)-\sigma {\lambda }_{k}Dh\left(s\right)\right){\xi }_k\hfill \\ & {\dot{\xi }}_i=\left(Df\left(s\right)-\sigma {\lambda }_{i}Dh\left(s\right)\right){\xi }_i.\hfill \end{array}$$

Thus, for each eigenmode ${\xi }_i$ of the perturbation $\delta X$, in the direction of the $i$-th Laplacian eigenvector $v_i$, the evolution depends uniquely on the corresponding eigenvalue ${\lambda }_i$.

2.4. The master stability function

For ${\lambda }_1=0$, the variational equation (10) evolves along the synchronous solution of the dynamics for $i=1$, while for $i>1$ it describes perturbations along the transverse directions. We now define the scalar parameter

(11) $$r=\sigma {\lambda }_i$$

and rewrite the equation (10) as the generic $n$-dimensional linear parametric equation,

(12) $$\dot{\xi }=\left(Df\left(s\right)-rDh\left(s\right)\right)\xi .$$

Therefore, if we know the stability of the trivial solution $\xi =0$ for any reasonable value of $r$, then we can deduce the stability for every perturbation mode, with $r=r_i=\sigma {\lambda }_i$[¹⁴].

In this way, we have reduced the original set of $nN$ nonlinear differential equations given in equation (2) to just one equation that gives us the synchronization manifold, equation (4), solution $x_i=s$ of the equation (3), and a single parametric linear equation (12), from which we can compute the set of $n$ conditional Lyapunov exponents⁵[¹⁹, ¹³] for each value of the parameter given in equation (11).

The parametric behavior of the largest of these exponents, $\Lambda =\Lambda \left(r\right)$, is what we refer to as the Master Stability Function (MSF) [¹⁴, ¹³]. Indeed, formally, the MSF is the function that maps the parameter $r=r_i=\sigma {\lambda }_i$ to the largest Lyapunov exponent $\Lambda$, i.e.,

$$\text{MSF}:r\mapsto \Lambda =\Lambda \left(r\right).$$

In practice, for each fixed value of $r$, one computes the largest Lyapunov exponent associated with the transverse variational equation. The collection of these values as a function of $r$ defines the MSF.

Transverse stability has a very intuitive physical interpretation. The synchronization manifold represents the physical situation in which all dynamical systems evolve identically, as if they were isolated. Meanwhile, transverse directions correspond to differences that cause some dynamical systems to differ from others. The implications summarized in Table 1 can be interpreted in simple dynamical terms: when the MSF is negative, these transverse perturbations die out with time [⁴³, ⁴⁵] and, if two dynamical systems are almost synchronized (i.e., the network is in the vicinity of the synchronization manifold), the coupling prevents them from drifting apart again, effectively acting as a restoring mechanism [⁴⁶] that continuously corrects state differences. In this way, the mathematical condition $\Lambda \left(r\right)<0$ simply expresses that network diffusive interaction is able to counterbalance the natural tendency of chaotic dynamics to amplify small perturbations.

Table 1
Key implications of the MSF.

Key implications of the MSF
✓	If $\Lambda \left(r\right)$ is negative, a small perturbation around the synchronization state will decrease exponentially, so that the synchronous solution is stable, at least when the dynamical systems are initialized in its vicinity (linearly stable).
✗	The synchronous solution is unstable and cannot be physically obtained or observed if $\Lambda \left(r\right)$ is positive because small perturbations of the synchronous state will lead to trajectories that diverge from the state.
☞	For a network of coupled dynamical systems, according to equation (3), a necessary condition for synchronization is that all generic coupling parameters $r_i$, for all $i=\mathrm{2,\; 3},\dots ,N$, fall within an interval on the $r$-axis where $\Lambda \left(r\right)$ is negative.

Recall that, according to our assumptions, the network considered is simple and undirected, and therefore the Laplacian is a symmetric, positive semidefinite matrix with real eigenvalues satisfying $0={\lambda }_1<{\lambda }_2\le \cdots \le {\lambda }_N$. Consequently, the parameter $r=r_i=\sigma {\lambda }_i\ge 0$ is also a non-negative real number.

Furthermore, note that the value of the MSF at $r=0$, i.e., $\Lambda \left(r=0\right)=\Lambda \left(\lambda ={\lambda }_1=0\right)$ can be either zero or positive, depending on whether the decoupled dynamics $\dot{x}=f\left(x\right)$, given in equation (1), are periodic or chaotic, respectively. In this work, we are interested in the synchronization of chaotic systems. Hence, in what follows, we assume that the largest Lyapunov exponent of the decoupled system to be positive, i.e., $\Lambda \left(r=0\right)>0$.

For $r>0$ (i.e., for ${\lambda }_2,{\lambda }_3,\dots ,{\lambda }_N$), there are three possible scenarios [¹²], illustrated in Figure 2. In order to study the stability of the synchronization manifold for given $f\left(x\right)$ and $h\left(x\right)$, we first compute and plot the MSF and then evaluate it at each eigenvalue of the Laplacian (scaled by $\sigma$).

Figure 2
Representative of MSF’s possible types for coupled chaotic systems. Type I (blue) is monotonically increasing, Type II (red) is monotonically decreasing, and Type III (green) exhibits a finite interval where $\Lambda \left(r\right)<0$ allowing stable synchronization for appropriate values of $r$. These curves illustrate different qualitative MSF behaviors that determine the coupling ranges for which the synchronous manifold is stable (i.e., synchronizable).

This approach elegantly separates dynamics from structure⁶: given a dynamical system $f\left(x\right)$ and a coupling function $h\left(x\right)$, we can calculate the MSF, universal and valid for any network topology, and immediately establish the structural conditions (i.e., what should the Laplacian spectrum be like?) that guarantee the stability of the synchronization manifold, i.e., identify spectral constraints on the Laplacian, with the requirement that all $r=\sigma {\lambda }_i$ must lie in the so-called “stability region”, where $\Lambda \left(r\right)<0$.

2.4.1. Type I: $\Lambda =\Lambda \left(r\right)$ is a monotonically increasing function

In this case, for a given dynamical system $f\left(x\right)$ and coupling function $h\left(x\right)$, no network topology can achieve stable synchronization. The synchronization manifold, equation (4), always exists, but it is inherently unstable: regardless of how the nodes are arranged or connected, the MSF remains positive for all values of $r$.

Consequently, for type I MSFs, neither adjusting the coupling strength $\sigma$ nor modifying the Laplacian spectrum can stabilize synchronization [⁴⁷, ⁴¹]. For every eigenmode, the product $r=\sigma {\lambda }_i$ yields a positive Lyapunov exponent, rendering the synchronization manifold transversely unstable.

From a physical perspective, type I MSFs describe coupled systems whose intrinsic dynamics are so sensitive to initial conditions that, regardless of the coupling strength, small perturbations around the synchronization manifold cannot be suppressed. As a result, synchronization is physically unattainable, regardless of network topology or coupling strength.

2.4.2. Type II: $\Lambda =\Lambda \left(r\right)$ is a monotonically decreasing function

In this case, the decreasing function $\Lambda =\Lambda \left(r\right)$ crosses zero at any point $r^{\ast }$, i.e., $\exists r^{\ast }\ge 0$ such that $\Lambda \left(r^{\ast }\right)=0$. Then, if the first nonzero Laplacian eigenvalue (times $\sigma$) is greater than $r^{\ast }$, by the Laplacian’s spectral properties, all the other eigenvalues (also multiplied by $\sigma$) will also exceed $r^{\ast }$. Consequently, the MSF is negative for all transverse modes and, therefore, the synchronization manifold becomes stable.

This MSF type characterizes systems that are always synchronizable, provided the network is configured so that $\sigma {\lambda }_2>r^{\ast }$. In other words, with a suitable choice of network topology and/or coupling parameter $\sigma$, synchronization can be achieved regardless of the details of the node dynamics or coupling function.

First, note that the analytical condition for synchronization can be written as $\sigma >r^{\ast }/{\lambda }_2$. Thus, the coupling strength must exceed this critical threshold for synchronization to occur. Once $f\left(x\right)$ and $h\left(x\right)$ (and therefore $r^{\ast }$) are fixed, and given a coupling parameter, the synchronization threshold in this type of MSF is determined solely by the second-smallest eigenvalue (${\lambda }_2$, for connected graphs) of the Laplacian [⁴⁸].

Conversely, if we want synchronization at a lower value of $\sigma$, then ${\lambda }_2$ must be increased⁷. Networks with larger ${\lambda }_2$ values therefore allow synchronization with weaker coupling, meaning that a higher ${\lambda }_2$ reduces the coupling effort needed to achieve synchrony. Consequently, for type II MSFs, larger values of ${\lambda }_2$ generally imply better synchronizability [⁴⁹].

This shows that, for type II MSFs, once the node dynamics $f\left(x\right)$ and coupling function $h\left(x\right)$ are fixed, synchronizability depends both on the coupling strength $\sigma$ and on the network topology [⁵⁰], as the latter determines the Laplacian eigenvalues. Crucially, the topology does not decide whether synchronization is possible (there always exists the synchronous state) but rather sets the threshold for achieving it. Specifically, ${\lambda }_2$ controls the critical value ${\sigma }_c=r^{\ast }/{\lambda }_2$[⁴¹].

In other words, any connected network can synchronize for a sufficiently large coupling, i.e., if $\sigma >{\sigma }_c$, regardless of the eigenvalue distribution: choosing a sufficiently large $\sigma$ guarantees that all transverse directions to $s$ have negative Lyapunov exponents, making the synchronization manifold transversely stable. Once $r^{\ast }$ is fixed by $f\left(x\right)$ and $h\left(x\right)$, the network’s role is purely to rescale this threshold through ${\lambda }_2$, dictating how much coupling effort is required for synchrony to emerge.

Physically, type II MSFs correspond to systems in which the coupling may eventually dominate internal chaotic tendencies. Once the normalized coupling strength parameter $r$ exceeds a critical threshold $r^{\ast }$, all transverse perturbations are damped and synchronizability becomes robust. In this case, the role of the network is not to determine whether synchronization is possible, but rather how easily it can be achieved.

2.4.3. Type III: $\Lambda =\Lambda \left(r\right)$ is a non-monotonic function with a local minimum

Here, $\Lambda \left(r\right)$ takes negative values within a finite interval $0⩽r_1<r<r_2$. Thus, synchronization is only possible if all eigenmodes $\sigma {\lambda }_i$ are within the interval $r_1<\sigma {\lambda }_i<r_2$[⁴¹]. This leads to two simultaneous constraints

$$\{\begin{array}{cc}\sigma {\lambda }_2>r_1,\hfill & \\ \sigma {\lambda }_N<r_2,\hfill & \end{array}$$

or equivalently, a restricted range for the coupling parameter:

$$\frac{r_1}{{\lambda }_2}<\sigma <\frac{r_2}{{\lambda }_N}.$$

If $\sigma >{\sigma }_{\text{max}}=r_2/{\lambda }_N$, at least the direction transverse to the synchronization manifold associated with the eigenvalue ${\lambda }_N$ becomes unstable, while if $\sigma <{\sigma }_{\text{min}}=r_1/{\lambda }_2$, in practice, the synchronization cannot be reached, since the smallest transverse eigenmode still has a positive Lyapunov exponent. Thus, synchronization is only stable when $\sigma$ is tuned within this “narrow” interval.

When $\sigma$ cannot be freely adjusted (or if synchronization is desired outside this critical range), the burden shifts entirely to network design: only certain network structures (those that shape ${\lambda }_2$ and ${\lambda }_N$ appropriately) can support stable synchronization.

In other words, not every network will admit synchronization in this case: only those whose Laplacian spectra satisfy the above constraints can do so. For MSFs of type III, this motivates the introduction of a network’s synchronizability [², ³⁸, ⁵¹], i.e., its structural capacity to support synchronization.

Here, synchronizability is quantified by the spectral ratio ${\lambda }_N/{\lambda }_2$ (which is always $>1$, since the Laplacian is positive semidefinite). A smaller ratio, i.e., closer to one, indicates a more compact spectrum and, therefore, better synchronizability.

More precisely, the condition for synchronizability can be written as:

(13) $$1\le \frac{{\lambda }_N}{{\lambda }_2}<\frac{r_2}{r_1},$$

which indicates that, for type III MSFs, synchronizability is fully determined by this spectral ratio [⁵², ⁵³]. Actually, in general, both the smallest nonzero eigenvalue and the eigenratio between the largest and the smallest nonzero eigenvalues of the Laplacian matrix can be used to characterize synchronizability [⁴⁹, ⁵⁴, ⁵⁵].

This means that synchronizability depends on how “bundled” the non-zero eigenvalues of the network’s Laplacian are: the closer they are to each other (i.e., the smaller the ratio ${\lambda }_N/{\lambda }_2$, or the closer to one), the greater the chance that all values $\sigma {\lambda }_i$ will fall within the stability interval $(r_1,r_2)$.

Nevertheless, it is worth emphasizing that these spectral indicators of synchronizability (i.e., ${\lambda }_2$ and ${\lambda }_N/{\lambda }_2$) are strictly defined within the MSF framework and that their relevance strongly depends on the qualitative type of the MSF. In particular, they provide necessary and convenient structural criteria for synchronizability only insofar as the transverse dynamics is well captured by the linear stability of the synchronization manifold. Outside the MSF framework, or in the presence of strong non-spectral dynamical effects, synchronizability cannot be fully characterized by Laplacian eigenvalues alone, especially when these non-spectral dynamical effects play a significant role [⁵⁵, ⁵⁶]. Recent studies have also proposed alternate frameworks to characterize network synchronizability that extend beyond spectral measures, highlighting structural or dynamical contributions not captured by Laplacian eigenvalues alone [⁵⁶].

Type III MSFs represent an intermediate physical scenario. Synchronization is possible, but only within a finite range of effective couplings [⁴¹]. If the coupling is too weak, the synchronous state cannot be stabilized; if it is too strong, the collective dynamics themselves can become unstable and synchronizability is lost [⁵⁷]. In this situation, synchronizability results from a delicate balance among dynamics, coupling strength, and network structure [⁴³], which explains why only certain network topologies are synchronizable in this case.

2.5. Scope and limitations of the MSF formalism

It is important to stress that the MSF formalism discussed in this work is inherently built upon the existence and invariance of the synchronization manifold, defined by the condition in equation (4). As a consequence, the MSF provides a rigorous criterion for the linear stability of this manifold with respect to transverse perturbations, but does not provide information about the global accessibility of the synchronized state from arbitrary initial conditions (i.e., it only guarantees local stability near the synchronization manifold), and it is not meant to describe synchronization phenomena that do not correspond to complete (i.e., identical) synchronization.

In particular, dynamical regimes such as partial or cluster synchronization, which are typically associated with network symmetries and the coexistence of distinct dynamical groups, fall outside the direct scope of the classical MSF approach [⁵⁸, ⁵⁹]. Likewise, phenomena such as intermittent synchronization and chimera states, which involve multistability or the absence of a single invariant synchronization manifold, cannot be captured within the standard MSF framework and require alternative theoretical descriptions [⁶⁰, ⁶¹]. Related symmetry-breaking collective behaviors in biological and social systems further illustrate the diversity of synchronization mechanisms beyond complete synchronization [⁶²].

In this sense, the MSF should be understood as a powerful and elegant tool for addressing a specific and well-defined synchronization problem (namely, the linear stability of complete synchronization by a set of identical systems) rather than as a general framework for all possible collective dynamical states in complexnetworks.

Examples: MSFs for Some Typical Nonlinear Dynamical Systems

In this section we will study three well-known nonlinear dynamical systems, which were coupled following a scheme in which the systems interact with each other through only one component. For example, for state vectors $⟨x^1,x^2,x^3⟩$, i.e., for three-dimensional systems described by the dynamical variables $(x,y,z)$, nine linear diffuse coupling configurations are possible: $x\to x$, $x\to y$, $x\to z$, $y\to x$, $y\to y$, $y\to z$, $z\to x$, $z\to y$ and $z\to z$, where the notation $x^i\to x^j$ represents a component-wise coupling scheme, from the $i$-th component of a dynamical system to the $j$-th component of another system (see Subsection 3.2 for an explicit example of the $x\to y$ coupling for the Lorenz system). The MSF was calculated and plotted for each of these possible configurations.

It is worth stressing that these linear coupling schemes are chosen here solely for illustrative and pedagogical purposes. The MSF formalism itself does not require the coupling function to be linear, nor to be applied to only one of the coordinates, provided it is at least continuously differentiable, so that the linearization and the variational analysis along the synchronous trajectory are well defined.

We focus on the Lorenz system [⁶³], the forced cubic Duffing oscillator (a paradigmatic forced, nonlinear oscillator) [⁶⁴], and the forced Van der Pol oscillator [⁶⁵]. In each case, the parameters were taken from the literature (see the corresponding reference) so that the resulting attractor is chaotic⁸, that is, with a positive Lyapunov exponent, i.e., $\Lambda \left(0\right)>0$.

Remark. The examples below are intentionally chosen within the standard continuous-time MSF setting. Nonetheless, it is useful to briefly comment on closely related extensions and modern developments.

First, it is worth recalling that, although we focus here on continuous-time chaotic systems described by ODEs, the MSF formalism is not restricted to flows. In fact, with minor adaptations [²⁴], it can also be applied to coupled discrete-time dynamical systems (i.e., chaotic maps [⁶⁶, ⁶⁷, ⁶⁸]). From a pedagogical standpoint, these examples may be particularly appealing since discrete-time chaos can exhibit strongly irregular (or pseudo-random-like) time series while remaining fully deterministic. An exploration of synchronization in coupled chaotic maps is left as a possible extension for interested readers.

Second, the MSF approach presented in this work corresponds to the classical setting of fixed coupling. Nevertheless, other extensions can naturally complement this introductory framework. In particular, adaptive synchronization strategies dynamically adjust coupling parameters in response to the synchronization error [⁶⁹, ⁷⁰, ⁷¹], improving robustness when the appropriate coupling strength is not known a priori. Such schemes may tune either an effective coupling gain [⁷²] or symmetry/structure-related aspects of the coupling [⁷³, ⁷⁴], aiming to keep the transverse dynamics within MSF-stable regimes (i.e., largest transverse Lyapunov exponent negative). A detailed implementation and comparison of adaptive controllers is left as a natural extension, since it would considerably expand the scope of the present introductory work.

In addition to adaptive strategies, another related topic that has recently attracted attention is time-reversible synchronization in chaotic systems [⁷⁵, ⁷⁶]. While the classical MSF framework provides a rigorous criterion for the linear stability of the synchronous manifold under fixed diffusive coupling, time-reversible synchronization may involve additional structural ingredients (e.g., nonstandard coupling architectures, symmetry constraints, or extended dynamical formulations) [⁷⁶, ⁷⁷, ⁷⁸]. A detailed treatment of this phenomenon would require a dedicated discussion beyond the present introductory scope, but we mention it here as an interesting modern development that may motivate further reading and student projects.

Finally, regarding discrete-time (digital) chaotic systems with a symmetry parameter, one may qualitatively expect such a parameter to influence the MSF through its effect on the synchronous trajectory $s\left(t\right)$ and, consequently, on the Jacobian $Df\left(s\right)$ (and possibly on $Dh$) [⁷⁹, ⁸⁰, ⁸¹, ⁸²]. In other words, changing a symmetry-related parameter typically alters the underlying local dynamics that enter the variational equation, and therefore modifies the resulting MSF profile $\Lambda \left(r\right)$. A parametric study of this dependence, however, would require introducing a specific map model and performing a dedicated numerical analysis, which we leave as a possible extension.

We now proceed to the case studies (i.e., numerical examples). The following subsection describes the computational steps used to obtain the MSF for each system and coupling configuration.

3.1. Methodology and computational details

The largest Lyapunov exponent, determined by the variational equation (12), was calculated as follows. We define the variational Jacobian matrix as $J\left(s\right)=Df\left(s\right)-rDh\left(s\right)$ and consider the time-dependent linear equation

$$\dot{\xi }\left(t\right)=J\left(s\right)\xi \left(t\right).$$

The trajectory $s=s\left(t\right)$ is previously obtained by numerical integration of the isolated system $\dot{s}=f\left(s\right)$, as in equation (1), using the odeint method from the SciPy library (adaptive step-size), with initial condition $s\left(0\right)=s_0=\left(\mathrm{1,\; 1,\; 1}\right)$, and in the interval $0\le t\le {10}^{\text{3}}$.

In general, for nonlinear dynamical systems, the Jacobian matrix $Df\left(s\right)$ depends explicitly on the trajectory $s\left(t\right)$. However, for a linear coupling function $h\left(x\right)$, the corresponding Jacobian matrix $Dh\left(s\right)$ is a matrix with constant entries.

In particular, in this case where we consider a single-component coupling scheme, the dynamical systems interact with each other through only one component. As a result, $Dh\left(s\right)$ will have only one non-zero entry, while all other entries will be equal to zero. For instance, for state vectors $x=⟨x^1,x^2,x^3⟩$ (i.e., 3D dynamical systems), in the situation where the second component of one system is coupled to the first component of another (coupling $x^1\to x^2$, i.e., from the first component of one system to the second component of another), the coupling function reads $h\left(x\right)=⟨0,x^1\mathrm{,\; 0}⟩$, and the only non-zero entry of $Dh\left(s\right)$ will be ${\left[Dh\left(s\right)\right]}_{21}=1$. In this case, the explicit coupled network equations read

$$\{\begin{array}{cc}{\dot{x}}_i^1=f^1(x_i^1,x_i^2,x_i^3),\hfill & \\ {\dot{x}}_i^2=f^2(x_i^1,x_i^2,x_i^3)-\sigma {\sum }_{j=1}^N{l}_{ij}{x}_j^1,\hfill & \\ {\dot{x}}_i^3=f^3(x_i^1,x_i^2,x_i^3),\hfill & \end{array}$$

for $i=1,\dots ,N$. Note that only the equation of the second component receives the coupling term. All other single-component coupling configurations considered in this work follow analogously from the choice of the coupling function $h$.

For each value of $r$, the initial perturbation ${\xi }_0$ was randomly sampled and normalized. Its evolution was calculated iteratively using the Euler method, with a fixed step size $\delta t={10}^{-3}$, over the interval $0\le t\le {10}^{\text{3}}$.

Because the systems studied in this work may operate in chaotic regimes (actually, the parameters were intentionally chosen so that each isolated system exhibits a chaotic attractor), it is important to note that numerical discretization may significantly affect the dynamics of chaotic systems, and that low-order explicit schemes (e.g., the explicit Euler method) can, in some situations, introduce spurious stabilization or even suppress chaotic behavior. This issue has been extensively discussed in the literature, particularly in the context of bifurcation diagrams and long-time integrations of chaotic flows [⁸³, ⁸⁴]. In this paper, however, the numerical integration is not intended to provide high-precision quantitative predictions, nor to compare different numerical schemes. Our focus is on the qualitative phenomenology of the MSF.

For this purpose, we used SciPy’s odeint (adaptive step-size solver) to obtain the synchronized trajectory $s\left(t\right)$, but used Euler’s method with a sufficiently small time step to integrate the linear variational equation for $\xi \left(t\right)$ for each value of $r$. More sophisticated schemes (e.g., semi-implicit or symplectic-type integrators, which are particularly natural in Hamiltonian/conservative settings [⁸⁵]) may certainly improve numerical fidelity and constitute an interesting extension activity, particularly from a didactic perspective, but a systematic comparison of numerical methods lies beyond the scope of this introductory study. In this sense, the chosen integrator serves primarily as an accessible computational tool to support the conceptual discussion of MSF-based stability.

After integrating the variational equation for $\xi \left(t\right)$, at each step $\delta t$, the norm of the perturbation was accumulated to calculate the Lyapunov exponent, and the perturbation was renormalized to avoid numerical exponential growth or disappearance,

$$\Lambda \left(r\right)=\underset{M\to \infty }{lim}\frac{1}{T}\sum _{k=1}^M\ln \Vert {\xi }_k\Vert ,$$

where $T=M\delta t$ is the total integration time.

The function $\Lambda \left(r\right)$ was obtained by repeating this procedure for values of $r$ uniformly distributed over a given interval, always using $300$ or more points. For each coupling scheme, the $r$-axis range was adjusted individually to include all intersection points with $\Lambda \left(r\right)=0$, allowing visualization of the complete MSF profile and all relevant transverse stability intervals. The zeros of the MSF were determined numerically using the bisection method, with a relative tolerance of ${10}^{-6}$.

All simulations were implemented in Python 3.11, using the NumPy, SciPy, and Matplotlib libraries. The code was developed and executed in the Jupyter Notebook 6.5.4 environment, with modest computational resources (single CPU, no parallelization). Python programs can be downloaded from GitHub at⁹https://github.com/javierg0mez/MSF_examples.git

3.2. The Lorenz system

The Lorenz system is described by the differential equations

(14) $$\{\begin{array}{cc}\dot{x}=\sigma \left(y-x\right),\hfill & \\ \dot{y}=x\left(\rho -z\right)-y,\hfill & \\ \dot{z}=xy-\beta z,\hfill & \end{array}$$

where $\sigma =10$, $\rho =28$ and $\beta =8/3$[⁶³].

The Jacobian matrix of the system is

(15) $$Df\left(s\right)=\left(\begin{array}{ccc}-\sigma & \sigma & 0\\ \rho -z^{\ast }& -1& -x^{\ast }\\ y^{\ast }& x^{\ast }& -\beta \end{array}\right),$$

where $x^{\ast }$, $y^{\ast }$ and $z^{\ast }$ are the time-dependent components of the solution $s$.

For pedagogical clarity, it is important to understand how the diffusive coupling introduced above translates into the network equations for the Lorenz system. Consider, e.g., a diffusive coupling acting only from the first component of one system to the second component of another (i.e., coupling $x\to y$). In this case, the network dynamics takes the form

$$\{\begin{array}{cc}{\dot{x}}_i=\sigma \left(y_i-x_i\right),\hfill & \\ {\dot{y}}_i=x_i\left(\rho -z_i\right)-y_i-\sigma {\sum }_{j=1}^N{l}_{ij}{x}_j,\hfill & \\ {\dot{z}}_i=x_i{y}_i-\beta z_i,\hfill & \end{array}$$

for $i=1,\dots ,N$. The explicit equations for the remaining coupling schemes and dynamical systems considered in this work follow the same construction for each local dynamics $f$ and the corresponding choice of the coupling function $h$, as detailed in Section 2.

The MSFs for each of the possible linear coupling configurations between the variables of the Lorenz system are shown in Figure 3. We note that, for some coupling schemes, the MSF remains strictly positive for all $r>0$, with no finite crossings of the $r$-axis, i.e., type I MSFs. This behavior is observed in the couplings $x\to z$, $y\to z$, $z\to x$, and $z\to y$. In these cases, no network of Lorenz systems coupled through these schemes will allow the emergence of the synchronized state, since all modes transverse to the synchronization manifold are unstable [⁴⁵] and therefore cannot suppress the intrinsic chaotic divergence of the Lorenz attractor [⁸⁶]. Consequently, small perturbations around the synchronization manifold will grow exponentially, driving the system away from the synchronized state. This lack of synchronizability is independent of both the coupling strength and the network topology.

Figure 3
Lorenz system, given by equation (14), for $\sigma =10$, $\rho =28$ and $\beta =8/3$[⁶³]. MSFs vs. the normalized coupling parameter $r$, under different coupling schemes. For each coupling scheme, the $r$-axis range was chosen to include all intersections with $\Lambda \left(r\right)=0$ and to capture the overall qualitative MSF profile (including the number of crossings and the existence of transverse stability regions). The same scaling criterion is adopted in Figures 4 and 5.

On the other hand, certain coupling schemes exhibit qualitatively different behavior, characterized by the presence of one or more finite crossing points of the MSF with the $r$-axis. At these points, the master function changes sign, indicating the existence of stability regions where the synchronization can occur, at least for initial conditions close to the synchronization manifold.

The first of these, e.g., is the coupling $x\to x$, which exhibits a single intersection point of the MSF with the $r$-axis, at $r_1=8.31$. Beyond this value, the MSF remains negative. In this case, every network of Lorenz systems coupled through their $x$ coordinates is potentially synchronizable whenever $\sigma {\lambda }_2>8.31$, i.e., requiring a minimum coupling strength $\sigma >{\sigma }_c=8.31/{\lambda }_2$, where ${\lambda }_2$ is the smallest non-trivial eigenvalue of the network’s Laplacian matrix. Physically, this indicates that coupling through the $x$ variable allows transverse instabilities to be effectively suppressed once the coupling exceeds a critical threshold [⁸⁶]. It is worth recalling that, for this type of MSF (i.e., type II MSF), given a fixed $\sigma$, we will always have the option of choosing networks with greater connectivity (greater ${\lambda }_2$), which allow synchronization for lower coupling strengths.

Similarly, for the coupling $x\to y$, the MSF crosses the $r$-axis at $r_1=7.221$, which indicates behavior similar to the previous case, with only a slight difference in the critical coupling strength required to ensure the transverse stability of the synchronization state. The coupling $y\to y$, on the other hand, shows stability starting at $r_1=2.904$, a value lower than the previous ones. This suggests that the coupling $y\to y$ is one of the most efficient ways (in terms of coupling “cost”) to achieve synchronizability in networks of coupled Lorenz systems.

In the case of the coupling $y\to x$, the MSF has two finite points of intersection with the $r$-axis, at $r_1=4.176$ and $r_2=23.517$, defining a single synchronizability band, which lies between two regions of instability. Within this range, the MSF takes negative values and, therefore, the transverse stability of the synchronization manifold is guaranteed. Outside this range, the stability of the synchronous state (and, consequently, the possibility of observing it) cannot be assured. In this scenario, the stability is more restrictive: for a fixed network topology and coupling strength, it is necessary that all products $\sigma {\lambda }_i$ over the entire Laplacian spectrum (for $i=\mathrm{1,2},\dots ,N$) fall within the interval $4.176<\sigma {\lambda }_i<23.517$ for the synchronized solution to be stable. This finite stability window is the “hallmark” of type III MSFs: for a network of Lorenz systems coupled via the $y\to x$ scheme, the synchronization manifold is transversely stable only within a limited range of coupling strengths.

The case of the coupling $z\to z$ is particularly interesting, as the master function falls outside the types I $\to$ III described above. Here, the MSF crosses the $r$-axis at three points: $r_1=1.473$, $r_2=6.065$, and $r_3=114.468$, creating two stability regions: the open interval $\left(\mathrm{1.473,\; 6.065}\right)$ and the range $r>114.468$. This implies that, in this configuration, the synchronization can only be observed if all products $\sigma {\lambda }_i$, for $i=\mathrm{2,\; 3},\dots ,N$, lie within one of these regions. However, the first region is quite narrow, which makes it difficult for many $\sigma {\lambda }_i$ values to fall inside it, especially in network topologies with broader Laplacian spectra. Conversely, reaching the second stability region, which would offer greater robustness, requires high coupling strength values $\sigma$, which may also be challenging in practice. While coupling through $z$ can stabilize the synchronization state in principle, it does so either within a very narrow parameter range or at the cost of extremely strong coupling strengths.

3.3. The forced Duffing oscillator

The forced cubic Duffing oscillator is governed by the differential equations

(16) $$\{\begin{array}{cc}\dot{x}=y,\hfill & \\ \dot{y}=-\delta y-\alpha x-\beta x^3+\gamma \sin \left(\omega t\right),\hfill & \end{array}$$

where $\omega =1$, $\delta =0.1$, $\alpha =0$, $\beta =1$ and $\gamma =7$[⁶⁴, ⁸⁷].

The Jacobian matrix of the system is

(17) $$Df\left(s\right)=\left(\begin{array}{cc}0& 1\\ -\alpha -3\beta {\left(x^{\ast }\right)}^2& -\delta \end{array}\right),$$

where $x^{\ast }$ denotes the first component of the trajectory $s\left(t\right)$.

There are three types of MSFs obtained for the forced Duffing oscillator under the four linear coupling configurations analyzed (Figure 4). Namely:

Figure 4
Forced Duffing oscillator, given by equation (16), with $\omega =1$, $\delta =0.1$, $\alpha =0$, $\beta =1$, and $\gamma =7$[⁶⁴]. MSFs as a function of the normalized coupling parameter $r$, for different coupling schemes.

In the couplings $x\to x$ and $y\to y$, the MSF crosses the $r$-axis once, at $r_1=0.161$ and $r_1=0.156$, respectively. Beyond these values, $\Lambda \left(r\right)<0$, meaning that any network satisfying $\sigma {\lambda }_2>r_1$ already guarantees the transverse stability of the synchronization manifold. This simple single-threshold behavior is characteristic of “diagonal” coupling schemes, which act directly on the corresponding variable and can efficiently suppress transverse perturbations [⁶⁴, ⁸⁷]. Note that these thresholds are relatively low (on the order of ${10}^{-1}$), so even modest coupling strengths, even in sparsely connected topologies, are sufficient to ensure stability of the synchronized state.

For the coupling $x\to y$, the MSF has three finite intersection points with the $r$-axis, at $r_1=0.95$, $r_2=2.767$ and $r_3=3.608$. In this case, we have two stability regions: $0.95<\sigma {\lambda }_i<2.767$ and $\sigma {\lambda }_i>3.608$. The appearance of disconnected stability regions already indicates that a non-diagonal coupling scheme requires greater control over both the coupling strength and network topology [⁸⁷]. The interaction does not act on the same variable, making the synchronizability slightly more sensitive to the normalized coupling strength [⁶⁴]. The first stability window is relatively narrow, which requires either networks with a tightly clustered Laplacian spectrum or a precise tuning of $\sigma$ to maintain stability of the synchronization manifold. The second region requires a little stronger coupling.

Finally, for the coupling $y\to x$, the MSF changes sign six times between $r_1=0.126$ and $r_6=1$, delimiting three very narrow stability bands, separated by regions of instability. This highly fragmented structure is a hallmark of “ragged synchronizability”, in which the synchronous state alternates between stability and instability as the coupling parameter is varied. Therefore, network synchronizability is only possible if all products $\sigma {\lambda }_i$ fall entirely within one of these bands. However, in networks with a dispersed Laplacian spectrum, this condition is hard to meet [⁶⁴], since even a single $r_i=\sigma {\lambda }_i$ outside the stability range is enough for the corresponding eigenmode of a small perturbation to grow exponentially and push the trajectory away from the synchronization manifold. Moreover, increasing $\sigma$ introduces another risk: some $\sigma {\lambda }_i$ may exceed $r_6=1$, completely ruling out synchronization because $\Lambda \left(r\right)$ remains strictly positive for $r>1$.

Once again, we see that, in the particular case of the forced Duffing oscillator, certain coupling schemes clearly favor synchronizability. Specifically, the couplings $x\to x$ and $y\to y$, i.e., when the chaotic Duffing oscillators are coupled either through their state variables (coordinates $x$) or through their velocities (coordinates $\dot{x}=y$), respectively. From a physical standpoint, these direct (diagonal) couplings correspond to adding spring- or damper-like forces, as reported in the literature [⁸⁸, ⁸⁷, ⁸⁹], which provide a simple dissipative mechanism that favors the stability of the synchronous state. These two coupling schemes can be regarded as relatively efficient “routes” for synchronizing networks of forced Duffing oscillators. This can be translated, in these two configurations, as a certain degree of robustness: once the minimum coupling threshold required to guarantee the transverse stability of the synchronization manifold in a given network is surpassed, further increases in coupling strength do not compromise synchronizability; on the contrary, they tend to reinforce it. By contrast, the $x\to y$ and especially $y\to x$ schemes require tighter control over $\sigma$ and/or over the network topology [⁹⁰]. This underlines the importance of carefully choosing the coupling scheme when attempting to synchronize networks of chaotic systems.

3.4. The forced Van der Pol oscillator

The forced Van der Pol oscillator is described by the differential equations

(18) $$\{\begin{array}{cc}\dot{x}=y,\hfill & \\ \dot{y}=-x+\mu \left(1-x^2\right)y+A\sin \left(\omega t\right),\hfill & \end{array}$$

where we use $\mu =3$, $A=15$ and $\omega =4.065$[⁶⁵].

The Jacobian matrix of the system (again, disregarding the explicit time dependence to focus only on the autonomous part) is

(19) $$Df\left(s\right)=\left(\begin{array}{cc}0& 1\\ -1-2\mu x^{\ast }{y}^{\ast }& \mu \left(1-{\left(x^{\ast }\right)}^2\right)\end{array}\right),$$

where $x^{\ast }$ and $y^{\ast }$ are the components of the solution $s$.

The four MSFs of the forced Van der Pol oscillator (Figure 5) show that all analyzed coupling configurations allow synchronization, as each curve crosses the $r$-axis at least once, but with varying degrees of “ease”. For instance, in the $x\to x$ and $y\to y$ couplings, the function $\Lambda \left(r\right)$ has a single intersection with the $r$-axis, which is typical of type II MSFs, at significantly low values: $r_1=0.146$ for $x\to x$ and $r_1=0.376$ for $y\to y$. These thresholds suggest relatively easy synchronization: for any network of chaotic Van der Pol oscillators coupled through their state variables ($x$ coordinates) or their velocities ($y=\dot{x}$ coordinates) that satisfies $\sigma {\lambda }_2>r_1$ (even with sparse connectivity), we will have the chance to observe the synchronization state. Furthermore, increasing $\sigma$ beyond this threshold only reinforces this chance, as the MSF remains negative. This robust stability aligns with the physical interpretation of these diagonal couplings as being analogous to adding simple, direct restoring (elastic) or damping (dissipative) forces between the oscillators, as reported in the literature [⁸⁹].

Figure 5
Forced Van der Pol oscillator, given by equation (18), for $\mu =3$, $A=15$, and $\omega =4.065$[⁶⁵]. MSFs vs. the normalized coupling parameter $r$, under different coupling schemes.

Likewise, when coupling the state variables to the velocities of the oscillators (i.e., $x\to y$), the MSF has a single crossing point at $r_1=1.372$, a threshold on the order of unity. This indicates that moderate coupling strengths or slightly higher network connectivity are enough to place all $\sigma {\lambda }_i$ products in the negative region of the MSF, ensuring the stability of the synchronization manifold.

However, when coupling velocities with dynamic variables (coupling $y\to x$), the picture changes and additional challenges appear. The MSF changes sign four times between $r_1=0.035$ and $r_4=1.027$, opening only two very narrow stability bands, separated by a broad region of instability. In this configuration, the synchronization state is stable if the entire rescaled Laplacian spectrum $\left\{\sigma {\lambda }_i\right\}$ of the network lies entirely within one of these two windows. This makes the scheme highly demanding, requiring either tightly compressed spectra or precise tuning of $\sigma$. Again, it is worth mentioning that increasing the coupling strength is not helpful (in fact, it can be detrimental), since if any $\sigma {\lambda }_i$ exceeds $r_4$, the function $\Lambda \left(r\right)$ becomes positive in that $i$-th eigen-direction, and the stability of the synchronization state is lost. Therefore, the $y\to x$ scheme requires strict control over both the network spectrum and coupling strength.

This extreme parameter sensitivity and fragmented stability are hallmarks of systems in which the coupling scheme (from velocity to position) interferes destructively with the intrinsic nonlinear damping of the Van der Pol oscillator, making the stability of the synchronous state fragile.

So, as in the Duffing oscillator, the diagonal couplings ($x\to x$ and $y\to y$) are simpler and more robust [⁸⁹], while the cross couplings are more delicate, especially $y\to x$.

Summary and Classification of Observed MSFs

Based on the numerical results from the previous section and recalling that each system studied was selected in its chaotic regime, meaning that all MSFs at $r=0$ coincide with the largest Lyapunov exponent of the isolated system, which is positive: $\Lambda \left(0\right)>0$. In this context, the synchronous solution is unstable and cannot be physically obtained whenever the MSF remains positive, $\Lambda \left(r\right)>0$, since small perturbations of the synchronous state lead to trajectories that diverge from it. So, for the synchronization manifold to be stable, and for a small perturbation around the synchronous state to decay exponentially, the MSF (which starts out positive) must change sign, i.e., it must cross the $r$-axis.

Under this criterion, we can classify each of the MSFs in terms of the behavior of $\Lambda \left(r\right)$ crossing the $r$-axis, i.e., according to the number of finite crossings with the $r$-axis. This provides an operational way to characterize the synchronization conditions in a network, for any given combination of local dynamics $+$ coupling scheme.

However, although the numerical results presented in the previous section focus on three representative chaotic systems (Lorenz, Duffing and Van der Pol), our goal in what follows is not to claim universality based solely on these examples. Rather, these systems were chosen because they illustrate distinct and historically relevant classes of nonlinear behavior (strange attractor, forced nonlinear oscillator, and relaxation dynamics), as well as for their authors’ contributions to nonlinear dynamics and the study of chaos. Moreover, their MSFs exhibit qualitative profiles that have already been reported in the literature for a wide range of coupled nonlinear systems under different coupling schemes (see, e.g., [²⁰, ⁴¹, ⁴³, ⁴⁵, ⁵², ⁵⁷, ⁶⁴, ⁸⁷]). Therefore, we emphasize that the classification proposed below should be understood as a phenomenological taxonomy of qualitative MSF profiles (e.g., absence of stable regions, stability beyond a single threshold, or finite stability windows), which reflects and synthesizes recurring qualitative MSF behaviors observed in coupled nonlinear dynamical systems, and is not restricted only to the specific examples discussed in this paper.

4.1. Classification of observed MSF profiles

Below, we outline the five MSF profiles that we found in the systems studied.

① $\Lambda \left(r\right)$ has no finite intersection points with the $r$-axis

When the MSF remains positive for all $r>0$, the synchronization manifold is transversely unstable for any network topology and any coupling strength. In this scenario, no choice of $\sigma$ or network structure achieves the stability of the synchronous state. This behavior corresponds to what we call type I MSF in Subsection 2.4. Examples include the $x\to z$ and $z\to y$ coupling schemes for the Lorenz system. In such cases, we say that these coupling schemes, as well as the $y\to z$ and $z\to x$ couplings, are inherently unfavorable to the stability of the synchronization manifold in a network of coupled Lorenz systems.

② $\Lambda \left(r\right)$ has only one finite intersection point, $r_1>0$

This MSF profile corresponds to the type II case discussed in Subsection 2.4, which characterizes networks that are potentially synchronizable.

When the MSF crosses the $r$-axis only once at $r_1$ and remains negative for all $r>r_1$, synchronizability is guaranteed as soon as $\sigma {\lambda }_2>r_1$, where ${\lambda }_2$ is the smallest non-trivial eigenvalue of the Laplacian matrix. Here, low thresholds, such as $r_1\approx 0.16$ in the $x\to x$ and $y\to y$ couplings of the Duffing oscillator, or $r_1=0.146$ in the $x\to x$ coupling of the Van der Pol oscillator, indicate that even weakly connected networks (i.e., with small ${\lambda }_2$ or a dispersed spectrum) may synchronize with modest couplings. While higher values, such as $r_1=8.31$ in the $x\to x$ coupling of the Lorenz system, require larger $\sigma$ or networks with higher connectivity (larger ${\lambda }_2$).

③ $\Lambda \left(r\right)$ has two finite intersection points, $r_1$ and $r_2$, where $0<r_1<r_2$

If the MSF changes sign at $r_1<r_2$ and becomes positive again for all $r>r_2$, the synchronous state is stable only within the finite interval $r_1<\sigma {\lambda }_i<r_2$, for all $i=\mathrm{2,\; 3},\dots ,N$. Therefore, as discussed in Subsection 2.4, the two simultaneous conditions that allow synchronization are $r_1<\sigma {\lambda }_2$ and $\sigma {\lambda }_N<r_2$, where ${\lambda }_N$ is the largest eigenvalue of the Laplacian matrix. For this MSF type, a network is synchronizable if its spectral ratio obeys $1\le {\lambda }_N/{\lambda }_2<r_2/r_1$. An example of this MSF profile is the coupling $y\to x$ of the Lorenz system, where $r_1=4.176$ and $r_2=23.571$.

In this case, wide stability windows provide more flexibility in choosing the coupling strength and tolerating dispersion in the Laplacian spectrum. In contrast, narrow windows require highly compressed spectra. However, even so, we have established that the more tightly packed the Laplacian spectrum is (i.e., the smaller the ratio ${\lambda }_N/{\lambda }_2$ or the closer it is to one), the greater the likelihood that all products $\sigma {\lambda }_i$ fall within the stability range and, therefore, ensuring the transverse stability of the synchronization manifold.

④ $\Lambda \left(r\right)$ has an odd number of finite intersections, where $0<r_1<r_2<\cdots <r_{2k+1}$, with $k\in {\mathbb{Z}}^{+}$

When the MSF has, e.g., three finite intersections $r_1$, $r_2$, and $r_3$ with the $r$-axis, as in the $z\to z$ coupling of the Lorenz system, the function $\Lambda \left(r\right)$ starts positive (ad hoc condition of the systems studied), crosses the axis at $r_1>0$ and changes sign, then crosses the axis again at $r_2>r_1$, returning to positive values. Unlike the previous case, however, it does not remain positive for all $r>r_2$: there is a third crossing at $r_3>r_2$, after which $\Lambda \left(r\right)$ becomes negative and stays below the axis for all larger $r$.

The same reasoning applies to cases with five, seven, nine, or any other odd number of crossings. In such cases, the MSF alternates between stable and unstable regions, i.e., it changes sign $2k+1$ times, where $k\in {\mathbb{Z}}^{+}$. Each interval where $\Lambda \left(r\right)<0$ defines a “window” of possible synchronization, while the intervals with $\Lambda \left(r\right)>0$ make it impossible to observe the synchronous state. In general, however, the asymptotic behavior of the function $\Lambda \left(r\right)$ is always negative as $r\to \infty$, i.e., there is a final crossing point $r_{2k+1}$ beyond which $\Lambda \left(r\right)<0$ for all larger $r$. These MSF profiles are particularly interesting in problems known as irregular synchronization [⁶⁴].

For pedagogical purposes, let us focus on the simpler case of only three intersection points. To ensure the stability of the synchronous state, we need all $\sigma {\lambda }_i$ must lie within the union of the two stability intervals, $(r_1,r_2)\cup (r_3,\infty )$, i.e., there must exist an index $j=\mathrm{3,\; 4},\dots ,N$ such that the following conditions hold: $r_1<\sigma {\lambda }_2$, $r_2>\sigma {\lambda }_j$ and $r_3<\sigma {\lambda }_{j+1}$. Alternatively, we may also have “extreme” cases in which all $\sigma {\lambda }_i$ lie entirely within a single stability interval. In this case, it is sufficient that the spectral ratio satisfies ${\lambda }_N/{\lambda }_2<r_2/r_1$, ensuring that all $\sigma {\lambda }_i$ fall within the first synchronization window $(r_1,r_2)$, or that $\sigma {\lambda }_2>r_3$, guaranteeing that all $\sigma {\lambda }_i$ fall within the second stability region.

In general, for MSFs with an odd number of crossings, a network is synchronizable if all $\sigma {\lambda }_i$ lie within the union of stability intervals $(r_1,r_2)\cup (r_3,r_4)\cup \cdots \cup (r_{2k+1},\infty )$. It is enough for a single transverse mode to fall into an interval $(r_{2l},r_{2l+1})$, for any $l\in {\mathbb{Z}}^{+}$, where $\Lambda \left(r\right)>0$, to destroy the stability of the synchronization manifold. Moreover, the larger the number of intersections, the narrower and more delicate these stability windows tend to be, requiring a more tightly packed Laplacian spectrum and a carefully tuned choice of $\sigma$.

However, as mentioned earlier, one common feature of these MSFs is that, given any topology, whenever the coupling is sufficiently strong (i.e., $\sigma >r_{2k+1}/{\lambda }_2$), the network will always be synchronizable. In the Duffing oscillator (with $x\to y$ coupling), we saw an example with three crossings, with $r_3=3.608$. In this situation, it would be enough, in principle, to push all products $\sigma {\lambda }_i$ into the final stability region by increasing $\sigma$, as long as the network supports such coupling strength. But, e.g., how practical (or even feasible) would this be for Lorenz systems coupled through their $z$ coordinates, where $r_3=114.468$?

⑤ $\Lambda \left(r\right)$ has an even number of finite intersections, where $0<r_1<r_2<\cdots <r_{2k}$, where $k\in {\mathbb{Z}}^{+}$

This is the most delicate case, requiring tighter control over both the network spectrum and the coupling strength, which makes synchronizability somewhat more difficult to achieve. MSFs that change sign $2k$ times (i.e., that have $2k$ crossing points with the $r$-axis) for $k\in {\mathbb{Z}}^{+}$ generate $k$ narrow stability windows interspersed between $k+1$ instability regions. Therefore, for MSFs of this type here, a network is synchronizable only if the entire rescaled Laplacian spectrum $\sigma {\lambda }_i$ lies entirely within the finite union of intervals $(r_1,r_2)\cup (r_3,r_4)\cup \cdots \cup (r_{2k-1},r_{2k})$, where the values of $\sigma {\lambda }_i$ may either occupy multiple intervals or be confined to just one of them.

Unlike the case where $\Lambda \left(r\right)$ has an odd number of finite crossings with the $r$-axis, here increasing $\sigma$ indiscriminately (e.g., taking $\sigma >r_{2k}/{\lambda }_N$) can push some modes out of the last window, rendering the system unsynchronizable regardless of the network topology.

An example of this behavior is found in the $y\to x$ coupling of the Duffing oscillator, where the MSF exhibits six crossings between $r_1=0.126$ and $r_6=1$. For a network of Duffing oscillators coupled through $y\to x$ to be synchronizable, its entire rescaled Laplacian spectrum $\sigma {\lambda }_i$ must fall within one of these narrow stability windows, each with a width of ≲ 0.1.

4.2. A step-by-step guide for setting up synchronizable networks

Each of the numerical results obtained and discussed in the previous section fits into one of the five classes defined above. From the perspective of designing synchronizable networks, we suggest the following practical and fairly straightforward strategy:

• ✐

  First, select the local dynamics and the coupling function that generate an MSF, preferably of class ② (one crossing) or, at most, of class ③ (two crossings). In these cases, simple spectral conditions and/or a moderate increase in $\sigma$ are usually sufficient to guarantee synchronizability.

• ✐

  If the MSF falls into class ④ or ⑤, greater caution is required: one must select topologies whose spectra are already naturally “compressed” and, at the same time, carefully tune $\sigma$ so as to keep all products $\sigma {\lambda }_i$ within the stability windows.

• ✐

  Finally, adjust the network topology (adjacency matrix, average degree, degree distribution, weights, etc.) until the required inequalities are satisfied. If there is still room for adjustment, $\sigma$ can then be increased, provided it remains within the stability windows.

Thus, the classification proposed here provides a practical “roadmap” for interpreting and comparing, in just a few minutes of calculation, all the cases detailed in the previous section and using only two or three parameters (such as the MSF class, the spectral ratio ${\lambda }_N/{\lambda }_2$ and the coupling strength $\sigma$). In this way, it becomes possible to identify which combinations of local dynamics, coupling scheme, and network structure are synchronizable.

4.3. Asymptotic MSF behavior and connection with classical types

Finally, it is worth noting that the asymptotic behavior of the MSF for large $r$ is locally robust¹⁰, a result already reported in the literature, see [²⁰]. The authors showed that, when an internal parameter of the local dynamics is slowly varied, the intersection points of the MSF with the $r$-axis are always created or annihilated in pairs. Consequently, the total number of crossings preserves its parity (even or odd). Thus, the MSF may “jump” from one class to another of the same parity (i.e., from even to even or from odd to odd), for instance, an MSF of class ③ (with two crossings) may shift to class ⑤ (four, six, …), but never to ② or ④. Likewise, an MSF of class ② can only jump to ④.

For example, in the Duffing oscillator (with coupling $x\to y$), when the parameter $\gamma$ is varied from $7$ to $11.3$, the master function changes from having three to five crossing points. This “parity principle” thus provides an additional guideline: moderate changes in the parameters do not alter the class (even or odd) of the MSF, which makes it easier to anticipate how a network will respond to small uncertainties in the internal parameters of the local dynamics.

It should be emphasized that this classification into classes (① to ⑤) complements, rather than replaces, the classic typology into types I, II, and III discussed in the theoretical framework (Subsection 2.4). While types I to III describe the asymptotic behavior of the MSF ($r\to \infty$), the classification into classes offers a more operational perspective, based on the number of finite intersections with the $r$-axis.

Therefore, each traditional type can be associated with one or more classes: e.g., type I corresponds to class ①, type II to class ②, and type III to class ③. The additional classes (④ and ⑤) naturally extend the original typology, accommodating more complex situations such as those observed in the Duffing oscillator ($y\to x$) or the Lorenz system ($z\to z$). In other words, every type (I, II, III) maps to a class, but not every new class fits into the original scheme. Still, the parity rule guarantees that the asymptotic behavior always remains consistent with the I $\to$ III logic.

What really matters is the limit ($r\to \infty$): it determines whether the MSF ends up being “type II” if $\Lambda \left(r\right)<0$ for all sufficiently large $r$ (equivalently, starting from $\Lambda \left(r\right)>0$, there is an odd number of finite zero-crossings), or “type III” if $\Lambda \left(r\right)>0$ for all sufficiently large $r$ (the negative region is bounded in $r$; equivalently, there is an even number of finite zero-crossings¹¹). The classification into classes simply details how many crossings/windows occur along the way, i.e., the “middle ground”.

This parity principle ensures compatibility between the two schemes: even if a curve “jumps” from one class to another as the parameters of the local dynamics vary, the asymptotic behavior always preserves the I $\to$ III logic. Thus, the two frameworks are not strictly equivalent but rather compatible, with the guarantee of compatibility coming precisely from the behavior of the master function in the limit $r\to \infty$.

5. Extending the MSF to General Diffusive Coupling

One of our assumptions was that the dynamics of the $i$-th system is governed by the uncoupled vector field $f\left(x_i\right)$ together with a diffusive coupling term of the specific form $h\left(x_j\right)-h\left(x_i\right)$, for every pair of connected nodes $i,j$, as described in the model of equation (2). However, this assumption can be relaxed, allowing for a generalization of the previous results [⁴³].

Let us consider a function $g:{ℝ}^n\times {ℝ}^n\to {ℝ}^n$, assumed to be sufficiently smooth (at least locally Lipschitz continuous). We say that $g$ is diffusive if

$$\{\begin{array}{cc}g(x,x)=0,\hfill & \\ g(x,y)=-g(y,x).\hfill & \end{array}$$

Clearly $g(x,y)=h\left(y\right)-h\left(x\right)$ is diffusive, but so is $g(x,y)=\sin \left(y-x\right)$, that appears in the Kuramoto model [²¹, ²²], and many other functions of interest in non-linear dynamics.

With this, the model in equation (2) can be extended to general diffusive couplings as:

(20) $${\dot{x}}_i=f\left(x_i\right)+\sigma \sum _{j=1}^N{a}_{ij}g(x_j,x_i).$$

Note that the existence and invariance of the synchronized solution $s=x_1=x_2=\cdots =x_N$, such that $\dot{s}=f\left(s\right)$, are guaranteed by the diffusive nature of the coupling function $g(x,y)$. Actually, under the assumed regularity of the coupling function, the synchronized solution is well defined. With this in place, we can now analyze its linear stability.

The procedure is analogous to that described in Section 2. We again consider a small perturbation around $s\left(t\right)$,

$$\begin{array}{cc}\hfill x_i=s+\delta x_i& ⟹{\dot{x}}_i=\dot{s}+{\dot{\delta x}}_i=f\left(s\right)+{\dot{\delta x}}_i\hfill \\ & ⟹{\dot{\delta x}}_i={\dot{x}}_i-f\left(s\right).\hfill \end{array}$$

Substituting ${\dot{x}}_i$ from equation (20),

$${\dot{\delta x}}_i=f\left(x_i\right)+\sigma \sum _{j=1}^N{a}_{ij}g(x_j,x_i)-f\left(s\right),$$

and linearizing around $s$,

$$\begin{array}{l}{\dot{\delta x}}_i={f\left(x_i\right)|}_s+{Df\left(x_i\right)|}_s\cdot \left(x_i-s\right)+𝒪\left({\Vert x_i-s\Vert }^2\right)\\ +\sigma \sum _{j=1}^N{a}_{ij}\{{g(x_j,x_i)|}_{(s,s)}+{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\\ \cdot \left(x_j-s\right)+𝒪\left({\Vert x_j-s\Vert }^2\right)+{\frac{\partial g(x_j,x_i)}{\partial x_i}|}_{(s,s)}\\ \cdot (x_i-s)+𝒪(\parallel x_i-s{\parallel }^2)\}-f(s)\\ \approx \cancel{f\left(s\right)}+Df\left(s\right)\delta x_i\\ +\sigma \sum _{j=1}^N{a}_{ij}\{\underset{\equiv 0}{\underbrace{g(s,s)}}+{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j\\ +{\frac{\partial g(x_j,x_i)}{\partial x_i}|}_{(s,s)}\delta x_i\}-\cancel{f\left(s\right)}\\ =Df\left(s\right)\delta x_i+\sigma \sum _{j=1}^N{a}_{ij}\{{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j\\ +{\frac{\partial g(x_j,x_i)}{\partial x_i}|}_{(s,s)}\delta x_i\},\end{array}$$

where $Df\left(s\right)$ denotes the $n\times n$ Jacobian matrix of $f$, evaluated at the synchronous state $s\left(t\right)$.

Now let us take an important conceptual step. Since the coupling function is diffusive, we have $g(x,x)=0$. Hence its value is constant (zero) on the synchronization manifold¹². It then immediately follows that its total derivative also vanishes, which in turn implies that

$${\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}+{\frac{\partial g(x_j,x_i)}{\partial x_i}|}_{(s,s)}=0,$$

which allows us to reorganize the linearized equation and write it in terms of the Laplacian,

$$\begin{array}{cc}\hfill {\dot{\delta x}}_i& =Df\left(s\right)\delta x_i+\sigma \sum _{j=1}^N{a}_{ij}\left\{{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j-{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_i\right\}\hfill \\ & =Df\left(s\right)\delta x_i+\sigma \left\{\sum _{j=1}^N{a}_{ij}{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j-\left(\sum _{j=1}^N{a}_{ij}\right){\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_i\right\}\hfill \\ & =Df\left(s\right)\delta x_i+\sigma \left\{\sum _{j=1}^N{a}_{ij}{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j-k_i{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_i\right\}\hfill \\ & =Df\left(s\right)\delta x_i+\sigma \left\{\sum _{j=1}^N{a}_{ij}{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j-\sum _{j=1}^N{\delta }_{ij}{k}_j{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j\right\}\hfill \\ & =Df\left(s\right)\delta x_i+\sigma \sum _{j=1}^N\left(a_{ij}-{\delta }_{ij}{k}_j\right){\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j\hfill \\ & =Df\left(s\right)\delta x_i-\sigma \sum _{j=1}^N\left({\delta }_{ij}{k}_j-a_{ij}\right){\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j\hfill \\ & =Df\left(s\right)\delta x_i-\sigma \sum _{j=1}^N{l}_{ij}{\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}\delta x_j.\hfill \end{array}$$

Furthermore, to simplify the notation, we introduce the $n\times n$ “pseudo-Jacobian” matrix of the function $g$, which we define as

$$Dg(s,s)=D_{x_j}g(s,s)={\frac{\partial g(x_j,x_i)}{\partial x_j}|}_{(s,s)}.$$

With this notation, we have:

(21) $${\dot{\delta x}}_i=Df\left(s\right)\delta x_i-\sigma \sum _{j=1}^N{l}_{ij}Dg(s,s)\delta x_j.$$

This last equation (21) corresponds to equation (6), with the substitution $Dh\left(s\right)\to Dg(s,s)$. Hence, we can follow the same algebraic steps as before. Indeed, equation (21) can be rewritten in block form using the Kronecker product. For this, we introduce the column vector of dimension $nN$, the perturbation vector $\delta X$ defined in equation (7). Stacking all $N$ coupled equations, but now in the form given in equation (21), each of dimension $n$, we obtain the following compact matrix equation:

(22) $$\dot{\delta X}=\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)-\sigma L\otimes Dg(s,s)\right)\delta X.$$

Once again, we can exploit both the spectral properties of the Laplacian matrix and the algebraic properties of the Kronecker product to decouple the system into independent equations along the Laplacian eigenvectors,

(23) $${\dot{\xi }}_i=\left(Df\left(s\right)-\sigma {\lambda }_{i}Dg(s,s)\right){\xi }_i,$$

where ${\xi }_i$ is the component of $\delta X$ along the $i$-th eigenvector. In other words, ${\xi }_i$ is the eigenmode previously defined, representing the projection of the perturbation $\delta X$ onto the direction of $v_i$, i.e.,

$$\delta X=\sum _{i=1}^N{v}_i\otimes {\xi }_i.$$

The MSF analysis proceeds in the same way. Each block in equation (23) has the same structure, differing only by the eigenvalue ${\lambda }_i$ of $L$. All of them can therefore be written in the parametric form

(24) $$\dot{\xi }=\left(Df\left(s\right)-rDg(s,s)\right)\xi .$$

As before, we can treat this generic equation and compute its largest Lyapunov exponent $\Lambda$. The dependence of $\Lambda$ on the generalized coupling strength $r=\sigma {\lambda }_i$ in equation (24) defines the MSF.

It is important to emphasize that, although the formulation obtained under general diffusive coupling reproduces the same results as the traditional approach for a coupling that is also diffusive, but in the specific form $h\left(x_j\right)-h\left(x_i\right)$, with the replacement $Dh\left(s\right)\to Dg(s,s)$ being sufficient, its significance lies precisely in extending the MSF framework. Recall that no assumption of linearity is made regarding the coupling function: the derivation only requires the coupling to be at least continuously differentiable, so that the linearization along the synchronization manifold is well defined. Indeed, this new approach extends the validity of the classical MSF formulation, allowing for the incorporation of more general interactions while preserving the analytical structure of the theory and thus covering a much broader class of coupling functions, e.g., the formal treatment of the Kuramoto model [²¹, ²², ¹⁸, ⁹¹], where $n=1$, in which each dynamical system or network unit $i$ is represented by the instantaneous phase ${\theta }_i$ of an oscillator, and the coupling between nodes $i$ and $j$ is given by the function $\sin \left({\theta }_j-{\theta }_i\right)$. This function does not have the specific classical form $h\left(x_j\right)-h\left(x_i\right)$, i.e., there is no function $h:ℝ\to ℝ$ such that $\sin \left({\theta }_j-{\theta }_i\right)=h\left({\theta }_j\right)-h\left({\theta }_i\right)$. Thus, this generalization not only confirms the validity of known results, but also extends the applicability of the MSF to a wider class of coupling schemes.

6. Conclusions

In this paper, we present the master stability function (MSF) approach to determine synchronizability in networks. The MSF formalism provides an elegant framework that separates local dynamics from network structure to evaluate the linear stability of the synchronization manifold. We first focus on the classic and original formulation by Pecora and Carroll for diffusive couplings in the specific algebraic form $h\left(x_j\right)-h\left(x_i\right)$. We have discussed how the MSF is obtained as the largest Lyapunov exponent transverse to the synchronization subspace, calculated from a simple set of linear equations as a function of the generic coupling parameter $r=\sigma {\lambda }_i$. A necessary (though not sufficient) condition for observing the synchronous state in practice is that the MSF must be negative. We also discuss how to characterize the possible MSF types for coupled chaotic systems, according to the regimes in which synchronization can emerge. The MSF is a fundamental and widely used tool for studying the synchronization of complex dynamical systems, e.g., chaotic systems.

Understanding the typical behaviors of the MSF for chaotic systems is essential for anticipating the collective dynamics of a network of such coupled systems. In this paper, we address this issue through numerical applications of the MSF formalism. We coupled representative nonlinear dynamical systems, widely studied in the chaos literature, under different coupling schemes between their components. From the results that we obtained, we present a generic classification of the MSF behaviors based on the number of finite crossings with the $r$-axis. Although our numerical applications focus on three representative chaotic systems, this classification is meant as a qualitative phenomenological taxonomy of MSF profiles, consistent with behaviors reported across a broad range of coupled nonlinear systems and coupling configurations in the literature. In this sense, the term “generic” refers to recurring qualitative scenarios in MSF curves, rather than to a universal claim based solely on the particular models analyzed.

We also show how this classification complements, and makes explicit, the characterization that has historically been implicitly assumed in the literature on network synchronization. The parity rule discussed provides an intuitive way to predict how small parameter changes can affect stability without altering the asymptotic behavior of the MSF. In this sense, the examples show that the MSF does not merely classify stability regimes, but also reflects how different coupling schemes act on the specific dynamical mechanisms of each system.

Beyond its pedagogical value and theoretical elegance, the MSF formalism also provides practical guidance for choosing and optimizing networks that are robustly synchronizable. In particular, it provides a reusable methodological “recipe” for setting up synchronizable networks in practice, once the node dynamics and the coupling scheme are specified. The classification of MSF profiles, based on the number of crossing points with the $r$-axis, provides a diagnostic tool for predicting whether a given combination of local dynamics, coupling scheme, and network structure will sustain stable synchronization. Once the local dynamics $f$ and the coupling function $h$ are given, the MSF determines the stability interval(s) of the synchronization manifold as a function of the effective coupling parameter $r$. Therefore, synchronizability can be promoted by shaping the Laplacian spectrum of the network so that all transverse modes $r=r_i=\sigma {\lambda }_i$ (for $i\ge 2$) fall within the region where $\Lambda \left(r\right)<0$. In particular, for type II MSFs, it is advantageous to increase ${\lambda }_2$, which reduces the critical coupling threshold ${\sigma }_c\sim r^{\ast }/{\lambda }_2$, while for type III MSFs, robustness typically requires a compact spectrum, often quantified by a small eigenratio ${\lambda }_N/{\lambda }_2$.

From an applied perspective, these results translate into a criterion for comparing topologies under fixed interaction restrictions: networks with higher values of ${\lambda }_2$ and/or lower values of ${\lambda }_N/{\lambda }_2$ favor the emergence of the synchronous state under weaker couplings and greater parameter uncertainty. Thus, the MSF formalism bridges mathematical stability analysis with concrete network design principles, showing how the “structure” must be adjusted to compensate for the intrinsic instability of nonlinear (e.g., chaotic) dynamics.

Finally, as a methodological extension, we explored a generalization of the classical MSF formulation for a general diffusive coupling, showing that simply replacing $Dh\left(s\right)$ with $Dg(s,s)$ preserves the analytical structure of the theory while allowing the incorporation of more general interactions. This approach encompasses a much broader class of coupling functions that do not fit the standard form $h\left(x_j\right)-h\left(x_i\right)$, confirming that the MSF remains a robust and versatile tool. The generalization of the MSF to general diffusive couplings further extends its applicability to systems with nonlinear or non-separable interactions, including phase-oscillator networks such as the Kuramoto model, a paradigm for collective synchronization in biological, chemical, and social systems. This flexibility allows us to assess synchronizability in settings where the traditional coupling assumptions of the classical MSF formulation do not hold.

It is important to emphasize, however, that the MSF formalism addresses the linear stability of the synchronization manifold and does not, by itself, guarantee that the synchronous state will be observed in practice. Even when the MSF is negative, indicating linear stability, the basin of attraction of the synchronized state may be very small, so that finite perturbations, parameter mismatch, or unfavorable initial conditions may prevent the system from reaching synchronization. In addition, the classical MSF framework assumes identical dynamical systems and the existence and invariance of the complete synchronization manifold, and therefore does not describe regimes such as partial synchronization, cluster states, or chimera states, as discussed in Subsection 2.5. These aspects do not reduce the value of the MSF as a stability analysis tool, but clarify its domain of validity and applicability and reinforce the importance of interpreting its predictions within the assumptions under which the formalism is constructed.

As possible continuations of the present work, one may explore adaptive synchronization schemes (adaptive gains or symmetry/structure control) as feedback strategies designed to steer the effective coupling into MSF-stable regions, going beyond the classical fixed-coupling setting considered here. Further continuations include extending the present discussion to discrete-time chaotic maps and other digital implementations, investigating how symmetry-related parameters may reshape the MSF profile, and exploring more recent synchronization paradigms such as time-reversible synchronization, which may require structural ingredients beyond the classical MSF setting.

We hope that this work will serve as both a useful reference and a starting point for future studies, whether applying MSF ideas to new coupling schemes, higher-order interactions, or emerging scenarios in the dynamics of complex networks. We also hope it can serve as a reusable methodological guide for setting up synchronizable networks in practice.

Acknowledgments

It is a pleasure to thank Leonardo L. Bosnardo and Sara F. J. Mion for carefully reading the manuscript and for helpful comments and suggestions. This work was partially supported by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), grant 88887.948035/2024-00 (CJG), the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), grant 2021/14335-0 (MAMA) and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant 303814/2023-3 (MAMA).

Supplementary Material

A. Some Basic Concepts of Spectral Graph Theory

In this appendix we review a few basic notions from spectral graph theory related to the adjacency and Laplacian matrices, as well as their main properties, which are important for understanding the stability analysis of synchronized states presented in the article. In particular, key synchronizability conditions discussed in Subsection 2.4 (§ 2.4.2 and § 2.4.3) are expressed in terms of the smallest nonzero and the largest Laplacian eigenvalues.

For this review, we mainly follow Van Mieghem’s book Graph Spectra for Complex Networks[⁹²] and the references contained therein as our primary reference and consultation source. However, for a broader overview of graph theory, we also recommend classic textbooks [⁹³, ⁹⁴, ⁹⁵] for learning the standard nomenclature and main results in this area. Furthermore, for an excellent didactic introduction to network science aimed at physicists, we suggest the reference [⁹⁶], as well as standard reference books in the field (e.g., [⁹⁷, ⁹⁸, ⁹⁹, ¹⁰⁰, ¹⁰¹]), to gain an insight into the importance and role of formal proofs and their computational analogues (i.e., the programs) as complementary tools for understanding.

Before presenting matrices and their properties, it is important to establish a basis of fundamental definitions. Thus, we will define a network as a graph $G=G(N,M)$, of order$N$ and size$M$, consisting of a set of $N$ nodes (or vertices) connected by a set of $M$ links (or edges). More generally, a graph, or a graphical representation of relationships within a set, is a mathematical structure that models pairwise relationships between objects. We usually refer to the topology of the graph (virtual “shape” or structure of the network) in terms of the pattern of interconnections between its elements.

In this context, we will also say that a graph is simple if its nodes have no self-connections, i.e., if there are no edges linking a node to itself. Furthermore, a graph is undirected if there is no distinction between the two vertices associated with each edge, i.e., the edges have no orientation. It should be noted that, in the MSF formalism, we restrict ourselves to simple, undirected networks (see assumptions in Subsection 2.1).

Finally, it is useful to introduce some additional concepts: we will call path the sequence of nodes (not repeated) connected in a graph, such that from each node in a path there is always a link to the next node in the sequence. The number of links in the path is called its length. The longest length among all shortest paths between any pair of nodes is called the diameter¹³. In other words, take all pairs of nodes, compute the shortest path length for each pair, and then take the maximum of these lengths. That maximum is the diameter. Thus, we say that a graph is connected if there is a path between any pair of nodes. In this case, the diameter is finite. And we will also say that a connected component of an undirected graph is defined as a subgraph, in which any pair of nodes is connected by a path.

Having established the above, we now introduce the matrix that encodes the topological information of a graph: the adjacency matrix. Indeed, a network can be represented in terms of its adjacency matrix $A$, a $\left(\mathrm{0,1}\right)$-valued (boolean) matrix, which describes the adjacency relationship between two vertices of a graph, i.e., whether there is ($=1$) or is not ($=0$) an edge connecting two vertices in that graph. Formally,

Definition 1 (Adjacency matrix)

Let $G=G(N,M)$ be a simple, undirected graph. The adjacency matrix of $G$ is the $N\times N$ matrix $A=\left(a_{ij}\right)$ defined by

(A.1) $$a_{ij}=\{\begin{array}{cc}1,\hfill & \mathit{\text{if vertices i and j are connected}},\hfill \\ 0,\hfill & \mathit{\text{otherwise}}.\hfill \end{array}$$

It is worth noting that several properties of a graph are directly reflected on the properties of its adjacency matrix and vice versa, which is a foundational premise of spectral graph theory [⁹²]. Some of these are immediate: for instance, $A$ is a real symmetric matrix (hence its eigenvalues are all real, and eigenvectors associated with distinct eigenvalues are orthogonal¹⁴). Moreover, all diagonal entries of $A$ are zero, reflecting the absence of self-connections, and a graph is uniquely determined by its adjacency matrix, which, as mentioned before, fully encodes the topological information (given an adjacency matrix $A$, one can reconstruct the graph $G$ associated with it and vice versa), among others. It should also be stressed that we restrict ourselves to simple, undirected networks.

Definition 2 (Degree of a vertex)

The degree $k_i$ of the $i$-th node (vertex) is the number of connections (edges) incident to it, i.e.,

(A.2) $$k_i=\sum _j{a}_{ij}.$$

New properties arise, e.g., if we sum the degrees of each node in a graph $G(N,M)$, each edge is counted twice, so ${\sum }_i{k}_i=2M$. Hence, the adjacency matrix provides not only information about the number of nodes or sites (by definition, this is the dimension of the matrix), but also about the number of links in a graph.

In addition to the adjacency matrix, another key matrix associated with a network is the combinatorial Laplacian, or simply the Laplacian, $L$. The eigenvalues and eigenvectors of matrices $A$ and $L$ provide a lot of information about the network structure [⁹²]. The eigenvalues of $L$, e.g., are closely related to how well the graph is connected [¹⁰², ¹⁰³, ¹⁰⁴] and how quickly a random walk spreads through it [¹⁰⁵]. In particular, the smallest nonzero eigenvalue of the Laplacian matrix (also known as the algebraic connectivity or Fiedler eigenvalue [¹⁰⁶]) plays a central role in shaping the synchronization properties of the network [¹⁰⁷, ¹⁰⁵, ¹⁰², ¹⁰⁸]. Moreover, since the graph is undirected, the matrix $L$ is symmetric (as we will see in the definition below), so its eigenvalues are real and it admits a complete set of orthonormal eigenvectors [¹⁰⁹].

Definition 3 (Laplacian matrix)

Given a simple, undirected graph $G(N,M)$, the combinatorial Laplacian matrix $L=\left(l_{ij}\right)$ is an $N\times N$ matrix with entries given by

(A.3) $$l_{ij}=\{\begin{array}{cc}{k}_i\mathit{\text{if i = j}},\hfill & \\ -1\mathit{\text{if i and j are connected}},\hfill & \\ 0\mathit{\text{otherwise}}.\hfill & \end{array}$$

Equivalently, the Laplacian can be written compactly as $L=D-A$, where $D=\mathrm{diag}(k_1,k_2,\dots ,k_N)$ is the degree matrix and $A$ is the adjacency matrix of the graph, i.e.,

(A.4) $$l_{ij}={\delta }_{ij}{k}_i-a_{ij},$$

where ${\delta }_{ij}$ is the Kronecker delta, which is ${\delta }_{ij}=1$ if $i=j$ and $0$ otherwise.

The following result summarizes some key properties of the Laplacian.

Definition 4 (spectral properties of the Laplacian)

Let $G(N,M)$ be a simple undirected graph and $L$ its associated Laplacian. Then:① $L$ has only real eigenvalues.② $\lambda =0$ is always an eigenvalue of $L$, with corresponding eigenvector $\left(\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right)$.③ $L$ is positive semidefinite, its eigenvalues listed in ascending order and repeated according to their multiplicity satisfy

(A.5) $$0={\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_N.$$

④ The multiplicity of the zero eigenvalue of $L$ is equal to the number of connected components of $G$.

For practical purposes, we will not provide a formal proof of the theorem (we suggest consulting Van Mieghem’s book Graph Spectra for Complex Networks[⁹²], § 2.1), but we briefly discuss its validity below.

Indeed, ① since $L$ is symmetric, all its eigenvalues are real and the respective eigenvectors (associated with different eigenvalues) are orthogonal, i.e., $L$ is orthogonally diagonalizable. ② Furthermore, it is clear that each row (or column) $i$ of $L$ has the degree of the $i$-th vertex $k_i$ on the diagonal and also has $-1$ repeated $k_i$ times (once for each neighbor), i.e., the sum of each row (or column) is zero. Consequently, $\lambda =0$ is an eigenvalue of $L$, associated with the eigenvector $v={(\mathrm{1,\; 1},\dots \mathrm{,\; 1})}^T$, satisfying $Lv=0$. This also implies that $L$ is singular, since $\det \left(L\right)=0$, and that there are linearly dependent columns, therefore, $\mathrm{rank}\left(L\right)<N$. ③ Note that the Laplacian matrix can be written as $L=R^{T}R$, where $R=\left(r_{li}\right)$ is the directed incidence matrix (even if the graph is undirected¹⁵), of dimension $M\times N$ with entries given by

$$r_{li}=\{\begin{array}{cc}+1\mathit{\text{if edge l enters vertex i}},\hfill & \\ -1\mathit{\text{if edge l leaves vertex i}},\hfill & \\ 0\mathit{\text{otherwise}}.\hfill & \end{array}$$

For each edge $l=(u,v)$, we assign an arbitrary orientation (e.g., from $u$ to $v$), such that $r_{lu}=-1$ (leaves from $u$) and $r_{lv}=+1$ (enters into $v$). Thus, given any eigenvector $v_i$ (already normalized) associated with the eigenvalue ${\lambda }_i$, we have

$${\lambda }_i=v_i^{T}L{v}_i=v_i^T{R}^{T}R{v}_i={\left(R{v}_i\right)}^T\left(R{v}_i\right)={\Vert R{v}_i\Vert }^2\ge 0,$$

i.e., the eigenvalues of $L$ are non-negative. Moreover, since we have already seen that $\lambda =0$ is always an eigenvalue, it is the smallest one, and furthermore, nothing prevents us from ordering them as in equation (A.5). ④ Finally, it is clear that in a graph with $n$ connected components ($n\le N$), the vertices can be reordered (relabeled) so that the adjacency matrix is block diagonal, where each block is the adjacency matrix of a connected component of the graph, with $A_m$ being the adjacency matrix of the $m$-th connected component among these $n$ components. The construction of the Laplacian preserves the block structure. Thus, each of the blocks of the Laplacian matrix is, in turn, a Laplacian matrix $L_m=D_m-A_m$ of a connected network and, therefore, each block has a zero eigenvalue with multiplicity one. Consequently, ${\lambda }_2>0$ whenever the network is connected.

On the other hand, the Laplacian spectrum is also related to other topological invariants. One of the most interesting is its relationship with the diameter, size, and degrees of the graph, as established in the following theorem.

Theorem 2

Let $G(N,M)$ be a simple, undirected graph with diameter $d$, and let $L$ be its associated Laplacian with eigenvalues given in equation (A.5). Then:

① [¹¹⁰]${\lambda }_2\ge \frac{4}{Nd}$,

② [¹⁰⁶]${\lambda }_2\le \frac{N}{N-1}{k}_{\min }$,

③ [¹⁰⁶]${\lambda }_N\ge \frac{N}{N-1}{k}_{\max }$,

④ [¹¹¹]${\lambda }_N\le \underset{(i,j)\mathrm{connected}}{\max }\left(k_i+k_j\right)$.

Here, $k_{\min }=\underset{i=1,\dots ,N}{\min }{k}_i$ and $k_{\max }=\underset{i=1,\dots ,N}{\max }{k}_i$ denote the smallest and largest degrees among all nodes of the graph, respectively. We will not present the proofs of this theorem here, but they can be found in the references we provide in it. Moreover, we also suggest consulting Section § 4.2 of reference [⁹²]. What we want to emphasize is that, for a fixed graph size, the magnitude of ${\lambda }_2$ reflects how well connected the graph is.

B. Kronecker Product

Let $A$ be a $p\times q$ matrix and $B$ an $r\times s$ matrix. The direct matrix product or Kronecker product (a special case of the tensor product of matrices) is the $pr\times qs$ block matrix given by:

$$A\otimes B=\left(\begin{array}{cccc}{a}_{11}B& a_{12}B& \cdots & a_{1q}B\\ a_{21}B& a_{22}B& \cdots & a_{2q}B\\ ⋮& ⋮& \ddots & ⋮\\ a_{p1}B& a_{p2}B& \cdots & a_{pq}B\end{array}\right).$$

For example,

$$\begin{array}{l}\left(\begin{array}{cc}1& 2\\ 3& 4\end{array}\right)\otimes \left(\begin{array}{cc}a& b\\ c& d\end{array}\right)=\left(\begin{array}{c}\begin{array}{cccc}a& b& 2a& 2b\\ c& d& 2c& 2d\\ & & & \\ 3a& 3b& 4a& 4b\\ 3c& 3d& 4c& 4d\end{array}\end{array}\right),\\ \left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\otimes \left(\begin{array}{cc}1& 2\\ 3& 4\end{array}\right)=\left(\begin{array}{c}\begin{array}{cccc}a& 2a& b& 2b\\ 3a& 4a& 3b& 4b\\ & & & \\ c& 2c& d& 2d\\ 3c& 4c& 3d& 4d\end{array}\end{array}\right).\end{array}$$

Note that in general, $A\otimes B\ne B\otimes A$.

The Kronecker product has several useful properties, among which we highlight the following: if $A$, $B$, $C$, and $D$ are matrices of such size that the products $AC$ and $BD$ are well defined, then

$$\left(A\otimes B\right)\left(C\otimes D\right)=\left(AC\right)\otimes \left(BD\right).$$

This relation is known as the mixed-product property because it links the ordinary matrix product with the Kronecker product in a single operation. A direct verification is straightforward: starting from the definition of the Kronecker product, one can expand both sides and check that the block structure matches. We sketch the argument below.

Suppose $A$ is $p\times q$ and $C$ is $q\times r$,

$$\begin{array}{cc}\hfill \left(A\otimes B\right)\left(C\otimes D\right)& =\left(\begin{array}{cccc}{a}_{11}B& a_{12}B& \cdots & a_{1q}B\\ a_{21}B& a_{22}B& \cdots & a_{2q}B\\ ⋮& ⋮& \ddots & ⋮\\ a_{p1}B& a_{p2}B& \cdots & a_{pq}B\end{array}\right)\left(\begin{array}{cccc}{c}_{11}D& c_{12}D& \cdots & c_{1r}D\\ c_{21}D& c_{22}D& \cdots & c_{2r}D\\ ⋮& ⋮& \ddots & ⋮\\ c_{q1}D& c_{q2}D& \cdots & c_{qr}D\end{array}\right)\hfill \\ & \stackrel{\underset{\text{①}}{}}{=}\left(\begin{array}{cccc}\left({\sum }_{s=1}^q{a}_{1s}{c}_{s1}\right)BD& \left({\sum }_{s=1}^q{a}_{1s}{c}_{s2}\right)BD& \cdots & \left({\sum }_{s=1}^q{a}_{1s}{c}_{sr}\right)BD\\ \left({\sum }_{s=1}^q{a}_{2s}{c}_{s1}\right)BD& \left({\sum }_{s=1}^q{a}_{2s}{c}_{s2}\right)BD& \cdots & \left({\sum }_{s=1}^q{a}_{2s}{c}_{sr}\right)BD\\ ⋮& ⋮& \ddots & ⋮\\ \left({\sum }_{s=1}^q{a}_{ps}{c}_{s1}\right)BD& \left({\sum }_{s=1}^q{a}_{ps}{c}_{s2}\right)BD& \cdots & \left({\sum }_{s=1}^q{a}_{ps}{c}_{sr}\right)BD\end{array}\right)\hfill \\ & \stackrel{\underset{\text{②}}{}}{=}\left(\begin{array}{cccc}{\left(AC\right)}_{11}BD& {\left(AC\right)}_{12}BD& \cdots & {\left(AC\right)}_{1r}BD\\ {\left(AC\right)}_{21}BD& {\left(AC\right)}_{22}BD& \cdots & {\left(AC\right)}_{2r}BD\\ ⋮& ⋮& \ddots & ⋮\\ {\left(AC\right)}_{p1}BD& {\left(AC\right)}_{p2}BD& \cdots & {\left(AC\right)}_{pr}BD\end{array}\right)\hfill \\ & =\left(AC\right)\otimes \left(BD\right),\hfill \end{array}$$

where, in step ①, we used the fact that the multiplication of two block matrices can be carried out as if their blocks were scalars, and in step ②, we applied the definition of matrix multiplication to deduce that

$${\left(AC\right)}_{uv}=\sum _{s=1}^q{a}_{us}{c}_{sv},$$

with ${\left(AC\right)}_{uv}$ denoting the $(u,v)$-th entry of $AC$.

An important consequence of this property is that if $A{v}_i={\lambda }_i{v}_i$ and $B{u}_j={\mu }_j{u}_j$, then

$$\left(A\otimes B\right)\left(v_i\otimes u_j\right)=\left(A{v}_i\right)\otimes \left(B{u}_j\right)=\left({\lambda }_i{v}_i\right)\otimes \left({\mu }_j{u}_j\right)={\lambda }_i{\mu }_j\left(v_i\otimes u_j\right).$$

In words, let $A$ and $B$ be square matrices of sizes $p\times p$ and $q\times q$, respectively. If ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p$ are the eigenvalues of $A$, associated with eigenvectors $v_1,v_2,\dots ,v_p$, and ${\mu }_1,{\mu }_2,\dots ,{\mu }_q$ are the eigenvalues of $B$, associated with $u_1,u_2,\dots ,u_q$ (listed according to their multiplicity), then, the eigenvalues of $A\otimes B$ are ${\lambda }_i{\mu }_j$, associated with the eigenvectors $v_i\otimes u_j$ of size $pq\times 1$, where $i=\mathrm{1,\; 2},\dots ,p$ and $j=\mathrm{1,\; 2},\dots ,q$.

For a more detailed and comprehensive review of the Kronecker product and its properties, we suggest the book Matrix Differential Calculus with Applications in Statistics and Econometrics by Jan R. Magnus and Heinz Neudecker [¹¹²], see § 2.2 and § 2.3, as well as the references contained therein. Moreover, we also suggest consulting the textbooks [¹¹³] and [¹¹⁴] for a practical and application-oriented presentation of the Kronecker product. Readers seeking a clear and self-contained exposition, including detailed proofs of standard properties, may also consult the online resource¹⁶[¹¹⁵]. Finally, a comprehensive list of properties can be found in the Ch. 4, § 4.2 of reference [¹¹⁶].

C. Correspondence Between Equations (6) and (8)

We start from equation (8) and explicitly expand its terms:

(A.6) $$\dot{\delta X}=\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)\right)\delta X-\sigma \left(L\otimes Dh\left(s\right)\right)\delta X.$$

We now analyze the two contributions on the right-hand side of equation (A.6) separately.

The first part corresponds to the term:

(A.7) $$\begin{array}{l}\left({\mathrm{𝟙}}_N\otimes Df\left(s\right)\right)\delta X=\left(\begin{array}{cccc}1& 0& \cdots & 0\\ 0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right)\otimes \stackrel{Df\left(s\right)}{\overbrace{\left(\begin{array}{cccc}{\partial }_1{f}^1\left(s\right)& {\partial }_2{f}^1\left(s\right)& \cdots & {\partial }_n{f}^1\left(s\right)\\ {\partial }_1{f}^2\left(s\right)& {\partial }_2{f}^2\left(s\right)& \cdots & {\partial }_n{f}^2\left(s\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\partial }_1{f}^n\left(s\right)& {\partial }_2{f}^n\left(s\right)& \cdots & {\partial }_n{f}^n\left(s\right)\end{array}\right)}}\left(\begin{array}{c}\delta x_1\\ \delta x_2\\ ⋮\\ \delta x_N\end{array}\right)\\ =\left(\begin{array}{c}\begin{array}{cccc}Df\left(s\right)& 0& \cdots & 0\\ & & & \\ 0& Df\left(s\right)& \cdots & 0\\ & & & \\ ⋮& ⋮& \ddots & ⋮\\ & & & \\ 0& 0& \cdots & Df\left(s\right)\end{array}\end{array}\right)\left(\begin{array}{c}\delta x_1\\ \delta x_2\\ ⋮\\ \delta x_N\end{array}\right)\\ =\left(\begin{array}{c}Df\left(s\right)\delta x_1\\ Df\left(s\right)\delta x_2\\ ⋮\\ Df\left(s\right)\delta x_N\end{array}\right).\end{array}$$

Here, ${\left(Df\left(s\right)\right)}_{ik}\equiv {\partial }_k{f}^i\left(s\right)={\frac{\partial f^i}{\partial x^k}|}_{x=s}$, where the superscript $i$ labels the component of the vector field and the subscript $k$ labels the component of the state vector. The same index convention applies to any Jacobian matrix appearing below, including $Dh\left(s\right)$.

The second part corresponds to the term:

(A.8) $$\begin{array}{cc}\hfill \left(L\otimes Dh\left(s\right)\right)\delta X& =\left(\begin{array}{cccc}{l}_{11}& l_{12}& \cdots & l_{1N}\\ l_{21}& l_{22}& \cdots & l_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ l_{N1}& l_{N2}& \cdots & l_{NN}\end{array}\right)\otimes \stackrel{Dh\left(s\right)}{\overbrace{\left(\begin{array}{cccc}{\partial }_1{h}^1\left(s\right)& {\partial }_2{h}^1\left(s\right)& \cdots & {\partial }_n{h}^1\left(s\right)\\ {\partial }_1{h}^2\left(s\right)& {\partial }_2{h}^2\left(s\right)& \cdots & {\partial }_n{h}^2\left(s\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\partial }_1{h}^n\left(s\right)& {\partial }_2{h}^n\left(s\right)& \cdots & {\partial }_n{h}^n\left(s\right)\end{array}\right)}}\left(\begin{array}{c}\delta x_1\\ \delta x_2\\ ⋮\\ \delta x_N\end{array}\right)\hfill \\ & =\left(\begin{array}{c}\begin{array}{cccc}{l}_{11}Dh\left(s\right)& l_{12}Dh\left(s\right)& \cdots & l_{1N}Dh\left(s\right)\\ & & & \\ l_{21}Dh\left(s\right)& l_{22}Dh\left(s\right)& \cdots & l_{2N}Dh\left(s\right)\\ & & & \\ ⋮& ⋮& \ddots & ⋮\\ & & & \\ l_{N1}Dh\left(s\right)& l_{N2}Dh\left(s\right)& \cdots & l_{NN}Dh\left(s\right)\end{array}\end{array}\right)\left(\begin{array}{c}\delta x_1\\ \delta x_2\\ ⋮\\ \delta x_N\end{array}\right)\hfill \\ & =\left(\begin{array}{c}{l}_{11}Dh\left(s\right)\delta x_1+l_{12}Dh\left(s\right)\delta x_2+\cdots +l_{1N}Dh\left(s\right)\delta x_N\\ l_{21}Dh\left(s\right)\delta x_1+l_{22}Dh\left(s\right)\delta x_2+\cdots +l_{2N}Dh\left(s\right)\delta x_N\\ ⋮\\ l_{N1}Dh\left(s\right)\delta x_1+l_{N2}Dh\left(s\right)\delta x_2+\cdots +l_{NN}Dh\left(s\right)\delta x_N\end{array}\right)\hfill \\ & =\left(\begin{array}{c}{\sum }_j{l}_{1j}Dh\left(s\right)\delta x_j\\ {\sum }_j{l}_{2j}Dh\left(s\right)\delta x_j\\ ⋮\\ {\sum }_j{l}_{Nj}Dh\left(s\right)\delta x_j\end{array}\right),\text{with}j=\mathrm{1,\; 2},\dots ,N\text{.}\hfill \end{array}$$

Substituting equations (A.7) and (A.8) into A.6, we obtain

$$\dot{\delta X}=\left(\begin{array}{c}Df\left(s\right)\delta x_1\\ Df\left(s\right)\delta x_2\\ ⋮\\ Df\left(s\right)\delta x_N\end{array}\right)-\sigma \left(\begin{array}{c}{\sum }_j{l}_{1j}Dh\left(s\right)\delta x_j\\ {\sum }_j{l}_{2j}Dh\left(s\right)\delta x_j\\ ⋮\\ {\sum }_j{l}_{Nj}Dh\left(s\right)\delta x_j\end{array}\right),$$

or, equivalently,

(A.9) $$\left(\begin{array}{c}{\dot{\delta x}}_1\\ {\dot{\delta x}}_2\\ ⋮\\ {\dot{\delta x}}_N\end{array}\right)=\left(\begin{array}{c}Df\left(s\right)\delta x_1-\sigma {\sum }_j{l}_{1j}Dh\left(s\right)\delta x_j\\ Df\left(s\right)\delta x_2-\sigma {\sum }_j{l}_{2j}Dh\left(s\right)\delta x_j\\ ⋮\\ Df\left(s\right)\delta x_N-\sigma {\sum }_j{l}_{Nj}Dh\left(s\right)\delta x_j\end{array}\right).$$

Thus, for the $i$-th component of equation (A.9), i.e., the $i$-th block of this vector, we recover exactly equation (6), corresponding to the $i$-th node. In other words, equation (8) already encodes, within its matrix structure, the full set of $N$ individual equations given by equation (6), completing the proof.

Data Availability

The computer codes and scripts used to generate the data and figures presented in this work are available in the GitHub repository at https://github.com/javierg0mez/MSF_examples.

Referências

• ^[1] A. Pikovsky, M. Rosenblum and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge University Press, Cambridge, 2001).
• ^[2] S. Boccaletti, A.N. Pisarchik, C.I. Del Genio and A. Amann, Synchronization: From Coupled Systems to Complex Networks (Cambridge University Press, Cambridge, 2018).
• ^[3] J. Pantaleone, American Journal of Physics 70, 992 (2002).
• ^[4] X. Wu, C. Zheng, Z. Lei, Y. Qian, Z. Di and X. Cui, Entropy 27, 908 (2025).
• ^[5] J.B. Buck, The Quarterly Review of Biology 13, 301 (1938).
• ^[6] J. Buck and E. Buck, Nature 211, 562 (1966).
• ^[7] J. Buck, The Quarterly Review of Biology 63, 265 (1988).
• ^[8] J. Peña Ramirez, L.A. Olvera, H. Nijmeijer and J. Alvarez, Scientific Reports 6, 23580 (2016).
• ^[9] C.M. Gray, Journal of Computational Neuroscience 1, 11 (1994).
• ^[10] E.M. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting (MIT Press, Cambridge, 2007).
• ^[11] S. Boccaletti, J. Kurths, G. Osipov, D.L. Valladares and C.S. Zhou, Physics Reports 366, 1 (2002).
• ^[12] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.U. Hwang, Physics Reports 424, 175 (2006).
• ^[13] L.M. Pecora, T.L. Carroll, G.A. Johnson, D.J. Mar and J.F. Heagy, Chaos: An Interdisciplinary Journal of Nonlinear Science 7, 520 (1997).
• ^[14] L.M. Pecora and T.L. Carroll, Physical Review Letters 80, 2109 (1998).
• ^[15] S.H. Strogatz and I. Stewart, Scientific American 269, 102 (1993).
• ^[16] L. Wijayasooriya, S. Saghafi, E. Khan and P. Sanaei, Physica D: Nonlinear Phenomena 471, 135071 (2025).
• ^[17] F. Dörfler, M. Chertkov and F. Bullo, Proceedings of the National Academy of Sciences of the United States of America 110, 2005 (2013).
• ^[18] J.A. Acebrón, L.L. Bonilla, C.J.P. Vicente, F. Ritort and R. Spigler, Reviews of Modern Physics 77, 137 (2005).
• ^[19] L. Barreira and Y.B. Pesin, Lyapunov Exponents and Smooth Ergodic Theory (American Mathematical Society, Providence, 2002).
• ^[20] L. Huang, Q. Chen, Y.C. Lai and L.M. Pecora, Physical Review E 80, 036204 (2009).
• ^[21] Y. Kuramoto, in: International Symposium on Mathematical Problems in Theoretical Physics, edited by H. Araki (Springer, Berlin, 1975), v. 39.
• ^[22] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer, Berlin, 1984).
• ^[23] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno and C. Zhou, Physics Reports 469, 93 (2008).
• ^[24] M. Ramasamy, S. Kumarasamy, S.K. Sampathkumar, A. Karthikeyan and K. Rajagopal, Complexity 2023, 6616560 (2023).
• ^[25] G. Strang, Introduction to Linear Algebra (SIAM, Philadelphia, 2022).
• ^[26] G. Strang, Linear Algebra and Its Applications (Cengage Learning, Boston, 2006).
• ^[27] S. Lang, Introduction to Linear Algebra (Springer, New York, 2012).
• ^[28] F.U. Coelho and M.L. Lourenço, Um Curso de Álgebra Linear (EdUSP, São Paulo, 2024).
• ^[29] G.F. Simmons, Differential Equations with Applications and Historical Notes (CRC Press, Boca Raton, 2016).
• ^[30] W.E. Boyce, R.C. DiPrima and D.B. Meade, Elementary Differential Equations and Boundary Value Problems (John Wiley & Sons, Hoboken, 2021).
• ^[31] S.H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering (Chapman and Hall/CRC, Boca Raton, 2024).
• ^[32] D. Jordan and P. Smith, Nonlinear Ordinary Differential Equations: An Introduction for Scientists and Engineers (Oxford University Press, Oxford, 2007).
• ^[33] T. Stankovski, T. Pereira, P.V. McClintock and A. Stefanovska, Reviews of Modern Physics 89, 045001 (2017).
• ^[34] A.A. de Castro Júnior, Curso de Equações Diferenciais Ordinárias (Instituto Nacional de Matemática Pura e Aplicada, Rio de Janeiro, 2009).
• ^[35] M. Viana and J. Espinar, Equações Diferenciais: Uma Abordagem de Sistemas Dinâmicos (Instituto Nacional de Matemática Pura e Aplicada, Rio de Janeiro, 2021).
• ^[36] J. Sotomayor, Lições de Equações Diferenciais Ordinárias (Instituto Nacional de Matemática Pura e Aplicada, Rio de Janeiro, 1979).
• ^[37] M.W. Hirsch, R.L. Devaney and S. Smale, Differential Equations, Dynamical Systems, and Linear Algebra (Academic Press, New York, 1974).
• ^[38] M. Barahona and L.M. Pecora, Physical Review Letters 89, 054101 (2002).
• ^[39] S. Boccaletti, D.U. Hwang, M. Chavez, A. Amann, J. Kurths and L.M. Pecora, Physical Review E 74, 016102 (2006).
• ^[40] M. Coraggio, P. De Lellis, S.J. Hogan and M. Di Bernardo, IEEE Control Systems Letters 2, 653(2018).
• ^[41] A. Bayani, F. Nazarimehr, S. Jafari, K. Kovalenko, G. Contreras-Aso, K. Alfaro-Bittner, R.J. Sánchez-García and S. Boccaletti, Nature Communications 15, 4955 (2024).
• ^[42] L.H.A. Monteiro, Sistemas Dinâmicos Complexos (LF Editorial, São Paulo, 2014).
• ^[43] S. Acharyya, P. Pradhan and C. Meena, arXiv:2412.19163 (2024).
• ^[44] M. Cattani, I.L. Caldas, S.L. Souza and K.C. Iarosz, Revista Brasileira de Ensino de Física 39, e1309 (2016).
• ^[45] L.V. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora and S. Boccaletti, Nature Communications 12, 1255 (2021).
• ^[46] P. Ji, T.K. Peron, F.A. Rodrigues and J. Kurths, Scientific Reports 4, 4783 (2014).
• ^[47] A. Buscarino, L.V. Gambuzza, M. Porfiri, L. Fortuna and M. Frasca, Scientific Reports 3, 2026 (2013).
• ^[48] S. Yu, J. Zhou and S. Guan, Chaos, Solitons & Fractals 164, 112607 (2022).
• ^[49] X. Wei, J. Emenheiser, X. Wu, J.A. Lu and R.M. D’Souza, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 013110 (2018).
• ^[50] X. Xi, S. Panahi, V.T. Pham, Z. Wang, S. Jafari and I. Hussain, Applied Mathematics and Computation 379, 125226 (2020).
• ^[51] F. Comellas and S. Gago, Journal of Physics A: Mathematical and Theoretical 40, 4483 (2007).
• ^[52] G. Chen, IEEE/CAA Journal of Automatica Sinica 9, 573 (2022).
• ^[53] Y. Song, X. Liu, D. Wang, P. Gao and M. Xue, Autonomous Intelligent Systems 4, 29 (2024).
• ^[54] Y. Deng, Z. Jia and F. Yang, Discrete Dynamics in Nature and Society 2020, 9143917 (2020).
• ^[55] K. Rajagopal, S. He, H. Natiq, A. Bayani, F. Nazarimehr and S. Jafari, Physics Letters A 516, 129637 (2024).
• ^[56] J.T. Lizier, F. Bauer, F.M. Atay and J. Jost, Proceedings of the National Academy of Sciences of the United States of America 120, e2303332120 (2023).
• ^[57] S. Jafari, A. Bayani, F. Parastesh, K. Rajagopal, C.I. del Genio, L. Minati and S. Boccaletti, Physical Review Research 6, 043105 (2024).
• ^[58] L.M. Pecora, F. Sorrentino, A.M. Hagerstrom, T.E. Murphy and R. Roy, Nature Communications 5, 4079 (2014).
• ^[59] Y. Wang, D. Zhang, L. Wang, Q. Li, H. Cao and X. Wang, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 093139 (2022).
• ^[60] S. Majhi, B.K. Bera, D. Ghosh and M. Perc, Physics of Life Reviews 28, 100 (2019).
• ^[61] Y. Zhang and A.E. Motter, SIAM Review 62, 817 (2020).
• ^[62] N. Zabzina, A. Dussutour, R.P. Mann, D.J.T. Sumpter and S.C. Nicolis, PLOS Computational Biology 10, e1003960 (2014).
• ^[63] E.N. Lorenz, Journal of the Atmospheric Sciences 20, 130 (1963).
• ^[64] A. Stefański, P. Perlikowski and T. Kapitaniak, Physical Review E 75, 016210 (2007).
• ^[65] R. Mettin, U. Parlitz and W. Lauterborn, International Journal of Bifurcation and Chaos 3, 1529 (1993).
• ^[66] M. Porfiri, Europhysics Letters 96, 40014 (2011).
• ^[67] M. Ramasamy, K. Rajagopal, B. Ramakrishnan and A. Karthikeyan, Physica A: Statistical Mechanics and Its Applications 625, 129032 (2023).
• ^[68] D. Joseph, S. Kumarasamy, S.A. Jose and K. Rajagopal, Cognitive Neurodynamics 18, 4089 (2024).
• ^[69] M. Lodi, S. Panahi, F. Sorrentino, A. Torcini and M. Storace, Communications Physics 7, 198 (2024).
• ^[70] P. Lellis, M. Bernardo and F. Garofalo, Automatica 45, 1312 (2009).
• ^[71] F. Sorrentino, Physical Review E 80, 056206 (2009).
• ^[72] Q. Ren and J. Zhao, Physical Review E 76, 016207 (2007).
• ^[73] A. Tutueva, L. Moysis, V. Rybin, A. Zubarev, C. Volos and D. Butusov, Chaos, Solitons & Fractals 159, 112181 (2022).
• ^[74] M.S. Anwar, S.N. Jenifer, P. Muruganandam, D. Ghosh and T. Carletti, Physical Review E 110, 064305(2024).
• ^[75] A. Karimov, V. Rybin, I. Babkin, T. Karimov, V. Ponomareva and D. Butusov, Mathematics 13, 1437 (2025).
• ^[76] D. Butusov, V. Rybin and A. Karimov, Physical Review E 111, 014213 (2025).
• ^[77] E. Steur, W. Michiels, H. Huijberts and H. Nijmeijer, Physica D: Nonlinear Phenomena 277, 22 (2014).
• ^[78] A. Tiwari, R. Nathasarma and B.K. Roy, Journal of the Franklin Institute 361, 106637 (2024).
• ^[79] G. Vidal and H. Mancini, Proceedings of the Net-Works 2008 Conference (University of Navarra, Pamplona, 2008).
• ^[80] J. Sun, E.M. Bollt and T. Nishikawa, Europhysics Letters 85, 60011 (2009).
• ^[81] D. Hu and H. Cao, Communications in Nonlinear Science and Numerical Simulation 35, 105 (2016).
• ^[82] A.V. Tutueva, L. Moysis, V.G. Rybin, E.E. Kopets, C. Volos and D.N. Butusov, Chaos, Solitons & Fractals 155, 111732 (2022).
• ^[83] R.M. Corless, C. Essex and M.A.H. Nerenberg, Physics Letters A 157, 27 (1991).
• ^[84] R.M. Corless, Computers & Mathematics with Applications 28, 107 (1994).
• ^[85] E. Hairer, G. Wanner and C. Lubich, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (Springer, Berlin, 2006).
• ^[86] F. Bagnoli and M. Baia, Algorithms 16, 213 (2023).
• ^[87] P. Perlikowski and A. Stefański, Mechanics and Mechanical Engineering 10, 110 (2006).
• ^[88] A. Dabrowski, K. Mnich and A. Stefański, Mechanics and Mechanical Engineering 19, 161 (2015).
• ^[89] U. Uriostegui-Legorreta and E.S. Tututi, Journal of Applied Research and Technology 21, 227 (2023).
• ^[90] Z. Dayani, F. Parastesh, S. Jafari, E. Schöll, J. Kurths and J.C. Sprott, International Journal of Bifurcation and Chaos 33, 2350122 (2023).
• ^[91] L.L. Bosnardo, Synchronization of Kuramoto Oscillators in Modular Networks Doctoral Thesis, Universidade Estadual de Campinas, Campinas (2025).
• ^[92] P.V. Mieghem, Graph Spectra for Complex Networks (Cambridge University Press, Cambridge, 2023).
• ^[93] R. Diestel, Graph Theory (Springer Nature, Berlin, 2025).
• ^[94] J.A. Bondy and U.S.R. Murty, Graph Theory with Applications (Macmillan, London, 1976).
• ^[95] D.B. West, Introduction to Graph Theory (Prentice Hall, Upper Saddle River, 2001).
• ^[96] P.F. Gomes, Revista Brasileira de Ensino de Física 46, e20240190 (2024).
• ^[97] J. Reichardt, Structure in Complex Networks (Springer, Berlin, 2008).
• ^[98] M. Newman, A.L. Barabási and D.J. Watts, The Structure and Dynamics of Networks (Princeton University Press, Princeton, 2011).
• ^[99] M. Newman, Networks (Oxford University Press, Oxford, 2018).
• ^[100] E. Estrada, The Structure of Complex Networks: Theory and Applications (American Chemical Society, Washington, DC, 2012).
• ^[101] A. Zinilli, Elements of Network Science: Theory, Methods and Applications in Stata, R and Python (Springer Nature, Cham, 2025).
• ^[102] T. Pereira, arXiv:1112.2297 (2011).
• ^[103] J.D. Cook, Connectivity and the Graph Laplacian, available in: https://www.johndcook.com/blog/2016/01/07/connectivity-graph-laplacian/
  » https://www.johndcook.com/blog/2016/01/07/connectivity-graph-laplacian/
• ^[104] Z. Lin, J. Wang and M. Cai, arXiv:2302.10491 (2023).
• ^[105] S. Hata and H. Nakao, Scientific Reports 7, 1121 (2017).
• ^[106] M. Fiedler, Czechoslovak Mathematical Journal 23, 298 (1973).
• ^[107] Z.S. Duan, W.X. Wang, L. Chao and G.R. Chen, Chinese Physics B 18, 3122 (2009).
• ^[108] M. Zou and W. Guo, IEEE Systems Journal 17, 2145 (2022).
• ^[109] J.C. Bronski, L. DeVille and T. Ferguson, SIAM Journal on Applied Mathematics 76, 1126 (2016).
• ^[110] B. Mohar, Graphs and Combinatorics 7, 53 (1991).
• ^[111] W.N. Anderson Jr. and T.D. Morley, Linear and Multilinear Algebra 18, 141 (1985).
• ^[112] J.R. Magnus and H. Neudecker, Matrix Differential Calculus with Applications in Statistics and Econometrics (John Wiley & Sons, Chichester, 2019).
• ^[113] A. Graham, Kronecker Products and Matrix Calculus with Applications (Courier Dover Publications, Mineola, 2018).
• ^[114] Y. Hardy and W.H. Steeb, Matrix Calculus, Kronecker Product and Tensor Product: A Practical Approach to Linear Algebra, Multilinear Algebra and Tensor Calculus with Software Implementations (World Scientific, Singapore, 2019).
• ^[115] M. Taboga, Properties of the Kronecker Product, available in: https://www.statlect.com/matri56tyx-algebra/Kronecker-product-properties
  » https://www.statlect.com/matri56tyx-algebra/Kronecker-product-properties
• ^[116] R.A. Horn and C.R. Johnson, Topics in Matrix Analysis (Cambridge University Press, Cambridge, 1994).

• 1
  https://demoweb.physics.ucla.edu/content/160-spontaneous-synchronization 160. Spontaneous Synchronization — UCLA Physics & Astronomy (University of California, 2010). Click the icon to watch a demonstration video.
• 2
  The matrix that encodes the topological information of the network.
• 3
  The local existence and uniqueness of solutions under local Lipschitz continuity follow from the classical Picard-Lindelöf theorem. A rigorous and accessible discussion can be found, e.g., in § 2.1 of the IMPA lecture notes [³⁴]. For complementary treatments from the perspective of dynamical systems, see also standard textbooks such as [³⁵, ³⁶, ³⁷]. A detailed analysis of these foundational results lies beyond the scope of the present work.
• 4
  See Appendix B for the definition of the Kronecker product, the mixed-product property, and its main algebraic consequences.
• 5
  For a pedagogical and accessible review of Lyapunov exponents and the foundations of deterministic chaos, see [⁴⁴], especially $\S$ 6 of that work.
• 6
  This is the genius of Pecora and Carroll’s MSF.
• 7
  See Appendix A (Theorem 2) for bounds relating ${\lambda }_2$ to topological features such as the diameter and degrees of the graph.
• 8
  We are interested precisely in the synchronization of chaotic systems.
• 9
  Click on the icon to access the GitHub repository.
• 10
  It cannot be globally robust, since the attractor itself may be completely destroyed if a parameter changes significantly.
• 11
  With “type I” being the zero-crossing case.
• 12
  Note that the antisymmetry condition $g(x,y)=-g(y,x)$ is no longer required, it is sufficient that $g(x,x)=0$.
• 13
  If there is an isolated node, i.e., a node with no connections, we say that the diameter is infinite, $d=\infty$.
• 14
  We know that a real symmetric matrix is simply a special case of a Hermitian matrix. Therefore, both $A$ and $L$ (to be introduced later) inherit all the properties of a Hermitian matrix.
• 15
  This means assigning an arbitrary orientation to each edge.
• 16
  https://www.statlect.com/matrix-algebra/Kronecker-product-properties Properties of the Kronecker product — Marco Taboga (Statlect collection, 2021) — Lectures on matrix algebra. Click the link icon to go to resource.

Edited by

• Editor-in-Chief:
  Marcello Ferreira https://orcid.org/0000-0003-4945-3169

Publication Dates

• Publication in this collection
  01 June 2026
• Date of issue
  2026

History

• Received
  26 Oct 2025
• Reviewed
  19 Feb 2026
• Accepted
  29 Mar 2026