# Dynamics of large oscillator populations with random interactions

Arkady Pikovsky  Institute of Physics and Astronomy, University of Potsdam,
Karl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm, Germany  Lev A. Smirnov
Department of Control Theory, Lobachevsky State University of Nizhny Novgorod,
23 Gagarin Avenue, Nizhny Novgorod 603022, Russia

###### Abstract

We explore large populations of phase oscillators interacting via random
coupling functions. Two types of coupling terms, the Kuramoto-Daido coupling
and the Winfree coupling, are considered. Under the assumption of statistical
independence of the phases and the couplings, we derive reduced averaged
equations with effective non-random coupling terms. As a particular example,
we study interactions that have the same shape but possess random coupling
strengths and random phase shifts. While randomness in coupling strengths just
renormalizes the interaction, a distribution of the phase shifts in coupling
reshapes the coupling function.

oscillatory ensembles, synchronization, disorder, random systems

> Populations of globally coupled oscillators appear in different fields of
> physics, engineering, and life sciences. In many situations, there is
> disorder in the coupling, and the coupling terms are not identical but vary,
> for example, due to different phase shifts and coupling strengths. We show
> that for large ensembles, one can effectively average over the disorder and
> describe the synchronization transition via simpler averaged equations.

##  I Introduction

Collective synchronization in oscillator populations attracted much interest
in the last decades [1, 2, 3, 4]. While the phenomenon is well-understood in a
regular situation, the influence of disorder remains a subject of intensive
current studies. The disordered case is relevant for many applications,
especially in neuroscience, where in the description of the correlated
activity of neurons [5, 6, 7], one can hardly assume the neurons themselves to
be identical and the coupling between them to be uniform.

Randomness in the form of quenched (time-independent) disorder can appear in
different forms. Already in the original paper by Kuramoto [2], a distribution
of the natural frequencies of oscillators was taken into account. In further
studies, effects of pairwise disorder in couplings have been formulated by a
random matrix of coupling strengths [8, 9, 10]. In this paper, we consider
another type of disorder in interactions, which is maximal in some sense: we
assume that all the pairwise coupling terms are different, taken from some
random distribution of random functions. In particular, the coupling functions
can have different shapes: some possess just one harmonic of the phases or of
the phase differences, and some can be more complex with many harmonics. Such
a setup accounts for a maximal heterogeneity of coupled units.

In numerical examples in this paper, we restrict to a subclass where the
coupling functions have the same shape but differ due to random coupling
strengths and to random phase shifts in the coupling [11]. Phase shifts in the
coupling lead to frustration, which can be understood as follows.
Synchronization is always possible for two coupled oscillators with different
phase shifts in the coupling terms. In this synchronous regime, the
frequencies of the oscillators coincide, but the phases are not equal; there
is some difference in the phases dependent on the phase shifts in the
coupling. If many oscillators interact with different phase shifts, optimal
phase relations in pairs may become conflicting, and the main problem is how
this frustration affects global synchrony. In a recent short communication
[12], we argued that in the thermodynamic limit of many coupled oscillators,
one can perform averaging over random phase shifts, which results in a new
effective coupling function. This paper extends this approach to different
setups and illustrates it with numerical simulations.

The paper is organized as follows. In Section II, we formulate different
setups for the phase dynamics of coupled oscillators and rotators. In Section
III we present our main theoretical result about reducing the disordered
situation to an effective regular averaged coupling function. A particular
situation where randomness is in the coupling strengths and in the phase
shifts is considered in Section IV. The validity of this reduction is
illustrated by numerical examples in Section V. In Section VI, we show how to
treat systems with partial disorder, where random phase shifts are attributed
to a driving or to a driven unit. We conclude with a discussion in Section
VII.

##  II Different regular setups for coupled phase oscillators

This section describes several popular models of interacting phase oscillators
without randomness in coupling. For our purposes, it is important to separate
the individual dynamics and the coupling terms, which we denote “c.t.” .

###  II.1 Winfree and Kuramoto-Daido coupling terms

The “standard” model describes the individual phase dynamics of an oscillator
as rotations with a natural frequency $\omega_{k}$, possibly with individual
noises $\sigma\xi_{k}(t)$:

$\dot{\varphi}_{k}=\omega_{k}+\sigma\xi_{k}(t)+\text{c.t.}\;.$ (1)

The model (1) is most popular because it can be directly derived for generic
coupled oscillators. However, in the literature, several similar particular
setups have been discussed. The model of coupled active rotators [13, 14, 15]
includes a non-uniform rotation of the variable $\varphi_{k}$ (which is thus
not interpreted as the oscillator phase):

$\dot{\varphi}_{k}=\omega_{k}-\nu\sin\varphi_{k}+\sigma\xi_{k}(t)+\text{c.t.}\,.$
(2)

Another case is where $\varphi_{k}$ is an angle variable (or the difference of
the phases of macroscopic wave functions for Josephson junctions) with a
second-order in time dynamics [16, 17, 18, 19] (recently, this model became
popular in the context of modeling power grids [20, 21]):

$\mu\ddot{\varphi}_{k}+\dot{\varphi}_{k}=\omega_{k}+\sigma\xi_{k}(t)+\text{c.t.}\,.$
(3)

Next, we specify the coupling terms “c.t.” according to the Winfree and the
Kuramoto-Daido approaches. In the Winfree model, the action on the oscillator
$k$ from the oscillator $j$ is proportional to the product
$S(\varphi_{k})Q(\varphi_{j})$, where $S(\varphi_{k})$ is the phase
sensitivity function of the unit $k$, and $Q(\varphi_{j})$ describes the force
with which the element $j$ is acting. The latter is usually renormalized
according to the number $N$ of interacting oscillators to have a well-defined
thermodynamic limit. The full Winfree-type coupling term reads [1, 22, 4]

$\frac{1}{N}S(\varphi_{k})\sum_{j}Q(\varphi_{j})\,.$ (4)

The Winfree model can be directly derived from the original equations
governing the oscillator dynamics, in the first order in the small parameter
describing the coupling strength [22, 4].

Further usage of this small parameter allows for a reduction of the Winfree
model (4) to the Kuramoto-Daido model. Because the coupling is weak, the phase
dynamics is represented by a fast rotation with frequency $\omega_{k}$ and a
slow variation due to the coupling. Because only the slow dynamics is
responsible for large deviations of the phases, one can perform averaging over
fast oscillations, keeping only the components with slow phase dependencies in
the coupling term. If the frequencies of the oscillators are nearly equal,
then the slow combinations of the phases are $\sim(\varphi_{j}-\varphi_{k})$.
Thus, keeping only slow terms yields the coupling term in the Kuramoto-Daido
form

$\frac{1}{N}\sum_{j}F(\varphi_{j}-\varphi_{k})\,.$ (5)

###  II.2 Formulation in terms of Kuramoto-Daido order parameters

Here, we show how the Kuramoto-Daido and the Winfree coupling functions can be
reformulated in terms of the Kuramoto-Daido order parameters. The latter are
defined as

$Z_{m}=\left\langle
e^{\mathrm{i}m\varphi_{j}}\right\rangle=\frac{1}{N}\sum_{j}e^{\mathrm{i}m\varphi_{j}}\,.$
(6)

We start with the Kuramoto-Daido-type coupling (5) and represent the
$2\pi$-periodic coupling function $F(x)$ as a Fourier series

$F(x)=\sum_{m}f_{m}e^{\mathrm{i}mx}\,,\qquad
f_{m}=\\!\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!dx\,F(x)e^{-\mathrm{i}mx}\,.$ (7)

Substituting this in the coupling term in (5) we obtain, using expression (6),

$\frac{1}{N}\sum_{j}F(\varphi_{j}-\varphi_{k})=\sum_{m}f_{m}e^{-\mathrm{i}m\varphi_{k}}\frac{1}{N}\sum_{j}e^{\mathrm{i}m\varphi_{j}}=\sum_{m}f_{m}e^{-\mathrm{i}m\varphi_{k}}Z_{m}\;.$
(8)

Similarly, for the Winfree case (4) we can represent the coupling functions
$Q(x)$ and $S(x)$ as Fourier series

$\begin{gathered}S(x)\\!=\\!\sum_{m}s_{m}e^{\mathrm{i}mx}\,,\qquad
s_{m}\\!=\\!\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!dx\,S(x)e^{-\mathrm{i}mx}\,,\\\
Q(x)=\sum_{l}q_{l}\,e^{\mathrm{i}lx}\,,\qquad
q_{l}=\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!dx\,Q(x)\,e^{-\mathrm{i}lx}\,.\end{gathered}$
(9)

Substitution in Eq. (4) allows for a representation of the coupling term as

$\frac{1}{N}\sum_{j}\sum_{m}s_{m}e^{\mathrm{i}m\varphi_{k}}\sum_{l}q_{l}e^{\mathrm{i}l\varphi_{j}}\\!=\\!\sum_{m}\sum_{l}s_{m}e^{\mathrm{i}m\varphi_{k}}\,q_{l}\left[\frac{1}{N}\sum_{j}e^{\mathrm{i}l\varphi_{j}}\right]\\!=\\!\sum_{m}s_{m}e^{\mathrm{i}m\varphi_{k}}\sum_{l}q_{l}Z_{l}\;.$
(10)

Below, we will use expressions (8) and (10) as “templates” for identifying the
effective coupling functions in the case of random interactions.

