DYNAMICS OF DISORDERED HETEROGENEOUS CHAINS OF PHASE OSCILLATORS

We study the problem of robustness of synchronous states to disorder in the chain of phase oscillators with local coupling. The study combines a numerical determination of the existence and stability of synchronous states with an analytical investigation of the role of the phase shift and the level of disorder in the natural frequencies in the destruction of synchrony. We show that the presence of the phase shift facilitates robustness of the synchronous regime, at least up to its certain threshold value.


Introduction
Synchronization is a basic concept of rhythms adjustment of self-sustained periodic oscillators due to their weak interaction. This adjustment can be described in terms of phase locking and frequency entrainment. Synchronization phenomena in large ensembles of coupled systems often manifest themselves as collective coherent regimes appearing via non-equilibrium phase transitions. Despite a recent progress in studies of this phenomenon in a wide range of systems [Dörfler and Bullo, 2014], there are still several less elaborated problems, such as synchronization in disordered chains and lattices. Recent studies [Kogan, 2009;Lee, 2009] have been restricted to the simplest case of pure sine-coupling of phase oscillators, which is strongly dissipative. Another approach to disordered lattices has been developed in paper [Lapteva et al., 2015], which is based on the reformulation of the problem in the basis of linear Anderson (localized) modes with the main focus on weakly nonlinear regimes.
Multistability [Tilles et al., 2011;Girnyk et al., 2012;Niu, 2017], multiple basins of attractions [Delabays et al., 2017;Ha and Kang, 2012], and traveling waves [Hong and Strogatz, 2011] are some of the fundamental phenomena directly related to equilibrium propeerties of lattices of phase oscillators with both attractive and repulsive phase couplings. Among other papers we would like to mention very recent investigations [Matheny et al., 2019], where quasi-sinusoidal oscillators with linear nearest-neighbour coupling have been shown to manifest exotic regimes, such as splay states, inhomogeneous synchronization, clusters, and weak chimeras. It is also important to mention exact analytical results by [Mihara and Nedrano, 2019] for stability and the dynamics close to the q-twisted states in a finite lattice of symmetrically coupled identical phase oscillators.
In this study we go beyond these investigations mainly by exploring disordered chains as the simplest case of disordered lattices. Namely, we characterize synchronous states in disordered one dimensional lattices of coupled phase oscillators. We considered disorder in natural frequencies of oscillators, taking the phase shift in the coupling, which determine whether interaction is attractive or neutral, as the main parameter.
Previous studies, particularly the fundamental paper by Ermentrout and Kopell [Ermentrout and Kopell, 1984], describe the dynamics and the bifurcations of the synchronous state in the case of zero phase shift in details. The main goal of this paper is to extent this study to the case of a non-zero phase shifts. In the framework of this study we numerically determine the existence and stability of synchronous states. The strategy is to start with a vanishing disorder, where the synchronous states have just a regular smooth profile, and to follow numerically these stationary solutions by increasing the disorder level, for different sets of the frequencies, checking for stability of the obtained solutions.
The paper is organized as follows. In the Section 2 we state the problem and present the mathematical model of the system under study. In Section 3 we first of all present results for a chain with vanishing disorder, including analytic description for synchronous regime in case of smooth profile of the solution. Our main analytical tool for the study of such regimes for different values of the phase shift is the quasicontinuous approximation. Then, applying an iterative procedure based on the Newton method to calculate an exact fixed point of lattice equations, we numerically find the domain of the existence of the stationary solutions associated with synchronous states in the disordered chains with the phase shift. Moreover, here we confirm obtained results by direct numerical simulations within the formulated model of the chain of phase oscillators. In Section 4 we present numerical study of the most interesting cases of the dynamics for non-zero disorder in natural frequencies in the region, where according to our analysis, the stationary solution does not exist. A summary of main results and plans for further studies can be found in Conclusion.

