# Complex Scalar fields in Scalar-Tensor and Scalar-Torsion theories

Andronikos Paliathanasis<EMAIL_ADDRESS>Institute of Systems Science,
Durban University of Technology, PO Box 1334, Durban 4000, South Africa
Instituto de Ciencias Físicas y Matemáticas, Universidad Austral de Chile,
Valdivia 5090000, Chile  Mathematical Physics and Computational Statistics
Research Laboratory, Department of Environment, Ionian University, Zakinthos
29100, Greece

###### Abstract

We investigate the cosmological dynamics in a spatially flat
Friedmann–Lemaître–Robertson–Walker geometry in scalar-tensor and scalar-
torsion theories where the nonminimally coupled scalar field is a complex
field. We derive the cosmological field equations and we make use of
dimensionless variables in order to determine the stationary points and
determine their stability properties. The physical properties of the
stationary points are discussed while we find that the two-different theories,
scalar-tensor and scalar-torsion theories, share many common features in terms
of the evolution of the physical variables in the background space.

Scalar field; Complex field; Scalar-tensor; Scalar-torsion; Dynamical
analysis.

###### pacs:

98.80.-k, 95.35.+d, 95.36.+x

##  I Introduction

Scalar fields in gravitational theory is a simple mechanism to introduce new
degrees of freedom which can play an important role in the description of
observable cosmological phenomena ratra . The early acceleration phase of the
universe is attributed to a scalar field known as inflaton in0 ; in1 ; in2 ;
in3 . In addition, scalar fields have been used to describe also the late-time
acceleration phase of the universe attributed to dark energy, or other matter
components such as dark matter, see for instance sf1 ; sf2 ; sf3 ; sf4 ; sf5 ;
sf6 ; sf7 ; sf8 ; sf9 ; sf10  and references therein.

Multi-scalar fields have been widely studied in the literature. Some well-
known two-scalar field models are the quintom qq1  or the Chiral model which
leads to hyperbolic inflation hy4 , while other proposed multi-scalar fields
theories can be found for instance in cc1 ; mm2 ; mm3 ; mm4 ; mm5 ; mm6  and
references therein. A simple mechanism to introduce a multi-scalar field
theory is to consider the existence of a complex scalar field, the real and
imaginary parts of which give the equivalent of a two scalar-field theory com1
; com2 .

An inflationary model with a complex scalar field was proposed in com3 .
Specifically it was found that, when inflation occurs, the imaginary component
of the complex scalar field does not contribute in the cosmological fluid,
that is, the phase of the complex scalar field is constant. The cosmological
perturbations with a complex scalar field were investigated in com4 .
Furthermore, in com5  the authors used cosmological observations to
reconstruct the quintessence potential for a complex scalar field. A
nonminimally coupled scalar field cosmological model has been studied in com6
, while for some recent studies of complex scalar field cosmological models we
refer the reader to com7 ; com8 ; com9 ; com10 ; com11 ; com12 ; com13 ; com14
.

In this study we consider the scalar-tensor and the scalar-torsion theories
with a complex scalar field sc1 . In scalar-tensor theory, the scalar field is
minimally coupled to gravity. The scalar field interacts with the
gravitational Action Integral of Einstein’s General Relativity, that is, the
Ricci-scalar of the Levi-Civita connection, the scalar-tensor theory satisfies
the Machian Principle Brans  and the theory is defined in the so-called Jordan
frame Jord . The Brans-Dicke Lagrangian Brans  is the most common scalar-
tensor theory. On the other hand, the scalar-torsion theory is the equivalent
of scalar-tensor model in teleparallelism. In the latter, the fundamental
geometric invariant is the torsion scalar determined by the curvatureless
Weitzenböck connection Hayashi79 ; Tsamp . There are many important results in
the literature on the cosmological studies of the scalar-tensor and scalar-
torsion theories, for an extended discussion we refer the reader to omegaBDGR
; st1 ; st2 ; st3 ; st4 ; st5 ; st6 ; st7 ; st8 ; te2 ; te3 ; te4  and
references therein.

The purpose of this study is to investigate the effects of a complex scalar
field in the evolution of the cosmological dynamics for the background space
for these two different gravitational theories and to compare the results.
Such analysis provides us with important results in order to understand the
differences between the use of the Ricci-scalar and of the torsion scalar in
the background geometry. We use dimensionless variables to perform a detailed
analysis of the dynamical systems which describe the evolution of the physical
variables. Such an approach has been widely studied before with many
interesting results about the viability of proposed gravitational theories,
see for instance dyn1 ; dyn2 ; dyn3 ; dyn4 ; dyn6 ; dyn5 . The plan of the
paper is as follows.

In Section II we consider a spatially flat Friedmann–Lemaître–Robertson–Walker
geometry in scalar-tensor gravitational theory with complex scalar field. We
derive the field equations and we write the point-like Lagrangian. Moreover,
we determined the stationary points for the field equations and we investigate
their stability properties. In Section III we perform a similar analysis but
now in the context of scalar-tensor theory. Finally, in Section IV we
summarize our results and compare the physical results provided by the two
theories and draw our conclusions.

##  II Scalar-Tensor Cosmology

Consider the complex scalar field $\psi$ in the case of scalar-tensor theory
sc1 , for which the gravitational Action Integral is defined as

$S_{Tensor}=\int
d^{4}x\sqrt{-g}\left[F\left(\left|\psi\right|\right)R+\frac{1}{2}g^{\mu\nu}\psi_{,\mu}\psi_{,\nu}^{\ast}-V\left(\left|\psi\right|\right)\right],$
(1)

where $R$ is the Ricci scalar related to the Levi-Civita connection for the
metric tensor $g_{\mu\nu}$; $\left|\psi\right|$ is the norm of the complex
field $\psi$, that is, $\left|\psi\right|^{2}=\psi\psi^{\ast}$,
$F(\left|\psi\right|)$ is the coupling function between the gravitational and
the scalar field $\psi$ and $V\left(\left|\psi\right|\right)$ is the potential
function which drives the dynamics. The Action Integral (1) admits the
$U\left(1\right)$ symmetry.

The Brans-Dicke theory with a complex field is recovered for
$F(\left|\psi\right|)=F_{0}\left|\psi\right|^{2}$, that is, the Action
Integral (1) is sc1

$S_{BD}=\int
d^{4}x\sqrt{-g}\left[F_{0}\left|\psi\right|^{2}R+\frac{1}{2}g^{\mu\nu}\psi_{,\mu}\psi_{,\nu}^{\ast}-V\left(\left|\psi\right|\right)\right].$
(2)

In the case of a spatially flat FLRW universe with line element,

$ds^{2}=-N^{2}\left(t\right)dt^{2}+a^{2}\left(t\right)\left(dx^{2}+dy^{2}+dz^{2}\right)~{},$
(3)

we derive for the Ricciscalar

$R=6\left(\frac{1}{N}\dot{H}+12H^{2}\right)~{},~{}$ (4)

where $H=\frac{1}{N}\frac{\dot{a}}{a}$, $\dot{a}=\frac{da}{dt}$, is the Hubble
function.

We substitute (4) into (1) which by integration by parts gives the point-like
Lagrangian

$L_{Tensor}\left(N,a,\dot{a},\psi,\dot{\psi}\right)=\frac{1}{N}\left(6F\left(\left|\psi\right|\right)a\dot{a}^{2}+6\dot{F}\left(\left|\psi\right|\right)a^{2}\dot{a}+\frac{1}{2}a^{3}\left(\dot{\psi}\dot{\psi}^{\ast}\right)\right)-a^{3}NV\left(\left|\psi\right|\right),$
(5)

or in the case of Brans-Dicke

$L_{BD}\left(N,a,\dot{a},\psi,\dot{\psi}\right)=\frac{1}{N}\left(6F_{0}\left|\psi\right|^{2}a\dot{a}^{2}+6F_{0}\left(\left|\psi\right|^{2}\right)^{\cdot}a^{2}\dot{a}+\frac{1}{2}a^{3}\left(\dot{\psi}\dot{\psi}^{\ast}\right)\right)-a^{3}NV\left(\left|\psi\right|\right).$
(6)

We focus now in the Brans-Dicke theory in which the field equations are
described by the point-like Lagrangian (6). Moreover, the complex scalar field
is written with the use of the polar form,
$\psi\left(t\right)=\phi\left(t\right)e^{i\theta\left(t\right)}$, such that
the Lagrangian (6) becomes

$L_{BD}\left(N,a,\dot{a},\phi,\dot{\phi},\theta,\dot{\theta}\right)=\frac{1}{N}\left(6F_{0}\phi^{2}a\dot{a}^{2}+12F\phi
a^{2}\dot{a}\dot{\phi}+\frac{1}{2}a^{3}\left(\dot{\phi}^{2}+\phi^{2}\dot{\theta}^{2}\right)\right)-a^{3}NV\left(\phi\right).$
(7)

It is obvious that the Lagrangian function (7) describes a multi-scalar field
cosmological model, where $\phi$ is the Brans-Dicke field and $\theta$ is a
second-scalar field minimally coupled to gravity but coupled to the Brans-
Dicke field $\phi$. The $U\left(1\right)$ for the point-like Lagrangian (7)
provides the invariant transformation $\theta=\theta+\varepsilon$, and the
conservation law $I_{0}=\frac{1}{N}a^{3}\phi^{2}\dot{\theta}$.

Variation with respect to the dynamical variables
$\left\\{N,a,\phi,\theta\right\\}$ of the Lagrangian (7) provides the
cosmological field equations which are

$0=6F_{0}\phi^{2}H^{2}+12F\phi
H\dot{\phi}+\frac{1}{2}\left(\dot{\phi}^{2}+\phi^{2}\dot{\theta}^{2}\right)+V\left(\phi\right)~{},$
(8)
$0=2F_{0}\phi^{2}\left(2\dot{H}+3H^{2}\right)+8F_{0}H\phi\dot{\phi}-\frac{1}{2}\dot{\phi}^{2}+4F_{0}\dot{\phi}^{2}-\frac{1}{2}\dot{\phi}^{2}\dot{\theta}^{2}+4F_{0}\phi\ddot{\phi}+V\left(\phi\right)~{},$
(9)
$0=\ddot{\phi}+\phi\left(12F_{0}\dot{H}-\dot{\theta}^{2}\right)+3\dot{\phi}\dot{H}+12F_{0}\phi
H^{2}+V_{,\phi}~{}$ (10)

and

$0=\phi\ddot{\theta}+\left(2\dot{\phi}+3H\phi\right)\dot{\theta}~{}\text{or
}I_{0}=a^{3}\phi^{2}\dot{\theta}\text{,}$ (11)

where without loss of generality we have assumed $N\left(t\right)=1$.

###  II.1 Cosmological dynamics

In order to reconstruct the cosmological history provided by this specific
complex scalar-tensor theory we make use of dimensionless variables in the
context of $H$-normalization and we investigate the dynamical evolution of the
field equations (8)-(11) by determining the stationary points and their
stability properties.

We consider the new dimensionless dependent variables dyn1

$\dot{\phi}=2\sqrt{3}H\phi
x~{},~{}V\left(\phi\right)=6H^{2}\phi^{2}y~{},~{}\dot{\psi}=2\sqrt{3}Hz~{},~{}V_{,\phi}=\lambda\frac{V\left(\phi\right)}{\phi}.$
(12)

Hence, for the new independent variable, $\tau=\ln a$, the field equations
(8)-(11) are written as the following dynamical system

