\section{Introduction}

The idea of self-organized criticality (SOC), introduced first time by Bak-Tang-Weinberg (BTW), is one of the most major concepts for understanding the physics of complex systems \cite{Bak}. A system displaying SOC organizes itself into a critical state with the statistical properties governed by power laws \cite{Tang, wat}. Unlike the standard critical systems, approaching the critical point occurs without tuning a control parameter. Avalanches with power law distribution emerge as a result of external driving and dissipation.

The BTW sandpile model is the simplest well-studied model which displays SOC behavior \cite{Bak, Dhar, Dhar2}. The model is simply defined on a square lattice. To each site, an integer height variable is assigned which represents the number of the sand grains of that site. A configuration in which all sites have heights less than a specific threshold is stable. The dynamics of the model is started from a stable configuration. At each step a grain of sand is added to a randomly chosen site so that its height is increased by $1$. If the height of that site reaches its threshold value, that node becomes unstable, topples and loses some grains of sand, redistributed among its nearest neighbours. As a result of the toppling, some of the neighbours may also become unstable. The toppling process continues until all sites become stable and an avalanche is formed. Avalanches that reach the boundaries allow sand to exit the boundaries. The balance between the external driving and dissipation tunes the system to the stationary state with  critical properties.

In a stochastic version of the BTW model, the so-called Manna model, the threshold value is set to $2$ \cite{manna}. If the height of a site reaches $2$, it topples and loses two grains. In this case the sand grains are not redistributed uniformly among the neighbors. Two neighbors are chosen randomly and each one receives a grain. The critical behavior of the Manna model and its universality class compared with the BTW model have been studied using many numerical methods \cite{Lubeck, Chessa, univ1, univ2, moghimi}.
\begin{figure*}[t]
\begin{center}
\scalebox{0.63}{\includegraphics[angle=0]{fig1n.eps}}
\end{center}
\caption{Time evolution of the density of active nodes, $\rho(t)$, for (a) $\epsilon < \epsilon_{c1}$, (b), (c) $\epsilon_{c1} <\epsilon < \epsilon_{c2}$, (d) $\epsilon_{c2} <\epsilon < \epsilon_{c3}$, (e) $\epsilon_{c3} <\epsilon < \epsilon_{c4}$ and (f) $\epsilon > \epsilon_{c4}$, on an ER network with $N=10^4$ and $\langle k\rangle=10$.}
 \label{f1}
\end{figure*}

To explain how the sandpile model self-organizes into a critical state, Dickman et. al. related the standard sandpile model to a conservative system, the so-called fixed energy sandpile (FES) model \cite{Dikman, FES, FES2}. In this model driving and dissipation are set to zero, such that the total number of sands (which can be interpreted as energy) is constant. A site at which the number of grains is equal to or greater than the threshold value is called active. The dynamics of the model is started with a few number of initial active sites. At each time step, the active sites topple and redistribute the sands to their neighbors. The process reaches either an absorbing state, in which no active site exists, or an active state with a constant fraction of active sites. Hence, the fixed energy sandpile model exhibits a second order phase transition, as the energy density grows, from an absorbing to an active phase \cite{PT1, PT2, montakhab}.
The SOC tunes the system to the edge of this transition through slow driving and dissipation. Hence, the scaling properties of the SOC system can be found from the critical properties of the absorbing phase transition in the FES model \cite{FES}. By extending the SOC to systems exhibiting bistability and phase coexistence, the concept of self-organized bistability (SOB) has been recently proposed as a mechanism through which a system self-organizes to the edge of bistability of a discontinuous phase transition \cite{munoz1}.

In general, the BTW sandpile dynamics can be defined on a random network with any arbitrary degree distribution \cite{Bonabeau, Goh1}. In this case, the threshold value for each node is different and is assumed to be equal to the degree of each node. The dissipation of sand grains occurs either through some sinks from which grains flow out of the system, or by assuming that at each time step a fraction of the grains is lost during the toppling. It is found that avalanche power-law distributions and critical exponents depend on the degree distribution and structure properties of the underlying network \cite{Goh2}.

As the structure of a network changes the avalanches distribution of the sandpile model, the sandpile dynamics can affect the network structure. The coevolution of the network structure and dynamics of the processes taking place in the network is studied in the framework of adaptive networks \cite{adaptive1, adaptive2}. In Ref \cite{Bianconi}, Bianconi et. al. proposed a simple model of network dynamics, which has the same dynamics as the sandpile model. The dynamics occurs in two parts: slow driving by adding new edges to the network and fast relaxation dynamics of avalanche of rewiring processes with dissipation. They showed that the network consequently self-organizes to a critical point in which avalanches present scaling behaviors. Also, in Ref. \cite{coevolution}, it was shown that how such an interplay between the sandpile dynamics and network structure leads to self-organization to a critical point in which the distribution of avalanches and degree distribution of the network similarly behave as power-law functions.

