Synchronization is a unique form of emergent behavior frequently observed in nature. Since Christiaan Huygens first noted this phenomenon in the 17th century, various types of synchronization among weakly coupled oscillators, including firing neurons and flashing fireflies, have been documented by numerous physicists and biologists [3, 12, 24]. Systematic studies using mathematical models were later initiated by Arthur Winfree [30] and Yoshiki Kuramoto [23]. Among them, the phase synchronization model proposed by Kuramoto is considered as one of the prototype models for synchronization due to its simple mathematical form, the emergence of various synchronization phenomena, and its applicability to real-world cases such as Josephson junction arrays. The model proposed in [23], which is now called the Kuramoto model, is given as follows. Let θi be the phase of the i-th oscillator in R. The dynamics of the N phase configuration Θ:=(θ1,…,θN) is governed by the system of ordinary differential equations:

(1) θ˙i=νi+κN∑j=1Nsin(θj-θi), t>0, θ0i=θini, i∈[N]:={1,2,…,N},

where the positive constant κ denotes the universal scale parameter of coupling strengths between each pair of particles, and νi∈R represents the natural frequency of the i-th oscillator in the absence of interaction, respectively.

As mentioned above, the Kuramoto model and its variations have been extensively studied in the literature as a prototype model for synchronization, particularly in relation to phase-transition behavior. For example, [1, 26] offer several numerical examples of the Kuramoto model from a physicist’s perspective. In [6, 9, 7, 8, 16, 22], the authors also studied variations of the Kuramoto model with additional structures, such as inertia or frustration. In [2, 18, 19], authors explored asymptotic phase and frequency synchronization for various initial conditions. Furthermore, [5, 10, 29, 28, 27] estimated the critical coupling strength needed to ensure the existence of phase-locked states under several conditions, including all-to-all and bipartite graph network topologies. Additionally, [13] provided sufficient conditions for the Kuramoto model with static network topology to exhibit asymptotic phase locking.

However, real-world systems exhibit inherent randomness and uncertainty that cannot be captured by deterministic models. Therefore, there have been many attempts to consider stochastically perturbed Kuramoto models to incorporate this variability, providing a more realistic representation. As far as we know, only stochastically perturbed Kuramoto models of the form

(2) dθti=νi+κN∑j=1Nsin(θtj-θti)dt+∑j=1Nηij(t,Θt)dBtj, t>0, θ0i=θini, i∈[N],

have been considered so far, where the functions B1,…,BN are independent and identically distributed (i.i.d.) standard one-dimensional Brownian motions, and ηij are bounded and Lipschitz continuous. For example, [11, 15] considered the additive white noise (ηij(t,Θt)=constant) to analyze its effect on the synchronization behavior of the Kuramoto model, and [17] used a multiplicative noise

(3)ηij(t,Θt)=2σN∑k=1Nsin(θtk-θti)δ1j, i,j∈[N],

to make the equilibria of (3) asymptotically stable. In [20, 21] sufficient framework was provided to obtain

Plimt→∞|θti-θtj|=0​(mod​π), i,j∈[N]=1,

when the phase configuration Θ satisfies

(4) dθti=κN∑j=1Nsin(θtj-θti)dt+2σN∑j=1Nsin(θtj-θti)dBti, t>0, θ0i=θini, i∈[N],

which is the model that was originally intended to be analyzed in [17] instead of (4). The main difference between additive and multiplicative noise models lies in their behavior around equilibrium points. In additive noise models, the phase configuration generally stays near an equilibrium point for a long time, but eventually moves to other equilibrium point. In contrast, the equilibrium points that were stable in the deterministic setting are also asymptotically stable in multiplicative noise models. This difference arises because, in multiplicative noise models, the diffusion coefficient near the equilibrium point is close to zero, whereas it is not in additive noise models. Among all these results, there are no conditions that are necessary and sufficient for (1) or (4) to produce specific types of synchronization phenomena. Therefore, it is still meaningful to properly discretize the Kuramoto model, as numerical experiments with nonlinear differential equations often reveal facts about their solutions that cannot be proven due to technical limitations.