##  III General random coupling functions

We start with the maximally disordered setup, where all pairwise coupling
functions are random. We consider first the Kuramoto-Daido-type coupling (5)
and rewrite it assuming that all the coupling terms are generally different

$\frac{1}{N}\sum_{j}F_{jk}(\varphi_{j}-\varphi_{k})\,.$ (11)

We use the Fourier representation like (7):

$F_{jk}(x)=\sum_{m}f_{m,jk}e^{\mathrm{i}mx}\;,\qquad
f_{m,jk}=\frac{1}{2\pi}\int_{0}^{2\pi}dxF_{jk}(x)e^{-\mathrm{i}mx}\,.$ (12)

Here, we treat the complex Fourier coefficients $f_{m,jk}$ as random numbers
with some distribution. The randomness in coupling is a quenched (i.e., time-
independent) disorder.

We now perform the averaging of the coupling (11). First, we substitute (12)
in (11)

$\frac{1}{N}\sum_{j}\sum_{m}f_{m,jk}e^{\mathrm{i}(m\varphi_{j}-m\varphi_{k})}=\sum_{m}e^{-\mathrm{i}m\varphi_{k}}\left[\frac{1}{N}\sum_{j}f_{m,jk}e^{\mathrm{i}m\varphi_{j}}\right].$
(13)

We identify the term in the squared brackets as the population average
$\left\langle f_{m,jk}e^{\mathrm{i}m\varphi_{j}}\right\rangle$. Next, we
assume statistical independence of the phases and the Fourier coefficients
$f_{m,jk}$. We expect this independence to be valid for a large population,
where many different couplings influence each phase. This assumption allows us
for representing this average as $\left\langle
f_{m,jk}\right\rangle\left\langle
e^{\mathrm{i}m\varphi_{j}}\right\rangle=\left\langle f_{m,jk}\right\rangle
Z_{m}$. As a result, we obtain the coupling term

$\sum_{m}\left\langle f_{m,jk}\right\rangle
e^{-\mathrm{i}m\varphi_{k}}Z_{m}\,.$ (14)

Comparing this expression with (8), we conclude that the coupling is described
with an effective deterministic coupling function, Fourier modes of which are
just $\left\langle f_{m,jk}\right\rangle$.

In the case of the Winfree-type coupling (4), the case of general randomness
is one with all different phase sensitivity functions $S_{jk}(x)$ and forcing
functions $Q_{jk}(x)$:

$\frac{1}{N}\sum_{j}S_{jk}(\varphi_{k})Q_{jk}(\varphi_{j})\;.$ (15)

Again, we represent these functions via random complex Fourier coefficients
$s_{m,jk}$ and $q_{l,jk}$

$\begin{gathered}S_{jk}(x)\\!=\\!\sum_{m}s_{m,jk}e^{\mathrm{i}mx}\,,\qquad
s_{m,jk}\\!=\\!\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!dx\,S_{jk}(x)e^{-\mathrm{i}mx}\,,\\\
Q_{jk}(x)=\sum_{l}q_{l,jk}\,e^{\mathrm{i}lx}\,,\qquad
q_{l,jk}=\frac{1}{2\pi}\int_{0}^{2\pi}\\!\\!dx\,Q_{jk}(x)\,e^{-\mathrm{i}lx}\,.\end{gathered}$
(16)

Substituting (16) in (15) and assuming statistical independence of the phases
and the Fourier coefficients, we arrive at an effective deterministic coupling

$\sum_{m}\left\langle s_{m,jk}\right\rangle
e^{\mathrm{i}m\varphi_{k}}\sum_{l}\left\langle q_{l,jk}\right\rangle Z_{l}\;.$
(17)

Comparison of this expression with (10) shows that the Fourier modes of the
effective phase sensitivity function and of the forcing are just the averaged
random Fourier modes.

The last remark is that because the Fourier transform is a linear operation,
averaging the Fourier coefficients is the same as averaging the functions.
Thus, the effective averaged coupling function in the Kuramoto-Daido case (11)
is

$\frac{1}{N}\sum_{j}F_{jk}(\varphi_{j}-\varphi_{k})\;\Rightarrow\;\frac{1}{N}\sum_{j}\mathcal{F}(\varphi_{j}-\varphi_{k})=\frac{1}{N}\sum_{j}\left\langle
F_{jk}(\varphi_{j}-\varphi_{k})\right\rangle\;.$ (18)

For the random Winfree coupling (15), we have

$\begin{gathered}\frac{1}{N}\sum_{j}S_{jk}(\varphi_{k})Q_{jk}(\varphi_{j})\;\Rightarrow\;\mathcal{S}(\varphi_{k})\frac{1}{N}\sum_{j}\mathcal{Q}(\varphi_{j})\;,\\\
\mathcal{S}(\varphi_{k})=\left\langle
S_{jk}(\varphi_{k})\right\rangle,\quad\mathcal{Q}(\varphi_{j})=\left\langle
Q_{jk}(\varphi_{j})\right\rangle\;.\end{gathered}$ (19)

##  IV Randomness in coupling strengths and in the phase shifts

The case of general randomness of interactions, presented in Section III,
includes a situation where different interactions have different shapes. For
example, some oscillators can be coupled via the first harmonic coupling
function, while others are coupled with the second harmonic coupling function.
A particular situation is one where all the shapes are the same, but the
interactions differ in their phase shifts and the coupling strengths. For the
Kuramoto-Daido couplings (11), this means that

$F_{jk}(\varphi_{j}-\varphi_{k})=A_{jk}F(\varphi_{j}-\varphi_{k}-\alpha_{jk})\,.$
(20)

Here, random numbers $A_{jk}$ describe different random coupling strengths,
and $\alpha_{jk}$ are random phase shifts. For simplicity of presentation, we
assume statistical independence of $A_{jk}$ and $\alpha_{jk}$. Then, the
random coupling (20) is described by two distributions. We denote the
distribution of the coupling strengths $U(A)$ and the distribution of the
phase shifts $G(\alpha)$. This latter distribution is defined on a circle and
can be represented by the Fourier series:

$G(\alpha)=\frac{1}{2\pi}\sum_{m}g_{m}e^{\mathrm{i}m\alpha}\,,\quad
g_{m}=\\!\int_{0}^{2\pi}\\!\\!d\alpha\,G(\alpha)e^{-\mathrm{i}m\alpha}\,.$
(21)

Such pairwise phase shifts naturally appear if there are time delays for the
pairwise couplings. As has been demonstrated in Ref. 23, a phase shift
$\alpha$ is determined in the leading order in coupling strength as follows:
$\alpha=\omega\tau$, where $\tau$ is the signal propagation time and $\omega$
is the characteristic frequency (e.g., a median frequency if the natural
frequencies are different). Remarkably, phase shifts like in (20) naturally
appear in the description of power grids (cf. Eq. (3) in Ref. 21).

Applying general expression (18) to this particular case, we obtain the
effective coupling function as

$\mathcal{F}(x)=\left\langle
A\right\rangle\int_{0}^{2\pi}F(x-\alpha)G(\alpha)\;d\alpha\;.$ (22)

where $\left\langle A\right\rangle=\int dA\,A\,U(A)$. One can see that the
randomness of coupling strengths renormalizes the total coupling strength, but
does not influence the shape of the coupling function. In contradistinction,
the randomness of the phase shifts changes the form of the coupling function
because of the convolution operator in (22). The best way to see this is to
return to the Fourier representation, where the effect of the randomness in
the phase shifts reduces to a factorization of the Fourier modes

$f_{m}\;\Rightarrow\;f_{m}g_{m}\;.$ (23)

One can see that some modes in the coupling function can even disappear if the
corresponding factors $g_{m}$ vanish. We will explore such cases in Section V
below.

In the context of the Winfree model (4), one can also restrict randomness to
coupling strengths and phase shifts so that the shapes of the phase
sensitivity function and of the forcing remain the same. We have, in general,
two phase shifts $\alpha_{jk}$ and $\beta_{jk}$, entering the phase
sensitivity function and the driving term, respectively, and coupling
strengths $A_{jk}$ (for simplicity of presentation we consider all the three
random quantities as statistically independent):

$\frac{1}{N}\sum_{j}A_{jk}S(\varphi_{k}-\alpha_{jk})Q(\varphi_{j}-\beta_{jk})\,.$
(24)

Of course, one of these phase shifts may be absent in a particular situation.
The distribution $B(\beta)$ of phase shifts $\beta$ is represented similarly
to (21):

$B(\beta)=\frac{1}{2\pi}\sum_{m}b_{m}e^{\mathrm{i}m\beta},\quad
b_{m}=\int_{0}^{2\pi}d\beta\,B(\beta)e^{-\mathrm{i}m\beta}\,.$ (25)

Now, we apply general expressions (19) to the case of randomness in phase
shifts and coupling strengths and obtain:

$\begin{gathered}\mathcal{S}(x)=\left\langle A\right\rangle\left\langle
S(x-\alpha_{jk})\right\rangle=\left\langle
A\right\rangle\int_{0}^{2\pi}S(x-\alpha)G(\alpha)\;d\alpha\;,\\\
\mathcal{Q}(y)=\left\langle
Q(y-\beta_{jk})\right\rangle=\int_{0}^{2\pi}Q(y-\beta)B(\beta)\;d\beta\;.\end{gathered}$
(26)

