Collective stochastic resonance behaviors of two coupled harmonic oscillators driven by dichotomous fluctuating frequency?

Lei Jiang(姜磊) Li Lai(賴莉) Tao Yu(蔚濤) Maokang Luo(羅懋康)

1College of Mathematics,Sichuan University,Chengdu 610064,China

2College of Aeronautics and Astronautics,Sichuan University,Chengdu 610064,China

Keywords: coupled harmonic oscillators, dichotomous fluctuating frequency, synchronization, stability,

1. Introduction

Stochastic resonance (SR) is one of several interesting phenomena, wherein the influence of noise is not restricted to destructive effects, but can have unanticipated ordered outcomes.[1–3]Counterintuitively, the output of such a noisy nonlinear dynamic system, driven by a weak periodic signal,can actually be amplified by a proper dose of noise. Three main ingredients(i.e.,nonlinearity,a weak periodic force,and a noise source)are required for the happening of SR.Early research of SR was confined to nonlinear systems driven by additive noise.[4–13]However, many recent investigations show that SR can also occur in linear systems with multiplicative colored noise.[14–24]

As a classical model,harmonic oscillators[25]are widely used to describe various natural phenomena. Thus, the research on harmonic oscillators have draw more and more attention. Especially, the SR of linear harmonic oscillators affected by complex environmental fluctuations has been widely investigated.[26–31]For example, Sauga[26]studied the resonant behavior of a fractional oscillator with fluctuating mass.Lang[27]investigated the SR of a linear system with random damping.Mankin[28]explored the SR of a harmonic oscillator with fluctuating frequency.

The above researches on SR mainly focused on uncoupled systems. However, in most stochastic dynamic systems,the particles are coupled in a complex network,[32]which means they may interact with each other. As the study of coupled system has continued to grow in importance and popularity, the SR investigation of coupled systems has become a hot research topic.[9–13,33]It was shown that the coupling can lead to the enrichment of the resonance behavior.[9,12,33]Pikovsky[12]demonstrated the existence of system size resonance in coupled noisy systems. The diversity-induced resonance in globally coupled bistable systems was researched in Ref. [9]. Nicolis[33]showed that nonlinear spatially coupled subsystems can induce coupling-enhanced SR.

In most previous studies, the impact of environmental fluctuations on each coupled particle is assumed to be identical,[17,19,31]which means the same noise will be added to each coupled oscillator. In realistic situations,although the coupled particles are in the same environment,the fluctuation noises acting on them are only independently and identically distributed,not exactly the same. Therefore,for coupled particles with independently and identically distributed noises,studying their collective SR behavior is an interesting problem.

As one of the most important noises, dichotomous noise has drawn great attention in theory and application. The widely application of dichotomous noise is due to its two properties. First,it is bounded,which makes it easy to be realized in physical and biological systems.[34–36]Second, dichotomous noise can be made by well-defined procedures to two limiting forms: Gaussian white noise and white shot noise,which is of great theoretical interest.[37]

In this paper, we study two coupled harmonic oscillators with different random frequency fluctuations, which are modeled as dichotomous noises. The two-coupled oscillators system can be regarded as a basic case of complex coupled systems,such as globally[13]and nearest-neighbor[33]coupled systems. Related studies[38–40]show that this kind of system can exhibit rich nonlinear dynamic behaviors and provide theoretical guidance for related investigations on complex systems. Through the Shapiro–Loginov formula[41]and the Routh–Hurwitz stability criterion,[42]we obtain the synchronization condition of the system. Based on this, the stability condition and the exact steady-state solution are derived.Finally,the collective behaviors of the coupled system are analyzed through theoretical and numerical methods, including stability,synchronization,and SR.

The remainder of our paper is organized as follows. In Section 2, we introduce the model of two coupled harmonic oscillators with fluctuating frequency. In Section 3,we derive the corresponding analytical results. The numerical simulation method is given in Section 4. In Section 5,the collective behaviors of the system are studied. The main conclusions of this paper are given in Section 6.

2. System model

2.1. The overdamped Langevin equation

In order to better explain the physical meaning of model(2),we first make a brief review of the related research on the overdamped Langevin equation

wherex(t)represents the position of the particle at timet.ωis the eigenfrequency, and its fluctuation is modeled as noise termξ(t).F(t)is the external force,usually taken as the sine(A0sin(?t))or cosine(A0cos(?t))periodic force.

