The mean field analysis of the Kuramoto model on graphs II.Asymptotic stability of the incoherent state, center manifold reduction, and bifurcations

The mean field analysis of the Kuramoto model on graphs II.
Asymptotic stability of the incoherent state, center manifold reduction, and bifurcations

Hayato Chiba  and Georgi S. Medvedev Institute of Mathematics for Industry, Kyushu University / JST PRESTO, Fukuoka, 819-0395, Japan, chiba@imi.kyushu-u.ac.jpDepartment of Mathematics, Drexel University, 3141 Chestnut Street, Philadelphia, PA 19104, medvedev@drexel.edu
Abstract

In our previous work [3], we initiated a mathematical investigation of the onset of synchronization in the Kuramoto model (KM) of coupled phase oscillators on convergent graph sequences. There, we derived and rigorously justified the mean field limit for the KM on graphs. Using linear stability analysis, we identified the critical values of the coupling strength, at which the incoherent state looses stability, thus, determining the onset of synchronization in this model.

In the present paper, we study the corresponding bifurcation. Specifically, we show that similar to the original KM with all-to-all coupling, the onset of synchronization in the KM on graphs is realized via a pitchfork bifurcation. The formula for the stable branch of the bifurcating equilibria involves the principal eigenvalue and the corresponding eigenfunctions of the kernel operator defined by the limit of the graph sequence used in the model. This establishes an explicit link between the network structure and the onset of synchronization in the KM on graphs. The results of this work are illustrated with the bifurcation analysis of the KM on Erdős-Rényi, small-world, as well as certain weighted graphs on a circle.

1 Introduction

In 1970s, a prominent Japanese physicist Yoshiki Kuramoto described a remarkable effect in collective dynamics of large systems of coupled oscillators [11]. He studied all-to-all coupled phase oscillators with randomly distributed intrinsic frequencies, the model which now bears his name. When the strength of coupling is small, the phases are distributed approximately uniformly around a unit circle, forming an incoherent state. Kuramoto identified a critical value of the coupling strength, at which the the incoherent state looses stability giving rise to a stable partially synchronized state. To describe the bifurcation corresponding to the onset of synchronization, Kuramoto introduced the order parameter - a scalar function, which measures the degree of coherence in the system. He further showed that the order parameter undergoes a pitchfork bifurcation. Thus, the qualitative changes in the statistical behavior of a large system of coupled phase oscillators near the onset to synchronization can be described in the language of the bifurcation theory.

Kuramoto’s discovery created a new area of research in nonlinear science [20, 21, 19]. The rigorous mathematical treatment of the pitchfork bifurcation in the KM was outlined in [4] and was presented with all details in [1]. The analysis in these papers is based on the generalized spectral theory for linear operators [2] and applies to the KM with intrinsic frequencies drawn from a distribution with analytic or rational probability density function. Under more general assumptions on the density, the onset of synchronization in the KM was analyzed in [6, 8, 5]. These papers use analytical methods for partial differential equations and build upon a recent breakthrough in the analysis of Landau damping [18].

In our previous work [3], we initiated a mathematical study of the onset of synchronization in the KM on graphs. Following [14, 15], we considered the KM on convergent families of deterministic and random graphs, including Erdős-Rényi and small-world graphs among many other graphs that come up in applications. For this model, we derived and rigorously justified the mean field limit. The latter is a partial differential equation approximating the dynamics of the coupled oscillator model in the continuum limit as the number of oscillators grows to infinity. In [3], we performed a linear stability analysis of the incoherent state and identified the boundaries of stability. Importantly, we related the stability region of the incoherent state to the structural properties of the network through the spectral properties of the kernel operator defined by the limit of the underlying graph sequence [12]. In the present paper, we continue the study initiated in [3]. Here, we analyze the bifurcations at the critical values of the coupling strength, where the incoherent state looses stability. In the remainder of this section, we review the model and the main outcomes of [3] and explain the main results of this work.

We begin with a brief explanation of the graph model that was used in [3] and will be used in this paper. In [3], we adapted the construction of W-random graphs [13] to define a convergent sequence of graphs. The flexible framework of W-random graphs allows us to deal with a broad class of networks that are of interest in applications. Let be a symmetric measurable function on the unit square with values in and let form a triangular array of points subject to the following condition