The Model
Let us consider the lattice consisting of N phase oscillators with a nearest-neighbour coupling. In this case the evolution of the phase ϕ n of each unit is given by the following equatioṅ ϕ n = σω n +sin(ϕ n+1 −ϕ n −α)+sin(ϕ n−1 − ϕ n −α), (1) where normalized natural frequencies ω n are chosen from a continuous uniform distribution over a line segment [0, 1], parameter σ defines the level of disorder in natural frequencies, and the phase shift α determines whether the interaction between elements is attractive, repulsive or neutral. Below we assume that the coupling is strongly attractive, i.e. α ≤ 0.3. It is natural to set the boundary conditions as follows: which corresponds to the free boundaries of the chain, i.e. there are no elements with indexes n = 0 and n = N + 1. Actually, the system (1) can be interpreted as the Kuramoto -Sakaguchi model, which is relevant to many physical, chemical and biological systems, e.g. lasers, biocircuits, electro-mechanical oscillators [Matheny et al., 2019]. This model is the paradigmatic and universal object of study, which allows one to describe in details the phenomenon of synchronization.
Previous studies [Ermentrout and Kopell, 1984], devoted to the case of zero coupling phase shift, reveal several nontrivial features of the dynamics of the system. First of all, here up to certain values of disorder in natural frequencies of oscillators, one can find synchronous states in the described system. Our main goal in this paper is to investigate, how these synchronous states and their stability are changing with the increase in the governing parameter α.
In the further exploration of the behaviour of the system (1), we first of all analyze the solutions of system (3). Our strategy is as follows. We present the analysis of stationary solutions for the phase differencesv n in the system (5)  Here we focus on the stable stationary solutions in the discrete system of nonlinear equations (3). We start our analysis with the case of σ = 0. It is well known that under this condition system (3) has a stable stationary solutionv n = 0 for α = 0. Taking this fact as a basis of our study, we describe modification of this solution for quite small values of the phase shift α, and for vanishing disorder σ = 0. In the considered case, it is instructive to start our analysis with the quasicontinuous approximation of the system (5). In order to do this, we associate a discrete index n with a continuous variable x, introducing a spatial step h. Then we approximate the discrete operators in the left-hand side of (5) with spatial derivatives, substitutingv(x) andv(x ± h) instead ofv n andv n±1 , and expanding sinv(x ± h) and cosv(x ± h) as the Taylor series expansion up to the second order. As a result, we arrive to the following ordinary differential equation for the functionv(x): with boundary conditionsv(0) = −α andv(L) = α, where parameter L is the length of the continuous medium. Double integration of (6) yields an implicit function for variablev of the following form where β = 2 sin α − Q and γ = 2 sin α + Q. Here a constant Q can be found using the boundary conditions as a solution of the transcendental equation The examples of the solution profilesv(x) are shown in Figure 1a. Formally, solution (7) of Eq. (6) exists for all α ∈ [0, π 2 ). Equation (6) is the quasicontinuous approximation of (5) for σ = 0. However, we stress here that it is not an actual asymptotic version of the lattice problem (5), because once the spatial length is normalized using the lattice length, the small parameter (which corresponds to h above) is eliminated, and higher-order terms are small to the extend that higher-order gradients are small. In other words, approximation (6) is not based on a small parameter, because real spacing between the sites of the considered chain is one. Hence, the higher-order terms are in general of the same order as the lower-order ones. Thus, the validity of this approximation may be supported only by a comparison of its predictions with the numerical solution of the discrete system of nonlinear equations (5).
According to our numerical calculations based on the Newton method, the stationary solution of lattice equation (6) with profile close to the kink solution (7) with h = 1 and L = N actually does not exist for values of the phase shift α > α * . We clearly see from our numerical analysis that the critical value α * is determined by the degree of smoothness of the solution (7) in its central part close to x = L/2. The possible reason for it is as follows. As it was mentioned above, the quasicontinuous approximation is valid for smooth profile ofv(x) defined by (7). The increase in α above the threshold value α * leads to the case, when the derivative dv/dx at x = L/2 becomes comparable with the value of α/h (Figure 1b). Therefore the critical value of α * can be calculated from the following equation: 3.2 Numerical Continuation of the Stationary Solution to the Case of the Non-zero Level of Disorder Now let us study an evolution of the stationary solution of the system (5) in the presence of disorder, i.e. when σ > 0. For this purpose we solve system (5) for the synchronous state numerically. Starting from a smooth profile of the phase differences v n obtained for the case of vanishing disorder as outlined above, we apply an iterative procedure based on the Newton method to find an exact fixed point of these equations. The strategy is to start from solutions of (5) with σ = 0 and different values of α, where the shape of the synchronous state is known, and to change parameter σ gradually for fixed values of α and ω n , to remain in the convergence domain of the Newton method. In this way, the solutions can be found in a large domain of parameters α, σ, and its borderlines can be identified for multiple sets of natural frequencies ω n .
Next, we analyze temporal stability of the stationary solutions under discussion, for each set of natural frequencies ω n and different values of parameters α and σ. For this goal we linearize the lattice equations (3) with boundary conditions (4), which are characterized by the phase shift α and the level of disorder σ, near one of such a synchronous statev n . Then, after applying some necessary transformations to it, we finally analyze the eigenvalues of a matrix with heterogeneous coefficients.
Here we present some details of basic calculations of time stability of stationary solutions in the presence of disorder. The main steps of this procedure are as follows. To study stability of the synchronous state with certain values of α and σ and corresponding set of natural frequencies ω n we apply a substitution v n (t) =v n +ṽ n (t), whereṽ n (t) describes a small deviation from the profilē v n , and accounts for all possible perturbations around the given stationary profilev n . Substituting (10) into (3) and linearizing the results with respect to the variatioñ v n (t), we get a system of linear evolutionary equations forṽ n (t) in the following form: (11) Solution of equations (11) is sought in the form of v n (t) = a n e λt (12) leading to the eigenvalue problem λa n = cos(v n−1 +α)a n−1 −2 cos α cos(v n )a n +cos(v n+1 −α) (13) Therefore the stability properties of synchronous state can now be investigated by analyzing the spectrum of the time-independent eigenvalue problem (13).
In order to characterize the borderlines of existence and stability of the stationary solution, we calculatē v n (α, σ) for different sets of individual frequencies ω n , and for different values of the governing parameters α > 0, σ > 0 as a C 0 -smooth (with respect to parameters) continuation of the zero solution v n (0, 0) = 0, n = 1, . . . , N in case σ = 0, α = 0 for different sets of natural frequencies ω n . Results of these calculations are shown in Figure 2, where one can see several realizations of the borderline of stability of the stationary solution (light blue lines), and their mean value (single dark blue line) on the parameter plane α, σ, where each realization of the described borderline corresponds to a particular set of natural frequencies ω n . For a particular set of ω n , one can observe stationary solution for any pair of parameters α, σ taken from the region below this borderline. Letters A − E describe sets of parameters for which we present the dynamical regimes below. At parameter set A one can observe the dynamics typical for the stationary solution, while at parameter sets B -E different types of more complex dynamics are observed.