Here, we, somewhat arbitrarily, attributed the average coupling strength to
the phase sensitivity function. Again, like in the case of the Kuramoto-Daido-
type coupling, the coupling strengths do not influence the shape of the
functions, while the distributions of the phase shifts do reshape them. In
terms of Fourier modes, one has a factorization by the modes of the
distributions:

$s_{m}\;\Rightarrow\;s_{m}g_{m}\;,\qquad q_{l}\;\Rightarrow q_{l}b_{l}\;.$
(27)

##  V Numerical examples

The theory above predicts that in the presence of random phase shifts in the
coupling, the effective coupling function is the convolution of the original
one with the phase shift distribution density. In terms of the Fourier modes,
one has a product of the modes of the original coupling function with the
Fourier modes of the distribution density of the phase shifts. The most
prominent effect appears if the original coupling function is rather complex
(has several harmonics), but the distribution of the phase shifts is simple,
possessing just one harmonic. For example, we take the density of the phase
shifts $G(\alpha)$ in the form ($M$ is an integer)

$G(\alpha)=\frac{1}{2\pi}(1+\cos M\alpha)\;.$ (28)

In the Fourier representation, this corresponds to one nontrivial Fourier mode
$g_{M}=0.5$. Correspondingly, only the $M$-th harmonic is present in the
effective coupling function with such a distribution of the phase shifts.

Figure 1: Estimations of the phase probability densities via histograms (200
bins are used) in model (29), (28) with $M=1$ (red lines), $M=2$ (green
lines), and $M=3$ (blue lines). These distributions are almost the same for
different system sizes $N$.

In Fig. 1, we show the distributions of the phases in globally coupled
populations of noisy identical phase oscillators ($\omega_{k}=0$) with the
Kuramoto-Daido coupling

$\dot{\varphi}_{k}=\sigma\xi_{k}(t)+\frac{1}{N}\sum_{j}F(\varphi_{j}-\varphi_{k}-\alpha_{jk}),\quad
F(x)=\sin(x)+2\sin(2x)+3\sin(3x)\;.$ (29)

We show the results for $\sigma=0.12$ and three distributions of type (28)
with $M=1,2,3$. One can see that in these cases, the densities possess the
corresponding $(2\pi/M)$-periodicities as functions of $\varphi$, confirming
that the phase shifts result in the effective one-mode coupling.

Figure 2: Dashed and dotted lines (order parameters
$\left\langle\left|Z_{1}\right|\right\rangle=\left\langle\left|\left\langle
e^{\mathrm{i}\varphi_{j}}\right\rangle\right|\right\rangle$ and
$\left\langle\left|Z_{2}\right|\right\rangle=\left\langle\left|\left\langle
e^{2\mathrm{i}\varphi_{j}}\right\rangle\right|\right\rangle$, respectively):
the behavior of an ensemble of noisy identical oscillators with Kuramoto-Daido
coupling in the form (31) with $K=1$ in dependence on the noise intensity
$\sigma^{2}$ in the absence of phase shifts in coupling. Markers: order
parameter $\left\langle\left|Z_{1}\right|\right\rangle$ in simulation of a
disordered population with the phase shifts $\alpha_{jk}$ sampled according to
(28) with $M=1$ ($N=500$: red pluses; $N=1000$: cyan circles, $N=2000$: brown
squares). Here and below, we also apply time averaging when calculating
complex mean fields to smooth out small-magnitude, non-regular fluctuations
caused by noise and finite-size effects. Solid green line: theoretical
prediction of the first order parameter $\left|Z_{1}\right|$ in the reduced
model according to Ref. 18. Figure 3: Panel (a): order parameter
$\left\langle\left|Z_{1}\right|\right\rangle$ for ensembles with random phase
shifts $\alpha_{jk}$ taken from the one-mode distribution (28) with $M=1$,
random natural frequencies $\omega_{k}$ drawn from the Cauchy distribution
(30) with $\Delta=1$, and coupling function (31). Dashed line: theoretical
prediction of the order parameter according to Ref. 2:
$|Z_{1}|=\sqrt{(K-4)\bigl{/}K\bigr{.}}$. Panels (b, c): averaged over time
absolute values of the second and third circular cumulants
$\kappa_{2}\\!=\\!\left\langle\left|Z_{2}-Z_{1}^{2}\right|\right\rangle$ and
$\kappa_{3}\\!=\\!\left\langle\left|Z_{3}-3Z_{2}Z_{1}+2Z_{1}^{3}\right|\right\rangle$,
respectively; their small values confirm that the distribution of the phases
is a wrapped Cauchy distribution.

A more quantitative test can be performed if an exact analytical solution
exists for the reduced system. Such solutions are known for the Kuramoto-Daido
system with only the first harmonic in the coupling function, in the cases of
noisy identical oscillators [18] and for deterministic oscillators with a
Cauchy distribution of natural frequencies [2]:

$W(\omega)=W_{C}(\omega-\Omega),\qquad
W_{C}(\varpi)=\Delta\Bigl{/}\pi\left(\varpi^{2}+\Delta^{2}\right)\Bigr{.}\,,$
(30)

where $\Omega$ is the mean value, and $\Delta$ is the parameter governing the
width of the distribution. In the next two examples, we fix $M=1$ in (28) and
take a two-harmonic original coupling function

$F(x)=K\bigl{(}\sin(x)+4\sin(2x)\bigr{)}\,.$ (31)

The effective coupling function is thus $\mathcal{F}(x)\\!=\\!0.5\,K\sin(x)$,
so that the analytical results mentioned are applicable. Figure 2 shows the
case of noisy oscillators. Figure 3 shows the case of deterministic
oscillators with the Cauchy distribution of the natural frequencies (the width
parameter of the Cauchy distribution is $\Delta\\!=\\!1$). In both cases, the
dynamics of a system with random phase shifts nicely corresponds to the
analytically predicted behavior of the reduced system. For the deterministic
case, there is an additional (to the comparison of the order parameter with
the theoretical prediction) indicator for the validity of the effective
coupling. As it follows from the Ott-Antonsen [24] solution of the problem,
the distribution of the phases in the thermodynamic limit $N\to\infty$ is a
wrapped Cauchy distribution. The most direct check of this prediction is the
calculation of the higher circular cumulants of the distribution [25], which
should vanish for the wrapped Cauchy distribution. We show the absolute values
of the second and the third cumulants in Fig. 3. Their values are indeed small
compared to the first cumulant (which is the Kuramoto order parameter
$Z_{1}$), and this smallness is even improved as the size of the population
grows.

We illustrate the case of disorder in coupling strengths in Fig. 4. As already
can be seen in the case of purely phase shifts disorder Fig. 3, one needs
large system sizes $N$ to become closer to the theoretical prediction. This
effect is even more pronounced when both the phase shifts and the coupling
strengths are random. Therefore, for the same setup as presented in Fig. 3, we
performed calculations for a selected coupling strength $K\\!=\\!0.625$, for
which theory predicts $|Z_{1}|\\!=\\!0.6$. We added disorder in coupling
strengths in such a way that $\left\langle A_{jk}\right\rangle\\!=\\!1$, and
followed deviations of the obtained values of
$\left\langle|Z_{1}|\right\rangle$ from the theoretical prediction. As it
follows from Fig. 4, the deviations decrease with $N$ roughly as
$\sim\\!N^{-1}$, although they are rather pronounced for small
$N\lesssim\\!1000$. We tested two distributions of $A_{jk}$, one the uniform
in the interval $0\leq A_{jk}\\!\leq\\!2$, and one bimodal, where $A_{jk}$
takes values $0$ and $2$ with probability $1/2$. In the latter case, finite-
size deviations from the theoretical prediction are stronger.

Figure 4: Deviations from the theoretical prediction of the main order
parameter vs. ensemble size $N$. The same coupling function and phase disorder
as in Fig. 3 were used but with an additional disorder in coupling strengths.
Red boxes: uniform distribution of $A_{jk}$ in the interval $0\leq
A_{jk}\\!\leq\\!2$; blue circles: $A_{jk}$ takes values $0$ or $2$ with
probability $1/2$. The black dotted line shows the law $\\!\sim\\!N^{-1}$.
Figure 5: Oscillatory partially synchronized regime the noiseless case the
Kuramoto-Daido model with two components of the symmetric bimodal distribution
having the Cauchy shapes:
$W(\omega)=\bigl{(}W_{C}(\omega+\Omega)+W_{C}(\omega-\Omega)\bigr{)}\bigl{/}2\bigr{.}$,
where $\Omega=0.25$ and $\Delta=0.125$. Bold blue solid lines: the dynamics of
the order parameter $\left\langle|Z_{1}|\right\rangle$ in the systems of
$N=64\times 10^{3}$ oscillators with three coupling functions: panel (a):
$F(x)=1.25\sin(x)+0.75\sin(3x)\bigr{)}$, panel (b):
$F(x)=1.25\sin(x)+2.25\sin(2x)\bigr{)}$, and panel (c):
$F(x)=1.25\sin(x)+2.25\sin(3x)\bigr{)}$. Thin red dashed curves: the dynamics
corresponding to a limit cycle in the ordinary differential equations for the
two subgroup order parameters, which can be obtained for the Kuramoto model
with the effective coupling $\mathcal{F}(x)=0.625\sin(x)$ using the Ott-
Antonsen approach in the thermodynamic limit.