Because of the wide application of Eq.(1)in many fields such as fluctuating barrier crossing,[43]enzimatic kinetics,[44]and nuclear magnetic resonance,[45]it has been extensively investigated. Berdichevskyet al.[46]firstly gave the analytical solution of Eq. (1) whenξ(t) is modeled as symmetric dichotomous noise, and analyzed the SR phenomenon of this model. After that, the authors in Refs. [16,47] extend the above results to the case of asymmetric dichotomous noise,and give some new applications of Eq.(1)in first-order linear circuit.

2.2. The overdamped coupled system

Most dynamic models in biology,physics,or engineering are complex systems composed of multiple interacting particles. And the coupling effect makes the system exhibit many new dynamic behaviors. However, due to nonlinearity and complexity, most complex coupled systems can only be analyzed by numerical simulation. As the simplest coupled case,the two coupled oscillators system can be regarded as a basic unit of complex coupled systems,such as globally and nearestneighbor coupled systems,and in many cases it can be solved theoretically. Thus,the study of such systems can provide theoretical guidance for the research of complex systems. Based on this,in this section we consider a coupled harmonic oscillators system consisting of two particles. The system is driven by a periodic cosine force. Random frequency fluctuation is modeled as the symmetric dichotomous noise.The system can be described as follows:

wherexi(i=1,2)represents the position of thei-th particle.ε≥0 is the coupling coefficient. Random frequency fluctuation is modeled as the symmetric dichotomous noiseξi(t)(i=1,2). Eachξi(t)(i=1,2)is a two-state Markov process with random values between±σ,which satisfies

where〈·〉represents statistical average,σis the noise amplitude andλrepresents the transition rate. Here,we assume thatξ1(t)andξ2(t)are independent of each other.

Due to the random fluctuation of noiseξi(t),the potential field of system(2)randomly switches between

This kind of fluctuating potential energy can be found in many systems. For example, in an ATP-ADP cycle system, ATP stores energy in its bond. By hydrolyze this bond,some electrons in the bond can enter into a lower energy state. That will result with some energy and ADP being put into the biological organism. Conversely, ADP can also generate ADT by using some energy.[48–50]Thus,system(2)can be viewed as the basic case of coupled protein molecular motor in the ATP-ADP potential field driven by a periodic cosine force.[31]Because in biological organism,ATP-ADP conversions do not all happen at the same time. Therefore,unlike Ref.[31],we assume that the fluctuationsξi(t)(i=1,2)of the potential field are independent of each other. The authors in Ref.[19]proposes some potential applications of this type of coupled system.[51–54]

3. Theoretical analysis and results

3.1. Synchronization

In the following theoretical derivation part, the Shapiro–Loginov formulas will paly an important role,so we list them here

In order to study the collective behaviors of the coupled system,it is necessary to analyze the synchronization betweenx1andx2,for which the solution of〈x1?x2〉is requisite. By subtracting Eq.(2b)from Eq.(2a),we obtain

There are two new variables〈ξ2x1?ξ1x2〉,〈ξ1ξ2(x1?x2)〉appearing in Eq.(9). In order to solve these two variables,multiplying Eq.(6)byξ1ξ2and taking average yield

Multiplying Eq.(2a)byξ1and Eq.(2b)byξ2,then subtracting the two results and averaging yield

Now, we have a set of ordinary differential equations[Eqs. (7), (9), (10), (11)], namely, synchronous equation set(SES),with four variables

Based on the Routh–Hurwitz stability criterion,we derive the stability condition of the SES as follows:

The detailed derivation of the inequality (13) is given in Appendix A.

Next, we will solve the SES when the condition (13) is satisfied, that is, find the steady-state solution of the SES. In this situation, the influence of the initial value will vanish ast →∞. So,without loss of generality,the initial value is set to zero. Taking Laplace transform to be the SES,we have

Applying the inverse Laplace transform to Eq.(16),we obtain

It can be seen from Eq.(17)that the average motions of the two particles are synchronized. So, inequality (13) is named the synchronization condition,that is to say,when inequality(13)is satisfied, the two particles are synchronized. The synchronization between particles makes the mean-field behavior consistent with any single particle behavior in the coupled system.So,in order to analyze the collective behaviors of the system,we only need to consider the first particle’s behaviors below.

3.2. The stability condition and output amplitude gain

In order to obtain the analytical expression of〈x1〉,averaging Eq.(2a),then using Eq.(17),we have