In this paper we consider the FES model with the BTW toppling rule (BTW-FES) on random networks. It is interesting to see the effect of the small-world property on the dynamics of this model. Using numerical simulations, we show that the active sites topple in synchrony in the stationary state if the density of sand grains exceeds a critical value.
We find the bifurcation points and discuss Hopf bifurcations which emerge when energy density of the system changes. To see the effect of symmetric toppling rules, we show that the oscillations disappear when we consider the stochastic Manna-FES model on random networks. At the following we use the interplay of structural effects and the BTW-FES dynamics which destroys the ordered pattern of topplings among the nodes. We rewire the network edges due to the state of nodes and show that the coevolution of structure and the BTW-FES dynamics makes the sustained oscillations to disappear. Also, as the result of rewiring, a coexistence region between the active and absorbing states emerge.

The article is organized as follows. In the next section, we define the dynamics of the BTW-FES on a random network and obtain the phase diagram and bifurcation points of the model. In Sec.~III we discuss the origin of emergence of oscillations in the dynamics of sandpile models.
The interplay between the network structure and the sandpile dynamics is discussed in Sec. IV. The paper is concluded in Sec. V.

\section{The model}

Let us consider an Erd\H{o}s--R\'enyi (ER) random network with $N$ nodes and mean degree $\langle k \rangle$, a paradigmatic example of a homogeneous network. We define the BTW-FES on this network as follows. A conserved quantity $E$ is considered as the total number of sand grains on the network. The quantity $E$ can be also interpreted as the total energy, such that the parameter $\epsilon=E/N$ is the energy density of the system. Each node has a threshold $h^c_{i}$, which is the capacity of the sand or energy that the node can hold. We set the threshold values equal to the degree of each node such that $h^c_{i}=k_i$. At the initial state, the total energy $E$ is divided randomly among $N$ nodes such that each node gets $h_i$ sands or energy. The nodes at which $h_i \geq h^c_i$ are called active, otherwise they are inactive. In other words at the initial state some nodes are active and the rest are inactive. Hence we can consider the model as a binary activation model with states $s_i=0$ for the inactive nodes and $s_i=1$ for the active nodes, such that the density of active nodes is defined as $\rho(t)= \frac{1}{N}\sum_i s_i(t)$ at each time step $t$. The activation dynamics occurs as follows. At each time step, each of the active nodes topples once following the same rule as the standard BTW model, i.e., $h_i= h_i-k_i$ and each of its neighbors gets one sand such that $h_j=h_j+1$. Active nodes can still have sands more than the threshold and therefore remain active. Also inactive nodes can become active as a result of toppings. All active nodes will topple at the next time step and the process continues until the system reaches a stationary state.

We carried out Monte Carlo simulations of the model on an ER network with $N=10^4$ and $\langle k \rangle=10$. Figure~\ref{f1} shows the evolution of the density of active nodes. The stationary state depends on the value of energy density $\epsilon$. For $\epsilon$ less than a critical value $\epsilon_{c1}$, the model shows an ``absorbing" phase, in which $\rho_{stat}=0$ (Fig.~\ref{f1} $(a)$). That is the same as what happens on a lattice. However for $\epsilon > \epsilon_{c1}$, the density of active nodes oscillates in the stationary state.
Close to the critical point $\epsilon_{c1}$, the oscillations are regular but by increasing $\epsilon$ the irregular oscillations appear (Fig.~\ref{f1} $(b),(c)$). A super-critical Hopf bifurcation takes place at $\epsilon_{c1}$
with the appearance of the limit cycle for the density of active nodes and the ``oscillatory" phase is observed in the range $\epsilon_{c1}< \epsilon < \epsilon_{c2}$. At the transition point $\epsilon_{c2}$, we observe that the active nodes start to oscillate between two states, up and down (Fig.~\ref{f1} $(d)$). We call this regime, which is in the range $\epsilon_{c2}< \epsilon < \epsilon_{c3}$, the ``binary" state. For $\epsilon > \epsilon_{c3}$ the density of active nodes increases to more than half while the density of inactive nodes decreases. Hence the model shows a symmetrical behavior such that the behavior of the inactive and active node densities are exchanged. Therefore for $\epsilon_{c3}< \epsilon < \epsilon_{c4}$, the oscillations appear again (Fig.~\ref{f1} $(e)$). With further increase $\epsilon$, the density of active nodes goes up while the amplitude of the oscillations decreases. Finally, at the transition point $\epsilon_{c4}$, the oscillations disappear and the density of active nodes starts approaching the trivial fixed point $\rho_{stat}=1$, in which each node is active. Hence, a sub-critical Hopf bifurcation occurs at $\epsilon_{c4}$ . We call the region $\epsilon > \epsilon_{c4}$ the ``active" phase (Fig.~\ref{f1} $(f)$). Figure~\ref{f2} shows the density of the active nodes as a function of $\epsilon$. The transition points depend on the initial state. For a uniformly random initial distribution of the sands, the transition points occur at $\epsilon_{c1}=5.51\pm0.01$, $\epsilon_{c2}=8.95\pm0.01$ , $\epsilon_{c3}=10.05\pm0.01$ and $\epsilon_{c4}=13.96\pm0.01$. We can also see the different phases of the model in Fig.~\ref{f3}. As we can see the transition points also depend on the mean degree of the network and shift to higher values for the networks with more connectivity.

