Symmetry broken states in an ensemble of globally coupled pendulums

TL;DR

This paper analyzes the rotational dynamics of globally coupled pendulums, revealing various symmetry-broken states, their stability, and bifurcations, extending the Kuramoto model to include inertia and nonlinearity.

Contribution

It introduces a generalized model of coupled pendulums with inertia, analyzing symmetry-breaking, stability, and bifurcations both analytically and numerically.

Findings

Identification of instability boundaries for in-phase rotations

Discovery of multiple out-of-phase regimes and their stability

Detailed bifurcation analysis revealing regular and chaotic states

Abstract

We consider the rotational dynamics in an ensemble of globally coupled identical pendulums. This model is essentially a generalization of the standard Kuramoto model, which takes into account the inertia and the intrinsic nonlinearity of the community elements. There exists the wide variety of in-phase and out-of-phase regimes. Many of these states appear due to broken symmetry. In the case of small dissipation our theoretical analysis allows one to find the boundaries of the instability domain of in-phase rotational mode for ensembles with arbitrary number of pendulums, describe all arising out-of-phase rotation modes and study in detail their stability. For the system of three elements parameter sets corresponding to the unstable in-phase rotations we find a number of out-of-phase regimes and investigate their stability and bifurcations both analytically and numerically. As a result, we obtain a sufficiently detailed picture of the symmetry breaking and existence of various regular and chaotic states.

Equations

$${\ddot{\phi }}_n+\lambda {\dot{\phi }}_n+\sin {\phi }_n=\gamma +\frac{K}{N}\sum _{\tilde{n}=1}^N\sin \left({\phi }_{\tilde{n}}-{\phi }_n\right),$$

$$\ddot{\varphi }+\lambda \dot{\varphi }+\sin \varphi =\gamma ,$$

$$\varphi (\tau )=\tau +\frac{{\lambda }^2}{{\gamma }^2}\sin \tau +o\left({\lambda }^4\right),\tau =\left(\frac{\gamma }{\lambda }-\frac{{\lambda }^3}{2{\gamma }^3}+o\left({\lambda }^7\right)\right)t.$$

$$\begin{array}{c}\delta {\ddot{\phi }}_n+\lambda \delta {\dot{\phi }}_n+\cos \varphi (t){\delta \phi }_n=\frac{K}{N}\sum _{\tilde{n}=1}^N({\delta \phi }_{\tilde{n}}-{\delta \phi }_n).\end{array}$$

$$\begin{array}{c}\eta =\frac{1}{N}\sum _{\tilde{n}=1}^N\delta {\phi }_{\tilde{n}},\\ {\xi }_n=\delta {\phi }_{n+1}-\delta {\phi }_n,n=1,\dots ,N-1,\end{array}$$

$$\begin{array}{c}\ddot{\eta }+\lambda \dot{\eta }+\cos \varphi (t)\eta =0,\end{array}$$

$$\begin{array}{c}{\ddot{\xi }}_n+\lambda {\dot{\xi }}_n+(K+\cos \varphi (t)){\xi }_n=0,\\ n=1,\dots ,N-1.\end{array}$$

$$\begin{array}{c}{\ddot{\xi }}_n+\lambda \omega {\dot{\xi }}_n+{\omega }^2\left[K+\cos (\varphi (\tau ))\right]{\xi }_n=0,\\ n=1,\dots ,N-1.\end{array}\text{vspace}-.1in$$

$$\begin{array}{c}{K}_{1,2}=\frac{1}{4}\left[\frac{{\gamma }^2}{{\lambda }^2}\mp 2\sqrt{1-{\gamma }^2}+\frac{1}{2}\frac{{\lambda }^2}{{\gamma }^2}\right]+O({\lambda }^4).\text{vspace}-.1in\end{array}$$

$${\ddot{\varphi }}_m+\lambda {\dot{\varphi }}_m+\sin {\varphi }_m=\gamma +\frac{K}{N}\sum _{\tilde{m}=1}^M{N}_{\tilde{m}}\sin \left({\varphi }_{\tilde{m}}-{\varphi }_m\right).\text{vspace}-.1in$$

$$\begin{array}{c}\delta {\ddot{\phi }}_{mk}+\lambda \delta {\dot{\phi }}_{mk}+\cos {\varphi }_m\delta {\phi }_{mk}\\ =\frac{K}{N}\sum _{\tilde{m}=1}^M\sum _{\tilde{k}=1}^{N_{\tilde{m}}}\cos \left({\varphi }_{\tilde{m}}-{\varphi }_m\right)\left(\delta {\phi }_{\tilde{m}\tilde{k}}-\delta {\phi }_{mk}\right).\text{vspace}-.1in\end{array}$$

$$\begin{array}{c}{\eta }_m=\frac{1}{N_m}\sum _{k=1}^{N_m}\delta {\phi }_{mk},m=1,\dots ,M,\\ {\xi }_{mk}=\delta {\phi }_{m,k+1}-\delta {\phi }_{mk},k=1,\dots ,N_m-1.\text{vspace}-.1in\end{array}$$

$$\begin{array}{c}{\ddot{\eta }}_m+\lambda {\dot{\eta }}_m+\cos {\varphi }_m{\eta }_m\\ =K\sum _{\tilde{m}=1}^M\frac{N_{\tilde{m}}}{N}\cos \left({\varphi }_{\tilde{m}}-{\varphi }_m\right)\left({\eta }_{\tilde{m}}-{\eta }_m\right),\text{vspace}-.1in\end{array}$$

$$\begin{gathered} {\ddot{\xi }}_{mk}+\lambda {\dot{\xi }}_{mk}+\left(\cos {\varphi }_m+K\sum _{\tilde{m}=1}^M\frac{N_{\tilde{m}}}{N}\cos \left({\varphi }_{\tilde{m}}-{\varphi }_m\right)\right){\xi }_{mk} \\ =0. \end{gathered}$$

$$\begin{array}{c}{\ddot{\tilde{\eta }}}_m+\left(2{\Lambda }^{in}+\lambda \right){\dot{\tilde{\eta }}}_m+\left({\left({\Lambda }^{in}\right)}^2+\lambda {\Lambda }^{in}+\cos {\varphi }_m\right){\tilde{\eta }}_m\\ =K\sum _{\tilde{m}=1}^M\frac{N_{\tilde{m}}}{N}\cos \left({\varphi }_{\tilde{m}}-{\varphi }_m\right)\left({\tilde{\eta }}_{\tilde{m}}-{\tilde{\eta }}_m\right).\text{vspace}-.1in\end{array}$$

$${\Lambda }^{in}=\sum _{j=1}^{+\infty }{\lambda }^j{\Lambda }^{(j)}.\text{vspace}-.08in$$

