Synchronization in Networks with Strongly Delayed Couplings

We investigate the stability of synchronization in networks of dynamical systems with strongly delayed connections. We obtain strict conditions for synchronization of periodic and equilibrium solutions. In particular, we show the existence of a critical coupling strength $\kappa_{c}$, depending only on the network structure, isolated dynamics and coupling function, such that for large delay and coupling strength $\kappa<\kappa_{c}$, the network possesses stable synchronization. The critical coupling $\kappa_{c}$ can be chosen independently of the delay for the case of equilibria, while for the periodic solution, $\kappa_{c}$ depends essentially on the delay and vanishes as the delay increases. We observe that, for random networks, the synchronization interval is maximal when the network is close to the connectivity threshold. We also derive scaling of the coupling parameter that allows for a synchronization of large networks for different network topologies.

1. Introduction. Coupled systems and networks with time-delayed interactions emerge in different fields including laser dynamics [21,32,11,10,36,40], neural networks [39,15,20,7], traffic systems [28], and others. Time delays in these systems are caused by finite signal propagation or finite reaction times. In many cases, especially in optoelectronics, delays are larger than the other time scales of the system, and they play a major role for system's dynamics [36,40]. For neural systems, large delays may emerge as a lump effect of a signal propagation along feed-forward chains of neurons [24].
When the interaction among oscillators is diffusive, the network possesses a synchronous subspace, where all subsystems constituting the network behave identically. The stability of the synchronous subspace plays a major role in applications. For instance, neural synchronization is known to be involved in brain functioning [34], synchronization of lasers is used for communication purposes [8,2].
Therefore, the stability of synchronization for strongly self-feedback delayed interactions has been addressed in various contexts and generalities in a mixture of analytical and numerical techniques [22,13].
Here, we aim at studying how stability of the synchronization manifold depends on the dynamics of the isolated nodes and network structure. We focus on the case where the isolated dynamics is either periodic or an equilibrium solution. This case is amenable to detailed understanding in the limit of large delays.
Model and discussion of results: We consider a network model with a selffeedback term in the coupling. This kind of coupling is widely studied, for example, in semiconductors lasers models [36,18], where the self-feedback interaction is generated by the optical feedback. We consider the following system of n identical oscillators with time-delayed interactions [18,13,22,37,9,7] x j (t) = f (x j (t)) + κ n l=1 A jl h(x l (t − τ ) − x j (t − τ )), (1.1) where j = 1, · · · , n, κ is the coupling strength, A is the adjacency matrix with the elements A jl = 1 if j = l and there is a link from l to j and A jl = 0 otherwise. The function f : R q → R q describes the uncoupled dynamics and the coupling function h : R q → R q , h(0) = 0, describes interactions between the nodes. More details and assumptions on the model (1.1) is given in Sec. 2. We consider the time-delay τ > 0 as a large parameter. Physically, this means that the interaction time is much larger than the typical time scale of an uncoupled system, which often occurs in, e.g. coupled laser systems [36,40].
The uncoupled units of the network have an attractor, which is either an equilibrium or periodic orbit. And, we obtain conditions for synchronization and desynchronization which depends on the network dynamics.
(i ) For the equilibrium, we show that with strong enough delay, there is a critical coupling parameter κ c , which depends on the dynamics, coupling function and network structure, such that the equilibrium solution in the network remains stable if and only if |κ| < κ c (see Theorem 3.1). For this case, κ c is independent of the delay. (ii ) Synchronization of stable periodic orbits is in sharp contrast to fixed point case. It is shown to be always desynchronized for long enough time-delay. However, for long but finite delay the synchronization of periodic orbits are attained for an interval that is shrinking as the delay grows to infinity, that is, for this case the critical coupling parameter depends also on the delay, κ c = κ(τ ), and κ(τ ) → 0 as τ → ∞. See Theorem 5.1 for the precise statement.
Fixing the vector field f , the coupling function h and the delay τ , we can study how the network structure affects the critical coupling κ c . Our results elucidate this dependence. For instance, we will show that -For Barabási-Albert scale free network, κ c ∼ O(1/ √ n). -For Erdös-Réniy (ER) random networks, κ c ∼ O(1/ log n).
This shows that having a large number of connections is detrimental, and close to percolation threshold is optimal for synchronization. This is in contrast to the non-delayed case [25,29] -best synchronization scenario happens for homogeneous networks and the large the degrees the better.
Here we use the theory of functional differential equations [17,35], elements of the spectral theory of graphs [12,3,6], as well as a spectral theory for delay differential equations with strong delays [38,23,33,40]. Some important ideas from Ref. [13] on the synchronization of strongly delayed networks are used as well.
The manuscript is organized as follows: In Sec. 2, we give some assumption on the considered model. In Sec. 3 we state the main results, in particular, a special case for the result of synchronization of periodic orbits. In Section 4 we give a description of the spectrum for the variational equation and prove the main results with illustrations given in Sec. 6. In Sec. 7 we discuss the relevance of the structure of the network to the synchronization condition and in Sec. 7.1 we discuss how scaling properties of the coupling parameter allow synchronization.

