Numerical Solution of Finite Kuramoto Model with Time-Dependent Coupling Strength: Addressing Synchronization Events of Nature

: The synchronization of an ensemble of oscillators is a phenomenon present in systems of different ﬁelds, ranging from social to physical and biological systems. This phenomenon is often described mathematically by the Kuramoto model, which assumes oscillators of ﬁxed natural frequencies connected by an equal and uniform coupling strength, with an analytical solution possible only for an inﬁnite number of oscillators. However, most real-life synchronization systems consist of a ﬁnite number of oscillators and are often perturbed by external ﬁelds. This paper accommodates the perturbation using a time-dependent coupling strength K ( t ) in the form of a sinusoidal function and a step function using 32 oscillators that serve as a representative of ﬁnite oscillators. The temporal evolution of order parameter r ( t ) and phases θ j ( t ) , key indicators of synchronization, are compared between the uniform and time-dependent cases. The identical trends observed in the two cases give an important indication that the synchrony persists even under the inﬂuence of external factors, something very plausible in the context of real-life synchronization events. The occasional boosting of coupling strength is also enough to keep the assembly of oscillators in a synchronized state persistently.


Introduction
The phenomenon by which a large population of elements oscillating with different natural frequencies lock to a common frequency is called synchronization.Such phenomena have been observed in many real-world systems, such as rhythmic response patterns of pacemaker cells in the heart [1], a synchronously flashing group of fireflies [2,3], crickets that chirp in unison [4], etc. Arrays of lasers and superconducting Josephson junctions are examples in physics that obey collective synchronization [5,6].
One of the most successful approaches for exploring and analyzing synchronization phenomenon for the population of elements as an ensemble of phase oscillators is given by Kuramoto [7], which treats each oscillator only by their phases [8].The oscillators proposed in his model are assumed to be non-linear, and there is a weak but uniform coupling between every oscillator (all-to-all coupling).
There is a threshold K c for the coupling strength K, which characterizes the interaction term in the Kuramoto mathematical model, above which a phase transition analogous to magnetic phase transition occurs [9].In other words, the system shifts from an incoherent state to a partially synchronized state and ultimately to the fully synchronized one with K >> K c .Kuramoto's original model can be considered a basic model in the sense that it describes oscillators of the fixed natural frequencies w j , connected via 'all-to-all' uniform K [10] and solvable only in the thermodynamic limit (N → ∞).However, the real-world systems that exhibit synchronous behavior consist of a finite number of oscillators, and the coupling is expected to be perturbed by time-dependent external factors.Some of the promising examples include the natural firing rate of neurons in the brain, which correspond to the natural frequencies, and the degree of excitation/inhibition from neurons to stimulate a neuron via a dendrite (modeled as time-dependent coupling strength) [10][11][12].
To address the inevitable influence of external fields, this paper incorporates a timedependent coupling strength of the form K(t) = K 1 (1 + A sin(Ωt)) as suggested in [10,11,13] and also mentioned in [14] with a slightly different scope.However, the term K 1 plays the same role as the uniform coupling strength of Kuramoto's original model, a different approach from [10,11,13].The influence of the external field can be oscillatory, so the choice of the term A sin(Ωt) is justifiable.
The above form of K(t) is tested here for an 32 oscillators system, an even power of 2 oscillators and a representative of a small number of finite oscillators such as synchronously clapping audiences at an orchestra.
This paper also intends to use coupling strength K(t) in the form of a step function.Unlike the trigonometric function, A sin(Ωt), which acts continuously in the system of oscillators, a different (higher) value of coupling strength K acts for only a certain time interval where the order parameter r(t), a key indicator that reflects the synchronous behavior of oscillators (Section 2), goes beyond its pre-defined threshold value.It is assumed that the system of such oscillators' assembly preserves its autonomy except at some intervals of time where r approaches some threshold value.K(t) serving as a stepfunction has been tested in various articles, notably in [15], which aims to implement this approach for two oscillators' assembly.However, the approach implemented here is unique in the sense that K(t) only takes a higher value if the order parameter r(t) goes beyond a certain threshold.Additionally, the step-function model of K(t) tested here is for 32 oscillators, which can serve as a representative for a wide range of finite oscillators exhibiting synchronous behavior in nature.
As mentioned earlier, there is a phase transition from an incoherent state to being partially synchronized in the coupling strength K ≥ K c and finally to the fully synchronized state for K >> K c .The testing of synchronized states at some value of K greatly depends on the critical value K c for the given oscillators' assembly.However, such critical value has only been found analytically in the infinite oscillators' assembly so far.Still, this paper picks the value of K c equal to the critical value associated with infinite oscillators that comes from Equation (8).The synchronization pattern, if observed in the vicinity of such analytical value, will help us to locate the respective K c of finite oscillators more precisely.
The following two methods will be followed to observe the impact of external fields,