$$\begin{array}{c}{\ddot{\tilde{\xi }}}_m+\left(2{\Lambda }^{ex}+\lambda \right){\dot{\tilde{\xi }}}_m+\left({\left({\Lambda }^{ex}\right)}^2+\lambda {\Lambda }^{ex}+\cos {\varphi }_m\right){\tilde{\xi }}_m\\ =-K\sum _{\tilde{m}=1}^M\frac{N_{\tilde{m}}}{N}\cos \left({\varphi }_{\tilde{m}}-{\varphi }_m\right){\tilde{\xi }}_m.\text{vspace}-.1in\end{array}$$

$$\begin{array}{c}{K}_c^{(2:1)}=\frac{{\gamma }^2}{4{\lambda }^2}+\frac{\sqrt{1-{\gamma }^2}}{3}+O(\lambda ).\end{array}$$

$$\begin{array}{c}{K}_c^{(1:1:1)}=\frac{{\gamma }^2}{4{\lambda }^2}+\frac{\sqrt{1-{\gamma }^2}}{3}+O(\lambda ).\end{array}$$

$$\begin{array}{c}{K}_{c1}^{(2:1:1)}=\frac{{\gamma }^2}{4{\lambda }^2}+\frac{\sqrt{1-{\gamma }^2}}{4}+O(\lambda ),\\ K_{c2}^{(2:1:1)}=\frac{{\gamma }^2}{4{\lambda }^2}+\frac{\sqrt{1-{\gamma }^2}}{2}+O(\lambda ).\end{array}$$

$${\varphi }_1^{\pm }(\tau )$$

$$+\frac{{\lambda }^2}{{\gamma }^2}\left[\frac{\gamma }{24{\gamma }_1}\left(36{\gamma }_1\mp {\left(A_2^{\pm }\right)}^2\right)-\gamma \cos 2\tau \pm {\gamma }_1\sin 2\tau \right]$$

$$\displaystyle+A_{2}^{\pm}\frac{\lambda^{3}}{\gamma^{3}}\Biggl{[}\pm\frac{16}{15}\gamma_{1}+B_{2}^{\pm}\cos\tau-\frac{\gamma}{2}\sin\tau$$

$$\displaystyle-\frac{1}{72}\left(\left(A_{2}^{\pm}\right)^{2}\mp 12\gamma_{1}\right)\cos 3\tau+\frac{\gamma}{6}\sin 3\tau\Biggr{]}+O(\lambda^{4}),$$

$${\varphi }_2^{\pm }(\tau )$$

$$+\frac{{\lambda }^2}{{\gamma }^2}\left[\frac{\gamma }{24{\gamma }_1}\left(36{\gamma }_1\mp {\left(A_2^{\pm }\right)}^2\right)-\gamma \cos 2\tau \pm {\gamma }_1\sin 2\tau \right]$$

$$\displaystyle+A_{2}^{\pm}\frac{\lambda^{3}}{\gamma^{3}}\Biggl{[}\pm\frac{16}{15}\gamma_{1}-2B_{2}^{\pm}\cos\tau+\gamma\sin\tau$$

$$\displaystyle+\frac{1}{36}\left(\left(A_{2}^{\pm}\right)^{2}\mp 12\gamma_{1}\right)\cos 3\tau-\frac{\gamma}{3}\sin 3\tau\Biggr{]}+O(\lambda^{4}),$$

$$\begin{array}{c}{A}_2^{\pm }=6\sqrt{\frac{2}{5}\left(\Delta K\pm \frac{{\gamma }_1}{2}\right)},\end{array}$$

$$\begin{array}{c}{B}_2^{\pm }=-\frac{3}{5{\left(A_2^{\pm }\right)}^2}-\frac{71}{720}{\left(A_2^{\pm }\right)}^2\pm \left(\frac{7}{15}{\gamma }_1-\frac{{\gamma }^2}{10{\gamma }_1}\right),\end{array}$$

$$\tau=\bigg{\{}\frac{\gamma}{2\lambda}-\frac{\lambda}{\gamma}\frac{\left(A_{2}^{\pm}\right)^{2}}{18}\\ -\frac{\lambda^{3}}{\gamma^{3}}\left[\frac{1}{4}+\frac{\left(A_{2}^{\pm}\right)^{2}}{18}\left(3B_{2}^{\pm}-\frac{\left(A_{2}^{\pm}\right)^{2}}{9}\right)\right]+O(\lambda^{5})\bigg{\}}t,$$

$${\gamma }_1=\sqrt{1-{\gamma }^2}.$$

$${\Lambda }_1^{in}$$

$${\Lambda }_2^{in}$$

Full text

I Introduction

The study of collective dynamics in networks of coupled oscillatory elements is one of the most attracting topics in modern nonlinear dynamics. It is important both for theoretical and practical points of view. A general phenomenon of collective behavior is synchronization Pikovsky ; Osipov ; Afraimovich ; Mosekilde ; Anishchenko ; Balanov . Synchronization is usually understood as a process of achieving the collective rhythm of functioning by coupled objects of different nature. Synchronization of two or three elements is possible, as well as of ensembles consisting of hundreds and thousands of elements Pikovsky ; Osipov ; Afraimovich ; Mosekilde ; Anishchenko ; Balanov . Even a weak attracting coupling can adjust phases and frequencies of oscillators, and they can synchronize. Now three types of synchronization in its networks are known: full (or global) synchronization, partial (or cluster) synchronization and chimera states.

The system of coupled pendulums is one of the widely used models in multiple fields of science and technology. Despite the simplicity of this model, it adequately describes not only mechanical objects Kecik , but also various processes that occur in semiconductor structures Barone . This model is often considered as the basis for the theoretical studies of coupled Josephson junctions Pikovsky ; Barone . Hence, it is very important to investigate the behavior of this system.

Cluster and solitary states in ensembles of coupled elements are of special interest for an investigation of the synchronization and symmetry breaking phenomena PikovskyRosenblum2015 ; Kurths2016 . The first type states consist of two or more groups in which an individual oscillators behave identically. These states have been well known for many years, but have still attracted great attention of investigators across different fields of science and engineering PikovskyRosenblum2015 ; Kurths2016 . Cluster states arise both in ensembles with a finite number of elements Kaneko1990 ; Okuda1993 ; Nakagawa1994 ; Schmidt2014 , and in distributed oscillatory media Vanag2000 ; Mikhailov2006 ; Lin2004 . In the article Kemeth2019 two-cluster regimes in small ensembles of Stuart-Landau oscillators are considered. Solitary state can be attributed to a special type of cluster state. Solitary state is formed when a single oscillator is separated from a synchronous cluster in an ensemble Maistrenko2014 . Examples of solitary states were found in various types of networks Mikhaylenko2019 ; Majhi2019 , systems of non-locally coupled elements Jaros2015 , systems of Stuart-Landau oscillators Schmidt2014 .