In this paper, we are interested in time discretization that preserves the asymptotic behavior of (4). More precisely, our interest lies in whether the already known asymptotic behavior of (4) can be accurately reproduced using the Euler-Maruyama method, one of the simplest methods for discretizing stochastic differential equations. Due to the similarity between drift and diffusion coefficients of (4), the approximated Euler-Maruyama solution for (4) can be obtained iteratively by using the recurrence relations

(5) θn+1i=θni+κh+2σhZniN∑j=1Nsin(θnj-θni), i∈[N], n∈N∪{0}, θ0i=θini, Zni∼i.i.d.N(0,1), i∈[N], n∈N∪{0}

where h:=Δt>0 is a fixed time step, and each Θn=(θn1,…,θnN) approximates the solution of (4) at time nh. In addition, each Zni is a standard normal random variable where Zni and Zmj are independent if and only if n≠m or i≠j. For the model (4), we will consider an order parameter

R(Θ)=1N∑j=1Neiθj

as an indicator of the phase synchronization, as it takes value in [0,1] and becomes 1 only if θ1=⋯=θN modulo 2π. Then, if there is a sufficient condition on κ,σ,h to satisfy

0≤E[R(Θ1)2]≤⋯≤E[R(Θn)2]≤⋯≤1,

we can say that (4) is synchronized in the sense that E[R2] increases monotonically. We will find such a condition by determining a function F=F(Θ) satisfying

E[R(Θn+1)2|Θn]-R(Θn)2≥F(Θn)≥0,

which make the series ∑nE[F(Θn)] converges to a finite number. Then, by using the usual Borel-Cantelli argument, we will obtain

limn→∞F(Θn)=0 almost surely,

which indicates a specific synchronization phenomenon that the phase configuration will asymptotically exhibit a special shape. Although the use of the order parameter was also an important idea in the results of [17, 21, 20], the difference is that the Itô calculus cannot be used for the discrete-time model (4). Because of the nonlinear structure of R2 and the presence of unbounded noises Zni, it is impossible to approximate R(Θn+1)2-R(Θn)2 with the desired accuracy with probability 1 by Taylor approximation. Therefore, the main novelty of this paper is to show that even though we cannot control the individual sample paths of R(Θn)2, we can compute its expectation accurately and achieve the desired result.

The rest of this paper is organized as follows. In Section 2, we review several concepts and properties on (1) and (4), such as the order parameter and phase synchronization. In Section 3, we present a sufficient framework leading to the monotonic increasing property of E[R2] for (5). In Section 4, we verify our analysis through numerical experiments and discuss the optimality of the results. Finally, Section 5 is devoted to a brief summary of our results.

In this section, we present several well-known results on the deterministic Kuramoto model (1) and the stochastic Kuramoto model (4). First, we introduce the definition and some properties of order parameters, which are frequently used throughout our discussion.

Definition 2.1. Let Θ=(θ1,⋯,θN) be a phase configuration in RN. Then, the order parameters (R(Θ),ϕ(Θ))∈R×R/2πZ are defined as

(1)R(Θ):=1N∑ j=1 Neiθj, eiϕ(Θ):=1NR(Θ)∑ j=1Neiθj(if R(Θ)≠0). 

Since the above order parameters are well defined by the phase configuration alone, they can be thought of in a general context, not just for the Kuramoto model. To be more precise, in the case of φ, it can be said that the values of cosϕ and sinϕ are determined rather than just φ. Therefore, the following lemma provides a useful way to incorporate information about φ.

Lemma 2.2. [25] Let Θ=(θ1,⋯,θN) be a phase configuration in RN. If R(Θ)>0, the order parameters (R(Θ),ϕ(Θ)) satisfy

(2)R(Θ)cos(ϕ(Θ)-θ)=1N∑k=1Ncos(θk-θ),​R(Θ)sin(ϕ(Θ)-θ)=1N∑k=1Nsin(θk-θ),

for every θ∈R. In particular, we have

(3) R(Θ)2=1N∑j=1Ncos(θj-θ)2+1N∑j=1Nsin(θj-θ)2 ∀​θ∈R, R(Θ)2=1N2∑i,j=1Ncos(θi-θj), R(Θ)=1N∑j=1Ncos(θj-ϕ(Θ)), 0=1N∑j=1Nsin(θj-ϕ(Θ)).