Our following example is a system with a bimodal distribution of natural
frequencies [26, 27, 28, 29]. Now, the corresponding dynamics may indeed be
more complicated than for an unimodal distribution, with a region of
bistability between asynchronous and synchronous states. A partially
synchronous state may be characterized either by a constant or oscillating
order parameter. We consider the noiseless case of the Kuramoto-Daido model
with several harmonics in the coupling function (7) and with two components of
the symmetric bimodal distribution $W(\omega)$ having the Cauchy shapes, i.e.
$W(\omega)=\bigl{(}W_{C}(\omega+\Omega)+W_{C}(\omega-\Omega)\bigr{)}\bigl{/}2\bigr{.}$.
According to our analysis, with the random phase shifts $\alpha_{jk}$ drawn
from the one-mode distribution (28) with $M=1$, the model under consideration
can be simplified to a system with an effective coupling function
$\mathcal{F}(x)\sim\sin(x)$ that contains only the first harmonic. The
dynamics of the latter oscillator systems with a symmetric bimodal
distribution was studied in detail [26, 27, 28, 29]. In particular, in Ref.
27, possible dynamical regimes and bifurcation between them were
comprehensively analyzed using the macroscopic approaches based on the Ott-
Antonsen ansatz.

In Fig. 5, we show the oscillatory regimes for the Kuramoto-Daido models with
three different coupling functions, each having two harmonics. However, the
effective averaged dynamics is the same, which is clearly seen in the Figure.
It is worth mentioning that for good quantitative (not just qualitative)
agreement, the number $N$ of ensemble units should be quite large; we used
$N=64\times 10^{3}$.

Figure 6:  Behavior of the first order parameter
$\left\langle\left|Z_{1}\right|\right\rangle$ an ensemble of
$N\\!=\\!12\\!\times\\!10^{3}$ noisy rotators (3) with equal natural
frequencies ($\omega_{k}=\Omega$) and coupling function (31) with $K\\!=\\!1$
in dependence on $\sqrt{\sigma}$ for different values of the moment of inertia
$\mu$: $(\mathrm{a})\,\mu=0.02$, $(\mathrm{b})\,\mu=0.1$, and
$(\mathrm{c})\,\mu=0.5$. Green squares and blue circles are simulations
without and with phase shifts (taken from a distribution (28) with $M=1$),
respectively. Solid red lines show the theoretical prediction for this order
parameter according to Eq. (32).

Next, we present simulations of the rotator model (3). We consider the case of
noisy oscillators with equal natural frequencies and coupling function (31).
We consider random phase shifts $\alpha_{jk}$ distributed according to (28)
with $M=1$. Thus, the effective coupling function is
$\mathcal{F}(x)=0.5\,K\sin(x)$. For such coupling, the analytical expression
for the order parameter in dependence on the noise intensity $\sigma^{2}$ can
be written in the parametric (parameter $R$) form [18]

$|Z_{1}|=\frac{2\pi RI_{0}^{2}(R)I_{1}(R)}{2\pi RI_{0}^{2}(R)+\mu
KI_{1}(R)},\quad\sigma^{2}=\frac{K|Z_{1}|}{2R}\;.$ (32)

Here, $I_{0}(R)$ and $I_{1}(R)$ are the principal branches of the modified
Bessel functions of the first kind with orders $0$ and $1$, respectively.
Relation (32) is valid for small $\mu$, in the first order in this parameter.
Note, for a vanishing $\mu$, we obtain the equations determining the solid
green line in Fig. 2. In Fig. 6, we compare the dynamics of a noisy ensemble
of rotators with and without random phase shifts. One can see that the effect
of the moment of inertia $\mu$ on the coherence level of the established state
is relatively small.

##  VI Reduction in the case of global random phase shifts

The theory above was based on the assumption of independence of the phases and
the phase shifts, which is plausible if there are many different phase shifts
for a given oscillator, like for random pairwise phase shifts. Here, we
consider a slightly different setup, where phase shifts $\alpha_{jk}$ are not
attributed to connections $j\to k$, but separately to the driven or to the
driving system; we call such a situation “global random phase shifts”. Thus,
these phase shifts have one index instead of two. In the case of global phase
shifts, the averaging requires additional justification. Such an analysis has
been performed in Ref. 12 for the Kuramoto-Daido setup; here, we present a
similar consideration for the Winfree model with noise.

We consider noisy oscillators with global random phase shifts in coupling

$\dot{\varphi}_{k}=\sigma\xi_{k}(t)+\frac{1}{N}S(\varphi_{k}-\alpha_{k})\sum_{j}Q(\varphi_{j}-\beta_{j})\;,$
(33)

where $\xi_{k}(t)$ are independent white Gaussian forces with zero mean
$\left\langle\xi\right\rangle=0$, and auto-correlation
$\left\langle\xi_{k}(t^{\prime})\,\xi_{j}(t+t^{\prime})\right\rangle\\!=\\!2\delta_{jk}\delta(t)$.
The phase shifts $\alpha_{k}$ are attributed to the driven systems; they
characterize the phase shifts in the phase sensitivity function independently
of the driving. On the other hand, the phase shifts $\beta_{j}$ are attributed
to the driving units; they characterize the force functions. Remarkably, by a
change of variables $\theta_{k}=\varphi_{k}-\alpha_{k}$, one can transfer both
shifts to the driving unit:

$\dot{\theta}_{k}=\sigma\xi_{k}(t)+\frac{1}{N}S(\theta_{k})\sum_{j}Q(\theta_{j}+\alpha_{j}-\beta_{j})$
(34)

We thus denote $\gamma_{j}=\beta_{j}-\alpha_{j}$ and consider it as a single
combined global random phase shit with probability density
$\mathit{\Gamma}(\gamma)$.

In the thermodynamic limit $N\to\infty$, we describe the population with the
probability density $P(\theta,t|\gamma)$, which generally can depend on the
phase shift $\gamma$. The Fokker-Planck equation for this density reads

$\frac{\partial P}{\partial
t}+\frac{\partial}{\partial\theta}\Bigl{(}S(\theta)\bar{Q}(t)\Bigr{)}=\sigma^{2}\frac{\partial^{2}P}{\partial\theta^{2}}\;,$
(35)

where

$\bar{Q}(t)=\\!\int_{0}^{2\pi}\\!\\!d\theta\int_{0}^{2\pi}\\!\\!d\gamma\,Q(\theta-\gamma)\mathit{\Gamma}(\gamma)P(\theta,t|\gamma)$
(36)

contains the averaging over the phases and the phase shifts. The crucial
observation is that, although $P(\theta,t|\gamma)$ can potentially depend on
$\gamma$, for example, an initial condition at time $t=0$ can contain a
$\gamma$-dependence, the equation for $P(\theta,t|\gamma)$ does not contain
$\gamma$. Thus, if this equation, for any prescribed time-dependent function
$\bar{Q}(t)$, has a unique attracting solution $P(\theta,t)$ which does not
depend on the initial conditions, then the statistical dependence of the
phases on the phase shits $\gamma$ disappears, even if it was presented
initially. The property to have a unique solution is sometimes called “Global
Asymptotic Stability”. It is rather natural for the dissipative parabolic
equation (35) on the circle. Although we found only proofs for time-
independent $\bar{Q}$ in the literature, see Refs. 30, 31, it appears that
such a proof can be extended to the nonstationary case too [32]. For a master
equation, the global asymptotic stability has been established in Ref. 33.

Taking into account that $P(\theta,t|\gamma)\to P(\theta,t)$ at large times,
we insert the $\gamma$-independent density in Eq. (36) and obtain

$\bar{Q}=\\!\int_{0}^{2\pi}\\!\\!d\theta\,P(\theta,t)\mathcal{Q}(\theta),\qquad\mathcal{Q}(\theta)=\\!\int_{0}^{2\pi}\\!\\!d\gamma\,\mathit{\Gamma}(\gamma)Q(\theta-\gamma)\,.$
(37)

We see that the forcing function reduces to an effective phase-shift-
independent forcing function $\mathcal{Q}(\theta)$, which is the convolution
of the original coupling function and the distribution of the phase shifts.
This result parallels expression (19). A similar expression holds for the
Kuramoto-Daido coupling (cf. Ref. 12).

Figure 7: The average values of the order parameter
$\left\langle|Z_{1}|\right\rangle$ vs the parameter $K$ of the coupling
function $F(x)=K\bigl{(}\sin(x)-1.5\sin(2x)+0.5\cos(2x)\bigr{)}$ for two
values of the width parameter $\Delta$ of the Cauchy distribution (30) of the
frequencies $\omega_{k}$: (a) $\Delta=0.2$ and (b) $\Delta=0.4$. Green boxes:
without phase shifts; blue circles: with phase shifts. The red lines show the
theoretical prediction for the order parameter $\left|Z_{1}\right|$ according
to Ref. 18. System size $N=48\times 10^{4}$.

In the example Fig. 7 we consider for the coupled rotators (3) (parameter of
inertia $\mu\\!=\\!0.1$), we include both a Cauchy distribution of natural
frequencies and noise with amplitude $\sigma\\!=\\!0.05$. Here, we take the
global phase shifts according to $F(\varphi_{j}-\varphi_{k}-\alpha_{j})$. In
Fig. 7, we choose the coupling function in the form
$F(x)=K\bigl{(}\sin(x)-1.5\sin(2x)+0.5\cos(2x)\bigr{)}$ and present the
average values of the order parameter $\left\langle|Z_{1}|\right\rangle$ for
simulations without phase shifts and for phase shifts distributed according to
(28) with $M=1$. One can see that the latter data nicely fits the theoretical
prediction of Ref. 18.

