Cluster synchronization of dry friction oscillators

Synchronization is a well known phenomenon in non-linear dynamics and is treated as correlation in time of at least two different processes. In scope of this article, we focus on complete and cluster synchronization in the systems of coupled dry friction oscillators, coupled by linear springs. The building block of the system is the classic stick-slip oscillator, which consists of mass, spring and belt-mass friction interface. The Stribeck friction itself is modelled using Stribeck friction model with exponential non-linearity. The oscillators in the systems are connected in nearest neighbour fashion, both in open and closed ring topology. We perform a numerical study of the properties of the dynamics of the systems in question, in two-parameter space (coupling coefficient vs. angular excitation frequency) and explore the possible configurations of cluster synchronization.


Introduction
Synchronization is well known phenomenon in the nonlinear dynamics. One of the first scientific paper on that topic dates back to the 17th century when Christian Huygens observed mutual synchronization of two pendulum clocks hanging on the common support [1]. Synchronization can be encountered in the dynamical systems in various fields of science including, mechanics, electric engineering, biology [2,3]. It is defined as adjustment of rhythms of oscillating objects due to their weak interaction [4].
There are many types of synchronization. In scope of this article we examine complete and cluster synchronization. Complete synchronization means that the trajectories of all oscillators in the system converge. It can happen only in the case of identical oscillators. Cluster synchronization is similar, however here the system is not in complete synchronization but respective oscillators form a subsets of oscillators, which are in sync with each other. More formal definition is presented in Section 2.
Synchronous state stability is a feature of synchronization, which enable us to analyse the synchronization thresholds. Master stability function (MSF) proposed by Pecora and Caroll [5] is a powerful mathematical tool for for that purpose. In this method one can separate the network topology from the local dynamics of the system. The synchronous state stability is estimated by means of Lyapunov exponents of variational equation representing local dynamics of the system, for parameters, which correspond to eigenvalues of the connectivity matrix (see Section 2).
Here, a more direct estimation of the master stability function is presented, which is based on two-oscillator probe [6]. Two-oscillator probe is suitable for system, pos-⋆ e-mail: michal.marszal@p.lodz.pl ⋆⋆ e-mail: andrzej.stefanski@p.lodz.pl sessing only real eigenvalues of the connectivity matrix. There is no need to calculate the Lyapunov exponents. MSF is estimated by synchronization error for reference probe of two coupled oscillators (see Section 2). Such a procedure is particularly convenient for non-smooth systems (e.g. stick-slip oscillator), where it is numerically difficult to calculate Lyapunov exponents.
Friction is an ubiquitous force in physics and is the main source of the energy dissipation. The phenomenon of friction was investigated by famous scientists including: da Vinci [7], Euler [8], Coulomb [9]. In many engineering applications the classic Coulomb model is sufficient. Nowadays various mathematical models of friction have been developed, such as: [10][11][12]. More advanced models include differential equation representing the internal states of the model to fully describe the friction phenomenon.
In this paper, we presents findings of our research on synchronization properties of systems consisting of coupled dry friction oscillators. Such kind of oscillators belong to the class of non-smooth systems, whose dynamics is influenced by the discontinuity plane, corresponding to sticking phase. We also show the application of master stability function based on two-oscillator probe in systems with discontinuity. This article is based on authors' previous works [13,14] and consists of following parts. Section 1 gives brief introduction to the topic of synchronization and friction. Section 2 discusses the methodology used in the research in more details. In Section 3 the mathematical model of the system in question is explained. Section 4 presents the results of numerical investigation. Finally, the conclusions are drawn in Section 5.

Methodology
In this Section we present the methodology used in the research. Subsection 2.1 gives basic notion of synchronization, while the 2.2 introduces the concept of master stability function (MSF).

Synchronization
Let us consider a a following dynamical system: where x = (x 1 , . . . , x N ) ∈ ℜ N is a state vector and ). The second term in (1) describes the coupling properties. G is a connectivity matrix, H : ℜ N → ℜ N linking function, H a linking matrix, σ coupling coefficient, ⊗ denotes the Kronecker product of two matrices. For the general case the properties of G and H matrices can be arbitrary. Let us narrow our investigations to only two types of synchronization, namely complete (global) synchronization and cluster synchronization. We assume that the oscillators in the network are identical and so are the linking functions between them. Complete synchronization can be observed when two trajectories converge to the same value and later hold that condition [15]. In [16] complete synchronization is defined for two dynamics systems represented with their phase plane trajectories x (t) and y (t), respectively, when for all t > 0, the following relation is fulfilled: lim Complete synchronization requires the synchronization of oscillators both in phase and in amplitude. However, there may be a situation when the whole system is not in complete synchronization, but one can distinguish a subset of oscillators which are in sync with each other and out of sync with the members of the other subset. Subsets of synchronized oscillators are called clusters. It is important to mention that we can talk about cluster synchronization when whole system is not in complete synchronization. The motion of different clusters may be uncorrelated or one can observe a shift phase between them. The topic of clusters can be found in literature in [17][18][19].