\begin{figure}[t]
\scalebox{0.45}{\includegraphics[angle=0]{fig2n.eps}}
\caption{Stationary values of the density of active nodes as a function of $\epsilon$, on an ER network with $N=10^4$ and $\langle k \rangle=10$. Transition points occur at $\epsilon_{c1}=5.51$, $\epsilon_{c2}=8.95$ , $\epsilon_{c3}=10.05$ and $\epsilon_{c4}=13.96$.}
\label{f2}
\end{figure}
\begin{figure}[t]
\begin{center}
\scalebox{0.45}{\includegraphics[angle=0]{fig3n.eps}}
\end{center}
\caption{Phase diagram of the BTW-FES model on an ER network with $N=10^4$ in the space $\epsilon-\langle k \rangle $.}
\label{f3}
\end{figure}

\section{Emergence of rhythmic behavior}
The process of activation of nodes in the BTW-FES model on ER random networks occurs as a rhythmic behavior, which is different from the well-known behavior of the model on lattices. The root of the oscillations can be in two factors: 1)  the small-world property of random networks and 2) symmetry in the dynamical rule of the BTW-FES model. We will discuss these two causes separately.

\subsection{small-world effect}

Let us consider the BTW-FES on a small-world network, produced following the Watts-Stogatz algorithm \cite{watts}. We start with a ring having $N$ nodes and node degree $k=2m$, such that each node is connected to its $m$ nearest neighbors on each side. Moving clockwise, for every node we select randomly an edge that connects that node to one of its neighbors, and rewire it with probability $p$. We continue this process until each edge in the original ring has been considered once. With increasing $p$ randomness increases such that at $p=1$ all edges are rewired and the resulting network corresponds with the ER random network.

We consider the dynamics of the BTW-FES on a small-world network with $N=10^4$ and $m=5$. At the initial time the fraction $\epsilon =5.7$ of sands are randomly distributed among the nodes of the ring. Figure \ref{f4} shows the evolution of the density of active nodes for some values of $p$.  As we can see, for $p=0$ (a regular ring), the density of active nodes decreases with the time and reaches a stationary fixed point. However with increasing the value of $p$, oscillations emerge in the stationary state. For higher values of $p$, the amplitude of the oscillations increases and we see a periodic behavior.

The oscillations correspond to a synchronized behavior for the active nodes. At any step, say at $t=0$, any node with index $j$ is active. These nodes become inactive altogether at the next step. After a period time $\tau_j$ they will be active again and this process is repeated. We can see this periodic pattern for the activation of the nodes in Fig.~\ref{f5}. Active nodes are represented by blue (dark) points, and white points show the inactive nodes. Because of periodic behavior, we can assign a phase $\phi_j$ to each node $j$, defined as $\phi_j= 2\pi(\tau_j-1) / T$, in which $T$ is the largest value of $\tau_j$. Hence, the order parameter for synchronization of phases $\phi_j$ at time step $t$ is defined as,
\begin{equation}
R(t) = \big|\frac{1}{N}\sum _j ^N e^{i \phi_j (t)}\big|.
\label{eq1}
\end{equation}

\begin{figure}[t]
\begin{center}
\scalebox{0.58}{\includegraphics[angle=0]{fig4n.eps}}
\end{center}
\caption{Time evolution of the density of active nodes for $\epsilon=5.7$ on small-world networks with different values of rewiring probability $p$.}
\label{f4}
\end{figure}
\begin{figure}[t]
\begin{center}
\scalebox{0.4}{\includegraphics[angle=0]{fig5n.eps}}
\end{center}
\caption{Raster plot for the activation of nodes in a small-world network with $m=5$ and $p=1$. The density of sand is $\epsilon=\epsilon_{c1}=5.51$.}
\label{f5}
\end{figure}