(1.1)

is a weighted graph on nodes labeled by the integers from , whose edge set is

Each edge is supplied with the weight In the theory of graph limits, is called a graphon [12]. It defines the asymptotic properties of for large .

Consider the initial value problem (IVP) for the KM on

(1.2)
(1.3)

The intrinsic frequencies are independent identically distributed random variables. The distribution of has density . For the spectral analysis in Section 3, we need to impose the following assumptions on : a) is an even unimodal function, and b) is real analytic function with finite moments of all orders: . For instance, the density of the Gaussian distribution satisfies these conditions. The KM on weighted graphs (1.2), (1.3) can be used to approximate the KM on a variety of random graphs (cf.  § 4.2 [3]).

Along with the discrete model (1.2) we consider the IVP for the following partial differential equation

(1.4)
(1.5)

where

(1.6)

Here, is the conditional density of the random vector given , and parametrized by , and is a circle. In particular,

(1.7)

It is shown in [3, Theorem 2.2] that

(1.8)

interpreted as a probability measure on Borel sets converges in the bounded Lipschitz distance [7] uniformly on bounded time intervals to the absolutely continuous measure

(1.9)

provided and are sufficiently close in the same distance. The latter can be achieved with the appropriate initial condition (1.3) and sufficiently large (see [3, Corollary 2.3]). Therefore, the IVP (1.4),(1.5) approximates the IVP (1.2),(1.3) on finite time intervals for sufficiently large .

An inspection of (1.4) shows that , the density of the uniform distribution on , is a steady state solution of (1.4). It corresponds to the incoherent (mixing) state of the KM. Numerics suggests that the incoherent state is stable for small The loss of stability of the incoherent state is interpreted as the onset of synchronization in the KM. This is the main focus of [3] and of the present paper. In [3], we identified the boundaries of the region of stability of the incoherent state in (1.4). Specifically, we showed that there exist such that is linearly stable for , and is unstable otherwise.

The critical values and depend on the network topology through the eigenvalues of the compact symmetric operator

(1.10)

The eigenvalues of are real with the only accumulation point at . Denote the largest positive and smallest negative eigenvalues of by and respectively. If all eigenvalues are nonnegative (nonpositive), we set (). The main stability result of [3] yields explicit expressions for the transition points

(1.11)

Thus, the region of linear stability of depends explicitly on the spectral properties of the limiting graphon . Recall that represents the graph limit of . Thus, (1.11) links network topology to synchronization in (1.2). For the classical KM (all-to-all coupling), and , which recovers the known Kuramoto’s formula.

In the present paper, we study the onset of synchronization in (1.2) in more detail. After some preliminaries and preparatory work in Sections 2 and 3, we revisit linear stability of the incoherent solution. This time, we show that despite the lack of eigenvalues with negative real part and the presence of the continuous spectrum on the imaginary axis, the incoherent state is an asymptotically stable solution of the linearized problem (cf. Theorem 4.1). This is a manifestation of the Landau damping in the KM.

In Section 5, we study the bifurcation at with a one-dimensional center manifold. To this end, we recall the order parameter

(1.12)

which was introduced in [3] as a measure of coherence in the KM on graphs. This is a continuous analog of the local order parameter

for the discrete model (1.2). The order parameter generalizes the original order parameter used by Kuramoto for the all-to-all coupled model. Note that (1.12) depends on and contains information about the structure of the network through . As will be clear below, the order parameter plays an important role in the analysis of the mean field equation. In particular, it can be used to locate nontrivial steady state solutions. To this end note that the velocity field (1.6) can be conveniently rewritten in terms of the order parameter

(1.13)

In particular, for a given steady state of the order parameter written in the polar form

(1.14)

the velocity field takes the following form

Setting , we find the corresponding steady state solution of (1.4)

(1.15)

where stands for the Dirac delta function. The stationary solution (1.15) has the following interpretation: the first line describes phase-locked oscillators, while the second line yields the distribution of the drifting oscillators. Thus, solutions of this form may combine phase-locked oscillators and those moving irregularly. Such solutions are called partially phase-locked or partially synchronized. The phase of a phase-locked oscillator at with a natural frequency , is given by

(1.16)