##  VII Conclusion

We first summarize the results of this paper. We have considered different
models of globally coupled phase oscillators and rotators. In the case of a
“maximal disorder”, all the coupling functions are distinct and random,
sampled from some distribution. Based on the assumption of independence of the
phases and the coupling functions in the thermodynamic limit, we derived the
averaged equations for the phases, where effective deterministic coupling
functions enter. A more detailed consideration was devoted to the case where
the shapes of the random coupling functions are the same, but the amplitudes
and the phase shifts are random. Then, the effective functions are
renormalized convolutions of the original coupling functions and the
distribution densities of the phase shifts. In the Fourier representation, the
Fourier modes of the coupling functions are multiplied by Fourier modes of the
distribution densities of the phase shifts. This means that the effective
averaged coupling is “simpler” than the original one. In particular, if the
distribution of the phase shifts possesses just one Fourier mode, the
effective coupling function will possess only this mode, too. This property
allows us to check the validity of the approach numerically because, for the
one-mode coupling function, there are theoretical predictions for the behavior
of the order parameters.

A special case is a maximally frustrating one, where the averaged coupling
function vanishes. For the coupling strengths disorder, this happens in the
Daido model [34]; in the case of random phase shifts, this occurs for a
uniform in the range $0\leq\alpha_{jk}\leq 2\pi$ distribution of the phase
shifts. Our theory predicts that in the thermodynamic limit the interaction
vanishes. However, certain synchronization phenomena can still be observed in
finite ensembles with random phase shifts, as demonstrated recently in Ref.
35. We also expect that for other distributions there can be pronounced
deviations from the averaged behavior for finite ensembles. This issue
definitely deserves further exploration.

###### Acknowledgements.

L.A.S. acknowledges support from the Russian Science Foundation (grant no.
22-12-00348) and the Ministry of Science and Education of Russian Federation
(project no. FSWR-2024-0005).

## Data availability

All numerical experiments are described in the paper and can be reproduced
without additional information.

## References

  * Winfree [1967] A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” J. Theor. Biol. 16, 15–42 (1967).
  * Kuramoto [1975] Y. Kuramoto, “Self-entrainment of a population of coupled nonlinear oscillators,” in  _International Symposium on Mathematical Problems in Theoretical Physics_ , edited by H. Araki (Springer Lecture Notes Phys., v. 39, New York, 1975) p. 420.
  * Strogatz and Stewart [1993] S. H. Strogatz and I. Stewart, “Coupled oscillators and biological synchronization,” Scientific American , 68–75 (1993).
  * Pikovsky, Rosenblum, and Kurths [2001] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization. A Universal Concept in Nonlinear Sciences._ (Cambridge University Press, Cambridge, 2001).
  * Glass and Mackey [1988] L. Glass and M. C. Mackey,  _From Clocks to Chaos: The Rhythms of Life._ (Princeton Univ. Press, Princeton, NJ, 1988).
  * Buzsáki [2006] G. Buzsáki,  _Rhythms of the brain_ (Oxford UP, Oxford, 2006).
  * Daffertshofer and Pietras [2020] A. Daffertshofer and B. Pietras, “Phase synchronization in neural systems,” Synergetics , 221–233 (2020).
  * Kalloniatis [2010] A. C. Kalloniatis, “From incoherence to synchronicity in the network Kuramoto model,” Physical Review E 82, 066202 (2010).
  * Chiba, Medvedev, and Mizuhara [2018] H. Chiba, G. S. Medvedev, and M. S. Mizuhara, “Bifurcations in the Kuramoto model on graphs,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (2018).
  * Juhász, Kelling, and Ódor [2019] R. Juhász, J. Kelling, and G. Ódor, “Critical dynamics of the Kuramoto model on sparse random networks,” Journal of Statistical Mechanics: Theory and Experiment 2019, 053403 (2019).
  * Park, Rhee, and Choi [1998] K. Park, S. W. Rhee, and M. Y. Choi, “Glass synchronization in the network of oscillators with random phase shifts,” Phys. Rev. E 57, 5030–5035 (1998).
  * Smirnov and Pikovsky [2024] L. A. Smirnov and A. Pikovsky, “Dynamics of oscillator populations globally coupled with distributed phase shifts,” Phys. Rev. Lett. 132, 107401 (2024).
  * Shinomoto and Kuramoto [1986] S. Shinomoto and Y. Kuramoto, “Phase transitions in active rotator systems,” Prog. Theor. Phys. 75, 1105–1110 (1986).
  * Sakaguchi, Shinomoto, and Kuramoto [1988] H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, “Phase transitions and their bifurcation analysis in a large population of active rotators with mean-field coupling,” Prog. Theor. Phys. 79, 600–607 (1988).
  * Klinshov  _et al._ [2021] V. V. Klinshov, S. Y. Kirillov, V. I. Nekorkin, and M. Wolfrum, “Noise-induced dynamical regimes in a system of globally coupled excitable units,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021).
  * Tanaka, Lichtenberg, and Oishi [1997] H. Tanaka, A. Lichtenberg, and S. Oishi, “First order phase transition resulting from finite inertia in coupled oscillator systems,” Phys. Rev. Lett. 78, 2104–2107 (1997).
  * Hong  _et al._ [1999] H. Hong, M. Y. Choi, J. Yi, and K.-S. Soh, “Inertia effects on periodic synchronization in a system of coupled oscillators,” Phys. Rev. E 59, 353–363 (1999).
  * Munyaev  _et al._ [2020] V. O. Munyaev, L. A. Smirnov, V. A. Kostin, G. V. Osipov, and A. Pikovsky, “Analytical approach to synchronous states of globally coupled noisy rotators,” New Journal of Physics 22, 023036 (2020).
  * Munyayev  _et al._ [2023] V. O. Munyayev, M. I. Bolotov, L. A. Smirnov, G. V. Osipov, and I. Belykh, “Cyclops states in repulsive kuramoto networks: The role of higher-order coupling,” Phys. Rev. Lett. 130, 107201 (2023).
  * Filatrella, Nielsen, and Pedersen [2008] G. Filatrella, A. H. Nielsen, and N. F. Pedersen, “Analysis of a power grid using a Kuramoto-like model,” Eur. Phys. J. B 61, 485–491 (2008).
  * Dorfler and Bullo [2012] F. Dorfler and F. Bullo, “Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators,” SIAM Journal on Control and Optimization 50, 1616–1642 (2012).
  * Kuramoto [1984] Y. Kuramoto,  _Chemical Oscillations, Waves and Turbulence_ (Springer, Berlin, 1984).
  * Izhikevich [1998] E. M. Izhikevich, “Phase models with explicit time delays,” Phys. Rev. E 58, 905–908 (1998).
  * Ott and Antonsen [2008] E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems of globally coupled oscillators,” CHAOS 18, 037113 (2008).
  * Tyulkina  _et al._ [2018] I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, and A. Pikovsky, “Dynamics of noisy oscillator populations beyond the Ott-Antonsen ansatz,” Phys. Rev. Lett. 120, 264101 (2018).
  * Bonilla, Vicente, and Spigler [1998] L. L. Bonilla, C. J. P. Vicente, and R. Spigler, “Time-periodic phases in populations of nonlinearly coupled oscillators with bimodal frequency distributions,” Physica D: Nonlinear Phenomena 113, 79 – 97 (1998).
  * Martens  _et al._ [2009] E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, and T. M. Antonsen, “Exact results for the Kuramoto model with a bimodal frequency distribution,” Phys. Rev. E 79, 026204 (2009).
  * Campa [2020] A. Campa, “Phase diagram of noisy systems of coupled oscillators with a bimodal frequency distribution,” Journal of Physics A: Mathematical and Theoretical 53, 154001 (2020).
  * Kostin  _et al._ [2023] V. A. Kostin, V. O. Munyaev, G. V. Osipov, and L. A. Smirnov, “Synchronization transitions and sensitivity to asymmetry in the bimodal Kuramoto systems with Cauchy noise,” Chaos: An Interdisciplinary Journal of Nonlinear Science 33, 083155 (2023).
  * Gardiner [1996] C. W. Gardiner,  _Handbook of Stochastic Methods_ (Springer, Berlin, 1996).
  * Calogero [2012] S. Calogero, “Exponential convergence to equilibrium for kinetic Fokker-Planck equations,” Comm. Part. Diff. Eqs. 37, 1357–1390 (2012).
  * [32] S. Zelik (private communication, 2023).
  * Earnshaw and Keener [2010] B. A. Earnshaw and J. P. Keener, “Global asymptotic stability of solutions of nonautonomous master equations,” SIAM Journal on Applied Dynamical Systems 9, 220–237 (2010).
  * Daido [1992] H. Daido, “Quasientrainment and slow relaxation in a population of oscillators with random and frustrated interactions,” Phys. Rev. Lett. 68, 1073–1076 (1992).
  * Pikovsky and Bagnoli [2024] A. Pikovsky and F. Bagnoli, “Dynamics of oscillator populations with disorder in the coupling phase shifts,” New Journal of Physics 26, 023054 (2024).


# Voronoi diagrams for the distributed sensor network system data processing
111Link to the journal Future Generation Computer Systems.