$\displaystyle 2F_{0}\left(12F_{0}-1\right)\frac{dx}{d\tau}$
$\displaystyle=2\sqrt{3}\left(7-48F_{0}\right)x^{2}+\left(3-24F_{0}\right)x^{3}$
(13)
$\displaystyle+2\sqrt{3}F_{0}\left(F_{0}+z^{2}+\left(\lambda-3\right)y\right)+3x\left(F_{0}+\left(4\lambda
F_{0}-1\right)y+\left(1-8F_{0}\right)z^{2}\right)~{},$ $\displaystyle
F_{0}\left(12F_{0}-1\right)\frac{1}{y}\frac{dy}{d\tau}$
$\displaystyle=2\sqrt{3}F_{0}x\left(4-\lambda+12F_{0}\left(\lambda-2\right)\right)+\left(3-24F_{0}\right)x^{2}$
(14) $\displaystyle+3\left(F_{0}\left(16F_{0}-1\right)+\left(4\lambda
F_{0}-1\right)y+\left(1-8F_{0}\right)z^{2}\right)~{},$
$\displaystyle-2F_{0}\left(12F_{0}-1\right)\frac{1}{3z}\frac{dz}{d\tau}$
$\displaystyle=4\sqrt{3}F_{0}\left(8F_{0}-1\right)x+\left(8F_{0}-1\right)x^{2}$
(15) $\displaystyle+\left(1-4\lambda
F_{0}\right)y+\left(8F_{0}-1\right)\left(F_{0}+z^{2}\right)~{},$
$\frac{d\lambda}{d\tau}=2\sqrt{3}\lambda
x\left(1-\lambda+\Gamma\left(\lambda\right)\right)~{},~{}\Gamma\left(\lambda\right)=\phi\frac{V_{,\phi\phi}}{V_{,\phi}}$
(16)

with algebraic constraint

$x^{2}+z^{2}+y+F_{0}\left(1+4\sqrt{3}x\right)=0~{}.$ (17)

For the scalar field potential we assume
$V\left(\phi\right)=V_{0}\phi^{\lambda_{0}}$; such that $\lambda=\lambda_{0}$
is always a constant and the dimension of the dynamical system is two, after
we apply the constraint equation (17). From (17) we substitute the parameter
$y$ and we end up with a two-dimensional system on the space of variables
$\left\\{x,z\right\\}$. Moreover, we observe that the dynamical system remains
invariant under the change of variable $z\rightarrow-z$. Thus, without loss of
generality, we can restrict our analysis to the region with $z\geq 0$.
Furthermore, for completeness of our analysis the scalar field potential can
be positive or negative such that $y\in\mathbb{R}$.

The stationary points of the reduced dynamical system
$A=\left(x\left(A\right),z\left(A\right)\right)$ are as follows

$A_{1}=\left(-\frac{\sqrt{3}F_{0}\left(\lambda-4\right)}{6F_{0}\left(\lambda+2\right)-3},0\right)~{}$

and

$A_{2}=\left(x_{2},\sqrt{-F_{0}-x_{2}\left(4\sqrt{3}F_{0}+x_{2}\right)}\right)\text{\thinspace}.$

Each stationary point describes an asymptotic solution for the field equations
in which the effective fluid has the equation of state parameter

$w_{eff}\left(x,y,z\right)=\frac{3\sqrt{3}F_{0}x+\left(3-34F_{0}\right)x^{2}+3\left(4F_{0}^{2}+\left(4F_{0}\lambda-1\right)y+\left(1-8F_{0}\right)z^{2}\right)}{3F_{0}\left(12F_{0}-1\right)}.$
(18)

We proceed with the discussion of the physical properties and the stability
properties for the admitted stationary points.

The stationary point $A_{1}$ describes a universe dominated only by the scalar
field $\phi$ and its potential function $V\left(\phi\right)$, because
$y\left(A_{1}\right)=\frac{F_{0}\left(12F_{0}-1\right)\left(3+F_{0}\left(\lambda-10\right)\left(2+\lambda\right)\right)}{3\left(1-2F_{0}\left(2+\lambda\right)\right)^{2}}$.
The effective equation of state parameter for the asymptotic solution is
$w_{eff}\left(A_{1}\right)=\frac{3+2F_{0}\left(2+\lambda\left(\lambda-9\right)\right)}{6F_{0}\left(2+\lambda\right)-3}$,
from which it follows that $A_{1}$ describes an accelerated universe as is
given in Fig. 1. Moreover, the eigenvalues of the linearized system around the
stationary point are derived,
$e_{1}\left(A_{1}\right)=\frac{3+F_{0}\left(\lambda-10\right)\left(2+\lambda\right)}{2F_{0}\left(2+\lambda\right)-1}~{}$,$~{}e_{2}\left(A_{1}\right)=\frac{3+F_{0}\left(\lambda-10\right)\left(2+\lambda\right)}{2F_{0}\left(2+\lambda\right)-1}$.In
Fig. 1 we present the region plot in the two-dimensional space of the free
variables $\left\\{F_{0},\lambda\right\\}$, where the eigenvalues have
negative real parts, that is, point $A_{1}$ is an attractor.

Figure 1: Region in the space of the free variables
$\left\\{F_{0},\lambda\right\\}$, where
$w_{eff}\left(A_{1}\right)<-\frac{1}{3}$ (Left Fig.) and where $A_{1}$ is an
attractor (Right Fig.). It is easy to observe that in this specific region
when $A_{1}$ describes an accelerated universe is also an attractor for the
dynamical system.

Points $A_{2}$ describe a family of stationary points which are real and
physical accepted when $-F_{0}-x_{2}\left(4\sqrt{3}F_{0}+x_{2}\right)>0$. For
the asymptotic solution we calculate
$w_{eff}\left(A_{2}\right)=1+\frac{8}{\sqrt{3}}x_{2}$, which means that the
family of points $A_{2}$ describe accelerated universes for
$x_{2}<-\frac{1}{2\sqrt{3}}$. Furthermore, the eigenvalues are calculated
$e_{1}\left(A_{2}\right)=6+2\sqrt{3}\left(2+\lambda\right)x_{2}$ and
$e_{2}\left(A_{2}\right)=0$. Then for
$6+2\sqrt{3}\left(2+\lambda\right)x_{2}<0$, where $e_{1}\left(A_{2}\right)<0$,
the Center Manifold Theorem (CMT) can be applied in order to investigate for a
possible stable submanifold and to infer about the stability. However, such an
analysis does not contribute to the physical discussion of the present theory
and we select to work numerically. In Fig. 2 we present the phase space
portrait in the two-dimensional place $\left(x,z\right)$ from which we observe
that the stationary points are always saddle points for
$e_{1}\left(A_{2}\right)<0$.

From the above analysis it is clear that various major eras of the
cosmological evolution can be described by the stationary points, $A_{1}$ and
$A_{2}$. For instance for
$F_{0}=-\frac{3}{2}\left(2-9\lambda+\lambda^{2}\right)^{-1}$,
$w_{eff}\left(A_{1}\right)=0$ which means that $A_{1}$ describes the matter
dominated era, and points $A_{2}$ can describe the early and late acceleration
phases of the universe. That is not the unique case, since, for
$x_{2}=-\frac{\sqrt{3}}{8},$ $w_{eff}\left(A_{2}\right)=0,$ and $A_{1}$ can
describe the future acceleration phase of the universe.

Figure 2: Phase-space portrait for the scalar-tensor theory on the two-
dimensional plane $\left(x,z\right)$ for different values of the free
parameters $F_{0}$ and $\lambda$. From the figure it is clear that the surface
of points described by $A_{2}$ are always saddle points or source points. With
solid line we present the family of points $A_{2}$ for
$e_{1}\left(A_{2}\right)<0$.

###  II.2 Analysis at Infinity

Because variables $\left(x,z\right)$ are not bounded, they can take values at
all the range of the real numbers, which means that they can take values at
infinity. Until now we have investigated the stationary points at the finite
regime. Hence in order to search for stationary points at infinity we consider
the Poincaré map

$x=\frac{X}{\sqrt{1-X^{2}-Z^{2}}}~{},~{}z=\frac{Z}{\sqrt{1-X^{2}-Z^{2}}}~{},~{}d\sigma=\sqrt{1-X^{2}-Z^{2}}d\tau.$
(19)

Therefore, at infinity, $1-X^{2}-Z^{2}=0$, i.e. $Z=\sqrt{1-X^{2}}$, the
dynamical system is reduced to the single ordinary differential equation

$\frac{dX}{d\sigma}=\frac{\sqrt{3}\left(\lambda-4\right)}{12F_{0}-1}\left(X^{2}-1\right).$
(20)

Consequently, the stationary points at the infinity
$B=\left(X\left(B\right),Z\left(B\right)\right)$ are

$B_{1}=\left(1,0\right)\text{ and }B_{2}=\left(-1,0\right)\text{.}$ (21)

Hence, at infinity only the scalar field $\phi$ contributes to the
cosmological solution and the physical properties of the points are similar to
those of $A_{1}$.

The eigenvalues of the two-dimensional dynamical system at the stationary
points at infinity are

$e_{1}\left(B_{1}\right)=\frac{\sqrt{3}\left(\lambda-4\right)}{12F_{0}-1}~{},~{}\operatorname{Re}\left(e_{2}\left(B_{1}\right)\right)=0$
(22)

and

$e_{1}\left(B_{2}\right)=-\frac{\sqrt{3}\left(\lambda-4\right)}{12F_{0}-1}~{},~{}\operatorname{Re}\left(e_{2}\left(B_{2}\right)\right)=0.$
(23)

From the phase-space portraits of Fig. 3 we can easily observe that stationary
point $B_{1}$ is an attractor when $e_{1}\left(B_{1}\right)<0$, while, when
$e_{1}\left(B_{1}\right)>0,$ point $B_{2}$ is an attractor.

Figure 3: Phase-space portrait for the scalar-tensor theory on the two-
dimensional plane $\left(X,Z\right)$ for $Z\geq 0,$ for different values of
the free parameters $F_{0}$ and $\lambda$. From the phase-space portaits we
observe that points $B_{1}$ and $B_{2}$ can be attractors at the infinity
regime for the dynamical system.

##  III Scalar-Torsion Cosmology

Scalar-torsion theory can be seen as the extension of scalar-tensor theory in
teleparallelism. The gravitational Action Integral is te2

$S_{Torsion}=\int
d^{4}x\sqrt{-e}\left[F\left(\left|\psi\right|\right)T+\frac{1}{2}g^{\mu\nu}\psi_{,\mu}\psi_{,\nu}^{\ast}-V\left(\left|\psi\right|\right)\right],$
(24)

where $T$ is the torsion scalar of the antisymmetric Weitzenböck connection,
$e=\sqrt{-g}$, $e_{i}$ describes the unholonomic frame, with
$g(e_{i},e_{j})=\mathbf{e}_{i}\cdot\mathbf{e}_{i}=\eta_{ij}$ or in terms of
coordinates, $e_{i}=h_{i}^{\mu}\left(x^{\kappa}\right)\partial_{i}$, where now
the Weitzenböck connection is expressed as revtel

$\hat{\Gamma}^{\lambda}{}_{\mu\nu}=h_{a}^{\lambda}\partial_{\mu}h_{\nu}^{a}~{}$
(25)

and $T={S_{\beta}}^{\mu\nu}{T^{\beta}}_{\mu\nu}$, with $T_{\mu\nu}^{\beta}$ to
be the torsion tensor
$T_{\mu\nu}^{\beta}=\hat{\Gamma}_{\nu\mu}^{\beta}-\hat{\Gamma}_{\mu\nu}^{\beta}$
and

${S_{\beta}}^{\mu\nu}=\frac{1}{2}({K^{\mu\nu}}_{\beta}+\delta_{\beta}^{\mu}{T^{\theta\nu}}_{\theta}-\delta_{\beta}^{\nu}{T^{\theta\mu}}_{\theta}),$
(26)

where ${K^{\mu\nu}}_{\beta}$ is the contorsion tensor

${K^{\mu\nu}}_{\beta}=-\frac{1}{2}({T^{\mu\nu}}_{\beta}-{T^{\nu\mu}}_{\beta}-{T_{\beta}}^{\mu\nu}).$
(27)

For the spatially flat FLRW spacetime with line element (3) we consider the
diagonal frame

$h_{~{}\mu}^{i}(t)=\mathrm{diag}(N\left(t\right),a(t),a(t),a(t)),$ (28)

from which we calculate revtel

$T=-6H^{2}.$ (29)

By replacing in the Action Integral (24) and substituting the complex scalar
field as $\psi=\phi e^{i\theta}$, we end up with the point-like Lagrangian
function

$L_{Torsion}\left(N,a,\dot{a},\phi,\dot{\phi},\psi,\dot{\psi}\right)=-\frac{6}{N}F\left(\phi\right)a\dot{a}^{2}+\frac{a^{3}}{2N}\left(\dot{\phi}^{2}+\phi\dot{\theta}^{2}\right)-a^{3}NV\left(\phi\right).$
(30)