Note that if R(Θ)=0, we can use ∑k=1Neiθk=0 to obtain

1N∑k=1Ncos(θk-θ)=Re1N∑k=1Nei(θk-θ)=Re1N∑k=1Neiθk·e-iθ=0,1N∑k=1Nsin(θk-θ)=Im1N∑k=1Nei(θk-θ)=Im1N∑k=1Neiθk·e-iθ=0.

Therefore, although ϕ(Θ) is not well-defined in this case, (2) and (3) holds by replacing the undetermined number ϕ(Θ) to arbitrary constant.

Now, given the dynamics of the phase configuration, the derivative of R(Θ) can be expressed in terms of Θ by differentiating (3) and combining the derivative of Θ with (2). If we do this for the Kuramoto model, we obtain the following equation.

Proposition 2.3. Let Θ(t)=(θ1(t),…,θN(t)) be a solution to the Kuramoto model (1). Then, R=R(Θ) satisfy

(4)R˙=-1N∑j=1Nνjsin(θj-ϕ)+κRN∑j=1Nsin2(θj-ϕ).

This is, in fact, a crucial estimate for the convergence of the Kuramoto model when all natural frequencies are identical. If νi≡ν, we can use Lemma 2.2 and Proposition 2.3 to simplify (4) as

dR2dt=2κR2N∑j=1Nsin2(θj-ϕ),

so that R is monotonically increasing. Therefore, if the initial configuration Θin satisfies R(Θin)>0, the phase ϕ(Θ)∈R/2πZ is well-defined. Then, since R is bounded above by 1, each sin(θj-ϕ) approaches zero as t tends to infinity, and we have the following result.

Proposition 2.4. [4] Let Θ(t)=(θ1(t),…,θN(t)) be a solution to (1), where the natural frequencies {νi}i=1N and initial configuration Θin=(θin1,…,θinN) satisfy

νi≡ν, R(Θin)>0.

Then, the limit

θ∞i:=limt→∞(θi(t)-νt)

exists for all i∈[N] and

sin(θ∞i-θ∞j)=0, ∀​i,j∈[N].

When we consider the stochastic Kuramoto model (4), we assume there exists a probability triple (Ω,F,P), and denote (Ft)t≥0 the natural filtration generated by {Bt1,…,BtN}. That is, Ft is the sub σ-algebra of F generated by {Bs1,…,BsN:s≤t}. Similar to Proposition 2.3, we can apply Itô's lemma to the stochastic Kuramoto model (4) as follows.

Proposition 2.5. [20, 21] Let Θt=(θt1,…,θtN) be a solution process to the stochastic Kuramoto model (4). Then, Rt=R(Θt) satisfy

(5)Ru2-Rs2=∫su[2κRt2Pt+2σRt2PtN-2σRt31N∑j=1Nsin2(ϕt-θtj)cos(ϕt-θtj)]dt  +22σN∑j=1N∫suRt2sin2(ϕt-θtj)dBtj,

where the process Pt is given by

Pt:=1N∑i=1Nsin2(θti-ϕt).

Therefore, by taking conditional expectation to both hand sides of (5), one can show that R(Θ)2 is a submartingale relative to ({Ft},P), provided that κ>1-1Nσ:

E[Ru2|Fs]=Rs2+∫suE[Rt2(2κPt+2σPtN-2σRt1N∑j=1Nsin2(ϕt-θtj)cos(ϕt-θtj))]dt ≥Rs2+∫suE[Rt2(2κPt+2σPtN-2σPt)]dt =Rs2+∫suE[Rt2Pt(2κ-2σ(1-1N))]dt ≥Rs2, 0≤s≤u<∞.

In Section 3, we will provide a direct proof of the monotone increasing property of the expectation of R(Θn)2, without reproducing Proposition 2.5 in the discrete model (5). This result will lead to the conclusion that Plimn→∞R(Θn)2P(Θn)=0=1.

In this section, we provide a synchronization estimate of (5). By using equation (2), we can rewrite the right hand side of (5) to

(1) θn+1i=θni-(κh+2σhZni)Rnsin(θni-ϕn), i∈[N], n∈N∪{0}, θ0i=θini, Zni∼i.i.d.N(0,1), i∈[N], n∈N∪{0}, Rn:=R(Θn), ϕn:=ϕ(Θn),