•
Step function: K < K c can be used except for the order parameter r(t) ≤ r th , where a higher K > K c is applied to elevate oscillators so that they maintain their synchronized state; • Trigonometric function: K 1 can be set in the vicinity of K c to observe whether the perturbation part A sin(Ωt) can make an impact to elevate the oscillators into the synchronized state.
The synchronization events observed in nature have displayed their synchronous behavior persistently, despite the presence of external fields, so the results from this proposed work are expected to find a close match between time-dependent and uniform parameters.In the case of a step function, its usefulness lies in whether the oscillators maintain the synchrony with the boost of a higher coupling strength K value for r ≤ r th .

Kuramoto Model with the Infinite Number of Oscillators
In a system of N phase oscillators, each oscillating with natural frequencies w j , the time evolution of phase θ j of the j th oscillator is given by following Kuramoto [7] equation: where K ≥ 0, is the homogeneous all-to-all coupling strength.The factor 1/N ensures that the system converges in the thermodynamic limit (as N → ∞).The frequencies w j are distributed according to some probability density g(w), where g(w) is assumed to be unimodal and symmetric about the mean frequency ν i.e., g(ν + ω) = g(ν − ω) [6,16].For example, g(w) can be a Gaussian distribution or a Cauchy-Lorentzian distribution.Nonidentical values of the frequency of each individual, especially in the context of biological oscillators such as fireflies, reflect their inevitable variability found in nature [17].Equation ( 1) can be related to a physically sensible quantity called order parameter r(t), which is simply an average (centroid vector of the phase distribution) of the complex number (Z) that represents the phase of the oscillators with each moving in a unit circle.
Multiplying Equation ( 2) by e −iθ j (t) , the following relation is obtained.
Now, comparing the imaginary part of the above relation and making the substitution in Equation ( 1), we find: In the thermodynamic limit, oscillators are expected to be distributed with a probability density ρ(θ, ω, t).Now, the arithmetic mean in Equation (2) becomes an average over phase and frequency: For K → 0, Equation (5) becomes Therefore, r(t) → 0 as t → ∞ and the oscillators are out of synchrony.However for These integrals are equal to 1 as they are the integral of densities over the entire range.It means r(t) → 1 and there is perfect synchrony among the oscillators [7,16].For the intermediate range of coupling strength, K c < K < ∞ (with K c , critical coupling strength; defined in Equation ( 8)), some of the oscillators are phase-locked ( θj = 0) and others rotate out of synchrony with the locked oscillators [8].
Kuramoto obtained the exact solution of K c after evaluating the contribution of both drifting and locked oscillators (Equation ( 5)) as where g(0) is the frequency distribution of interest, evaluated at ω = 0.
For the Gaussian distribution of frequencies that is used in this research, the approximate solution or r(t) near the bifurcating point K c is given by [8] r ≈ where g (0) is the second derivative of g(0) evaluated at ω = 0.

Synchronization Visualization
The exact solution for the critical coupling strength K c given in Equation ( 8) is achievable as N → ∞.The numerical solution of N = 32 oscillators for this work will pick the value of K c that comes from Equation (8).
Figure 1 compares the behavior of r(t) versus K for four different assemblies of oscillators, with numbers ranging from 32 to 256 near the onset of the synchronization.The natural frequencies w j are drawn from the Gaussian.The number of oscillators used in four cases replicates the number in various synchronized events of nature, such as fireflies displaying synchrony.As N increases, the bifurcation point approaches its counterpart for N → ∞.In the case of 32 oscillators, however, there is a larger discrepancy.Figure 2a displays the finite size effect on the temporal evolution of phase angle θ j (t) and the order parameter r(t) of the 32 oscillators.The use of K c = 0.32 for dealing with 32 oscillators cannot be justified; however, K = 0.5 > K c keeps the finite size effect on the sideline.The phases of oscillators start gathering coherently for K > K c .The fluctuation (≈ N − 1 2 ) of r(t) is significantly higher in the 32 oscillator case compared to Figure 2b, which uses 256 oscillators.
(a) (b) Figure 2. Synchronization pattern: 32 vs. 256 oscillators.The critical value K c is expected to reach closer to 0.32 (exact value in the thermodynamic limit, N → ∞) with the increase in the number of oscillators.Two partially synchronized groups are visible in the 32 oscillator case (a), while there are many of them in the case of 256 oscillators (b).Additionally, the fluctuation in the typical r(t) vs. t plot, which is proportional to N − 1 2 , is less prominent in the case of the 256 oscillator system compared to 32 oscillators.