In an analogue to the Brans-Dicke model we assume now that
$F\left(\phi\right)=F_{0}\phi^{2}$. This is not a random choice, the Dilaton
field, or the Brans-Dicke field for the specific potential function
$V\left(\phi\right)=V_{0}\phi^{2}$ admits a discrete symmetry known as the
Gasperini-Veneziano duality transformation. Similarly, the scalar-torsion
theory for $F\left(\phi\right)=F_{0}\phi^{2}$ and potential function
$V\left(\phi\right)=V_{0}\phi^{2}$ admits a discrete symmetry similar to the
Gasperini-Veneziano duality transformation. Thus, the
$F\left(\phi\right)=F_{0}\phi^{2}$ can be seen as the teleparallel Brans-Dicke
equivalent model.

Variation with respect to the dynamical variables in the Lagrangian function
(30) gives the field equations

$0=-6F_{0}\phi^{2}H^{2}+\frac{1}{2}\left(\dot{\phi}^{2}+\phi^{2}\dot{\theta}^{2}\right)+V\left(\phi\right)~{},$
(31) $4F_{0}\phi^{2}\left(2\dot{H}+3H^{2}\right)+16F_{0}\phi
H\dot{\phi}+\left(\dot{\phi}^{2}+\phi^{2}\dot{\theta}^{2}\right)-V\left(\phi\right)=0~{},$
(32)
$\ddot{\phi}-\phi\dot{\theta}^{2}+3H\left(4F_{0}H\phi+\dot{\phi}\right)+V_{,\phi}=0~{}$
(33)

and

$\phi\ddot{\theta}+\dot{\theta}\left(3H\phi+2\dot{\phi}\right)=0~{}.$ (34)

As in the case of the scalar-tensor theory we proceed with the analysis for
the dynamics of the latter system of differential equations.

###  III.1 Cosmological dynamics

We make use of the dimensionless variables (12) and we write the field
equations (31)-(34) in the equivalent dynamical system

$F_{0}\frac{dx}{d\tau}=\left(4\sqrt{3}F_{0}+3x\right)x^{2}-2\sqrt{3}F_{0}\left(2F_{0}+\lambda
y-2z^{2}\right)-3x\left(F_{0}+y-z^{2}\right)~{},~{}$ (35)
$F_{0}\frac{dy}{d\tau}=y\left(2\sqrt{3}F_{0}\left(2+\lambda\right)x+3x^{2}+3\left(F_{0}-y+z^{2}\right)\right)~{},$
(36) $F_{0}\frac{dz}{d\tau}=\frac{3}{2}z\left(x^{2}-y+z^{2}-F_{0}\right)~{},$
(37) $\frac{d\lambda}{d\tau}=2\sqrt{3}\lambda
x\left(1-\lambda+\Gamma\left(\lambda\right)\right)~{},~{}~{}\Gamma\left(\lambda\right)=\phi\frac{V_{,\phi\phi}}{V_{,\phi}}$
(38)

with constraint equation

$x^{2}+y+z^{2}-F_{0}=0,$ (39)

where now the parameter for the equation of state for the effective fluid is
expressed as follows

$w_{eff}\left(x,y,z\right)=\frac{8}{\sqrt{3}}x+\frac{1}{F_{0}}\left(x^{2}-y+z^{2}\right)\text{.}$
(40)

We observe that the field equations are invariant under the change of variable
$z\rightarrow-z$. Thus without loss of generality we select to work in the
region $z\geq 0$. Moreover, from the constraint equation (39) the dimension of
the dynamical system is reduced by one, while, for the power-law potential
function $V\left(\phi\right)=V_{0}\phi^{\lambda}$, the field dimension of the
dynamical system is two.

The stationary points on the two-dimensional plane $\left\\{x,z\right\\}$, are
of the form $C=\left(x\left(C\right),z\left(C\right)\right)$, they are

$C_{1}=\left(\frac{F_{0}}{\sqrt{3}}\left(2+\lambda\right),0\right)~{},$
$C_{2}=\left(x_{2},\sqrt{F_{0}-\left(x_{2}\right)^{2}}\right)~{}.$

Point $C_{1}$ describes a universe in which the field $\phi$ and the potential
function contribute to the cosmological solution. The physical parameter
$w_{eff}$ is calculated to be
$w_{eff}\left(C_{1}\right)=\frac{1}{3}\left(2F_{0}\left(\lambda^{2}-4\right)-3\right)$
from which it easily follows that $C_{1}$ describes an accelerated asymptotic
cosmological solution for
$\left\\{\left|\lambda\right|<2,F_{0}>\frac{1}{\lambda^{2}-4}\right\\}$ or
$\left\\{\left|\lambda\right|>2,F_{0}<\frac{1}{\lambda^{2}-4}\right\\}$ while,
for $\left|\lambda\right|=2$, the de Sitter Universe is recovered. The
eigenvalues of the linearized system near to the stationary point are

$e_{1}\left(C_{1}\right)=-3+F_{0}\left(\lambda^{2}+2\right)~{},~{}e_{2}\left(C_{1}\right)=-3+F_{0}\left(\lambda^{2}+2\right)\text{~{}}.$
(41)

Thus $C_{1}$ is an attractor when $-3+F_{0}\left(\lambda^{2}+2\right)<0$.

The family of points $C_{2}$ exists when$~{}F_{0}>0$ and describes asymptotic
solutions with $w_{eff}\left(C_{2}\right)=1+\frac{8}{\sqrt{3}}x_{2}$ and
eigenvalues

$e_{1}\left(C_{2}\right)=6+2\sqrt{3}\left(2+\lambda\right)x_{2}~{},~{}e_{2}\left(C_{2}\right)=0.$

Because one of the eigenvalues is zero, in Fig. 4 we present the phase-space
portraits for the dynamical system in the two-dimensional plane
$\left(x,z\right)$ for positive values of $F_{0}$ and different parameters of
$\lambda$. From the diagrams we observe that when $e_{1}\left(C_{2}\right)<0$
the family of points $C_{2}$ describes attractors, otherwise the points are
sources.

Figure 4: Phase-space portrait for the scalar-torsion theory on the two-
dimensional plane $\left(x,z\right),$ for different values of the free
parameters $F_{0}$ and $\lambda$ in order to investigate the stability
properties of the family of points $C_{2}$ (solid lines). From the plots we
observe that, when $e_{1}\left(C_{2}\right)<0$, the points ae attractors,
otherwise they are source points

###  III.2 Analysis at Infinity

We consider now the Poincaré map (19) and we write the two-dimensional system
in the equivalent form

$F_{0}\frac{dX}{d\sigma}=\left(F_{0}-\left(1+F_{0}\right)\left(X^{2}+Z^{2}\right)\right)\left(\sqrt{3}F_{0}\left(2+\lambda\right)\left(X^{2}-1\right)-3X\sqrt{1-X^{2}-Z^{2}}\right)~{},$
(42)
$F_{0}\frac{dZ}{d\sigma}=\left(F_{0}-\left(1+F_{0}\right)\left(X^{2}+Z^{2}\right)\right)\left(\sqrt{3}F_{0}\left(2+\lambda\right)X-3\sqrt{1-X^{2}-Z^{2}}\right).$
(43)

The stationary points at infinity are
$D=\left(X\left(D\right),Z\left(D\right)\right)$ with
$1-X\left(D\right)^{2}-Z\left(D\right)^{2}=0$. We find the points

$D_{1}=\left(1,0\right)\text{ and }D_{2}=\left(-1,0\right)\text{,}$ (44)

from which it is clear that only the scalar field $\phi$ contributes in the
cosmological solution.

The eigenvalues of the linearized system are

$e_{1}\left(D_{1}\right)=-\sqrt{3}\left(2+\lambda\right)~{},~{}\operatorname{Re}\left(e_{2}\left(D_{1}\right)\right)\simeq-
sign\left(F_{0}\right)$ (45)

and

$e_{1}\left(D_{2}\right)=\sqrt{3}\left(2+\lambda\right)~{},~{}\operatorname{Re}\left(e_{2}\left(D_{2}\right)\right)\simeq-
sign\left(F_{0}\right).$ (46)

Therefore point $D_{1}$ is an attractor when $F_{0}>0$ and $2+\lambda>0$,
while $D_{2}$ is attractor for $F_{0}>0$ and $2+\lambda<0$.

In Fig. 5 we present phase-space portraits for the dynamical system of scalar-
torsion theory on the Poincaré variables.

Figure 5: Phase-space portraits for the scalar-torsion theory on the two-
dimensional plane $\left(X,Z\right),$ for different values of the free
parameters $F_{0}$ and $\lambda$.

##  IV Conclusions

In this piece of work we considered a complex scalar field in the context of
scalar-tensor and scalar-torsion theories in a spatially flat FLRW background
geometry. For these two gravitational models we derived the field equations
and we investigated the dynamical evolution of the physical quantities by
using dimensional variables. With the use of the latter variables the field
equations for both theories are reduced to a two-dimensional dynamical system
of first-order ordinary differential equations. In order to construct the
cosmological history as provided by the given models we determined the
stationary points and we investigated their stability properties.

For the two different dynamical systems which correspond to the scalar-tensor
and scalar-torsion theories respectively, we found that the stationary points
for the two dynamical systems are four. They are two points at the finite
regime and two points at infinity. The stationary points can describe
accelerated asymptotic solutions which can describe the early and late-time
acceleration phases of the universe. The stationary points and their stability
properties are summarized in Table 1.

We can easily conclude that the two theories, that is, the scalar-torsion
theory and the scalar-tensor theory have the same number of stationary points
and similar stability properties. The physical variables of the background
geometry at the stationary points have similar properties with a different
functional dependence on the free variables $F_{0}$ and $\lambda$. The
dynamical equivalence of the two theories is an expected result. Torsion is
not dynamical and there are no new degrees of freedom in scalar-torsion
cosmology with respect to scalar-tensor cosmology.

From the above analysis it is clear that the two theories cannot be
distinguished from the evolution of the background geometry and the analysis
of the perturbations should be performed. However, that overpasses the scopus
of this work and will be published elsewhere.

Finally, as far as the complex scalar field is concerned, we remark that there
only at the (families) stationary points $A_{2}$ and $C_{2}$ the imaginary
part of the complex scalar field contributes in the physical quantities of the
asymptotic solutions. However, these families of points can describe important
eras of the cosmological history, for instance, the early inflationary era,
the radiation era or the matter dominated era. What will be of special
interest is to investigate the existence of a scalar field potential which
does not depend only on the norm of the complex scalar field.

Table 1: Stationary points and their stability properties for the cosmological models of our consideration. Theory | Point | Finite/Infinity | Stability
---|---|---|---
Scalar-Tensor |  |  |
| $A_{1}$ | Finite | Attractor see Fig. 1
| $A_{2}$ | Finite | Saddle
| $B_{1}$ | Infinity | Attractor $\frac{\sqrt{3}\left(\lambda-4\right)}{12F_{0}-1}<0$
| $B_{2}$ | Infinity | Attractor $-\frac{\sqrt{3}\left(\lambda-4\right)}{12F_{0}-1}<0$
Scalar-Torsion |  |  |
| $C_{1}$ | Finite | Attractor $-3+F_{0}\left(\lambda^{2}+2\right)<0$
| $C_{2}$ | Finite | Attractor $6+2\sqrt{3}\left(2+\lambda\right)x_{2}<0$
| $D_{1}$ | Infinity | Attractor $F_{0}>0$, $2+\lambda>0$
| $D_{2}$ | Infinity | Attractor $F_{0}>0$, $2+\lambda>0$

###### Acknowledgements.

This work was partially supported by the National Research Foundation of South
Africa (Grant Numbers 131604). The author thanks Dr. G. Anargirou for the
hospitality provided while part of this work carried out.