where we use the following relation instead of (1) when Rn=0:

Rn=0⟹Rnsin(θni-ϕn)=0 ∀​i∈[N].

Then, the idea is calculating the conditional expectation of R(Θn+1)2 for a fixed Θn. Once Θn is given, (1) shows that each θn+1i is a normal random variable, or just equal to θni if Rnsin(θni-ϕn). If we apply (3)₂ in Lemma 2.2, Rn+12 can be considered as an average of resultant probability distributions obtained by taking the cosine of normal random variables. Therefore, we can use the following lemma to get E[Rn+12|Θn] explicity.

Lemma 3.1. For every λ,μ∈R and a standard normal random variable Z, we have

E[cos(λZ+μ)]=e-λ22cosμ.

Proof. We use the fact that the characteristic function of the Gassian distribution is given by

E[eiλZ]=e-λ22

(see [14], p.210). Therefore, we have

E[cos(λZ+μ)]=E[Re(eiλZ+iμ)] =ReeiμE[eiλZ] =Ree-λ22(cosμ+isinμ) =e-λ22cosμ.

Now, using Lemma 3.1, (3)₂ and introducing a simple notation, we can compute the conditional expectation as follows.

Lemma 3.2. Let (Θn)n≥0=(θn1,…,θnN)n≥0 be a solution process to (1). Then, we have

E[Rn+12|Θn] =1N2N+∑j≠ke-σhRn2sin2θ^nj+sin2θ^nkcosθ^nj-θ^nk-κhRn(sinθ^nj-sinθ^nk),

where the augmented phase θ^ni is defined as

(2)θ^ni:=θni-ϕn, i∈[N], n∈N∪{0}.

Proof. First, let's fix j,k∈[N] and apply Lemma 3.1 to compute E[cos(θn+1j-θn+1k)|Θn]. From (1), we can represent θn+1j-θn+1k with respect to Θn and (Zn1,…,ZnN) as follows:

(3)θn+1j-θn+1k=θnj-(κh+2σhZnj)Rnsin(θnj-ϕn) -θnk-(κh+2σhZnk)Rnsin(θnk-ϕn) =θnj-θnk-κhRn(sin(θnj-ϕn)-sin(θnk-ϕn)) +2σhRn-sin(θnj-ϕn)Znj+sin(θnk-ϕn)Znk.

Since θnj-θnk=(θnj-ϕn)-(θnj-ϕk), we may use the augmented phases (2) to rewrite (3) as

(4)θn+1j-θn+1k=θ^nj-θ^nk-κhRn(sinθ^nj-sinθ^nk) +2σhRn-sinθ^njZnj+sinθ^nkZnk.

If j≠k, the latter term in the right hand side of (4) is also a normal random variable with variance

(λnjk)2:=2σhRn2sin2θ^nj+sin2θ^nk,

as Znj and Znk are independent. Therefore, we denote

μnjk:=θ^nj-θ^nk-κhRn(sinθ^nj-sinθ^nk)

and apply Lemma 3.1 to obtain

E[cos(θn+1j-θn+1k)|Θn]=e-(λnjk)22cosμnjk =e-σhRn2sin2θ^nj+sin2θ^nkcosθ^nj-θ^nk-κhRn(sinθ^nj-sinθ^nk),

for every j≠k. Finally, using (3)₂ in Lemma 2.2 yields

E[Rn+12|Θn] =E[1N2∑j,kcos(θn+1j-θn+1k)|Θn] =1N2(N+∑j≠kE[cos(θn+1j-θn+1k)|Θn]) =1N2N+∑j≠ke-σhRn2sin2θ^nj+sin2θ^nkcosθ^nj-θ^nk-κhRn(sinθ^nj-sinθ^nk).

Next, we need to find the condition that makes E[Rn+12|Θn]≥Rn2. Since κ,σ, and N are constants that do not change and are variables used in the continuous-time model, we find a condition for the newly added variable h in the discrete-time model (1). In addition, since Θn varies depending on n, the desired condition should not depend on Θn. If we consider the explicit formula for E[Rn+12|Θn] in Lemma 3.2 as a function fn=fn(h) that depends only on h, we have