Numerical Studies of Dynamical Properties in the Disordered Chain
Now let us describe numerically typical dynamical regimes in the lattice within the framework of the initial model (3), for several values of the parameters α and σ inside and outside the region of existence of synchronous solution ( Figure 2). All the results presented here are obtained for one fixed set of ω n ; the dynamics of the system under study for other possible sets of natural frequencies is qualitatively similar. First of all, we note that direct numerical simulations confirm the presence of synchronous solution for vanish- Figure 4. The dynamics of a complex non-stationary solution of the system (1) for σ = 0.6, α = 0.1 (point B in Figure 2). (a) Time series for phases ϕ n and initial distribution of natural frequencies ω n , (b) snapshots for phase differences v n . Here in coloured regions, one may observe phase slips that are markers of loss of the stability of synchronous state and the occurence of more complex dynamics in the system. In non-coloured regions, the behaviour of oscillators is close to a synchronous uniform rotation with insufficient quasiperiodic deviations from given synchronous state, which are caused by phase slips in neighbouring regions.   ing disorder and the phase shift from the region α < α * . In Figure 3 one can see, after some transient process, that the phases of all units are locked together and the phase differences v n do not change in time (see Figure 3b).
Outside of the predescribed region, i.e. in the case of a sufficient disorder in the oscillators' natural frequencies, one can observe more complicated dynamical regimes in the system. In Figure 4 one can see, for parameter values α = 0.1 and σ = 0.6, a complex non-stationary regime with phase slips, which means the breakups in the synchronous dynamics of the chain (1). Far beyond the borderline of the existence of the stationary synchronous regime, e.g. for parameter values α = 0.1 and σ = 1, these phase slips become even more evident and rough (see Figure 5). With the increase in the phase shift parameter up to the value α = 0.3, we observe that the threshold value of σ, which is related to the first occurrence of the phase slips, increases. This fact is quite ev-ident from the comparison of pairs of Figures 4-6 and Figures 5-7.

Conclusion
In this paper we studied the robustness of the stationary solution in a disordered lattice of phase oscillators against the phase shift parameter and the level of disorder in the natural frequencies. We focused on the case of a smooth profile of the stationary solution. Under this assumption, we have shown that the phase shifts in a certain range make the regions of existence of the stable synchronous regime even wider in comparison to the case of zero phase shift. This result may be counterintuitive, because a phase shift is believed to be harmful for synchronization in this case.
Also, in this system we observe different patterns for different disorder realization. Of the main interest were stationary synchronous solutions. As we expected to see, these clusters dominate weakly asynchronous states for the phase shifts of the interaction close to zero. Dynamically, the simplest cluster states are periodic or quasiperiodic ones, while for a stronger disorder and a larger number of clusters, chaotic states are expected.
This study is the first part of intensive investigations, where the main goal is to find statistical properties of the interplay between structural (characterized via spatial configurations, e.g., clusters) and dynamical (characterized via Lyapunov exponents, correlations, diffusion properties) complexity. More specifically, we plan to study complex states that appear in different nonhomogeneous chains of phase elements. A detailed analytical and numerical study of such ensembles should reveal basic features, stability conditions and bifurcation scenarios of birth and death of complex states in such systems. Here the interdisciplinary approach, where a physical statistical analysis should be complemented by a mathematical bifurcation analysis, is essential.