Time-Dependent Parameters
The inherent time-variability of oscillators can be addressed by incorporating timedependent coupling strength K(t) in Kuramoto's original model, which assumes all-to-all homogeneous coupling K.
In this research work, the following form of K(t) has been used: where K 1 is the uniform coupling strength with the value chosen in the vicinity of K c ).Other parameters, A and Ω, represent the amplitude and frequency of the perturbation term, respectively.Now, we introduce K(t) for the system of N = 32 oscillators as: The natural frequencies w j are chosen randomly from the standard Gaussian density function.The amplitude A and frequency Ω of the perturbation term have also been treated as variables separately to see their impact on the synchronization pattern in the r(t) vs. t graph using Equation (2).

Results
As discussed in Section 2, the main objective of this work is to check for the similarity in r(t) and θ(t) vs. t graphs using the time-dependent and uniform coupling strengths in the vicinity of K c .
Numerical integration of Equation ( 11) was carried out using the Runge-Kutta 4th order method, with a time step of dt = 0.01 throughout 100 (some unit of time).The values of K 1 = K have been chosen in the vicinity of critical coupling strength K c based on Equation (8).The natural frequencies of the oscillators were randomly chosen to be from standard Gaussian distribution with mean µ = 0 and standard deviation σ = 0.2.Initial phases θ j (0) were also randomly chosen in the range of [0,2π].A default setting of the random number generation was used to ensure the same numbers were being used in the successive iteration.Now, the standard Gaussian distribution of natural frequencies given by From the above equation g(0) ≈ 1.995 for µ = 0 and σ = 0.2.Critical coupling strength K c based on Equation ( 8) equals to

Using Periodic Function as Perturbation
The comparisons were made in the r(t) vs. t graphs using the uniform and timevarying coupling strengths.Additionally, the response of K(t) was observed more closely in r(t) vs. t by varying amplitude A and frequency Ω of the perturbation term separately.
The temporal evolution of r(t) in Figure 3a for the K = 0.5 & K(t) = K 1 + A sin(Ωt) with K 1 = 0.5 are similar from the perspectives of their synchronization pattern.The extra oscillatory behavior in the graph is simply due to the sinusoidal term embedded in K(t).Similar trends are observed even for the 3b, except for the fact that they no longer exhibit the coherent behavior.
(a): Temporal evolution of r(t) for K = 0.5 = K 1 > K c .Results are identical for the time-dependent and uniform coupling strength.(b): r(t) vs. t comparison using K = 0.2 & K(t) with K 1 = 0.2.These graphs are fairly similar.

Figure 4a,c display the time evolution of θ(t) and r(t) with variable A for
The oscillators remain out of synchrony except for some periodic peaks at which K(t) → K c .It is usually seen if amplitude A goes on increasing.
Similarly, Figure 4b,d show θ j (t) and r(t) vs. t plots for K(t) = K 1 + A sin(Ωt) with K 1 = 0.5, using different A values.As K 1 = 0.5, which is well above K c , the oscillators display synchronous behavior.However, oscillators start losing their rhythmic pattern when A increases significantly.The impact of various values Ω in K(t) was not significant in both regions, mainly from the perspective of reshaping the synchronization pattern.

Using Step Function as a Perturbation
An alternative approach that essentially replaces the sinusoidal function as a perturbation term is to provide a uniform but higher value of coupling strength K(t) for the time at which order parameter r(t) goes below some threshold value r th .A step function may fulfill such requirements and thus keep oscillators in their synchronized state persistently.A system of such oscillators displays autonomous behavior, except for r(t) ≤ r th , where a uniform but higher value of K(t) is applied to maintain the synchronized state of oscillators.
One important assumption is made in the numerical solution using the additional coupling strength K(t) in the form of a step function.

•
The threshold value r th of the order parameter selected in the r(t) vs. t graph, below which an additional value of coupling strength (K add ) is applied, is a tentative midvalue (0 ≤ r(t) ≤ 1).Oscillations observed in the r(t) vs. t graph above r th are assumed to be holding the synchronous behavior of oscillators.
The order parameter r(t) vs. t graph (Figure 5a) looks very random in the case of 32 oscillators.Additionally, the lack of precise information regarding the critical coupling strength K c gives the real difficulty in the numerical solution.The extrapolation technique by which special time points t 2 , which correspond to the r th , are estimated, fails to reflect the closer values.There is a discrepancy between actual time points corresponding to r th and the estimated t 2 (represented by red asterisk in Figure 5a).Figure 5c (corresponding K vs. t plot in Figure 5b) reflects such discrepancy, as r(t) occasionally goes below its threshold value r th = 0.2.However, such occasional drifting can be eliminated if the tolerance value is increased.The higher the tolerance, the more of these r(t) ≤ r th points are boosted with the additional coupling strength K add , thus ultimately, preserving the synchronous behavior.This will also increase the average value of the coupling strength K ave , as K add will be applied for a longer interval of time.Figure 5d, which is generated using t tol = t 2 ± 4, does not include such occasional drifting of the oscillators or, to be more precise, the order parameter r(t) at least does not go beyond r th = 0.25.(a): Estimating the t 2 values at r = r th for a 32 oscillator system.These points lack periodicity and are located at random distances in the r(t) vs. t curve.(b): K vs. t graph for 32 oscillator system after giving K add in the vicinity of r = r th = 0.2.The higher value of K is applied in a narrow interval, causing no noticeable increase in the K ave .(c): Respective r(t) vs. t graph for a 32 oscillator system after giving K add in the vicinity of r = r th = 0.2.Finally, the higher tolerance value ±4 is used in (d) of r(t) vs. t graphs.Now, every point in the graph lies above the r th due to the increased value of tolerance.