provided In Sections 5 and 6, we will identify branches of stable equilibria bifurcating from in terms of the corresponding values of the order parameter. Then Equation (1.16) will be used to describe the corresponding stable phase-locked solutions.

In Section 5, assuming that is a simple eigenvalue of we show that the coupled system (1.2) undergoes a supercritical pitchfork bifurcation at . Specifically, we derive an ordinary differential equation for the order parameter and show that the trivial solution of this equation looses stability at and gives rise to a stable branch of (nontrivial) equilibria, corresponding to partially synchronized state (cf. (5.1)). In Section 6, we consider the onset of synchronization in networks with certain symmetries (cf. (6.1)). This leads to the bifurcation with a two-dimensional center manifold. The bifurcation analysis in Sections 5 and 6 is illustrated with the analysis of the KM on Erdős-Rényi, small-world graphs, and to a class of weighted graphs on a circle.

2 Preliminaries

In the remainder of this paper, we will assume that . The case of negative is reduced to that above by switching to and . Furthermore, without loss of generality we assume that .

2.1 Fourier transform

We rewrite (1.4) in terms of the complex Fourier coefficients

(2.1)

Applying the Fourier transform to (1.4) and using integration by parts, we obtain

(2.2)

where

(2.3)

From (2.2), we further obtain

(2.4)

By (1.7), . Further, because is real. Thus, in (2.4) we can restrict to .

Let

(2.5)

Combining these observations, we rewrite (2.4):

(2.6)
(2.7)

Note that the trivial solution is a steady state solution of (2.6), (2.7). It corresponds to the uniform distribution , a constant steady state solution of (1.4). Linearizing around we arrive at

(2.8)
(2.9)

where is a linear operator on

(2.10)

2.2 The eigenvalue problem

The multiplication operator defined by

(2.11)

is a closed operator. The continuous spectrum of fills the imaginary axis

(2.12)

Since is compact (as a Hilbert-Schmidt operator), is closed and .

Next we turn to the eigenvalue problem

(2.13)

where and are operators on (cf. (2.10) and (2.5)).

We will locate the eigenvalues of through the eigenvalues of (cf. (1.10)). Since is a compact symmetric operator on , it has a countable set of real eigenvalues with the only accumulation point at zero. All nonzero eigenvalues have finite multiplicity.

Suppose is an eigenvalue of and is the corresponding eigenfunction. Then a simple calculation yields (cf. [3])

(2.14)

where

(2.15)
(2.16)

Equation (2.14) yields the equation for eigenvalues of

(2.17)

where is a nonzero eigenvalue of .

Using (2.17), we establish a one-to-one correspondence between the eigenvalues of and those of . Specifically, for every positive eigenvalue of , , there is a branch of eigenvalues of ,

(2.18)

such that

(2.19)

Recall that stands for the largest positive eigenvalue of . Then for there are no eigenvalues with positive real part. Furthermore, for small and there is a unique positive eigenvalue of , which vanishes as (see [3] for more details).

3 The generalized spectral theory

The major obstacle in studying stability and bifurcations of the incoherent state is the continuous spectrum of the linearized problem on the imaginary axis (cf. (2.12)). To deal with this difficulty, we develop the generalized spectral theory following the treatment of the classical KM in [1]. Below, we outline the key steps in the analysis of the generalized eigenvalue problem referring the interested reader to [1] for missing proofs and further details.

3.1 The rigged Hilbert space

Here, we define the rigged Hilbert space which will be used in the spectral analysis below. Let be the set of holomorphic functions on the region such that the norm

(3.1)

is finite. With this norm, is a Banach space [1]. Let be their inductive limit with respect to and

(3.2)

With the inductive limit topology, is a complete Montel space 111A locally convex topological vector space is called Montel if every bounded set in it is relatively compact.. The properties of are described in detail in [1].

Next, we define as a topological tensor product of and . is a complete Montel space. It is a dense subspace of The topology of is stronger than that of . For every we have for each . In addition, is holomorphic in on the upper half plane, where it can grow at most exponentially.

Let be the dual space of . It is the space of continuous antilinear functionals on . Let denote the pairing between and i.e., for and , stands for the corresponding antilinear functional.

