# Coherence resonance in influencer networks

Ralf Tönjes1, Carlos E. Fiore2 and Tiago Pereira3,4

  * 1

Institute for Physics and Astronomy, University of Potsdam, Karl-Liebknecht-
Str. 24, 14476 Potsdam, Germany

  * 2

Instituto de Física, Universidade de São Paulo, São Paulo, Brazil

  * 3

Instituto de Ciências Matemáticas e Computação, Universidade de São Paulo, São
Carlos, São Paulo, Brazil

  * 4

Department of Mathematics, Imperial College London, SW7 2AZ, London, UK

## Abstract

Complex networks are abundant in nature and many share an important structural
property: they contain a few nodes that are abnormally highly connected
(hubs). Some of these hubs are called influencers because they couple strongly
to the network and play fundamental dynamical and structural roles.
Strikingly, despite the abundance of networks with influencers, little is
known about their response to stochastic forcing. Here, for oscillatory
dynamics on influencer networks, we show that subjecting influencers to an
optimal intensity of noise can result in enhanced network synchronization.
This new network dynamical effect, which we call coherence resonance in
influencer networks, emerges from a synergy between network structure and
stochasticity and is highly nonlinear, vanishing when the noise is too weak or
too strong. Our results reveal that the influencer backbone can sharply
increase the dynamical response in complex systems of coupled oscillators.

## Introduction

A central discovery in network science is that a small group of highly
connected hubs can couple to the network more strongly than their peers and
greatly influence the network behavior [1, 2, 3, 4, 5, 6]. Examples of network
influencers can be found in neuroscience (e.g., normal and aberrant synaptic
connectivity [7, 8, 9, 10]), political opinions (e.g., election blogging [11]
or social networks), and man-made scale-free networks (e.g., the internet
[1]). Surprisingly, the presence of such influencers makes synchronization of
deterministic network dynamics more difficult because networks with
influencers require stronger coupling than homogenous networks [12, 13];
indeed, in many situations, synchronization of influencer networks cannot be
achieved at all [14, 15]. This observation is all the more remarkable because
synchronization plays a fundamental role in regulating network function [16,
17] and is mediated predominantly through influencers [7, 18, 19]. This raises
a crucial question: why have many real-world networks evolved to contain
influencers when they appear to be detrimental to the network dynamics, at
least at face value?

While strong random fluctuations usually have a negative effect in complex
systems it has long been recognized that a small amount of noise can actually
improve the system response and its ability to process information. Known
mechanisms for such a constructive influence of noise are stochastic
resonance, coherence resonance and noise induced synchronization [20, 21, 22,
23, 24, 25, 26, 27, 28, 29, 30]. The term coherence resonance is used to
describe an optimal response of noise-induced oscillations without external
stimmulus in excitable cells [22]. It was observed in globally coupled systems
[23], in homogeneous networks [24, 25], in non-excitable systems near a Hopf
bifurcation [26] and two coupled oscillators [27]. The effects of coherence
resonance and its role in heterogeneous networks such as influencer networks
remains elusive.

In this work, we show that stochastic forcing of influencers can lead to an
optimal collective network response. Strikingly, introduction of noise
synergizes with the network structure to create collective oscillations that
become optimal at a given noise strength in the influencers. This phenomenon
emerges in two steps. First, the network acts as a nonlinear filter for the
stochastic influencer dynamics, and at an optimal noise strength, the
influencers induce synchronization in the nodes directly connected to them.
Second, different parts of the network develop macroscopic dynamics and
interact indirectly through the influencers. We develop an adiabatic theory to
uncover this macroscopic interaction law and show that it mediates the
emergence of global collective oscillations. When the noise in the influencers
is either too weak or too strong, the coupling vanishes. Interestingly, at a
macroscopic level, the interaction between different parts of the network can
be described by a hyper-graph.

We refer to a network where most nodes couple predominantly to a small number
of influencers as an influencer network, and refer to the remaining nodes as
followers (Figure 1). As generic oscillatory dynamics, we consider a network
of phase oscillators

$\dot{\vartheta}_{n}=\omega_{n}+\frac{\lambda_{n}}{\mu_{n}}\sum_{m=1}^{N}W_{nm}g(\vartheta_{m},\vartheta_{n})+\sqrt{2D_{n}}\xi_{n}.$
(1)

Here, $\omega_{n}$ is the natural frequency of node $n$, which couples with
strength $\lambda_{n}$ to the weighted mean of the coupling functions $g$ to
neighboring nodes $m$. Given a network weight coupling matrix $W_{nm}\geq 0$,
which is nonzero if node $n$ receives a link from node $m$, the intensity
$\mu_{n}$ is the total coupling weight received by the $n$th node. A table of
parameters and their function is provided in Methods. The coupling

$g(\vartheta_{m},\vartheta_{n})=\sin\left(\vartheta_{m}-\vartheta_{n}-\alpha\right)+c_{0}$
(2)

is generic for weakly coupled, nearly identical oscillators [31, 32]. The
parameter $\alpha$ is called phase frustration and the bias $c_{0}$ is due to
shear, an amplitude dependence of the frequency [33, 34]. The effect of shear
is a shift in the average frequency proportional to the coupling strength.
Phase equations with this form of coupling $g$ are known as the Kuramoto-
Sakaguchi model [35, 36] and are widely applied across scientific disciplines
[31, 32, 33, 34, 35, 36]. In addition, each term $\sqrt{2D_{n}}\xi_{n}$
denotes uncorrelated Gaussian white noise of strengths $D_{n}$. In many
studies, the coupling strength $\lambda_{n}$ to the local mean-field is chosen
to be uniform, in which case the coupling is called normalized. In real-world
and experimental systems, though, coupling may be heterogenous and hubs can
couple more strongly to the network [5, 18]. We model this coupling as