Close to the critical point $\epsilon_{c1}$, a large number of nodes are synchronized and are activated together, such that the value of $R$ is large. The synchronization is complete, $R=1$, if all nodes become active at the same time. However as we move away from the critical point, the nodes start to be  activated irregularly, such that the corresponding phases become different and are distributed between $[0, 2\pi]$. In this case the order parameter $R$ gets close to zero.

Figure.~\ref{f6} shows the behavior of $R$ in terms of the parameter $p$. To see a better synchronized behavior we consider a larger network with $N=10^6$ nodes. Since close to the critical point $\epsilon_{c1}$, the activated nodes are more synchronized, we consider $\epsilon \sim \epsilon_{c1}$ for each network with a definite value of the mean degree. We can see a phase transition occurs from a non-synchronized to a synchronized state at a definite value $p_c$. The value of $p_c$ depends on the number of neighbors $2m$. For a network with larger mean degree, the synchronization emerges at a smaller value of $p_c$.

\begin{figure}[t]
\begin{center}
\scalebox{0.4}{\includegraphics[angle=0]{fig6n.eps}}
\end{center}
\caption{Synchronization order parameter $R$ as a function of $p$ on small-world networks with $N=10^6$. The curves are plotted at $\epsilon_{c1}$ for different mean-degree values $\langle k\rangle = 2m$.}
\label{f6}
\end{figure}


\subsection{symmetric dynamical rule}

In the BTW-FES model, the redistribution rule of sand grains is symmetric. It means that when a node topples it gives exactly one sand to each of its neighbors. This symmetric rule makes an order among the nodes such that each node reaches its threshold at the same time as a number of other nodes. If the symmetric redistribution of sand grains is broken, there is no synchrony among the nodes. To see that, we consider the Manna-FES model with stochastic and asymmetrical dynamics. The model is defined on the ER network with $N$ nodes and mean degree $\langle k \rangle$. At first, $\epsilon N$ grains are randomly distributed among the nodes such that each node $i$ receives $h_i$ grains. It is assumed that the threshold value is the same for all nodes, i.e., $h_{ic}=2$. Hence, at each time $t$, the nodes with $h_i \geq 2$ are active and topple. When each node topples it loses two grains and gives one grain to two of its randomly selected neighbors. The dynamics continues until the system approaches a stationary state. The Manna-FES model also undergoes an absorbing phase transition at a definite value $\epsilon_c$. For uniformly random initial distribution of the sand grains on the ER network with $\langle k\rangle =4$, the critical sand density value is $\epsilon_c= 0.74\pm 0.01$. Figure.~\ref{f7} shows the evolution of the density of active nodes $\rho$ for two values of $\epsilon$ in the absorbing and active states. As we can see for $\epsilon > \epsilon_c$, the density of active nodes decreases and reaches a nonzero stationary fixed point. In other words, the FES models with stochastic rules do not show oscillations. Also for $\epsilon < \epsilon_c$, the same behavior is observed with the difference that $\rho$ approaches zero in the stationary state.

\begin{figure}[t]
\begin{center}
\scalebox{0.4}{\includegraphics[angle=0]{fig7an.eps}}
\scalebox{0.4}{\includegraphics[angle=0]{fig7bn.eps}}
\end{center}
\caption{Time evolution of the density of active nodes for the Manna-FES model on an ER network with $N=10^4$ and $\langle k\rangle = 4$. In $(a)$ $\epsilon = 0.5 < \epsilon_{c}$ and in $(b)$ $\epsilon = 0.9 > \epsilon_{c}$. }
\label{f7}
\end{figure}

\section{Activity-dependent rewiring}

So far, we have considered the dynamics of the fixed energy sandpile models on static networks. In this case the degree of nodes and thus the value of thresholds are constant. Here we allow the network structure changes during the dynamics. We consider a simple rewiring rule such that during the dynamics of the BTW-FES model, active nodes avoid to connect with each other. In this case, the toppling of the active nodes makes it possible for more inactive nodes to be activated. Also the threshold values change during the rewiring and then the ordered pattern for the toppling of the nodes breaks down.