The dual space is equipped with the weak dual topology (cf. [1])222 converges to if tends to for every .. is a dense subspace of , and the topology of is stronger than that of . Hence,

form a rigged Hilbert space (a.k.a. Gelfand triple) [10].

Note that if then

for , where Thus, and the rigged Hilbert space defined above satisfy all assumptions of the generalized spectral theory in [1].

3.2 The generalized eigenvalue problem

In this subsection we calculate the resolvent of and spectral projections. With the rigged Hilbert space defined above, we will view the resolvent as an operator from to .

Below, we will need to construct analytic continuation for certain functions involving integrals of Cauchy type. For this, we are going to use an implication of the Sokhotski formulas, which we formulate as a separate statement for convenience.

Lemma 3.1.

(Sokhotski, cf. [9]) Let be a complex valued function on . Suppose has at most a finite number of integrable discontinuities. Then

(3.3)

is an analytic function in the right and left open half-planes of . Furthermore, for , the following formulas determine the limits of as :

(3.4)

where stands for the principal value in the sense of Cauchy and denotes the Hilbert transform of .

Corollary 3.2.

Suppose is holomorphic on the real axis and admits the analytic continuation to the upper half-plane. Then

(3.5)

is an entire function.

3.3 The generalized resolvent

Our next goal is to compute the resolvent of

(3.6)

To this end, we first compute for and extend it analytically to the left half-plane as an operator from to .

In the right half-plane can be rewritten as follows

(3.7)

where

(3.8)

and stands for the identity operator. Note that ceases to exist as the multiplication operator on as (recall that the imaginary axis is the continuous spectrum of ). However, it can be extended to the left half-plane as as an operator defined as follows

(3.9)

By Corollary 3.2, is an entire function in for all . This suggests an appropriate generalization of defined by

(3.10)

where is the dual operator of .

For each , is an -valued meromorphic function. For , it coincides with the restriction of to . Thus, is a meromorphic continuation of from the right half-plane to the left-half plane as an -valued operator. Note that since has the continuous spectrum on the imaginary axis, can not be continued to the left-half plane as an operator on .

We define the generalized eigenvalues of as the singularities of the generalized resolvent

Definition 3.3.

is called a generalized eigenvalue of if there is a nonzero such that

(3.11)

In this case, is called a generalized eigenfunction.

It turns out that the generalized eigenvalues and the corresponding eigenfunctions of are, in fact, the eigenvalues and eigenfunctions of the dual of .

Theorem 3.4.

(cf. [2]) Let be a generalized eigenvalue of and is the corresponding eigenfunction. Then .

Remark 3.5.

Using (3.9) and (3.11), one can see that the generalized eigenvalues of are the roots of the following equation

(3.12)

where is a nonzero eigenvalue of and

(3.13)

The right hand side of (3.13) is an entire function (cf. Corollary 3.2). For , (3.12) is reduced to the equation for the eigenvalues of (cf. (2.17)). In this case, the corresponding generalized eigenfunction is included in i.e., is an eigenvalue of . On the other hand, for , the generalized eigenfunction is not in but is an element of the dual space .

3.4 The generalized Riesz projection

Let be a positive eigenvalue of and be the corresponding eigenfunction. The largest positive eigenvalue of and the corresponding eigenfunction are denoted by and respectively. For every there is a real positive eigenvalue of . The corresponding eigenfunction is given by

(3.14)

As approaches the critical value from above, the eigenvalue converges to along the real axis and at it hits the continuous spectrum on the imaginary axis. The corresponding eigenfunction approaches the critical vector

(3.15)

where stands for the limit in with respect to the weak dual topology, i.e., the action of on is given by

(3.16)

Let be a generalized eigenvalue of . Then the generalized Riesz projection is defined by

(3.17)

where is a simple closed curve around oriented counterclockwise that does not encircle or intersect the rest of the spectrum. Below, we shall refer to such curves as contours. The image of gives the generalized eigenspace of  [1].

Theorem 3.6.

Suppose the algebraic and geometric multiplicities of coincide. Then the generalized Riesz projection of the generalized eigenvalue of for , has the following form

(3.18)

where is a positive constant, was defined in (3.8), and stands for the Riesz projection onto the eigenspace of corresponding to the eigenvalue . The operator on is defined by

(3.19)