$\lambda_{n}=\left\\{\begin{array}[]{cl}\beta_{n}\lambda_{0}&\mbox{~{}for~{}influencers}\\\
\lambda_{0}&\mbox{~{}for~{}followers.}\end{array}\right.$ (3)

For simplicity, throughout this exposition we consider the coupling intensity
$\beta_{n}=\beta$, the noise strength $D_{n}=D$, and the frequency
$\omega_{n}=\omega$, to be identical for all influencers. For the followers we
assume a Lorentzian frequency distribution with mean frequency $\omega_{0}$
and width $\gamma_{0}$. Noise is of identical strength $D_{n}=D_{0}$ in all
followers. We denote $\Delta\omega=\omega-\omega_{0}$ the average gap in
natural frequencies between influencers and followers.

In Methods we show how Eq. (1) can be recast in terms of dimensionless
effective parameters shown in Table 1. These effective parameters, and in
particular the influencer effective noise strength $q=D/\Delta\Omega$, play
key roles in the collective dynamics of the system. The dynamical frequency
gap $\Delta\Omega/\lambda_{0}$ leads to a time scale separation between the
dynamics of the followers and the influencers. A coupling intensity $\beta$ of
comparable but smaller magnitude leads to an effective coupling strength
$\Lambda$ close to one for which the effect of coherence resonance is most
pronounced. We note that the dynamical frequency gap needs to be large in
units of $\lambda_{0}$, but it can be small in natural time units. In
Supplementary Note 1, we present an example of the transformation for
realistic parameters in Eq. (1) to effective parameters.

Parameter | Meaning | range
---|---|---
$\Delta\Omega/\lambda_{0}=\Delta\omega/\lambda_{0}+(\beta-1)c_{0}$ | dynamical frequency gap | $\Delta\Omega/\lambda_{0}\gg 1$
$\Lambda=\beta\lambda_{0}/\Delta\Omega$ | dimensionless coupling strength | $\Lambda<1$
$q=D/\Delta\Omega$ | influencer effective noise strength | $q=O(1)$
$D_{0}/\lambda_{0}$ | followers effective noise strength | $D_{0}/\lambda_{0}\ll 1$
$\gamma_{0}/\lambda_{0}$ | followers frequency heterogeneity | $\gamma_{0}/\lambda_{0}\ll 1$
Table 1: Effective dynamical parameters in influencer networks. We obtain
these parameters as described in Methods. These are key parameters in the
description of coherence resonance and the optimal noise strength in the
influencers.

We divide the followers into partitions $P_{\sigma}$ of nodes connected to the
same set of influencers. In Figure 1, we show an influencer network with two
influencers ($a$ and $b$) and three partitions of followers ($\sigma$=1,2,3),
which are connected to influencers $a$, $b$, or both (see additional examples
in Supplementary Note 2). To capture the collective dynamics in each partition
$\sigma$, we introduce the complex mean-fields

$Z_{\sigma}(t)=\frac{1}{|P_{\sigma}|}\sum_{n\in
P_{\sigma}}e^{i\vartheta_{n}(t)}$ (4)

The modulo of the complex mean-field $R_{\sigma}=|Z_{\sigma}|$ is the
partition order parameter; that is, $R_{\sigma}=0$ for incoherent, uniformly
distributed phases and $R_{\sigma}=1$ in full synchrony. Similarly, the global
mean-field $Z$ and order parameter $R$ are defined by summing over all
followers in the network.

Figure 1:  Coherence resonance in an influencer network.  Distribution of the
order parameter $R$ versus the effective diffusion $q$ in the influencers. (a)
Influencers $a$ and $b$ are hubs that couple strongly to the network, and all
other nodes are regarded as followers. Three distinct partitions of followers
are shown in red, blue, and green, which connect to influencer $a$, $b$, and
both $a$ and $b$, respectively. In our simulation, each partition has 300
followers. (b) Mean-field theory predicts that the interactions of the
partition mean-fields take place in a hyper-graph mediated by the coupling
functions $G_{a}$ and $G_{b}$ see Methods. (c,d) For each value of effective
noise strength $q$ in the influencers, we plot the density of the global order
parameter $R$ on a color scale from $0$ (white) to the maximum value (dark
blue). At an optimal noise strength, the mean value of the global order
parameter reaches a maximum, revealing the coherence resonance effect. In (c)
the dynamical frequency gap between influencers and followers
$\Delta\Omega/\lambda_{0}=18$ is moderate, whereas in (d)
$\Delta\Omega/\lambda_{0}=198$ is large. The solid red line in (c) is our
analytical prediction.

## Results

With deterministic influencers where $q=0$, and when $|\Lambda|<1$, the
influencers cannot frequency lock to the followers. Synchronization of the
followers through the influencer backbone is poor and counteracted by noise
and frequency heterogeneity in the followers. Our results show that by setting
a weak noise strength or frequency heterogeneity in the followers and by
changing the effective noise $q$ in the influencers, synchronization of the
whole network increases, reaches a maximum, and then decreases.

We numerically integrate our model Eq. (1) in dimensionless units (Table 1)
for the network with two influencers (as shown in Figure 1) with $300$
identical followers in each partition, and a small fixed noise strength in the
followers. By changing the noise strength in the influencers, we then obtain
the distribution of the order parameter $R$ as a function of $q$. After a
transient, the order parameter is independent of the initial conditions. At an
optimal noise strength, $R$ reaches its maximum (in expected value), as shown
in Figure 1 for $\Delta\Omega/\lambda_{0}=18$ (panel c) and for
$\Delta\Omega/\lambda_{0}=198$ (panel d). The solid line is a theoretical
prediction in the thermodynamic limit for heterogeneous followers using a
slow-fast approximation. Frequency heterogeneity and noise in the followers
have qualitatively and quantitatively the same desynchronizing effect. Optimal
synchronization of the whole network is predicted theoretically and achieved
in all simulations for an effective noise strength $q\approx 1$ in the
influencers, see details in Methods. Our mean-field analysis predicts that the
effect of coherence resonance is only observed for very small frequency
heterogeneity or noise in the followers, below a threshold that depends on
$\Lambda$ (Supplementary Note 3).

In Figure 2, we show the time series of the order parameter $R$ for two
complex and real-world networks. The upper row represents a scale-free network
and the lower row the directed neural network in the model organism
Caenorhabditis elegans. We assign the role of influencers to the $K$ most
strongly connected nodes and use a weighted connectivity matrix $W_{nm}=1$ for
all connections from or to an influencer, and $W_{nm}=0.01$ for all other
connections. For small effective noise $q$ in the influencers (Figure 2, $q$
weak), the order parameter fluctuates at a low level. When $q=1$ ($q$
optimal), the order parameter fluctuates around a value close to $1$,
revealing coherent collective oscillations. Finally, when $q$ is large ($q$
strong), the order parameter decreases again, revealing the loss of synchrony.
All parameters for the simulation and numerical scheme can be found in
Methods. In Supplementary Note 4, we show three additional examples of
coherence resonance in influencer networks; with 3 influencers, a random
network with $100$ influences, and a network of linked political blogs.

Figure 2:  Coherence resonance of the order parameter in different complex
networks. (a-c, e-g) show the time series of the order parameter $R$ for three
values of noise strength in the influencers for weak $q=0.1$ (a,e), optimal
$q=1$ (b,f), and strong q=10 (c,g). (d,h) show the corresponding networks with
(d) a scale-free network with exponent $2$ [1] and (h) C. elegans directed
neural network [3]. See Methods for further details. Additional examples can
be found in Supplementary Note 5.

### Stochastic forcing by a single influencer

Let us consider a single influencer. When its followers are asynchronous, the
sinusoidal contributions in the sum of the coupling functions for that
influencer average out, and the influencer phase is effectively decoupled from
the followers. The influencer is independent of the network and acts as a
common stochastic force on the followers connected to it. The additive noise
in the influencer enters the dynamics of its followers multiplicatively
through the coupling function. That is, the network acts as a nonlinear filter
for the noise in the influencers. We can show that the effective diffusion
constant of the integrated stochastic forcing in the followers, a proxy for
the noise strength, attains a maximum at an optimal effective noise strength

$q_{\rm opt}=1\,.$ (5)

We present the calculations in Methods and further details in Supplementary
Note 5. When oscillations are driven by common multiplicative noise, the
effect of noise-induced synchronization can be observed [28, 29, 30]. As the
common noise intensity is increased the oscillators synchronize faster. This
suggests that at this optimal value of $q$ the incoherent state will be most
unstable. However, behavior of the noise transfer does not explain the
synchronization between different partition mean-fields. Because this
synchronization requires studying macroscopic dynamics of $R$ far from zero,
our next goal is to uncover the interaction function between the mean-fields
in different partitions.

### Mean-field dynamics of partitions takes place in a hyper-network

For simplicity, we provide an analysis of the influencer network shown in
Figure 1. Mean-field equations for mixed repulsive and attractive coupling or
intra- and inter-partition interactions in the followers can be generalized
from these results. In Methods, we show that assuming large partition sizes
$|P_{\sigma}|\gg 1$, large dynamical frequency gap
$\Delta\Omega/\lambda_{0}\gg 1$, and noise free followers $D_{0}=0$ with
frequency heterogeneity $\gamma_{0}$, it is possible to derive averaged
dynamics of the partition mean-fields $Z_{\sigma}$ in an adiabatic
approximation. The resulting deterministic equations have the structure of a
hyper-network. In our example, the governing equations are

$\dot{Z}_{1}=F(Z_{1},G_{a}),\quad\dot{Z}_{2}=F\left(Z_{2},\frac{1}{2}G_{a}+\frac{1}{2}G_{b}\right)\quad\mbox{and}\quad\dot{Z}_{3}=F(Z_{3},G_{b})$
(6)

where the coupling functions are

$G_{a}=G(w_{a1}Z_{1}+w_{a2}Z_{2})\mbox{~{}~{} and
~{}~{}}G_{b}=G(w_{b2}Z_{2}+w_{b3}Z_{3}).$ (7)

$F$ describes a Riccati force (see Methods). The weights $w_{k\sigma}$ with
$\sum_{\sigma}w_{k\sigma}=1$ denote the relative size of partition $\sigma$
among all followers of influencer $k$. Note that, while in the microscopic
description the connections between nodes are pairwise, at the level of mean-
fields, edges represented by a coupling function $G_{k}$ for each influencer
can connect multiple partitions of followers. Thus, the mean-field interaction
between different parts of the network is described by a hyper-graph.

The interaction functions $G_{a}$ and $G_{b}$ can be determined analytically;
they depend on $\Lambda$ and $q$, are maximal at an optimal noise strength,
and vanish at critical values of $q$. That is, at weak or strong noise in the
influencers, the hyper-network interactions vanish, revealing the highly
nonlinear nature of the phenomenon. In particular, this means that the
macroscopic fields will not interact in the strong noise limit. We derive the
analytic expressions for the coupling functions $G_{k}(Z;\Lambda,q)$ in
Methods.

### Global synchronization and resonance

When influencers have equal parameters $q$ and $\Lambda$ the synchronization
manifold $Z_{\sigma}=Z$ is invariant under (6) and (7) and stable for phase-
attractive coupling. Hence, the macroscopic fields synchronize and we can
explain the global coherence resonance by restricting the analysis to this
invariant subspace. Our mean-field theory predicts both effects: the coherence
resonance of the partition order parameters and phase synchronization of the
partition mean-fields as shown in Methods. The solid line in Figure 1 (right)
is the stationary average order parameter predicted by our theory in the
infinite time-scale separation limit and with frequency heterogeneity
$\gamma_{0}/\lambda_{0}=0.02$ in the followers. The predicted values agree
with simulations of the finite size network, large dynamical frequency gap
$\Delta\Omega/\lambda_{0}=198$ and identical followers with noise
$D_{0}/\lambda_{0}=0.02$. In Supplementary Information, we provide two short
movies displaying synchronization of the network in Fig.1 at an optimal noise
strength in the two influencers.

## Discussion

We have found a new effect induced by a synergy between noise and network
structure to generate a transition towards a synchronization that would not be
possible in the absence of noise. The key element for this effect is the
existence of influencers – a group of hubs that couple strongly and connect
different parts of a network. Although deterministic network parameters
prevent synchronization, we show that an optimal noise strength in the
influencers can induce and mediate synchronization. The mechanism for this
coherence resonance in influencer networks is different from the known effect
of coherence resonance in homogeneous networks with excitatory dynamics, where
noise simply excites oscillations [23, 24, 25]. At the macroscopic level, the
interaction between different parts of the network is indirect and takes place
on an emerging hyper-network, thus changing the interaction structure from the
microscopic level. Such higher order interactions have previously been
conjectured and reported in neuronal data recordings [37]. Our findings
suggest that the emergent order in complex systems could be controlled by
regulating the noise in only a few key nodes.

## Methods

### Canonical form

To bring the Eq. (1) into a dimensionless form with effective parameters given
in Table 1, we change the time scale to units of $1/\lambda_{0}$ and add the
frequency shift from the bias $c_{0}$ in the coupling function to the natural
frequencies of the oscillators, i.e.
$\omega_{n}\mapsto\omega_{n}+\lambda_{n}c_{0}$ and
$g(\vartheta_{m},\vartheta_{n})\to\sin(\vartheta_{m}-\vartheta_{n}-\alpha)$.
The difference between the average follower frequency and the frequency of an
influencer, both including the frequency shift from the coupling bias, is the
dynamical frequency gap $\Delta\Omega/\lambda_{0}$ (in units of
$\lambda_{0}$). Observing the invariance of the phase equations under a global
phase shift, i.e. $\vartheta\to\vartheta-\omega_{0}t$, we can go into a co-
rotating reference frame where the average follower frequency is zero. The
deviations of the follower frequencies from their mean frequency $\omega_{0}$
may be written as $\gamma_{0}\nu_{n}$ where the $\nu_{n}$ are taken from some
standard distribution with mean zero and the factor $\gamma_{0}\geq 0$
characterizes the frequency heterogeneity. Then the phase equations for the
followers in the new time units and co-rotating reference frame are

$\dot{\vartheta}_{n}=\frac{\gamma_{0}}{\lambda_{0}}\nu_{n}+\frac{\lambda_{n}}{\lambda_{0}}\frac{1}{\mu_{n}}\sum_{k}W_{nk}\sin(\vartheta_{k}-\vartheta_{n}-\alpha)+\sqrt{2D_{n}/\lambda_{0}}\tilde{\xi}_{n}$
(8)

and for the influencers with phases $\psi_{k}$

$\dot{\psi}_{k}=\frac{\Delta\Omega_{k}}{\lambda_{0}}\left(1+\Lambda_{k}\frac{1}{\mu_{k}}\sum_{m}W_{km}\sin(\vartheta_{m}-\psi_{k}-\alpha)\right)+\sqrt{2q_{k}\frac{\Delta\Omega_{k}}{\lambda_{0}}}\tilde{\xi}_{k}.$
(9)

Here, $\Lambda_{k}=\lambda_{k}/\Delta\Omega_{k}$ is the ratio between the
coupling strength and the frequency of the influencer. In the noise free case,
phase locking is only possible for $|\Lambda_{k}|>1$. Changing $\Lambda$ can
lead to a discontinuous, explosive synchronization [12, 38]. The terms
$\tilde{\xi}_{m}$ are independent white noise with
$\langle\tilde{\xi}_{m}(t)\tilde{\xi}_{n}(t^{\prime})\rangle=\delta_{mn}\delta(t-t^{\prime})$
in the new units of time and $q_{k}$ is the effective noise strength in the
influencers on the fast time scale $\Delta\Omega_{k}/\lambda_{0}$. In
Supplementary Note 1, we provide examples of such rescaling.

### Parameters and their meaning

In Table 2, we present the main parameters that naturally appear in the phase
model Eq. (1) and give rise to effective parameters, as shown in Table 1 in
the main text. The main parameters in our mean-field analysis are shown in
Table 3.

Parameter | Meaning
---|---
$\omega_{n}$ | isolated frequency of the $n$th oscillator;
| set as $\omega_{n}=\omega_{0}+\gamma_{0}\nu_{n}$ for followers and $\omega$
for influencers
$\omega_{0}$ | mean frequency of the followers
$\gamma_{0}\nu_{n}$ | frequency deviation $\omega_{n}-\omega_{0}$ of the $n$th follower
$\gamma_{0}$ | scale parameter of follower frequency distribution
$\Delta\omega$ | gap $(\omega-\omega_{0})$ between influencer and average follower frequency
$W_{nm}$ | nonnegative matrix of connection weights
$\mu_{n}$ | connection intensity ( $\mu_{n}=\sum_{m}W_{nm}$)
$\lambda_{n}$ | coupling strength of the $n$th oscillator;
|  $\lambda_{0}$ for followers and $\beta\lambda_{0}$ for influencers
$\beta$ | coupling intensity for influencers
$D_{n}$ | noise strength set as $D$ for influencers and $D_{0}$ for followers
$\alpha$ | phase frustration in the coupling function $g$
$c_{0}$ | shear parameter in the coupling function $g$
Table 2: Parameters in the model presented in Eq. (1) of the main manuscript. Parameter | Meaning
---|---
$P_{\sigma}$ | follower partitions according to the influencers they connect to
$I_{\sigma}$ | set of influencers of a partition $\sigma$
$Z_{\sigma}$ | complex mean-field of partition $\sigma$ (order parameter $R_{\sigma}=|Z_{\sigma}|$)
$G_{k}$ | coupling function between mean-fields mediated by influencer $k$
$w_{k\sigma}$ | relative size of partition $\sigma$ among the followers of influencer $k$
$F$ | Ricatti vector field see Eq. (15)
$h_{\sigma}$, $h_{k}$ | forces on oscillators in partition $\sigma$ and on influencer $k$
$H_{\sigma}$ | average of $h_{\sigma}$ obtained from adiabatic mean-field approximation
Table 3: Parameters of the mean-field analysis presented in Eq. (6).

### Simulations and parameter values

In our analysis and our simulations, we use the transformed, dimensionless
canonical form (8) and (9) of the phase equations (1). Existing connections in
the network from and to the influencers are given the weight $W_{mn}=1$ and
connections between followers $W_{nm}=0.01$. The followers couple to their
neighbors with strength $\lambda_{n}=\lambda_{0}$ and influencers with
strength $\lambda_{k}=\beta\lambda_{0}$. The phase frustration in the coupling
function is set to $\alpha=-0.1$. The frequency deviations $\nu_{n}$ of the
followers are drawn from a Cauchy distribution $p(\nu)=1/\pi(\nu^{2}+1)$ and
multiplied by $\gamma_{0}/\lambda_{0}$. Thus, frequency heterogeneity and
noise strength in the followers are given by $\gamma_{0}/\lambda_{0}$ and
$D_{0}/\lambda_{0}$, respectively. We chose $\beta_{k}=\beta$,
$\lambda_{k}=\beta\lambda_{0}$ and $D_{k}=D$ for all influencers, so that
$\Delta\Omega/\lambda_{0}$, $\Lambda$ and $q$ are identical for all
influencers. We integrate the Langevin equations of the phases with an Euler-
Maruyama scheme and small time steps $dt=5\cdot 10^{-4}$ because of the large
time scale separation $\Delta\Omega/\lambda_{0}\gg 1$.

The parameters of Figure 1 are as follows: The network structure is a pure
influencer network without connections between followers or between
influencers. We simulate 300 identical followers $\gamma_{0}/\lambda_{0}=0$ in
each of the three partitions with small independent noise
$D_{0}/\lambda_{0}=0.02$. In the lower left panel we have $\beta=10$, a
dynamical frequency gap of $\Delta\Omega/\lambda_{0}=18$ and an effective
influencer coupling strength $\Lambda=10/18$. In the lower right panel
$\beta=100$, $\Delta\Omega/\lambda_{0}=198$ and $\Lambda=100/198$. For each
value of the effective influencer noise strength $q=D/\Delta\Omega$ we record
a histogram of the order parameter over $T=10^{4}$ time units, which is much
longer than the relaxation time of $R$. The theoretical prediction, the solid
line in the lower right panel, is for noiseless followers $D_{0}=0$ and
$\gamma_{0}/\lambda_{0}=0.02$.

Parameters of Figure 2 are as follows: For the C. elegans directed neuronal
network [3], we choose the top $K=15$ out-degree nodes as influencers. All
nodes with zero in-degree have been removed, resulting in a network with
$N=268$ nodes. Connections between followers are given the weight
$W_{nm}=0.01$. We simulate identical followers $\gamma_{0}/\lambda_{0}=0$ with
small independent noise $D_{0}/\lambda_{0}=0.02$. The dynamical frequency gap
between followers and influencers is $\Delta\Omega/\lambda_{0}=18$ and the
effective coupling strength in the influencers is $\Lambda=10/18$. Shown are
three time series of the network order parameter for small ($q=0.1$), optimal
($q=1$), and large ($q=10$) noise strength in the influencers. The undirected
scale-free network with exponent 2 is the largest connected component of a
network generated via a configurational algorithm [1] without self loops or
double edges. We chose the top 5 degree nodes as influencers. The other
parameters are the same as in the C. elegans neuronal network.

### Mean-field dynamics in influencer networks

We have developed a mean-field theory for undirected influencer networks with
connections exclusively between influencers and followers, as shown in Figure
1. This theory can be generalized to more complex configurations, heterogenous
influencers, directed, attractive, or repulsive coupling between followers and
influencers, within partitions or between different partitions. While these
generalizations may lead to more complex dynamic behavior, the mechanism for
the coherence resonance is apparent in the simplest model.

We consider the network as a union of a set $P$ of followers and a set $I$ of
influencers. The nodes $n$ connected to an influencer $k$ are elements $n\in
P_{k}$ of the periphery of the influencer $k$. Intersections of the sets
$P_{k}$ form equivalence classes or partitions $P_{\sigma}$ of followers that
are connected to the same subsets $I_{\sigma}$ of influencers such as in Fig.1
all followers connected to influencer $a$ or $b$ or to both influencers. The
phases of the oscillators are encoded as complex variables
$z_{n}=\exp(i\vartheta_{n})$ for the followers and $z_{k}=\exp(i\psi_{k})$ for
the influencers. The dynamics can be formulated in terms of partition averages
and averages over the influencers of these partitions

$\displaystyle Z_{\sigma}$ $\displaystyle=$
$\displaystyle\frac{1}{|P_{\sigma}|}\sum_{n\in P_{\sigma}}z_{n}$ (10)
$\displaystyle h_{\sigma}$ $\displaystyle=$
$\displaystyle\frac{e^{-i\alpha}}{2i}\frac{1}{|I_{\sigma}|}\sum_{k\in
I_{\sigma}}z_{k}$ (11) $\displaystyle h_{k}$ $\displaystyle=$
$\displaystyle\frac{e^{-i\alpha}}{2i}\frac{1}{|P_{k}|}\sum_{n\in
P_{k}}z_{n}=\frac{e^{-i\alpha}}{2i}\sum_{\sigma}w_{k\sigma}Z_{\sigma}.$ (12)

Here, $h_{\sigma}$ and $h_{k}$ are the forces acting on the followers in
partition $\sigma$ and on the influencer $k$. The weight $w_{k\sigma}$ is the
relative size of partition $\sigma$ in the periphery of an influencer $k$;
that is, $w_{k\sigma}=|P_{\sigma}|/|P_{k}|$ when $P_{\sigma}\subseteq P_{k}$
or $w_{k\sigma}=0$ otherwise. Thus, the phase dynamics (8) and (9) can be
written in complex form as

$\displaystyle\dot{z}_{n}$ $\displaystyle=$ $\displaystyle
iz_{n}\left(\bar{h}_{\sigma}z_{n}+\frac{\gamma_{0}}{\lambda_{0}}\nu_{n}+h_{\sigma}\bar{z}_{n}\right)+iz_{n}\sqrt{2D_{0}/\lambda_{0}}\xi_{n}(t),\qquad
n\in P_{\sigma}$ (13) $\displaystyle\dot{z}_{k}$ $\displaystyle=$
$\displaystyle
iz_{k}\frac{\Delta\Omega_{k}}{\lambda_{0}}\left(\Lambda_{k}\bar{h}_{k}z_{k}+1+\Lambda_{k}h_{k}\bar{z}_{k}\right)+iz_{k}\sqrt{2D_{k}/\lambda_{0}}\xi_{k}(t)$
(14)

The first reduction of model complexity is via the Ott-Antonsen approach [39]
for followers without Gaussian white noise but frequency heterogeneity
$\gamma_{0}/\lambda_{0}$ with Cauchy-distributed frequency deviations
$\nu_{n}$. In the thermodynamic limit $|P_{\sigma}|\to\infty$ (keeping the
ratios $w_{k\sigma}$ of the partition sizes constant) there exists a globally
attractive invariant manifold on which the mean-fields $Z_{\sigma}$ evolve by
a complex Riccati equation as

$\dot{Z}_{\sigma}=i\left(\bar{h}_{\sigma}Z_{\sigma}^{2}+i\frac{\gamma_{0}}{\lambda_{0}}Z_{\sigma}+h_{\sigma}\right)=F(Z_{\sigma},h_{\sigma}).$
(15)

For large partition sizes, Eqs. (10-15) provide a good description of the
system dynamics, including an accurate description of the fluctuations of the
mean-fields (see Supplementary Note 6). The effect of small noise
$D_{0}/\lambda_{0}$ in the followers is comparable to the effect of frequency
heterogeneity $\gamma_{0}/\lambda_{0}$. For small white noise, the Ott-
Antonsen manifold is no longer invariant but one can derive a hierarchy of
corrections to the dynamics (15) in increasing orders of the noise strength.
To the zeroth order the effects of frequency heterogeneity and noise are
identical [40]. In fact, if the noise in the followers is white Cauchy noise,
the equivalence of noise and frequency heterogeneity is exact [41].

### Slow-fast dynamics

If there is a large dynamical frequency gap $\Delta\Omega_{k}/\lambda_{0}\gg
1$ between the followers and an influencer, oscillators in the follower group
experience an average force from the fast influencer. Conversely, if the
followers are desynchronized, the mean-field of the followers vanishes and the
influencer phases perform a drift diffusion process on the circle

$\dot{\psi}_{k}=\frac{\Delta\Omega_{k}}{\lambda_{0}}+\sqrt{2D_{k}/\lambda_{0}}\tilde{\xi}_{k}(t)$
(16)

whereas the followers connected to only that influencer experience a
stochastic forcing by the influencer phase

$\dot{\vartheta}_{n}=\frac{\gamma_{0}}{\lambda_{0}}\nu_{n}+\sin\left(\psi_{k}-\vartheta_{n}-\alpha\right)+\sqrt{2D_{0}/\lambda_{0}}\tilde{\xi}_{n}(t).$
(17)

This forcing is multiplicative since
$\sin\left(\psi_{k}-\vartheta_{n}-\alpha\right)=s_{k}\cos\vartheta_{n}-c_{k}\sin\vartheta_{n}$
with two uncorrelated but not independent random processes
$s_{k}(t)=\sin(\psi_{k}-\alpha)$ and $c_{k}(t)=\cos(\psi_{k}-\alpha)$. The
diffusion constants $D_{s}$ and $D_{c}$ for the integrated stochastic forces
quantify an effective noise strength and can be calculated as the integral of
the respective autocorrelation functions [42]

$\frac{D_{s}}{\lambda_{0}}=\frac{D_{c}}{\lambda_{0}}=\frac{\lambda_{0}}{2}\frac{D_{k}}{\Delta\Omega_{k}^{2}+D_{k}^{2}}=\frac{\lambda_{0}}{2\Delta\Omega_{k}}\frac{q_{k}}{1+q_{k}^{2}}.$
(18)

We present the details in Supplementary Note 5. By changing the noise strength
$D_{k}$, the effective noise strengths $D_{s}$ and $D_{c}$ have a maximum at
$D_{k}=\Delta\Omega_{k}$ or $q_{k}=D_{k}/\Delta\Omega_{k}=1$. At this noise
value, and for incoherent followers, the effect of noise-induced
synchronization [29] is expected to be strongest. As the amplitudes of $s_{k}$
and $c_{k}$ are bounded, when $D_{k}$ or the time scale separation
$\Delta\Omega/\lambda_{0}$ are further increased the effective noise strengths
go to zero.

For $\Delta\Omega/\lambda_{0}\gg 1$, the system has slow and fast dynamics and
we can replace the influencer phases $z_{k}$ contributing to the force fields
$h_{\sigma}$ (11) in each partition $\sigma$ by the expected values $G_{k}$ of
$z_{k}$ subject to Langevin equation (14). On the fast time scale of the
influencers, the fields $h_{k}$ are changing very slowly and can assumed to be
constant for the calculation of the $G_{k}$. In this averaged dynamics, the
influencers create an average force $H_{\sigma}$ that follows the partition
mean-fields adiabatically. The slow dynamics of the partition mean-fields is
thus given as

$\displaystyle\dot{Z}_{\sigma}$ $\displaystyle=$ $\displaystyle
i\left(\bar{H}_{\sigma}Z_{\sigma}^{2}+i\frac{\gamma_{0}}{\lambda_{0}}Z_{\sigma}+H_{\sigma}\right)=F(Z_{\sigma},H_{\sigma})$
(19) $\displaystyle H_{\sigma}$ $\displaystyle=$
$\displaystyle\frac{e^{-i\alpha}}{2i}\frac{1}{|I_{\sigma}|}\sum_{k\in
I_{\sigma}}G_{k}\left(\sum_{\sigma^{\prime}}w_{k\sigma^{\prime}}Z_{\sigma^{\prime}};~{}\Lambda_{k},q_{k}\right).$
(20)

This corresponds to a hyper-graph with partitions $\sigma$ as nodes and
coupling functions $G_{k}$ for each edge $k$ of the hyper-graph. General
setups can be considered as well, with intra and inter-partition coupling and
connections between influencers. The absence of such connections shows that
the synchronization is indeed a noise-induced effect.

### Mean-field of the fast influencers

The Langevin equation (14) for $z_{k}$ with constant fields $h_{k}$ is indeed
a complex formulation of the noisy Adler equation [42]

$\dot{\psi}=\frac{\Delta\Omega}{\lambda_{0}}\left(1+2\Lambda\textrm{Im}\left[he^{-i\psi}\right]\right)+\sqrt{2D/\lambda_{0}}\xi(t).$
(21)

The expected value $G$ of $z=\exp(i\psi)$ is the first circular moment of the
stationary distribution which has an expression as a continued fraction [42]
and evaluates to a ratio of confluent hypergeometric limit functions
${}_{0}F_{1}(o,x)$ [43]. It can be derived from the Fokker-Planck equation
noting that the Fourier modes $p_{k}=\langle\exp(ik\psi)\rangle$ of the
stationary distribution $p^{st}(\psi)$ are in a tridiagonal recurrence
relation

$0=ik\Delta\Omega\left(\Lambda\bar{h}p_{k+1}+p_{k}+\Lambda
hp_{k-1}\right)-Dk^{2}p_{k}$ (22)

which is solved by a continued fraction. Defining $q=D/\Delta\Omega$ and

$s=\frac{|h|}{ih},\quad o=1-i\frac{1}{q},\quad\textrm{and}\quad
x=\frac{q}{\Lambda|h|}$ (23)

we have $G=p_{1}$ and

$G=\frac{1}{o+\frac{x^{-2}}{(o+1)+\frac{x^{-2}}{(o+2)+\dots}}}\frac{1}{sx}=\frac{{}_{0}F_{1}(o+1,x^{-2})}{{}_{0}F_{1}(o,x^{-2})}\frac{1}{sox}.$
(24)

### Synchronization manifold and prediction of order parameter

If all influencers have the same effective noise strength
$q=q_{k}=\frac{D_{k}}{\Delta\Omega_{k}}$ and the same effective coupling
strength $\Lambda=\Lambda_{k}=\frac{\beta\lambda_{0}}{\Delta\Omega_{k}}$, the
synchronization manifold where all partitions have identical mean-fields
$Z_{\sigma}=Z=R^{i\Theta}$ is invariant under the averaged dynamics (19) and
(20) on the hyper-graph and we can write

$\displaystyle\dot{Z}$ $\displaystyle=$ $\displaystyle
i\left(\bar{H}Z^{2}+i\frac{\gamma_{0}}{\lambda_{0}}Z+H\right)$ (25)
$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{e^{-i\alpha}}{2i}G\left(Z;~{}\Lambda,q\right)$ (26)

where $G(Z;~{}\Lambda,q)$ is defined as (24) with

$h=\frac{e^{-i\alpha}}{2i}Z.$ (27)

In particular, because of rotational symmetry, the dynamics of the amplitude
$R=|Z|$ does not depend on the angle $\Theta$ of the mean-field

$\dot{R}=\textrm{Re}\left[\frac{e^{-i\alpha}}{2}G(R;~{}\Lambda,q)\right]\left(1-R^{2}\right)-\frac{\gamma_{0}}{\lambda_{0}}R.$
(28)

If the synchronization manifold is stable, the stable fixed points of this
dynamics where $\dot{R}=0$ approximate the average order parameter over all
followers. From (28) we find that the level sets of the right-hand side of

$\frac{\gamma_{0}}{\lambda_{0}}=\textrm{Re}\left[\frac{e^{-i\alpha}}{2}G(R;~{}\Lambda,q)\right]\frac{1-R^{2}}{R}$
(29)

determine this average order parameter $R$ for any given
$\gamma_{0}/\lambda_{0}$ implicitly. We show this prediction for
$\gamma_{0}/\lambda_{0}=0.02$ and $\Lambda=0.51$ as a solid line in the lower
right panel in Fig.1. Further resonance curves and maxima of $R$ for different
heterogeneities $\gamma_{0}$ and different $\Lambda$ can be found in
Supplementary Note 3.

Acknowledgements We would like to thank C. Sagastizábal, D. Eroglu, S. van
Strien, D. Turaev, J. Lamb, and A. Pikovsky for enlightening discussions. This
work was supported in parts by the DFG and FAPESP through the IRTG 1740/TRP
2015/50122-0, by the Center for Research in Mathematics Applied to Industry
(FAPESP Cemeai grant 2013/07375-0) and grants 2015/04451-2, by the Royal
Society London, CNPq grant 302836/2018-7, and by the Serrapilheira Institute
(Grant No. Serra-1709-16124).

## Publication

This work was published under Creative Commons CC BY license at

https://doi.org/10.1038/s41467-020-20441-4

together with Supplementary Notes, Movies and referee reports, and may be
cited as

R. Tönjes,C.E. Fiore and T. Pereira,Coherence resonance in influencer
networks, Nat Commun 12, 72 (2021)

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Code Availability

Input files or sets of input parameters for Fortran as well as self-developed
Python codes are available from the corresponding author upon request. Python
code for full oscillator network simulations for an undirected scale free
network and for the C. elegans directed neuronal network can be found at
https://github.com/rahleph/CR-on-Influencer-Networks.git

## Author contributions

RT and TP wrote the text and developed the theory. RT made the slow fast
approximation and numerical solutions of the special functions. CEF made
simulations of Figure 1. TP and RT made the Figures. 5

## Competing interests.

Authors declare no competing interests.

## References

  * [1] Newman, M. E. J.  _Networks: An Introducton_ (Oxford University Press, Oxford, 2010).
  * [2] Xu, X.-K., Zhang, J. & Small, M.  Rich-club connectivity dominates assortativity and transitivity of complex networks.  _Phys. Rev. E_ 82, 046117 (2010).
  * [3] Watts, D. J. & Strogatz, S. H.  Collective dynamics of ‘small-world’networks.  _Nature_ 393, 440–442 (1998).
  * [4] van den Heuvel, M. P. & Sporns, O.  Network hubs in the human brain.  _Trends Cogn. Sci._ 17, 683–696 (2013).
  * [5] Song, S., Sjöström, P. J., Reigl, M., Nelson, S. & Chklovskii, D. B.  Highly nonrandom features of synaptic connectivity in local cortical circuits.  _PLoS Biol._ 3, e68 (2005).
  * [6] Fulcher, B. D. & Fornito, A.  A transcriptional signature of hub connectivity in the mouse connectome.  _Proc. Natl. Acad. Sci. USA_ 113, 1435–1440 (2016).
  * [7] Bonifazi, P. _et al._ Gabaergic hub neurons orchestrate synchrony in developing hippocampal networks.  _Science_ 326, 1419–1424 (2009).
  * [8] Schneidman, E., Berry, M. J., Segev, R. & Bialek, W.  Weak pairwise correlations imply strongly correlated network states in a neural population.  _Nature_ 440, 1007–1012 (2006).
  * [9] Buzsáki, G. & Mizuseki, K.  The log-dynamic brain: how skewed distributions affect network operations.  _Nat. Rev. Neurosci._ 15, 264–278 (2014).
  * [10] Gal, E. _et al._ Rich cell-type-specific network topology in neocortical microcircuitry.  _Nat. Neurosci._ 20, 1004 (2017).
  * [11] Adamic, L. A. & Glance, N.  The political blogosphere and the 2004 us election: divided they blog.  In _Proceedings of the 3rd international workshop on Link discovery_ , 36–43 (2005).
  * [12] Zou, Y., Pereira, T., Small, M., Liu, Z. & Kurths, J.  Basin of attraction determines hysteresis in explosive synchronization.  _Phys. Rev. Lett._ 112, 114102 (2014).
  * [13] Rodrigues, F. A., Peron, T. K. D., Ji, P. & Kurths, J.  The kuramoto model in complex networks.  _Phys. Rep._ 610, 1–98 (2016).
  * [14] Pereira, T., van Strien, S. & Tanzi, M.  Heterogeneously coupled maps: hub dynamics and emergence across connectivity layers.  _J. Eur. Math. Soc._ 22, 2183–2252 (2020).
  * [15] Nishikawa, T., Motter, A. E., Lai, Y.-C. & Hoppensteadt, F. C.  Heterogeneity in oscillator networks: Are smaller worlds easier to synchronize?  _Phys. Rev. Lett._ 91, 014101 (2003).
  * [16] Fries, P., Reynolds, J. H., Rorie, A. E. & Desimone, R.  Modulation of oscillatory neuronal synchronization by selective visual attention.  _Science_ 291, 1560–1563 (2001).
  * [17] Bruno, R. M. & Sakmann, B.  Cortex is driven by weak but synchronously active thalamocortical synapses.  _Science_ 312, 1622–1627 (2006).
  * [18] Morgan, R. J. & Soltesz, I.  Nonrandom connectivity of the epileptic dentate gyrus predicts a major role for neuronal hubs in seizures.  _Proc. Natl. Acad. Sci. USA_ 105, 6179–6184 (2008).
  * [19] Zamora-López, G., Zhou, C. & Kurths, J.  Cortical hubs form a module for multisensory integration on top of the hierarchy of cortical networks.  _Front. Neuroinf._ 4, 1 (2010).
  * [20] Lindner, B., Garcıa-Ojalvo, J., Neiman, A. & Schimansky-Geier, L.  Effects of noise in excitable systems.  _Phys. Rep._ 392, 321–424 (2004).
  * [21] Gammaitoni, L., Hänggi, P., Jung, P. & Marchesoni, F.  Stochastic resonance.  _Rev. Mod. Phys._ 70, 223–287 (1998).
  * [22] Pikovsky, A. S. & Kurths, J.  Coherence resonance in a noise-driven excitable system.  _Phys. Rev. Lett._ 78, 775 (1997).
  * [23] Wang, Y., Chik, D. T. & Wang, Z.  Coherence resonance and noise-induced synchronization in globally coupled hodgkin-huxley neurons.  _Phys. Rev. E_ 61, 740 (2000).
  * [24] Semenova, N., Zakharova, A., Anishchenko, V. & Schöll, E.  Coherence-resonance chimeras in a network of excitable elements.  _Phys. Rev. Lett._ 117, 014102 (2016).
  * [25] Masoliver, M., Malik, N., Schöll, E. & Zakharova, A.  Coherence resonance in a network of fitzhugh-nagumo systems: Interplay of noise, time-delay, and topology.  _Chaos_ 27, 101102 (2017).
  * [26] Ushakov, O. _et al._ Coherence resonance near a hopf bifurcation.  _Phys. Rev. Lett._ 95, 123903 (2005).
  * [27] Amro, R. M., Lindner, B. & Neiman, A. B.  Phase diffusion in unequally noisy coupled oscillators.  _Phys. Rev. Lett._ 115, 034101 (2015).
  * [28] Nakao, H., Arai, K. & Kawamura, Y.  Noise-induced synchronization and clustering in ensembles of uncoupled limit-cycle oscillators.  _Phys. Rev. Lett._ 98, 184101 (2007).
  * [29] Pimenova, A. V., Goldobin, D. S., Rosenblum, M. & Pikovsky, A.  Interplay of coupling and common noise at the transition to synchrony in oscillator populations.  _Sci. Rep._ 6, 38518 (2016).
  * [30] Dolmatova, A. V., Goldobin, D. S. & Pikovsky, A.  Synchronization of coupled active rotators by common noise.  _Physical Review E_ 96, 062204 (2017).
  * [31] Acebrón, J. A., Bonilla, L. L., Vicente, C. J. P., Ritort, F. & Spigler, R.  The kuramoto model: A simple paradigm for synchronization phenomena.  _Rev. Mod. Phys._ 77, 137 (2005).
  * [32] Stankovski, T., Pereira, T., McClintock, P. V. & Stefanovska, A.  Coupling functions: universal insights into dynamical interaction mechanisms.  _Rev. Mod. Phys._ 89, 045001 (2017).
  * [33] Montbrió, E. & Pazó, D.  Shear diversity prevents collective synchronization.  _Phys. Rev. Lett._ 106, 254101 (2011).
  * [34] Ko, T.-W. & Ermentrout, G. B.  Partially locked states in coupled oscillators due to inhomogeneous coupling.  _Phys. Rev. E_ 78, 016203 (2008).
  * [35] Sakaguchi, H. & Kuramoto, Y.  A soluble active rotater model showing phase transitions via mutual entertainment.  _Progress of Theoretical Physics_ 76, 576–581 (1986).
  * [36] Omel’chenko, E. & Wolfrum, M.  Nonuniversal transitions to synchrony in the sakaguchi-kuramoto model.  _Physical Review Letters_ 109, 164101 (2012).
  * [37] Giusti, C., Ghrist, R. & Bassett, D. S.  Twoś company, three (or more) is a simplex.  _J. Comp. Neurosci._ 41, 1–14 (2016).
  * [38] Vlasov, V., Zou, Y. & Pereira, T.  Explosive synchronization is discontinuous.  _Phys. Rev. E_ 92, 012904 (2015).
  * [39] Ott, E. & Antonsen, T. M.  Low dimensional behavior of large systems of globally coupled oscillators.  _Chaos_ 18, 037113 (2008).
  * [40] Tyulkina, I. V., Goldobin, D. S., Klimenko, L. S. & Pikovsky, A.  Dynamics of noisy oscillator populations beyond the ott-antonsen ansatz.  _Phys. Rev. Lett._ 120, 264101 (2018).
  * [41] Tönjes, R. & Pikovsky, A.  Low-dimensional description for ensembles of identical phase oscillators subject to cauchy noise.  _arXiv preprint arXiv:2009.04725_ (2020).
  * [42] Risken, H.  _Fokker-Planck Equation_ (Springer-Verlag, Berlin, 1989).
  * [43] Stratonovich, R. L.  _Topics in the theory of random noise_ , vol. 1 (Gordon and Breach, New York, 1963).

## Supplementary Note 1: Example of transformation into canonical form

We start with weakly non-identical ($\omega_{0}=2.0$,
$\Delta\omega/\omega_{0}<10\%$), weakly coupled ($\lambda_{0}=0.01$,
$\beta=10$) phase oscillators with moderately large shear ($c_{0}=0.5$) strong
noise ($D=1.4$) in the hubs and weak noise (or weak frequency heterogeneity)
in the followers ($D_{0}=2\times 10^{-4}$).

$\displaystyle\dot{\vartheta}_{n}$ $\displaystyle=$ $\displaystyle
2.0+\frac{0.01}{\mu_{n}}\sum_{m}W_{nm}\left[\sin(\vartheta_{m}-\vartheta_{n}-\alpha)+0.5\right]+\sqrt{4\times
10^{-4}}\xi_{n}$ (30) $\displaystyle\dot{\vartheta}_{\rm hub}$
$\displaystyle=$ $\displaystyle 2.135+\frac{0.1}{\mu_{\rm hub}}\sum_{m}W_{{\rm
hub}~{}m}\left[\sin(\vartheta_{m}-\vartheta_{\rm
hub}-\alpha)+0.5\right]+\sqrt{2\times 1.4}\xi_{\rm hub}$

After going into a co-rotating reference frame, where the follower mean
natural frequency is zero
$\vartheta\mapsto\vartheta-(\omega_{0}+\lambda_{0}c_{0})t$ and a change of
time scale $t\mapsto\lambda_{0}t$, we obtain

$\displaystyle\dot{\vartheta}_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{\mu_{n}}\sum_{m}W_{nm}\sin(\vartheta_{m}-\vartheta_{n}-\alpha)+\sqrt{2\times
0.02}\tilde{\zeta}_{n}$ (32) $\displaystyle\dot{\vartheta}_{\rm hub}$
$\displaystyle=$ $\displaystyle 18\left[1+\frac{10}{18}\frac{1}{\mu_{\rm
hub}}\sum_{m}W_{{\rm hub}~{}m}\sin(\vartheta_{m}-\vartheta_{\rm
hub}-\alpha)\right]+\sqrt{2\times\frac{140}{18}18}\tilde{\zeta}_{\rm hub}$

The effective dimensionless parameters for this system are
$\Delta\Omega/\lambda_{0}=18$, $\Lambda=10/18$ and $q=140/18\approx 7$. The
large noise strength $D/\lambda_{0}=140$ is only large in units of
$\lambda_{0}=0.01$. Here, $D=1.4$ is already much larger than the optimal
noise strength $D_{\rm opt}=0.18$.

## Supplementary Note 2: Hyper-graph

Figure 3: Three-influencer network and hyper-network. (a) network with three
influencers. (b) a hyper-network of seven partitions connected by three hyper-
edges.

We discuss the hyper-graph structure emerging at the macroscopic interaction
of mean-fields. We will provide illustration of an influencer network with
seven partitions of equal size connected through three influencers. After the
Ott-Antonsen reduction of the partition mean-field dynamics and after
averaging over the fast influencers, the resulting mean-field equations have
the following form

$\displaystyle\dot{Z}_{1}$ $\displaystyle=$ $\displaystyle
F\left(Z_{1},G_{a}\right),\quad\dot{Z}_{2}=F\left(Z_{2},G_{b}\right),\quad\dot{Z}_{3}=F\left(Z_{3},G_{c}\right),$
(34) $\displaystyle\dot{Z}_{4}$ $\displaystyle=$ $\displaystyle
F\left(Z_{4},\frac{1}{2}G_{a}+\frac{1}{2}G_{b}\right),\quad\dot{Z}_{5}=F\left(Z_{5},\frac{1}{2}G_{b}+\frac{1}{2}G_{c}\right),$
(35) $\displaystyle\dot{Z}_{6}$ $\displaystyle=$ $\displaystyle
F\left(Z_{6},\frac{1}{2}G_{c}+\frac{1}{2}G_{a}\right),\quad\dot{Z}_{7}=F\left(Z_{7},\frac{1}{3}G_{a}+\frac{1}{3}G_{b}+\frac{1}{3}G_{c}\right)$
(36)

with coupling functions

$\displaystyle G_{a}$ $\displaystyle=$ $\displaystyle
G\left(\frac{1}{4}Z_{1}+\frac{1}{4}Z_{4}+\frac{1}{4}Z_{6}+\frac{1}{4}Z_{7};\Lambda_{a},q_{a}\right)$
(37) $\displaystyle G_{b}$ $\displaystyle=$ $\displaystyle
G\left(\frac{1}{4}Z_{2}+\frac{1}{4}Z_{4}+\frac{1}{4}Z_{5}+\frac{1}{4}Z_{7};\Lambda_{b},q_{b}\right)$
(38) $\displaystyle G_{c}$ $\displaystyle=$ $\displaystyle
G\left(\frac{1}{4}Z_{3}+\frac{1}{4}Z_{5}+\frac{1}{4}Z_{6}+\frac{1}{4}Z_{7};\Lambda_{c},q_{c}\right)$
(39)

each connecting four partition mean-fields, and with Riccati dynamics

$F(Z,G)=\left(\frac{e^{-i\alpha}}{2}G-\frac{\gamma_{0}}{\lambda_{0}}Z-\frac{e^{i\alpha}}{2}\bar{G}Z^{2}\right)$
(40)

## Supplementary Note 3: Impact of follower heterogeneity and noise

In the main body of the manuscript we presented predictions and simulations
with fixed values for frequency heterogeneity $\gamma_{0}/\lambda_{0}=0.02$ or
noise strength $D_{0}/\lambda_{0}=0.02$ in the followers. The influence of
these quantities on the coherence resonance can be seen in the mean-field
analysis with slow-fast approximation. The frequency heterogeneity can be
expressed as a function of all other system parameters, in particular the
effective noise strength $q$ and the global order parameter $R$. The contour
plot of $\gamma_{0}/\lambda_{0}$ as a function of $q$ and $R$ shows the
resonance curves for constant follower heterogeneity (see Supplementary Figure
4).

Figure 4:  Effect of follower frequency heterogeneity on the coherence
resonance. (a) shows resonance curves of the order parameter as a function of
effective noise $q$ in the influencers, with effective coupling
$\Lambda=0.56$. Shown are the level sets of the right-handside. of Eq. (26),
which correspond to different frequency heterogeneities
$\gamma_{0}/\lambda_{0}$ in the followers, as indicated by the values in the
labels. Using Gaussian white noise of equal strengths $D_{0}/\lambda_{0}$
instead of Lorentzian frequency distribution for follower frequencies of width
$\gamma_{0}/\lambda_{0}$ results in similar resonance curves. (b) shows the
maximum value of the order parameter as a function of the frequency
heterogeneity for different values of effective coupling strength $\Lambda$.

## Supplementary Note 4: Examples of coherence resonance

We show three examples of coherence resonance in Supplementary Figure 5. For
the three-influencer network motif, we choose each partition with $100$ nodes
and undirected connections. We show an influencer network of $3000$ nodes with
$100$ influencers. Every node (also the influencers) is connected
symmetrically to $k+1$ randomly selected influencers where $k$ is
geometrically distributed with mean $10/6$. Since there are no direct links
between followers, all connections have the weight $W_{mn}=1$. For the
directed hyperlink network of 1033 political weblogs after the 2004 US
election, we choose the top 8 in-degree nodes as influencers. All nodes with
zero out-degree have been removed. Phase coupling in the oscillator dynamics
is realized in the opposite direction of the hyperlinks Ref. [11] of the main
manuscript. The time scale separation is $\Delta\Omega/\lambda_{0}=18$,
$\beta=10$ ($\Lambda=10/18$). The noise strength in the followers is
$D_{0}/\lambda_{0}=0.02$ in the three-influencer network and
$D_{0}/\lambda_{0}=0.04$ in the election blog network.

Figure 5: Coherence resonance of the order parameter for different influencer
networks. We show the time series of the order parameter $R$ for three values
of noise intensity in the influencers and the network in the corresponding
row. Parameters are $\Delta\Omega/\lambda_{0}=18$ and $\Lambda=10/18$. In
(a-d) we show the influencer network with three influencers and seven
partitions as shown in SI Note 2, with $D_{0}/\lambda_{0}=0.02$. In (e-h) we
show an influencer network with $3000$ nodes including $100$ influencers. In
(i-l), we show the directed hyperlink network of $1033$ political weblogs
after the 2004 US election (Ref. [11] of main manuscript) and the top 8-out
degree nodes as influencers. Arrows indicate coupling direction which is in
the opposite direction of the hyperlinks. The noise strength in the followers
is $D_{0}/\lambda_{0}=0.04$. The noise intensity in the influencers for the
simulations in (a,e,i) is weak $q=0.1$, in (b,f,j) optimal $q=1.0$, and in
(c,g,k) strong $q=10$.

## Supplementary Note 5: Effective diffusion

Here, we give an estimate of the effective noise strength for random forces
$c(t)=\lambda_{0}\cos(\psi-\alpha)$ and $s(t)=\lambda_{0}\sin(\psi-\alpha)$
when $\psi$ is a drift-diffusion process

$\psi(t+\tau)=\psi(t)+\Omega\tau+W_{D}(\tau)\quad(\textrm{mod}~{}2\pi)$ (41)

on the circle with constant velocity $\Omega$, Brownian diffusion
$W_{D}(\tau)$, and diffusion constant $D$. Such forces act on the follower
phases coupled to influencers with phase $\psi$ when the global order
parameter is zero, that is, in an incoherent state.

The conditional probability density $p(\psi,t+\tau|\psi_{0},t)$ is a wrapped
normal distribution $p_{WN}(\psi;\mu,\sigma^{2})$ with mean
$\mu=\psi_{0}+\Omega\tau$ and variance $\sigma^{2}=2D\tau$ , which in the
limit $\tau\to\infty$ becomes a stationary, uniform distribution on the
circle. The complex autocorrelation function

$C(\tau)=\lambda_{0}^{2}\left\langle e^{-i\psi(t)}\cdot
e^{i\psi(t+\tau)}\right\rangle_{t}$ (42)

is the sum of the autocorrelation functions of the two forces
$c(t)=\lambda_{0}\cos(\psi-\alpha)$ and $s(t)=\lambda_{0}\sin(\psi-\alpha)$ in
the real part and the cross-correlation function in the imaginary part.
Replacing the time average by the average with respect to the conditional
probability density and the average over the initial conditions $\psi_{0}$
with respect to the stationary probability density, the complex
autocorrelation function has the form of the first circular moment of the
wrapped normal distribution

$C(\tau)=\lambda_{0}^{2}\left\langle
e^{i\left(\Omega\tau+W_{D}(\tau)\right)}\right\rangle_{p_{{}_{WN}}}=\lambda_{0}^{2}e^{(i\Omega-D)\tau}.$
(43)

Because of the uniform stationary distribution, the autocorrelation functions
of $c(t)$ and $s(t)$ are identical and both forces are uncorrelated; that is
$\left\langle c(t)s(t)\right\rangle_{t}=0$. The time integrals of the random
forces are random variables that, according to the central limit theorem, have
a variance that grows asymptotically linearly in time. The effective diffusion
constant, which is half of the asymptotic speed of this growth, is the
integral of the autocorrelation functions

$D_{\rm
eff}=\frac{\lambda_{0}^{2}}{2}\int_{0}^{\infty}\textrm{Re}\left[e^{(i\Omega-D)\tau}\right]d\tau=\frac{\lambda_{0}^{2}}{2}\frac{D}{D^{2}+\Omega^{2}}$
(44)

Measuring time in units of $1/\lambda_{0}$, we obtain the expression in the
Methods section of the manuscript.

Figure 6: Time series of the order parameters $R_{\sigma}$ ($R_{1}$ red,
$R_{2}$ green, $R_{3}$ blue) for the network in Figure 1 of the main
manuscript with two influencers for (a-c) moderate time scale separation
$\Delta\Omega/\lambda_{0}=18$, noise $q=20/18$ and reduced coupling strength
$\Lambda=10/18$ and (d-f) large time scale separation
$\Delta\Omega/\lambda_{0}=198$ with noise $q=200/198$ and $\Lambda=100/198$.
Full network simulations with noise $D_{0}/\lambda_{0}=0.02$ in the followers
are shown in (a) and (d). Full network simulations with frequency
heterogeneity $\gamma_{0}/\lambda_{0}=0.02$ in the followers are shown in (b)
and (e). In (c) and (f) we show simulations of the Ott-Antonsen mean-field
equations under stochastic forcing by the influencers and of the averaged
dynamics (smooth curves) obtained under the assumption of infinite time scale
separation. The initial conditions are $Z_{1}=0.9\exp(i\pi/2)$,
$Z_{3}=0.9\exp(-i\pi/4)$, $Z_{3}=0$, $\psi_{a}=\pi$ and $\psi_{b}=-\pi/2$. The
drop in the order parameter $R_{3}$ is predicted by the averaged dynamics and
can be observed in the full network simulations with noise $D_{0}$ or
frequency heterogeneity $\gamma_{0}$.

## Supplementary Note 6: Comparison of reduced models

The Ott-Antonsen ansatz for ensembles of phase oscillators with Lorenzian
frequency distribution and under common forcing in the zeroth and in the first
harmonics results is an exact expression for the dynamics of the mean-field.
The forcing can be stochastic if interpreted as Stratonovich stochastic
differential equation. In the thermodynamic limit $N\to\infty$ but for finite
dynamical frequency gap $\Delta\Omega/\lambda_{0}$, the fluctuations in the
influencer phases lead to fluctuations in the partition mean-fields that are
comparable to the partition mean-field fluctuations in the finite size network
dynamics. For a large dynamical frequency gap, the effective noise strength of
the stochastic forcing vanishes. In the slow-fast approximation, the followers
are only subject to the average forces from the influencers, which in turn
depend adiabatically on the partition mean-fields. While fluctuations are
completely absent in the averaged dynamics, the average forces depend
continuously on the influencer effective noise strength $q=D/\Delta\Omega$.
Indeed, while for $\Delta\Omega/\lambda_{0}=18$ the averaged dynamics
describes the equilibrium order parameters near the optimal noise strength
reasonably well, for much larger time scale separation
$\Delta\Omega/\lambda_{0}=198$ even the transient to equilibrium is reproduced
robustly Supplementary Figure 6.


# On relation between renormalized frequency and heat capacity for particles
in an anharmonic potential

Y. T. Liu  School for Theoretical Physics, School of Physics and Electronics,
Hunan University, Changsha 410082, China  Y. H. Zhao  School for Theoretical
Physics, School of Physics and Electronics, Hunan University, Changsha 410082,
China  Y. Zhong<EMAIL_ADDRESS>School for Theoretical Physics, School of
Physics and Electronics, Hunan University, Changsha 410082, China  J. M. Shen
<EMAIL_ADDRESS>School for Theoretical Physics, School of Physics and
Electronics, Hunan University, Changsha 410082, China  J. H. Zhang  College of
Material Sciences, Hunan University, Changsha 410082, China  Q. H. Liu
<EMAIL_ADDRESS>School for Theoretical Physics, School of Physics and
Electronics, Hunan University, Changsha 410082, China

###### Abstract

For free particles in a simple harmonic potential plus a weak anharmonicity,
characterized by a set of anharmonic parameters, Newtonian mechanics asserts
that there is a renormalization of the natural frequency of the periodic
motion; and statistical mechanics claims that the anharmonicity causes a
correction to the heat capacity of an ideal gas in the anharmonic potential.
The orbital motion and thermal motion depend on the same anharmonic
parameters, but in different combinations. These two manners of combinations
are fundamentally different, demonstrating that statistical law can not emerge
from the many-body limit __ of deterministic law for one-body.

renormalization, anharmonicity, heat capacity, Poincare–Lindstedt method,
classical orbits, statistical law.

##  I Introduction

Anharmonicity plays a crucial role in modern physics, for instance, classical
anharmonic $\phi^{4}$ model in quantum field theory, qft  and various
anharmonic effects in condensed matter physics, solid1 ; solid2 ; solid3 ;
solid4 ; solid5 ; solid6 ; solid7 ; solid8 ; solid9  and secular evolution in
planetary orbits. liu19  In the perturbation theory, the anharmonicity is
usually associated with the renormalization of the natural frequency to remove
the superficial divergence, i.e., to eliminate the secular term in the naive
solution; liu19 ; RONALD ; PLRE ; liu  and the relationship between the
renormalization in quantum field theory and removal of the secular term in
perturbation expansion has been a subject under intensive investigations. QFT1
; QFT2 ; QFT3 ; QFT4 ; QFT5 ; QFT6 ; QFT7 ; QFT8  However, there are still
problems not yet fully understood. On one hand, in Newtonian mechanics the
anharmonicity usually leads to a renormalization of the natural frequency; and
on the other, in statistical mechanics the simplest situation is that the
anharmonicity causes a correction to the heat capacity for an ideal gas in the
anharmonic potential. We have therefore two different mechanics both start
from the same Hamiltonian to deal with the same many-body system: one is the
Newtonian mechanics from which every particle has its own orbit; and another
is the statistical mechanics from which every particle situates at a
microstate with a definite probability. An immediate question then arises: Can
the statistical law emerge from the _many-body limit_ of deterministic law for
one-body? This question may be of _fundamental importance_ , and in present
paper, an exactly solvable one-dimensional system is used to understand this
question in some depth.

Assume that there is a particle of mass $m$ moving in a potential field
$U\left(x\right)$ $\left(x\in\left(-\infty,\infty\right)\right)$ given by,

$U\left(x\right)=U\left(0\right)+\frac{1}{2}m\omega_{0}^{2}\ell^{2}\left(\left(\frac{x}{\ell}\right)^{2}-a\left(\frac{x}{\ell}\right)^{3}+b\left(\frac{x}{\ell}\right)^{4}+c\left(\frac{x}{\ell}\right)^{5}+d\left(\frac{x}{\ell}\right)^{6}\right),$
(1)

where $\omega_{0}$ is the natural frequency of the unperturbed potential
$m\omega_{0}^{2}x^{2}/2$, and $a$, $b$, $c$, and $d$ are four small
dimensionless parameters accounting for various orders of the anharmonicity,
and $\ell\left(\neq 0\right)$ is a characteristic length accounting for, e.g.
the anharmonicity and we usually set $U\left(0\right)=0$. We will call these
four constants $a$, $b$, $c$, and $d$ as anharmonic parameters. Once the
anharmonicity happens at infinity $\ell\rightarrow\infty$, $U\left(x\right)$
reproduces the usual harmonic one. The characteristic length $\ell$ can be
conveniently chosen to be the amplitude of the initial position, and can in
fact be freely specified because our conclusion is independent of its specific
value. When $c=d=0$, the potential (1) becomes,

$U\left(x\right)=\frac{1}{2}m\omega_{0}^{2}\ell^{2}\left(\left(\frac{x}{\ell}\right)^{2}-a\left(\frac{x}{\ell}\right)^{3}+b\left(\frac{x}{\ell}\right)^{4}\right).$
(2)

Landau used such a form of potential (2) describing the anharmonicity of the
vibrations and their interaction with the rotation within a diatomic molecule
(see Eq. (49.11) in Ref. Landau ), and we follow Landau’s convention Landau
to take a negative sign before first order anharmonicity
$a\left(x/\ell\right)^{3}$ in (1) and (2) though the interval of $x$ in Landau
model Landau  is half space $x\in\left(0,\infty\right)$ but what we are
interested in is the full one.

Every particle in the potential $U\left(x\right)$ (1) moves along an exclusive
trajectory, no matter what energy it has. However, each trajectory has its own
frequency provided that it takes an exclusive value of energy. As we show
shortly (c.f. Eqs. (11a) and (57)), we can expand the renormalized frequency
up to order $\left(x/\ell\right)^{4}$ in the following form,

$\omega\approx\omega_{0}\left(1+\chi^{\left(1\right)}\mu+\chi^{\left(2\right)}\mu^{2}+\chi^{\left(3\right)}\mu^{3}+\chi^{\left(4\right)}\mu^{4}\right),$
(3)

where parameter $\mu$ is dimensionless parameter, defined by,

$\mu\equiv\frac{A}{\ell}\succ 0,$ (4)

where $A\equiv x\left(t=0\right)\succ 0$ is initial position of the particle,
which can also be used to characterize the value of the energy the particle
takes, and once letting $\ell=A$, we have $\mu=1$. To note that the parameter
$\ell$ must be the same in both Newtonian mechanics and statistical physics,
otherwise the comparison of their results is meaningless. The common features
all trajectories in the potential (1) share are from (3)
$\chi^{\left(i\right)}=$ $\chi^{\left(i\right)}\left(a,b,c,d\right)$
($i=1,2,3,4$), and we call $\chi^{\left(i\right)}$ the $i$-th order _orbital
anharmonicity_ (OA). It is worth stressing that OAs are independent of initial
conditions $x$ and $dx/dt$, whose possible uncertainty or stochasticity does
not effect OAs. In analogue, we will introduce $i$-th order _thermal
anharmonicity_ $\zeta^{\left(i\right)}=$
$\zeta^{\left(i\right)}\left(a,b,c,d\right)$ in the similar expansion of the
heat capacity,

$C\approx
C_{0}\left(1+\alpha_{1}\zeta^{\left(1\right)}+\alpha_{2}^{\left(2\right)}\zeta^{\left(2\right)}+\alpha_{3}\zeta^{\left(3\right)}+\alpha_{4}\zeta^{\left(4\right)}\right),$
(5)

where $\alpha_{i}$ are some expansion coefficients, and $C_{0}=Nk_{B}$ with
$N$ the number of the particle and $k_{B}$ the Boltzmann constant. The key
finding of the present study is,

$\zeta^{\left(i\right)}\text{ (}i=1,2,3,4\text{) is linearly independent of
OAs
(}\chi^{\left(1\right)},\chi^{\left(2\right)},\chi^{\left(3\right)},\chi^{14}\text{).}$
(6)

We are confident that it is a completely novel and physically significant
result for OA and _thermal anharmonicity_ offer proper and faithful
characterization of the anharmonicity from the point of Newtonian dynamics and
thermodynamics, respectively. The importance of the inequivalence between OA
and _thermal anharmonicity_ can be understood from the opposite but untrue
limit: If $\chi^{\left(i\right)}=\zeta^{\left(i\right)}$, we could safely say
that statistical law for the many-body system and the many-body limit of
Newtonian mechanics for every particle in it are at least heavily overlapped
and are even of same origin in nature. Otherwise, Eq. (6) strongly suggests
that the statistical law can not emerge from the many-body limit __ of
deterministic law for one-body. In very rough terms, the molecular dynamics
can not exactly and completely reproduce all thermodynamic results.

This paper is organized as follows. Section II and III give the detailed steps
of calculations of $\omega$ and $C$ with $c=d=0$, and Section IV presents only
the final results of both $\omega$ and $C$ with nonvanishing $c$ and $d$, and
all calculational steps are omitted. Explicitly, in section II, we utilize the
Poincare–Lindstedt method to solve the equation of motion of position $x$ in
terms of time $t$, from which we see in detail how the natural frequency is
renormalized. In section III, the anharmonicity induced correction of the heat
capacity is calculated and an order-by-order comparison between the
renormalized frequency and the heat capacity is made. In section IV, the
potential containing higher order anharmonicities with $c\neq 0$ and $d\neq 0$
in (1) is studied, and we see clearly that heat capacity $C$ not only involves
these OAs $\chi^{\left(1\right)},\chi^{\left(2\right)},\chi^{\left(3\right)},$
and $\chi^{\left(4\right)}$ but also the anharmonic parameters $a,b,c$, and
$d$. In final section V, a brief conclusion is given.

##  II Renormalized frequencies for second order anharmonic oscillator

The equation of motion for one particle in the potential (2) is,

$m\frac{d^{2}x}{dt^{2}}=-\frac{dU\left(x\right)}{dx}=-m\omega_{0}^{2}\left(x-\frac{3}{2}a\ell\left(\frac{x}{\ell}\right)^{2}+2b\ell\left(\frac{x}{\ell}\right)^{3}\right),$
(7)

where the initial conditions at instant $t=0$ are,

$x\left(0\right)=A,\frac{dx(0)}{dt}=0.$ (8)

Making a variable transform,

$x(t)\rightarrow\ell\varphi(t),$ (9)

we have from (4), (7) and (8),

$\frac{d^{2}\varphi}{dt^{2}}=-\omega_{0}^{2}\left(\varphi-\frac{3}{2}a\varphi^{2}+2b\varphi^{3}\right),\varphi\left(0\right)=\mu,\frac{d\varphi(0)}{dt}=0.$
(10)

To this equation, no exact solution is possible due to the nonlinearity in
$\varphi$, and even worse, the regular perturbation approaches fail for they
lead to the secular term in the solutions of $\varphi=\varphi(t)$. Instead,
the Poincare–Lindstedt method gives uniformly valid asymptotic expansions for
the periodic solutions of weakly nonlinear oscillations. RONALD ; PLRE ; liu
By the method, we mean that following three transformations must be done
simultaneously,

$\displaystyle\omega_{0}$ $\displaystyle\rightarrow$
$\displaystyle\omega=\omega_{0}+\omega_{1}+\omega_{2}+...,$ (11a)
$\displaystyle t$ $\displaystyle\rightarrow$
$\displaystyle\tau=\frac{\omega_{0}}{\omega}t,\text{ }$ (11b) $\displaystyle
x\left(t\right)$ $\displaystyle\rightarrow$
$\displaystyle\ell\xi\left(\tau\right)=x\left(t\left(\tau\right)\right),$
(11c) where $\omega_{1}\sim$ $O(a)$ and $\omega_{2}\sim$ $O(a^{2})\sim O(b)$
are the first and second order renormalization of the frequency, and so forth.
In the same time, we have,
$\xi(\tau)\approx\xi_{0}(\tau)+\xi_{1}(\tau)+\xi_{2}(\tau)+...$ (12)

in which $\xi_{0}(\tau)$ is the equation of motion for the unperturbed
oscillator satisfying the initial conditions,

$\xi_{0}(0)=\mu,\frac{d\xi_{0}(0)}{d\tau}=0,$ (13)

and $\xi_{1}\sim$ $O(a)$ and $\xi_{2}\sim$ $O(a^{2})$ are the first and second
order correction of the position $\xi(\tau)$, with the initial conditions,
respectively,

$\xi_{i}(0)=0,\frac{d\xi_{i}(0)}{d\tau}=0,(i=1,2).$ (14)

The correct equation of motion takes the following form, accurate up to $O(b)$
or $O(a^{2})$,

$\frac{d^{2}\xi}{d\tau^{2}}\approx-\left(\omega_{0}+\omega_{1}+\omega_{2}\right)^{2}\left(\xi-\frac{3}{2}a\xi^{2}+2b\xi^{3}\right),$
(15)

The zeroth, first, and second order equations of motion of Eq. (15) are,
respectively,

$\displaystyle\frac{d^{2}\xi_{0}}{d\tau^{2}}+\omega_{0}^{2}\xi_{0}=0,$ (16)
$\displaystyle\frac{d^{2}\xi_{1}}{d\tau^{2}}+\omega_{0}^{2}\xi_{1}+\left(-\frac{3}{2}a\omega_{0}^{2}\xi_{0}^{2}+2\omega_{1}\omega_{0}\xi_{0}\right)=0,$
(17)
$\displaystyle\frac{d^{2}\xi_{2}}{d\tau^{2}}+\omega_{0}^{2}\xi_{2}+\left(2\omega_{1}\omega_{0}-3a\omega_{0}^{2}\xi_{0}\right)\xi_{1}+2\omega_{0}^{2}b\xi_{0}^{3}-3a\omega_{1}\omega_{0}\xi_{0}^{2}+\left(2\omega_{2}\omega_{0}+\omega_{1}^{2}\right)\xi_{0}=0.$
(18)

The zeroth order equation of motion gives the usual harmonic oscillatory
solution,

$\xi_{0}(\tau)=\mu\cos\left(\omega_{0}\tau\right).$ (19)

The naive solution of the first order equation of motion is then,

$\xi_{1}(\tau)=-\omega_{1}\tau\mu\sin\left(\omega_{0}\tau\right)+\frac{1}{4}a\mu^{2}\left(3-2\cos\left(\omega_{0}\tau\right)-\cos\left(2\omega_{0}\tau\right)\right).$
(20)

The first term in the right-hand side gives the divergent oscillatory
amplitude $\omega_{1}\tau\mu$ as time $\tau\rightarrow\infty$ with
$\omega_{1}\neq 0$. To remove the divergence, we have to choose,

$\omega_{1}=0.$ (21)

The correct first order solution of equation of motion is thus,

$\xi_{1}(\tau)=\frac{1}{4}a\mu^{2}\left(3-2\cos\left(\omega_{0}\tau\right)-\cos\left(2\omega_{0}\tau\right)\right).$
(22)

The naive solution of the second order equation of motion is,

$\displaystyle\xi_{2}(\tau)$ $\displaystyle=$
$\displaystyle\frac{\mu\tau}{16}\left(3\mu^{2}\left(5a^{2}-4b\right)\omega_{0}-16\omega_{2}\right)\sin(\omega_{0}\tau)$
(23)
$\displaystyle+\frac{\mu^{3}}{16}\left(-12a^{2}+\left(\frac{29}{4}a^{2}-b\right)\cos\left(\omega_{0}\tau\right)+4a^{2}\cos\left(2\omega_{0}\tau\right)+\left(b+\frac{3}{4}a^{2}\right)\cos\left(3\omega_{0}\tau\right)\right).$

The first term in the right-hand side gives also the divergent oscillatory
amplitude as time $\tau\rightarrow\infty$, and this divergence can simply be
removed with $\omega_{2}$ being selected to satisfy,

$3\mu^{2}\left(5a^{2}-4b\right)\omega_{0}-16\omega_{2}=0.$ (24)

We have the second order correction of the frequency $\omega_{2}$,

$\omega_{2}=\frac{3}{16}\left(5a^{2}-4b\right)\mu^{2}\omega_{0}=\chi^{\left(2\right)}\mu^{2}\omega_{0}.$
(25)

where $\chi^{\left(2\right)}$ is the second order OF __ which is a combination
of second order anharmonic parameters $b$ and $a^{2}$,

$\chi^{\left(2\right)}\equiv\frac{3}{16}\left(5a^{2}-4b\right).$ (26)

The correct second order solution of equation of motion (18) is,

$\xi_{2}(\tau)=\frac{\mu^{3}}{16}\left(-12a^{2}+\left(\frac{29}{4}a^{2}-b\right)\cos\left(\omega_{0}\tau\right)+4a^{2}\cos\left(2\omega_{0}\tau\right)+\left(\frac{3}{4}a^{2}+b\right)\cos\left(3\omega_{0}\tau\right)\right).$
(27)

The important result is then that the natural frequency $\omega_{0}$ is
renormalized to be, up to accuracy of second order anharmonicity $O(a^{2})\sim
O(b)$,

$\omega=\left(1+\chi^{\left(2\right)}\mu^{2}\right)\omega_{0}.$ (28)

The oscillation is composed of a single prime frequency $\omega_{0}$ and its
higher order harmonics,

$\displaystyle\xi(\tau)$ $\displaystyle\approx$
$\displaystyle\mu\cos\left(\omega_{0}\tau\right)+\frac{1}{4}a\mu^{2}\left(3-2\cos\left(\omega_{0}\tau\right)-\cos\left(2\omega_{0}\tau\right)\right)$
(29)
$\displaystyle+\frac{\mu^{3}}{16}\left(-12a^{2}+\left(\frac{29}{4}a^{2}-b\right)\cos\left(\omega_{0}\tau\right)+4a^{2}\cos\left(2\omega_{0}\tau\right)+\left(b+\frac{3}{4}a^{2}\right)\cos\left(3\omega_{0}\tau\right)\right).$

It is easily to verify that the energy is conserved for we have,

$E\left(t\right)=\frac{1}{2}m\left(\frac{dx}{dt}\right)^{2}+U\left(x\right)=\frac{1}{2}m\omega_{0}^{2}A^{2}\left(1-a\mu+b\mu^{2}\right)=E\left(t=0\right).$
(30)

The anharmonicity induced correction of the energy is,

$\Delta E=\frac{1}{2}m\omega_{0}^{2}A^{2}\left(-a\mu+b\mu^{2}\right).$ (31)

Requiring that the energy shift is small, we have,

$\frac{\left|\Delta
E\right|}{\frac{1}{2}m\left(\omega_{0}A\right)^{2}}=\left|-a\mu+b\mu^{2}\right|\ll
1.$ (32)

The sufficient conditions for this equation are,

$\left|a\right|\ll 1,\left|b\right|\ll 1.$ (33)

This is what small constants $a$ and $b$ mean in Newtonian mechanics. Once
these conditions break, the perturbation method (15) does not apply.

##  III Second order anharmonicity induced correction of heat capacity

For our purpose, we need to compute the partition function in Boltzmann
statistical mechanics, with $H=p^{2}/2m+U(x)$,

$Z\equiv\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\exp(-\beta
H)\frac{dxdp}{h}=Z_{T}Z_{U},$ (34)

where $h$ is the Planck’s constant, and $Z_{T}$ is the momentum factor of the
partition function divided by $h$,

$Z_{T}\equiv\int_{-\infty}^{\infty}\exp(-\beta\frac{p^{2}}{2m})\frac{dp}{h}=\frac{\sqrt{2\pi}}{h}\sqrt{\frac{m}{\beta}},$
(35)

and $Z_{U}$ is the configurational factor of the partition function, with
transform $x\rightarrow\ell\xi$,

$\displaystyle Z_{U}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\exp(-\beta U(x))dx$ (36)
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\exp\left(-\beta\left(\frac{1}{2}m\omega_{0}^{2}\ell^{2}\left(\left(\frac{x}{\ell}\right)^{2}-a\left(\frac{x}{\ell}\right)^{3}+b\left(\frac{x}{\ell}\right)^{4}\right)\right)\right)dx$
$\displaystyle=$
$\displaystyle\ell\int_{-\infty}^{\infty}\exp\left(-\frac{\xi^{2}}{2\eta^{2}}\right)\exp\left(-\frac{-a\xi^{3}+b\xi^{4}}{2\eta^{2}}\right)d\xi$
$\displaystyle\approx$
$\displaystyle\ell\int_{-\infty}^{\infty}\exp\left(-\frac{\xi^{2}}{2\eta^{2}}\right)\left(1+\frac{a\xi^{3}-b\xi^{4}}{2\eta^{2}}+\frac{1}{8}\left(\frac{a\xi^{3}}{\eta^{2}}\right)^{2}\right)d\xi$
$\displaystyle=$
$\displaystyle\sqrt{2\pi}\ell\eta\left(1+\frac{3}{8}\left(5a^{2}-4b\right)\eta^{2}\right),$

where $\eta$ is the dimensionless parameter, defined by,

$\eta\equiv\sqrt{\frac{1}{\beta
m\omega_{0}^{2}\ell^{2}}}=\sqrt{\frac{k_{B}T}{m\omega_{0}^{2}\ell^{2}}}.$ (37)

The partition function is then,

$Z\equiv Z_{T}Z_{U}\approx\frac{2\pi}{\beta
h\omega_{0}}\left(1+\frac{2\chi^{\left(2\right)}}{\beta
m\omega_{0}^{2}\ell^{2}}\right)=\frac{2\pi}{\beta
h\omega_{0}}\left(1+2\chi^{\left(2\right)}\eta^{2}\right).$ (38)

Once $\chi^{\left(2\right)}$ is negligible, the partition function reduces to
be,

$Z\approx\frac{2\pi}{\beta h\omega_{0}},$ (39)

which leads to the energy equipartition result for the heat capacity of the
oscillatory degree of freedom,

$C_{0}=Nk_{B}.$ (40)

The anharmonicity gives rise to the correction to the internal energy, up to
accuracy of second order anharmonicity $O(a^{2})\sim O(b)$,

$\Delta
U^{(2)}=-N\frac{\partial}{\partial\beta}\ln\left(1+\frac{2\chi^{\left(2\right)}}{\beta
m\omega_{0}^{2}\ell^{2}}\right)\approx-\frac{2\chi^{\left(2\right)}N}{m\omega_{0}^{2}\ell^{2}}\frac{\partial}{\partial\beta}\frac{1}{\beta}=\frac{2\chi^{\left(2\right)}N}{m\omega_{0}^{2}\ell^{2}\beta^{2}}.$
(41)

The corresponding correction to heat capacity is proportional to the first
power of the temperature via $\eta^{2}$ (37),

$\Delta C^{(2)}=\frac{\partial\Delta U}{\partial T}\approx
4\eta^{2}\chi^{\left(2\right)}Nk_{B},$ (42)

which is also compatible with the Landau’s result. Landau

Evidently, in statistical mechanics, $\Delta C^{(2)}$ is the second order
quantity which linearly depends on $\chi^{\left(2\right)}$. Thus, we have the
heat capacity,

$C\approx\left(1+4\eta^{2}\chi^{\left(2\right)}\right)Nk_{B}.$ (43)

It takes the form (5) which suggests in general a  _linear_ dependence of $C$
on $\chi^{\left(i\right)}$. Whether such a form (5) persists for the heat
capacity with higher order anharmonicities is an interesting problem. In next
section, we show that this _linear_ dependence on the OAs breaks.

##  IV High order anharmonicities: renormalization of the natural frequency
and heat capacity correction

When $c\neq 0$ and $d\neq 0$, the equation of motion (15) becomes,

$\frac{d^{2}\xi}{d\tau^{2}}\approx-\left(\omega_{0}+\omega_{1}+\omega_{2}+\omega_{3}+\omega_{4}\right)^{2}\left(\xi-\frac{3}{2}a\xi^{2}+2b\xi^{3}+\frac{5}{2}c\xi^{4}+3d\xi^{5}\right).$
(44)

Utilization of the Poincare–Lindstedt method to solve this equation of motion
of position $\xi$ in terms of time $\tau$, Eq. (44) gives results for each
order in the following. The first three solutions
($\xi_{0}(\tau),\xi_{1}(\tau),\xi_{2}(\tau)$) are already given in (29), and
the third order solution $\xi_{3}(\tau)$ is,

$\xi_{3}(\tau)=\frac{\eta^{4}}{128}(\Lambda_{0}+\Lambda_{1}\cos(\omega_{0}\tau)+\Lambda_{2}\cos(2\omega_{0}\tau)+\Lambda_{3}\cos(3\omega_{0}\tau)+\Lambda_{4}\cos(4\omega_{0}\tau)),$
(45)

where,

$\displaystyle\Lambda_{0}$ $\displaystyle=$ $\displaystyle
3\left(75a^{3}-84ab-40c\right),$ (46a) $\displaystyle\Lambda_{1}$
$\displaystyle=$ $\displaystyle-119a^{3}+140ab+64c,$ (46b)
$\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{32}{3}\left(-9a^{3}+12ab+5c\right),$ (46c)
$\displaystyle\Lambda_{3}$ $\displaystyle=$
$\displaystyle-3a\left(3a^{2}+4b\right),$ (46d) $\displaystyle\Lambda_{4}$
$\displaystyle=$ $\displaystyle\frac{1}{3}\left(-3a^{3}-12ab+8c\right).$ (46e)
The fourth order solution $\xi_{4}(\tau)$ is,
$\xi_{4}(\tau)=\frac{3\eta^{5}}{64}(\Omega_{0}+\Omega_{1}\cos(\omega_{0}\tau)+\Omega_{2}\cos(2\omega_{0}\tau)+\Omega_{3}\cos(3\omega_{0}\tau)+\Omega_{4}\cos(4\omega_{0}\tau)+\Omega_{5}\cos(5\omega_{0}\tau)),$
(47)

where,

$\displaystyle\Omega_{0}$ $\displaystyle=$
$\displaystyle-a\left(75a^{3}-116ab-56c\right),$ (48a)
$\displaystyle\Omega_{1}$ $\displaystyle=$
$\displaystyle\frac{2357a^{4}}{64}+\frac{23b^{2}}{12}-\frac{1475a^{2}b}{24}-\frac{292ac}{9}-\frac{8d}{3},$
(48b) $\displaystyle\Omega_{2}$ $\displaystyle=$
$\displaystyle\frac{16}{9}a\left(18a^{3}-30ab-13c\right),$ (48c)
$\displaystyle\Omega_{3}$ $\displaystyle=$
$\displaystyle\frac{93a^{4}}{16}-2b^{2}-\frac{11a^{2}b}{4}+\frac{3ac}{4}+\frac{5d}{2},$
(48d) $\displaystyle\Omega_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{9}a\left(3a^{3}+12ab-8c\right),$ (48e)
$\displaystyle\Omega_{5}$ $\displaystyle=$
$\displaystyle\frac{5a^{4}}{192}+\frac{b^{2}}{12}+\frac{5a^{2}b}{24}-\frac{11ac}{36}+\frac{d}{6}.$
(48f) The third order and fourth order normalized frequencies are,
respectively,
$\text{$\omega$}_{3}=-a\chi^{\left(2\right)}\eta^{3}\omega_{0}=\chi^{\left(3\right)}\eta^{3}\omega_{0},$
(49)

where $\chi^{\left(3\right)}$ is the third order OA, defined by,__

$\chi^{\left(3\right)}\equiv-a\chi^{\left(2\right)},$ (50)

and,

$\text{$\omega$}_{4}=\frac{3\eta^{4}}{1024}\left(1155a^{4}-2200a^{2}b-1120ac+304b^{2}-320d\right)\omega_{0}=\chi^{\left(4\right)}\eta^{4}\omega_{0},$
(51)

where $\chi^{\left(4\right)}$ is the fourth order OA _,_ formed by a non-
trivial combination of all fourth order parameters
$(a^{4},a^{2}b,ac,b^{2},d)$, defined by,

$\chi^{\left(4\right)}\equiv\frac{3\eta^{4}}{1024}\left(1155a^{4}-2200a^{2}b-1120ac+304b^{2}-320d\right).$
(52)

We anticipate that $\chi^{\left(3\right)}$ and $\chi^{\left(4\right)}$ will
appear in the heat capacity.

To see how two quantities $\chi^{\left(3\right)}$ and $\chi^{\left(4\right)}$
may appear in the higher order corrections to heat capacity, let us compute
the partition function $Z\equiv Z_{T}Z_{U}$ with the full form of potential
(1). The result is, with calculations similar to (36)-(38),

$Z\equiv Z_{T}Z_{U}\approx\frac{2\pi}{\beta
h\omega_{0}}\left(1+2\chi^{\left(2\right)}\eta^{2}+\gamma\eta^{4}\right).$
(53)

where,

$\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{15}{128}\left(7\left(33a^{4}-72a^{2}b-32ac+16b^{2}\right)-64d\right)$
(54) $\displaystyle=$
$\displaystyle-8b\chi^{\left(2\right)}+8\chi^{\left(4\right)}.$

The anharmonicity gives rise to the correction to the internal energy in the
following,

$\Delta
U=-N\frac{\partial}{\partial\beta}\ln\left(1+\frac{2\chi^{\left(2\right)}}{\beta
m\omega_{0}^{2}\ell^{2}}+\frac{\gamma}{\left(\beta
m\omega_{0}^{2}\ell^{2}\right)^{2}}\right)\approx
Nk_{B}T\left(\frac{2\chi^{\left(2\right)}}{\beta
m\omega_{0}^{2}\ell^{2}}+2\frac{\gamma-2\left(\chi^{\left(2\right)}\right)^{2}}{\left(\beta
m\omega_{0}^{2}\ell^{2}\right)^{2}\beta}\right)$ (55)

The corresponding correction to the heat capacity is,

$\Delta C=\frac{\partial\Delta U}{\partial
T}=Nk_{B}\left(4\chi^{\left(2\right)}\eta^{2}-12\left(\left(\chi^{\left(2\right)}\right)^{2}+4b\chi^{\left(2\right)}-4\chi^{\left(4\right)}\right)\eta^{4}\right).$
(56)

Thus fourth order correction to the heat capacity does not depend on the OAs
__ alone because of presence of a term $b\chi^{\left(2\right)}$. Collecting
all results together, we have the renormalized frequency $\omega$ and heat
capacity $C$, respectively,

$\displaystyle\omega$ $\displaystyle\approx$
$\displaystyle\omega_{0}\left(1+\chi^{\left(2\right)}\mu^{2}-a\chi^{\left(2\right)}\mu^{3}+\chi^{\left(4\right)}\mu^{4}\right),\text{
and }$ (57) $\displaystyle C$ $\displaystyle\approx$ $\displaystyle
Nk_{B}\left(1+4\chi^{\left(2\right)}\eta^{2}-12\left(\left(\chi^{\left(2\right)}\right)^{2}+4b\chi^{\left(2\right)}-4\chi^{\left(4\right)}\right)\eta^{4}\right).$
(58)

In terms of the thermal anharmonicity, we have two nonzero elements, the 2nd
and fourth order anharmonicity $\chi^{\left(2\right)}$ and
$\left(\chi^{\left(2\right)}\right)^{2}+4b\chi^{\left(2\right)}-4\chi^{\left(4\right)}$.
It is evidently that both OA and the thermal anharmonicity have only two
common elements, the first and second order anharmonicities
$\chi^{\left(1\right)}\left(=0\right)$ and $\chi^{\left(2\right)}$, showing
that from particle orbits statistical law can be largely re-constructed.
However, in essence, the thermal properties are independent of the orbits for
statistical mechanics dictate its own manner of dependence on the anharmonic
parameters.

Before closing this section, we discuss an interesting situation which a
little bit deviates the theme of present study. In statistical mechanics, we
can theoretically assume the characteristic length to be a ”thermal one”
$\ell=\sqrt{k_{B}T/m\omega_{0}^{2}}$. Then, we have $C=Nk_{B}$ which is
irrelevant to the anharmonicity. The demonstration is straightforward, because
the configurational factor of the partition function $Z_{U}$ is, with
transform $x\rightarrow\ell\xi$,

$\displaystyle Z_{U}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\exp(-\beta U(x))dx$ (59)
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\exp\left(-\frac{\beta}{2}m\omega_{0}^{2}\ell^{2}\left(\left(\frac{x}{\ell}\right)^{2}+\sum_{j=3}a_{j}\left(\frac{x}{\ell}\right)^{j}\right)\right)dx$
$\displaystyle=$
$\displaystyle\ell\int_{-\infty}^{\infty}\exp\left(-\frac{1}{2}\sum_{j=2}a_{j}\xi^{j}\right)d\xi$
$\displaystyle=$ $\displaystyle\ell f\left(\left\\{a_{j}\right\\}\right),$

where $U(x)$ contains anharmonicity of arbitrarily high orders, and $a_{2}=1$
and $a_{j}$ $(j\geq 3)$ are anharmonic parameters, and
$f\left(\left\\{a_{j}\right\\}\right)\equiv\int_{-\infty}^{\infty}\exp\left(-\frac{1}{2}\sum_{j=2}a_{j}\xi^{j}\right)d\xi$
is independent of the temperature. Since the dependence of $Z_{U}$ on the
temperature is via $\ell=\sqrt{k_{B}T/m\omega_{0}^{2}}$ only, the partition
function becomes
$Z=\left(k_{B}T/\hbar\omega_{0}\right)f\left(\left\\{a_{j}\right\\}\right)$.
We have immediately $C=Nk_{B}$. To note that the anharmonic potential of form
(1) with $\ell=\sqrt{k_{B}T/m\omega_{0}^{2}}$ is in fact problematic because
by definition, Hamiltonian must be temperature-indepedent. However, from the
pure theoretical consideration, such a potential may be considered as an
effective one, which may enrich our understanding of the energy equipartition
theorem.

##  V Conclusions

A particle in a simple harmonic potential is fully understood. However, once
the potential is added by some weakly anharmonic terms, the problem becomes
highly non-trivial as we learn from Kolmogorov–Arnold–Moser theorem and
Fermi–Pasta–Ulam–Tsingou nonlinear lattice oscillations. Once a particle moves
in the anharmonic potential field, the natural frequency must be renormalized
and the OAs can then be introduced; and the dependence of the OAs on the
anharmonic parameters is dictated by the Newtonian mechanics. For a classical
ideal gas in the same anharmonic potential field, the internal energy and heat
capacity have their own ways of dependence on the anharmonic parameters,
determined by the statistical mechanics, and the corresponding so-called
thermal anharmonicities are introduced. The OAs and thermal anharmonicities
reflect the anharmonicity in potential in mechanics and thermodynamics,
respectively. The first two order anharmonicities in the mechanics and the
thermodynamics are the same, whereas the third and fourth order
anharmonicities are different. Though no higher order anharmonicities are
calculated, we can conclude that the statistical law can not emerge from the
_many-body limit_ of deterministic law for few-body.

It seems to us that the clear difference between the OA and the thermal
anharmonicity is useful in exploring the relation between a single particle
that obeys the Newtonian mechanics and many particles that follows the
statistical mechanics. The applications to other problems are under
exploration.

###### Acknowledgements.

QHL is grateful to the members, especially to Professor Hong Qian at
University of Washington and Hong Zhao at Xiamen University and Professor
Zhigang Zheng at Huaqiao University, of Online Club Nanothermodynamica
(Founded in June 2020) for extensive discussions of various problems in
statistical physics. We are indebted to Dr. Xinyuan Ai for participation of
the early stage of the present work. This work is financially supported by
National Natural Science Foundation of China under Grant No. 11675051 and No.
11905056.

## References

  * (1) Claude Itzykson, Jean-Bernard Zuber, Quantum Field Theory, (New York: Dover, 2006)
  * (2) R. Mankowsky, A. Subedi, M. Först, et al. Nonlinear lattice dynamics as a basis for enhanced superconductivity in $YBa_{2}Cu_{3}O_{6.5}$. Nature 516(2014)71–73.
  * (3) Yue Chen, Xinyuan Ai, and C. A. Marianetti, First-Principles Approach to Nonlinear Lattice Dynamics: Anomalous Spectra in $PbTe$, Phys. Rev. Lett. 113(2014)105501.
  * (4) C. W. Li, J. Ma, H. B. Cao, A. F. May, D. L. Abernathy, G. Ehlers, C. Hoffmann, X. Wang, T. Hong, A. Huq, O. Gourdon, and O. Delaire, Anharmonicity and atomic distribution of $SnTe$ and $PbTe$ thermoelectrics, Phys. Rev. B 90(2014)214303.
  * (5) Yao Tian, Shuang Jia, R. J. Cava, Ruidan Zhong, John Schneeloch, Genda Gu, and Kenneth S. Burch, Understanding the evolution of anomalous anharmonicity in $Bi_{\mathit{2}}Te_{\mathit{3-x}}Se_{\mathit{x}}$, Phys. Rev. B 95(2017)094104.
  * (6) R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Maur, Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation, Phys. Rev. B 96(2017)014111.
  * (7) A. S. T. Nguetcho, G. M. Nkeumaleu, and J M Bilbault, Behavior of gap solitons in anharmonic lattices, Phys. Rev. E, 96(2017)022207.
  * (8) Yi Xia, Revisiting lattice thermal transport in PbTe: The crucial role of quartic anharmonicity, Appl. Phys. Lett., 113(2018)073901.
  * (9) Qiya Liu, Ruihui Shao, Ning Li, Weizheng Liang, Xinsheng Yang, S. N. Luo, and Yong Zhao, Anharmonicity of $Bi_{2}Se_{3}$ revealed by fs transient optical spectroscopy, Appl. Phys. Lett. 115(2019)201902.
  * (10) M. E. Manley, O. Hellman, N. Shulumba, A. F. May, P. J. Stonaha, J. W. Lynn, V. O. Garlea, A. Alatas, R. P. Hermann, J. D. Budai, H. Wang, B. C. Sales, and A. J. Minnich, Intrinsic anharmonic localization in thermoelectric $PbSe$, Nat. Commun. 10(2019)1928.
  * (11) Q. H. Liu, Q. Li, T. G. Liu and X. Wang, Gravitational major-axis contraction of Mercury’s elliptical orbit, Mod. Phys. Lett. A, 34(2019)1950159.
  * (12) Ronald E. Mickens, An introduction to nonlinear oscillations, Cambridge university press, (London: Cambridge, 1981). Chapter 2\.
  * (13) L.-Y. Chen, N. Goldenfeld and Y. Oono, The renormalization group and singular perturbations: multiple scales, boundary layers and reductive perturbation theory, Phys. Rev. E 54(1996)376.
  * (14) Q. H. Liu, Z. Li, M. N. Zhang, Q. Li, B .J. Chen, Exponential functions of perturbative series and elimination of secular divergences in time-dependent perturbation theory in quantum mechanics, Results in Physics, 7(2017)890.
  * (15) B. Craps, O. Evnin, and J. Vanhoof, Renormalization group, secular term resummation and AdS (in)stability. J. High Energ. Phys. 2014(2014)48.
  * (16) O. J. C. Dias, G. T. Horowitz and J. E. Santos, Gravitational turbulent instability of anti-de Sitter space, Class. Quant. Grav. 29(2012)194002.
  * (17) O. J. C. Dias, G. T. Horowitz, D. Marolf and J. E. Santos, On the nonlinear stability of asymptotically anti-de Sitter solutions, Class. Quant. Grav. 29(2012)235019.
  * (18) M. Maliborski and A. Rostworowski, Time-periodic solutions in an Einstein AdS-massless-scalar-field system, Phys. Rev. Lett. 111(2013)051102.
  * (19) A. Buchel, S.L. Liebling and L. Lehner, Boson stars in AdS spacetime, Phys. Rev. D 87(2013)123006.
  * (20) J. Abajo-Arrastia, E. da Silva, E. Lopez, J. Mas and A. Serantes, Holographic relaxation of finite size isolated quantum systems, J. High Energ. Phys. 05(2014)126.
  * (21) M. Maliborski and A. Rostworowski, What drives AdS unstable?, Phys. Rev. D 89(2014)124006.
  * (22) V. Balasubramanian, A. Buchel, S.R. Green, L. Lehner and S.L. Liebling, Holographic thermalization, stability of AdS and the Fermi-Pasta-Ulam-Tsingou paradox, Phys. Rev. Lett. 113(2014)071601.
  * (23) L D. Landau and E. M. Lifshitz, Statistical Physics (I) [M], 3rd ed., (Oxford: Pergamon, 1980). p. 146.