In this paper we consider features found in rotational dynamics of ensemble of globally coupled pendulums. Note that one can interpret our basic model as a generalization of the standard Kuramoto model, which takes into account the inertia and the intrinsic nonlinearity of the elements of the discussed population. The influence of the first type effects have been intensevely studied in recent years (e.g., see Ji ; Ha ; Belykh2016 ). Actually, the modifications of the second type have been also considered in the literature Komin ; Daido ; Lafuerza . In this work, we assume, that both effects are present and play an important role. Hence, the Kuramoto model describing the evolution of a group of phase oscillators with global coupling transforms to the system, consisting of a number of pendulums with meanfield interaction. We are interested in in-phase rotations and nontrivial out-of-phase ones. In the present article, we restrict ourselves to considering ensembles with a small number $N$ of elements and focus on all the possible limit motions of the system.

This paper is organized as follows. In Section II.1, we describe the model, state the problem. In Section II.2 we report on the numerically observed effect: in-phase periodic motion instability. Using the assumption of small dissipation we developed an asymptotic theory, which explains instability of the in-phase periodic motion (essentially presented limit cycle on the cylinder) of the pendulums. Here we also find analytical expression for the boundaries of the in-phase limit rotation mode instability interval regarding the coupling strength. During the nonlinear stage of this instability a periodic out-of-phase rotation emerges, in particular, a solitary state for which the phases of some pendulums coincide, while the phases of the other pendulums differ from the rest. Section II.3 contains a detailed description of the general analytical approach for stability analysis of out-of-phase periodic limit cycles existing in the ensemble of globally coupled pendulums. In Section III theoretical and numerical results are presented. In particular,in Section III.1 the main results for stability analysis of out-of-phase rotation modes within the framework of the considered model are presented. In Section III.2 bifurcations that lead to the appearance and disappearance of the out-of-phase limit rotation modes are analyzed. Bistability of the in-phase and out-of-phase limit periodic modes is established for the system under study. In Section III.3, scenario of chaotic rotational dynamics emergence is described. In Section IV, examples of out-of-phase regimes, existing in larger ensembles, are considered. A summary of the main results can be found in Conclusion. Appendix A contains the details about the linear stability analysis for out-of-phase rotation modes. In Appendix B we present a brief description of the numerical methods used for calculating any possible periodic modes and their linear stability within the framework of the considered model.

II Nature of symmetry breaking in a system of coupled pendulums

II.1 Model under study

We study the rotational dynamics of an ensemble of $N$ globally coupled identical pendulum-like oscillators (index $n$). The setup we employ can be considered as a generalization of the standard Kuramoto model, which takes into account the inertia and the intrinsic nonlinearity of the structural elements. The nonlinearity in the form of sinus was chosen because it is rather typical for many real systems indicated above. In this case, the evolution of the generalized coordinate $\varphi_{n}(t)$ of the $n$-th unit belonging to such a community is given by

[TABLE]

where $\lambda$ is the damping coefficient responsible for all the dissipative processes in the system, $\gamma$ is a constant external force identical for all pendulums, $K$ characterizes the strength of global coupling between the oscillators in the discussed group.

Since we are interested in the rotational dynamics of system (1), we begin with considering the simplest in-phase rotational mode and its stability. When at any given time $t$ coordinates $\varphi_{n}(t)$ and velocities $\dot{\varphi}_{n}(t)$ $(n=1,\dots,N)$ coincide, the system demonstrates in-phase dynamics, i.e. at any time $t$: $\varphi_{1}\left(t\right)=\dots=\varphi_{N}\left(t\right)=\phi\left(t\right)$. We denote such a regime as $(N:0)$. In this regime all pendulums move synchronously and their dynamics is described by the following single equation:

[TABLE]

which essentially represents a nonlinear damped pendulum equation with a constant external force. The dynamical behaviour, all possible equilibrium states and limit motions of a particle that is governed by Eq. (2) have been well studied and thoroughly analysed in many publications (e.g., see Andronov ; Tricomi1933 ; Belykh1977 ; Peng for details). Based on the previous results, here we briefly outline the main specific features of the model (2), that make the formulation of the problem investigated in the presented paper clear and accessible to the broad readership. The parameter plane $\lambda$, $\gamma$ is divided into three domains Tricomi1933 ; Belykh1977 ; Peng ; Strogatz , which correspond to different structurally stable cylindrical phase spaces of the system (2). For $\lambda$, $\gamma$ from the first domain, only two steady states, namely, a saddle and a focus (node), exist in the phase space of the corresponding pendulum. When $\lambda$, $\gamma$ are from the second domain in addition to these steady states, there is also a stable $2\pi$-periodic in $\phi$ limit cycle. Noteworthy, in this case, the attraction region of a focus (node) is small enough and is delimited by separatrices of a saddle. For $\lambda$, $\gamma$ from the third domain, equilibrium states disappear from the system, and there remains only an attractive rotation mode. The third domain is separated from the others by the straight line $\gamma\!=\!1$ in the parameter plane $\lambda$, $\gamma$. In turn, the first and second domains are separated by the so-called Tricomi bifurcations curve Tricomi1933 ; Belykh1977 ; Peng ; Strogatz .

In the case of small dissipation (i.e. when $\lambda\!\ll\!\gamma$), it is possible to analytically describe a rotational limit motion of one pendulum, using an asymptotic approach based on, e.g. the Lindstedt-Poincaré method Nayfeh . It has been shown in Bolotov2019 that an approximation for limit rotation solution of the system (2) has the following form:

[TABLE]

We note, the system (2) has a continuum of oscillatory (not rotatory) type solutions at $\lambda=0$ and $\gamma\in[0,1)$. If $\lambda$ is not equal to 0, the stable steady state (center) is transformed to the stable focus and all these solutions do not exist and only one periodic rotatory solution can appear.

II.2 Stability of in-phase rotation mode

The in-phase rotation mode is stable in the wide range of the parameters $\lambda$, $\gamma$ and $K$, according to our numerical simulations directly in the framework of the basic model (1). However, for certain values of $\lambda$, $\gamma$ and $K$ the system (1) demonstrates non-trivial behaviour: the instability of in-phase rotation limit motion is developed and, as a result, one of a variety of out-of-phase states appears, that leads to partially broken symmetry in the considered ensemble of globally coupled pendulums. It is worth mentioning that similar in nature effect takes place for the system of only two elements Smirnov2016 . This effect is also observed in the case of a small chain of nearest-neighbour interacting identical pendulums Bolotov2019 .