A new variable〈ξ1x1〉appears in Eq. (18). In order to solve this variable, multiplying Eq. (2a) byξ1and averaging, we

Now, we get a closed set of equations [Eqs. (18)–(21)],namely,OAG equation set(OAGES),for four variables:

Based on the Routh–Hurwitz stability criterion, the stability condition for the OAGES can be written as follows:

We call inequality (23) the stability condition. The concrete derivation is given in Appendix B.

Next,we will solve the OAGES with the stability condition(23)satisfied. Without loss of generality, the initial conditions are assumed to be zero. Using the Laplace transform technique,we obtain

Applying the inverse Laplace transform technique to Eq.(26),we obtain the steady state-solution of Eq.(2a)as follows:

whereh1(t)is the inverse Laplace transform ofH1(s).

On the other hand,we can also solve〈x1〉through the linear response theory. Since the input of Eq. (26) is a cosine signal,the system response should also be a cosine signal that differs from the input signal only in amplitude and phase angle. So we have

whereAandφare the amplitude and the phase angle of〈x1〉,respectively. Taking Laplace transform to Eq. (29) and comparing to Eq.(26)withs=j?,we have

According to its definition, the output amplitude gain (OAG)can be written as follows:

From inequality (23) and Eq. (31), we can see that the coupling strength will affect the stability region and the OAG of the system. This is mainly caused by the different random fluctuations. When the noises are identical, the effect of coupling strength disappears.[31]Therefore,compared with some existing works,we think our results can reflect better the influence of the interaction between particles,i.e., coupling strength,on the collective behavior of the system. This is also one of the main results of this article.

3.3. Stability analysis

As can be seen from the inequality (23), the stability of the system is mainly determined by the value ofσ. Whenσis smaller than a certain critical valueσ??, the system is stable. This is because the potential field of system(2)randomly switches betweenV±(xi). Note that (i) for 0＜σ ＜ω, bothV+(xi)andV?(xi)are in steady states,and(ii)for 0＜ω ＜σ,V+(xi)is stable andV?(xi)is unstable. Direct calculations indicate that the right side of the inequality(23)is greater thanω. So the stability of the system is not obvious. According to the physical meaning of (i), whenσ ＜ω, affected by the potential field force,particles will have a tendency to slide to the bottom of the potential field. This tendency leads to the stability of system (2). Meanwhile, it can be seen from inequality (23) thatσ??increases with the increase ofω, and whenω=0,σ??=0. Through the above analysis, it can be obtained that among all the parameters,the stability range(inequality(23))of the system(2)is mainly affected byω.

Sinceλ,ω,εare all positive numbers, direct calculation shows that the stability condition is stronger than the synchronization condition,that is,for any positive numbersλ,ω,ε,

The concrete derivation is given in Appendix C. Inequality(32)shows that if system(2)is stable then the two particles must be synchronized. Due to the inclusion relationship represented by inequality(32),below we mainly discuss the two particles’behaviors when the stability condition(23)is satisfied.

4. Numerical simulation method

The numerical simulation method used in this article is stochastic Taylor expansion. For sufficiently small time step?t, the discrete time representation of the system (2) can be written as follows:

and ?Xiis simulated by the method described in Ref.[55].

In this study, we focus on the average displacement〈xi(t)〉, which can be approximated by Monte Carlo method as follows:

wherexki(t)denotes displacement of thei-th particle in thek-th realization andKis the number of the realizations.

To verify the effectiveness of this algorithm, we need to compare the simulation result with theoretical value. According to Eq. (29),〈xi(t)〉is a cosine signal with the same frequency as the input signalA0cos?t. The theoretical value ofGis given in Eq. (31). Meanwhile, we can obtain the simulation output amplitudeAby performing Fourier transform on the numerical result obtained by Eq.(35),which corresponds to the first panel of Fig. 1. The first panel of Fig. 1 is a concrete realization of〈xi(t)〉. As shown in the bottom panel of Fig.1, the numerical result=A/A0=0.4388 is consistent with the theoretical resultG=0.43906. The relative error is only?5.9217×10?4,which is mainly affected by the number(K)of Monte Carlo simulations and the time step ?t.

Fig.1.The average realization and the corresponding frequency domain representation with ε=1,A0=1,σ=2,ω=2,λ =0.9,? =π/2,?t=10?3,N=104. Different colors represent the two different particles.