## References

  * (1) B. Ratra and P.J.E Peebles, Phys. Rev. D 37 3406 (1988)
  * (2) J. Rubio, Front. Astron. Space Sci. 5, 50 (2019)
  * (3) L.N. Granda and D.F. Jimenez, EPJC 79, 772 (2019)
  * (4) P. Parsons and J.D. Barrow, Phys. Rev. D 51, 6757 (1995)
  * (5) J.D. Barrow and A. Paliathanasis, Gen. Rel. Grav. 50, 82 (2018)
  * (6) E.V Linder, Phys. Rev. D. 70 023511 (2004)
  * (7) S. Barshay and G. Kreyerhoff, EPJC 5, 369 (1998)
  * (8) T. Matos, J.-R. Luevano, I. Quiros, L.A. Urena-Lopez and J.A. Vazquez, Phys. Rev. D 80, 123521 (2009)
  * (9) A.R. Liddle, C. Pahud and L.A. Urena-Lopez, Phys. Rev. D 77, 121301 (2008)
  * (10) L.P. Chimento, V. Mendez and N. Zuccala, Class. Quantum Grav. 16, 3749 (1999)
  * (11) C. Deffayet, S. Deser and G. Esposito-Farese, Phys. Rev. D 80, 064015 (2009)
  * (12) D. Bertacca, N. Bartolo and S. Matarrese, Adv. Astron. 2010, 904379 (2010)
  * (13) D. Bertacca, A. Raccanelli, O.F. Piatella, D. Pietrobon, N. Bartolo, S. Matarrese and T. Giannantonio, JCAP 03, 039 (2011)
  * (14) A. Suárez and T. Matos, MNRAS 416, 87 (2011)
  * (15) C. Cosme, J.G. Rosa and O. Bertolami, Phys. Lett. B 781, 639 (2018)
  * (16) W. Hu, Phys. Rev. D 71, 047301 (2005)
  * (17) D. Wands, Lect. Notes Phys. 738, 275 (2008)
  * (18) A.A. Coley and R.J. van den Hoogen, Phys. Rev. D 62, 023517 (2000)
  * (19) A. Arbey, Phys. Rev. D 74, 043516 (2006)
  * (20) C.-J. Gao and Y.-G. Shen, Phys. Lett. B 541, 1 (2002)
  * (21) D.S.M. Alves and G.M. Kremer, JCAP 10, 009 (2004)
  * (22) V. Sivanesan, Phys. Rev. D 90, 104006 (2014)
  * (23) A. Paliathanasis, Universe 8, 325 (2022)
  * (24) I.M. Khalatnikov and A. Mezhlumian, Phys. Lett. A 169, 308 (1992)
  * (25) I.M. Khalatnikov, Lect. Notes Phys. 455, 343 (1995)
  * (26) D. Scialom, Helv. Phys. Acta 69, 190 (1996)
  * (27) P. Jetzer and D. Scialom, Phys. Rev. D 55, 7440 (1997)
  * (28) J.-A. Gu and W.-Y. Hwang, Phys. Lett. B 517, 1 (2001)
  * (29) A.Yu.Kamenshchik, I.M.Khalatnikov and A.V.Toporensky, Phys. Lett. B 357, 36 (1995)
  * (30) M. Arik, M. Calik and N. Katirci, Central Eur. J. Phys. 9, 1465 (2011)
  * (31) G. Rosen, EPL 89, 19002 (2010)
  * (32) Y.-G. Shen and X.-H. Ge, Int. J. Theor. Phys. 45, 17 (2006)
  * (33) A. Arbey, Phys. Rev. D 74, 043516 (2006)
  * (34) R.C.G. Landim, EPJC 76, 1 (2016)
  * (35) B. Li, P.R. Shapiro and T. Rindler-Daller, PoS BASH2015, 028 (2016)
  * (36) A. Bevilacqua, J. Kowalski-Glikman and W. Wislicki, Phys. Rev. D 105, 105004 (2022)
  * (37) H. Foidl and T. Rindler-Daller, Phys. Rev. D 105, 123534 (2022)
  * (38) V. Faraoni, Cosmology in Scalar-Tensor Gravity, Fundamental Theories of Physics, Springer Dordrecht (2004)
  * (39) C. Brans and R.H. Dicke, Phys. Rev. 124, 195 (1961)
  * (40) P. Jordan, Nature 164, 637 (1937) and Schwerkraft und Weltfall, 2nd ed., Vieweg und Sohn, Braunschweig, (1955)
  * (41) K. Hayashi and T. Shirafuji, Phys. Rev. D 19, 3524 (1979)
  * (42) M. Tsamparlis, Phys. Lett. A 75, 27 (1979)
  * (43) V. Faraoni, Phys. Rev. D 59, 084021, (1999)
  * (44) J. Garcia-Bellido, D. Wands, Phys. Rev. D 52, 6739 (1995)
  * (45) S.Sen and A.A. Sen, Late time acceleration in Brans Dicke Cosmology, Phys. Rev. D 63, 124006 (2001)
  * (46) M. Artymowski, Z. Lalak and M. Lewicki, Inflation and dark energy from the Brans-Dicke theory, JCAP 06, 031 (2015)
  * (47) M. Gasperini and G. Veneziano, Astropart. Phys. 1, 317 (1993)
  * (48) L. Jarv and J. Lember, Universe 7, 178 (2021)
  * (49) M. Hohmann, Phys. Rev. D 98, 064002 (2018)
  * (50) M. Hohmann and C. Pfeifer, Phys. Rev. D 98, 064003 (2018)
  * (51) H. Wei, Phys. Dynamics of Teleparallel Dark Energy, Lett. B 712, 430 (2012)
  * (52) C.-Q. Geng, C.-C. Lee, E.N. Saridakis and Y.-P. Wu, ”Teleparallel” Dark Energy, Phys. Lett. B 704, 384 (2011)
  * (53) C. Xu, E.N. Saridakis and G. Leon, Phase-Space analysis of Teleparallel Dark Energy, JCAP 07, 005 (2012)
  * (54) A. Paliathanasis, Eur. Phys. J. Plus 136, 674 (2021)
  * (55) E.J. Copeland, A.R. Liddle and D. Wands, Phys. Rev. D 57, 4686 (1998)
  * (56) . Wainwright and G.F.R Ellis, Cambridge University Press, Cambridge (1997)
  * (57) A.A. Coley, Dynamical Systems and Cosmology, Springer, Dordrecht (2003)
  * (58) L. Amendola and S. Tsujikawa, Dark Energy, Cambridge University Press, Cambridge (2010)
  * (59) C.R. Fadragas and G. Leon, Class. Quantum Grav. 31, 195011 (2014)
  * (60) O. Hrycyna and M. Szydlowski, JCAP 1312, 016 (2013)
  * (61) S. Bahamonte, K.F. Dialektopoulos, C. Escamilla-Rivera, V. Gakis, M. Hendry, J.L. Said, J. Mifsud and E. Di Valentino, Teleparallel Gravity: From Theory to Cosmology, [arXiv:2106.13793] (2021)


# Data-driven stochastic modeling of coarse-grained dynamics with finite-size
effects using Langevin regression

Jordan Snyder<EMAIL_ADDRESS>Jared L. Callaham  Steven L. Brunton  J. Nathan
Kutz  Department of Applied Mathematics, University of Washington  Department
of Mechanical Engineering, University of Washington

###### Abstract

Obtaining coarse-grained models that accurately incorporate finite-size
effects is an important open challenge in the study of complex, multi-scale
systems. We apply Langevin regression, a recently developed method for finding
stochastic differential equation (SDE) descriptions of realistically-sampled
time series data, to understand finite-size effects in the Kuramoto model of
coupled oscillators. We find that across the entire bifurcation diagram, the
dynamics of the Kuramoto order parameter are statistically consistent with an
SDE whose drift term has the form predicted by the Ott-Antonsen ansatz in the
$N\to\infty$ limit. We find that the diffusion term is nearly independent of
the bifurcation parameter, and has a magnitude decaying as $N^{-1/2}$,
consistent with the central limit theorem. This shows that the diverging
fluctuations of the order parameter near the critical point are driven by a
bifurcation in the underlying drift term, rather than increased stochastic
forcing.

##  1 Introduction

Collective behavior is an important and fascinating set of phenomena that
occurs in a wide array of both naturally occurring and engineered systems.
Examples include flocking and collective foraging of groups of animals,
synchronous flashing of fireflies, correlated activity of neurons in sensory
systems, and rolling blackouts in electric power grids. Characterizing such
phenomena relies on a systematic construction of coarse-grained, reduced-order
models that parsimoniously capture the physics of system-wide (macro-scale)
observables. For instance, the Kuramoto model of coupled oscillators, with its
synchronization phase transition, is a paradigmatic model of collective
behavior, with the construction of reduced-order models having been
successfully developed over the many years since its introduction. In
addition, much is known about the statistics of fluctuations near the phase
transition for finite system sizes, and their scaling with system size has
been shown to be anomalous relative to typical statistical physics problems.
However, comparatively little work has been done to systematically incorporate
finite-size effects in reduced-order models of coupled oscillators. In this
paper, we apply Langevin regression, a recently introduced method for
extracting SDE descriptions of realistic time-series data [8], to the problem
of obtaining physically interpretable stochastic models of finite-size effects
in the Kuramoto model.

The determination of collective behavior is related to the construction of
closure models which aim to characterize the effect of unresolved degrees of
freedom upon observed coarse-grained dynamics. The unresolved and unmeasured
degrees of freedom lead to the manifestation of randomness and memory effects,
even when the underlying laws of motion are deterministic and memoryless [57,
58, 31]. Thus it is natural to seek a stochastic differential equation (SDE)
as a coarse-grained model. We may further desire that our coarse-grained model
be interpretable, suggesting that it may be appropriate to regress to a given
family of SDEs.

Coarse-grained model reductions are often considered with the canonical
Kuramoto model of coupled oscillators [28], which has long been studied as an
example of collective behavior [3, 48, 17, 1, 30]. The Kuramoto model features
a composition of many coupled, nonlinear oscillatory units with distinct
natural frequencies, where the pairwise coupling tends to drive phases to
synchronicity. Importantly, the Kuramoto model provides tractable coarse-
grained approximations in certain special limits [47, 34] despite its rich
behavioral phenomenology [36, 5, 55]. In addition, the Kuramoto model and its
variants describe a number of synchronization phenomena in domains of diverse
nature such as coupled Josephson junctions [53], neuroscience [50], chemical
oscillators [27, 56], and the power grid [15]. The synchronization observed in
the Kuramoto model allows for asymptotic approximations of a closure model. In
the seminal work of Ott and Antonsen [34], they showed that in the
$N\to\infty$ limit and under certain conditions on the distribution of natural
frequencies and initial phases, the center of mass of a population of Kuramoto
oscillators remarkably obeys an autonomous ODE which undergoes a pitchfork
bifurcation as the coupling strength is varied, i.e. an exact closure model
can be constructed. Other studies have focused on complex coupling topologies,
proposing techniques using spectral information to merge nodes together [20]
or otherwise systematically discarding irrelevant degrees of freedom [26].
Still others take a more strictly data-driven approach and seek closed
equations of motion for low-order moments of the distribution of phases [33,
39] or to identify good coarse-grained variables via manifold learning
techniques [49]. Finally there are approaches that employ a “collective
coordinate” ansatz governing the phase of each oscillator within a phase-
locked cluster, and thereby arrive at a closed equation of motion [21, 22, 44,
45, 54].

Our aim is to construct an order-parameter description for finite-sized
Kuramoto systems, and more broadly provide mathematical methods for
characterizing finite-sized coarse-grained models. This has been studied
previously using nonequilibrium statistical mechanics [7], whose results are
indirect and difficult to interpret. Others have considered the ensemble-
variance in the steady state order parameter magnitude, and dynamical
fluctuations under quenched disorder [10, 11, 12, 13, 14, 24, 23]. While
precise finite-size scaling results can be obtained for ensemble statistics of
the steady-state order parameter, investigations of dynamical fluctuations
were limited to numerical simulation. Most importantly, investigations of
dynamical fluctuations did not take the step of interpreting the order
parameter dynamics as governed by an SDE. In this paper, we build on prior
work concerning coarse-grained modeling of Kuramoto oscillators by applying
Langevin regression to simulated trajectories. We obtain physically meaningful
SDE models for finite-size effects, demonstrating that the exotic finite-size
scaling behavior near the critical point is consistent with white noise
forcing with a magnitude that scales as $N^{-1/2}$. We conclude that critical
finite-size scaling behavior is driven mainly by instabilities in the
underlying dynamics rather than anomalously large forcing.