Before proceeding to a detail description of the variety of out-of-phase rotation modes observed in direct numerical simulations within the model (1), we propose an analytical approach, which permits one to explain an instability of periodic limit cycles existing in the system of globally coupled pendulums. As a starting point, we consider a possibility of developing of such an instability process for the in-phase state discussed above. This allows us to demonstrate the basic idea of our theoretical analysis while avoiding cumbersome calculations.

Let us find the stability conditions for the in-phase rotational mode. First we linearize the system (1) around $\phi(t)$, then ${\varphi}_{n}(t)=\phi(t)+\delta{\varphi}_{n}(t)$ and $\dot{\varphi}_{n}(t)=\dot{\phi}(t)+\delta\dot{\varphi}_{n}(t)$, where perturbations $\delta{\varphi}_{n}(t)$, $\delta\dot{\varphi}_{n}(t)$ are small in magnitude. Next we get the corresponding equations for variations $\delta{\varphi}_{n}(t)$:

[TABLE]

To continue with Eqs. (4), we introduce new variables

[TABLE]

for which Eqs. (4), corresponding to the in-phase $T$-periodic rotational cycle of the system (1), take the following form

[TABLE]

[TABLE]

The Floquet multipliers of Eq. (6) are equal to $1$ and $e^{-\lambda T}$. Indeed, differentiation of Eq. (2) with respect to $t$ shows that periodic function $\dot{\phi}\left(t\right)$ is the solution to Eq. (6), so the corresponding multiplier is equal to 1. The second multiplier can be found using Liouville’s formula. Therefore $\eta$ mode is always stable, so instability in the system can arise only due to the excitation of modes $\xi_{n}$ $(n=1,\dots,N-1)$. Eq. (7) belongs to the Mathieu-type equation. Hence, the parametric instability effects can be observed for some values of the parameter $K$ depending on $\lambda$ and $\gamma$ McLachlan . To find the boundaries of the instability domain of the in-phase rotation mode, we determine the coupling parameter $K$ values for which the Eq. (7) admits a solution with period $2T$ or, equivalently, with frequency $\omega/2$. Passing to a dimensionless time $t=\omega\tau$, we reduce the system (7) to the following form

[TABLE]

Using the perturbation theory, taking result (3) and searching for a solution to Eq. (8) with $\omega/2$ frequency, we get boundaries $K_{1,2}$ for the instability domain

[TABLE]

When the critical values of $K_{1,2}$ are reached, the instability of the $\xi_{n}$, $n=1,\dots,N-1$ modes simultaneously develops, which is a common property at destabilization of the synchronous states in globally coupled systems of $N$ identical units. Thus, in the system of globally coupled pendulums (1) only one instability region ($K_{1}<K<K_{2}$) of the in-phase rotational motion $\phi(t)$ exists, which boundaries $K_{1}$ and $K_{2}$ are determined by the expression (9).

Fig. 1 shows the width of the range of values of coupling strength $K$, when in-phase mode instability is observed, depending on the values of parameters $\gamma$ and $\lambda$ as a function $[K_{2}-K_{1}](\lambda,\gamma)$. The parameter plane $(\lambda,\gamma)$ can be divided into 3 areas: $L_{1}$, $L_{2}$ and $L_{3}$. In the area $L_{1}$, the in-phase mode is stable for any $K$. In the area $L_{2}$, bounded by the Tricomi curve and line $\gamma=1$, this regime does not exist. The Tricomi curve corresponds to the values of the parameters $(\lambda,\gamma)$, for which the rotational motion in the system (2) disappears as a result of the saddle-node homoclinic bifurcation. In the color region $L_{3}$ there is a range of values of $K$ for which the in-phase rotation is unstable ($K_{1}<K<K_{2}$). The boundaries $K_{1}$ and $K_{2}$ of this range are determined by expressions (9). Black color means, that the length of the instability interval of $(N:0)$ regime is close to one. Orange color means, that length of the instability interval of $(N:0)$ regime is close to zero. By the fixed $\lambda$ with increase of $\gamma$ the interval $K_{2}-K_{1}$ becomes smaller. The region of instability arises at values $\gamma$ close to unity and becomes wider as $\gamma$ decreases. When the parameters $(\gamma,\lambda)$ reach the Tricomi curve, the in-phase rotational mode disappears.

II.3 General approach for stability of analysis of symmetry broken states

Let us proceed to a detailed description of the general analytical approach for stability analysis of out-of-phase periodic limit cycles existing in the model (1) of $N$ globally coupled pendulums. Just for definiteness assume that pendulums form $M$ synchronous clusters with $N_{m}$ pendulums in $m$-th cluster ($N_{1}+\ldots+N_{M}=N$). Here and below such regimes will be denoted as $\left(N_{1}:N_{2}:\ldots:N_{M}\right)$. Note that a cluster also may consist of only one element. Assuming what every pendulum in the $m$-th cluster is moving with the same phase $\phi_{m}\left(t\right)$ and velocity $\dot{\phi}_{m}\left(t\right)$ the system (1) of $N$ second order equations reduces to the set of $M$ second order differential equations which defines dynamics of the unknown phases $\phi_{m}\left(t\right)$:

[TABLE]

Applying the same approach described in Section II.1 for the in-phase solution of (1) to the current problem of finding $T$-periodic limit cycles of the Eq. (10) the phases $\phi_{m}\left(t\right)$ can be determined.

After phases $\phi_{m}\left(t\right)$ are found it is possible to proceed with stability analysis of the regime $\left(N_{1}:N_{2}:\ldots:N_{M}\right)$. For that one imposes a linear perturbation $\varphi_{mk}\left(t\right)=\phi_{m}\left(t\right)+\delta\varphi_{mk}\left(t\right)$ and $\dot{\varphi}_{mk}\left(t\right)=\dot{\phi}_{m}\left(t\right)+\delta\dot{\varphi}_{mk}\left(t\right)$, where $\varphi_{mk}\left(t\right)$, $\dot{\varphi}_{mk}\left(t\right)$ are the perturbed phase and velocity of $k$-th pendulum from the $m$-th cluster ($k=1,\ldots,N_{m}$). The equations for variations $\delta\varphi_{mk}\left(t\right)$ thus read

[TABLE]

The next key step is the linear change of variables which breaks system (11) into a set of independent differential equations systems of smaller dimensions:

[TABLE]