The proof of Theorem 3.6 relies on three technical lemmas. Below we state and prove these lemmas first and then prove the theorem.

Lemma 3.7.

Let then

(3.20)
Proof.

By definition of (3.6), for any , we have

and, thus,

(3.21)

Using Fubini theorem, from (2.5) we have

(3.22)

On the other hand, integrating both sides of (3.21) against , we obtain

(3.23)

and

(3.24)

Plugging (3.24) into (3.22), we have

(3.25)

The combination of (3.21) and (3.25) yields (3.20). ∎

Lemma 3.8.

Let be an eigenvalue of corresponding to the positive eigenvalue of , , and , and suppose that the geometric and algebraic multiplicities of coincide.

Then

(3.26)

where is the Riesz projection onto the eigenspace of corresponding to .

Proof.

As before, let denote a contour around . From (3.20), we have

(3.27)

We change variable in the integral on the right–hand side to . By deforming the contour if necessary, we can always achieve for so that this change of variable is well defined. Under this transformation, is mapped to , a contour around . Thus, we have

(3.28)

Since the algebraic and geometric multiplicities of are equal, the singularity of at is a simple pole, and the other factor in the integrand of the above is regular at . Therefore, the right–hand side of (3.28) simplifies to

(3.29)

By multiplying both sides of (3.29) by , we have

Finally, since is the projection on the eigensubspace of

Thus,

(3.30)

Lemma 3.9.
(3.31)
Proof.

Using integration by parts times, we obtain

The application of Lemma 3.1 to the integral on the right-hand side yields (3.31). ∎

Below will need the following implications of Lemma 3.9.

Corollary 3.10.
(3.32)
(3.33)
Proof.

Differentiating and using (3.31), for off the imaginary axis we have

(3.34)

The integral on the right–hand side is of Cauchy type and Lemma 3.1 applies. By (3.4),

(3.35)

Since is even, and is odd. Because is also nonnegative and unimodal Thus,

(3.36)

The combination of (3.4) and (3.36) yields (3.32).

Likewise, (3.33) follows from Lemma 3.9 for and Lemma 3.1. ∎

Proof of Theorem 3.6.

Theorem 3.6 follows from (3.30) and (3.32). ∎

4 Asymptotic stability of the incoherent state

We now return to the problem of stability of the incoherent state. Recall that in the Fourier space the incoherent state corresponds to the trivial solution . The linearization about shows that it is a neutrally stable equilibrium of (2.8), (2.9) for . There are no eigenvalues of for these values of and the continuous spectrum fills out the imaginary axis. Nonetheless, we show that the incoherent state is asymptotically stable with respect to the weak dual topology.

Theorem 4.1.

For the trivial solution of (2.8), (2.9) is an asymptotically stable equilibrium for initial data from with respect to the weak dual topology on .

Remark 4.2.

The stability with respect to the weak dual topology is weaker than that with respect to the topology of the Hilbert space . Still it is a natural topology for the problem at hand. In particular, Theorem 4.1 implies that the order parameter evaluated on the trajectories of the linearized problem tends to as .

Figure 1: Deformation of the integral path for the Laplace inversion formula.
Proof.

Integrating (2.9) subject to , we have

By the Riemann-Lebesgue lemma,

(4.1)

We now turn to (2.8). Consider the semigroup generated by the operator on given by the Laplace inversion formula

(4.2)

where is arbitrary. Thus, the (continuous) spectrum of lies to the left of the integration path along (see Fig. 1a).

For arbitrary , we have

(4.3)

For , is an analytic function in the right half–plane, which can be extended to the entire complex plane as a meromorphic function . Thus,

(4.4)

Let be fixed. Next we claim that one can choose such that there are no generalized eigenvalues of on or inside the contour

for every . To construct with the desired property, we first fix . Then we recall that generalized eigenvalues of satisfy (3.12). From (3.13), under our assumptions on , there exists such that there are no roots of (3.12) in the region

because (3.12) can be reduced to in for . On the other hand, is holomorphic. Thus, the set of roots of (3.12) (i.e., the set of generalized eigenvalues) does not have accumulation points in

Thus, we can choose such that there are no generalized eigenvalues in . This completes the construction of .

By the Cauchy Integral theorem,

(4.5)

for any , and