Figure 1: Emergence of synchronization in a system of $N=32$ coupled Kuramoto
oscillators, and the associated order parameter, at the level of individual
phases (left) and probability density (right). On the left, each colored line
is the phase of one oscillator through time, plotted on a circle in the plane
perpendicular to the $t$-axis. Line color corresponds to the oscillator’s
natural frequency $\omega_{i}$. On the right, color indicates the number
density of oscillators as a function of phase $\theta$ and time $t$. In both
plots, the solid black curve depicts the order parameter
$z=\langle\exp(\mathrm{i}\theta_{j})\rangle_{j=1}^{N}$ over time, and the
dashed line is the origin. The trajectory of the order parameter is also
projected onto the back plane. Initially, the oscillators’ phases are spread
evenly around the circle and $|z|\sim 0$. Due to the attractive coupling, the
oscillators’ phases quickly approach one another, forming a synchronized
cluster of oscillators that precesses with a single frequency. The magnitude
of the order parameter quantifies how synchronized the system is as a whole,
while its phase tracks the phase of the dominant synchronized cluster. Note
that the oscillators participating in the synchronized cluster are those with
natural frequencies near the middle of the distribution, while oscillators
with large (positive or negative) natural frequencies continue to drift.

##  2 Background

Networks of coupled oscillators are an extremely popular object of study for
applied mathematicians interested in collective behavior, owing to their rich
phenomenology, relative tractability, and broad abstract representation of
dynamic behavior. A paradigmatic model of coupled oscillators is the Kuramoto
model:

$\dot{\theta_{i}}=\omega_{i}+{1\over
N}\sum_{j=1}^{N}K_{ij}\sin(\theta_{j}-\theta_{i}),\hskip
28.45274pti=1,2,\dots,N$ (1)

where $\theta_{i}$ is the _phase_ of oscillator $i$, $\omega_{i}$ is its
_natural frequency_ , and $(K_{ij})$ is the _coupling matrix_. A fundamental
result due to Kuramoto [28] is that if coupling is mean-field (i.e.
$K_{ij}\equiv K$ for all $i,j$) and natural frequencies are distributed
according to a symmetric, unimodal probability density, then there is a
critical coupling strength $K_{c}$ such that for $K>K_{c}$ the system is
(partially) synchronized while for $K<K_{c}$ the system is incoherent. In
particular, synchrony can be measured by the _order parameter_
$z\coloneqq\langle\exp(\mathrm{i}\theta_{j})\rangle_{j=1}^{N}\in\mathbb{C}$,
whose magnitude $|z|$ is near zero if oscillators’ phases are spread evenly
around the unit circle, and is nonzero if oscillators’ phases break rotational
symmetry and cluster around a preferred phase. See Fig. 1 for a visualization.

Next, we review two major threads of research concerning the Kuramoto model.
First is research that seeks equations of motion for reduced or coarse-grained
degrees of freedom, towards obtaining a mechanistic understanding of
collective behavior, and is generally carried out in the thermodynamic
($N\to\infty$) limit. Second is research that seeks to understand
synchronization as a phase transition by investigating how the properties of
finite systems behave as $N\to\infty$. Our work bridges these threads by
learning stochastic reduced-order models for finite systems and studying their
scaling with $N$.

###  2.1 Dimension Reduction Methods

There has been extensive work studying the conditions under which the Kuramoto
model admits description as a lower-dimensional dynamical system [47, 2, 41,
4]. In the case that all natural frequencies are the equal, it is known that
for any $N\geq 3$, the full $N$-dimensional phase space is foliated by
invariant two-dimensional tori [52], and the dynamics on those tori are
generated by the action of the Möbius group [32]. In effect, this means that
the dynamics of the Kuramoto model with identical natural frequencies is two-
dimensional, regardless of the number of oscillators $N$.

Concerning different natural frequencies, Ott and Antonsen [34] showed that if
coupling is through a mean-field (i.e. $K_{ij}\equiv K$ for all $i,j$) and if
natural frequencies $\omega_{i}$ are drawn from a Cauchy distribution with
location $\Omega$ and scale $\delta$, then there is a closure for the order
parameter
${z\coloneqq\langle\exp(\mathrm{i}\theta_{j})\rangle_{j=1}^{N}\in\mathbb{C}}$,
namely

$\dot{z}=\mathrm{i}(\Omega+\mathrm{i}\delta)z+{K\over 2}(z-|z|^{2}z).$ (2)

This equation is valid in the $N\to\infty$ limit and technically only for
certain (distributional) initial conditions; a follow-on paper [35]
established under mild assumptions that the manifold in question is attractive
and so the above equation should hold after some initial transient, while
recent work has shown that the Ott-Antonsen (OA) manifold is not attracting
for finite $N$ [16]. Equivalent results hold for the case that the oscillators
are organized into modules, where within each module the oscillators’ natural
frequencies are drawn from a Cauchy distribution and the coupling strength
between oscillators $i$ and $j$ depends only on the modules they belong to. If
the distribution of natural frequencies is not Cauchy but a different rational
function of some order, then a similar reduction is possible, but with one
variable for each pair of poles of the natural frequency distribution [43].

Another analytical tool for deriving reduced models of coupled oscillator
systems is the “collective coordinate” method introduced by Gottwald [21].
This method represents the state of the population of oscillators by a shape
function together with a time-varying constant of proportionality, the
collective coordinate, such that the residual phase of any phase-locked
oscillator is given by the collective coordinate multiplied by the shape
function evaluated at its natural frequency. The resulting equations offer
remarkably accurate predictions for finite networks, including the scaling
behavior of the position of the critical point [21], have been applied to more
complex network structures [45, 44], and have been used as inductive biases to
guide data-driven discovery of coarse-grained models [46].

Finally, there are data-driven approaches to discovering coarse-grained models
for coupled oscillator systems. Prominent among these is the equation-free
(EF) framework, which estimates the time derivative of a coarse variable based
on projections of numerically integrated fine-grained dynamics. The EF
approach has been used to obtain accurate coarse-grained numerical models of
oscillator systems with a collective coordinate-type approach [33] and for
networks with spectral gaps [39]. More recently, this approach has been
combined with manifold learning techniques to learn the coarse-graining map
itself in addition to the dynamics of the coarse variables [49].

###  2.2 Finite-Size Scaling

While reduced-order modeling has proven to offer accurate and physically
meaningful descriptions of populations of coupled oscillators, little of this
work deals directly with the manner in which models for finite systems
converge as $N\to\infty$. On the other hand, much is known about the finite-
size scaling of certain summary statistics in the Kuramoto model from a
statistical physics standpoint.

The synchronization transition in the Kuramoto model is characterized by a
continuous yet abrupt change in synchrony and a diverging susceptibility, or
sensitivity to external forcing, closely related to temporal variance of the
order parameter [11, 12, 10, 13]. Both the abruptness of the change in
synchrony and the divergence of susceptibility are smoothed out in finite
systems, and that smoothing is described by finite-size scaling relations.
Progress on this problem for the Kuramoto model has focused on applying
scaling ansätze to the self-consistent equation that constrains the stationary
value order parameter, as well as to averages and variances of the order
parameter computed via numerical simulation [24, 23, 9].

While fundamental, the results described above still leave questions
unanswered regarding the origin of the particular scaling relations that are
observed. Through stochastic modeling over a range of system sizes, we
demonstrate that the scaling relations present in the Kuramoto model are
consistent with an SDE whose drift term has the form predicted by the
$N\to\infty$ theory due to Ott-Antonsen, and whose diffusion term, capturing
finite-size fluctuations, scales as $N^{-1/2}$, consistent with the central
limit theorem.

##  3 Methods

###  3.1 Model System and Analytics

We now compute the form of the SDE that we will fit to the dynamics of the
order parameter $z$. According to the Ott-Antonsen ansatz, $z$ should evolve
according to an ODE of the form

$\dot{z}=\lambda z+\mu|z|^{2}z$ (3)

where $\lambda=\mathrm{i}\Omega-\delta+K/2$ and $\mu=-K/2$. The solution to
this ODE either approaches $z=0$, or a fixed point or limit cycle with
$|z|^{2}=-\lambda/\mu$. However, it is clear from simulations that when $N$ is
finite, the value of $z$ exhibits fluctuations that persist through time (see
Fig. 2). The reason for these fluctuations is that when $N$ is finite, the
distribution of oscillators’ phases is necessarily discrete, and so cannot
approach the smooth distribution described by the OA anstaz that allows $z$ to
evolve autonomously.

\begin{overpic}[width=433.62pt]{figs/density-evolution-with-order-parameter-
annotated.pdf} \put(30.0,15.0){increasing synchrony} \put(88.0,22.5){
\parbox{50.58878pt}{finite-size effects}} \end{overpic} Figure 2:
Visualization of the emergency of synchronization in a system of $N=200$
Kuramoto oscillators in terms of the probability density (top) and order
parameter (bottom). The density plot shows that the oscillators’ phases
rapidly align with one another, forming an obvious peak in density starting at
time $t\gtrsim 5$. The degree of alignment can be quantified by the magnitude
of the order parameter $z$, depicted in orange in the bottom plot. Rather than
settling down to a constant value, the magnitude of the order parameter
fluctuates due to the system’s finite size. The angle of the order parameter
meanwhile closely tracks the location of the peak in oscillator density.
Oscillators’ natural frequencies were drawn from a Cauchy distribution with
width $\delta=1/2$ and coupled with coupling strength $K=2$.

We therefore suppose that when $N$ is finite, the statistics of
$z(t)=r(t)\exp(\mathrm{i}\Theta(t))$ will be well described by an SDE whose
drift term has the form of Eq. (3), and whose noise term quantifies finite-
size effects. Namely, we suppose that for finite $N$, the order parameter is
governed by an SDE of the form

$\mathrm{d}z=\left(\lambda
z+\mu|z|^{2}z\right)\mathrm{d}t+\sigma\exp(\mathrm{i}\Theta)\mathrm{d}W$ (4)

where $\sigma\in\mathbb{R}$ depends on $N$, and $W=W_{1}+\mathrm{i}W_{2}$ is a
complex-valued Brownian motion with independent real and imaginary parts
$W_{1}$ and $W_{2}$ respectively. Because $W_{1}$ and $W_{2}$ are independent,
adjusting the phase of the noise term by $\exp(\mathrm{i}\Theta)$ does not
affect the distribution of the solution, while it does simplify the next step
of the analysis. We emphasize that modeling finite-size effects as Brownian
motion reflects an assumption that fluctuations have a finite correlation
time, which may or may not be satisfied in practice, and that our goal is not
to exactly model a given trajectory $z(t)$, but rather its statistics.

In addition to fluctuations in time, finite Kuramoto systems may also exhibit
uncertainty in their parameter values. For instance,
$\lambda=\mathrm{i}\Omega-\delta+K/2$ depends on population statistics: the
location, $\Omega$, and scale, $\delta$, of the distribution of natural
frequencies. However, the realized set of natural frequencies is a finite
sample, and thus may not be precisely representative of the infinite
population. We therefore treat $\lambda$ and $\mu$ as unknowns to be fit to
observed data, in addition to $\sigma$.

For simplicity, we recast the complex-valued SDE (4) into polar coordinates
$z=r\exp(\mathrm{i}\Theta)$ so that

$\displaystyle\mathrm{d}r$
$\displaystyle=\left(\lambda_{r}r+\mu_{r}r^{3}+{\sigma^{2}\over
2r}\right)\mathrm{d}t+\sigma\mathrm{d}W_{1}$ (5)
$\displaystyle\mathrm{d}\Theta$
$\displaystyle=\left(\lambda_{i}+\mu_{i}r^{2}\right)\mathrm{d}t+{\sigma\over
r}\mathrm{d}W_{2}$ (6)

where subscripts $i$ and $r$ denote real and imaginary values of the
respective variables. Note the $\sigma^{2}/2r$ term in the drift function for
$r$, which arises due to the Itô formula for the change of variables [19] and
ensures that $r$ never becomes negative despite being subject to additive
noise. Since $r$ is the variable associated directly with the symmetry-
breaking synchronization transition, and since the dynamics of $r$ do not
depend on $\Theta$, we hereafter focus exclusively on the dynamics of $r$.