Excitation of mode $\eta_{m}$ in the case $\xi_{\tilde{m}\tilde{k}}=0$ can be interpreted as perturbation of $m$-th cluster like whole unit, because in this case all pendulums forming cluster $m$ are identically perturbed. Stability associated with such kind of perturbations are often called internal stability LyapunovExponents . On the contrary excitation of modes $\xi_{mk}$ for fixed $m$ in the case $\eta_{\tilde{m}}=0$ means arbitrary perturbation of pendulums only inside the $m$-th cluster and associated stability is called external stability LyapunovExponents .

For new variables Eqs. (11), corresponding to the $T$-periodic rotational cycle of the system (1), take the following form

[TABLE]

[TABLE]

So actually it is necessary to analyze the system of $M$ second order equations (13) and $\tilde{M}$ independent equations of second order (14) (because of their similarity for different $k$, so index $k$ will be omitted below), where $\tilde{M}=\left|\left\{m:N_{m}>1\right\}\right|$ is the number of clusters in which the number of elements is greater than one.

Floquet’s Theorem implies that to each characteristic exponent $\Lambda$ of the linear differential equations system with periodic coefficients, such as Eqs. (13) is, there corresponds a solution of the form $\eta_{m}(t)=e^{\Lambda^{in}t}\tilde{\eta}_{m}(t)$ ($m=1,\ldots,M$), where $\tilde{\eta}_{m}(t)$ is periodic in $t$: $\tilde{\eta}_{m}\left(t+T\right)=\tilde{\eta}_{m}\left(t\right)$. New variables $\tilde{\eta}_{m}(t)$ satisfy the equations:

[TABLE]

To determine the characteristic exponent $\Lambda$, we use the Lindstedt-Poincaré method. First, since our theory is being developed near the conservative limit, i.e. for small dissipation $\lambda$, the characteristic exponent $\Lambda$ is represented in the form of the asymptotic expansion

[TABLE]

Coefficients $\Lambda^{(j)}$ should be chosen so as to avoid secular terms in (15), i.e. so that $\tilde{\eta}_{1},\ldots,\tilde{\eta}_{M}$ are periodic.

Equivalently, using Floquet’s Theorem we are searching for the solution of (14) in the form $\xi_{m}(t)=e^{\Lambda^{ex}t}\tilde{\xi}_{m}(t)$, where $\tilde{\xi}_{m}(t)$ is periodic, i.e. $\tilde{\xi}_{m}\left(t+T\right)=\tilde{\xi}_{m}\left(t\right)$:

[TABLE]

Using asymptotic expansion (16) and cancelling the secular terms in expression for $\tilde{\xi}_{m}$ unknown coefficients $\Lambda^{(j)}$ can be found.

The signs $\Lambda^{in}$ and $\Lambda^{ex}$ determine the stability of the modes $\eta_{m}$ and $\xi_{m}$ respectively. Thus obtaining the expression for $\phi_{m}(t)$ as an asymptotic expansion, similar to (3) in the case of in-phase rotations, we can get the decomposition in the form of a series for the Lyapunov exponents $\Lambda^{in}$ and $\Lambda^{ex}$ for the systems (15) and (17). In the end we have $2M$ exponents $\Lambda^{in}$ determining internal stability and $2(N-M)$ exponents $\Lambda^{ex}$ determining external stability. An example of symmetry broken states stability analysis is presented in Sec. III.1.

III Theoretical and numerical results for the system of three elements

III.1 Regimes (2:1) and (1:1:1) stability analysis

In this section, we consider the application of the approach developed in the Sec. II.3 for analyzing the stability of cluster rotational motions resulting from the development of in-phase mode $\phi(t)$ instability for $N=3$. In this case, there are two types of cluster modes: (2:1) solitary state and (1:1:1) regime. Detailed stability analysis of these regimes is presented in Appendix A.

It is shown in A.1 that for the solitary state (2:1) solutions $\phi_{m}^{-}$ arising on the right side $K_{2}$ of the instability interval of parameter $K$ defined by Eq. (9) is always internally unstable for all values of $K$ (see Eq. (27)). Conversely, solution $\phi_{m}^{+}$ arising on the left side $K_{1}$ of the instability interval is always internally stable, but changes its external stability with $K=K_{c}^{(2:1)}$, and it is stable for $K>K_{c}^{(2:1)}$, where

[TABLE]

Similarly for (1:1:1) regime from (35) the instability of the solution $\phi_{m}^{-}$ for any value of $K$ follows. Solution $\phi_{m}^{+}$ changes its stability with $K=K_{c}^{(1:1:1)}$, and it is stable for $K<K_{c}^{(1:1:1)}$, where

[TABLE]

Thus, for certain values of the coupling strength $K$ the synchronous in-phase rotational mode in the ensemble of three coupled pendulums becomes parametrically unstable and, depending on the initial conditions, regimes (2:1) or (1:1:1) can be realized. The values of coupled strength parameters $K_{c}^{(2:1)}$ and $K_{c}^{(1:1:1)}$ at which the change of stability of cluster modes occurs were obtained. From numerical simulations it is known that $K_{c}^{(2:1)}\leq K_{c}^{(1:1:1)}$ (see Section III.2).

III.2 Regular dynamics of in-phase and out-of-phase rotational modes

In this and following section, we present the results of the numerical simulations (for details see Appendix B) which are performed directly within the framework of the discussed model (1) for a wide range of the parameters $\lambda$, $\gamma$ and $K$ in the case $N=3$. First of all, we consider in detail the development of the self-induced parametric instability of the in-phase synchronous regime and focus our attention on the nonlinear stage of this process and the resulting movements that can be set over long time. Our numerical calculations employed a commonly used fifth-order Runge-Kutta scheme (with fixed time step $dt=0.001$) to integrate the system (1).

Let us consider the case $\gamma=0.97$. The bifurcation diagrams of evolution of periodical motions of the system (1) are shown in Fig. 2. The diagram shows the dependence of synchronism characteristics $\Xi$ from magnitude of the coupling strength $K$. First let us describe the diagram for $\gamma=0.97$, $\lambda=0.2$ (Fig. 2a). The horizontal segments $A_{1}$, $A_{3}$ correspond to the stable synchronous in-phase regime ($\Xi=0$). There is a region $A_{2}$ of the parameter $K$ values, when this regime becomes unstable. As shown above, in the course of the asymptotic consideration (the expression (9)), it is for the values of the coupling parameter $K$ that the parametric instability of the in-phase periodic motion develops from these intervals.