Conclusions and Discussions
This paper was focused on the synchronization of finite phase oscillators using timedependent coupling strength K(t).In contrast, the original Kuramoto model [7] was formulated for infinite oscillators using a homogeneous all-to-all coupling strength, K.
Keeping synchronization phenomena exhibited by different real-world examples, such as periodic flashing of fireflies, in consideration, the choice of using finite oscillators is a necessity.Additionally, considering the influence of external fields on the behavior of each oscillator in the form of a perturbation term makes the original Kuramoto model more realistic.Finding a good match between time-dependent and uniform cases is a required outcome, as real-life synchronization events have been in existence despite the influence of external factors.
A similar pattern observed in r(t) and θ(t) vs. t graphs for uniform and timedependent coupling strengths have a broad indication in nature: time variability is expected but well-adjusted by the oscillators displaying the synchronous behavior.A and Ω of the perturbation part had little impact on the r(t) vs. t graphs.However, the drifting of oscillators even with K = 0.5 > K c was noticeable for higher values of A (Figure 4b).
The numerical solution that involved the step function (Section 3.2) was carried out with the assumption that the persistence of such synchronized states can be achieved by supplying higher coupling strength K(t) = K initial + K add in the system of oscillators if the order parameter r(t) goes beyond some threshold r th .Such thresholds are userdefined, which means a limit can be set for the oscillator's phase difference from each other that can still be included in the group of oscillators that are displaying the synchronous behavior.For example, a firefly (oscillator) flashing with some period plus some margin can still be included in the main group of fireflies that are exhibiting synchronous behavior.Mathematically, they can be kept in one group if their frequencies are within some value of σ.The synchronous behavior of oscillators is always ensured for K(t) ≥ K c if K c is precisely known.Thus, Section 3.2 is simply trying to figure out whether such synchronized states are achievable for the coupling strength K(t) which is, in general, kept less than K c , except for the occasional cases (r(t) → r th ), where a higher K(t) acts for a brief period.
This approach, however, does not provide a conclusive picture.This is mainly due to the uncertainties associated with K c for 32 oscillators.The other parameters such as K initial , K add , r th , and t tol also come with some degree of uncertainties because of K c .However, such a limitation was essentially overcome using the higher tolerance values as in Figure 5d.

Figure 3 .
Comparison of r(t) vs. t graphs for K

Figure 4 .
θ(t) and r(t) vs. t graphs for different values of A, Ω = 1.Two regions, K 1 < K c and K 1 > K c , are tested.(a): K 1 = 0.2 < K c forces the oscillators to remain out of synchrony except for some with higher values of A. (b): Drifting of some oscillators are associated with the higher values of A (=1.2) even for K 1 = 0.5 > K c .It gives the important indication that the higher A values are destabilizing the synchronization pattern.(c,d) are the respective r(t) vs. t graphs for K 1 < K c and K 1 > K c .

Figure 5 .
Figure5.(a): Estimating the t 2 values at r = r th for a 32 oscillator system.These points lack periodicity and are located at random distances in the r(t) vs. t curve.(b): K vs. t graph for 32 oscillator system after giving K add in the vicinity of r = r th = 0.2.The higher value of K is applied in a narrow interval, causing no noticeable increase in the K ave .(c): Respective r(t) vs. t graph for a 32 oscillator system after giving K add in the vicinity of r = r th = 0.2.Finally, the higher tolerance value ±4 is used in (d) of r(t) vs. t graphs.Now, every point in the graph lies above the r th due to the increased value of tolerance.

Funding:
National Science Foundation (NSF)-HBCU-UP-Implementation Grant: Preparing the Pipeline of Next Generation STEM Professionals.(Award # HRD 2011917).Data Availability Statement: Not applicable.