Fig.2. The mean absolute percentage error with the system parameters being the same as Fig. 1, (a) as a function of ?t with K =104 and (b) as a function of K with ?t=10?3.

To better reflect the effect of simulation parameters on error,the following mean absolute percentage error is considered:

Under the same system parameters and simulation parameters,the above simulation process is repeatedNtimes,and thenNsimulation valuesj, (j=1,2,...,N)are obtained.Gis the theoretical value given by Eq.(31). Below,we takeN=10 for simulation analysis.

Figure 2 shows the trend ofGMAPEwith the time step ?tand the simulation numberK. From Fig.2(a),we can see thatGMAPEdecreases with the decreasing of ?t. And the simulation results will not depend on ?twhen ?t ≤10?2.6. From Fig. 2(b), we can see thatGMAPEdecreases with the increasing ofK. WhenK=104,GMAPE=5.147×10?4,which has meet the accuracy requirements. So,in all subsequent simulations,the time step ?tand the simulation numberKare set to be 10?3and 104,respectively.

5. Particle behavior analysis

5.1. The effect of stable region

In Section 3, the theoretical solution is derived on the premise that the stability condition (23) is satisfied. So we need to analyze the influence of the stability condition on the particle’s motion behavior. To this end, figure 3 provides the stable region in theσ–εplane. From Fig. 3, we can see that theσ–εplane is divided into two regions. The black region is the stable region which means the stability condition (23)holds,while the white area represents the unstable region.

Fig.3. The stable region in σ–ε plane with ω =2,λ =0.9.

To intuitively show the effect of the stability condition on system (2), we choose one pointA=(2,1) in the stable region and one pointB=(3.5,1)in the unstable region. Figure 4 illustrates the average displacements of the two particles with different noise intensitiesσ. The initial positions of the two particles are random variables obeying normal distribution with mean 0 and standard deviation 0.5.

A stable system is a dynamic system with a bounded response to a bounded input.[42]From the simulations in Fig. 4(a), we observe that the average particle displacements are bound and significantly smaller than that in Fig. 4(b). In Fig. 4(b), the system response is unbounded, which means the system is unstable. This is consistent with the conclusion given in Fig.3.

Fig. 4. Average displacements of the two particles with ω =2, λ =0.9,A0=1,ε =1,? =π/2 for(a)σ =2 and(b)σ =3.5.

5.2. Synchronization

From the inequality(32),it can be seen that when the particles meet the stability condition(23),the average motions of the two particles tend to be synchronized after a period of time.In order to analyze the effect of coupling on particle synchronization, we select two points (A= (2,1),C= (2,4)) with different coupling strengths in the stable region of Fig.3. The corresponding numerical results are given in Fig.5.

Figure 5 shows the average displacements of the two particles and the corresponding variances. In the first panel of Figs.5(a)and 5(b),the average displacement values of the two particles are depicted. It can be seen that, for differentε, the two particles will be synchronized after a period of time. And the particles in Fig.5(a)enter the synchronization state faster than that in Fig. 5(b). This indicates that increasing the coupling strength will reduce the synchronization time. To further analyze the synchronization time,we define the variance ?(t)of the average displacements as follows:

where

The second panel of Figs.5(a)and 5(b)depicts ?(t)as a function oft. We can intuitively see that ?(t) falls faster in Fig. 5(b) than in Fig. 5(a). Whent=2,?=1.1×10?5in Fig.5(a),while?=1.571×10?7in Fig.5(b). To further clarify the difference,we use log-log plot to redraw the two figures in the third panel. In the two new figures,?(t)monotonically decreasing at the beginning,then it began to randomly fluctuate after falling below a certain value. This phenomenon can be explained as follows: The potential energy of thei-th particle is(ω+ξi)x2i/2, and the potential field force it receives is(ω+ξi)xi.Since the noiseξi(i=1,2)are independent of each other and the number(K)of realizations is limited,the difference between the two random potential field forces cannot be completely eliminated,thereby hindering the further reduce of the?. In fact,we have verified through simulation that whenξ1=ξ2,even ifK=1,the displacements of the two particles can be completely synchronized(?＜10?30).In order to illustrate the synchronization time in quantity, we set a threshold??=10?4. When?＜??, we assume that the two particles are synchronized. We define

as the synchronization time. Then the synchronization times in Figs.5(a)and 5(b)are 1.276 and 0.26,respectively.