Kondybayeva Almagul  Giovanna Di Marzo  Carouge 1212, Geneve  Centre
Universitaire d’Informatique, Universite de Geneve<EMAIL_ADDRESS>Carouge, Universite de Geneve, Switzerland

###### Abstract

This article represents the computational model for spacial addresation of the
sensors in the dynamically changing real-time internet of things system. The
model bases on the Voronoi diagrams as a basic data structure. Problem - the
correct data addresation without time delays in real-time processing and
database indexation in distributed storages. Goal - to develop the real-time
processing model of the object location identification in the Voronoi diagram
data structure. Relevance - the research and development on the contemporary
issues on the convergence (the N limit up to which the model presents the
algorithm convergence), time-delay, correct indexation in the database
transactions and adressation throughout wireless radio frequencies in the
distribiuted sensor systems. Solution proposes the Voronoi diagram
computational geometry data structure and the sweeping curve algorithm.
Methods represents the following steps: simulation of the dynamically changing
agent system using set of points that are based on a contemporary paths of the
public transport routes, bykes, vehicles circulations; 3D map rendering on the
SITG Canton of Geneva map; Voronoi diagrams calculation with the Fortune’s
algorithm. Results - this work presents 2D static and dynamic Fortune’s
realization of the Voronoi diagrams, the architecture of the model for
constructing a distributed sensor system based on Voronoi diagrams is
described. Scope - geographic charts showing the distribution of parameters,
determination of the nearest points from the object, spacial arrangement of
objects in nature, solving the problem of finding all nearest neighbors to the
object, database indexing and search transactions based on quandrant trees.
Conclusions - the research shows the great potential of the new data
structures based on a class of computational geometry data structures and
sweeping curves in spacial distributed dimensions.

###### keywords:

sensor networks, computational geometry, internet of things , Voronoi diagrams
, spacial join , spacial index, sweeping curves, Morton code, Lebesgue curve ,
search operations , distributed database , publisher subscriber

###### MSC:

[2010] 68-06, 68P05

††journal: Future Generation Computer Systems

##  1 Introduction

Historically, noteworthy applications of the Voronoi diagram or Tissen’s
polygons emerged in the writings of Jon Snow, widely regarded as the father of
modern epidemiology, who studied the cholera that plagued London in 1854 At
that time, neither the etiology (note. Cause of occurrence), nor the mode of
transmission of this disease were precisely known, there were disputes about
two possibilities: infection through contact with the patient, his clothes and
/ or property, and the theory that the disease spread through atmospheric
phenomena such as wind. Snow, using a geographical method, revealed that the
cause of the disease was the consumption of water contaminated with feces. To
do this, he mapped the distribution of cholera deaths. Then he studied the
location of drinking water sources in the city and identified the Voronoi
regions for each of them. After calculating the distance between the residence
of each victim and the nearest pumping station, he concluded that the area
most severely affected by cholera corresponds to the Voronoi area associated
with the pumping station on Broad street, as 73 of 83 people died here. After
removing the handle of this pumping station, the outbreak of cholera was
extinguished. The Figure 1 represents that map.

Figure 1: The area most severely affected by cholera

### Systematization of the main problems

IT architecture includes two seemingly incompatible things: on the one hand,
it is a large number of peripheral devices with low computing power, low power
consumption, high speed of reaction to events, with delay in signal
transmission, and on the other hand, cloud servers with high computing power
for processing large amounts of data, storing and classifying them, often with
elements of machine intelligence and analytics. These two worlds use
completely different principles of construction and internal architecture [1].

### Problem identification

There are several levels of integration at which the Internet of Things and
edge technologies pose challenges. The first level is the integration of IoT
and Edge with the underlying IT systems deployed in manufacturing, financial,
engineering and other areas. Many of these basic systems are obsolete. If they
do not have APIs that enable IoT integration, then batch ETL (extract,
transform, load) software may be required to load data into these systems.

The second area of challenge is the IoT itself. Many IoT devices are built
by independent vendors. They use their own proprietary OS. This makes it
difficult to ”mix and match” different IoT devices in a single Edge
architecture. The area of IoT security and compliance is still evolving.
Meanwhile, organizations’ IT professionals can now ask potential IoT vendors
about what is already available for heterogeneous integration and whether they
plan to provide interoperability in future products.

When it comes to integrating legacy systems, ETL is the one integration method
that can help if APIs for the systems are not available. Another alternative
is to write an API, but it takes a long time. The good news is that most
legacy vendors are aware of the upcoming IoT wave and are already developing
their own APIs if they haven’t already. IT departments should check with major
system vendors to find out what their plans for IoT integration are. Once the
data volume grows significantly, it turns out that building the spacial index
itself can be problematic. Either the algorithm goes beyond the definition, or
the build time becomes unacceptably long, or the search is ineffective … As
soon as the spacial data begins to change intensively, the spacial index can
begin to degrade. Rebuilding the index helps, but not for long. Etl layer -
the level of collection, processing and storage of data. The back-end ETL
(extract, transform, and load) is the third ETL operation. The first was in
the peripheral, the second was in the gateway. The back-end ETL collects data
from all peripherals and gateways and is responsible for the following
operations:

  * 1.

Collection of information,

  * 2.

Bringing information to standard form,

  * 3.

Saving information for future use,

  * 4.

Information lifecycle management including archiving and destruction,

  * 5.

Notifying other services when new data arrives.

The Figure 2 represents The typical ETL layer in the IoT sensor network
systems.

Figure 2: The typical ETL layer in the IoT sensor network systems

To organize spacial search in multidimensional indices, this work proposes to
find a logarithmic way of organizing a computational model. In connection with
the problems that arose with R trees [2], quadrant trees [3], this work
settled on the development of algorithms based on self-similar functions [4],
isomorphic to binary trees[5]. Self-similar functions [6] are used to organize
multidimensional indexing and searching.

### Basic mathematical model

Basic mathematical model represents the architechture model for the sensor
network system based on publisher/subscriber architecture with the Voronoi
diagrams data structure [7] data circulation. Storage operations (load) are
intended for storing, sorting and subsequent retrieval of information.
Different tools are used depending on the type of information and its use
cases. If the data does not have a strict schema (table columns), then it is
stored in NoSQL databases. However, if the data can be systematized with a
fixed schema, then SQL database types are used. The latter, in turn, have 2
types - OLTP (Online Transactional Processing) and OLAP (Online Analytic
Processing). As the name suggests, the first type is more suitable for the ETL
process itself - writing new values to the database, while the second is
more convenient for searching and analyzing data. Therefore, often after
loading the OLTP database, in the background, the data is copied to OLAP.
There are situations when it is not convenient or possible to store data in
databases, for example, in the form of a record. This data is written to the
data bucket, and the metadata of records is stored in databases. To reduce
storage costs, obsolete data is archived or deleted. And the last component of
this layer is the internal notification of the presence of new stored data for
presentation to clients and for analysis services [8]. This work proposes the
NoSQL approach for data processing based on Voronoi diagrams representation of
the set of the objects in the system within the construction of the spacial
indexation in NoSQL databases. Every sensor has its own unique identifier ID
correspondenting to his IP address. The appriximate model of the traffic
exchanges from the sensor to the server proposes the following (1):

$\gamma=\frac{\rho\bar{t}}{2(1-\rho)}\frac{\sigma_{a}^{2}+\sigma_{s}^{2}}{\bar{t}^{2}}\frac{\bar{t}^{2}+\sigma_{s}^{2}}{\bar{a}^{2}+\sigma_{s}^{2}}+\bar{t}$
(1)

where $\sigma_{a}^{2},\sigma_{s}^{2}$ \- variance values of time interval
between packets and time service, $\bar{a}$ \- average value of the interval
between packets, $\bar{t}$ \- average service time, $\rho$ \- system’s
loading.

The flow of messages entering the messaging queue from publishers in the
system (2):

$\varrho(x)=\frac{\exp^{-\lambda}\lambda^{x}}{x!}$ (2)

where $\varrho(x)$ \- possibility of obtaining the $x$ income signals in the
time unit, $x$ \- number of requests per unit of time, $\lambda$ \- the
average number of applications per unit of time (rate of receipt of
applications), $exp$ = 2.7182 - base of the natural logarithm.

##  2 Methodology

The main specificity of the Lebesgue’s curve [9] is localized in the splitting
of the subquery.

  * 1.

first we find the starting extent, which is the minimum rectangle that
includes the search extent and contains one contiguous range of key values,

  * 2.

calculate the key values,

  * 3.

for the lower left and upper right points of the search extent (KMin, KMax),

  * 4.

find a common binary prefix (from high order to low order) for KMin, KMax,

  * 5.

zeroing all the digits behind the prefix, we get SMin, filling them with ones,
we get SMax,

  * 6.

do the inverse transformation of keys to coordinates and find the corners of
the starting extent.

In the case of the Hilbert curve, by the way, the lower left corner of the
starting extent does not necessarily come from SMin, one needs to choose the
minimum value. The same with the upper right corner. The starting extent can
be noticeably larger than the search extent, if one is not unlucky, it will
turn out to be the maximum extent of the layer (empty prefix). For the
Z-curve, optimization can be done and the starting extent can be equated to
the search extent. This is possible due to the peculiarity of the z-curve -
for any rectangle, its lower left corner gives the minimum key value, and the
upper right corner gives the maximum (in the rectangle). Moreover, such a
starting extent can contain more than one range of values, but this is further
removed by filtering. If one pushes the starting extent onto the subquery
stack until the subquery stack is empty one will get the top element, if it
does not overlap with the search extent, one should discard it and skip
iteration. This is in the case the starting extent is larger than the search
extent and one needs to ignore the excess. Furthermore, if the minimum point
of the subquery is within the search extent, one queries the index by the
value of this point. As a result, one gets two values - ”the first key” higher
or equals the desired one and the ”last key” on the physical (leaf) page,
where the ”first key” lies:

  * 1.

