Correlations, fluctuations and stability of a finite-size network of coupled oscillators
Michael A. Buice and Carson C. Chow
Laboratory of Biological Modeling, NIDDK, NIH, Bethesda, MD
(Dated: November 19, 2018)
The incoherent state of the Kuramoto model of coupled oscillators exhibits marginal modes in
mean field theory. We demonstrate that corrections due to finite size effects render these modes
stable in the subcritical case, i.e. when the population is not synchronous. This demonstration is
facilitated by the construction of a non-equilibrium statistical field theoretic formulation of a generic
model of coupled oscillators. This theory is consistent with previous results. In the all-to-all case,
the fluctuations in this theory are due completely to finite size corrections, which can be calculated
in an expansion in 1/N , where N is the number of oscillators. The N → ∞ limit of this theory is
what is traditionally called mean field theory for the Kuramoto model.
I. INTRODUCTION
Systems of coupled oscillators have been used to de-
scribe the dynamics of an extraordinary range of phe-
nomena [1], including networks of neurons [2, 3], syn-
chronization of blinking fireflies [4, 5], chorusing of chirp-
ing crickets [6], neutrino flavor oscillations [7], arrays of
lasers [8], and coupled Josephson junctions [9]. A com-
mon model of coupled oscillators is the Kuramoto model
[10], which describes the evolution of N coupled oscilla-
tors. A generalized form is given by
θ̇i = ωi +
f(θj − θi) (1)
where i labels the oscillators, θi is the phase of oscillator
i, f(θ) is the phase dependent coupling, and the intrin-
sic driving frequencies ωi are distributed according to
some distribution g(ω). In the original Kuramoto model,
f(θ) = sin(θ). Here, we consider f to be any smooth
odd function. The system can be characterized by the
complex order parameter
Z(t) =
eiθj(t) ≡ r(t)eiΨ(t) (2)
where the magnitude r gives a measure of synchrony in
the system.
In the limit of an infinite oscillator system, Kuramoto
showed that there is a bifurcation or continuous phase
transition as the coupling K is increased beyond some
critical value, Kc [10]. Below the critical point the steady
state solution has r = 0 (the “incoherent” state). Be-
yond the critical point, a new steady state solution with
r > 0 emerges. Strogatz and Mirollo analyzed the lin-
ear stability of the incoherent state of this system us-
ing a Fokker-Planck formalism [11]. In the absence of
external noise, the system displays marginal modes as-
sociated with the driving frequencies of the oscillators.
However, numerical simulations of the Kuramoto model
for a large but finite number of oscillators show that the
oscillators quickly settle into the incoherent state below
the critical point. The paradox of why the marginally
stable incoherent state seemed to be an attractor in sim-
ulations was partially resolved by Strogatz, Mirollo and
Matthews [12] who demonstrated (within the context of
the N → ∞ limit) that there was a dephasing effect akin
to Landau damping in plasma physics which brought r to
zero with a time constant that is inversely proportional to
the width of the frequency distribution. Recently, Stro-
gatz and Mirollo have shown that the fully locked state
r = 1 is stable [13] but the partially locked state is again
marginally stable [14]. Although dephasing can explain
how the order parameter can go to zero, the question of
whether the incoherent state is truly stable for a finite
number of oscillators remains unknown. Even with de-
phasing, in the infinite oscillator limit the system still has
an infinite memory of the initial state so there may be
classes of initial conditions for which the order parameter
or the density exhibits oscillations.
The applicability of the results for the infinite size
Kuramoto model to a finite-size network of oscillators
is largely unknown. The intractability of the finite size
case suggests a statistical approach to understanding the
dynamics. Accordingly, the infinite oscillator theories
should be the limits of some averaging process for a finite
system. While the behavior of a finite system is expected
to converge to the “infinite” oscillator behavior, for a fi-
nite number of oscillators the dynamics of the system will
exhibit fluctuations. For example, Daido [15, 16] consid-
ered his analytical treatments of the Kuramoto model
using time averages and he was able to compute an ana-
lytical estimate of the variance. In contrast, we will pur-
sue ensemble averages over oscillator phases and driving
frequencies. As the Kuramoto dynamics are determin-
istic, this is equivalent to an average over initial phases
and driving frequencies. Furthermore, the averaging pro-
cess imparts a distinction between the order parameter
Z and its magnitude r. Namely, do we consider 〈Z〉 or
〈r〉 = 〈|Z|〉 to be the order parameter? This is important
as the two are not equal. In keeping with the density as
the proper degree of freedom for the system (as in the in-
finite oscillator theories mentioned above), we assert that
〈Z〉 is the natural order parameter, as it is obtained via
a linear transformation applied to the density.
Recently, Hildebrand et al. [17] produced a kinetic
theory inspired by plasma physics to describe the fluctu-
ations within the system. They produced a Bogoliubov-
Born-Green-Kirkwood-Yvon (BBGKY) moment hierar-
http://arxiv.org/abs/0704.1650v3
chy and truncation at second order in the hierarchy
yielded analytical results for the two point correlation
function from which the fluctuations in the order param-
eter could be computed. At this order, the system still
manifested marginal modes. Going beyond second or-
der was impractical within the kinetic theory formalism.
Thus, it remained an open question as to whether going
to higher order would show that finite size fluctuations
could stabilize the marginal modes.
Here, we introduce a statistical field theory approach
to calculate the moments of the distribution function gov-
erning the Kuramoto model. The formalism is equivalent
to the Doi-Peliti path integral method used to derive sta-
tistical field theories for Markov processes, even though
our model is fully deterministic [18, 19, 20, 21]. The
field theoretic action we derive produces exactly the same
BBGKY hierarchy of the kinetic theory approach [17].
The advantages of the field theory approach are that 1)
difficult calculations are easily visualized and performed
through Feynman graph analysis, 2) the theory is eas-
ily extendable and generalizable (e.g. to local coupling),
3) the field theoretic formalism permits the use of the
renormalization group (which will be necessary near the
critical point), and 4) in the case of the all-to-all ho-
mogeneous coupling of the Kuramoto model proper, the
formalism results in an expansion in 1/N and verifies
that mean field theory is exact in the N → ∞ limit. We
will demonstrate that this theory predicts that finite size
corrections will stabilize the marginal modes of the infi-
nite oscillator theory. Readers unfamiliar with the tools
of field theory are directed to one of the standard texts
[22]. A review of field theory for non-equilibrium dynam-
ics is [23].
In section II, we present the derivation of the the-
ory and elaborate on this theory’s relationship to the
BBGKY hierarchy. In section III, we describe the com-
putation of correlation functions in this theory and, in
particular, describe the tree level linear response. This
will connect the present work directly with what was
computed using the kinetic theory approach [17]. In sec-
tion IV, after describing two example perturbations, we
calculate the one loop correction to the linear response
and demonstrate that the modes which are marginal at
mean field level are rendered stable by finite size effects.
In addition, we demonstrate how generalized Landau
damping arises quite naturally within our formalism. We
compare these results to simulations in section V.
II. FIELD THEORY FOR THE KURAMOTO
MODEL
The Kuramoto model (1) can be described in terms of
a density of oscillators in θ, ω space
η(θ, ω, t) =
δ(θ − θi(t))δ(ω − ωi) (3)
that obeys the continuity equation
C(η) ≡
f(θ′ − θ)η(θ′, ω′, t)η(θ, ω, t)dθ′dω′
= 0 (4)
Equation (4) remains an exact description of the dynam-
ics of the Kuramoto model (1) [17]. Although equa-
tion (4) has the same form as that used by Strogatz
and Mirollo [24], it is fundamentally different because
solutions need not be smooth. Rather, the solutions of
Eq. (4) are treated in the sense of distributions as defined
by Eq. (3). As we will show, imposing smooth solutions
is equivalent to mean field theory (the infinite oscilla-
tor limit). Drawing an analogy to the kinetic theory of
plasmas, Eq. (4) is equivalent to the Klimontovich equa-
tion while the mean field equation used by Strogatz and
Mirollo [24] is equivalent to the Vlasov equation [25, 26].
Our goal is to construct a field theory to calculate the
response functions and moments of the density η. Even-
tually we will construct a theory akin to a Doi-Peliti field
theory [18, 19, 20, 21], a standard approach for reaction-
diffusion systems. We will arrive at this through a con-
struction using a Martin-Siggia-Rose response field [27].
Since the model is deterministic, the time evolution of
η(θ, ω, t) serves to map the initial distribution forward in
time. We can therefore represent the functional proba-
bility measure P [η(θ, ω, t)] for the density η(θ, ω, t) as a
delta functional which enforces the deterministic evolu-
tion from equation (4) along with an expectation taken
over the distribution P0[η0] of the initial configuration
η0(θ, ω) = η(θ, ω, t0). We emphasize that no external
noise is added to our system. Any statistical uncertainty
completely derives from the distribution of the initial
state. Hence we arrive at the following path integral,
P [η(θ, ω, t)]
Dη0P0[η0]δ [N {C(η(θ, ω, t))− δ(t− t0)η0(θ, ω)}]
The definition of the delta functional contains an ar-
bitrary scaling factor, which we have taken to be N .
We will show later that this choice is necessary for the
field theory to correctly describe the statistics of η. The
probability measure obeys the normalization condition
DηP [η].
We first write the generalized Fourier decomposition of
the delta functional.
P [η(θ, ω, t)]
Dη̃Dη0P0[η0]
× exp
dθdωdt η̃ [C(η)− δ(t− t0)η0(θ, ω)]
where η̃(θ, ω, t) is usually called the “response field” after
Martin-Siggia-Rose [27] and the integral is taken along
the imaginary axis. It is more convenient to work with
the generalized probability including the response field:
P̃[η, η̃] =
Dη0P0[η0] exp
dθdωdt η̃C(η)
dθdωη̃(θ, ω, t0)η0(θ, ω)
which obeys the normalization
DηDη̃P̃[η, η̃; η0] (8)
We now compute the integral over η0 in equation (7),
which is an ensemble average over initial phases and driv-
ing frequencies. We assume that the initial phases and
driving frequencies for each of the N oscillators are inde-
pendent and obey the distribution ρ0(θ, ω). P [η0] repre-
sents the functional probability distribution for the ini-
tial number density of oscillators. Noting that η0(θ, ω) is
given by
η0(θ, ω) =
δ(θ − θi(0))δ(ω − ωi) (9)
Dη0P0[η0] =
dθidωiρ0(θi, ωi), one can show
that the distribution from equation (7) is given by
P̃ [η, η̃] = exp
dθdωdt η̃C(η)
+ N ln
eη̃(θ,ω,t0) − 1
ρ0(θ, ω)
In deriving Eq. (10), we have used the fact that
Dη0P0[η0] exp
dθdωη̃(θ, ω, t0)η0(θ, ω)
dθidωiρ0(θi, ωi) exp
η̃(θi, ωi, t0)
dθdωρ0(θ, ω)e
η̃(θ,ω,t0)
= exp
dθdωρ0(θ, ω)
eη̃(θ,ω,t0) − 1
We see that the fluctuations (i.e. terms non-linear in η̃)
appear only in the initial condition of (10), which is to
be expected since the Kuramoto system is deterministic.
In this form the continuity equation (4) appears as a
Langevin equation sourced by the noise from the initial
state.
Although the noise is contained entirely within the ini-
tial conditions, it is still relatively complicated. We can
simplify the structure of the noise in (10) by performing
the following transformation [21]:
ϕ(θ, ω, t) = η exp(−η̃)
ϕ̃(θ, ω, t) + 1 = exp(η̃) (12)
Under the transformation (12), P̃ [n, ñ] becomes P̃[ϕ, ϕ̃],
which is given by
P̃ [ϕ, ϕ̃] = exp (−NS[ϕ, ϕ̃]) (13)
where the action S[ϕ, ϕ̃] is
S[ϕ, ϕ̃] =
dωdθdt
dω′dθ′ (ϕ̃′ϕ̃+ ϕ̃)
{f(θ′ − θ)ϕ′ϕ}
dθdωϕ̃(θ, ω, t0)ρ0(θ, ω)
The form (14) is obtained from the transformation (12)
only after several integrations by parts. In most cases,
these integrations do not yield boundary terms because
of the periodic domain (i.e. those in θ). In the case of the
∂t operator, however, we are left with boundary terms of
the form [ln(ϕ̃ + 1) − 1]ϕ̃ϕ. These terms will not affect
computations of the moments because of the causality of
the propagator (see section III A).
We are interested in fluctuations around a smooth so-
lution ρ(θ, ω, t) of the continuity equation (4) with initial
condition ρ(θ, ω, t0) = ρ0(θ, ω). We transform the field
variables via ψ = ϕ−ρ and ψ̃ = ϕ̃ in (14) and obtain the
following action:
S[ψ, ψ̃] =
dωdθdtψ̃
dω′dθ′
{f(θ′ − θ) (ψ′ρ+ ρ′ψ + ψ′ψ)}
dω′dθ′ψ̃′
{f(θ′ − θ) (ψ′ + ρ′) (ψ + ρ)}
(−1)k+1
dθdωψ̃(θ, ω, t0)ρ0(θ, ω)
For fluctuations about the incoherent state: ρ(θ, ω, t) =
ρ0(θ, ω) = g(ω)/2π, where g(ω) is a fixed frequency dis-
tribution. The incoherent state is an exact solution of
the continuity equation (4). Due to the homogeneity in
θ and the derivative couplings, there are no corrections to
it at any order in 1/N . The action (15) with ρ = g(ω)/2π
therefore describes fluctuations about the true mean dis-
tribution of the theory, i.e. 〈η(θ, ω, t)〉 = g(ω)/2π. We
can evaluate the moments of the probability distribution
(13) with (15) using the method of steepest descents,
which treats 1/N as an expansion parameter. This is a
standard method in field theory which produces the loop
expansion [22]. We first separate the action (15) into
“free” and “interacting” terms.
S[ψ, ψ̃] = SF [ψ, ψ̃] + SI [ψ, ψ̃] (16)
where
SF [ψ, ψ̃]
dωdθdtψ̃
dω′dθ′
{f(θ′ − θ) (ψ′ρ+ ρ′ψ)}
dxdtdx′dt′ψ̃(x′, t′)Γ0(x, t;x
′, t′)ψ(x, t)
SI [ψ, ψ̃]
dωdθdtψ̃
dω′dθ′
{f(θ′ − θ) (ψ′ψ)}
dω′dθ′ψ̃′
{f(θ′ − θ) (ψ′ + ρ′) (ψ + ρ)}
(−1)k+1
dθdωψ̃(θ, ω, t0)ρ0(θ, ω)
In deriving the loop expansion, the action is expanded
around a saddle point, resulting in an asymptotic series
whose terms consist of moments of the Gaussian func-
tional defined by the terms in the action (15) which are
bilinear in ψ and ψ̃, i.e. SF [ψ, ψ̃]. Hence, the loop expan-
sion terms consist of various combinations of the inverse
of the operator Γ0, defined by SF , called the bare prop-
agator, with the higher order terms in the action, called
the vertices. Vertices are given by the terms in SI .
The terms in the loop expansion are conveniently rep-
resented by diagrams. The bare propagator is repre-
sented diagrammatically by a line, and should be com-
pared to the variance of a gaussian distribution. Each
term in the action (other than the bilinear term) with n
powers of ψ and m powers of ψ̃ is represented by a vertex
with n incoming lines and m outgoing lines. The initial
state vertices produce only outgoing lines and, like the
non-initial state or “bulk” vertices, are integrated over
θ, ω, and t for each point at which the operators are de-
fined. The bulk vertices are represented by a solid black
dot (or square, see Figure 1) and initial state vertices
by an open circle. The bare propagator and vertices are
shown in Figure 1. Unlike conventional Feynman dia-
grams used in field theory, the vertices in Figure 1 rep-
resent nonlocal operators defined at multiple points. In
particular, the initial state terms involve operators at a
different point for each outgoing line. Although uncon-
ventional, this is the natural way of characterizing the
1/N expansion.
Adopting the shorthand notation of x ≡ {θ, ω}, each
arm of a vertex must be connected to a line (propagator)
at either x or x′ and lines connect outgoing arms in one
vertex to incoming arms in another. The moment with
n powers of ψ and m powers of ψ̃ is calculated by sum-
ming all diagrams with n outgoing lines and m incoming
{ f (θ′−θ) . . .}
{ f (θ′−θ) . . .}
g(ω′)
{ f (θ′−θ) . . .}
{ f (θ′−θ)}
(2π)2
g(ω)g(ω′)
{ f (θ′−θ)}
P0(x, t|x
′)x x
(−1)k+1
ρ0(θ j,ω j)
FIG. 1: Diagrammatic (Feynman) rules for the fluctuations
about the mean. Time moves from right to left, as indicated
by the arrow. The bare propagator P0(x, t|x
′, t′) (see Eq. (31))
connects points at x′ to x, where x ≡ {θ, ω}. Each branch of
a vertex is labeled by x and x′ and is connected to a factor of
the propagator at x or x′. Each vertex represents an operator
given to the right of that vertex. The “. . . ” on which the
derivatives act only include the incoming propagators, but
not the outgoing ones. There are integrations over θ, θ′, ω, ω′
and t at each vertex.
lines. This means that each diagram will stand for sev-
eral terms which are equivalent by permutations of the
vertices or the edges in the graph, equivalently permuta-
tions of the factors of ψ̃ and ψ in the terms in the series
expansion. In a typical field theory, this results in combi-
natoric factors. In the present case, diagrams which are
not topologically distinct can produce different contribu-
tions to a given moment. Nonetheless, we will designate
the sum of these terms with a single graph. Generically,
the combinatoric factors we expect are due to the ex-
change of equivalent vertices, which typically cancel the
factorial in the series expansion. Additionally, each line
in a diagram contributes a factor of 1/N and each vertex
contributes a factor of N . Hence each loop in a diagram
carries a factor of 1/N . The terms in the expansion with-
out loops are called “tree level”. The bare propagator
is the tree level expansion of 〈ψ(θ, ω, t)ψ̃(θ′, ω′, t′)〉. The
tree level expansion of each moment beyond the first car-
ries an additional factor of 1/N , i.e. the propagator and
two-point correlators are each O(1/N).
Mean field theory is defined as the N → ∞ limit of this
field theory. In the infinite size limit, all moments higher
than the first are zero (provided the terms in the series
are not at a singular point, i.e. the onset of synchrony).
Hence, the only surviving terms in the action (15) are
those which contribute to the mean of the field at tree
level. These terms sum to give solutions to the continuity
equation (4). If the initial conditions are smooth, then
mean field theory is given by the relevant smooth solution
of (4). In most of the previous work (e.g. [12, 24]),
smooth solutions to (4) were taken as the starting point
and hence automatically assumed mean field theory.
We can now validate our choice of N as the correct
scaling factor in the delta functional of (5) by considering
the equal time two-point correlator. Using the definition
of η from (3) we get
〈η(x, t)η(x′, t)〉 = C(x, x′, t) + ρ(x, t)ρ(x′, t)
ρ(x, t)ρ(x′, t) +
δ(x − x′)ρ(x′, t)
where ρ(x, t) = 〈η(x, t)〉 = 〈ϕ(x, t)〉. Using the fields ϕ
and ϕ̃ (defined in (12)) and taking η at different times
gives
〈η(x, t)η(x′, t′)〉
= 〈[ϕ̃(x, t)ϕ(x, t) + ϕ(x, t)] [ϕ̃(x′, t′)ϕ(x′, t′) + ϕ(x′, t′)]〉
for t > t′. The response field has the property that ex-
pectation values containing ϕ̃(x, t) are zero unless an-
other field insertion of ϕ(x, t) is also present but at a
later time (this is because the propagator is causal; it is
zero for t− t′ ≤ 0). Therefore
〈η(x, t)η(x′, t′)〉 = 〈ϕ(x, t)ϕ(x′, t′)〉
+ 〈ϕ(x, t)ϕ̃(x′, t′)〉〈ϕ(x′, t′)〉 (21)
As we will show later when we discuss the propagator in
more detail, we have
〈ϕ(x, t)ϕ̃(x′, t′)〉 =
δ(x − x′) (22)
Comparing (21) in the limit t → t′ with (19) allows the
immediate identification of
〈ϕ(x, t)ϕ(x′, t)〉 = C(x, x′, t) + ρ(x, t)ρ(x′, t)
ρ(x, t)ρ(x′, t) (23)
C(x, x′, t) is the two oscillator “connected” correlation
or moment function. This is consistent with (15) which
gives
〈ϕ(x, t0)ϕ(x
′, t0)〉 = ρ(x, t0)ρ(x
′, t0)−
ρ(x, t0)ρ(x
′, t0)
as the initial condition. Thus, comparing the second mo-
ment of η using the Doi-Peliti fields (20) with the expres-
sion from the direct computation given by (19) shows
that the factor of N in the delta functional of (5) was
necessary to obtain the correct scaling for the moments.
The Doi-Peliti action (15) can also be derived by con-
sidering an effective Markov process on a circular lattice
representing the angle θ where the probability of an oscil-
lator moving to a new point on the lattice is determined
by its native driving frequency ω and the relative phases
of the other oscillators (see Appendix C). This is the
primary reason we refer to the theory as being of Doi-
Peliti type. The continuum limit of this process yields a
theory described by the action (15). The Markov picture
provides an intuitive description and underscores the fun-
damental idea that we have produced a statistical theory
obeyed by a deterministic process.
Although our formalism is statistical, we emphasize
that no approximations have been introduced. The sta-
tistical uncertainty is inherited from averaging over the
initial phases and driving frequencies. This formalism
could be applied to a wide variety of deterministic dy-
namical systems that can be represented by a distribu-
tional continuity equation like Eq. (4). In general, a solu-
tion for the moment generating functional for our action
(15) is as difficult to obtain as solving the original sys-
tem. The advantage of formulating the system as a field
theory is that a controlled perturbation expansion with
the inverse system size as the small parameter is possible.
A. Relation to Kinetic Theory and Moment
Hierarchies
The theory defined by the action (15) is equivalently
expressed as a Born-Bogoliubov-Green-Kirkwood-Yvon
(BBGKY) moment hierarchy starting with the continuity
equation (4). To construct the moment hierarchy, one
takes expectation values of the continuity equation with
products of the density, η(θ, ω, t). This results in coupled
equations of motion for the various moments of η(θ, ω, t),
where each equation depends upon one higher moment in
the hierarchy. The Dyson-Schwinger equations derivable
from (15) with 〈ϕ̃〉 = 0 are exactly the BBGKY hierarchy
derived in Ref. [17]. Thus the kinetic theory and field
theory approaches are entirely equivalent.
The moments of η can be computed in the BBGKY hi-
erarchy by truncating at some level. In Ref. [17], this was
done at Gaussian order by assuming that the connected
three-point function was zero. The first two equations
of the hierarchy then form a closed system which can by
solved. The first two equations of the BBGKY hierar-
chy [17] are
f(θ′ − θ)ρ(x, t)ρ(x′, t)dθ′dω′
f(θ′ − θ)C(x, x′, t)dθ′dω′ (25)
where ρ(x, t) = 〈η(x, t)〉, and the connected (equal-time)
correlation function
C(x, x′, t) = 〈η(x, t)η(x′, t)〉 − ρ(x, t)ρ(x′, t)
ρ(x, t)ρ(x′, t)−
δ(x− x′)ρ(x′, t)
obeys
f(θ3 − θ1) +
f(θ3 − θ2)]ρ(x3, t)dθ3dω3
C(x1, x2, t)
f(θ3 − θ1)ρ(x1, t)C(x2, x3, t)dθ3dω3 +K
f(θ3 − θ2)ρ(x2, t)C(x3, x1, t)dθ3dω3}
f(θ2 − θ1) +
f(θ1 − θ2)]ρ(x1, t)ρ(x2, t).
In the field theoretic approach, instead of truncating
the BBGKY hierarchy, one instead truncates the loop
expansion. Truncating the moment hierarchy at the mth
order is equivalent to truncating the loop expansion for
the lth moment at the (m − l)th order. Thus the solu-
tion to the moment equations (25) and (27) is the one
loop expression for the first moment and the tree level
expression for the second moment. The advantage of us-
ing the action (15) is that the terms in the perturbation
expansion are given automatically by the relevant dia-
grams at any level of the hierarchy. Ref. [17] suggested
that a higher order in the hierarchy would be necessary to
check whether the mean field marginal modes are stable
for finite N . We demonstrate below that the field the-
ory facilitates the calculation of the linearization of equa-
tion (25) to higher order in 1/N and show that marginal
modes are stabilized by finite size fluctuations.
One can compare this approach to the maximum en-
tropy approach of Rangan and Cai [28] for developing
consistent moment closures for such hierarchies. In the
moment hierarchy approach of Ref. [17], moment closure
is obtained via the somewhat ad hoc approach of setting
the nth cumulant to zero. In contrast, Rangan and Cai
maximize the entropy of the distribution subject to cer-
tain normalization constraints. The moment closure is
facilitated by constraining higher moments from the hi-
erarchy. However, one still must solve the resulting equa-
tions. In the loop expansion, moment closure is obtained
implicitly via truncating the loop expansion. The loop
expansion approach offers the advantage of providing a
natural means for determining when the approximation,
thus the implicit closure, breaks down and avoids deal-
ing with the moment hierarchy explicitly. In fact, Ran-
gan and Cai’s procedure has a natural interpretation in
field theory, namely the minimization of a generalized
effective action in terms of various moments. The sim-
plest and most common is the effective action in terms
of the mean field, which is the generating functional of
one partical irreducible (1PI) graphs [22]. The next level
of approximation is a generalized effective action (the
“effective action for composite operators” [29]) in terms
of the mean and the two-point function (or functions),
which is the generating functional of two particle irre-
ducible (2PI) graphs. One can continue in this way. The
equations of motion of these effective actions will produce
a closure of the moment hierarchy implicit in the action
for the theory. At tree level these equations will be equiv-
alent to those produced by Rangan and Cai’s maximum
entropy approach. The loop expansion allows for sys-
tematic corrections to these equations without explicitly
invoking higher equations in the hierarchy.
III. TREE LEVEL LINEAR RESPONSE,
CORRELATIONS AND FLUCTUATIONS
As a first example, we reproduce the calculation of the
variation of the order parameter Z, which was calculated
previously using the BBGKY moment hierarchy [17]. To
do so requires the calculation of the tree level linear re-
sponse or bare propagator and the tree level connected
two-oscillator correlation function.
A. The Propagator
The propagator P (θ, ω, t|θ′, ω′, t′) is given by the ex-
pectation value
P (θ, ω, t|θ′, ω′, t′) = 〈ϕ(θ, ω, t)ϕ̃(θ′, ω′, t′)〉 (28)
It is the linear response of (4). This can be shown by
considering a small perturbation δρ0 to the initial state
ρ in the action (15). Expanding to first order then yields:
δρ(θ, ω, t) = N
dθ′dω′ 〈ϕ(θ, ω, t)ϕ̃(θ′, ω′, t′)〉 δρ0(θ
′, ω′)
The tree level linear response or bare propagator
P0(θ, ω, t|θ
′, ω′, t′) ≡ P0(x, t;x
′, t′) is the functional in-
verse of the operator Γ0 defined by the free part of the
action (17). The bare propagator is therefore given by
Γ0 · P0 ≡
dx′′dt′′Γ0(x, t;x
′′, t′′)P0(x
′′, t′′;x′, t′)
δ(x− x′)δ(t− t′) (30)
Using the action (15) with Eq. (30) gives
Γ0 · P0 ≡
f(θ1 − θ)ρ(x1, t)dθ1dω1
P0(x, x
′, t− t′)
f(θ1 − θ)ρ(x, t)P0(x1, x
′, t− t′)dθ1dω1
δ(θ − θ′)δ(ω − ω′)δ(t− t′) (31)
Due to the rotational invariance in θ of f(θ),
P (x, t;x′, t′) ≡ P (θ − θ′, ω, ω′, t− t′).
In the incoherent state, ρ(θ, ω, t) = g(ω)/2π. Thus, for
f(θ) odd, Eq. (31) becomes
P0(x, x
′, t− t′)
f(θ1 − θ)P0(x1, x
′, t− t′)dθ1dω1
δ(θ − θ′)δ(ω − ω′)δ(t− t′) (32)
We can invert this equation using Fourier and Laplace
transforms.
Taking the Fourier transform of Eq. (32) (with respect
to θ and θ′) yields
+ inω
P0(n, ω;m,ω
′, t− t′)
+ inKg(ω)f(−n)
P0(n, ω1;m,ω
′, t− t′)dω1
δ(t− t′)δ(ω − ω′)δn+m (33)
where we use the following convention for the Fourier
transform
f(n) =
f(θ)e−inθdθ
f(θ) =
f(n)einθ (34)
Hereon, we will suppress the index m since the propaga-
tor must be diagonal (i.e. P0(n,m)) ∝ δm+n).
We Laplace transform in τ = t− t′ to get
[s+ inω] P̃0(n, ω, ω
′, s)
+ inKg(ω)f(−n)
P̃0(n, ω1, ω
′, s)dω1
δ(ω − ω′) (35)
using the convention
f̃(s) =
f(t)e−sτdτ
f(τ) =
f̃(s)esτds (36)
where the contour L is to the right of all poles in f̃(s).
We can solve for P̃0(n, ω, ω
′, s) using a self-consistency
condition. Integrate (35) over ω after dividing by s+ inω
to get
dωP̃0(n, ω, ω
′, s)
inKg(ω)f(−n)
s+ inω
dω1P̃0(n, ω1, ω
′, s)
s+ inω′
which we can solve to obtain
dωP̃0(n, ω, ω
′, s) =
s+ inω′
Λn(s)
where
Λn(s) = 1 + inKf(−n)
s+ inω
Λn(s) is defined for Re(s) ≤ 0 via analytic continuation.
In the kinetic theory context of an oscillator density obey-
ing the continuity equation (4), Λn(s) is analogous to a
plasma dielectric function [17]. If we assume that g(ω) is
even and f(θ) is odd, then there is a single real number
sn such that Λn(sn) = 0. (Mirollo and Strogatz proved
that there is at most one single, real root of (39) and that
it must satisfy Re(s) ≥ 0 [11, 30]. In our case, Λn(s) is
defined for Re(s) < 0 not by (39), but rather via ana-
lytic continuation.) Using (39) in (35) and solving for
P̃0(n, ω, ω
′, s) gives
P̃0(n, ω, ω
′, s) =
δ(ω − ω′)
s+ inω
inKg(ω)f(−n)
(s+ inω) (s+ inω′)
Λn(s)
Here we identify the spectrum with the zeroes of Γ0 or,
equivalently, the poles of the propagator, as these will de-
termine the time evolution of perturbations. Analogous
to the analysis of Strogatz and Mirollo [24] we define the
operator O by
O[bn(ω, t)] ≡ inωbn(ω, t)
+ inKf−n
bn(ω1, t)g(ω1)dω1 (41)
FIG. 2: Diagrams for the connected two-point function at tree
level a) and to one loop b).
(cf. equation (62) below). The continuous spectrum of
O consists of the frequencies inω whereas the discrete
spectrum (according to Ref. [24]) only exists for K >
Kc. Consistent with that approach (i.e. linear operator
theory), we identify the poles in P due to s+ inω as the
continuous spectrum and those due to the zeroes of Λn(s)
as the discrete spectrum. If Λn(s) is not analytically
continued for Re(s) < 0, it will not have zeroes for that
domain, as in Ref. [24]. However, zeros can exist for
Re(s) < 0 when Λn(s) is analytically continued and this
is why analytic continuation is of such crucial importance
to the conclusions of Ref. [12].
B. Correlation Function
The connected correlation function (cumulant func-
tion) is given by
C(x1, t1;x2, t2) = 〈ψ(x1, t1)ψ(x2, t2)〉 (42)
This is equivalent to (26), which was computed in
Ref. [17] when t1 = t2. If the initial phases are un-
correlated (i.e. C(x1, 0;x2, 0) = 0) then at tree level
C(x1, t1;x2, t2) is given by the diagram shown in Fig-
ure 2 a). It is comprised of vertex III (see Figure 1) com-
bined with a bare propagator on each arm. For general
(i.e. not odd) f(θ), the diagram in Figure 2 a) actually
corresponds to two different terms because the arms of
vertex III can be interchanged, giving two terms in equa-
tion (43) below. Unlike, conventional field theory, these
interchanges are not symmetric. These two terms are
equal when f(θ) is odd. More generally other vertices
do not exhibit the symmetries typical of Feynman dia-
grams even for odd f(θ). Applying the Feynman rules
then gives at tree level
C(x1, t1;x2, t2)
(2π)2
dω′1dω
dθ′1dθ
× P0(x1, x
1, t1 − t
′)P0(x2, x
2, t2 − t
f(θ′2 − θ
f(θ′1 − θ
2)]g(ω
1)g(ω
2)(43)
where t = min(t1, t2). This is essentially identical to the
ansatz used in Ref. [17] for the solution of the second
moment equation in the BBGKY hierarchy (27). For
f(θ) = sin θ and t1 > t2, Fourier transforming Eq. (43)
and inserting the bare propagator from Eq. (40) (after
an inverse Laplace transformation) gives
C1(ω1, t1, ω2, t2) =
g(ω1)g(ω2)
(iω1 +
)(−iω2 +
(iω1 +
)(−iω2 +
e(i(ω1−ω2)t2)
i(ω1 − ω2)
i(ω1 − ω2)
eiω1(t1−t2)
(iω2 −
− iω1)((
)2 + (ω2)2)
+iω2)(t2) − 1
)(t1−t2)
(−iω1 −
+ iω2)((
)2 + (ω1)2)
−iω1)t2 − 1
e−iω1(t1−t2)
4(K − 2Kc
− iω1)(
+ iω2)
e−(Kc−K)t2 − 1
e−(Kc−K))(t1−t2)
The equal time correlator is given by (t1 = t2 = τ):
C1(ω1, ω2, τ) =
g(ω1)g(ω2)
(iω1 +
)(−iω2 +
(iω1 +
)(−iω2 +
e(i(ω1−ω2)τ)
i(ω1 − ω2)
i(ω1 − ω2)
(iω2 −
− iω1)((
)2 + (ω2)2)
+iω2)τ − 1
(−iω1 −
+ iω2)((
)2 + (ω1)2)
−iω1)τ − 1
4(K − 2Kc
− iω1)(
+ iω2)
e−(Kc−K)τ − 1
C−1 = C
1 and the other modes vanish. Note that the
initial condition C1(ω, ω
′, 0) = 0 is satisfied and that the
time constants and frequencies which appear are every
possible way of pairing those from the tree level propa-
gator.
For illustrative purposes, Figure 2 b) shows the one
loop diagrams which contribute to C; these diagrams are
O(1/N2). We should note here the special role played by
the diagrams with initial state terms, in particular the
vertex proportional to ψ̃(θ, ω, t0)
2. This diagram evalu-
ates to exactly the same result as (44), with an additional
factor of −1/N . It serves to provide the proper normal-
ization for the two point function, which should go as
(N − 1)/N since the self-interaction (diagonal) terms are
not included. The other diagrams diverge faster as one
approaches criticality (K = Kc). They are of negligible
importance at small coupling but become increasingly
important near the onset of synchrony.
C. Order Parameter Fluctuations
We now compute the fluctuations in the order parame-
ter Z given in Eq. (2). The variance of Z (second moment
〈ZZ̄〉) is given by
〈ZZ̄〉 = 〈r2(t)〉
dωdω′dθdθ′〈η(ω, θ, t)η(ω′, θ′, t)〉ei(θ−θ
Using equation (19) in equation (46) gives
〈r2(t)〉 =
dωdω′dθdθ′C(x, x′, t)ei(θ−θ
since in the incoherent state ρ(x, t) = g(ω)/2π is inde-
pendent of θ, so that 〈Z〉 = 0. Hence
〈r2(τ)〉 = 4π2
dωdω′C−1(ω, ω
′, τ) +
which evaluates to [17]
〈r2(τ)〉 =
Λ1(s− s0)− 1
Λ1(s− s0)
× Res
Λ1(s)
esτ +
The time evolution of 〈r2〉 is then determined by the
poles of Λ1(s). As an example, consider f(θ) = sin θ and
g(ω) a Lorentz distribution (see Appendix B); we have
the result
〈r2(τ)〉 =
Kc −K
Kc −K
e−(Kc−K)τ (50)
where Kc = 2γ. Note that this diverges as τ → ∞ for
K = Kc. In the mean field limit N → ∞, 〈r
2〉 = 0 as
expected. As was shown in Ref [17], the tree level calcu-
lation adequately captures the fluctuations except near
the onset of synchrony (K = Kc). An advantage of the
field theoretic formalism is that it allows us to approach
even higher moments without needing to worry about
the moment hierarchy. In particular, for f(θ) = sin θ it
is straightforward to show that higher cumulants, such
as 〈(ZZ̄)2〉 − 〈ZZ̄〉2, must be zero at tree level because
of rotational invariance (more precisely, any cumulant of
Z higher than quadratic). The cumulants are given by
graphs which are connected. Vertex III produces two
lines with wave numbers n = ±1. Additionally, the IV
and V vertices impose a shift in wave number, whereas
Z and Z̄ project onto ±1. In order to calculate these
higher fluctuations it is necessary to go to the one loop
level. Note that this does not imply that the higher cu-
mulants of η are zero. Figure 2 b) gives the diagrams
for the correlation function at one loop. The one loop
calculation would also give a better estimate for 〈r2〉, es-
pecially nearer to criticality. The non-interacting distri-
bution is Gaussian, with non-Gaussian behavior growing
as one approaches criticality.
IV. LINEAR STABILITY AND MARGINAL
MODES
We analyze linear stability by convolving the linear
response or propagator with an initial perturbation. Al-
though we are free in our formalism to consider arbitrary
perturbations, we will consider two specific kinds for il-
lustrative purposes. In the first case we perturb only the
angular distribution. In the second case we consider a
perturbation which fixes one oscillator to be at a given
angle θ and frequency ω at a given time t. We first calcu-
late the results at tree level which reproduces the mean
field theory results of Ref. [24]. In particular, we arrive
at the same spectrum and Landau damping results of
Refs. [24] and [12]. We then define and calculate the op-
erator Γ (an extension of Γ0) to one loop order which
ultimately allows us to calculate the corrections to the
spectrum to order 1/N .
A. Mean field theory
The bare propagator which is the full linear response
for mean field theory is given by Eq. (40). The zeroes
of the operator Γ with respect to s specify the spectrum
of the linear response. There is a set of marginal modes
(continuous spectrum) along the imaginary axis spanning
an interval given by the support of g(ω). There is also a
set of discrete modes given by the zeros of the dielectric
function Λn(s) as was found in Ref. [24], aside from the
issue of the analytic continuation of Λn(s).
However, even though there are marginally stable
modes, the order parameter Z can still decay to zero due
to a generalized Landau damping effect as was shown in
Ref. [12]. Consider a generalization of the order param-
Zn(t) =
einθj , (51)
which represents Fourier modes of the density integrated
over all frequencies:
Zn(t) =
dθdωη(θ, ω, t)einθ (52)
and hence
〈Zn(t)〉 =
dθdωρ(θ, ω, t)einθ
dωρ(−n, ω, t). (53)
In the incoherent state, ρ(θ, ω, t) is independent of θ so
〈Zn〉 is zero for n > 0. The density response δρ(θ, ω, t)
to an initial perturbation δρ(θ, ω, 0) is given by
δρ(θ, ω, t) = N
P0(x, x
′, t)δρ(θ′, ω′, 0)dθ′dω′
Recall from the definition of the action (15), the prop-
agator operates on an initial condition defined by Nρ0.
The perturbed order parameter thus obeys
δ〈Zn(t)〉 =
dθdωδρ(θ, ω, t)einθ. (55)
We will show that for any initial condition involving a
smooth distribution in frequency and angle, δ〈Zn(t)〉 will
decay to zero. However, for non-smooth initial perturba-
tions, δ〈Zn(t)〉 will not decay to zero but will oscillate.
We first consider an initial perturbation of the form
δρ(θ, ω, 0) = g(ω)c(θ) (56)
where
c(θ)dθ = 0 (57)
Inserting into (54) yields
δρ(θ, ω, t) = N
P0(x, x
′, t)c(θ′)g(ω′)dθ′dω′
which is consistent with the perturbation considered in
Ref. [24]. Taking the Laplace transform of (58) gives
δρ̃n(ω, s) = Ncn
P̃0(n, ω, ω
′, s)g(ω′)dω′ (59)
Using the tree level propagator (40), we can show that
P̃0(n, ω, ω
′, s)g(ω′)dω′ =
s+ inω
Λn(s)
where Λn(s) is given in (39). Hence
δρ̃(n, ω, s) =
s+ inω
Λn(s)
From Eq. (61), we see that the continuous spectrum is
given by inω and the discrete spectrum by the zeros of
Λn(s). If we define bn(ω, t) = δρ(n, ω, t)/(Ncng(ω)) then
using (58) and (40) we can show
+ inω
bn(ω, t)
+ inKf(−n)
bn(ω1, t)g(ω1)dω1
δ(t) (62)
which is equivalent to the linearized perturbation equa-
tion derived by Strogatz and Mirollo [24], with the ex-
ception that (62) includes the effects of the initial config-
uration through the source term proportional to δ(t).
Inserting (61) into the the Laplace transform of (55)
yields
δ〈Z̃n(s)〉 = c−n
s− inω
Λ−n(s)
= c−n
Λ−n(s)− 1
−inKf(n)
Λ−n(s)
. (63)
For f(θ) = sin θ, f(±1) = ∓i/2 which leads to
δ〈Z̃1(s)〉 = c−1
Λ−1(s)− 1
Λ−1(s)
. (64)
We note that δZ1(s) is identical to what was calculated
in Ref. [12] in which it was shown that δZ1(t) → 0 as
t → ∞. Even in the presence of marginal modes, the
order parameter decays to zero through dephasing of the
oscillators. This dephasing effect is similar to Landau
damping in plasma physics.
We can see this explicitly for the case of the Lorentz
distribution
g(ω) =
γ2 + ω2
From (39) we can calculate
Λ±1(s) =
s+ γ − K
and Λn(s) = 1 for n 6= ±1. The zero of Λ±1(s) is at
s±1 = −(γ −K/2), which provides a critical coupling
Kc = 2γ (67)
above which the system begins to synchronize. The in-
coherent state is reached when K < Kc, which gives
s±1 < 0. Thus
δ〈Z±1〉 = c∓1e
−(γ−K
δ〈Zn6=±1〉 = c−ne
−|n|γτ (68)
Hence angular perturbations decay away in the order pa-
rameter.
Landau damping due to dephasing is sufficient to de-
scribe the relaxation of Zn(t) to zero for a smooth per-
turbation. However, for non-smooth perturbations, this
may not be true. Consider the linear response to a stim-
ulus consisting of perturbing a single oscillator to have
initial position θ0 and frequency ω0:
ρ0(θ, ω) =
(N − 1)
+ δ(θ − θ0)δ(ω − ω0)
and so the initial perturbation is
δρ0(θ, ω) =
+ δ(θ − θ0)δ(ω − ω0)
Inserting into (54) gives the time evolution of this initial
perturbation
δρ(θ, ω, t) = −
+ P (θ, ω, t; θ0, ω0, t
′) (71)
Substituting into (55) and taking the Lapace transform
gives
δ〈Z̃n(s)〉 = 2π
dωδρ̃−n(ω, s)
dωP̃ (−n, ω, ω0, s)
s− inω0
Λ−n(s)
There are therefore two modes in δZn(t), one which de-
cays due to dephasing (determined by the zero of Λ−n(s))
and one which oscillates at frequency ω0. Thus, for this
perturbation, the tree level prediction is that δ〈Z〉 is not
zero but oscillates. Inverse Laplace transforming (72),
gives the time dependence of the order parameter
δ〈Z1(t)〉 =
ω20 +
γ − K
+ ω20 −
e−iω0t
−iω0 + γ −
e−(γ−
The perturbed oscillator has phase θ0 − ω0t. It can al-
ways be located; no information is lost as the time evo-
lution progresses. Hence, for a single oscillator pertur-
bation, the tree level calculation predicts that the order
parameter will not decay to zero. In the next section, we
show that to the next order in the loop expansion, which
accounts for finite size effects, the marginal modes are
moved off of the imaginary axis stabilizing the incoher-
ent state, and the order parameter for a single oscillator
perturbation decays to zero.
B. Finite size effects
Let us define a generalization of the operator Γ0.
Γ · P =
δ(x− x′)δ(t− t′) (74)
where P without a subscript denotes the full propagator
and the operator Γ is the functional inverse of P . We
can estimate the effect of finite size on the stability of the
incoherent state by calculating the one loop correction to
the operator Γ. We will see that the one loop correction
FIG. 3: The diagrams contributing to the propagator at one
loop order, organized by topology. We consider d) to be of
different topology than c) because it is equivalent to a tree
level diagram with an additional factor of 1/N , due to the
initial state vertex.
FIG. 4: Diagrammatic equation for the propagator. The dou-
ble lines represent the summation of the entire series in 1/N
for the propagator.
produces the effect, among others, of adding a diffusion
operator to Γ, which is enough to stabilize the continuum
of marginal modes because the continuous spectrum is
pushed off the imaginary axis by an amount proportional
to the diffusion coefficient.
We calculate the correction to Γ to one loop order. The
propagator is represented by diagrams with one incom-
ing line and one outgoing line. There are four groups of
diagrams which contribute to the propagator at one loop
order. They are shown in Figure 3 and are labeled by
a), b), c), and d). Using these graphs to calculate the
propagator to order 1/N is not sufficient to demonstrate
the behavior of the spectrum to order 1/N . However,
we can use these graphs to construct an approximation
of Γ to order 1/N and derive the spectrum from this.
If we denote the full propagator (i.e. the entire series
in 1/N for P ) by a double line, we can approximate the
full propagator recursively by the diagrammatic equation
shown in Figure 4. The only terms which are neglected
in this relation are those which are from two loop and
higher graphs and therefore would contribute O(1/N2)
to Γ. Readers familiar with field theory will note that we
are simply calculating the two point proper vertex, which
is the inverse of the full propagator, to one loop order.
If we act on both sides of this equation with Γ0 (the
operator whose inverse is the tree level propagator), we
arrive at an equation of the form:
Γ0 · P =
δ(x− x′)δ(t− t′)− Γ1 · P (75)
where we have implicitly defined the one loop correction
to Γ, which we label Γ1. The action of Γ0 converts the
leftmost propagator in each diagram into a delta func-
tion, so that the delta function term in (75) arises from
the tree level propagator line and Γ1 is then comprised of
the loop portions of the remaining diagrams (the “ampu-
tated” graphs). We denote the contribution to Γ1 from
each group of one loop diagrams by Γ1r(θ, ω;φ, η; t− t
where r represents a, b, c, or d indicating the group of
diagrams in question. Γ1b = 0 because the derivative
coupling acts on the incoherent state, which is homoge-
neous in θ. The equation of motion for the one loop
propagator P1(x, x
′, t) then has the form
Γ0 · P1 + Γ1 · P1 =
δ(θ − θ′)δ(ω − ω′)δ(t− t′) (76)
where Γ0 · P1 is given by
Γ0 · P1 =
f(θ1 − θ)ρ(x1, t)dθ1dω1
P1(x, x
′, t− t′)
f(θ1 − θ)ρ(x, t)P1(x1, x
′, t− t′)dθ1dω1 (77)
Γ1 · P1 =
dt′′Γ1a(θ, ω;φ, η; t− t
′′)P1(φ, η, t
′′; θ′, ω′; t′)
dt′Γ1c(θ, ω;φ, η; t− t
′′)P1(φ, η, t
′′; θ′, ω′; t′)
dt′Γ1d(θ, ω;φ, η; t− t
′′)P1(φ, η, t
′′; θ′, ω′; t′) (78)
is the one loop contribution. The kernels Γ1a, Γ1c, and
Γ1d are explicitly computed in Appendix A.
The expressions for the one loop contribution to Γ are
rather complicated but several key features can be ex-
tracted, namely 1) the introduction of a diffusion opera-
tor, 2) a shift in the driving frequency, and 3) the addi-
tion of higher order harmonics to the coupling function
f . The diffusion operator has the effect of shifting the
marginal spectrum from the imaginary axis into the left
hand plane. The effect is that the finite size fluctuations
to order 1/N stabilize the incoherent state.
We can see these effects more easily by considering the
special case of f(θ) = sin θ and g(ω) being a Lorentz dis-
tribution (see Appendix B). The Fourier-Laplace trans-
formed equation of motion for the one loop propagator
has the form
α±1(s;ω)P̃1(n, ω, ω
′, s)
g(ω)(1 + β1(s;ω))
dνP̃1(n, ν, ω
′, s)
δ(ω − ω′), n = ±1, (79)
α±2(s;ω)P̃1(n, ω, ω
′, s)
dνβ±2(s;ω, ν)P̃1(n, ν, ω
′, s)
δ(ω − ω′), n = ±2, (80)
αn(s;ω)P̃1(n, ω, ω
′, s) =
δ(ω − ω′), |n| > 2
where
αn(s;ω) =
s+ inω +
s+ γ − K
+ i(m+ n)ω
n(2m+ n) + n(m+ n)
2γ −K
β±1(s;ω) = −
s+ γ − K
± 2iω
γ − K
2γ − K
2γ −K
s+ γ − K
β±2(s;ω, η) = −
s∓ iω
s∓ iω + γ − K
s± i(η − ω)
±iω + γ
γ − K
s+ γ − K
s+ 2γ −K
s+ γ − K
2γ −K
−γ + K
2γ −K
s+ γ − K
s+ γ − K
s+ 2γ − K
s+ 2γ −K
. (82)
We can solve for P̃1 using the same kind of self-
consistency computation that we used for P̃0. This pro-
duces an analogously defined dialectric function, Λ1n(s).
Λ1n(s) = 1−
g(ω)(1 + βn(s;ω))
αn(s;ω)
Stability of the incoherent state is determined by the
spectrum of the operator Γ. Analogous to tree level, the
continuous spectrum is given by the zeros of αn(s;ω) and
the discrete spectrum by the zeros of Λ1n(s) given by (83).
However, at one loop order, the expressions for P̃1 rep-
resent solutions to a coupled system of equations. Thus,
the poles of tree level are shifted, and there are also new
poles reflecting the interaction of the mean density with
the two-point correlation function. These poles will have
residues of O(1/N2) owing to their higher order nature.
For n = ±1,±2, we cannot solve for the poles exactly
but we can approximate the shift in the tree level spec-
trum by evaluating the loop correction at the value of the
tree level pole, s = ∓inω, which is equivalent to using
the “on-shell” condition in field theory. Since the higher
order modes will decay faster than the tree level modes,
this essentially amounts to ignoring short time scales and
is similar to the Bogoliubov approximation [25, 26]. The
remaining effective equation for P̃1 is now first order and,
consequently, we can consider the spectrum of the implic-
itly defined operator analogous to (41).
The continuous spectrum consists of all the zeros of the
function αn. We expect a term of the form inω+O(1/N)
because of the tree level continuous spectrum. This will
govern the behavior at large times. In this case, the “on-
shell” condition is equivalent to Taylor expanding the
loop correction via s = inω + O(1/N) and keeping only
terms which are O(1/N). This yields:
αn(s;ω) = s+ in(ω + δω) + n
2D (84)
where
δω = −
γ − K
4γ −K
2γ −K
is a frequency shift and
γ − K
is a diffusion coefficient. The frequency shift, which is
negative, serves to tighten the distribution around the
average frequency. The diffusion operator serves to damp
the modes which are marginal at tree level.
The discrete spectrum arises from the zeroes of Λ1n(s).
We can again approximate the shift in the tree level zero
by using the on-shell condition. This gives
Λ1n(s) = 1−
s+ in(ω + δω) + n2D
2γ − K
2γ −K
γ − K
+ inω
γ − K
− inω
We assume the shift will be small which allows us to write
the zero of Λ1n(s) as
s = −
δΛ1n(s0)
) (88)
where s0 = −
γ − K
and δΛ1n(s) is the O(1/N) correc-
tion to Λ1n(s). This results in
sn = −
2γ −K
γ − K
6γ −K
2γ −K
Away from criticality (K ≪ 2γ) or for large N , this cor-
rection is small.
We conclude this section by writing down an effective
equation of motion for the density function which incor-
porates the effect of fluctuations. Recall the first equa-
tion of the BBGKY hierarchy is
f(θ′ − θ)ρ(x, t)ρ(x′, t)dθ′dω′
f(θ′ − θ)C(x;x′, t)dθ′dω′(90)
This equation has a term (the “collision” integral) involv-
ing the 2-point correlation function, C, on the right hand
side. The equation determining the tree level propaga-
tor (31) is the linearization of the first BBGKY equation
with C considered to be zero. The one loop correction
to this equation (78) incorporates the effect of the cor-
relations on the linearization. The diagrams in Figure 3
provide the linearization of the collision term, where C
is considered as a functional of ρ. Using our one loop
calculation, we can propose an effective density equation
at one loop order
sin(θ′ − θ)ρ(θ′,Ω′, t)ρ(θ,Ω, t)dθ′dω′
= −K2(ω)
sin(2θ′ − 2θ)ρ(θ′,Ω′, t)ρ(θ,Ω, t)dθ′dω′ (91)
where Ω = ω + δω and D is given by Eq. (86). The field
ρ(θ,Ω, t) is now defined in terms of the shifted frequency
distribution
G(Ω) ≈ g(ω)(1−
) (92)
The new coupling constant K2(ω) is O(1/N) and is due
solely to the fluctuations. It arises from the term β±2
in the equation for P̃1. In fact, given the structure of
the diagrams, it is clear that for O(1/Nn) there will be a
new coupling, Kn+1, which corresponds to a sin[(n+1)θ]
term. In the language of field theory, all odd couplings
are generated under renormalization.
The generation of higher order couplings is especially
interesting in light of the results of Crawford and Davies
concerning the scaling of the density η beyond the onset
of synchronization, i.e. η − ρ0 ∼ (K − Kc)
β [31, 32].
Although our calculations pertain to the incoherent state,
the fact that the loop corrections generate higher order
couplings is a general feature of the bulk theory defined
by the action of Eq. (14) as well. Thus, we expect a
crossover from β = 1/2 to β = 1 behavior to occur as
N gets smaller. This is consistent with [31], wherein a
crossover manifested as the rate constant became smaller
than the externally applied diffusion. In our case, the
magnitude of the diffusion is governed by the distance to
criticality and the number of oscillators.
Our proposed effective equation (91) is not self-
consistent because we use the propagator to infer the
form of the mean field equation. Thus, we neglect non-
linear terms which may arise due to the loop corrections.
In addition, our calculation applies specifically to per-
turbations in the incoherent state. There are likely other
terms we are neglecting for both of these reasons. The
consistent approach would be to calculate the effective ac-
tion to one loop order and derive the equation for ρ from
that. This would involve essentially the same calculation
we have performed here, but for arbitrary ρ(θ, ω, t) (i.e.
we would need to solve (31) for the propagator in the
presence of an arbitrary mean).
V. NUMERICAL SIMULATIONS
We compare our analytical results to simulations of sin-
gle oscillator perturbations, since this provides a direct
measurement of the propagator per equation (71). We
perform simulations of N oscillators with f(θ) = sin θ.
We fix 2% of the oscillators at a specific angle (θ0 = 0)
and driving frequency (unless N = 10, in which case we
fix a single oscillator; the plots with N = 10 have been
rescaled to match the other data). The remaining oscilla-
tors are initially uniformly distributed over angle θ with
driving frequencies drawn from a Lorentz distribution.
We measure the real part of Z1(t). This measurement
allows us to observe the behavior of the modes which are
marginal at tree level.
Equation (73) gives the behavior of δZ1(t) with a single
oscillator fixed at θ0 and ω0 at time t = 0. Recall that
Z1 = 0 in the incoherent state, so that we expect δZ1 ≈
Z1. To tree level
Z1(t) =
ω20 +
γ − K
(γ(γ −
) + ω20 −
iω0)e
−iω0t − (−iω0 + γ −
e−(γ−
In other words, the initially fixed oscillator has phase θ0−ω0t; no information is lost as the time evolution progresses.
Incorporating the one loop computation gives
Z1(t) =
ω20 +
γ − K
(γ(γ −
) + ω20 −
iω0)e
i(ω+δω)t−Dt − (−iω0 + γ −
where we have ignored a term of amplitude O(1/N2); we are only considering the contributions coming from the poles
described in the previous section. With the one loop corrections taken into account, we see that Z1(t) relaxes back
to zero as t→ ∞. In the simulations, we compute the real part of Z1(t) with θ0 = 0. This gives
Re(Z1(t)) =
ω20 +
γ − K
γ(γ −
) + ω20
cos ((ω0 + δω)t) e
sin ((ω0 + δω)t) e
−Dt −
)es1t
The special case of ω0 = 0 and θ0 = 0 gives:
Z1(t) =
γ − K
γe−Dt −
The imaginary part vanishes so that Re(Z1(t)) = Z1(t).
We first compare our estimate of the diffusion coef-
ficient D given by (86) with the simulations. We plot
the measured decay constant D of Z1 compared to the
theoretical estimate of (96) for the long time behavior in
Figure 5. These data only include values of K = 0.3Kc
and K = 0.5Kc (Kc = 2γ = 0.1). Higher values of K did
not yield good fits due to the neglected contributions to
Z1. These decay constants are obtained via fitting the
time evolution of Z1 to an exponential for t > 200s. In
both cases they behave as 1/N for large N as predicted.
There is a consistent discrepancy likely due to round-
ing error after simulating for such a long period of time
(≈ 30000 time steps). This error appears as a small de-
gree of noise which further damps the response, hence the
decay constants appear slightly larger in Figure 5. This
effect can be seen in Figures 6 and 7 as well. For large
times, the simulation data consistently fall slightly under
10 100
0.00001
0.00010
0.00100
0.01000
K = 0.03
K = 0.05
FIG. 5: The large time (> 200s) decay constants with zero
driving frequency for the perturbed oscillator. Lines are the
predictions given by Eq. 86). Solid line and circles represent
K = 0.03. Dashed line and boxes represent K = 0.05.
the analytic prediction. Similarly, the data is noisier at
large times.
Figures 6 and 7 show the evolution of Z1(t) over time
along with the analytical predictions and the tree level
result for K = 0.3Kc and K = 0.5Kc respectively. For
K = 0.3Kc, the prediction works quite well, with perhaps
the beginning of a systematic deviation appearing atN =
10 and N = 50 (there is a slight initial overshoot followed
by an undershoot at larger times). This same deviation
is more pronounced for K = 0.5Kc, although the data
follow the prediction quite well nonetheless.
Consistent with our expectations from the loop expan-
sion, as we move closer to criticality, i.e. the onset of syn-
chronization, the results for K = 0.7Kc and K = 0.9Kc
do not fare as well. Figure 8 demonstrates a marked devi-
ation from the prediction. We have not shown analytical
results for the lower values of N because the deviation
is so severe. The same holds true for all the results for
K = 0.9Kc, so that we have just plotted the simula-
tion data in Figure 9. The general trend of approaching
the mean field result still holds. The primary feature to
take from these plots is that the fluctuations increase the
decay constant. The closer to criticality, the more impor-
tant the fluctuations and the faster the decay, hence the
systematic undershoot which grows as one nears critical-
ity. The fastest relaxation to the incoherent state appears
at high K for a given N and at low N for a given K, ei-
ther limit results in increased effects from fluctuations.
It would be necessary to carry the loop expansion to two
or more loops in order to obtain good matches with these
data.
In Figures 10 and 11, we plot the time evolution of
Z1(t) given that the favored oscillator has a driving fre-
quency of ω0 = 0.05. Note first that Z1(t) approaches
the tree level calculation as N → ∞. The amplitude
of the oscillation also shows the same deviation as the
ω0 = 0 data, namely that of a slight initial overshoot of
the one loop prediction followed by an undershoot. In
0 500 1000 1500
Tree Level
N = 10
N = 50
N=100
N = 500
N = 1000
FIG. 6: Z1(t) vs. t for various values of N and K = 0.3Kc.
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as
N → ∞ the curve approaches the tree level value. From top
to bottom: Black line represents tree level. X’s and violet line
represent N = 1000. Triangles and blue line represent N =
500. Diamonds and green line represent N = 100. Boxes and
red line represent N = 50. Circles and purple line represent
N = 10.
0 500 1000 1500
Tree Level
N = 10
N = 50
N = 100
N = 500
N = 1000
FIG. 7: Z1(t) vs. t for various values of N and K = 0.5Kc.
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as
N → ∞ the curve approaches the tree level value. Symbols
as in Figure 6.
addition to this, we can see an increasing frequency shift
as N → 0. The data, prediction, and mean field results
eventually become out of phase. For intermediate values
of N , one can see that the one loop correction follows
this shift, while for N = 10, the mean field, data, and
one loop results each have a different phase.
In the case of n > 2 it is easier to write down a complete
analytical solution for the time evolution. With ω0 = 0,
θ0 = 0 we have,
Zn(t) =
γ − K
γ − K
)e−(γ−
)t+Dn2t
0 500 1000 1500
Tree Level
N = 50
N = 100
N = 500
N = 1000
0 500 1000 1500
Tree Level
N = 10
N = 50
N = 100
N = 500
N = 1000
FIG. 8: Z1(t) vs. t for various values of N and K = 0.7Kc.
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as
N → ∞ the curve approaches the tree level value. Symbols
as in Figure 6.
This is compared with a simulation result in Figure 12.
We see the same general trends as the previous graphs.
At large N , the simulation follows the prediction quite
well. For small N , the simulation seems consistently
higher than the one loop prediction with K = 0.3Kc.
For K = 0.5Kc, the prediction is again sufficiently sin-
gular that we have not plotted N = 10. The deviation is
already apparent for N = 50.
VI. DISCUSSION
Using techniques from field theory, we have produced
a theory which captures the fluctuations and correlations
of the Kuramoto model of coupled oscillators. Although
we have used the Kuramoto model as an example system,
the methodology is readily extendible to other systems of
coupled oscillators, even those which are not interacting
via all-to-all couplings. Moreover, the methodology can
be readily applied to any system which obeys a conti-
nuity equation. We derive an action that describes the
dynamics of the Kuramoto model. The path integral de-
fined by this action constitutes an ensemble average over
the configurations of the system, i.e. the phases and driv-
ing frequencies of the oscillators. Because the dynamics
0 500 1000 1500
Tree Level
N = 10
N = 50
N = 100
N = 500
N = 1000
FIG. 9: Z1(t) vs. t for various values of N and K = 0.9Kc.
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as
N → ∞ the curve approaches the tree level value. Symbols
as in Figure 6.
of the model are deterministic, this is equivalent to an
ensemble average over initial phases and driving frequen-
cies. Using the loop expansion, we can compute moments
of the oscillator density function perturbatively with the
inverse system size as an expansion parameter. However,
it is important to point out that the loop expansion is
equivalent to an expansion in 1/N only because of the
all-to-all coupling. A local coupling will produce fluctu-
ations which do not vanish in the thermodynamic limit.
Our previous work in this direction developed a mo-
ment hierarchy analogous to the BBGKY hierarchy in
plasma physics. This paper fully encompasses that earlier
work. The equations of motion for the multi-oscillator
density functions derivable from the action are in fact the
equations of that BBGKY hierarchy. The calculation in
Ref. [17] is in the present context the tree level calculation
of the 2-point correlation function, given by the Feyn-
man graph in Figure 2a. With the BBGKY hierarchy,
the calculational approach involves arbitrary truncation
at some order, with no a priori knowledge of how this
approximation is related to the system size, N . Here we
show that this approximation is entirely equivalent to the
loop expansion approximation. Truncating the hierarchy
at the nth moment is equivalent to truncating the loop
expansion at the (n− l)th loop for the lth moment. The
one loop calculation is performed in the BBGKY context
by considering the linear response in the presence of the
2-point correlation function. This would produce a more
roundabout manner of arriving at our one loop linear
response. One should also compare our one loop calcula-
tion with the Direct-Interaction-Approximation of fluid
dynamics; the path integral approach in that context is
the Martin-Siggia-Rose formalism [27]. Another possible
equivalent means of approaching this problem is through
the Ito Calculus, treating the density as a stochastic vari-
able and developing a stochastic differential equation for
An important aspect of our theory is that it is di-
0 500 1000 1500
-0.03
-0.02
-0.01
Tree Level
One Loop
N = 10
0 500 1000 1500
-0.03
-0.02
-0.01
Tree Level
One Loop
N = 50
0 500 1000 1500
-0.03
-0.02
-0.01
Tree Level
One Loop
N = 100
0 500 1000 1500
-0.03
-0.02
-0.01
N = 500
One Loop
Tree Level
FIG. 10: As the previous figures, but with K = 0.3Kc and
ω0 = 0.05. Solid line represents tree level. Dashed blue line
represents the one loop calculation. Circles represent the sim-
ulation data.
0 500 1000 1500
Tree Level
One Loop
N = 10
0 500 1000 1500
-0.03
-0.02
-0.01
Tree Level
One Loop
N = 50
0 500 1000 1500
-0.03
-0.02
-0.01
Tree Level
One Loop
N = 100
0 500 1000 1500
-0.03
-0.02
-0.01
Tree Level
One Loop
N = 500
FIG. 11: As the previous figures, but with K = 0.5Kc and
ω0 = 0.05. Symbols as in Figure 10.
0 500
Tree Level
N = 10
N = 50
N = 100
N = 500
N = 1000
0 100 200 300 400 500
0.005
0.015
Tree Level
N = 50
N = 100
N = 500
N = 1000
FIG. 12: δZ3(t) vs. t for K = 0.3Kc (top) and K = 0.5Kc
(bottom). Symbols as in Figure 6.
rectly related to a Markov process derivable from the
Kuramoto equation. One can employ the standard Doi-
Peliti method for deriving an action from a Markov pro-
cess to arrive at the same theory. Although the Ku-
ramoto model is deterministic, the probability distribu-
tion evolves in a manner indistinguishable from a funda-
mentally random process. The stochasticity of the effec-
tive Markov process is due to the distribution of phases
and driving frequencies. In other words, it is a state-
ment about information available to us about the state
of the system. The incoherent state is a state of high en-
tropy. The single oscillator perturbation is one in which
we have gained a small amount of information about the
system and we ask a question concerning our knowledge
about future states. In the mean field limit for the sin-
gle oscillator perturbation, we always know where to find
the perturbed oscillator given a prescription of its initial
state. In the finite case, our ignorance of the positions
and driving frequencies of the other oscillators makes a
determination of its future location difficult. Eventually,
we lose all ability to locate the perturbed oscillator as it
interacts with the “heat bath” of the population. Fur-
thermore, this result should be time reversal invariant.
Just as we have no way of determining with accuracy
where to find the oscillator in the future, likewise we
have no means of determining where it has been at some
time in the past. To prove this statement in the con-
text of our theory would require an analysis of the “time
reversed” theory, obtained essentially by switching the
roles of ϕ̃ and ϕ. The relevant propagator for this time
reversed theory will be the solution of the linearization
of the adjoint of the mean field equation. Accordingly,
this adjoint theory will have loop corrections which will
damp the time reversed propagator as well.
It is important to point out that our formulation ac-
counts for the local stability of the incoherent state
ρ(θ, ω, t) = g(ω)/2π to linear perturbations along with
demonstrating the order parameters Zn approach zero.
In mean field theory, there is the possibility of quasiperi-
odic oscillations so that Zn = 0, i.e. the modes de-phase,
while the incoherent state is marginally stable and infor-
mation of the initial state is retained. Our work shows
that in a finite size system, this is does not happen; the
incoherent state is linearly stable.
We have considered exclusively the case of fluctua-
tions about the incoherent state below the critical point.
Above criticality, a fraction of the population synchro-
nizes. In this case, to analyze the fluctuations one may
need to employ a “low temperature” expansion in con-
trast to our “high temperature” treatment. In essence,
one separates the populations into locked and unlocked
oscillators and derives a perturbation expansion from the
locked action. At criticality, each term in the loop expan-
sion diverges. This is an indication that fluctuations at
all scales become relevant near the transition and thus
a renormalization group approach is suggested. Our for-
malism provides a natural basis for this approach.
In summary, we have provided a method for deriving
the statistics of theories defined via a Klimontovich, or
continuity, equation for a number density. This method
produces a consistent means for approximating arbitrary
multi-point functions. In the case of all-to-all coupling,
this approximation becomes a system size expansion. We
have demonstrated further that the system size correc-
tions are sufficient to render the incoherent state of the
Kuramoto model stable to perturbations.
ACKNOWLEDGMENTS
This research was supported by the Intramural Re-
search Program of NIH/NIDDK.
APPENDIX A: ONE LOOP CALCULATION OF THE PROPAGATOR
The loop correction Γ1a applied to P is given by
dt′Γ1a(θ, ω;φ, ν; t− t
′′)P (φ, ν, t′′; θ′, ω′; t′) = (A1)
dθ1dω1dθ
[f(θ′2 − θ) {P0(θ
2, t; θ
1, t1)P0(θ, ω, t; θ1, ω1, t1)
+ P0(θ
2, t; θ1, ω1, t1)P0(θ, ω, t; θ
1, t1)}]
[f(θ′1 − θ1) {ρ(θ
1, t1)P (θ1, ω1, t1; θ
′, ω′, t′) + ρ(θ1, ω1, t1)P (θ
1, t1; θ
′, ω′, t′)}]
This term arises from one vertex with a single incoming line and two outgoing lines and one vertex with two incoming
lines and a single outgoing line, hence the product of two tree level propagators, P0. We can represent Γ1a in
Fourier/Laplace space as
Γ1a(n, ω, ν, t− t
′) = (A2)
dω′1dω
2(2π)
n(m+ n)f(−m)g(ω′1)P
0(−m,ω′2, t;ω
′)P 0(m+ n, ω, t; ν, t′)
(−nm)f(m)g(ω′1)P
0(−m,ω′2, t;ω
′)P 0(m+ n, ω, t; ν, t′)
(−mn)f(m+ n)g(ω′1)P
0(−m,ω′2, t; ν, t
′)P 0(m+ n, ω, t;ω′1, t
n(m+ n)f(−n−m)g(ω′1)P
0(−m,ω′2, t; ν, t
′)P 0(m+ n, ω, t;ω′1, t
If f(θ) is odd then f(m) = −f(−m). Therefore,
Γ1a(n, ω, ν, t− t
′) = (A3)
dω′1dω
2(2π)
n(2m+ n)f(−m)g(ω′1)P
0(−m,ω′2, t;ω
′)P 0(m+ n, ω, t; ν, t′)
n(2m+ n)f(−n−m)g(ω′1)P
0(−m,ω′2, t; ν, t
′)P 0(m+ n, ω, t;ω′1, t
In this form it is easy to see the different channels which appear in the correction.
Evaluating this expression using the tree level propagator (40) gives us
Γ̃1a(n, ω, ν, s) = (A4)
(2π)2
n(2m+ n)f(−m)
imKf(m)
Λ−m(s1)
s1=sn
NP̃ 0(m+ n, ω; ν, s− sn)
n(2m+ n)f(−m− n)
Λ−m(s1)
s1=sn
(sn − imν)
(s− sn + i(m+ n)ω)
Λm+n(s− sn)
n(2m+ n)f(−m− n)
Λ−m(imν)
(s− imν + i(m+ n)ω)
Λm+n(s− imν)
The diagram Γ1c is given by
dt′Γ1c(n, θ, ω;φ, ν; t− t
′)P (φ, ν; θ′, ω′; t′) = (A6)
(2π)2
dθ′3dω
dθ′idω
idθidωidti
f(θ′3 − θ)
f(θ′1 − θ1)g(ω1)g(ω
{P0(θ
3, t; θ2, ω2, t2)P0(θ, ω, t; θ1, ω1, t1)
+ P0(θ
3, t; θ1, ω1, t1)P0(θ, ω, t; θ2, ω2, t2)}
f(θ′2 − θ2) {P0(θ
2, t2; θ
1, t1)P (θ2, ω2, t2; θ
′, ω′, t′)
+ P (θ′2, ω
2, t2; θ
′, ω′, t′)P0(θ2, ω2, t2; θ
1, t1)}] (A7)
In Fourier space, we can write:
Γ1c(n, ω, ν, t− t
′) = 2N2K3(2π)3
dω′3dω1dω
1dt1g(ω1)g(ω
dω′2(imn)(m+ n)f(−m− n)f(m)f(−m)P (n+m,ω
3, t; ν, t
′)P (−m,ω, t;ω1, t1)P (m,ω
′;ω′1, t1)
dω′2inm(m+ n)f(m)f(m)f(−m)P (−m,ω
3, t;ω1, t1)P (n+m,ω, t; ν, t
′)P (m,ω′2, t
′;ω′1, t1)
dω2inm(m+ n)f(−n−m)f(m)f(−n)P (n+m,ω
3, t;ω2, t
′)P (−m,ω, t;ω1, t1)P (m,ω2, t
′;ω′1, t1)
dω2inm(m+ n)f(m)f(m)f(−n)P (−m,ω
3, t;ω1, t1)P (n+m,ω, t;ω2, t
′)P (m,ω2, t
′;ω′1, t1)
and taking the Laplace transform:
Γ̃1c(n, ω, ν, s) =
2N2K3(2π)3
dω′3dω1dω
1ds1g(ω1)g(ω
dω′2(imn)(m+ n)f(−m− n)f(m)f(−m)P̃ (n+m,ω
3; ν, s− s1)P (−m,ω;ω1, s1)P (m)ω
1,−s1)
dω′2inm(m+ n)f(m)f(m)f(−m)P (−m,ω
3;ω1, s1)P (n+m,ω; ν, s− s1)P (m,ω
1,−s1)
dω2inm(m+ n)f(−n−m)f(m)f(−n)P (n+m,ω
3;ω2, s− s1)P (−m,ω;ω1, s1)P (m,ω2;ω
1,−s1)
dω2inm(m+ n)f(m)f(m)f(−n)P (−m,ω
3;ω1, s1)P (n+m,ω;ω2, s− s1)P (m,ω2;ω
1,−s1)
where the contour for s1 lies between 0 and 0 < s < −Re(sn). Performing the integrals, we have:
Γ̃1c(n, ω, ν, s) =
(imn)(m+ n)
f(−m− n)f(m)f(−m)
Λn+m(s− s1)
s− s1 + i(m+ n)ν
s1 − imω
Λ−m(s1)
imKf(−m)
Λm(−s1)
+ f(m)f(m)f(−m)
imKf(m)
Λ−m(s1)
(2πN)P (n+m,ω; ν, s− s1)
imKf(−m)
Λm(−s1)
dω2f(−n−m)f(m)f(−n)
Λn+m(s− s1)
s− s1 + i(m+ n)ω2
s1 − imω
Λ−m(s1)
g(ω2)
−s1 + imω2
Λm(−s1)
dω2f(m)f(m)f(−n)
imKf(m)
Λ−m(s1)
(2πN)P (n+m,ω;ω2, s− s1)
g(ω2)
−s1 + imω2
Λm(−s1)
(A10)
Finally, the diagram Γ1d is given by
dt′Γ1d(θ, ω;φ, ν; t− t
′)P (φ, ν; θ′, ω′; t′) = (A11)
(2π)2
dθ′3dω
dθ′idω
idθidωi
[f(θ′3 − θ)g(ω1)g(ω
{P0(θ
3, t; θ2, ω2, t2)P0(θ, ω, t; θ1, ω1, t0)
+ P0(θ
3, t; θ1, ω1, t0)P0(θ, ω, t; θ2, ω2, t2)}
f(θ′2 − θ2) {P0(θ
2, t2; θ
1, t0)P (θ2, ω2, t2; θ
′, ω′, t′)
+ P (θ′2, ω
2, t2; θ
′, ω′, t′)P0(θ2, ω2, t2; θ
1, t0)}] (A12)
Using the fact that
dω′dθ′P (θ, ω, t; θ′, ω′, t′)g(ω′) = g(ω)/N and
dθf(θ) = 0 we have
dt′Γ1d(θ, ω;φ, ν; t− t
′)P (φ, ν; θ′, ω′; t′) = (A13)
(2π)2
dθ′3dω
dθ′idω
idθidωi
[f(θ′3 − θ)g(ω1)g(ω
1)P0(θ
3, t; θ2, ω2, t2)
f(θ′2 − θ2) P (θ
2, t2; θ
′, ω′, t′)] (A14)
In Fourier-Laplace we have
Γ̃1d(n;ω, ν; s) =
inKf(−n)g(ω)
Λn(s)
(A15)
APPENDIX B: f(θ) = sin θ AND g(ω) LORENTZ
In order to both simplify the correction and to provide
a concrete example, we will specialize to the case that
g(ω) is a Lorentz distribution and f(θ) = sin θ. f(θ) =
sin θ is the traditional coupling for the Kuramoto model
(and has the advantage of being bounded in Fourier space
so that we avoid “ultraviolet” singularities) and using a
Lorentz frequency distribution yields analytical results.
The Lorentz frequency distribution is given by
g(ω) =
γ2 + ω2
From this we can calculate
Λ±1(s) =
s+ γ − K
and Λn(s) = 1 for n 6= ±1. The residue of the function
Λ±1(s) is
Λ−m(s1)
s1=sn
and the pole is at s±1 = −(γ − K/2). This provides a
critical coupling
Kc = 2γ (B4)
above which the system begins to synchronize. The in-
coherent state is reached when K < Kc, which gives
s± < 0. From f(θ) = sin θ we also have f(±1) = ∓i/2.
In this case, diagram a evaluates to
Γ̃1a(n, ω, ν, s) = (B5)
(2π)2
n(2m+ n)f(−m)
imKf(m)
NP̃ 0(m+ n, ω; ν, s+ γ −
n(2m+ n)f(−m− n)
− γ − imν
s+ γ − K
+ i(m+ n)ω
s+ 2γ − K
s+ 2γ −K
n(2m+ n)f(−m− n)
Λ−m(imν)
(s− imν + i(m+ n)ω)
s+ γ − imν
s+ γ − K
− imν
There is no contribution for n = 0 which we expect
from probability conservation. The terms proportional
to n(2m + n)f(−n − m) will always evaluate to 0, be-
cause f(n 6= ±1) = 0. Since the tree level propagator
contains such a term as well, we have the simplification
Γ̃1a(n, ω, ν, s) = (B7)
n(2m+ n)
δ(ω − ν)
s+ γ − K
+ i(m+ n)ω
For diagram c, we see that we can immediately ignore
the third term because nmf(−n)f(−m− n) is always 0.
We also see that the first term is only non-zero for n =
±2, and the last term only for n = ±1. After performing
the s1 integration, this leaves us with
Γ̃1c(n, ω, ν, s) =
n(m+ n)
δ(ω − ν)
s+ γ − K
+ i(n+m)ω
2γ −K
+ δn±1
s+ γ − K
+ 2inω
γ − K
+ inω
2γ − K
2γ −K
+ δn±2(−g(ω))
iω + γ − K
i(ν − ω)
iω + γ
γ − K
s+ γ − K
s+ 2γ −K
s+ γ − K
2γ −K
−γ + K
2γ −K
s+ γ − K
s+ γ − K
s+ 2γ − K
s+ 2γ −K
Γ1d is given simply by
Γ̃1d(±1;ω, ν; s) = −
s+ γ − K
Γ̃1d(n 6= ±1;ω, ν; s) = 0 (B9)
APPENDIX C: EQUIVALENT MARKOV
PROCESS
The action (15) can be derived by applying the Doi-
Peliti method to a Markov process equivalent to the Ku-
ramoto dynamics. Consider a two dimensional lattice L,
periodic in one dimension, with lattice constants aθ in
the periodic direction and aω in the other. (The radius
of the eventual cylinder is R.) The indices i and j will
be used for the frequency and periodic domains, respec-
tively. The oscillators obey an equation of the form:
θ̇i = v
~θ, ~ω
The indices on ~θ and ~ω run over the lattice points of the
periodic and frequency variables. The state of the system
is described by the number of oscillators ni,j at each site.
Given this, the fraction of oscillators found on the lattice
sites is governed by the following Master equation:
dP (~n, t)
ni,jP (~n, t)
vi,j−1
(ni,j−1 + 1)P
~nj , t
where the indices of the vector ~n run over the lattice
points, L, and ~nj (note the superscript) is equal to ~n
except for the jth and j − 1st components. At those
points we have n
i,j = ni,j − 1 and n
i,j−1 = ni,j−1 + 1.
The first term on the RHS represents the outward flux
of oscillators from the state with nj oscillators at each
periodic lattice point while the second term is the inward
flux due to oscillators “hopping” from j − 1 to j. There
is no flux in the other direction (ω); this lattice variable
simply serves to label each oscillator by its fundamental
frequency.
Consider a generalization of the Kuramoto model of
the form:
θ̇i = ωi +
f (θj − θi) (C3)
N is the total number of oscillators and we impose f(0) =
0. The velocity in equation (C1) now has the form:
vij = iaω +
i′,j′
f ([j′ − j]aθ)ni′,j′ (C4)
In the limit aω → 0 we have iaω = ω. Similarly, iaθ → θ.
The factor of ni′,j has been added because the sum must
cover all oscillators, and this factor describes the number
at each site. We also sum over all frequency sites i′. The
master equation (C2) now takes the form:
dP (~n, t)
i′,j′
f ([j′ − j]aθ)ni′,j′
×ni,jP (~n, t)
i′,j′
f ([j′ − j + 1]aθ)ni′,j′
×(ni,j−1 + 1)P
~nj , t
= −HP (~n, t) (C5)
The matrix H is the Hamiltonian. From this point, one
can develop an operator representation as in Doi-Peliti.
Using coherent states and taking the continuum and ther-
modynamic limits results in the action (15), after “shift-
ing” the field ϕ̃.
[1] A. T. Winfree, Journal of Theoretical Biology 16, 15
(1967).
[2] C. Liu, D. Weaver, S. Strogatz, and S. Reppert, Cell 91,
855 (1987).
[3] D. Golomb and D. Hansel, Neural Computation 12, 1095
(2000).
[4] G. B. Ermentrout and J. Rinzel, Am. J. Physiol 246,
R102 (1984).
[5] G. B. Ermentrout, Journal of Mathematical Biology 29,
571 (1991).
[6] T. J. Walker, Science 166, 891 (1969).
[7] J. Pantaleone, Physical Review D 58 (1998).
[8] S. Y. Kourtchatov, V. V. Likhanskii, and A. P. Npar-
tovich, Physical Review A 52 (1995).
[9] K. Wiesenfeld, P. Colet, and S. H. Strogatz, Physical
Review Letters 76 (1996).
[10] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu-
lence (Springer-Verlag, 1984).
[11] S. H. Strogatz, Physica D 143, 1 (2000).
[12] S. H. Strogatz, R. E. Mirollo, and P. C. Matthews, Phys-
ical Review Letters 68, 2730 (1992).
[13] R. E. Mirollo and S. H. Strogatz, Physica D 205, 249
(2005).
[14] R. E. Mirollo and S. H. Strogatz, The spectrum of the
partially locked state for the kuramoto model, eprint
arXiv:nlin/0702043.
[15] H. Daido, Journal of Statistical Physics 60, 753 (1990).
[16] H. Daido, Progress of Theoretical Physics 75, 1460
(1986).
[17] E. J. Hildebrand, M. A. Buice, and C. C. Chow, Physical
Review Letters 98 (2007).
[18] L. Peliti, Journal de Physique 46, 1469 (1985).
[19] M. Doi, Journal of Physics A: Mathematical and General
9, 1465 (1976).
[20] M. Doi, Journal of Physics A: Mathematical and General
9, 1479 (1976).
[21] H.-K. Janssen and U. C. Tauber, Annals of Physics 315,
147 (2005).
[22] J. Zinn-Justin, Quantum Field Theory and Critical Phe-
nomena (Oxford Science Publications, 2002), 4th ed.
[23] U. C. Tauber, Lecture Notes in Physics 716, 295 (2006).
[24] S. H. Strogatz and R. E. Mirollo, Journal of Statistical
Physics 63, 613 (1991).
[25] S. Ichimaru, Basic principles of Plasma Physics: A Sta-
tistical Approach (W.A. Benjamin Advanced Book Pro-
gram, 1973).
[26] D. R. Nicholson, Introduction to Plasma Theory (Krieger
Publishing Company, 1992).
[27] P. C. Martin, E. D. Siggia, and H. A. Rose, Physical
Review A 8, 423 (1973).
[28] A. V. Rangan and D. Cai, Physical Review Letters 96
(2006).
[29] J. M. Cornwall, R. Jackiw, and E. Tomboulis, Physical
Review D 10, 2428 (1974).
[30] R. E. Mirollo and S. H. Strogatz, Journal of Statistical
Physics 60, 245 (1990).
[31] J. D. Crawford, Physical Review Letters 74 (1995).
[32] J. D. Crawford and K. Davies, Physica D 125, 1 (1999).

  The incoherent state of the Kuramoto model of coupled oscillators exhibits
marginal modes in mean field theory. We demonstrate that corrections due to
finite size effects render these modes stable in the subcritical case, i.e.
when the population is not synchronous. This demonstration is facilitated by
the construction of a non-equilibrium statistical field theoretic formulation
of a generic model of coupled oscillators. This theory is consistent with
previous results. In the all-to-all case, the fluctuations in this theory are
due completely to finite size corrections, which can be calculated in an
expansion in 1/N, where N is the number of oscillators. The N -> infinity limit
of this theory is what is traditionally called mean field theory for the
Kuramoto model.

Introduction to Plasma Theory (Krieger
Publishing Company, 1992).
[27] P. C. Martin, E. D. Siggia, and H. A. Rose, Physical
Review A 8, 423 (1973).
[28] A. V. Rangan and D. Cai, Physical Review Letters 96
(2006).
[29] J. M. Cornwall, R. Jackiw, and E. Tomboulis, Physical
Review D 10, 2428 (1974).
[30] R. E. Mirollo and S. H. Strogatz, Journal of Statistical
Physics 60, 245 (1990).
[31] J. D. Crawford, Physical Review Letters 74 (1995).
[32] J. D. Crawford and K. Davies, Physica D 125, 1 (1999).