Let us consider processes occurring in an ensemble when $K$ takes values from the $A_{2}$, and values outside it. As the parameter $K$ increases, the in-phase periodic motion undergoes period doubling bifurcation ($K\approx 5.764$), while from the stable $2\pi$-periodic in $\varphi=(\varphi_{1},\varphi_{2},\varphi_{3})^{T}$ rotation the stable (1:1:1) $4\pi$-periodic regime (branch $B_{1}$) and unstable $4\pi$-periodic motion (2:1) (branch $B_{4}$) are generated, and $2\pi$-periodic regime loses stability. In addition to the two above described $4\pi$-periodic motions, there are also two unstable out-of-phase $4\pi$-periodic (1:1:1) and (2:1) rotations in $\varphi$ (branches $B_{3}$ and $B_{6}$, respectively), which are generated from an unstable $2\pi$-periodic motion as a result of a subcritical period doubling bifurcation ($K\approx 6.008$) with increasing $K$. Corresponding values of the parameter $K$ for which doubling bifurcation occurs calculated analytically using Eq. (9) are $K\approx 5.764$ and $K\approx 6.007$.

The (1:1:1) regime (Fig. 3a) corresponds to the branches $B_{1}$, $B_{2}$, $B_{3}$, and the (2:1) solitary regime (Fig. 3b) corresponds to the branches $B_{4}$, $B_{5}$, $B_{6}$. As the parameter $K$ increases the stability change of the (1:1:1) and (2:1) regimes is to see. The first one loses stability at $K_{c}^{(1:1:1)}\approx 5.9699$, and the second one becomes stable at $K_{c}^{(2:1)}\approx 5.9687$ (according to theoretical predictions, see Eqs. (18) and (19), $K_{c}^{(1:1:1)}\approx K_{c}^{(2:1)}\approx 5.962$). So there is a bistability region of (1:1:1) and (2:1) regimes. The stability change occurs through the pitchfork bifurcation, and herewith the unstable (1:1:1) regime appears (brunch $B_{7}$). Further with the increasing $K$ the periodic motions, corresponding to (1:1:1) regime (branches $B_{2}$ and $B_{3}$) and (2:1) regime (branches $B_{5}$ and $B_{6}$), merge and disappear resulting from the saddle-node bifurcation at $K\approx 7.023$ and $K\approx 7.042$ respectively. For given values of the parameters, the obtained numerical results are in good agreement with the analytical theory described in Sec. III.1.

On Fig. 2 it is to see that with the increasing values of $\lambda$ the stability region length of (1:1:1) regime becomes longer and the length of stability region for (2:1) regime becomes shorter.

III.3 Chaotic dynamics

Next let us consider the case with a lower value of the external force $\gamma=0.8$. Then at lower values of $\lambda$ the bifurcations of $2\pi$- and $4\pi$-periodic regimes qualitatively coincide with the above described $\gamma=0.97$. With the further increasing damping parameter from branch $B_{5}$ of the solitary regimes motions with bigger periods appear resulting from doubling period bifurcations, and it leads to the appearance of chaotic dynamics. Let us describe the same scenario for the case $\lambda=0.3$ (Fig. 4). Here at the increasing $K$ the (2:1) solitary regime (branch $B_{5}$) loses stability (branch $B_{8}$) at $K\approx 2.507$ resulting from period doubling bifurcation. Further after cascade of period doubling bifurcations the (2:1) solitary regime becomes chaotic.

Branch $B_{9}$ corresponds to the stable (2:1) regime. Similarly at a lower $K$ parameter at $K\approx 4.063$ the regime loses stability due to period-doubling bifurcation, what leads to the appearance of chaotic attractor. That exists at $2.7\lesssim K<2.754$. At $K>2.754$ the in-phase rotation becomes stable again.

Next at the increasing $\lambda$ the length of the stable branch of the (2:1) solitary state $B_{5}$ decreases, and it disappears. Fig. 5 shows the bifurcation diagrams of synchronous regimes for $\lambda=0.35$. In this case chaotic solitary regime (Fig. 6) is realized resulting from the cascade of period doubling bifurcations at $1.7<K<1.77$ and $3.0<K<3.01$. However, at $1.77<K<3.01$ the chaotic motion is unstable because of the chaotic attractor crisis, and in-phase rotation is to see.

The next case is for $\lambda=0.5$. At the increasing $K$ (1:1:1) motion (branch $B_{1}$) loses stability ($K\approx 0.618$) resulting from the Neimark-Sacker bifurcation, herewith quasiperiodic (1:1:1) motion appears. At $K\approx 0.656$ this quasiperiodic regime becomes chaotic after torus destruction bifurcation (Fig. 8a). At $0.89<K<1.272$ there is branch $B_{10}$ of $4\pi$-periodic motions, which merges with branch $B_{2}$ at $K\approx 0.956$ and $K\approx 1.272$ resulting from pitchfork bifurcations. At $0.89<K<0.907$ this branch has a stable region which corresponds to "stability window" when the chaos in system (1) is not observed (Fig. 7b).

Further at the increasing $K$ the regime is realized, when time intervals, at which pendulums’ phases $\varphi_{1}\approx\varphi_{2}\approx\varphi_{3}$ alternate with intervals, where $\varphi_{1}$, $\varphi_{2}$ and $\varphi_{3}$ do not coincide (Fig. 8b), i.e. there is an intermittency of chaotic oscillations (3:0) and (1:1:1). If $K\approx 1.007$, the in-phase regime becomes stable and chaotic oscillations are not realized, and herewith the system shows the in-phase rotations.

IV Rotational states in larger ensembles

In this chapter we present a variant of the development of the in-phase rotational motion’s instability for the system (1) at $N>3$. Let us consider the case $N=4$. The boundaries of the stability region of the in-phase rotation are defined by expressions (9) according to asymptotic theory developed in Section II for $\lambda\ll\gamma$. Next applying our analytical approach it can be shown there are three out-of-phase rotation motions: $(2:2)$ and $(3:1)$ two-clusters regimes and also $(2:1:1)$ cluster state. What is more according to developed theory this cluster state is stable for $K_{c1}^{(2:1:1)}<K<K_{c2}^{(2:1:1)}$, where

[TABLE]

For example let us numerically and analytically describe system (1) with parameters $\gamma=0.97$, $\lambda=0.2$ (Fig. 9). Direct numerical simulations show that for this values of coupling parameter $K$ in-phase periodic motion undergoes period doubling bifurcation. While out-of-phase (2:2), (2:1:1) and (3:1) regimes are generated (branches $B_{1}$, $B_{2}$ and $B_{3}$, respectively). Fig. 9(b) shows (2:1:1) regime, when two pendulums form an in-phase cluster and the others demonstrate out-of-phase dynamics. We note that in this case the range of parameter $K$ exists where out-of-phase (2:2), (2:1:1) and (3:1) regimes coexist and are stable, i.e. multistability is observed. In particular it is found from numerical calculations that for values $K_{c1}^{(2:1:1)}<K<K_{c2}^{(2:1:1)}$, where $K_{c1}^{(2:1:1)}\approx 5.9859$ and $K_{c2}^{(2:1:1)}\approx 6.0443$, (2:1:1) regime is stable. Corresponding theoretical values for the boundaries $K_{c1}^{(2:1:1)}$, $K_{c2}^{(2:1:1)}$ of the stability region following from the Eqs. (20) are $K_{c1}^{(2:1:1)}\approx 5.9414$ and $K_{c2}^{(2:1:1)}\approx 6.0022$.