We wish to investigate whether or not the dynamics of $r$ for finite Kuramoto
systems are consistent with a model of the form (5), and if so, what are the
properties of the coefficients $\lambda,\mu$, and $\sigma$ as a function of
system size $N$ and coupling strength $K$?

###  3.2 Data Generation

To investigate these issues, we performed numerical simulations of the
Kuramoto model for different values of the number of oscillators $N$, and the
coupling strength $K$. For all experiments, natural frequencies were drawn
i.i.d. from a Cauchy distribution with location $\Omega=0$ and scale
$\delta=1/2$, meaning that the critical coupling strength is $K_{c}=2/\pi
g(0)=1$ [28]. We define the bifurcation parameter $\epsilon\coloneqq K-K_{c}$,
such that the symmetry-breaking transition occurs at $\epsilon=0.$ The number
of oscillators was taken to be each of $\\{64,256,1024,4096\\}$, $\epsilon$
was set at 20 evenly spaced values between -1 and 1. For each value of
$(N,\epsilon)$, we simulate an ensemble of 100 systems, each with a different
sample of natural frequencies $\\{\omega_{i}\\}$ and initial phases
$\\{\theta_{i}(0)\\}$. For each sample we integrate the ODE (1) to a final
time $T=10^{3}$ with a time resolution of $\Delta t=10^{-3}$, using a 5th-
order Tsitouras scheme [37]. Finally we coarse-grained the resulting time
series $\\{\theta_{i}(t)\\}$ to obtain

$r(t)\exp(\mathrm{i}\Theta(t))\coloneqq{1\over
N}\sum_{j=1}^{N}\exp(\mathrm{i}\theta_{j}(t)).$ (7)

Our goal is to obtain an SDE model that is statistically consistent with the
observed time series $r(t)$, after discarding an initial transient.
Specifically, we seek an SDE of the form (5).

###  3.3 Inference Procedure

Based on the discussion in Sec. 3.1, we expect that the coarse-grained SDE
will be of the form

$\mathrm{d}r=\left(\xi_{0}r+\xi_{1}r^{3}+{\xi_{2}^{2}\over
2r}\right)\mathrm{d}t+\xi_{2}\mathrm{d}W_{t},$ (8)

where the parameters $\xi=(\xi_{0},\xi_{1},\xi_{2})$ are to be identified with
the recently proposed Langevin regression method for identifying stochastic
differential equations from data [8]. This method solves both the forward and
adjoint Fokker-Planck equations to enforce consistency with both the finite-
time conditional moments and the steady-state probability distribution (PDF),
both computed after discarding an initial transient. A schematic of the
overall procedure performed in this paper is given in Fig. 3.

Figure 3: Schematic of our analysis pipeline. We first sample $N$ i.i.d.
natural frequencies $\\{\omega_{i}\\}$ from a Cauchy distribution with
location $\Omega$ and scale $\delta$, then integrate the Kuramoto ODE (1)
forward in time from randomly assigned initial phases. We then coarse-grain
the trajectory $\\{\theta_{i}(t)\\}$ to obtain the complex order parameter
$z(t)=r(t)\exp(\mathrm{i}\Theta(t))$, and compute the joint distribution
$p(r(t),r(t+\tau))$. From this joint distribution we compute the first two
finite-time Kramers-Moyal averages (9), and finally use Langevin regression to
infer SDE coefficients $\xi$ that are consistent with the measured finite-time
Kramers-Moyal averages.

For each of the replicate experiments described above, we computed the
empirical PDF $\hat{p}(r)$ and the first two ($n=1,2$) empirical finite-time
conditional moments $\hat{m}^{(n)}_{\tau}(r)$, defined as

$\displaystyle\hat{m}^{(n)}_{\tau}(r)$
$\displaystyle=\left(\mathbb{E}\left[(r(t+\tau)-r(t))^{n}|r(t)=r\right]\right)^{{1\over
n}}.$ (9)

The power of $1/n$ ensures that the conditional moments of different orders
are dimensionally consistent with one another, so that they can be
appropriately weighted in the objective function that we define below. For a
stationary process near thermal equilibrium with a sampling rate much faster
than any natural time scale of the system, these moments can be directly
related to the drift and diffusion functions [40, 42, 18, 6]. However, for
many closure problems the neglected degrees of freedom have unresolved
dynamics that lead to non-Markovian memory effects upon truncation.

Sampling the system too quickly will then invalidate the assumption of white-
in-time forcing underlying the convergence of Eq. (9) to a useful estimator of
the drift and diffusion. On the other hand, coarse sampling rates lead to
significant finite-time distortion of the conditional moments [38]. For
example, the evolution of the order parameter is smooth at short times; Eq.
(9) would significantly underestimate the drift and diffusion coefficients at
fast sampling rates. We therefore compute the conditional moments using
$\tau=0.1$. Langevin regression uses the adjoint Fokker-Planck equation to
correct for finite-time effects [29, 25]. The empirical conditional moments
$\hat{m}_{\tau}^{(n)}(r)$ are then compared against the finite-time
conditional moments $m_{\tau}^{(n)}(r,\xi)$ that would be observed given a set
of parameters $\xi$ (see Ref. [8] for details). Because the computed
trajectories are finite, we estimate $\hat{m}_{\tau}^{(n)}(r)$ by computing a
histogram with equally-sized bins centered at locations
$\\{r_{i}|i=1,\dots,M\\}$. In the following, we use $M=40$ equally-sized bins
that cover full range of observed $r$ values.

Langevin regression solves the following optimization problem:

$\xi=\operatorname*{arg\rm{}min}_{\xi}\sum_{n=1}^{2}\sum_{i=1}^{M}w_{i}^{(n)}\left[m_{\tau}^{(n)}(r_{i},\xi)-\hat{m}_{\tau}^{(n)}(r_{i})\right]^{2}+\eta_{\text{KL}}\mathcal{D}_{\text{KL}}\left(\hat{p}(r)\|p(r,\xi)\right)+\eta_{\ell_{1}}\|s\odot\xi\|_{1},$
(10)

where $\odot$ denotes the entrywise (Hadamard) product of vectors, the weight
variables $w^{(n)}_{i}$ reflect pointwise uncertainty in the empirical
estimate of the moments, and the $\eta$ variables control the relative
contributions of each of the terms in the cost function. We have added two
regularizing terms. The Kullback-Leibler (KL) divergence
$\mathcal{D}_{\text{KL}}$ is a statistical measure of difference between
probability distributions, and is defined as

$\mathcal{D}_{\text{KL}}(p\|q)=\int p(x)\log\left({p(x)\over q(x)}\right).$
(11)

The KL regularization ensures that the inferred model has a steady-state PDF
comparable to the empirical one. Finally, we normalize the weights
$w^{(n)}_{i}$ so that they integrate to unity. Because the domain is discrete,
this means enforcing the condition

$\sum_{i}(r_{i+1}-r_{i})w^{(n)}_{i}=1,\quad n=1,2$ (12)

The final term in the optimization problem, the $\ell_{1}$ norm of the
rescaled coefficient vector $s\odot\xi$, is new relative to the original
method. The subcritical systems ($\epsilon<0$) are linearly stable about the
origin, introducing a redundancy between the linear and cubic terms. The
$\ell_{1}$ regularization promotes sparse coefficient vectors, encouraging the
optimization problem to select zero or near-zero coefficients where these are
consistent with the data. The rescaling vector $s$ allows us to adjust the
magnitudes of the entries of $\xi$ so that they are comparable to each other.
As we expect the noise coefficient $\xi_{2}$ to decay as $N^{-1/2}$, we use
$s=(1,1,N^{1/2})$. We used regularization parameters
$\eta_{\text{KL}}=10^{-3}$ and $\eta_{\ell_{1}}=10^{-4}$, as we found these
values to prevent both blowup and premature dropout of inferred coefficients.
Optimization was carried out using the Nelder-Mead method as implemented in
scipy.optimize [51].

##  4 Results

Figure 4: Left: Violin plots of the distribution of the temporal average of
$r$ as a function of bifurcation parameter $\epsilon=K-K_{c}$ and $N$. For
each value of $(N,\epsilon)$, we sample 100 Kuramoto systems and compute the
time-average of $r$ (after discarding the transient), and depict the
distribution of these 100 values as a vertical violin. Vertical lines with
tick-ends depict the max and min of each set of 100 values. Dashed line is
$r_{OA}$, the value of the fixed point of the ODE model derived from the Ott-
Antonsen ansatz, Eq. (2). Note that to the right of the critical point
($\epsilon>0$), values cluster around the $r_{OA}$ with a variance that
shrinks with $N$. To the left of the critical point, values also approach
$r_{OA}$ as $N$ increases, but are offset from zero due to the singularity at
the origin (see the term $\sigma^{2}/2r$ in the drift term of Eq. (5)). Right:
Violin plots of the distribution of the temporal variance in $r$, scaled by
$N$. As on the left, for each $(N,\epsilon)$ we observe 100 Kuramoto systems,
compute the temporal variance of $r$, and depict the distribution of these 100
values as a vertical violin. The ensure comparability across $N$, we scale the
variance by $N$. Compare to Hong et al., Figure 1 [24]. Note that while far
from the critical point ($\epsilon=0$) variance appears to be inversely
proportional to $N$ (since $\chi\coloneqq N(\langle r^{2}\rangle-\langle
r\rangle^{2})\sim\mathcal{O}(1)$), as would be predicted by the central limit
theorem, near the critical point $\chi$ exhibits a peak whose height grows
with $N$, meaning that the variance decays more slowly than $N^{-1}$

We first report summary statistics about the trajectories under study. Fig. 4
(left) shows the time-average of the order parameter, relative to the value
predicted by the Ott-Antonsen ansatz,
$r_{OA}=\max(0,\sqrt{\epsilon/(\epsilon+1)})$. To the right of the bifurcation
($\epsilon>0$), the distribution of the time-average of $r(t)$ has its peak
near the Ott-Antonsen value, and has a width that shrinks with $N$. On the
left, however, the presence of noise makes it impossible for $r$ to remain
near zero, leading to a distribution of time-average $r(t)$ values with a
nonzero peak whose location decreases with $N$.

Next, Fig. 4 (right) depicts the variance of $r(t)$ over time, as a function
of $\epsilon$ and $N$. In particular we are visualizing the variance scaled by
$N$, which we refer to as $\chi=N\langle(r-\langle r\rangle)^{2}\rangle$. If
at any given time $t$ the complex phases $\\{\exp(\mathrm{i}\theta_{i}(t))\\}$
were independent, then we would expect their mean, $z$, to exhibit finite-
sample fluctuations scaling as $N^{-1/2}$ by virtue of the central limit
theorem (CLT). That is, we would expect $\chi\sim\mathcal{O}(1)$, and this
appears to be the case away from the critical point. Near the critical point,
however, substantial correlations build up among the individual phases, and
CLT scaling is no longer a good approximation.

Based on extensive numerical study, Hong et al. report that the peak
$\chi_{\text{max}}$ of $\chi$ with respect to $\epsilon$ scales as
$\chi_{\text{max}}(N)\sim N^{0.40}$. Finally, Fig. 1 of [24] reports the
ensemble average of $\chi$ over 100-1000 trials, obtaining highly precise
estimates of the population average of $\chi$; we remark that the full
distribution of $\chi$ has nontrivial scaling properties with $N$. Namely, in
the far subcritical regime ($\epsilon<0$) the distribution appears to be
independent of $N$, while in the supercritical regime ($\epsilon>0$) the
distribution of $\chi$ appears to narrow with increasing $N$. We explore
explanations for this fact in Appendix A.

An overview of parameter values inferred via Langevin regression based on
trajectories of the Kuramoto order parameter are given in Fig. 5. In the top
two panels, the brown line depicts the value predicted by the theory derived
by Ott-Antonsen, namely $\xi_{0}=\epsilon/2$ and $\xi_{1}=-(\epsilon+1)/2$. In
the bottom-left panel, the horizontal line is at height $1/\sqrt{2}$.