2.
Assumptions on the network model. The initial conditions of (1.1) are history functions φ j : [−τ, 0] → R q and we assume that φ j is continuous for all j = 1, · · · , n. The state space where the solutions x j and initial conditions φ j lives is the space of continuous functions The function f : R q → R q is of class C r , r ≥ 2. Assume that the solutions of the system are uniformly bounded, i.e., the overall system is dissipative. So the solution operator is eventually compact and thus we can apply the spectral mapping theorem. Thus, we can obtain detailed information by studying the linearized system.
The following assumptions for the vector field f restrict the model (1.1) to the synchronization scenarios that we are interested in.
Assumption 2.1. The uncoupled system possesses an exponentially locally stable equilibrium s(t) = s 0 .
Assumption 2.2. The uncoupled system possesses a locally exponentially stable periodic solution s(t).
The coupling function h : R q → R q , is assumed to be bounded and differentiable, H = Dh(0) is the Jacobian matrix at 0 and det(H) = 0.
Due to the diffusive nature of the coupling, the coupling term vanishes identically if the states of all oscillators are identical. This ensures that the synchronized state x j (t) = s(t) for all j = 1, 2, . . . , n persists for all coupling strengths κ. Let s(t) be the solution specified in Assumptions 2.1 or 2.2 and U ⊂ R q a tubular neighborhood which belongs to the basin of attraction of s(t). The set S := {x j ∈ U ⊂ R q : j = 1, · · · , n and x 1 = · · · = x n } is called the synchronization manifold. The local stability of S determines the stability of the synchronization, and it depends, in particular, on spectral properties of Laplacian matrix L defined as L = D − A, where A is the adjacency matrix and D is the diagonal matrix of the network's degrees. We assume that the matrix L is diagonalizable.
3. Main Results. If the uncoupled system has stable equilibrium, then the coupled system possesses the stable equilibrium (s 0 , . . . , s 0 ) at κ = 0, and the equilibrium solution remains stable for small κ and small delay by the principle of linearization.
Our result concerns the limit of large delay. In this setting, the principle of linearization does not apply, and the spectral theory of DDE offers the tools to tackle this problem.
Theorem 3.1, gives the interval of the coupling parameter κ where the equilibrium solution of (1.1) is locally exponentially stable for large delays. More specifically, it provides a critical coupling parameter κ c such that the destabilization occurs for κ > κ c . Moreover, κ c does not depend on the delay. and τ > τ 0 , the equilibrium solution remains locally exponentially stable. Here, ρ L is the spectral radius of the Laplacian matrix L.
If the condition κ > κ c is fulfilled, then there exists such τ 0 > 0 that the equilibrium solution is unstable for τ > τ 0 .
As follows from Theorem 3.1, the condition κ = κ c is the strict destabilization value for τ → ∞, while the value κ c does not depend on time delay τ . The proof of Theorem 3.1 is presented in Sec. 4. The interval (0, r 0 /ρ L ) of the coupling strength will be called the synchronization window.
We also show, in Sec. 6.3, that the characteristic time t tr in which the trajectories of (1.1) approach the synchronous (equilibrium) solution scales as for κ κ c and large delay. In Sec. 6.3 we derive the estimate (3.2) and illustrate it with an example.
In the case of the synchronous periodic solutions, the situation is subtle. In particular, one can show that the periodic orbits will be desynchronized with an increasing delay, however, it is possible to synchronize them for any finite time-delay.
The generic condition for synchronization of periodic orbits is a co-dimension one condition. This condition is technical and it will be discussed at length later in Sec. 5.2. To avoid overly technical statements we first give our result for the special case when the isolated dynamics is given by a Stuart-Landau (SL) system. The network model readsż where z j ∈ C. The uncoupled SL system is retrieved if κ = 0. Note that it has one equilibrium at the origin, which is asymptotically stable for α < 0, and a stable periodic orbit for α > 0, which emerges from the Hopf bifurcation at α = 0. This model is illustrative since it represents a normal form for a Hopf bifurcation.
be fulfilled. Then, there is τ 0 > 0 such that for each τ > τ 0 there is κ c = κ c (τ ) such that the synchronous periodic solution is locally exponentially stable The Fig. 6.5 will give the numerically computed stability region from Theorem 3.2. The full statement will be given Sec. 5.2, Theorem 5.1.
The main consequences of the theorems 3.1 and 3.2 for different networks are presented as corollaries below. From the point of view of the network structure, the stability of the synchronization manifold is related to the spectral radius of the Laplacian matrix ρ L . Therefore, fixing f , h and τ > 0 large enough, one can study the effect of the changes in the network to the synchronization.
We give the asymptotic behavior of the synchronization window for two important examples of complex networks, namely, heterogeneous (e.g. Barabási-Albert (BA) scale-free network), and homogeneous network (e.g., Erdös-Réniy (ER) random graph). More details about BA and ER networks and the proofs of the corollaries that follow are given in Sec. 7.
Corollary 3.1 (Synchronization window for large BA networks). Consider the model (1.1) and its assumptions 3.1 with a BA network with n nodes. Then, for sufficiently large n, the length of the synchronization window scales as 1/ √ n.
Corollary 3.2 (Synchronization window for large ER networks). Consider the model (1.1) and its assumptions with a connected ER network with n nodes. Then, for sufficiently large n, the length of the synchronization window scales as 1/ ln n. Furthermore, the synchronization window is maximal near the percolation threshold.