To further illuminate the effect of coupling on particle synchronization,we depicts the variance?as a function of the timetand the coupling strengthεin Fig.6(a)and depicts the corresponding synchronization timet0as the function of the coupling strengthεin Fig.6(b). From this we can see that,although the particles will synchronize eventually,the coupling strengthεwill affect the synchronization timet0. The synchronization time decreases with the increasing ofε. But the rate of decline gets slower and slower.

Fig.5. The average displacements of two particles and the corresponding variances with ω=2,λ =0.9,σ =2,A0=1,? =π/2 for(a)ε=1 and(b)ε =4.

Fig. 6. Dependence of synchronization time on the coupling strength with ω =2,λ =0.9,A0=1,σ =2,? =π/2. (a)The variance ?as a function of t and ε;(b)the synchronization time t0 as a function of ε.

5.3. Stochastic resonance

SR is the result of the cooperation between the randomness of noise and the orderliness of driving forces under nonlinear action. In our model, the randomness is provided by noise,and the orderliness is yielded by three parts:cosine driving force,coupling force and potential field force generated by frequencyω,i.e.,ωxi. Since the effect of cosine driving force has been extensively discussed in the previous literature, this article mainly analyzes the cooperation phenomenon of noise,coupling force. The mechanism of their cooperation to produce SR is as follows: The noisesξi, (i=1,2) provide randomness to the system through the random fluctuation of the potential field. The system potential field switches randomly betweenV+(xi)andV?(xi).V+(xi)makes thei-th particle have a tendency to slide to the bottom of the potential field, andV?(xi)will slow down this trend and even make thei-th particle slide to infinity(σ ＞ω). Therefore,frequencyωreduces the randomness caused by noise by weakening the amplification effect ofV?(xi). The coupling forces±ε(x1?x2) make the two particles attract each other and move synchronously to provide orderliness for the system. The greater the coupling strengthε, the stronger the orderliness. Under the nonlinear action mode of“multiplicativeness”,orderliness and randomness compete with each other. When the two reach a certain optimal match,the SR occurs,which maximizes the output of the system. Next,we analyze the effect of coupling and noise on the resonance behavior of OAG obtained in Eq.(31).

Firstly, we analyze the influence of coupling strengthεon OAG under different noise amplitudeσ. Figure 7(a)plots the curves of OAG as a function ofεwhenσtakes different values. It can be seen that whenσ=0,Gtakes a constant value, which means that the coupling has no effect on OAG.This is because, whenσ=0, the noise disappears, and system (2) becomes a deterministic system. The statistical synchronization(〈x1〉=〈x2〉)between two particles becomes traditional synchronization (x1=x2). Therefore, the coupling forcesε(x1?x2) andε(x2?x1) disappear, and the effect of the coupling disappears. In fact, whenσ=0, direct calculation can getH1(j?)=1/(ω+j?). So,G=|H1(j?)| has no relationship with the coupling coefficientεin this situation. Whenσis greater than 0,the curvesG(ε)show a singlepeak resonance. This indicates that there is an optimal coupling strength which maximizes the OAG. Asεcontinues to increase,Gwill converge to a limit value. This is because,in strong coupling region, the two particles are bounded as a whole,so that the effect of coupling gradually disappears. Asσincreases, the peak values first increase and then decrease.This indicates that the resonance intensity can be controlled by adjusting the noise amplitude,and there is an optimal noise amplitude to maximize the resonance intensity.

Fig. 7. The OAG versus noise amplitude σ and coupling strength ε with? =π/2, ω =2, λ =0.9, A0 =1. (a)Parameter-induced stochastic resonance (PSR) versus coupling strength ε; (b) SR versus noise amplitude σ,different colored curves represent different coupling strength ε from 0 to 20.

Secondly, we analyze the effect of the noise amplitudeσon OAG under different coupling strengthε. Figure 7(b)plots the curves ofGas a function ofσwhenεtakes different values. It can be seen that single-peak resonance occurs. Asεincreases, the peak values increase first and then decrease.This indicates that there is an optimal coupling strength which maximizes the resonance intensity. In addition, whenεincreases,the curves gradually coincides,that is to say,whenεis large enough,continuing to increaseεwill have no changes.This is consistent with the conclusion given in Fig.7(a).