if “first key” lower than the maximum value of the current subquery, one
should ignore the current subquery, then skip iteration,

  * 2.

if the entire subquery is inside the search extent, one subtracts it by the
traverse of the index from the minimum to the maximum value, then end of
iteration,

  * 3.

if the “last key” higher than the maximum value of the current subquery, then
all the data of the current subquery is on this page, one needs to read the
page to the end and filter out the unnecessary, then end of iteration,

  * 4.

the separate case ”last key” equials to the maximum value of the current
subquery, then one processes separately and traverse forward splitting the
current subquery,

  * 5.

the one adds 0 and 1 to its prefix - to get two new prefixes,

  * 6.

then one fills in the remainder of the key 0 or 1 - to get the minimum and
maximum values of the new subqueries,

  * 7.

one pushes them onto the stack, first the one that added 1, then 0 (this is
for unidirectional reading of the index).

The Figure 3 represents the proposed transactional model for the ETL level.

Figure 3: The proposed transactional model for the ETL level

### Results

The results represents the first step realizations of the developed
transactional ETL model based on Voronoi Diagrams and Lebesgue curves and
consists of the realization of the Voronoi diagrams using the C++, Ubuntu OS
and OpenGL in 2D dimension in two performances: static version and dynamically
recalculating version in time. The Figure 4 represents the static realization
of the Voronoi diagrams.

Figure 4: The static realization of the Voronoi diagrams in 2D

The Figures 5 and 6 represents the dynamic realization of the Voronoi
diagrams.

Figure 5: The dynamic realization of the Voronoi diagrams in 2D Figure 6: The
static realization of the Voronoi diagrams in 2D

Moreover, the results represents the transactional model based on the
publisher/subscriber architecture. The Figure 3 represents the proposed
transactional model for the ETL level. The static Voronoi Diagrams situates in
this repository "the static Voronoi Diagrams realization". The dynamic Voronoi
Diagrams situates in this repository "the dynamic Voronoi Diagrams
realization".

##  3 Research

The research potential of the work lies in the possibility of using algorithms
for building maps with the need to plot the density of parameters. The
potential lies in the displaying the density maps for the necessary parameters
(for example, the dynamics of animals migration, the assessment of the
geological resources in different regions, the dynamics of the extinction of
wild animal species, the dynamics of the disappearance of forests, the traffic
congestion, the consumption of housing and communal resources, the sensitivity
of crystals on the surface, the chemical properties of materials, the physical
properties of materials, etc…) in loci / zones of the same density.

The research potential of the work lies in the possibility of applying
algorithms for constructing multidimensional spacial indexes to optimize
search queries in a growing volume of data.

### Applications

The applications of the proposed model can be realized in multivarious
applications in the sensor network systems:

  * 1.

for search operations in the database optimization of the multidimentional
spacial geo systems,

  * 2.

for search operations in the distibuted transactions such as blockchain,

  * 3.

for identification of the closest neighbor objects to the object.

### Implication examples

The class of the Lebesgue curves (Morton code)[10], [11] can be found in the
following database management systems:

  * 1.

AWS database Amazon Aurora [12],

  * 2.

Amazon RDS for MySQL [12].

##  4 Discussion

Due to the recursive structure of the Voronoi diagrams and Lebesgue curves
there exists a possibility to use the recursive defragmentation of the
iterations in the both algorithms. Due to the spacial adaptiveness of the
Voronoi diagrams this work proposes to explore the properties of the
adaptiveness and time costs of these algorithms in real-time systems. The
estimated time for Voronoi diagrams calculation is $N\log{N}$, where $N$ \- is
the number of iterations. The estimated time of the Lebesque curve indexation
also is $N\log{N}$, where $N$ \- is the number of iterations.

## References

## References

  * [1] D. Nimratz, Iot architecture - a first look under the hood, Electronic resource (Aug. 2018).
URL https://habr.com/ru/post/420173/

  * [2] A. Guttman, R-trees: A dynamic index structure for spatial searching, SIGMOD Rec. 14 (2) (1984) 47–57.  doi:10.1145/971697.602266.
URL https://doi.org/10.1145/971697.602266

  * [3] B. Muratshin, Hilbert, lebesgue … and emptiness, Electronic resource (Aug. 2019).
URL https://habr.com/ru/users/zzeng/

  * [4] B. Muratshin, Hilbert curve vs z-order, Electronic resource (Oct. 2017).
URL https://habr.com/ru/post/340100/

  * [5] B. Muratshin, About z-order and r-tree, Electronic resource (Jan. 2017).
URL https://habr.com/ru/post/319096/

  * [6] H. C. J. Moon, BongkiJagadish, Analysis of the clustering properties of hilbert space-filling curve, IEEE Transactions on Knowledge and Data Engineering.
URL https://drum.lib.umd.edu/handle/1903/804

  * [7] L. Guibas, J. Stolfi, Primitives for the manipulation of general subdivisions and the computation of voronoi, ACM Trans. Graph. 4 (2) (1985) 74–123.  doi:10.1145/282918.282923.
URL https://doi.org/10.1145/282918.282923

  * [8] D. Nimratz, Iot architecture, Electronic resource (Jun. 2019).
URL https://habr.com/ru/post/455377/

  * [9] A. Esculier, Lebesgue curve, lebesguesche kurve, Electronic resource (2006).
URL https://mathcurve.com/fractals/lebesgue/lebesgue.shtml

  * [10] G. M. Morton, A computer Oriented Geodetic Data Base; and a New Technique in File Sequencing, Technical report, IBM Ltd., 1966.
  * [11] Wikipedia, Z-order curve, Electronic public library article (Jan. 2022).
URL https://en.wikipedia.org/wiki/Z-order_curve

  * [12] S. Chandrasekaran, Amazon aurora under the hood: indexing geospatial data using z-order curves, AWS Database Blog (Jan. 2018).
URL https://aws.amazon.com/blogs/database/amazon-aurora-under-the-hood-
indexing-geospatial-data-using-z-order-curves/


# First-principles molecular dynamics simulation of liquid indium

Yu. D. Fomin 111Corresponding author<EMAIL_ADDRESS>Vereshchagin Institute
of High Pressure Physics, Russian Academy of Sciences, Kaluzhskoe shosse, 14,
Troitsk, Moscow, 108840 Russia  Moscow Institute of Physics and Technology
(National Research University), 9 Institutskiy Lane, Dolgoprudny, Moscow
region, 141701, Russia

###### Abstract

We report an ab-initio simulation of liquid Indium in a wide range of
pressures and temperatures. We calculate equation of state, thermal expansion
and compressibility coefficients. The structure of the system is analyzed by
radial distribution functions and structure factors. The results are compared
with available experimental data.

###### pacs:

61.20.Gy, 61.20.Ne, 64.60.Kw

##  I Introduction

Liquid indium is an important metal in many fields of industry. It can be used
as coverege to reduce friction coefficient, as a solder, as an additive to
glasses to modify their optical properties, etc. It also attracted attention
of researchers for many years. Since 196-th there were several important works
which studied the structure of liquid indium at different temperatures str-1 ;
str-2 ; str-3 . Most of them, however, considered the system at ambient
pressure only. At the same time the melting curve of indium is measured up to
rather high pressure of 10.5 GPa dudley  (see also mc  for the melting curve
up to about 1 GPa).

Basing on the results of Refs. dudley  and mc  the authors of Ref. shen
studied the structure and the density of liquid indium up to the melting line
at $T=710$ K. Basing on X-ray techniques the authors measured the density of
molten indium at several pressures up to the melting point and the structure
factors. They also determined the coordination number by integrating of radial
distribution functions $g(r)$. From their results one can see that the
coordination number monotonically increases with pressure.

Another X-ray study of liquid indium is reported in Ref. mudry . In this study
experiments were performed at ambient pressure and a set of temperatures. The
authors report that the coordination nember behaves non-monotonically and
conclude that a liquid-liquid phase transition takes place.

In Ref. shen1  the molar volume of liquid indium at $T=710$ K was measured up
to $P=8.5$ GPa. Later on these data were fitted to Birch-Murnaghan equation in
Ref. liq-in . Basing on these fits the authors calculated the bulk modulus and
thermal expansion coefficient of liquid indium along the $T=710$ K isotherm.

In the present work we perform ab-initio molecular dynamics calculations of
liquid indium. We calculate the equation of state of the system, its bulk
modulus and thermal expansion coefficient. We also compute the structure
factors. All quantities are compared with the experimental data where it is
possible.

##  II System and methods

In the present study we investigate a system of 100 atoms of indium in a cubic
box with periodic boundary conditions. The system is simulated by means of ab-
inition molecular dynamics method as realized in VASP simulation package.
Projector augmented wave method is used for the treatment of the electronic
structure. The energy cut-off is set to 96 eV. The time step is set to 3 fs.
Only gamma point was taken into account in the ${\bf k}$ space calculations.