Let us consider the BTW-FES model on the ER random network with $N=10^4$ and $\langle k\rangle=10$. We take $\epsilon=\epsilon_{c1}=5.51$. The sands are randomly distributed among the nodes.
Hence, the initial state of the nodes is $(s_1(0),\ldots s_N(0))$, such that some nodes are active $(s_i(0)=1)$ and the rest are inactive $(s_i(0)=0)$ initially. At each time step, we first apply the following rewiring rule for $f$ fraction of the active edges (edges which connect two active nodes). An active edge is randomly chosen and from one of its active end nodes, it is disconnected and rewired to a randomly chosen inactive node in the network. Then we apply the dynamics of the BTW-FES model and allow all active nodes to topple. At the next time step, we repeat the rewiring rule and synchronously update the state of the nodes. The process is continued until the system reaches a stationary sate.

Figure~\ref{f8} shows the evolution of the active node density during the rewiring process. As we can see, during rewiring the density of active nodes increases. Whenever an active node topples, it loses as much sand as its degree. However, the node degree and therefore the activation threshold of the nodes change due to rewiring such that the ordered activation pattern of the nodes is broken. That is why oscillations become smaller and after a while the dynamics evolves to a stationary fixed point. At this point, we stop the rewiring and allow the dynamics of the FES model to proceed. After stopping the rewiring, we still have a random network that has the same mean degree as the network before the rewiring. Hence, as we expect after a while the oscillations emerge again due to the small-world effect and the order of topplings among the nodes. However, the rewiring causes that the arrangement of the active and inactive nodes changes and therefore the density of the active nodes oscillates around a higher mean value.
\begin{figure}[t]
\begin{center
\scalebox{0.45}{\includegraphics[angle=0]{fig8n.eps}}
\end{center}
\caption{Time evolution of the density of the active nodes for the BTW-FES model on an adaptive network. The initial network is an ER random network with $N=10^4$ and $\langle k\rangle=10$. At each time step edges are rewired with probability $f=0.1$. The dashed line indicates the time at which the rewiring stops.}
\label{f8}
\end{figure}
\begin{figure}[t]
\begin{center}
\scalebox{0.42}{\includegraphics[angle=0]{fig9n.eps}}
\end{center}
\caption{Phase transitions of the BTW-FES model on an adaptive network. The initial network is an ER random network with $N=10^4$ and $\langle k\rangle=10$.}
\label{f9}
\end{figure}

Fig.~\ref{f9} shows the steady state density of the active nodes $\rho$
as a function of the density of sand grains $\epsilon$. As we already saw, in the case of static networks a continuous absorbing phase transition occurs. However by including rewiring, a bistable region is observed in which the active and absorbing phases coexist. The coevolution of structure and dynamics makes the system abruptly jumps from an absorbing phase to an active phase and the absorbing phase transition shifts to a higher value of $\epsilon$. On the other hand with decreasing of sand density, the system tends to remain in the active phase and finally collapses in the absorbing phase in a lower value of $\epsilon$. This observation can be related to the concept of self-organized-bistability, such that the system self-organizes to a bistable region instead of a critical point.



\section{Conclusion}

In this work, we studied the dynamics of the FES models using numerical simulations on random networks. The BTW-FES model with deterministic dynamics undergoes an absorbing phase transition from a fixed point with zero density of active nodes to an oscillatory active phase. We discussed that the small-world effect and symmetric redistribution of sand grains lead to the emergence of sustained oscillations. By increasing the sand density, the model undergoes two other subsequent phase transitions and binary phase and active fixed point emerge. We obtained the phase diagram of the model and found the transition points for random networks with different mean degrees. The values of transition points depend on the mean degree of the network and the initial distribution of sand (energy).

If the toppling rule is not symmetric, as it is in the Manna-FES model, there is not an ordered toppling pattern among the nodes. In this case, the oscillations do not emerge and similar to what happens on lattices, the transition occurs from an absorbing state to an active fixed point. The ordered toppling pattern can be also destroyed by rewiring the network edges during the BTW-FES dynamics. In the rewiring process, we assumed that the active nodes avoid to connect with each other. In this case, the density of active nodes grows up and the oscillations will gradually disappear during the rewiring. Moreover, the coevolution of the BTW-FES dynamics and the structure of the network leads to the coexistence of absorbing and active phases and bistability emerges.

The sustained oscillations also appear in some cyclically three-state systems such as the susceptible-infected-recovered-susceptible (SIRS) epidemic model \cite{SIRS} and the rock-paper-scissor (RPS) game \cite{RPS}. For these systems, a phase transition occurs from an absorbing state to a self-sustained oscillation state if the network randomness exceeds a threshold value. Hence, the small-world effect results in qualitatively similar behaviors between the BTW-FES model with a symmetrical toppling rule and the cyclic three-state models. We propose the BTW-FES model has an inherent cyclic three-state dynamics that results in periodic oscillations. However finding an exact equivalence of the FES models and the cyclic three-state models is a challenge for future work.