fn(0)=1N2N+∑j≠kcos(θ^nj-θ^nk) =1N2∑1≤j,k≤Ncos(θnj-θnk) =Rn2.

Then, we want to find the condition of h such that fn(h)≥fn(0) for every n, and for this we can use Taylor's theorem

fn(h)=fn(0)+hfn′(0)+h22fn′′(h*), ∃​h*∈[0,h].

Therefore, in the following lemma, we find a positive constant h0 satisfying

• h0=h0(κ,σ,N),

• independent of n and Θn,

• h0≤minh*∈[0,h0]2fn′(0)|fn′′(h*)| for all n,

so that the relation E[Rn+12|Θn]≥Rn2 holds for every h≤h0.

Lemma 3.3. Let (Θn)n≥0=(θn1,…,θnN)n≥0 be a solution process to (1), where κ,σ,N and h satisfy

κ>1-1Nσ, 0<h≤2κ-21-1Nσκ2+21-1Nσ2.

Then, we have the following inequality: for every n∈N∪{0},

(5)E[Rn+12|Θn]≥Rn2+h2κ-21-1Nσ-hκ2+21-1Nσ2Rn2Pn ≥Rn2.

Proof. From now, we use the following notation for simplicity:

cosθ^ni=c^ni, sinθ^ni=s^ni, i∈[N], n∈N∪{0}.

Then, we can write fn(h) as

fn(h)=1N2N+∑i≠je-σhRn2(s^ni)2+(s^nj)2cosθ^ni-θ^nj-κhRn(s^ni-s^nj).

By differentiating fn with respect to h, one can obtain

fn′(h)=1N2∑i≠je-σhRn2((s^ni)2+(s^nj)2)   [-σRn2((s^ni)2+(s^nj)2)cos(θ^ni-θ^nj-κhRn(s^ni-s^nj)) +κRn(s^ni-s^nj)sin(θ^ni-θ^nj-κhRn(s^ni-s^nj))] =1N2∑i≠je-σhRn2((s^ni)2+(s^nj)2)gnij(h),

where gnij(h) is defined as

(6)gnij(h):=-σRn2((s^ni)2+(s^nj)2)cos(θ^ni-θ^nj-κhRn(s^ni-s^nj)) +κRn(s^ni-s^nj)sin(θ^ni-θ^nj-κhRn(s^ni-s^nj)).

Now, we are ready to estimate f′n(0) and fn′′(h).

(1) For fn′(0), the direct computation yields

fn′(0)=1N2∑i≠jgnij(0) =1N2∑i≠j[-σRn2((s^ni)2+(s^nj)2)cos(θ^ni-θ^nj)+κRn(s^ni-s^nj)sin(θ^ni-θ^nj)] =1N2∑i≠j[-σRn2((s^ni)2+(s^nj)2)(c^nic^nj+s^nis^nj)+κRn(s^ni-s^nj)(s^nic^nj-c^nis^nj)] =1N2∑1≤i,j≤N[-σRn2((s^ni)2+(s^nj)2)(c^nic^nj+s^nis^nj)+κRn(s^ni-s^nj)(s^nic^nj-c^nis^nj)] +2σRn2N2∑i=1N(s^ni)2.

Then, we use (3)₃ of Lemma 2.2, i.e.,

1N∑i=1Nc^ni=Rn, 1N∑i=1Ns^ni=0,

to obtain

fn′(0)=1N2∑1≤i,j≤N[-σRn2((s^ni)2+(s^nj)2)(c^nic^nj)+κRn((s^ni)2c^nj+(s^nj)2c^ni)] +2σRn2N2∑i=1N(s^ni)2 =2κRn2N∑i=1N(s^ni)2-2σRn3N∑i=1N(s^ni)2c^ni+2σRn2N2∑i=1N(s^ni)2.

This exactly coincides with the drift part of (4), confirming that all calculations up to this point are correct. Therefore, if we use the notation

Pn:=1N∑i=1N(s^ni)2,

we have

(7)fn′(0)≥2Rn2Pnκ-σRn+1Nσ.

(2) For fn′′(h), we first represent fn′′(h) in terms of gnij and (gnij)′ as follows:

(8)fn′′(h)=1N2∑i≠je-σhRn2((s^ni)2+(s^nj)2)[gnij(h)·(-σRn2((s^ni)2+(s^nj)2))+(gnij)′(h)],

where each (gnij)′(h) is given by

(9)(gnij)′(h)=-κσRn3((s^ni)2+(s^nj)2)(s^ni-s^nj)sin(θ^ni-θ^nj-κhRn(s^ni-s^nj)) -κ2Rn2(s^ni-s^nj)2cos(θ^ni-θ^nj-κhRn(s^ni-s^nj)).

To compute an upper bound of (8), we again introduce the following notation:

Cnij(h):=cos(θ^ni-θ^nj-κhRn(s^ni-s^nj)), Snij(h):=sin(θ^ni-θ^nj-κhRn(s^ni-s^nj)).

Then, we combine (4) and (9) to obtain

gnij(h)·(-σRn2((s^ni)2+(s^nj)2))+(gnij)′(h) =(-σRn2((s^ni)2+(s^nj)2)Cnij(h)+κRn(s^ni-s^nj)Snij(h))·(-σRn2((s^ni)2+(s^nj)2)) -κσRn3((s^ni)2+(s^nj)2)(s^ni-s^nj)Snij(h)-κ2Rn2(s^ni-s^nj)2Cnij(h) =σ2Rn4((s^ni)2+(s^nj)2)2-κ2Rn2(s^ni-s^nj)2Cnij(h) -[2κσRn3((s^ni)2+(s^nj)2)(s^ni-s^nj)]Snij(h).

Since |ACnij(h)+BSnij(h)|≤A2+B2 for every A,B∈R, we have

|gnij(h)·(-σRn2((s^ni)2+(s^nj)2))+(gnij)′(h)| ≤σ2Rn4((s^ni)2+(s^nj)2)2-κ2Rn2(s^ni-s^nj)22+[2κσRn3((s^ni)2+(s^nj)2)(s^ni-s^nj)]2 =σ2Rn4((s^ni)2+(s^nj)2)2+κ2Rn2(s^ni-s^nj)2.

Therefore, we apply this inequality to (8) to obtain

(10)|fn′′(h)|≤1N2∑i≠j[gnij(h)·(-σRn2((s^ni)2+(s^nj)2))+(gnij)′(h)] ≤1N2∑i≠j[σ2Rn4((s^ni)2+(s^nj)2)2+κ2Rn2(s^ni-s^nj)2] =1N2∑1≤i,j≤N[σ2Rn4((s^ni)2+(s^nj)2)2+κ2Rn2(s^ni-s^nj)2]-1N2∑i=1N4σ2Rn4(s^ni)4 =σ2Rn41N2∑1≤i,j≤N(s^ni)4+2(s^ni)2(s^nj)2+(s^nj)4-4N2∑i=1N(s^ni)4 +κ2Rn21N2∑1≤i,j≤N(s^ni)2-2s^nis^nj+(s^nj)2 =σ2Rn42-4N1N∑i=1N(s^ni)4+2Pn2+κ2Rn22Pn,

where we used

1N∑i=1Ns^ni=0,

in the last equality.Then, we combine (7) and (10) to obtain

(11)2fn′(0)|fn′′(h*)|≥4Rn2Pnκ-σRn+1Nσσ2Rn42-4N1N∑i=1N(s^ni)4+2Pn2+2κ2Rn2Pn ≥4Rn2Pnκ-σRn+1Nσσ2Rn42-4NPn+2Pn2+2κ2Rn2Pn =2κ-σRn+1Nσκ2+σ2Rn21-2N+Pn ≥2κ-21-1Nσκ2+21-1Nσ2,

for every h*∈R. Finally, we apply (7) and (11) to obtain

E[Rn+12|Θn]=fn(h) =fn(0)+hfn′(0)+h22fn′′(h*) ≥fn(0)+h1-κ2+21-1Nσ22κ-21-1Nσ·hfn′(0) ≥fn(0)+h2κ-21-1Nσ-hκ2+21-1Nσ2Rn2Pn.

for all n∈N∪{0}, which is our desired result.

Remark 3.4. In [25], the author proved that the discrete deterministic Kuramoto model