Master stability function
Master stability function is a power mathematical tool used in determining the stability of the synchronous state [5]. The classic algorithm to estimate the MSF is to calculate transversal Lyapunov exponents (TLE) of the Eq. (6) derived below. Let us start with deriving a variational equation of Eq. (1): where ξ i is the variation of the i-th node, ξ = (ξ 1 , ξ 2 , ...ξ N ) is variation vector, DF is the Jacobian of any node, DH is the Jacobian of the linking function. Diagonalization of Eq. (3) leads to uncoupling the variational Eq.(3) into N block having a form of: where γ k is the k-th eigenvalue of the G, i = 0, 1, 2, ..., N − 1, ξ k is transverse mode of perturbation from the synchronous state. For k = 0, eigenvalue is γ 0 = 0, and consequently Eq. (4) can be reduced toξ associated with the longitudinal direction located within the synchronization manifold. The other k-th eigenvalues correspond to transverse eigenvectors [5]. The tendency to synchronization is a function of eigenvalues γ k . Let us substitute σγ = α + iβ in Eq. (4), where α and β are respective real and imaginary part of eigenvalues.
where ξ is an arbitrary transverse mode. Invariant synchronization manifold exists when the connectivity matrix G has zero sum rows [5]. All the real parts of eigenvalues, which correspond to transversal modes, are negative (Re(γ k 0 ) < 0). The spectrum of eigenvalues has the descending form, i.e., γ 0 ≥ γ 1 ≥ ... ≥ γ N−1 . In general case, [5] defines MSF as the largest transversal Lyapunov exponent λ T surface, computed basing on Eq. (6), on a complex numbers plane (α, β). If the interaction between the nodes be mutual (e.g. mechanical systems), then the eigenvalues have only real part and then MSF is represented only by a curve describing the largest TLE as a function of real number α, defined as The synchronous state of dynamical system is stable when all eigenmodes of the discrete eigenvalue spectrum σγ k lay in ranges of the largest negative TLE. Even if only one eigenvalue is in the range, where λ T > 0, the global synchronization is unstable, however, cluster synchronization is still possible.
In case of non-smooth dynamical systems, the computation of TLE requires is numerically troublesome and requires algorithms. Solutions to that problem are techniques called three-oscillator probe and two-oscillator probe [6,20]. Oscillator probe is based on estimating the MSF on the complex plane by direct detection of the complete synchronization by means of numerical integration or by experimental measurement. The methodology is simple and efficient. When calculating MSF in threeoscillator probe for system containing N oscillators, one initially investigates the reference probe of three oscillators. The area on the complex plane (α, β), where the complete synchronization or imperfect complete synchronization occurs, is the equivalent of the area of negative transversal Lyapunov exponents. The two-oscillator probe can be applied for mechanical systems, where due to mutual interaction between the nodes, the eigenvalues of the connectivity matrix G are real and the MSF is reduced to a curve. Should all non-zero eigenvalues of the connectivity matrix of the system of N oscillators be in the surface or range of complete (imperfect complete) synchronization for the reference probe, then the complete synchronization (imperfect complete

Mathematical model
Let us consider a system of coupled dry friction, selfinduced oscillators with additional kinematic excitation (see Fig. 1). The equation of motion in a non-dimensional form can be as follows: (8) where χ is a non-dimensional displacement, σ is a coupling coefficient, corresponding the strength of the coupling, G N is the connectivity matrix, ǫ normal load coefficient, ϑ r is a relative tangential velocity between contacting interfaces, u amplitude of excitation, ω angular excitation frequency, τ non-dimensional time.
Stribeck friction model with the exponential nonlinearity is applied as the basis for friction modelling in this work. The other models of friction may also be applied, however, they do not make significant difference from the point of view of the synchronization properties. The Stribeck friction model with the exponential nonlinearity is formulated as follows: where µ s -static friction coefficient, µ k -kinetic friction coefficient, a is constant defining the shape of the friction -relative velocity curve.

Results
In this Section let us show the result of two parameter study in (σ, ω) space, as well an example of MSF with two-oscillator probe.
The numerical simulations are based on author's own program written in C++. Following non-dimensional parameters are used in simulations: ϑ = 0.1, µ s = 0.35, µ k = 0.2, a = 2.5, u = 0.1. If different values are used, information is placed in figure caption or legend respectively. A transient time equal to 1000 excitation periods (∆τ t = 1000 · 2π/ω) is applied, which is followed by measurement of average synchronization error for time interval corresponding to 200 excitation periods. The investigated systems are started from the initial conditions, when the complete synchronization is slightly perturbed. Should the synchronous state be stable, the trajectories return to synchronous state after time ∆τ t . Fig.s 2 and 3 depict different behaviour of the investigated system for particular set of parameters. Black colour marks the complete synchronization areas. Red colour in Fig. 2 denotes the cluster synchronization zone. The cluster synchronization layout is presented on the right hand side of each sub-figure. In Fig. 3, which shows result for the closed ring topology, the variety of cluster synchronization configurations is richer. The cluster layouts are marked with respective colour. Both for open and close ring topology the complete synchronization region is located in the weak coupling zone. For stronger coupling the complete synchronization transforms to cluster synchronization eventually to disappear.
More detailed analysis of global (CS) and cluster synchronization thresholds is performed for certain selection of ω, including verification of the obtained results by means of MSF (see Fig.s 4, 5). The MSF is defined as average synchronization error for two-oscillator probe �e II �(α) as a function of real number α (see Equation (7)). Next, MSF �e II �(α) is projected via eigenvalues of connectivity matrix � G N onto bifurcation diagrams of average synchronization error for networks consisting of N oscillators with σ as a bifurcation parameter. The complete synchronization for network of N oscillators occurs, provided all eigenvalues spectra of connectivity matrix � G N lie within zero �e II (α)� function. The areas of zero MSF within the eigenvalues spectrum, as well as complete synchronization regions in networks of N oscillators are marked with grey colour respectively.
The method described above is robust for predicting the global CS thresholds and along with the analysis of eigenvectors can be used to explain the cluster synchronizability. However, due to additional coupling factors (i.e., excitation, friction or coexistence of the system attractors) the MSF method indicates only the tendencies of the oscillators to synchronize and might be not always verified in given network configuration.
for the networks consisting from three to six oscillators. For the case of three oscillators ǫ = 1.5, N = 3, ω = 1.6 ( Fig. 4) the CS region is in range σ ∈ [0.12, 0.27]. It is worth mentioning that the first and last oscillators in the network are in cluster synchronization for almost all investigated range of σ. Similar behaviour can be seen in Fig. 2a, marked as red region, wherein cluster synchronization occupies large area of the investigated parameter space. For the case of four oscillators (see Fig. 5     In the presented systems it is possible to observe the so called ragged synchronizability phenomenon [21]. In Fig.s 4

Conclusions
To sum up, we have performed a study of synchronization properties in arrays of coupled dry friction oscillators. Oscillators are connected in nearest neighbour fashion. The goal was to find synchronization thresholds in various networks of oscillators. The used methodology is based on master stability function. However, MSF is not estimated in a traditional way by TLE but by means of more direct approach, namely, two-oscillator probe.
One needs to bear in mind that MSF describes tendencies of the system to synchronize. It cannot be treated as the final condition for the synchronization. The oscillators in question are coupled also by common excitation and the friction force itself, which may contribute also to the synchronizability of the system. Velocity of the conveyor belt, which is equal for all oscillators as well as the same friction model, provide identity of the parameters, which is necessary condition for the occurrence of CS. Moreover, the common harmonic excitation correlates in time with the driving components of all oscillators and as a consequence facilitates the synchronization.
On the contrary, important factor leading to the desynchronization or appearance of cluster is the coexistence of attractors, which is characteristic and often encountered for the systems with friction and impact oscillators. Coexistence of attractors is a property of non-linear systems, which can occur also in smooth, time-continuous dynamical systems. Therefore, the MSF concept and eigenvectors analysis can be treated only as a tool for estimating the overall, global predisposition of the system of coupled oscillators to synchronize or to cluster. This may explain fact, that for some configuration in closed ring topology, some of the cluster layouts cannot be explained by the eigenvectors or eigenvalues interpretation. The results presented in this paper also show ragged synchronization phenomenon (i.e., complete synchronization windows).
In general, the synchronization stability criterion given by the MSF does not provide for proper detection of global network synchronization state even in the case of smooth systems described by continuous differential equations. The more this problem occurs in non-smooth systems where the structure of attractors coexistence and their basins of attraction is usually more complex than in smooth systems. Hence, on the basis of our research, we can conclude that for non-smooth dynamical systems the MSF estimated with two oscillators probe can be even more effective than one calculated with use of the TLE, because then we can be sure that the synchronous region was really detected and it is not only a projection of an interval of the negative TLE.