Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 Contents lists available at ScienceDirect Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns Research paper Stability of small-amplitude periodic solutions near Hopf bifurcations in time-delayed fully-connected PLL networks Diego P. Ferruzzo Correa a , ∗, Átila M. Bueno b , José R. Castilho Piqueira c a Universidade Federal do ABC, São Bernardo, São Paulo, Brazil b Universidade Estadual Paulista Júlio de Mesquita Filho, Campus Experimental de Sorocaba, Sorocaba, São Paulo, Brazil c Escola Politécnica da Universidade de São Paulo, São Paulo, Brazil a r t i c l e i n f o Article history: Received 13 August 2015 Revised 8 September 2016 Accepted 2 October 2016 Available online 3 October 2016 Keywords: Time-delay differential equations Networks Synchronization Hopf bifurcation Stability a b s t r a c t In this paper we investigate stability conditions for small-amplitude periodic solutions emerging near symmetry-preserving Hopf bifurcations in a time-delayed fully-connected N-node PLL network. The study of this type of systems which includes the time delay be- tween connections has attracted much attention among researchers mainly because the delayed coupling between nodes emerges almost naturally in mathematical modeling in many areas of science such as neurobiology, population dynamics, physiology and engi- neering. In a previous work it has been shown that symmetry breaking and symmetry preserving Hopf bifurcations can emerge in the parameter space. We analyze the stabil- ity along branches of periodic solutions near fully-synchronized Hopf bifurcations in the fixed-point space, based on the reduction of the infinite-dimensional space onto a two- dimensional center manifold in normal form. Numerical results are also presented in order to confirm our analytical results. © 2016 Elsevier B.V. All rights reserved. 1. Introduction Networks of oscillators have been studied for decades because their models can represent dynamics in a very wide range of fields as astronomy, biology, neurology, economics, population dynamics, and the stock market (see [1–7] and references therein). Much of the research has focused into understand the influence that changes in the parameter space have over the dynamics. We are interested in studying the influence of the lag between nodes in the stability of the network syn- chronization. There is a considerable body of literature on time-delayed networks, for example: in [3] and [8] are explored conditions for the global exponential stability; in [9] is addressed an statistical approach including analysis of noise influ- ence in linearly-coupled oscillators, in [10] are studied Hopf and Bogdanov-Takens bifurcations in a small neural network; in [11] different kind of solutions including amplitude death, spatiotemporal, phase-locked, standing-waves and synchro- nized oscillations are studied considering the time-delay along with a distance-dependent coupling, in [12] is addressed a comparative study of three different models for a fully-connected N-node network and it is also shown the existence of multiple eigenvalues forced by symmetry, in [13] is presented a stability criterion for the synchronization in a two-node network. ∗ Corresponding author. E-mail addresses: diego.ferruzzo@ufabc.edu.br (D.P. Ferruzzo Correa), atila@sorocaba.unesp.br (Á.M. Bueno), piqueira@lac.usp.br (J.R. Castilho Piqueira). http://dx.doi.org/10.1016/j.cnsns.2016.10.001 1007-5704/© 2016 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cnsns.2016.10.001 http://www.ScienceDirect.com http://www.elsevier.com/locate/cnsns http://crossmark.crossref.org/dialog/?doi=10.1016/j.cnsns.2016.10.001&domain=pdf mailto:diego.ferruzzo@ufabc.edu.br mailto:atila@sorocaba.unesp.br mailto:piqueira@lac.usp.br http://dx.doi.org/10.1016/j.cnsns.2016.10.001 D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 67 In this contribution we delve in the study of the full-fase model presented in [12] , specifically we focus on the stabil- ity along branches of small-amplitude periodic solutions near symmetry-preserving Hopf bifurcations for a N-node Phase- Locked-Loop (PLL) network using the time-delay between nodes as bifurcation parameter (the existence of these solutions was proved in [12] ). We analyze the stability of the branches in the parameter space ( μ, τ ) for the non degenerative case ( K > 1) applying the center manifold theorem extended to functional differential equations and the normal form in the center space to reduce the order of the infinite-dimensional system. It is important to note that when the lag between the nodes in a network is considered the ordinary differential equation (ode) system which describe the network dynamics becomes a delayed differential equation (dde) system. Solutions for a dde system lie in the function space and its characteristic equation has infinitely many roots, this particular kind of functional differential equations appear in many engineering problems [14–16] . The main approach used for our analysis is the decomposition of the infinite-dimensional space into to two subspaces, a two-dimensional center space spanned by the eigenvectors corresponding to the simple imaginary eigenvalues λ = ±i ω, ω > 0, and an infinite-dimensional space “orthogonal” to the first one (the orthogonality condition will be defined below). We will follow closely the theory and procedures presented in [14,17–20] . 2. The full phase model The general model for a N -node, fully-connected, second-order oscillator network in terms of the i -th node output phase φi ( t ), is: φ̈i (t) + μ ˙ φi (t) − μ − Kμ N − 1 N ∑ j � = i j=1 f (φi , φ j ) = 0 , i = 1 , . . . , N, (1) where the sumatory term represents the coupling function and the remaining terms represent the local inner second-order dynamics in each node. The coupling function for a PLL oscillator is given by: f (φi , φ j ) = sin (φ j (t − τ ) − φi (t)) + sin (φ j (t − τ ) + φi (t)) , (2) we consider that all signals coming from the other N − 1 nodes are affected by a lag τ ; thus, f : R × R → R ; μ, K, τ ∈ R + are parameters, and N ∈ N − { 1 } . The equilibria φ±, in Eq. (1) , are: φ+ (n ) = 1 2 ( arcsin ( − 1 K ) + 2 nπ ) φ−(n ) = 1 2 ( π − arcsin ( − 1 K ) + 2 nπ ) , (3) n ∈ Z , K ≥ 1 . For our analysis, we consider three main assumptions: (a) The critical eigenvalue λ of the linearization of (1) at equilibria crosses the imaginary axis with non vanishing velocity, i.e. Re( λ′ ( φ±)) � = 0. (b) The purely imaginary eigenvalue λ = i ω is simple. (c) The linearization of (1) at equilibria has no eigenvalues of the form i k ω, k ∈ Z − { 1 , −1 } . The Taylor expansion of (1) at equilibria is: δφ̈i + μδ ˙ φi − Kμ N − 1 N ∑ j � = i j=1 ∞ ∑ r=1 { 1 r! ( δφi ∂ ∂φ′ i + δφ jτ ∂ ∂φ′ jτ )r f (φi , φ jτ ) } φ′ i = φ± φ′ jτ = φ± = 0 . (4) where φ jτ := φ j (t − τ ) . Truncate the series up to the third-order term: φ̈i + μ ˙ φi = Kμ N − 1 N ∑ j � = i j=1 { ( φ jτ − φi ) + ( φ jτ + φi ) cos 2 φ± −1 2 ( φ jτ + φi )2 sin 2 φ± − 1 6 [ ( φ jτ − φi )3 + ( φ jτ + φi )3 cos 2 φ± ] } , (5) i = 1 , . . . , N, here for the sake of notation we changed δφ → φ . i i 68 D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 The vector field form ˙ x = G (x t , x ;η) , G : R 2 N × R 2 N × R 3 → R 2 N , can be obtained by choosing x (i ) 1 = φi and x (i ) 2 = ˙ φi , then the restriction G | i or ˙ x (i ) = G (i ) (x t , x (i ) ;η) gives: ˙ x (i ) 1 = x (i ) 2 ˙ x (i ) 2 = −μx (i ) 2 + Kμ N − 1 N ∑ j � = i j=1 { (−1 + cos 2 φ±) x (i ) 1 + (1 + cos 2 φ±) x ( j) 1 τ − 1 2 ( x ( j) 1 τ + x (i ) 1 )2 sin 2 φ± −1 6 [ ( x ( j) 1 τ − x (i ) 1 )3 + ( x ( j) 1 τ + x (i ) 1 )3 cos 2 φ± ] } , i = 1 , . . . , N. (6) Following [17,21] , we can represent the dynamics in (6) by the abstract differential equation: d dt x t (φ) = A (η) x t (φ) + F(x t (φ) , η) . (7) We define X := C([ −τ, 0] , R 2 N ) the Banach space of continuous functions from [ −τ, 0] into R 2 N equipped with the usual norm ‖ ϑ‖ = sup −τ≤θ≤0 | ϑ(θ ) | , ϑ ∈ C([ −τ, 0]) , x t , in (7) , lies in X and satisfies (T (t) φ)(θ ) = (x t (φ))(θ ) = x (t + θ ) , T ( t ) is a semigroup of family of operators, θ ∈ [ −τ, 0] , and η is a vector of parameters. The linear operator A ( η) ∈ Mat(2 N ) is: (A (η) ϑ) = { ∂ϑ ∂θ (θ ) , −τ < θ ≤ 0 A 0 (η) ϑ(0) + A τ (η) ϑ(−τ ) , θ = 0 , (8) where A 0 (η) := ∂G ∂x ∣∣ φ± , A τ (η) := ∂G ∂x τ ∣∣ φ± and, ( F(x ) ) (θ ) = { ∂x ∂θ (θ ) , −τ ≤ θ < 0 F (x (0) , x (−τ ) , η) , θ = 0 , (9) F = ( f (1) , . . . , f (N) ) T , f (i ) = ( f (i ) 1 , f (i ) 2 ) , f (i ) 1 = 0 , and f (i ) 2 = Kμ N − 1 N ∑ j � = i j=1 { − 1 2 ( x ( j) 1 τ + x (i ) 1 )2 sin 2 φ± − 1 6 [ ( x ( j) 1 τ − x (i ) 1 )3 + ( x ( j) 1 τ + x (i ) 1 )3 cos 2 φ± ] } . In order to build the decomposition of the infinite-dimensional space we need two tools: the adjoint operator associated to the linear part of the linearization, and an inner product via a bilinear form. Associated to the linear part of (7) the formal adjoint equation is dy dt (t, η) = A T 0 (η) y (t, η) + A T τ (η) y (t + τ, η) . (10) The strongly continuous semigroup (T ∗(t) ψ)(θ ) = (y t (ψ))(θ ) = y (t + θ ) , defines the infinitesimal generator: (A ∗(η) ψ) = { ∂ψ ∂θ (θ ) , 0 < θ ≤ τ A 0 (η) T ψ(0) + A τ (η) T ψ(τ ) , θ = 0 , (11) such that d dt T ∗(t) ψ = A ∗T ∗(t) ψ, ψ ∈ X ∗ := C([0 , τ ] , R 2 N ) . The natural inner product has the form [22] : 〈 x, y 〉 = x̄ T (0) y (0) + ∫ 0 −τ x̄ T (s + τ ) A τ (η) y (s ) ds, x ∈ X and y ∈ X ∗; thus, we have [17] : 1. λ is an eigenvalue of A ( η) if and only if λ̄ is and eigenvalue of A ∗( η). 2. If ϕ 1 , . . . , ϕ d is a basis for the eigenspace of A ( η) and ψ 1 , . . . , ψ d is a basis for the eigenspace of A ∗( η), construct the matrices � = (ϕ 1 , . . . ϕ d ) and � = (ψ 1 , . . . , ψ d ) . Define the bilinear form: 〈 �, �〉 = I. (12) D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 69 3. The fixed point space S N Due to the S N -symmetry of (1) the space where solutions φi lie can be decomposed into the fixed point subspace where symmetry-preserving solutions emerge and a subspace with symmetry-breaking solutions, this was shown in [12] . We an- alyze stability of the small-amplitude periodic solutions near Hopf bifurcations in the fixed point space, these bifurcations satisfy assumptions (a)–(c) for K > 1. In this subspace Eq. (6) has the form: ˙ x 1 = x 2 ˙ x 2 = −μx 2 + Kμ(−1 + cos 2 φ±) x 1 + Kμ { (1 + cos 2 φ±) x 1 τ − 1 2 (x 1 τ + x 1 ) 2 sin 2 φ± −1 6 [ (x 1 τ − x 1 ) 3 + (x 1 τ + x 1 ) 3 cos 2 φ±]} , (13) then matrices A 0 ( η) and A τ ( η) in (11) become: A 0 (η) = ( 0 1 Kμ(−1 + cos 2 φ±) −μ ) , (14) A τ (η) = ( 0 0 Kμ(1 + cos 2 φ±) 0 ) , (15) and F in (9) takes the form F = ( f 1 f 2 ) T with f 1 = 0 and f 2 : f 2 (x t , η) = Kμ { − 1 2 (x 1 τ + x 1 ) 2 sin 2 φ± − 1 6 [ (x 1 τ − x 1 ) 3 + (x 1 τ + x 1 ) 3 cos 2 φ±]} . (16) We need the complex eigenfunctions As (ϑ) = i ωs (ϑ) , A ∗n (θ ) = i ωn (θ ) , associated to the critical eigenvalues λ = i ω, and λ̄ = −i ω with s (ϑ) = s 1 (ϑ) + i s 2 (ϑ) and n (θ ) = n 1 (θ ) + i n 2 (θ ) . These eigenfunctions can be computed solving the boundary value problem d dϑ s 1 , 2 = ∓ωs 2 , 1 (ϑ) , and d dθ n 1 , 2 = ±ωn 2 , 1 (θ ) , which after substituting the operator A ( η) becomes: A 0 (η) s 1 (0) + A τ (η) s 1 (−τ ) = −ωs 2 (0) A 0 (η) s 2 (0) + A τ (η) s 2 (−τ ) = ωs 1 (0) (17) and A T 0 (η) n 1 (0) + A T τ (η) n 1 (−τ ) = ωn 2 (0) A T 0 (η) n 2 (0) + A T τ (η) n 2 (−τ ) = −ωn 1 (0) , (18) with general solutions: s 1 (ϑ) = cos (ωϑ) c 1 − sin (ωϑ) c 2 s 2 (ϑ) = sin (ωϑ) c 1 + cos (ωϑ) c 2 n 1 (θ ) = cos (ωθ ) d 1 − sin (ωθ ) d 2 n 2 (θ ) = sin (ωθ ) d 1 + cos (ωθ ) d 2 . (19) The coefficients c 1 = [ c 11 c 12 ] T , c 2 = [ c 21 c 22 ] T , d 1 = [ d 11 d 12 ] T , d 2 = [ d 21 d 22 ] T can be obtained by considering the boundary conditions ( A 0 (η) + cos (ω τ ) A τ (η) ω I + sin (ω τ ) A τ (η) )T ( c 1 c 2 ) = 0 ( A T 0 (η) + cos (ω τ ) A T τ (η) −ω I − sin (ω τ ) A T τ (η) )T ( d 1 d 2 ) = 0 , (20) the “orthonormality” condition 〈 s, n 〉 = I, and setting c 11 = 1 and c 21 = 0 , see [14,23] for more details. It is also possible to decompose the solution x t ( ϑ) to Eq. (7) into x t (ϑ) = y 1 (t) s 1 (ϑ) + y 2 (t) s 2 (ϑ) + w t (ϑ) , where y 1 and y 2 lie in the center subspace, such that y 1 , 2 (t) = 〈 n 1 , 2 (0) , x t (0) 〉 , and w t in the infinite-dimensional component subspace, thus, we have ˙ y 1 = ωy 2 + n T 1 (0) F ˙ y 2 = −ωy 1 + n T 2 (0) F (21) ˙ w = A (η) w t + F(x t , η) − n T 1 (0) F s 1 − n T 2 (0) F s 2 , (22) where F = { 0 , ϑ ∈ [ −τ, 0) (23) F (y 1 (t) s 1 (0) + y 2 (t) s 2 (0) + w (t)(0)) , ϑ = 0 . 70 D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 3.1. The center manifold Following [14,18,24] we know that w t ( ϑ) can be approximated by the second-order expansion: w (y 1 , y 2 )(ϑ) = 1 2 (h 1 (ϑ) y 2 1 + 2 h 2 (ϑ) y 1 y 2 + h 3 (ϑ) y 2 2 ) , (24) thus, by differentiating and substituting Eq. (22) keeping up to second order terms, we obtain: ˙ w = −ω h 2 y 2 1 + ω (h 1 − h 3 ) y 1 y 2 + ωh 2 y 2 2 + O (y 3 ) , (25) and from Eq. (22) , d w dt = A (η) w t + F( w + y 1 s 1 + y 2 s 2 ) − (d 12 s 1 + d 22 s 2 ) f 2 . (26) From the definition of A ( η), equivalent to (11) , we see that A (η) w = { 1 2 ( ̇ h 1 y 2 1 + 2 ̇ h 2 y 1 y 2 + ˙ h 3 y 2 2 ) , ϑ ∈ [ −τ, 0) A 0 ( η) w ( 0) + A τ (η) w (−τ ) , ϑ = 0 , (27) then, from Eqs. (24) to (27) we can obtain the unknown coefficients h 1 , h 2 , and h 3 solving: ˙ h 1 = 2(−ωh 2 + f 20 2 (d 12 s 1 (ϑ) + d 22 s 2 (ϑ))) , ˙ h 2 = ω(h 1 − h 3 ) + f 11 2 (d 12 s 1 (ϑ) + d 22 s 2 (ϑ)) , ˙ h 3 = 2(ωh 2 + f 02 2 (d 12 s 1 (ϑ) + d 22 s 2 (ϑ))) , (28) and, A 0 (η) h 1 (0) + A τ (η) h 1 (−τ ) = 2(−ωh 2 (0) + f 20 2 (d 12 s 1 (0) + d 22 s 2 (0))) , A 0 (η) h 2 (0) + A τ h 2 (−τ ) = ω(h 1 (0) − h 3 (0)) + f 11 2 (d 12 s 1 (0) + d 22 s 2 (0))) , A 0 (η) h 3 (0) + A τ (η) h 3 (−τ ) = 2(ωh 2 (0) + f 02 2 (d 12 s 1 (0) + d 22 s 2 (0))) , (29) where f 20 = 1 2 ∂ 2 f ∂y 2 1 ∣∣∣ 0 , f 11 = ∂ 2 f ∂ y 1 ∂ y 2 ∣∣∣ 0 , and f 02 = 1 2 ∂ 2 f ∂y 2 2 ∣∣∣ 0 . Eq. (28) is written as the inhomogeneous differential equation: dh dϑ = Ch + p cos (ωϑ) + q sin (ωϑ) (30) where h := ⎛ ⎝ h 1 h 2 h 3 ⎞ ⎠ , C := ω ⎛ ⎝ 0 −2 I 0 I 0 −I 0 2 I 0 ⎞ ⎠ 6 ×6 p := ⎛ ⎜ ⎝ f 20 2 p 0 f 11 2 p 0 f 02 2 p 0 ⎞ ⎟ ⎠ , q := ⎛ ⎜ ⎝ f 20 2 q 0 f 11 2 q 0 f 02 2 q 0 ⎞ ⎟ ⎠ , p 0 := ( d 12 c 22 d 22 ) , q 0 := ( d 22 −c 22 d 12 ) , with general solution: h (ϑ) = e Cϑ K + M cos (ωϑ) + N sin (ωϑ) . (31) After substituting the general solution into (30) we solve for M and N , and then from the boundary value problem we solving for K , ( C −ωI ωI C )( M N ) = − ( p q ) (32) P h (0) + Qh (−τ ) = p − r, (33) D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 71 where P := ( A 0 0 0 0 A 0 0 0 0 A 0 ) − C, Q := ( A τ 0 0 0 A τ 0 0 0 A τ ) , (34) and r := ( 0 f 20 2 0 f 11 2 0 f 02 2 )T . The expressions for w 1 (0) and w 1 (−τ ) , necessary in (23) , are: w 1 (0) = 1 2 ( (M 1 + K 1 ) y 2 1 + 2(M 3 + K 3 ) y 1 y 2 + (M 5 + K 5 ) y 2 2 ) , w 1 (−τ ) = 1 2 ( (e −Cτ K| 1 + M 1 cos (ωτ ) − N 1 sin (ωτ )) y 2 1 +2(e −Cτ K| 3 + M 3 cos (ωτ ) − N 3 sin (ωτ )) y 1 y 2 +(e −Cτ K| 5 + M 5 cos (ωτ ) − N 5 sin (ωτ )) y 2 2 ) , (35) note that we only need w 1 ( ϑ) since the nonlinear function in (16) only depends on x 1 ; then by substituting (35) into (21) we obtain: ˙ y 1 = ωy 2 + g 1 (y 1 , y 2 ;η) ˙ y 2 = −ωy 1 + g 2 (y 1 , y 2 ;η) , (36) or ˙ y 1 = ωy 2 + a 20 y 2 1 + a 11 y 1 y 2 + a 02 y 2 2 + a 30 y 3 1 + a 21 y 1 y 2 + a 12 y 1 y 2 2 + a 03 y 3 2 , ˙ y 2 = −ωy 1 + b 20 y 2 1 + b 11 y 1 y 2 + b 02 y 2 2 + b 30 y 3 1 + b 21 y 1 y 2 + b 12 y 1 y 2 2 + b 03 y 3 2 . (37) In [20] is computed the coefficient a which determines stability of the normal form (37) , a = 1 16 [ g 03 2 + g 21 2 + g 12 1 + g 30 1 ] + 1 16 ω [ g 11 2 ( g 02 2 + g 20 2 ) − g 11 1 ( g 02 1 + g 20 1 ) − g 02 2 g 02 1 + g 20 2 g 20 1 ] , (38) where g i j r = ∂ i + j ∂ i y 1 ∂ j y 2 g r (0 , 0) . Periodic orbits near Hopf bifurcation at the critical eigenvalue λ = i ω will be stable if a < 0 and unstable if a > 0. 4. Numerical results We reproduce some of the computations for the Hopf bifurcations curves in the fixed point space for the case K > 1 presented in [12] in order to compute the coefficient a for small-amplitude periodic solutions near these bifurcation curves using results obtained in the previous section. Fig. 1 shows the symmetry-preserving Hopf bifurcation curves in the parameter space ( μ, τ ) for K = 1 . 05 for both cases: bifurcations with Re( λ′ ) > 0 (black curves), and with Re( λ′ ) < 0 (red Fig. 1. Symmetry-preserving bifurcation curves in Fix( S N ) for K = 1 . 05 . In black Hopf bifurcations with Re( λ′ ) > 0, and in red Hopf bifurcations with Re( λ′ ) < 0. 72 D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 Fig. 2. Coefficient a ( Eq. (38) ) computed for the Hopf bifurcation curves in Fix( S N ) for K = 1 . 05 (see Fig. 1 ). Figure (a): curves in plane ( μ, a ), figure (b): curves for a in the parameter space ( μ, τ ). Fig. 3. (a) Branch of periodic solutions emerging from point A = (μ, τ ) = (0 . 15 , 7 . 46197) . (b) Periodic solution profile at μ = 0 . 15 , τ = 7 . 5315 , T = 12 . 0364 seg, (point psol ). (c) Floquet multipliers for the periodic solution psol . curves). These curves correspond to the equilibrium φ−(n ) in Eq. (3) , and each lobe correspond to a different value of n ∈ N . We choose three testing point for numerical simulation A = (μ, τ ) = (0 . 15 , 7 . 46) , B = (0 . 3 , 11) , and C = (0 . 421 , 7 . 10) . Fig. 2 a shows the coefficient a computed using Eq. (38) in the parameter space ( μ, τ ) for K = 1 . 05 related to the Hopf bifurcations curves shown in Fig. 1 ; the black curves correspond to stability of periodic orbits near Hopf bifurcations with Re( λ′ ) > 0, as we can see also in Fig. 2 b all these curves are under the plane a = 0 , therefore, all these periodic solutions are stable; the red curves correspond to stability of periodic orbits near Hopf bifurcations with Re( λ′ ) < 0, these periodic orbits are unstable for μ < μc ( n ), and stable for μ > μc . Fig. 2 a also shows points A, B , and C ; small amplitude periodic orbits are stable at points A and C whilst at point B they are unstable. In order to confirm our results we computed branches of periodic solutions near the Hopf bifurcations points A, B , and C using DDE-BIFTOOL [25,26] along with the Floquet multipliers for a specific periodic solution chosen in the branch ( Figs. 3–5 ). D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 73 Fig. 4. (a) Branch of periodic solutions emerging from point B = (μ, τ ) = (0 . 3 , 11 . 001518) . (b) Periodic solution profile at μ = 0 . 3 , τ = 11 . 3744 , T = 12 . 8506 seg, (point psol ). (c) Floquet multipliers for the periodic solution psol . Fig. 5. (a) Branch of periodic solutions emerging from point C = (μ, τ ) = (0 . 421 , 7 . 101329) . (b) Periodic solution profile at μ = 0 . 421 , τ = 7 . 00 , T = 8 . 8704 seg, (point psol ). (c) Floquet multipliers for the periodic solution psol . Fig. 3 -(a) shows a branch of periodic solutions with small amplitude emerging from the Hopf bifurcation point A = (μ, τ ) = (0 . 15 , 7 . 46) ; Fig. 3 -(b) shows the periodic solution profile psol at τ = 7 . 5315 ; Fig. 3 -(c) shows the Floquet mul- tipliers related to psol . It is clear that this periodic solution is stable since there is no Floquet multiplier outside the unity circle. 74 D.P. Ferruzzo Correa et al. / Commun Nonlinear Sci Numer Simulat 45 (2017) 66–74 For the point B = (μ, τ ) = (0 . 3 , 11) , the branch of periodic solutions is shown in Fig. 4 -(a), Fig. 4 -(b) shows the profile of the periodic solution psol chosen at τ = 11 . 3744 , this solution is unstable because there is a Floquet multiplier outside the unity circle, see Fig. 4 -(c). Finally, the branch of periodic solutions near the Hopf bifurcation point C = (μ, τ ) = (0 . 421 , 7 . 10) is shown in Fig. 5 -(a); the periodic solution chosen in the branch is at τ = 7 . 00 , its profile is shown in Fig. 5 -(b); all the Floquet multipliers shown in Fig. 5 -(c) are within the unity circle, therefore the solution is stable. 5. Conclusions The reduction of the infinite-dimensional space onto the center manifold in normal form was applied to the fixed point space for the full phase model in order to analyze the stability of small-amplitude periodic orbits near simple Hopf bifurca- tions. For the case Re( λ′ ) > 0 we found that stable ( a < 0) periodic orbit can emerge, and for the case Re( λ′ ) < 0 unstable ( a > 0) periodic orbits can emerge for μ < μc ( n ), and stable ( a < 0) periodic orbits for μ > μc ( n ). The numerics show that our analytical results are correct. Although we computed the coefficient a for a specific value of K the procedure shown is valid for all the parameter space where simple Hopf bifurcations emerge. Finally, it is important to spotlight some points for further research: First, analyze the nature of the solutions at the special point μ = μc (n ) at which the coefficient a changes sign. Second, analyze the case K = 1 , degenerate Hopf bifurca- tions codimension 2 may appear (pure imaginary and zero eigenvalue), where fold-Hopf and Bautin bifurcations (generalized Hopf bifurcations) could emerge; and third, the stability of the symmetry-breaking degenerate Hopf bifurcations which have multiplicity N − 1 . Acknowledgment We would like to thank the Escola Politécnica da Universidade de São Paulo and FAPESP for their support. References [1] Gy ̋ori I , Ladas G . Oscillation theory of delay differential equations: with applications. New York: Oxford Mathematical Monographs. The Clarendon Press Oxford University Press; 1991 . [2] Martins A , Monteiro LHA . Frequency transitions in synchronized neural networks. Commun Nonlinear Sci Numer Simul Jul 2013;18(7):1786–91 . [3] Li W , Yang H , Wen L , Wang K . Global exponential stability for coupled retarded systems on networks: a graph-theoretic approach. Commun Nonlinear Sci Numer Simul Jun 2014;19(6):1651–60 . [4] Earl M , Strogatz S . Synchronization in oscillator networks with delayed coupling: a stability criterion. Phys Rev E Jan 2003;67(3) . [5] Ponzi A , Aizawa Y . Self-organized criticality and partial synchronization in an evolving network. Chaos Solitons Fractals Jan 20 0 0;11(7):1077–86 . [6] Zhang C , Li W , Wang K . Graph-theoretic approach to stability of multi-group models with dispersal. Discrete Continuous Dyn Syst - Ser B Jan 2015a;20(1):259–80 . [7] Zhang C , Li W , Wang K . Graph theory-based approach for stability analysis of stochastic coupled systems with lvy noise on networks. IEEE Trans Neural Netw Learn Syst Aug 2015b;26(8):1698–709 . [8] Yuan Y , Zhao X-Q . Global stability for non-monotone delay equations (with application to a model of blood cell production). J Differ Equ Feb 2012;252(3):2189–209 . [9] Vuellings A , Schoell E , Lindner B . Spectra of delay-coupled heterogeneous noisy nonlinear oscillators. Eur Phys J B Feb 3 2014;87(2) . [10] Giannakopoulos F , Zapp A . Bifurcations in a planar system of differential delay equations modeling neural activity. Physica D Nov 2001;159(3–4):215–32 . [11] Guo Y , Niu B . Amplitude death and spatiotemporal bifurcations in nonlocally delay-coupled oscillators. Nonlinearity May 2015;28(6):1841 . [12] Correa DPF , Wulff C , Piqueira JRC . Symmetric bifurcation analysis of synchronous states of time-delayed coupled phase-locked loop oscillators. Com- mun Nonlinear Sci Numer Simul Aug 2015;22(13):793–820 . [13] Correa DPF , Piqueira JRC . Synchronous states in time-delay coupled periodic oscillators: a stability criterion. Commun Nonlinear Sci Numer Simul Aug 2013;18(8):2142–52 . [14] Kalmár-Nagy T , Stépán G , Moon FC . Subcritical Hopf bifurcation in the delay equation model for machine tool vibrations. Nonlinear Dyn Oct 2001;26(2):121–42 . [15] Stone E , Campbell SA . Stability and bifurcation analysis of a nonlinear DDE model for drilling. J Nonlinear Sci Jan-Feb 2004;14(1):27–57 . [16] Gilsinn DE . Estimating critical Hopf bifurcation parameters for a second-order delay differential equation with application to machine tool chatter. Nonlinear Dyn Oct 2002;30(2):103–54 . [17] Gilsinn DE . Bifurcations, center manifolds, and periodic solutions. In: Delay differential equations. New York: Springer; 2009. p. 155–202 . [18] Zhao S , Kalmár-Nagy T . Center manifold analysis of the delayed Lienard equation. In: Delay differential equations. US: Springer; 2009. p. 1–17 . [19] Campbell S . Calculating centre manifolds for delay differential equations using maple. In: Delay differential equations. US: Springer; 2009. p. 1–24 . [20] Guckenheimer J , Holmes P . Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. New York: Springer Verlag; 1983 . [21] Hale JK . Functional differential equations. New York: Springer-Verlag; 1971 . [22] Hale JK , Lunel SMV . Introduction to functional differential equations. London: Springer-Verlag; 1993 . [23] Hale JK . Theory of functional differential equations (applied mathematical sciences). Springer; 1977 . [24] Hassard BD , Kazarinoff ND , Wan YH . Theory and applications of Hopf bifurcation. London mathematical society lecture note series, volume 41. Cam- bridge: Cambridge University Press; 1981 . [25] Engelborghs K , Luzyanina T , Samaey G . DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differential equations. Numerical analysis and applied mathematics section. K.U.Leuven, Leuven, Belgium: Department of Computer Science; Oct 2001 . [26] Engelborghs K , Luzyanina T , Roose D . Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL. ACM Trans Math Softw Mar 2002;28(1):1–21 . http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0001 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0001 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0001 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0002 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0002 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0002 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0003 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0003 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0003 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0003 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0003 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0004 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0004 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0004 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0005 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0005 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0005 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0006 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0006 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0006 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0006 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0007 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0007 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0007 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0007 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0008 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0008 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0008 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0009 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0009 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0009 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0009 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0010 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0010 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0010 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0011 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0011 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0011 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0012 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0012 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0012 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0012 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0013 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0013 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0013 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0014 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0014 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0014 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0014 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0015 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0015 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0015 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0016 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0016 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0017 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0017 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0018 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0018 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0018 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0019 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0019 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0020 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0020 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0020 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0021 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0021 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0022 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0022 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0022 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0023 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0023 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0024 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0024 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0024 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0024 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0025 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0025 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0025 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0025 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0026 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0026 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0026 http://refhub.elsevier.com/S1007-5704(16)30329-X/sbref0026 Stability of small-amplitude periodic solutions near Hopf bifurcations in time-delayed fully-connected PLL networks 1 Introduction 2 The full phase model 3 The fixed point space SN 3.1 The center manifold 4 Numerical results 5 Conclusions Acknowledgment References