Figure 5: Overview of $\xi$ parameters inferred by Langevin regression from
Kuramoto order parameter trajectories. As in Fig. 4, for each value of
$(N,\epsilon)$ we sample 100 Kuramoto systems, and for each system obtain an
estimate of SDE parameters $\xi=(\xi_{0},\xi_{1},\xi_{2})$. Vertical violins
depict the distribution of the corresponding 100 quantities. “Theory” lines in
the top panels are based on the Ott-Antonsen ansatz. $\xi_{2}$ is scaled by
$N^{1/2}$ so that values for all $N$ are comparable, and we can see that
$\xi_{2}$ is approximately consistent with $\xi_{2}=\sigma N^{-1/2}$ for a
global constant $\sigma\approx 1/\sqrt{2}$.

The key result in Fig. 5 is in the bottom-left panel, which depicts
distributions of inferred values of the noise intensity, $N^{1/2}\xi_{2}$.
Langevin regression clearly shows that a noise intensity scaling like
$N^{-1/2}$ is consistent with the data. Moreover, the distributions of
inferred values of $N^{1/2}\xi_{2}$ are centered on $1/\sqrt{2}$, for all $N$
and for all $\epsilon$. This remarkable regularity shows that when we regard
the Kuramoto order parameter as following an SDE, its critical fluctuations
can be explained by a simple (white) noise process applied to a drift term
that undergoes a bifurcation.

On the subcritical side of the bifurcation, inferred values for $\xi_{0}$
cluster around the theoretical value, with a spread that shrinks as $N$
increases. For $\xi_{1}$, however, inferred values show much larger spread.
This is because in the subcritical phase, $r$ is small, and so $r^{3}$ is very
small. Thus its coefficient, $\xi_{1}$, is very weakly constrained and highly
susceptible to blowup. The $\ell_{1}$ term in the cost function substantially
mitigates this problem, but does not eliminate it entirely.

On the supercritical side, we also see noticeable trends with $\epsilon$ and
$N$. Interestingly, inferred coefficient values do not cluster around the
theoretical line. Rather, they are systematically farther from zero (i.e.
$\xi_{0}$ is more positive and $\xi_{1}$ is more negative). This indicates
that finite Kuramoto systems exhibit drift behavior that is nontrivially
different from what is predicted in the $N\to\infty$ limit. Alternatively, it
is possible that a different noise model is more appropriate, and that with
such a noise model in place the drift coefficients would align with what we
expect from Ott-Antonsen theory. In particular, it is possible that the
fluctuations affecting the order parameter are correlated on a timescale long
enough that white noise is not a good approximation, meaning that colored
noise would be needed. Further research is warranted to derive possible
higher-order corrections to the Ott-Antonsen theory.

##  5 Discussion and Conclusions

We have shown that it is possible to use Langevin regression to regress data
from finite-size Kuramoto oscillator systems onto a simple class of SDEs
consisting of linear and cubic drift terms and additive white noise. Since in
the absence of noise, this model class contains the equation predicted by the
Ott-Antonsen ansatz for the $N\to\infty$ limit, we can evaluate the limiting
behavior of our learned models. We find that while the linear and drift terms
exhibit the correct direction of dependence on $K$, there are substantial
quantitative discrepancies even after accounting for variance across samples
of natural frequencies.

Remarkably, however, we find that the fluctuations associated with finite
system size are consistent with forcing by white noise at a magnitude that
scales with $N$ as $N^{-1/2}$. This demonstrates that finite-size effects in
complex systems can be effectively modeled as stochastic forcing that follows
CLT-like scaling, even as response fluctuations follow highly nontrivial
scaling laws, especially near critical transitions.

## Acknowledgements

JS, SLB, and JNK acknowledge funding support from the Air Force Office of
Scientific Research (AFOSR FA9550-19-1-0386). JLC acknowledges funding support
from the Department of Defense (DoD) through the National Defense Science &
Engineering Graduate (NDSEG) Fellowship Program.

##  Appendix A Scaling of $\chi$ distributions

We now discuss the scaling of the distributions of $\chi$ with $N$, for both
large and small values of the bifurcation parameter, $\epsilon$. In
particular, we see in Fig. 4 (right) that the distributions of $\chi$ are
approximately constant with respect to $N$ when $\epsilon$ is large and
negative, while they become more narrow with increasing $N$ when $\epsilon$ is
large and positive.

This phenomenon can be explained by the correlation between $\langle r\rangle$
and $\chi$ in the sub- and super-critical regimes, respectively, as visualized
in Fig. 6.

Figure 6: Relationship between $\langle r\rangle$ and $\chi$ in the
subcritical (left, $\epsilon\leq-0.5$) and supercritical (right, $\epsilon\geq
0.5$) regimes. Each dot corresponds a Kuramoto system, i.e. Eq. (1) with a
sample of natural frequencies $\\{\omega_{i}\\}$. In the subcritical regime,
$\chi$ is positively correlated with $\langle r\rangle$, while the opposite is
true in the supercritical regime.

We can explain these relationships by appealing to approximate SDE models that
hold when $r\approx 0$ and when $r\approx r_{0}\neq 0$. We begin with the SDE
(5):

$\mathrm{d}r=\left(\lambda r+\mu r^{3}+{\sigma^{2}\over
2r}\right)\mathrm{d}t+\sigma\mathrm{d}W$ (13)

First, when $r\approx 0$, as in the subcritical regime, the $r^{3}$ term is
small and we omit it, to obtain

$\mathrm{d}r=\left(\lambda r+{\sigma^{2}\over
2r}\right)\mathrm{d}t+\sigma\mathrm{d}W.$ (14)

In this case the steady-state PDF is

$p(r)={-2\lambda\over\sigma^{2}}r\exp\left({\lambda
r^{2}\over\sigma^{2}}\right)$ (15)

and its first two moments are

$\displaystyle\langle r\rangle$ $\displaystyle={\sqrt{\pi}\over
2}{\sigma\over\sqrt{-\lambda}}$ (16) $\displaystyle\langle r^{2}\rangle$
$\displaystyle={-\sigma^{2}\over\lambda}={4\over\pi}\langle r\rangle^{2}$ (17)

which implies that $\chi=N(\langle r^{2}\rangle-\langle
r\rangle^{2})=N({4\over\pi}-1)\langle r\rangle^{2}$. In particular, the
distribution of $\chi$ should scale with $N$ in the same manner as $N\langle
r\rangle^{2}$. Since in the subcritical regime, correlations between phases
are weak, it is reasonable to view $\langle r\rangle$ as the mean of $N$
i.i.d. samples, and therefore subject to the central limit theorem. The CLT
would predict that $\langle r\rangle$ is approximately normally distributed
with standard deviation $\mathcal{O}(1/\sqrt{N})$, so $N\langle r\rangle^{2}$
should be $\mathcal{O}(1)$ for all $N$.

In the supercritical regime, on the other hand, $r$ is not small. Rather, it
falls in a narrow band around some nonzero value, $r_{0}$. In this case we
suppose that the drift term of Eq. (5) can be approximated linearly:

$\mathrm{d}r=\overline{\lambda}(r-r_{0})\mathrm{d}t+\sigma\mathrm{d}W.$ (18)

where $\overline{\lambda}$ depends on $\lambda,\mu,\sigma$, and $r_{0}$. In
this case the steady-state PDF is Gaussian:

$p(r)=\sqrt{{\overline{\lambda}\over\pi\sigma^{2}}}\exp\left({\overline{\lambda}\over\sigma^{2}}(r-r_{0})^{2}\right)$
(19)

which (neglecting the bounds $0<r<1$) has moments

$\displaystyle\langle r\rangle$ $\displaystyle=r_{0}$ (20)
$\displaystyle\langle r^{2}\rangle$ $\displaystyle={-\sigma^{2}\over
2\overline{\lambda}}+r_{0}^{2}$ (21)

and thus $\chi=-N\sigma^{2}/\overline{\lambda}$.

Next, we can examine the dependence of $\chi$ on $r_{0}$ by partial
differentiation:

$\displaystyle{\partial\chi\over\partial r_{0}}$
$\displaystyle={\partial\chi\over\partial\overline{\lambda}}{\partial\overline{\lambda}\over\partial
r_{0}}$ (22) $\displaystyle={N\sigma^{2}\over
2\overline{\lambda}^{2}}\left[6\mu r_{0}+{\sigma^{2}\over r_{0}^{3}}\right].$
(23)

Clearly $N\sigma^{2}/2\overline{\lambda}^{2}$ is positive; we now show that
the term in square brackets is negative. To do this we multiply it by
$r_{0}^{3}$ and obtain

$\displaystyle r_{0}^{3}\left[6\mu r_{0}+{\sigma^{2}\over r_{0}^{3}}\right]$
$\displaystyle=6\mu r_{0}^{4}+\sigma^{2}$ (24)
$\displaystyle=6\left({-\sigma^{2}\over 2}-\lambda
r_{0}^{2}\right)+\sigma^{2}$ (25) $\displaystyle=-2\sigma^{2}-6\lambda
r_{0}^{2}<0$ (26)

where we have used the condition that $r_{0}$ is a root of the drift term,
namely that $\lambda r_{0}^{2}+\mu r_{0}^{4}+\sigma^{2}/2=0$, and that
$\lambda>0$ since we consider the supercritical regime. Thus,
$\partial\chi/\partial r_{0}<0$, and indeed we see that $\chi$ and $\langle
r\rangle$ are negatively related in the supercritical phase. Thus we should
expect the distribution of $\chi$ to be comparable to the distribution of
$\langle r\rangle$, i.e. it should become narrower with larger $N$.

These results can be understood in terms of the population of oscillators. In
the subcritical regime, when coupling is weak, both $\langle r\rangle$ and
$\chi$ are driven by the amount of correlation among the oscillators’ phases,
and are therefore positively related. In the supercritical regime, on the
other hand, fluctuations are driven by the drifting subpopulation of
oscillators, i.e. those whose frequency does not lock to the global frequency.
The larger the value of $r$, the fewer oscillators remain drifting, and
therefore the fewer oscillators there are to contribute to fluctuations in
$r$.

##  Appendix B Normally distributed natural frequencies

Here we consider stochastic modeling of the order parameter for Kuramoto
systems where the natural frequencies are drawn not from a Cauchy
distribution, but a normal (Gaussian) distribution. For comparability with the
Cauchy-distributed case, we take the standard deviation $\sigma$ to be
$\sqrt{\pi/8}$, so that the critical coupling strength is $K_{c}=1$. As
before, we take 20 values of $\epsilon=K-K_{c}$ evenly spaced between $-1$ and
$1$, and for each value of $\epsilon$ perform 100 replicate experiments
consisting of sampling natural frequencies and initial conditions, and
integrating the equations of motion. We do this for the number of oscillators
$N=\\{64,256,1024,4096\\}$.

Figure 7: Distributions of the temporal average $\langle r\rangle$ and scaled
temporal variance $\chi=N(\langle r^{2}\rangle-\langle r\rangle^{2})$ of the
order parameter as a function of number of oscillators $N$ and bifurcation
parameter $\epsilon$. Trends agree qualitatively with those for the Cauchy
case depicted in Fig. 4, including the width of the distribution of $\chi$
above vs. below the bifurcation. Note that for small $N$ and large $\epsilon$,
many sampled systems achieve complete synchronization, leading to nearly-
constant $r(t)$ and vanishingly small $\chi$. We omit these from the
statistics here, and therefore the remaining cases have comparatively large
values of $\chi$ (see the $N=2^{6}$ cases for $\epsilon\approx 1$).

In Fig. 7 we depict distributions of the temporal mean $\langle r\rangle$ and
temporal variance $\langle(r-\langle r\rangle)^{2}\rangle$ of the order
parameter as a function of $N$ and $\epsilon$ (cf. Fig. 4 for the Cauchy-
distributed case). As there is no known closed-form expression for the value
of the order parameter as $N\to\infty$ in the case that natural frequencies
are normally distributed, we omit the $N\to\infty$ theory line in the left
panel. Note that the value of $\chi$ decays with $\epsilon$ more quickly than
it does in the case of Cauchy-distributed natural frequencies. This is because
when natural frequencies are normally distributed, much less coupling is
required to synchronize a given fraction of the population of oscillators than
when natural frequencies are Cauchy-distributed.

Figure 8: Distributions of inferred drift and diffusion coefficients for
Kuramoto systems with normally-distributed natural frequencies (cf. Fig. 5).