Finally, we have to point out that the reasons for SR in Figs.7(a)and 7(b)are slightly different. In Fig.7(a),two independent noises introduce different randomness to the trajectories of the two particles. This difference keeps the effect of coupling strength from disappearing over time. And from Eq. (31) we can also see that although the difference in the trajectory of the two particles can be eliminated by averaging,theεhas always had an impact on the particles trajectories.This indicates that the SRversus εis mainly caused by the difference between the two noises. When the two noises are identical,the SR will disappear.[31]However,we can see from Fig. 7(b) that even if the coupling between the two particles disappears (ε=0), which is equivalent to the two noises being identical,the system still exhibits SR.This shows that no matter whether the two noises are identical, the randomness introduced by the two noises will cause the system to produce SRversus σ.

Finally, some numerical simulations are given to verify the above theoretical results. Figure 8 is a comparison diagram of the theoretical results and the simulation results. It can be seen that the theoretical results are basically consistent with the simulations. Due to the greater noise amplitude will bring more randomness,the simulation error increases slightly with the growth ofσin Fig. 8(b). All the numerical simulations provide a good insight into the theoretical results.

Fig.8. Comparison of theoretical and simulation results,with(a)σ =2 and the other system parameters being the same as Fig.7(a); (b)ε =1 and the other system parameters being the same as Fig.7(b).

6. Discussion and conclusion

In this paper, the two coupled harmonic oscillators with dichotomous fluctuating frequency were investigated. First,the synchronization between the two particles and the synchronization condition were derived. Then the stability condition and the OAG of the coupled system were obtained. Based on these theoretical results,we analyzed the collective behaviors of the system,including stability,synchronization,and SR.

The synchronization means that, the mean field behavior is consistent with the behavior of any single particle in this coupled system. Therefore, to study the mean-field behavior of the system, we can only study the behavior of the first particle. Moreover, synchronization allows us to obtain exact steady-state solutions of the system, which provides a theoretical basis for the follow-up research on the collective behaviors.

Comparative analysis shows that, the stability condition is stronger than the synchronization condition, that is to say,when the stability condition is satisfied,the system is simultaneously synchronous and stable. In this situation, the system output is bounded,and through simulation analysis,it is found that the coupling strength has an important effect on the synchronization behavior. The bigger the coupling strength is,the shorter the synchronization time is. This indicates that increasing the coupling strength will reduce the synchronization time.

In the weak coupling region,there is an optimal coupling strength that maximizes the OAG, thus the coupling-induced SR behavior occurs. In the strong coupling region, the two particles are bounded as a whole,so that the coupling’s effect gradually disappears. Coupling will also affect the system’s traditional SR behavior. As the coupling strength increases,the resonance intensity first increases and then decreases.This indicates that there is an optimal coupling strength that maximizes the traditional SR behavior.

The model of this paper is just a basis case with linear coupling form and minimum system size. The collective behaviors of more complex coupled systems remain problems that require further research. We will extend the results of this paper to more complex situations in our future research.

Appendix A: Derivation of the synchronization condition(13)

In this section,we derive the synchronization condition(13). Using Eq.(12),we rewrite the SES as follow:

From Eq.(A1), we get the characteristic equation of the SES as follows:

with

Based on the Routh–Hurwitz stability criterion, we can obtain the stability conditions of the SES in the following inequations:

Then,by simplifying Eq.(A5),the stability conditions of the SES as the function ofσcan be written as follows:

where

Sinceε,σ,λ,ωare all positive numbers,through direct calculation,we have

In conclusion,the stability condition of the SES can be written as follows:

which is the synchronization condition(13).

Appendix B:The derivation of stability condition(23)

In this section,the stability condition(23)is derived. Using Eq.(22),we rewrite the OAGES as follows:

Based on Eq.(B1),the characteristic equation is

with

Using the Routh–Hurwitz stability criterion,the stability condition for the OAGES is as follows:

By using Eq.(B4)to simplify Eq.(B5),the stability conditions for the OAGES as the function ofσcan be written as follows:

where

Sinceε,σ,λ,ωare all positive numbers,through direct calculation,we have

Let

In conclusion,the stability conditions of OAGES can be written as follows:

Appendix C: Comparison of stability condition and synchronization condition

In this section, we will derive the following conclusion:the stability condition(23)is stronger than the synchronization condition(13).

Proof We only need to proveσ??＜σ?(Eqs. (A9) and(B10)). We abbreviateσ?,σ??(Eqs.(A8)and(B9))as

where

By direct calculation,we have

Thus,the conclusion is proved.