θn+1i=θni-κhRnsin(θni-ϕn), i∈[N], n∈N∪{0}, θ0i=θini, i∈[N],

has a monotone increasing property for Rn2 when κ>0 and 0<κh≤2. Thus, the result of Lemma 3.3 is consistent with [25] as σ tends to 0.

Consequently, we can show the convergence of the series ∑nE[Rn2Pn]. This not only implies that Rn2Pn will converge to zero almost surely, but also leads to all of the following results.

Theorem 3.5. Let (Θn)n≥0=(θn1,…,θnN)n≥0 be a solution process to (1), where κ,σ,N and h satisfy

κ>1-1Nσ, 0<h<2κ-21-1Nσκ2+21-1Nσ2.

Then, we have

Plimn→∞Rns^ni=0=Plimn→∞Rns^niZnj=0=Plimn→∞|θn+1i-θni|=0=1, ∀​i,j∈[N].

In addition, there exists a random variable R∞:Ω→[0,1] satisfying

Plimn→∞Rn=R∞=1.

Proof. First, we set a constant independent of n:

C(κ,σ,N,h):=h2κ-21-1Nσ-hκ2+21-1Nσ2>0.

Then, we apply the tower property of the conditional expectation to (4) to obtain

E[Rn+12]=EE[Rn+12|Θn] ≥ERn2+C(κ,σ,N,h)Rn2Pn =E[Rn2]+C(κ,σ,N,h)E[Rn2Pn], ∀​n∈N∪{0}.

Since the order parameter Rn is always containted in [0,1], this yields

C(κ,σ,N,h)∑n=0∞E[Rn2Pn]≤∑n=0∞(E[Rn+12]-E[Rn2])≤1-Rin2<∞.

Note that Θn and Zni are independent for every i∈[N] and n∈N∪{0}. Thus, we can also obtain

∑n=0∞E[Rn2Pn]=∑n=0∞E[Rn2Pn(Zni)2]<∞, ∀​i∈[N].

Then, we use the Markov inequality to get

(12) ∑n=1∞PRn2Pn≥ε≤1ε∑n=1∞E[Rn2Pn]<∞, ∀​ε>0, ∑n=1∞PRn2Pn(Zni)2≥ε≤1ε∑n=1∞E[Rn2Pn(Zni)2]<∞, ∀​ε>0.

By using the Borel-Cantelli lemma, (12) implies that

P∩n=1∞∪k=n∞{Rk2Pk≥ε}=0, P∩n=1∞∪k=n∞{Rk2Pk(Zki)2≥ε}=0, ∀​i∈[N], ε>0,

which shows the P-almost sure convergence of Rk2Pk and Rk2Pk(Zki)2 to 0. Since Rn2Pn and Rn2Pn(Znj)2 can be written as

Rn2Pn=1N∑i=1N(Rns^ni)2, Rn2Pn(Znj)2=1N∑i=1N(Rns^niZnj)2,

we have the desired P-almost sure convergence of Rns^ni and Rns^niZnj as n→∞. In addtion, we apply these almost sure convergence results to (1) to obtain

Plimn→∞|θn+1i-θni|=0=Plimn→∞Rns^ni|κh+2σhZni|=0=1.

Finally, since Rn2 is a bounded submartingale (relative to natural filtration generated by Θn), we apply Doob's supermartingale convergence theorem to show the existence of finite random variable R∞:=limn→∞Rn.

Remark 3.6. Suppose the system parameters κ,σ,N,h satisfy the assumption of Theorem 3.5. If R(Θin)>0, the expectation of R∞2 must be positive due to the monotone increasing property of E[Rn2]. Therefore, there are positive measures of ω∈Ω such that

limn→∞Rn(ω)=R∞(ω)>0, limn→∞R(Θn(ω))s^ni(ω)=0, i∈[N].

If all s^ni(ω) converges to zero as n tends to infinity, we have

limn→∞sin(θni(ω)-ϕn(ω))=0, i∈[N],

and this leads to the following conclusion: for every i,j∈[N], we have

limn→∞sin(θni-θnj)=0.

This result corresponds to the synchronization of deterministic Kuramoto model, for instance the result in Proposition 2.4.