Preliminaries: variational equation and spectral theory.
In this section, we introduce the relevant elements to proof our main results.

4.1.
Dynamics near the synchronization manifold. In order to study the local stability of the synchronization manifold we linearize (1. where A lj = g l δ lj − L lj with g l the degree of the node l and δ lj = 1 if l = j and 0 otherwise. Here, J(t) = Df (s(t)) is the Jacobian matrix of the vector field f along s(t) and H = Dh(0) is the Jacobian of h at 0. Using η = η T 1 , · · · , η T n we rewrite (4.1) aṡ where ⊗ is the Kronecker product and I is the identity matrix. As L is diagonalizable we use L = RΛR −1 , with Λ = diag(µ 1 , · · · , µ n ), and the change of coordinates η = (R −1 ⊗ I)ξ to get the block diagonal equatioṅ Then, each block of (4.2) is given bẏ Note that µ 1 = 0, and the variational equation corresponding to µ 1 iṡ which describes the perturbations within the synchronization manifold. So, the interest concentrates now on the spectrum of (4.3) for µ 2 , · · · , µ n . It is shown in [42,33,23,40] that the spectrum of (4.3) with large delay consists of two parts which we introduce here following Ref. [33]. We consider the following two cases.
-Case 1: the synchronized solution s(t) = s 0 is an equilibrium of the uncoupled system and then J(t) = J does not depend on t.
-Case 2: s(t) is a T −periodic solution of the uncoupled system, implying that

Spectral Theory.
Omitting the index j and denoting σ = −κµ, Eq. (4.3) can be rewritten asξ In a general case, the Floquet theory can be used to study the stability of (4.4) [17]. Using the following Floquet-like ansatz where y(t) is a non-trivial T -periodic function, λ and e λT are Floquet exponent and Floquet multipliers respectively, one can obtain a delay differential equation As the function y(t) is T periodic we can rewrite y(t − τ ) in terms of a fraction of the period T so that the large parameter τ appears only as a parameter in e −λτ in Eq. (4.6): where η = τ mod T . The system (4.7) is said to have a non-trivial periodic solution y(t) = y(t + T ) when it posses a Floquet multiplier equals to 1. Let U be the monodromy operator of (4.7). Then, U = U (λ, σe −λτ ) and y(t) is an eigenvector of U associated with the eigenvalue 1 (trivial Floquet multiplier), that is, [U (λ, σe −λτ ) − I]y(t) = 0. This leads to a characteristic equation [33] H(λ, σe −λτ ) = 0. (4.8) The solutions λ of H forms the whole spectrum of (4.4). In the case of s(t) = s 0 (equilibrium solution) then We will define now certain objects called "instantaneous spectrum", "strongly unstable spectrum" as well as "asymptotic continuous spectrum". Strictly speaking, these objects do not belong to the spectrum of the synchronous states, i.e., they are not the Lyapunov exponents of system (4.4). However, they play an important role since the spectrum will be well approximated by the "strongly unstable spectrum" and "asymptotic continuous spectrum" as the delay becomes large. Another advantage of these limiting spectra is that they can be much more easily found or computed in comparison to the actual spectrum of system (4.4), see more details in [42,41,38,23,33,40].
has a non-trivial periodic solution y(t) = y(t+T ) for some φ ∈ R. Here η = τ mod T .
Remark 4.2. If J(t) does not depend on time, the asymptotic continuous spectrum can be determined from the following equation As follows from [42,41,38,19,23,33,40], the spectrum of generic linear delay system (4.4) converges to either the instantaneous spectrum or to the curves in the complex plane γ(ω)/τ + iω, ω ∈ R, where γ(ω) is defined from the asymptotic continuous spectrum (note the division by τ ). The second part of the spectraconsisting of eigenvalues with asymptotically vanishing real parts -approaches a continuous curves asymptotically, while being still discrete for any finite τ ; for this reason, it is called pseudo-continuous spectrum.
An example of a typical spectrum of the system with long delay is shown in Fig. 6.2. In particular, the distances between neighboring eigenvalues within one curve of the pseudo-continuous spectrum scale as 2π/τ for large delay, and in the limit of τ → ∞ they vanish and the eigenvalues fill the curve [33].
The union of the strongly unstable spectrum, Γ SU , with the asymptotic continuous spectrum Γ A forms the approximation of the whole spectrum of Eq. (4.3), given that the instantaneous spectrum does not contain eigenvalues with zero real parts and some non-degeneracy conditions are fulfilled [23,33]. Moreover, if the set of the strongly unstable spectrum is empty (Γ SU = ∅), and the asymptotic continuous spectrum is entirely contained on the left side of the complex plane (Γ A ⊂ C − ) then the trivial solution of Eq. (4.3) is exponentially asymptotically stable.

5.
Proof of the main theorems. In this section we use the given description of the spectrum of a general linear DDE to prove Theorems 3.1 (persistence of equilibrium solution), and 5.1 (synchronization of periodic orbit).
The result obtained in Theorem 3.1 could also be proved by using Lyapunov-Kravoskii functionals [16]. However, we use the approach of studying the spectrum for large delay not only because it is powerful, but it also serves as a preparation to prove the result of synchronization of periodic solutions. 5.1. Persistence of equilibrium solution: Theorem 3.1. In this section we specify the stability condition for the case when system (1.1) possesses a synchronous solution s 0 , i.e. f (s 0 ) = 0, and the Jacobian J = Df (s 0 ) is a constant q × q matrix. The corresponding characteristic equation is The real parts of the eigenvalues λ determine the stability, and, as it was discussed above, the whole spectrum converges to either the strongly unstable spectrum or pseudo-continuous spectrum. Lemma 5.1 gives explicit dependence of the asymptotic continuous spectrum on the coupling strength σ.
Lemma 5.1 (Ref. [42]). Consider the linear delay differential equatioṅ with J, H ∈ R q×q and det(H) = 0. Then, the asymptotic continuous spectrum for (5.2) is given by the branches where l = 1, · · · , q and g l are complex roots of the polynomial The application of Lemma 5.1 to system (4.3) leads to the proof of Theorem 3.1.
Proof of Theorem 3.1 (Persistence of equilibrium solution). The basic idea of the proof is to show that the strongly unstable spectrum is empty and the asymptotic continuous spectrum is stable if and only if the condition (3.1) is satisfied. The rest of the proof elaborates this idea more carefully in details. The instantaneous spectrum Γ I coincides with the spectrum of J and, by the Assumption 2.1, the equilibrium solution of the isolated system is asymptotically stable, thus there are no eigenvalues of J with positive real parts, hence Γ SU = ∅.
Using Lemma 5.1 for the variational Eq. (4.3). Then, the real part of the asymptotic continuous spectrum is The condition γ l,j (ω, σ) < 0 is fulfilled if and only if The asymptotic continuous spectrum lies strictly on the left side of the complex plane if the inequality (5.6) holds for all j = 2, · · · , n, l = 1, · · · , q, and ω ∈ R. This leads to the condition where r 0 = min l=1,··· ,q inf ω∈R |g l (ω)| exists and is always bounded from zero. Indeed, if we assume the opposite, i.e., r 0 = 0, then it means that there exist such l 0 SYNCHRONIZATION IN NETWORKS WITH STRONGLY DELAYED COUPLINGS 9 and ω 0 that g l0 (ω 0 ) = 0, and hence, det(−iω 0 + J) = 0 implying that iω 0 belongs to the spectrum of J. In this way we arrive at the contradiction to the assumption that the spectrum of J has strictly negative real parts. Moreover, r 0 = r 0 (J, H) since the functions g l (ω) are the roots of the characteristic equation (5.4). The number ρ L = max j=2,··· ,n |µ j | stands for the spectral radius of the Laplacian matrix L. Therefore, under the condition (3.1), the asymptotic continuous spectrum Γ A ⊂ C − and therefore the zero solution of (4.3) is stable.
Having both Re Γ A < 0 and Γ SU = ∅, Theorem 5.2 implies that there exists such τ 0 > 0 that for all τ > τ 0 the equilibrium s 0 is locally exponentially stable under the condition (5.7).
The statement about the instability of the equilibrium follows from the fact that for |κ| > κ c there exist such l ∈ {1, . . . , q}, j ∈ {2, . . . , n} and ω ∈ R that |κµ j | > |g l (ω)|, and, hence the real part of the asymptotic continuous spectrum is positive γ l,j (ω, σ) > 0. Therefore, for sufficiently large time-delays and for |κ| > κ c , there will be eigenvalues from the pseudo-continuous spectrum with positive real parts.
With the result of persistence of equilibrium solution proven, we move to the case of periodic synchronous solution.

5.2.
Synchronization of periodic orbits. We have previously stated our result on the synchronization of periodic orbits for the SL oscillators (Theorem 3.2), so, we could avoid the discussion about the characteristic equation H. Now we present a generalization of the announced result and we prove it in this section.
Remark 5.1. The co-dimension 1 condition ∂ 2 H(0, 0) = 0 may lead to isolated points, where the synchronization of the synchronous periodic orbit cannot be achieved (see example in Sec. 6.2).
The proof of Theorem 5.1 rely on the following lemma, adapted from [33] for the case of the stability of a periodic synchronized solution s(t) with respect to desynchronized perturbations (transverse to the synchronization manifold).
Lemma 5.2. The synchronous periodic orbit s(t) of system (1.1) with period T is exponentially stable with respect to perturbations transverse to the synchronization manifold for all sufficiently large τ if all of the following conditions hold: S-I (No strong instability) All elements of the instantaneous spectrum Γ I have negative real parts (this implies in particular that the strongly unstable spectrum Γ SU is empty), S-II (Simple trivial multiplier) The trivial Floquet exponent λ = 0 is simple.
And, it is exponentially unstable with respect to perturbations transverse to the synchronization manifold for sufficiently large τ if one of the following conditions hold: U-I (Strong instability) the strongly unstable spectrum Γ SU is non-empty, or U-II (Weak instability) a non-empty subset of the asymptotic continuous spectrum Γ A (−κµ j ) has positive real part for nontrivial eigenvalue µ j of the Laplacian matrix L.
Proof. This theorem follows almost directly from Theorem 6 of [33], and we highlight the main features for our setup. First of all, the existence of the periodic solution s(t) does not depend on time delay τ . As a result, τ can be considered as a continuous parameter here, instead of the parameter N which is an integer part of τ /T used in [33]. Hence, the asymptotic statements with respect to N are equivalently formulated for τ in the present theorem.
Secondly, the stability of the periodic solution is determined by the variational equation (4.3) that splits into the part along the synchronization manifold with µ 1 = 0 and the remaining part for the transverse perturbations with nonzero µ j . Hence, the transverse stability is determined by the variational equation (4.6) with σ ∈ {−κµ 2 , · · · − κµ n }. The conditions for the stability S-I -S-III guarantee that the spectrum of the corresponding variational equations is stable for large enough τ .
The linear stability of a synchronous periodic orbit is governed by the linearized equation (4.4) with periodic J(t). In Ref. [33] the authors proved results about the stability of this equation for large τ . In particular, the existence of a holomorphic function H : Ω 1 × Ω 2 ⊆ C × C → C is shown, such that λ is a Floquet exponent of the linear DDE (4.4) if and only if H λ, σe −λτ = 0 (see Eq. (4.8)).
The function H in this case is analogous to the characteristic equation (5.1) for the case of equilibrium. The main difference of the periodic case is that the function H is not given explicitly. The asymptotic continuous spectrum is determined from the equation H iω, σe −γ e −iφ = 0, (5.8) (compare (5.8) and (4.10)). Here we should emphasize that, in contrast to Ref. [33], we explicitly write the parameter σ in the argument of the function H since the dependence on σ is of interest for this study. Moreover, this parameter can be eliminated from Eq.
In particular, this expression shows that for α P > 0, the periodic solution will be destabilized for positive σ, and stabilized for negative σ. For α P < 0, the stabilization occurs for positive σ and destabilization for negative. It is important to notice, that the stabilization (or destabilization) occurs for all transverse modes simultaneously, since the eigenvalues µ 2 , . . . , µ n of the Laplacian matrix are positive and σ can admit values −κµ j , which have the same sign for all µ j , j = 2, . . . , n.

6.
Example: Coupled Stuart-Landau systems. Let us consider an example of the ring of coupled Stuart-Landau (SL) oscillators shown in Fig. 6.1. The dynamics of a node j is given by (3.3). Throughout this section, we use the identity as the coupling function. We rewrite Eq.  with identity coupling function and parameters α = 1, β = π, τ = 20, and κ = 0.7. The points approaching the curve on the left side of the figure belongs to the pseudo-continuous spectrum and the isolated points on the right belongs to the strongly unstable spectrum. Solid lines show the re-scaled asymptotic continuous spectrum Γ A . The gray strip represents a break on the figure, which is necessary due to the different scales of the two parts of the spectrum.
6.1. Persistence of equilibrium. Firstly, we consider the case α < 0 (equilibrium solution) so that the eigenvalues of the uncoupled system have negative real parts. The Jacobian matrix at zero z = 0 has the eigenvalues α±iβ. Therefore, the strongly unstable spectrum is empty, and the whole spectrum of the zero equilibrium for the network system (4.3) consists only of the asymptotic-continuous spectrum for long delays.
The laplacian matrix of the network in Fig. 6.1 has spectral radius ρ L = 2. In order to compute the value r 0 from the condition (3.1) of Theorem 3.1 we compute the functions g l (ω), l = 1, 2 using Eq. (5.4):

SYNCHRONIZATION IN NETWORKS WITH STRONGLY DELAYED COUPLINGS 13
Hence, we find r 0 = min l inf ω |g l (ω)| = |α|. Therefore, the synchronization manifold for the considered network is locally exponentially stable if and only if 0 < κ < κ c = |α| 2 (6.2) for sufficiently long time delay τ (and unstable if κ > κ c ). Fig. 6.3 illustrates the convergence of the trajectories to the equilibrium (left panel) for the case when the condition (6.2) is fulfilled and the absence of convergence in the opposite case. The detailed parameter values are given in the figure's caption. Moreover, in order to compute the synchronization error we compare each solution to a fixed one taking the norm of maximal difference, that is, we measure the synchronization error by max j z 1 (t) − z j (t) . 6.2. Periodic orbit synchronization . Let us consider the case α > 0 when the synchronous periodic solution √ αe iβt exists. In this section we prove Theorem 3.2. Consider the transformation z(t) = r(t)e iβt , then Eq. (6.1) reads, in terms of r(t),ṙ Note that the solution z(t) = √ αe iβt is transformed into a family of equilibria. The variational equation of (6.3) (considering real and imaginary parts) around the equilibrium solution r = √ α, equivalent to (4.3), iṡ ξ j (t) = −2α 0 0 0 ξ j (t) − κµ j cos(βτ ) − sin(βτ ) sin(βτ ) cos(βτ ) ξ j (t − τ ). (6.4) where µ j ≥ 0 are the eigenvalues of L. For short, we write Eq. (6.4) asξ j (t) = J 0 ξ j (t) − κµ j T ξ j (t − τ ).
The instantaneous spectrum can be computed using Definition 4.1. It consists of the eigenvalues of the non-delayed part of (6.4) and reads Γ I = {−2α, 0}. The eigenvalue 0 in the instantaneous spectrum is associated with the trivial Floquet multiplier producing a singularity in the asymptotic spectrum. Then, as predicted by Theorem 5.1, the synchronization cannot be achieved for an arbitrary large delay. However, Theorem 5.1 guarantee the existence of a synchronization interval (0, κ c ) (or (κ c , 0) with κ c < 0) with the decreasing length with increasing τ .
We can also compute the spectrum of (6.4) by using Eq. (5.1). So, we get the transcendental equation with σ = −κµ j . Eq. (6.7) can be numerically solved. We compare the solutions of (6.7) to the analytical approximation. The asymptotic spectrum (given in terms of Eq. (6.6)) and the spectrum points (solutions of Eq. (6.7)) are given in Fig. 6.4 with parameters given in the figure's caption.
From the data of Fig. 6.4 we notice that for α = 1, β = π and τ = 20, stale synchronization is attained for 0 < κ 0.04. Now, in order to have a complete view of the stability domain, we produce a synchronization map in the parameter space σ × τ . The Fig. 6.5 shows such a color map presenting the values of max m (λ m ) < 0, where λ m , m ∈ Z, is a solution of Eq. 6.7, in the σ × τ parameter space. Note that σ = −κρ L and, as predicted by Theorem 5.1, the stable synchronization either occurs for 0 < κ < κ c or for −κ c < κ c < 0 with κ c shrinking as τ grows.
The parameter values −1 ≤ σ ≤ 1 of the color map in Fig. 6.5 can be related to more complicated connected network for coupling parameter values in the range −1/ρ L ≤ κ < 1/ρ L . The interchange between stability domains occurs at the values  with η > 0 to be the largest real part of the eigenvalues. Hence, we define the characteristic time t tr as t tr = 1/η.
The characteristic time measures how fast the slowest solution ξ(t) approaches 0. One can write η = −γ max /τ where γ max is the maximum of the real part of the asymptotic continuous spectrum given in Eq. (5.5). It can be computed as Hence, the characteristic time is Clearly, t tr → ∞ as κ → κ c . Fig. 6.6 shows a test for characteristic time given by Eq. (6.9) using two coupled Stuart-Landau oscillators with parameters α = −1, β = π, and τ = 20. In this example, Eq. (6.9) reads t tr (κ) = −20 ln −1 (2κ).  The red curve is t tr (κ) = 20 ln −1 (2κ). The blue dots were obtained by fixing κ and computing η, which stands for the angular coefficient of Eq. (6.8) in log scale in which ||ξ(t)|| = ||x 1 (t) − x 2 (t)||, and then taking t tr = 1/η. The parameters used were α = −1, β = π and τ = 20. The history functions were taken as constant and non-zero. 7. Synchronization loss versus network structure. Here, we explore the relationship between growing networks (with strongly delayed connection) and its synchronization window in the case of equilibrium. The main results are the corollaries 3.1 and 3.2.
We will concentrate on the behavior of large networks such as Barabási-Albert scale-free network, Erdös-Réniy random network, and some regular graphs, for which the expressions for the Laplacian spectral radius ρ L are known.
First, let us take a look at how regular graphs respond to the synchronization (using the network model (1.1) with long delay) in the limit of large network size. Table 7.1 lists spectral radius of the Laplacian matrix of the main regular graphs: complete, ring, star, and path.
Using Theorem 3.1, and Table 7.1 we observe that: For strong delay and a large network size n → ∞, the synchronization manifold tends to be always unstable in networks of the types Complete or Star provided that the coupling strength κ does not scale with n.
For strong delay and a large network size n → ∞, the synchronization manifold tends to be stable for a certain interval of the coupling strength κ ∈ (0, r/4) in simple networks of the types Ring or Path. Now, let us consider some complex networks. We relate the synchronization condition to the graph structure for two important examples of complex networks.
Homogeneous networks, characterized by a small disparity in the node degrees. A canonical example is the Erdös-Réniy (ER) random network: Starting with n nodes the graph is constructed by connecting nodes randomly. Each edge is included in the graph with probability 0 < p < 1 independent from every other edge. If p = p 0 ln n/n with p 0 > 1 then the ER network is almost surely connected [5]. We will consider that p is chosen so that the ER network is connected.
Heterogeneous networks, where some nodes (called hubs) are highly connected whereas most of the nodes have only a few connections. A canonical example of such networks is the Barabási-Albert (BA) scale-free network. One way to construct it is to start with two nodes and a single edge that connect them. Then at each step, a new node is created and it is connected with one of the preexisting nodes with a probability proportional to its degree. This process is called preferential attachment.
More details about the construction, structure, and dynamics of such networks can be found in [27,4]. Illustrations of Erdös-Réniy (ER) random networks (homogeneous) and Barabási-Albert (BA) Scale-Free networks (heterogeneous) can be seen in Figure 7  The response of the BA and ER networks to synchronization under the considered delayed model already stated in Sec. 3 and encoded in corollaries 3.1 and 3.2. The both results says that ER and BA networks tend to not have stable synchronization in the limit of large n. But, ER networks, losses stable synchronization at a slow rate compared to BA networks.
Here we prove the both cited corollaries. The proof of Corollary 3.1 is based on the following well-known results.
Lemma 7.1 (Ref. [3]). Let g max denote the largest degree of a undirected network G of size n. Then the spectral radius ρ L of the Laplacian matrix L of G has the following estimates: n n − 1 g max ≤ ρ L ≤ 2g max . (7.1) Lemma 7.2 (Ref. [26]). Consider a BA undirected network of size n and g max its largest degree. With probability 1 we have the limit is almost surely positive and finite.
Proof of Corollary 3.1. The Lemma 7.2 implies that the maximum degree of a BA network grows as √ n when we let n → ∞. Therefore, using Eq. (7.1) and (7.2) we see that ρ L ≥ µ √ n with probability 1 in the limit of large n which implies that the synchronization window, that is, the interval (0, r/ρ L ) (see Theorem 3.1) shrinks and vanishes with large n.
Remark 7.1. Corollary 3.1 shows that the synchronization manifold tends to be unstable in large and strongly delayed BA (heterogeneous) networks if the coupling strength is not scaled with n or τ .
For the proof of Corollary 3.2, we use Lemma 7.1 and the following result. Lemma 7.3 (Ref. [30]). Consider an Erdös-Réniy network with n nodes, connection probability 0 < p < 1 and q = 1 − p. Then, the probability that the maximum degree g max being at most np + b √ npq, for some constant b > 0, tends to 1 as n → ∞.
Proof of Corollary 3.2. For large n and using the probability p = p 0 ln n/n with p 0 > 1 one get g max ∼ p 0 ln n + O √ ln n .
So, using Lemma 7.1, we end up with ρ L ≥ n n − 1 p 0 ln n + O √ ln n which means that ρ L → ∞ when n → ∞ with a rate of order ln n. The spectral radius ρ L is the lowest possible when p 0 → 1 + maintaining the network connected, or, ρ L is the lowest possible when the network crosses the connectivity threshold becoming connected.
Although Corollary 3.2 says that the ER random network does not allow the synchronization manifold to be stable in the limit of n → ∞, the rate in which it happens is slow, making it possible to synchronize very large ER networks with long time delays.
The properties observed for the stability of the synchronization manifold for BA and ER networks with strong delay, that is, large BA networks doesn't support strong delay interaction and relatively large ER networks supports strong delay interaction, are similar to the persistence of the synchronization when the nondelayed coupling functions are nonidentical [25]. 7.1. Synchronization of BA and ER networks when coupling strength is scaled with n. As we have seen from the previous subsection, not only simple but also complex networks tend to have an always unstable synchronization manifold with strong delay τ and large network size n.
In many network dynamical models, it is common to normalize the coupling parameter, in our case κ, in order to preserve certain behaviors and properties of the network, such as synchronization, mean field oscillation, community clustering, etc [1,31,14]. When the network size is dynamical, for instance, the size of the network is growing with time, then this normalization becomes even more important.
With this notion in mind, we can see that in the regime of large network size, the coupling parameter κ should have some natural scaling depending on the network structure. If those scaling are taken into account in the network model then the stability of the synchronization manifold can always be preserved for n → ∞.
For example, if we know beforehand that we have a large simple network of the type Complete (or Star) then the natural scaling of the coupling parameter would be κ → κ n .
We then state here further consequences of the corollaries 3.1 and 3.2 when considering the scaling of the coupling parameter. Then for any large enough size n and long enough delay, there is always a nonempty interval I = (0, κ max (n)) ⊂ R such that any synchronous steady state is locally exponentially stable for all κ s ∈ I and unstable if κ s ∈ R \ I (with probability 1). Moreover, the length of this interval converges to a nonzero constant with n → ∞ with probability 1: lim n→∞ κ max (n) = κ ∞ , 0 < κ ∞ < ∞.
Proof. Any equilibrium solution is locally exponentially stable if 0 < κ < r/ρ L and unstable for κ > r/ρ L for some r > 0. For BA scale-free networks we know that ρ L ∼ √ n (see Lemma 7.2). Then, if we re-scale the coupling parameter as κ → κ s / √ n the new synchronization condition is 0 < κ s / √ n < r/ρ L ≈ r/ √ n, which leads to 0 < κ s < r. Then for large enough network size n and long enough time delay τ , there is an interval I = (0, κ max (n)) ⊂ R such that any equilibrium solution is locally exponentially stable if κ s ∈ I and unstable if κ s ∈ R \ I (with probability 1). Moreover, the length of the interval I converges to a nonzero constant with n → ∞ with probability 1: lim n→∞ κ max (n) = κ ∞ , 0 < κ ∞ < ∞.
The proof of the Corollary 7.2 follows the same steps of the Corollary 7.1 using the Lemma 7.3.