We have analyzed the dynamics of ensembles with other number of elements. It was found that in the process of development of the in-phase ($N:0$) mode instability with an increase in the coupling strength $K$, a stable two-cluster state $(N_{1}:N_{2})$ is first observed. Noteworthy, in such case the number of elements $N_{1}$ and $N_{2}$ in each cluster is approximately equal to each other. This state becomes unstable with further increase of the parameter $K$, and other cluster modes are realized. In all cases the sequence of cluster mode type is finished by solitary (two-cluster) state $(N-1:1)$.

We demonstrate this observation for the case $N=7$ pendulums with parameters $\gamma=0.97$ and $\lambda=0.2$. Due to in-phase rotation mode instability development the stable (4:3) cluster regime arises (Fig. 10a). Further there is a sequence of stable cluster modes of other types, for example (3:2:2) (Fig. 10b). Solitary state (6:1) is observed last in a sequence of stable cluster modes (Fig. 10c).

V Conclusion

We have studied the dynamics of an ensemble of globally coupled identical pendulums. A relatively simple model demonstrates a big variety of regular and chaotic in-phase and out-of-phase cluster regimes. We found and theoretically approved self-induced parametric instability of symmetric in-phase rotation regime. It is shown, that in the system with the growth of coupling strength the generation of out-of-phase rotation periodic motions occurs resulting from period doubling bifurcation.

Note that there is one instability region, where different out-of-phase regimes can be realized. It describes an approach that allows to build asymptotic expansion for rotational modes and to investigate their stability in low dissipation limit. The developed method was applied in the case of a system consisting of three elements. The results of direct numerical modeling are well correlated with theoretical estimates of the boundaries of the existence and stability of rotational regimes, including solitary states, in the case of small dissipation. Bistability of in-phase and out-of-phase rotational periodic regimes can also be observed.

It is also demonstrated that at larger values of dissipation parameter chaotic regimes and intermittency can be realized. The several examples of cluster modes, as well as analysis of their stability, were demonstrated in the case of $N=4$ and $N=7$.

Acknowledgements.

Authors acknowledge A.O. Kazakov, G.V. Osipov and A.S. Pikovsky for valuable advices and fruitful discussion. Results presented in Section II were supported by the RFBR grant No. 17-32-50096. Results presented in Section III were supported by the RSF grant No. 19-12-00367. Results presented in IV were supported by the RFBR grant No. 19-52-12053.

Appendix A Regimes (2:1) and (1:1:1) stability analysis

A.1 Regime 2:1

Let us consider solitary state $(2:1)$, when two clusters exist: $\varphi_{1}(t)=\varphi_{2}(t)=\phi_{1}(t)$, $\varphi_{3}(t)=\phi_{2}(t)$. We construct an asymptotic expansion for $\phi_{1}(t)$ and $\phi_{1}(t)$ in the case of $\lambda\!\ll\!\gamma$. Applying coupling strength $K=\dfrac{\gamma^{2}}{4\lambda^{2}}+\Delta K$, where the first summand $\dfrac{\gamma^{2}}{4\lambda^{2}}$, to the first two leading orders, determines the middle of the coupling strength $K$ values range, for which the in-phase periodic motion becomes unstable (see Eq. 9), $\Delta K$ – deviation from this value. Considering that $\phi_{1}^{(0)}=\phi_{2}^{(0)}$ and assuming for simplicity initial condition $\left(\phi_{1}-\phi_{2}\right)^{\prime}|_{t=0}=0$, we construct two asymptotic solutions $\phi_{m}^{+}$ and $\phi_{m}^{-}$ using methods described above for the case $\lambda\ll\gamma$:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus there are two $(2:1)$ solitary states. Solution $\phi_{m}^{+}$ corresponds to the $4\pi$-periodic regime, arising as a result of direct period doubling bifurcation from in-phase rotation $\phi(t)$ for $K=K_{1}$. Similarly, solution $\phi_{m}^{-}$ corresponds to $4\pi$-periodic regime, arising as a result of inverse period doubling bifurcation for $K=K_{2}$.

Further, using the method described in Section II.3, we study the stability of the obtained solutions $\phi_{m}^{\pm}$. For the internal stability (see Eq. (15)) of the first solution $\phi_{m}^{+}$ we obtain

[TABLE]

Similarly, for the second one $\phi_{m}^{-}$:

[TABLE]

For the two-pendulums cluster external stability (see Eq. (17)) of the motion $\phi_{m}^{+}$ we have the following characteristic exponents:

[TABLE]

And for $\phi_{m}^{-}$:

[TABLE]

A.2 Regime 1:1:1

The methods used to investigate stability of the regime (1:1:1) are completely analogous. Therefore, this section contains only brief description of the main results about the stability of the (1:1:1) motion.

Expressing $K$ in term of $\Delta K$ again, considering that $\phi_{1}^{(0)}=\phi_{2}^{(0)}=\phi_{3}^{(0)}$ and assuming initial conditions $\left(\phi_{1}-\phi_{3}\right)^{\prime}|_{t=0}=~{}0$, we get two solutions $\phi_{m}^{+}$ and $\phi_{m}^{-}$:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

In the case of (1:1:1) regime there are only one-pendulum clusters. For this $\phi_{m}^{+}$ solution characteristic exponents are

[TABLE]

Respectively for $\phi_{m}^{-}$:

[TABLE]

Appendix B Methods for calculation of periodic motions and their stability

The theoretical analysis above allows us to describe the initial stage of the discussed instability of the synchronous rotation mode. One also can find all intervals of values of the coupling coefficient $K$, for which the development of the set-induced parametric instability is possible, and estimate the boundaries of these ranges with rather good accuracy. The direct numerical simulations of an initial value problem for the dynamical system (1) give us general ideas about the evolution in time of the ensemble of coupled pendulums and the nonlinear stage of the developed instability. In order to connect and complete these two pictures, we also identify periodic rotations and explore their parametric continuation within the framework of the model (1). To this end, taking into account that the $\varphi_{n}\left(t\right)$ is determined in the range from $-\pi$ to $\pi$ and using the property of closure of the considered trajectories in the phase space $\left\{\varphi_{n}\left(t\right),\dot{\varphi}_{n}\left(t\right)\right\}$, we construct the Poincaré map and employ the Newton-Raphson algorithm to find a fixed point there and a period $T$ of motion along a corresponding trajectory for each given set of parameters $\lambda$, $\gamma$ and $K$ BifurcationTheory .