Large set of densities and temperatures is considered. The densities varied
from $\rho_{min}=3.88$ to $\rho_{max}=19.06$ $g/cm^{3}$. At each density we
firstly simulated the high temperature system with $T_{max}=1000$ K. Then the
last configuration of this simulation was used as an initial one for
simulation of the system at the temperatures from $T=300$ up to $T=900$ K with
the step of $100$ K. Each simulation lasted for 10000 steps, i.e. 30 ps. The
set of simulated points is shown on the phase diagram in Fig. 1

Figure 1:  Points where simulations were performed on the phase diagram. The
curve ’Dudley et. al.’ shows the results from Ref. dudley , the curve
’McDaniel’ - results from Ref. mc , ’Dudley et. al. S’ - fitting the results
of Ref. dudley  to Simon equation (Eq. (2) of Ref. dudley ), TW - the points
of the present work.

##  III Results and discussion

At the first step we calculate the equation of state on liquid indium and
compare it to the experimental data from ref. liq-in .

Fig. 2 presents a comparison of the data of this work and the experimental
results from liq-in . The experimental data are taken at $710$ K, while in the
present work $T=700$ K is used. We believe that this difference sufficiently
small to not affect the comparison.

The difference between the experimental and simulation data does not exceed
$1.12\%$. Basing on this we conclude that the results of simulations are
reliable.

Figure 2:  Comparison of experimental data at $T=710$ K from Ref. liq-in  and
the data from ab-initio simulation of this work at $T=700$ K.

In order to evaluate the equation of state of liquid indium on more solid
basis we fit it to a polynomial function

$P(\rho,T)=\sum_{n,m}a_{n,m}\rho^{n}T^{m},$ (1)

where $\rho$ is density in $g/cm^{3}$, T is temperature in Kelvins, $a_{n,m}$
are fitting coefficients and the exponents n and m are choosen in such a way
that $n+m<=5$.

Fig. 3 shows a comparison of the data for equation of state from ab-inition
simulation and from fitting to Eq. (1). One can see that the quality of the
fitting is perfect. Below we use this fitted equation of state for calculation
of the response functions $B_{T}=\rho\left(\frac{\partial
P}{\partial\rho}\right)_{T}$ and
$\alpha_{P}=-\frac{1}{T}\left(\frac{\partial\rho}{\partial T}\right)_{T}$.
Similarly, fitting of the energy is used for calculation of the heat capacity
of the system.

n | m | $a_{n,m}$
---|---|---
0 | 0 | -305.57333657145500
1 | 0 | 238.70598753111230
2 | 0 | -63.833351593789722
3 | 0 | 6.5433130551773946
4 | 0 | -0.24008639734032150
5 | 0 | 3.9155839840238878E-003
0 | 1 | -0.17373675573617220
1 | 1 | 1.4237276981756433E-002
2 | 1 | 8.0437749691983397E-003
3 | 1 | -5.7392554358567341E-004
4 | 1 | 1.5599905130353485E-005
2 | 2 | 5.8408548738952959E-004
1 | 2 | -1.3516570420262682E-004
2 | 2 | -1.5393806829344470E-006
3 | 2 | -3.4804293612014242E-008
0 | 3 | -5.1700880554506057E-007
1 | 3 | 1.7197816344291688E-007
2 | 3 | 1.3138966034853920E-009
0 | 4 | -6.8034449601794833E-011
1 | 4 | -7.7188558042335661E-011
0 | 5 | 1.8898204983398251E-013
Table 1: Coefficients of Eq. (1). The resulting pressure is in kbars (0.1
GPa). The dimensionality of coefficients $a_{n,m}$ are choosen in such a way
that $a_{n,m}\rho^{n}T^{m}$ is in kbars.  Figure 3:  Comparison of the data
from ab-inition molecular dymanics and fitting to Eq. LABEL:fit-pr. The data
for $T=1000$ K are shifted to 100 GPa upwards in order to avoid overlaps.

Fig. 4 (a) and (b) show a comparison of the data of the present calculations
with the experimental data from Ref. liq-in  for the bulk modulus and thermal
expansion coefficient. In the case of bulk modulus the curves are very close.
However, the slope of the calculated curve is slightly larger. The largest
descreepancy takes place at low pressure and reaches about $20\%$.

The agreement between simulation and experement becomes worse in the case of
thermal expansion coefficient. Experemental values are smaller then the
calculated ones. The descreepancy increases with pressure reaching about
$40\%$ at $P=8$ GPa.

Figure 4: (a) Comparison of bulk modulus $B_{T}$ at $T=700$ K from ab-initio
simulation of the present work (AMD) and from experiments of Ref. liq-in . (b)
The same for the thermal expansion coefficient. (c) The same for isochoric
heat capacity $c_{V}$ (in units of $k_{B}$). The straight line shows
$c_{V}=2k_{B}$ level.

In order to characterize the structure of the liquid we calculated the radial
distribution functions $g(r)$. Fig. 5 shows $g(r)$ for a set of pressures
along the isotherm $T=700$ K. According to Ref. liq-in  the melting point of
indium at $T=710$ K is $P_{m}=7.77$ GPa. From Fig. 5 (a) one can see that even
at the pressure as high as $P=26.70$ GPa the system preserves liquid-like
structure which should be related to formation of metastable liquid in ab-
initio simulations. Some signs of crystallinity can be observed at pressures
as high as $P=50.75$ GPa whish is much higher then experimental melting
pressure.

Figure 5: (a) Radial distribution functions of liquid indium at $T=700$ K
isotherm for the pressures from $P=-0.96$ GPa up to 26.60 GPa. (b) The same
for higher pressures.

X-ray experiments measure structure factor $S({\bf k})$ which is a Fourier
transform of radial distribution function. Fig. 6 shows the structure factor
of the system obtained in the present work in comparison to the experimental
ones from Refs. shen ; mudry . The structure factors of this work are obtained
by numerical Fourier transform of $g(r)$. The results for $k<1.6$ $\AA^{-1}$
should be addressed to numerical errors and should be discarded from
consideration.

One can see that experimental structure factors demonstrate very slow decay of
the right branch of the first peak. Unlike this the first peak of $S({\bf k})$
of our calculations look more symmetric and ’Gauss-like’. However, the overall
aggreement is sufficiently good.

Figure 6: Comparison of the structure factors obtained in the present work
with the experimental ones. The curves are shifted upward to avoid overlaps.
Pressures are given next to the curves. The curve at $P=0$ GPa and $T=690$ K
is from Ref. mudry . The curves at $P=1.0$, $2.0$, $3.1$, $5.1$ and $6.3$ GPa
are from Ref. shen . The curves with label ’MD’ after the pressure are our
calculations.

Although the crystallization is not visible in radial distribution functions,
it is more pronounced in the dynamical properties of the system characterized
by means square displacement. Fig. 7 (a) shows mean square displacement long
the $T=700$ K isotherm. One can see that some mobility of atoms preserves up
to the pressure $=12.28$ GPa and disappearsh at higher pressures. Fig. 7 (b)
demonstrates the diffusion coefficient along the same isotherm. One can see
that it monotonously decreases with pressure and vanishes at $P=16.82$ GPa.

Figure 7: (a) Mean square displacement of liquid indium at a set of pressures
at $T=700$ K. (b) Diffusion coefficient of liquid indium along the same
isotherm.

This work was carried out using computing resources of the federal collective
usage center ”Complex for simulation and data processing for mega-science
facilities” at NRC ”Kurchatov Institute”, http://ckp.nrcki.ru, and
supercomputers at Joint Supercomputer Center of the Russian Academy of
Sciences (JSCC RAS). The work was supported by the Council of the President of
the Russian Federation for State Support of Young Scientist (Grant
MD-6103.2021.1.2).

## References

  * (1) H. Ocken and C. N. J. Wagner, Temperature Dependence of the Structure of Liquid Indium, Phys. Rev. 149, 122 (1966)
  * (2) H. Ruppersberg and K. H. Winterberg, Structure factor and resistivity of liquid indium at temperature between 165 and 665 o C
  * (3) M. L. Harris, R. L. Collier, R. W. Gruebel and T. O. Callaway, The direct pair correlation function of liquid indium, Physics Letters 65A, 244 (1978).
  * (4) J. D. Dudleyj and H. T. Hall, Experimental Fusion Curves of Indium and Tin to 105 000 Atmospheres, Phys. Rev. 118, 1211 (1960)
  * (5) M. L. McDaniel, S. E. Babb Jr., and G. J. Scott, Melting Curves of Five Metals under High Pressure, J. Chem. Phys. 37, 822 (1962)
  * (6) G. Shen, N. Sata, N. Taberlet, M. Newville, M. L. Rivers, and St. R. Sutton, Melting studies of indium: determination of the structure and density of melts at high pressures and high temperature, J. Phys.: Condens. Matter 14, 10533 10540 (2002)
  * (7) S. Mudry, I. Shtablavyi, and U. Liudkevych, The relation between structure changes and thermal expansion in liquid indium, Phys. and Chem. of Liquids 55, 254-263 (2017)
  * (8) G. Shen, N. Sata, M. Newville, M. L. Rivers, and S. R. Sutton, Molar volumes of molten indium at high pressures measured in a diamond anvil cell, Applied Physics Letters 81, 1411 (2002)
  * (9) H. Li, Yo. Sun, and M. Li, Equation of state of liquid Indium under high pressure, AIP Advances 5, 097163 (2015).