In Fig. 8 we show the distribution of inferred drift and diffusion
coefficients for Kuramoto systems with normally distributed natural
frequencies, obtained via Langevin regression. Several key differences from
the Cauchy-distributed case (Fig. 5) are apparent. First, the absolute size of
the inferred coefficients is smaller (compare vertical scales). Next, the
distributions of the (rescaled) noise strength $N^{1/2}\xi_{2}$ are not
constant with $\epsilon$. Rather, their centers show mild increase with
$\epsilon$ up to the critical point, and decay towards zero as $\epsilon$
grows above the critical point. This indicates that when natural frequencies
are normally-distributed, the effective strength of noise decreases as more of
the population becomes synchronized. Perhaps surprisingly, this is not the
case when natural frequencies are Cauchy-distributed, owing to the extreme
heterogeneity of the Cauchy distribution.

Similarly to the Cauchy-distributed case, we also see significant noise in the
recovery of the cubic coefficient in the subcritical regime. On the other
hand, we also observe irregularity at the far-supercritical side of the
bifurcation. This is due to the fact that when $\epsilon$ is large, the system
can easily become completely synchronized, and the order parameter becomes
constant with time. In these cases a stochastic model is clearly inappropriate
and the inferred coefficients lose their meaning.

## References

  * [1] Alexandre Arenas, Albert Díaz-Guilera, and Conrad J. Pérez-Vicente.  Synchronization reveals topological scales in complex networks.  Physical Review Letters, 96(11):1–4, 2006.
  * [2] Alexandre Arenas, Jürgen Kurths, Yamir Moreno, Albert Díaz-Guilera, and Changsong Zhou.  Synchronization in complex networks.  Physics Reports, 469(3):1–56, 2008.
  * [3] Jürgen Aschoff and Rütger Wever.  The circadian system of man.  In Biological rhythms, pages 311–331. Springer, 1981.
  * [4] Christian Bick, Marc Goodfellow, Carlo R Laing, and Erik Andreas Martens.  Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review.  The Journal of Mathematical Neuroscience, 10(1):9, 12 2020.
  * [5] Christian Bick, Mark J Panaggio, and Erik A. Martens.  Chaos in Kuramoto oscillator networks.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7):071102, 7 2018.
  * [6] Lorenzo Boninsegna, Ralf Banisch, and Cecilia Clementi.  A Data-Driven Perspective on the Hierarchical Assembly of Molecular Structures.  Journal of Chemical Theory and Computation, 14(1):453–460, 2018\.
  * [7] Michael A. Buice and Carson C. Chow.  Correlations, fluctuations, and stability of a finite-size network of coupled oscillators.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 76(3), 2007.
  * [8] Jared L. Callaham, Jean-Christophe Loiseau, Georgios Rigas, and Steven L. Brunton.  Nonlinear stochastic modeling with Langevin regression.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 477(2250):20210092, 2021.
  * [9] Tommaso Coletta, Robin Delabays, and Philippe Jacquod.  Finite-size scaling in the Kuramoto model.  Physical Review E, 95(4):1–8, 2017.
  * [10] H. Daido.  Scaling behaviour at the onset of mutual entrainment in a population of interacting oscillators.  Journal of Physics A: Mathematical and General, 20(10), 1987.
  * [11] Hiroaki Daido.  Lower critical dimension for populations of oscillators with randomly distributed frequencies: A renormalization-group analysis.  Physical Review Letters, 61(2):231–234, 1988.
  * [12] Hiroaki Daido.  Intrinsic Fluctuation and Its Critical Scaling in a Class of Populations of Oscillators with Distributed Frequencies.  Progress of Theoretical Physics, 81(4):727–731, 4 1989.
  * [13] Hiroaki Daido.  Intrinsic fluctuations and a phase transition in a class of large populations of interacting oscillators.  Journal of Statistical Physics, 60(5-6):753–800, 1990.
  * [14] Hiroaki Daido.  Susceptibility of large populations of coupled oscillators.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 91(1):1–9, 2015.
  * [15] Florian Dörfler, Michael Chertkov, and Francesco Bullo.  Synchronization in complex oscillator networks and smart grids.  Proceedings of the National Academy of Sciences, 110(6):2005–2010, 2013.
  * [16] Jan R. Engelbrecht and Renato Mirollo.  Is the Ott-Antonsen manifold attracting?  Physical Review Research, 2(2):023057, 4 2020.
  * [17] T.D. Frank, A. Daffertshofer, C.E. Peper, P.J. Beek, and H. Haken.  Towards a comprehensive theory of brain activity: coupled oscillator systems under external forces.  Physica D: Nonlinear Phenomena, 144(1-2):62–86, 9 2000.
  * [18] Rudolf Friedrich, Joachim Peinke, Muhammad Sahimi, and M. Reza Rahimi Tabar.  Approaching complexity by stochastic methods: From biological systems to turbulence.  Physics Reports, 506(5):87–162, 2011.
  * [19] Crispin W. Gardiner.  Handbook of stochastic methods: For Physics, Chemistry and the Natural Sciences, 1996.
  * [20] David Gfeller and Paolo De Los Rios.  Spectral coarse graining and synchronization in oscillator networks.  Physical Review Letters, 100(17):1–4, 2008.
  * [21] Georg A. Gottwald.  Model reduction for networks of coupled oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(5):053111, 5 2015.
  * [22] Edward J. Hancock and Georg A. Gottwald.  Model reduction for Kuramoto models with complex topologies.  Physical Review E, 98(1):012307, 7 2018.
  * [23] Hyunsuk Hong.  Finite-size scaling in the system of coupled oscillators with heterogeneity in coupling strength.  Physical Review E, 96(1):1–7, 2017.
  * [24] Hyunsuk Hong, Hugues Chaté, Lei Han Tang, and Hyunggyu Park.  Finite-size scaling, dynamic fluctuations, and hyperscaling relation in the Kuramoto model.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 92(2):1–8, 2015.
  * [25] Christoph Honisch and Rudolf Friedrich.  Estimation of Kramers-Moyal coefficients at low sampling rates.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 83(6):1–6, 2011.
  * [26] Yuki Izumida and Hiroshi Kori.  Coarse-grained description of general oscillator networks.  arXiv preprint arXiv:1311.0917, pages 1–8, 2013.
  * [27] István Z Kiss, Yumei Zhai, and John L Hudson.  Emerging coherence in a population of chemical oscillators.  Science, 296(5573):1676–1678, 2002.
  * [28] Yoshiki Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In International symposium on mathematical problems in theoretical physics, volume 39, pages 420–422, Berlin/Heidelberg, 1975. Springer, Springer-Verlag.
  * [29] S. J. Lade.  Finite sampling interval effects in Kramers-Moyal analysis.  Physics Letters A, 373(41), 2009.
  * [30] D. Li, I. Leyva, J. A. Almendral, I. Sendiña-Nadal, J. M. Buldú, S. Havlin, and S. Boccaletti.  Synchronization interfaces and overlapping communities in complex networks.  Physical Review Letters, 101(16):2–5, 2008.
  * [31] Kevin K. Lin and Fei Lu.  Data-driven model reduction, Wiener projections, and the Koopman-Mori-Zwanzig formalism.  Journal of Computational Physics, 424:109864, 2021.
  * [32] Seth A. Marvel, Renato E. Mirollo, and Steven H. Strogatz.  Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action.  Chaos, 19(4), 2009.
  * [33] Sung Joon Moon, R. Ghanem, and Ioannis G. Kevrekidis.  Coarse graining the dynamics of coupled oscillators.  Physical Review Letters, 96(14):1–4, 2006.
  * [34] Edward Ott and Thomas M. Antonsen.  Low dimensional behavior of large systems of globally coupled oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 18(3):037113, 9 2008.
  * [35] Edward Ott and Thomas M. Antonsen.  Long time evolution of phase oscillator systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 19(2):023117, 6 2009.
  * [36] Mark J Panaggio and Daniel M. Abrams.  Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators.  Nonlinearity, 28(3):R67–R87, 2015.
  * [37] Christopher Rackauckas and Qing Nie.  Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia.  Journal of Open Research Software, 5(1), 2017.
  * [38] Mario Ragwitz and Holger Kantz.  Indispensable finite time corrections for Fokker-Planck equations from time series data.  Physical Review Letters, 87(25):254501–1, 2001.
  * [39] Karthikeyan Rajendran and Ioannis G. Kevrekidis.  Coarse graining the dynamics of heterogeneous oscillators in networks with spectral gaps.  Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 84(3):1–11, 2011.
  * [40] H. Risken.  The Fokker-Planck equation: methods of solution and applications.  Springer-Verlag Berlin Heidelberg, 1996.
  * [41] Francisco A. Rodrigues, Thomas K Dm Peron, Peng Ji, and Jürgen Kurths.  The Kuramoto model in complex networks.  Physics Reports, 610:1–98, 2016.
  * [42] S. Siegert, R. Friedrich, and J. Peinke.  Analysis of data sets of stochastic systems.  Physics Letters, Section A: General, Atomic and Solid State Physics, 243(5-6):275–280, 1998.
  * [43] Per Sebastian Skardal.  Low-dimensional dynamics of the Kuramoto model with rational frequency distributions.  Physical Review E, 98(2):1–11, 2018.
  * [44] Lachlan D. Smith and Georg A. Gottwald.  Chaos in networks of coupled oscillators with multimodal natural frequency distributions.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(9):093127, 2019.
  * [45] Lachlan D. Smith and Georg A. Gottwald.  Model reduction for the collective dynamics of globally coupled oscillators: From finite networks to the thermodynamic limit.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(9):093107, 9 2020.
  * [46] Jordan Snyder, Anatoly Zlotnik, and Andrey Y Lokhov.  Data-driven selection of coarse-grained models of coupled oscillators.  Physical Review Research, 2(4):043402, 12 2020.
  * [47] Steven H. Strogatz.  From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 143(1):1–20, 9 2000.
  * [48] Steven H. Strogatz, Charles M. Marcus, Robert M. Westervelt, and Renato E. Mirollo.  Collective dynamics of coupled oscillators with random pinning.  Physica D: Nonlinear Phenomena, 36(1-2):23–50, 1989.
  * [49] Thomas N. Thiem, Mahdi Kooshkbaghi, Tom Bertalan, Carlo R. Laing, and Ioannis G. Kevrekidis.  Emergent Spaces for Coupled Oscillators.  Frontiers in Computational Neuroscience, 14:1–31, 5 2020.
  * [50] Francisco Varela, Jean-Philippe Lachaux, Eugenio Rodriguez, and Jacques Martinerie.  The brainweb: phase synchronization and large-scale integration.  Nature reviews neuroscience, 2(4):229–239, 2001.
  * [51] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, and Paul van Mulbregt.  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature Methods, 17(3):261–272, 3 2020.
  * [52] Shinya Watanabe and Steven H. Strogatz.  Integrability of a globally coupled oscillator array.  Physical Review Letters, 70(16):2391–2394, 1993.
  * [53] Kurt Wiesenfeld, Pere Colet, and Steven H Strogatz.  Frequency locking in josephson arrays: Connection with the kuramoto model.  Physical Review E, 57(2):1563, 1998.
  * [54] Wenqi Yue, Lachlan D. Smith, and Georg A. Gottwald.  Model reduction for the Kuramoto-Sakaguchi model: The importance of nonentrained rogue oscillators.  Physical Review E, 101(6):062213, 6 2020.
  * [55] Yuanzhao Zhang, Zachary G. Nicolaou, Joseph D. Hart, Rajarshi Roy, and Adilson E. Motter.  Critical Switching in Globally Attractive Chimeras.  Physical Review X, 10(1):011044, 2 2020.
  * [56] Anatoly Zlotnik, Raphael Nagao, István Z. Kiss, and Jr-Shin Li.  Phase-selective entrainment of nonlinear oscillator ensembles.  Nature communications, 7(1):1–7, 2016.
  * [57] Robert Zwanzig.  Memory effects in irreversible thermodynamics.  Physical Review, 124(4):983–992, 1961.
  * [58] Robert Zwanzig.  Nonlinear generalized Langevin equations.  Journal of Statistical Physics, 9(3):215–220, 1973.