The main idea of this method is as follows. Each of the solutions $\phi_{n}\!\left(t\right)$ we are interested in is primarily characterized by its period $T$ (which is strictly speaking unknown and is to be defined at the end of numerical computations) and the number $m$ of changes of phases $\varphi_{n}\!\left(t\right)$ by $2\pi$ during the period $T$. Hence, the Poincaré map \bigl{\{}\varphi_{n}\left(0\right),\dot{\varphi}_{n}\left(0\right)\bigr{\}}\to\bigl{\{}\varphi_{n}\left(T\right)-2\pi{m},\dot{\varphi}_{n}\left(T\right)\bigr{\}} has a fixed point corresponding to a trajectory \bigl{\{}\phi_{n}\left(t\right),\dot{\phi}_{n}\left(t\right)\bigr{\}}. Using this fact that $\phi_{n}\left(T\right)=\phi_{n}\left(0\right)+2\pi{m}$ and $\dot{\phi}_{n}\left(T\right)=\dot{\phi}_{n}\left(0\right)$, we construct the following system of equations

[TABLE]

where \bigl{\{}\varphi_{n}\left(t\right),\dot{\varphi}_{n}\left(t\right)\bigr{\}} is the solution to Eqs. (1) with initial conditions \bigl{\{}{\varphi_{0}}_{n},\dot{\varphi_{0}}_{n}\bigr{\}}, i.e. \bigl{\{}\varphi_{n}\left(0\right),\dot{\varphi}_{n}\left(0\right)\bigr{\}}=\bigl{\{}{\varphi_{0}}_{n},\dot{\varphi_{0}}_{n}\bigr{\}}. Therefore, a periodic solution with period $T$ of Eqs. (1) will be a root to (36). Because of the translational invariance symmetry (in time), we note that one value from the set \bigl{\{}{\varphi_{0}}_{n}\bigr{\}} can always be taken to be zero without loss of generality. We use the Newton-Raphson algorithm NumericalRecipes to approximate the roots of \mathbf{P}\left(T,\bigl{\{}{\varphi_{0}}_{n},\dot{\varphi_{0}}_{n}\bigr{\}}\right). It is also noteworthy that the Jacobian is $\widehat{\mathbf{J}}=\widehat{\mathbf{I}}-\widehat{\mathbf{Q}}\left(T\right)$, where $\widehat{\mathbf{I}}$ is the identical matrix and $\widehat{\mathbf{Q}}\left(T\right)$ is matrix obtained from the monodromy matrix $\widehat{\mathbf{M}}\left(T\right)$ (see its definition below) by replacing one of the columns by the vector of values of the right-hand sides of Eqs. (1) at the time $t=T$. As a result, we numerically obtain, with high precision, stable (above dynamically generated) and unstable rotational modes as exact time-periodic solutions of Eqs. (1). Continuing these solutions in value of the coupling strength $K$ within the interval of instability of in-phase rotational mode allows us to trace the entire family of nontrivial periodic motions and to analyze their bifurcations (see, e.g., Fig. 2) and study in detail their bifurcations and a process of transition to chaos. This is one of the main goals of the paper.

The linear (spectral) stability of the arbitrary ($2\pi$-, $4\pi$-, $8\pi$- and etc.) periodic motions (on the cylinder) of the dynamical system (1) is investigated by means of a Floquet analysis, chiefly relying on numerical calculations (see, e.g., BifurcationTheory ). To this end, we add a small perturbation $\delta\varphi_{n}\!\left(t\right)$ to a given periodic solution $\phi_{n}\!\left(t\right)$ of Eqs. (1). The linearized equations satisfied to the first order in $\delta\varphi_{n}\!\left(t\right)$ read:

[TABLE]

The Floquet analysis of Eqs. (37) can be performed due to the property of periodicity of the trajectory \bigl{\{}\phi_{n}\left(t\right),\dot{\phi}_{n}\left(t\right)\bigr{\}}. Therefore, the stability of the considered motion is defined by the spectrum of the Floquet operator (monodromy matrix) $\widehat{\mathbf{M}}\left(T\right)$ given by

[TABLE]

The eigenvalues $\mu_{n^{\prime}}=\exp\left(iq_{n^{\prime}}\right)$ (here and below $n^{\prime}=1,\ldots,2N$) of the matrix $\widehat{\mathbf{M}}\left(T\right)$ are dubbed the Floquet multipliers with being Floquet exponents $q_{n^{\prime}}$ of the periodic solution $\phi_{n}\!\left(t\right)$. Because of the damping and the external force, only one property (implied by the real character of the monodromy) can be extracted for $\mu_{n^{\prime}}$, which is that they are or appear in complex conjugated pairs. To examine the stability of each of the rotation motion under discussion, we compute their Floquet multipliers. If $\left|\mu_{n^{\prime}}\right|\leq{1}$ for all $n^{\prime}$, then the rotation mode is linearly stable. It is worth mentioning that one of the eigenvalues $\mu_{n^{\prime}}$ must be strictly equal to one, because we investigate the stability of a periodic motion. Hence, using this fact it is possible to check that the trajectory \bigl{\{}\phi_{n}\left(t\right),\dot{\phi}_{n}\left(t\right)\bigr{\}} we find numerically belongs to the family of periodic rotations (on the cylinder). If at least one of Floquet multipliers $\mu_{n^{\prime}}$ locates outside the unit circle in the complex plane, then the rotation mode is linearly unstable. Noteworthy, our calculations show that the most of the discussed periodic solutions $\varphi_{n}\left(t\right)$ exhibit a relatively strong instability due to a real Floquet multiplier $\mu_{n^{\prime}}$ with an absolute value greater than one, i.e. $\left|\mu_{n^{\prime}}\right|>1$.

As a characteristic of the degree of synchronization, we consider the value $\Xi$, which is the frequency lag of pendulums:

[TABLE]

where $T$ is the period of rotational mode. It follows from the definition (39) that $\Xi$ takes non-negative values, and $\Xi~{}=~{}0$ only in the case of in-phase mode. In the case of an out-of-phase regime, when there exists such a pair of pendulums that $\dot{\varphi}_{n_{1}}~{}\neq~{}\dot{\varphi}_{n_{2}}$, where ${n_{1}}$ and ${n_{2}}$ are the numbers of pendulums, $\Xi~{}>0$.